Spontaneous pain dynamics characterized by stochasticity in awake human LFP with chronic pain

11 Chronic pain involves persistent fluctuations lasting seconds to minutes, yet there are limited studies on 12 spontaneous pain fluctuations utilizing high-temporal-resolution electrophysiological signals in humans. 13 This study addresses the gap, capturing data during awake deep brain stimulation (DBS) surgery in five 14 chronic pain patients. Patients continuously reported pain levels using the visual analog scale (VAS), and 15 local field potentials (LFP) from key pain-processing structures (ventral parietal medial of the thalamus, 16 VPM; subgenual cingulate cortex, SCC; periaqueductal gray, PVG) were recorded. Our novel AMI 17 analysis revealed that regular spike-like events in the theta/alpha band was associated with higher pain; 18 and regular events in the gamma band was associated with opioid effects. We demonstrate a novel 19 methodology that successfully characterizes spontaneous pain dynamics with human electrophysiological 20 signals, holding potential for advancing closed-loop DBS treatments for chronic pain. 21 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint Chronic pain is a pervasive condition, impacting approximately 20% of the global population1 22 and contributing to 15-20% of medical appointments2. Despite its prevalence, there is a limited 23 understanding of the spontaneous fluctuations in pain experienced by patients with chronic pain. These 24 fluctuations, a key symptom for these patients, involve continuous changes in pain intensity and sensation 25 lasting seconds to minutes, and are observed in conditions like irritable bowel syndrome (IBS)3, chronic 26 back pain4, postherpetic neuropathy5, chronic pelvic pain syndrome6, and chronic knee osteoarthritis7. 27 Earlier investigations on pain have aimed at elucidating nociceptive neurons8,9 and central pain 28 pathways primarily through animal studies10. Human neuroimaging studies have advanced this 29 understanding, revealing that chronic pain is associated with abnormalities in network connectivity11,12. 30 These studies have recognized the broader involvement of brain networks, extending beyond sensation to 31 include affect and cognition13,14, along with the salience network and default mode network15. 32 Among studies focused on spontaneous pain, many have examined the static resting-state 33 characteristics in chronic pain patients, revealing distinct network and connectivity strength with specific 34 structures 16,17. However, these endeavors merely offer a snapshot and fail to capture the ongoing pain 35 fluctuations reported by chronic pain patients in a dynamic manner18. Given the intricate link between 36 dynamic brain network processes and pain intensity19,20, it is important to understand the brain's dynamic 37 response to pain on a moment-to-moment basis. A small subset of studies sought to dynamically 38 characterize spontaneous pain fluctuations but were limited by neuroimaging techniques 5,21,22, which lack 39 temporal resolution necessary to dynamically capture neural activity. For that reason, examining the 40 dynamic mechanisms behind spontaneous pain would greatly improve our understanding of its emergence 41 and, ultimately, advance treatment options for chronic pain. 42 To achieve this objective, we captured local field potentials (LFP) from structures engaged in 43 pain processing from five awake patients with chronic pain undergoing awake deep brain stimulation 44 (DBS) implantation surgery. Specifically, we recorded from brain areas implicated in pain processing, 45 conceptualized as pain processing hubs - the ventral posteromedial nucleus (VPM) of the thalamus, the 46 periaqueductal gray region (PVG) 23,24, and the subgenual cingulate cortex (SCC) 25 (Figure 1A). During 47 awake DBS surgery, patients continuously rated their pain levels by moving a red circle along the visual 48 analog scale (VAS) line displayed on a screen, using a trackpad or buttons for about 10 minutes (Figure 49 1B). The continuously recorded pain levels were used to segment the data based on varying pain intensity. 50 These segments were analyzed by comparing two consecutive segments and examining the changes in 51 neural activity as pain intensity changed (Figure 1C). 52 To characterize pain dynamics, we present a novel computational method in which we quantify 53 the stochasticity in the LFP time series. We conceptualize irregularity to signify diverse (pain- and non-54 pain-related) information flow from nearby structures, while greater regularity to indicate a more focused 55 flow of pain-related information 11,26,27 (Figure 1D). We therefore hypothesize that irregularity in the LFP 56 would characterize a lower level of pain, and regularity to reflect a higher level of pain. Our method 57 involves extracting the time between two consecutive peaks of the theta/alpha bandpass filtered LFP, 58 termed inter-event-interval (IEI) (Figure 1E), which is essentially an instantaneous theta/alpha frequency, 59 as it is an inverse of the corresponding frequency cycle. Here, we quantify the level of regularity by 60 assessing the IEI stochasticity with auto-mutual information (AMI), a nonlinear analogue to auto-61 correlation. Here, a high AMI would suggest greater regularity, while a low AMI would suggest 62 irregularity. We demonstrate this novel methodology to effectively capture the dynamic changes in 63 response to shifts in pain experience. We further compare the performance of this characterization method 64 with conventional electrophysiology metrics to gain insights into the neural patterns and identify essential 65 features in characterizing spontaneous pain dynamics. 66 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 67 Figure 1. LFP recording, behavioral task, and computational method. 68 A. (Left) The DBS lead of a representative patient (P5) is targeting the left and right of the subgenual cingulate 69 cortex (SCC), reconstructed using 28. (Right) The schematic diagram of the DBS lead depicts a full-ring contact at 70 both the distal end (channel 1) and the proximal end (channel 8) of the cable. In between, there are contacts with 71 equally spaced directional contacts (channel 2-7). B. The patient continuously rated their pain level on a VAS from 72 0 to 10 by moving a red circle using either a trackpad or a button box. C. (Top) A representative patient (P5) 73 continuously rated the pain level (blue). Because this patient had low back pain, pain was induced by raising the leg 74 29 , and leg raises were recorded with accelerometers attached at the ankle (red). Instances of brief leg raises were 75 excluded from analysis. (Middle, Bottom) The spectrograms of the recorded LFP signal were normalized by z-score. 76 The pain VAS trajectory and their corresponding spectrograms for all other patients are shown in Figure S1. D. We 77 hypothesize that lower pain levels are characterized by signal irregularity, while higher pain levels exhibit regularity. 78 We conceptualize that irregularity reflects diverse information flow from both pain and non-pain-related sources, 79 while greater regularity indicates focused pain-related information. The schematic diagram illustrates this concept, 80 where black circles indicate pain-related information sources, white circles indicate non-pain-related sources, and 81 line thickness reflects the amount of information flow. E. We demonstrate a novel computational method to 82 characterize pain dynamics, which involves extracting the inter-event-intervals (IEI) - time between peaks in the 83 LFP time series - and assessing its stochastic patterns at different pain levels. 84 Abbreviation: VAS, visual analog scale; SCC, subgenual cingulate cortex; V pain rating number on a VAS (e.g., 85 V1 denotes pain rating 1 on a VAS) 86

