Model-Based Estimation of Inspiratory Effort Using Surface EMG

<italic>Objective:</italic> The quantification of inspiratory patient effort in assisted mechanical ventilation is essential for the adjustment of ventilatory assistance and for assessing patient-ventilator interaction. The inspiratory effort is usually measured via the respiratory muscle pressure (<italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{mus}}$</tex-math></inline-formula>) derived from esophageal pressure (<italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{es}}$</tex-math></inline-formula>) measurements. As yet, no reliable non-invasive and unobtrusive alternatives exist to continuously quantify <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{mus}}$</tex-math></inline-formula>. <italic>Methods:</italic> We propose a model-based approach to estimate <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{mus}}$</tex-math></inline-formula> non-invasively during assisted ventilation using surface electromyographic (sEMG) measurements. The method combines the sEMG and ventilator signals to determine the lung elastance and resistance as well as the neuromechanical coupling of the respiratory muscles via a novel regression technique. Using the equation of motion, an estimate for <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{mus}}$</tex-math></inline-formula> can then be calculated directly from the lung mechanical parameters and the pneumatic ventilator signals. <italic>Results:</italic> The method was applied to data recorded from a total of 43 ventilated patients and validated against <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{es}}$</tex-math></inline-formula>-derived <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{mus}}$</tex-math></inline-formula>. Patient effort was quantified via the <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{mus}}$</tex-math></inline-formula> pressure-time-product (PTP). The sEMG-derived PTP estimated using the proposed method was highly correlated to <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{es}}$</tex-math></inline-formula>-derived PTP (<inline-formula><tex-math notation="LaTeX">$\mathsf{r}=\text{0.95}\pm \text{0.04}$</tex-math></inline-formula>), and the breath-wise deviation between the two quantities was <inline-formula><tex-math notation="LaTeX">$-\text{0.83}\pm \text{1.73}\,\text{cmH}_\text{2}\text{O}\text{s}$</tex-math></inline-formula>. <italic>Conclusion:</italic> The estimated, sEMG-derived <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{mus}}$</tex-math></inline-formula> is closely related to the <italic>P</italic><inline-formula><tex-math notation="LaTeX">$_{\mathsf{es}}$</tex-math></inline-formula>-based reference and allows to reliably quantify inspiratory effort. <italic>Significance:</italic> The proposed technique provides a valuable tool for physicians to assess patients undergoing assisted mechanical ventilation and, thus, may support clinical decision making.


I. INTRODUCTION
A SSISTED mechanical ventilation modes are based on the idea that the respiratory workload is shared between patient and machine, thus relieving the patient while encouraging spontaneous breathing. Recently, it has been shown that poor patient-ventilator interaction is very common during assisted mechanical ventilation and highly detrimental to the patient [1]. Inadequate assistance levels, leading to diaphragm atrophy or injury, as well as excessive driving pressures, leading to selfinflicted lung injury, have to be recognized and avoided. Thus, monitoring the patient's own inspiratory effort is crucial for titrating a proper level of assistance and improving lung and diaphragm health [2]- [4].
The current clinical standard for the quantification of inspiratory effort is based on esophageal pressure (P es ), measured with an invasive esophageal balloon catheter, which allows derivation of the respiratory muscle pressure (P mus ) [5]. Esophageal pressure measurements, however, have not yet been fully adopted in most clinics because of the technical intricacies involved in the placement and calibration of the esophageal balloon [4], [6], [7]. Also, the measured pressure is often heavily perturbed by cardiogenic artifacts and has to be corrected for the elastance of the patient's chest-wall to retrieve P mus [5].
As an alternative to P es , the electrical activity of the respiratory muscles has been proposed for assessing patient activity [8]- [11]. The electrical activity of the diaphragm (EAdi) can be measured using a special nasogastric catheter equipped with an electrode array. EAdi measures the electromyogram (EMG) of the crural part of the diaphragm and is assumed to reflect the neural respiratory drive of a patient. EAdi has been shown to be tightly correlated with transdiaphragmatic pressure (P di ) and total respiratory muscle pressure (P mus ) [12], [13]. Bellani et al. [13] have shown that P mus can be closely approximated by multiplying EAdi with a patient-specific neuromechanical coupling index, relating the measured EMG to the generated muscle force/pressure. The authors also found that this index is stable within one patient over time and that it can be estimated via airway occlusion maneuvers. However, EAdi shares the disadvantage of requiring the correct placement of an invasive nasogastric catheter with P es . The respiratory surface EMG (sEMG) represents a promising, non-invasive alternative This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ to EAdi and can be measured in a straightforward manner with electrodes placed on the thorax [14]. Compared to EAdi, sEMG offers the advantage that the electrical activity of accessory respiratory muscles, such as the intercostal muscles, can also be recorded. This is particularly useful in patients who exhibit a predominantly thoracic breathing pattern. Similar to EAdi, recent studies [15], [16] have demonstrated an occlusion-based technique for estimating P mus from sEMG. A patient-specific neuromechanical coupling index was determined during occlusions and an approximate P mus signal was derived by scaling the respiratory sEMG channels.
The occlusion-based estimation of P mus from EAdi and sEMG has some practical drawbacks. Firstly, occlusion manuevers necessarily interfere with the normal breathing of the patient and, therefore, their use should be minimized. Secondly, the occlusion-based approach only uses EMG data to form a P mus estimate and does not incorporate the pneumatic measurements available from the ventilator during assisted mechanical ventilation. Therefore, the quality of the estimated P mus depends strongly on the validity of the airway occlusion data and on the signal-to-noise ratio of the measured EMG.
The model-based estimation of respiratory effort was proposed by Yamada et al. [17], who determined lung mechanical parameters, namely the respiratory elastance (E) and resistance (R), on the passive patient and then calculated P mus during assisted ventilation using non-invasive measurements of airway pressure and flow. Similar approaches were also investigated by other authors [18]- [21]. Fewer attempts have been made to estimate lung mechanical parameters and breathing activity simultaneously during assisted mechanical ventilation, which is a challenging task due to the underdetermined nature of the problem [22]. Researchers have proposed different ways to ensure identifiability during spontaneous breathing: Younes et al. [23]- [25] proposed to use different ventilatory maneuvers to identify lung parameters and derive a P mus signal. This approach is promising but the maneuvers necessarily interrupt the normal ventilation, potentially disturbing the patient. An alternative approach for achieving identifiability during spontaneous breathing uses physiological knowledge about the shape of respiratory efforts to constrain the estimated P mus waveform. This was demonstrated by incorporating a P mus template waveform or physiologically motivated constraints into the parameter estimation [22], [26], [27]. These approaches have so far only been demonstrated on small patient cohorts and it is still unclear if the assumptions about the shape of P mus will hold in real ventilation scenarios.
In this work, we pursue a novel model-based approach that incorporates respiratory sEMG signals as additional measurements to gain further information about P mus . We have recently proposed a concept for neurally adjusted ventilation based on incorporating electromyographic measurements in [28]-here, we aim to substantially extend this technique and validate it on a large cohort of 43 patients. The herein proposed method allows to non-invasively estimate the lung parameters and inspiratory effort by combining ventilator signals with respiratory sEMG during assisted breathing. The estimation employs a sensor fusion approach, which uses both the EMG data and pneumatic measurements over the full course of the ventilation to form a P mus estimate that is only weakly affected by noise in the individual signals; see [29] for a comprehensive overview of sensor fusion methods. The performance of the proposed method is evaluated by comparison against the invasive P es -based measurement, and compared to the performance of the previously proposed occlusion method [16]. One of the main advantages of a modelbased approach lies in the improvement of signal quality when all measurements are combined, thus, the novel method attains a better estimation accuracy than the occlusion-based method. Also, by using a model-based approach, P mus can be estimated without using highly obtrusive maneuvers such as airway occlusions. Rather, the novel approach allows for P mus estimation with less obtrusive maneuvers that aim to increase breathing pattern variability, such as changes in pressure support level.

