Seventeen international wheelchair court sport players split into two groups according to SCI level—9 TETRA (1 woman, 8 men) and 8 PARA (1 woman, 7 men)—volunteered for this study (Table 1). All athletes used a wheelchair for daily ambulation. Most athletes self-reported to have a sensory- and motor-complete SCI. 2 TETRA had a sensory-incomplete and 1 a sensory- and motor-incomplete injury. 3 PARA had a motor-incomplete and 1 a sensory-incomplete injury. The study was approved by the Ethical Committee of Loughborough University and conducted in accordance with the Declaration of Helsinki. All participants provided written informed consent prior to the data collection and further completed separate health, training and impairment questionnaires. Note that the data were obtained as part of a previously published study by Leicht et al. (2014). However, we removed the data of one PARA and one TETRA due to incomplete data that would have been required for the model fitting approach chosen in the present study.

For a detailed description of the experimental design, see Leicht et al. (2014). In brief, participants performed a continuous graded exercise test to exhaustion in their sports wheelchair on a motorized treadmill (HP Cosmos, Traunstein, Germany) at a constant 1.0% gradient. Wheelchairs were secured to a safety rail located on the side of the treadmill that kept the wheelchair centered but allowed free back and forward movement within the treadmill dimension. The starting speed of the test was set between 1.2 and 2.0 m⋅s^–1 with 0.2–0.4 m⋅s^–1 increments every 3 min, which were individually adjusted based on the differences in impairment and fitness levels. The test was terminated when participants were unable to maintain the speed of the treadmill, touching the “backstop” of the safety rail for the third time. Respiratory data (V̇O₂, V̇CO₂, V̇E) were recorded continuously using an online gas analysis system in breath-by-breath mode (Meta-Lyzer 3B, Cortex Biophysik GmbH, Leipzig, Germany).

The data of TETRA and PARA was previously analyzed in the study of Leicht et al. (2014), where blood lactate and respiratory thresholds between TETRA and PARA were determined by expert opinion and then compared. This is a subjective approach, which we extend on in the current study by determining thresholds using objective, automated procedures.

In the current study, a breakpoint model (i.e., two regression lines, Equation 1) and a continuous curvilinear (no-breakpoint) model (i.e., an exponential curve, Equation 2) were fitted to the V̇O₂-V̇CO₂ (identification of GET) and the V̇CO₂-V̇E (identification of RCP) data. For the V̇E/V̇O₂-time (identification of GET) and V̇E/V̇CO₂-time (identification of RCP) data, the continuous curvilinear (no-breakpoint) model was a second order polynomial (Equation 3). We considered the identification of the GET and RCP by the ventilatory equivalent method valid if the breakpoint in the V̇E/V̇O₂-time occurred before the breakpoint in the V̇E/V̇CO₂-time data (Powers et al., 1984; Gaskill et al., 2001). Matlab (R2021a; Mathworks Inc., Natick, MA) was used to perform the model fitting procedures. The function slmengine of the Shape Language Modeling (SLM) toolkit (D’Errico, 2021) was used to fit two regression lines to the data, including the function fmincon with the interior-point algorithm to find the best fitting model. A custom-made function was written to fit the exponential curve to the data, and the function fminsearch using the Nelder-Mead approach (Nelder and Mead, 1965) was used to determine the best fit with a maximum of 40,000 evaluations and 1,500 iterations. The function polyfit was used to fit the second-order polynomial to the data.

where y is the variable of interest, a the y-axis offset, b to f the slope coefficients, and k represents the intersection point of the first and the second regression lines of the piecewise function. Prior to fitting the breakpoint and no-breakpoint models to the breath-by-breath data an outlier detection was performed. Data points were considered outliers and removed if they exceeded 3 standard deviations from a third-order polynomial that was fitted to the data of each individual participant (removal of mean percentage ± SD of data points in V̇O₂-V̇CO₂ data: 1.4 ± 0.6%, V̇E/V̇O₂-time data: 1.2 ± 0.4%, V̇E/V̇CO₂-time data: 1.9 ± 0.8%, V̇CO₂-V̇E data: 1.4 ± 0.6%).

The participant characteristics provided in Table 1 were compared by employing independent samples T-tests. An α level of 0.05 was used to indicate statistical significance. The adjusted R² (Equation 4) was used to compare the fit of the breakpoint model (Equation 1) with the fit of the continuous curvilinear model (Equations 2 and 3) for each individual participant.

