Selection of stimulus parameters for enhancing slow wave sleep events with a neural-field theory thalamocortical 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 with 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 random 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.


Introduction
There is evidence supporting that sensory manipulation of SO leads to memory improvement [3] or disruption [20], with significant effects in enhancing the declarative memory and no impact on procedural memory [20,21]. Moreover, the number of SPs during SWS, and their power, is also related with best results in word-pairs memory tests [2]. However, the knowledge and optimization of stimulation protocols have been based on correlation analysis rather than causal evidence [22].
Even if auditory closed-loop stimulation has been used to improve the memory consolidation process, there is no evidence of tuning the stimulation protocol to maximize this effect. Specifically, sensory stimulation looks to maximize the correct timing between the evoked SP and SO events associated with the memory consolidation process. The manipulation of SP in auditory closed-loop [23], and the use of memory cues (sounds related to the learning episode), show better results when the external stimulus is applied at a specific phase of the SO [24].
An alternative approach to improving the stimulation paradigm that attempts to increase SP and SO events during SWS is the use of mathematical models of brain activity [25,26]. In this way, mechanistic insights about the generation of SP and SO activity can guide the search for a better stimulation that maximizes its efficacy. In this work, we present a systematic study of possible stimulation protocols and the more efficient features of the stimuli in promoting these SWS events' co-occurrence.
Specifically, in this work we use a model based on the Neural Field Theory (NFT) [27,28] with corticothalamic loops between two cortical neural populations: excitatory or pyramidal (e) and inhibitory or interneurons(i), and two thalamic populations: reticular nucleus(r) and relay nuclei(s). Previously, this model has been tuned to reproduce the electrical brain activity in different states: from arousal to slow-wave sleep and including the different stages of NREM sleep [29,30]. In this NFT model, sensory stimulation can be introduced as spatial-localized perturbations on the noisy input neural population, obtaining an EEG-like signal as the model's output. From this signal, we can assess the efficacy of the stimulation paradigms and implement a closed-loop feedback mechanism that can fine-tune the stimuli delivery timing.
With the purpose of simulating a single-modality sensory stimulation, we searched for a one-channel impulse pattern that alters the dynamics of the entire cortex. We assessed the efficacy of different pulse shapes, delivered with different frequency, intensity, and periodicity, including a closed-loop algorithm that intended to phase-lock the stimuli with the SO. The stimulation goals were: 1. To increase the power at both frequency bands of interest (SO and SP), measured by the scalogram (wavelet spectrogram).
2. To increase the quantity of SO and SP events and detecting their occurrence on the time series.
3. To rise the co-occurrences between SO and SP events, over the level of co-occurrences obtained by chance. 4. To increment the occurrence of SP on the up-peak of the SO, measured by the delay time between the down-peak of the slow oscillation and the center of the coincident sleep spindle.
Our results propose a sensory stimulation protocol, maximizing the SP and SO co-occurrence associated with the memory consolidation process, for its use and verification in experimental setups of memory consolidation enhancement. Specifically, the SP and SO events' best co-occurrence is obtained when 50 ms decreasing ramp pulses were delivered at the zero phases (sine reference) of the ongoing activity filtered at a specific central frequency on the SO band. Nevertheless, the behavioral effect of this stimulation protocol on memory performance should be evaluated with clinical experiments.

Simulation of EEG recordings of NREM sleep
In EEG recordings, the classification as NREM sleep comes from the presence of SWS events on epochs of 30 seconds. The prominence of SO indicates the N3 stage, and the occurrence of SPs indicates the N2 stage [31]. Furthermore, the evidence suggests that SP nested in the up phase of the SO may play a role in memory consolidation [15]. A SWS model, including SOs and SPs, presents advantages to characterize the conditions promoting their co-occurrence. We used the neural-field theory model with thalamic and cortical populations, eirs-NFT [27], with the connections shown in Fig 1A to obtain simulated EEG-like recordings. The spatiotemporal propagation of the firing rate in each population (e,i,r,s) is used as input of the synaptic-dendritic dynamics of the soma voltage in other populations ( Fig 1B). The connections strengths and the steady-state of the input population � ð0Þ n determines the behavior of the 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' spatial distribution. The location of stimuli is highlighted by the circles in color, representing the standard deviations of the Difference of Gaussians spatial filter.
https://doi.org/10.1371/journal.pcbi.1008758.g001 model (see Eq (7) in the Materials and methods section) by changing the value of the connection gains G ab , being a the receiving population and b the delivering one. Our model includes a stage for sensory stimuli input (left block of Fig 1C), consisting on Gaussian white noise along all the spatial extension, v(r,t), to which the stimuli signal u 0 (r,t) is added when applied. The stimuli signal location comes from applying a Difference of Gaussians filter to the stimuli signal u(t). A spatial filter acts as a receptive field, with excited and inhibited areas, on the relay nuclei of the thalamus. Finally, the output is the spatio-temporal activity of the excitatory population ϕ e (r, t), which is averaged in the spatial domain to generate a simulated EEG-like signal (right block of Fig 1C).