A. Signal Processing
Respiratory sEMG signals were measured via two pairs of electrodes located on the thorax at the following positions: (1) bilaterally at the lower costal margin on the midclavicular line and (2) bilaterally in the second intercostal space on the parasternal line. The two channels were measured differentially between the electrode pairs and sampled at a rate of 1000 Hz. A common electrode was placed above the sternum. As cardiac artifacts are very prominent in thoracic sEMG, a simple gating technique was used to suppress the ECG signal. To this end, the R peak of the ECG was detected and segments affected by the QRS complex were discarded. Finally, a smooth EMG amplitude signal was formed using a moving 250 ms root mean square filter; the filtered signal is hereafter referred to as the envelope signal, refer to [30]. These processing steps were performed on all available channels, resulting in one sEMG envelope signal for each channel. We have previously shown that respiratory sEMG channels often do not carry the same amount of information and that a channel selection method, i.e., using the channel with higher signal-to-noise ratio (SNR), can lead to a higher correlation to P mus [16]. Therefore, we employ our previously described automatic selection technique and only use the EMG channel with the highest SNR for the EMG-derived measure of inspiratory effort. To this end, for each channel, the SNR was approximated by forming the ratio between the maximum EMG amplitudes reached during tidal breathing and the level of the measurement noise in between patient activities (both noise and signal power were estimated across full recordings), see [16] for further details. The selected envelope with higher estimated SNR is denoted by EMG sel (t). The sEMG envelopes are often subject to substantial offsets in the order of several µV due to measurement noise, which can be observed during the pauses between patient efforts. This was accounted for by calculating a baseline EMG sel,0 (t) and subtracting it from the selected envelope. To compensate for possible changes of the offset (e.g., due to fluctuations in noise level) an adaptive, time-varying baseline was calculated using the first tercile of envelope values within a moving 5 s window. Airway pressure (P aw ) and flow (V ) were collected from the ventilator and a volume signal V (t) was calculated directly fromV ; see Section III-A for further details. is measured with an air-filled esophageal balloon and the respiratory muscle pressure is calculated as the difference between P es and the elastic recoil pressure of the chest-wall E cw V (orange line). (b) Occlusion-based estimation: The surface EMG is measured at multiple positions on the thorax and the best channel EMG sel is automatically selected. Using multiple subsequent occlusion maneuvers, a neuromechanical coupling indexK (occl) EMG is determined, which is then multiplied with the EMG to retrieve an estimate for P mus , denoted asP (occl) mus,EMG . (c) Model-based estimation: All available ventilator signals (P aw ,V, V) and the selected EMG channel (EMG sel ) are incorporated into a parameter estimation algorithm, which is based on a model of the respiratory system. Two variants are investigated in this study. In the first variant (mdl), all parameters are estimated during tidal breathing on multiple pressure support levels. In the second variant (mdl+occl) the neuromechanical coupling index K EMG is determined during occlusion maneuvers and the remaining parameters are estimated during tidal ventilation. In both versions, a P mus signal is calculated using the pneumatic signals and the estimated set of parameters. The two estimates are denoted byP (mdl) mus,EMG andP (mdl+occl) mus,EMG , respectively.

