WO2024092280A1 - A system and method for measuring analgesia and responses to noxious or painful stimuli - Google Patents

Abstract

The present disclosure provides systems and methods for providing a quantitative assessment of pain and analgesia. The systems and methods provide, direct, objective, assessments of pain and the reduction of pain from different drugs or therapeutics. As a result, such systems and methods provide feedback for adjusting drug or therapeutics to a desired endpoint and could be used in the process of discovering or developing novel drugs or therapeutics. The system and methods include a novel state space evoked response potential (SS-ERP) model for extracting ERP signals from electroencephalography (EEG) data from a patient. The disclosed systems and methods enable the measurement of pain and pain therapeutic effects in an individualized manner.

Description

A System and Method for Measuring Analgesia and Responses to Noxious or Painful Stimuli Cross Reference to Related Applications [0001] The present application is based on, claims priority to, and incorporates herein by reference in its entirety for all purposes, US Provisional Application Serial No. 63/381,312, filed October 28, 2022. Statement of Government Support [0002] This invention was made with government support under R01DA056593 awarded by the National Institutes of Health. The U.S. government has certain rights in the invention. Background [0003] Each year more than 60 million surgeries are performed in the US, and more than 250 million worldwide. A considerable proportion of these patients may experience acute postoperative pain. Large-scale studies indicate that 26% of them experience moderate to severe acute pain after surgery., and nearly 10% will go on to develop chronic pain. Postoperative pain may also lead to post-operative opioid abuse, contributing to the opioid crisis. Accordingly, a critical priority in perioperative pain management is to avoid or minimize the risk of both acute and chronic postoperative pain. Postoperative pain also negatively impacts hospital and health-care systems by increasing returns to emergency department and readmissions, hospital length of stay and cost of care. [0004] Currently there are multiple strategies to provide analgesia during surgical procedures, including opioid-based strategies during general anesthesia, multimodal anesthesia, and regional anesthesia. Unfortunately, a major limitation underlying these strategies is the absence of monitors that can directly and objectively measure pain and/or nociception, and that can also assess whether an adequate level of analgesia is being maintained. Currently, anesthesiologists rely on changes in heart rate (HR) and blood pressure (BP) as indirect measures of nociception and titrate opioids and other drugs to minimize these changes. Unfortunately, HR and BP are influenced by numerous intraoperative factors, such as blood loss, anesthetic drugs, and anti-hypertensive medications, making them unreliable indicators of nociception. On the other hand, a significant fraction of surgeries rely on regional anesthesia techniques such as epidurals or regional blocks for nociception and post-operative pain management. However, epidurals and blocks are not foolproof, and may fail during surgery, leading to potentially uncontrolled 1  Q   B\125141.04438\85541799.2 postoperative pain that is apparent only when the patient recovers consciousness after surgery in intense pain. At present anesthesiologists have no way to directly assess the ongoing efficacy of their regional blocks, leading to unreliable nociception and pain management. Technologies have been proposed to monitor nociception during surgery using either HR or BP measurements or both, as well as 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) that can be measured at the scalp by averaging waveforms from repeated stimuli generated by a variety of methods including electrical stimulation. These ERPs are challenging to use in clinical settings because the responses are very small, on the order of between ~1 to 5 microvolts, compared to background EEG noise and oscillations that may be an order of magnitude larger in amplitude. To overcome the small effect size (or low signal-to-noise ratio) of these ERPS, a large number of stimuli must be provided, on the order of many 100’s or more, to average out the much larger background noise. Even then, the residual noise may be so large that it is difficult to obtain reliable ERP responses in individual subjects or patients. Accordingly, many research studies employing ERPs only show results in the aggregate across a population of subjects and cannot provide information on individual subjects. Differentiating between different conditions, e.g., the presence or absence or degree of a pain therapeutic, would be even more challenging. Even if all these difficulties could be overcome with a sufficiently high number of stimuli, practical clinical compliance challenges would remain: Could patients endure 10’s of minutes of repeated painful stimuli? Would the response to such stimuli ultimately saturate due to an underlying process of adaptation, or tolerance? Could the ERP stimuli, delivered over such a long time-frame, track changes in the nociception or pain during a procedure? Finally, these ERPs would be even more challenging if not impossible to use during anesthesia or when analgesic drugs such as opioids are administered, since the very small, between ~1 to 5 microvolts, and are overshadowed by anesthesia- or analgesia-induced electroencephalogram oscillations that are 10- to 100-fold larger in amplitude. Thus, although ERPs do have a relevant fundamental physiological basis for objective pain and analgesia measurement, in their present form the systems and methods available to compute ERPs make the ERPs completely unworkable in any practical setting, clinical or otherwise, particularly if individualized responses are required. 2  Q   B\125141.04438\85541799.2 Summary [0006] The present disclosure provides systems and methods that overcome the aforementioned drawbacks by providing a quantitative assessment of pain and analgesia, and in doing so offering a significant advantage over products that employ indirect measurements based on HR, BP, and other autonomic signals. Furthermore, the systems and methods provide, direct, objective, assessments of pain and the reduction of pain from different drugs or therapeutics. Such systems and methods provide feedback for adjusting drug or therapeutics to a desired endpoint, are useful in the process of discovering or developing novel drugs or therapeutics. The disclosed systems and methods enable the measurement of pain and pain therapeutic effects in an individualized manner. The disclosed systems also make it possible to evaluate pain in populations where classical assessment is difficult or not possible, such as in infants, non-verbal patients, or intubated and sedated patients. [0007] In one aspect of the present disclosure, a method for measuring intraoperative stimuli is described. The method comprises using a patient monitor to receive electroencephalography (EEG) data from a patient while one or more stimuli are applied to the patient. The method further comprises, using a computer processor, extracting one or more event-related potential (ERP) signals from the EEG data related to one or more discrete events. The method further comprises, using the computer processor, determining the one or more stimuli from the extracted one or more ERP signals and, using a display, indicating an adjustment of one or more intraoperative parameters correlated to counteract the one or more stimuli in a future one or more discrete events. [0008] In another aspect of the present disclosure, a system for measuring intraoperative stimuli is described. The system comprises one or more sensors configured to measure electroencephalogram (EEG) signals of a patient and a processor. The processor is configured to receive EEG signals from a patient while one or more stimuli are applied to the patient, extract one or more event-related potential (ERP) signals from the EEG data related to one or more discrete events, determine the 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 comprises a display to communicate the adjustment to the one or more intraoperative parameters. [0009] These aspects are nonlimiting. Other aspects and features of the systems and methods described herein will be provided below. 3  Q   B\125141.04438\85541799.2 Brief Description of the Drawings [0010] The foregoing features of embodiments will be more readily understood by reference to the following detailed description, taken with reference to the accompanying drawings, in which: [0011] FIG. 1 is a schematic illustration of a system for measuring pain and analgesia of a patient under anesthesia in accordance with the present invention. [0012] FIG. 2 is a is a flow chart setting forth the steps of a method for measuring pain and analgesia of a patient under anesthesia in accordance with the present invention [0013] FIG. 3 is a plot of a P300 oddball response. SS-ERP outperforms average ERP in extracting P300 in an oddball paradigm. The results are from only 52 trials (42 frequent + 10 oddball). The shaded regions show 90% CI constructed from posterior variances in SS-ERP and sample variances in avg. ERP. [0014] FIG.4A is a diagram of the experimental design for pain ERPs under analgesia. [0015] FIG. 4B is a plot of electrode Cz comparison of ERPs elicited by 4 different noxious stimulus intensity levels before administering Remifentanil. [0016] FIG. 4C is a plot of electrode Cz comparison of ERPs elicited by 4 different noxious stimulus intensity levels after administering Remifentanil. [0017] FIG. 4D is a plot of Fp2-Fp1 bipolar electrode comparison of ERPs elicited by 4 different noxious stimulus intensity levels before administering Remifentanil. [0018] FIG. 4E is a plot of Fp2-Fp1 bipolar electrode comparison of ERPs elicited by 4 different noxious stimulus intensity levels after administering Remifentanil. [0019] FIG. 5 is a receiver operating characteristics curve showing discrimination of low vs high electrical stimuli using SS-ERP P100 amplitude. AUC = area under curve. [0020] FIG.6A is a plot of group average results (n=10) for SS-ERPs under different stimulus intensity levels at baseline. [0021] FIG.6B is a plot of group average results (n=10) for SS-ERPs under different stimulus intensity levels and a mild analgesia regime. [0022] FIG.6C is a plot of group average results (n=10) for SS-ERPs under different stimulus intensity levels and a moderate analgesia regime. [0023] FIG.7A shows SS-ERP under general anesthesia. Simulated EEG time series using the oscillator parameters estimated from EEG recording from a sevoflurane induced anesthesia. [0024] FIG. 7B shows SS-ERP under general anesthesia. SS-ERP makes extraction of ERPs, dwarfed by the overlying slow and alpha, possible.7C 4  Q   B\125141.04438\85541799.2 [0025] FIG. 7C shows SS-ERP under general anesthesia. SS-ERP (MSE: 0.658 a.u) follows ground truth better than average ERP (MSE: 11.51). The shaded regions show 90% CI constructed from posterior variances in SS-ERP and sample variances in average ERP. [0026] FIG.8A shows a simulation study, where SS-ERP captures the time-locked responses while simultaneously extracting the background oscillations and thus removing their effects on the ERP waveforms in a data-driven manner. [0027] FIG.8B shows a simulation study, where SS-ERP performance remains unchanged as number of trials drops, unlike classical average ERPs. The background oscillations produce spurious peaks in the average ERPs with small number of trials while SS-ERP is relatively immune to such interference due to the explicit modeling of background oscillations [0028] FIG.8C shows a simulation study, where KE-ERP captures the time-locked responses while simultaneously extracting the background oscillations and thus removing their effects on the ERP waveforms in a data-driven manner. [0029] FIG.8D shows a simulation study, where KE-ERP performance remains unchanged as number of trials drops, unlike conventional average ERPs employing kernel expansions. The background oscillations produce spurious peaks in the average ERPs with small number of trials while KE-ERP is relatively immune to such interference due to the explicit modeling of background oscillations. [0030] FIG. 9A is a plot of oddball P300 response in healthy vs. Amyloid positive individual from 42 frequent/10 oddball trials. Traditional average ERP is confounded by the background oscillation. [0031] FIG. 9B is a plot of oddball P300 response in healthy vs. Amyloid positive individual from 42 frequent/10 oddball trials. SS-ERP shows the difference in ERPs, possibly due to underlying pathology. The shaded regions show 90% confidence interval constructed from posterior variances in SS-ERP and sample variances in averaging ERP. [0032] FIG. 9C is a plot of oddball P300 response in healthy vs. Amyloid positive individual from 42 frequent/10 oddball trials. A conventional average ERP employing a kernel expansion is confounded by the background oscillation. [0033] FIG. 9D is a plot of oddball P300 response in healthy vs. Amyloid positive individual from 42 frequent/10 oddball trials. KE-ERP shows the difference in ERPs, possibly due to underlying pathology. The shaded regions show 90% confidence interval constructed from posterior variances in KE-ERP and sample variances in averaging ERP. 5  Q   B\125141.04438\85541799.2 Detailed Description [0034] Disclosed herein is a novel technology for processing ERPs that increases their precision ~150-fold even in the presence of orders-of-magnitude larger background oscillations. The 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 to study psychophysiological correlates of sensory and cognitive processes with components that are informative of the course of sensory (‘exogenous’) and cognitive (‘endogenous’) processes with millisecond temporal resolution. ERPs are tiny ~1 μV signals that are embedded in background spontaneous oscillations that may be 10 to 100 times larger. Classically, combining a large number single- trial waveforms into an averaged ERP waveform is considered to be a “gold standard” for averaging out the event-independent activities. However, it is often difficult to obtain enough trials to extract these time-locked elicited responses, since a large number of trials are required to effectively average or cancel out the background of persistent oscillations. Also, excessive averaging could easily obscure fine structures of the ERPs, since intrinsic factors like mental workload, fatigue, attention are known to modulate time locked responses. To reduce this reliance on large number of trials, the method and system described herein models the background oscillations using a novel linear state-space model, so that the persistent, spontaneous activities are effectively removed. [0036] Herein, the human brain is modeled as a linear time-invariant system, consistent with classical ERP literature, so the response evoked by stimulus presentation manifested in EEG is written as a convolution between ERP waveforms, ^ ^^_^^^{^} ^, and an impulse train pertaining to discrete events of stimulus presentation, ^^_௧: ^ ^^{^^evoked^ ൌ ^ ^^ ∗ ൌ} _^ ^ℎ _^ [0037] However, unlike