Selection of model parameters for SWS stage
Previous works using the eirs-NFT model were specific for the N2 stage (only SP) [29], or for the N3 stage only (without SP) [30]. Then, our first goal will be to tune the model parameters into a regime that exhibits both events.
In the eirs-NFT model, the sleep (and awake) stages can be represented as points in a XYZ space of model parameters, related to physiological connections. Fig 2A shows this space where the X-axis represents the strength of the intra-cortical connections, higher as sleep deepens. The Y-axis represents the internal thalamocortical connections that have a negative value during sleep. The more pronounced variations between sleep stages are in the Z-axis, representing the intra-thalamic connections [30]. In terms of the loop gains, the axes are defined as: where α and β are the rising and decaying rate of the synaptic-dendritic potential, respectively. Then, the position in the XYZ space depends on the model parameters and � ð0Þ n . Furthermore, Fig 2A shows excerpts of the cortical excitatory population activity for N2 and N3 stages, and we traced a trajectory between these two stages using linear interpolation on the connection strengths. On top of the cortical activity, the root mean square (RMS) value of the bands shows the higher presence of low-frequency SO activity (orange line) when the point is at a low Z-value. The SP activity (green dashed line) increases as the point moves towards more negative values on the Y-axis. Dominant spindle activity is appreciable using the N2-stage parameters indicated in [29]. The spindle activity decreases as the simulation point slides along the interpolated trajectory towards the N3 stage, which shows slow-wave sleep dynamics but not noticeable spindle activity.
We are interested in an intermediate state between N2 and N3. For this, we concentrated on having energy both in N2 and N3 frequency bands. We obtained this condition in the second half of the linear interpolation between N2 and N3. Within this segment, we arbitrarily selected the point at 2/3 on the N2 to N3 trajectory. Fig 2B shows the linearized model spectrums of N2 and N3 stages and three extra points in the interpolation trajectory, including the 2/3 position in the blue-line.
The increase of the steady-state input value (� ð0Þ n ) raises the model response in the Z-axis (Fig 2A), indicating lightness of NREM sleep [30,32]. Besides, it also increases the power and frequency of the SP-band peak (Fig 2C), then a higher steady-state input value facilitates SP occurrence.
We selected � ð0Þ n ¼ 1 s -1 as steady-state input in order to have similar peak amplitudes for frequency band lobes of both SWS events. The black line in Fig 2C shows the corresponding spectrum of the chosen simulation parameters. The peak frequency is close to the lower SP frequency with the selected parameters for the connections gains and the steady-state input (see Materials and methods).
The variation of the input, expected with the addition of stimuli signal, has a high effect in the relay nuclei activity � ð0Þ s , and a less considerable effect in the reticular nucleus activity � ð0Þ r , and the cortical activity � ð0Þ e . Fig 2D displays the relative changes of the loop-gains and the steady-state of the neural populations for the changes of the steady-state input activity.

Enhancement of SWS events on the eirs-NFT model
We searched for changes in the occurrence of the slow-wave sleep events of interest (SO and SP) when we added stimulation u(t) to the noisy input propagation field ϕ n (r, t). Slow oscillations (SO) were defined to occur every time the signal, representing the cortical activity at the band 0.5-1.25 Hz, had a negative peak below -40 μV, and a peak-to-peak amplitude higher than 75 μV. We selected the higher amplitude oscillations at the SO band (see Materials and methods). On the other hand, sleep spindles (SP) were defined as bursts of oscillations in the 9-16 Hz band. Fig 3A shows the simulated cortical activity for three different stimulation types: random (STIM-R), periodic (STIM-P), and closed-loop with a target-phase (STIM-CL). The top panel shows the simulated cortical activity for the SHAM condition, b(t). STIM-R stimulation delivers stimuli pulses at random times from a Poisson stochastic process with a mean inter-stimulus time. With STIM-P, the pulses occur at a fixed inter-stimulus period. STIM-CL follows the ongoing activity at a narrow band delivering a stimulus every time the signal arrives at a target-phase (see below, section SWS events activity variation with closed-loop stimulation). Fig  3B shows the correspondent scalograms for each time series in Fig 3A, where δ-band activity (0.5-4 Hz) -that includes SO-band (0.5-1.25 Hz)-is present through all the time. On the other hand, the SP-band activity (9-16 Hz) occurs at discrete periods, some of them with enough duration and amplitude to be classified as spindles.
Fig 3C focuses on a smaller time window when two SPs and one SO occur. The vertical marks in Fig 3A and 3C show the start time point of each stimulus pulse that, in this case, consists of decreasing ramps of 0.1 s duration and 40 a.u. of energy. The horizontal marks in Fig  3A and 3C represent detected SWS events. Their duration and peak-to-peak amplitude have the distributions shown in Fig 3E and 3F, respectively.
Using the time-averaged scalograms, we calculated the spectrums shown in Fig 3D. We defined the power difference index, I (SO, SP) , as the area between the corresponding STIM condition and SHAM spectrums inside the frequency band of interest, normalized by the sum of the power of both cases. Then, I (SO, SP) quantifies the differences with respect to the SHAM case for the frequency bands of interest (see Materials and methods).

