Anal. Chem. 2005, 77, 840-847
General Approach to High-Efficiency Simulation of Affinity Capillary Electrophoresis Ning Fang and David D. Y. Chen*
Department of Chemistry, University of British Columbia, Vancouver, British Columbia, Canada V6T 1Z1
The differential equation describing electrophoretic migration can be evaluated with various finite difference schemes (FDSs). However, the accuracy and efficiency can be dramatically different depending on the FDS chosen and the way the algorithm is implemented in a computer simulation program. The monotonic transport scheme is used as the algorithm for the hyperbolic part of the differential equation, and the first-order fully explicit scheme is used for the parabolic part of the equation. The combination of these algorithms minimizes the errors and maintains high efficiency. A circular arrangement of the cells in the computer’s memory is used in the implementation of the algorithms, and the use of concentration thresholds to enable and disable cells along the capillary makes the new algorithm highly efficient. Either thermodynamic or kinetic constants can be used in this program to simulate binding interactions between two species for equilibrium and nonequilibrium affinity CE. Simulation results with various parameters are presented. The simulated peak with proper parameters for an equilibrium affinity CE experiment has shape and position similar to that of the experimental peak. The simulated electropherograms for a nonequilibrium affinity CE experiment also show characteristics of the experimental electropherograms. As the knowledge of analyte behavior in capillary electrophoresis (CE) accumulates, it becomes obvious that the theory based on the observation of retention time and peak shape at the detector is inadequate for describing the analyte migration and bandbroadening processes. The behavior of analyte in a separation column can be described in detail by appropriate differential equations.1 However, it can be extremely challenging to expand the equations into solvable forms and design the experiment to meet the requirement of the simplified equations in order to solve for useful parameters such as the binding constant and complex mobilities. Computer simulation provides a unique opportunity for accurately describing the behavior of analytes in a column when the differential equations are translated into finite difference schemes. However, even with the same algorithm, different implementations could result in significantly different computing times. Because it usually requires hundreds to thousands of simulated runs to solve for thermodynamic or kinetic parameters for the interaction of a single pair of compounds, computation * To whom correspondence should be addressed. Tel. +01 (604) 822-0878. Fax +01 (604) 822-2847. E-mail:
[email protected]. (1) Galbusera, C.; Thachuk, M.; De Lorenzi, E.; Chen, D. D. Y. Anal. Chem. 2002, 74, 1903-1914.
840 Analytical Chemistry, Vol. 77, No. 3, February 1, 2005
time becomes crucial for the successful application of computer simulation in high-throughput applications where tens of thousands of compounds need to be characterized in a short time. With appropriate implementations of the algorithm used for describing the analyte migration behavior in CE, the computing time can be significantly reduced. Because CE has a well-controlled physical and chemical environment, it can be used to determine binding or dissociation constants. When CE is used for these studies, it is often referred to as affinity CE or ACE.2 There are two strategies used in ACE: the equilibrium method and the nonequilibrium method. In the equilibrium method, a capillary is filled with a background electrolyte (BGE) with an additive at various concentrations, and the analyte is injected into the capillary to form a narrow plug. When a high voltage is applied on the capillary, the analyte plug migrates toward the detector at a speed determined by the interaction between the analyte and the additive. This method is suitable for studying fast protein-ligand interactions. In the nonequilibrium method, the capillary is filled with a BGE that does not have mobility altering additives, and the sample injected is a mixture of analyte and additive. When a high voltage is applied, the complex is separated from the free analyte and additive as they migrate at different velocities. Eventually, several portions of the analyte band will pass through the detector window rather than only one peak as in equilibrium CE. For traditional nonequilibrium CE, the key to a successful binding study is that the dissociation rate constant (k-1) must be small enough for the complex to be held together before it is completely separated from the starting materials. Berezovski et al. introduced nonequilibrium capillary electrophoresis of equilibrium mixtures (NECEEM) as a method to determine kinetic and equilibrium parameters of protein-DNA interactions.3,4 In NECEEM, k-1 is determined from the exponential part of the electropherogram by fitting the experimental data with a single-exponential function. Since Bier et al. used partial differential equations to describe the electrophoretic migration process,5 several computer simulation models of CE have been developed.6-14 The finite difference (2) Chu, Y.-H.; Avila, L. Z.; Gao, J.; Whitesides, G. M. Acc. Chem. Res. 1995, 28, 461-468. (3) Berezovski, M.; Krylov, S. J. Am. Chem. Soc. 2002, 124, 13674-13675. (4) Berezovski, M.; Nutiu, R.; Li, Y.; Krylov, S. Anal. Chem. 2003, 75, 13821386. (5) Bier, M.; Palusinski, O. A.; Mosher, R. A.; Saville, D. A. Science 1983, 219, 1281-1287. (6) Saville, D. A.; Palusinski, O. A. AIChE J. 1986, 32, 207-214. (7) Palusinski, O. A.; Graham, A.; Mosher, R. A.; Bier, M.; Saville, D. A. AIChE J. 1986, 32, 215-223. (8) Dose, E. V.; Guiochon, G. A. Anal. Chem. 1991, 63, 1063-1072. (9) Ikuta, N.; Hirokawa, T. J. Chromatogr., A 1998, 802, 49-57. 10.1021/ac048938o CCC: $30.25
© 2005 American Chemical Society Published on Web 12/23/2004
techniques used in these models can be divided into two major categories: the Euler/Runge-Kutta integration methods6-10,13,14 and the finite difference scheme (FDS) with artificial dispersion.11,12 The Euler/Runge-Kutta methods are unstable due to spurious oscillations and numerical instability as demonstrated by Ermakov et al.11 The FDS with artificial dispersion is stable but could produce undesired cusps at sharp gradients. In this paper, various FDSs are compared, and a more robust and accurate one is implemented. Most simulation models of CE carefully account for electroneutrality and mass balance. All of the ions (H+, OH-, conjugate acid and base forms of buffer ions, and the ions of analytes and additives) in the system are considered in the calculation, and a large number of differential equations had to be used. As a consequence, the computing time is long, and the memory required to store simulation data is large. In practice, H+, OH-, and buffer ion concentrations are often kept near constant and are not the point of interest in a well-buffered system where the analyte and additive concentrations are much less than the BGE. If only the ions of analytes and additives are monitored, as in the model proposed in this paper, the local electric field strength can be assumed constant throughout the capillary. We will demonstrate that the simulation results agree well with the experimental data even with the assumption of constant local electric field strength. The computing time and the memory requirement can be significantly reduced. Dose and Guiochon implemented a simulation model, in which the calculations were carried out only for the active cells that contained ion concentrations above a set threshold.8 Theoretically, this method could reduce the computing time by up to 90% because only a small portion of the capillary was monitored. However, their implementation of activation/deactivation of cells relied on complicated column segmentation, which required frequent dynamic operations on the data, such as relocating, creating, merging, deleting, etc. These operations increased the complexity of the simulation program, as well as the computing time. Thus, the performance boost was only 3-30%.8 A computer simulation model of ACE was developed by Andreev et al.15 In their model, the rate constants, instead of binding constants, were included in the differential equations to describe the electrophoretic migration and the binding interaction at the same time. For many systems, the rate constants are often not available, but the binding constants are. In equilibrium CE, the binding constants are sufficient for the simulation of fast binding interactions. Recently, Okhonin et al. developed a mathematical model for NECEEM.16 Their key assumptions are that interacting species are separated immediately and that there is no new complex formed during the CE process. These assumptions made it (10) Ikuta, N.; Sakamoto, H.; Yamada, Y.; Hirokawa, T. J. Chromatogr., A 1999, 838, 19-29. (11) Ermakov, S. V.; Bello, M. S.; Righetti, P. G. J. Chromatogr., A 1994, 661, 265-278. (12) Ermakov, S.; Mazhorova, O.; Popov, Y. Informatica 1992, 3, 173-197. (13) Dubroa´kova´, E.; Gasˇ, B.; J., V.; Smolkova´-Keulemansova´, E. J. Chromatogr., A 1992, 623, 337-344. (14) Gasˇ, B.; Vacı´k, J.; Zelensk×c6, I. J. Chromatogr., A 1991, 545, 225-237. (15) Andreev, V. P.; Pliss, N. S.; Righetti, P. G. Electrophoresis 2002, 23, 889895. (16) Okhonin, V.; Krylova, S. M.; Krylov, S. N. Anal. Chem. 2004, 76, 15071512.
possible to find solutions for the differential equations describing the mass-transfer process with the Green’s function. However, the application of this model is limited to NECEEM. On the other hand, FDSs are still needed to obtain the precise solutions of the integrals describing the migration of the dissociated reactants. In this paper, a new algorithm is proposed to simulate both equilibrium and nonequilibrium CE. A number of finite difference schemes are considered, and the one with minimal errors was implemented to evaluate the differential equations describing the analyte behavior in ACE. The electrophoretic migration and the binding interaction are considered in two separate steps. The binding constants and rate constants are used for the simulation of equilibrium CE and nonequilibrium CE, respectively. To eliminate the laborious dynamic operations on the data, a circular arrangement of cells is used to implement the algorithm. As a result, more than 90% decrease in the computing time is achieved. EXPERIMENTAL SECTION A well-characterized 1:1 interaction between p-nitrophenol and β-cyclodextrin (β-CD) is chosen to compare the experimental and simulation results. The experiments were carried out on a Beckman Coulter ProteomeLab PA800 (Beckman Coulter, Inc., Fullerton, CA) with a built-in PDA detector. A 64.5-cm-long (54.3 cm to detector), 50µm-i.d. fused-silica capillary (Polymicro Technologies, Phoenix, AZ) was used. β-CD (Sigma, St. Louis, MO) was dissolved in 160 mM borate buffer (pH 9.1) at 5 mM as the BGE, and p-nitrophenol (Fisher, Fair Lawn, NJ) was dissolved in the same borate buffer at 2.0 mM. The capillary was filled with the BGE with β-CD, and pnitrophenol was injected into the capillary to form a very narrow analyte plug. The analyte was introduced into the capillary by applying a pressure of 0.5 psi (3447 Pa) for 3 s, and the volume was calculated to be 3.6 nL, with a plug length of 0.18 cm. This calculation was based on Beckman P/ACE injection parameters provided by the manufacturer. The CE process was operated at 10 kV, and a temperature of 20.0 °C was maintained throughout the experiments. Methanol was used as the electroosmotic flow (EOF) marker. β-CD is a neutral molecule and migrates at the same speed as the EOF. The mobility of p-nitrophenol in neat buffer (µA) has to be measured separately. In this measurement, the sample was injected into the capillary filled with neat borate buffer with no additive under a pressure of 0.5 psi (3447 Pa) for 3 s, and the migration was driven by a voltage of 15 kV. The measurements were repeated multiple times before carrying out the affinity CE experiment. β-CD has no UV absorbance in the range of 190-500 nm, and the complex is assumed to have a UV absorbance similar to the analyte. The binding constant of 533 M-1 for this 1:1 interaction is determined using nonlinear regression or linear regression methods.17 This number was used as one of the inputs of the simulation program. The results of the simulation, including peak shapes and migration times, were compared with the experimental data. Nonequilibrium CE experiments, similar to the NECEEM experiments reported by Berezovski et al.,4 are simulated with the same program. No real nonequilibrium CE experiments were (17) Bowser, M. T.; Chen, D. D. Y. J. Phys. Chem. A 1998, 102, 8063-8071.
Analytical Chemistry, Vol. 77, No. 3, February 1, 2005
841
carried out. Instead, the conditions emulating NECEEM experiments are input into the simulation programs to demonstrate the ability of this program to simulate all ACE processes. ALGORITHM AND IMPLEMENTATION The capillary is divided into thousands of narrow cells. Between adjacent cells, there is migrational flux of all species. At the same time, there are interactions between the additive and the analyte. The electrophoretic migration and the association/dissociation processes are considered separately: the changes in concentrations due to the electrophoretic migration are first calculated; then the association/dissociation processes of the analyte (A), the additive (P), and the complex (AP) are considered in order to calculate the final concentrations of all species in each cell at this moment before the next step of electrophoretic migration. Electrophoretic Migration. The electrophoretic migration process for both equilibrium and nonequilibrium CE is identical and can be described by the following equation:6, 7
∂Cz,t,i ∂Cz,t,i ∂2Cz,t,i ) - µiEz + Di ∂t ∂z ∂z2
(1)
where Cz,t,i is the concentration of species i at position z, and time t, Ez is the total local electric field strength, µi is the apparent mobility of the ion i (µi ) µep,i + µeo, where µep,i is the electrophoretic mobility of the ion i and µeo is the EOF), and Di is the diffusion coefficient of the ion. The local electric field strength is assumed constant throughout the capillary. With a selected time increment ∆t, which defines the step, and a selected space increment ∆z, which is the length of a cell, eq 1 can be evaluated with FDSs. Finite Difference Analysis of Hyperbolic Partial Differential Equations. In the absence of diffusion, eq 1 can be simplified into a hyperbolic partial differential equation:
∂Cz,t,i/∂t ) - µiEz∂Cz,t,i/∂z
(2)
Equation 2 can be evaluated with one of the following finite difference schemes:
Cz,t + ∆t,i - Cz,t,i Cz,t,i - Cz-∆z,t,i ) - µiEz ∆t ∆z Cz,t+∆t,i - Cz,t,i Cz + ∆z,t,i - Cz,t,i ) - µiEz ∆t ∆z Cz,t+∆t,i - Cz,t,i Cz + ∆z,t,i - Cz - ∆z,t,i ) - µiEz ∆t 2∆z Cz,t+∆t,i - (1/2)(Cz+∆z,t,i + Cz-∆z,t,i) ) ∆t Cz + ∆z,t,i - Cz-∆z,t,i - µiEz 2∆z Cz,t + ∆t,i - Cz,t-∆t,i Cz + ∆z,t,i - Cz-∆z,t,i ) - µiEz 2∆t 2∆z
(3)
Analytical Chemistry, Vol. 77, No. 3, February 1, 2005
|µiEz|
∆t e1 ∆z
(8)
The solution of a hyperbolic problem at a point depends on information within some domain of dependency to the past. Equation 8 states that the space increment (∆z) must not be smaller than the distance that ion i travels during ∆t. In other words, ∆z defines the boundary of the spatial region that is allowed to calculateCz,t + ∆t,i by the differencing schemes. Mathematically, the Lax-Friedrichs scheme (eq 6) and the staggered leapfrog scheme (eq 7) are more accurate than eqs 3 and 4, because associated errors for eqs 6 and 7 (proportional to ∆z2) are smaller than those for eqs 3 and 4 (proportional to ∆z), as ∆z approaches 0. Using higher orders of centered differencing on the right-hand side (RHS) of eq 6 or 7, greater accuracy can be achieved. However, if there are sharp gradients in the system, the use of eq 6 or 7 will lead to numerical dispersion or phase errors. Sharp gradients cannot be avoided either in the operation or the simulation of CE. For a rectangular injection profile, there are sharp gradients on both edges of the injection plug; while for a parabolic profile, there is a sharp gradient at the inlet of the capillary.20 In addition to the sharp gradients resulting from the injection, the difference in concentrations between two adjacent cells does not approach 0 if the length of cell does not approach 0; i.e., there is always a small gradient between two adjacent cells. Because numerical dispersion is a propagation error, it becomes more significant as the simulation proceeds. On the other hand, the smaller the length of the cells (∆z), the slower numerical dispersion becomes significant. Numerical dispersion in various FDSs has been demonstrated by Ermakov et al.11,12 In addition, eqs 6 and 7 do not describe the physical properties of CE systems correctly, because the equations describing CE systems are advective equations, not wave equations. Without diffusion, any species migrates in only one direction due to the applied electric field. The migrational flux (F), which is the amount of a species migrating out of a cell, is defined as
(4) (5)
(6)
Fz,t,i ) Cz,t,iµiEz
(9)
If a species migrates toward the detector (µiEz > 0), only the migrational fluxes of the previous cell and the current cell determine the new concentration in the current cell. Therefore, the forward-time backward-space scheme (eq 3) should be used in this case. Substituting eq 9 into eq 3, we have
(7)
Equations 3-7 are the forward-time backward-space scheme, forward-time forward-space scheme, forward-time centered-space scheme, Lax-Friedrichs scheme (an enhanced scheme over eq 5), and staggered leapfrog (centered-space centered-time) scheme, respectively.18,19 842
Stability and convergence are required of FDSs. The Courant condition, which is described by the following equation, must be satisfied in all five schemes to avoid amplitude errors.19
Cz,t + ∆t,i - Cz,t,i Fz,t,i - Fz-∆z,t,i )∆t ∆z
(10)
(18) Strikwerda, J. C. Finite Difference Schemes and Partial Differential Equations; Wadsworth & Brooks/Cole Advanced Books & Software: Pacific Grove, CA, 1989. (19) Press: W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in C, 2nd ed.; Cambridge University Press: New York, 1992.
If a species migrates away from the detector (µiEz < 0), eq 4 should be used to give
Fz + ∆z,t,i - Fz,t,i Cz,t + ∆t,i - Cz,t,i )∆t ∆z
(11)
Equations 3, 4, 10, and 11 are also called the first-order upwind difference schemes because the spatial differences are skewed in the “upwind” direction of µiEz. In most CE systems, because of the existence of a strong electroosmotic flow, all species migrate toward the detector. Thus, in the following discussion, only the equations describing the movements toward the detector (µiEz > 0), such as eq 10, are given. The first-order upwind difference schemes not only provide the “fidelity” to CE systems, which is lacking in other schemes, but also eliminate numerical dispersion. But the cost for stability and simplicity is a new type of errorsnumerical diffusion. In eqs 10 and 11, the ions always arrive only from one cell (∆z) away. But in reality, the ions from a distance µiEz∆t away arrive after a time interval ∆t. The Courant condition (eq 8) states the requirement of ∆z and ∆t. If |µiEz|∆t ) ∆z, the numerical diffusion error is eliminated. However, if |µiEz|∆t < < ∆z, numerical diffusion can cause large errors. This type of error can be reduced by using higher-order upwind difference schemes. The monotonic transport scheme is an efficient second-order upwind difference scheme, which is derived by considering the analytic transport equation.21 The equation for the case of µiEz > 0 is
(
[
)
Cz,t + ∆t,i - Cz,t,i µiEz µiEz∆t 1 )Cz,t,i - Cz-∆z,t,i + 1 ∆t ∆z 2 ∆z
]
(dz,t,i - dz-∆z,t,i) (12)
where dz,t,i (defined by eq 13) and dz-∆z,t,i are the concentration gradients at position z and z-∆z, respectively.
dz,t,i )
2(Cz,t,i - Cz-∆z,t,i)(Cz + ∆z,t,i - Cz,t,i) , (Cz + ∆z,t,i - Cz-∆z,t,i) if (Cz,t,i - Cz-∆z,t,i)(Cz + ∆z,t,i - Cz,t,i) > 0 (13)
dz,t,i ) 0, otherwise. With the consideration of both accuracy and efficiency, eq 12 is chosen as the algorithm for the hyperbolic part of the differential equation. Finite Difference Analysis of Parabolic Partial Differential Equations. The following diffusion equation is a parabolic partial differential equation.
∂Cz,t,i/∂t ) Di(∂2Cz,t,i/∂z2)
(14)
Equation 14 can be evaluated with the following three schemes: the first-order fully explicit scheme (eq 15), the first(20) Peng, X.; Chen, D. D. Y. J. Chromatogr., A 1997, 767, 205-216. (21) Hawley, J. F.; Smarr, L. L. Astrophys. J. Suppl. 1984, 55, 221-246.
order fully implicit scheme (eq 16), and the second-order CrankNicholson scheme (eq 17).
Cz,t + ∆t,i - Cz,t,i Cz + ∆z,t,i - 2Cz,t,i + Cz-∆z,t,i ) Di ∆t (∆z)2
(15)
Cz,t + ∆t,i - Cz,t,i Cz + ∆z,t+∆t,i - 2Cz,t + ∆t,i + Cz-∆z,t+∆t,i ) Di ∆t (∆z)2 (16) Cz,t + ∆t,i - Cz,t,i ) ∆t Di
(Cz + ∆z,t+∆t,i - 2Cz,t + ∆t,i+ Cz-∆z,t+∆t,i) + (Cz + ∆z,t,i - 2Cz,t,i + Cz-∆z,t,i) 2(∆z)
(17)
2
Equation 15 (used in the simulation model of Dose el al.8) is stable only for sufficiently small time increments (∆t), and the stability criterion is
2D∆t/(∆z)2 e 1
(18)
Both eqs 16 and 17 are unconditionally stable; that is, arbitrarily large time increments can be used. Equation 17 is more accurate than eq 16 because of its second-order expression on RHS. Although eqs 16 and 17 are better choices, eq 15 is used in this model for the following reasons. First, for most large molecules such as proteins, the diffusion coefficients are usually quite small (< 10-6cm2/s). The simulation results are already very good using eq 15 as shown later in this paper. Second, eq 15 can be solved in a way similar to eq 12, so that they can be easily incorporated together in a program evaluating both the hyperbolic and parabolic parts of the function for analyte migration. But in order to solve eqs 16 and 17, a set of simultaneous linear equations must be solved at each time step for Cz,t + ∆t,i. When µiEz > 0, the overall finite difference scheme for the evaluation of the partial differential equation for electrophoretic migration in CE is the combination of eqs 10 and 15 or the combination of eqs 12 and 15 as the following:
Cz,t + ∆t,i - Cz,t,i Fz,t,i - Fz-∆z,t,i )+ ∆t ∆z Cz + ∆z,t,i - 2Cz,t,i + Cz-∆z,t,i (19) Di (∆z)2
[
(
)
µiEz∆t µiEz Cz,t + ∆t,i - Cz,t,i 1 Cz,t,i - Cz-∆z,t,i + 1 )∆t ∆z 2 ∆z Cz + ∆z,t,i - 2Cz,t,i + Cz-∆z,t,i (dz,t,i - dz-∆z,t,i) + Di (20) (∆z)2
]
To achieve sufficient accuracy with eq 19, very small ∆z and ∆t have to be used, which requires a long computing time to finish a single simulation run. But with the more accurate eq 20, it is not necessary to use extremely small ∆z and ∆t, so that the Analytical Chemistry, Vol. 77, No. 3, February 1, 2005
843
simulation becomes much faster. The simulation results with both equations will be compared in the Results and Discussion section. Association and Dissociation Processes. At this point, the electrophoretic migration process in this step is finished, and the new concentrations of all species in each cell are obtained. In the next step, the association/dissociation processes are taken into account, and the concentrations at the newly established equilibrium are calculated. The concentrations of free A, free P, and AP after the electrophoretic migration process of the current step are denoted as [Afree]j, [Pfree]j, and [AP]j. And the concentrations after the association/dissociation process are denoted as [Afree] /j , [Pfree] /j , and [AP] /j . The equations for equilibrium CE and nonequilibrium CE are different. Equilibrium CE. In this case, the binding interaction is assumed fast, and the equilibrium can be reached within ∆t. Thus, the binding constant (K) is used in the calculations:
K ) [AP] /j /[Afree] /j [Pfree] /j
(21)
Rearrange eq 21 by replacing [AP] /j and [Pfree] /j with the total concentrations of A and P ([Atotal]j and [Ptotal]j), which are obtained from the calculations for the electrophoretic migration process and remain constant at the equilibrium stage.
K[Afree] /j 2 + (K[Ptotal]j - K[Atotal]j + 1)[Afree] /j [Atotal]j ) 0 (22) [Afree] /j can be easily solved from eq 22. Nonequilibrium CE. Because slower interactions are studied with this technique, we have to use rate constants instead of binding constants. Different binding interactions may have different mechanisms, and each mechanism needs its own set of equations. For a simple 1:1 interaction, A+P h AP, forward rate constant and reverse rate constant are k + 1 and k-1, respectively. The final concentrations of the additive and the analyte in cell j are
[Pfree] /j ) [Pfree]j - (k + 1[Afree]j[Pfree]j - k-1[AP]j)∆t (23) [Afree] /j ) [Afree]j - (k + 1[Afree]j[Pfree]j - k-1[AP]j)∆t (24) Once we have [Afree] /j and [Pfree] /j , the calculation cycle is finished for the current step (time t). We then move on to the next step (time t + ∆t) and do the same calculations again. Implementation of Circular Cell Arrangement. An efficient implementation is as important as developing appropriate mathematical algorithms in computer simulation. Conventionally, the entire capillary is equally divided into m very narrow cells. A linear array of m elements is allocated in the computer’s memory, and each element in this array stores the data for one cell. The way that the data for the divided cells are stored and manipulated in the computer’s memory is the key to the efficiency improvement of the simulation programs. In this work, a circular cell arrange844 Analytical Chemistry, Vol. 77, No. 3, February 1, 2005
Figure 1. Circular arrangement of cells in the memory of the PC. Cell 0 is the starting point. Cell j is the last cell that contains an appreciable amount of the analyte. Cell n is the last cell allocated in the memory, but when the analyte migrates pass this cell, it automatically enters and reuses cell 0.
ment and data manipulation algorithm is used to significantly reduce the computing time and the data storage. In equilibrium CE, the concentration in most cells does not change, because only the cells within or near the analyte plug have a noticeable amount of the analyte. Therefore, there is no need to track the changes in concentration for every cell during the electrophoretic migration of the analyte. Our program also divides the capillary into m cells. The cells within and near the analyte plug are labeled as active cells, in which there is a noticeable amount of analyte. All other cells are called inactive cells, in which the analyte concentration is less than a threshold value. Only the concentration changes in active cells are tracked and stored. At the beginning of an equilibrium CE experiments, the analyte exists only in the analyte plug. Due to the applied electric field, the analyte plug moves toward the detector; due to diffusion, the analyte spreads to the adjacent cells in both directions. To determine whether a cell is active, a userdefined threshold concentration of the analyte is used. When the concentration of the analyte in a cell is greater than the threshold, the cell is activated; when the concentration of the analyte in a cell is smaller than the threshold, the cell is deactivated. An array with n elements is allocated in the memory. Although n can be much smaller than m, it has to be large enough to store the data for all active cells. This array is arranged as a circle as shown in Figure 1, which makes it possible to reuse the cells. Two pointers (p1 and p2) are used with p1 pointing at the starting point of the active segment and p2 pointing at the end of the active segment. Only the elements within p1 and p2 are updated in the current step. At the beginning of a simulation, p1 points at element 0, and p2 points at element j (j ) Lplug/∆z), which is the last cell of the analyte plug. Any cell outside this range (from j + 1 to n) has the same default conditions: the concentration of the additive is [P]0 and the concentration of the analyte is 0. In the first step, the changes in the concentrations of all species for the cell from p1 (0) to p2 (j) are calculated using eqs 19-26. Some analyte moves into cell j + 1, so the pointer p2 is moved forward to j + 1 as well. At the same time, the concentration of analyte in the cell 0 decreases.
In every following step, the pointer p2 is moved one cell forward. But the pointer p1 still points at the first cell until the concentration of the analyte in the first cell is lower than the predefined threshold value. As the pointers p1 and p2 keep being moved forward, p2 may finally bypass the last physical element in the computer’s memory, n of the array. However, because the elements are arranged as a circle, logically there is no “last” element. The pointer p2 is moved to 0 again after n. This operation is valid also because the elements are released on the side of pointer p1 as the analyte moves forward. This new algorithm is very efficient for the simulation of equilibrium CE, because there is always just one analyte plug needed to be monitored. In nonequilibrium CE, the analyte band often becomes quite broad. Therefore, two or more significant portions of the analyte band are monitored in the electropherogram, and more cells are active at any given time, which leads to a longer computing time for running a simulation. Similar implementation strategy as described for the simulation of equilibrium CE is used to simulate nonequilibrium CE, which is an advantage of this new simulation model because it can be used on nearly any kind of CE application with increased speed. However, the initial conditions for the simulation of nonequilibrium CE are different: the equilibrated mixture of the additive and the analyte is injected into a capillary filled with a plain buffer. In CE applications for separation purposes, there are gaps between two adjacent, but completely separated peaks. If the cells within these gaps, having no analyte inside, are treated as inactive cells, the speed of the simulation can also be increased. In that case, more pointers (p3, p4, ...) are needed to track the start and the end of each active segment. RESULTS AND DISCUSSION Based on eq 19 and eq 20, two CE simulation programs were written in C++. Most of the discussion is based on the 1:1 equilibrium CE experiments. However, these properties are applicable to both equilibrium and nonequilibrium CE. Time and Space Increments. The time increment ∆t and the space increment (the length of cell) ∆z are the most important simulation parameters. The implementation of eq 19 is first tested with different ∆t and ∆z. Figure 2A shows the simulated peaks for the equilibrium CE experiment described in the Experimental Section (initial conditions: [P]0 ) 5 mM, [A]0 ) 2 mM) with the same ∆z (0.0005 cm) but different ∆t (0.0005-0.002 s). The three simulated curves only have small differences that can be ignored. Thus, there is no need to use extremely small values for ∆t. The space increment ∆z plays a more important role than the time increment ∆t in predicting peak shapes with eq 19. Figure 2B shows the simulated peaks with the same ∆t (0.0005 s) but different ∆z (0.0001-0.003 cm). With a larger ∆z, the numerical diffusion is so significant that the simulated peaks are broad and very different from the real peak. As ∆z becomes smaller, the simulated peaks become narrower and closer to the real one. However, even if we can afford to spend longer computing time to integrate using extremely small ∆t and ∆z, round-off errors, originated from machine errors, will accumulate and eventually become large enough to invalidate the solution. The implementation of eq 20 is then tested in a similar way. Because the numerical diffusion is greatly reduced with eq 20,
Figure 2. Simulated peaks with the implementation of eq 19. (A) Same ∆z and different ∆t. (B) Different ∆z and same ∆t.
Figure 3. Simulated peaks with the implementation of eq 20. Even though different simulation parameters are used, the generated peaks have a very similar shape.
no extremely small ∆t and ∆z are required. The simulated peaks shown in Figure 3 are all very similar even though very different values of ∆t and ∆z are used. The simulated peaks generated using eq 19 and eq 20 are plotted together with the real peak in Figure 4. All three peaks have very similar shapes: a steep edge on the right and a gentle Analytical Chemistry, Vol. 77, No. 3, February 1, 2005
845
Figure 4. Experimental and simulated peaks. The y-axis for the simulated peaks is the concentration of the analyte, and the y-axis for the experimental peak is the absorption (mAu). In this figure, all peaks are normalized to the same scale. The experimental conditions for generating this experimental peak are described in the Experimental Section.
slope on the left. The cause of this kind of shape has been discussed in another paper.22 However, the real peak is a little broader (having larger variance) than the simulated peaks. The total peak variance is the sum of the variances caused by injection, longitudinal diffusion, and other factors.20 The variance caused by diffusion has been considered in the simulation program. Treating the parabolic profile caused by pressure injection as a rectangle and the inability to measure the length of the injection plug directly lead to an uncertainty in the variance caused by injection. In addition, other variances, such as band broadening caused by the length of the detector window and wall adsorption, are not considered. Therefore, the simulated peaks usually have smaller variances. The simulated peak positions (1389 s from eq 19 and 1390 s from eq 20) are slightly different from the real peak position (1386 s), ∼0.3%, mainly due to the errors associated with the measurements of the electrophoretic mobility of the free analyte (µep,A) and the viscosity correction factor (ν). The speed of the simulation is determined by ∆z and ∆t. The computing time for generating the peak with ∆t ) 0.0001 s and ∆z ) 0.0001 cm using eq 19 is more than 8 h with a laptop PC powered by an Intel Pentium III 866. On the other hand, the computing time for generating the peak with∆t ) 0.001 s and ∆z ) 0.001 cm using eq 20 is only 7 min with the same PC. It is clear that eq 20 is superior to eq 19, because rather large increments can be used with eq 20 to achieve similar accuracy compared to much smaller increments used with eq 19. The speed of the simulation is improved close to 2 orders of magnitude. Active/Inactive Cells. Another important step of the simulation is to determine which cell should be activated or deactivated. This is achieved by using threshold values on both inlet and outlet sides of the capillary as described in the previous section. These threshold values play important roles in the simulation. If the threshold values are too large, a cell with a significant amount analyte could be deactivated if its total concentration of analyte is (22) Fang, N.; Ting, E.; Chen, D. D. Y. Anal. Chem. 2004, 76, 1708-1714.
846 Analytical Chemistry, Vol. 77, No. 3, February 1, 2005
Figure 5. Simulated peaks for nonequilibrium CE experiments. All three figures are generated with identical experimental conditions and simulation parameters except the dissociation constant (k-1). (A) k-1 ) 0.01 s-1. Peak 1 corresponds to the free analyte under the initial equilibrium conditions. The exponential part 2 corresponds to the dissociation of the complex, and the small peak 3 corresponds to the remaining complex. (B) k-1 ) 0.05 s-1. In this case, there is no complex peak because the complex has been completely dissociated. (C) k- ) 0.001 s-1. In this case, the dissociation process is so slow that the exponential part 2 is not visible.
smaller than the threshold values, which leads to large errors in the simulation. If the threshold values are too small, a large
number of cells remain active for a long time, which decreases the efficiency of the simulation program. The acceptable choices of threshold values are usually 1/10 000-1/1 000 000 of the concentration of injected analyte solution. Diffusion Coefficients. The longitudinal diffusion of all species in the capillary contributes to the peak shape in CE experiments. The more significant the diffusion is, the broader the final peak is, which leads to more active cells and longer computing time. According to the Stokes-Einstein equation for isolated hard spheres, the diffusion coefficient depends on the size and the shape of the molecule, the interaction with the solvent, and the viscosity of the solvent. The diffusion coefficients for proteins in aqueous solution are usually on the order of 10-6 cm2 s-1 or smaller. The influence of the longitudinal diffusion on the peak shape is not significant for studying protein binding interactions. For all systems presented in this paper, the diffusion coefficients are not measured experimentally, instead, estimated values are given to the simulation program. Simulation of Nonequilibrium CE. The simulation model can also be used to simulate nonequilibrium CE. However, no real experimental data are compared to the simulation results in this section. Instead, the characteristics of Berezovski et al.’s NECEEM experiments 4 will be demonstrated by simulation. The conditions given to the program are as follows: the binding constant (Ka ) 2.0 × 107 M-1), the dissociation rate constant (k-1 ) 0.01 s-1 in Figure 5A; k-1 ) 0.05 s-1 in Figure 5B; k-1 ) 0.001 s-1 in Figure 5C), the initial concentration of the additive (2 µM), the initial concentration of the analyte (61 nM), the length of the capillary (33 cm), the length to detector (33 cm), the length of the injected plug (0.036 cm), the applied voltage (13200 V), the diffusion coefficients (10-5 cm2 s-1 for all species), the mobility of free analyte (0.000 357 cm2 V-1 s-1), and the complex mobility (0.000 146 cm2 V-1 s-1). The association rate is set to 0, because the re-formation of the complex is completely prevented as in Berezovski et al.’s NECEEM experiments.4 In addition to the experimental conditions, the simulation parameters are set as follows: the time increment (∆t ) 0.003
s), the length of the cell (∆z ) 0.003 cm), and the concentration threshold (10-14 M). Figure 5 shows the simulated peaks under the given experimental conditions and simulation parameters. The peaks 1 and 3 are sharper than normal experimental peaks due to the incomplete consideration of various causes of peak variance. Because the entire peak region in nonequilibrium CE is much longer than that in equilibrium CE, the number of active cells in nonequilibrium CE is also much larger. The computing time for the simulation of nonequilibrium CE is usually longer than the simulation of equilibrium CE, if the same ∆z and ∆t values are used. CONCLUSION The monotonic transport scheme is used to solve the hyperbolic part of the partial differential equation. With this implementation, numerical dispersion is avoided, and numerical diffusion is greatly reduced. With the implementation of this finite difference scheme, the time and space increments do not need to be extremely small to achieve higher accuracy. At the same time, a circular cell arrangement and manipulation algorithm is implemented to avoid unnecessary calculations and provide an efficient use of the computer’s memory. This new CE simulation model can simulate both equilibrium and nonequilibrium CE accurately in a reasonably short time. The simulation can be finished in a few minutes to a few hours depending mainly on the simulation parameters. The high efficiency of this model makes it possible to bridge the gap between theoretical modeling and real CE experiments and to predict the peak shape and peak position at any given time of a CE process. ACKNOWLEDGMENT The authors thank Professors Mark Thachuk, Lionel G. Harrision, and Y. Alexander Wang of Chemistry Department, UBC, for helpful discussions. The Beckman Coulter ProteomeLab PA800 automated CE system was kindly lent to us by Beckman Coulter Inc., Fullerton, CA. Received for review July 20, 2004. Accepted October 27, 2004. AC048938O
Analytical Chemistry, Vol. 77, No. 3, February 1, 2005
847
