Interoception of breathing and its relationship with anxiety

1We perceive the world through our body. Although questions regarding how we sense and interpret our external environment (exteroception) have been highly prominent across centuries of research, the importance and cognitive mechanisms of monitoring our internal environment have only more recently gained traction within the neuroscience community. “Interoception,” the perception of our body and inner physiological condition, constitutes a fundamental component of cerebral processes for maintaining bodily homeostasis. However, it has also been suggested to play a wider role within the systems governing emotion, social cognition, and decision making. An impaired ability to monitor bodily signals has also been hypothesized to exist across a host of psychiatric illnesses, in particular anxiety. As sympathetic arousal is a reflexive response to a perceived threat, many symptoms associated with anxiety manifest themselves in the body (such as a racing heart or shortness of breath). Conversely, perceiving bodily states compatible with sympathetic arousal in the absence of external triggers can itself induce anxiety. Miscommunications between the brain and body are thus thought to represent a key component of anxiety, where bodily sensations may be under-, over-, or mis-interpreted, which can initiate and perpetuate symptoms of anxiety.

2Studying interoception is not without significant challenges, as bodily signals are both noisy and difficult to safely manipulate. Controlled manipulations of respiratory processes represent a promising way to address these challenges: suitable experimental setups allow dynamic yet safe changes in visceral aspects of respiration as one interoceptive modality. Furthermore, given the vitally important role of breathing for survival, respiratory changes are highly salient. Indeed, labored, unsatisfied, unexpected, or uncontrolled breathing can itself be perceived as a dangerous and debilitating interoceptive threat. Beyond respiratory diseases, aversive breathing symptoms have been noted to be particularly prevalent in individuals suffering from psychiatric conditions such as anxiety and panic disorder.

3Work toward conceptualizing interoceptive dimensions has provided us with a framework to integrate the growing body of interoception research. Instead of treating interoception as a single entity, studies now consider both different sensory channels (e.g., organ-specific and humoral signals) and cognitive layers of interoceptive processing. These layers encompass multiple levels, ranging from metrics of afferent signal strength at “lower” levels (using techniques such as heartbeat evoked potentials;; ) and psychophysical properties (such as measuring perceptual sensitivity;;; ) to psychological and cognitive components at “higher” levels. Notable domains within these higher levels include attention toward bodily signals, static and dynamic beliefs and models of body state, and insight into both our interoceptive abilities and the accuracy of our interoceptive beliefs (“metacognition”). Importantly, research into dynamic models of body state has also connected the interoceptive literature to that of learning, where influential (Bayesian) theories of inference about the external world, such as predictive coding, have been extended to interoception and used to propose how the brain may build models of the changing internal environment. Although theoretical models have been proposed, realistic synthetic data have been produced, and initial learning models have now been fit to empirical cardiac data, concurrent measures of dynamic brain processes during interoceptive learning have not yet been demonstrated.

4Here, we build on these conceptual advances and assess the relationship between anxiety and breathing-related interoception across the multiple hierarchical levels of processing. Importantly, although many theoretical proposals have been put forward as to how anxiety may interrupt the brain’s processing of dynamic (trial-by-trial) interoceptive predictions and/or prediction errors, these are as yet untested. Therefore, within a rigorous assessment profile, we included neuroimaging of a novel breathing-related interoceptive learning paradigm, providing the first empirical insight into the brain activity associated with interoceptive predictions and prediction errors. Furthermore, it is not yet known how alterations in interoceptive learning may relate to previously identified relationships between anxiety and lower level breathing sensitivity, higher level beliefs, or metacognition. Therefore, we aimed to both assess how anxiety is related to dynamic interoceptive learning and additionally provide a unifying perspective on anxiety and breathing-related interoception across the hierarchical levels of interoceptive processing. We adopted a multi-modal experimental approach to investigate multiple levels of breathing-related interoceptive processing, including low-level perceptual sensitivity and related higher level metacognition via the filter detection task (FDT), subjective interoceptive beliefs via questionnaires, and trial-by-trial interoceptive learning and related brain activity in a novel breathing learning task (BLT). Both the FDT and trial-by-trial behavioral and functional magnetic resonance imaging (fMRI) data from the BLT were analyzed with separate computational models. All tasks were performed by two matched groups of low- and moderate-anxiety individuals, allowing us to evaluate the relationship between anxiety and each level of breathing-related interoceptive processing across the hierarchy, from sensitivity to metacognition.

5Below we present the results from each of our task modalities: questionnaires, a breathing perception and metacognition task (the FDT), and a novel interoceptive learning task (the BLT), in which group-wise comparisons between each of the measures of interest were conducted. The results from the questionnaires and FDT are contextualized by previous findings related to anxiety, while the results from the BLT were validated against an additional unseen dataset and the relationship with anxiety was assessed. We then present the results of a combined multi-modal analysis, in which we compared our measures both within and across task modalities. The principal components (PCs) of these measures were identified to formalize and asses any shared variance between measures, which spanned multiple dimensions of breathing-related interoceptive processing.

6The group summaries and comparisons for each of the affective and interoceptive questionnaires (excluding the trait anxiety score that was used for group allocation) are displayed in Figure 1. The group summary values and statistics presented in text are either mean ± SE when values were normally distributed and thus compared using unpaired t tests or median ± inter-quartile range (IQR) when values were not normally distributed and thus compared using Wilcoxon rank-sum tests. Scores from all questionnaires of affective symptoms used were found to be highly significantly different between groups with low and moderate trait anxiety: individuals with moderate levels of trait anxiety demonstrated higher state anxiety (Spielberger State Anxiety Inventory [STAI-S] mean ± SE: low anxiety, 25.7 ± 0.7; moderate anxiety, 34.1 ± 1.2; t = −6.1, p < 0.01), higher levels of anxiety disorder symptoms (Generalized Anxiety Disorder Questionnaire [GAD-7] median ± IQR: low anxiety, 1.0 ± 2.0; moderate anxiety, 4.0 ± 3.0; Z = −5.9, p < 0.01), greater anxiety sensitivity (Anxiety Sensitivity Index [ASI] mean ± SE: low anxiety, 6.8 ± 0.8; moderate anxiety, 18.4 ± 1.5; t = −6.9, p < 0.01), and higher levels of depression symptoms (Center for Epidemiologic Studies Depression Scale [CES-D] median ± IQR: low anxiety, 6.5 ± 3.0; moderate anxiety, 14.0 ± 6.0; Z = −6.0, p < 0.01).