87 Participants 88 Three patients with trigeminal neuropathic pain (participant ID - P2, P3, P4) and two patients 89 with chronic low back pain (participant ID - P1, P5), undergoing a bilateral or unilateral deep brain 90 stimulation (DBS) implantation surgery, participated in this study. The target(s) of the DBS lead included 91 the ventral-parietal-medial (VPM) of the thalamus, periventricular gray (PVG), and/or subgenual 92 cingulate cortex (SCC) (Figure 1A). All subjects provided written informed consent to participate in a 93 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint research study approved by the institutional review board (IRB) at the University of California, Los 94 Angeles. Details of the patient demographics along with their relevant clinical information can be found 95 in Table 1. 96 Table 1. Participant demographic, relevant clinical information, and change in pain events. 97 Participant ID P1 P2 P3 P4 P5 Diagnosis CLBP TNP TNP TNP CLBP Sex M F M M M Age 68 64 67 67 50 Weight (kg) 79 97 96 103 135 Recording Site1 B SCC, OFC, PFC R VPM L VPM, L PVG L VPM, OFC, PFC B SCC, OFC, PFC Bipolar Reference Channels2 1 - μ (2,3,4) μ (2,3,4)- μ (5,6,7) VPM:1- μ (2,3,4) PVG:1-μ (2,3) μ (2,3,4)- μ (5,6,7) μ (2,3,4)- μ (5,6,7) Sampling Rate (Hz) 5500 5500 5500 4800 5500 Condition3 Spontaneous (N) 1 1 0 4 7 Δ VAS4 1 → 0 9 → 3 - 1→ 0, 8→ 0, 8→ 4, 6→ 4 2→ 1, 3→ 2, 4→ 3, 6→ 5, 5→ 4, 7→ 6, 8→ 7 Fentanyl (N) 0 2 1 0 1 Bolus (mcg) n/a 200 (alfentanil5), 50 12.5 n/a 25 Propofol (N) 0 1 1 1 0 Bolus (mg) (infusion rate ug/kg/min) n/a 50 mg (125) 30mg (50) 30mg (25) n/a 1Recording sites are listed, but only areas that are underlined are analyzed for the purpose of the study. 2Channels 98 were bipolar referenced to capture the local activity of the structure. μ denotes average of the channels listed in 99 parenthesis. Channel numbers correspond to the location along the DBS lead schematically depicted in Figure 1A. 100 3Condition shows the number of comparisons (N) in sequential events, where we assess lower pain states by 101 comparing them to the previous or subsequent higher pain state. For spontaneous lower pain events (Spontaneous), 102 changes in VAS (visual analog scale) pain ratings are listed. For pain events via fentanyl, the bolus for each event is 103 listed. Events via propofol induction show the bolus and infusion rate. 4Change in VAS lists events where the rating 104 either increased or decreased in time. We assessed lower pain events by comparing them against the immediately 105 preceding or subsequently higher pain instances. 5Initially alfentanil (faster-acting but less potent than fentanyl) was 106 administered, then followed by fentanyl. 107 Abbreviation. CLBP, chronic low back pain; TNP, trigeminal neuropathic pain; B SCC, bilateral subgenual 108 cingulate cortex; OFC, orbitofrontal cortex; PFC, prefrontal cortex; R VPM, right ventral posteromedial nucleus of 109 the thalamus; L VPM, left ventral posteromedial nucleus of the thalamus; L PVG, left periventricular gray; Δ VAS, 110 change in visual analog scale pain rating. 111 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint Behavioral Task 112 The patient used a visual analog scale (VAS) displayed on a monitor, along with a trackpad (P1-113 P4) or button box (P5), to continuously indicate their pain level. Each patient moved a red dot along the 114 VAS scale, which ranged from 0 to 10. The trackpad enabled continuous ratings between 0 and 10, while 115 the button box facilitated discrete integer rating. Patients continuously assessed their pain levels for a 116 duration of 5 to 15 minutes. Afterwards, they were given either fentanyl or propofol based on the clinical 117 protocol reflecting their specific needs. The recording time varied depending on how much the patient's 118 pain level fluctuated. If a patient had no pain (e.g., patient P3), the pain rating ended after 5 minutes. 119 However, if a patient's pain rating fluctuated significantly (e.g., patient P4, shown in Figure S1-D), they 120 continued rating their pain for up to 15 minutes. For the two chronic back pain patients (P1, P5), we also 121 utilized intermittent leg raising to induce low back pain similar to the endogenous pain experienced by 122 these patients 29. To record the time and extent of the leg raise, we strapped 2 opal inertial measurement 123 unit (IMU) movement sensors (APDM, USA) to the patients’ ankles. The pain rating trajectory of each 124 patient, along with the timing of leg raises for the two chronic low back pain patients are shown in Figure 125 1A and S1. 126 Surgery and data acquisition 127 During an awake DBS surgery, intraoperative recordings were conducted for research purposes. 128 For all patients, local field potentials (LFP) at the DBS target (P1: bilateral SCC; P2: right VPM; P3: left 129 VPM and PVG; P4: left VPM; P5: bilateral SCC) were recorded from the DBS lead in its final implant 130 position. Patients P1 and P5 had DBS leads manufactured by Abbott (Model 6173ANS, 1.27 mm lead 131 body diameter, inter-contact distance 1.5mm; Abbott Laboratories, Chicago, IL, USA), while the other 132 patients had DBS leads manufactured by Medtronic (Model B3301542 or B3301542M, 1.36 mm lead 133 body diameter, inter-contact distance 2.2mm; Medtronic, Inc., Minneapolis, MN, USA). All DBS leads 134 have a total of 8 electrode contacts and are designed as those shown in Figure 1A. Based on the channel 135 labels in the diagram of Figure 1A, channels 1 and 8 are the most distal and proximal ones respectively, 136 and have full contact. Between them, there are equally spaced 3 directional contacts each at 2 levels, 137 ranging from channels 2 to 4 and channels 5 to 7. LFPs were recorded using these amplifiers: 138 NeuroOmega (AlphaOmega Engineering, Nazareth Israel) for patients P1, P2, P3, P5, Matlab/Simulink 139 software connected to an amplifier (g.Tec, g.USBamp 2.0) for patient P4. The signals obtained from the 140 former amplifier were sampled at 5500Hz, and the latter at 4800Hz. The ground and reference signals 141 were registered using alligator clip electrodes attached to the cannula that was used for the DBS 142 implantation. 143 For patients with chronic back pain (P1, P5), their legs were intermittently raised to simulate their 144 low back pain. We used two Opal IMU movement sensors (APDM, USA) strapped on their ankles to 145 capture the kinematics, specifically the angular velocity, and measure the timing and extent of the leg 146 raise. To synchronize between the electrophysiological and kinematic signals, we utilized an external 147 synchronization equipment (“sync box” provided by APDM, USA). This equipment sent a digital output 148 trigger to the electrophysiological signal amplifier, precisely indicating the start and end times of the IMU 149 sensor recording. 150 For patients who received either fentanyl or propofol, we manually recorded the timing of the 151 medication administration. Due to potential human error, the precise timing of the administration may 152 have a deviation less than 5 seconds. However, since both medications typically take around less than a 153 minute to take effect, such minor imprecision is not expected to significantly impact the overall results 154 presented in this study. 155 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint Analysis 156 Preprocessing 157 LFPs recorded from each of the 8 channels of the DBS lead were applied with a 60 Hz notch 158 filter (band-stop IIR, filter order 80) and then visually examined for artifacts and excess noise. To reduce 159 noise from the 8 recorded channels, we performed the following data processing steps. First, we averaged 160 the signals from the 2nd level (channels 2-4) and 3rd level (channels 5-7) of the DBS lead (Figure 1A). 161 Then, we applied bipolar referencing between adjacent levels as follows: Channel 1 was bipolar 162 referenced against the average of channels 2-4; the average of channels 2-4 was bipolar referenced against 163 the average of channels 5-7; and the average of channels 5-7 was bipolar referenced against channel 8. 164 Given that the distal part of the DBS lead (channel 1) is typically implanted closest to the targeted 165 structure, we prioritized examining the bipolar referenced signal of channel 1 against the average of 166 channels 2-4. However, in certain cases where channel 1 exhibited excess noise (i.e., lacking 1/f structure 167 in the power spectrum), we resorted to examining the bipolar referenced signal of the average of channels 168 2-4 against the average of channels 5-7. Additionally, for patient P3, the PVG LFPs recorded showed 169 excess noise in channels 2-5. As a result, we opted to bipolar reference the signal between channel 1 and 170 the average of channels 6-7. The exact channel numbers used for the bipolar referencing are provided in 171 Table 1. Once the channels were identified for further analysis, we conducted visual inspections to 172 exclude segments with electrical artifacts. These artifacts were characterized by a sudden and distinct 173 change in amplitude (>500uV), lasting for more than 1 second. We removed such segments, and the total 174 length of removed segments did not exceed more than 5 seconds in duration for any of the patients. 175 To segment the data by epochs based on the patient's pain levels, we considered both the pain 176 rating and the verbal reports conveyed during the recording. As an example, for patient P1, who 177 expressed feeling minimal pain throughout the recording, we excluded data associated with leg raises to 178 avoid any interference from the effects of leg raise, and only selected times when the patient verbally 179 expressed change in pain. We also considered the acclimation time for patients who used the trackpad or a 180 button box to indicate their pain rating. Patients using the trackpad took around 1-2 minutes to get 181 acclimated, while patients using the button box took less than 1 minute. For all patients, we excluded the 182 initial acclimation period from the analysis. Furthermore, we addressed the issue of excessive fluctuations 183 in VAS ratings registered through the trackpad (due to patients' unintentional motions) by applying a 184 moving average method with a smoothing factor 0.25 ("smoothdata.m" function in Matlab). This allowed 185 us to visually identify time segments with different pain levels. Lastly, for patient P5, whose leg was 186 raised intermittently to induce pain, we excluded times when the leg was in motion. Specifically, we 187 selected time segments when there were no leg movements (i.e., angular speed of the angle was less than 188 0.01 degrees/second) for at least 10 seconds. After carefully considering these criteria, which are 1) 189 ensuring that the pain ratings aligned with the patient's verbal reports, 2) accounting for their familiarity 190 with the trackpad or button box, 3) identifying stable VAS ratings, and 4) excluding data with leg 191 movements, we proceeded to segment the data into epochs of different pain levels obtained from this 192 process. To capture the impact of pain medication, we created an epoch by extracting data starting from 1 193 minute after the administration of the medication and spanning approximately 45 to 60 seconds. The 194 segmentation is visually depicted by green boxes in Figure 1C for patient P4 and in Figure S1 for all other 195 patients. At the end, these resulted in 33 epochs across all patients (median epoch duration 42 seconds; 196 IQR = 17, 60 seconds). 197 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 198 199 Figure 2. Auto mutual information (AMI) computation pipeline 200 A. LFP from each recorded structure was divided into epochs. Within each epoch, two types of windows were 201 chosen for analysis. The first window had a duration of 5 seconds, and it was paired with another window delayed 202 by 175ms (shown in green). These windows were used for theta/alpha bandwidth analysis (4-12hz). Another type of 203 window had a duration of 1 second and was paired with another window delayed by 30ms (shown in orange). These 204 windows were used for gamma bandwidth analysis (30-80hz). For both 5-second and 1-second paired windows, 205 these were slid forward with 50% overlap across the entire epoch time frame. B. To analyze within the theta/alpha 206 bandwidth, the 5-second time windows were bandpass filtered at 4-12hz. The local peaks (denoted in red) within 207 these windows were marked to compute the inter-event intervals (IEI), which is defined as the time between peaks. 208 C. IEI time series were created for the paired time windows. D. Within the 5-second time window, probability 209 distributions of the inter-event intervals (IEIs) were computed for the paired window, along with the joint 210 probability distribution of paired-IEIs. This process yielded a single AMI at each slid time point, and the AMI 211 values across the entire epoch were collected. The mean of the AMI (shown in red line) was extracted for later 212 comparison between epochs. E. To analyze within the gamma bandwidth, 1-second time windows were bandpass 213 filtered at 30-80hz, and local peaks were extracted (shown in red dots). F. Similar to C., time series data for IEI were 214 generated for the corresponding windows. G. Probability distributions and a joint probability distribution of the IEIs 215 were computed and used to calculate the AMI at each time point. The AMIs from all time points were collected, and 216 the mean of the AMIs (shown as a red line) was extracted for later comparisons. 217 Auto mutual information (AMI) based on inter-event-interval (IEI) 218 After dividing the dataset into epochs, we proceeded to calculate the auto mutual information 219 (AMI) for each epoch and observed the changes in AMI when pain was lower. AMI is considered a 220 nonlinear analog of the auto-correlation function, and employs information theory 30to quantify the 221 predictability of future points in a time series based on past points. The relationship between past and 222 predicted future time points is reflected in the diminishing AMI as the time delay increases, providing 223 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint insights into the system's complexity 31. Notably, research has demonstrated that higher AMI values 224 (indicative of reduced complexity) are linked to with various conditions such as schizophrenia32, 225 idiopathic dilated cardiomyopathy 33, Alzheimer’s disease 34,35, multiple organ dysfunction syndrome and 226 heart failure 36. In this context, we computed the epoch-specific AMI to characterize the dynamics of pain, 227 and to gain insight into the obtained results in relation to complexity. 228 The formula to compute AMI is identical to mutual information (MI), a measure that quantifies 229 the probabilistic dependency between signals from two different sources. In the context of AMI, the 230 dependency is computed within a single source, where one signal is temporally shifted by a specified time 231 delay in relation to the other signal. To elaborate, let’s consider X to represent the IEI series from the first 232 time window (i.e., future points) and Y to represent the IEI series from the time window with a time delay 233 (i.e., past points). AMI, denoted as /g1835 /g3025/g3026 in the formula below, quantifies the extent to which the 234 uncertainty of X can be reduced by having knowledge of Y. This concept is defined in equation (1), 235 wherein /g1834 represents an entropy function that quantifies the degree of uncertainty. 236 /g1835 /g3025/g3026 /g3404 /g1834 /g4666 /g1850 /g4667 /g3398/g1834 /g4666 /g1850 | /g1851 /g4667 ( 1 ) 237 Here, /g1834 /g4666 /g1850 /g4667 is the entropy of X, which signifies the inherent uncertainty within X, as defined in equation 238 (2). The term /g1834 /g4666 /g1850 | /g1851 /g4667 represents the conditional entropy of X given Y, depicting the remaining 239 uncertainty in X when Y is known, as defined in equations (3-5). 240 /g1834 /g4666 /g1850 /g4667 /g3404 /g3398 ∑ /g1842/g4666/g1876 /g3036 /g4667/g1864/g1867/g1859 /g2870 /g1842/g4666/g1876 /g3036 /g4667/g3051 /g3284 ( 2 ) 241 /g1834 /g4666 /g1850 | /g1851 /g4667 /g3404 ∑ /g1842 /g3026 /g3435/g1877 /g3037 /g3439/g1834/g3435/g1850/g3627/g1851 /g3404 /g1877 /g3037 /g3439/g3052 /g3285 ( 3 ) 242 /g3404 /g3398 ∑ /g1842 /g3025/g3026 /g3435/g1876 /g3036 ,/g1877 /g3037 /g3439/g1864/g1867/g1859 /g2870 /g3017 /g3273/g3274/g3435/g3051 /g3284,/g3052 /g3285/g3439 /g3017 /g3274/g3435/g3052 /g3285/g3439/g3051 /g3284/g3052 /g3285 (4) 243 /g3404/g1834 /g4666 /g1850 , /g1851 /g4667/g3398/g1834 /g4666 /g1851 /g4667 ( 5 ) 244 Based on equations (2) and (5), /g1835 /g3025/g3026 can be simplified as equation (6), which was used to calculate the 245 AMI values in this study. 246 /g1835 /g3025/g3026 /g3404/g1834 /g4666 /g1850 /g4667 /g3397/g1834 /g4666 /g1851 /g4667 /g3398/g1834 /g4666 /g1850, /g1851 /g4667 ( 6 ) 247 To do this, first, the LFP was initially down-sampled to 400Hz. Subsequently, we computed two 248 distinct sets of AMIs, each based on different bandwidths. Note, we chose to initially down-sample the 249 data to 400hz, because the minimum IEI was found to be 2.5ms based on the raw data that was filtered at 250 30-80hz (gamma), and the duration of a single frame within a 400hz data equals 2.5ms. 251 The first AMI was computed using data that underwent bandpass filtering within the theta and 252 alpha frequency range (4-12Hz; FIR, 200th order zero-phase, transition width 0.2). These filtered time 253 series were then partitioned into 5-second intervals, which were paired with a 5-second window delayed 254 by 175ms (Figure 2A). For each of these time windows, local peaks were identified (marked as red dots in 255 Figure 2B). The duration between consecutive local peaks (derived with “findpeaks.m” function in 256 Matlab; minimum peak prominence 0, threshold 0, minimum peak distance 0), referred to as inter-event-257 intervals (IEIs), was extracted for every time window. Given the sequence of IEIs, these were represented 258 as a time series aligned with the 400hz frame size (Figure 2C). For each paired time windows, probability 259 distributions were constructed for the IEIs of the preceding time window (referred to as P(IEI) in Figure 260 2D or P(X) in the equations above) and for the IEIs of the subsequent time window (referred to as P(IEI’) 261 in Figure 2D or P(Y) in the equations above). Additionally, a joint probability distribution (P(IEI,IEI’) in 262 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint Figure 2D or P(X,Y) in the equations above) was constructed based on the IEI and IEI’ occurrences. 263 These distributions enabled us to compute the components /g1834 /g4666 /g1850 /g4667 ,/g1834 /g4666 /g1851 /g4667 ,/g1834 /g4666 /g1850, /g1851 /g4667 , which ultimately 264 produced a single AMI value. Paired time windows were shifted forward with 50% overlap across the 265 entire epoch, resulting in a series of AMI values (rightmost graph in Figure 2D). The mean AMI was 266 selected for comparison between epochs. Note, the variables /g1861 and /g1862 in the equations (2-4) correspond to 267 the bins encompassing the range of IEIs, spanning from 10 to 130 frames with increments of 10 frames 268 (i.e., 25ms to 325ms in 25ms increments). We chose the time window length of 5-second to gather around 269 50 IEI datapoints, enabling a meaningful frequency distribution, while capturing the dynamical changes 270 in AMI over the epoch duration (about 45 seconds). We chose the shift size of 70 frames (175ms), as it 271 effectively represents the IEI’s mean duration, and optimally differentiates the epochs with differing pain 272 intensity states (See Figure S2 for different lengths of shift size). 273 The second AMI was computed using data that was bandpass filtered within the gamma 274 frequency range (30-80Hz; FIR, 200th order zero-phase, transition width 0.2). These filtered time series 275 were then partitioned into 1-second intervals, which were paired with a 1-second windows delayed by 276 30ms. Then, local peaks were identified (marked as red dots in Figure 2E) to compute the IEIs for every 277 time window. Given the sequence of IEIs, these were represented as a time series (Figure 2F), and the 278 paired time window’s IEI time series were used to construct the necessary probability distributions to 279 compute the AMI at a given time point (Figure 2G). Paired time windows, overlapping by 50%, were 280 shifted across the entire epoch. Mean AMI within each epoch was then compared to another (rightmost 281 graph in Figure 2G). Note, the distribution bins (variables /g1861 and /g1862 in the equations (2-4)) spanned from 5 282 to 14 frames with increments of 1 frame (i.e., 12.5ms to 35ms in 2.5ms increments). We chose the time 283 window length of 1-second to gather around 50 IEI datapoints, enabling a meaningful frequency 284 distribution, while capturing the dynamical changes in AMI over the epoch duration. We chose the shift 285 size of 12 frames (30ms), as it effectively represents the IEI’s mean duration, and optimally differentiates 286 the epochs with differing pain intensity states (See Figure S2 for different shift size). 287 We highlight that while the AMI metric has found utility in prior neural studies (e.g.,33–37), 288 variations exist in defining the constituent variables (i.e., X and Y in equations 1-6) for constructing 289 probability distributions. Indeed, converting a continuous signal like electrophysiological data into an 290 appropriate discrete probability space is a challenging task. To that end, we introduce an innovative 291 approach rooted in kinematics studies 38,39 to define constituent variables. Diverging from conventional 292 AMI computations involving amplitude binning, our method concentrates solely on the timing 293 information, disregarding amplitude details. By reducing the raw time series to a binary spike train, this 294 approach effectively simplifies the computation while presenting a versatile solution suited to diverse 295 signal modalities. Notably, prior studies (e.g., electrocorticography40 , EEG41, voice42, magnetometer43) 296 have successfully harnessed the inter-event interval (IEI) information to characterize a range of 297 biophysical signals, showing its adaptability across different modalities. 298 Statistical comparisons of AMI between low and high pain levels 299 To assess the AMI changes when pain is lower, we analyzed the data in two ways, by: 1) 300 combining all lower pain instances across all brain structures (VPM, PVG, SCC) and 2) separating the 301 lower pain instances by structure. When these data were separated by structure, we are left with a small 302 sample size. For that reason, when we combined the structures, we performed the Wilcoxen signed rank 303 test for zero median, and when we separated the data by structures, we performed a series of permutation 304 tests to see whether there were significant changes in the signal’s AMI. 305 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint Wilcoxon signed rank test for zero median. Using each instance of lower pain, we conducted a 306 Wilcoxon signed rank test for zero median (non-parametric analog of a one-sample t-test) to determine if 307 the difference between matched samples comes from a distribution with a median of zero. To do this, we 308 collected the mean AMI values from each lower pain epoch and its adjacent (subsequent or prior) higher 309 pain epoch (from here on, we will refer these as the paired epochs). Then, we subtracted the lower pain 310 epoch’s mean AMI value from the higher pain epoch’s mean AMI value for each paired epochs (low pain 311 AMI – high pain AMI). This non-parametric test identifies statistical significance when the paired epochs 312 exhibit an AMI difference departing from zero. Note, based on a power analysis considering our available 313 data points, achieving significant results requires a total sample size of 16 for a two-sided one-sample t-314 test with a power of 0.8 and an alpha of 0.05. While the test applied to AMI comparisons during 315 spontaneous lower pain is adequate for reaching statistical significance (N=21), the sample size is 316 insufficient for AMI comparisons during lower pain induced by fentanyl (N=8) and propofol (N=4). We 317 complement these findings with permutation tests performed separately for each structure below. 318 Permutation test. To address the challenge posed by a relatively small sample size when 319 analyzed per structure, we performed a two-sided permutation test involving 1000 permutations. The 320 permutation test aimed to assess whether the mean difference between the lower and higher pain states 321 within each paired epochs significantly differed from zero. For the theta/alpha range AMI comparisons, 322 we randomly selected 3 AMI values without replacement from each epoch. Note, within a single epoch, 323 there were 115 AMI values on average (ranging from 19 to 409). We then shuffled the AMI values 324 between the epoch pairs (higher and lower), and computed the difference between mean AMI from the 325 shuffled-lower and shuffled-higher epochs. We repeated this 1000 times, to generate a null distribution of 326 5000 datapoints (5 paired epoch x 1000 repetition) for the spontaneous lower pain in VPM, 4000 (4 327 paired epoch x 1000 repetitions) for fentanyl-induced in VPM, 3000 (3 paired epoch x 1000 repetitions) 328 for propofol-induced in VPM, 16000 (16 paired epoch x 1000 repetition) for spontaneous lower pain in 329 SCC, 4000 (4 paired epoch x 1000 repetition) for fentanyl-induced in SCC, and 1000 (1 paired epoch x 330 1000 repetition) for propofol-induced in PVG. The observed mean difference between the paired epochs 331 were then compared to this null distribution, allowing us to determine statistical significance. For the 332 gamma range AMI comparisons, we repeated the above permutation test methods, but by randomly 333 selecting 25 AMI values without replacement from each epoch. Within a single epoch for gamma band-334 passed signal, there were 125 AMI values on average (ranging from 27 to 125). 335 Kruskal-Wallis nonparametric one-way analysis of variance (ANOVA). To confirm that the 336 AMI values do not differ between different types of lower pain, we performed the Kruskal-Wallis test on 337 the AMI values across spontaneous, fentanyl-induced, and propofol-induced conditions at both higher and 338 lower pain states. Given the skewed non-normal distribution of the AMI, and unequal size and variance 339 between conditions, we chose this non-parametric test that allowed to address these data features. 340 Comparison metrics 341 To illustrate the utility of AMI and assess its critical features in characterizing pain dynamics, we 342 also applied established metrics commonly used to analyze neural activity. These metrics capture 343 different signal attributes such as non-stationarity, non-linearity, and dynamic past dependency. By 344 comparing the efficacy of pain characterization between AMI and these metrics, we identified the critical 345 attributes to capture the neurophysiological patterns of pain. 346 Power Spectrum. Using the preprocessed LFP, we conducted power spectrum analysis utilizing 347 the BOSC algorithm44. This examination covered individual frequencies ranging from 1 to 100Hz with a 348 step width of 1Hz and employed 6th order Morlet wavelets. To standardize the power time series, we 349 applied z-scoring to each frequency by normalizing across the entire recording. We then extracted the z-350 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint scores from epochs and compared these across different epochs. By analyzing the changes in z-scores 351 derived from the Morlet wavelet power spectrum, we aimed to determine whether a metric capable of 352 mitigating the stationarity assumption proves effective in characterizing pain dynamics. 353 As an alternative to the BOSC method, we also computed the power spectrum using the Hilbert-354 Huang Transform (HHT)45 with frequency limits from 1 to 100hz. Intrinsic mode functions (IMF), that 355 are necessary to compute the HHT power spectrum, were set with default parameters specified in the 356 Matlab functions emd.m and hht.m. After obtaining the Hilbert spectrum for the entire recording, we 357 applied z-scoring to each frequency for normalization and compared the z-scores corresponding to 358 different epochs. Using the HHT method, we aimed to determine whether a metric that relaxes both 359 linearity and stationarity constraints is effective in characterizing pain dynamics. 360 To statistically compare across frequency bands ranging from 1 to 100 Hz, we applied the 361 aforementioned permutation test at each individual frequency level. Additionally, we accounted for 362 multiple comparisons by employing the Bonferroni correction. This involved dividing the significance 363 level (alpha 0.05) by the total number of frequency bins (100) to maintain an adjusted threshold for 364 statistical significance. This enabled us to assess whether a non-dynamic metric with relaxed assumptions 365 of stationarity and/or linearity could effectively characterize pain dynamics. 366 Auto-Correlation. Using the bandpass-filtered and epoch-segmented LFP, we computed auto-367 correlation for each paired epochs, adhering to the identical time window lengths and sliding percentage 368 as specified in the AMI approach. We performed the Wilcoxen signed rank test based on the mean auto-369 correlation values from the paired epochs, as was done with AMI analyses. Using this metric, we aimed to 370 evaluate whether a linear and dynamic metric that captures past dependency such as auto-correlation 371 exhibited a positive correlation with the pain VAS. 372 Entropy. Using the spike trains generated for AMI calculation, we extracted the inter-event 373 interval (IEI) data points within each epoch. We then computed entropy using the same binning 374 parameters as in the AMI method. We computed the entropy from paired epochs and compared the 375 change in entropy using the Wilcoxen signed rank test. Notably, this approach captured the static 376 "uncertainty" observed throughout the epoch and did not capture dynamic past dependency. The objective 377 of this comparison was to determine whether a static information-theory metric, such as entropy, captures 378 pain more effectively than the dynamic AMI metric. 379 Statistical correlation with pain VAS 380 To determine if there's a correlation between the pain VAS and the aforementioned metrics, we 381 calculated the Kendall rank correlation coefficient (referred to as Kendell’s τ coefficient) between pain 382 VAS and each metric. Kendall’s coefficient is a suitable choice since it's a non-parametric test commonly 383 used to gauge ordinal associations. 384 In addition, we assessed whether the direction of change in pain (increasing vs. decreasing in time) 385 makes a difference in the AMI. To do this, in contrast to how we assessed lower pain by comparing any 386 two sequential epochs, we separated by those that are increasing pain in time versus those that are 387 decreasing in time. After separating the AMI values by increasing pain and decreasing pain, we 388 conducted a Kendall’s correlation between the corresponding pain ratings and the AMI for each set. 389

