The sample used for this post hoc analysis was selected from the INTREPIBD-TIMEBASE cohort [23,24]. Inclusion criteria were broad to have a representative sample of the population: (1) aged 18 to 75 years; (2) patients with a current or past diagnosis of BD according to DSM-5 criteria and confirmed with a semistructured interview (SCID-5-RV [Structured Clinical Interview for DSM-5 Disorders, Research Version]); (3) sustained Young Mania Rating Scale (YMRS) and Hamilton Depression Rating Scale (HDRS) scores of ≤7 for at least 8 weeks; and (4) willing and able to give consent. The exclusion criteria were as follows: (1) presence of severe cardiovascular or neurological medical conditions that may lead to autonomic dysfunction, ongoing cardiovascular arrhythmia, or pacemaker use; (2) current comorbid substance use disorder as per DSM-5 criteria, excluding nicotine substance use disorder; (3) current comorbid psychiatric disorder causing significant interference with symptoms; (4) current pharmacological treatment involving beta-blockers; and (5) ongoing pregnancy.

The sample consisted of 80 patients; 49 (61.3%) were female participants, with a mean age of 46.9 (SD 14.1) years (Table 1). The mean number of years since illness onset was 15.3 (SD 12.1), most had been between 2 and 5 lifetime episodes (37/75, 49.3%), and 65% (52/80) had BD type I. All participants scored very low on the YMRS and HDRS, as they needed to be euthymic for inclusion in this analysis. Of the 80 participants, 51 (63.8%) were currently receiving antipsychotic treatment and had a high level of functioning (Table 1).

The sample was recruited at the Bipolar and Depressive Disorder Unit of Hospital Clínic de Barcelona between 2021 and 2024 [24]. Clinicians involved in the study invited eligible patients to participate during routine clinical visits. Thus, a nonprobability convenience sampling approach was used, based on clinician referral and patient willingness to participate. The number of individuals who declined participation was not systematically recorded. Participation was voluntary, and no financial compensation was provided. All participants were aware of the nature of the data collection, and no blinding was implemented. No a priori sample size calculation was performed, as this was a post hoc exploratory analysis of an existing cohort. From the total cohort of 362 participants, we selected 80 for this analysis. The remaining 282 were excluded due to the presence of acute symptoms (YMRS or HDRS scores >7) or did not meet the diagnostic criteria for BD.

This study was reviewed and approved by the Hospital Clinic of Barcelona Ethics and Research Board (HCB/2021/104). The study was also conducted in accordance with the ethical standards established in the Declaration of Helsinki and Good Clinical Practice. All participants provided written informed consent prior to inclusion in the study. All data were stored anonymously and encrypted on servers complying with the General Data Protection Regulation and the Health Insurance Portability and Accountability Act regulations. No monetary compensation was given for their participation in the study. No identifiable individual patient information is presented.

During an initial interview, the diagnosis was confirmed through the Structured Clinical Interview (SCID) for DSM-5, administered by trained mental health professionals [25]. The following clinical data were obtained from the interview and the patients’ medical records: age, sex assigned at birth, medical comorbidities, pharmacological treatment, years since BD onset, number of lifetime episodes, predominant polarity, diagnosis (BD type I, BD type II, unspecified BD, or schizoaffective disorder), and functionality (Social and Occupational Functioning Assessment Scale [SOFAS]). Pharmacological treatment was binarily categorized to reflect the presence or absence of each major pharmacological drug group, for example, current antipsychotic use. All antipsychotic medications used in the sample are reported in Table S1 in Multimedia Appendix 1. The number of lifetime episodes was grouped into clinically interpretable categories (first episode, 2‐5, 6‐9, and ≥10 episodes) to reduce recall bias and enhance data reliability during acquisition. Although initially reported by participants, episode counts were supplemented with information extracted from the electronic health record. Predominant polarity was defined as manic predominant polarity (MPP) if two-thirds or more of lifetime episodes were manic, depressive predominant polarity (DPP) if two-thirds or more of lifetime episodes were depressive, and undetermined predominant polarity (UPP) if it did not meet either of these criteria [26]. Euthymic state was defined as a score of 7 or fewer points in the Young Mania Rating Scale (YMRS) [27] and the Hamilton Depression Rating Scale (HDRS), which have well-established psychometric properties [28,29].

