Annals of Biomedical Engineering, Vol. 36, No. 11, November 2008 (
2008) pp. 1844–1855 DOI: 10.1007/s10439-008-9550-8

Detection of One-Lung Intubation Incidents LIOR WEIZMAN, JOSEPH TABRIKIAN, and ARNON COHEN Department of Electrical and Computer Engineering, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel (Received 20 October 2007; accepted 4 August 2008; published online 22 August 2008)

There are currently several means for OLI detection:

Abstract—Lung sounds are very common source for monitoring and diagnosis of pulmonary function. This approach can be used for detecting one lung intubation (OLI) during anesthesia or intensive care. In this paper, an algorithm for detecting OLI from lung sounds is presented. The algorithm assumes a multiple-input-multiple-output system, in which a multi-dimensional auto-regressive model relates the input (lungs) and the output (recorded sounds). An OLI detector is developed based on the generalized likelihood ratio test (GLRT), assuming coherent distributed sources for each lung. This method exhibited reliable results also when the lungs were modeled by incoherent distributed sources, which is a more accurate model for lung sources. The algorithm was tested using real breathing sounds recorded in an operating room, and it achieved an OLI detection rate of more than 95%, for each breathing cycle.


The stethoscope is the simplest mean, but it requires high attention and its reliability is low.30
The Capnograph, which measures the exhaled amount of CO2, was proved to be an unreliable method for detection of OLI.26
Oxygen saturation, measured by pulse oximeter is one of the most widely used methods today19,27 and the most reliable one, but its results are provided with latency of 2–5 min, after the patient becomes symptomatic without indicating the reason for desaturation. This long latency also may be too long to prevent damage.22 Although the above methods are commonly used to detect OLI during a surgery, none of them are useful for the detection of this incident prior to the development of symptoms (desaturation, CO2 changes, ventilation pressures changes etc.) when there is already an imminent life threatening. There is a need to develop a technique, by which the Anesthesiologist will be alerted once the patient is at risk but still stable. Detection of OLI is requested before any clinical signs appear, in order to enable the Anesthesiologist the opportunity to check the patient and correct the position of the tube prior to the standard practice today. The proposed method for OLI detection is based on analyzing the breathing sound signals. An algorithm for detection of the number of ventilated lungs from the recorded breathing sounds is developed and tested. The proposed algorithm is developed under the assumption that the chest cavity is a linear time-invariant (LTI) system, modeled by a multiple-input-multiple-output auto-regressive (MIMO-AR) system. The decision on the number of ventilated lungs is based on the generalized likelihood ratio test (GLRT). Analysis of respiratory sound signals, in application to diagnosis and monitoring of pulmonological malfunctioning, is well known in the literature.4,5,8,9,12–16,21,23,28–31 Several papers have addressed the problem of modeling of the acoustic transmission of the respiratory system.4,8,16,31 Others dealt with analysis of the

Keywords—MIMO-AR, GLRT, Ventilation, OLI, TRI.

INTRODUCTION Ventilation of patients under general anesthesia or in intensive care is performed by an endotracheal tube, which is placed in the trachea. The location of the tip of tube is critical and it should be placed and maintained above the bifurcation of the trachea. A correct position of the tube, in which both lungs are ventilated, is called tracheal intubation (TRI). If the tube is misplaced or shifted due to patient movements, cases of one lung intubation (OLI) may occur. Prolonged cases of OLI should be avoided since it may cause insufficient oxygenation and may damage the ventilated or the non-ventilated lung. The Australian Incident Monitoring Analysis22 indicates that OLI is the most prevalent incident in anesthesia today. It is more prevalent in obese and pediatric patients (who have relatively short span between the vocal cord and the Carina).

Address correspondence to Joseph Tabrikian, Department of Electrical and Computer Engineering, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel. Electronic mail: [email protected]. ac.il

1844 0090-6964/08/1100-1844/0


2008 Biomedical Engineering Society

Detection of One-Lung Intubation Incidents

breathing sound signals for certain detection or classification purposes.5,13,15 As for detection of OLI from breathing sounds signals, several works have been published.9,12,14,21,23,29 Sod-Moriah et al.21 proposed an OLI algorithm based on linear predictive coding (LPC)11 and partial correlation (PARCOR)3 features. In their study dogs were ventilated while 2 microphones were attached to their chest. In addition, a strain gauge was used to monitor inspiration/expiration/pause periods. The recorded signals were first pre-processed in order to cancel the phonocardiogram (PCG) signal and then classified by extracting LPC and PARCOR features. These features were proved to be unreliable for OLI detection. In addition, classification by energy ratio was also tested. In their study, 89% of OLI detection rate was reported using an energy ratio based classifier. As further explained in this section, the energy ratio method provides poor results when applied to our database, which is based on humans. A pioneering study which uses a sound source to monitor the tube location was published by Juan et al.9 The method is based on analyzing the reflections on an acoustic pulse emitted into the tube in order to estimate the location and the position of the tube. The method was tested on rabbits, and it has the potential to exhibit reliable results while tested on humans. In this paper, we propose a passive method which does not require an external energy source. Another study performed by Tejman-Yarden et al.23 was also based on detection of OLI by monitoring the lungs sounds. The patients’ lungs were ventilated with a double lumen tube, which is a special tube used for lungs surgery. Two microphones were attached to each patient, one on each side of the chest. Another microphone was attached to the patient’s arm, in order to record the environment sound signals for adaptive noise cancellation. A simple algorithm, which was based on energy envelope calculation, breathing sound signals segmentation and energy ratio between samples collected from each microphone, was tested. Additional papers on detection of OLI using energy ratios between left and right breathing sound signals were published.12,14 This energy ratio method is not expected to exhibit reasonable results under every circumstances during a surgery in an operating room, due to mutual coupling between the lungs. This conclusion is also supported by Tejman-Yarden et al.23 In addition, the energy ratio method was tested on our database and exhibited poor results as shown in section ‘‘Experimental Results’’ of this paper. Therefore, a new and reliable approach, based on the technology of monitoring the lungs’ sound signals to detect OLI is presented in this paper. The clinical consequences of

