$$p(y_{(k)}|v_{(k)}, \zeta) = \left( \frac{v_{(k)}^{\zeta}}{v_{(k)}^{\zeta} + (1 - v_{(k)})^{\zeta}} \right)^{y_{(k)}} \left( \frac{(1 - v_{(k)})^{\zeta}}{v_{(k)}^{\zeta} + (1 - v_{(k)})^{\zeta}} \right)^{1-y_{(k)}} \quad (\text{Equation 2})$$

Here, $y_{(k)}$ represents the expressed decision of a subject given the cue (contingency pairs) on trial k. The parameter $\zeta$ captures how deterministically y is associated with v. The higher $\zeta$, 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 (Mathys et al., 2014), given the trajectories of trial-wise cues and responses (see Figure 3).

In 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 (Iglesias et al., 2013). 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).

### Comparison of fitted model parameters

Group differences in model parameter estimates of learning rate ($\alpha$) and inverse decision temperature ($\zeta$), 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.

### Model validation

Following random effects Bayesian model selection (BMS (Rigoux et al., 2014; Stephan et al., 2009)), 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) (Daw, 2011) 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).

To 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.

### Computationally informed regressors

The 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:

$$v_{(k)}^{\text{stim}} \stackrel{\text{def}}{=} \begin{cases} v_{(k)}, & \text{if cue type} = 1 \\ 1 - v_{(k)}, & \text{if cue type} = 2 \end{cases} \quad (\text{Equation 3})$$

$$\delta_{(k)}^{\text{stim}} \stackrel{\text{def}}{=} \begin{cases} \delta_{(k)}, & \text{if cue type} = 1 \\ -\delta_{(k)}, & \text{if cue type} = 2 \end{cases} \quad (\text{Equation 4})$$

Here, $v_{(k)}^{\text{stim}}$ and $\delta_{(k)}^{\text{stim}}$ now represent the prediction and prediction error values in stimulus space, with $v_{(k)}^{\text{stim}} = 1$ representing maximal predictions of no resistance and $v_{(k)}^{\text{stim}} = 0$ maximal predictions of resistance. Similarly, $\delta_{(k)}^{\text{stim}} = 1$ represents maximal prediction errors of no resistance and $\delta_{(k)}^{\text{stim}} = -1$ maximal prediction errors of resistance (see Figure S4 for details).

Second, 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
$$p(y_{(k)}|v_{(k)}, \zeta) = \left( \frac{v_{(k)}^{\zeta}}{v_{(k)}^{\zeta} + (1 - v_{(k)})^{\zeta}} \right)^{y_{(k)}} \left( \frac{(1 - v_{(k)})^{\zeta}}{v_{(k)}^{\zeta} + (1 - v_{(k)})^{\zeta}} \right)^{1-y_{(k)}} \quad (\text{Equation 2})$$

Here, $y_{(k)}$ represents the expressed decision of a subject given the cue (contingency pairs) on trial k. The parameter $\zeta$ captures how deterministically y is associated with v. The higher $\zeta$, 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 (Mathys et al., 2014), given the trajectories of trial-wise cues and responses (see Figure 3).

In 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 (Iglesias et al., 2013). 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).

### Comparison of fitted model parameters

Group differences in model parameter estimates of learning rate ($\alpha$) and inverse decision temperature ($\zeta$), 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.

### Model validation

Following random effects Bayesian model selection (BMS (Rigoux et al., 2014; Stephan et al., 2009)), 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) (Daw, 2011) 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).

To 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.

### Computationally informed regressors

The 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:

$$v_{(k)}^{\text{stim}} \stackrel{\text{def}}{=} \begin{cases} v_{(k)}, & \text{if cue type} = 1 \\ 1 - v_{(k)}, & \text{if cue type} = 2 \end{cases} \quad (\text{Equation 3})$$

$$\delta_{(k)}^{\text{stim}} \stackrel{\text{def}}{=} \begin{cases} \delta_{(k)}, & \text{if cue type} = 1 \\ -\delta_{(k)}, & \text{if cue type} = 2 \end{cases} \quad (\text{Equation 4})$$