Participants wore an Empatica E4 wristband (Empatica Inc) for 48 hours. The E4 is a research-grade device specifically designed and validated for ambulatory assessment of the ANS in daily life. Participants followed their daily activities and removed the wristband only when there was a risk of water contact (eg, showering and doing dishes) to avoid harming the hardware. The E4 continuously records EDA at a 4 Hz sampling rate through dry sensors and HR through a photoplethysmography sensor with a 64 Hz sampling rate. In addition, the device records triaxial acceleration (32 Hz) and skin temperature (4 Hz), which were not used in this analysis. As mentioned previously, EDA provides an established index of SNS activity and is widely used in stress and arousal research [5,30]. Cardiac measures derived from blood volume pulse include HR and HRV (and additional parameters) as an index of parasympathetic activity or balance between both branches (Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, 1996). Both HR and EDA are prone to artifacts due to sensor displacement or movement, among other reasons. A motion-artifact removal algorithm is applied by default to the HR signal; therefore, we did not add any additional filtering after extracting the data from the E4 Manager platform. Conversely, EDA is not preprocessed by Empatica and, consequently, the signal might be subject to motion artifacts and poor wrist contact, which can affect its quality. To detect artifacts and obtain clean data, we submitted the raw EDA signal to the E4 dashboard preprocessing pipeline [30-32]. The pipeline uses several quality indicators to extract clean and artifact-free data. For the current analyses, we used HR, average EDA (SCL), and nonspecific SCR (expressed as EDA peaks per minute) [33]. We specifically computed the EDA peaks per minute variable (EDA peaks) by counting the number of peaks for each minute according to the timestamp. Participants with more than 75% of artifact-contaminated signal were excluded from the analyses. Additionally, for each model, participants with missing data for the corresponding illness course variable were not included. Wearable missing data were evaluated through visual inspection, and analyses were performed using available valid data without imputation.

Statistical outliers were not explicitly defined, as extreme values were considered to have been largely addressed during the preprocessing stage through artifact removal procedures. Data distributions were visually inspected to assess plausibility and skewness. No additional data transformations were applied.

This study follows an observational design. To describe the sample, we present the mean (SD) for continuous variables and the number and percentage for categorical variables. To avoid collapsing EDA and HR high-granularity data into a daily mean, or forcing them into a linear model, we applied generalized additive models (GAMs) to high-level time-series data. These types of models adapt to the data by combining splines and capturing their natural form. Two types of models were performed for each outcome (mean EDA, peaks, or HR) and each illness course variable (years since onset, type of BD, number of lifetime episodes, and predominant polarity). The years since onset variable was categorized into 4 groups based on the quartile distribution observed in our sample (1‐4, 5‐13, 14‐21, and ≥22 y). This data-driven categorization allowed us to model the effect of illness duration independently of the participant-specific effects, given that each participant had a unique value for years since BD onset. By grouping this continuous variable into quartiles, we reduced sparsity and ensured sufficient sample size within each category for stable estimation. For BD type, only BD type I and II were included in the models to avoid small-group overfitting due to the small sample of nonspecified BD and schizoaffective disorder. For the predominant polarity model, participants with schizoaffective disorder were discarded as there is no official definition of this clinical course variable in this population.

In the first type of model, the dependent variable (outcome) was determined by the illness course variable, the minute of the day, and a random effect. The minute of the day was included, as both HR [34-36] and EDA [37,38] vary systematically over the course of the day, following circadian-related patterns, and the random effect represented interindividual differences. The second model included the interaction between the clinical variable and the minute of the day in order to test whether the clinical variables disrupted the circadian patterns of HR or EDA. All models were adjusted for 3 covariates: age [39], sex [40,41], and antipsychotic medication [42], because of its reported effects on HR and EDA.

