Systems and methods for measuring analgesia and response to harmful or pain stimuli

Abstract

The present disclosure provides systems and methods for providing a quantitative assessment of pain and analgesia. The systems and methods provide a direct, objective assessment of pain as well as a direct, objective assessment of pain alleviation by different drugs or treatments. Thus, such systems and methods provide feedback for adjusting a drug or treatment to a desired endpoint, and may be used in the discovery or development of novel drugs or treatments. The system and method include a novel state space evoked response potential (SS-ERP) model for extracting electroencephalogram (EEG) signals from electroencephalogram (EEG) data of a patient. The disclosed systems and methods enable the measurement of pain and pain treatment effects in an individualized manner.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application is based upon and claims priority to U.S. Provisional Application Serial No. 63 / 381,312, filed on October 28, 2022, which is incorporated herein by reference in its entirety for all purposes. Government Support Statement

[0002] This invention was made with government support under Grant R01DA056593 from the National Institutes of Health. The U.S. Government has certain rights in this invention. Background Art

[0003] Over 60 million surgeries are performed annually in the United States and over 250 million worldwide. A significant proportion of these patients may experience acute postoperative pain. Large-scale studies have shown that 26% of patients experience moderate to severe acute pain after surgery, and nearly 10% will go on to develop chronic pain. Postoperative pain can also lead to postoperative opioid misuse, contributing to the opioid crisis. Therefore, a key priority in perioperative pain management is to avoid or minimize the risk of both acute and chronic postoperative pain. Postoperative pain also negatively impacts hospitals and healthcare systems by increasing emergency department visits and readmissions, length of stay, and costs of care.

[0004] Currently, there are multiple strategies for providing analgesia during surgery, including opioid-based strategies during general anesthesia, multimodal anesthesia, and regional anesthesia. Unfortunately, a major potential limitation of these strategies is the lack of monitors that can directly and objectively measure pain and / or pain perception and assess whether adequate analgesia levels are maintained. Currently, anesthesiologists rely on changes in heart rate (HR) and blood pressure (BP) as indirect measures of pain perception and titrate opioids and other medications to minimize these changes. Unfortunately, HR and BP are affected by numerous intraoperative factors, such as blood loss, anesthetics, and antihypertensive medications, making them unreliable indicators of pain perception. On the other hand, a large proportion of surgeries rely on regional anesthetic techniques, such as epidurals or regional blocks, for pain perception and postoperative pain management. However, epidurals and blocks are not foolproof and can fail during surgery, resulting in potentially uncontrolled postoperative pain that is only apparent when the patient regains consciousness after surgery and is in severe pain. Currently, anesthesiologists cannot directly assess the continued efficacy of their regional blocks, resulting in unreliable pain perception and pain management. Techniques have been proposed to monitor pain sensation during surgery using either or both HR or BP measurements along with EEG. Unfortunately, the clinical evidence supporting the use of these devices is mixed. Use of these devices does not appear to reduce postoperative pain scores.

[0005] It is well known that painful stimuli generate event-related potentials (ERPs), which can be measured on the scalp by averaging the waveforms of repetitive stimuli generated by various methods, including electrical stimulation. These ERPs are challenging to use in clinical settings because their responses are very small, ranging from approximately 1 to 5 microvolts, compared to background EEG noise and oscillations, which can be an order of magnitude larger in amplitude. To overcome the small effect size (or low signal-to-noise ratio) of these ERPs, large numbers of stimuli, on the order of hundreds or more, must be delivered to average out the much larger background noise. Even then, the residual noise can be so large that obtaining reliable ERP responses in individual subjects or patients is difficult. Consequently, many studies employing ERPs present only aggregated results across groups of subjects and lack information on individual subjects. Distinguishing between different conditions, such as the presence or degree of pain treatment, is even more challenging. Even if all these difficulties can be overcome with sufficiently high stimulation numbers, practical clinical compliance challenges remain: Can patients tolerate repeated painful stimulation for tens of minutes? Will responses to such stimulation eventually saturate due to underlying processes of adaptation or tolerance? Can ERP stimulation delivered over such long timeframes track changes in pain perception or pain during surgery? Finally, these ERPs are even more challenging to use, if not impossible, during anesthesia or when taking analgesic drugs (such as opioids) because they are very small, on the order of 1 to 5 microvolts, and are masked by anesthesia-induced or analgesia-induced EEG oscillations that are 10 to 100 times larger in amplitude. Thus, although ERPs do have a relevant underlying physiological basis for objective pain and analgesia measurements, in their current form, the systems and methods available for calculating ERPs make them completely impractical in any practical setting, clinical or otherwise, especially where individualized responses are desired. Summary of the Invention

[0006] The present disclosure provides a system and method for overcoming the above-mentioned shortcomings by providing a quantitative assessment of pain and analgesia, and in the process of doing so provides significant advantages over products that use indirect measurements based on HR, BP and other autonomic signals. In addition, the system and method provide a direct, objective assessment of pain and pain relief for different drugs or treatments. Such systems and methods provide feedback for adjusting drugs or treatments to desired endpoints, which is useful in the process of discovering or developing novel drugs or treatments. The disclosed system and method can measure pain and pain treatment effects in an individualized manner. The disclosed system also makes it possible to assess pain in groups that are difficult or impossible to perform classical assessments (such as infants, non-verbal patients or intubated and sedated patients).

[0007] In one aspect of the present disclosure, a method for measuring intraoperative stimulation is described. The method includes receiving electroencephalogram (EEG) data from a patient using a patient monitor while applying one or more stimuli to the patient. The method further includes extracting one or more event-related potential (ERP) signals from the EEG data associated with one or more discrete events using a computer processor. The method further includes determining one or more stimuli from the extracted one or more ERP signals using the computer processor, and indicating, using a display, an adjustment of one or more intraoperative parameters associated with offsetting the one or more stimuli in the future one or more discrete events.

[0008] In another aspect of the present disclosure, a system for measuring intraoperative stimulation is described. The system includes one or more sensors and a processor, wherein the one or more sensors are configured to measure an electroencephalogram (EEG) signal of a patient. The processor is configured to receive an EEG signal from the patient when one or more stimuli are applied to the patient, extract one or more event-related potential (ERP) signals from the EEG data associated with one or more discrete events, determine one or more stimuli from the extracted one or more ERP signals, and determine an adjustment to one or more intraoperative parameters based on the one or more stimuli. The system further includes a display for communicating the adjustment to the one or more intraoperative parameters.

[0009] These aspects are non-limiting.Other aspects and features of the systems and methods described herein are provided below. BRIEF DESCRIPTION OF THE DRAWINGS

[0010] The foregoing features of the embodiments will be more readily understood from the following detailed description with reference to the accompanying drawings, in which:

[0011] Figure 1 is a schematic diagram of a system for measuring pain and analgesia in a patient under anesthesia according to the present invention.

[0012] Figure 2 is a flow chart setting forth the steps of a method according to the present invention for measuring pain and analgesia in a patient under anesthesia.

[0013] Figure 3 This is a graph of the P300 oddball response. In the oddball paradigm, the SS-ERP outperformed the average ERP in extracting the P300. These results are from only 52 trials (42 frequent trials + 10 oddball trials). The shaded area shows the 90% CI constructed from the posterior variance of the SS-ERP and the sample variance of the mean ERP.

[0014] Figure 4A This is the experimental design diagram of pain ERP under analgesia.

[0015] Figure 4B is a graph comparing electrode Cz ERPs evoked by four different levels of noxious stimulus intensity before remifentanil administration.

[0016] Figure 4C is a graph comparing electrode Cz ERPs evoked by four different levels of noxious stimulus intensity after remifentanil administration.