Here, $v_{(k)}^{\text{stim}}$ and $\delta_{(k)}^{\text{stim}}$ now represent the prediction and prediction error values in stimulus space, with $v_{(k)}^{\text{stim}} = 1$ representing maximal predictions of no resistance and $v_{(k)}^{\text{stim}} = 0$ maximal predictions of resistance. Similarly, $\delta_{(k)}^{\text{stim}} = 1$ represents maximal prediction errors of no resistance and $\delta_{(k)}^{\text{stim}} = -1$ maximal prediction errors of resistance (see Figure S4 for details).

Second, 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 (Paulus and Stein, 2010; Paulus et al., 2019) 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):

$$ \text{If } v_{(k)}^{\text{stim}} > 0.5 v_{(k)}^{\text{pos}} \underline{\underline{\text{def}}} v_{(k)}^{\text{stim}} - 0.5 \qquad (\text{Equation 5}) $$

$$ \text{If } v_{(k)}^{\text{stim}} < 0.5 v_{(k)}^{\text{neg}} \underline{\underline{\text{def}}} 0.5 - v_{(k)}^{\text{stim}} \qquad (\text{Equation 6}) $$

Here, both $v_{(k)}^{\text{pos}}$ and $v_{(k)}^{\text{neg}}$ exist between 0 and 0.5, with values closer to zero indicating less certain predictions.

Like 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):

$$ \text{If } \delta_{(k)} > 0 \delta_{(k)}^{\text{pos}} \underline{\underline{\text{def}}} \delta_{(k)} \qquad (\text{Equation 7}) $$

$$ \text{If } \delta_{(k)} < 0 \delta_{(k)}^{\text{neg}} \underline{\underline{\text{def}}} - \delta_{(k)} \qquad (\text{Equation 8}) $$

Here, both $\delta_{(k)}^{\text{pos}}$ and $\delta_{(k)}^{\text{neg}}$ 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 (Iglesias et al., 2013)).

## Physiological data processing

Physiological data were recorded at a sampling rate of 1000 Hz, and included heart rate, chest distension, pressure of expired carbon dioxide ($P_{ETCO_2}$) and oxygen ($P_{ETO_2}$), and pressure at the mouth (for equipment details see (Rieger et al., 2020)). In addition to the task, small boluses of a $CO_2$ gas mixture (20% $CO_2$; 21% $O_2$; balance $N_2$) were administered during some rest periods, allowing for de-correlation of any changes in $P_{ETCO_2}$ from task-related neural activity, as previously described (Faull and Pattinson, 2017; Faull et al., 2015, 2016, 2018).

Physiological noise regressors were prepared for inclusion into single-subject general linear models (GLMs, described below). Linear interpolation between $P_{ETCO_2}$ peaks was used to form an additional $CO_2$ 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 (Chang and Glover, 2009). Temporal and dispersion derivatives of this $CO_2$ noise regressor were also included. An additional three cardiac- and four respiratory-related waveforms (plus one interaction term) were created using PhysIO (Kasper et al., 2017). Four respiratory volume per unit time (RVT) regressors (delays: -5, 0, 5, 10) were created using the Hilbert-transform estimator in PhysIO (Harrison et al., 2021b), and convolution with a respiratory response function (Kasper et al., 2017).

## Magnetic resonance imaging

MRI 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.

## MRI preprocessing

MRI 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) (Jenkinson et al., 2012) and SPM12 (Statistical Parametric Mapping software, London, UK) (Friston et al., 2011) 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 (Jenkinson and Smith, 2001)), removal of non-brain structures (BET (Smith, 2002)), and high-pass temporal filtering (Gaussian-weighted least-squares straight line fitting; 100 s cut-off period)(Woolrich et al., 2001). Independent component analysis (ICA) was used to identify noise due to motion, scanner and cerebrospinal fluid artifacts (Griffanti et al., 2017), 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
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 (Paulus and Stein, 2010; Paulus et al., 2019) 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):

$$ \text{If } v_{(k)}^{\text{stim}} > 0.5 v_{(k)}^{\text{pos}} \underline{\underline{\text{def}}} v_{(k)}^{\text{stim}} - 0.5 \quad (\text{Equation 5}) $$