Selection of the stimulus pulse characteristics with open-loop stimulation
We characterized the efficiency of different stimulation signals to enhance the occurrence of slow-wave events. In a first step, we tested stimulus pulse features with the open-loop stimulation (STIM-R and STIM-P). We looked at power differences at the frequency bands of interest and the number of occurrences of SWS events while changing the shape, the energy, and the duration of the stimulus pulse. In the following, we will evaluate each of these parameters' impact on the co-occurrence of SP and SOs.
Changes on SWS events by pulse shape. With STIM-P and STIM-R, the shape, duration, and pulse energy change the temporal and frequency domain characteristics of u(t). The stimulation frequency also modifies them in the case of STIM-P. We tested for the six shapes in Table 1, and shown in Fig 4A and 4B, with the same duration (0.1 seconds), keeping constant the stimulus pulse energy (E p = 40 a. u.) and the stimulation frequency (0.2 pulses per second).  The co-occurrence of SP and SO was evaluated using: (i) the conditional probability of counted co-occurrences given the number of detected spindles, P(C|SP); and (ii) the probability of slow oscillations, P(SO). In Fig 4E, a positive result means that P(C|SP) is higher than P(SO); in other words, that the probability of spindles co-occurring with SOs is higher than chance. Given this, both ramp shapes (decreasing and increasing) have positive results for the co-occurrence of events.
Changes on SWS events by duration and energy of stimulus pulse. Other stimulus pulse characteristics could modify the SWS events measurements. The pulse duration modifies the influence of stimulation on the frequency domain (narrower stimulus, wider spectral  The minimum pulse energy was selected to overpass the input-noise power of 9.67 s -2 per second in each node. I (SO) and I (SP) increase as the pulse energy rises with statistical differences at each increase step (p < 0.01), but N SO /min. and N SP /min. detected by stimulation with pulse energy of 100 (a. u.) are not different from the detected events with pulse energy of 80 (a. u.) (p > 0.08). These measurements are also not different between the minimum tested pulse energy and the SHAM condition (p > 0.1). Interestingly, Fig 5D shows a peak of N SO /min. at 60 (a. u.), but without statistical difference from their neighboring pulse energies. The decrease of N SO /min. with higher pulse energies arises from the impossibility of evoking more SO events without overlapping them. Regarding the co-occurrence of events, the rectangular pulse shape of 0.05 s, and the decreasing ramp pulses of 0.1 s with energies above 40 (a. u.) show positive results in Fig 5C and 5F, respectively.
The overall results of Fig 5 suggest that the best choice is the 0.05 s duration pulse, always considering that lower duration translates to higher amplitude in order to keep the same pulse energy, which in turn must be enough to overpass the intrinsic activity. Nevertheless, optimal results are found without need to rise over a six-fold value of the intrinsic power per second.
When exploring the stimulation frequency using a decreasing ramp of 0.1 seconds, we found no statistical differences between frequencies in the 0.5-1.25 Hz range (see S2 Fig). Stimulation with pulses at random intervals (STIM-R) has very similar effects of periodic pulses (STIM-P), except for a higher N SO /min. with STIM-P than STIM-R (see the first row of Fig 5 and the second row of S3 Fig). Nevertheless, the results keep the same order between the characteristics variation. Thus, the difference of effects in Fig 5 comes from the pulse characteristics and not from the stimulation periodicity.

SWS events activity variation with closed-loop stimulation
Next, we asked whether the phase of the ongoing activity could influence the power and the number of slow-wave sleep events at the moment of the stimulus onset. To answer this, we implemented a closed-loop stimulation protocol. Here, the stimulus pulse is applied when the inverse-notch filtered signal arrives at a particular target phase (see Fig 6A and Materials and methods section). Our phase detector uses a fixed frequency of the ongoing activity, and 0.85 Hz has been used in other closed-loop implementation as the central frequency of a phaselocked loop (PLL) stage [33,34].
We applied the decreasing ramp stimulation pulses of 0.1 s duration and E p = 40 (a. u.) in the three stimulation cases: STIM-P, STIM-R, and STIM-CL. STIM-CL was further divided into STIM-CL 0 (stimulus onset when the phase of slow oscillations arrives at 0 degrees; the upside cross-zero point), STIM-CL 45 (45 degrees) and STIM-CL 90 (90 degrees; up-peak of the slow oscillation).
The STIM-CL 0 case outperforms the other two closed-loop stimulations in all the measurements as shown in Fig 6B-6D, with statistical significance in I (SP) and P(SO) (p < 0.01 Welch's independent t-test if both cases pass Shapiro, otherwise Wilcoxon test, see S1 and S2 Tables for all t-values and p-values). The STIM-CL 0 case is also statistically different from STIM-CL 90 for I (SO) .
In relation to the open-loop stimulation cases, STIM-CL 0 outperforms STIM-P and STIM-R in I (SO) (p < 0.01 vs STIM-P), I (SP) (p < 0.01 vs STIM-R) and N SP /min. Interestingly, STIM-P outperforms all other stimulations in N SO /min. but is only statistically different with STIM-R. These results probably come from the zero chance in STIM-P of two pulses occurring with an inter-pulse interval less than the stimulation frequency period.
The co-occurrence of events shows positive results in Fig 6D for all stimulation cases. STIM-P and STIM-CL 0 are very similar in the P(SO) versus P(C|SP) plane, but this plot does not say anything about the timing of occurrence of these temporal coincidences.

Timing of SWS events occurrence
A higher occurrence, or higher power of slow-wave events, may not be enough to enhance memory consolidation [35]. The relative timing of occurrence between SPs and SOs also play a role, with SPs being expected to be nested in the up-phase of a slow oscillation.  Table 2). SPs' timing distribution from STIM-CL 45 is also significantly different from the SHAM distribution (p < 0.01). Other significant differences (p < 0.05) are STIM-R vs. SHAM, and STIM-R vs. STIM-CL 0. The STIM-CL 45 and STIM-CL 90 show better SPs' timing, but the co-occurred events and the other measurements are higher with STIM-CL 0. Then, STIM-CL 0 is the best among the tested cases of stimulation. Further, inter-stimuli interval distribution is different between the simulation cases (see S5 Fig). A stimulus modifies the ongoing activity. The activity's slope increases, even when the pulse is delivered nearer to the up-peak. Then, STIM-CL 45 and STIM-CL 90 delay the peak occurrence, needing more time to detect the same target phase again. Moreover, the

Discussion
In this study, we used a neural field model with thalamocortical populations to represent the SWS activity. We then explored the parameter space to maximize the occurrence of the SWS events (SOs and SPs), maximizing their frequency power and co-occurrence. In the openloop, between different pulse shapes, we propose the decreasing ramp pulse, with a pulse duration of 50ms, like the one maximizing the SWS features. On the other hand, the closed-loop stimulation, with stimulus applied at zero degrees of the ongoing activity (filtered at a specific frequency), increases the power and the number of both SWS events, together with their cooccurrences. Here, we analyze in detail the advantages and limitations of our contribution.

Selection of stimulation parameters
This research asks which stimulation pattern is effective in raising the power of slow-wave sleep events, consequently increasing the SP and SO occurrences. Moreover, we were also interested in maximizing the co-occurrence of SP with the up-peak of SO.
Our most significant result is that a pulse of 0.05 seconds with decreasing ramp shape maximizes the power and number of both SWS events. The features of the stimulation pulses delivered with STIM-P that result in higher I (SO) , I (SP) , N SO /min., and N SP /min. are also above the diagonal of chance in the plane P(C|SP) vs P(SO). It is interesting to note that in experimental acoustic stimulation studies, the pulse duration was 50 ms [3,33,36], a right choice following our results for pulse duration.
The pulse energy must be higher than the intrinsic activity. Following our results, the preferred energy value is 60 (a. u.). However, the translation of this energy value to a physical measurement of sensory stimulation is difficult. Thus, the useful recommendation is to use an energy value around six times the intrinsic EEG activity power per second.
To achieve SPs' precise timing at the up-peak of SOs, STIM-CL 45 and STIM-CL 90 conditions significantly differ in the timing distribution compared with the SHAM condition. The evoked response's effect (see S1 Fig) shifts the signal's peak, mainly because the activity slopewhich was about to reach zero-increases up to the following maximum. On the other hand, the stimulus generates a burst activity at the thalamus that could last enough to cause cortical SPs. These responses on both brain regions could justify the change of the timing distribution. However, the STIM-CL 0 condition provides a better result in power and number of events than STIM-P except in N SO /min. as shown in Fig 6B, and with the higher number of co-occurrences indicated at Table 2. Morevoer, the STIM-CL 0 condition actually shortens the inter-stimuli interval.
The pulse shape and its onset amplitude reshape the SWS events response. The shape of the stimulation pulses had small relevance in the previous neuro-stimulation works [37]. One probable cause is the technical difficulty of performing the stimulation following a particular shape-wave with precision in short duration times. For example, the triangular shape in Fig 4A and the Gaussian shape in Fig 4B may be indistinguishable from each other if the stimulation device has a small sampling rate. However, it is essential to consider that the mean firing rate from all neural populations does not represent the modulated amplitude of the sensory pulse linearly.
In any case, and regarless of the amplitude modulation, the pulse shapes in this work may represent variations in the firing response from different populations; each shape producing a different outcome from the thalamocortical system (see S4 Fig). Then, the rectangular shape relates to constant firing response. The rectangular trapezoid and the decreasing ramp are responses with short-term adaptation [38]. The Gaussian and triangular relate to a set of heterogeneous populations with different time-lag for the response's start. Finally, the rising ramp, less naturalistic, could be seen as an exaggeration of the sequential response of time-lag populations.
The presented shape of firing response includes the external sensory noise, but it excludes the variation by other brain internal sources, which is assumed the input noise ϕ n (r, t).
At the cellular level, the work of Pyragas et al. [39] demonstrates that using Pontryagin's maximum principle (a control theory formulation), a bipolar bang-off-bang waveform is optimal for the entrainment of a neuron to a specific frequency. This waveform type is not suitable for sensory stimulation, because a sensory stimulus affects excitatory and inhibitory neuronal populations at the same time [40]. With basis on the receptive fields of the visual perception system and with the assumption of the stimulus independence from the constant variance Gaussian noise input [41], we use a difference of Gaussians as a spatial filter of the simulated sensory input.
Different pulse shapes could generate benefits for stimulating during sleep, like soft slopes for avoiding abrupt sensory inputs of high prediction error and contrast-dependent temporal dynamics. A sudden stimulus is detrimental in sleep by evoking non-desired arousal. However, the results for inducing SPs and SOs were better when each pulse began with a high amplitude (see also S3 Table). The decreasing ramp shapes produce more beneficial and notable changes for both SWS events. In the same direction, small duration pulses with higher initial amplitude for keeping the same pulse energy show a higher increase in the power of both SWS events bands and a higher N SP /min. in Fig 5B. Furthermore, the rising ramp (time-mirrored version of the decreasing ramp) results in Fig 4D have the lowest performance of the tested shapes. These results suggest high relevance of the onset amplitude and the intrinsic phase of the stimulus pulse. Secondary effects of these stimuli, such as non-desired arousal, will have to be experimentally determined.

Results-driven by the model dynamic linear and non-linear effects
Mean-field models of the brain activity relay in two assumptions: The average population dynamics are similar to one neuron dynamics, and the average is enough statistic to characterize the population activity [27,42]. Further, the output activity of a mean-field neural population depends on the dynamics (see Eq (1) in Materials and methods) and statistics of the incoming signals (Eq 4). In addition to the u(t) dynamics, the input activity's statistics could play a role on the output dynamics of the model.
As in previous works of the Neural Field Theory [28][29][30]32], we use Gaussian white noise as input, which provides the linear impulse response of the model. The noisy input activity represents the net contribution from other sources as the basal ganglia, brainstem, and peripheral nervous system [30]. At the minute scale, variations from the average input value come from these region's activity rather than modifications to the connection's strengths. In our simulations, we used a permanent average value for the noisy input.
We superimpose a pulse train in this approach, modifying the input mean value, and its variance. Still, the steady-state results in Fig 2C do not predict all the dynamic effects produced by adding a pulse train at the model's input.
The most notable difference from the expected results with the steady-state approximated response appears when the tested stimulation parameter is the pulse energy. Fig 5D and 5E show that with the rise of pulses energy, there is an increase in the power of both SWS events bands. The increase of I (SO) is opposite to the expected result from Fig 2C when the steadystate input rises. This increase in power does not come only from the periodicity of stimulation. The random stimulation also increases the event's bands' power similarly (see S3 Fig). Moreover, in conjunction with results in S2 Fig, a specific stimulation frequency may not be too relevant for enhancing the power of both SWS events. It is also important to note that there are no significant differences in the results between STIM-R simulations.
Another possible way to change the input signal's statistics, and consequently its steadystate value, is by modifying their spatial characteristics. The coupling with sensory regions could matter to get results and predictions with sensory stimulation. However, the nodes' position in our model does not have a particular relation with a cortical region. Results from harmonics studies [43,44] make it possible to assign a physiologically-related topography in future works. Consequently, the results reported here are not related to a particular sensory type; they are general in how the increase of activity in a delimited region could modify the whole system's output.
More work is needed to consider non-linear effects that surpass the expected response using the linear approximation to explain the obtained results with stimulation. The use of modern control techniques [45], and the integration of information theory measurements [46] could provide a better analytical approximation of the results than traditional linear systems analysis.