[0017] Figure 4D is a graph comparing the Fp2-Fp1 bipolar electrodes of the ERP evoked by four different levels of noxious stimulus intensity before remifentanil administration.

[0018] Figure 4E is a graph of the Fp2-Fp1 bipolar electrode comparison of ERPs evoked by four different levels of noxious stimulus intensity after remifentanil administration.

[0019] Figure 5 Receiver operating characteristic curve showing the use of SS-ERP P100 amplitude to discriminate between low and high electrical stimulations. AUC = Area Under the Curve.

[0020] Figure 6A is a graph of the group average results (n=10) of SS-ERP at different stimulation intensity levels at baseline.

[0021] Figure 6B is a graph showing the group mean results (n=10) of SS ERP under different stimulation intensity levels and mild analgesia regime.

[0022] Figure 6C is a graph showing the group average results (n=10) of SS-ERP under different stimulation intensity levels and moderate analgesia regime.

[0023] Figure 7A Figure 5. SS-ERP under general anesthesia. EEG time series were simulated using oscillator parameters estimated from EEG recordings of sevoflurane-induced anesthesia.

[0024] Figure 7B SS-ERP under general anesthesia is shown. SS-ERP makes it possible to extract ERPs that appear weak due to being covered by slow waves and alpha waves.

[0025] Figure 7CSS-ERP under general anesthesia is shown. SS-ERP (MSE: 0.658 au) follows the ground truth better than mean ERP (MSE: 11.51). The shaded area shows the 90% CI constructed from the posterior variance of SS-ERP and the sample variance of mean ERP.

[0026] Figure 8A A simulation study is shown, in which SS-ERP captured time-locked responses while extracting background oscillations, thereby eliminating the influence of background oscillations on the ERP waveform in a data-driven manner.

[0027] Figure 8B A simulation study is presented in which SS-ERP performance remains constant as the number of trials decreases, unlike the classic average ERP. Background oscillations produce spurious peaks in the average ERP over a smaller number of trials, whereas SS-ERP is relatively immune to such perturbations due to the explicit modeling of background oscillations.

[0028] Figure 8C A simulation study is shown, in which KE-ERP captured time-locked responses while extracting background oscillations, thereby eliminating the influence of background oscillations on the ERP waveform in a data-driven manner.

[0029] Figure 8D We present a simulation study in which KE-ERP performance remains constant as the number of trials decreases, unlike conventional ERP averaging using kernel expansion. Background oscillations can produce spurious peaks in ERP averaged over a smaller number of trials, whereas KE-ERP is relatively immune to such perturbations due to its explicit modeling of background oscillations.

[0030] Figure 9A Figure 2 is a plot of the oddball P300 responses from 42 frequent / 10 oddball trials in healthy individuals and amyloid-positive individuals. The traditional average ERP is contaminated by background oscillations.

[0031] Figure 9B Figure 2 shows the oddball P300 responses from 42 frequent / 10 oddball trials in healthy individuals and amyloid-positive individuals. The SS-ERP shows differences in the ERP, likely due to underlying pathology. The shaded area shows the 90% confidence interval constructed from the posterior variance of the SS-ERP and the sample variance of the ERP mean.

[0032] Figure 9C Figure 2 shows the oddball P300 responses from 42 frequent / 10 oddball trials in healthy individuals and amyloid-positive individuals. Conventional ERP averaging using kernel expansion is contaminated by background oscillations.

[0033] Figure 9D Figure 2 shows the oddball P300 responses from 42 frequent / 10 oddball trials in healthy individuals and amyloid-positive individuals. The KE-ERP shows differences in the ERP, likely due to underlying pathology. The shaded area shows the 90% confidence interval constructed from the posterior variance of the KE-ERP and the sample variance of the ERP mean. DETAILED DESCRIPTION

[0034] This paper discloses a novel technique for processing ERPs that improves their accuracy by a factor of approximately 150, even in the presence of orders of magnitude larger background oscillations. These methods and systems make it possible, for the first time, to use ERPs to actively assess pain and analgesia during general and regional anesthesia.

[0035] ERPs provide a non-invasive method for studying the psychophysiological correlates of sensory and cognitive processes, with their components providing information about sensory (“exogenous”) and cognitive (“endogenous”) processes with millisecond temporal resolution. ERPs are tiny, approximately 1 μV signals embedded in background spontaneous oscillations, which can be 10 to 100 times larger than ERPs. Traditionally, combining a large number of individual trial waveforms into an averaged ERP waveform has been considered the “gold standard” for averaging event-independent activity. However, it is often difficult to obtain enough trials to extract these time-locked evoked responses, as a large number of trials is required to effectively average out or cancel out the ongoing oscillatory background. Furthermore, over-averaging can easily obscure the fine structure of ERPs, as intrinsic factors such as mental workload, fatigue, and attention are known to modulate time-locked responses. To reduce this reliance on a large number of trials, the methods and systems described herein use a novel linear state-space model to model background oscillations, effectively removing the ongoing spontaneous activity.

[0036] Here, the human brain is modeled as a linear time-invariant system, consistent with the classic ERP literature, so the response evoked by stimulus presentation shown in the EEG is written as an ERP waveform and pulse trains associated with discrete events of stimulus presentation t Convolution between:

[0037] However, unlike conventional analyses, spontaneous stimulus-independent activity in the EEG is modeled as the superposition of multiple oscillators, modeled as state-space oscillators: Among them, assuming is independent, |f (m) |<f s , f s is the sampling frequency, and a (m)<1. In other words, the oscillatory temporal dynamics are explicitly imposed on the background activity, rather than treating it as white noise or other unstructured noise. In summary, oscillator models are used in conjunction with ERP models to interpret EEG recordings, effectively removing contamination from strong oscillations:

[0038] We further include priors on the ERP waveforms to ensure robust recovery and improve the variability of the estimates. Furthermore, we impose a temporal continuity constraint in the form of the following state-space model: This results in a multi-level state-space framework called state-space ERP (SS-ERP).

[0039] Alternatively, assuming that the ERP, h allows different priors for the following expansion on the dictionary consisting of the given basis functions, in terms of the expansion coefficients has a zero-mean Gaussian prior on , with a diagonal ∧ covariance matrix: This results in a state-space framework called kernel-extended ERP (KE-ERP).

[0040] Reference Figure 1 , a system 100 configured for use in accordance with the present invention includes a patient monitoring device 102, such as a physiological monitoring device, Figure 1 , which is an electroencephalogram (EEG) electrode array. However, it is contemplated that the patient monitoring device may also include a mechanism for monitoring galvanic skin response (GSR), for example, to measure arousal to external stimuli. One specific implementation of this design utilizes a frontal Laplace EEG electrode layout with additional electrodes to measure the GSR. Another implementation of this design combines a frontal array of electrodes that can be combined in post-processing to obtain any combination of electrodes found to optimize detection of the previously described EEG features, also with separate GSR electrodes. Another implementation of this design utilizes a high-density layout, using between 64 and 256 sensors to sample the entire scalp surface for the purpose of source localization, also with separate GSR electrodes.

[0041] The patient monitoring device 102 is connected via a cable 104 to communicate with a monitoring system 106. Furthermore, the cable 104 and similar connections may be replaced by wireless connections between the components. As shown, the monitoring system 106 may be further connected to a dedicated analysis system 108. Furthermore, the monitoring system 106 and the analysis system 108 may be integrated. The patient monitoring device 102 and the monitoring system 16 may each further include a user interface 110. The user interface 110 may be, but is not limited to, a keyboard, mouse, touch screen, or buttons or switches integrated into the monitoring system 106.