where p is the number of coefficients in the model and N the number of data points. The adjusted R² was used to account for the different number of coefficients, with four coefficients for the two regression lines, and three coefficients for the exponential curve and second-order polynomial. The overlapping blocked bootstrap technique with 10-s intervals was chosen to account for dependency in the data (Davison and Hinkley, 1997; Bühlmann, 2002), and used to attain one hundred random samples of the underlying respiratory data for each participant. An adjusted R² value was then obtained for the fit of each model to each of the resampled V̇O₂-V̇CO₂, V̇E/V̇O₂-time, V̇CO₂-V̇E, and V̇E/V̇CO₂-time data. This resulted in 100 adjusted R² values for each model for every individual participant. The difference between the adjusted R² values, and thus fit of the breakpoint model compared with the non-breakpoint model, for each individual and each of the approaches (i.e., two approaches for the GET and two approaches for the RCP) was assessed by comparing notched boxplots. The notches are 95% CIs that are constructed around the median (see Supplementary Figure 1). No overlap in the notches indicates 95% confidence of a difference between medians (Chambers et al., 1983). Our rationale was that a better fit of the continuous curvilinear model (no-breakpoint model) as compared with the two regression lines (breakpoint model), would challenge the presence of a distinct breakpoint. Notched boxplots were also used to assess whether the adjusted R² differed between breakpoint and no-breakpoint models with the data of TETRA and PARA pooled, as well as to investigate whether there was a difference in model fit within and between those two groups.

Whilst the employed breakpoint models to determine the GET and RCP do not improve the goodness of fit to the data when compared with no-breakpoint models, breakpoints may occur after all. These breakpoints may be too subtle for the two regression lines to result in a better fit compared to the continuous no-breakpoint models. For example, breakpoints may be identified by only a handful of data points that do not follow a general curvilinear trend, too few to significantly alter any regression metrics. Automatic procedures may require improved fitting algorithms that more subtly address if functions, which best explain the data, are differentiable or not (the mathematical description of a breakpoint). Searching for more suitable differentiable continuous functions is only one option. Applying spectral density analysis procedures may be another. Machine learning, pattern recognition and artificial intelligence approaches may be a third. The success for any of such approaches depends heavily on the signal-variance ratio of the data, particularly if true breakpoints are subtle (as may be the case in this study).

In favor of the argument that breakpoints exist, GET and RCP methods are relatively widely used to identify intensity zones and, in some studies, high to very high inter- and intra-rater reliability has been reported (Gaskill et al., 2001; Amann et al., 2004; Pallarés et al., 2016; Au et al., 2018; Kouwijzer et al., 2019). There are few publications that use automatic determination of thresholds (Zignoli et al., 2019, 2020), with the majority opting for a fully subjective approach of visual determination (Aunola and Rusko, 1984; Gladden et al., 1985; Bhambhani et al., 1993; Holland et al., 1994; Kouwijzer et al., 2019). The inherent inaccuracy of this approach is reflected by a statement by Hopker et al. (2011), who noted that the GET and RCP methods “very crudely identify a time-point somewhere near to when the non-linearity in respiratory responses becomes more pronounced.” The investigated data set of the current study confirms this deviation from non-linearity, observed most strikingly in the ventilatory equivalent plots (Figure 3). However, our data also cast further doubt on the practical use of GET and RCP. Indeed, we showed that in the majority of cases, automatic threshold determination using the ventilatory equivalent method resulted in the RCP to be determined at lower exercise intensities than the GET, which is in contradiction to their definition, in which the RCP occurs at higher exercise intensity than the GET (Gaskill et al., 2001; Meyer et al., 2005). Given that the relationship in the V̇O₂-V̇CO₂ data is rather linear, it appears that it is mainly the rapid increase in V̇E that leads to the disproportionate rise in both the V̇E/V̇O₂-time and the V̇E/V̇O₂-time data. As such, it is not surprising that also for the remainder of cases, GET and RCP were determined in such proximity that their use in exercise prescription would be limited, i.e., it would not be possible to create a large enough moderate intensity zone within their boundaries. It is possible that subjective approaches are prone to subconscious bias (Hopker et al., 2011), leading to threshold identification that is “in line with underlying theory,” and/or being sufficiently far apart to allow creation of an intensity zone. Threshold identification based on traditional, subjective approaches may hence be limited to the approximate identification of a change in data patterns at best, or suffer from such bias at worst.

