Estimation of Work of Breathing from Respiratory Muscle Activity In Spontaneous Ventilation: A Pilot Study

Abstract

Work of breathing (WOB) offers information that may be relevant to determine the patient’s status under spontaneous mechanical ventilation in Intensive Care Unit (ICU). Nowadays, the most reliable technique to measure WOB is based on the use of invasive catheters, but the use of qualitative observations such as the level of dyspnea is preferred as a possible indicator of WOB level. In this pilot study, the activity of three respiratory muscles were recorded on healthy subjects through surface electromyography while they were under non-invasive mechanical ventilation, using restrictive and obstructive maneuvers to obtain different WOB levels. The respiratory pattern between restrictive and obstructive maneuvers was classified with the Nearest Neighbor Algorithm with a 91% accuracy and a neural network model helped classify the samples into three WOB levels with a 89% accuracy, Low: [0.3–0.8) J/L, Medium: [0.8–1.3] J/L and Elevated: (1.3–1.8] J/L, demonstrating the relationship between the respiratory muscle activity and WOB. This technique is a promising tool for the healthcare staff in the decision-making process when selecting the best ventilation settings to maintain a low WOB. This study identified a model to estimate the WOB in different ventilatory patterns, being an alternative to invasive conventional techniques.

1. Introduction

The work of breathing (WOB) has been defined as the energy required to accomplish the body´s ventilatory demand, which in spontaneous breathing depends on the effort performed by respiratory muscles [1]. Monitoring WOB levels provides several advantages for healthcare staff in the decision-making process, making it possible to adjust the mechanical ventilator’s settings, and the pharmacological treatment; furthermore, it helps understand the mechanism of the pathologies and assess the patient’s evolution post-extubation [1,2,3].

Some authors have identified WOB values for different patients, both with and without mechanical ventilation. WOB values in the range of 0.3 to 0.8 J/L are considered indicators of low respiratory demand in patients under spontaneous mechanical ventilation or low WOB [4,5,6]. Sharp et al. [6] found that 96% of patients with a WOB lower than 0.8 J/L had a successful weaning process. On the other hand, values in the range of 0.8–1.3 J/L are considered a population with high respiratory demand, such as people with airway and pulmonary parenchymal abnormalities, intubated patients, people with dynamic hyperinflation, airflow obstruction (Chronic Obstructive Pulmonary Disease or COPD patients) or obese people [4,7,8]. In addition, WOB values in the range of 1.3–1.8 J/L are considered elevated, making it impossible to maintain an adequate spontaneous ventilation for a long period of time [5,8].

According to the state of the art, the use of an esophageal balloon catheter is required to accurately measure the work of breathing. However, this technique is invasive and uncomfortable for the patient, which restricts its use in clinical practice [9,10].

Several clinical strategies to detect an increase in work of breathing include distress or dyspnea, vital signs, physical signs of increased breathing effort, breathing rate, tidal volume, and frequency/tidal volume ratio [3,5,11]. However, these strategies are not reliable in all cases, which hampers its generalization [3].

Dyspnea is defined as “breathing discomfort”, which occurs when there is an imbalance between the load imposed upon the respiratory muscles and the capacity of these to overcome it. For that reason, many authors have successfully used surface electromyography (sEMG) signals of respiratory muscles to quantify the level of dyspnea [12,13,14]. Therefore, a relationship between respiratory muscle activity and work of breathing can be inferred.

Some authors have proposed indexes from sEMG signals of respiratory muscles, in the following domains: time, frequency, and time-frequency [12,15,16,17]. However, the relationship between respiratory muscle activity and work of breathing has not been proven yet, because those authors had no WOB measurements. Furthermore, considering that multiple respiratory muscles are involved, patients may have different ventilation patterns to achieve the same work, and to our knowledge this has not been studied.

Some researchers have studied muscle activity patterns from sEMG signals and used different machine learning algorithms to recognize activation patterns in several muscles, mainly of the limbs. Chowdhury et al. [18] presented a review describing the use of sEMG, showing selected features, used muscles, classification algorithms, objective, and the accuracy of the studies. A common objective was the classification of patients with neuromuscular diseases and muscular fatigue using sEMG signals of the Biceps brachii muscle. Consequently, classifiers like the Simple Fisher Linear Discriminant Analysis (LDA), statistical K-Nearest Neighbor (KNN), Support Vector Machine (SVM), Neural Networks (NN), among others, were used.

In mechanical ventilation, in order to achieve an optimal configuration of the ventilator, some authors have used Artificial Neural Networks to monitor parameters like total positive end-expiratory pressure (PEEP) [19] and power of breathing or work of breathing per minute [11,20] from ventilatory signals.

In previous studies, our team used unsupervised learning techniques to find indexes that were related to the patient’s condition or identify muscle respiratory patterns [15,17,21]. Based on the hierarchical cluster, an unsupervised learning technique, Hernández et al. [21] characterized different muscular responses attributed to changes in the mechanical ventilator settings. Using similar techniques, Salazar et al. [15] identified muscle weakness in ventilated patients, and in a complementary study, patients who presented failures in the weaning process were separated from patients with a successful weaning [17].