1845

this research are available.24 This paper presents the signal processing aspects of the research, and it includes rigorous derivation of the algorithm which does not appear in our previous publication.24 A discussion about the differences between this paper and our previous publication is further presented in section ‘‘Discussion and Conclusions.’’

OLI DETECTION Model Formulation In the proposed method, the breathing sound signals are recorded by 4 microphones attached to the patient’s back. Each one of the microphones records reverberated sounds generated by both lungs. Accordingly, a convolutive mixture model approach is adopted in this paper. In the proposed method, an auto-regressive (AR) model that relates the lungs and the microphones was assumed. The AR model was chosen because it is commonly used in applications of speech and audio processing and its computational complexity is relatively low. In this model, each ventilated lung represents a distributed source. Our goal is to detect the situation in which only one lung is ventilated, from the received signals by the sensors. It is assumed that the signals generated by the ventilated lungs are independent. Figure 1 shows a block diagram of the proposed MIMO-AR model, in which x[n] represents the sources (lungs), and y[n] represents the sensor (microphones) measurements. Let K and L denote the number of sources (lungs) and sensors (microphones), respectively (K < L). Therefore, the vector of source signals, x[n], is defined as a K 9 1 vector as follows x½n ¼ ½x1 ½n x2 ½n



xK ½nT :

ð1Þ

The L 9 1 measurement vector is defined as y½n ¼ ½y1 ½n

y2 ½n



yL ½nT :

ð2Þ

The relation between the source signals and the measurements is assumed to be given by a MIMO-AR model y½n ¼ AyðMÞ ½n þ Cx½n þ e½n; where yðMÞ ½n is an ML 9 1 vector defined as h iT ðMÞT ðMÞT ðMÞT yðMÞ ½n ¼ y1 ½n y2 ½n    yL ½n

ð3Þ

ð4Þ

ðMÞ

and yi ½n is an M 9 1 vector which contains the past values of the i-th sensor, yi[n]: ðMÞ

yi

½n ¼ ½yi ½n  1 yi ½n  2



yi ½n  MT : ð5Þ

1846

WEIZMAN et al.

FIGURE 1. Block diagram of the MIMO-AR model.

A is an L 9 ML matrix defined as 3 2 T a11    aT1L 6   7 7 6 6 A¼6   7 7; 4   5 aTL1    aTLL

to distinguish between the cases of OLI (K = 1) and TRI (K = 2). The ML Estimator

ð6Þ

where aij is an M 9 1 vector, which relates the samples of the i-th sensor, yi[n], with the past values of the j-th ð MÞ sensor, yi ½n ¼ ½yi ½n  1; . . . ; yi ½n  MT : C is an L 9 K matrix whose i,j-th element relates the samples of source j and sensor i. Finally, e[n] is an L 9 1 vector representing additive white noise. It is assumed that the noise and the source signals are independent, zero-mean, Gaussian with covariance matrices r2I and I, respectively. These assumptions may cause improper modeling of the problem and lead to suboptimal results, but they enable us to deal with a simple formulation which leads to reasonable computation time for future real-time implementation. The last assumption can be employed with no loss of generality, because the covariance of the sources is determined by the unknown matrix C, as it can clearly be seen from (3). As a result, it
is obtained that the conditional distribution of y½n
yðMÞ ½n is Gaussian:  
y½n
yðMÞ ½n  N AyðMÞ ½n; R ; where R is defined as

Estimation of the unknown matrices, A and R, from the N samples of the data: y[1], …, y[N], shall be performed prior to determination of the number of sources, K. For this purpose, the maximum-likelihood (ML) estimator is used. The ML estimator of the matrices A and R, is obtained by maximizing the logarithm of the conditional probability density function (PDF) of the output samples given the unknown matrices, which is26 NL N logð2pÞ  logjRj log fðy½1; . . . ; y½NjR; AÞ ¼  2 2 N     X T 1  y½n  AyðMÞ ½n R1 y½n  AyðMÞ ½n : 2 n¼1 ð8Þ The ML estimator of the matrix A is obtained by maximization of (8) with respect to A: ^ ML ¼ arg max log fðy½1; . . . ; y½NjR; AÞ: A A

In Appendix A it is shown that ! !1 N N X X ðMÞT ð MÞ ðMÞT ^ y½ny ½n y ½ny ½n : AML ¼ n¼1

R ¼ CCT þ r2 I:

n¼1

ð7Þ

ð9Þ

Note that the unknown parameters: A, R, M and K must be estimated from a set of N measurements, y[1], …, y[N]. It is also assumed that all the initial conditions are zero, i.e., y[n] = 0 for n < 0, and that the input and noise signals are stationary. Our goal is