stimulus-agnostic activities in EEG are modeled as a superimposition multiple oscillators, modeled after state-space oscillators: 6  Q   B\125141.04438\85541799.2 ^ ^{2 ^^ ^^^^^ 2 ^^^} ^ ^{^^^} ௧_,^ ^{é ^^ ^^ cos െ sin ù ^ ^ ^^} ^^ ^{^^ ê ^^ ^^ ^^} _௧,^ ^^{^^ ^^^^^ ^ ^ ú} _^ ^{௧ି^,^} _^ ^^_^:sampling imposed

on as or noise. In

summary, the models are used to work in tandem with ERP model to explain the EEG recording, effectively removing the contamination coming from strong oscillations: ^_{^ ൌ ^^} ^{^evoked^ ^^^^^௧.^} ௧ _{௧ ^ ^^௧ .} [0038] A prior is further included on the ERP waveform to ensure robust recovery and improve the variability of the estimates. Further, a temporal continuity constraint is imposed in form of following state-space model: ℎ_^ ൌ ℎ_^ି^ ^ ^^_^ , ^^_^ ∼ ^^^0, ^^^ଶ^, resulting in a multi-level state-space ERP (SS-ERP).

[0039] Alternatively, a different prior that assumes that the ERP, ^^ admits following expansion on a dictionary consisted of a given basis function, with a zero mean gaussian prior on the _{expansion coefficients,} ^{^} _{^^ with a diagonal Λ covariance matrix: ^^ ൌ ^^} ^{^} _{^^ ,} ^{^} _{^^ ~ ^^^ ^^,Λ^,} resulting in a state-space framework referred to as Kernel Expansion ERP (KE-ERP). [0040] Referring to FIG. 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, illustrated in FIG. 1 as an electroencephalography (EEG) electrode array. However, it is contemplated that the patient monitoring device may also include mechanisms for monitoring galvanic skin response (GSR), for example, to measure arousal to external stimuli. One specific realization of this design utilizes a frontal Laplacian EEG electrode layout with additional electrodes to measure GSR. Another realization of this design incorporates a frontal array of electrodes that could be combined in post-processing to obtain any combination of electrodes found to optimally detect the EEG signatures described earlier, also with separate GSR electrodes. Another realization of this design utilizes a high-density layout sampling the entire 7  Q   B\125141.04438\85541799.2 scalp surface using between 64 to 256 sensors 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. Also, cable 104 and similar connections can be replaced by wireless connections between components. As illustrated, the monitoring system 106 may be further connected to a dedicated analysis system 108. Also, the monitoring system 106 and analysis system 108 may be integrated. The patient monitoring device 102 and monitoring system 16 may also each include a user interface 110. The user interface 110 may be, but is not limited to, a keyboard, mouse, touch screen, or buttons, switches integrated in the monitoring system 106. [0042] For example, as noted above, it is contemplated that the patient monitoring device 102 may be an EEG electrode array, for example, a 64-lead EEG electrode array. However, as will be apparent below, greater spatial accuracy can be achieved by increasing the number of electrodes from 64 to 128, 256, or even higher. Similarly, the present invention can be implemented with substantially fewer electrodes. In any case, the monitoring system 106 may be configured to receive raw signals acquired by the EEG electrode array and assemble, and even display, the raw signals as EEG waveforms. Accordingly, the analysis system 108 may receive the EEG waveforms from the monitoring system 106 and, as will be described, analyze the EEG waveforms and signatures therein based on a selected anesthesia compound, determine a state of the patient based on the analyzed EEG waveforms and signatures, and generate a report, for example, as a printed report or, preferably, a real-time display of signature information and determined state. However, it is also contemplated that the functions of monitoring system 106 and analysis system 108 may be combined into a common system. [0043] Referring to FIG.2, a method 200 of measuring pain and analgesia of a patient under anesthesia is shown. The method may be performed by a computer processor associated with the monitoring system 106 of FIG.1. 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 a non-limiting example, the one or more stimuli may include an electrical shock, or an audible or visual stimulus. An electrical shock may be provided by one or more electrodes attached to the patient. The electrical shock may be variable to adjust the voltage of the electrical shock presented to the patient. 8  Q   B\125141.04438\85541799.2 [0045] In a non-limiting example, the EEG data includes a signature of one or more analgesic agents. Furthermore, the EEG signals are acquired by one or more sensors, such as patient monitoring device 102. [0046] Next, at step 204, one or more ERP signals is extracted from the EEG data related to one or more discrete events. Further, the one or more stimuli are determined from the one or more ERP signals at step 206. In a non-limiting example, the one or more stimuli include noxious stimuli 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, the processor of the monitoring system decomposes the EEG data using a plurality of state-space oscillators derived from background signal for the EEG data, the one or more ERP signals, and the one or more discrete events. Includes a state pace ERP (SS-ERP) model for extracting the one or more ERP signals from the EEG data. Further, the SS-ERP model includes a plurality of state space oscillators derived from background signal for the EEG signals, the one or more ERP signals, and the one or more discrete events. [0048] In a non-limiting example, the one or more ERP signals are separated from oscillatory background using following composite model: ^^{^^^^ cos 2 ^^ ^^^^^/ ^^ െ si ^^^ ^^^^^ ^^^^^} ^ ^{௧,^} ^_{^^൩ ൌ ^^} ^{^^^} _^ ^{^ n 2 ^^ ^^ / ^^^ ௧ି^,^} ^_{^^ ^^^ ^ ^ ^ ൩ ^ ^} ^{௧,^} _{൩ , ^^ ൌ 1, 2,⋯ , ^^ ^^ ^^ ^^^ ௧,ଶ}