This paper presents a strategy to estimate work of breathing exerted by the subject during the ventilation process using indexes from sEMG signals based on machine learning techniques to identify and discriminate different respiratory patterns in spontaneous ventilation.

2. Materials and Methods

2.1. Experimental Design

The study was conducted according to the Helsinki Declaration with subsequent revisions. The experimental design and written consent were approved by the Ethics Committee of the University of Antioquia, approval report 15-59-664.

Two experiments were designed to achieve different levels of work of breathing. Experiment one (E1) involved restrictive maneuvers to reduce the respiratory system compliance and experiment two (E2) involved obstructive maneuvers to increase airway resistance.

In both experiments, the participants were connected to a mechanical ventilator Hamilton G5 (Hamilton Medical, Bonaduz, Switzerland) through a non-invasive oronasal mask. The ventilator was set in the spontaneous mode, with a PEEP of 0 cmH2O, pressure support of 0 cmH2O, and FiO2 of 21%.

In E1, a cuirass [22] was placed on the chest and abdomen of the volunteers changing the pressure in the cuirass with an external controlled pump, which allowed the reduction of the respiratory system compliance. Four samples were recorded for each subject. The first sample was recorded without the cuirass, followed by three samples sequentially with pressures of 0, 10, and 20 cmH2O inside the cuirass.

In E2, two different external airflow resistors (Model 7100R; Hans Rudolph, Inc., Kansas City, MO) were used between the oronasal mask and the flow sensor. Three samples were recorded in each subject sequentially with resistive loads of 0, 5, and 20 cmH2O/L/s (see Figure 1).

Each sample consisted of a 3-minute recording, where the first two minutes were considered an adaption period after a change of stimulus because there is high variability among subjects during the transient. Therefore, only the last minute associated with stable state was analyzed in this study. As shown in Figure 1, the samples obtained per volunteer are recorded sequentially, e.g., in E2, the volunteer began with an external resistor of 0 cmH2O for 3 minutes, and subsequently, the external resistor was changed to 5 cmH2O; the 2 minutes following the change in resistance is considered the adaption period or time needed to stabilize the ventilatory process.

2.2. Subjects

The inclusion criterion was: Adult male subjects. The exclusion criteria were: subjects with a body mass index higher than 30, those undergoing medical treatment, with implanted electronic devices or thoracic trauma, those who had ingested alcohol 48 h prior, those who frequently use hallucinogens and who practiced yoga or pilates due to atypical breathing patterns [23].

A public and independent call for volunteers was made to participate in each experiment, looking for at least twenty-two participants per experiment in accordance with Equation (1) [24].

$$\mathrm{N}={\left(\frac{\sqrt{\mathsf{\theta }(1-\mathsf{\theta })}{\mathrm{z}}_{\mathsf{\beta }}+0{.5\mathrm{z}}_{\mathsf{\alpha }}}{0.5-\mathsf{\theta }}\right)}^2$$

(1)

where α is the significance value, β the statistical test power and θ the probability that the alternative hypothesis is met. With a reported variability of respiratory pattern of 0.15 (θ = 0.15) [25], α = 0.05 (Z_α = 1.645), and β = 0.99 (Z_β = 2.326), the sample size N should be at least 22.3. According to this result, twenty-two volunteers (age of 25.95 ± 6.54 years, height of 173.18 ± 5.11 cm and weight of 72.90 ± 6.54 Kg) were enrolled to the restrictive maneuver experiment (E1) and thirty-two different volunteers (age of 26.62 ± 5.79 years, height of 172.45 ± 4.48 cm and weight of 72.82 ± 9.43 Kg) were enrolled to the obstructive maneuver experiment (E2) (see Figure 1).

2.3. Measurements and Signal Acquisition

2.3.1. Respiratory and Surface Electromyography Signals

Airflow, airway pressure and tidal volume signals were recorded using a datalogger attached to the mechanical ventilator (Hamilton Medical, Bonaduz, Switzerland) and sEMG of respiratory muscles was acquired with the biopotential amplifier Bagnoli-8 (Bagnoli EMG system, Delsys, Boston, US, 2003). A sampling frequency of 1024 Hz was used in both equipment for each single channel.

Ag-AgCl electrodes in bipolar configuration were placed on all subjects to record three respiratory muscles, diaphragm (Dia), intercostal (Int) and sternocleidomastoid (Strn). A pair of electrodes was used to record Dia, placed between the seventh and eighth intercostal spaces, located in the middle of the mid axillary line and the external clavicular line [26]. Electrodes placed in the second and third intercostal spaces, at the midclavicular line, were used to record Int [27], and electrodes between the mastoid process and mid-point of the muscle were used to record Strn [26].

2.3.2. Respiratory Mechanics

To estimate the respiratory system compliance, an occlusion for 3 seconds was set in the ventilator at the end of inspiration, the occlusion maneuver was performed by the mechanical ventilator during the last minute in the first respiratory cycle. Once the respiratory system compliance was estimated, an optimization algorithm was applied to calculate the airway resistance (see technical details in [28,29]).

After estimating the airway resistance and respiratory system compliance, the muscular pressure (Pmus) was calculated with equation of motion (2) and WOB with Equation (3) for each inspiratory cycle.