Note that this estimation is independent of R. In order to maximize the log-likelihood function with respect to R, we shall define ^ ML yðMÞ ½n; z½n ¼ y½n  A

ð10Þ

Detection of One-Lung Intubation Incidents

and (8) turns into   ^ ML log f y½1; . . . ; y½NjR; A ¼

N NL N 1X logð2pÞ  logjRj  zT ½nR1 z½n: 2 2 2 n¼1

ð11Þ The sum term in (11) is a scalar, and therefore N X

zT ½nR1 z½n ¼

n¼1

N X     tr R1 z½nzT ½n ¼ N tr R1 S ; n¼1

ð12Þ where tr(Æ) denotes the trace operation, and S denotes the sample covariance matrix of the residual process z[n]: N N 1X 1X S¼ z½nzT ½n ¼ y½nyT ½n N n¼1 N n¼1 ! !1 N N X 1 X ðMÞT ðMÞ ðMÞT  y½ny ½n y ½ny ½n N n¼1 n¼1 ! N X ðMÞ T  y ½ny ½n : ð13Þ n¼1

1847

^ R ¼ US maximized with respect to UR by choosing U ^ T US ¼ UT U ^ R ¼ I6 : Accordingly, maximiand thus U R S zation of (16) with respect to UR gives   ^ ML max log f y½1; . . . ; y½NjR; A UR " # L X  1  N L logð2pÞ þ log kRi þ tr KR KS ¼ 2 i¼1 " # L L X X N kSi L logð2pÞ þ ¼ log kRi þ 2 k i¼1 i¼1 Ri " K L X X N ¼ log kRi þ log r2 L logð2pÞ þ 2 i¼1 i¼Kþ1 # K L X X kS i kS i þ þ ; ð17Þ k r2 i¼1 Ri i¼Kþ1 where fkRi gLi¼Kþ1 are constrained to be r2. The ML estimator of fkRi gK i¼1 can be derived by equating the partial derivatives of (17) with respect to kRi ; 1 £ i £ K, to zero. As a result, it is obtained that  ^kR ¼ kSi ; i ¼ 1; . . . ; K; i r2 ; i ¼ K þ 1; . . . ; L: Therefore, the ML estimator of R under the constraint that its lower L-K eigenvalues are equal to r2, is

The singular value decomposition (SVD) of the matrix S is

^ RU ^ RK ^ T; ^ ML ¼ U R R

S ¼ US KS UTS ;

^ R is a diagonal matrix whose ^ R ¼ US ; and K where U diagonal elements are given by kS1 ; . . . ; kSK ; r2 ; . . . ; r2 : |fflfflfflfflfflffl{zfflfflfflfflfflffl}

ð14Þ

in which KS is a diagonal matrix with diagonal elements kS1  kS2     kSL representing the eigenvalues of S, and US is the matrix of eigenvectors. Similarly, the SVD of R is given by R ¼ UR KR UTR ;

ð15Þ

where KR is a diagonal matrix with diagonal elements kR1 kR2      kRL representing the eigenvalues of R, and UR is the matrix of eigenvectors. As a result, (11) can be rewritten as   ^ ML  NL logð2pÞ log f y½1; . . . ; y½NjR; A 2 ! L Y   N N T T  log kRi  tr UR K1 R UR US KS US 2 2 i¼1 L NL NX logð2pÞ  log kRi 2 2 i¼1  T   T T  N   tr K1 : R UR US KS UR US 2

¼

ð16Þ

The matrix R, defined in (7), has K eigenvalues greater than r2, and L-K eigenvalues equal to r2, because the rank of CCT is equal to K. Equation (16) is

ð18Þ

LK

In the presence of white noise with variance, r2, the eigenvalues of the matrix R are greater or equal to r2. Harmanci et al.6 showed that the constrained ML estimator of the matrix R is also given by (18) with ^ R ¼ US and eigenvalues U  2 ^kR ¼ maxðkSi ; r Þ; i ¼ 1; . . . ; K; ð19Þ i 2 r; i ¼ K þ 1; . . . ; L:

GLRT Model order selection methods based on information theoretic criteria1–3,17,20 can be a possible approach in order to determine the AR model order, M, and the number of sources, K. This method was developed and tested during this work, but due to modeling mismatches its performance was found to be unreliable when applied to real breathing sound signals.28 Therefore, a GLRT-based method10 was developed and tested as presented in this section. In the OLI detection problem, the lungs represent the source signals and thus, the number of sources can

1848

WEIZMAN et al.

be either one or two. Thus, for the purpose of decision between TRI and OLI cases, the GLRT is used. This test is based on the ratio between the PDF’s under each hypothesis. Let us denote the following hypotheses: H 1:

and since from (19) we have kR^2 ¼ maxðkS2 ; r2 Þ  r2 ; then the logarithm term in (24) can be removed and the GLRT can be rewritten as: H2

kS2

Only one source exists—OLI case, K(H1) = 1,

H2: There are two independent sources—TRI case, K(H2) = 2. Then, the GLRT is given by:   ^ ML ; H1 H1 ^ ð1Þ ; A f y½1; . . . ; y½NjR ML   > c; ð20Þ ð2Þ ^ ^ f y½1; . . . ; y½NjRML ; AML ; H2 < H2

^ ð1Þ ; R ML

^ ð2Þ R ML