B. Occlusion-Based Effort Estimation
It was previously proposed to convert between electrophysiological signals (e.g., EAdi or sEMG) and muscle pressure (e.g., P di or P mus ) using a linear neuromechanical coupling index, defined as K EMG = ΔP mus /ΔEMG (or K EMG = ΔP di /ΔEMG) and measured in cmH 2 O/µV, cf. [13], [15], [16]. Beck et al. [12] showed that the assumption of linearity between EAdi and P di is valid over a wide range of physiological activities and different levels of pressure support. These findings were confirmed in later studies for both EAdi and sEMG [13], [15], [16]. This would allow to produce an estimate for the muscle pressure by multiplication of the EMG signal with K EMG . When P es is available, the value of K EMG can be determined directly by fitting the available sEMG channels to the muscle pressure P mus , cf. [16] for details.
Several groups have investigated whether the neuromechanical coupling index K EMG can be accurately identified without P es using occlusion maneuvers for EAdi [13], [31] and sEMG [15], [16]. During end-expiratory occlusions, when flow and volume are zero, the inspiratory muscle pressure P mus is directly visible in P aw , such that the EMG can be related to the generated pressure. It was previously shown that the neuromechanical coupling index measured during occlusions is well correlated with but higher than the coupling index observed during normal breathing, due to the isometric configuration of the diaphragm. A constant correction factor of 0.8 was therefore proposed to account for the systematic deviation [15], [16]. In our previous work [16], K EMG was determined over the course of multiple subsequent occlusions by means of linear regression between the selected sEMG channel EMG sel and P aw . The resulting index was corrected with the factor 0.8 and denoted asK (occl) EMG , with the hat symbol indicating estimated parameters. An sEMG-derived estimate of P mus was then calculated viâ In this work we compare the performance of the occlusion-based estimate from our previous study [16] to the herein proposed model-based approaches. The occlusion-based estimation is attractive due to its simplicity but strongly relies on a good EMG quality to form a reliable P mus estimate. Moreover, it does not use all available signals, i.e., the pneumatic signals during tidal breathing are not incorporated. A visual overview of all considered methods is provided in Fig. 1 including the occlusion-based method described in this section, and the model-based methods described in the following.

C. Model-Based Effort Estimation
The occlusion-based method can already be understood as a first, simple model-based method because a model about the relation of different measurements (P aw and EMG) is combined with the occlusion data to derive parameters of the physiological system (K EMG ) and determine P mus . In the following, we further pursue this model-based path and attempt to incorporate all measurements available during tidal ventilation to derive respiratory parameters and to estimate a reliable muscle pressure P mus with a high signal-to-noise ratio.