$${\mathrm{P}}_{\mathrm{mus}}=\frac{1}{\mathrm{C}}{\ast \mathrm{V}}_{\mathrm{T}}+\mathrm{R}\ast \mathrm{Q}+{\mathrm{PEEP}}_{\mathrm{T}}-{\mathrm{P}}_{\mathrm{aw}}$$

(2)

$${\mathrm{WOB}}_{\mathrm{r}}=\frac{1}{{\mathrm{V}}_{\mathrm{T}}}{{\displaystyle \int }}_0^{\mathrm{t}}{\mathrm{P}}_{\mathrm{mus}}\times \mathrm{Q}\mathrm{dt}$$

(3)

where P_aw, Q, and V_T represent pressure measured at the mouth, airflow, and tidal volume respectively. PEEP_T is total positive pressure at the end of expiration, C is respiratory system compliance, R is airway resistance, and t is the inspiratory time [30].

According to the literature, depending of the subject’s effort and respiratory load, the WOB takes values in different ranges that can be divides into three groups [4,5,6,7,8], G1: [0.3–0.8) J/L, G2: [0.8–1.3] J/L, and G3: (1.3–1.8] J/L. The determined group was considered as the true value in the classification process.

2.4. Data Analysis

2.4.1. Indexes from Surface Electromyography of Respiratory Muscles

From the last minute, variable epochs of at least thirty seconds of sEMG signals were selected, avoiding motion artifacts per epoch. Each signal was filtered with a 4th-order Butterworth 60 Hz bandstop filter and a FIR bandpass filter between 10 Hz and 500 Hz [31,32]. Furthermore, a 5th-order RLS adaptive filter was used to reduce cardiac artifacts [33].

The ventilatory and sEMG signals were synchronized through an auxiliary signal recorded by both the ventilator and EMG amplifier [34] in order to identify the respiratory cycles. All indexes were calculated for each respiratory cycle during the epoch of thirty seconds, and the average was obtained for each sample.

Time Domain Indexes

The amplitude of sEMG signals was obtained through the Hilbert transformation and its maximum was calculated as the difference between maximum and minimum values during the inspiratory phase (see Figure 2) [16]. Information concerning muscle synchronization was obtained by calculating the Pearson’s correlation coefficient between sEMG amplitude signals [15] (see

Table 1. sEMG indexes in time and frequency domains for each one of the recorded muscles, diaphragm (Dia), intercostal (Int) and sternocleidomastoid (Strn).

Index	Time Domain	Index	Frequency Domain
MaxAmpj	Maximum amplitude of the sEMG signal of muscle j: Dia, Int or Strn	Fcj	The central frequency of power spectral density of sEMG signal of muscle j
Rj-k	Pearson correlation between pairs of amplitudes of sEMG signals of the muscles j and k	RHLj	The ratio between high and low frequencies of power spectral density of sEMG signal of muscle j
		Cj-k	Spectral coherence between pairs of sEMG signals of the muscles j and k

).

Frequency Domain Indexes

The power spectral density (PSD) of sEMG signals was estimated using the Burg method with an 8th-order (see Figure 2). Using the PSD, the central frequency and ratio between high and low frequencies, RHL (LF: 20–40 Hz and HF: 138–240 Hz) were calculated to quantify activation level of different types of muscle fibers [17,21]. Finally, the coherence between sEMG signals of different muscles was obtained to identify the muscle engagement in the frequency domain (see Table 1).

2.4.2. Work of Breathing Estimation

Figure 3 summarizes the methodology to estimate the level of WOB from sEMG signals. The first part consists of determining which indexes from the sEMG signals are relevant in describing the breathing pattern and the work of breathing. The second part includes two classification steps. The first one identifies whether the breathing pattern corresponds to a restrictive or obstructive one. Once the breathing pattern is identified, the subjects are classified into one of the three WOB levels, as established in the literature [4,5,6,7,8]. Lastly, the work of breathing is estimated using information from respiratory muscle sEMG signals.

• Relevant Indexes Selection

Before the first classification, a dimension reduction was performed to choose the most relevant indexes (see Section 2.3.1). The following two algorithms of relevant indexes selection were used: the stepwise regression and the mutual information criterion. Statistically significant indexes were examined using Kruskal-Wallis test, followed by the post hoc Tukey-Kramer test for multiple comparisons. Additionally, a Mann–Whitney test on each pair of groups was done to determine which groups differed from each other.

1. Stepwise Regression

The stepwise regression is a systematic method for adding and removing indexes based on their statistical significance in a regression [35]. In this study, the p-value for determining the use of an index was 0.05. The summarized sequential algorithm steps are [35]:

• Calculus of correlation coefficients of n-indexes with the response variable Y.
• Sort the n-indexes from highest to lowest correlation coefficient X = [x1, x2…..xn].
• Add index xi to the P vector.
• Calculus of regression and partial F_j-values between response variable Y and P, where j is the index in the P vector.
• Selection of the lowest F_j-value and comparison between this and the F-value at α = 0.05 in the F-distribution.
• If the F_j-value is higher than the F-value at α = 0.05, i is incremented by one, the index P_j is removed from the P vector.
• Steps 3–6 are repeated until the last index enters the P vector.