where ^^^{^^^} ௧_,∗ ~ ^^^0, ^^_^ ^ଶ ^, ^^_௧~ ^^^0, ^^^ଶ^ are assumed to be independent, | ^^^{^^^}| ^ ^^_^, ^^_^ , fs is the ^ 1, and ^^^{^^^} ௧_,^, … , ^^^{^^^} ௧_,^ are the state space oscillators, and where

and imaginary

phasor representation of constituent oscillation, ^^_௧, ^^ is the oscillation frequency, ^^_^ is the sampling frequency, ^^ is the damping factor, ^^ is oscillation state noise, ^^^ଶ is the state noise covariance, ^^_௧is the measurement noise and ^^^ଶ is noise variance The evoked responses are represented as a convolution between ERP waveforms, ^ ^^_^^^{^} ^, and one or more stimulus presentations, ^^_௧ : ^ . The ERP models are

[0049] Alternatively, the evoked component can be represented as an expansion on a dictionary consisting of a given basis function. As a representative of such basis functions, a 9  Q   B\125141.04438\85541799.2 Gabor dictionary, ^^, is considered whose atoms are given by Gaussian kernels of fixed variance, normalized to have maximum 1. ^_{^ ൌ ^^} ^{^} _{^^ ,} ^{^} _{^^ ~ ^^^ ^^,Λ^ A zero mean gaussian prior is considered on the expansion coefficients,} ^{^} _{^^ with a diagonal Λ} covariance matrix. The ERPs are referred to as Kernel Expansion ERP (KE-ERP). [0050] At step 208, an adjustment of one or more intraoperative parameters correlated to counteract the one or more stimuli in a 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 an administration of one or more analgesic agents. [0051] As used in this specification and the claims, the singular forms “a,” “an,” and “the” include plural forms unless the context clearly dictates otherwise. [0052] As used herein, “about”, “approximately,” “substantially,” and “significantly” will be understood by persons of ordinary skill in the art and will vary to some extent on the context in which they are used. If there are uses of the term which are not clear to persons of ordinary skill in the art given the context in which it is used, “about” 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 meaning as the terms “comprise” and “comprising.” The terms “comprise” and “comprising” should be interpreted as being “open” transitional terms that permit the inclusion of additional components further to those components recited in the claims. The terms “consist” and “consisting of” should be interpreted as being “closed” transitional terms that do not permit the inclusion of additional components other than the components recited in the claims. The term “consisting essentially of” should be interpreted to be partially closed and allowing the inclusion only of additional components that do not fundamentally alter 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 analogous to “at least one of A, B and C, etc.” is used, in general such a construction is intended in the sense of one having 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 not be limited to systems that have A alone, B alone, C alone, 10  Q   B\125141.04438\85541799.2 A and B together, A and C together, B and C together, and/or A, B, and C together.). It will be further understood by those within the art that virtually any disjunctive word and/or phrase presenting two or more alternative terms, whether in the description or figures, should be understood to contemplate the possibilities of including one of the terms, either of the terms, or both terms. For example, the phrase “A or B” will be understood to include the possibilities of “A” or “B” or “A and B.” [0056] All language such as “up to,” “at least,” “greater than,” “less than,” and the like, include the number recited and refer to ranges which can subsequently be broken down into ranges and subranges. A range includes 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, and so forth. [0057] The modal verb “may” refers to the preferred use or selection of one or more options or choices among the several described embodiments or features contained within the same. Where no options or choices are disclosed regarding a particular embodiment or feature contained in the same, the modal verb “may” refers to an affirmative act regarding how to make or use an aspect of a described embodiment or feature contained in the same, or a definitive decision to use a specific skill regarding a described embodiment or feature contained in the same. In this latter context, the modal verb “may” has the same meaning and connotation as the auxiliary verb “can.” Example [0058] The state space event-related potential (SS-ERP) and kernel expanded event-related potential (KE-ERP) methods can separate ERPs from background oscillations with up to ~150x greater precision 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 waveforms between ~1 and ~10 microvolts in size that occur over 10’s to 100’s of milliseconds. The major drawback of current approaches is that they rely on averaging over hundreds of trials to extract this time-locked elicited responses from a backdrop of persistent slow and/or alpha wave oscillations. Consequently, the structure of these conventional averaged ERPs is easily obscured by these background oscillations, particularly if they 10- or 100-fold larger as they are during general anesthesia and sedation. However, as described herein, an “SS-ERP” method that uses multi-level state space models, and an “KE-ERP” method that uses sparse representation on a pre-defined dictionary of bases or kernels to improve ERPs by explicitly modeling and effectively removing persistent, spontaneous background activity. The method further improves the 11  Q   B\125141.04438\85541799.2 precision of ERPs by imposing a temporal continuity prior on the ERP waveform to recover a denoised ERP. [0059] The performance of the SS-ERP was benchmarked by analyzing the P300 auditory “oddball” response, an ERP that is generated by comparing responses to frequently presented auditory tones (e.g., 2000 Hz) with those of less frequent “oddball” tones (1000 Hz). Several hundred trials are typically required to resolve ERP waveforms. FIG.3 shows the results from the SS-ERP algorithm, working with only 52 trials (42 frequent, 10 oddball) in a single subject. Whereas the difference between the frequent and oddball waveforms is difficult to appreciate with traditional averaging (Cohen’s d: 3.57, precision (inverse of the estimation variance): 0.56, minimally-detectable event (MDE, 2.485*std): 3.75 µV), the SS-ERP- derived waveforms are readily appreciated with clear differences between frequent and oddball conditions, and with narrow confidence intervals (Cohen’s d: 21.27, precision: 85.33, MDE: 0.303 µV). This analysis showcases an impressive improvement in effect size (Cohen’s d) of 5.96-fold, a reduction in the MDE of 12.37-fold, and an enormous 152.4-fold increase in precision using the SS-ERP in a single subject from only 52 trials. [0060] We ran the same benchmark for KE-ERP, but with different numbers of trials (118 frequent, 30 oddball) in a single subject. Like SS-ERP, KE-ERP (cohen’s d: 167.68, precision: 178.98, MDE: 0.2µV) shows impressive improvements over classical averaging ERP (cohen’s d: 5.10, precision: 0.112, MDE: 7µV). [0061] ERP features correlated with pain stimulus intensity are diminished by opioids to a level equivalent to non-painful stimuli; these features can be estimated from as few as 30 trials from parietal or frontal electrodes using the SS-ERP method. ERPs elicited by painful somatosensory stimuli can be clearly seen in EEG channels directly placed over somatosensory cortex, such as Cz, but suffer from the same problems described above for the P300 oddball response: the pain ERPs are orders of magnitude smaller than overlying background EEG signals. [0062] Further, EEG recordings were analyzed from an individual healthy volunteer who received electrical stimulation at four different intensities calibrated to the subject’s pain perception so that only the highest stimulus intensity, “Level 4,” was perceived to be painful. The stimulation consisted of 10 train-of-3 electrical pulses for each intensity level, shuffled randomly (FIG.4A). Remifentanil (3ng/mL) was then administered and then 5 train-of-3 electrical pulses (shuffled randomly) were presented at each pain intensity level. Similar to the P300, the SS-ERP method significantly improved the precision of pain ERP estimates compared to traditional averaging (FIGS.4B-4E). The waveforms from the Cz electrode 12  Q   B\125141.04438\85541799.2 (FIGS.4B-4C) show visible differences between stimulus intensities, but with far higher precision when using the SS-ERP method. With the SS-ERP method it is easy to see that ERP waveforms diminish after remifentanil administration (FIG 4C). In contrast, traditional averaging shows considerable residual structure in the ERP with large overlapping confidence limits (FIG.4B). [0063] Detection of pain ERPs were considered in this individual subject over prefrontal electrodes Fp1 and Fp2, which could be more accessible for clinical applications as they sit on the forehead, and which might be more informative about affective perception of pain due to the proximity of the prefrontal cortical structures that mediate that perception. Using traditional averaging (FIG.4D), ERPs cannot be effectively estimated from Fp1-Fp2 and no differences can be discerned between stimulus intensity levels or when the opioid is administered. In contrast, the SS-ERP method can detect clear differences between stimulus intensity levels at baseline (FIG.4E). Furthermore, using SS-ERP, when remifentanil is administered the responses to the painful stimulus “Level 4” becomes indistinguishable from the lower intensity stimuli (FIG.4E). [0064] FIG.5 shows the receiver operating characteristics (ROC) curve discriminating low vs high electrical stimulation using the P100 ERP peak during the baseline period prior to analgesia. The AUC for discriminating low vs. high electrical stimulation is 0.86, indicating highly accurate discrimination in an individual subject. This analysis suggests that the methods described here can be used for highly accurate diagnostic inferences in individual patients with only a small number of stimulus events. [0065] Next it was verified whether, on average, these same results would hold over a larger group of subjects. FIGS.6A-6C show the group average results for n = 10 subjects who underwent the same protocol presenter earlier, in which four different electrical stimulus intensities were presented, tuned so that Intensity level 4 was subjectively rated by the subject as being painful. The subjects were studied under a baseline condition (FIG.6A), as well as mild (FIG.6B) and moderate analgesia (FIG.6C) conditions provided by intravenous infusions of remifentanil. The group average results replicate those observed in the individual subject: the SS-ERP response increases in amplitude with increasing stimulus intensity and decreases progressively with administration of mild and moderate analgesia, respectively. [0066] The diagnostic value of the SS-ERP signal for predicting the level of pain was further considered. To summarize the SS-ERP response, in one example, the amplitude between the initial trough of the response to the subsequent peak at the central EEG channel Cz are 13  Q   B\125141.04438\85541799.2 analyzed, which are denoted as a variable “SS-ERP” in the statistical models. A mixed-effect model is built to evaluate the association between the SS-ERP amplitude and the pain reported by each subject: Pain ~ SS-ERP + (SS-ERP | Subject). [0067] It was found that the SS-ERP amplitude was positively correlated with the pain reported by each subject at every single 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 could influence the pain reported, a second mixed-effect model was built, where now two covariates (level of analgesia administered and the intensity of the stimuli) are included: Pain ~ SS-ERP + Analgesia + Stim + (SS-ERP | Subject). [0069] The fixed variables were the amplitude of the ERP obtained with the method (SS- ERP), the amount of opioid administered (Analgesia) and the intensity of the stimuli (Stim). The random variable was the subject. Chi² p-value SS-ERP amplitude 69 0008