1) Model:
The standard linear respiratory model reads and is commonly referred to as the equation of motion (EqM). It assumes a linear resistance R as well as a linear elastance E (with compliance C = 1/E), P 0 = PEEP, and V 0 is a volume offset due to intrinsic PEEP (iPEEP). It can be easily observed that this equation is underdetermined: with no further knowledge about P mus , the respiratory parameters E and R cannot be determined from the available pneumatic measurements P aw ,V and V [22]. As in the occlusion-based methods, we model the relation between P mus and EMG through a neuromechanical coupling index via P mus (t) = K EMG · EMG(t). Combining this with (2), we obtain which now is in the form of a multivariate linear regression model and, given linearly independent measurements, can be solved for the parameters with appropriate regression methods. As described in Section II-A, we employ the channel selection strategy and only use the most informative sEMG channel, i.e., The approach uses the EMG channel with higher SNR and thereby improves correlation to P mus [16]. Given some estimate for the respiratory mechanical parameters (R,Ê,K EMG andV 0 ), the model-based estimation of the muscular pressure is done usinĝ which is a fully pneumatic reconstruction and coincides with the P mus calculation that is used, e.g., in PAV+ [32], and other model-based approaches for inspiratory effort estimation [17], [18]. In our model, the estimation of parameters becomes possible through the incorporation of EMG. The proposed P mus reconstruction does not directly use the EMG signal and therefore offers an improvement in the SNR over using the original EMG directly. Notice that the regression model in (3) can be solved on data collected during normal assisted ventilation, even when no special maneuvers are performed.
2) Parameter and Effort Estimation: The determination of 'good' parameters for (3) proves to be a difficult task, and the naïve solution using ordinary least squares (OLS) regression on all measured data points during spontaneous breathing appears to be brittle [21], [28], [33], [34]. In fact, the samplewise OLS solution in some cases fails to capture the 'big picture' of the ventilation data and tends to overfit samples that a user would deem less important, such as uninformative time phases, identical repetitive breaths, or artifacts. In this work, we put a special focus on the estimation of respiratory parameters and propose a novel method that manages to fit clinically useful models to the data. Two key novelties are introduced for parameter estimation: (1) the determination of the respiratory time constant τ = R/E (cf. [35]) is done separately by first analyzing expirations only, and (2) an integrated breath-wise form of the equation of motion is introduced to determine all remaining parameters. In both steps, we use robust regression methods, which represent an alternative to OLS when data are contaminated by outliers, refer to [36] for an overview of robust regression.
For the estimation of the respiratory time constant τ , we exploit the fact that under assisted spontaneous ventilation, most patients tend to be (more) passive during expirations. This can be recognized by the typical exponential flow curve and would allow a direct estimation of τ from the slope of the flow-volume loop and has already been exploited in previous works [37], [38]. Unfortunately, in many cases, expirations are confounded by residual patient activity, e.g., missed efforts or active expirations, which complicate the automatic determination of the time constant. The approach taken here is to fit the linear model simultaneously to many expirations via robust regression. The model in (5) can be directly derived from the general EqM (2) when P aw = P 0 and P mus = 0 are substituted. To account for potential changes in absolute lung volume, we use a parallel slopes regression model, i.e., the slope τ is shared between all included expiratory datapoints while an individual volume offset V 0,k is used in each expiration; k indicates breaths. Thus, the parameter V 0,k accounts for breath-wise volume offsets, for example due to intrinsic PEEP. We solve the regression using iteratively reweighted least squares with Tukey's loss function [36] and denote the estimated parameters asτ and V 0,k . As a final step, a continuous signalV 0 (t) is generated from the individual, breath-wise offsetsV 0,k . This is achieved by smoothly interpolating the individual offsets between the expirations; details are provided in Appendix A. The resulting signalV 0 (t) provides a smooth and continuous estimate for the offset V 0 over the full length of the signal. The herein used robust regression procedure rejects rare outliers, which may, e.g., correspond to occasional expiratory patient activity. These outliers thus do not (negatively) influence the estimation of the time constant τ . With OLS regression, on the other hand, any outliers have an out-sized effect on the estimation result. Robustness with respect to expiratory activity can be further increased through a selection strategy, i.e., through additional criteria for selecting 'good' expirations: here, we use very simplistic criteria based on the length and depth of breaths. Beyond that, more advanced selection strategies are possible for selecting mostly passive segments, e.g., based on the shape of the expiratory flow curve similar to [39], or based directly on the EMG.
The determination of the time constant τ solves an important part of the full parameter identification problem, because we can now fix the ratio of R and E.
where the artificially generated signal is substituted. Given our previous estimatesτ andV 0 (t),V rs (t) is calculated according to (8) and only two model parameters (E and K EMG ) thus remain to be identified in the simplified EqM (7). It should be noted that V rs (t) is an interesting quantity in its own, as it forms a natural counterpart to the total pressure P rs (t) = P mus (t) + P aw (t) − P 0 acting on the respiratory system. Notably, the equation of motion can be rewritten as P mus (t) + P aw (t) − P 0 = E · V rs (t). Thus, V rs (t) is the respiratory system pressure P rs (t) scaled by the elastance E, refer to Fig. 2 for a visualization.  (2), where a superposition of P aw and P mus is used as the input. In such cases, the flowV and volume V tend to be difficult to interpret. Once an estimate for the time constant τ is available, the flow and volume can be aggregated via V rs = V + τ ·V. The resulting signal is proportional to the pressure P rs = E · V rs = P aw + P mus − P 0 acting on the respiratory system and therefore provides an intuitive insight into the timing/amplitude of breathing efforts from the ventilator and patient.
The second novelty for improving parameter estimation is based on a reformulation of the equation of motion into a breathwise form. This idea is based on the observation that both sides of the EqM can be integrated in a breath-by-breath fashion over inspirations, leading to The integral of a pressure signal as in (9) is known as pressuretime product [5], defined as PTP = P (t) dt, and is often employed as a metric for total pneumatic energy expenditure. A very commonly used PTP metric is the muscle pressure-time product (PTP mus ), quantifying the patient's breathing effort [5].
In the following, we denote the integral over volume and EMG in (9) as volume-and EMG-time products (VTP = V (t) dt and ETP = EMG(t) dt). We have recently shown that the integrated surface EMG curve, i.e., ETP, is highly informative as a metric for patient activity [16]. The breath-wise evaluation of (9) results in an integrated equation of motion that is now fully defined in terms of time product values, with k indicating the breath index. In (11), all time products are known: PTP aw is the integrated ventilatory support, VTP rs is estimated by integration over the variableV rs (t) (which in turn was calculated usingτ andV 0 (t)), ETP is the EMG-time product [16], and, finally, PTP 0 is just calculated by integrating the PEEP. One set of time product values is calculated from each inspiration (detected directly fromV ). The remaining parameters E and K EMG can now be simply estimated from equation (9) using a robust regression over these data points, i.e., linearly independent time product values calculated from multiple inspirations. We solve the regression using iteratively reweighted least squares with Huber's loss function [36]. Given the estimated parameterÊ, we can estimate the muscle pressure viaP mus,EMG (t) =Ê ·V rs (t) − P aw (t) + P 0 (12) which corresponds to the pneumatic P mus reconstruction in (4), andV rs (t) is calculated from estimatesτ ,V 0 (t). An estimate for the resistance R can be obtained viaR =Ê ·τ . The pneumatic calculation of P mus can lead to high-frequency artifacts within the reconstructed signal, which may originate from model mismatch or pneumatic measurement errors. These artifacts can be easily removed by means of lowpass filtering because they have much higher frequencies than the physiological P mus signal. This was implemented using a 7th order Butterworth lowpass filter with a cutoff frequency of 3.5 Hz. While this cutoff may seem rather restrictive, our main goal was reconstructing a smooth signal, and in our experience, a 3.5 Hz cutoff did not negatively affect the P mus waveform. In the following, this model-based, filtered estimate for P mus will be denoted byP (mdl) mus,EMG (t), and corresponding parameters byK (mdl) EMG ,Ê (mdl) andR (mdl) . Rather than fitting the respiratory model to single samples, the proposed reformulation in (11) implies a breath-wise interpretation, which means that a model is fitted to explain complete breaths. Thus, the model relates a change in flow and volume over the course of a full inspiration, which is aggregated in VTP rs , to some equivalent aggregated pressure counterpart in either PTP aw or PTP mus . One of the main advantages of this reformulation is that PTP mus (k) = K EMG · ETP(k) is a much weaker assumption than P mus (t) = K EMG · EMG(t), and the measured EMG data only ever enter the model in their integrated form. Our method therefore does not depend on the assumption that the EMG waveform is directly proportional to P mus . Instead, we only require a high degree of correlation between ETP and PTP mus , which we have previously observed experimentally [16]. Our formulation is thus robust to (zero-mean) noise in either of the involved signals (pneumatic or electromyographic), as well as to small time delays between the signals. Such delays have been investigated by, e.g., [33], and pose serious difficulties for a sample-wise regression approach. It can be easily seen in the integrated equation (11), that at least two linearly independent operating points are needed to solve the regression for the two parameters E and K EMG . More precisely, sufficient breath-wise variation of PTP mus (and ETP) against VTP rs must be available within the ventilation data to ensure identifiability of all parameters. This can be achieved either by relying on the natural variation of the patient's spontaneous breathing or by inducing some external variation, e.g., by varying the ventilator pressure support level. However, highly obtrusive maneuvers such as airway occlusions used in [15], [16] are not at all needed for the herein proposed estimation procedure, and parameter estimation can, e.g., be performed during normal tidal breathing at multiple pressure support levels. It should be emphasized that our approach presented in equations (10) and (11) can be applied very generally to determine respiratory parameters by using any effort metric correlated with PTP mus , and is not restricted to using ETP. Finally, notice that a very similar reformulation of the EqM can be done by integrating over volume (instead of time), which leads to a formulation in terms of work of breathing (instead of time products).