The final P vector contains the relevant indexes.

2. Mutual Information

The mutual information gives an idea of the general statistical dependence between two variables [36]. In this study, the mutual information was calculated between each one of the sEMG indexes and the target group through Equation (4).

$$I\left(X;Y\right)=H\left(X\right)-H\left(X,Y\right)$$

(4)

where X was each one of the indexes, Y the target group, and H the entropy of a random variable defined by Equation (5).

$$H\left(X\right)=-{{\displaystyle \sum }}^{}p\left(X\right)lo{g}_{2}p\left(x\right)$$

(5)

where p(X) is the probability density function of index X.

The indexes were sorted from higher to lower mutual information, and each one was examined with a multiple comparisons test.

3. Kruskal-Wallis and Post Hoc Analysis

The indexes selected with each one of the algorithms were examined with the Kruskal-Wallis test followed by the post hoc Tukey-Kramer multiple comparisons correction and Mann–Whitney tests on each pair of groups. Only the indexes that separate statistically (p-value < 0.05) in any of the groups were considered in the machine learning process.

• Machine Learning

Steps 2.1 and 2.2 in Figure 3 identify the respiratory muscular pattern and the WOB level, respectively.

Due to the fact that in clinical practice, knowing the ventilatory pattern is not common, its identification is necessary before the WOB level classification. In this study, the ventilatory pattern was defined according to the respiratory muscular activity, in order to know what muscles, participate in the ventilatory process.

Several supervised learning techniques were selected according to their wide and successful use in classifying biomedical data [18]. The techniques tested were: LDA, SVM, KNN, an ensemble algorithm, and NN. The parameters of the three first classifiers (LDA, SVM, and KNN) were optimized using the Bayesian optimization through the option of the MATLAB algorithms. In the case of the Neural Network, three hidden layer sizes, with default values of 10, 5 and 20, were tested.

Table 2. Information about configuration and validation of each of the classification algorithms. The algorithms used in each of the cases, classification of the respiratory muscular pattern and WOB level, are specified. LDA: Linear Discriminant Analysis, SVM: Support Vector Machine, KNN: K-Nearest Neighbor, NN: Neural Network.

Respiratory Muscular Pattern
Classifier Algorithm	Configuration		Validation
LDA	Gamma = 0 Delta = 0 Discriminant type: linear Indexes: logarithmic transformed		K-fold 20%
SVM	Box Constraint = 962.47 Kernel Scale = 8.4 Kernel Function: Gaussian Indexes: Standardized		K-fold 20%
KNN	Distance: Euclidean Number of Neighbors: 2 Indexes: Standardized		K-fold 20%
WOB LEVEL
Classifier Algorithm	Experiment 1	Experiment 2	Validation
	Configuration	Configuration
LDA	Gamma = 0.49 Delta = 0.06 Discriminant type: linear Indexes: logarithmic transformed		K-fold 20%
SVM	Box Constraint = 1 Kernel Scale = 1 Kernel Function: Gaussian Indexes: Standardized	Box Constraint = 0.7416 Kernel Scale = 110.77 Kernel Function: Gaussian Indexes: Standardized	K-fold 20%
KNN	Distance: Euclidean Number of Neighbors: 2 Indexes: Standardized		K-fold 20%
META-CLASSIFIER KNN	Distance: Euclidean Number of Neighbors: 2 Indexes: Standardized		K-fold 20%
NN	Training function: Levenberg-Marquardt hidden layer: 1 hidden layer sizes: 10 input layer sizes: 4 output layer sizes: 3 Indexes: Standardized		K-fold 20%

shows the parameters obtained for greater accuracy (see Table 2).

For linear discriminant analysis (LDA), the indexes of sEMG were logarithmically transformed based on the fundamental assumption of normally distributed variables for this method. Indexes were standardized for the algorithms where the assumption is not applicable [37] (see configuration details in Table 2).

An ensemble technique was implemented to enhance the power of single classification techniques, looking to combine single classifiers through different strategies. In this study, the stacking algorithm was selected to combine three single classifiers or supervised methods (LDA, SVM, and KNN). The stacking algorithm has two levels (see Figure 4), the single classifiers (C_k) on the first level where the indexes are the inputs (F_n) and the correct classes are the outputs (G_n); at this level the single classifiers predict the classes (φ_k), which will be the inputs for the meta-classifier (supervised method) on the second level, where the outputs are also the correct classes [38]. The meta-classifier implemented in this study was a KNN model and the predicted classes (Φ) were compared with the correct classes.

Once the hidden layer size was determined for the fourth implemented classification algorithm, NN, one thousand neural networks were tested using a different randomized segmentation for training and validation dataset (80% and 20% respectively). Also, the initial weights were randomized for each one. The best NN was selected based on its relation to training and validation loss function.

The accuracy, specificity and sensitivity were calculated for the best algorithms according to Subasi et al. [39]. In the case of the WOB level classification, the calculus was as follows:

• Accuracy: number of samples correctly classified/number of total samples
• Specificity: number correctly classified as normal WOB/number of total normal WOB.
• Sensitivity (Group 2 only): number correctly classified as medium WOB/number of total medium WOB.
• Sensitivity (Group 3 only): number correctly classified as elevated WOB/number of total elevated WOB.
• Sensitivity (Group 2 and Group 3): number correctly classified as high WOB/number of total high WOB.

• Modeling to Estimate the WOB

From previous patient observations [15,17], our research group state the following hypothesis: “high work of breathing implies high respiratory muscle activity, differentiating the diaphragm activity from accessory muscles”. Equation (6) is proposed in this article as a model to quantify work of breathing from the muscular activity of two accessory muscles and diaphragm.

$$WO{B}_e= \frac{{\alpha }_1 \times MaxAmpDia +{\alpha }_2 \times \left(MaxAmpInt + MaxAmpStrn \right) }{MaxAmpDia + MaxAmpInt + MaxAmpStrn}$$

(6)

where MaxAmpDia, MaxAmpInt, and MaxAmpStrn are the maximum amplitudes of the surface electromyographic signals of the diaphragm, intercostal and sternocleidomastoid muscles respectively. α1 and α2 are coefficients that express the contribution of the main and accessory muscles respectively to achieve the necessary work of an inspiratory cycle.

Surface electromyography is a technique sensitive to size and type of tissue over the muscle, which influences the level of attenuation of the signal [18]. The standardization of Equation (6) by the sum of the maximum amplitudes was done to know the contribution of each muscle to the ventilatory process, reducing the effect of different anatomy structures in the surface electromyography recording.

Optimization of Model Parameters

An optimization algorithm based on sequential quadratic programming (SQP) was applied to find α₁ and α₂ values in each WOB group (G). Equation (7) was the function to be minimized.

$$\mathrm{J}=\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}|{\mathrm{WOB}}_{\mathrm{r}}-{\mathrm{WOB}}_{\mathrm{e}}|$$

(7)

where N is the number of samples in each group, WOBr is the real work of breathing represented by Equation (3) and WOBe is the estimated work of breathing from muscular activity, using the proposed Equation (6).

In the optimization algorithm, two nonlinear constraints were proposed to find parameters with a physical meaning. The first constraint, Equation (8), looks for values of α₁ and α₂ that express the contribution of the diaphragm and accessory muscles, the sum of α₁ and α₂ must be equal to the proportion between the average work of breathing ($\overline{{\mathrm{WOB}}_{\mathrm{rG}}}$) and the average maximum amplitude of sEMG signals of all muscles ($\overline{{\mathrm{MaxAmp}}_{\mathrm{G}}}$) for each group (G).

$$\begin{array}{c}{\alpha }_{1_{\mathrm{G}}}+{\alpha }_{2_{\mathrm{G}}}={\mathsf{\beta }}_{\mathrm{G}};\mathrm{where}\mathrm{G}=1,2,3\\ {\mathsf{\beta }}_{\mathrm{G}}=\frac{\overline{{\mathrm{WOB}}_{\mathrm{rG}}}}{\overline{{\mathrm{MaxAmp}}_{\mathrm{G}}}}\end{array}$$

(8)

The second constraint, Equation (9), ensures that the highest coefficient is equivalent to the highest muscular activity. For example, if the maximum diaphragm amplitude is higher than the other two maximum amplitudes (intercostal and sternocleidomastoid), α₁ will be higher than α₂. In addition, if any of the maximum amplitudes of accessory muscles is higher than the maximum diaphragm amplitude, α₂ will be higher than α₁.

$$\begin{array}{c}{\mathrm{MaxAmpDia}}_{\mathrm{G}}{>\mathrm{MaxAmpInt}}_{\mathrm{G}}{\&\mathrm{MaxAmpDia}}_{\mathrm{G}}{>\mathrm{MaxAmpStrn}}_{\mathrm{G}}⟹{\mathsf{\alpha }}_{1_{\mathrm{G}}}{>\mathsf{\alpha }}_{2_{\mathrm{G}}}\\ {\mathrm{MaxAmpInt}}_{\mathrm{G}}{>\mathrm{MaxAmpDia}}_{\mathrm{G}}{\left|\right|\mathrm{MaxAmpStrn}}_{\mathrm{G}}{>\mathrm{MaxAmpDia}}_{\mathrm{G}}⟹{\mathsf{\alpha }}_{2_{\mathrm{G}}}{>\mathsf{\alpha }}_{1_{\mathrm{G}}}\end{array}$$

(9)

The model implementation and the parameter optimization were performed for each data set obtained from the Experiments 1 and 2. The validation technique was mainly K-fold 20%.

2.5. Statistical Analysis

For the normal distribution data, the mean and standard deviation are presented, and all other data are characterized by the median and interquartile range. Normal distribution was tested with the Lilliefors test. To compare two datasets, ANOVA and Kruskal-Wallis analyses were applied to the normal and non-normal distribution respectively. To compare more than two data sets, a Tukey-Kramer test was implemented. For all tests, statistical significance was assumed at p-value < 0.05.

3. Results

Seventy-six and eighty-three samples with work of breathing lower than 1.8 J/L (superior limit according to literature) were obtained from E1 and E2, respectively. Therefore, one hundred and fifty-nine samples were used in the training and testing of the algorithms. Figure 5 shows a scheme of the best model proposed.

