Anal. Chem. 1994,66, 3066-3072
Deconvolution of Overlapping Skewed Peaks Attila Fellnger Department of Analytical Chemistry, University of Veszprhm, Veszprhm, P.O. Box 158, H-820 1, Hungary
Asymmetrical chromatographicpeaks can be well modeled by the exponentially modified Gaussian model (EMG). In this study, an analytical expression is derived to give the model for the sharpened overlapping peaks when the deconvolution is performed by another EMG peak as a sharpening function. By means of this model, it can easily be explained what kind of peak shapes can be expected after deconvolution. The validity of the model derived was tested by numerical deconvolution methods-the relaxation and the Fourier deconvolution. The model derived here proves that it is impossible to conserve the peak shape during deconvolution. The integration of the peak areas is almost always burdened with error. The model gives an explanation of why side lobes and spurious peaks emerge during the deconvolution procedure. Fused peaks cause severe difficulties in chromatography and other separations. When analyzing a multicomponent mixture, we have to face that fact that only a small fraction of the components will appear as isolated in spite of having extremely efficient columns. For this reason, different deconvolution techniques have been elaborated with a diverse sophistication. Besides these algorithms, routine chromatographic integrators still apply the simplest graphical methods as perpendicular drop or tangent skim,4 which are known to fail at highly skewed peaks. Numerical deconvolution methods for resolution of overlapping peaks that do not assume the a priori knowledge of the peakshape have been around since the advent of computers. The two major kinds of these deconvolutions are the direct deconvolution in Fourier domain and the iterative relaxation method. The deconvolution of band broadening effects was extensively studied during the course of the early 1970s utilizing the fast Fourier transform (FFT) algorithm. The theoretical basis of frequency domain peak sharpening was developed by Kirmse and Westerberg5for symmetrical peaks. They showed that overlapping Gaussian peaks conserve their shape during deconvolution, provided that the transfer function is also Gaussian. The variance of the sharpened peak is the difference between the original variance and that of the transfer function uf2 = a* - Ud2. However, the deconvolution of skewed peaks is still unexplained, while it is obvious that a Gaussian sharpening function fails to eliminate peak tailing5 The deconvolution of skewed peaks has only been studied numerically up to now6 because of the lack of an analytical closedform model for sharpened skewed peaks. (1) Davis, J. M.; Giddings, J. C. Anal. Chem. 1983,55, 418.
(2) Pietrogrande, M. C.; Dondi, F.; Felinger, A.; Davis, J. M.; Chemom. Intell. Lab. Syst., submitted for publication. (3) Felinger, A.; Pasti, L.;Dondi, F. Anal. Chem. 1990,62, 1846. (4) Papas, A. N.; Tougas, T. P. Anal. Chem. 1990, 62, 234. (5) Kirmse, D. W.; Westerberg, A. W. Anal. Chem. 1971,43, 1035.
3066 Analytical Cbmlstty, Vol. 66, No. 19, October 1, 1994
One of the most critical steps of frequency domain deconvolution is to determine the transfer function. It can be a nonretained peak that is affected by most of the band broadening effect^.^ The total extracolumn broadening effects can be eliminated by removing the column and connecting the injector directly to the detector.* This response function can also be used as a transfer function. Major disadvantages of Fourier domain deconvolution are the sensitivity to noise, the difficulties arising from the determination of the cutoff frequency, and the smoothing window .9J0 A more stable and powerful deconvolution method, the constrained iterative relaxation method-or Jansson's method"-was unknown to chromatographers until the late 198Os.l2J3 This method can easily use such constraints as maximum peak height and nonnegativity. Moreover, it is able to use a sharpening function whose width depends on location; in this way, more band broadening effect can be deconvoluted at higher retention times.14
THEORY EMG Function. In analytical chromatography, peakshapes are usually assumed to be symmetrical. In fact, the theory of linear chromatography demonstrates that the elution profile is very close to the Gaussian function, even if the column is ineffi~ient.'~Unfortunately, many disturbing effects accompany chromatographic separations, which leads to distorted peak profiles. A great number of extracolumn distortion effects can becharacterized by the timeconstant of the transfer function of a first-order linear system, assuming dead volumelike behavior. Surface heterogeneity can also result in highly asymmetric peak shape.15J6 A popular model for describing asymmetrical chromatographic peaks-and maybe the only one that has physical justification-is the exponentially modified Gaussian (EMG) function. Its use arises from the fact that the ideal Gaussian peak is distorted by first-order delays caused by, for instance, dead volumes of the injector or detector volume, etc. It was (6) Johansson, M.;Berglund. M.; Baxter, D. C. Spectrochim. Acta 1993,48B,
1393.
(7) Maldacker, T. A.; Davis, J. E.; Rogers, L. B. Anal. Chem. 1974,46,637. (8) Wright, N. A.; Villalanti, D. C.; Burke, M . F. Anal. Chem. 1982,54,1735. (9) Kiillik, E.; Kaljurand, M.; Ess, L. J . Chromatogr. 1976,118,313. (10) Felinger, A.; Pap, T.L.; InczCdy, J. Anal. Chim. Acta 1991,248,441. (1 1) Jansson, P. A. Deconuolution With Applicationr in Spectroscopy; Academic Press, Inc.: Orlando, FL, 1984. (1 2) Crilly, P. B. J . Chemom. 1987,1, 79. (13) Schure, M. R. J. Chromatogr. 1991,550, 51. (14) Schure, M. R.; Barman, B. N.; Giddings, J. C. Anal. Chem. 1989,61,2735. (1 5 ) Golshan-Shirazi, S.;Guiochon, G. In Theoretical Aduancement in Chromatography and Related Separation Techniques; Dondi, F.,Guiochon, G., Eds.; Kluwer Academic Publishers: Dordrecht, the Netherlands, 1992; pp 61-92. (16) Dondi, F.; Blo, G.; Remelli, M.;Rcschiglian, P. In Dondi, F., Guiochon, G., Eds.; Theoretical Aduancement in Chromatography and Relared Separarion Techniques; Kluwer Academic Publishers: Dordrecht, the Netherlands, 1992; pp 173-210.
0003-2700/94/0366-3066$04.50/0
0 1994 American Chemical Society
found to fit very well LC peaks except for some early eluted, highly asymmetrical peaks.” This model is derived when convolving the density function of the normal distribution and an exponential decay function: y ( t ) = h(t)*f(t) = c h ( t ’ ) f ( t - t’) dt’
&a
(2)
and f ( t ) is the exponential decay function,
f ( t ) =;e1
-t/r
, t10
(3)
The convolution of eqs 2 and 3 results in
Since y ( t ) contains the error function (erf), different approximations have been proposed. l 8 9 l 9 Mathematical subroutines are available for its computation,20 and many calculations with EMG peaks can easily be accomplished in frequency domains2’ Deconvolution of EMG Peaks. In the following, we will derive the function describing sharpened EMG peaks. This includes the deconvolution of symmetrical peaks, too, when substituting zero for all the 7 parameters. This derivation can easily be done in frequency domain utilizing some Fourier theorems. First, we calculate the Fourier transform of an EMG peak. According to the convolution theorem of Fourier transform,22 the convolution performed in time domain is equivalent to a multiplication in frequency domain. The Fourier transforms of h ( t ) andf(t) are
H ( ~=) ~ ~ - o ~ s ~ / 2 - j w m
(5)
and 1 1 +jar respectively. Applying the theorem above recalled, the Fourier transform of the EMG function is F(w) =
and the exponentially decaying part,
(1)
where the asterisk operator denotes the convolution process, h ( t ) is the Gaussian function, A e-(t-m)2/2s2 h(t) = -
both the Gaussian component,
(9)
are of unit area, the area of the original peak to be deconvoluted is conserved during the procedure. Moreover, as hd(f) is centered around the origin, the deconvolution performed by yd(f) = hd(f)*fd(f) does not affect the location of the original peak. Using yd(f) as a sharpening function, a wide variety of peak broadening effects can be deconvoluted. When Ud = 0, the Gaussian part of the peak is untouched; only the asymmetry is decreased, deconvoluting the total asymmetry if Td = 7 . Conversely, when 7d = 0, only the peak width is modified; the time constant is not affected. This latter case is, however, an unfavorable solution, since when the peak width is decreased, the asymmetry increases when the time constant is unchanged. The Fourier transform of yd(f) is e-U2Ud2/2
Yd(w) = 1 + jw7,
(10)
The deconvolution by yd(f) is identical to dividing Y(w) by yd(w):
X ( w ) can be divided into two parts, and a simple relationship exists between the two, recalling the derivative theorem of Fourier transform,22 which says that if f ( t ) has the Fourier transform F(w), thenf’(t) has the Fourier transform jwF(w):
-e*2d2/2-jwm
+
X ( w ) = A 1 jwr
e-du’2/2-jwm
+ AJU7d 1 + j w r
(12)
X , (a) X2(w) where d = (a2 - 6d2)1/2. Because of the theorem above recalled, the connection between the inverse Fourier transforms of the two parts is
~
Accordingly, the deconvoluted EMG peak can be written as
e-u2s2/2-jwm
Y(w) = H ( w ) F ( w ) = A
(7) 1+jar The deconvolution is the inverse procedure of convolution, i.e., we peel off the band broadening effects that lead to wide or asymmetrical band shape. For these reasons, it seems obvious to choose another EMG peak as a sharpening function. The sharpening function is the convolution of a unit area Gaussian peak centered around the origin (m= 0) with a unit area exponential decay function. In this way, we are able to decrease both the asymmetry and the width of the peak. When
The derivative of the above function is
(17) Naish, P. J.; Hartwell, S.Chromutographia 1988, 26, 285. (18) Foley, J. P.; Dorsey, J. G. J . Chromatogr. Sci. 1984, 22, 40. (19) Berthod, A. Anal. Chem. 1991, 63, 1879. (20) Press, W. H.; Teukolsky, S.A.; Vcttcrling, W. T.; Flanncry, B. P. Numerical Recipesin FORTRAN. The ArtofScientiJicComputing,2ndd;Cambridge University Press: Cambridge, U.K., 1992. (21) Felinger, A,; Pap, T. L.;InczCdy, J. Talanta, in press. (22) Bracewell, R. N . The Fourier Transform and Its Applications, 2nd 4.; McGraw-Hill Book Company: New York, 1986.
From the above equations we can get the shape of the sharpened
Comparing Xl(w) and eq 7, we can see that xl(t) is also an EMG peak, whose Gaussian component is of a’width, and the exponential decay part is unchanged:
AnalyticalChemistry, Vol. 66,No. 19, October 1, 1994
3067
peak:
T-Td[l -erf(L-’-m)] v5-7
e x p ( 5 - e )
v5-uf
(17)
This peak shape is a linear combination of a Gaussian peak of u’ width and an EMG function of parameters d and T :
RESULTS AND DISCUSSION The peak model given in eq 18 is unique for several reasons. Looking at x ( t ) , we can see that its EMG part is being suppressed as Td increases, but the tailing is always present. This is because the tailing is caused by 7 , and this parameter in the EMG part of x ( t ) is unchanged. It is rather surprising that the EMG part of the above function depends on the original T parameter. Since this parameter does not decrease during the deconvolution, we always have to pay attention to the tailing of the deconvoluted peaks. Of course, as Td approaches 7 , the tailing is less enhanced, but always present, except for the case when r = Td. When Td is increased from 0 up to 7 , the tailing part is linearly suppressed. The rear part of the peak does not get narrower as we might expect, but the tailing is gradually suppressed into the base line. The moments of x ( t ) can be calculated by using the wellknown relationship between the nth moment and the nth differential coefficient of F(o)at the origin:22
The skew and the excess are the measures of the peak shape. They can be calculated from the higher order central statistical moments:
The moments of x ( t ) show some logical consequences of the deconvolution procedure. As the deconvolution of the tailing reduces the asymmetry, the first moment decreases; the peak maximum approaches the value m, where the maximum of the Gaussian peak is located: 1.1, =
m+T
- T ~
(21)
The second central moment is a linear function of the squares of the peak shape parameters: jL,’
+
= d2
T2
- 72
(22)
The third, jL3’
= 273 - 2Td3
(23)
and the fourth central moments,
j~l= 3 d 4 + 6d2(r2- 72)+ 9r4 - 37: 3068
- 6T2Td2 (24)
AnalyticalChemistfy, Vol. 66,No. 19, October 1, 1994
ud
/.
Figure 1. Relative values of the peak shape parameters of the sharpening function neededto conserve the original skew (dashed line) and excess (solid line).
are needed to calculate the skew, 2T3 - 2Td3
S= ((7’2
+ r 2- Td*)3/2
and the excess,
We can calculate the shape parameters of the sharpening function (ad, Td) needed in order to conserve the original shape of the overlapping peaks. This is interesting when the peak shape holds some important information, but the individual peak shapes cannot be determined because of overlap. We can calculate the skew and the excess for both the original and the deconvoluted peaks (see eqs 25 and 26) and calculate the Ud and Td needed to keep the original S and E values. In Figure 1, we show the desired values of parameters to reach this situation. The dashed line shows the combination of the Ud and Td parameters at which the skew is not changed during the deconvolution. The solid line presents the values needed to conserve the excess. It can beseen that, besides the singular point where the peak is narrowed into a spike at Ud = u and Td = 7 , which makes no sense at all, it is impossible to conserve the peak shape of asymmetrical profiles during the deconvolution. In Figure 2, a series ofdeconvolutions with different transfer functions is shown. The original skewed peak was generated with u = 50 and T = 50 parameters. This is shown with solid lines in Figures 2a-f. When a transfer function of Ud = 25, Td = 25 is used, both the Gaussian part (dashed line) and the EMG part (dotted line) of the sharpened peak are comparable. The sharpened peak-the chain line, which is the sum of the dashed and dotted lines-does not really seem to be narrower, although the location of its maximum is lower and the peak is higher. When we choose a transfer function with Ud = 0 and Td = 50, the asymmetry is eliminated, but the Gaussian part is untouched, see Figure 2b. When the transfer function is close to the original peak, Ud = 40, Td = 40, only the Gaussian part of the deconvoluted peak is dominant, the EMG part is almost completely suppressed, see Figure 2c.
60
a
-
-
50
.t" 40 30
4 20 -
5
10
0
90 0
Flgure3. Absolute value of the relativeerror of peak areadetermlnation. = 50, T~ = T~ = 50; Rs = 0.6. First peak, surface of solid lines, dashed contour ilnes. Second peak, surface of dotted lines, chain contour Ilnes.
AI = 10, u1 = u2
times the peak standard deviation+
R,= .:.::.;. Flgure 2. Deconvolutionperformedwith different Sharpeningfunctions. The origlnai peak parameters were u = 50, 7 = 50 (solid lines). Chain line, sharpened peak; dashed line, its Gaussian part: dotted line, its EMG part. The peak shape parameters of the sharpening function Were (a) Ud = 25, 7 d = 25; (b) u d = 0,7 d = 50; (c) u d = 40, 7 d = 40; (d) ud = 0,Td = 100;(e) u d = 40, 7 4 = 100; and (f) = 50, 7 d
= 0.
Something very strange happens when the peak asymmetry is overcompensated. The parameters of the transfer function are Ud = 0 and 74 = 100 in Figure 2d, whereas Ud = 40 and 7 d = 100 in Figure 2e. To understand this effect, we have to recall eq 18. The EMG part of this equation is weighted by parameter 1 - 7d/7,which is negative if 7 d > 7. This fact will cause a negative side lobe on the right side of the peak. Usually it has been thought that during Fourier domain deconvolution, the cutoff of the domain of the noise causes spurious oscillations observed in the deconvoluted signal. That effect is, of course, present. But from the above result it is obvious that when we use a skewed transfer function to sharpen overlapping peaks, there is a theoretical explanation of why side lobes emerge independently of whether Fourier deconvolution or another technique was applied. Finally, in Figure 2f, all the Gaussian contribution is deconvoluted (Ud = 50), while the exponential decay part is not affected (74 = 0). This combination of the parameters is rather unfortunate, as the tailing of the peak is almost unaffected during the deconvolution procedure. A systematic grid search has been carried out in order to map the error of peak integration based on deconvoluted peaks. First, a symmetrical case, in which two identical peaks overlap, was looked at. The peaks were generated as 01 = u2 = 50, 71 = 72 = 50; the ratio of the two areas was Al/A2 = 1; the peak resolution was R,= 0.6. The two peak shape parameters of the transfer function were changed in the intervals: 0 5 Ud I 50, 0 I 74 I 100. The peak resolution was calculated as the difference of the first moments over the average base line width, where the base line peak width is estimated as 4
m2 + r2- m, - 71 2[(u,2
+ T12)1'2 + (a? + T22)1'2]
(27)
Figure 3 shows the absolute value of the relative error of peak area determination. In the front part of the figure, the two surfaces are identical. In this case, the absolute and the relative errors are equal as the peak areas are the same. Because of the tailing, the area of the first peak is underestimated when the skewed peaks still overlap (solid line). This causes the overestimation of the second peak area (dotted line). The contour lines show where the peak area is determined with minimum error for the first peak (dashed line) and for the second peak (chain line). In the area between the contour lines, the error is less than 2%. We can see that this happens exactly when the total asymmetry is deconvoluted, resulting in symmetrical peaks, whatever the Gaussian contribution is. This result is not surprising. On the other hand, it is astonishing to observe that no other combination of Ud and 7 d parameters gives error-free estimation of the peak area. When 7 d is increased over 71 = 72 = 50, the negative side lobes appear, which increases the error of the area calculation, as the peaks are not integrated for the domain of the negative side lobes. During the deconvolution, the total area of the peak is conserved. However, the emergence of the negative side lobes increases the area of the positive parts, which leads to the overestimation of the peak area. The contour line shows that extremely high 7 d also can give correct integration for the second peak. This contour line, along the configuration for base line separation, is caused by the deletion of the individual negative and positive errors. In Figure 4, the ratio of the peak areas is changed to AlIA2 = 2. The outline of the error surfaces is very similar to the previous situation, except that the values of the optimum Ud and 7 d parameters depend on how much peak broadening effect is deconvoluted. When the peak width is deconvoluted only partially, the optimum configuration depends on the ratio of the peak areas and the extent of skew. For this reason, it cannot be used as a rule; only the exact base line separation gives satisfying results. Deconvoluting a general case, when all peak parameters differ (A1 = 10, A2 = 5, u1 = 25, u2 = 75, 71 = 75, 72 = 25), although the variance of the peaks is identical, even at the AnaIyticalChemistry, Vol. 66,No. 19, October 1, 1994
3069
60
-
2c0
-
50 150
- 40
-
9 30
4
0
2130
\
7
20
a 50
0 10
-
0
2 0
d
90
00
0
Figure 4. Same as Figure 3, except 25
-
= IO,
Figure 7. Same as Figure 3, except Al = 10, A2 = 25, 7 1 = 100, 7 2 = 25.
= 5.
=
1, u1 = 25, u2
r------
-
20
15
1
40 lo 5 2 90
0
Flgure 5. Same as Flgure 3, except A, = 10, A2 = 5, u1 = 25. u2 = 75, 71 = 75, 7 2 = 25. 200
-
-
150
100 200 300 400
9 100 T
a 50 5
500 600 7Ci
Flgure 8. Deconvolutlon of the case shown in Flgure 7, wlth dlfferent methods. The parameters of the sharpening functlon are = 10, r d = 45.
1
0 90
0
Figure 8. Same as Figure 3, except AI = 1, A2 = 10, u1 = 25, o2 = 25, 71 = 25, 72 = 100.
maximum possible deconvolution of the Gaussian contribution and at total deconvolution of the asymmetry, the integration for the two peaks requires different sharpening functions; see the contour lines in Figure 5. Strangely, this means that whatever sharpening function we use, the integration will not be faultless. The further results are still more disappointing when the peak shape parameters are very different. Fortunately, this seldom happens in practice. In Figure 6, we present the errors for the two overlapping peaks: A1 = 1, A2 = 10, u1 = 25, 6 2 = 25, 71 = 25, 72 = 100. In this case, the integration of the first peak is error-free when exactly the asymmetry of the first peak is deconvoluted. When the time constant of the transfer function is larger than this threshold, the error increasesvery fast, because of the small relative concentration ofthis component. Becauseof the big concentrationdifference, the integration of the second peak gives error smaller than 2% in a wide range, but the integration for this component is never error-free. 3070
Analytical Chemistty, Vol. 66, No. 19, October 1, 1994
When we simply swap the two peaks, as results show in Figure 7, no error-free peak area can be determined for either peak. Moreover, the second peak is minor, slightly skewed, and hidden in the tail of the main peak. Because of the enhanced tailing of the main peak, both peaks are seemingly very asymmetric, see Figure 8. But when the asymmetry is deconvoluted, the negative side lobe of the second component and the still present tailing of the first peak interfere. Consequently, a spurious peak appears, mistakenly indicating the presence of another component. It is very important to emphasize that the model derived in this study explains why this happens. These kinds of effects-the emergence of side lobes or spurious peaks-have usually been attributed to the numerical method used for deconvolution itself. Figure 8 demonstrates that independent of the method used, the same effect is forecast by the theoretical model, although this does not mean that the numerical methods cannot contribute to this effect. In the numerical calculations, the direct Fourier deconvolution was carried out on the basis of eq 11 by an FFT subroutine. The iterative relaxation method was used with the following algorithm:11J2
a
b
c
d
numerically calculated signal oscillates around the one calculated by our model (see eq 18).
e
Figure B. Overlapplng peaks analyzed during the grid search. For the parameters of peaks a-e, see Figures 3-7, respectively.
Figure 10. Deconvolutionof noisy chromatograms by Fourler method. The original peaks (solid Ilne) and the sharpened peaks calculated by Fourier method (dotted llne) and by eq 18 (dashed line) are shown. A1/A2= 2.5, ml = 110, m2 = 150, u1 = a2 = 10, r1 = r2 = 15, = 3.5,7 4 = 20.
where the r [ ~ ( relaxation ~)] function, r[x(Q]= r,(l - 2/cIdk)- c/21)
(29)
constrains the estimated between values 0 and c. The different overlapping peaks analyzed during the grid search are plotted in Figure 9. Finally, we have toemphasize again that during the analysis of noisy chromatograms, the presence of noise (via the apodization function), of course, causes oscillations in the deconvoluted signal, because in this case the cutoff is at low frequencies. This is demonstrated in Figure 10: two noisy chromatograms ( S I N = 50 and 20) are deconvoluted by the Fourier method. It can be seen that the amplitude of this oscillation depends on the signal-to-noise ratio; moreover, the
CONCLUSIONS The analytical expression derived for sharpened, overlapping peaks explores many aspects of peak deconvolution. Although the EMG peak model has been studied, the qualitative results can be extended to other skewed peaks. The most important factor is that we have to know exactly the amount of skew we want to get rid of. For this reason, great care has to be taken when selecting the sharpening function. When the peak asymmetry is overestimated or overcompensated, negative side lobes will appear. Depending on the configuration, i.e., the asymmetry, or the area of the individual peaks, this can lead to the emergence of false peaks. False peaks can appear particularly when, for instance, a constrained deconvolution technique is applied that does not allow negative values. A rather surprising result is that the sharpened peak is a linear combination of a Gaussian and an EMG function. Moreover, the shape of the tailing part of the sharpened peak does not depend on 7d; for this reason, the tailing does not narrow during the deconvolution, but it is gradually suppressed into the base line. Since we have to be very careful about the value of Td, precise deconvolution can be expected only when the total skew is deconvoluted. This makes the deconvolution rather problematic. Consequently, when we suspect that the shapes of the individual peaks differ, deconvolution should not be used. In this case, only base line separation can give reliable results. This can often happen, since the dispersion increases linearly with retention time, and it has been shown that the time constant 7 also depends on retention." GLOSSARY peak area maximum peak height allowed by the iterative deconvolution excess exponential decay function exponential decay part of the sharpening function Fourier transform off Gaussian function Gaussian part of the sharpening function Fourier transform of h imaginary unit index for iteration location of the maximum of the Gaussian peak relaxation function maximum value of the relaxation function peak resolution skew sharpened peak parts of the sharpened peak Fourier transform of x Fourier transforms of X I and x2, respectively original EMG peak sharpening function Fourier transform of y Fourier transform of yd time Analytical Chemisw, Vol. 66, No. 19, October 1, 1994
3071
t' P P' a U'
ad 7
8072
dummy variable statistical moment statistical moment centered around the mean Gaussian contribution to the peak width Gaussiancontribution to thewidth of the sharpened peak Gaussian contribution to the width of the sharpening function time constant off
Analytical Chemistry, Vol. 66, No. 19, October 1, 1994
time constant of fd frequency
Td
w
~ for ~review~February b 10,d 1994. Accepted June 17,
~ 1994,"
*Abstract published in Aduunce ACS Absrrucrs. August 1, 1994.