where are the ML estimators of R under hypotheses H1 and H2, respectively, and c is a threshold which is determined according to the trade-off between false and miss detections. The log-likelihood function under the i-th hypothesis can be obtained by substituting the L-K lowest eigenvalues from (19) into (17):

i

Thus, the likelihood ratio for decision between H1 and H2 is as follows:   ^ ML ; H1 ^ ð1Þ ; A f y½1; . . . ; y½NjR ML  log  ð2Þ ^ ^ f y½1; . . . ; y½NjRML ; A ML ; H2 " # kR^1 kR^2 r2ðL2Þ kR^2 N log ¼  2 þ1 2 r kR^1 r2ðL1Þ   k^ k^ N ¼ log R22  R22 þ 1 : ð22Þ 2 r r By taking the logarithm of the GLRT in (20) and using (22), one obtains H1

k^ k^ log R22  R22 r r

> <

2 log c  1; N

ð23Þ

H2

or

<

c00 :

ð25Þ

H1

It can be observed from (25), that the algorithm for OLI detection is simple to implement and it is given by thresholding the second highest eigenvalue of S defined in (13). The threshold c¢¢ is set according to the allowed false alarm rate. Implementation of the algorithm involves computation of the matrix S, defined in (13), and its eigenvalues. The matrix S can be computed as follows: T S ¼ C1 þ C2 C1 3 C2 where: C1 ¼

N 1X y½nyT ½n; N n¼1

and C3 ¼

N X

C2 ¼

N 1X y½nyðMÞT ½n; N n¼1

yðMÞ ½nyðMÞT ½n:

n¼1

ðiÞ

^ ML ; Hi Þ ^ ;A log fðy½1; . . . ; y½NjR ML " ! ! KðH Yi Þ N N 2ðLKðHi ÞÞ kR^i r ¼ L logð2pÞ þ log 2 2 l¼1 # L X kR^l þ KðHi Þ þ : ð21Þ r2 l¼KðH Þþ1

>

The orders of the computational complexity of C1 ; are: C1 : NLðLþ1Þ ; C2 : NL2 M; C2 ; C3 and C1 3 2 3 3 C3 : NLMðLMþ1Þ and C1 3 : L M and the computation 2 T is of order: complexity of S ¼ C1 þ C2 C1 3 C2 2 3 3 2 2M L þ 2ML þ L : |fflfflffl{zfflfflffl} |fflffl{zfflffl} |{z} C2 C1 3

T C2 C1 3 C2

T C1 þC2 C1 3 C2

In addition, the eigenvalues computation has the complexity of L3. Under the assumption that M  L, the computational complexity of S, considering  one frame with N samples is: O L2 M2 LM þ N2 : Considering the parameters in our real data experiments described in section ‘‘Experimental Results’’ (i.e., L = 4, M = 15 and N = 2000), and a frames overlapping value of 80%, the time complexity would be ~19 MIPS (Million Instructions Per Second), which requires a modest processing power for real data implementations. EXPERIMENTAL RESULTS In this section, the performance of the proposed OLI detection algorithm is tested with simulated and real recorded data.

H2

kR^2 k^  log R22 2 r r

> <

c0 ;

ð24Þ

H1

where c0 ¼  N2 log c þ 1: Since the function f(x) = xlog x is a monotonically increasing function for x > 1,

Simulation Results The OLI algorithm derived in this paper was developed under the assumption that each lung is modeled by a coherently distributed source. In practice, this model is not accurate, and a more accurate

Detection of One-Lung Intubation Incidents

model is the incoherently distributed source model. In this section, the performance of the proposed algorithm was tested using simulated data with both coherently and incoherently distributed sources.

0.5

1849

Pe of estimating number of sources, K , σ 2=1

0.45 0.4 0.35

Coherent Sources In this test, the MIMO-AR system defined in (3), was simulated and the ML estimators A of and R were calculated according to (9) and (18). In the simulations the system parameters were chosen as follows (unless otherwise is indicated): the number of sources was K = 2, the number of sensors was L = 4 with AR order M = 5, and the noise variance was r2 = 1. In order to control the angle between the columns of the matrix C, we chose the columns of C as follows: 2 1 ffi3 0 pffiffiffiffiffiffiffi c2 þ3 6 0 pffiffiffiffiffiffiffi 1 ffi7 6 c2 þ3 7 C¼6 : 1 ffi7 4 0 pffiffiffiffiffiffiffi 5 c2 þ3 c ffi 1 pffiffiffiffiffiffiffi c2 þ3

Pe

0.3 0.25 0.2 0.15 0.1 0.05 0 2 10

3

10

Number of samples - N

FIGURE 2. The probability of error, Pe, of the proposed OLI detector with simulated coherent distributed source, with u 5 90.

It can clearly be seen that the columns of C have unit norm. In addition, the angle between the columns  of C,  c ffi u, is controlled by c, according to: u ¼ arccos pffiffiffiffiffiffiffi : 2 c þ3