3.1. Ventilatory Pattern Classification

The first step was classifying the samples into a respiratory pattern. Four out of the fifteen indexes from dimension reduction (stepwise regression and mutual information) were relevant in differentiating the samples between the restrictive and obstructive patterns. The indexes were coherence between the diaphragm and the accessory muscles (C_Dia-Int and C_Dia-Strn), the central frequency, and ratio between high and low frequencies of the sternocleidomastoid muscle (FC_Strn and RHL_Strn, respectively).

Table 3. Restrictive and obstructive breathing pattern classification. The median values and interquartile ranges of each index are shown.

Index	Experiment		P-Value
	Restrictive (E1)	Obstructive (E2)
CDia-Int	0.67 [0.64–0.73]	0.75 [0.72–0.80]	<0.01
CDia-Strn	0.72 [0.70–0.76]	0.65 [0.62–0.70]	<0.01
FCStrn (Hz)	77.84 [68.24–86.71]	62.37 [58.46–67.88]	<0.01
RHLStrn	169.76 [107.71–229.17]	92.07 [55.00–131.52]	<0.01

shows the median values and interquartile ranges of these four indexes. Furthermore, the p-value that indicates significant differences between E1 and E2 is presented.

The intercostal accessory muscle supports the diaphragm during obstructive maneuvers more than the sternocleidomastoid, with the highest coherence between the diaphragm and intercostal, which indicates synchronization of these muscles in the ventilatory process. In contrast, the sternocleidomastoid supports the diaphragm during restrictive maneuvers, with the highest coherence between diaphragm and sternocleidomastoid, and highest central frequency and ratio between high and low frequencies for the sternocleidomastoid, which indicates that this muscle has a higher engagement in the ventilatory process of restrictive maneuvers.

The classification capabilities of the machine learning algorithms to distinguish restrictive and obstructive patterns were 86.16%, 90.57% and 90.85% for LDA, SVM, and KNN, respectively, with KNN being the best algorithm; therefore, the latter was selected for the final model (see Figure 5).

3.2. WOB Level Classification

The second step was classifying the samples in one level of work of breathing. After applying the stepwise regression and mutual information tests, the MaxAmpDia and MaxAmpInt were relevant in identifying the level of WOB for E1, the MaxAmpDia and MaxAmpStrn for E2, and C_Dia-Int was relevant for both experiments.

Table 4. Median and standard deviation of real work of breathing values for each group. G1: Low WOB; G2: medium-high WOB; G3: elevated WOB. Furthermore, the parameters of each model and the selected index values for classification model are shown.

Parameter/Index	Restrictive (E1)			P-Value
	G1	G2	G3	G1- G2	G1- G3	G2- G3
WOBr (J/L)	0.61 ± 0.13	1.02 ± 0.14	1.52 ± 0.13	<0.05	<0.05	<0.05
MaxAmpDia (mV)	0.024 [0.022–0.030]	0.032 [0.022–0.049]	0.054 [0.040–0.068]	*	<0.05	<0.05
MaxAmpInt (mV)	0.020 [0.016–0.027]	0.032 [0.020–0.043]	0.034 [0.020–0.049]	<0.05	<0.05	*
MaxAmpStrn (mV)	0.018 [0.010–0.024]	0.021 [0.014–0.036]	0.021 [0.017–0.041]	*	*	*
CDia-Int	0.66 [0.62–0.68]	0.67 [0.65–0.75]	0.71 [0.67–0.83]	*	<0.05	*
α1	0.55 ± 10-5	1.07 ± 0.04	1.55 ± 10–5	-	-	-
α2	0.45 ± 10-5	0.93 ± 0.04	1.45 ± 10–5	-	-	-
	Obstructive (E2)
WOBr (J/L)	0.63 ± 0.12	1.03 ± 0.14	1.51 ± 0.12	<0.05	<0.05	<0.05
MaxAmpDia (mV)	0.028 [0.021–0.038]	0.031 [0.023–0.047]	0.050 [0.038–0.083]	*	<0.05	<0.05
MaxAmpInt (mV)	0.016 [0.010–0.029]	0.019 [0.013–0.036]	0.023 [0.015–0.045]	*	*	*
MaxAmpStrn (mV)	0.016 [0.010–0.023]	0.017 [0.013–0.022]	0.030 [0.019–0.048]	*	<0.05	<0.05
CDia-Int	0.79 [0.75–0.84]	0.75 [0.72–0.78]	0.73 [0.69–0.75]	<0.05	<0.05	*
α1	0.55 ± 10–5	1.05 ± 10–5	1.55 ± 10–5	-	-	-
α2	0.45 ± 10–5	0.95 ± 10–5	1.45 ± 10–5	-	-	-

* There are no statistically significant differences between groups.

shows the median values of each of the indexes, with the p-value for each pair of groups.

It should be noted that MaxAmpDia, MaxAmpInt and MaxAmpStrn median values increase when the work of breathing increases. MaxAmpDia separates G₁ and G₂ from G₃ in both experiments. MaxAmpInt separates G₁ from G₂ and G₃ in E1. MaxAmpStrn separates G₁ and G₂ from G₃ in E2. In addition, C_Dia-Int separates G₁ from G₃ in E1 and G₁ form G₂ and G₃ in E2.