Relation to similar computational works
The work of Schellenberger et al. [10] tries to reproduce in a neural-mass model the experimental results from Ngo et al. [3] with closed-loop stimulation. They recover the experimental results using two external calcium-dependent currents and an anomalous rectifier current responsible for the waxing and waning structure of the spindle oscillations. The authors employed the same neural populations of eirs-NFT with one additional connection, the autoconnection of the thalamic reticular nucleus, and two additional noise inputs to the cortical populations. Similar to our approach, the sensory stimulation (80 ms square pulses with amplitude of 70 s -1 ) is added to the noise input of the relay nuclei. The simulation results show high similarity with the SO band's experimental ones on amplitude and timing, but differences in the SP band activity's timing.
Wei et al. [26,47,48] tested several ways of closed-loop stimulation, including the protocol of Ngo et al. [3], in a cellular level model. Their results showed that the correct spatio-temporal localization of stimuli could facilitate the replay of a cortical sequence associated with the stimulation site, pointing to the need for peaks nesting for strengthening the neural connections. For this purpose, they propose stimuli application at the end of the SO's down-state, near the start of the SO up-state (0˚by sine reference), almost matching with our stimulation targetphase proposal.

Relation with experimental works and further work
Non-invasive methods of neural system stimulation operate at large spatial scales, and they can be by direct application of electrical or magnetic fields or by the use of sensory stimuli to modify neural activity [37,49]. The non-invasive techniques are under active-research because of the potential clinical applications in humans. Several works [17,18] tried to enhance the SWS activity to achieve memory improvements. There were applications of olfactory [50,51], auditory [3,33], and tactile stimuli [52] with similar or even better results than direct stimulation [53][54][55][56][57][58][59]. In counterpart with in-vivo works, one limitation of this work is the assumption of a steady-state of the sleeping brain during the entire simulation time, ignoring the changes known to occur within the same sleep stage.
Following the discoveries from previous works of Ngo et al., the application of stimulus at a precise phase of the ongoing cortical activity [3] produces better behavioral results than rhythmic stimulation [36,60]. Another result of the same authors shows that the application of more than two consecutive pulses does not improve physiological outcomes [61]. These results suggest the presence of mechanisms that prevent an over-driving of SO activity, which could be present in our case given a similar increase of I (SO) but lesser N SO /min. with STIM-R (see S1 Fig). A recent work [35] shows that this closed-loop stimulation pattern does not always enhance memory consolidation even with the increase of the power of the SWS events. Furthermore, based on the stimulation protocol of Ngo et al. [3], the recent work [62] identifies the phase at which the amplitude of SO increases, also observing an increase in the SP quantity and amplitude. They encountered the SO up-peak, or 90˚, as optimal timing for stimuli application.
Phase detection for closed-loop stimulation depends on the selected SO frequency band. Pulse stimulation at the SO rising slope increments its frequency power [3,33], which was also observed in our results for the three closed-loop conditions studied (see I (SO) ). However, precisely determine the phase signal to stimulate neither at zero-cross, nor peak highly depends on the selected SO frequency band.
Our results show that the shape, energy, and application phase of the stimulation pulses modifies the power and occurrence of SWS events. On the framework of EEG, these stimuli parameters are often neglected in sensory event-related potential (ERP) analysis, given the confounds with attention in awake experimental studies [63] and the incomplete knowledge about the codification pathways of sensory information. Despite this, the stimuli used in this work represent incoming firing responses to the thalamus, not the sensory stimuli properly. However, these stimulation characteristics are under research for neural implants, and direct stimulation [64]. Finally, the selected stimulation parameters need experimental verification to correlate with memory task performance.