The matrix A was randomly chosen, while ensuring it is a stable matrix, representing system poles inside the unit circle in the Z plane. Simulation results of the ML estimators of the matrices A and R are available.27 Note that the values for the matrices A and R were chosen for simulation purposes. In the actual experiments, which will follow, these matrices are estimated as explained in the previous section. The performance of the GLRT as a function of the number of independent samples, N, is examined. The second highest eigenvalue, kS2 ; was extracted and compared to a threshold value, c¢¢, which was set to the noise level r2. Under each one of the hypotheses (K = 1, 2), and for each N, J = 500 iterations were performed. The probability of error for determining K is defined as: number of events of incorrectly estimated K : Pe ¼ J ð26Þ Figure 2 shows the probability of error, Pe, as a function of the number of samples, N, while the columns of C are orthogonal (u = 90). It can be seen that the probability of error decreases as the number of samples grows. Figure 3 shows Pe as a function of the angle between the columns of the matrix C, while u varies between 9 and 90. Based on Fig. 2, the number of samples, N, was set to 250, in order to exhibit a meaningful examination of the behavior of Pe as a function of u. It can clearly be seen that Pe decreases as

FIGURE 3. The probability of error, Pe, of the proposed OLI detector with simulated coherent distributed source, as a function of the angle between the columns of C, with N 5 250.

the angle between the columns of C grows. Since the angle between the columns of C represents the distance between the spatial signatures of the sources, the estimation of the number of sources is improved as they become orthogonal. Incoherent Distributed Sources Each lung can be modeled by a distributed source, which can be approximated by several independent point sources. In order to evaluate the performance of the algorithm under this assumption, two spatially distributed sources representing the two lungs were synthesized. Each distributed source was composed of four independent point sources with close spatial

1850

WEIZMAN et al.

FIGURE 4. The eigenvalues of S vs. the scattering level, e.

signatures. Therefore, the columns of C were chosen to be: " C ¼ c1 c1 þ eDc11 c1 þ 2eDc12 c1 þ 3eDc13 ; |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} source 1

#T

c2 c2 þ 2eDc21 c2 þ 3eDc22 c2 þ eDc23 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} source 2

ð27Þ where fci Dci1 Dci2 Dci3 g is an orthonormal group. e is a constant which controls the spatial distribution width of each source. The product Cx[n] represents the contribution of the two distributed sources, where each one is composed of 4 independent point sources. In order to examine the performance of the proposed GLRT under the incoherently distributed sources condition, the eigenvalues of S for the cases of one and two distributed sources are computed. In Fig. 4 the eigenvalues of S are drawn as a function of e for K = 1 and K = 2. It can be seen from Fig. 4 that even when the sources are widely distributed, the eigenvalues of the source signal subspace are well separated from the eigenvalues of the noise subspace. Therefore, despite of the fact that each lung cannot be modeled by a point source or coherently distributed sources, the second highest eigenvalue of S was chosen as a

detector for OLI situation in real breathing sound signals. In addition, it can also be seen that the threshold value in the case where the sources are widely distributed should be higher than the threshold in the case of point sources. Real Data Results In order to examine the proposed algorithm for OLI detection, a database of recorded breathings was established. The database was composed of 24 adult surgical patients (ASA I and II), which were recorded in an operating room in both situations: during correct ventilation, when the tip of the tube is placed above the carina, and during a situation of OLI when the tip of the tube is below the carina and only one lung is ventilated. The study was approved by the ethics committee of the Soroka Medical Center in BeerSheva, Israel. The patients signed informed consent forms prior to participation in the study. All of the patients were scheduled for a routine surgical abdominal or orthopedic procedure that required the insertion of a single lumen endotracheal tube. No patients of open chest surgeries were participated in this study. Prior to anesthesia, four Tommyscope piezoelectric acoustic sensors (KOL Medical Ltd., Innovative Medical Equipment, Ramat-Gan, Israel) were attached to the patient’s back, as shown in Fig. 5. After induction of anesthesia with thiopental and succinylcholine,7 the

Detection of One-Lung Intubation Incidents

tube was inserted and advanced down the airway into the right or left mainstem bronchus under fiberoptic guidance, and six breaths were then delivered via positive-pressure bag ventilation. The tube was then withdrawn stepwise until a correct position of the tube is reached, and six similar breaths were delivered again. The experiments were performed in the main operating room of Soroka University Medical Center—Israel. Because of the fact that most of the information in the breathing sound signals exist in low frequencies, the data

FIGURE 5. Experimental set-up—view from back (TRI case).

1851

were recorded at 4 kHz sampling rate, using a 12 bit A/D (model: KPCMCIA-12AI-C, Keithley, Cleveland, Ohio, USA). The experimental setting is identical to those published in our clinical paper.24 The data were filtered with bandpass Butterworth filter type with cut-off frequencies29 of 100 and 600, downsampled with factor of 0.3, and divided into windows of 2000 samples each, with 80% of overlapping. An arbitrary AR order of 15 was set considering the computation complexity and the available processing time. The recorded breathing signals contain both situations of OLI and TRI. Figure 6 shows a few breathing cycles of both OLI and TRI situations, recorded by the four microphones after preprocessing. As it can be seen from this figure, the OLI and TRI cases cannot be identified by using only the amplitude of the recorded sounds. This statement is supported further in this paper when we compare the performance of our proposed method with the energy ratio method. Figure 7 shows the GLRT stated in (25), as a function of time using the measurements shown in Fig. 6. Figure 8 shows the histograms of the second highest eigenvalue of S for every breathing cycle, in both OLI and TRI cases, along the 24 experiments. As it can clearly be seen from Figs. 7 to 8, OLI and TRI

FIGURE 6. Examples of the recorded signals—four microphones, OLI and TRI cases.

1852

WEIZMAN et al.