[00 0] emar aby, a ter adjust ng by t e ntens ty o t e st mu , and the amount of analgesia received, the amplitude of the SS-ERP is still associated with the pain reported by subjects, indicating that the method is able to capture the individual variation between subjects in their response to stimulus intensity and analgesia, i.e., the SS-ERP is measuring pain information individualized to each specific subject. When repeated, the analysis using the difference between the frontal EEG electrodes, Fp2-Fp1, the statistical significance for the association between the SS-ERP and pain is even higher (p-value = 8.9e^-5), consistent with the notion that frontal electrodes may provide more information about affective components of pain mediated by prefrontal cortical structures. [0071] These data show how the SS-ERP method can detect differences in pain stimulus intensity as well as diminished ERP waveforms in response to administration of a pain therapy (opioid analgesia). Moreover, these data illustrate how these differences in pain intensity and analgesia may be detected using only 15 to 30 stimulus events, implying that this methodology could be applied practically in clinical scenarios to objectively measure pain in a short duration of time. 14  Q   B\125141.04438\85541799.2 [0072] The SS-ERP method can recover ERPs accurately against a background of ~10-fold larger anesthesia-induced oscillations. During general anesthesia and sedation, the EEG shows large, stereotyped oscillations at specific frequency bands. Propofol and sevoflurane, the two drugs most frequently used to maintain unconsciousness during general anesthesia, induce slow and delta oscillations (~1-4 Hz) alongside alpha oscillations (~9-12 Hz). These slow/delta and alpha oscillations are much larger than the ~1 to ~10 microvolt ERP changes expected during somatosensory stimulation. In this scenario, extracting ERP responses using traditional averaging methods would require a massive number of stimulus events to overcome the substantially larger noise. This would make ERPs for pain assessment impractical in a clinical setting due to an unacceptably long latency in assessing the ERP measurement. In addition, the overly long train of stimulus events could cause discomfort to some patients. To verify that SS-ERP can effectively extract ERPs in the presence of large anesthesia-induced background oscillations, a realistic simulation study was performed featuring simulated ERP waveforms obscured by much larger simulated anesthesia-induced oscillations. The background oscillations were generated from state space oscillator models whose parameters were estimated from a patient undergoing sevoflurane general anesthesia. As shown in FIG.7A, the simulated data are similar in form, both in time and frequency domain, to the real data. An event related signal, generated by replicating a train of 30 realistic ERP waveforms, was added to the simulated oscillations with an SNR of ~ 1/10. FIG.7B shows the decomposition of that signal into oscillatory and event related components. The SS-ERP method is clearly able to extract the tiny ERP signals from the background of strong oscillations. Compared to the averaged ERP (mean squared error, MSE: 11.51 a.u.), the SS-ERP (MSE: 0.658 a.u.) follows the ground truth more faithfully (17.5-fold reduction in MSE) and does so with a much tighter confidence interval (FIG.7C). These results illustrate how the SS-ERP method can estimate ERP waveforms even in the presence of much larger anesthesia-induced EEG oscillations. [0073] Methods Overview [0074] Novel ERP analysis method. [0075] State-space ERP. The SS-ERP algorithm employs a multi-level model that represents persistent oscillatory background signals distinct from evoked responses, making it possible to separate these two components and dramatically reduce the uncertainty in the ERP estimate. The background oscillations, ^^^{^^^ ^^^} ௧_,^, … , ^^_௧,^ , are represented by state space oscillators. The event-related time series are modeled as a convolution between ERP 15  Q   B\125141.04438\85541799.2 waveforms, ℎ, and discrete events, ^^ (modeled as impulses), and the observed EEG signal is then represented as the sum of background and ERP components with an observation noise term: ^^{^^^^ ^^^ ^^^ ^ ^^^ ^^^} ^ ^{௧,^} ^_{^^൩ ൌ ^^} ^{^^^} _^ ^{cos 2 ^^ ^^ / ^^^ െ sin 2 ^^ ^^ / ^^^} s_{in ^^^ ^^^ ^ ^} ^{^௧ି^,^} ^_{^^ ൩ ^ ^} ^{^^௧,^} ^_{^^൩ , ^^ ൌ 1, 2,⋯ , ^^ (A) ^^ 2 ^^ ^^ / ^^^ cos 2 ^^ ^^ / ^^^}

௧_,∗ ~ are frequency, and ^^^{^^^} ^ 1. This multi-level framework also allows us to include a prior on the ERP waveform to ensure robust recovery and reduce the variability of the estimates, implemented here in the form of a random walk prior: ℎ_^ ൌ ℎ_^ି^ ^ ^^_^, ^^_^ ∼ ^^^0, ^^^ଶ) (C) [0076] ^^^{^^^} ^^^ଶ,

and smoothness parameter, ^^^ଶ are learned using an Maximization