Our findings support the general notion that respiratory changes in response to a gradual increase in exercise intensity are curvilinear. We here discuss this line of thought, especially as this may appear to be in contradiction to published research. Indeed, in contrast to the current study, distinct breakpoints in V̇E or related variables have been previously presented convincingly in figure format (Beaver et al., 1986; Wasserman, 1987). It is worth noting, however, that similar to the point made regarding bias in threshold identification, researchers have a tendency to present their “clearest dataset,” i.e., often the dataset with the least variability, and/or the one that matches the underlying theory (van der Steen et al., 2019). An example to support this notion in the present study are the data of PARA₈ in Figure 1, which display one of the most distinct breakpoints in our study sample, even though the box plots show no significantly better fit of the breakpoint compared to the no-breakpoint model. To avoid such a selection bias, we here present all individual datasets (see Supplementary Figures 2–4) in addition to the data of example participants provided in Figures 1–3. These data indeed support the notion that V̇E is highly variable and increases in a curvilinear rather than a breakpoint fashion in response to an increase in exercise intensity. This pattern is likely reflected in V̇CO₂ data, as it has previously been suggested that increases in V̇CO₂ relate to increases in V̇E (Hopker et al., 2011). As V̇O₂ and increases in submaximal exercise intensity are linearly related (Arts and Kuipers, 1994), it would hence not be expected that any relationships between exercise intensity, V̇O₂, and/or V̇CO₂ (derived) data should be characterized by any distinct breakpoints. Indeed, there is evidence that the transition from aerobic to anaerobic contributions with increasing exercise intensity occurs gradually, rather than at a specific intensity (Helal et al., 1987; Chamari and Padulo, 2015). Further, even though some studies report high intra-rater reliability as outlined above, there is also a considerable number of studies that report lower reliability (Shimizu et al., 1991; Bhambhani et al., 1993; Holland et al., 1994). This implies that the occurrence of a distinct point at which thresholds occur may either not be clear cut or, indeed, be non-existent. Altogether, this raises doubts on the concept of identifying breakpoints in respiratory data to distinguish low, moderate and high intensity exercise.

In addition, the absence of distinct breakpoints may be linked to the specific population tested. Breath-by-breath data are relatively variable, and when identified within a limited range of respiratory responses, breakpoints may not be visible. It is worth noting that respiratory responses are lower during upper-body compared to lower-body exercise, even in able-bodied participants (Nagle et al., 1984). Importantly, these responses are further reduced in people with more severe impairments such as a cervical SCI (Coutts et al., 1983), and accordingly, the ratio of the variability to the range of all respiratory variables was larger in TETRA in the present study (Supplementary Figures 2–4). This may explain to some extent that the goodness of fit of both, the breakpoint and no-breakpoint models, was lower in TETRA when compared to PARA, even though the fit between the breakpoint and no-breakpoint models was not significantly different between subgroups. Overall, we question how distinct these breakpoints are in the respiratory data of wheelchair athletes.

In the current study, we selected two very commonly used analytical methods to identify the GET and the RCP. There are other methods to identify the GET and RCP (Gaskill et al., 2001), and it remains speculative whether breakpoints are also not that distinct when applying these other methods. In addition to respiratory methods for identifying breakpoints that distinguish low, moderate and high intensity, a range of methods that set out to do the same are based on blood lactate concentration data (Faude et al., 2009). Using the concept of fitting breakpoint vs no-breakpoint curvilinear models to such data was beyond the scope of this study but is likely subject to the same inherent problems of absence of distinct breakpoints. Further, while the exponential or polynomial models have a strong grounding in existing literature (Fairshter et al., 1987; Dennis et al., 1992; Su et al., 2007), it remains to be investigated if there are curvilinear models with an even better fit to the data. A superior model fit of the curvilinear models may further strengthen our argument of the absence of distinct breakpoints in the respiratory data.

We focused on the respiratory responses of wheelchair court sport players with SCI and the interpretation of our results is therefore specific to this population. In intermittent sports, a greater endurance capacity is linked to better repeated sprint ability as well as faster recovery between games (McMahon and Wenger, 1998; Baumgart and Sandbakk, 2016), pointing out the relevance to determine endurance capacity markers such as GET or RCP in this cohort. Our results apply to highly trained athletes, who are characterized by a larger range of respiratory responses or markers of endurance performance than are seen in untrained wheelchair users (Janssen et al., 2002). Indeed, for untrained individuals, threshold determination may be further complicated by the lower range of respiratory variables. While we did not record information on indicators of reduced respiratory function, such as autonomic dysfunction and early onset sleep disordered breathing, the smaller range of the respiratory data of TETRA compared to PARA is in support of this.

Given that the present study was not designed to distinguish effects of disability and exercise mode, replicating it by including able-bodied upper-body trained participants would help to isolate the effect of the upper-body exercise mode from the disability. Finally, future investigations are needed to look into whether our findings of subtle or absent breakpoints in the respiratory responses during upper-body exercise extend to lower- and whole-body exercise, which elicit higher and less variable respiratory responses (Sawka et al., 1982; Castro et al., 2011).

We here show that breakpoints in the respiratory data of wheelchair athletes obtained during upper-body exercise testing are either very subtle, or they do not exist. This information is of value to both athletes and coaches, as it indicates that the use of the methods employed in the current study should be discouraged to identify low-, moderate- and high-intensity zones. As this implies that intensity zone boundaries cannot be accurately determined based on ventilatory data from graded exercise tests to exhaustion, such intensity zones need to be verified using alternative methods. This may include performing steady state exercise tests to establish the maximum lactate steady state or the highest exercise intensity that does not raise blood lactate above resting concentrations (“the lactate threshold”) (Beneke et al., 2011). Alternatively, critical power may be determined using a series of supramaximal sprints (Jones et al., 2010). If verified as such, thresholds may remain a valuable tool to guide training. If applicable, training zones should be assessed separately for different exercise modalities.