Thalamic deep brain stimulation modulates cycles of seizure risk in epilepsy

Chronic brain recordings suggest that seizure risk is not uniform, but rather varies systematically relative to daily (circadian) and multiday (multidien) cycles. Here, one human and seven dogs with naturally occurring epilepsy had continuous intracranial EEG (median 298 days) using novel implantable sensing and stimulation devices. Two pet dogs and the human subject received concurrent thalamic deep brain stimulation (DBS) over multiple months. All subjects had circadian and multiday cycles in the rate of interictal epileptiform spikes (IES). There was seizure phase locking to circadian and multiday IES cycles in five and seven out of eight subjects, respectively. Thalamic DBS modified circadian (all 3 subjects) and multiday (analysis limited to the human participant) IES cycles. DBS modified seizure clustering and circadian phase locking in the human subject. Multiscale cycles in brain excitability and seizure risk are features of human and canine epilepsy and are modifiable by thalamic DBS.


Results
Subjects. Seven dogs (D1-7) and one human (H1) met inclusion criteria. Subjects were monitored with the NV device (n = 3 dogs), NV transitioned to RC + S™ (n = 2 dogs), or RC + S™ device (n = 2 dogs; n = 1 human); both devices provide continuous intracranial EEG recordings (iEEG). The RC + S™ device provides both sensing and electrical stimulation. Subjects implanted with the NV device had bilateral subdural strip electrodes spanning the frontoparietal cortices; subjects transitioned from NV to RC + S had custom connectors to interface the RC + S with the NV strip electrodes. Subjects who were implanted exclusively with the RC + S™ device had penetrating depth electrodes targeting bilateral anterior thalamic nuclei (Medtronic 3387 or 3389 electrodes) and bilateral hippocampi (Hc) (Medtronic 3391 electrodes). The median recording duration was 298 days (range 51-395), and median seizure burden was 45 seizures (range 11-338) (Supplementary Table 1). D1-5 were monitored in a research kennel, and D6, D7, and H1 (all with RC + S™ system) were monitored in their home environments. Three dogs died during monitoring. D5 died secondary to a device related infection; D1 died from a traumatic head injury; D2 died as a complication of surgery during implantation of the RC + S system (subsequent to NV monitoring). Chronic brain recordings in human and canine epilepsy. Figure 1A provides a schematic of implanted chronic brain recording devices and shows a representative spontaneous seizure. The running hourly average IES rate is shown for H1 and D1 (Fig. 1B,C). Circadian fluctuations in the IES rate showing more frequent spiking during sleep/night are more clearly visualized with a magnified view (Fig. 1D,F), and seen as the high www.nature.com/scientificreports/ frequency oscillation in the hourly average IES rate tracing. Multiday fluctuations in IES rate are evident in the two-day moving average of the spike rate (Fig. 1B,C) and amplitude spectral density (Fig. 1E,G).
Circadian spike rate cycles and seizures. Circadian IES rate cycles were present in all subjects, and are seen as peaks in the wavelet transform of IES rate above the 95% threshold of normally distributed white noise ( Fig. 2; Supplementary Fig. 1 and methods for more detail). Figure 2B shows four of seven dogs and the human subject had 12-h IES rate ultradian cycles, and two of those dogs additionally had 8-h cycles. Interestingly, subject specific multiday cycles were observed in all subjects (further discussion below). Circadian IES rate cycles and seizure timing was evaluated relative to the hour of the day. Figure 3A shows the average normalized daily IES rate, with superimposed seizure onset relative to clock time. Two periods of the averaged cycle are displayed to aid visualization. Twelve-hour IES rate cycles are evident as minor daytime peaks and are most clearly seen in Dogs 6 and 7.  Circadian IES rate and seizure cycles. (A) Average normalized daily IES rate, with superimposed seizure onset (black circles) relative to clock time. Six hundred and twenty seven total seizures. To aid visualization two periods of the averaged cycle are displayed. Day/night periods are indicated by white/black bars at the bottom of each plot. Ultradian IES rate cycles are evident as minor daytime peaks in the canines and are most clearly seen in Dogs 6 and 7. The human subject has a pronounced nocturnal IES peak. (B) Seizure phase locking to the circadian IES rate cycle. Five subjects had significant seizure phase locking to the circadian spike rate cycle (marked by '*' next to the subject identifier in (A)). There was no clear group level phase preference of seizure timing relative to the circadian spike rate cycle. www.nature.com/scientificreports/ Seizure phase locking to the circadian spike rate cycle is shown in Fig. 3B. Five subjects had significant seizure phase locking to the circadian spike rate cycle (marked by '*' next to the subject identifier in Fig. 3A). There was no clear group level phase preference of seizure timing relative to the circadian spike rate cycle.
Multiday spike rate cycles and seizures. All subjects had significant multiday spike rate cycles (Figs. 1E,F, 2), which were also evident in the 2-day moving average trend of the hourly spike rate (Fig. 1B,C). Prior studies have shown inter-individual epilepsy chronotypes, including approximately weekly, 15-day, and 20-day cycle durations 10 . In our cohort, five subjects had significant weekly IES rate cycles, and five subjects had significant 2-3 weeklong cycles (Fig. 2). The human subject showed circadian and multiday IES cycles.
Significant IES rate cycles identified in Fig. 2 were further evaluated using the filter-Hilbert method, and the filtered multiday cycles reflect slow changes seen in the IES rate 2-day moving average (Fig. 1B,C). The combined circadian and multiday cycle data provide better estimates of the raw hourly spike rate trend when compared to an isolated cycle period, as seen in Fig. 1D,F.
Seizure phase locking to multiday spike rate cycles was evaluated in subjects with approximately weekly or 2-3 weeklong spike rate cycles. Representative circular histograms are shown in Fig. 4B, and the group level resultant vector is represented by a black arrow. Figure 4C shows the resultant vector of seizure timing relative to spike rate cycles for all subjects with weekly and 2-3 weeklong spike rate cycles. The group average vector is shown as a black dashed arrow. Four subjects had significant seizure phase locking to weekly spike rate cycles, and 5 subjects had phase locking to 2-3 weeklong spike rate cycles (P < 0.05); there was a group level seizure phase preference for the rising phase of multiday spike rate cycles.
Multiscale spike rate and seizure cycles in epilepsy. Seizure timing was evaluated relative to multiscale spike rate cycles (circadian and multiday). Phase-phase plots show seizure timing relative to the phase of circadian and multiday spike rate cycles for the four subjects with seizure phase locking to both circadian and multiday cycles (Fig. 4). Phase-phase plots show that seizure risk is maximal with circadian and multiday cycles in their respective high risk phase, with declining risk when moving away from the high risk phase of either cycle. Peak risk using phase-phase analyses was 5.8 fold greater relative to multiday cycle-only analysis, and 5.3 fold greater relative to circadian cycle-only analysis for D1; 3.1 and 5.1 for D3; 2.7 and 4.0 for D6; and 2.0 and 1.6 for H1 ( Supplementary Fig. 2).
Thalamic DBS and rhythms in epilepsy. Three community dwelling subjects (H1, D6, D7) living in their natural environment received ANT DBS during the study, in addition to a baseline recording period without DBS. Figure 5A shows the time-frequency spectrogram (CWT scalogram) of hourly average spike rate span- The resultant vector of seizure timing relative to IES rate cycles for subjects with weekly and 2-3 weeklong IES rate cycles (black dashed arrow is group average). Four subjects had significant seizure phase locking to weekly IES rate cycles, and 5 subjects had phase locking to 2-3 weeklong IES rate cycles (P < 0.05). There was a group level seizure phase preference for the rising phase of multiday IES rate cycles. (D) Phase-phase plots of seizure phase locking to circadian and multiday spike rate cycles shows the interplay of circadian and multiday cycles. As might be expected, multiday cycles are associated with seizures occurring during existing circadian patterns. D1 dog 1, H1 human 1. www.nature.com/scientificreports/ ning more than 50 weeks of Hc recordings from H1. Changes to thalamic stimulation were accompanied by clear changes in the amplitude of circadian and multiday (13-21 day periods) spike rate cycles. Following device placement, H1 had 112 days of baseline monitoring with no stimulation, 116 days of low frequency stimulation, and 110 days of high frequency stimulation. The two pet canines living with their owners, Dogs 6 and 7, had 206 and 339 days of baseline recording, and 52 and 0 days, and 0 days and 127 days of high frequency stimulation or low frequency stimulation, respectively (includes recording downtime). A decline in recording adherence during the stimulation-on periods for Dogs 6 and 7 prevented full assessment of multiday cycles with long cycle periods.
The spectral amplitude across stimulation epochs is shown in Fig. 5B. There was significant modulation of circadian (all subjects) and multiday (H1) IES rate cycle amplitude across thalamus stimulation epochs ( Fig. 5C; P < 10 -5 for all subjects; statistical significance level was maintained while controlling for the median IES rate for each stimulation state). Interruptions in chronic recordings during thalamic DBS of D6 and D7 limited the duration of multiday cycles that could be assessed, and seizure counts were insufficient for stimulation state subgroup analysis.
Thalamic stimulation had a significant impact on the relative burden of clustered seizures (seizures that occur within 24 h of a preceding seizure) relative to lead seizures for H1 (Fig. 5D)-there was an increased burden of clustered seizures during high frequency stimulation relative to baseline (Fisher's exact test P = 0.03) and low Seizure phase locking to the circadian spike rate cycle was anticorrelated with the amplitude of the multiday spike rate cycle: with low frequency stimulation (and low amplitude multiday spike rate cycle) the circadian seizure phase locking R value = 0.35, P value 3.7e−6; no stimulation (middle amplitude multiday spike rate cycle) R value = 0.29, P value 1.2e−4; high frequency stimulation (and high amplitude multiday spike rate cycle) R value = 0.22, P value = 1.3e−3 (phase doubling was used to account for the bimodal circular distribution). Figure 6 presents a model of seizure risk in which changes in the amplitude of multiscale IES rate cycles are positively correlated with seizure network excitability or seizure risk. In this model IES rate cycle amplitude and seizure risk cycle amplitude are positively correlated, but there can be phase shifts between these cycles (peak multiday seizure risk occurs at the rising phase of IES rate cycle; peak circadian seizure risk has subject specific phase preference (or may fall within several circadian chronotypes)). This model reconciles the changes in multiday IES rate cycle amplitude, seizure clustering, and circadian seizure phase locking seen in H1 across thalamic stimulation states.

Discussion
Using long-term brain recordings from novel implantable sensing and stimulation devices in one human and two pet dogs living in natural environments, and five kennel dwelling research dogs, we demonstrate that circadian and multiday cycles of seizure risk are features of human and canine epilepsy and can be modified by thalamic DBS. Furthermore, thalamic DBS was shown to modulate epilepsy cycles in a stimulation frequency dependent manner (high vs. low frequency stimulation). To the best of our knowledge this work provides the first evidence that neuromodulation can modulate cycles of IES rates and seizure risk, suggests that thalamocortical circuits have a role in regulating circadian and multiday rhythms in epilepsy, and expands the evidence that naturally occurring canine epilepsy can provide valuable insights into the multiscale temporal dynamics of seizures and epilepsy 7 .
Circadian IES rate cycles were present in all subjects. Behavioral state dependent changes in IES rate are well established in epilepsy, and sleep/wake related changes likely have a primary role in these cycles. However, beyond behavioral state, there is evidence for circadian cycles in epilepsy dependent on the endogenous circadian pacemaker 38,39 . Five of eight subjects had seizure phase locking to the circadian IES rate cycle. Recent work based on chronic brain recordings has further found that seizure likelihood varies with day-to-day changes in sleep duration 40 . There was not a consistent inter-subject pattern of the preferred phase of seizure occurrence, consistent with prior human work 10 . We are engaged in ongoing work to sleep stage these chronic recordings to further elucidate the impact of behavioral state on circadian cycles in epilepsy. Four dogs had 12-h IES rate cycles, which may correspond to daytime naps-naps are a common feature of canine sleep 22,41 . Interestingly, both community dwelling pet dogs (in contrast kennel dwelling research dogs) had particularly robust ultradian spike rate cycles, which might reflect more stable nap routines. , and in combination (center panels). The right panel phase-phase plots mark the phase angle above a theoretical high seizure risk threshold for both the short and long cycles. As the amplitude of the long seizure risk cycle increases (lower panel), a greater proportion of the short cycle falls above the threshold line, which would result in reduced seizure phase locking to the short cycle. During the high risk phase of the long cycle there is an increased risk of seizure clustering as a greater proportion of each short cycle falls above the high seizure risk threshold. In this model, the amplitude of seizure risk cycles is positively correlated with the amplitude of IES rate cycles, however, there are phase shifts between these cycles (for multiday IES rate cycles peak seizure risk occurs during the rising phase of cycles; for circadian IES rate cycles the phase of peak seizure risk is subject specific). IES interictal epileptiform spikes. www.nature.com/scientificreports/ Multiday IES rate cycles are a well-established feature of human epilepsy 5,6,8,10,11 , and this work provides the first evidence of multiday spike rate cycles in a cohort of dogs with naturally occurring epilepsy (data from D2 was included in a recent review paper on cycles in epilepsy 11 ). Multiday chronotypes with approximately weekly and 2-3 weeklong cycles have been described in a large cohort (222 subjects) of people with epilepsy undergoing monitoring with a clinical responsive neurostimulation system 10 . All subjects in our study had weekly and/ or 2-3 weeklong spike rate cycles, with seizure phase locking to these cycles in four and five subjects, respectively. Across subjects, seizure occurrences were increased during the rising phase of multiday spike rate cycles, consistent with prior human and rat studies 4,8 . As others have noted 4 , the long period lengths of these cycles preclude characterization during inpatient monitoring stays and may contribute to conflicting descriptions of the relationship between seizures and spike rates.
An important issue to emphasize is that in the absence of a Zeitgeber (from German, 'time giver'-environmental signals that synchronize or entrain biological rhythms to the external world), endogenous rhythms drift relative to clock time 42 . The gradual desynchronization of a biological rhythm (seizure risk) relative to clock time would limit the usefulness of seizure risk assessments based on a fixed clock time-based cycle. By evaluating seizure occurrences relative to a biomarker of endogenous seizure risk cycles with fine temporal resolution (spike rate), the problem of desynchronization from clock time is avoided. Additionally, by using a composite of neighboring spike rate cycle periods (center cycle period +/− 25%), this method accommodates for fluctuations in the cycle period around a central tendency 43 . These issues highlight the importance of temporally resolved biomarkers of seizure risk, i.e. chronic iEEG. Circadian seizure risk cycles are less disposed to desynchronization from clock time given the presence of powerful daily Zeitgebers, however, desynchronization from clock time could result from shifting behaviors relative to clock time.
Our spike rate and seizure data were generated from a particularly robust dataset of continuous full bandwidth iEEG recordings, using validated spike and seizure classifiers with epileptologist visual review for confirmation. The continuous full bandwidth research recordings support prior human work based on constrained chronic brain recordings from a clinical system 4,10,44,45 .
The combination of short (circadian) and long (multiday) timescale cycles can improve the characterization of seizure risk (Fig. 4D) 4 . Such multiscale analyses may hone assessments of seizure risk, reduce the unpredictability of epilepsy, and potentially prompt temporally informed therapeutic and behavioral interventions.
Little is known about the impact of electrical brain stimulation on cycles in epilepsy. In contrast to prior chronic brain recording studies (NV sense only system; Neuropace RNS ® system with constrained clinician defined-detectors that are commonly adjusted over time limiting comparisons between detector settings), the RC + S™ device is uniquely suited to evaluate the longitudinal impact of neuromodulation on cycles in epilepsy.
The RC + S™ device requires significant patient and caregiver support to maintain continuous iEEG recordings and streaming telemetry, and several extended gaps of non-recording time during DBS for D6 and D7 required that the data be analyzed in separate chunks, limiting the multiday cycle analysis. H1 maintained reliable iEEG recordings throughout the study, allowing full multiday analysis.
In this study H1, D6, and D7 had bilateral anterior nucleus of thalamus DBS and concurrent bilateral Hc iEEG recording. Circadian spike rate cycles were modulated in a DBS frequency dependent manner-low frequency DBS was accompanied by a reduction in the circadian spike rate cycle amplitude, while high frequency stimulation was accompanied by an increase in circadian spike rate cycle amplitude, in human and canine subjects. Prior work shows that ANT DBS can increase arousals during sleep and more work is needed to assess the potential role of sleep disruption on cycles in epilepsy 46 . Similar to the modulation of circadian spike rate cycles by thalamic DBS, there was DBS frequency dependent modulation of multiday spike rate cycle amplitude (analysis limited to H1). There was a marked increase in the multiday (17-day) spike rate cycle amplitude during high frequency thalamic DBS, and a reduction in amplitude during low frequency stimulation, when compared to baseline. Additionally, seizure clustering increased during high frequency DBS compared to baseline and low frequency DBS. And finally, seizure phase locking to the circadian spike rate cycle was anti-correlated with multiday spike rate cycles amplitude, with increased seizure circadian phase locking during low frequency DBS (and low multiday spike rate cycle amplitude), and reduced phase locking during high frequency DBS (and high multiday spike rate cycle amplitude), relative to baseline.
These findings are consistent with a model of seizure risk in which changes in the amplitude of multiscale spike rate cycles is positively correlated with cycles of seizure network excitability or seizure risk (Fig. 6). As the multiday seizure risk cycle amplitude increases (Fig. 6 lower panels), there is an increase in the proportion of the circadian cycle that falls above the seizure threshold. This model suggests that the increased amplitude of multiday spike rate cycles, resulting in a larger circadian phase angle above the seizure threshold, would increase the risk of seizure clustering and reduce circadian seizure phase locking, as seen in H1. These findings suggest that beyond seizure phase locking to multiday and circadian spike rate cycles, changes in spike rate cycle amplitude may influence seizure risk and be a biomarker to guide epilepsy management. This extends prior work on multiscale cycles of seizure risk 47 and proposes a framework to understand the relationship between changes in multiscale cycle amplitude, seizure clusters, and seizure phase locking.
There is evidence that the thalamus plays an important role in focal epilepsy 27 , with evidence of enhanced structural and functional thalamocortical connectivity in individuals with temporal lobe epilepsy [48][49][50] , thalamus involvement in the initiation and propagation of seizures 28,29 , and numerous studies showing alterations of thalamocortical circuits in individuals with epilepsy 26,51,52 . Here, we provide the first evidence that thalamic DBS can modulates circadian and multiday cycles in epilepsy. This finding raises interesting questions about the role of thalamocortical circuits in the regulation of long timescale dynamics of epilepsy.
DBS is an established therapy for the management of focal epilepsy 27 and various neurological and psychiatric disorders 53 , however, DBS mechanisms of action are unresolved 54 www.nature.com/scientificreports/ oscillations) 31 , increased seizure threshold to chemotoxin 56 , and reduced hippocampal synchronization 31,57 . These investigations relied on acute stimulation protocols, and little is known about the electrophysiological changes associated with chronic stimulation. The pivotal study for ANT DBS 27 delivered high frequency stimulation on a duty cycle (1 min on, 5 min off) without clear supporting evidence; DBS for movement disorders is typically provided using continuous high frequency stimulation. Less is known about low frequency DBS, however there is empirical evidence to support low frequency cortical and subcortical stimulation for the treatment of epilepsy 31,34,35,58,59 , and evidence that noninvasive low frequency stimulation depresses cortical excitability 32 . In our study, subjects received continuous (not duty cycle) high or low frequency stimulation. Continuous stimulation was used in part due to device programming constraints, and was considered acceptable given the established safety of continuous high frequency stimulation based movement disorders protocols, and the lack of a theoretical framework for duty cycle stimulation.
The increased amplitude of IES rate cycles in H1 and D6 and increased seizure clustering in H1 with high frequency stimulation was unexpected. In accordance with the model above (Fig. 6), the apparent accentuation of epilepsy cycles could result in a greater proportion of the circadian seizure risk cycle falling above a seizure risk threshold, with increased seizure clustering and reduced circadian seizure phase locking, as was seen in H1. We did find desynchronization of hippocampal activity and reduced spike rate with acute high frequency stimulation similar to prior work 31 , however, this effect dissipated over time with chronic continuous high frequency stimulation (Supplementary Fig. 3), in contrast to chronic duty cycle stimulation 60 . This may suggest that epileptic networks are able to adapt to network desynchronization from intermittent vs. continuous high frequency stimulation. Additionally, data from transcranial magnetic stimulation (TMS) 61 suggests that high frequency stimulation can be excitatory, in contrast to the more inhibitory influence of low frequency stimulation. Further work is needed to study the different mechanisms and effects of continuous vs. duty cycle and low frequency vs. high frequency stimulation. Given these data, continuous high frequency stimulation should be used with caution.
Little is known about the mechanisms of multiday cycles. Modulation of multiday cycles in epilepsy by thalamic stimulation implicates thalamocortical circuits and intrinsic brain network dynamics in the regulation of these long cycles. This does not exclude a role for other cyclical endogenous or exogenous factors. It does, however, suggest that brain stimulation can influence the neurophysiology of cycles in epilepsy.
This work is limited by the number of subjects with chronic iEEG and concurrent thalamic DBS. Despite the subject count, these chronic recordings form a large dataset with 627 seizures recorded over 2146 days. Still, this work will benefit from continued iEEG monitoring and additional trials of thalamic DBS. The investigational devices used in the study require daily charging and device management by participants to maintain continuous full bandwidth recording and subject data did have gaps in recordings (Supplementary Table 1). Nonrecording segments were filled by normally distributed noise as previously described 16 which could impact assessments of IES rate cycles. For extended gaps, data were analyzed in independent epochs, and this did preclude analysis of the impact of thalamic stimulation on multiday cycles in D6 and D7. H1 maintained reliable iEEG recordings throughout the study without extended nonrecording periods, allowing full multiday analysis. Seizure cycles were characterized using non-causal filters, consistent with prior studies 4,6,10,44 , and future efforts will use causal methods compatible with prospective cycle characterization. We used objective electrographic iEEG seizures in this study, given the unreliability of patient reported seizure diaries and the challenge of tracking clinical seizures in canines. We did not track the behavioral state of subjects in the study so we cannot characterize the role of sleep in cycles in epilepsy, and our group's automated iEEG based sleep classification efforts are ongoing.
Our results demonstrate that cycles of seizure risk are common features of canine and human epilepsy, which can be modified by thalamic DBS in a frequency dependent manner. This work motivates further study of thalamocortical circuits in the regulation of brain excitability and seizure risk. Naturally occurring canine epilepsy can provide important insights into the temporal dynamics of seizure risk, and next generation implantable sensing and stimulation devices may improve our understanding of cycles in epilepsy and enable new treatment paradigms that leverage the temporal profile of seizure risk and adaptively adjust stimulation to prevent seizures.