FIGURE 7. Second highest eigenvalue of S as a function of time.

FIGURE 8. Histograms of the second highest eigenvalue of S over 24 experiments for OLI (left curve) and TRI (right curve) cases.

cases can clearly be discriminated by the GLRT given by the second highest eigenvalue of S in every breathing cycle. Evaluation of the performance of the system was carried out using the ‘‘leave some out’’ method, which is an extension of the ‘‘leave one out’’ method25 as follows. Twenty different experiments were used to extract the data for training the system. The rest of the 4 experiments were used to validate the system and were tested according to the extracted statistics in the training process. This process was repeated 6 times. Each time a different group of four validation experiment was used and the results were validated using a total number of 24 experiments, in a ‘‘patient-independent’’ mode. There are two types of errors in OLI detection: Pmiss, the probability of a true OLI to be wrongly detected as TRI, and PFA, the probability of TRI to be

FIGURE 9. The DET curve of energy ratio OLI detector.

FIGURE 10. The DET curve of the proposed OLI detector.

detected as OLI. The detection error tradeoff (DET) curve is a common mean to display the trade-off between these errors. The DET curve provides information about the device’s performance, where each point on the curve shows the PFA and Pmiss for a given threshold. The threshold of a real monitoring system should be calculated according to the required sensitivity of the system, while taking into consideration the allowed Pmiss of the system. The DET curve in our experiments computed using a total number of 24 experiments, according to the ‘‘leave some out’’ method described above, based on a single breathing cycle. In order to examine the results of the energy ratio method on our database, an OLI detector based on this method12,14 was also applied to our database. Figure 9 presents the DET curve of this detector, and Fig. 10 shows the DET curve of our proposed decision system. The equal error rate (EER) point is defined as the point on the DET curve where Pmiss = PFA.

Detection of One-Lung Intubation Incidents

The improvement of our method with relation to the energy ratio method can clearly be noticed by comparing the EER point for the energy ratio detector, which is 34%, with the EER of our method, which is 4.8%. Naturally, more importance should be given to Pmiss rather than to PFA. Therefore, it is assumed that in a practical system the selected activity point on a DET curve based on our method will be at Pmiss = 2% and PFA = 9%.

DISCUSSION AND CONCLUSIONS In this paper, a new concept for automatic OLI detection system is presented. An algorithm for detection of OLI by monitoring lung sounds is developed, based on assumption of Gaussian distribution of the signals. The Gaussian assumption may lead to suboptimal results, since it ignores the higher-order statistics. Therefore, a relatively simple detector was developed, which is more suitable for developing a real-time monitoring system. Avoiding the Gaussian assumption for the breathing sound signals will probably lead to a better system performance,18 but the computation complexity may not meet the real-time system requirements. Note that the OLI detector of the proposed system, i.e., the second highest eigenvalue of the residual covariance matrix, has also an acoustic meaning. The eigenvalues of the residual covariance matrix represent the energy of the input sources. Under the assumption that the lungs signals are stronger than the noise signals, in the TRI case the energy of the lungs is spread more or less equally over the two highest eigenvalues. During cases of OLI, there is only one major energy source (the ventilated lung), and as a result, the second highest eigenvalue decreases dramatically. Therefore, the second highest eigenvalue was found to be an indicator to the existence of an additional acoustic source. In our experiment, the method was implemented using four microphones. From our point of view, applying the proposed method with less than four microphones may degrade the performance due to improper estimation of the noise level. In addition, the usage of four microphones enables to exhibit a reliable method even in a presence of background noise. In such cases, most of the energy that is originated from the noise source will be expressed in the third and fourth highest eigenvalues of R, and the decision about OLI, which is based on the second highest eigenvalue, will not be affected. Implementation of the presented method into a real-time monitoring system will bring to a major improvement in the ability to detect a ventilation incident as it occurs, while the patient is still well oxygenated and stable. However, few points should be

1853

considered prior to developing a real-time system. The issue of setting the threshold level for the decision system also needs to be taken into consideration while designing a real-time monitoring system. Due to the physiological variability between patients, we suggest that the threshold level will be automatically set for every patient after a short training over a validated short TRI period in the beginning of the operation. From the clinical point of view, since this method analyses the patient’s ventilation status rather than endotracheal tube position, it may be possible to use the method for detection of other surgical life threatening events that lead to unequal breath sounds, such as esophageal intubation, massive unilateral atelectasis or pneumothorax. Moreover, since the system is attached to the patients’ back in parallel position to the ECG monitoring system, it could be designed in the future, to combine these two monitors into one system, with both modalities sampled by the same sensors.23 Note that the proposed method is not meant to replace the anesthesiologist during the intubation process. However, the main contribution of the proposed method is that it can monitor the tube’s location during the surgery (while the anesthesiologist may be busy with other tasks). In addition, it reduces the anesthesiologist need to rely on other means, which are not directly meant for OLI purposes. Another clinical issue which should be considered is the effect of microphones locations. After performing several experiments, we concluded that the microphones should be attached to the patient’s back (as illustrated in Fig. 5), rather than the patient’s chest, in order to attenuate the background noise. However, this setting of the microphones might lead to pressure necrosis, particularly during long procedures. Although this complication has never occurred in our experiments, new lightweight and small microphones (S.L.P, Tel-Aviv, Israel), designed especially for this purpose, will be used in future experiments in order to prevent such a phenomenon. However, the results of our method should be limited to the conditions of the experiments that were considered in our work, i.e., detection of OLI under operating room conditions with no open chest surgeries. Detection of OLI in Intensive Care Unit (ICU) or during open chest surgeries may not meet the requirements of our method. As mentioned above, the clinical aspects of our work have been published.24 The major contribution of this paper is from the signal processing point of view of the problem. This paper includes mathematical derivations of the ML estimator and the GLRT detector, along with a computation of the complexity of the problem, for real time implementation purposes. In addition, an analysis of the performance of the model with respect to the