[0042] For example, as described above, it is contemplated that the patient monitoring device 102 may be an EEG electrode array, e.g., a 64-lead EEG electrode array. However, as will be apparent below, increasing the number of electrodes from 64 to 128, 256, or even more can achieve greater spatial accuracy. Similarly, the present invention can be implemented with substantially fewer electrodes. In any case, the monitoring system 106 can be configured to receive the raw signals collected by the EEG electrode array and assemble and even display the raw signals as EEG waveforms. Thus, the analysis system 108 can receive the EEG waveforms from the monitoring system 106 and, as will be described, analyze the EEG waveforms and features therein based on the selected anesthetic compound, determine the patient's status based on the analyzed EEG waveforms and features, and generate a report, e.g., as a printed report or, preferably, as a real-time display of the feature information and the determined status. However, it is also contemplated that the functionality of the monitoring system 106 and the analysis system 108 can be combined into a common system.

[0043] Reference Figure 2 , a method 200 for measuring pain and analgesia in a patient under anesthesia is shown. The method may be performed by Figure 1 The method is executed by a computer processor associated with the monitoring system 106. At step 202, EEG data is received from a patient using a patient monitor while one or more stimuli are applied to the patient.

[0044] In non-limiting examples, the one or more stimuli may include an electric shock, or an auditory or visual stimulus. The electric shock may be provided by one or more electrodes attached to the patient. The electric shock may be variable to adjust the voltage of the shock presented to the patient.

[0045] In a non-limiting example, the EEG data includes a signature of one or more analgesics.Additionally, the EEG signal is collected by one or more sensors (such as the patient monitoring device 102).

[0046] Then, at step 204, one or more ERP signals are extracted from the EEG data associated with the one or more discrete events. Furthermore, at step 206, one or more stimuli are determined from the one or more ERP signals. In a non-limiting example, the one or more stimuli include noxious or painful stimuli. In a further non-limiting example, the ERP signals are associated with the one or more stimuli.

[0047] In a non-limiting example, a processor of a monitoring system decomposes EEG data using a plurality of state-space oscillators derived from a background signal of the EEG data, one or more ERP signals, and one or more discrete events. A state-space ERP (SS-ERP) model is included to extract the one or more ERP signals from the EEG data. Furthermore, the SS-ERP model includes a plurality of state-space oscillators derived from the EEG signal, the one or more ERP signals, and the background signal of the one or more discrete events.

[0048] In a non-limiting example, one or more ERP signals are separated from an oscillatory background using the following composite model: Among them, assuming is independent, |f (m) |<f s , f s is the sampling frequency, a (m) <1, and is a state-space oscillator, and where x t,1 and x t,2 is the generalized phasor representation of the component oscillations x t The real and imaginary components, f is the oscillation frequency, f s is the sampling frequency, α is the damping factor, ∈ is the oscillation state noise, σ 2 is the state noise variance, n t is the measurement noise, and τ 2 is the noise variance. The evoked response is represented as an ERP waveform Presented with one or more stimuli t Convolution between: The following examples describe the ERP model in more detail.

[0049] Alternatively, the induced component can be expressed as an extension of a dictionary consisting of given basis functions. As a representative of such basis functions, consider a Gabor dictionary G, whose atoms are given by Gaussian kernels of fixed variance, normalized to have a maximum value of 1. In the expansion coefficient with diagonal ∧ covariance matrix A zero-mean Gaussian prior is considered above. This ERP is called kernel-extended ERP (KE-ERP).

[0050] At step 208, adjusting one or more intraoperative parameters associated with counteracting one or more stimuli in the future one or more discrete events is indicated using a display of the patient monitor. In a non-limiting example, adjusting the one or more intraoperative parameters includes adjusting administration of one or more analgesics.

[0051] As used in this application and the claims, the singular forms "a," "an," and "the" include plural referents unless the context clearly dictates otherwise.

[0052] As used herein, "approximately," "approximately," "substantially," and "significantly" will be understood by those of ordinary skill in the art and will vary to some extent depending on the context in which they are used. If the usage of a term is unclear to one of ordinary skill in the art given the context in which the term is used, "approximately" and "approximately" will mean up to plus or minus 10% of the particular term, and "substantially" and "significantly" will mean more than plus or minus 10% of the particular term.

[0053] As used herein, the terms "include" and "including" have the same meanings as the terms "comprise" and "comprising." The terms "comprising" and "comprising" should be interpreted as "open" transition terms that permit the inclusion of further additional components beyond the components recited in the claims. The terms "consisting of" and "consisting of" should be interpreted as "closed" transition terms that do not permit the inclusion of additional components beyond those recited in the claims. The term "consisting essentially of" should be interpreted as partially closed and only permits the inclusion of additional components that do not fundamentally change the nature of the claimed subject matter.

[0054] The phrase "such as" should be interpreted as "for example, including." Moreover, the use of any and all exemplary language, including but not limited to "such as," is intended merely to better illuminate the invention and does not pose a limitation on the scope of the invention unless otherwise claimed.

[0055] Furthermore, in those instances where a convention similar to "at least one of A, B, and C, etc." is used, generally speaking, such construction is intended to be in the sense that one of ordinary skill in the art would understand the convention (e.g., "a system having at least one of A, B, and C" would include, but is not limited to, systems having A alone, B alone, C alone, A and B together, A and C together, B and C together, and / or A, B, C together). Those skilled in the art will further understand that virtually any disjunctive word and / or phrase presenting two or more alternative terms, whether in the specification or the drawings, should be understood to contemplate the possibility of including either, either, or both of the terms. For example, the phrase "A or B" would be understood to include the possibility of "A" or "B" or "A and B."

[0056] All language such as "up to," "at least," "greater than," "less than," etc., includes the recited numbers and refers to ranges that can be subsequently subdivided into ranges and subranges. Ranges include each individual member. Thus, for example, a group having 1-3 members refers to groups having 1, 2, or 3 members. Similarly, a group having 6 members refers to groups having 1, 2, 3, 4, or 6 members, etc.

[0057] The modal verb "may" refers to a preferred use or selection of one or more options or choices among several described embodiments or features included in such embodiments. If there is no disclosure regarding the options or choices included in a particular embodiment or feature, the modal verb "may" refers to an affirmative action regarding how to make or use an aspect of the described embodiment or feature, or a definite decision regarding the use of specific techniques associated with the described embodiment or feature. In this latter context, the modal verb "may" has the same meaning and connotation as the auxiliary verb "may." Example

[0058] State-space event-related potentials (SS-ERPs) and kernel-extended event-related potentials (KE-ERPs) can separate ERPs from background oscillations with approximately 150 times greater accuracy than conventional averaging methods. Event-related potentials (ERPs) are time-locked EEG responses to sensory, cognitive, or motor events, characterized by prominent peaks and troughs in a waveform ranging from approximately 1 microvolt to approximately 10 microvolts in size, occurring within 10 to 100 milliseconds. A major drawback of current methods is that they rely on averaging over hundreds of trials to extract these time-locked evoked responses from a background of ongoing slow and / or alpha wave oscillations. Consequently, the structure of these conventionally averaged ERPs can be easily obscured by these background oscillations, particularly if these background oscillations are 10 or 100 times larger than the conventionally averaged ERPs during general anesthesia and sedation. However, as described here, the "SS-ERP" method, which uses a multilevel state-space model, and the "KE-ERP" method, which employs a sparse representation on a predefined basis, or kernel, dictionary, improve ERPs by explicitly modeling and effectively removing ongoing spontaneous background activity. This method further improves the accuracy of ERP by imposing a temporal continuity prior on the ERP waveform to recover the denoised ERP.