Methods
Subjects. Nineteen dogs with naturally occurring epilepsy and one human with epilepsy were implanted and monitored with the RC + S™ or NV device (Supplementary Table 1), described previously 7 . Three dogs and the human participant lived in their home environments. All other dogs were housed in a research kennel. Two dogs underwent recording with both the NV and RC + S™ devices in sequence; the lead positions were kept however the leads themselves were exchanged due to cross-device incompatibility. Inclusion criteria required recording of at least 10 electrographic seizures, and greater than 50 days of iEEG. This study was completed in accordance with relevant guidelines and regulations. Experimental protocols pertaining to the human participant were approved by the Food and Drug Administration (IDE G180224) and the Mayo Clinic Institutional Review Board. This study and experimental protocols pertaining to the canine subjects were approved by the Mayo Clinic, University of California, Davis, and the University of Minnesota Institutional Animal Care and Use Committees. Methods were carried out in accordance with the ARRIVE guidelines. All participants and/or legal guardians provided written informed consent.
Interictal spike and electrographic seizure classification. Data consisted of continuous full bandwidth (250-500 Hz sampling rate) iEEG recorded by the NV or RC + S™ systems. Interictal spikes were identified using a validated automated spike classifier 62 . Electrographic seizures were identified by a validated, highsensitivity automated seizure classifier 63 , followed by board certified epileptologist review(G.W.), as previously described 30  www.nature.com/scientificreports/ rate data were excluded for hours in which data drops accounted for more than 50% of the block. Contacts with inconsistent recordings, artifacts, or lack of seizure involvement were excluded.