We checked for concurvity between variables, especially between age and years since BD onset, as we suspected a possible correlation between them. Results showed little to no concurvity between age and years since onset, but full concordance between the random effect and age (Table S2 in Multimedia Appendix 1). All models were fitted with the bam() function from the mgcv (v 1.9‐1) package in RStudio (v. 4.4.1). Smooth terms were fitted using penalized regression splines, with smoothing parameters selected via fREML. Effective degrees of freedom were estimated during model fitting and penalized to prevent overfitting. Temporal autocorrelation was addressed using an AR(1) structure (ρ=0.9), and between-participant variability was modeled using random effects. The predicted effects of the models were illustrated with the ggplot2 (v 4.0.0) package. All analyses were conducted considering a 2-tailed type I error rate of 5% and 95% CIs. When predictors were statistically significant in the models, we performed pairwise comparisons, adjusting for multiple comparisons using the Tukey method with the emmeans package (v 1.10.5).

Lastly, we performed a sensitivity analysis excluding participants with schizoaffective disorder, as it is usually considered a different clinical entity from BD [25]. The effects of years since BD onset, number of lifetime episodes, and diagnosis were tested with GAMs against mean EDA, peaks, and HR in this sample subset.

When the analysis was repeated after excluding patients with schizoaffective disorder, the effect of years since illness onset on EDA peaks gained strength and precision through narrow CIs (β=−1.06, 95% CI −1.71 to −0.41; P=.001; Table S5A in Multimedia Appendix 1 and Figure 2). Pair-to-pair comparison revealed that the group with more than 22 years of illness had significantly fewer EDA peaks per minute than the 2 younger illness groups, even after multiple comparison correction (Table S6 in Multimedia Appendix 1). Additionally, predominant polarity significantly predicted EDA peaks, with the UPP group having significantly fewer EDA peaks than the DPP group, even after adjusting for multiple comparisons (Table S7 in Multimedia Appendix 1 and Figure 3).

This exploratory, post hoc, cross-sectional study investigated whether BD illness course variables (years since illness onset, diagnosis, number of lifetime episodes, and predominant polarity) predicted EDA and HR in euthymic patients. Using GAMs, the main results were negative; however, a trend emerged suggesting that longer illness duration was linked to fewer EDA peaks per minute. This trend became statistically significant when excluding patients with schizoaffective disorder, pointing to a potential cumulative dampening of nonspecific SCR over time. Lastly, the absence of significant effects of the 4 analyzed clinical variables on HR and the circadian-related patterns of EDA and HR further supports the limited influence of these illness course variables on digital biosignals. Altogether, our findings suggest that the influence of BD illness course on EDA and HR in naturalistic settings is negligible, supporting the view that these measures could function primarily as state-dependent biomarkers that capture transient physiological fluctuations rather than stable illness traits.

To our knowledge, no previous studies have examined EDA in euthymic patients with BD under naturalistic conditions. In contrast, there is prior evidence showing associations between illness duration and HRV [14,43]. Although HR and HRV are closely related, their relationship is nonlinear; it has been described that in conditions of repeated exposure to stress, HR rises and HRV decreases [44,45]. Prior orthostatic testing observed that older age at onset correlated with lower very-low-frequency and high-frequency bands, reflecting reduced ANS activation [46]. However, this correlation arose only during the third position of an orthostatic test, which corresponds to the moment when the patient lies down after standing for 5 minutes. Thus, it arose under artificial testing conditions rather than in naturalistic environments [46]. In naturalistic settings, Ortiz et al [14] also reported a negative association between HRV, specifically root mean square of successive differences and the SD of all the normal R-R intervals, and age at BD onset [14]. Importantly, while previous work has focused exclusively on HRV, our study is the first to examine whether HR itself varies with illness course indicators in everyday conditions, providing novel evidence that HR does not show these associations in euthymic patients with BD.