WEIZMAN et al.

1854

mixing matrix (C) and the scattering level of the sources is presented. We also present a comparison between the performances of our method to the energy ratio detection method. In addition, the analysis presented in this paper enables the reader to continue the research in this area, in order to further develop a reliable real-time monitoring system. In conclusion, we established a database of recorded breathing sound signals of patients during OLI and TRI situations. It has been shown that a MIMO-AR model and thresholding the second highest eigenvalue of the residual covariance matrix provides a reliable method for OLI detection on the recorded breathing sound signals. Additional work for developing a realtime monitoring system in a future research based on the propose method is currently in progress.

APPENDIX A: PROOF OF (9) Maximization of (8) with respect to the unknown matrix, A, is achieved via equating the corresponding partial derivative to zero. The last term of (8) is relevant to calculate the derivative of the log-likelihood with respect to A: N  T   1X  y½n  AyðMÞ ½n R1 y½n  AyðMÞ ½n 2 n¼1 " N 1X ¼ yT ½nR1 y½n  yðMÞT ½nAT R1 y½n 2 n¼1 #  T ð MÞ ðMÞ ðMÞ T 1 1  y ½nR Ay ½n þ Ay ½n R Ay ½n : ð28Þ The derivative of a scalar c with respect to a matrix A is as a matrix whose ij-th element is given by  @c defined @c ¼ : We shall use the following identities (for @A ij @Aij proof of non-trivial identities, see Makhoul11): @  T @  T T x Ay ¼ y A x ¼ xyT ; @A @A

ð29Þ

i @ h ðAxÞT CðAxÞ ¼ ðC þ CT ÞðAxÞxT : @A

ð30Þ

The derivative of (28) w.r.t. A is calculated using (29) and (30), and thus, the derivative of (8) can be written as @ logðfðy½1; . . . ; y½NjR; AÞ @A N h 1X ¼ R1 y½nyðMÞT ½n  RT y½nyðMÞT ½n 2 n¼1 i   þ R1 þ RT AyðMÞ ½nyðMÞT ½n ð31Þ

In order to find the ML estimator of A, the above derivative should be equated to zero. Since R is a covariance matrix, then R ¼ RT and R1 ¼ RT : Therefore we obtain the following equation: 

N h i 1X ^ ML yðMÞ ½nyðMÞT ½n 2R1 y½nyðMÞT ½n þ 2R1 A 2 n¼1

¼ 0:

ð32Þ