(EM) algorithm. EM alternates between optimizing the distribution over the hidden oscillator states and ERP waveform given the current parameters (the E-step) and updating the parameters given the distribution of hidden states (the M-step). The E-step was further simplified by constraining the posterior distributions of ERP waveform, ℎ and oscillation states to be statistically independent: ^^^ ^^, ^ ^^_௧^ ∣ ^^, ^^^ଶ, ^^^ଶ^ ൌ ^^^ ^^ ∣ ^^^ଶ, ^^^ଶ^ ^^^^ ^^_௧^ ∣ ^^, ^^^ଶ^, to be able to cyclically update the oscillator state and the ERP distribution. The different model parameters can then be re-estimated in closed from in the M-step. This approach is known as Variational Bayes approximation, and the modified algorithm is known as generalized EM. Specifically, the ERP update step resembles ridge regression, with the estimated oscillations removed from EEG signal. It is noted that in contrast to point estimates in classical averaging technique, the framework provides the posterior distribution of the ERP. The updates are stopped when the value of the likelihood of the observation stabilizes. [0077] Kernel Expanded ERP. The model is utilized in equation (A) with the assumption that the ERP admits following expansion on a dictionary formed of a given basis function. As a representative of such basis functions, Gabor dictionary, ^^, was considered, whose atoms are given by Gaussian kernels of fixed variance, normalized to have maximum 1. ^_{^ ൌ ^^} ^{^} _{^^ ,} ^{^} _{^^ ~ ^^^ ^^,Λ^} 16  Q   B\125141.04438\85541799.2 _{We again consider a zero mean gaussian prior on the expansion coefficients,} ^{^} _{^^ with a} diagonal Λ covariance matrix. ^{[0078] Learning the ERP. Similar to SS-ERP, Oscillator parameters, ^^^^^} _{:ൌ ൫ ^^^^^, ^^^^^, ^^^ଶ ൯, noise variance, ^^ଶ, and expansion coefficient prior variance, Λ are learned} an instance of Expectation Maximization (EM) algorithm. EM alternates between

over the hidden oscillator states and ERP waveform given the current parameters (the E-step) and updating the parameters given the distribution of hidden states (the M-step). The E-step was further simplified by constraining the posterior distributions of ERP waveform expansion coefficients,^{^} ^^ and oscillation states to be statistically independent: ^_{^൫ ^} ^{^} _{^, ^ ^^௧^ ∣ ^^, ^^} ^ଶ _{,Λ൯ ൌ ^^൫ ^} ^{^} _{^ ∣ ^^} ^ଶ _{,Λ൯ ^^^^ ^^௧^ ∣ ^^, ^^} ^ଶ _^, to be able to The different

model parameters can re- This learning of Λ, in fact, encourages the coefficients of irrelevant lags to drive to zero, thereby, performing an automatic relevant detection (ARD) in a completely data-driven manner. [0079] In summary, SS-ERP and KE-ERP readily capture the time-locked responses while simultaneously extracting the background oscillations and thus removing their effects on the waveforms in a data-driven manner. Also, in contrast to point estimates in the traditional averaging technique, both the frameworks provide the posterior distribution of the ERP. [0080] Simulation Results [0081] FIGS.8A-8B illustrate the SS-ERP extraction from simulated auditory tone responses, contaminated by strong slow and alpha oscillations. FIG.8A depicts how the oscillatory components are effectively removed to retrieve the auditory evoked responses by explicit modeling of the oscillatory activities. FIG.8B compares the SS-ERP to the traditional average ERP alongside the ground truth for increasing number trials (n=25,50,100). Clearly, SS-ERP follows the ground truth very closely in all three cases, with the narrow confidence intervals deteriorating slightly as n decreases. On the contrary, average ERP is dominated by the large oscillatory components present in the background: meaningful results are only obtained after averaging 100 trials. This demonstrates that the SS-ERP is superior at capturing the true latencies and amplitudes of the ERP peaks from small number of trials. 17  Q   B\125141.04438\85541799.2 [0082] FIGS.8C-8D illustrate the KE-ERP extraction from the similar example. They correspond to the FIGS 8A-B respectively. Clearly, the observations made for SS-ERP also applies here. [0083] Application on real EEG data [0084] SS-ERP performance was also benchmarked in detecting P300 using EEG data (single channel, Parietal lobe) recorded from a single subject: the subject (eyes closed) was asked to respond to a target tone (2000 Hz) that occurs 20% of the time intermixed with a regular tone (1000 Hz). Only the first 52 out of a total of 200 trials were included in the analysis for to demonstrate the improved ERP extraction with SS-ERP method: even with 52 trials which barely allow traditional averaging methods to identify the P300 response (effect size: 3.57, Cohen’s d, precision: 0.56), SS-ERP captures the response with high confidence by removing the slow wave and persistent alpha from the background (effect size: 21.27 Cohen’s d, precision: 85.33). This showcases an impressive improvement in effect size of 5.96 folds, and in precision by 12.37 with only 52 trials. In addition, SS-ERP clearly shows diminished P300 response in the amyloid positive individual compared to healthy individual (FIG.9B), while the strong background oscillation activity confounds such effect in traditional averaging ERPs (FIG.9A). [0085] We ran the same benchmark for KE-ERP, but with different numbers of trials (118 frequent, 30 oddball) in a single subject. Like SS-ERP, KE-ERP (Cohen’s d: 167.68, precision: 178.98, MDE: 0.2µV) shows impressive improvements over classical averaging ERP (Cohen’s d: 5.10, precision: 0.112, MDE: 7µV). KE-ERP again shows diminished P300 response in the amyloid positive individual compared to healthy individual (FIG.9D), while the strong background oscillation activity confounds such effect in traditional averaging ERPs (FIG.9C) [0086] Thus, while significantly reducing the number of trials requirements, the proposed method allows tracking of short-term changes in ERP due to various intrinsic and extrinsic reasons, which has significant implications for neuroscience studies and clinical applications. [0087] Learning Parameters for State Space Oscillator Models [0088] Given a univariate time series data y from neural signal recordings, a key objective is to determine what, if any, neural oscillations are present in the data y. A related question is to characterize the properties or parameters of those neural oscillations. The iterative oscillator method (iOsc) addresses these questions. [0089] iOsc is a greedy search algorithm that attempts to represent ^^ by adding neural oscillation components one at a time until a pre-specified stopping number. Model selection 18  Q   B\125141.04438\85541799.2 is then performed based on some metric of how well different numbers of oscillations can represent the original data. The output of the iOsc algorithm is this selected set of neural oscillations, inclusive of the number of oscillations and all of their parameters, which are regarded as the most representative oscillator components underlying ^^. [0090] ^^ is represented in iOsc is via time domain modeling using the following class of parametric oscillator state-space models: _^ ^{^^௧,^ co ^^ ^^ ^^ ଶ} ^_{^௧,ଶ^ ൌ ^^ ^} ^{s2 ^^ ^^/ ^^^ െsin2 ^^ ^^/ ^^^} s_{in2 ^^ ^^/ ^^ ൨ ^} ^{௧ି^,^} _{^ ^ ^} ^{௧,^} _{^ , ^} ^{௧,^} _{^ ∼ ^^ ^ ^^, ^} ^{^^ 0} _{, (1); ^ cos2 ^^ ^^/ ^^^ ^^௧ି^,ଶ ^^௧,ଶ ^^௧,ଶ 0 ^^ଶ൨^} . it