Our negative findings regarding the associations between the analyzed biosignals and diagnosis and number of episodes cannot be compared with previous studies because, as far as we know, no other work has tested these variables. However, the ANS typically responds to stress and arousal, which are states, and is not as sensitive to traits, not even to sex in some studies [12,47]. Consequently, it is plausible that diagnosis, which reflects clinical phenotypes that influence behavior and symptomatology, does not substantially affect ANS functioning per se. In addition, the number of lifetime episodes could be considered a reasonable proxy for neuroprogression, as this theory is based on the presumed toxicity of episodes [22]. However, it does not capture episode intensity or duration and relies primarily on patient recall; therefore, it is not very robust, and in addition, extreme groups are underrepresented in our sample. By contrast, years since BD onset is more objective and captures the physiological load accumulated over time. Overall, our predominantly negative findings support the notion that digital biosignals are state-dependent rather than trait-dependent, following the trend observed in other peripheral biomarkers, such as inflammatory markers and neurotrophic factors, which tend to normalize during euthymic phases in BD, reaching levels similar to those of healthy controls [48-51]. From a clinical perspective, this reinforces the potential role of wearable-derived digital biosignals as tools for monitoring dynamic changes in mental states and relapse risk, rather than cross-sectional indicators of illness severity [11,12,52,53].

We also did not find any interaction between illness course variables and the biosignals’ circadian patterns.. HR and tonic EDA circadian patterns are highly stable due to their dependence on core biological factors such as the sleep-wake cycle, hormonal rhythms, and thermoregulation [54]. As these remain stable, the influence of a long-term clinical factor could become negligible. This would explain why, in our models, the effect of minute of day on EDA and HR is large, while the effects of clinical variables, if present at all, were small (Tables S4-S7 in Multimedia Appendix 1). Moreover, previous studies have only reported alterations in the mean or basal tone of HR or EDA, but not in their circadian patterns [10,14]. Thus, circadian patterns of sympathetic arousal, indexed by EDA and cardiac activity, remain largely preserved and primarily driven by sleep-wake cycles and central circadian mechanisms, which likely explains the absence of significant clinical and time-of-day interactions [55].

From the sensitivity analyses, we observed that excluding participants with schizoaffective disorder strengthened the group differences. Specifically, the group with more than 22 years of BD onset showed significantly fewer EDA peaks per minute than the group with less than 13 years since onset. The consistency of findings across sensitivity analyses suggests a modest effect of clinical heterogeneity on EDA. In addition, the predominant polarity model reached statistical significance, with the UPP group displaying fewer EDA peaks than the DPP group. However, we do not believe these findings reflect cross-sectional diagnostic differences but rather the clinical characteristics of these specific participants. Moreover, the predominant polarity model had the lowest statistical power, as this variable was available for only 35 participants. Thus, although the UPP group exhibited fewer EDA peaks than the DPP group, even after adjusting for multiple comparisons, these results warrant replication in a larger sample. Predominant polarity was obtained primarily relying on patient recall, making it particularly vulnerable to recall and cognitive biases and, therefore, the least reliable of the 4 illness course variables explored. Lastly, the predominant polarity groups did not differ significantly in sex, years since BD onset, or antipsychotic use (Table S8 in Multimedia Appendix 1).

The association we observed between longer illness duration and fewer EDA peaks per minute could be explained by several underlying mechanisms. First, habituation refers to the gradual reduction in physiological responsiveness that occurs when repeated or prolonged sympathetic activation leads the system to produce weaker reactions over time [43]. If we understand BD as a possible cause of repeated and prolonged sympathetic activation [44], habituation could explain the reduction in skin conductance reactivity and, therefore, the fewer EDA peaks reported in our sample. This mechanism is known as allostatic load, whereby prolonged exposure to stress leads to a lower phasic response and a higher basal tone [45,46]. Second, long-term antipsychotic use or chronic benzodiazepine intake can attenuate sympathetic activity [47,48]. Although we controlled for current antipsychotic use, we were unable to control for cumulative antipsychotic exposure across the lifetime. Third, the findings could also be attributed to neuroprogression hypotheses within BD staging frameworks [49]. Focusing on Berk’s staging proposal, which emphasizes episode number, as opposed to Kapczinski’s model, which prioritizes functional impairment, the repetition of episodes may induce biological changes in frontolimbic networks, oxidative stress, and inflammation [22,50]. However, this neuroprogression effect should be more clearly detected through the number of lifetime episodes rather than years since onset, as it is postulated that the number of episodes is what drives these changes, but this effect was not observed in our data.