390 391 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 392 393 Figure 3. Results of Theta/Alpha Frequency Band (A-D) and Gamma Band (E-H) Dynamics. 394 Theta/Alpha Band Dynamics: A. AMI values did not show differences between S, F, and P conditions under high 395 (left) and low (middle) pain states, based on Kruskal-Wallis non-parametric ANOVA test. When changes in AMI 396 are examined for each lower pain instances (i.e., difference in paired epochs) using the Wilcoxon signed-rank test 397 for a zero median (right), AMI is reduced under S and F conditions, but not under P condition. B. Mean AMI values 398 are plotted for each paired epochs - low and high pain states - per patient (differentiated by marker color) and by 399 conditions (differentiated by marker shape). Each data point corresponds to a single paired epoch. C. The change in 400 AMI from each condition is assessed using permutation tests, showing significant reduction for each. The red line 401 denotes the observed AMI change, while the distribution represents the shuffled null distribution. D. When the 402

of AMI change are assessed separately by structures, the decrease in AMI is consistent for S and F conditions 403 across all structures. Gamma Band Dynamics: E. AMI values do not show differences between S, F, and P 404 conditions under high (left) and low (middle) pain states, based on the Kruskal-Wallis non-parametric ANOA test. 405 However, when changes are examined per paired epochs using the Wilcoxon signed-rank test for a zero median 406 (right), AMI values increase only under F condition. F. Similar to B, single data points are plotted for each paired 407 epochs differentiated by marker color (patient) and shape (condition), where F conditions are denoted by a larger 408 diamond shape. G. Similar to C, permutation test results are performed for each condition, revealing an increase in 409 AMI under F condition, and a reduction in AMI under S and P conditions. H. Similar to D, the Wilcoxon signed-410 rank test is performed separately by structure, revealing an increase in AMI values under F conditions for all 411 structures. Abbreviation: S, spontaneous lower pain condition; F, fentanyl based lower pain condition; P propofol 412 administered condition; S-Sc, spontaneous lower pain condition in SCC; F-Sc, fentanyl-based lower pain condition 413 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint in SCC; P-P, propofol-induced condition in PVG; L-H, Low-High; Δ AMI, change in AMI from lower pain (i.e., 414 AMI from low pain – AMI from high pain); **, significance p<0.01; *, significance p<0.05; ns, non-significant. 415 Lower pain is characterized by irregularity (reduced AMI) in the theta/alpha frequency 416 band 417 The AMI derived from the LFP within the theta/alpha frequency bands did not differ across the 418 three conditions – spontaneous vs. fentanyl-induced vs. propofol-induced (during high pain, χ 2(2,30)=3.1, 419 p=0.21; during low pain, χ 2(2,30)=3.6, p=0.16). However, when we compared between the paired epochs, 420 AMI was significantly less than zero when pain was lower (Wilcoxon signed-rank test for zero median). 421 This AMI reduction was observed both when pain was lower spontaneously (n=21; p=0.04) and when 422 fentanyl was administered (n=8; p=0.02) but not under propofol administration (n=4; ns) (Figure 3A). 423 The individual data points illustrating the mean AMI for the paired epochs separated by recording site per 424 patient and by conditions are shown in Figure 3B. A majority of datapoints were below the unity line, 425 indicating that AMI is generally lower when pain is lower. 426 To address the limitation posed by a small sample size, permutation testing was also performed 427 on these AMI data points. These tests revealed a consistent AMI reduction when pain was lower across all 428 three conditions (p<0.01) (Figure 3C). Given that these data points encompass the VPM (n=5 for S; n=4 429 for F; n=3 for P), SCC (n=16 for S; n=4 for F; n=1 for P), and PVG (n=1 for P) structures, we also 430 segregated these structures, and confirmed that the directional shift in AMI remains uniform across 431 structures, with the exception of the case involving propofol administration in the VPM (Figure 3D) . 432 It is worth noting that we also conducted tests for the temporal shift between time windows at 433 various lengths (ranging from 25ms to 275ms). The results indicate that a temporal shift of 175ms (as 434 shown in the main results) provides the optimal distinction between high and low pain states (Figure S2 435 A). Furthermore, we explored different frequency bands, including theta only, alpha only, alpha/low beta, 436 and theta/alpha/low beta, and found that the theta/alpha bands offer the most effective characterization of 437 pain dynamics (Figure S2 C). 438 Effect of opioid pain medication is associated with regularity (increased AMI) in the 439 gamma frequency band 440 To examine the effect of opioid pain medications on LFPs within the VPM, PVG and SCC, we 441 performed a similar analysis as the above, but limited the bandwidth to the gamma frequency range. Here, 442 the AMIs did not exhibit any differences between the three conditions, whether it was during high pain 443 (χ 2(2,30)=2.19, p=0.3) or low pain (χ 2(2,30)=0.76, p=0.7) states. However, we observed an increase in 444 AMI upon fentanyl-induced lower pain, (Wilcoxon signed-rank test for a zero median) (n=8; p=0.02), and 445 no change in spontaneous (n=21; p=0.24) and propofol-induced (n=4; p=0.38) conditions (see Figure 3E). 446 Figure 3F presents individual data points depicting the mean AMI for each paired epochs denoted by 447 recording site, per patient, and condition. 448 We also conducted permutation testing on the change in AMI datapoints per paired epochs. These 449 tests consistently revealed an increase in AMI under fentanyl-induced lower pain (p<0.01), and AMI 450 reduction for the other two conditions (p<0.01), as illustrated in Figure 3G. Similar to the above analysis, 451 we also segregated these by recorded structures and confirmed that the directional shift in AMI under 452 fentanyl-induced lower pain remained the same for VPM (n=4) and SCC (n=4) (see Figure 3H). 453 We also conducted tests for temporal shifts between time windows of various lengths, ranging 454 from 15ms to 35ms, and found that 30ms temporal shift (as shown in the main results) provided the best 455 distinction between high and low pain states (see Figure S2 B). Furthermore, we explored whether 456 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint restricting the bandwidth to low gamma (30-60hz) or high gamma (60-80hz) would be more effective in 457 characterizing pain dynamics, but found the 30-80hz range to be the most optimal (see Figure S2 D). 458 Dynamical past dependency captures pain dynamics 459 AMI based on the IEI parameter is a novel method that captures dynamical past dependency and 460 nonlinearity of a single LFP signal, and we demonstrated how this successfully characterized pain 461 dynamics. Next, we aimed to understand which features of this method contributed to an effective 462 characterization of pain dynamics. To explore this further, we applied a variety of conventional metrics 463 using datasets obtained from the spontaneous condition (as these included pain VAS scores) and aimed to 464 identify key features to capture the pain dynamics. 465 466 Figure 4. Comparison of metrics in characterizing pain dynamics. 467 A. AMI positively correlated with pain VAS (left), and was indifferent to whether the pain was increasing (middle) 468 or decreasing (right) in time. B. AMI positively correlated against pain VAS individually for patient P4 (left) and P5 469 (middle, right), and left structures showed a stronger correlation than the right structure. For ease of reading, for P4 470 the x-axis (pain VAS) is ordered by pain scores; whereas for P5, the x-axis follows a temporal order. The red line 471 marks the highest pain score. C. Changes in PSD (z-score) using the Morlet wavelet (left) and the HHT method 472 (right) following spontaneous lower pain are plotted across all frequency bands. The shaded area represents the 473 standard error of the mean. In PSD-HHT at 2 Hz, significance was observed at p < 0.05 without multiple corrections. 474 No significant correlation was found between PSD and pain VAS (right). D. Changes in entropy following 475 spontaneous lower pain were not evident (left), nor did they exhibit any correlation with the pain VAS (right). E. 476 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint Changes in auto-correlation following spontaneous lower pain were not observed across all conditions, although the 477 auto-correlation value exhibited a positive correlation with pain VAS scores. 478 Abbreviation: L VPM, left ventral-parietal-medial of the thalamus; L SCC, left subgenual cingulate cortex; R SCC, 479 right subgenual cingulate cortex; PSD, power spectrum density; HHT, Hilbert-Huang Transform; R, auto-480 correlation coefficient; τ , Kendall’s tau correlation coefficient; L-H, Low-High; S, spontaneous lower pain 481 condition; F, fentanyl based lower pain condition; P propofol administered condition; **, significance p<0.01; *, 482 significance p<0.05; ns, non-significant. 483 First, we found that AMI doesn't merely characterize pain dynamics (i.e., AMI reduction 484 corresponds to lower pain) in categorical manner (as shown in Figure 3A), but that it is correlated with 485 the pain VAS (τ =0.50; p<0.01). We further validated this for both rising (τ =0.64; p<0.01) and falling pain 486 trends (τ =0.79; p<0.01) and confirmed that AMI mirrors the overall pain intensity level (Figure 4A). We 487 also compared this per patient P4 and P5 – the two patients who experienced a wide range of pain VAS 488 during the recording – and confirmed that their AMI values positively correlated with their pain VAS at a 489 marginal (p<0.1) level on the left side of the recorded structure for P4 and P5 (P4: τ = 0.80, p=0.08; P5: τ 490 = 0.69, p=0.02). P5’s right side structure did not show significant positive correlation (τ = 0.50, p=0.14) 491 (Figure 4B). 492 To compare with other metrics, we also examined whether power spectrum density (PSD) could 493 effectively capture pain dynamics. This assessment was carried out using both the Morlet wavelet and the 494 Hilbert-Huang Transform (HHT) methods. It is noteworthy that while neither of these methods reflects 495 past dependency, the latter does capture nonlinearity. Our findings revealed that neither of the PSD 496