$$ \text{If } v_{(k)}^{\text{stim}} < 0.5 v_{(k)}^{\text{neg}} \underline{\underline{\text{def}}} 0.5 - v_{(k)}^{\text{stim}} \quad (\text{Equation 6}) $$

Here, both $v_{(k)}^{\text{pos}}$ and $v_{(k)}^{\text{neg}}$ exist between 0 and 0.5, with values closer to zero indicating less certain predictions.

Like 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):

$$ \text{If } \delta_{(k)} > 0 \delta_{(k)}^{\text{pos}} \underline{\underline{\text{def}}} \delta_{(k)} \quad (\text{Equation 7}) $$

$$ \text{If } \delta_{(k)} < 0 \delta_{(k)}^{\text{neg}} \underline{\underline{\text{def}}} - \delta_{(k)} \quad (\text{Equation 8}) $$

Here, both $\delta_{(k)}^{\text{pos}}$ and $\delta_{(k)}^{\text{neg}}$ 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 (Iglesias et al., 2013)).

## Physiological data processing

Physiological data were recorded at a sampling rate of 1000 Hz, and included heart rate, chest distension, pressure of expired carbon dioxide ($P_{\text{ETCO}_2}$) and oxygen ($P_{\text{ETO}_2}$), and pressure at the mouth (for equipment details see (Rieger et al., 2020)). In addition to the task, small boluses of a $CO_2$ gas mixture (20% $CO_2$; 21% $O_2$; balance $N_2$) were administered during some rest periods, allowing for de-correlation of any changes in $P_{\text{ETCO}_2}$ from task-related neural activity, as previously described (Faull and Pattinson, 2017; Faull et al., 2015, 2016, 2018).

Physiological noise regressors were prepared for inclusion into single-subject general linear models (GLMs, described below). Linear interpolation between $P_{\text{ETCO}_2}$ peaks was used to form an additional $CO_2$ 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 (Chang and Glover, 2009). Temporal and dispersion derivatives of this $CO_2$ noise regressor were also included. An additional three cardiac- and four respiratory-related waveforms (plus one interaction term) were created using PhysIO (Kasper et al., 2017). Four respiratory volume per unit time (RVT) regressors (delays: -5, 0, 5, 10) were created using the Hilbert-transform estimator in PhysIO (Harrison et al., 2021b), and convolution with a respiratory response function (Kasper et al., 2017).

## Magnetic resonance imaging

MRI 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.

## MRI preprocessing

MRI 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) (Jenkinson et al., 2012) and SPM12 (Statistical Parametric Mapping software, London, UK) (Friston et al., 2011) 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 (Jenkinson and Smith, 2001)), removal of non-brain structures (BET (Smith, 2002)), and high-pass temporal filtering (Gaussian-weighted least-squares straight line fitting; 100 s cut-off period)(Woolrich et al., 2001). Independent component analysis (ICA) was used to identify noise due to motion, scanner and cerebrospinal fluid artifacts (Griffanti et al., 2017), 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 (Jenkinson et al., 2002); 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 (Andersson et al., 2007) 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.

## Single-subject general linear model

Single-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 CO₂ 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.

The regressors of interest in the design matrix were as follows (compare Figure 3A and Figure S5):

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

2) 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 $v_{(k)}^{\text{pos}}$ in Equation 5;

3) 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 $v_{(k)}^{\text{neg}}$ in Equation 6;

4) 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;

5) 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;

6) 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 $\delta_{(k)}^{\text{pos}}$ in Equation 7;

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 $\delta_{(k)}^{\text{neg}}$ in Equation 8;

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.

Regressors 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.

## Group fMRI analysis

First, 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 (Smith and Nichols, 2009) within a combined mask of the anterior insula and periaque ductal 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 (Berner et al., 2018; Faull and Pattinson, 2017; Faull et al., 2016, 2018; Paulus et al., 2012; Walter et al., 2020) as well as prediction errors (Roy et al., 2014) and precision (Grahl et al., 2018) 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 (Fan et al., 2016) (bilateral ventral and dorsal anterior insula regions), and the PAG incorporated an anatomically-defined mask that has been used in previous fMRI publications (Faull and Pattinson, 2017; Faull et al., 2016).