Frequency domain analyses.
Frequency domain analyses of IES rate data were performed using the continuous wavelet transform (CWT) and filter-Hilbert transform method. Both of these methods provide instantaneous phase, frequency and amplitude information for time-frequency analyses. The CWT was implemented using analytic Morlet wavelets with 10 voices per octave, L2 normalization, and minimum/maximum scales determined by the energy spread of the wavelet in frequency and time, as described 64 . The filter-Hilbert transform method involved narrowband filtering of the EEG signal prior to finding instantaneous phases and amplitudes from the analytic signal generated by the Hilbert transform. Filtering was implemented using least-squares finite impulse response (FIR) bandpass filters of order-3 with zero-phase shift (non-causal). The bandpass center period was incremented by 0.1 days between 0.3 and 2 days, 0.5 days between 2.5 and 10 days, and by 1 day for periods longer than 10 days. The maximum period duration was defined by the recording duration and temporal spread of FIR filters. When computing the CWT and Hilbert transform, missing hourly IES rate data were filled with random noise drawn from the normal distribution, with mean and standard deviation equal to the distribution of all recorded data from the hour of the day being filled, as previously described 16 . For long non-recording periods (> 15 days) the adjacent intact data segments were analyzed independently. For analyses using instantaneous phase/power data from the CWT and Hilbert transform, data within the cone of influence (COI) of boundary effects were excluded (COI defined as √ 2s where s is the CWT scale or FIR passband center period-COI shown in Fig. 5A. Circadian and multiday spike rate cycles. The CWT was used to assess circadian and multiday cycles in IES rate. Significant cycles were defined as relative maxima in the IES rate CWT amplitude spectral density (ASD) that were above the 95th percentile of white noise. ASD was calculated as the time averaged amplitude of the CWT of the spike rate. To define the 95th percentile threshold of white noise simulations of the CWT ASD (n = 1000) were run on data drawn from the standard normal distribution, and the 95th percentile value was identified for each period duration ( Supplementary Fig. 1).
Seizure phase locking to circadian and multiday spike rate cycles. The CWT ASD was used to identify significant circadian and multiday IES rate cycles. Identified IES rate cycles were further evaluated using the filter-Hilbert method, and seizure timing was evaluated relative to the phase of IES rate cycles. For each cycle period a composite of FIR filtered IES rate data (spanning period lengths of center frequency +/− 25%) were used to calculate the Hilbert transform analytic signal. This composite cycles data accommodates for drift of endogenous free running rhythms around a central tendency as previously described 4 -for example the precise duration of a "monthly" cycle may vary over time around a central trend. Seizure timing was evaluated with respect to the phase of these circadian and multiday oscillations in IES rate, and circular statistics methods characterized the preferred phase and circular nonuniformity of seizure occurrences relative to IES rate cycles (Fig. 4A), analogous to prior work 4 . Prior work on human epilepsy has demonstrated group level epilepsy chronotypes, including approximately weekly, 15 day, and 20 day cycle durations 10 . Here we emphasized analyses of significant multiday rhythms with weekly and 2-3 weeklong cycle periods. For instances with two significant cycles within a period range of interest, the cycle with greater amplitude was selected.
Circadian seizure and spike rate cycle. Circadian IES rate and seizure cycles were evaluated using two methods. Method 1 is the same as outlined above, in which a composite of data from neighboring IES rate cycles (centered around the 24-h period) are used to define the analytic signal-here, events are not evaluated relative to fixed hours of the day, but rather the phase is defined by the circadian (approximately 24-h) oscillation in the IES rate (the exact duration and phase of this physiological circadian rhythm may shift relative to clock time from day to day). Circadian cycles are entrained by a powerful exogenous Zeitgeber (24-h rotation of the earth). Given this synchronizing Zeitgeber, method 2 evaluated changes in IES rate and seizures relative to daily clock time, and were not defined by the circadian IES rate cycles itself, as in method 1. Method 2 was used for the evaluation of circadian rhythms shown in Fig. 3A, otherwise analyses used method 1.
Circular statistics. Circular statistics analyze angular data-in our case the phase of oscillations in the rate of IES. Circular statistics account for the phase similarity between 359° and 1° (Fig. 4A) 65,66 . To perform phase analyses of seizures with respect to IES rate cycles, seizure timing is evaluated with respect to the instantaneous phase of the underlying IES rate cycle of interest (from filter-Hilbert analytic signal), and the phase is assigned to the 360° arc of a circle. Circular histograms provide a representation of the phase preference and spread of seizures relative to circadian and multiday IES rate cycles (Fig. 4B). The R value (phase locking value) quantifies the circular nonuniformity of events. The R value is equal to the mean resultant vector of seizure events (uniform arbitrary amplitude; phase corresponding to the phase of the cycle of interest). Randomly distributed events will have an R value of zero, perfectly phase-locked events will have an R value of 1 4,65,66 . For evaluations of circadian seizure timing, the analysis was implemented as outlined above (relative to the circadian oscillation in IES rate), and additionally, seizures were evaluated relative to daily clock time. Statistical significance was assessed by the Rayleigh test 65 , which determines the R value required to reject the null hypothesis that events have a uniform circular distribution, assuming von Mises distributed data. www.nature.com/scientificreports/ Phase-phase plots. For individuals with significant phase locking of seizures to circadian and weekly or 2-3 weeklong IES rate cycles, phase-phase plots were used to further evaluate how combined IES rate cycles inform seizure risk. The circadian and multiday phase data of each seizure was used to generate 3-dimensional linearly interpolated histograms.
Anterior thalamus stimulation and rhythms in epilepsy. The human participant and two dogs (Dog 6 and Dog 7) underwent chronic ANT therapeutic stimulation. These subjects completed extended baseline "stim-off " recording, followed by extended "stim-on" monitoring. Subjects received bipolar stimulation to bilateral thalami, with the cathodal electrode selected to be positioned within the ANT-stimulation was delivered continuously at low-frequency (< 10 Hz) and high frequency (> 100 Hz) settings. Monopolar stimulation can introduce stimulation artifact into recording channels and bipolar stimulation was used to limit stimulation artifact. RC + S™ delivers current clamped stimulation, and stimulation amplitudes ranged between 2 and 6 mA. For the human subject, low amplitude stimulation ≤ 3 mA was considered to be subtherapeutic and was classified as baseline along with stimulation off time periods. The impact of ANT DBS on circadian and multiday IES rate cycles were assessed using frequency domain and circular statistics methods described above. The impact of ANT DBS on seizure rate and seizure clustering was evaluated using statistics described below.
Statistics. Significant IES rate cycles were defined as relative maxima in the IES rate ASD with amplitudes above the 95th percentile of normally distributed white noise. A one-thousand trial Monte-Carlo simulation was used to define the white noise 95th percentile ( Supplementary Fig. 1). For assessments of seizure phase locking to IES rate cycles, circular nonuniformity was assessed by the Rayleigh test 65 . The Rayleigh test determines the amplitude of the resultant vector needed to reject the null hypothesis that events have a uniform circular distribution 65 . One-way analysis of variance (ANOVA) was used to assess for changes in the spectral amplitude of circadian and multiday IES rate cycles relative to thalamic DBS stimulation state. Fisher's exact test assessed the burden of clustered vs. lead seizures during baseline monitoring and during thalamic DBS (low frequency and high frequency stimulation). All statistics were evaluated at 0.05 significance level. Analyses were performed using MATLAB (version 2020b, Mathworks Inc, Natick, MA, USA). www.nature.com/scientificreports/