random frequency modulation, which is unanimously observed in human EEG recordings, distinct from pure sinusoidal waves encountered in physical science problems such as radar. When there are multiple neural oscillations present at the same time, multiple independent oscillators are easily concatenated in block-diagonal matrix form in equation (1) and arrive at an equally interpretable model with simple summations in equation (2) to produce ^^. An additional benefit of time domain modeling is that log likelihoods provide a natural metric to evaluate how well different models with different numbers of oscillators can represent ^^, allowing for a straightforward model selection scheme. [0092] At each iteration, iOsc performs time domain modeling of ^^ using one or more oscillator components in the form of equation (1). The unknown parameters are ^ ^^, ^^, ^^^ଶ^ for each oscillator and a single observation noise variance ^^^ଶ. If the values of these parameters are known, state estimation (also called exact inference or belief propagation) can be done efficiently using Kalman filtering and fixed-interval smoothing algorithms to estimate the hidden variables ^^, which represent the time series of the underlying oscillation components. Conversely, given values of ^^, one can update the parameter values in this linear Gaussian setting using closed-form equations. Alternating between these two steps is an instance of expectation-maximization (EM) algorithm with oscillators ^^ and data ^^. Detailed equations for exact inference and parameters are outlined below. [0093] While each iteration of iOsc simply carries out EM learning, and the computations are known and efficiently implemented, a crucial challenge is how to initialize EM in iOsc. The complete data log likelihood landscape in this linear Gaussian case is concave with a single global maximum; however, there are non-trivial numbers of unknown parameters and only 19  Q   B\125141.04438\85541799.2 limited numbers of time points in ^^. Given poor initialization, EM may require infinite iterations or data length to converge to the appropriate values and render poor results of iOsc. Thus, a crucial step to successfully identify underlying oscillator components using iOsc is to choose informed starting points throughout the greedy search. [0094] Oscillator initialization [0095] The algorithm design of iOsc adds one oscillator at a time; therefore, a scheme to initialize model parameters is constructed for a single oscillator. Consider that data ^^ was observed under sampling rate Fs, and the pre-specified stopping number of oscillations is ^^, meaning there are at most ^^ oscillations present in the recording. For example, ^^ ൌ 7 could be a reasonable but non-limiting choice in human EEG that typically contains 2-4 distinct oscillations. The following steps to initialize an oscillator. [0096] Step 1: fit an autoregressive (AR) model [0097] An AR model of order 2 ^^ െ 1 is first fit to ^^. The order is chosen as 2 ^^ െ 1 because the complex roots of AR models appear as pairs in frequency. For example, when ^^ ൌ 7, a zero at DC and 6 frequencies above 0 Hz for oscillatory poles using an AR(13) model are identified. This AR fitting is done using the Yule-Walker method or the Burg algorithm, for example, and the fitted AR coefficients are denoted as ^^ ൌ ^ ^^_^, ^^_ଶ,⋯ ^^_ଶ^ି^^: ^^_௧ ൌ ∑^ଶ ^^{^} ୀ^ି ^^{^} ^^_^ ^^_௧ି^ ^ ^^_௧ , ^^_௧~ ^^^0, ^^^ଶ^ (3). [0098] Step 2: select the largest pole/zero [0099] The complex roots of the fitted AR model can be identified by solving for the roots of the polynomial with coefficients ^1,െ ^^^. The theoretical PSD spectrum ^^ of the AR process can also be obtained. To initialize the oscillator parameters ^ ^^, ^^^, the root corresponding to the highest power is first found. The oscillator strength at a pole is then estimated by weighting the theoretical PSD ^^ at the pole frequency ^^_^ by the pole radius ^^_^: ^^̂ ൌ arg max_^ ^^൫2 ^^ ^^_^^/ ^^_^൯ ^^_^, (4); where ^^ comes from the set of all poles/zeros with unique oscillatory frequencies using the complex roots. Once the pole/zero with the largest power is identified, the oscillator parameters ^ ^^, ^^^ can be initialized as: ^^ ൌ ^^_{^^ ^^} When a zero is selected as the largest

0 Hz is not a desired choice because these oscillator parameters are initialized for EM learning. Kalman filtering starting from this value prevents updates of the oscillator frequency to be different from 0 Hz. Thus, a lower 20  Q   B\125141.04438\85541799.2 bound is placed on ^^ at 0.1 Hz during initialization, although subsequent EM learning can update it to be closer to 0 Hz. [0100] Step 3: initialize oscillator state noise covariance [0101] After fitting an AR model, there is not an obviously apparent way to initialize the state noise variance parameter ^^^ଶ for a single oscillator. This is because the estimated white noise process variance ^^^ଶ excites all poles and zeros in the AR process at the same time, instead of being specific to one or multiple poles at the frequency ^^ where one might wish to place an oscillator. This technical challenge can be solved by transforming the AR process and the noise variance ^^^ଶ into distinct eigenmodes and read off noise variances only from the poles that locate at the selected frequency. [0102] The key idea here is that an informed estimate of ^^^ଶ is derived using spectral decomposition of the fitted AR process. The AR(2 ^^ െ 1) model is written in the following state-space form: ^^{^௧ ^^^ ⋯ ^^ଶ^ିଶ ^^ଶ^ି^ ^^௧ି^} ^^ ^^ ^ _൦ ^ଶ ൦ _௧ି^ ൪ ൌ ^^{1 ⋯ 0 0} _௧ିଶ ൪ ^{^ ∼ ^^ ^ ^ ^^ 0 ;} .

^^{^} _௧ ^{^^} ^^ _ି^ ^{௧ି^} Therefore, the AR(2

form: ^^_௧ ൌ ^^ ^^ _^^െ1 ^ ^^ _^^ (7). ^^_௧ ൌ ^^ ^^ _^^ (8).

^^ ൌ ^^ ^^ ^^^{ି^} (9). Equation (7), top line can be rearranged by multiplying ^^^{ି^}, an orthogonal matrix, to both sides: ^^^{ି^} ^^_௧ ൌ ^^ ^^^{ି^} ^^_௧ି^ ^ ^^^{ି^} ^^_௧ (10) 21  Q   B\125141.04438\85541799.2 [0103] The diagonal matrix ^^ contains complex eigenvalues, which are actually the same poles and zeros of the characteristic function with AR coefficients ^ ^^^. This equivalence relation exists because ^^ is the companion matrix of the polynomial equation with coefficients ^1,െ ^^^. In fact, modern computing software solves the polynomial via eigendecomposition of the companion matrix. Therefore, not only the complex eigenvalues in ^^ match the complex roots in step 2, they are in the exact same order. [0104] Notice that ^^ can be viewed as the transition matrix of a set of parametric oscillators. Each complex number can be written in a complex exponential form, which coincides with a rotation matrix with frequency ^^ multiplied by a damping parameter ^^ that is the magnitude of the complex eigenvalue (or equivalently, the radius of the complex root). In other words, except a slight difference in the state noise processes not being diagonal, equation (10) is also describing a set of block-diagonally concatenated oscillators in the space rotated by ^^^{ି^} just like equation (1). The only mathematical distinction lies in that the noise processes driving the real and imaginary components of each oscillator are now correlated, as well as among different oscillators. While this oscillator interpretation of eigenmodes is not further utilized, it provides the motivation for initializing the oscillator parameter ^^^ଶ using variances of the noise processes driving these eigenmodes. [0105] Observation equation in (8) is also modified to be in the same rotate space: ^^_௧ ൌ ^^ ^^ ^^^{ି^} ^^_௧ , with the new observation vector ^^′ ൌ ^^ ^^, which is simply the first row of ^^, or ^^_^,⋅ in short: ^^_௧ ൌ ^^_^,⋅ ^^^{ି^} ^^_௧. (11). [0106] Since ^^_^,⋅ is taking a linear combination of the state values from the different oscillators

and zeros) in the rotated space to form the observed data ^^_௧, all entries of ^^_^,⋅ are masked out as 0 except for the selected pole by equation (4) : ^_{^௧ ൌ ^} ^{0 ⋯ ^^} _^,^^ ^{⋯ 0} _{^ ^^} ^{ି^} _{^^௧. (12).} [0107] Observe that the only noise process in this rotated space is ^^^{ି^} ^^_௧. Multiplying by the _{masked observation vector ^} ^{0 ⋯ ^^} _^,^^ ^{⋯ 0} _{^ scales the noise variance back in scale} with the observed data ^^ and keeps only the contribution from a single pole. Thus, an estimate of the pole-specific noise variance is derived by applying the masking observation vector to the rotated noise covariance matrix ^^′ ൌ ^^^{ି^} ^^^{^} ^^^{ି^^∗}. Because ^^ is only non-zero in the first diagonal entry, the rotated noise covariance matrix is rewritten as: ^^′ ൌ ^^^ଶ ^^_⋅,^ ^{^}൫ ^^_⋅,^ ^{^}൯^∗ (13); where ^^_⋅,^ ^{^} is the first column of the inverse of ^^, and ^^^∗ is the conjugate transpose of ^^. 22  Q   B\125141.04438\85541799.2 [0108] Recall that the complex roots (or equivalently, the complex eigenvalues) appear as pairs in frequency except at DC. This means two poles in the AR process contribute to the spectral decomposition at any frequency greater than 0 Hz. This can be taken into account by summing variances across the poles with the same rotation frequency, which ignores the cross-correlation between them. Finally, the desired estimate of the oscillator state noise variance parameter ^^^ଶ is obtained: ^^{^ଶ ൌ ∑} _{^∈^^:^^ୀ^^^^} ^{^^ଶ} _ห ^{^^} _^,^ ^{^^} _^ ^ି ,_^ ^{^} _ห ^(14). [0109] in iOsc are: 1)