[0059] The performance of the SS-ERP is benchmarked by analyzing the P300 auditory "oddball" response, which is generated by comparing the response to a frequently occurring auditory tone (e.g., 2000 Hz) with the response to a less frequent "oddball" tone (1000 Hz). Typically, several hundred trials are required to resolve the ERP waveform. Figure 3 Results from the SS-ERP algorithm are presented for a single subject using only 52 trials (42 frequent trials, 10 oddball trials). Although the differences between the frequent and oddball waveforms are difficult to interpret using traditional averaging (Cohen's d: 3.57, precision (inverse of the estimated variance): 0.56, minimum detectable event (MDE, 2.485*std): 3.75 μV), the waveforms derived from the SS-ERP are easily interpreted, as there are clear differences between the frequent and oddball conditions, with narrow confidence intervals (Cohen's d: 21.27, precision: 85.33, MDE: 0.303 μV). This analysis revealed a significant 5.96-fold improvement in effect size (Cohen's d), a 12.37-fold reduction in MDE, and a substantial 152.4-fold improvement in precision, all within a single subject's 52 trials using the SS-ERP.

[0060] We performed the same benchmark on KE-ERP, but with a different number of trials in a single subject (118 frequent and 30 oddball trials). Like SS-ERP, KE-ERP (Cohen's d: 167.68, precision: 178.98, MDE: 0.2 μV) showed significant improvement over the classic average ERP (Cohen's d: 5.10, precision: 0.112, MDE: 7 μV).

[0061] Opioids reduce ERP features associated with the intensity of painful stimuli to levels comparable to those of non-painful stimuli; these features can be estimated using SS-ERP methods from as few as 30 trials at parietal or frontal electrodes. ERPs evoked by painful somatosensory stimulation can be clearly seen in EEG channels placed directly over the somatosensory cortex (such as Cz), but the same problem applies to the P300 oddball response: the pain ERP is orders of magnitude smaller than the overlying background EEG signal.

[0062] In addition, EEG recordings were analyzed from healthy individual volunteers who received electrical stimulation at four different intensities calibrated to the subjects' pain perception so that only the highest stimulation intensity, "level 4," was perceived as painful. The stimulation consisted of 10 trains of 3 electrical pulses at each intensity level, randomly shuffled ( Figure 4A ). Remifentanil (3 ng / mL) was then administered, followed by 5 trains of 3 electrical pulses at each pain intensity level (randomly shuffled). Similar to P300, the SS-ERP method significantly improved the accuracy of pain ERP estimation compared to traditional averaging ( Figures 4B-4E ). The waveform from the Cz electrode ( Figures 4B-4C ) shows visible differences between stimulus intensities, but with much higher precision when the SS-ERP method is used. Using the SS-ERP method, it is easy to see that the ERP waveform is reduced after remifentanil administration ( Figure 4C ). In contrast, conventional averaging revealed considerable residual structure in the ERPs, with large overlapping confidence bounds ( Figure 4B ).

[0063] In this individual subject, considering the detection of pain ERP by the prefrontal electrodes Fp1 and Fp2, the prefrontal electrodes Fp1 and Fp2 can be more easily used for clinical applications when they are located on the forehead, and due to their proximity to the prefrontal cortical structures that mediate this perception, the prefrontal electrodes Fp1 and Fp2 can provide more information about the emotional perception of pain. Figure 4D ), cannot effectively estimate ERP from Fp1-Fp2, and cannot discern differences between stimulus intensity levels or when opioids are administered. In contrast, the SS-ERP method can detect clear differences between stimulus intensity levels at baseline ( Figure 4E Furthermore, using SS-ERP, responses to painful stimuli of “level 4” became indistinguishable from stimuli of lower intensities when remifentanil was administered ( Figure 4E ).

[0064] Figure 5 The receiver operating characteristic (ROC) curve for distinguishing low from high electrical stimulation using the P100 ERP peak during the baseline period before analgesia is shown. The AUC for distinguishing low from high electrical stimulation was 0.86, indicating that individual subjects had highly accurate discrimination. This analysis shows that the method described herein can be used for highly accurate diagnostic inference in individual patients with only a small number of stimulation events.

[0065] Next, we verified whether these same results hold true, on average, for a larger group of subjects. Figures 6A-6C Shown are group average results for n=10 subjects who underwent the same protocol presented previously, where four different electrical stimulation intensities were presented, tuned so that intensity level 4 was subjectively rated as painful by the subjects. Figure 6A ) and mild analgesia provided by intravenous infusion of remifentanil ( Figure 6B ) and moderate analgesia ( Figure 6C ) conditions. Group average results were consistent with those observed in individual subjects: the amplitude of the SS-ERP response increased with increasing stimulus intensity and gradually decreased with the administration of mild and moderate analgesia, respectively.

[0066] We further considered the diagnostic value of the SS-ERP signal for predicting pain levels. To summarize the SS-ERP response, in one example, we analyzed the amplitude between the initial trough and the subsequent peak of the response at the central EEG channel Cz, which was represented in the statistical model as the variable "SS-ERP." A mixed-effects model was constructed to assess the association between the SS-ERP amplitude and each subject's reported pain: pain ~ SS-ERP + (SS-ERP | subject).

[0067] The SS-ERP amplitude was found to be positively correlated with the pain reported by each subject on each trial (estimate 0.18, SE 0.02, t-value 6.88, p-value 6.2e -12 ).

[0068] Given that other variables inherent to the experimental protocol might influence reported pain, a second mixed-effects model was constructed, which now included two covariates (the level of analgesia administered and the intensity of stimulation): pain ~ SS-ERP + analgesia + stimulation + (SS-ERP | subject).

[0069] The fixed variables were the amplitude of the ERP obtained using this method (SS-ERP), the amount of opioid administered (analgesia), and the intensity of the stimulation (stimulation). The random variable was the subject. <![CDATA[Chi-square (Chi 2 )]]> p-value SS-ERP amplitude 6.9 0.008 analgesia 62.1 <![CDATA[2.12e -13 ]]> Stimulus intensity 73.7 <![CDATA[7.0e -16 ]]>

[0070] Notably, after adjusting for stimulus intensity and the amount of analgesia received, the magnitude of the SS-ERP remained correlated with subject-reported pain, suggesting that the method is able to capture individual differences in responses to stimulus intensity and analgesia between subjects, i.e., the SS-ERP measures pain information that is personalized to each specific subject. When the analysis was repeated using differences between frontal EEG electrodes (Fp2-Fp1), the statistical significance of the association between the SS-ERP and pain was even higher (p-value = 8.9e -5 ), which is consistent with the notion that frontal electrodes may provide more information about the affective component of pain mediated by prefrontal cortical structures.

[0071] These data demonstrate how the SS-ERP method can detect differences in painful stimulus intensity and a reduction in the ERP waveform in response to pain treatment (opioid analgesia). Furthermore, these data demonstrate how these differences in pain intensity and analgesia can be detected using only 15 to 30 stimulus events, suggesting that this method could be practically applied in clinical settings to objectively measure pain over a short period of time.