It is noted that the features are not related in a linear manner, since the accuracy, specificity, and sensitivity were lower for the LDA model. For the case of the NN model, it is noted that the volunteers with elevated WOB (group 3) were different from the other groups in the restrictive database (E1), with a sensitivity of 100%, in contrast to 80% for the obstructive database (E2), such as group 1 or low WOB. However, in medium WOB, a sensitivity of 100% was achieved in the obstructive database, in contrast to 91.67% in the restrictive database. The above could indicate that during restrictive maneuvers, it was possible to achieve a greater number of samples for level 3 of WOB.

The selection of the neural network was according to the training and validation loss function relation, which is shown in Figure 6 for the restrictive and obstructive dataset. It is possible to observe that depending on the dataset obtained for validation (20%) and training (80%), the model could be overfitted, as is the case with the points located in the fourth quadrant in the figure, where it is possible to observe a training loss function of zero but a validation loss function of almost one, while in the third quadrant, there is the NN model, where the loss function of both datasets can be considered low. Therefore, this NN was selected.

The neural network selected for the final model (see Figure 4 and

Table 5. Classification accuracies, specificities, and sensitivities obtained for each of the machine learning algorithms for both experiments. In the case of neural networks, the accuracy is determined using 20% of the dataset for validation and the between-brackets accuracy is determined using 100% of the dataset.

Algorithm	Accuracy (%)		Specificity (%)		Sensitivity (%) (Group 2)		Sensitivity (%) (Group 3)
	E1	E2	E1	E2	E1	E2	E1	E2
LDA	56.58	55.42	57.69	17.39	63.89	75.56	35.71	53.33
SVM	47.37	56.63	0	0	100	100	0	0
KNN	39.47	54.22	34.62	52.17	41.67	64.44	42.86	26.67
STACKING (LDA-SVM-KNN)	35.53	46.99	53.85	17.39	25.00	77.78	21.43	0
STACKING (LDA-KNN-KNN)	40.79	48.19	42.31	17.39	44.44	75.56	28.57	13.33
STACKING (SVM-KNN-KNN)	47.37	37.35	42.31	65.22	66.67	24.44	7.14	33.33
STACKING (LDA-SVM-KNN-KNN)	50.00	43.37	34.62	34.78	66.67	62.22	35.71	0
NN	80.00 [93.42]	82.35 [91.57]	92.31	82.61	91.67	100	100	80

) obtained an accuracy of 89%, a specificity of 83.67%, a sensitivity of 96.36% (group 2 and group 3), a sensitivity of 90.12% (group 2), and a sensitivity of 89.66% (group 3) (see

Table 6. Confusion matrix for the final classification model of Figure 4.

		Real Values
		Group 1	Group 2	Group 3	Total
Predicted Values	Group 1	41	4	0	45
	Group 2	5	73	3	81
	Group 3	3	4	26	33
	Total	49	81	29	159

).

3.3. WOB Estimation

The last step was estimating the work of breathing value. Figure 7 shows the relation between the real work of breathing (WOBr) and the estimated work of breathing (WOBe) with coefficients α₁ and α₂ (see Table 4) of the groups found with the NN model. The real WOB and the estimated WOB have a Pearson correlation coefficient higher than 80%.

Figure 6 shows the dispersion of the data, where it is possible to observe that this is smaller for G3, offering higher precision in the prediction for cases that have a risky Work of Breathing (near the superior level). It is also possible to observe the samples that were classified with an incorrect WOB level, being points placed near the next or before the neighboring level, exactly in the overlapping region.

4. Discussion

It is widely known that the use of mechanical ventilation is a treatment that causes complications when used for prolonged periods, increasing patient mortality [40]. The work of breathing (WOB) is an accepted index to know the state of the ventilatory mechanics and respiratory muscles, nevertheless, the monitoring of WOB during spontaneous ventilation is not easy to perform [10].

This pilot study seeks to estimate the WOB in both restrictive and obstructive situations, taking advantage of its relationship with the respiratory muscle activity. Indexes from respiratory muscle sEMG signals were used to classify each one of the muscular patterns (restrictive and obstructive maneuvers) and each one of the levels of WOB.

Three indexes from the frequency domain, FC, RHL of the sternocleidomastoid and coherence (C) between the diaphragm and the accessory muscles discriminated the patterns between restrictive and obstructive (see Table 3). The FC and RHL of the sternocleidomastoid were higher in the restrictive pattern than in the obstructive one, which indicates that the activity of the sternocleidomastoid is located in the high frequency band, therefore, it has a high activity [17,21]. The coherence between the diaphragm and any accessory muscles makes it possible to infer the level of engagement of the involved accessory muscle, the sternocleidomastoid in the restrictive case. On the contrary, in the obstructive case, the muscle that supports the ventilatory process is the intercostal muscle, as observed in the coherence value between the diaphragm and intercostal, which is higher in this pattern than in the restrictive one.