D. Combined Approach
In the fully model-based approach described above, all model parameters (R, E, K EMG , V 0 ) are estimated simultaneously during continuous, tidal breathing on multiple pressure support levels. As a final approach, we attempt to combine the benefits of the model-based and the occlusion-based approach. The former uses all available signals during tidal breathing to form a P mus estimate while the latter gives an instantaneous and simple estimate for the neuromechanical coupling index K EMG . In the combined approach, we therefore estimate K EMG during occlusions only, and then assume it to be known during the model-based estimation of the remaining parameters and P mus . The parameter estimation follows the same steps as before, but K (occl) EMG is used for K EMG and not estimated in the final robust regression (11). We denote the results of this second model-based variant byP (mdl+occl) mus,EMG and the corresponding parameters bŷ K (mdl+occl) EMG ,Ê (mdl+occl) andR (mdl+occl) . In this variant, K EMG is determined in exactly the same way as in the occlusion-based method-the only difference between the two approaches being howP mus,EMG is calculated from K EMG .

A. Study Data
To assess the validity of the proposed technique, the novel P mus estimation was tested on study data of 43 intubated patients undergoing bronchoscopy with assisted mechanical ventilation; details on demographic characteristics are given in [16]. The study was registered in the German Clinical Trials Register (DRKS00021524) and data were collected at the department of pneumology, cardiology and intensive care medicine of the Klinikum Konstanz (Konstanz, Germany). The protocol was approved by the ethics committee of Witten/Herdecke University on July 12, 2018 (protocol number 137/2017) and conducted in adherence to the ethical standards laid down in the Declaration of Helsinki in its most current form. The two proposed modelbased estimatesP (mdl) mus,EMG andP (mdl+occl) mus,EMG are validated against P mus derived from esophageal pressure-the occlusion-based estimateP (occl) mus,EMG serves as a benchmark for the two modelbased methods. A visual overview of all methods is provided in Fig. 1.
After patients were enrolled, a nasogastric catheter equipped with an esophageal balloon (Bösch, Gottenheim, Germany) was filled according to [7] and connected to a pressure transducer at the proximal end of the catheter. The EMG signals were amplified and recorded at a sampling rate of 1000 Hz using a dedicated acquisition system and software by Dräger (Drägerwerk AG & Co KGaA, Lübeck, Germany), which also provided proprietary digital filters for the removal of baseline wander from the raw signals. The same acquisition system was used to record time traces of the esophageal pressure transducer at 200 Hz. In addition, airway pressure P aw and flowV were collected from the ventilator, sampled with a rate of 100 Hz and then synchronized with the remaining data. A continuous volume signal was generated by means of running integration overV . Prior to parameter estimation, long-term drifts in the volume signal due to leakages and flow-sensor inaccuracies were reduced using a trend removal algorithm. Briefly, end-expiratory volumes were filtered with a running median and the filtered points were interpolated to form a continuous baseline that was then subtracted from the volume signal. The baseline-corrected volume signal is denoted by V (t) and used in the proposed estimation methods.
At the beginning of the study protocol, the correct positioning of the esophageal balloon was validated via the Baydur maneuver [40]. Following the initial positioning, a series of spontaneous efforts under airway occlusions was recorded. Patients were ventilated with continuous positive airway pressure (CPAP) on a Dräger V500 ventilator (Drägerwerk AG & Co KGaA, Lübeck, Germany) and PEEP = 5 cmH 2 O. Next, patients were switched to pressure support ventilation (PSV), and three support levels (5, 10 and 15 cmH 2 O) were applied in random order, again with PEEP = 5 cmH 2 O.
As a first step towards generating a clean P mus reference signal, cardiogenic oscillations in P es were removed using our previously proposed template subtraction method [41]. Next, we checked and corrected for possible scaling errors in P es which may for example arise from inaccurate balloon positioning. During airway occlusion maneuvers, the pressure drop in P aw can be assumed to be equal to P es and, thus, the ratio of P es and P aw should ideally be close to one, however, even after proper positioning, small scaling errors usually remain. To mitigate this effect and account for positioning errors, P es was fitted to P aw via linear regression across the occlusion data and the resulting linear correction factor was used to scale the full P es signal.
The total respiratory muscle pressure was calculated via where E cw denotes the elastance of the chest-wall. As fully passive ventilation was not reached throughout the protocol, E cw was approximated in a region with overall low patient activity by annotating points in the P es -V loop (Campbell diagram) that were not subject to patient efforts. The value of E cw was calculated as the slope between annotated end-inspiratory and end-expiratory points. Finally, as the absolute value of esophageal pressure often contains positive offsets, we only ever evaluate relative pressure swings of P mus against a baseline calculated during passive expirations.