7The interoceptive questionnaires we used measured “positively minded” interoceptive awareness, overall body awareness, breathing symptom catastrophizing, and breathing symptom vigilance. All except breathing-related vigilance were also found to be significantly different between groups. Individuals with moderate levels of trait anxiety demonstrated reduced “positively minded” interoceptive awareness (Multidimensional Assessment of Interoceptive Awareness [MAIA] mean ± SE: low anxiety, 109.1 ± 4.6; moderate anxiety, 84.6 ± 3.7; t = 4.2, p < 0.01) and greater reports of overall body awareness (Body Perception Questionnaire [BPQ] median ± IQR: low anxiety, 66.0 ± 68.0; moderate anxiety, 104.0 ± 52.0; Z = −2.5, p = 0.01) in line with previous research. Additionally, elevated levels of breathing-related catastrophizing were observed in the moderate anxiety group (Pain Catastrophizing Scale [PCS-B] median ± IQR: low anxiety, 3.5 ± 11.0; moderate anxiety, 14.0 ± 17.0; Z = −3.3, p < 0.01), while no statistically significant difference was observed for breathing-related vigilance (Pain Vigilance Awareness Questionnaire [PVQ-B] mean ± SE: low anxiety, 16.3 ± 2.2; moderate anxiety, 21.4 ± 2.5; t = −1.5, p = 0.13). Results for sub-component scores and additional questionnaires can be found in Figure S1.

8The group summaries and comparisons for each of the FDT measures are displayed in Figure 2. The FDT output includes the number of filters at perceptual threshold (indicative of perceptual sensitivity, where a greater number of filters indicates lower perceptual sensitivity), decision bias (with c < 0 indicating a tendency to report the presence of a resistance), metacognitive bias (calculated from average confidence scores), and metacognitive performance (reflecting the congruence between confidence scores and performance accuracy). Individuals with moderate levels of trait anxiety demonstrated both lower perceptual sensitivity (in line with previous findings;; ) (filter number median ± IQR; low anxiety, 3.0 ± 2.0; moderate anxiety, 4.0 ± 2.0; Z = −2.4, p = 0.01) and lower metacognitive bias (average confidence score median ± IQR: low anxiety, 6.7 ± 2.2; moderate anxiety, 6.2 ± 2.1; Z = 2.0, p = 0.02) than those with low levels of anxiety, with a similar level of metacognitive performance (Mratio median ± IQR: low anxiety, 0.8 ± 0.0; moderate anxiety, 0.8 ± 0.1; Z = 0.7, p = 0.23). Decision bias was not found to be different between the groups (decision bias c parameter mean ± SE; low anxiety, 0.15 ± 0.06; moderate anxiety, 0.05 ± 0.06; t = 1.1, p < 0.14). The relationship between greater anxiety and reduced confidence is consistent with results previously observed in the exteroceptive (visual) domain, where decreased confidence related to individual levels of both anxiety and depression.