Conclusion
The neural field model can reproduce SWS activity with the occurrence of slow oscillations and sleep spindles. Rhythmic, random, and closed-loop stimulation was applied using this model to enhance the SWS events. Closed-loop stimulation, with stimulus applied at zero degrees as target phase of the ongoing activity filtered at a specific frequency, increases the power and the number of both SWS events, together with their co-occurrences. Also, closedloop stimulation at other target-phases changes spindles' timing in respect to the down-peak of slow oscillations. Our results indicate that the number of slow oscillations decreases when the stimulation pulse energy surpasses more than six times the background noise power per second. On the other hand, increasing the stimulus energy increases the SPs' power and their number of occurrences. The stimulus pulse's promoted parameters to get a higher occurrence of both events is the decreasing ramp shape with a duration of 50 ms applied at target-phase zero of closed-loop stimulation.

Large-scale brain model: Neural field theory
Neural field theory represents the the activity of a whole population of neurons as its mean firing rate. Spatial propagation can then be calculated to obtain the spatio-temporal evolution of neuronal activity [27].
Thalamocortical interactions can be modeled using Neural Field Theory, where both cortical and thalamic regions are each represented by a pair of excitatory and inhibitory neural populations. The populations relate to known structures. In the cortex, the excitatory population e corresponds to the pyramidal neurons, while inhibitory population i represents the inter-neurons. In the thalamus, the excitatory population represents the relay nucleus s of incoming sensory information. Similarly, the thalamic inhibitory population represents the reticular nucleus r. Fig 1A (see Results) schematizes the connections between the four populations in the eirs-NFT model.
The dynamics of the mean soma voltage, V a , of the neural population a, with incoming activity from another population b (with a, b being elements of the set of neural populations p = {e, i, r, s}) is represented as (see also Fig 1B) 1 ab where α is the inverse of the decaying time of the impulse response, and β is the inverse of the rising time of the impulse response. Both parameters shape the post-synaptic-dendritic response. ν ab is the strength of the connection from population b to a, and its sign determines whether the connection is excitatory (positive connection strength) or inhibitory (negative connection strength). τ ab is the propagation delay between the populations a and b. Then, the model uses a sigmoid function to convert the soma voltage into the firing rate population response, Q a , as follows Q a ðr; tÞ ¼ S½V a ðr; tÞ� ¼ Q max 1 þ exp ðÀ ðV a ðr; tÞ À yÞ=s r Þ ; ð2Þ where Q max is the maximum firing response, and the sigmoid function approximates the normal cumulative probability distribution of the firing threshold voltage with mean θ and standard deviation σ where s r ¼ s ffi ffi ffi 3 p =p. The spatial propagation of the firing response only applies to the excitatory population e, the one with axons long enough for making it relevant [28]. It is expressed as 1 g 2 e @ 2 � e ðr; tÞ @t 2 þ 2 g e @� e ðr; tÞ @t þ � e ðr; tÞ À r 2 e r 2 � e ðr; tÞ ¼ Q e ðr; tÞ; ð3Þ where γ e = v e /r e is the temporal damping coefficient for pulses propagating at a velocity v e through axons with a characteristic range r e . For the other populations, ϕ a (r, t) = Q a (r, t).
The expressions Eqs (1)-(3) govern the dynamics of each neural population. Their parameters are physiologically related values [28] and are listed in Table 3. Note that the subscripts of the parameters follow the order: 'input to'-'output from'. Fig 1B (see Results) represents the connections inside a population and relationship between the variables. Table 3 presents the used parameters in our simulations of SWS activity on the eirs-NFT model.
Linear transfer function and loop-gains. To calculate the spectrum of the system's linear approximation with diverse sets of connections strengths, we used the steady-state solutions of the Eqs (1) and (3). The steady-state solution of the corticothalamic model comes from setting all temporal and spatial derivatives to zero. Defining � ð0Þ e and � ð0Þ n as the steady-state value of the cortical neural population and the input, respectively, the steady-state solution [29] is S À 1 ½� ð0Þ e � À ðn ee þ n ei Þ� ð0Þ e ¼ n es S½n se � ð0Þ e þ n sr S½n re � ð0Þ e þ ðn rs =n es ÞðS À 1 ½� ð0Þ e � À ðn ee þ n ei Þ� ð0Þ e Þ� þ n sn � ð0Þ n �: ð4Þ We searched the steady-state solutions that accomplish Eq (4) using the function fsolve of Matlab with an initial value for � ð0Þ e ¼ 10 s -1 . The linear approximation of the sigmoid function in Eq (2) uses its Taylor's serie expansion where ρ a is the slope of the sigmoid function evaluated at V ð0Þ a , the steady-state soma voltage of the population a. Note from Eq (3) that, in the steady-state, � ð0Þ a ¼ Q ð0Þ a and we can treat each variable of the model as the perturbation from their steady-state value, reducing the above expression to � a ðr; tÞ ¼ Q a ðr; tÞ ¼ r a V a ðr; tÞ: ð6Þ Parameters of the model used in all the simulations. The first three parameters correspond to the synaptic-dendritic and soma-voltage dynamics Eq (1)). Q max , σ ρ , and θ are parameters related to the activation function of the firing response (Eq (2)). r e and γ are parameters of the spatiotemporal propagation (Eq (3)). The following eight parameters are the connection strengths between the neural populations. Finally, the two last parameters are the mean and standard deviation (SD) of the Gaussian white noise input. https://doi.org/10.1371/journal.pcbi.1008758.t003 Taking the Fourier transform of Eqs (1), (3) and (6), and rearranging and eliminating ϕ r and ϕ s , leads to the system's spectrum with output ϕ e , and input ϕ n [29,30].
� e ðk; oÞ � n ðk; oÞ ¼ L 2 G es G sn exp ðiot 0 =2Þ ð1 À LG ei Þð1 À L 2 G sr G rs Þðd À k 2 r 2 e Þ ; ð7Þ where G ab = ρ a ν ab is the gain of the connection ab, and L is an abbreviation for L(ω), the reciprocal of the Fourier transform of the operators applied to V a (r, t) in Eq (1) and expressed as The linear approximation in Eq (7) uses a reduction of the model parameters with the same incoming connection strengths for the excitatory and inhibitory populations of the cortex, i. e. the connection strength ν ee = ν ie , ν ei = ν ii , and ν es = ν is . The gray arrows in Fig 1A represents the gain loops in Eq (7); G ei becomes an additional loop gain with the reduction of parameters. Further, the non-zero propagation delays considered in Eq (7) only correspond to the ones between the corticothalamic and thalamocortical populations, τ es = τ se = τ re = t 0 /2. Spatial coupling of stimuli input. The eirs-NFT model uses spatio-temporal variables, but our stimuli input is a single-channel time series. To give the stimuli a spatial representation, the signal goes through the sensory input module (the block at the left of Fig 1C in Results). We formed a receptive field 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 256 nodes with positive standard deviation σ E = 1 and negative standard deviation σ I = 2. The chosen center and standard deviations values avoid the grid boundaries.
Before entering the thalamocortical model, we added Gaussian white noise, v(r, t), forming the potential field of the nth population, ϕ n (r, t).
Simulation procedure. We solved the model dynamics in Eqs (1)-(3) using the Euler's integration method with a time step h = 10 −4 seconds, and a uniform sheet of N = 256 equidistant nodes in a square cortex of 0.5 meters per side using periodic boundary conditions (toroidal boundary conditions). The model was implemented on Matlab, and it is available on https://github.com/mjescobar/SWSNeuralField. The simulations were performed on a computer running Ubuntu Server 18.04 with a double Xeon E5-2630 v4 processor.
Each simulation lasts 910 seconds. We performed five simulations for every tested stimulation pattern u(t), each with a different random number generator seed, giving a total of 4550 seconds per stimulation case. The stimulation started after the first 5 seconds and stopped before the last 5 seconds in open loop stimulation. The initial conditions of neural populations variables come from a previous 6-second simulation with initial conditions equal to zero, and without stimulation.
We stored the output cortical activity ϕ e (r, t)), and the input signal ϕ n (r, t) with a sampling frequency of 100 samples per second.
We obtained an EEG-like signal, x(t), spatially averaging the output activity for all time points, and then we subtracted its temporal mean value where r is the grid position, N the number of nodes, and � � e ðr; tÞ the average of the cortical activity (see Fig 1C, module at the right).