Concerning WOB classification, two of the calculated indexes presented a statistically significant contribution to discriminate the samples: the amplitude and the coherence (see Table 4). The median of the maximum amplitude of sEMG signals of the diaphragm increased with increments in work of breathing in both experiments, confirming that this is the principal muscle in the ventilation process. A differentiated pattern was found related to the accessory muscle activity, the more significant accessory muscle for the restrictive maneuver (E1) was the intercostal muscle (maximum amplitude of sEMG) and in the case of the obstructive maneuver (E2), the significant changes were due to the sternocleidomastoid muscle. For the same reason, the coherence between the diaphragm and the intercostal increased in E1 and decreased in E2.

The findings show that during the ventilatory process of a subject with a restrictive pattern, the muscles with permanent activation are the diaphragm as the main muscle and the sternocleidomastoid as the accessory muscle engaged. However, as the work of breathing increases, there is also an increase in the intercostal muscle activity. Conversely, in a subject with an obstructive pattern, where the muscles are permanently activated are the diaphragm and the intercostal, but when the work of breathing increases, the muscle that showed increments was the sternocleidomastoid.

To identify the muscular pattern between restrictive or obstructive, the best technique was the KNN, obtaining an accuracy of 91%. Meanwhile, to identify the level of work of breathing, LDA, SVM and KNN showed accuracies lower than 60%. Combining the properties of some of these classifiers, the accuracy reached is lower than 80%. However, the neural network was the algorithm with the greatest accuracy, achieving a mean of 89% for the all dataset (see Table 5). It should be noted that the accuracy of the classification of WOB level depends on the real WOB (WOBr), and this was calculated from values of compliance, resistance, and muscular pressure, obtained with an occlusion maneuver and optimization algorithms, thus WOBr may have an estimation error [28,29]. In future work, the classifiers should be validated to increase the accuracy, recording the real WOB with invasive methods in patients with invasive mechanical ventilation.

The proposed algorithm might be useful in critical care medicine, due to its ability to identify the level of work of breathing in a non-invasive manner. With a specificity of 83.67% and a sensitivity (Group 2 and Group 3) of 96.36% (see Figure 5 and Table 6), it is possible to predict who could have a successful weaning, which in terms of low WOB is in line with the study of Kirton et al. [8], where 96% of patients with a work of breathing lower than 0.8 J/L had a successful weaning process.

The model proposed in this study relates work of breathing with the maximum amplitudes of sEMG signals of the diaphragm and two accessory muscles, the intercostal and sternocleidomastoid. However, due to the normalization process, it was necessary to find different coefficients according to the level of WOB. The coefficients indicate the contribution to the ventilatory process of each of the muscles. In all cases, the muscle with the most contribution was the diaphragm with 55% for the first group and 53% for the other groups; additionally, the remaining contribution was distributed between the other respiratory muscles.

The coefficients (α1 y α2) in both experiments were similar, because the minimum error for all the dataset was in the median value of the group. In that case, it was possible to consider only one model for both cases, restrictive and obstructive patterns.

The proposed approach estimates the value of work of breathing with an accuracy of 80.43%. Other authors have shown non-invasive techniques to estimate the work of breathing. Vicario et al. [41] presented a model that estimates the work of breathing of mechanically ventilated patients, using ventilatory signals such as airway pressure, volume, and airflow. Their model offers favorable results in cases in which the patients have normal and controlled respiratory efforts, but it has limitations for uncontrolled efforts. Banner et al. [11,20] presented an Artificial Neuronal Network to estimate the power of breathing or work of breathing per minute from parameters like spontaneous minute ventilation, inspiratory trigger pressure, inspiratory flow rise time, intrinsic PEEP and respiratory muscle pressure, where the last two parameters depend on the respiratory system compliance and airway resistance calculated from controlled ventilations. The model presented in this study can predict the WOB at different levels during total spontaneous ventilation and has the advantage that it does not need ventilatory signals, i.e., it is possible to know the work of breathing in situations that do not require mechanical ventilation or devices to record ventilatory signals that can affect the respiratory mechanics.

This pilot study suggests that the work of breathing can be estimated by measuring the activity of respiratory muscles using surface electromyography. In accordance with other authors who inferred it in previous works, like Salazar M, et al. [15], who found a promising index to assess the participation of respiratory muscle on spontaneous breathing test in poisoned patients, Schmidt M, et al. [14], who showed that measurements of accessory inspiratory muscles by sEMG can contribute to a reliable non-invasive indicator of respiratory effort in mechanically ventilated patients, and Hernandez et al. [21], who showed that changes in indexes from sEMG of the diaphragm and sternocleidomastoid are related to changes in the level of PEEP.

Furthermore, the work of breathing was estimated in different ventilatory patterns in spontaneous ventilation, which corroborates the close relationship, not only of the diaphragm, but of the accessory muscles with the ventilatory process. This offers a promising technique in different areas, mainly in critical care medicine, because the information obtained with this model will not only indicate the state of the patient during mechanical ventilation, but it will also help medical staff decide whether or not a patient needs treatment with mechanical ventilation, determine the right extubation moment and optimize the mechanical ventilator settings. However, for future work, a database of patients in ICU would be recorded to validate the model and its possible use in the daily routine of healthcare professionals.