## Multi-modal analysis

### Multi-modal data

The 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:
scan to the whole-brain EPI image (Jenkinson et al., 2002); 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 (Andersson et al., 2007) 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.

## Single-subject general linear model

Single-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 CO₂ 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.

The regressors of interest in the design matrix were as follows (compare Figure 3A and Figure S5):

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

2) 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 $v_{(k)}^{\text{pos}}$ in Equation 5;

3) 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 $v_{(k)}^{\text{neg}}$ in Equation 6;

4) 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;

5) 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;

6) 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 $\delta_{(k)}^{\text{pos}}$ in Equation 7;

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 $\delta_{(k)}^{\text{neg}}$ in Equation 8;

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.

Regressors 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.

## Group fMRI analysis

First, 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 (Smith and Nichols, 2009) within a combined mask of the anterior insula and periaque ductal 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 (Berner et al., 2018; Faull and Pattinson, 2017; Faull et al., 2016, 2018; Paulus et al., 2012; Walter et al., 2020) as well as prediction errors (Roy et al., 2014) and precision (Grahl et al., 2018) 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 (Fan et al., 2016) (bilateral ventral and dorsal anterior insula regions), and the PAG incorporated an anatomically-defined mask that has been used in previous fMRI publications (Faull and Pattinson, 2017; Faull et al., 2016).

## Multi-modal analysis

### Multi-modal data

The 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:
1) The scores from the four main affective questionnaires that were not used to pre-screen the participants (STAI-S (Spielberger et al., 1970), GAD-7 (Spitzer et al., 2006), ASI-3 (Taylor et al., 2007) and CES-D (Radloff, 1977));

2) The four interoceptive questionnaires (BPQ (Porges, 1995), MAIA (Mehling et al., 2012), PCS-B (Sullivan et al., 1995) and PVQ-B (McCracken, 1997));

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

4) The 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 (Fan et al., 2016), using the first eigenvariate of the data.

## Multi-modal correlations and shared variance

A 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.

To 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 $n \times m$ data matrix $X$ (participants × measures) into a new matrix $P$, 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 $w$, corresponding to the contribution of each measure $m$ to each component. The PCA also transforms the original $n \times m$ data matrix $X$ to map each row (participant) vector $x_i$ of $X$ onto a new vector of principal component scores $t_i$, given by:

$$t_{k(i)} = x_i \times w_k \quad \text{for } i = 1, \dots, n; k = 1, \dots, m \qquad (\text{Equation 9})$$

where $t_{k(i)}$ is the score for each participant $i$ within each component $k$. 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.

To 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.

An 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).
1) The scores from the four main affective questionnaires that were not used to pre-screen the participants (STAI-S (Spielberger et al., 1970), GAD-7 (Spitzer et al., 2006), ASI-3 (Taylor et al., 2007) and CES-D (Radloff, 1977));

2) The four interoceptive questionnaires (BPQ (Porges, 1995), MAIA (Mehling et al., 2012), PCS-B (Sullivan et al., 1995) and PVQ-B (McCracken, 1997));

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

4) The 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 (Fan et al., 2016), using the first eigenvariate of the data.

## Multi-modal correlations and shared variance

A 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.

To 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 $n \times m$ data matrix $\mathbf{X}$ (participants × measures) into a new matrix $\mathbf{P}$, 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 $\mathbf{w}$, corresponding to the contribution of each measure $m$ to each component. The PCA also transforms the original $n \times m$ data matrix $\mathbf{X}$ to map each row (participant) vector $\mathbf{x}_i$ of $\mathbf{X}$ onto a new vector of principal component scores $\mathbf{t}_i$, given by:

$$t_{k(i)} = \mathbf{x}_i \times \mathbf{w}_k \quad \text{for } i = 1, \dots, n; k = 1, \dots, m \qquad (\text{Equation 9})$$

where $t_{k(i)}$ is the score for each participant $i$ within each component $k$. 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.

To 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.

An 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).