Event detection in the time domain
We detected the slow-wave sleep events (slow oscillations and sleep spindles) from the output signal of the complete model. In the following, we refer as x(t) to the output signal of the stimulation cases (STIM), while b(t) represents the output signal of the non-stimulation cases or SHAM condition.
Detection of slow-oscillations events. As our EEG-like signal, x(t) is not in voltage units; then, we cannot use the classical method to detect SOs in EEG recordings. The classical method establishes SOs every time the signal has a negative peak below -40 μV, and a peak-topeak amplitude higher to 75 μV. As an alternative, we searched for events with time-lengths between 0.8 and 2 seconds (frequency: 0.5-1.25 Hz) [3], and we kept the higher amplitude oscillations.
The detection of SOs is based on the search for high amplitude oscillations in a delta-band sub-region. First, we filtered the signal in the SO frequency band (0.5-1.25 Hz, Chebyshev type I filter, fourth-order, 10 −6 dB allowed ripple, zero-phase lag). Second, we used the z-score of the filtered signal for searching the zero-crossing points. Later, we applied a threshold to the filtered signal, and we searched for negative peaks below the arbitrary low value of −10 −6 s -1 . Using the zero-crossing points, we calculated the peak-to-peak amplitude for the single oscillations that underpass the negative peak threshold. Finally, slow oscillations are confirmed if the peak-to-peak amplitude is higher than 1.25 times the baseline peak-to-peak amplitude average.
Detection of sleep spindles. Spindles are bursts of oscillations in the 9-16 Hz band. To detect them, we normalized each simulated output x(t) by the mean and variance of their corresponding b(t) signal. The results are single-channel time series expressed as z-scored EEG registers for each simulation. Then, we computed the RMS value of the z-score x(t) filtered in the SPs frequency band (9-16 Hz, Chebyshev type I filter, fourth-order, 10 −6 dB allowed ripple, zero-phase lag) and applied an amplitude threshold for the automatic detection of SPs [23,65]. The RMS value calculated at each sample uses a 0.2 s window length, and we also smoothed it with a Hamming window of the same time length [66]. SPs are then detected every time the signal exceeds a detection threshold of 1.25 standard deviations of the RMS value of the filtered baseline signal.
Counting of co-occurrent events. We defined co-occurrence of events every time SOs and SPs overlap at least 250ms. The counting process of co-occurrences uses tagged events to avoid duplication. We also extracted the time percentage of occurrence of SOs, SPs, and cooccurrences. From here, we use the identifier N SO for the number of slow oscillations, N SP for spindles, and N C for co-occurrence of events.

