Anal. Chem. 2001, 73, 4862-4872
Predicting Peak Shape in Capillary Zone Electrophoresis: a Generic Approach to Parametrizing Peaks Using the Haarhoff-Van der Linde (HVL) Function Guillaume L. Erny, Edmund T. Bergstro 1 m, and David M. Goodall*
Department of Chemistry, University of York, York YO10 5DD, U.K. Sally Grieb
Pfizer Global Research and Development, Pfizer Limited, Kent CT13 9NJ, U.K.
We have found that the Haarhoff-Van der Linde (HVL) peak function provides excellent fitting to the shapes of CZE peaks. Initially designed for overloaded peaks in gas chromatography, this function describes a Gaussian peak when there is no peak distortion, and a triangular peak when there is no diffusional peak broadening. As such, it is ideal for CZE peaks distorted by electromigration dispersion (EMD). Fitting peaks with this function gives four parameters: three of them can be related to the Gaussian peak that would have been obtained in case of no EMD; the last one is a measure of the peak distortion. Using moving boundary theory, this peak distortion parameter may readily be expressed in terms of analyte and background electrolyte mobilities and concentrations, electric field, and sample injection length. The variance of an HVL peak is shown to be described by a universal function, and a master equation is presented. The region where EMD adds less than 10% to the Gaussian variance is shown to be very narrowly spread around the mobility matching condition. Under typical CZE operating conditions with an analyte at 1% of the BGE concentration, significant peak distortion is always present. Because the total peak variance is not an addition of the Gaussian and triangular contributions, the HVL model and the methodology introduced here should always be used to correctly combine variances. Electromigration dispersion (EMD) is a problem in efforts to extend the linear dynamic range of capillary electrophoresis (CE). It can drastically spoil a separation by increasing peak variance. Because the EMD is due to a difference of electric field strength between the analyte zone and the background electrolyte (BGE),1-4 it can be reduced to an acceptable level only by matching the * To whom correspondence should be addressed. (1) Mikkers, F. E. P.; Everarts, F. M.; Verheggen, Th. P. E. M. J. Chromatogr. 1979, 169, 1-10. (2) Sustacek, V.; Foret, F.; Bocek, P. J. Chromatogr. 1991, 545, 239-248. (3) Gas, B.; Stedry, M.; Kenndler, E. Electrophoresis 1997, 18, 2123-2133. (4) Poppe, H. Anal. Chem. 1992, 64, 1908-1919.
4862 Analytical Chemistry, Vol. 73, No. 20, October 15, 2001
co-ion mobility to the analyte mobility.1,3,5-8 Some work has been done to try to match as precisely as possible the two ions.6,8-10 A number of different approaches have been used to study EMD. In computer simulation, the capillary is divided into a series of cells, each of which is considered to be homogeneous. Concentrations in each cell are calculated at small intervals of time using a numerical algorithm to solve the one-dimensional migrationdiffusion differential equation.11-13 This approach gives very good results, but it is time-consuming, and problems of numerical diffusion tend to degrade the precision as the overall run time increases.6 In another approach developed by Mikkers, EMD is studied assuming that diffusional processes can be neglected. In this approach, the concentration distribution of the analyte due to EMD can be predicted using the initial condition.14 However, no link has been made between the triangular distribution that is predicted and the peak shape where EMD and diffusion are both taken into account. In this paper we will show how the Haarhoff-Van der Linde (HVL) function15,16 provides a perfect fit to a CE peak distorted only by EMD. More importantly, the HVL peak model can be viewed in a first approximation as convoluting a Gaussian peak due to diffusion inside the column and a triangular peak due to EMD. Using the theory developed by Mikkers,1,14 we have been able to parametrize the triangular peak in terms of concentrations (5) Williams, R. L.; Childs, B.; Dose, E. V.; Guiochon, G.; Vigh, G. Anal. Chem. 1997, 69, 1347-1354. (6) Dose, E.; Guiochon, G. Anal. Chem. 1991, 63, 1063-1072. (7) Rawjee, Y. Y.; Williams, R. L.; Vigh, Gy. Anal. Chem. 1994, 66, 37773781. (8) Stalberg, O.; Hedeland, M.; Petterson, C.; Westerlund, D. Chromatographia 1998, 48, 415-421. (9) Williams, R. L.; Vigh, Gy.; J. Chromatogr. A 1996, 730, 273-278. (10) Williams, R. L.; Vigh, Gy. J. Liq. Chromatogr. 1995, 18, 18-19. (11) Martens, J. H. P. A.; Reijenga, J. C.; Boonkkamp, J. H. M. T.; Mattheij, R. M. M.; Everaerts, F. M. J. Chromatogr. A 1997, 772, 49-62. (12) Reijenga, J. C.; Kenndler, E. J. Chromatogr. A 1994, 659, 403-415; 417426. (13) Mosher, R. A.; Saville, D. A.; Thormann, W. The Dynamics of Electrophoresis; VCH: Weinheim, 1992. (14) Mikkers, F. E. P. Anal. Chem. 1999, 71, 522-523. (15) Haarhoff, P. H.; Van der Linde, H. J. Anal. Chem. 1966, 38, 573. (16) PeakFit 4.0 for Windows, User’s Manual; SPSS Science: Chicago, 1999. 10.1021/ac010758g CCC: $20.00
© 2001 American Chemical Society Published on Web 09/14/2001
and mobilities of the different ions present in the system. Master equations will be presented that allow calculation of the extra variance added to the CE peak due to EMD. METHOD A dynamic simulation program written by Reijenga was used to simulate CE peaks.11,17 The column parameters were column i.d., 75 µm; length to the detector, 10 mm. The pKa values of all cations used as co-ion or analyte were set at 14; in this software, this value is the default value for fully charged cations. The buffer was an acid HA with pKa 4.0; its conjugate base, A-; and its pH, fixed at 4.7. The mobility of A-, the counterion, was selected to be -5.0 × 10-8 m2 V-1 s-1, giving an effective mobility at pH 4.7 of -4.17 × 10-8 m2 V-1 s-1. The current was fixed at 60 µA. The program offers a choice of algorithms; for our simulations, we used the default algorithm, dime. The electric field in the BGE was calculated by the software for the different buffers simulated. The HVL function in the program PeakFit v4 (SPSS Inc., Chicago)16 was used to fit simulated peaks. RESULTS AND DISCUSSION HVL Function. The HVL function was developed for nonideal gas chromatography and can describe both tailing and fronting peaks.16 The HVL function is defined by four parameters as shown in eq 1
[ ( )]
1 t - a1 2 exp 2 a2 a1a3x2π f(t) ) t - a1 1 1 + 1 + erf a1a3 2 x2a2 exp -1 a22 a0a2
( )
[ ( )]
Figure 1. Fitting of simulated data with the HVL function.
Case 1: a22 . a1|a3|
f(t) ≈
2
(2)
When the distortion parameter a3 is very small in comparison to a22/a1, the HVL function is equal to a Gaussian function with area a0, peak center a1, and standard deviation a2. Case 2: a22 , a1|a3|. If a3 is positive,
f(t) ) (1)
[ ( )]
1 t - a1 exp 2 a2 a2x2π a0
a0 a0t for a1 - x2a1a3 < t < a1 a3 a1a3
f(t) ≈ 0 for t < a1 - x2a1a3 and t > a1
(3)
If a3 is negative,
f(t) )
a0 a0t for a1 < t < a1 + x2a1|a3| a3 a1a3
where a0 is the peak area; a1, the peak center; a2, the standard deviation for a nondistorted peak; and a3, a measure of the peak distortion. For a Gaussian peak, a1 corresponds to the retention time in chromatography or the migration time in CE. Eq 1 is given with time, t, as the variable. If it is required to give a spatial rather than a temporal description of the HVL function, an equation analogous to eq 1 may be used with t replaced by x, where x is distance. As shown in Figure 1, the quality of the fit with the HVL function in the case of a CE peak distorted only by EMD is close to perfect. The analyte and BGE conditions were 3, 5, and -4.17 × 10-8 m2 V-1 s-1 for analyte, co-ion mobilities, and counterion effective mobility, respectively; the analyte concentration was 5 mM; the BGE concentration 50 mM; and the electric field was 28.11 kV m-1. Limit Functions. Two limit functions can be derived from the HVL function: the first when a22 . a1|a3|, the second when a22 , a1|a3|. Derivation of these equations is given in the appendix.
When the diffusional process can be neglected, that is, a22/a1 is small in comparison with a3, the HVL function is seen to be a triangular function, with origin at t ) a1, slope equal to - a0/ a1a3, and width at base equal to x2a1|a3|. A summary of these results is given in Figure 2. The triangular function is independent of a2 and describes EMD when diffusion is negligible. The Gaussian function is independent of a3, and describes diffusional peak broadening when EMD is negligible. Statistical Moments. When the peaks are non-Gaussian, the different parameters characterizing peaks in chromatography and electrophoresis can be calculated using statistical moments.18-20 Definitions and significance of the statistical moments are given in Table 1. Values of the moments for the Gaussian and triangular function are given in Table 2. An analytical solution for the integral
(17) Dynamic simulation program available from Reijenga, J. C., Dept. of Chemical Engineering and Chemistry, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands; E-mail:
[email protected].
(18) Foley, J. P.; Dorsey, J. G. Anal. Chem. 1983, 55, 730-737. (19) Dyson, N. Chromatographic Integration Methods; Royal Society of Chemistry: Cambridge, U.K. 1990. (20) Morton, D. W.; Young, C. L. J. Chromatogr. Sci. 1995, 33, 514-524.
f(t) ≈ 0 for t < a1 and t > a1 + x2a1|a3|
Analytical Chemistry, Vol. 73, No. 20, October 15, 2001
(4)
4863
As can be seen in Table 3, the true second moment of the HVL function, M2T, is not given by the additivity of variances, as in eq 7. The error in the second moment, ξM2, calculated using eqs 6 and 7, is shown in Figure 3 as a function of G. It is evident that the error is systematic and follows a universal function that applies for all combinations of a0, a1, a2, and a3. The basis for a universal function for evolution of an HVL peak in CE lies in the interplay between diffusion and electromigration dispersion. Viewed in a time frame in which the origin moves with the peak center, any given time slice of the triangular function subsequently broadens as a result of diffusion but also shifts because of EMD. The universal function for the error in M2 was fitted with an r2 of 0.9994 by the following polynomial:
ξM2(G) )( -6.44 × 10-3)G2 - (2.80 × 10-4)G3 + Figure 2. Triangular, Gaussian, and HVL peak shapes for a0 ) 1, a1 ) 5, a2 ) 0.5, a3 ) 0.3, a′3 ) -0.3: (a) tri (a0, a1, a3); (b) tri (a0, a1, a′3); (c) Gau (a0, a1, a2); (d) HVL (a0, a1, a2, a3); (e) HVL (a0, a1, a2, a′3). Table 1. Definition of Statistical Moments and of Some Chromatographic Parameters moment
symbol
general expression
zeroth first second central nth central
M0 M1 M2 Mn γS γE N
+∞ ∫-∞ f(t) dt +∞ ∫-∞ t × f(t) dt/M0 +∞ ∫-∞ (t - M1)2 f(t) dt/M0 +∞ ∫-∞ (t - M1)n f(t) dt/M0 M3/M23/2 (M4/M22) - 3 M12/M2
chromatographic parameter peak area peak center peak variance peak skew peak excess efficiency
of an HVL peak could not be found. Instead, we postulate that the value of the different moments will be a combination of the moments of both limit functions. HVL peaks were calculated using the PeakFit program. Parameters for a0, a1, a2, and a3 were chosen randomly, and true values of M1 and M2, M1T and M2T, were calculated by the software. Some values and calculation for the different parameters and moments are shown in Table 3. We define the percentage of Gaussian variance, G, as
(6.56 × 10-6)G4 - (3.11 × 10-8)G5 (8) Because the variance due to diffusion is normally easy to calculate from parameters of the analyte in CE, it is helpful to consider how the total variance of an HVL peak relates to M2(Gau). Combination of eqs 5, 6, and 8 gives
M2T 100 100 ) M2(Gau) 100 - ξM (G) G
( )
(9)
2
Figure 4 shows plots of M2/M2(Gau) over the range G ) 50100%. Curve a is calculated using eq 7 to determine M2; curve b uses eq 9 to determine M2, which in this case is the true variance, M2T. It is evident that the assumption of no correlation of variances (curve a) will always provide an overestimate of the true variance (curve b). For example, at G ) 50%, M2/M2(Gau) ) 2 from summation of second moments, whereas the true value is M2T/ M2(Gau) ) 1.67. The plot of eq 9 in Figure 4 was successfully fitted (r2 ) 0.99994), for G ) 50%, with a polynomial of degree 3,
M2T ) 6.115 - (1.523 × 10-1)G + (1.524 × 10-3)G2 M2(Gau) (5.12 × 10-6)G3 (10)
If the HVL were simply a convolution between a Gaussian and a triangular peak, and the variances were uncorrelated, then the second moment of the HVL function would be the summation of the second moments of the Gaussian and triangular peaks.
When 90% e G e 100%, M2T/M2(Gau) is close to unity; the peak starts to be distorted, but the peak width and, thus, the variance stay almost constant. At G ) 90% (M2(tri) ) 0.1 M2(Gau)), M2T/M2(Gau) ≈ 1.015, so the peak variance is only 1.5% more than it would be if the peak were Gaussian. Using the same approach as that used for setting the limit for extracolumn dispersion in LC, which is that in order not to spoil resolution, the extracolumn variance should be no more than 10% of the column variance,21 we evaluate the point in Figure 4 where EMD adds 10% to the Gaussian variance. At G ) 75% (3M2(tri) ) M2(Gau)), M2T/M2(Gau)) 1.1. Thus EMD would not be expected to significantly degrade resolution until its variance is more than one-third of the Gaussian variance. This finding will be used in a subsequent section when modeling CE peaks in more detail. This section has shown that the HVL function is not a simple convolution of a Gaussian function and a triangular function, but,
M2 ) M2(tri) + M2(Gau)
(21) Klinkenberg, A. Gas Chromatography 1960; Scott, R. P. W., Ed.; Butterworths: London, 1960, p 194.
G)
M2(Gau) M2(Gau) + M2(tri)
× 100
(5)
and the error in the ith moment, ξMi, as
ξMi )
4864
MiT - Mi MiT
× 100
Analytical Chemistry, Vol. 73, No. 20, October 15, 2001
(6)
(7)
Table 2. Statistical Moments of the HVL, Triangular and Gaussian Functions
Table 3. Calculations of Moments by PeakFit Program, and Using Eqs 5, 6, 7, and Equations Given in Table 2a a0 au
a1 au
a2 au
a3 au
M1T au
M2T au2
M1(tri) au
M2(Gau) au2
M2(tri) au2
G%
M2(Gau) + M2(tri) au2
ξM2 %
1 1 1 1 1 1 1 1 1 1 1 1 0.1 0.1 6 5 0.1 0.1 1 1 2 5 1
6 6 6 6 6 6 6 6 7 8 9 10 3 2 3 10 3 2 6 6 6 6 6
0 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.2 0.05 0.2 2 0.2 0.05 0.3 0.3 0.3 0.3 0.3
0.1 0.1 0.08 0.06 0.04 0 0.2 0.5 0.2 0.2 0.2 0.2 0.3 0.01 0.3 0.3 0.03 0.001 0.09 0.08 0.08 0.09 0.07
5.27 5.60 5.67 5.75 5.83 6 5.28 4.62 6.19 7.10 8.02 8.95 2.23 1.91 2.23 9.56 2.87 1.99 5.55 5.59 5.59 5.55 5.63
0.067 0.18 0.18 0.17 0.16 0.16 0.24 0.44 0.26 0.28 0.31 0.33 0.13 0.0039 0.13 4.03 0.042 0.0025 0.12 0.12 0.12 0.12 0.11
5.27 5.27 5.35 5.43 5.54 6 4.97 4.37 5.88 6.81 7.74 8.67 2.11 1.87 2.11 8.37 2.72 1.96 5.31 5.35 5.35 5.31 5.39
0 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.16 0.04 0.0025 0.04 4 0.04 0.0025 0.09 0.09 0.09 0.09 0.09
0.067 0.067 0.053 0.04 0.027 0 0.13 0.33 0.16 0.18 0.2 0.22 0.1 0.0022 0.1 0.33 0.010 0.0002 0.06 0.053 0.053 0.06 0.047
100 29.41 25 20 14.29 0 45.46 67.58 49.30 52.63 55.56 58.14 71.43 47.06 71.43 7.69 20 8.16 40 37.21 37.21 40 34.15
0.067 0.23 0.21 0.2 0.19 0.16 0.29 0.49 0.32 0.34 0.36 0.38 0.14 0.0047 0.14 4.33 0.05 0.0027 0.15 0.14 0.14 0.15 0.14
0 -22.49 -20.82 -17.93 -13.58 0 -21.49 -10.38 -19.97 -18.41 -16.92 -15.54 -8.36 -20.90 -8.36 -7.59 -17.92 -8.02 -23.15 -23.46 -23.46 -23.15 -23.46
a
au ) arbitrary units.
because it is a universal function, the true variance can be calculated for any value of G using eq 8 or eq 10. Simulations and HVL Fitting over a Range of Analyte Concentrations. Equations 7 and 9 were applied to CE peaks distorted by EMD for the case of a buffer with a single co-ion, and a fully charged cationic analyte. For this set of experiments, the BGE
conditions were held constant, and the peak shape was studied as a function of analyte concentration. The peak shapes were simulated using a dynamic CE simulation, and values of the HVL parameters were found by fitting the CE peaks with PeakFit. BGE conditions and results are shown in Table 4, and peak shapes as a function of the analyte concentration are presented in Figure 5. Analytical Chemistry, Vol. 73, No. 20, October 15, 2001
4865
Figure 3. Error in calculation of M2 as a function of percentage of Gaussian variance.
Figure 4. Ratio of total to Gaussian peak variance for an HVL peak as a function of G, with total variance calculated (a) assuming additivity of variances and (b) using the universal function for the error in M2 to give the true variance.
The value of the electric field in the BGE was calculated by the program to equal 25.56 kV m-1. The HVL model fits the different peak shapes perfectly, with r2 values always >0.998. The value of the peak center, a1, extracted from the fit should be independent of analyte concentration and equal to the theoretical value of 19.56 s. Extracted values lie in the range 19.54-19.87 s, and all differ by 15% for highly distorted peaks (G < 30%). In contrast, mobilities calculated using a1 give errors in µA of 0 (A4) 2
e-t dt ≈ 2
e-x xπx 2
1 + erf(x) ≈ 2 +
for x large, and x < 0 (A5)
1 ≈ exp(-x) exp(x) - 1
1 ≈ -[1 + exp(-x)] exp(x) - 1
[ ( )]
1 t - a1 2 exp 2 a2 a1a3x2π f(t) ) t - a1 1 1 + 1 + erf a1a3 2 x2a2 exp -1 a22 a0a2
[ ( )]
( )
( )
exp(x) ≈ 1 + x for small x, x positive or negative
a1a3
exp -
Case 1: a2 . a1|a3|. The HVL function can be simplified under this condition, using the following approximation:
[ ( )]
a1a3 2
>
a2
f(t) ≈
a22 a1a3
[ ( )]
2
a1a3x2π a2
f(t) )
)
a0t a0 a3 a1a3
x2π(a1 - t)
a1a3 2
a2
a1a3
0 and a3 < 0 shown in Figure 2.
)
[ ( )] 2
(A3)
When the distortion parameter, a3, is very much less than a1/ a22, the HVL function is equal to a Gaussian function with area, a0; peak center, a1; and standard deviation, a2. 4872
a0a2
( )
1 a1 - t 2 a2
a22
1 t - a1 exp 2 a2 a2x2π a0
(A8)
When
equation A2 can be approximated by
2
2
When
[ ( )]
1 t - a1 exp 2 a2
1 a1 - t exp 2 a2 x2π(a1 - t) a2
(A2)
a22 t - a1 1 + 1 + erf a1a3 2 x2a2
a22 x - a1 1 . 1 + erf a 1a 3 2 x2a2
a1a3x2π
[ ( )] [ ( )] 2
When
2
[ ( )]
+
a22
By definition, the values of the error function lie between 1 and -1. Because
a0a2
(A7)
(A1)
Using eq A1, the HVL function is equivalent to
f(t) )
for x large, and x < 0
1 t - a1 exp 2 a2 a1a3x2π a 0a 2
2
1 t - a1 exp 2 a2 a1a3x2π
(A6)
For example, when a3 > 0 and t < a1, using eqs A4 and A7, the HVL function becomes
f(t) ≈
a0a2
for x large, and x > 0
Analytical Chemistry, Vol. 73, No. 20, October 15, 2001
ACKNOWLEDGMENT Support from the University of York and Pfizer Global Research and Development for a graduate research studentship for G. L. Erny is gratefully acknowledged. Received for review July 9, 2001. Accepted July 9, 2001. AC010758G