B. Data Analysis
Invalid segments subject to irregular patient behaviour such as ventilator fighting or coughing were excluded from the analysis. Parameters of the model-based method were determined by fitting the model to all pressure support levels (i.e., CPAP and PSV at 5, 10 and 15 cmH 2 O) andP mus,EMG estimates were calculated as described in II-C. To evaluate the performance of the proposed methods, inspiratory effort was quantified via the pressure-time product, see [5], leading to PTP mus for the P mus reference, PTP (mdl) mus,EMG and PTP (mdl+occl) mus,EMG for the two model-based estimates and PTP (occl) mus,EMG for the simple occlusion-based approach. For the calculation of PTP values, recordings were segmented into inspirations and expirations using a threshold-based detector [16] applied directly to the P mus signal, and time product values were calculated over inspirations. Compared to using a flow-based segmentation, this offers the advantage that missed efforts are also recognized and included in the analysis. As a second effort metric, we calculated aggregated PTP/min values (expressed as cmH 2 Os/min), see [5]. To this end, all detected efforts within each support level were summed and divided by the length of that segment, leading to four PTP/min values for each patient. A reference value for the neuromechanical coupling index K EMG was determined by fitting the EMG directly to P mus by means of linear regression as described in [16].
Results are expressed as mean ± standard deviation. Correlation between variables was quantified by means of Pearson's correlation coefficient r. To analyze slopes and intercepts between estimated PTP mus,EMG and measured PTP mus , we used the linear regression model where K is the slope and P bias is a constant bias term. Here, as we solve the regression over time products, the regressor T k denotes the length of detected efforts and accounts for the integration of bias over each effort. Deviations between estimated PTP mus,EMG and measured PTP mus values were quantified by the mean absolute deviation (MAD) and using the Bland-Altman method for repeated measurements [42]. Differences between estimated lung parameters were assessed by the two-tailed Wilcoxon signed-rank test.

IV. RESULTS
Of the 43 analyzed patients, 16 were diagnosed with chronic obstructive pulmonary disease (COPD), 5 with asthma, 7 with interstitial lung disease and another 4 patients were diagnosed with obstructive sleep apnea syndrome. The remaining patients were diagnosed with rheumatic/infectious diseases or lung cancer and in many patients comorbidities existed. We have analyzed 23.3 ± 4.0 min of ventilation data per patient and a total of 19 540 efforts (454 ± 137 efforts per patient with a length of 0.89 ± 0.31 s). The number of breaths detected by the flowbased segmentation and used during parameter identification was 399 ± 137 per patient. In two patients, no sufficiently long (> 0.35 s) occlusion maneuvers were available, which is why these two patients were excluded from any analysis relying on occlusions. An illustrative excerpt of patient data and the result of the first model-based estimateP (mdl) mus,EMG are shown in Fig. 3.

A. Parameter Estimates
The respiratory parameters estimated by the two model-based methods across all patients are provided in Table I. The estimated resistance and compliance values were physiologically plausible and within the range of values reported for normal and COPD lungs in [43]. Passively measured references for lung mechanical parameters were not available, which is why we focus on evaluating the validity of estimated PTP values and neuromechanical coupling. We compared the neuromechanical coupling index K (mdl) EMG estimated by the model-based approach against the reference K EMG calculated according to [16]. We foundK (mdl) EMG to be highly correlated to K EMG (r = 0.91, p < .001, slope = 0.91), see Fig. 4, providing strong evidence for the validity of a model-based approach for determining the neuromechanical  Table I. In turn, and in accordance with the EqM, the resistance and elastance values in the first model-based estimator (mdl) were larger than those in the combined approach (mdl+occl) (p < .001).