Our study has several limitations. First, although restricting the sample to euthymic patients is justified, it also introduces methodological constraints. Selecting only euthymic patients reduces the variability of HR and EDA signals, as these mainly fluctuate during acute states. This reduced variability may limit the ability to detect associations with illness course variables, particularly if their effects are small. Additionally, if the impact of illness course variables is primarily expressed during acute states, as previous literature suggests [16], analyzing only euthymic patients may further reduce the capacity to detect such effects. Finally, our design precludes examining potential interactions between illness course variables and mood states.

Another limitation was the categorization of years since illness onset into quartiles, which may result in loss of information and reduced variability. However, this approach allowed more stable group comparisons and improved clinical interpretability, and we were still able to detect an effect of this variable. Additionally, the post hoc exploratory design may increase the probability of false-positive findings because the hypothesis is not predefined [56]. Although models were adjusted for current antipsychotic treatment, other potential confounding factors, such as lifetime pharmacological exposure, cumulative medication effects, and comorbid medical conditions, were not systematically accounted for. These factors may independently influence autonomic nervous system activity and could therefore partially explain the observed associations. The relatively small sample size in some subgroup analyses and the short recording period (48 h) may further limit the robustness of the findings. In particular, short monitoring durations may not fully capture intraindividual variability in physiological signals, especially in BD, where mood episodes evolve over longer timescales, such as 2 weeks for depressive episodes. Consequently, these factors may reduce the ability to detect subtle physiological patterns and limit the generalizability of the results. Our analysis focused on the tonic and phasic components of EDA and HR, but including HRV metrics and other physiological or behavioral digital biosignals may provide complementary information [57]. Finally, the lack of external validation of our models and the absence of comparable methodologies in the literature hinder direct comparison with other studies.

On the other hand, our study has several strengths, such as the use of high-resolution, real-world physiological data collected via wearable devices. Wearables allow a granular and objective assessment of physiological digital biosignals that is easy to implement in real clinical settings [52,58]. Also, the application of advanced statistical methods, such as GAMs, enabled us to capture nonlinear relationships and circadian-related patterns in the data. Additionally, the inclusion of a diverse set of clinical variables (number of episodes, diagnosis, predominant polarity, and years since onset) allowed for a comprehensive exploration of illness course and its impact on digital biosignals.

On the other hand, our study has several strengths, such as the use of high-resolution, real-world physiological data collected via wearable devices. Wearables allow a granular and objective assessment of physiological digital biosignals that is easy to implement in real clinical settings [52,58]. Also, the application of advanced statistical methods, such as GAMs, enabled us to capture nonlinear relationships and circadian-related patterns in the data. Additionally, the inclusion of a diverse set of clinical variables (number of episodes, diagnosis, predominant polarity, and years since onset) allowed for a comprehensive exploration of illness course and its impact on digital biosignals.

Future research should aim to replicate and extend our findings to determine whether BD illness course variables truly have no meaningful impact on ANS functioning, as our data suggest. The relationship between HRV and illness course variables should also be further explored. Additionally, although illness duration was associated with physiological signatures compatible with greater allostatic load, this finding neither confirms nor refutes the neuroprogression hypothesis, and rigorously designed longitudinal studies explicitly targeting neuroprogressive mechanisms are still needed to address this question. In parallel, efforts should shift toward more precise characterization of the specific states to which digital biosignals such as EDA and HR are most sensitive, so they can be effectively integrated into clinical management of BD and other conditions where autonomic arousal is relevant. A state-focused approach may ultimately enable the development of transdiagnostic, stress-tailored interventions that improve functioning and long-term outcomes.

Our findings support a predominantly state-dependent role of digital biosignals in BD and suggest that, in clinically stable patients, EDA and HR alone may have limited value as cross-sectional indicators of illness course, although subtle autonomic alterations associated with years since illness onset cannot be excluded and warrant confirmation in larger, longitudinal studies [2,3,52,59]. Despite the limitations, our study demonstrates the potential of wearable technologies for objective, real-world monitoring in BD, specifically as state markers. Finally, we emphasize the importance of further research to rigorously validate these digital biosignals and support their translation into clinically relevant biomarkers [52].