Patent Application
20220338786
Deep brain stimulation (DBS) has proven to be effective in treating several neurological disorders, including Parkinson's Disease, epilepsy, and Obsessive Compulsive Disorder. Advances in device development over the last decade have enabled concurrent stimulation and recording at adjacent locations in the brain. These chronic neural recordings present an opportunity for identification of neural signatures of behaviors relevant to disease states. However, DBS and other electrical stimulations creates high amplitude, high frequency stimulation artifacts that contaminate the electrical recordings of the brain. As a result, there is a need for effective periodic artifact removal systems and methods for recovering the underlying neural signal and thus, for understanding the biological mechanisms underlying the qualitative successes of DBS.
In DBS applications, the artifact removal problem is complicated by a variety of factors. For example, while the stimulation frequency can be set by the DBS device, the device setting provides an estimate of the stimulation frequency (and period) that is not accurate enough for successful artifact removal. Furthermore, the amplitude of DBS artifacts are often orders of magnitude larger than the underlying neural activity, and the sampling rate of the neural recordings is often low compared to the high frequency of therapeutic DBS. For example, power constraints of DBS devices often require low sampling rates (e.g., 200-250 Hz), which in turn alias the stimulation frequency near the frequencies of underlying signals of interest. In addition, DBS recordings are often broken into many time segments with unknown phases, e.g., due to modulation of device settings or missing data of unknown length, the latter of which is a common limitation of the wireless data transmission methods used by some DBS devices.
There is a desire and need for computationally efficient artifact removal systems and methods that are effective in removing periodic artifacts from DBS and other electrical stimulation recordings (e.g., sampled at less than twice the frequency of DBS). Moreover, as periodic artifacts in DBS recordings may often cause multiple unknown phase shifts, there is a desire and need for computationally efficient estimation of periodic artifacts in the presence of such phase shifts. Furthermore, there is a desire and need for systems and methods of recovering and/or assessing the scale of lost packets during transmission of waveform data produced by DBS and other electrical stimulation. This assessment of lost packets may facilitate data alignment in the artifact removal process and in the recovery of neural signals necessary for, e.g. control of therapy.
Various embodiments of the present disclosure describe systems and methods that address one or more of the above described shortcomings.
The present disclosure presents new and innovative systems and methods for identifying and removing DBS artifacts from electrical measurements of brain activity. The system may perform several steps or processes to identify and remove DBS artifacts, including (1) generating an initial guess for a period of artifacts created by a DBS therapy signal, (2) determining the true period for the artifacts created by the DBS therapy signal based on the initial guess, and removing artifacts based on the true period.
In at least one embodiment, a system for stimulation period-based artifact removal and reconstruction is disclosed. The system may include intracranial electroencephalography (iEEG) device; one or more processors; and memory. The memory may store instructions that, when executed by the one or more processors, cause the one or more processors to: receive, in real-time from the iEEG device, waveform data caused by a deep brain stimulation of a patient-specific area of interest; and determine a stimulation period relative to a sampling rate. The stimulation period may be determined using an iterative (e.g., stage-wise) process, and may utilize regulation of the waveform data. The instructions may further cause the one or more processor to identify, by applying Nadaraya-Watson kernel regression and a linear filter to the waveform data and the stimulation period, a stimulation artifact; subtract the identified stimulation artifact from the received waveform data; and generate, in real-time, a filtered waveform data indicating an absence of the stimulation artifact.
In some embodiments, the instructions, when executed, may further cause the one or more processors to determine the stimulation period by identifying, between two contiguous time segments of the waveform data, one or more phase shifts in the waveform data; estimating, based on the one or more phase shifts, a number of time points in one or more missing packets of unknown duration; and realigning, based on the number of time points in the one or more missing packets, the waveform data. The stimulation period may be based on the realigned waveform data.
In another embodiment, a method of period-based estimation of electrical stimulation artifacts in the presence of phase shifts is disclosed. The method may include: receiving, in real-time by a computing device having a processor, waveform data caused by a deep brain stimulation of a patient-specific area of interest. The waveform data may comprise: a plurality of runs, a plurality of phase shifts, and a periodic artifact. Each run may indicate a contiguous portion of the waveform data separated by a packet loss. The method may further include generating a model of the received waveform data based on: a neural signal of interest and a model of the periodic artifact. The model of the periodic artifact may be based on a stimulation period of the periodic artifact and a phase shift of the plurality of phase shifts. For example, a model of the periodic artifact may be defined as:
where δi* comprises the phase shift, of the plurality of phase shifts, between the 0th and i-th run of the plurality of runs, A is the model of the periodic artifact based on a stimulation period
of the periodic artifact and the phase shift (δi*) of the plurality of phase shifts, Bi(t) is the neural signal of interest at the i-th run, and ηi(t) is noise at the i-th run. The method may further include applying harmonic regression to optimize the model of the periodic artifact. In some embodiments, the model of the periodic artifact to which harmonic regression is applied may be characterized as a(t|ξ,δ,α0,αk,βk,K)=α0+Σk=1K(αk cos(2πk(ξt+δ))+βk sin(2πk(ξt+δ))). wherein K is a finite number of harmonics, k. The method may further include determining simultaneously, using an optimization of a loss model, the plurality of phase shifts and the neural signal of interest, wherein the loss model is based on the model of the received waveform data and an optimized model of the periodic artifact. For example, the loss function model may be characterized as L(ξ,δi,θ)=Σi=0nΣtϵT
In another embodiment, the present disclosure presents new and innovative systems and methods for improved packet loss estimation. In one aspect a method is provided that includes dividing time series data into a plurality of runs. A period of stimulation may be estimated from data in the plurality of runs. A harmonic regression model may be fit to a longest run of the plurality of runs and packet loss may be estimated based on the harmonic regression mode.
In another aspect, a method for periodic estimation of lost packets included: receiving, in real-time by a computing device having a processor, waveform data caused by a deep brain stimulation of a patient-specific area of interest. The waveform data comprises a plurality of runs. Each run may indicate a contiguous portion of the waveform data separated by a packet loss. The method may further include: determining, by the computing device, a stimulation period relative to a sampling rate. Furthermore, the computing device may apply, based on the stimulation period, a harmonic regression model to a longest run of the plurality of runs; and determine, in real-time and based on the harmonic regression model, an optimal size of packet loss for the waveform data.
The features and advantages described herein are not all-inclusive and, in particular, many additional features and advantages will be apparent to one of ordinary skill in the art in view of the figures and description. Moreover, it should be noted that the language used in the specification has been principally selected for readability and instructional purposes, and not to limit the scope of the disclosed subject matter.
Today, state-of-the-art artifact removal methods from electrical stimulation therapies (e.g., deep brain stimulation (DBS) may often rely on template subtraction of the artifact at each electrical stimulation (e.g., DBS) pulse. However, existing methods for identifying individual pulses (e.g. thresholding) are not robust for low sampling rates or in the presence of other spurious high amplitude artifacts. Various embodiments of the present disclosure describe a novel and nonobvious systems and methods that identify each instance of the artifact by windowing the signal based on the exact period of stimulation. Finding the true period is not straightforward due to discrepancies between reported sampling rate and stimulation frequency of existing DBS and other electrical stimulation devices. There is a desire and need for systems and methods of finding the true period of stimulation. The presently disclosed systems and methods effectively and efficiently find the true period of any periodic signal, including that of DBS artifact. After the true period is found, a template of the artifact can be reconstructed and subtracted at each instance of the DBS artifact without contaminating the underlying signal of interest. The presently disclosed system and method is robust for non-periodic spurious high amplitude artifacts and for low sampling rates. The recovered underlying neural signal of interest is not contaminated even when the frequency content of the aliased artifact and underlying neural signal overlaps. Additionally, the presently disclosed systems and methods are computationally efficient and could be implemented onboard existing or future DBS devices to remove DBS artifacts in real time. The present disclosure also provides artifact removal methods that are effective in removing periodic artifacts from local field potential (LFP) and other electrical stimulation recordings sampled at less than twice the frequency of the electrical field stimulation (e.g., DBS).
The presently disclosed systems and methods can be used by DBS device companies and researchers to implement offline or onboard artifact subtraction to enable biomarker detection during concurrent sensing and DBS therapy. Additionally, the presently disclosed systems and methods could be used in other fields to recover any type of signal with a period artifact. These techniques are discussed in further detail below.
High amplitude and high frequency DBS artifacts complicate biomarker detection necessary for the development of adaptive DBS therapy. DBS artifacts can be difficult to remove, especially in low sampling rate recordings where components of the artifact are aliased. Existing methods rely on finding individual DBS pulses via thresholding which may be ineffective for the low sampling rates and long artifacts characteristic of electrical stimulation recordings (e.g., LFP)
Various embodiments of the present disclosure describe a novel and nonobvious method (referred to herein as period-based artifact reconstruction and removal method (PARRM)), and systems for implementing the same, which improve on existing methodologies by accurately finding each instance of the artifact. Experiments and results described herein demonstrate that PARRM is effective in removing high frequency DBS artifact from relatively low sample rate electrical stimulation recordings (e.g., DBS recordings, LFP recordings, etc.), without introducing contamination to the underlying neural signal.
Biomarker detection requires a hardware platform capable of concurrent sensing and stimulation. The amplitude of DBS therapy is often orders of magnitude greater than underlying neural signals causing the signal of interest to be heavily contaminated by artifact. In high resolution recordings, DBS artifact is typically removed using a simple low-pass filter. However, in low sampling rate recordings, the artifact is aliased into biologically significant frequency bands requiring more complex approaches to artifact removal. Stimulation artifacts may need to be removed from electrical stimulation recordings (e.g., LFP recordings) in order to identify neural biomarkers of disease symptoms and side effects of DBS. The period may be guessed by dividing the sampling rate by the stimulation frequency. Slight differences in external factors may make the sampling rate inexact. The true period may need to be found in order to produce accurate templates. A grid search may be done around the initial guess to find the best period. Periods may be evaluated via linear regression with a sum of sinusoidal harmonics of the period. The period which minimizes mean squared error with the raw data is chosen. Experimental results demonstrate that false periods can mimic harmonics of the artifact over short timespans. In some embodiments, to avoid such ‘distractor’ solutions, a penalty may be included for regression coefficients for high frequency sinusoids which are less significant for the true solution. After an optimal stimulation period is found, a region around each sample can be used to construct a local template. The value of the template corresponding to the phase of the sample can be subtracted in order to recover the artifact free recording.
Furthermore, it is expected that periodic artifacts in received waveform data from DBS and other electrical stimulations may cause multiple undesired and unknown phase shifts. Computationally efficient systems and methods are disclosed for period-based estimation of electrical stimulation artifacts in the presence of these multiple phase shifts.
Furthermore, as waveform data caused by deep brain stimulation is processed (e.g., to remove stimulation artifacts via PARRM), it can be expected that portions of the original waveform data may be lost in the transmission process. Various embodiments of the present disclosure describe novel and nonobvious systems and methods for estimating packet loss using the stimulation period (the systems and methods referred to herein as periodic estimation of lost packets (PELP)).
For example, during wireless transmission, neural data samples can be grouped into formatted units called ‘packets’. Packets contain a series of subsequent samples of a particular length as well as timing information and other relevant metadata. When transmitted, it is possible for packets to fail to reach the receiver leading to missing samples. These missing samples need to be properly accounted for when the time series is reconstructed. The timing information contained in each packet aids in this process but may be frequently inexact resulting in uncertainty in the number and location of the samples missing from a recording.
Particularly in less controlled environments away from the clinic, recordings can be especially prone to packet losses. Lower sampling rates may reduce the number of dropped packets and increase transmission ranges. However, it may still be typical for as much as 5% of the data to be lost even with such adjustments. Failure to accurately account for packet losses leads to the introduction of timing inaccuracies, artifacts during filtering, and reduced ability to identify meaningful neural signals.
Various embodiments for Periodic Estimation of Lost Packets (PELP) are disclosed herein, which include systems and methods for precisely estimating packet losses in bidirectional recordings from implanted devices where stimulation is active. PELP leverages a data driven procedure for determining the period of stimulation and the knowledge that stimulation continues identically during regions where data are missing to accurately account for the number of samples missing due to each dropped packet. As will be discussed herein, PELP may be robust across a range of amplitude ratios between stimulation and signal, pulse to pulse variations in stimulation amplitude, drift in stimulation frequency, and uncertainties in loss size estimates.
Streaming of intracranial electrophysiology data in the clinic and at home in ecologically valid environments may be essential for biomarker discovery in a variety of neurological disorders. Bidirectional implanted devices have enabled the acquisition of such datasets, however, data losses during wireless streaming hinder accurate analyses of neural signals. To address these issues, systems and methods for PELP are disclosed to precisely estimate and account for data losses from implanted recordings where stimulation is on.
Successful applications of PELP to recorded data are described herein that to using MEDTRONIC's SUMMIT RC+S from a human participant performing a behavioral task both in the clinic and at home to exactly estimate every occurrence of packet loss. PELP is widely applicable to other stimulating devices capable of wireless data streaming.
This disclosure contemplates any suitable number of computer systems 100. This disclosure contemplates the computer system 100 taking any suitable physical form. As example and not by way of limitation, the computer system 100 may be an embedded computer system, a system-on-chip (SOC), a single-board computer system (SBC) (such as, for example, a computer-on-module (COM) or system-on-module (SOM)), a desktop computer system, a laptop or notebook computer system, an interactive kiosk, a mainframe, a mesh of computer systems, a mobile telephone, a personal digital assistant (PDA), a server, a tablet computer system, an augmented/virtual reality device, or a combination of two or more of these. Where appropriate, the computer system 100 may include one or more computer systems 100; be unitary or distributed; span multiple locations; span multiple machines; span multiple data centers; or reside in a cloud, which may include one or more cloud components in one or more networks. Where appropriate, one or more computer systems 100 may perform without substantial spatial or temporal limitation one or more steps or processes of one or more methods described or illustrated herein. As an example and not by way of limitation, one or more computer systems 100 may perform in real time or in batch mode one or more steps, processes, or blocks of one or more methods described or illustrated herein. One or more computer systems 100 may perform at different times or at different locations one or more steps, processes, or blocks of one or more methods described or illustrated herein, where appropriate.
In particular embodiments, computer system 100 includes a processor 106, memory 104, storage 108, an input/output (I/O) interface 110, and a communication interface 112. Although this disclosure describes and illustrates a particular computer system having a particular number of particular components in a particular arrangement, this disclosure contemplates any suitable computer system having any suitable number of any suitable components in any suitable arrangement.
In particular embodiments, the processor 106 includes hardware for executing instructions, such as those making up a computer program. As an example and not by way of limitation, to execute instructions, the processor 106 may retrieve (or fetch) the instructions from an internal register, an internal cache, memory 104, or storage 108; decode and execute the instructions; and then write one or more results to an internal register, internal cache, memory 104, or storage 108. In particular embodiments, the processor 106 may include one or more internal caches for data, instructions, or addresses. This disclosure contemplates the processor 106 including any suitable number of any suitable internal caches, where appropriate. As an example and not by way of limitation, the processor 106 may include one or more instruction caches, one or more data caches, and one or more translation lookaside buffers (TLBs). Instructions in the instruction caches may be copies of instructions in memory 104 or storage 108, and the instruction caches may speed up retrieval of those instructions by the processor 106. Data in the data caches may be copies of data in memory 104 or storage 108 that are to be operated on by computer instructions; the results of previous instructions executed by the processor 106 that are accessible to subsequent instructions or for writing to memory 104 or storage 108; or any other suitable data. The data caches may speed up read or write operations by the processor 106. The TLBs may speed up virtual-address translation for the processor 106. In particular embodiments, processor 106 may include one or more internal registers for data, instructions, or addresses. This disclosure contemplates the processor 106 including any suitable number of any suitable internal registers, where appropriate. Where appropriate, the processor 106 may include one or more arithmetic logic units (ALUs), be a multi-core processor, or include one or more processors 106. Although this disclosure describes and illustrates a particular processor, this disclosure contemplates any suitable processor.
In particular embodiments, the memory 104 includes main memory for storing instructions for the processor 106 to execute or data for processor 106 to operate on. As an example, and not by way of limitation, computer system 100 may load instructions from storage 108 or another source (such as another computer system 100) to the memory 104. The processor 106 may then load the instructions from the memory 104 to an internal register or internal cache. To execute the instructions, the processor 106 may retrieve the instructions from the internal register or internal cache and decode them. During or after execution of the instructions, the processor 106 may write one or more results (which may be intermediate or final results) to the internal register or internal cache. The processor 106 may then write one or more of those results to the memory 104. In particular embodiments, the processor 106 executes only instructions in one or more internal registers or internal caches or in memory 104 (as opposed to storage 108 or elsewhere) and operates only on data in one or more internal registers or internal caches or in memory 104 (as opposed to storage 108 or elsewhere). One or more memory buses (which may each include an address bus and a data bus) may couple the processor 106 to the memory 104. The bus may include one or more memory buses, as described in further detail below. In particular embodiments, one or more memory management units (MMUs) reside between the processor 106 and memory 104 and facilitate accesses to the memory 104 requested by the processor 106. In particular embodiments, the memory 104 includes random access memory (RAM). This RAM may be volatile memory, where appropriate. Where appropriate, this RAM may be dynamic RAM (DRAM) or static RAM (SRAM). Moreover, where appropriate, this RAM may be single-ported or multi-ported RAM. This disclosure contemplates any suitable RAM. Memory 104 may include one or more memories 104, where appropriate. Although this disclosure describes and illustrates particular memory implementations, this disclosure contemplates any suitable memory implementation.
In particular embodiments, the storage 108 includes mass storage for data or instructions. As an example and not by way of limitation, the storage 108 may include a hard disk drive (HDD), a floppy disk drive, flash memory, an optical disc, a magneto-optical disc, magnetic tape, or a Universal Serial Bus (USB) drive or a combination of two or more of these. The storage 108 may include removable or non-removable (or fixed) media, where appropriate. The storage 108 may be internal or external to computer system 100, where appropriate. In particular embodiments, the storage 108 is non-volatile, solid-state memory. In particular embodiments, the storage 108 includes read-only memory (ROM). Where appropriate, this ROM may be mask-programmed ROM, programmable ROM (PROM), erasable PROM (EPROM), electrically erasable PROM (EEPROM), electrically alterable ROM (EAROM), or flash memory or a combination of two or more of these. This disclosure contemplates mass storage 108 taking any suitable physical form. The storage 108 may include one or more storage control units facilitating communication between processor 106 and storage 108, where appropriate. Where appropriate, the storage 108 may include one or more storages 108. Although this disclosure describes and illustrates particular storage, this disclosure contemplates any suitable storage.
In particular embodiments, the I/O Interface 110 includes hardware, software, or both, providing one or more interfaces for communication between computer system 100 and one or more I/O devices. The computer system 100 may include one or more of these I/O devices, where appropriate. One or more of these I/O devices may enable communication between a person (i.e., a user) and computer system 100. As an example and not by way of limitation, an I/O device may include a keyboard, keypad, microphone, monitor, screen, display panel, mouse, printer, scanner, speaker, still camera, stylus, tablet, touch screen, trackball, video camera, another suitable I/O device or a combination of two or more of these. An I/O device may include one or more sensors. Where appropriate, the I/O Interface 110 may include one or more device or software drivers enabling processor 106 to drive one or more of these I/O devices. The I/O interface 110 may include one or more I/O interfaces 110, where appropriate. Although this disclosure describes and illustrates a particular I/O interface, this disclosure contemplates any suitable I/O interface or combination of I/O interfaces.
In particular embodiments, communication interface 112 includes hardware, software, or both providing one or more interfaces for communication (such as, for example, packet-based communication) between computer system 100 and one or more other computer systems 100 or one or more networks 114. As an example and not by way of limitation, communication interface 112 may include a network interface controller (NIC) or network adapter for communicating with an Ethernet or any other wire-based network or a wireless NIC (WNIC) or wireless adapter for communicating with a wireless network, such as a Wi-Fi network. This disclosure contemplates any suitable network 114 and any suitable communication interface 112 for the network 114. As an example and not by way of limitation, the network 114 may include one or more of an ad hoc network, a personal area network (PAN), a local area network (LAN), a wide area network (WAN), a metropolitan area network (MAN), or one or more portions of the Internet or a combination of two or more of these. One or more portions of one or more of these networks may be wired or wireless. As an example, computer system 100 may communicate with a wireless PAN (WPAN) (such as, for example, a Bluetooth® WPAN), a WI-FI network, a WI-MAX network, a cellular telephone network (such as, for example, a Global System for Mobile Communications (GSM) network), or any other suitable wireless network or a combination of two or more of these. Computer system 100 may include any suitable communication interface 112 for any of these networks, where appropriate. Communication interface 112 may include one or more communication interfaces 112, where appropriate. Although this disclosure describes and illustrates a particular communication interface implementations, this disclosure contemplates any suitable communication interface implementation.
The computer system 102 may also include a bus. The bus may include hardware, software, or both and may communicatively couple the components of the computer system 100 to each other. As an example and not by way of limitation, the bus may include an Accelerated Graphics Port (AGP) or any other graphics bus, an Enhanced Industry Standard Architecture (EISA) bus, a front-side bus (FSB), a HYPERTRANSPORT (HT) interconnect, an Industry Standard Architecture (ISA) bus, an INFINIBAND interconnect, a low-PIN-count (LPC) bus, a memory bus, a Micro Channel Architecture (MCA) bus, a Peripheral Component Interconnect (PCI) bus, a PCI-Express (PCIe) bus, a serial advanced technology attachment (SATA) bus, a Video Electronics Standards Association local bus (VLB), or another suitable bus or a combination of two or more of these buses. The bus may include one or more buses, where appropriate. Although this disclosure describes and illustrates a particular bus, this disclosure contemplates any suitable bus or interconnect.
The development of closed-loop electrical neuromodulation therapies, for example adaptive Deep Brain Stimulation (aDBS) and adaptive Spinal Cord Stimulation (aSCS) can revolutionize the efficacy of neurostimulation therapies for treatment of many disorders, including Parkinson's Disease (PD), Epilepsy, Essential Tremor, Obsessive Compulsive Disorder (OCD), Treatment Resistant Depression (TRD), and chronic pain. While significant advances have been made in biomarker identification in PD and Epilepsy, there may not definitive biomarker for a single mental disorder, including OCD and TRD, or chronic pain. Biomarker identification and development of an adaptive neurostimulation system may require a hardware platform capable of simultaneous sensing and stimulation. This is particularly challenging when the neural signal of interest originates in or nearby the stimulation target as the amplitude of stimulation therapy is typically several orders of magnitude greater than the amplitude of signals of interest in the brain and spinal cord. Therefore, recordings for adaptive control may be heavily contaminated by high amplitude, high frequency stimulation artifact. In order to extract the underlying neural signatures of disease state, it may be necessary to remove the stimulation artifact.
Typically, high frequency artifacts are removed using a lowpass filter, however, limited sampling rates of existing implantable DBS and SCS devices and aliasing of stimulation pulses into low frequencies render lowpass filters ineffective. Stimulation artifact removal methods robust to aliasing may include: signal reconstruction via deletion and interpolation, decomposing and subtracting components of the signal related to the artifact, and subtracting a template of the artifact at each stimulation pulse. Deletion and interpolation may not be effective when artifact duration is long, and may not be ideal due to signal loss over the duration of each artifact. Signal decomposition methods have shown some success, but may require a large set of recording channels to be effective. While template subtraction methods have proven to be successful, such methods rely on accurate detection of each stimulation pulse. Existing methods for identifying individual stimulation pulses in recorded data (e.g. thresholding) are not robust to low sampling rates, the presence of other spurious high amplitude artifacts, or stimulation artifacts with broad peaks. Thus, there is a desire and need for a system and method that is effective at removing stimulation artifacts from local field potential (LFP) and other electrical stimulation recordings sampled at less than twice the frequency of stimulation without contaminating the underlying neural signal. The absence of such systems and methods thus greatly hinders control signal identification.
To overcome the challenges in removing periodic stimulation from neural recordings, various embodiments of the present disclosure describe a novel artifact removal method (e.g., a Period-Based Artifact Reconstruction and Removal Method (PARRM), to remove high frequency stimulation artifact in low and high-resolution LFP and other electrical stimulation recordings. As will be discussed herein, PARRM has superior performance to existing, state-of-the-art filters in saline experiments, computer simulations, and five unique in vivo recording paradigms. As will be discussed herein, PARMM can allow the recovery of a previously obscured biomarker in Parkinson's Disease participants and could be implemented online to perform real-time biomarker detection (e.g., as illustrated in
An example design of PARMM may involve subtracting an estimate of the stimulation artifact at each time bin from the recorded signal at that time bin. The artifact estimate may be formed by averaging the recorded signal at other time bins that are close to the current time bin in both time and stimulation phase. The artifact may bepresumed to be roughly identical for all of these time bins. Averaging may reduce the influence of brain signals and additional sources of noise, so that the estimate is primarily artifact. This process can be implemented as a linear filter (e.g., asaa weighted average using a sliding window). In some embodiments, PARRM may need a precise estimate of the stimulation period relative to the sampling rate. Slight inaccuracies in device system clocks can necessitate using a data-driven method to determine this period, which is done by finding the period that, when the data are divided into epochs the length of one period and overlapped, the samples will consolidate around the shape of the high-resolution artifact waveform. The complete process of data-driven period finding, artifact estimation, and signal reconstruction is illustrated in
Method 300 may begin with the application of a deep brain stimulation at a patient-specific area of interest (block 302). For example, as previously discussed in
The computing device may thus receive waveform data causes by the deep brain stimulation (block 304). In some aspects, the waveform data may be received, and one or more subsequent steps, processes, or blocks may be performed, in real-time or near real-time. At block 306, the computing device may determine a stimulation period relative to a sampling rate. As will be discussed, method 400 shown in
The computing device may identify, a stimulation artifact in the waveform data (block 310). In some aspects, the stimulation artifact may be identified by averaging recorded image data of the waveform data at time bins proximal to the time and/or proximal to the simulation phase of time bin (t) of the waveform data (block 308). The artifact may be presumed to be roughly identical for all of these time bins, including time bin t. Averaging may reduce the influence of brain signals and additional sources of noise, so that the subtracted signal may be deemed as primarily stimulation artifact. After the stimulation artifact is identified, the computing device may remove the stimulation artifact from the waveform data. Thus, for each time bin (t), the computing device may subtract the identified stimulation artifact at time bin (t) from the received waveform data at time bin (t) (e.g., the recorded signal at time bin (t)) (block 312).
For example, in some embodiments, blocks 308 through 314 may be shown through the following operations. Let T denote the stimulation period relative to the sampling rate (in units of sampling time bins) (e.g., as determined in block 306). The time bins included in the average (e.g., determined in block 308) may be those times bins s such that Nskip<|s−t|≤Nbins and such that |s−t| (mod T)≤Dperiod or |s−t| (mod T)≤T−dperiod, where α (mod T) denotes α modulo T, and where 0≤Nskip<Nbins and 0≤Dperiod≤T are user-chosen design parameters. In some aspects, the additional criterion s−t<0 can be included so that only past observations are used to estimate the stimulation artifact. Let Bt denote the collection of those times bins s that are used for averaging and let |Bt| denote the number of such time bins. Using rt to denote the recorded signal at time bin t, the corrected signal is ct defined by
where wi is a list of weights defined by w0=1, and wi=−1/B0| if −i ϵ B0, and wi=0 otherwise. The final expression shows that the PARRM correction can be implemented by a fixed linear filter (with the filter weights denoted by w′), making it fast and simple to implement. (e.g., if the additional criterion s−t<0 is used, then the final summation would begin at i=0. The final methodology may be modified near the start and end of stimulation.
In some embodiments, the design parameters for the PARRM filter are Nbins, Nskip, and Dperiod. Larger choices of Nbins allow more data to be averaged in order to estimate the artifact, reducing estimation variability. But, larger choices of Nbins also lengthen the temporal window used to estimate the artifact, perhaps introducing estimation bias if the artifact shape is changing in time. Because neural signals have temporal autocorrelation, it is important to avoid averaging data too close to time bin t or the neural signal itself could be subtracted during artifact removal. Larger choices of Nskip help to mitigate this danger, but also reduce the amount of data used to estimate the artifact. Similar to Nbins, larger choices of Dperiod allow more data to be averaged, but also introduce more estimation bias by temporally smoothing the estimated artifact. The optimal choices for these design parameters may vary depending on the situation.
Method 400 is an automated, data-driven method to determine the simulation period, T, by searching for a period that creates a strongly resolved template. Method 400 begins with selecting a candidate period, δ>0 (block 402). For each candidate period δ>0, the method estimates a waveform template with this period (block 404). The computing device may thus quantify a deviation from the estimated template (block 406). For example, for each candidate period, fB,δwhere the potential period is δ>0, and β=(β=, . . . , β>2m+1) is a parameter vector, the computing device may determine whether a deviation between the candidate period, fB,δand an estimated waveform template (e.g., its amplitude, yk) is the smallest (block 408). If not the smallest, another candidate period may be selected. The candidate period with the smallest deviation may be identified as the final estimate T of the period that is used by PARRM (block 410).
In at least one embodiment, the aforementioned blocks may be described by the following operations. Let m≥0 be an integer. For each potential period δ>0 and each parameter vector β=(β=, . . . , β>2m+1) define the functions
The function fB,δis a periodic function with period δ. Each fB,δis a candidate artifact waveform. The parameter vector β controls the strength of the different frequencies that define fB,δ, and m controls the number of allowed frequencies. Let ((tk, yk): k=1, . . . , n) be a collection of (time, value) pairs. The yk value used here is the change in recorded electrical stimulation (e.g., LFP) amplitude at time tk with some preprocessing to obtain standardized units, reduce the influence of outliers, and reduce the size of the dataset. Mean squared error is used to measure how well the function fB,δ, fits these pairs:
In some embodiments, the minimization over δ may complicated by many local minima, spurious ‘distractor’ solutions that mimic the harmonics of the true waveform, and high sensitivity to small changes in δ. In some aspects, a penalized, stagewise search that begins with smaller intervals of data (to reduce the sensitivity to δ), smaller m (to reduce the number of local minima), and a penalty for higher frequency solutions (to help avoid distractor solutions) may be used to overcome these issues. After the stimulation period, T is found, it may be fixed for PARRM (e.g., as previously discussed in method 300). Furthermore, as will be discussed in
The similarity between the artifact free and filtered signals is visually apparent in 0.2 seconds of data for both the 10 and 50 Hz injected signals (e.g., as shown in graphs 614 and 616). Afterwards, filter performance can be quantified by comparing the distribution of absolute errors between artifact free signals and unfiltered, MAS filtered, notch filtered, and PARRM filtered signals to a baseline noise recording (e.g., no stimulation and no injected signal) (e.g., as shown in graphs 618 and 620). Filtering using MAS may not reduce the error to the level of baseline for either injected signal (p<0.0005). While effective at 10 Hz (p>0.05), notch filtering removed the injected signal along with the artifact leading to a large reduction in error yet still significantly different from baseline (p<0.0005). For both the 10 and 50 Hz injected signals, PARRM outperformed the other methods with no significant difference (p>0.05) from baseline, indicating that the remaining errors were expected due to noise in saline.
Having shown that PARRM is effective for recovering simple sinusoidal signals recorded in saline, a comparison of the method's performance to conventional filters in recovering more complex, injected, chirp signals for simulated data is discussed.
Next, a parameter sweep was performed to test the effect of varying chirp length (1-10 s), amplitude (0.5-5 V), pulse width (30-180 μV), and frequency (80-180 Hz) on PARRM performance, measured by RRMSE and Relative R Ratio (S.
As waveform data caused by deep brain stimulation is received by the computing device 100, it can be expected that portions of the original waveform data (e.g., packets) may be lost in the transmission process. For example,
For example,
In some embodiments, before subsequent steps, processes, or blocks of PELP can be applied to the received recording of waveform data, the locations of packet losses and their sizes may initially be determined (block 904). In some aspects (e.g., in recordings by the SUMMIT RC+S device), each packet may have has three integer timing variables of note for this purpose: (1) ‘dataTypeSequence’, which may indicate the packet number that rolls over every 256 packets, (2) ‘systemTick’ time of the last sample in a packet with 0.1 ms resolution that rolls over every 6.5536 s, and (3) ‘timestamp’ with 1 s resolution and no rollover. A packet loss may occur when the dataTypeSequence between subsequent packets skips an index or the timestamps is more than 6 seconds apart corresponding to cases of large losses where the dataTypeSequence may roll over. For sampling rate fs, the number of samples lost, N, can be estimated according to the following equation when there is no systemTick rollover: N=((s2−s1−104 fs−1 (n−1)) mod 65536)×10−4, here s1 and s2 are the system ticks of the packets preceding and following the loss respectively and n is the number of samples in the packet after the loss. If systemTick rollover occurs, estimation may no longer be acceptably accurate due to timing drift between the systemTick and timestamp, and may need to be reset using a coarser metric, such as the UNIX timestamp (±50 ms vs. ±3 ms) corresponding to the time when the packet was received or generated. These estimates may not be sufficiently accurate to ensure the exact reconstruction of the timing between received packets down to sample resolution. PELP leverages the presence of regular stimulation in both the received and missing data to ensure precise estimates of data losses.
PELP may proceed by dividing the time series into a set of R continuous ‘runs’ (block 906). Each run may comprise a continuguous packet that is separated from another run by losses (block 906). Thus, computing device 100 may recognize runs (e.g., based on non-contiguos nature of the received waveform) and identify such runs accordingly.
At block 908, the computing device may determine, for each run, a period of stimulation (6) (e.g., from the waveform data of each respective run) relative to the sampling rate. The computing device 100 may use the period estimation component of PARRM, as previously discussed (e.g., block 304 of
A harmonic regression model, fβ,δ, is then fit to the longest run to estimate the coefficients (β) and the optimal number of harmonics (m) (block 910). The optimal number of harmonics (m) may be chosen using Akaike Information Criterion (AIC). In some embodiments, additional harmonic regression models may be fit, e.g., for other runs of the received waveform data. The computing device 100 may thus use the harmonic regression model, fβ,δ, to determine the optimal size of the packets loss in the received waveform data (block 912). An example process for performing blocks 910 and 912 (e.g., training and applying a harmonic regression model and determining the optimal size of packet loss) is explained in more detail in method 1000 shown in
Method 1000 may begin with the computing device receiving a plurality of runs, including the longest run (block 1002), of the received waveform from the deep brain stimulation of method 900. At block 1004, the computing device may generate a candidate harmonics function, fB,δ(t) for the training process. An example candidate harmonics function, fB,δ(t), and the variables and coefficients to be determined for the candidate harmonics function is explained below:
Let ((tk,r,yk,r): k=1, . . . , nr: r=1, . . . , R) be a collection of (time, value) pairs.
At block 1006, the computing device may determine a root mean squared error (RMSE) between the waveform corresponding to the run (e.g., the longest run) and a candidate harmonics function. In some aspects, the waveform corresponding to the run may be represented by a function that outputs the waveform's amplitude. For example, a value, yk,r may be used to represent the change in recorded electrical stimulation (e.g., LFP) amplitude of a run, r, at a time relative to the start of the recording tk,r for run r. Root mean squared error (RMSE) may be used to measure how well the candidate harmonics function fB,δ(t) fits these pairs for a packet loss size Δ:
Since, an objective may be to find a candidate harmonics function for which the RMSE between the candidate harmonics function and the waveform is at a minimum, t the optimal size of the packets loss for run r, Δr, may therefore be: Δr=argminΔ(rmse(Δ,r)).
As different candidate harmonic functions are tested, the computing device may determine, at each iteration, whether a given harmonic function satisfies a threshold (block 1008). For example, the threshold may be based on whether a previously held record of a minimum RMSE has been beaten.
After an RMSE satisfies a threshold (e.g., a minimum RMSE is reached between a candidate harmonics function and the waveform corresponding to the run), the computing device 100 may identify, at block 1010 the loss size associated with the run (run-specific loss size). Furthermore, the computing device may use the candidate harmonics function (for which the RMSE satisfied the threshold) as the harmonics function model for the run (e.g., at block 1012).
Blocks 1002 through 1012 were thus explained for determining a harmonics function model based on the longest run. The harmonics model may be applied to other runs of the waveform data received in method 900. In some embodiments, other harmonics model may be determined for other runs of the waveform data. Thus, if the computing device determines that there are additional runs (block 1014—Yes), the computing device 100 may receive, at block 1016, the waveform for the longest run of the additional runs (i.e., the next longest run of the entire waveform received in method 300). Blocks 1006 through 1014 may thus be repeated for each next longest run until all runs of the entire waveform data has been assessed. At each iteration, the run-specific loss sizes may be identified and stored.
At block 1018, the computing device 100 may aggregate run-specific loss sizes (e.g., identified at block 1010 during one or more iterations) to determine the size of the packets loss in the received waveform data. The method of estimating packet loss is further illustrated in
Various experiments were performed using waveforms received from deep brain stimulation that show the effectiveness of PELP as a method of determining packet loss. In at least one experiment, neural data was recorded from one participant undergoing DBS surgery for treatment of obsessive compulsive disorder (OCD). Electroencephalography (EEG) data were recorded both with and without stimulation when the participant visited the clinic for DBS programming. LFP data were recorded both in the clinic and when the patient was at home. Electrodes on both the right and left sides were implanted in the VC/VS according to standard stereotactic procedures using computed tomography for target determination. The location of electrode placement was made entirely on clinical grounds. Bilateral 150.6 Hz stimulation with a pulse width of 90 ps and amplitude of 4 mA for the left side and 4.5 mA for the right side was used for all recordings where stimulation was turned on. The participant gave informed consent and data presented were collected in accordance with recommendations of the federal human subjects regulations and under protocol H-40255 approved by the Baylor College of Medicine Institutional Review Board.
Continuous electroencephalography (EEG) was recorded using a 64-channel ACTICAP BRAIN VISION SYSTEM (BRAIN VISION, Morrisville, N.C., USA). A common mode sense electrode was located at FCz. The EEG was band-pass filtered online between 0.1 and 1000 Hz and digitized at 5 kHz. The EEG was downsampled offline to 1000 Hz with an anti-aliasing filter prior to analysis. The continuous LFP was recorded using the SUMMIT RC+S (MEDTRONIC, Minneapolis, Minn., USA) via wireless data streaming from implanted electrodes to the device running the task. Each DBS probe (Model 3387, one per hemisphere) included four electrode contacts two of which were used per side to conduct bipolar recordings. LFP recordings were sampled at 1 kHz in the clinic and 250 Hz at home in order to minimize data losses. Signal processing and analysis were performed in MATLAB (MATHWORKS, Natick, Mass., USA) using in-house code.
In order to simulate stimulation in the DBS recordings, DBS artifacts were modeled as a sum of sinusoidal harmonics of the stimulation frequency (e.g., as previously discussed in
Furthermore, three sets of experiments were performed to simulate the accuracy of loss estimation while varying different parameters in the stimulation model. For each experiment, Monte Carlo analyses were used to simulate a large number of experiments by randomly sampling subsets of 50 sample “packets” to remove from the recording. These simulations were applied while varying one of amplitude ratio, amplitude variability, or drift as a function of the loss uncertainty. For each simulation, the uncertainty ranged from 0-50 samples in 1 sample increments while the dependent parameters ranged from 0-4 in increments of 0.1 for the amplitude ratio, 0-10% in increments of 1% for the amplitude variability, and 0-0.6% in increments of 0.015% per 1000 cycles for the drift. PELP was applied to each simulation to determine the proportion of losses it was able to estimate correctly depending on the stimulation parameters.
The histograms 1304-1308 in
Thus, streaming of intracranial electrophysiology data in the clinic and at home in ecologically valid environments may be essential for biomarker discovery in a variety of neurological disorders. Bidirectional implanted devices have enabled the acquisition of such datasets. However, data losses during wireless streaming hinder accurate analyses of neural signals. To address these issues, various aforementioned embodiments of the present disclosure describe novel and nonobvious systems and methods (e.g., PELP) to precisely estimate and account for data losses from implanted recordings where stimulation is on.
In some embodiments, PELP may be widely applicable to other stimulating devices capable of wireless data streaming. The stimulation model disclosed herein accurately accounts for the range of amplitude and variability parameters that could be expected for other implanted devices. For example, in recordings where sensing and stimulation occur on nearby contacts, stimulation amplitude can exceed the underlying neural signal by a factor of 10. For recordings where sensing and stimulation occur far apart or the stimulation harmonics fall within the transition band of an online low-pass filter, the amplitude ratio will be closer to 1.
Our recordings using the SUMMIT RC+S had amplitude ratios of 27 and 1.2 and estimate accuracy was consistent with predictions from the simulation. Pulse to pulse amplitude variability for the SUMMIT RC+S is well within the range of values where PELP was most accurate. Fluctuations in battery or the surrounding medium could influence amplitude on longer timescales. Instead of using the longest run for modeling stimulation in the entire recording, neighboring runs could be used for each prediction to account for slower drift in parameters or longer recordings. Similar considerations would also be effective if stimulation frequency drift or errors in period estimation occur.
Since PELP requires stimulation artifacts to be present in order to model the signal during data losses, the method may not necessarily be applicable for recordings where stimulation is off or significantly attenuated by online filters. In such circumstances, PELP may be complemented or replaced with less accurate methods utilizing packet timing metadata for loss estimation. In theory, stimulation could be applied below therapeutic amplitudes and can still be used for reconstruction using PELP although such modifications would be reasonable if the inevitable stimulation artifacts did not obscure neural signals of interest.
In some embodiments, PELP could be adapted to determine phase differences between discontinuous bursts of stimulation by minimizing mean squared error as a function of the stimulation phase offset between bursts. Such an approach could be useful for increasing the amount of data available to artifact removal methods relying on template subtraction.
By leveraging the periodicity of stimulation artifacts, PELP can produce highly accurate reconstructions of timing information from wirelessly transmitted neural data. PELP effectively accounts for the timing of missing data enabling the analysis of complex, naturalistic neural time series data from next-generation bidirectional implanted devices aiding in the development of novel therapeutic approaches.
As previously discussed, various embodiments of the present disclosure describe systems and methods for period-based estimation of electrical stimulation artifacts. Such systems and methods may include the determination of a stimulation period of the artifact.
In some embodiments, estimation of electrical stimulation artifacts from waveform data may be made difficult or unfeasible due to a variety of factors. Such factors may include but are not limited to missing data of unknown duration (e.g., due to packet losses in wireless transmission combined with inaccuracies in device timing capabilities), device modulation between time segments of constant stimulation (e.g., caused by trials of interleaved stimulation-on and stimulation-off, and/or by the adjustment of device settings for clinical purposes), and device irregularities where the period of stimulation and/or sampling is momentarily different. Such factors may cause one or more phase shifts in the waveform data, further confounding the ability to estimate electrical stimulation artifacts from waveform data. Various embodiments of systems and methods for period-based estimation of electrical stimulation artifacts in the presence of phase shifts are thus described below. In some embodiments, such systems and methods may estimate these multiple phase shifts simultaneously with the estimation of the stimulation artifacts. Furthermore, such systems and methods may also or alternatively estimate the stimulation period of artifacts simultaneously with the estimation of the phase shifts from the waveform data.
Method 1600 may begin with receiving waveform data caused by deep neural stimulation of a patient-specific area of interest (block 1602). For example, the patient-specific area of interest may comprise an area or section of a brain of a patient being treated, on which electrodes may be applied. The received waveform data may be characterized by a plurality of runs (e.g., segments). Each run (e.g., segment) may comprise a contiguous portion of the waveform data separated by a packet loss. Furthermore, the received waveform data may include (e.g., in addition to the underlying neural signal of interest), one or more periodic artifacts. Such artifacts may be estimated (e.g., based on their stimulation periods) using methods presented herein. Furthermore, as will be discussed in method 1600, the waveform data may be characterized as having multiple phase shifts. In some aspects, the waveform data may be received in real-time or near real-time. The received waveform data (e.g., observed signal) may be incorporated into a model (model of the received waveform data) in order to estimate and determine a periodic artifact, as will be discussed herein. An example of the model of the received waveform data, Si, may be as follows. In at least one embodiment, given n+1 segments of data, an i-th run (e.g., i-th segment) of the received waveform data (e.g., observed signal), Si, may be modeled as:
where δi* represents the phase shift between the 0th and i-th run of the received waveform data, A represents a model of the periodic artifact based on a stimulation frequency ξ* (or stimulation period
of the periodic artifact and the phase shift (δi*) of the plurality of phase shifts, Bi(t) represents the neural signal of interest at the i-th run, and ηi(t) represents noise at the i-th run. Since, the periodic artifact and the phase shifts may not be known when the waveform data is received, initial estimates for the periodic artifact and phase shifts may be generated (block 1604) as part of an initialization process described further below.
An objective of method 1600 may be to estimate and remove the model of the periodic artifact, A, from the model of the received waveform data at each run, Si. As a result, the underlying signal of interest that may be desired to be recovered may include: Bi+ηi, where i=0, . . . , n. In some embodiments, Si may be sampled at some collection of discrete times Ti. Then, the following loss function can be used to reconstruct and remove the artifact: L(ξ,δi,θ)=Σi=0nΣtϵT
of the periodic artifact, a phase shift δi of the plurality of phase shifts of the received waveform data, and a set (possibly empty) of parameters θ for the periodic artifact, and where Si(t) is the model of the received waveform data at time t, and sampled at discrete times, Ti. As shown in the loss function model, the model of the periodic artifact may be need to initially be determined. In some embodiments, the periodic artifact may be modeled by optimizing a parametric model of the periodic artifact using harmonic regression. In some aspects, a(·|ξ,δi,θ) approximates, represents, and/or is approximated by the previously discussed
Thus, at block 1606, a periodic artifact model may be optimized (e.g., via harmonic regression) to estimate the periodic artifact and phase shifts from the received waveform data. For example, the following parametric equation may be used for periodic artifact model: a(t|ξ,δ,α0,αk,βk,K)=α0+Σk=1K(αk cos(2πk(ξt+δ)+βk sin(2πk(ξt+δ))), where K may represent a specific number of harmonics, k. In some embodiments, the number of harmonics, K, applied during the harmonic regression must be chosen to be finite. Since the periodic artifacts model takes phase shifts into account, an optimization of the periodic artifacts model to determine or estimate the periodic artifact may also simultaneously determine or estimate the multiple phase shifts exhibited in the waveform data. Using the above described parametric model may allow for fast computations, thereby freeing up computing resources of computing device 100. Based on the model of the received waveform data, Si, and the aforementioned periodic artifact model to be optimized via harmonic regression, a loss model may be generated (block 1608). The loss function model (e.g., L(ξ,δi,θ)=Σi=0nΣtϵT
As previously discussed, in relation to
In some embodiments, the aforementioned parametric model for the periodic artifact, a(t|ξ,δ,α0,αk,βk,K) (“periodic artifact model”) may be used to remove the periodic artifact from the received waveform data (modeled as Si) using harmonic regression. For example, harmonic regression may be applied to optimize the periodic artifact model a(t|ξ,δ,α0,αk,βk,K) to estimate one or more periodic artifacts from the received waveform data (block 1606). For fixed frequencies of artifacts and phase shifts (ξ,δ), the harmonic regression may be equivalent to minimizing the loss function model, L(ξ,δi,θ) with respect to θ={α0,αk,βk,k=1, . . . , K}. Thereafter, the previously discussed loss function model, L(ξ,δi,θ)=Σi=0nΣtϵT
relative RMSE=√{square root over (Σi=1N(f(ti)−{circumflex over (f)}(ti))2/Σi=1Nf(ti)2)}, where f is the true signal of interest, {circumflex over (f)} is the estimated signal, and ti are the sample times. It is to be appreciated that even somewhat small errors in the frequency may lead to large reconstruction errors.
Many periodic artifact removal methods rely on the discrete Fourier transform (DFT) to estimate the frequency. However, DFT-based methods generally cannot handle phase shifts and, even segment-wise, are limited in their accuracy for frequency estimation. For example, consider the following simple DFT-based frequency detection method, where the artifact frequency is estimated as the frequency that maximizes the energy (magnitude squared of the DFT) of the observed signal. Graph 1704 of
which may be the case for most real-life applications. Specifically, even with 105 samples, this method may achieve, on average, greater than 0:001% relative error in the frequency estimate, which, as discussed above, is insufficiently accurate for harmonic regression.
In contrast, graph 1704 also shows that the method discussed herein (e.g., in
In some embodiments, the iterative periodic artifact removal methods and systems described herein can accurately estimate the artifact frequency in the presence of unknown phase shifts and may involve few tuning parameters and relatively simple computations. An example system and method is further described below. First, the artifact removal process is further described, which removes the artifact by estimating the frequency and phase shifts while jointly fitting the artifact using harmonic regression. Second, an initialization process is further described, which finds an initial estimate of the frequency and phase shifts via an energy maximization method.
For n+1 runs (e.g., segments of waveform data Si), where i=0; . . . , n, a model for the waveform data, Si, may be defined as follows:
(e.g., as previously discussed), and δi* may represent the (unknown) true phase shift between the 0-th and i-th runs (e.g., segments) (by convention, δ0*=0). The artifact may be assumed to be represented by a(t|ξ,δ,α0,αk,βk,K) (e.g., as previously discussed), where the true fundamental frequency ξ*, the true amplitudes α0*,αk*,βk*, and the true number of harmonics K are unknown. It may be assumed for this process that the energy of the artifact at its fundamental frequency and harmonics can be sufficiently larger than the energy of the brain at those frequencies (or at the frequencies that they are aliased to). In some embodiments, there may not be access to the continuous signals and there may only be access to each Si at some discrete sample times. While these sample times may be irregularly spaced, the process considers only uniformly spaced sample times as most DBS devices typically use uniformly spaced sampling based on relatively stable system clocks. Thus, it is assumed that the i-th segment is sampled at the times
where fs is the sampling rate.
Based on the assumed form of the artifact a(t|ξ,δ,α0,αk,βk,K), harmonic regression may be is a natural choice for removing the artifact given an estimate of the frequency. However, as shown in
where the objective function g is defined by:
where the i-th segment of the recovered signal can be modeled as: Ri(t, α0,αk,βk)=Si(t)−a(t|ω, δi, α0, αk, βk, {circumflex over (K)}), where a(t|ω, δi, α0, αk, βk, {circumflex over (K)}) is defined by a(t|ξ,δ,α0,αk,βk,K). For each fixed (ω, δi, . . . , δn), g(ω, δi, . . . , δn) is the result of using harmonic regression to remove a periodic artifact of the form a(t|ξ,δ,α0,αk,βk,K) with K={circumflex over (K)} harmonics, fundamental frequency ξ*=ω, and phase shifts δi from the 0-th observed segment S0. In some embodiments, the number of harmonics to fit {circumflex over (K)} may effectively be the only parameter that needs to be tuned. Minimizing α0 in
may mean that the recovered signal will have a mean of 0. Including this amplitude may be optional, and if α0 is omitted in Ri(t,α0,αk,βk), then the recovered signal may have a nonzero mean, but the reconstructed artifact may necessarily have a mean of 0. Furthermore, each fixed (ω, δi, . . . , δn), g(ωm δi, . . . , δn) can then be computed up to numerical precision errors. Specifically, the gradient of the objective function in
may be set to zero to solve the resulting linear system using an appropriate numerical linear algebra solver.
For example, Newton's descent method may be used to solve for the optimization problem
in some aspects, the numerical complexity of the above described artifact removal process may be dominated by the computation of the descent direction, which may requires solving for an (n+1)×(n+1) and a {circumflex over (K)}×{circumflex over (K)} linear system. Thus, the numerical complexity of the algorithm is O(n3+{circumflex over (K)}3). In practice, both n and {circumflex over (K)} may be relatively small since the data can always be segmented so as not to consider too many gaps at a time and many DBS devices have built-in low-pass filters that effectively band limit the recorded signals. Newton's descent can also be applied to this problem since g can be twice differentiable. Also while g may not be a convex function, g may be locally convex at the global minimizer. Hence, the initialization for the artifact removal process may be crucial for convergence. An initialization process for the artifact removal process is thus described. Based on numerical experiments, the initialization process may typically have relatively short computational runtime and may provide sufficiently accurate initializations for the artifact removal process.
The initialization process locally maximizes the energy of the observed signal using discrete samples but without being constrained to frequencies on the DFT grid. A Fourier transform for the initialization process may be defined as follows:
The above Fourier transform aligns all of the segments Sl (l=0, 1, . . . , n) in time to be in phase with a single wave e−2πiωt. Furthermore, the energy may be defined as:
E(ω,δ1, . . . ,δn=|(ω,δ1, . . . ,δn)|2,
By straightforward calculation, one can easily verify that E is twice differentiable and nonconcave. The True frequency ξ* and true phase shifts δi* may be recovered by solving the following optimization problem:
Note that for fixed ω,δj,j≠i, the function δi→(ω,δ1, . . . ,δn) (and hence, the function δi→E(ω, δ1, . . . , δn)) is 1-periodic. Thus, solving
may only recover the phase shifts i up to a multiple of 1. Moreover, since E is nonconcave, the maximizers of (11) may not necessarily be unique.
In the initialization process,
can be solved numerically by applying a modified Newton's ascent method with backtracking line search and random initialization.
Since there may only be access to discrete samples (and not the continuous signals), the integrals in the definitions of the energy E, its gradient ∇E, and its Hessian ∇E2 may be approximated using trapezoidal rule. The computational effort for each iteration of the initialization process may be proportional to (n+1)3. Typically, the data can always be segmented so that n is relatively small. Additionally, the implementation of the initialization process can be made more efficient by parallelizing the computations for each random initialization.
Thus, a periodic artifact removal process (and a corresponding initialization process) has been described herein that is able to accurately recover underlying signals in the presence of unknown phase shifts (e.g., missing data of unknown length), and in some cases where the stimulation frequency exceeds the Nyquist frequency. The disclosed processes effectively require only one tuning parameter—the number of harmonics {circumflex over (K)} to fit. Choosing k too small may cause an underfitting of the artifact, while choosing a {circumflex over (K)} that is too large may cause an overfitting of the artifact. Thus, the parameter {circumflex over (K)} can be easily interpreted and could be roughly estimated by counting the number of bands corresponding to the harmonics of the artifact that are visible in a spectrogram of the observed signal.
Herein, a computer-readable non-transitory storage medium or media may include one or more semiconductor-based or other types of integrated circuits (ICs) (e.g., field-programmable gate arrays (FPGAs) or application-specific ICs (ASICs)), hard disk drives (HDDs), hybrid hard drives (HHDs), optical discs, optical disc drives (ODDs), magneto-optical discs, magneto-optical drives, floppy diskettes, floppy disk drives (FDDs), magnetic tapes, solid-state drives (SSDs), RAM-drives, SECURE DIGITAL cards or drives, any other suitable computer-readable non-transitory storage media, or any suitable combination of two or more of these, where appropriate. A computer-readable non-transitory storage medium may be volatile, non-volatile, or a combination of volatile and non-volatile, where appropriate.
Herein, “or” is inclusive and not exclusive, unless expressly indicated otherwise or indicated otherwise by context. Therefore, herein, “A or B” means “A, B, or both,” unless expressly indicated otherwise or indicated otherwise by context. Moreover, “and” is both joint and several, unless expressly indicated otherwise or indicated otherwise by context. Therefore, herein, “A and B” means “A and B, jointly or severally,” unless expressly indicated otherwise or indicated otherwise by context.
The scope of this disclosure encompasses all changes, substitutions, variations, alterations, and modifications to the example embodiments described or illustrated herein that a person having ordinary skill in the art would comprehend. The scope of this disclosure is not limited to the example embodiments described or illustrated herein. Moreover, although this disclosure describes and illustrates respective embodiments herein as including particular components, elements, features, functions, operations, blocks, or steps, any of these embodiments may include any combination or permutation of any of the components, elements, features, functions, operations, blocks, or steps described or illustrated anywhere herein that a person having ordinary skill in the art would comprehend. Furthermore, reference in the appended claims to an apparatus or system or a component of an apparatus or system being adapted to, arranged to, capable of, configured to, enabled to, operable to, or operative to perform a particular function encompasses that apparatus, system, component, whether or not it or that particular function is activated, turned on, or unlocked, as long as that apparatus, system, or component is so adapted, arranged, capable, configured, enabled, operable, or operative. Additionally, although this disclosure describes or illustrates particular embodiments as providing particular advantages, particular embodiments may provide none, some, or all of these advantages.
All of the disclosed methods and procedures described in this disclosure can be implemented using one or more computer programs or components. These components may be provided as a series of computer instructions on any conventional computer readable medium or machine readable medium, including volatile and non-volatile memory, such as RAM, ROM, flash memory, magnetic or optical disks, optical memory, or other storage media. The instructions may be provided as software or firmware, and may be implemented in whole or in part in hardware components such as ASICs, FPGAs, DSPs, or any other similar devices. The instructions may be configured to be executed by one or more processors, which when executing the series of computer instructions, performs or facilitates the performance of all or part of the disclosed methods and procedures.
It should be understood that various changes and modifications to the examples described here will be apparent to those skilled in the art. Such changes and modifications can be made without departing from the spirit and scope of the present subject matter and without diminishing its intended advantages. It is therefore intended that such changes and modifications be covered by the appended claims.
The present application claims priority to and the benefit of U.S. Provisional Application 63/175,838, filed Apr. 16, 2021, and U.S. Provisional Application 63/184,891, filed May 6, 2021, the entireties of which are hereby incorporated herein by reference.
This invention was made with government support under Grant Number NS100549 awarded by the National Institutes of Health. The government has certain rights in the invention.
| Number | Date | Country | |
|---|---|---|---|
| 63175838 | Apr 2021 | US | |
| 63184891 | May 2021 | US |