B. Validation Against PTP mus
Next, we validate the estimated PTP mus,EMG effort metrics against P es -derived PTP mus and report results in Table II. The results of the occlusion-based estimate PTP (occl) mus,EMG from [16] serve as a benchmark for the results of the model-based estimates PTP (mdl) mus,EMG and PTP (mdl+occl) mus,EMG . The breath-wise correlation in individual patients between the simple occlusion-based estimate PTP (occl) mus,EMG and PTP mus was r = 0.87 ± 0.09. For the herein proposed model-based estimators, we found a very high correlation to PTP mus (r = 0.95 ± 0.04 for both approaches), which was an improvement over the occlusion-based method (p < .001). Thus, both model-based methods manage to substantially reduce random scatter that is associated with using a scaled sEMG signal as in the simpler occlusion-based approach. The slope K and bias P bias of EMG-derived effort against PTP mus were calculated using the linear regression defined by equation (14). Results are reported in Table II. The bias was small for all three methods. For the model-based estimates PTP (mdl) mus,EMG and PTP (mdl+occl) mus,EMG , we found a systematic overestimation as the slope K was smaller than one (p < .001).
The breath-wise deviation between estimated and measured PTP mus was calculated across all efforts from all patients with sufficiently long occlusions (n = 41 patients and m = 18 341 efforts). As the data contain multiple efforts from each patient, the total breath-wise deviation across all patients was calculated using the approach by Bland and Altman for repeated measurements [42] and results are reported in Table II. The smallest deviation and MAD was found for the combined estimator (PTP (mdl+occl) mus,EMG ). The deviation of the two model-based approaches was also visualized by means of the Bland-Altman plots, see Fig. 5, which show that the deviation is within clinically acceptable bounds (−0.23 ± 1.44 cmH 2 Os) and would allow to assess patient efforts on a breath-by-breath level. As a last step, we evaluated the possibility to quantify aggregated efforts over a longer period of time. To this end, we calculated PTP/min values from each pressure support level and report deviations between PTP mus,EMG /min and PTP mus /min in Table II. The deviation of PTP/min values was also visualized by means of Bland-Altman plots, see Fig. 6. The deviation of the combined estimator (mdl+occl) was −2.37 ± 30.6 cmH 2 Os/min, which is small compared to a typical value range of 0-250 cmH 2 Os/min (and a clinical target range of 100-150 cmH 2 Os/min). We provide additional performance metrics for several baseline estimators (including OLS-based regression) and intermediary variants of the proposed methodology in Appendix B.