[0072] The SS-ERP method can accurately recover the ERP in the presence of anesthesia-induced oscillations that are approximately 10 times larger than the ERP. During general anesthesia and sedation, the EEG exhibits large, patterned oscillations in specific frequency bands. Propofol and sevoflurane, two of the drugs most commonly used to maintain unconsciousness during general anesthesia, induce slow and delta oscillations (approximately 1-4 Hz) as well as alpha oscillations (approximately 9-12 Hz). These slow / delta and alpha oscillations are much larger than the expected ERP variations of approximately 1 to 10 microvolts during somatosensory stimulation. In this case, extracting the ERP response using traditional averaging methods would require a large number of stimulation events to overcome the substantially greater noise. This would make ERPs for pain assessment impractical in clinical settings due to unacceptably long latency times for evaluating ERP measurements. Furthermore, excessively long stimulation event trains may cause discomfort to some patients. To verify that SS-ERP can effectively extract the ERP in the presence of large anesthesia-induced background oscillations, a realistic simulation study was conducted, characterized by the masking of the simulated ERP waveform by the much larger simulated anesthesia-induced oscillations. Background oscillations were generated from a state-space oscillation model whose parameters were estimated from patients receiving general anesthesia with sevoflurane. Figure 7AAs shown in , the simulated data are similar in form to the real data in both the time and frequency domains. An event-related signal generated by replicating a sequence of 30 real ERP waveforms was added to the simulated oscillations, which had an SNR of approximately 1 / 10. Figure 7B The decomposition of the signal into oscillatory and event-related components is shown. The SS-ERP method is clearly able to extract the small ERP signal from the background of strong oscillations. Compared with the average ERP (mean square error, MSE: 11.51 au), the SS-ERP (MSE: 6.658 au) follows the ground truth more faithfully (MSE is reduced by 17.5 times) and follows the ground truth with tighter confidence intervals ( Figure 7C These results demonstrate that the SS-ERP method can estimate ERP waveforms even in the presence of much larger anesthesia-induced EEG oscillations.

[0073] Method Overview

[0074] Novel ERP analysis method.

[0075] State-space ERP. The SS-ERP algorithm employs a multilevel model that represents the ongoing oscillatory background signal distinct from the evoked response, making it possible to separate these two components and significantly reduce the uncertainty of the ERP estimate. Background oscillations The event-related time series is modeled as the convolution between the ERP waveform h and the discrete events e (modeled as impulses), and then the observed EEG signal is represented as the sum of the background and ERP components with an observation noise term: Among them, assuming is independent, |f (m) | <f s , f s is the sampling frequency, and a (m) < 1. This multi-level framework also allows us to include priors on the ERP waveforms to ensure robust recovery and reduce the variability of the estimates, implemented here in the form of a random walk prior: h k =h k-1 +v k , v k ~N(0,ξ 2 ) (C)

[0076] Study ERP. Oscillator parameters Noise variance τ 2 and smoothness parameter ξ 2Learning is performed using an instance of the expectation-maximization (EM) algorithm. EM alternates between optimizing the distribution of hidden oscillator states and ERP waveforms given the current parameters (E-step) and updating the parameters given the hidden state distribution (M-step). The E-step is further simplified by constraining the posterior distribution of the ERP waveform h and the oscillatory state to be statistically independent: q(h,{x t}|θ,τ 2 ,ξ 2 )=q(h|τ 2 ,ξ 2 )q({x t}|θ,τ 2 ) This allows for periodic updates of the oscillator states and ERP distributions. The various model parameters can then be re-estimated in closed form in an M-step. This approach is known as variational Bayesian approximation, and the modified algorithm is known as generalized EM. Specifically, the ERP update step resembles ridge regression, removing the estimated oscillations from the EEG signal. Notably, this framework provides a posterior distribution of the ERP, compared to point estimates in classical averaging techniques. The update stops when the observed likelihood stabilizes.

[0077] Kernel-extended ERP. This model is applied to Equation (A), assuming that the ERP allows the following extension on a dictionary formed by a given basis function. As a representative of such basis functions, consider a Gabor dictionary G, whose atoms are given by Gaussian kernels of fixed variance, normalized to have a maximum value of 1. We again consider the expansion coefficients with a diagonal ∧ covariance matrix A zero-mean Gaussian prior on .

[0078] Learning ERP. Similar to SS-ERP, the oscillator parameters Noise variance τ 2 The prior variance of the expansion coefficient Λ is learned using an instance of the expectation maximization (EM) algorithm. EM alternates between optimizing the distribution of hidden oscillator states and ERP waveforms given the current parameters (E step) and updating the parameters given the hidden state distribution (M step). ～ The E-step is further simplified by constraining the posterior distribution of the oscillatory state to be statistically independent: This allows for periodic updates of the oscillator state and the ERP distribution. The different model parameters can then be re-estimated in closed form in the M-step. In fact, this learning of Χ encourages the coefficients of irrelevant lags to approach zero, thus performing automatic correlation detection (ARD) in a fully data-driven manner.

[0079] In summary, SS-ERP and KE-ERP easily capture time-locked responses while simultaneously extracting background oscillations, thereby eliminating the influence of background oscillations on the waveform in a data-driven manner. Furthermore, compared to point estimates in traditional averaging techniques, both SS-ERP and KE-ERP frameworks provide posterior distributions of ERPs.

[0080] Simulation results

[0081] Figures 8A-8B Shown is the extraction of SS-ERPs from simulated auditory tone responses contaminated by strong slow and alpha oscillations. Figure 8A We depict how oscillatory components can be effectively removed to retrieve auditory evoked responses through explicit modeling of oscillatory activity. Figure 8B The SS-ERP was compared with the traditional averaged ERP and the ground truth for increasing numbers of trials (n = 25, 50, 100). Clearly, the SS-ERP is very close to the ground truth in all three cases, with narrow confidence intervals slightly deteriorating as n decreases. In contrast, the averaged ERP is dominated by large oscillatory components present in the background: meaningful results were only obtained after averaging over 100 trials. This suggests that the SS-ERP has an advantage in capturing the true latency and amplitude of the ERP peak from a small number of trials.

[0082] Figures 8C-8D KE-ERP extractions from similar examples are shown. They correspond to Figure 8A - B. Obviously, the observations made for SS-ERP also apply here.

[0083] Application of real EEG data

[0084] SS-ERP performance was also benchmarked when detecting the P300 using EEG data recorded from a single subject (single-channel, parietal lobe): the subject was asked (with their eyes closed) to respond to a target tone (2000 Hz) that was intermixed with a regular tone (1000 Hz) 20% of the time. Only the first 52 trials out of a total of 200 trials were included in the analysis to demonstrate the improved ERP extraction using the SS-ERP approach: while 52 trials barely allowed the traditional averaging method to identify the P300 response (effect size: 3.57, Cohen's d, precision: 0.56), SS-ERP captured the response with high confidence by removing slow waves and sustained alpha from the background (effect size: 21.27 Cohen's d, precision: 85.33). This shows a significant improvement in effect size, with a 5.96-fold increase in effect size and a 12.37-fold increase in precision using only 52 trials. Furthermore, SS-ERP clearly showed a reduction in P300 responses in amyloid-positive individuals compared to healthy individuals ( Figure 9B ), while strong background oscillatory activity interferes with such effects in traditional averaged ERPs ( Figure 9A ).

[0085] We tested the KE-ERP on the same benchmark, but with a different number of trials in a single subject (118 frequentist trials, 30 oddball trials). Like the SS-ERP, the KE-ERP (Cohen's d: 167.68, precision: 178.98, MDE: 0.2 μV) showed a significant improvement over the classical average ERP (Cohen's d: 5.10, precision: 0.112, MDE: 7 μV). The KE-ERP again showed a reduction in the P300 response in amyloid-positive individuals compared to healthy individuals ( Figure 9D ), while strong background oscillatory activity interferes with such effects in traditional averaged ERPs ( Figure 9C ).

[0086] Thus, while significantly reducing the number of experimental requirements, the proposed method allows tracking short-term changes in ERPs due to various intrinsic and extrinsic causes, which has important implications for neuroscience research and clinical applications.

[0087] Learning parameters for a state-space oscillator model

[0088] Given univariate time series data y from neural signal recordings, a key goal is to determine which neural oscillations, if any, are present in the data y. A related question is to characterize the properties or parameters of these neural oscillations. The iterative oscillator method (iOsc) addresses these issues.

[0089] iOsc is a greedy search algorithm that attempts to represent y by adding one neural oscillatory component at a time until a pre-specified number of stops is reached. Model selection is then performed based on some measure of how well different numbers of oscillations represent the original data. The output of the iOsc algorithm is the set of selected neural oscillations, including the number of oscillations and all parameters of the oscillations, which are considered to be the most representative oscillator components of y.

[0090] y is represented in iOsc via time-domain modeling using a parametric oscillator state-space model of the following class:

[0091] This model class was chosen as a building block for single oscillations because it provides an interpretable and effective model to represent random oscillatory processes with random frequency modulations that are consistently observed in human EEG recordings, unlike pure sine waves encountered in physical science problems such as radar. When multiple neural oscillations are present simultaneously, multiple independent oscillators are easily cascaded in block diagonal matrix form in Eq. (1) and yield an equally interpretable model by simple summation in Eq. (2) to produce y. An additional benefit of time-domain modeling is that the log-likelihood provides a natural metric to assess how well different models with different numbers of oscillators can represent y, allowing a straightforward model selection scheme.

[0092] At each iteration, iOsc performs time-domain modeling of y using one or more oscillator components in the form of equation (1). The unknown parameters of each oscillator are {a, f, σ 2}, and the single observation noise variance τ 2 . If the values of these parameters are known, state estimation (also known as exact inference or belief propagation) can be performed efficiently using Kalman filtering and fixed-interval smoothing to estimate the hidden variable x (which represents the time series of the underlying oscillatory component). Conversely, given the value of x, the parameter values in this linear Gaussian setting can be updated using closed-form equations. Alternating between these two steps is an instance of the expectation-maximization (EM) algorithm with the oscillator x and the data y. The detailed equations for exact inference and the parameters are outlined below.

[0093] While each iteration of iOsc simply performs EM learning, and the computations are known and efficiently implemented, a key challenge is how to initialize EM in iOsc. In this linear Gaussian case, the complete data log-likelihood landscape is concave, with a single global maximum; however, there are a large number of unknown parameters in y and only a finite number of time points. In the case of poor initialization, EM may require infinite iterations or data length to converge to an appropriate value, resulting in inferior iOsc results. Therefore, a key step in successfully identifying potential oscillator components using iOsc is to choose an informed starting point throughout the greedy search process.

[0094] Oscillator initialization

[0095] iOsc's algorithm is designed to add one oscillator at a time; therefore, a scheme for initializing model parameters is constructed for a single oscillator. Consider data y observed at a sampling rate Fs with a pre-specified number of oscillation stops, n, implying a maximum of n oscillations in the recording. For example, in human EEG, which typically includes 2-4 distinct oscillations, n = 7 is a reasonable but non-limiting choice. The following steps initialize the oscillator.

[0096] Step 1: Fit an autoregressive (AR) model

[0097] First, an AR model of order 2n-1 is fitted to y. The order 2n-1 is chosen because the complex roots of the AR model appear in pairs in frequency. For example, when n = 7, the AR (13) model is used to identify the oscillation pole at DC and six frequencies above 0 Hz. The AR fit is performed, for example, using the Yule-Worker method or the Burg algorithm, and the fitted AR coefficients are expressed as a = {a1, a2, ... a 2n-1}:

[0098] Step 2: Select the largest pole / zero

[0099] The complex roots of the fitted AR model can be identified by solving the roots of a polynomial with coefficients {1, -a}. The theoretical PSD spectrum S of the AR process can also be obtained. To initialize the oscillator parameters {a, f}, first find the root corresponding to the highest power. Then, by replacing the pole frequency f p Theoretical PSD S weighted pole radius r at p To estimate the oscillator strength at the extremes: where p comes from the set of all poles / zeros with unique oscillation frequencies using complex roots. Once the pole / zero with maximum power is identified, the oscillator parameters {a,f} can be initialized as: When zero is chosen as the largest root, a frequency of 0 Hz is not ideal because these oscillator parameters are initialized for EM learning. Kalman filtering starting at this value prevents the oscillator frequency from being updated to anything different than 0 Hz. Therefore, during initialization, a lower limit is placed on f at 0.1 Hz, although subsequent EM learning can update it closer to 0 Hz.

[0100] Step 3: Initialize the oscillator state noise covariance

[0101] After fitting the AR model, there is no obvious way to initialize the state noise variance parameter σ of a single oscillator 2 This is because the estimated white noise process variance s 2 Simultaneously excite all poles and zeros in the AR process, rather than specifically one or more poles at a frequency f where one might wish to place an oscillator. This technical challenge can be overcome by combining the AR process and the noise variance s 2 It is solved by converting to different eigenmodes and reading the noise variance only from the poles located at the selected frequencies.

[0102] The key idea here is to use the spectral decomposition of the fitted AR process to derive σ 2 The AR(2n-1) model is written in the following state-space form: The following shorthand notations are defined: as well as G=[1 0 … 0]. Therefore, the AR(2n-1) equation is rewritten in the following matrix form: x t =Fx t-1 +w t (7). y t =Gx t (8). Consider the eigendecomposition of the transfer matrix F: F=VΛV -1 (9). Equation (7), the top row can be obtained by multiplying both sides of the equation by the orthogonal matrix V -1 To rearrange: V -1 x t=ΛV -1 x t-1 +V -1 w t (10)

[0103] The diagonal matrix ∧ contains complex eigenvalues, which are actually the same poles and zeros of the characteristic function with AR coefficients {a}. This equivalence exists because F is the adjoint matrix of the polynomial equation with coefficients {1, -a}. In fact, modern computing software solves polynomials via eigendecomposition of the adjoint matrix. Therefore, not only do the complex eigenvalues in ∧ match the complex roots in step 2, but they are also in exactly the same order.

[0104] Note that ∧ can be viewed as a transfer matrix for a set of parametric oscillators. Each complex number can be written in complex exponential form that corresponds to a rotation matrix of the frequency f times the damping parameter a, where the damping parameter a is the magnitude of the complex eigenvalue (or equivalently, the radius of the complex root). In other words, except that the small differences in the state noise process are not diagonal, equation (10) also describes the state of the oscillator after rotating V. -1 A set of block diagonal cascaded oscillators in the space of , just like in equation (1). The only mathematical difference is that the noise processes driving the real and imaginary components of each oscillator are now correlated and also correlated between different oscillators. Although this interpretation of oscillators is not further exploited for the eigenmodes, it provides the possibility of initializing the oscillator parameters σ using the variance of the noise processes driving these eigenmodes. 2 motivation.

[0105] The observation equation in (8) is also modified to be in the same rotation space: t =GVV -1 x t , with a new observation vector G'=GV, which is just the first row of V, or V for short 1,· : y t =V1,V -1 x t (11).

[0106] Since V 1,· The observation data y is formed by taking a linear combination of the state values from different oscillators (poles and zeros) in the rotation space t , all entries of V1,· are masked to 0 except for the poles selected by equation (4):

[0107] Observe that the only noise process in this rotation space is V -1 w t . Multiply by the mask observation vector The noise variance is scaled back to the scale of the observation data y, and only the contribution of a single extreme point is kept. Therefore, by applying the mask observation vector to the rotated noise covariance matrix Q' = V -1 Q(V -1 ) * , an estimate of the pole-specific noise variance can be derived. Since Q is nonzero only in the first diagonal entries, the rotation noise covariance matrix is rewritten as: in is the first column of the inverse of V, M * is the conjugate transpose of M.

[0108] Recall that complex roots (or equivalently, complex eigenvalues) behave as pairs in frequency, except at DC. This means that both poles in the AR process contribute to the spectral decomposition at any frequency greater than 0 Hz. This can be accounted for by summing the variances between poles with the same rotation frequency, which ignores cross-correlations between poles. Finally, we obtain the oscillator state noise variance parameter σ 2 The expected estimated value of :

[0109] In summary, the key insights of using this eigendecomposition in iOsc are: 1) relating the complex eigenvalues to the poles and zeros of the AR process; 2) identifying the complex eigenvalues under spectral decomposition as the rotated oscillator; and 3) applying a masked observation vector in the rotated space to combine the variances of the associated poles and zeros to initialize the oscillator parameters σ 2 .