9A visual depiction of the BLT and example fitted trajectories for prediction certainty and prediction error magnitude are provided in Figure 3. When comparing the plausibility of the three alternative models (a Rescorla-Wagner [RW] model; a two-level hierarchical Gaussian filter [HGF2], and a three-level hierarchical Gaussian filter [HGF3]) using random-effects Bayesian model selection, no single model was found to have a protected exceedance probability (PXP) greater than 90% (RW:HGF2:HGF3 PXP = 0.30:0.40:0.30; Table S4). Therefore, as specified in our analysis plan (https://gitlab.ethz.ch/tnu/analysis-plans/harrison_breathing_anxiety), we conducted our model-based analysis using the conceptually most simple model (the RW model), in accordance with Occam’s razor. Importantly, the finding that none of the models demonstrated a PXP greater than 90% does not provide any absolute statement about the quality of these models. Rather, this finding indicates that none of the chosen models is conclusively superior to the others in explaining the data. To ensure that the chosen model (RW) provided an adequate explanation of the data, we compared it to a “null model” (i.e., where the choices were due to chance and not related to any associative learning mechanism) using a likelihood ratio test. We found that in 58 of the 60 participants, the RW model fit the behavior significantly better than the null model, demonstrating that the chosen model captured important aspects of their behavior. The two participants (one from each anxiety group) who did not show model fits above chance were excluded from any further model-based analyses and comparisons.

10To further establish the adequacy of our chosen model to explain learning behavior in this novel interoceptive learning task, we completed a model validation on 15 additional held-out datasets for the BLT. These participants were not pre-selected for any particular level of anxiety (see STAR Methods for details). A logistic regression was conducted to assess whether the model prediction trajectory from the original data was able to significantly explain the prediction decisions made by the 15 unseen participants in the validation sample. A representative prediction trajectory from the original 60 participant model fits (the trajectory from the participant with the closest learning rate to the mean) was used in this regression, as well as an intercept term. The beta estimate for the original prediction trajectory was 3.1 ± 0.3 (t = 12.1, p = 1.0 × 10−33), denoting a highly significant ability of the trajectory to predict unseen data. The beta estimate for the intercept term was −0.2 ± 0.1 (t = −3.1, p = 1.7 × 10−3). For a qualitative representation of model fits for both the original and validation data, see Figure S3.

11Both model-based and behavioral parameter comparisons are presented in Table 1. For the estimated model parameters, no difference was observed between the groups for either learning rate (α) or inverse decision temperature (ζ) (Table 1). Results from parameter comparisons between groups including the excluded participants can be found in Table S5. For the subjective measures, no difference was observed between the groups for breathing difficulty ratings, while the task-induced anxiety ratings were significantly greater in those with moderate anxiety (Table 1). Additionally, no difference in any physiological measures were observed (Table S1), nor in relative head motion during the task (average root-mean-square displacement ± SD: low anxiety, 0.17 ± 0.10 mm; moderate anxiety, 0.18 ± 0.07 mm; Wilcoxon rank-sum p = 0.91).

12The overall and between-group BLT brain activity analysis results are displayed in Figures 4 and 5. In the analysis for the entire field of view, activations associated with breathing-related prediction certainty and prediction errors across all participants are shown in Figure 4. Dorsolateral prefrontal cortex (dlPFC), anterior insula (aIns), anterior cingulate cortex (ACC), and middle frontal gyrus (MFG) all demonstrated significant deactivations with overall prediction certainty (i.e., averaged across trials with positive and negative prediction certainty; Figure 4A). In contrast, aIns, ACC, MFG, and the periaqueductal gray (PAG) demonstrated significant activations with overall prediction error values (i.e., averaged across trials with positive and negative prediction errors; Figure 4B). A small number of differences due to valence (differences between positive and negative outcomes) were found for prediction errors but not prediction certainty, with negative prediction errors associated with deactivations of left dlPFC and activations of left posterior insula (Figure 4B). Although no main effect of anxiety group was observed, an interaction effect was found using the region-of-interest (ROI) analysis between valence and groups for predictions in the bilateral aIns (Figure 5). In contrast, no group or interaction effects were found for prediction errors. Brain activity associated with inspiratory resistance is provided in Figure S6 for comparison with previously published results.

13First, the key measures from each of the different modalities were combined into a multi-modal correlation matrix. This analysis allowed us to assess the relationships both within and across task modalities as well as across levels of breathing-related interoceptive processing. The full correlation matrix of all 16 included measures is displayed in Figure 6A and Table S7. To briefly summarize, the strongest across-task modality correlations were found between all affective and interoceptive questionnaires (Figure 6A). Concerning affective questionnaires and the FDT measures, state anxiety was weakly correlated with the FDT perceptual threshold, decision bias, and metacognitive bias, while anxiety sensitivity was additionally weakly related to metacognitive bias. Depression scores were also weakly related to the FDT perceptual threshold. Between the interoceptive and the FDT measures, breathing-related catastrophizing was weakly related to metacognitive performance on the FDT. Last, between the FDT and aIns activity, metacognitive performance was strongly related to the peak aIns activity associated with negative (i.e., resistance-related) prediction errors, while metacognitive bias was weakly related to aIns activity associated with negative (i.e., resistance-related) prediction certainty. Non-parametric correlations (using Spearman’s rho values) produced highly consistent results and are presented in Table S7.

14Finally, to assess the extent of shared variance across interoceptive measures, the multi-modal data matrix was then subjected to a principal-component analysis (PCA). This analysis allowed us to delineate how many underlying dimensions may exist within the data, as well as which measures were most strongly associated with anxiety. Two PCs were found to be significant, where the variance explained with each component was above the 95% confidence interval of its null distribution. The properties of each of these significant components are displayed in Figures 6B and 6C. The first PC demonstrated a highly significant (p < 1 × 10−11) difference in scores between the anxiety groups. Correspondingly, the greatest coefficients within the first PC were from the affective measures of depression scores, state anxiety, anxiety sensitivity, and anxiety disorder scores. Additionally, breathing-related catastrophizing and negative interoceptive awareness also had strong coefficient values, followed by negative metacognitive bias (i.e., lower confidence scores), body perception scores (from the BPQ), and negative metacognitive performance (i.e., lower metacognitive performance). In contrast, the second PC demonstrated a weak difference (p = 0.05) in scores between the anxiety groups. This component had the highest coefficient scores from the peak aIns activity related to positive and negative prediction certainty, as well as negative coefficients for negative prediction errors, metacognitive performance, and positive prediction errors.

15Interoceptive abilities are thought to be tightly linked to affective properties such as anxiety. Here we have provided a unifying analysis by characterizing this relationship across multiple interoceptive levels in the breathing domain, including novel findings of altered brain activity within the aIns when processing dynamic breathing predictions. This study is the first to demonstrate brain activity related to dynamic interoceptive learning and, specifically, activity in the insula that is related to breathing-related prediction error and prediction certainty in a trial-by-trial fashion. Notably, aIns activity related to the certainty of predictions about breathing resistance was found to differ between trait anxiety levels. Furthermore, this study is also among the first to simultaneously tackle multiple levels of interoceptive processing, using breathing as a salient and accessible channel of body perception. The tasks used reflect the broad range of targeted processes; not only were questionnaires used that spanned affect and body perceptions, but behavioral data from two different tasks were assessed by separate computational models. These analyses allowed formal assessments of both breathing-related interoceptive learning and metacognition, including the first computational assessment of trial-by-trial learning in the interoceptive domain, as well as applying state-of-the-art models of metacognition to interoception of breathing. Our multi-modal approach revealed that not only is the relationship between breathing-related interoception and trait anxiety broad, it is most strongly detected (i.e., greatest PCA weights; Figure 6) at the higher levels of interoceptive processing, which includes specific subjective measures of interoceptive beliefs (often termed “interoceptive sensibility”; ) followed by metacognitive aspects of breathing perceptions. Notably, the peak aIns brain activity associated with breathing-related interoceptive learning appeared to be largely independent of other interoceptive measures, with the exception of negative prediction error-related brain activity and metacognitive performance.

16Beyond consequences at single levels of interoceptive processing, here we aimed to assess how the relationship with trait anxiety may cross multiple interoceptive levels related to breathing. Using PCA (with permutation testing) allowed us to identify any components that share common variance within our multi-modal dataset and additionally assess the relative contribution of our measures to each dimension (Figure 6B). Here we found that all affective qualities loaded strongly onto the first principal component, with the greatest additional contributions from subjective measures of negative interoceptive awareness and breathing-related catastrophizing. General body awareness and the breathing-related metacognitive measures (bias and performance) were the next largest contributors to this shared variance, followed by the perceptual sensitivity and decision bias parameters, and last peak aIns activity from the BLT. These results suggest that the relationship with anxiety is particularly prominent at the level of subjective interoceptive beliefs in the breathing domain, which are thought to exist at the higher levels of interoceptive space, followed by metacognitive insight into breathing perception. In comparison, the relationship of trait anxiety to lower level properties such as interoceptive sensitivity appear to be present but less prominent in the breathing domain. However, it must be noted that quantifying these higher interoceptive levels may be less noisy in comparison with measuring psychophysical properties such as sensitivity, and thus the relationship with anxiety might be most easily detected rather than being inherently stronger.

17Although strong relationships were observed between affective qualities and many of our interoceptive breathing measures, a sparse number of correlations were found between interoceptive measures themselves, and in particular across task modalities (Figure 6A). These findings support the idea that there are potentially separable levels of breathing-related interoception, as proposed. The only notably strong cross-modal relationship was found between metacognitive performance and aIns activity, for which greater insight into breathing sensitivity correlated with greater aIns activity during negative prediction errors. This relationship may reflect a previously proposed contribution of error processing toward metacognitive awareness, whereby deviations between actual and predicted bodily inputs are propagated to metacognitive areas via interoceptive brain structures such as the aIns.

18Many theories surrounding anxiety have hypothesized an important role of altered predictions regarding upcoming threat and in particular interoceptive threat in the aIns. Although numerous studies have used inspiratory resistive loads to evoke threatening interoceptive stimuli, the BLT approach presented here is, to our knowledge, the first investigation of dynamic (trial-by-trial) brain activity associated with model-based interoceptive predictions and prediction errors for respiration. By fitting an associative learning model to each participant’s behavioral responses, we could quantify both the certainty of predictions and magnitude of prediction errors on each trial. We could then compare both the parameter estimates and the brain activity associated with these computational quantities, with a particular focus on the aIns and PAG (Figure 4). Here, we observed evidence for a relationship between anxiety and aIns reactivity to threat valence in the prediction domain (Figure 5). Specifically, while individuals with low trait anxiety demonstrated greater aIns deactivation that scaled with predictions of breathing resistance compared with no resistance, the opposite was true in individuals with moderate trait anxiety (creating an interaction effect). This demonstrates a shift in the aIns processing of threat valence with different levels of anxiety, in line with hypothesized differences in brain prediction processing. In comparison, no anxiety group differences or interactions were found in the breathing prediction error domain, contrasting with some previously proposed hypotheses.

19Beyond the aIns and independent of anxiety, multiple (and largely consistent) proposals have been made regarding which brain networks might be involved in processing interoceptive predictions and prediction errors. These proposed networks are loosely hierarchical in structure and typically assign metacognitive processes to higher cortical areas (e.g., prefrontal cortex [PFC]), while interoceptive predictions are thought to originate from regions that may serve as an interface between interoceptive and visceromotor function (e.g., aIns and ACC). In these concepts, prediction errors have two different roles: on one hand, they are thought to be sent up the cortical hierarchy of interoceptive regions in order to update predictions in aIns and ACC; on the other hand, they are thought to determine regulatory signals, sent from visceromotor regions and brainstem structures such as the PAG to the autonomic nervous system and bodily organs.

20Although these theories have received considerable attention, there has been little empirical evidence thus far. In particular, we are not aware of any studies that have demonstrated, using a concrete computational model, trial-by-trial prediction and prediction error activity in interoceptive areas. Here, we report evidence of relevant computational quantities being reflected by activity within several areas of a putative interoceptive breathing network. Although activity related to trial-wise prediction certainty was found in higher structures such as dorsolateral PFC, ACC, and aIns, prominent prediction error responses were not only found in aIns and ACC but also in the midbrain PAG (Figure 4). Importantly, concerning predictions, widespread brain activity was found to be related mainly to prediction uncertainty, where blood-oxygen-level-dependent (BOLD) activity was decreased in proportion to increases in the certainty of predictions. Furthermore, it is notable that the aIns displayed both deactivation for more certain predictions and activation for greater prediction errors. This might reflect the proposed critical role of aIns in the representation and updating of models of bodily states, given that greater certainty (precision of beliefs) reduces and greater prediction errors increase belief (model) updating.

21Our PAG findings are of particular interest. Although the PAG has been previously noted during anticipation of certain breathing resistance stimuli and has been related to the precision of prior beliefs about placebo-induced reductions in pain intensity, here we observed that PAG activity did not appear to be related to the extent of prediction certainty toward upcoming breathing stimuli (Figure 4). Concerning prediction error activity in the PAG, this has previously been demonstrated in relation to pain; here, we found PAG activity in relation to the magnitudes of trial-wise interoceptive prediction errors (Figure 4B), consistent with a role of PAG in homeostatic control.

22Finally, overall prediction and prediction error-related activity did not appear to be dissociated between anterior and posterior insula cortices (respectively), as has been previously hypothesized. However, a small valence difference in prediction errors was observed, with negative prediction errors (the unexpected presence of inspiratory resistance stimuli) producing greater activity in the right posterior insula than positive prediction errors (the unexpected absence of inspiratory resistance stimuli; Figure 4). It is therefore possible that the representation of homeostatically relevant inputs in the posterior insula is enhanced for events that may negatively affect homeostasis. However, these results are the first demonstration of brain activity related to dynamic interoceptive prediction and prediction errors; furthermore, the functional images from this study do not have the resolution required for layer-specific identification of prediction and prediction error processing in the insula. Additionally, the represented prediction error-related activity may be specific to the breathing domain within interoceptive processing. This latter caveat of course also applies to the wider results presented here; as only one interoceptive channel (i.e., inspiratory resistances within the breathing domain) was tested, we cannot assume these results would translate to other interoceptive processes (e.g., related to cardiac or gastric states).

23The relationship between anxiety and breathing crosses multiple levels of the interoceptive hierarchy. In particular, anxiety and associated affective dimensions appear to be most strongly related to subjective negative body awareness and catastrophizing about breathing symptoms, followed by metacognitive measures related to breathing perception. Furthermore, a novel interaction between trait anxiety and valence was found within the aIns, associated with dynamic prediction certainty (but not prediction errors) of breathing-related interoceptive processing. More generally, this study provides the first empirical demonstration of brain activity associated with dynamic (trial-by-trial) interoceptive learning. These results provide new and comprehensive insights into how anxiety is related to levels of interoceptive processing in the breathing domain and provide evidence of brain activity associated with trial-wise predictions and prediction errors about bodily states in interoceptive breathing networks.

24Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Olivia Harrison (née Faull; faull@biomed.ee.ethz.ch).

25This study neither used any reagent nor generated new materials.

26Thirty individuals (pre-screened online for MRI compatibility, right handedness, non-smoking status, and no history of major somatic or psychological conditions) were recruited into each of two groups, either with very low anxiety (score of 20-25 on the Spielberger State-Trait Anxiety Inventory; STAI-T), or moderate anxiety (score >= 35 STAI-T). The resulting mean (±std) trait anxiety score for the low anxiety group was 23.2 ± 1.8 and for the moderate anxiety group 38.6 ± 4.6. Groups were matched for age and sex (15 females in each group), with mean (±std) ages of 25.4 ± 3.9 and 24.2 ± 5.0 years for low and moderate anxiety groups, respectively. Study numbers were based on a power calculation for a two-sided two-sample t test with an α-level of 5% (moderate effect size d = 0.5), where a power of 90% is achieved with 30 participants in each group. Behavioral data (not used in any other analyses) from an additional 8 participants served to determine model priors, with four participants (two females and two males) from each of the low and moderate anxiety groups. In this way, prior values could be drawn from a comparable group of participants as the main study sample. All participants signed a written, informed consent, and the study was approved by the Cantonal Ethics Committee Zurich (Ethics approval BASEC-No. 2017-02330).

27Behavioral data for an additional group of 15 individuals (12 female) were collected for model validation purposes. These participants were not pre-selected based on anxiety values (mean ± std trait anxiety = 38.9 ± 12.5), but were screened for non-smoking status and no history of major somatic or psychological conditions. Participants were aged (mean ± std) 23.1 ± 5.6 years. All participants signed a written, informed consent, and the study was approved by the New Zealand Health and Disability Ethics Committee (HDEC) (Ethics approval 20/CEN/168).

28Each participant completed three tasks over two testing sessions: a behavioral session that included questionnaires and a task probing interoceptive sensitivity and metacognition (the Filter Detection Task, or FDT), and a brain imaging session where the Breathing Learning Task (BLT) was paired with fMRI. Each of these tasks and analyses are described below, and all analyses were pre-specified in time-stamped analysis plans (https://gitlab.ethz.ch/tnu/analysis-plans/harrison_breathing_anxiety). The length of time between testing sessions (mean ± std) was 4 ± 3 days for all participants. Participants in the validation group completed the BLT in a behavioral session only.

29The main questionnaire set employed was designed to first capture subjective affective measures, and second both general and breathing-specific subjective interoceptive beliefs. The assignment of participants to groups was based on the Spielberger Trait Anxiety Inventory (STAI-T). Affective qualities that were additionally assessed included state anxiety (Spielberger State Anxiety Inventory; STAI-S), symptoms that are part of anxiety disorder (Generalized Anxiety Disorder Questionnaire; GAD-7), anxiety sensitivity (anxiety regarding the symptoms of anxiety; Anxiety Sensitivity Index; ASI-3), and symptoms of depression (Centre for Epidemiologic Studies Depression Scale; CES-D). To obtain self-reports of body awareness we used the Body Perception Questionnaire (BPQ), while the Multidimensional Assessment of Interoceptive Awareness Questionnaire (MAIA) was used to measure positive and ‘mindful’ attention toward body symptoms. We also measured breathing-related catastrophising using the Pain Catastrophising Scale (PCS-B), and breathing-related vigilance using the Pain Vigilance Awareness Questionnaire (PVQ-B) (in both questionnaires, the word ‘breathless’ or ‘breathlessness’ was substituted for ‘pain’). Finally, the following supplementary questionnaires were included to explore possible contributing factors (e.g., general positive and negative affect, resilience, self-efficacy and fatigue): Positive Affect Negative Affect Schedule (PANAS-T), Connor-Davidson Resilience Scale, General Self-Efficacy Scale, Fatigue Severity Scale (FSS). The STAI-T and CES-D were completed online as part of the pre-screening process; all other questionnaires were completed in the behavioral session at the laboratory.

30To systematically test properties of breathing perception and related metacognition, we utilized a perceptual threshold breathing task (the Filter Detection Task; FDT), described in detail elsewhere. The FDT was used to determine interoceptive perceptual sensitivity, decision bias, metacognitive bias (self-reported confidence) and metacognitive performance (congruency between performance and confidence scores) regarding detection of very small variations in an inspiratory load. In this task (outlined in Figure 2A), following three baseline breaths either an inspiratory load was created via the replacement of an empty filter with combinations of clinical breathing filters, or the empty filter was removed and restored onto the system (sham condition) for three further breaths. All filter changes were performed behind participants, out of their field of view. After each trial of six breaths, participants were asked to decide whether or not a load had been added, as well as reporting their confidence in their decision on a scale of 1-10 (1 = not at all confident in decision, 10 = extremely confident in decision). An adapted staircase algorithm was utilized to alter task difficulty until participants were between 60%–85% accuracy, and 60 trials were completed at the corresponding level of filter load once the threshold level had been identified (using a ‘constant’ staircase procedure, as described by. Respiratory threshold detection, metacognitive bias and interoceptive metacognitive performance have previously been linked to anxiety symptomology.

31To measure behavior and brain activity concerning the dynamic updating of interoceptive beliefs or expectations under uncertainty, a novel associative learning task was developed and employed during functional magnetic resonance imaging (fMRI). In this Breathing Learning Task (BLT), 80 trials were performed where on each trial two visual cues were paired with either an 80% or 20% chance of a subsequent inspiratory resistive load. Participants were explicitly told the probabilities that were being used, as well as that the cues were paired together – if one cue indicated an 80% chance of resistance then the other must indicate a 20% chance. Participants were also told that the cues could only swap their contingencies throughout the task, and could not act independently of each other. The number of trials was limited to 80 to ensure feasibility for participants. The number and structure of the cue contingency swaps (n = 4) within the 80 trials was chosen as simulated data demonstrated both parameter recovery and model identifiability of the three candidate learning models (see Figure S2 and ‘Quantification and Statistical Analysis’ details below).

32The visual information for the task was presented through the VisualStim system (Resonance Technology, Northridge, CA, US). As outlined in Figure 3, participants were required to predict (via button press) whether they would experience a breathing resistance following the presentation of one of the cues. The visual cues were counter-balanced such that each was first matched with an 80% chance of resistance for half the participants, as well as counter-balancing of whether the answer ‘yes’ to the prediction question was presented on the left or right of the screen. Following this prediction and a short (2.5 s) pause, a circle appeared on the screen to indicate the stimulus period (5 s), where participants either experienced inspiratory resistance (70% of their maximal inspiratory resistance, measured in the laboratory, delivered via a PowerBreathe KH2; PowerBreathe International Ltd, Warwickshire, UK) or no resistance was applied. Rest periods of 7-9 s were pseudo-randomized between trials.

33For the inspiratory resistances we used a mechanical breathing system that allows for remote administration and monitoring of inspiratory resistive loads (for technical details on resistance administration see previous work). The cue presentations were balanced such that half of all trials delivered the inspiratory resistance. Following an initial stable period of 30 trials, the stimulus-association pairing was swapped four times during the remainder of the 80 trials (i.e., repeated reversals; Figure 3). The trial sequence was pseudorandom and fixed across subjects to ensure comparability of the induced learning process. Following every stimulus, participants were asked to rate ‘How difficult was it to breathe”?, on a visual analog scale (VAS) from “Not at all difficult” to “Extremely difficult.” Immediately following the final trial of the task, participants were also asked to rate “How anxious were you about your breathing” on a VAS from “Not at all anxious” to “Extremely anxious.”

34Two representations of trial-wise quantities were employed for subsequent analyses of data from this task. First, a computational model (see below) provided dynamic estimates of both predictions and prediction errors on each trial. Second, a standard categorical approach represented trial-by-trial whether the subjects’ prediction decisions indicated the anticipated presence or absence of an upcoming inspiratory resistance, as well as unsurprising (i.e., following correct predictions) and surprising (i.e., following incorrect predictions) respiratory stimuli. The latter results are presented in the Supplemental information (Figures S4–S7).

35All analyses and hypotheses were pre-specified in an analysis plan (https://gitlab.ethz.ch/tnu/analysis-plans/harrison_breathing_anxiety). Each measure within each task modality was first compared between groups, and multiple comparison correction was applied within modalities. Pre-selection and comparisons across low or moderate anxiety scores allowed us to control for important factors such as age and sex, with equal numbers of men and women recruited into each group. Finally, a cross-modal analysis was performed on the key measures from each task. As the trait anxiety measure that was used to recruit participants into the separate groups was not included in this final analysis, we employed a correlation-based approach as the full spectrum of scores for each task were available.

36Group differences were tested individually for the 13 scores resulting from the 12 questionnaires, with all questionnaires included except the trait anxiety score that was used to screen participants and assign them to groups. The data that was used for group comparisons were first tested for normality (Anderson-Darling test, with p < 0.05 rejecting the null hypothesis of normally distributed data), and group differences were determined using either two-tailed independent t tests or Wilcoxon rank sum tests. For the questionnaires, Bonferroni correction for the 13 tests was applied, requiring p < 0.004 for a corrected significant group difference. Results with p < 0.05 not surviving correction are reported as exploratory for questionnaires as well as all other data. In a secondary exploratory step, group difference analyses were then conducted on the questionnaires’ subcomponent scores (22 scores); please see Figure S1.

37Breathing-related interoceptive sensitivity (i.e., perceptual threshold) was taken as the number of filters required to keep task performance between ∼60%–85% accuracy. Both decision bias and metacognitive performance from the FDT were analyzed using the hierarchical HMeta-d statistical model, as previously described. This model first utilizes signal detection theory to provide single subject parameter estimates for task difficulty (d’; not analyzed as performance is fixed between 60%–85% by design) and decision bias (c, akin to over- or under-reporting the presence of resistance with values below and above zero, respectively), as well as using a hierarchical Bayesian formulation of metacognitive performance (Mratio, calculated by fitting metacognitive sensitivity meta-d’, then normalizing by single subject values for d’). Finally, metacognitive bias was calculated as the average confidence scores across all analyzed trials.

38The four FDT parameters that were used for group comparison analyses were: sensitivity (filter number at perceptual threshold), decision bias (c), metacognitive bias (average confidence scores over threshold trials) and metacognitive performance (Mratio). These data were first tested for normality (Anderson-Darling test, with p < 0.05 rejecting the null hypothesis of normally distributed data), and group differences were determined using either two-tailed independent t tests or Wilcoxon rank sum tests. Importantly, as the single subject Mratio parameters were fit using a whole-group hierarchical model (with equal group numbers), standard statistical comparisons between groups for these parameters can also be performed. Hypothesized group differences were based on previous findings, where respiratory threshold detection level was hypothesized to be higher, perceptual decisions biased toward ‘yes’ (or deciding the resistance was present; denoted by more negative values for c), metacognitive bias to be lower and interoceptive metacognitive performance to be lower with greater anxiety. Bonferroni correction for the four tests was applied, requiring p < 0.013 for a corrected significant group difference.

39For the trial-by-trial analysis of behavioral data from the BLT, we considered three computational models that are routinely used for associative learning tasks. This included a Rescorla Wagner (RW) model (Equation 1) and 2 variants of the Hierarchical Gaussian Filter (HGF) with 2 or 3 levels (HGF2 and HGF3). While the RW model assumes a fixed learning rate, the HGF allows for online adaption of learning rate as a function of volatility. All learning models were paired with a unit-square sigmoid response model (Equation 2) and were implemented using the Hierarchical Gaussian Filter Toolbox (version 5.3) from the open-source TAPAS software (https://www.translationalneuromodeling.org/tapas/).

40Prior means and variances were determined using the distribution of maximum likelihood estimates fit across a holdout dataset (consisting of 8 participants who were distinct from the participants of our study). Individual fits and estimated prior densities of the free parameters are given in Table S2, as well as values for all remaining parameters of the three models. By adopting this procedure, prior densities were in a regime of the parameter space that is representative of the actual behavioral responses observed when participants performed the task. At the same time, the arbitrariness inherent to the specification of prior densities in non-hierarchical inference is reduced to a minimum.

41To demonstrate face validity of the three models considered in our model space (each with one parameter free: in RW, in HGF2, in HGF3), we assessed both parameter recovery and model identifiability for each. Data for 60 synthetic subjects were generated for each of the candidate models by randomly sampling values from the prior densities that were placed over the parameters of the perceptual model. This synthetic data was generated for different noise levels. Subsequently, maximum a posteriori (MAP) parameter estimates were obtained using the Brayden-Fletcher-Goldfarb-Shanno algorithm, as implemented in the HGF toolbox, to fit the synthetic datasets. This allowed us to quantify parameter recovery and model identifiability across three different noise levels for each of the candidate models. Parameter recovery of the perceptual parameters was assessed using Pearson’s correlation coefficient (PCC) and by visual inspection of simulated and recovered parameter values (Figure S2). Mean and standard deviation of estimated values (from the response model) were computed for every noise level. Model identifiability was quantified by calculating the proportion of correctly identified models using approximate log model evidence (LME) scores, and assessing whether the former was greater than the upper bound of the 90% confidence interval when assuming every model is equally likely a priori. For the resulting confusion matrices (Figure S2), we additionally computed the mean proportion of correctly identified models (balanced accuracy scores).

42The outlined procedure for assessing parameter recovery and model identifiability was repeated over 10 iterations with different seed values, to ensure robustness against any particular setting of the random number generator. The final results (PCCs for the perceptual parameters and values) for every given level of noise were calculated as the average over all iterations, and are presented in Table S3.

43Following MAP estimates of the 60 empirical datasets using each of the candidate models, the models were formally compared using random effects Bayesian model selection (BMS) as implemented in SPM12. BMS utilizes the LME to determine the most likely among a set of competing hypotheses (i.e., models) that may have generated observed data, and is robust to outliers. Our analysis plan had specified that a model would be chosen as the ‘winning’ model if it demonstrated a protected exceedance probability (PXP) greater than 90%. As explained in the Results section, none of our models reached this criterion (although simulations indicated that the proposed models could in principle be differentiated; see Figure S2 for details). We therefore applied the simplest of the models considered (i.e., the RW model), as pre-specified in our analysis plan (https://gitlab.ethz.ch/tnu/analysis-plans/harrison_breathing_anxiety).

44In our application of the RW model as a perceptual model, the update equation corresponded to a simple delta-learning rule with a single free parameter, the learning rate:where is the predicted probability for a specific outcome (encoded as 0 or 1) on trial, is the estimated outcome probability on the trial, α ∈ [0, 1] is a constant learning rate parameter, and is the prediction error magnitude at trial.

45The above perceptual model was paired with a unit-square sigmoid response model. This response model accounts for decision noise by mapping the predicted probability that the next outcome will be 1 onto the probabilities and that the agent will choose response 1 or 0, respectively:Here, represents the expressed decision of a subject given the cue (contingency pairs) on trial. The parameter captures how deterministically y is associated with. The higher, the more likely the agent is to choose the option that is more in line with its current prediction. The decision model uses the perceptual model indirectly via its inversion, given the trajectories of trial-wise cues and responses (see Figure 3).

46In our paradigm, trial-wise outcomes are categorical (resistance versus no resistance), which raises the question of how outcomes should be coded in the computational model. One way would be to model two trajectories, separately for resistance and no resistance outcomes, and indicate on any given trial whether the respective outcome has occurred (1) or not (0). However, due to the fixed coupling of contingencies in our paradigm (see above) – which the participants were explicitly instructed about – a computationally more efficient approach that requires only a single model is to code the outcome in relation to the cue. Here, we adopted this coding in “contingency space,” following the same procedure as in the supplementary material of Iglesias and colleagues. Specifically, due to the fixed coupling of contingencies in our paradigm (see above), we represented the occurrence of “no resistance” given one cue and the occurrence of “resistance” given the other cue as 1, and both other cue-outcome combinations as 0 (note that under the subsequent transformations described below, the resulting trajectories of predictions and prediction errors would remain identical if the opposite choice had been made).

47Group differences in model parameter estimates of learning rate and inverse decision temperature, as well as perception measures of stimulus intensity (averaged across all trials), breathing-related anxiety (rated immediately following the task) and prediction response times were compared following tests for data normality. Bonferroni correction for five tests was applied, requiring p < 0.01 for a corrected significant group difference. Results from additional exploratory models encompassing anxiety, depression and gender are reported in Table S6.

48Following random effects Bayesian model selection (BMS), the chosen model was examined in each participant with regard to whether it demonstrated an adequate fit. To this end, model fit in each individual was compared to the likelihood of obtaining the data by a ‘null model’ (i.e., due to chance) using the likelihood ratio test (lratiotest function) provided in MATLAB. The final behavioral and brain imaging analyses presented here were run without two subjects in which non-significant (p > 0.05) differences to randomness were encountered. To demonstrate the extent to which the chosen model (the Rescorla Wagner) captured important aspects of participant performance, the proportion of incorrect responses across participants at each trial was compared to the mean prediction error trajectory (Figure S3B).

49To further validate the application of the chosen model to the BLT, we compared the fitted model trajectory to unseen behavioral data from an additional 15 participants. For qualitative assessment, the mean prediction error trajectory from the original dataset was compared to the proportion of incorrect responses across these held-out participants at each trial (Figure S3C). For quantitative assessment, a logistic regression was conducted to assess whether the model prediction trajectory from the original data was able to significantly explain the prediction decisions made by the 15 participants in the validation sample.

50The trajectories of predictions and prediction errors estimated by the RW model were used to construct regressors representing computational trial-by-trial quantities of interest for subsequent GLM analyses. In order to investigate the salient effects of inspiratory resistance as an interoceptive stimulus, we separated trials into “negative” (occurrence of resistance) and “positive” (no resistance) events and represented these events by separate regressors in the GLM (see Figure 4). To achieve this, we first transformed both the original prediction and prediction error values (estimated in contingency space) back into the stimulus space, according to the cue presented at each trial:Here, and now represent the prediction and prediction error values in stimulus space, with representing maximal predictions of no resistance and maximal predictions of resistance. Similarly, represents maximal prediction errors of no resistance and maximal prediction errors of resistance (see Figure S4 for details).

51Second, trial-wise prediction values were then transformed to represent the deviation from maximally uninformed predictions (i.e., guessing), by taking the distance from 0.5 (see Equations 5 and 6). In the RW model prediction values are probabilities bounded by 0 and 1, hence the distance from ‘guessing’ (at 0.5) reflects the ‘certainty’ by which the absence or presence of respiratory resistance was predicted. This simple transformation enabled us to take into account the role of (un)certainty of predictions – which plays a crucial role in interoception-oriented theories of anxiety but, in contrast to Bayesian models, is not represented explicitly in the RW model. Specifically, separately for the two event types, we defined certainty of positive predictions (no resistance) and of negative predictions (resistance) as the absolute deviation from a prediction with maximum uncertainty (i.e., 0.5):Here, both and exist between 0 and 0.5, with values closer to zero indicating less certain predictions.

52Like predictions, prediction errors were also divided between positive (no resistance) and negative events (resistance) values. This was again determined as the absolute deviation from the mid-point of the prediction errors (i.e., 0):Here, both and exist between 0 and 1, with values closer to zero indicating smaller prediction errors. Note that this derivation gives prediction error values identical to those that would have been obtained by modeling two separate trajectories for resistance and no resistance outcomes (see above and).

53Physiological data were recorded at a sampling rate of 1000 Hz, and included heart rate, chest distension, pressure of expired carbon dioxide (PETCO2) and oxygen (PETO2), and pressure at the mouth (for equipment details see). In addition to the task, small boluses of a CO2 gas mixture (20% CO2; 21% O2; balance N2) were administered during some rest periods, allowing for de-correlation of any changes in PETCO2 from task-related neural activity, as previously described.

54Physiological noise regressors were prepared for inclusion into single-subject general linear models (GLMs, described below). Linear interpolation between PETCO2 peaks was used to form an additional CO2 noise regressor, which was convolved using a response function based on the haemodynamic response function (HRF) provided by SPM with delays of 10 s and 20 s for the overshoot and undershoot, respectively. Temporal and dispersion derivatives of this CO2 noise regressor were also included. An additional three cardiac- and four respiratory-related waveforms (plus one interaction term) were created using PhysIO. Four respiratory volume per unit time (RVT) regressors (delays: −5, 0, 5, 10) were created using the Hilbert-transform estimator in PhysIO, and convolution with a respiratory response function.

55MRI was performed using a 7 Tesla scanner (Philips Medical Systems: Achieva, Philips Healthcare, Amsterdam, the Netherlands) and a 32 channel Head Coil (Nova Medical, Wilmington, Massachusetts, United States of America). A T2∗-weighted, gradient echo EPI was used for functional scanning, using a reduced field of view (FOV) with an axial-oblique volume centered over the insula and midbrain structures. The FOV comprised 32 slices (sequence parameters: TE 30ms; TR 2.3 s; flip angle 75°; voxel size 1.5x1.5x1.5mm; slice gap 0.15mm; SENSE factor 3; ascending slice acquisition), with 860 volumes (scan duration 33 mins 9 s). A matched whole-brain EPI scan (96 slices) was immediately acquired following the task scan for registration purposes. Additionally, a whole-brain T1-weighted structural scan with 200 slices was acquired (MPRAGE, sequence parameters: TE 4.6ms; TR 10ms; segment-TR 3000ms; TI 1000ms; flip angle 8°; voxel size 0.8x0.8x0.8mm; bandwidth; 153.1Hz/Px; sagittal slice orientation). Finally, a task-free (resting-state) functional scan (250 volumes) was obtained, with participants instructed to keep their eyes open and fixating a white fixation cross on a black screen.

56MRI data analysis was performed using a combination of FSL version 6.0.1 (the Oxford Centre for Functional Magnetic Resonance Imaging of the Brain Software Library, Oxford, UK) and SPM12 (Statistical Parametric Mapping software, London, UK) as prespecified in our analysis plan (https://gitlab.ethz.ch/tnu/analysis-plans/harrison_breathing_anxiety). Image preprocessing was performed using FSL, including motion correction (MCFLIRT), removal of non-brain structures (BET), and high-pass temporal filtering (Gaussian-weighted least-squares straight line fitting; 100 s cut-off period). Independent component analysis (ICA) was used to identify noise due to motion, scanner and cerebrospinal fluid artifacts, and the timeseries of these noise components were entered into single-subject GLMs (described below) as nuisance regressors. The functional scans were registered to the MNI152 (1x1x1mm) standard space using a three-step process: 1) Linear registration (FLIRT) with 6 degrees of freedom (DOF) to align the partial FOV scan to the whole-brain EPI image; 2) Boundary-based registration (BBR; part of the FMRI Expert Analysis Tool, FEAT) with 12 DOF and a weighting mask of the midbrain and insula cortex to align the whole-brain EPI to T1 structural image; and 3) Non-linear registration using a combination of FLIRT and FNIRT to align the T1 structural scan to 1mm standard space. Functional MRI scans were resampled once into standard space with a concatenated warp from all three registration steps, and then spatial smoothing in standard space was performed using a Gaussian kernel with 3mm full-width half-maximum using the fslmaths tool.

57Single-subject estimates of the general linear model (GLM) were performed using SPM. A GLM was constructed for each of the participants, with a design matrix informed by trial-wise estimates from the RW model of each participant (see above). An additional analysis, using a more classical (non-computational model-based) design matrix, is presented in the Supplemental information (Figures S4–S6). Alongside the task regressors described below, rigorous de-noising was performed by the inclusion of the following regressors (all described above): the convolved end-tidal CO2 regressor plus temporal and dispersion derivatives, six motion regressors trajectories plus their first-order derivatives, physiological noise regressors (provided by the PhysIO toolbox) and ICA components identified as noise.

58The regressors of interest in the design matrix were as follows (compare Figure 3A and Figure S5):A ‘Cue’ regressor (80 repeats), with onsets and durations (2.5 s) determined by the presentations of visual cues and a magnitude of 1;A ‘Positive prediction’ regressor, with onsets given by the presentation of each corresponding visual cue (when no-resistance was predicted), durations of 0.5 s and magnitudes given by in Equation 5;A ‘Negative prediction’ regressor, with onsets given by the presentation of each corresponding visual cue (when resistance was predicted), durations of 0.5 s and magnitudes given by in Equation 6;A ‘No resistance’ stimulus regressor, with onset timings according to the first inspiration that occurred after the presentation of the visual cue, and durations as the remaining time of the potential resistance period (circle in Figure 3), with a magnitude of 1;A ‘Resistance’ stimulus regressor, with onset timings according to the initiation of the inspiratory resistance (identified from the downward inflection of the inspiratory pressure trace) after the presentation of the visual cue, and durations as the remaining time of the resistance period (circle in Figure 3), with a magnitude of 1;A ‘Positive prediction error’ regressor, with onsets given by the start of each corresponding no resistance period, durations of 0.5 s and magnitudes given by in Equation 7;A ‘Negative prediction error’ regressor, with onsets given by the start of each corresponding resistance period, durations of 0.5 s and magnitudes given by in Equation 8;A ‘Rating period’ noise regressor, with onsets and durations covering the period where participants were asked to rate the difficulty of the previous stimulus, and with a magnitude of 1.

59A ‘Cue’ regressor (80 repeats), with onsets and durations (2.5 s) determined by the presentations of visual cues and a magnitude of 1;

60A ‘Positive prediction’ regressor, with onsets given by the presentation of each corresponding visual cue (when no-resistance was predicted), durations of 0.5 s and magnitudes given by in Equation 5;

61A ‘Negative prediction’ regressor, with onsets given by the presentation of each corresponding visual cue (when resistance was predicted), durations of 0.5 s and magnitudes given by in Equation 6;

62A ‘No resistance’ stimulus regressor, with onset timings according to the first inspiration that occurred after the presentation of the visual cue, and durations as the remaining time of the potential resistance period (circle in Figure 3), with a magnitude of 1;

63A ‘Resistance’ stimulus regressor, with onset timings according to the initiation of the inspiratory resistance (identified from the downward inflection of the inspiratory pressure trace) after the presentation of the visual cue, and durations as the remaining time of the resistance period (circle in Figure 3), with a magnitude of 1;

64A ‘Positive prediction error’ regressor, with onsets given by the start of each corresponding no resistance period, durations of 0.5 s and magnitudes given by in Equation 7;

65A ‘Negative prediction error’ regressor, with onsets given by the start of each corresponding resistance period, durations of 0.5 s and magnitudes given by in Equation 8;

66A ‘Rating period’ noise regressor, with onsets and durations covering the period where participants were asked to rate the difficulty of the previous stimulus, and with a magnitude of 1.

67Regressors 1-8 were included in the design matrix after convolution with a standard HRF in SPM12, together with their temporal and dispersion derivatives. Contrasts of interest from this design examined brain activity associated with the average across positive and negative valence for both predictions and prediction errors, as well as the difference due to valence (i.e., positive versus negative) for both predictions and prediction errors.

68First, for the analysis of our entire field of view, contrasts of interest were assessed using random effects group-level GLM analyses based on the summary statistics approach in SPM12. The group-level GLM consisted of a factorial design with both a group mean and group difference regressor. The analyses used a significance level of p < 0.05 with family-wise error (FWE) correction at the cluster-level, with a cluster-defining threshold of p < 0.001. Second, for our region of interest (ROI) analysis, we used FSL’s non-parametric threshold-free cluster enhancement within a combined mask of the anterior insula and periaqueductal gray (PAG), as pre-specified in our analysis plan (section 7.5.5). This analysis employed a significance level of p < 0.05, with FWE correction across the joint mask. While the anterior insula and PAG have previously been shown to be involved in both conditioned anticipation and perception of inspiratory resistances as well as prediction errors and precision toward pain perception, our current analysis considers computational trial-by-trial estimates of interoceptive predictions and prediction errors for the first time. The mask of the anterior insula was taken from the Brainnetome atlas (bilateral ventral and dorsal anterior insula regions), and the PAG incorporated an anatomically-defined mask that has been used in previous fMRI publications.

69The different task modalities were then combined into a multi-modal analysis to assess both the relationships between and shared variance among measures. The data entered into this analysis consisted of:The scores from the four main affective questionnaires that were not used to pre-screen the participants (STAI-S, GAD-7, ASI-3 and CES-D);The four interoceptive questionnaires (BPQ, MAIA, PCS-B and PVQ-B);The four FDT measures (breathing sensitivity, decision bias c, metacognitive bias, metacognitive performance Mratio); andThe individual peak anterior insula activity associated with both positive and negative predictions, as well as positive and negative prediction errors. Activity was extracted from a 4mm sphere, centered on each participant’s maximal contrast estimate within a Brainnetome atlas mask of the anterior insula, using the first eigenvariate of the data.

70The scores from the four main affective questionnaires that were not used to pre-screen the participants (STAI-S, GAD-7, ASI-3 and CES-D);

71The four interoceptive questionnaires (BPQ, MAIA, PCS-B and PVQ-B);

72The four FDT measures (breathing sensitivity, decision bias c, metacognitive bias, metacognitive performance Mratio); and

73The individual peak anterior insula activity associated with both positive and negative predictions, as well as positive and negative prediction errors. Activity was extracted from a 4mm sphere, centered on each participant’s maximal contrast estimate within a Brainnetome atlas mask of the anterior insula, using the first eigenvariate of the data.

74A Pearson’s correlation matrix of all 16 included measures was calculated in order to visualize the relationships between all variables. The significance values of the correlation coefficients were taken as p < 0.05 (exploratory), and a false discovery rate (FDR) correction for multiple comparisons was applied (using the mafdr function in MATLAB). A supplementary non-parametric correlation matrix was additionally calculated using Spearman’s rho values, and these results are presented in Table S7.

75To assess the shared variance across measures and delineate which measures were most strongly associated with affective qualities, we entered all specified data into a principal component analysis (PCA), following normalization using z-scoring within each variable. PCA is an orthogonal linear transformation that transforms the data matrix (participants measures) into a new matrix, where the dimensions of the variance explained in the data are projected onto the new ‘principal components’ in descending order. Each principal component consists of a vector of coefficients or weights, corresponding to the contribution of each measure to each component. The PCA also transforms the original data matrix to map each row (participant) vector of onto a new vector of principal component scores, given by:where is the score for each participant within each component. The number of significant components were then determined by comparing the variance explained of each component to a null distribution, created by randomly shuffling (n = 1000) the measures from each variable across participants. Components were considered significant if the variance explained was above the 95% confidence interval of the corresponding component’s null distribution.

76To assess the relationship between each of the significant components and anxiety, the component scores for low and moderate anxiety were compared using either independent t tests or Wilcoxon rank sum tests (following Anderson-Darling tests for normality). The significance values of the group differences in component scores were taken as p < 0.05 (exploratory), and a false discovery rate (FDR) correction for multiple comparisons (number of significant components) was applied.

77An independent code review was performed on all data analysis procedures, and the analysis code is available on GitLab (https://gitlab.ethz.ch/tnu/code/harrison_breathing_anxiety_code; https://doi.org/10.5281/zenodo.5523258).