to recognize the complex eigenvalues under spectral decomposition as oscillators after rotation; and 3) to apply masking observation vectors in the rotated space to combine variances from relevant poles and zeros to initialize the oscillator parameter ^^^ଶ. [0110] Step 4: repeat for additional oscillators [0111] Steps 1-3 provide informed starting points of parameters ^ ^^, ^^, ^^^ଶ^ for the first oscillator when ^^ is the observed data. This first oscillator added is most often a slow oscillator (< 1 Hz) in human EEG recordings, since such slow oscillations tend to be large in amplitude and thus power. Recall that the search algorithm of iOsc will now attempt to add an additional oscillator on top of the first oscillator that has gone through EM learning. The outcomes of EM learning can be used for existing oscillators, but the parameters ^ ^^, ^^, ^^^ଶ^ need to be initialized for the new oscillator to be added. To achieve this, the AR fitting procedure is reapplied to the innovations, i.e., one-step prediction error (OSPE): ^^^_௧ ൌ ^^_௧ െ ^^ ^^ ^^_{௧ି^|௧ି^}. (15); [0112] ^^{^}^_௧ is a reasonable time series to repeat steps 1-3 because the current model only includes existing oscillators in the transition matrix ^^, therefore it only predicts oscillatory activity characterized by existing oscillators but not additional oscillations present in the observed data ^^ that are yet to be characterized. Therefore, the next oscillator to be added by iOsc can be identified by selecting the strongest pole/zero of an AR model fitted to this OSPE ^^^_௧. [0113] In contrast, if the residual error based on the filtered estimate ^^_௧|௧ ൌ ^^ ^^_௧|௧ were employed to discover additional oscillations, i.e., ^^^_௧ ൌ ^^_௧ െ ^^ ^^_௧|௧, those additional oscillations would not be visible since the Kalman filter adjusts the state estimates based on the most recent observations, absorbing the error from any model mis-specification into the state estimate. The same would be true for smoothers such as the fixed interval smoother. 23  Q   B\125141.04438\85541799.2 Thus the use of the OSPE for discovery of additional oscillators as described above is an essential feature of the disclosed approach. [0114] This algorithm design offers an efficient implementation because the same routines in steps 1-3 can be repeated identically after replacing ^^ with ^^^. But, there are two additional technical details that require some attention. First, when there are only a small number of oscillators present in a model, strong oscillations in the low frequency range may not be explained fully because an EM learned oscillator has skewed frequency and/or damping parameter to accommodate spectral content in higher frequency. As a result, the next strongest pole identified by an AR process fitted to the residual ^^^ might be still in the low frequency. This is problematic because at the next iteration there are two oscillators competing with each other in this frequency range. To circumvent this issue, a minimal frequency resolution (typically 1 Hz) is imposed when initializing new oscillators, meaning poles and zeros are only considered if they are more than 1 Hz away from existing oscillators. This is a mild constraint only applied during initialization, and if there are indeed oscillations closer together, EM learning is able to and will update oscillator frequencies to be less than 1 Hz apart. This frequency resolution is only used to prevent erroneous duplicated oscillators during search. [0115] Second, the choice of AR model order as 2 ^^ െ 1 in step 1 assumes at most ^^ oscillators. During subsequent search iterations after the first oscillator has been added, fitting AR(2 ^^ െ 1) will again identify ^^ additional frequencies. Since there are already oscillators present in the model, a fixed AR model order does not strictly respect the upper bound of numbers of oscillators. The AR model is kept order constant in order to ensure comparable initialization of oscillator parameters ^ ^^, ^^, ^^^ଶ^ across iterations. When there are no further oscillations remaining in the residual during later search iterations, fitting AR(2 ^^ െ 1) is similar to fitting a high order AR model to a white noise process. Diffuse poles with small radii will be identified, which are still compatible with iOsc and the use of log likelihoods for model selection, since only a single oscillator is added at a time. [0116] Decomposed oscillator modeling [0117] A non-iterative variant of the iOsc algorithm can be easily derived from the above routine with slight modifications. Since the spectral decomposition using complex eigenvalues in step 3 is not limited to just the largest pole/zero, the desired state noise variance parameter ^^^ଶ is initialized for any of the rotated oscillators after AR modeling and applying separate masking observation vectors. Similarly there is access to oscillator 24  Q   B\125141.04438\85541799.2 parameters ^ ^^, ^^^ for all poles and zeros in step 2. Therefore, a non-iterative approach to model ^^ can be constructed by initializing all of ^^ oscillators at once after an initial fit of an AR model of order 2 ^^ െ 1 to ^^. This removes the need for step 4 that examines OSPE because all oscillators to be added have been initialized appropriately. [0118] We refer to this approach as decomposed oscillator modeling (dOsc) because it relies on AR-based decomposition of a time series into complex eigenmodes, from which oscillators are initialized. The rest of processing in dOsc is identical to iOsc and will add all initialized oscillators one at a time in descending order of ^^^ଶ. The only difference is that each oscillator is no longer initialized after preceding oscillators have been added and learned with EM; initial parameters of all added oscillators are determined in one step with a single AR model fitting. [0119] As expected, iOsc and dOsc perform similarly in many cases given the overlap in the underlying computations. However, a key distinction lies in how the two approaches extract oscillations underlying recordings. iOsc is better geared to uncover masked oscillations that are only clearly visible after other oscillations have been explained away since it operates on innovation time series except the first oscillator. This is at the cost of potentially biasing subsequent added oscillators: this can be understood from the maximal likelihood nature of Kalman filtering and smoothing. When an insufficient number of oscillators than the true number of underlying oscillations are employed in a model, the oscillators at lower frequencies tend to get skewed away from their central frequencies in order to account for unexplained activity in other frequencies. The innovation spectrum may then contain more power in frequency ranges with existing oscillators and less power in frequency ranges without existing oscillators. Since oscillator parameters in iOsc are initialized from AR fitting of the innovation spectrum, this results in biased initial estimates and parameter values after EM learning. [0120] In comparison, dOsc does not suffer from biases due to oscillator fitting in previous iterations because oscillator parameters are initialized at once on the original data ^^. Thus when there are clear and separable oscillations in a recording, dOsc often recovers more precise parameter estimates and state estimation results. However, dOsc is not capable of distinguishing between oscillations that are very close by, especially if there is a weaker oscillation that is nested under a stronger oscillation that can only be modelled 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 out eigenmodes; dOsc is still much superior 25  Q   B\125141.04438\85541799.2 than AR modeling due to employing interpretable and efficient model structures, allowing more accurate representation of neural oscillations, and adding oscillators sequentially for model selection. Nevertheless, if no oscillator gets initialized at a value after AR fitting to ^^, dOsc does not have a mechanism to adapt and add an oscillator in frequency regions not well captured by previous oscillators (which requires examining the OSPE in iOsc), since it simply adds the set of oscillators initialized at the beginning in descending order of initialized values of ^^^ଶ. [0121] When it comes to studying real recordings of neural activity with a variable number of oscillations as well as uncertain oscillation properties, dOsc and iOsc nicely complement each other. Clearer understanding of neural oscillations is often achieved by applying both approaches to the recording and contrasting their modeling results. [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 ^^. In different recording modalities, there are different ways to estimate the observation noise. For instance, MEG can utilize empty room recording to directly measure the observation noise process. In EEG however, there is not an equivalent empty room recording available. Thus, an empirical estimate is derived of the assumed white observation noise after high-pass filtering above a selected frequency. This heuristic estimate is motivated by with the idea that EEG recordings above a certain frequency is predominated by noise processes. For example, in observed data ^^ sampled at Fs ൌ 100 ^^ ^^, the variance of ^^ high-pass filtered above 30 Hz, scaled by ^ହ ଷ to account for full frequency range of white noise, is used to obtain an initial estimate of

[0124] We employ an informed inverse gamma prior when updating the observation noise variance ^^^ଶ. The motivation is that oscillator state noise processes are additive with the observation noise in the formulation of equation (2). Thus, during Kalman filtering and fixed- interval smoothing, some of the observation noise gets explained away by the hidden oscillator processes, resulting in a biased estimate of ^^ during parameter learning. This can also be understood by noting that theoretical spectra of oscillators have non-zero values above 30 Hz, and the additive effect of multiple oscillators in this frequency range can account for a portion of a white observation noise process during state estimation, resulting in underestimated ^^. Thus, parameter updates of ^^ are constrained during iOsc and dOsc search by imposing a prior whose mode is placed at the empirically estimated ^^ described above using high-pass filtered data variance. 26  Q   B\125141.04438\85541799.2 [0125] To ensure comparable evaluation of log likelihoods for model selection, the same ^{prior must be used throughout iOsc iterations since the log likelihood expression involves ^^ଶ:} l_{og ^^ ൌ െ^ ଶ∑ ୀ^ log |2 ^^ ^^| െ ^ ଶ∑ ୀ^ ^ ^^௧ െ ^^ ^^௧ି^ ௧ ^ୃ ^^ି^^ ^^ ௧ି^ ௧ െ ^^ ^^௧ ^ (16)} where ^^ ൌ ^^ ^^_௧ ^{௧ି^} ^^^ୃ ^

with ^^_௧ ^{௧ି^} and ^^_௧ ^{௧ି^} being the conditional mean and variance of the ^^_௧ given ^ ^^_^,⋯ , ^^_௧ି^^. [0126] Model selection [0127] Now that the oscillators can be iteratively added to model the observed

proper initialization, each iteration of iOsc and dOsc simply performs EM learning for some number of iterations to update model parameters and estimate hidden oscillator states. Model selection is based on which model, i.e., which combination of oscillators, provides the best representation of ^^ using log likelihoods. Since more oscillators provide more degrees of freedom and can over-fit ^^, the knee point on the log likelihood curve spanning 1 to ^^ oscillators is selected to account for diminishing returns of additional model complexity (Satopaa et al.2011). The model and its learned oscillators at the knee point are the outputs of iOsc and dOsc as a set of oscillatory components underlying observed data ^^. 27  Q   B\125141.04438\85541799.2

Claims

Claims 1. A method for measuring intraoperative stimuli, the method comprising: using a patient monitor to receive electroencephalography (EEG) data from a patient while one or more stimuli are applied to the patient; using a computer processor, extracting one or more event-related potential (ERP) signals from the EEG data related to one or more discrete events; using the computer processor, determining the one or more stimuli from the extracted one or more ERP signals; and using a display, indicate an adjustment of one or more intraoperative parameters correlated to counteract the one or more stimuli in a 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 background signal for the EEG data, the one or more ERP signals, and the one or more discrete events.

3. The method of claim 2, wherein waveforms from the one or more ERPs signals are extracted using a state-space ERP (SS-ERP) model defined by the function: ^^{^^ ^^^^^ ௧,^ cos} ଶగ^ ^{ଶగ^^^^} ^^^ ^{^^^ ^ೞ െ sin ^ೞ ^^௧ି^,^ ^^௧,^ ,}

^ ^^_^, ^^ the

^ 1, and ^^^{^^^} ௧_,^, … , ^^_௧,^ are the state space oscillators, and where ^^ is the real component and ^^_௧,ଶ is the

component of phasor ^^_௧, ^^ is frequency, ^^_^ is the sampling frequency, ^^ is the damping factor, ^^ is system state noise, ^^^ଶ is the state noise variance, and ^^_௧ is measurement noise, ^^^ଶ is measurement noise variance.

4. The method of claim 3, wherein an evoked response is represented as a convolution between the one or more ERP waveforms, ^{^} ^^_^ ^{^^} ^, and one or more stimulus presentations, ^^_௧: ^ ^^{evoked^} , Q   B\125141.04438\85541799.2

with elements of ^^ following a random walk model with gaussian driving noise, ^^_^ with variance ^^^ଶ. 5. The method of claim 2, wherein the one or more ERP waveforms are extracted using a kernel expansion ERP (KE-ERP) model defined by the function: ଶగ^^{^^^ ଶగ^^^^ ^^^^^ cos െ sin ^^} ^^^ ^{^^^} ^ ^{௧,^} ^_{^^൩ ൌ ^^} ^{^^^} _൦ ^{^ೞ ^ೞ ௧ି^,^ ^^} ^_{^^ ^^^} ^{൪ ^} _^^^ ^{൩ ^ ^ ௧,^} ^_^^ ^{൩ , ^^ ൌ 1, 2,⋯ , ^^ ,} , fs

the ^^^{^^^} ௧_,^, … , ^^_௧,^ are the state space oscillators, and where

^^ is the real

^^_௧,ଶ is the imaginary component of phasor ^^_௧, ^^ is frequency, ^^_^ is the sampling frequency, ^^ is the damping factor, ^^ is system state noise, ^^^ଶ is the state noise variance, and ^^_௧ is measurement noise, ^^^ଶ is measurement noise variance.

5. The method of claim 4, wherein an evoked response is represented as a convolution between the one or more ERP waveforms, ^{^} ^^_^ ^{^^} ^, and one or more stimulus presentations, ^^_௧: ^ ^^^{^evoked^} ௧ ൌ ^ ^^ ∗ ^^^_௧ ൌ ^ℎ_^ ^^_௧ି^ , with its expansion or similar dictionary of overlapping kernels/bases.

6. The method of claim 1, wherein the one or more stimuli include noxious stimuli or painful stimuli.

7. The method of claim 1, wherein the one or more stimuli includes one or more electrical stimuli.

8. The method of claim 1, wherein the ERP signals are associated with the one or more stimuli. 29  Q   B\125141.04438\85541799.2

9. The method of claim 1, wherein adjusting the one or more intraoperative parameters includes adjusting an administration of one or more analgesic agents.

10. The method of claim 1, wherein the EEG data includes a signature of one or more analgesic agents.

11. A system for measuring intraoperative stimuli, the system comprising: one or more sensors configured to measure electroencephalogram (EEG) signals of a patient; and a processor configured to: receive EEG signals from a patient while one or more stimuli are applied to the patient; extract one or more event-related potential (ERP) signals from the EEG data related to one or more discrete events; determine the one or more stimuli from the extracted one or more ERP signals; determine an adjustment to one or more intraoperative parameters based on the one or more stimuli; and a display to communicate the adjustment to the one or more intraoperative parameters.

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

13. The system of claim 11, wherein waveforms from the one or more ERP signals are extracted using a state-space ERP (SS-ERP) model defined by the function: ^ ^cos ଶగ ^{^^^ ^^^ ^^} ^ െ ^{sin ଶగ^} ^^^ ^{^^^ ^ ^^ ^^ ,}  

where ^^^{^^^} ௧_,∗ ~ ^^^0, ^^_^ ^ଶ ^, ^^_௧~ ^^^0, ^^^ଶ^ are assumed to be independent, | ^^^{^^^}| ^ ^^_^, ^^_^ is the ^^^{^^^} ^ 1, and ^^^{^^^} ௧_,^, … , ^^^{^^^} ௧_,^ are the state space oscillators, and where and ^^_௧,ଶ is the component of phasor ^^_௧, ^^ is the

^^_^ frequency, ^^

damping factor, ^^ is system state noise, ^^^ଶ is the state noise variance, and ^^_௧ is measurement noise, ^^^ଶ is measurement noise variance.

14. The system of claim 13, wherein an evoked response is represented as a convolution between the one or more ERP waveforms, ^ ^^_^^^{^} ^, and one or more stimulus presentations, ^^_௧: ^ ^_^ ^{^evoked^} ௧ _ൌ ^{^} _{^^ ∗ ^^} ^{^} _{௧ ൌ ^ℎ^ ^^௧ି^} ^, with elements of ^^ following a driving noise, ^^_^ with variance ^^^ଶ.

15. The system of claim 11, wherein the one or more ERP waveforms are extracted using a kernel expansion ERP (KE-ERP) model defined by the function: ^^{^^ ^^^ ^^^^^ cos} ଶగ^ ^{െ sin ଶగ^ ^^} ^^^ ^{^^^^^ ,} Where

is ^ 1, and ^^^{^^^}, … , ^^_௧,^ are the state

and

and ^^_௧,ଶ is the imaginary component of phasor ^^_௧, ^^ is frequency, ^^_^ is the sampling frequency, ^^ is the damping factor, ^^ is system state noise, ^^^ଶ is the state noise variance, and ^^_௧ is measurement noise, ^^^ଶ is measurement noise variance.

16. The system of claim 15, where an evoked response is represented as a convolution between the one or more ERP waveforms, ^ ^^_^^^{^} ^, and one or more stimulus presentations, ^^_௧: ^ ^^{^^evoked^ ^ ^^ ^^^} ^^ , with its expansion or similar dictionary of

overlapping kernels/bases. 31  Q   B\125141.04438\85541799.2

17. The system of claim 11, wherein the one or more stimuli include noxious stimuli or painful stimuli.

18. The system of claim 11, wherein the one or more stimuli includes one or more electrical stimuli.

19. The system of claim 11, wherein the ERP signals are associated with the one or more stimuli.

20. The system of claim 11, wherein adjusting the one or more intraoperative parameters includes adjusting an administration of one or more analgesic agents.

21. The system of claim 11, wherein the EEG data includes a signature of one or more analgesic agents. 32  Q   B\125141.04438\85541799.2