demonstrated effective characterization of pain dynamics. Furthermore, we conducted a 497 correlation between the normalized (z-scored) Morlet wavelet PSD values and the pain VAS scores but 498 did not find any significant patterns (τ = -0.23, p=0.11) (Figure 4C). This suggests that the inclusion of 499 dynamic past dependency may be a critical factor in characterizing pain dynamics. 500 To assess the importance of nonlinearity, we compared entropy values between paired epochs, as 501 the entropy metric captures the nonlinearity feature but does not account for past dependency. Here, our 502 findings showed that entropy did not effectively capture pain dynamics across conditions (S, p=0.57; F, 503 p=0.25; P, p=0.50), nor did it correlate with the pain VAS (τ = 0.03, p=0.9), This indicates that 504 nonlinearity is not a critical feature (see Figure 4D) to characterize pain dynamics. To investigate the 505 importance of past dependency, we compared the auto-correlation values between paired epochs. This 506 metric does not reflect nonlinearity features but does capture past dependency within a specific time 507 window (in this case, 5 seconds). As a result, we did not observe auto-correlation to effectively capture 508 pain dynamics under all conditions (S, p=0.86; F, p=0.20; P, p=0.38). However, we did find a positive 509 correlation between auto-correlation values and the pain VAS (τ = 0.33, p=0.02) (Figure 4E). This finding 510 suggests that past dependency may play a critical role in characterizing pain dynamics, although the 511 presence of linearity might diminish its effectiveness when compared to the current AMI method. 512 Table 2. Comparison of metrics for characterizing pain. 513 Metrics Reflect non-linearity? Reflect past dependency? Characterize pain dynamics? Correlate with pain VAS? AMI o o o * o ** PSD-Morlet x x x x PSD-HHT o x x nc Entropy o x x x Auto-Correlation x o x o * (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint Abbreviation: AMI, auto-mutual information; PSD-Morlet, power spectrum density using Morlet wavelet; PSD-514 HHT, power spectrum density using Hilbert-Huang Transform method; Notation: o, yes; x, no; nc, not computed; 515 **, p<0.01; *, p<0.05. 516

517 Our goal was to apply a novel method - AMI based on IEIs – on the LFP signals at the VPM, 518 PVG, and SCC to characterize pain dynamics. As we conceptualized AMI as a measure of regularity in 519 the signal (indicated by increased past dependency), we discovered a link between pain and regularity 520 (high AMI) within the theta/alpha band. Additionally, we found that opioid-based pain medication led to 521 regularity (high AMI) within the gamma band. This discovery marks the first to dynamically characterize 522 spontaneously occurring pain fluctuations with human electrophysiological signals. 523 In addition, our novel computational approach has several merits that potentially make it a 524 valuable pain biomarker for application in closed-loop DBS treatment. First, we demonstrate the 525 sufficiency of using signals from a single site to characterize pain. This was possible because we recorded 526 from a neural hub that processes diverse pain-related information 23–25. Notably, we could not achieve this 527 when calculating AMI in the prefrontal cortical sites (shown in Figure S2 D). Therefore, our findings 528 differ from viewpoints favoring the use of signals from multiple sites as an effective way to characterize 529 pain 49. Additionally, the AMI method we introduced is computationally more efficient than typical 530 machine-learning methods 46,50, which require extensive data and training (lasting hours to days). Note, 531 we merely needed about 10 seconds worth of LFP time series to compute the AMI. Our method also 532 surpasses the efficiency of prior neural studies using the AMI metric (e.g., 34,35,51 ), which typically 533 depend on signal amplitude values. Here, we reduce the LFP time series to binary spike trains for 534 computing the inter-event intervals (IEIs) thus simplifying the computation. Considering these strengths, 535 we believe this novel methodology holds significant potential to serve as a robust biomarker for the 536 development of closed-loop DBS for chronic pain. 537 We further explored which features of this AMI method most contributed to its effective 538 characterization in pain dynamics. To do this, we employed various metrics that capture different features 539 of the signal and assessed their effectiveness in characterizing pain dynamics. As a result, we found the 540 auto-correlation metric to be the most suitable alternative for characterizing pain dynamics. As auto-541 correlation measures the extent of past dependency in the signal, we consider that a decrease in past 542 dependency (i.e., irregularity) in the signal plays a critical role in effectively characterizing lower pain. 543 Moreover, considering our current AMI method's superior performance compared to auto-correlation, our 544 speculation is that its effectiveness stems from its nonparametric and nonlinear approach in quantifying 545 past dependency. 546 We interpret these results with the conceptualization that irregularity in the LFP signal represents 547 diverse (pain- and non-pain-related) information flow from nearby structures, whereas greater regularity 548 reflects a more focused flow of pain-related information. Specifically, we consider that in the absence of 549 pain, the signal becomes more irregular and less reliant on past data (i.e., lower AMI), as diverse 550 information about the external world is constantly communicated, updating the signal moment-to-moment. 551 Conversely, during the experience of pain, analogous to how an individual in pain concentrates on the 552 painful sensation, the communication primarily involves pain-related information originating from within 553 the internal body. In this case, the signal tends to hold onto its internal pain information, impeding the 554 input about the external world or other non-pain-related factors. With this viewpoint, we hypothesized 555 that higher pain would be associated with more signal regularity and lower pain with signal irregularity, 556 and we found this to be the case within the theta/alpha frequency band. Our results align with prior 557 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint research, where lower variability was found in fMRI BOLD signals from somatosensory cortex and PVG 558 during pain states 20, and in EEG signals among sensory-over responsive individuals 41. Interestingly, we 559 find an inverse trend between the spectral power and AMI, as we note an increase in theta/alpha power 560 coincided with lower AMI (irregularity) upon pain reduction, and a decrease in gamma band power 561 (Figure S3 A, B) coincided with higher AMI (regularity) following fentanyl administration. However, we 562 also note that these two metrics (power spectrum and AMI) do not exhibit a significant negative 563 correlation when directly compared against each other (Figure S3 C). From this, we conjecture that when 564 new information flows into a site within a specific oscillatory bandwidth, the signal becomes more 565 irregular (lower AMI), while more oscillatory power emerges within that frequency band. To gain a 566 deeper understanding of this hypothesis, investigating the neural connectivity between the pain processing 567 hub structures and neighboring sites would provide a more comprehensive view of this concept of 568 information flow. 569 We also interpret our findings on past dependency within the context of complexity. Although 570 infrequent, AMI has been applied in electrophysiological studies as a measure of complexity, representing 571 rich interconnectedness and information content 31. Specifically, higher AMI would indicate lower 572 complexity, and has been linked to various diseases (e.g., Alzheimer’s disease 34,35, schizophrenia 32, 573 idiopathic dilated cardiomyopathy 33, multiple organ dysfunction syndrome36). Regarding our results, we 574 interpret that higher AMI (found in higher pain instances) is associated with a less complex brain state, 575 indicating sparsely interconnected information flow within the theta/alpha band. Interestingly, we 576 observed that pain relief induced by opioid medication increases AMI within the gamma band, implying 577 reduced complexity. This suggests that opioids may introduce an artificial information flow leading to 578 less complex interactions within the gamma band. While drawing definitive conclusions is premature, we 579 argue that our findings provide insights into the frequency bandwidths worth exploring in future studies 580 involving opioid pain medication, to deepen our understanding of the complexity of brain dynamics in 581 such contexts. 582 We recognize that the study's small sample size, consisting of 5 participants, limits the 583 generalizability of our findings. However, it's crucial to underscore that each individual data point 584 originates from an extensive dataset (ranging approximately from 500 to 2500 datapoints) acquired from 585 each patient. Consequently, these individual data points possess notable statistical power and serve as 586 meaningful summary statistics on their own. Furthermore, we acknowledge the variability in epoch 587 lengths, both among patients and within paired epochs, which results in an uneven sample size from 588 which the mean AMIs are derived. This approach was unavoidable as we aimed to maintain an 589 observational approach within a naturalistic context. Still, we assert that the AMI values, rooted in the 590 mean of a stochastic process with more than 100 datapoints, would not be too sensitive to the sample size 591 (compared to statistical parameters like variance). 592 In future studies, it would also be valuable to conduct this research within a more controlled 593 setting, especially when investigating the impact of pain medication. Additionally, the calculation of AMI 594 in our study involved selecting parameters such as frequency bandwidth, bin size, and time window 595 overlap, which would affect the resulting AMI values. We used a data-driven approach to derive these 596 selections, going through multiple trial-and-error iterations (depicted in Figure S3 and S4 B) to discover 597 optimal parameters that effectively differentiated the pain states. For subsequent research, we suggest 598 applying these selected parameters to a larger group of participants and confirm whether these parameter 599 choices hold true for a wider range of individuals. We also recognize that the patients in this study 600 monitored pain in an intraoperative setting, which is different from a more naturalistic everyday setting. 601 To validate our findings, assessing LFP in a naturalistic environment, especially in patients with 602 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint externally recorded DBS implants, would enhance the authenticity of pain experiences in chronic pain 603 patients. Lastly, gaining a deeper understanding of the underlying neural activity could be achieved by 604 simultaneously recording single-unit spike activity alongside population-based activity. This approach 605 could eventually lead to the development of a model-based understanding of this phenomenon. 606 In summary, our study presents a novel AMI computational method for characterizing pain 607 dynamics through the analysis of electrophysiological signals in the VPM, SCC, and PVG regions. This 608 discovery offers significant potential for advancing the development of closed-loop DBS treatments for 609 individuals with chronic pain, ultimately leading to enhanced care and improved outcomes for a wide 610 range of chronic pain patients. 611 Conflict of interest statement 612 The authors have no conflicts of interest to declare. 613 This research was funded by UCLA Training in Neurotechnology Translation (TNT) NIH T32 fellowship 614 (JR). 615 616 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint

617 1. Gureje, O. et al. The relation between multiple pains and mental disorders: Results from the World 618 Mental Health Surveys. doi:10.1016/j.pain.2007.05.005. 619 2. Mantyselka, P. et al. Pain as a reason to visit the doctor: a study in Finnish primary health care. 620 3. Kwan, C. L. et al. Abnormal forebrain activity in functional bowel disorder patients with chronic 621 pain. Neurology 65, 1268–1277 (2005). 622 4. Hashmi, J. A. et al. Shape shifting pain: chronification of back pain shifts brain representation 623 from nociceptive to emotional circuits. A JOURNAL OF NEUROLOGY doi:10.1093/brain/awt211. 624 5. Foss, J. M., Apkarian, A. V. & Chialvo, D. R. Dynamics of pain: Fractal dimension of temporal 625 variability of spontaneous pain differentiates between pain states. J Neurophysiol 95, 730–736 626 (2006). 627 6. Farmer, M. A. et al. Brain Functional and Anatomical Changes in Chronic Prostatitis/ Chronic 628 Pelvic Pain Syndrome HHS Public Access. J Urol 186, 117–124 (2011). 629 7. Parks, E. L. et al. Brain activity for chronic knee osteoarthritis: Dissociating evoked pain from 630 spontaneous pain. European Journal of Pain 15, 843.e1-843.e14 (2011). 631 8. Zhang, Z. H., Dougherty, P. M. & Oppenheimer, S. M. Monkey insular cortex neurons respond to 632 baroreceptive and somatosensory convergent inputs. Neuroscience 94, 351–360 (1999). 633 9. Hutchison, W. D., Davis, K. D., Lozano, A. M., Tasker, R. R. & Dostrovsky, J. O. Pain-related 634 neurons in the human cingulate cortex. (1999). 635 10. Davis, K. D., Kucyi, A. & Moayedi, M. The pain switch: an “ouch” detector. 156, 2164–2166 636 (2015). 637 11. Meeker, T. J. et al. Tonic pain alters functional connectivity of the descending pain modulatory 638 network involving amygdala, periaqueductal gray, parabrachial nucleus and anterior cingulate 639 cortex. Neuroimage 256, 119278 (2022). 640 12. Baliki, M. N., Mansour, A. R., Baria, A. T. & Apkarian, A. V. Functional reorganization of the 641 default mode network across chronic pain conditions. PLoS One 9, (2014). 642 13. Melzack, R. & Casey, K. L. Sensory, motivational, and central control determinants of pain: a new 643 conceptual model. The skin senses 1, 423–443 (1968). 644 14. Melzack, R. PAIN From the gate to the neuromatrix. Pain Suppl 6, 121–126 (1999). 645 15. Kucyi, A. et al. Enhanced medial prefrontal-default mode network functional connectivity in 646 chronic pain and its association with pain rumination. Journal of Neuroscience 34, 3969–3975 647 (2014). 648 16. May, E. S. et al. Dynamics of brain function in patients with chronic pain assessed by microstate 649 analysis of resting-state electroencephalography. Pain 162, 2894–2908 (2021). 650 17. Baliki, M. N., Schnitzer, T. J., Bauer, W. R. & Apkarian, A. V. Brain Morphological Signatures 651 for Chronic Pain. PLoS One 6, 26010 (2011). 652 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 18. Baliki, M. N. et al. Chronic Pain and the Emotional Brain/i6 : Specific Brain Activity Associated 653 with Spontaneous Fluctuations of Intensity of Chronic Back Pain. 26, 12165–12173 (2006). 654 19. Kucyi, A. & Davis, K. D. The Neural Code for Pain: From Single-Cell Electrophysiology to the 655 Dynamic Pain Connectome. Neuroscientist 23, 397–414 (2017). 656 20. Rogachov, A. et al. Regional brain signal variability: A novel indicator of pain sensitivity and 657 coping. Pain 157, 2483–2492 (2016). 658 21. Kucyi, A., Salomons, T. V. & Davis, K. D. Mind wandering away from pain dynamically engages 659 antinociceptive and default mode brain networks. Proc Natl Acad Sci U S A 110, 18692–18697 660 (2013). 661 22. Apkarian, A. V., Krauss, B. R., Fredrickson, B. E. & Szeverenyi, N. M. Imaging the pain of low 662 back pain: Functional magnetic resonance imaging in combination with monitoring subjective pain 663 perception allows the study of clinical pain states. Neurosci Lett 299, 57–60 (2001). 664 23. Bittar, R. G. et al. Deep brain stimulation for pain relief: A meta-analysis. 665 doi:10.1016/j.jocn.2004.10.005. 666 24. Wu, D., Wang, S., Stein, J. F., Aziz, T. Z. & Green, A. L. Reciprocal interactions between the 667 human thalamus and periaqueductal gray may be important for pain perception. Exp Brain Res 232, 668 527–534 (2014). 669 25. Kashanian, A., Tsolaki, E., Pouratian, N. & Bari, A. A. Deep Brain Stimulation of the Subgenual 670 Cingulate Cortex for the Treatment of Chronic Low Back Pain. Neuromodulation: Technology at 671 the Neural Interface 2021, (2021). 672 26. Valet, M. et al. Distraction modulates connectivity of the cingulo-frontal cortex and the midbrain 673 during pain—an fMRI analysis. Pain 109, 399–408 (2004). 674 27. Park, B. yong, Lee, J. J., Kim, H. J., Woo, C. W. & Park, H. A neuroimaging marker for predicting 675 longitudinal changes in pain intensity of subacute back pain based on large-scale brain network 676 interactions. Scientific Reports 2020 10:1 10, 1–11 (2020). 677 28. Horn, A. et al. Lead-DBS v2: Towards a comprehensive pipeline for deep brain stimulation 678 imaging. (2018) doi:10.1016/j.neuroimage.2018.08.068. 679 29. Rebain, R., David Baxter, G. & McDonough, S. A Systematic Review of the Passive Straight Leg 680 Raising Test as a Diagnostic Aid for Low Back Pain (1989 to 2000). Spine (Phila Pa 1976) 27, 681 388–395. 682 30. Shannon, C. E. A mathematical theory of communication. The Bell system technical journal 27, 683 379–423 (1948). 684 31. Fraser, A. M. & Swinney, H. L. Independent coordinates for strange attractors from mutual 685 information. Physical Review 33, 1134 (1986). 686 32. Na, S. H., Jin, S. H., Kim, S. Y. & Ham, B. J. EEG in schizophrenic patients: mutual information 687 analysis. Clinical Neurophysiology 113, 1954–1960 (2002). 688 33. Palacios, M. et al. Changes of autonomic information flow due to idiopathic dilated 689 cardiomyopathy. Physiol Meas 28, 677–688 (2007). 690 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 34. Jeong, J., Gore, J. C., & Peterson, B. S. Mutual information analysis of the EEG in patients with 691 Alzheimer’s disease. Clinical neurophysiology 112, 827–835 (2001). 692 35. Abasolo, D., Escudero, J., Hornero, R., Gomez, C. & Espino, P. Approximate entropy and auto 693 mutual information analysis of the electroencephalogram in Alzheimer’s disease patients. Med 694 Biol Eng Comput 46, 1019–1028 (2008). 695 36. Hoyer, D. et al. Interactions between short-term and long-term cardiovascular control mechanisms. 696 Chaos 17, (2007). 697 37. Melia, U. et al. Auto-Mutual Information Function for Predicting Pain Responses in EEG Signals 698 during Sedation. in XIII Mediterranean Conference on Medical and Biological Engineering and 699 Computing 2013 623–626 (2013). doi:10.1007/978-3-319-00846-2. 700 38. Torres, E. B. Two classes of movements in motor control. 269–283 (2011) doi:10.1007/s00221-701 011-2892-8. 702 39. Torres, E. B., Cole, J. & Poizner, H. Motor output variability, deafferentation, and putative deficits 703 in kinesthetic reafference in Parkinson’s disease. Front Hum Neurosci 8, 823 (2014). 704 40. Ryu, J., Choi, J.-W., Niketeghad, S., Torres, E. B. & Pouratian, N. Irregularity of instantaneous 705 gamma frequency in the motor control network characterize visuomotor and proprioceptive 706 information processing. doi:10.1101/2023.07.28.551050. 707 41. Ryu, J., Bar-Shalita, T., Granovsky, Y., Weissman-Fogel, I. & Torres, E. B. Personalized 708 Medicine Personalized Biometrics of Physical Pain Agree with Psychophysics by Participants with 709 Sensory over Responsivity. J. Pers. Med 11, 93 (2021). 710 42. Ryu, J. & Torres, E. B. Motor Signatures in Digitized Cognitive and Memory Tests Enhances 711 Characterization of Parkinson’s Disease. Sensors 22, (2022). 712 43. Ryu, J. & Torres, E. Toward interpretable digital biomarkers of walking and reaching in 713 Parkinson’s disease. Wearable Technologies 3, e21 (2022). 714 44. Whitten, T. A., Hughes, A. M., Dickson, C. T. & Caplan, J. B. A better oscillation detection 715

robustly extracts EEG rhythms across brain state changes: The human alpha rhythm as a 716 test case. Neuroimage 54, 860–874 (2011). 717 45. Huang, N. E. et al. The empirical mode decomposition and the Hubert spectrum for nonlinear and 718 non-stationary time series analysis. Proceedings of the Royal Society A: Mathematical, Physical 719 and Engineering Sciences 454, 903–995 (1998). 720 46. Chen, Z. S. Decoding pain from brain activity Decoding pain from brain activity. 721 47. Ploner, M. & May, E. S. Electroencephalography and magnetoencephalography in pain research - 722 Current state and future perspectives. Pain 159, 206–211 (2018). 723 48. Hemington, K. S., Wu, Q., Kucyi, A., Inman, R. D. & Davis, K. D. Abnormal cross-network 724 functional connectivity in chronic pain and its association with clinical symptoms. Brain Struct 725 Funct 221, 4203–4219 (2016). 726 49. Shirvalkar, P., Veuthey, T. L., Dawes, H. E. & Chang, E. F. Closed-Loop Deep Brain Stimulation 727 for Refractory Chronic Pain. 12, 1–13 (2018). 728 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 50. Shirvalkar, P. et al. nature neuroscience First-in-human prediction of chronic pain state using 729 intracranial neural biomarkers. Nature Neuroscience | 26, 1090 (2023). 730 51. Erpelding, N. & Davis, K. D. Neural underpinnings of behavioural strategies that prioritize either 731 cognitive task performance or pain. Pain 154, 2060–2071 (2013). 732 733 734 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint Supplemental Figures 735 736 737 Figure S1. Pain VAS alongside the corresponding spectrogram and epoch specifications for each patient. 738 A. Patient P1. (top) Pain VAS is denoted by blue, while leg angular velocity along the Y-axis is represented by red. 739 Initially, despite the VAS being marked at 4, we labeled it as "V0" because the patient verbally reported no pain. A 740 subsequent increment of +1 in pain VAS was noted as "V1”. B. Patient P2. The patient's trackpad trajectory is 741 depicted in cyan, and the smoothed VAS is shown in blue. C. Patient P3. This patient verbally indicated the 742 absence of pain and was accordingly administered propofol, followed by fentanyl as per clinical protocol. D. Patient 743 P5. The patient’s incremental fluctuations in pain are shown in blue and angular velocity along the Y-axis are 744 shown in red. Subsequently the patient performed an unrelated experiment for approximately 15 and then was 745 administered with fentanyl due to experienced pain. 746 Abbreviation. V, pain VAS; F, fentanyl administration; P, propofol administration; R SCC, right subgenual 747 cingulate cortex; L SCC, left subgenual cingulate cortex; R VPM, right ventral parietal medial of the thalamus; L 748 VPM, left ventral parietal medial of the thalamus; L PVG, left periventricular gray. 749 750 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 751 Figure S2. Change in AMI under varying parameters. 752 A. When applying a theta/alpha bandpass filter to the LFP, we computed the change in AMI resulting from pain 753 relief across different time-delayed windows (ranging from 25ms to 275ms). It was determined that the 175ms time 754 delay yielded the most optimal characterization of pain relief. B. For the LFP filtered within the gamma frequency 755 band, the change in AMI stemming from fentanyl intake was best characterized using a time delay of 30ms. C. The 756 computation of change in AMI from pain relief involved exploring different frequency bandwidths within the lower 757 frequency range. Notably, the theta/alpha band emerged as the most suitable for characterizing pain relief from both 758 spontaneous and fentanyl-induced instances. D. Similarly, change in AMI from pain relief was calculated across 759 different frequency bandwidths within the higher frequency range. Results indicated that the 30-80Hz band provided 760 the most effective characterization of pain relief from fentanyl. E. We extended the AMI comparison to the 761 prefrontal cortex. No significant change in AMI was detected for either the theta/alpha or gamma bands upon 762 experiencing pain relief. 763 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 764 Figure S3. Change in power spectrum from pain relief and its relation to AMI. 765 A. The change in power spectrum (z-score) from pain relief is depicted using Morlet wavelets for all types of pain 766 relief, and categorized into spontaneous, fentanyl-induced, and propofol-induced. The shaded region represents the 767 standard error of the mean. Instances of statistical significance at p<0.05, without corrections for multiple 768 comparisons, are indicated in green. However, no statistical significance was observed across all frequency ranges 769 when corrections were applied. B. The change in power spectrum (z-score) from pain relief is presented using the 770 Hilbert Huang Transform method. This method exhibits a comparable pattern to the Morlet wavelet approach, 771 showing an observable increase in power spectrum within the lower frequency range (around the theta band) and a 772 decrease in the gamma band range. Despite this consistency, no statistical significance was detected upon applying 773 corrections for multiple comparisons. C. When comparing the power spectrum density (PSD) and AMI values, no 774 direct correlation is evident between the two metrics. 775 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint 776 Figure S4. Miscellaneous analyses on AMI computation. 777 A. The magnitude of change in AMI does not exhibit correlation with VAS measurements. Specifically, the change 778 in AMI resulting from pain relief does not correlate with the lower pain VAS, as evidenced by the relationship 779 between Δ AMI and the corresponding lower pain VAS (τ = -0.20, p=0.24), nor does it correlate with higher pain 780 VAS (τ = -0.20, p=0.24). Furthermore, the change in AMI does not show correlation with the magnitude of change 781 in pain VAS, as indicated by the correlation between Δ AMI and Δ VAS (τ = -0.24, p=0.18). B. Between two time 782 windows with time delay, observations revealed that when the two windows experienced no time delay (indicated by 783 100% overlap), the mean AMI was the highest. Nonetheless, when the overlap reduced, alterations in time delay did 784 not result in significant changes to overall AMI values. Thus, while we identified the optimal time delay for AMI 785 computation in this study, we believe that variations in time delay would not substantially impact the overall 786 outcomes. C. In nonlinear dynamical systems, a steeper decline in the rate of AMI with time delay serves as a 787 normalized complexity measure 31. In our present study, we observed an exponential decrease in the AMI function 788 within the lag range of up to 250ms. Consequently, a steeper decline would correspond to lower AMI values. 789 790 791 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted February 27, 2024. ; https://doi.org/10.1101/2024.02.22.581655doi: bioRxiv preprint