[0110] Step 4: Repeat for Additional Oscillators

[0111] When y is the observed data, steps 1-3 provide the first oscillator with the parameters {a, f, σ 2} is an informed starting point. In human EEG recordings, the first oscillators added are most often slow oscillators (<1 Hz), because the amplitudes of such slow oscillations tend to be large, and therefore the power is also large. Recall that the search algorithm of iOsc will now try to add additional oscillators on top of the first oscillator learned by EM. The results of EM learning can be used for the existing oscillators, but the parameters {a, f, σ 2 To achieve this, the AR fitting process is reapplied to the new information, i.e., the one-step forecast error (OSPE):

[0112] is a reasonable time series to repeat steps 1-3, because the current model only includes the existing oscillators in the transfer matrix F, so it only predicts the oscillatory activity characterized by the existing oscillators, not the additional oscillations present in the observed data y that have yet to be characterized. Therefore, the next oscillator to be added by iOsc can be fitted to this OSPE by choosing The strongest poles / zeros of the AR model are identified.

[0113] In contrast, if we use the filter-based estimate y t|t =Gx t|t The residual error is used to find the additional oscillator, that is, These additional oscillations will not be visible because the Kalman filter adjusts the state estimate based on the most recent observations, absorbing any model misspecification errors into the state estimate. This is also true for smoothers such as fixed-interval smoothers. Therefore, using OSPE as described above to discover additional oscillators is an essential feature of the disclosed method.

[0114] The algorithm design provides an efficient implementation because After replacing y, the same routine in steps 1-3 can generally be repeated identically. However, there are two additional technical details to note. First, when there are only a small number of oscillators in the model, strong oscillations in the low-frequency range may not be fully explained because the oscillators learned by EM have tilted frequencies and / or damping parameters to accommodate the spectral content in the higher frequencies. Therefore, by fitting the residuals The next strongest pole identified by the AR process may still be at a low frequency. This is problematic because in the next iteration, there are two oscillators competing with each other in this frequency range. To circumvent this problem, when initializing new oscillators, a minimum frequency resolution (typically 1 Hz) is imposed, meaning that poles and zeros are only considered if they are more than 1 Hz away from existing oscillators. This is a moderate constraint applied only during initialization, and if closer oscillations do exist, EM learning can and will update the oscillator frequencies to be less than 1 Hz apart. This frequency resolution is only used to prevent erroneous duplicate oscillators during the search.

[0115] Second, when the AR model order is chosen as 2n-1 in step 1, a maximum of n oscillators is assumed. During subsequent search iterations after the first oscillator is added, fitting AR(2n-1) will again identify n additional frequencies. Since oscillators already exist in the model, a fixed AR model order does not strictly adhere to the upper limit on the number of oscillators. The AR model maintains a constant order to ensure that the oscillator parameters {a, f, σ 2}. When no further oscillations remain in the residual during subsequent search iterations, fitting an AR(2n-1) is similar to fitting a higher-order AR model to a white noise process. Scattering poles with small radii will be identified, which are still compatible with iOsc and model selection using log-likelihood because only a single oscillator is added at a time.

[0116] Decomposition of Oscillator Modeling

[0117] The non-iterative variant of the iOsc algorithm can be easily derived from the above routine with only minor modifications. Since the spectral decomposition using complex eigenvalues in step 3 is not limited to the largest poles / zeros, the required state noise variance parameter σ is initialized for any rotating oscillator after AR modeling and applying an independent mask observation vector. 2 Likewise, in step 2, the oscillator parameters {a,f} for all poles and zeros are accessible. Therefore, a non-iterative approach to model y can be constructed by initializing all n oscillators immediately after an initial fit of an AR model of order 2n-1 to y. This eliminates the need to check the OSPE in step 4, since all oscillators to be added have already been appropriately initialized.

[0118] We refer to this approach as decomposed oscillator modeling (dOsc) because it relies on the AR-based decomposition of the time series into complex eigenmodes from which the oscillator is initialized. The rest of the processing in dOsc is the same as in iOsc, with the oscillator initialized by σ 2 All initialized oscillators are added at once in descending order of . The only difference is that each oscillator is no longer initialized after the previous oscillator has been added and learned with EM; the initial parameters of all added oscillators are determined in one step by a single AR model fit.

[0119] As expected, given the overlap in underlying computations, iOsc and dOsc perform similarly in many cases. However, a key difference lies in how the two methods extract the underlying oscillations of the recording. iOsc is better suited to revealing masked oscillations that become clear only after accounting for other oscillations, because it operates on the innovation time series except for the first oscillator. This comes at the expense of potentially biasing subsequently added oscillators: this can be understood from the maximum likelihood properties of Kalman filtering and smoothing. When the number of oscillators used in the model is insufficient to account for the true number of underlying oscillators, lower-frequency oscillators tend to shift their center frequencies to account for unexplained activity at other frequencies. Consequently, the innovation spectrum can include more power in the frequency range in the presence of existing oscillators and less power in the absence of existing oscillators. Because the oscillator parameters in iOsc are initialized from the AR fit of the innovation spectrum, this can lead to biased initial estimates and parameter values after EM learning.

[0120] In contrast, because the oscillator parameters are initialized immediately on the raw data y, dOsc is not biased by oscillator fitting in previous iterations. Consequently, when clear and separable oscillations are present in the recording, dOsc generally recovers more accurate parameter and state estimates. However, dOsc cannot distinguish between closely spaced oscillations, particularly if there is a weaker oscillation nested beneath a stronger oscillation, which can only be modeled after the stronger oscillation is captured by its oscillator. In other words, the performance of dOsc is limited by how well AR modeling can separate the eigenmodes; dOsc still significantly outperforms AR modeling because it employs an interpretable and efficient model structure, allows for a more accurate representation of neural oscillations, and sequentially adds oscillators for model selection. However, if no oscillators are initialized at a certain value after AR fitting to y, dOsc does not have a mechanism to adapt and add oscillators in frequency regions that were not well captured by the previous oscillators (this requires examining the OSPE in iOsc), as it simply adds oscillators according to the initialization value σ. 2 Adds a set of oscillators initialized at the beginning in descending order.

[0121] When it comes to studying real recordings of neural activity with variable numbers of oscillations and uncertain oscillatory properties, dOsc and iOsc complement each other well. By applying both methods to recordings and comparing their modeling results, a clearer understanding of neural oscillations can often be achieved.

[0122] Observation noise initialization

[0123] The only remaining parameter that needs to be initialized before iOsc and dOsc can proceed is the observation noise variance R. In different recording modes, there are different ways to estimate the observation noise. For example, MEG can make use of empty room recordings to directly measure the observation noise process. However, in EEG, there is no equivalent empty room recording. Therefore, an empirical estimate is derived from the assumed white observation noise after high-pass filtering above the chosen frequency. The motivation for this heuristic estimate is that EEG recordings above a certain frequency are dominated by noise processes. For example, in the observation data y sampled at Fs = 100 Hz, the variance of y high-pass filtered above 30 Hz (by scaled to account for the full frequency range of the white noise) is used to obtain an initial estimate of R.

[0124] In updating the observation noise variance τ 2 We employ an informed inverse γ prior when . The motivation is that in the formulation of Equation (2), the oscillatory state noise process is additive with the observation noise. Therefore, during Kalman filtering and fixed-interval smoothing, some of the observation noise is explained by the hidden oscillator process, leading to a biased estimate of R during parameter learning. This can also be understood by noting that the theoretical spectrum of the oscillator has non-zero values above 30 Hz, and the additive effect of multiple oscillators in this frequency range can explain part of the white observation noise process during state estimation, leading to an underestimate of R. Therefore, during the iOsc and dOsc searches, the parameter updates of R are constrained by imposing a prior that places the mode at the above empirically estimated R using the variance of the high-pass filtered data.

[0125] To ensure comparable evaluation of the log-likelihood for model selection, since the log-likelihood expression involves τ 2 , so the same prior must be used throughout the iOsc iterations: in Given {y1,…,y t-1}, and is x t The conditional mean and variance of .

[0126] Model selection

[0127] Since oscillators can be iteratively added to model the observed data y given appropriate initialization, each iteration of iOsc and dOsc simply performs a certain number of iterations of EM learning to update the model parameters and estimate the hidden oscillator states. Model selection is based on the model, i.e., the combination of oscillators, that provides the best representation of y using the log-likelihood. Because more oscillators provide more degrees of freedom and may overfit y, the inflection point on the log-likelihood curve spanning from 1 to n oscillators is chosen to account for diminishing returns on additional model complexity (Satopaa et al., 2011). The model and its learned oscillators at the inflection point are the outputs of iOsc and dOsc as a set of oscillatory components underlying the observed data y.

Claims

1. A method for measuring intraoperative stimulation, the method comprising: receiving electroencephalogram (EEG) data from the patient using a patient monitor while applying one or more stimuli to the patient; extracting, using a computer processor, one or more event-related potential (ERP) signals from the EEG data associated with one or more discrete events; determining, using the computer processor, the one or more stimuli from the extracted one or more ERP signals; as well as The display is used to indicate adjustment of one or more intraoperative parameters associated with counteracting the one or more stimuli in the future one or more discrete events.

2. The method of claim 1 , further comprising decomposing the EEG data using a plurality of state-space oscillators derived from a background signal of the EEG data, the one or more ERP signals, and the one or more discrete events.

3. The method according to claim 2, characterized in that The waveform from the one or more ERP signals is extracted using a state-space ERP SS-ERP model defined by the following function: Among them, assuming n t ~N(0,τ 2 ) is independent, |f (m) |<f s , f s is the sampling frequency, a (m) <1, and is a state-space oscillator, and where x t,1 is the phasor x t The real component of x t,2 is the phasor x t The imaginary component, f is the frequency, f s is the sampling frequency, α is the damping factor, ∈ is the system state noise, σ 2 is the state noise variance, and n t is the measurement noise, τ 2 is the measurement noise variance.

4. The method according to claim 3, characterized in that The evoked response is represented as one or more ERP waveforms Presented with one or more stimuli t Convolution between: The elements of h follow a random walk model with a Gaussian driving noise v l The variance of 2 .

5. The method according to claim 2, characterized in that , the one or more ERP waveforms are extracted using the kernel-extended ERP KE-ERP model defined by the following function: Among them, assuming n t ~N(0,τ 2 ) is independent, |f (m) | <f s , f s is the sampling frequency, a (m) <1, and is a state-space oscillator, and where x t,1 is the phasor x t The real component of x t,2 is the phasor x t The imaginary component, f is the frequency, f s is the sampling frequency, α is the damping factor, ∈ is the system state noise, σ 2 is the state noise variance, and n t is the measurement noise, τ 2 is the measurement noise variance.

6. The method according to claim 4, characterized in that The evoked response is represented as one or more ERP waveforms Presented with one or more stimuli t Convolution between: Its expansion coefficient is determined by the Given, G is a Gabor dictionary or similar overlapping kernel / basis dictionary.

7. The method according to claim 1, characterized in that The one or more stimuli include a noxious stimulus or a painful stimulus.

8. The method according to claim 1, characterized in that The one or more stimuli include one or more electrical stimuli.

9. The method according to claim 1, characterized in that The ERP signal is associated with the one or more stimuli.

10. The method according to claim 1, characterized in that Adjusting the one or more intraoperative parameters includes adjusting the administration of one or more analgesics.

11. The method according to claim 1, wherein The EEG data includes one or more signatures of analgesics.

12. A system for measuring intraoperative stimulation, the system comprising: one or more sensors configured to measure electroencephalogram (EEG) signals of the patient; as well as a processor configured to: receiving an EEG signal from the patient while applying one or more stimuli to the patient; extracting one or more event-related potential (ERP) signals from the EEG data associated with one or more discrete events; determining the one or more stimuli from the extracted one or more ERP signals; determining an adjustment to one or more intraoperative parameters based on the one or more stimuli; and a display for communicating adjustments to the one or more intraoperative parameters.

13. The system according to claim 11, wherein: The processor is further configured to decompose the EEG data using a plurality of state-space oscillators derived from a background signal of the EEG data, the one or more ERP signals, and the one or more discrete events.

14. The system according to claim 11, wherein: A waveform from the one or more ERP signals is extracted using a state-space ERPSS-ERP model defined by the following function: Among them, assuming n t ~N(0,τ 2 ) is independent, |f (m) |<f s , f s is the sampling frequency, a (m) <1, and is a state-space oscillator, and where x t,1 is the phasor x t The real component of x t,2 is the phasor x t The imaginary component, f is the frequency, f s is the sampling frequency, α is the damping factor, ∈ is the system state noise, σ 2 is the state noise variance, and n t is the measurement noise, τ 2 is the measurement noise variance.

15. The system according to claim 13, wherein: The evoked response is represented as one or more ERP waveforms Presented with one or more stimuli t Convolution between: The elements of h follow a random walk model with a Gaussian driving noise v l The variance of 2 .

16. The system according to claim 11, characterized in that , the one or more ERP waveforms are extracted using the kernel-extended ERP KE-ERP model defined by the following function: Among them, assuming n t ~N(0,τ 2 ) is independent, |f (m) | <f s , f s is the sampling frequency, a (m) <1, and is a state-space oscillator, and where x t,1 is the phasor x t The real component of x t,2 is the phasor x t The imaginary component, f is the frequency, f s is the sampling frequency, α is the damping factor, ∈ is the system state noise, σ 2 is the state noise variance, and n t is the measurement noise, τ 2 is the measurement noise variance.

17. The system according to claim 15, wherein: The evoked response is represented as one or more ERP waveforms Presented with one or more stimuli t Convolution between: Its expansion coefficient is determined by the Given, G is a Gabor dictionary or similar overlapping kernel / basis dictionary.

18. The system according to claim 11, wherein: The one or more stimuli include a noxious stimulus or a painful stimulus.

19. The system according to claim 11, wherein: The one or more stimuli include one or more electrical stimuli.

20. The system according to claim 11, wherein The ERP signal is associated with the one or more stimuli.

21. The system according to claim 11, wherein: Adjusting the one or more intraoperative parameters includes adjusting the administration of one or more analgesics.

22. The system according to claim 11, wherein: The EEG data includes one or more signatures of analgesics.

Images