Power increases in the event's frequency bands
We define the scalogram power difference index (I (j) ), which relates the power of the oscillations in the jth frequency band as Eq (8) t2T jS x ðs; tÞj À jS b ðs; tÞj P s2s j P t2T jS x ðs; tÞj þ jS b ðs; tÞj where S x (s, t) is the STIM register scalogram, S b (s, t) is the baseline or SHAM register scalogram. s SO comprises the scales inside the 0.5-1.25 Hz spectrum related to slow oscillations and the s SP the range of scales correspondent to the frequency range of the sleep spindles of 9-16 Hz. T is the total simulation time.
We used the Morlet wavelet in the scalogram computation. This wavelet shape is a Gaussian-modulated complex sinusoid with template C(t) = A exp (−t 2 /(2σ 2 )) exp (jω 0 t). The mother wavelet has central frequency ω 0 = 15, using scaling factors, that result in a frequency resolution of 300 scales related to a starting frequency at 0.1 and ending at 30 Hz. The scales are equally distanced in a logarithmic range [67].

Closed-loop phase-lock driver
The closed-loop stimulation requires the on-line detection of instantaneous phase of the ongoing activity, to apply stimuli at the target phase. For this, we filtered x(t) with an IIR inversenotch filter for ω 0 = 2πf 0 h with poles z = {0.9999exp(±jω 0 ), 0}, and zeros z = {−1, −1, 1}, considering a scaling factor of 10 −4 . In addition, we extracted the envelope of the filtered signal at each sample to normalize itself. The envelope detector (see Fig 6A in Results) has three stages: (i) We detected the peaks of the rectified signal using its absolute value. (ii) We exponentially decreased the peak value, with a decay rate of 12.5 s -1 , until the detection of a new peak or the reaching of the limit of a sample counter. (iii) We normalized the filtered signal with the current envelope value at each sample.
On the normalized signal, z(t), an additional IIR filter gets its π/2 shifted-phase version y(t). This shift-phase filter has the poles z = {0.9999exp(±jω 0 ), 0.999}, and zeros z = {−0.5 exp (±0.495π), −0.9}. The filter scaling factor is the half quotient of the sum of the denominator coefficients over the sum of the initial numerator coefficients.
The stimulus pulse generator is triggered when the phase arrives inside a hysteresis window of ε = ±9 degrees (±5%) around the desired target-phase. The trigger waits for the phase 1.9 radians or 342 degrees to be reached before resetting for the next pulse application.
We applied the stimuli at a specific target-phase of a particular frequency in the closed-loop condition. The signal goes through a narrow band-pass filter and we detected the instantaneous phase with the arc-tangent with an online π/2-shifted phase filtered version of the narrow band signal. We normalized the amplitude of the signal before the arc-tangent by the envelope amplitude. The envelope value came from a peak detector, where the peak value decreases exponentially with a decay rate of 12.5 s -1 until the detection of a new peak.

Statistical tests
We tested for normality of each set of results for each stimulation pattern with the Shapiro-Wilk test.
To compare between pairs of stimulation patterns following a normal distribution, we used the Welch's t-test. On the contrary, when at least one pattern does not follow a normal distribution, we used the Wilcoxon signed-rank test. We used the Kolmogorov-Smirnoff test for the similarity between the distributions of time delays between the down-peak of the slow oscillation and the center of the coincident spindle.
All the statistical analysis were implemented using the library scipy.stats 1.5.2 in Python 3.7.

Supporting information
The supporting information contains results by other stimulation frequencies, random stimulation results not plotted in Figs 4 and 5, tables with values from statistical tests, and plots of the obtained results versus the onset amplitude of the stimulation pulses.