V. DISCUSSION
We have shown that a model-based approach for estimating respiratory muscle pressure P mus from sEMG is reliable and allows for non-invasive quantification of inspiratory patient efforts. The proposed model-based estimates PTP mus,EMG were highly correlated to the reference PTP mus obtained from esophageal pressure (P es ) measurements (r = 0.95 ± 0.04 for both PTP In the previously introduced occlusion-based method [15], [16], P mus was approximated by a scaled EMG envelope, which is why its performance strongly depends on the quality of the measured sEMG signal and of the data obtained during the occlusion. The herein proposed model-based method takes all available measurements into account, i.e., both the ventilator and the sEMG signals, and thus reduces measurement noise and does not rely on a single occlusion, thereby improving the quality of the P mus estimate. We have shown that this leads to a significantly improved correlation to the P es -derived P mus . The   [42] for the difference between estimated PTP mus,EMG /min from the two model-based variants and invasively measured PTP mus /min (aggregated over the available pressure support levels). The plots depict k = 161 ventilation phases from n = 41 patients and the coloring corresponds to the four pressure support levels (CPAP and PSV at 5, 10 and 15 cmH 2 O). The plot shows the mean across all values (solid grey line) and 95% intervals (dashed grey lines). advantage of our model-based approach is best motivated by patient recordings with a very low signal-to-noise ratio in the sEMG data. One such example is shown in Fig. 7(b), where the pneumatic P mus estimation substantially improves the signal quality and restores a high-fidelity P mus signal despite the high sEMG noise level. Hence, the proposed method mitigates a potential drawback of the respiratory sEMG, namely its low SNR, and may therefore promote its clinical application.
Another important advantage of the model-based approach is the possibility to gain insights into lung mechanical parameters during spontaneous breathing. We have shown that the estimated neuromechanical coupling indexK (mdl) EMG is highly correlated to the reference value (r = 0.91, p < .001). Passively measured reference values for R and E were not available; therefore, these estimates have to be further validated in future studies. We have investigated two variants for the model-based estimator: in the first variant (mdl) no maneuvers were used, i.e., all parameters were estimated during normal breathing, and in the second variant (mdl+occl), the neuromechanical scaling factor K EMG was taken from occlusions and the remaining parameters were estimated during tidal breathing. Both approaches led to an equally high correlation to P mus and would allow to reliably quantify patient efforts. This demonstrates that a completely occlusion-free estimation procedure (mdl) is feasible. The combined approach (mdl+occl) showed the smallest breath-wise deviation to P mus among all considered methods. The slightly higher breath-wise deviation of the occlusion-free method (mdl) cannot be explained through scatter (because of the high correlation) but must be attributed to an overestimation of P mus , visible in the larger neuromechanical scaling factorK (mdl) EMG and larger lung mechanical parametersÊ (mdl) andR (mdl) as well as in the slope K against P es -derived PTP mus (indicating that P mus is ∼ 21% smaller than the estimateP (mdl) mus,EMG ). Here, it must be taken into account that the reference P es itself can be  Fig. 7. Excerpt of signals from two patients. In both cases, the modelbased estimators manage to reconstruct an estimate for P mus which is strongly correlated to the P es -based measurement (see correlation plots ofP mus,EMG against P es -derived P mus ). (a) Patient response to a changing pressure support level. (b) An example of a patient with low EMG signalto-noise ratio. subject to substantial scaling errors. The overestimation of P mus by the fully model-based estimator (mdl) should therefore not be quickly dismissed as an error against P es but might be considered as an independent result, contributing to our understanding of a patient.
From a technical viewpoint, some important progress with respect to respiratory parameter estimation has been achieved. It was already suggested by several authors that ordinary least squares (OLS) regression for respiratory parameter estimation appears to be brittle: e.g., [33] reported a strong sensitivity of OLS to small time delays in the pneumatic signals. Several authors [21], [28], [34] discussed potential shortcomings of the sample-wise OLS regression approach and presented several cases were it failed to identify good model parameters for the data. Additional results given in Appendix B confirm the potential shortcomings of using the direct OLS-based regression approach. In this work, we introduced several new techniques for robust respiratory parameter estimation. Firstly, a robust regression method was proposed to determine the expiratory time constant over many breaths, while rejecting outliers such as expirations with residual patient activity, as long as they are rare enough. And secondly, an integrated form of the equation of motion was introduced, which allows to estimate all remaining lung parameters, again, by using robust regression instead of OLS. The integrated equation of motion implies a breath-wise interpretation of the ventilation data, i.e., the model is fitted to explain full breaths and not individual samples. Very importantly, the pneumatic and electromygraphic signals are only used in an integrated form, i.e., as the time product values, which reduces the sensitiviy to artifacts, time delays, and the shape of the EMG envelope. Appendix B provides detailed results on the performance improvements achieved by the herein proposed novelties as well as a comparison to simpler OLS-based estimators. One of the strengths of the present study is the large cohort of patients and the large number of efforts on which the method's performance was evaluated.
Some limitations of the study and method need to be discussed. A total of 16 patients were diagnosed with COPD, but only four patients showed a strong airflow limitation (with iPEEP > 4 cmH 2 O). Therefore, we do not have evidence yet for the feasibility of a model-based approach in such cases. We hypothesize that the model will have to be adapted to better represent patients with severe airflow limitation, e.g., by using a switching model with two time constants. In regard of this, other authors have also expressed their concerns on the use of the simple equation of motion in severe COPD [44].
Another potential limitation of the study is the composition of the input data, which do not reflect normal ventilation. Here, the models were fitted to all four available support/activity levels, which contain a variety of workload shares between patient and ventilator and therefore support the identifiability of parameters. More precisely, breath-wise variation in V rs induces identifiability of the elastance E in equation (11) and breath-wise variation in patient activity, i.e., in the sEMG, induces identifiability of the neuromechanical index K EMG . An important future research question is how identifiability, i.e., sufficient variation within the  III  METRICS QUANTIFYING THE RELATION TO PTP MUS FOR BASELINE ESTIMATORS AND DIFFERENT VARIANTS OF THE  ventilation, can be reached with a less obtrusive method. We hypothesize that in most cases, two pressure support levels would be sufficient to solve the regression problem and even relying on the natural variation of the patient's spontaneous breathing might be feasible. Other strategies, for inducing variation, would also be possible, e.g., changing cycling-off criteria or using noisy PSV. These strategies should be investigated in future studies.
As a last important remark, we want to emphasize the usefulness of the proposed estimator beyond the herein discussed effort quantification. Firstly, the pneumatically calculatedP mus,EMG signal might be highly useful for assessing the timing of efforts, e.g., to non-invasively recognize patient-ventilator asynchrony, cf. [32]. Secondly, the described method naturally lends itself to proportional ventilation modes [25], [28], where inspiratory assistance is to be provided in proportion to patient efforts. To this end, the calculatedP mus,EMG signal, when multiplied with a patient-specific proportionality factor, could be used as a control signal for the airway pressure. Such a ventilation strategy would position itself between the already existing approaches for proportional ventilation such as neurally adjusted ventilatory assist (NAVA), where support is given in proportion to the EMG, and proportional assist ventilation (PAV/PAV+) as well as proportional pressure support (PPS), where support is given in proportion to a P mus signal calculated from the equation of motion (using either a maneuver-based parameter estimation or manually chosen parameters).

VI. CONCLUSION
We have shown that patient efforts can be quantified noninvasively using surface EMG measurements with a modelbased estimator. The calculatedP mus,EMG signal is closely related to the P mus signal derived from esophageal pressure measurements and might facilitate the monitoring of inspiratory efforts in patients under assisted mechanical ventilation.

APPENDIX A
The calculation of the continuous volume offsetV 0 (t) comprises the following steps: (A) The estimateτ from equation (5) is used to calculate the artificial signalṼ rs (t) := −V (t) −τV (t) (here, the previous breath-wise estimatesV 0,k are not incorporated, therefore,Ṽ rs (t) still contains expiratory offsets). (B) A running baseline valueV 0 (t) is calculated across all expiratory segments ofṼ rs (t) providing smooth estimates of expiratory offsets. (C) As a last step, these expiratory segments have to be interpolated to generate a continuousV 0 (t) signal covering the full domain. To this end, the missing inspiratory segments ofV 0 (t) are recovered by means of linear interpolation between the adjacent expiratory values.

APPENDIX B
We report performance metrics for additional variants of the proposed methodology in Table III. To this end, several intermediary estimators were applied to the study data described in Section IV, and performance metrics from Section III-B were used to quantify the relation of the estimated PTP mus,EMG values to the esophageal pressure reference PTP mus . As before, all parameters (K EMG , E, R) were estimated on full recordings using all pressure support levels, and muscle pressureP mus,EMG was reconstructed using equation (4) or (12). The first two of the additional variants provide a baseline performance, i.e., OLS was applied directly to EqM (3) using a constant and a timevarying offset parameter, respectively; see Table III for details. The remaining estimators constitute different combinations of the novelties proposed in this paper, namely (1) the usage of robust regression instead of OLS, (2) the separate estimation of τ and (3) the integrated breath-wise formulation of the EqM. The highest correlation and smallest bias to PTP mus was reached by combining all aforementioned techniques (corresponding to theP (mdl) mus,EMG estimate, see last row in Table III). It also attained the smallest standard deviation of breath-wise deviation to PTP mus . As discussed before, the estimators reaching the highest correlation also showed a slight overestimation of PTP mus (visible from slope K < 1 throughout the proposed variants in Table III). thank Thomas Handzsuj and Marcus Eger (Drägerwerk AG & Co. KGaA) for many fruitful discussions and technical support.