Equation (9) can be obtained by solving (32) for ^ ML under the assumption that the matrix A PN ðMÞ ½nyðMÞT ½n is invertible. ( n¼1 y

ACKNOWLEDGMENTS The authors wish to thank Prof. G. Gurman, Dr. S. Teiman and Dr. A. Zlotnik from the Soroka University Medical Center—Israel, for their contribution with collecting the real data and supporting the medical part of the paper. While completing this paper, Professor Arnon Cohen, head of the signal processing laboratory in Ben-Gurion University had passed away. By publishing this paper we continue Prof. Cohen’s legacy of honest and fair research.

REFERENCES 1

Akiake, H. Information theory and an extension of the maximum likelihood principle. In: Proc. 2nd Int. Symp. Inform. Theory, 1973, pp. 267–281. 2 Akiake, H. A new look at the statistical model identification. IEEE Trans. Automat. Contr. 19(6):716–723, 1974. doi:10.1109/TAC.1974.1100705. 3 Cohen, A. Biomedical Signal Processing. FL: CRC Press, 1986. 4 Cohen, A., and A. Berstein. Acoustic transmission of the respiratory system using speech stimulation. IEEE Trans. Bio-Med. Eng. 38(2):126–132, 1991. doi:10.1109/10.76377. 5 Cohen, A., and D. Landsberg. Analysis and automatic classification of breathing sounds. IEEE Trans. Bio-Med. Eng. 31(9):585–590, 1984. doi:10.1109/TBME.1984.325301. 6 Harmanci, K., J. Tabrikian, and J. L. Krolik. Relationships between adaptive minimum-variance beamforming and optimal source localization. IEEE Trans. Signal Proces. 48(1):1–13, 2000. doi:10.1109/78.815474. 7 Heinonen, J., S. Takki, and T. Tammisto. Effect of the trendelenburg tilt and other procedures on the position of endotracheal tubes. Lancet 1:850–853, 1969. doi:10.1016/ S0140-6736(69)91898-4. 8 Iyer, V. K., P. A. Rammoorthy, and Y. Ploysongsang. Autoregressive modeling of lung sounds: characterization of source and transmission. IEEE Trans. Bio-Med. Eng. 36(11):1133–1137, 1989. doi:10.1109/10.40821. 9 Juan, E. J., J. P. Mansfield, and G. R. Wodica. Miniature acoustic guidance system for endotracheal tubes. IEEE Trans. Bio-Med. Eng. 49(6):584–596, 2002. doi:10.1109/ TBME.2002.1001973.

Detection of One-Lung Intubation Incidents 10

Kay, S. Fundamentals of Statistical Signal Processing, Vol. II. Prentice Hall, 1998. 11 Makhoul, J. Linear prediction: a tutorial review. Proc. IEEE 63(4):561–580, 1975. doi:10.1109/PROC.1975.9792. 12 Mansy, H. A., C. J. O’Connor, R. A. Balk, and R. H. Sandler. Breath sound changes associated with malpositioned endotracheal tubes. Med. Bio. Eng. Comp. 43:206– 211, 2005. doi:10.1007/BF02345956. 13 Mansy, H. A., T. J. Royston, R. A. Balk, and R. H. Sandler. Pneumothorax detection using computerized analysis of breath sounds. Med. Bio. Eng. Comp. 40:529–532, 2002. 14 O’Connor, C. J., H. A. Mansy, R. A. Balk, K. J. Tuman, and R. H. Sandler. Identification of endotracheal tube malpositions using computerized analysis of breath sounds via electronic stethoscopes. Anesth. Analg. 101:735–739, 2005. doi:10.1213/01.ane.0000167068.71601.e4. 15 Ploy-Song-Sang, Y., P. T. Macklen, and W. R. D. Ross. Distribution of regional ventilation measured by breath sounds. Am. Rev. Respir. Dis. 117:657–664, 1978. 16 Ra¨sa¨nen, J. O., G. Rosenhouse, and N. Gavrieli. Effects of diameter, length and circuit pressure on sound conductance through endotracheal tube. IEEE Trans. Bio-Med. Eng. 53(7):1255–1264, 2006. doi:10.1109/TBME.2006.875667. 17 Rissanen, J. Modeling by shortest data description. Automatica 14:465–471, 1978. doi:10.1016/0005-1098(78)90005-5. 18 Routtenberg, T., and J. Tabrikain. MIMO-AR system identification and blind source separation using GMM. To appear in: Proc. of 32nd Intl. Conf. on Acoustics, Speech and Signal Processing (ICASSP), 2007. 19 Runciman, W. B., R. K. Webb, L. Barker, and M. Currie. The pulse oximeter: applications and limitations—An analysis of 2000 incident reports. Anesth. Intensive Care 21:543–550, 1993. 20 Schwartz, G. Estimating the dimension of a model. Ann. Stat. 6:461–464, 1978. doi:10.1214/aos/1176344136. 21 Sod-Moriah, G., O. Gelber, G. Gurman, and A. Cohen. Monitoring of separate lung ventilation in anesthesia and intensive care. In: 18th Convention of Electrical and Electronics Engineers in Israel, 1995 (Available at: http:// ieeexplore.ieee.org/servlet/opac?punumber=3899).

22

1855

Szekely, S. M., R. K. Webb, J. A. Williamson, and W. J. Russell. The Australian Incident Monitoring Study Problems related to the endotracheal tube: an analysis of 2000 incident reports. Anesth. Intensive Care 21:611–616, 1993. 23 Tejman-Yarden, S., D. Lederman, I. Eilig, A. Zlotnik, N. Weksler, A. Cohen, and G. M. Gurman. Acoustic monitoring of double-lumen ventilated lungs for the detection of selective unilateral lung ventilation. Anesth. Analg. 103:1489–1493, 2006. doi:10.1213/01.ane.0000240909. 48774.49. 24 Tejman-Yarden, S., A. Zlotnik, L. Weizman, J. Tabrikian, A. Cohen, N. Weksler, and G. M. Gurman. Acoustic monitoring of lung sounds for the detection of one lung intubation. Anesth. Analg. 105:397–404, 2007. doi:10.1213/ 01.ane.0000267263.40505.c0. 25 Wahba, G., and S. Wold. A completely automatic French curve: Fitting spline by functions by cross-validation. Commun. Stat. 4:1–17, 1975. 26 Webb, R. K., J. H. van der Walt, W. B. Runciman, J. A. Williamson, et al. Which monitor? An analysis of 2000 incident reports. Anesth. Intensive Care 21:529–542, 1993. 27 Webster, T. A. Now that we have pulse oximeters and capnographs, we don’t need precordial and esophageal stethoscope. J. Clin. Monitor. 3:191–192, 1987. doi:10.1007/ BF01695942. 28 Weizman, L. Detection of One-Lung Intubation Incidents, M.Sc. Thesis. Ben-Gurion University of the Negev, Dept. of Electrical and Computer Engineering, Beer-Sheva, Israel, 2004. (Available at: http://www.ee.bgu.ac.il/~spl/ publication/theses.htm). 29 Weizman L., J. Tabrikian, and A. Cohen. Detection of one lung intubation by monitoring luns sounds. In: Proc. of the 24th Annual International Conference of the IEEE EMBC. San-Francisco, CA, USA, pp. 917–920, 2004. 30 Wodicka, G. R., P. D. DeFrain, and S. S. Kraman. Bilateral asymmetry of respiratory acoustic transmission. Med. Bio. Eng. Comp. 32(5):489–494, 1994. doi:10.1007/BF02515306. 31 Wodicka, G. R., and H. L. Golub. A model of acoustic transmission in the respiratory system. IEEE Trans. BioMed. Eng. 36(9):925–934, 1989. doi:10.1109/10.35301.