Anal. Chem. 1998, 70, 1066-1075
Three-Dimensional Molecular Simulation of Electrophoretic Separations Daniel L. Hopkins and Victoria L. McGuffin*
Department of Chemistry and Center for Fundamental Materials Research, Michigan State University, East Lansing, Michigan 48824-1322
A three-dimensional stochastic simulation has been developed for electrophoretic separations. This simulation follows the trajectories of individual molecules by applying the fundamental equations of motion for diffusion, electroosmotic flow, and electrophoretic migration. The molecular zone profile may then be examined and characterized at any specified time or distance. The stochastic simulation has been validated by comparison with classical transport models, having average relative errors for the mean zone distance and variance of 0.01-0.04% and 2.67-4.36%, respectively. The simulation also shows excellent agreement with experimental results for the electrophoretic separation of the nucleotide monophosphates in phosphate buffer solution, having average relative errors for the mean zone distance and variance of 1.65% and 12.99%, respectively. This stochastic simulation provides an accurate and versatile means to examine and characterize transport processes in complex electrophoretic systems. Electrophoresis is an important analytical technique for the separation of strong and weak electrolytes as well as large polyelectrolytes of biological and industrial importance. Our understanding of the fundamental processes that comprise electrophoresis has been achieved through comparison of experimental results with the predictions of theoretical models. The classical models used to describe mass transport in electrophoretic separations are usually developed to address the conditions of a specific mode of separation, such as zone electrophoresis (ZE),1-6 moving boundary electrophoresis (MBE),7 isotachophoresis (ITP),5,8-12 and isoelectric focusing (IEF).13-17 However, unified models that can be used to describe all of the electrophoretic
separation modes have also been developed.18-21 Descriptions of these and other electrophoretic models can be found in several reviews.22-24 These classical transport models require the formulation of detailed differential equations of mass and charge balance, together with the associated boundary conditions. The solution of these equations yields the concentration profile as a function of time and/or the steady-state concentration profile in the cases of ITP10,11 and IEF.15-17 A comprehensive model must include all of the physical and chemical processes that contribute to the migration and dispersion of the solute zone in the electrophoretic system. In practice, however, these comprehensive differential equations cannot be solved analytically in closed form and are generally solved by numerical methods after simplifying assumptions are invoked.25 For example, some models may consider electrophoretic migration alone,6,26 whereas others may include additional processes such as molecular diffusion and electrodiffusion,27 convection,4,5,28 chemical equilibria,12 and various combinations thereof.1,2,29 As a consequence of these simplifying assumptions, the classical models do not provide a completely general description of transport phenomena in electrophoretic systems. Furthermore, because the differential equations are derived from macroscopic principles and properties, these models are inadequate to describe transport in systems that are heterogeneous at the microscopic and molecular levels.
(1) Reijenga, J. C.; Kenndler, E. J. Chromatogr. A 1994, 659, 403-415. (2) Reijenga, J. C.; Kenndler, E. J. Chromatogr. A 1994, 659, 417-426. (3) Cann, J. R.; Goad, B. W. J. Biol. Chem. 1965, 240, 1162-1164. (4) Ermakov, S. V.; Righetti, P. G. J. Chromatogr. A 1994, 667, 257-270. (5) Dose, E. V.; Guiochon, G. A. Anal. Chem. 1991, 63, 1063-1072. (6) Mikkers, F. E. P.; Everaerts, F. M.; Verheggen, T. P. E. M. J. Chromatogr. 1979, 169, 1-20. (7) Cann, J. R.; Goad, B. W. J. Biol. Chem. 1965, 240, 148-155. (8) Vacik, J.; Fidler, V. In Analytical Isotachophoresis; Everaerts, F. M., Ed.; Elsevier: Amsterdam, The Netherlands, 1981; pp 19-24. (9) Moore, G. T. J. Chromatogr. 1975, 106, 1-16. (10) Shimao, K. Electrophoresis 1986, 7, 121-128. (11) Shimao, K. Electrophoresis 1986, 7, 297-303. (12) Gebauer, P.; Bocek, P. J. Chromatogr. 1984, 299, 321-330. (13) Cann, J. R.; Stimpson, D. I. Biophys. Chem. 1977, 7, 103-114. (14) Stimpson, D. I.; Cann, J. R. Biophys. Chem. 1977, 7, 115-119.
(15) Palusinski, O. A.; Allgyer, T. T.; Mosher, R. A.; Bier, M.; Saville, D. Biophys. Chem. 1981, 13, 193-202. (16) Bier, M.; Mosher, R. A.; Palusinski, O. A. J. Chromatogr. 1981, 211, 313335. (17) Shimao, K. Electrophoresis 1987, 8, 14-19. (18) Bier, M.; Palusinski, O. A.; Mosher, R. A.; Graham, A.; Saville, D. A. In Electrophoresis ‘82; Stathokos, D., Ed.; Walter de Gruyter & Co.: Berlin, Germany, 1983; pp 51-60. (19) Bier, M.; Palusinski, O. A.; Mosher, R. A.; Saville, D. A. Science 1983, 219, 1281-1287. (20) Saville, D. A.; Palusinski, O. A. AIChE J. 1986, 32, 207-214. (21) Palusinski, O. A.; Graham, A.; Mosher, R. A.; Bier, M.; Saville, D. A. AIChE J. 1986, 32, 215-223. (22) Thormann, W.; Mosher, R. A. In Advances in Electrophoresis, Vol. 2; Chrambach, A., Dunn, M. J., Radola, B. J., Eds.; VCH Publishers: Weinheim, Germany, 1988; pp 47-108. (23) Mosher, R. A.; Saville, D. A.; Thormann, W. The Dynamics of Electrophoresis; VCH Publishers: Weinheim, Germany, 1992. (24) Hjerten, S. Electrophoresis 1990, 11, 665-690. (25) Ermakov, S. V.; Bello, M. S.; Righetti, P. G. J. Chromatogr. A 1994, 661, 265-278. (26) Thormann, W.; Mosher, R. A. Electrophoresis 1985, 6, 413-418. (27) Boyack, J. R.; Giddings, J. C. J. Biol. Chem. 1960, 235, 1970-1972. (28) Biscans, B.; Bertrand, J. Chem. Eng. Res. Des. 1987, 65, 224-230. (29) Andreev, V. P.; Lisin, E. E. Chromatographia 1993, 37, 202-210.
1066 Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
S0003-2700(97)00802-0 CCC: $15.00
© 1998 American Chemical Society Published on Web 02/03/1998
A completely different approach, which may be called stochastic (or Monte Carlo) molecular simulations, can be utilized to describe mass transport in separation systems. These simulations are based on the migration of individual molecules by the sequential application of independently defined processes (Markov chain). These simulations range in complexity from onedimensional models with fixed step size30-32 to three-dimensional models with step size distributions.33-37 In the most rigorous and comprehensive treatments, the fundamental equations of motion for diffusion, convection, and other transport processes are applied to the individual molecules. In contrast to the complex differential equations of mass and charge balance, these equations of motion are relatively simple and require few, if any, simplifying assumptions. However, extremely powerful computers are required in order to model the motion of a statistically significant number of molecules through meaningful intervals of time or distance. This problem has been overcome with the routine availability of highspeed, high-memory computers, particularly those with parallel processors.38 In recent years, three-dimensional molecular simulations have been applied to the study of flow injection analysis,33,34 field-flow fractionation,35,39 gas and liquid chromatography,40,41 and capillary electrophoresis.36 In the present study, we describe the development of a threedimensional molecular simulation of electrophoretic separations. Algorithms describing diffusion, electroosmotic flow, and electrophoretic migration have been implemented in a sequential manner, such that each may be independently validated by comparison with classical transport models. All physical and chemical parameters that are routinely used to control these transport phenomena are invoked as intrinsic and independent variables. Hence, this simulation provides a powerful and versatile means to examine and characterize the complex behavior of solute molecules in electrophoretic separations. COMPUTER SIMULATION Simulation Program. The computer simulation is written in the FORTRAN 90 programming language and optimized for execution on an IBM RS/6000 model 580 computer. A simplified flowchart of the program is shown in Figure 1. The program incorporates algorithms for diffusion, electroosmotic flow, and electrophoretic migration, which are applied to each molecule at each time increment until the total simulation time is reached. The simulation may be performed in Cartesian global coordinates (X, Y, Z), which are most appropriate for electrophoretic separa(30) Weber, S. G. Anal. Chem. 1984, 56, 2104-2109. (31) Giddings, J. C. Dynamics of Chromatography; Marcel Dekker: New York, 1965. (32) Giddings, J. C. J. Chem. Educ. 1958, 35, 588-591. (33) Betteridge, D.; Marczewski, C. Z.; Wade, A. P. Anal. Chim. Acta 1984, 165, 227-236. (34) Wentzell, P. D.; Bowridge, M. R.; Taylor, E. L.; MacDonald, C. Anal. Chim. Acta 1993, 278, 293-306. (35) Schure, M. R. Anal. Chem. 1988, 60, 1109-1119. (36) Schure, M. R.; Lenhoff, A. M. Anal. Chem. 1993, 65, 3024-3037. (37) McGuffin, V. L.; Wu, P.; Hopkins, D. L. In Proceedings of the International Symposium on Chromatography; Hatano, H., Hanai, T., Eds.; World Scientific Publishing: River Edge, NJ, 1995; pp 45-69. (38) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, England, 1989. (39) Schure, M. R.; Weeratunga, S. K. Anal. Chem. 1991, 63, 2614-2626. (40) Guell, O. A.; Holcombe, J. A. Anal. Chem. 1990, 62, 529A-542A. (41) McGuffin, V. L.; Wu, P. J. Chromatogr. A 1996, 722, 3-17.
tions in planar media, or, alternatively, in cylindrical global coordinates (R, Θ, Z) for separations in capillary tubes or fibers. Similarly, the local frame for each molecule may be expressed in Cartesian coordinates (x, y, z) or spherical coordinates (F, φ, θ). Because of their mathematical simplicity, the cylindrical global frame and the spherical molecular frame will be described in detail here. Simulation Input. The input parameters for the simulation are divided into three classes, as shown in Table 1. The system parameters include the properties of the fluid and surface phases, the electrical conditions, and the spatial dimensions of the system. The molecular parameters describe the physical and chemical attributes of the solute molecules, and the computational parameters describe certain constraints that are required for the simulation. On the basis of these input parameters, arrays are created to contain the properties and coordinates of the molecules. To initialize the simulation, the molecules are distributed with a delta, rectangular, or Gaussian axial profile of specified variance at a specified mean distance in the global coordinate frame. Simulation Output. The program allows the molecular zone profile to be examined as the distance distribution at fixed time(s) or as the time distribution at fixed distance(s). The latter case may also be used to simulate the signal from a detector at a fixed distance with a window or transducer of finite length. The molecular distributions in distance or time are characterized by means of the statistical moments.42 For example, the first statistical moment, or mean distance, Zh , is calculated as N
∑Z
Z h ) N-1
(1)
i
i)1
where Zi is the axial coordinate of an individual molecule and N is the total number of molecules. The second statistical moment, or variance, σ2, is calculated as N
σ2 ) N-1
∑(Z - Zh )
2
i
(2)
i)1
and the higher-order statistical moments are determined in a similar manner. These statistical moments and other quantitative figures of merit for the electrophoretic separation are stored in data files at specified intervals during the simulation. In addition, the molecular population is summed in discrete segments of distance or time and then smoothed by Fourier transform methods43 to provide a continuous zone profile for graphical display. Examples of the segmented and smoothed profiles are shown in Figure 2. Because the molecular distributions can be examined at any distance or time, these output routines provide an extensive numerical and visual record of the transport processes throughout the simulation. (42) Grushka, E.; Myers, M. N.; Schettler, P. D.; Giddings, J. C. Anal. Chem. 1969, 41, 889-892. (43) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: Cambridge, England, 1989.
Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
1067
Figure 1. Flowchart of the computer simulation program.
DESCRIPTION AND VALIDATION OF TRANSPORT ALGORITHMS Assumptions. To facilitate the development of algorithms for the transport processes, the following assumptions have been invoked: (1) The molecules are represented as point masses and point charges. (2) The molecules are independent and do not interact with one another. (3) The properties of the fluid and surface phases are assumed to be constant in time and constant in all spatial dimensions. (4) The electric field is applied along the axial dimension only. The electric field gradient is assumed to be constant in time and constant in the axial spatial dimension. 1068 Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
(5) The passage of current through the fluid and surface phases does not cause significant Joule heating. The temperature is assumed to be constant in time and constant in all spatial dimensions. (6) The transport processes are assumed to be independent and fully developed at all times (no induction period). (7) The chemical reactions (e.g., acid-base and complexation) are assumed to be rapid and to attain equilibrium within the time increment of the simulation. The models of the individual transport processes may contain other assumptions that are inherent in their mathematical development (vide infra). Many of these assumptions are not intrinsic limitations of the stochastic simulation but, rather, are necessary
Table 1. Input Parameters for the Computer Simulation Program parameter
symbol
System Parameters radius of the system length of the system position of the detection zone length of the detection zone position of the injection zone length (variance) of the injection zone zeta potential of the surface pH of the fluid ionic strength of the fluid viscosity of the fluid dielectric constant of the fluid temperature voltage
R0 L Ldet ldet Linj linj (σ2inj) ζ pH I η T0 V
Molecular Parameters diffusion coefficient of species i equilibrium constant of species i electrophoretic mobility of species i charge of species i
Di Ki µi z(i
Figure 3. Comparison of the radial probability distribution PF for simulated three-dimensional diffusion with theoretical prediction (s) according to eq 4.
N t T
a constant step size in one or three dimensions.31,32,36,39,41,44,45 This approach has been extended by varying the step size about a nonzero average value.30,40 Although this extension allows for some statistical variation, it does not provide an accurate distribution of step sizes. Betteridge et al.33 and Wentzell et al.34 used a three-dimensional approach where the diffusion step in each direction of the Cartesian molecular frame (x, y, z) was chosen randomly from a uniform distribution in the range (2(2Dt)1/2. As in the previous case, this approach provides the correct average step size but does not provide an accurate distribution of step sizes. The molecular diffusion algorithm incorporated in the present simulation is based on a Gaussian (normal) distribution integrated through three-dimensional space in spherical coordinates (F, φ, θ)
Computational Parameters number of molecules time increment total simulation time molecular coordinate systems spherical coordinates Cartesian coordinates global coordinate systems cylindrical coordinates Cartesian coordinates
F, φ, θ x, y, z R, Θ, Z X, Y, Z
∫ ∫ ∫ (2π)1 2π
0
π
0
∞
0
1/2
σ
exp
( )
-F2 F sin θ dF dφ dθ 2σ2
(3)
By integration of the radial component of eq 3, the probability distribution PF is expressed as
PF )
Figure 2. Segmented and smoothed molecular distributions in time at fixed distance (top) and in distance at fixed time (bottom). Simulation conditions: N ) 500; t ) 1.0 × 10-4 s; T ) 25.0 s; R0 ) 2.5 × 10-3 cm; L ) 2.5 × 101 cm; D ) 1.0 × 10-4 cm2 s-1; v0 ) 1.0 cm s-1.
in the present study for validation and comparison to classical transport theory. Diffusion. Diffusion arises from the random Brownian motion of molecules and is often modeled as a random-walk process with
( )
F2 -F2 exp 1/2 (2π) σ 2σ2
(4)
This theoretical probability distribution is illustrated in Figure 3. The variance is calculated from the Einstein-Smoluchowski equation44 as σ2 ) 2Dt, where D is the diffusion coefficient and t is the time increment. Under these conditions, the mean, rootmean-square, and most probable values of F from this distribution are (16Dt/π)1/2, (6Dt)1/2, and (4Dt)1/2, respectively. It is noteworthy that this root-mean-square value is equal to the constant step size that has been used successfully in previous threedimensional stochastic simulations.41 To implement the diffusion algorithm, the radial distance F traveled by each molecule during the time increment t is randomly (44) Einstein, A. Ann. Phys. 1905, 17, 549. (45) Chandrasekhar, S. Rev. Mod. Phys. 1943, 15, 1-89.
Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
1069
Electroosmotic Flow. Electroosmotic flow is induced when an electric field is applied tangential to a charged surface in contact with a fluid phase. The axial distance z traveled in the time increment t is given by
z ) vt
(5)
where v is the velocity of the molecule at a specified radial position. The radial profile of electroosmotic flow has been modeled as flat in previous stochastic simulations,36 which is a reasonable approximation for many experimental conditions. However, the assumption of a flat profile may lead to error when the fluid phase has low ionic strength or when the molecules have small diffusion coefficients. The electroosmotic flow algorithm incorporated in the present simulation is based on the radially dependent flow profile given by the Rice-Whitehead equation,46
Figure 4. Validation of the diffusion algorithm by comparison of the mean zone distance and variance with theoretical predictions. Simulation conditions: N ) 500; t ) 1.0 × 10-4 s; R0 ) 2.5 × 10-3 cm; D ) 1.0 × 10-4 (3), 1.0 × 10-5 (]), 1.0 × 10-6 (4), 1.0 × 10-7 (0), and 1.0 × 10-8 cm2 s-1 (O). Theory (s): Z h ) 0; σ2 ) 2DT.
selected from the probability distribution given in eq 4, and the direction is subsequently randomized through the spherical coordinate angles (φ, θ). For simplicity and computational speed, the radial integration limits are taken as 0-4σ rather than 0-∞. This includes >99.9% of the integrated Gaussian distribution, which suggests that less than 1 in 1000 diffusion steps would be expected statistically to be greater than the allowed range. As shown in Figure 3, the distribution of F for 50 000 molecules that undergo simulated three-dimensional diffusion agrees well with the theoretical distribution. The randomly selected coordinate increments in the spherical molecular frame are used to calculate the new position of each molecule in the cylindrical global frame. Because diffusion provides radial in addition to axial transport, the possibility must be considered that a molecule will encounter the surface. Other stochastic simulations have addressed this situation by returning the molecule to its original position,33 by repositioning the molecule at the point of intersection with the surface,39 or by reflecting the molecule along the incident trajectory.34 Another approach, called infinite radial diffusion, completely eliminated encounters with the surface by randomizing the radial position of each molecule at each time increment.34 In the present simulation, the molecule undergoes an elastic collision when it intersects the surface.41 This approach preserves both the distance and the angle of the diffusion step and is a physically meaningful description of the molecule-surface encounter. To verify the accuracy of the diffusion algorithm, the mean zone distance and variance for an ensemble of 500 molecules were monitored as a function of the simulation time. These results were compared to theoretical predictions based on the EinsteinSmoluchowski equation.44 As shown in Figure 4, excellent agreement is observed for diffusion coefficients in the range from 10-4 to 10-8 cm2 s-1. The average relative errors for the zone distance and variance are 0.01% and 2.73%, respectively. 1070 Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
(
v ) v0 1 -
I0(κR) I0(κR0)
)
(6)
where R is the radial coordinate in the cylindrical global frame, R0 is the radius, v0 is the maximum velocity, κ-1 is the Debye length, and I0 is the zero-order modified Bessel function of the first kind. The maximum velocity may be specified as an input parameter or may be calculated from the Helmholtz-Smoluchowski equation47,48 as
v0 )
0ζV ηL
(7)
where and 0 are the permittivity of the fluid phase and vacuum, respectively, η is the viscosity of the fluid phase, and V is the applied voltage. The zeta potential ζ can be measured experimentally or estimated by means of a semiempirical relationship.49,50 The Debye length is a measure of the depth of the electrical double layer at the surface and is given by
κ-1 )
0kT0 2Ce2z(2
(8)
where e is the electron charge, k is the Boltzmann constant, T0 is the absolute temperature, and C and z( are the concentration (number density) and charge, respectively, of the electrolyte in the fluid phase. The electrokinetic radius, κR0, is the ratio of the capillary radius to the Debye length and indicates the extent of curvature in the electroosmotic flow profile. The flow profiles derived from simulations of electroosmotic flow for electrokinetic radii of 5, 10, 20, 50, and 100 are compared to the theoretical profiles from the Rice-Whitehead equation in Figure 5. To verify the accuracy of the electroosmotic flow algorithm, the mean zone distance and variance for an ensemble of 500 (46) Rice, C. L.; Whitehead, R. J. Phys. Chem. 1965, 69, 4017-4024. (47) von Smoluchowski, M. Bull. Int. Acad. Sci. Cracovic 1903, 1903, 184. (48) Huckel, E. Phys. Z. 1924, 25, 204. (49) Tavares, M. F. M.; McGuffin, V. L. Anal. Chem. 1995, 67, 3687-3696. (50) Hayes, M. A.; Ewing, A. G. Anal. Chem. 1992, 64, 512-516.
V L
v)µ
Figure 5. Comparison of the simulated electroosmotic flow profiles with theoretical prediction (s) according to McEldoon and Datta51 for electrokinetic radii κR0 ) 5 (O), 10 (0), 20 (4), 50 (]), and 100 (3).
where µ is the electrophoretic mobility. The mobilities are usually reported in the literature at infinite dilution; hence, they are corrected by means of the modified Onsager equation53 to the specified ionic strength of the fluid phase in the simulation. If the molecule exists as a single species (e.g., strong electrolyte), then the mobility is constant. In this case, the algorithm provides constant axial displacement for all molecules during each time increment according to eqs 5 and 9. This axial coordinate increment in the molecular frame is used to calculate the new position for each molecule in the global frame. To verify the accuracy of this algorithm, the mean zone distance and variance for an ensemble of 500 molecules were monitored as a function of the simulation time. As shown in Figure 7, excellent agreement was observed between the computer simulation and theoretical predictions for electrophoretic mobilities in the range from +10-3 to -10-3 cm2 V-1 s-1. The average relative errors for the zone distance and variance are 0.04% and 2.67%, respectively. If the molecule exists as n multiple species in dynamic equilibrium (e.g., weak electrolyte), then the mobility of an individual molecule is determined from statistical probability at each time increment. The fraction Ri of each species i is calculated from the appropriate equilibrium constants,54,55 which are corrected for ionic strength by means of the Davies equation.56 The identity of a molecule is determined by selecting a random number, ξ, between 0 and 1 to establish the value of i that satisfies the relationship n-i
1-
∑ j)1
Figure 6. Validation of the electroosmotic flow algorithm by comparison of the mean zone distance and variance with theoretical predictions. Simulation conditions: N ) 500; t ) 1.0 × 10-4 s; R0 ) 2.5 × 10-3 cm; L ) 1.0 × 102 cm; V ) 2.0 × 104 V; D ) 1.0 × 10-5 cm2 s-1; pH ) 7.0; I ) 1.0 × 10-7 (]), 1.0 × 10-6 (4), 1.0 × 10-4 (0), and 1.0 × 10-2 M (O). Theory (s): Z h ) (0ζVT/ηL)(1 - C1); σ2 ) 2DT + C2R02Z h 2/DT, where C1 ) 2I1(κR0)/κR0I0(κR0) and C2 ) (2C12/ (1 - C1)2)(3/8 + 2/(κR0)2 - 1/C1(κR0)2 - 1/C12(κR0)2).
molecules were monitored as a function of the simulation time. These results were compared to theoretical predictions based on the analytical solution of the Rice-Whitehead equation by McEldoon and Datta.51 As shown in Figure 6, excellent agreement is observed for fluid phases with ionic strength in the range from 10-2 to 10-7 M. The average relative errors for the zone distance and variance are 0.01% and 4.36%, respectively. Electrophoretic Migration. For molecules under the influence of an applied electric field,52 the velocity of electrophoretic migration is given by (51) McEldoon, J. P.; Datta, R. Anal. Chem. 1992, 64, 227-230. (52) Atkins, P. W. Physical Chemistry; Oxford University Press: Oxford, England, 1978.
(9)
i+1
Rn-j < ξ e
∑R
j-1
(10)
j)1
The molecule is then assigned the mobility µi corresponding to species i during that time increment, and its electrophoretic migration is subsequently calculated via eqs 5 and 9. This axial coordinate increment in the molecular frame is used to calculate the new position of the molecule in the global frame. To verify the accuracy of this algorithm, phosphate was selected as a well-characterized model system. The equilibrium constants and electrophoretic mobilities of the individual species H3PO4, H2PO4-, HPO42-, and PO43- are summarized in Table 2.57 The mean zone distance and variance for an ensemble of 500 molecules were monitored as a function of the simulation time. As shown in Figure 8, excellent agreement is observed between the computer simulation and theoretical predictions for phosphate in the pH range from 3.0 to 9.0. The average relative errors for the zone distance and variance are 0.01% and 3.38%, respectively. (53) Robinson, R. A.; Stokes, R. H. Electrolyte SolutionssThe Measurement and Interpretation of Conductance, Chemical Potential and Diffusion in Solutions of Simple Electrolytes; Butterworth: London, England, 1959. (54) Butler, J. N. Ionic Equilibrium: A Mathematical Approach; AddisonWesley: Reading, MA, 1964. (55) Laitinen, H. A.; Harris, W. E. Chemical Analysis; McGraw-Hill: New York, 1975. (56) Davies, C. W. J. Chem. Soc. 1938, 1938, 2093-2098. (57) Hirokawa, T.; Kobayashi, S.; Kiso, Y. J. Chromatogr. 1985, 318, 195-210.
Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
1071
Figure 7. Validation of the electrophoretic migration algorithm for single species by comparison of the mean zone distance and variance with theoretical predictions. Simulation conditions: N ) 500; t ) 1.0 × 10-4 s; R0 ) 2.5 × 10-3 cm; L ) 1.0 × 102 cm; V ) 2.0 × 104 V; v0 ) 0.0 cm s-1; D ) 1.0 × 10-5 cm2 s-1; µ ) +1.0 × 10-3 (]), +1.0 × 10-4 (4), -1.0 × 10-4 (0), and -1.0 × 10-3 cm2 V-1 s-1 (O). Theory (s): Z h ) µVT/L; σ2 ) 2DT. Table 2. Electrophoretic Mobilities and Equilibrium Constants for the Individual Phosphate Species at Infinite Dilution57 species
charge
mobility (cm2 V-1 s-1)
equilibrium constant
H3PO4 H2PO4HPO42PO43-
0 -1 -2 -3
0.00 × 10-4 -3.42 × 10-4 -5.91 × 10-4 -7.15 × 10-4
6.90 × 10-3 6.21 × 10-8 4.73 × 10-13
Effect of Simulation Parameters. Among the input parameters shown in Table 1, the computational parameters serve as fundamental constraints of the simulation. The most important of these parameters are the number of molecules, the time increment, and the total simulation time. These parameters control the computational time required for the simulation but also influence the accuracy and precision of the results. To illustrate the effect of these parameters, a series of simulations was performed under standardized conditions, and the results are summarized in Table 3. As the number of molecules is increased from 100 to 10 000, the computational time increases linearly. The accuracy, expressed as the relative error of the mean zone distance and variance, improves with increasing number of molecules. It is noteworthy that the relative error of the mean zone distance improves approximately as the square root of the number of molecules, which is consistent with statistical sampling theory.58 The precision, expressed as the relative standard deviation of the mean zone distance and variance, improves in a similar manner. As the time increment is decreased from 1.0 × 10-3 to 1.0 × 10-5 s, the computational time increases inversely, and the accuracy (58) Meyer, S. L. Data Analysis for Scientists and Engineers; Wiley: New York, 1975.
1072 Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
Figure 8. Validation of the electrophoretic migration algorithm for multiple species (phosphate) at different values of pH by comparison of the mean zone distance and variance with theoretical predictions. Simulation conditions: N ) 500; t ) 1.0 × 10-4 s; R0 ) 2.5 × 10-3 cm; L ) 1.0 × 102 cm; V ) 2.0 × 104 V; v0 ) 0.0 cm s-1; D ) 1.0 × 10-5 cm2 s-1; phosphate species at pH 3.0 (O), 5.0 (0), 7.0 (4), and 9.0 (]). Equilibrium constants and electrophoretic mobilities of individual phosphate species are given in Table 2. Theory (s): Z h) (∑Riµi)VT/L; σ2 ) 2DT.
and precision are improved. Finally, as the total simulation time is increased from 1 to 100 s, the computational time increases linearly. The accuracy and precision of the mean zone distance improve with increasing simulation time, whereas the accuracy and precision of the variance become worse. These results provide guidance in the selection of computational parameters to achieve a satisfactory level of accuracy and precision within a reasonable computational time. CORRELATION OF EXPERIMENTAL AND SIMULATED SEPARATIONS In the studies described above, the algorithms for diffusion, electroosmotic flow, and electrophoretic migration have been validated by comparison to classical transport models. It is equally important, however, to establish the extent of correlation between the stochastic simulation and experimental results. For this purpose, the 5′-monophosphates of adenosine (AMP), guanosine (GMP), cytidine (CMP), and uridine (UMP) were chosen as simple, well-characterized solutes. The separation of these nucleotides by capillary zone electrophoresis has been experimentally studied and optimized previously in our laboratory.59 The optimum separation was achieved by using a phosphate buffer solution of pH 9.915 and ionic strength 12.5 mM in a fused-silica capillary of 75.5 µm inner diameter and 112.15 cm total length that is maintained under constant current conditions of 12.5 µA. To simulate this separation, it is necessary to provide the input parameters related to the properties of the nucleotides and the electrophoretic system. The equilibrium constants and electro(59) McGuffin, V. L.; Tavares, M. F. M. Anal. Chem. 1997, 69, 152-164.
Table 3. Performance of the Computer Simulation under Standardized Conditionsa as a Function of the Number of Molecules, Time Increment, and Total Simulation Time distance (cm) RSD (%)b error (%)c
no. of molecules
time increment (s)
total time (s)
mean
100 1000 10 000 1000 1000 1000 1000 1000 1000
1.0 × 10-4 1.0 × 10-4 1.0 × 10-4 1.0 × 10-3 1.0 × 10-4 1.0 × 10-5 1.0 × 10-4 1.0 × 10-4 1.0 × 10-4
10 10 10 10 10 10 1 10 100
2.973 2.974 2.974 2.974 2.974 2.973 0.2976 2.974 29.74
0.034 0.032 0.009 0.026 0.032 0.013 0.043 0.032 0.001
-0.046 -0.011 -0.005 -0.013 -0.011 -0.028 0.054 -0.011 -0.001
mean 3.18 × 10-4 3.23 × 10-4 3.38 × 10-4 3.47 × 10-4 3.23 × 10-4 3.41 × 10-4 3.25 × 10-5 3.23 × 10-4 3.25 × 10-3
variance (cm2) RSD (%)b 5.33 4.71 1.72 9.07 4.71 3.26 3.85 4.71 9.29
error (%)c
CPU timed (min)
-4.30 -2.73 1.81 4.44 -2.73 2.61 -2.06 -2.73 -2.12
4.67 46.4 461 4.78 46.4 465 4.74 46.4 466
a Simulation conditions: N ) variable; t ) variable; T ) variable; R ) 2.5 × 10-3 cm; L ) 1.0 × 102 cm; V ) 2.0 × 104 V; T ) 25 °C; ζ ) 2.553 0 0 × 102 mV; ) 78.5; η ) 9.0 × 10-1 cP; pH ) 7.0; I ) 1.0 × 10-6 M; D ) 1.0 × 10-5 cm2 s-1. Equilibrium constants and electrophoretic mobilities b c of individual phosphate species are given in Table 2. Relative standard deviation (%) ) (standard deviation) × 100/mean. Relative error (%) ) (mean - theory) × 100/theory. d Central processing unit time required for IBM RS/6000 model 580 computer.
Table 4. Electrophoretic Mobilities and Equilibrium Constants for the Individual Nucleotide Species at Infinite Dilution60 solute AMP
GMP
CMP
UMP
charge
mobility (cm2 V-1 s-1)
0
0.00 × 10-4
-1
-2.39 × 10-4
-2 0
-4.00 × 10-4 0.00 × 10-4
-1
-2.46 × 10-4
-2
-3.96 × 10-4
-3 0
-5.59 × 10-4 0.00 × 10-4
-1
-2.22 × 10-4
-2 0
-4.12 × 10-4 0.00 × 10-4
-1
-2.60 × 10-4
-2
-4.18 × 10-4
-3
-5.86 × 10-4
equilibrium constant 1.04 × 10-4 3.60 × 10-7 1.43 × 10-3 2.56 × 10-7 1.35 × 10-10 3.40 × 10-5 5.56 × 10-7 3.17 × 10-3 3.08 × 10-7 1.34 × 10-10
phoretic mobilities for the individual species of each nucleotide at infinite dilution are summarized in Table 4. These values were determined from experimental measurements of the effective mobility as a function of pH by using a nonlinear regression method.59,60 The parameters of the electrophoretic system have been measured or derived directly from the experimental system, whenever possible, to ensure that the simulation conditions are as representative as possible of the true experimental conditions. On the basis of these input parameters, the electrophoretic separation of the nucleotides was simulated under the optimal conditions specified above. The simulated evolution of the nucleotide separation is shown as a function of time and distance in Figure 9. In addition, the simulated separation is shown at a fixed distance of 43.4 cm (Figure 10, top), which corresponds to the position of the detector in the experimental system (Figure 10, bottom). The qualitative features of the simulated and experimental separations appear to (60) Tavares, M. F. M. Ph.D. Dissertation, Department of Chemistry, Michigan State University, East Lansing, MI, 1993.
Figure 9. Simulated evolution of the electrophoretic separation of nucleotides as a function of distance and time. The initial solute zone is shown at zero time. Simulation conditions: N ) 500; t ) 1.0 × 10-4 s; T ) 600 s; R0 ) 3.775 × 10-3 cm; L ) 1.1215 × 102 cm; Ldet ) 4.34 × 101 cm; ldet ) 5.0 × 10-1 cm; Linj ) 1.065 × 10-1 cm; linj ) 2.13 × 10-1 cm; V ) 2.533 × 104 V; T ) 25 °C; ζ ) 1.0613 × 102 mV; ) 78.5; η ) 9.0 × 10-1 cP; pH ) 9.915; I ) 1.25 × 10-2 M; D ) 1.0 × 10-5 cm2 s-1. Solutes: 1, AMP; 2, CMP; 3, GMP; 4, UMP. Equilibrium constants and electrophoretic mobilities of individual nucleotide species are given in Table 3.
be very similar; thus, the stochastic simulation has correctly predicted the elution order of the nucleotides and their resolution. A more quantitative comparison of the zone profiles is provided Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
1073
Table 5. Comparison of the Mean Migration Time and Variance for the Nucleotides from Experimental and Simulation Results variance (σ2)
migration time (s) solute AMP CMP GMP UMP
experiment
simulation
390.2 ( 4.1 400.6 ( 4.1 449.4 ( 4.2 477.1 ( 5.9
382.67 ( 0.04 391.63 ( 0.01 444.70 ( 0.04 470.50 ( 0.04
error
average a
(%)a
-1.93 -2.24 -1.05 -1.38 -1.65
experiment
simulation
error (%)a
2.34 ( 0.13 2.61 ( 0.14 3.04 ( 0.40 3.19 ( 0.30
2.52 ( 0.01 2.68 ( 0.06 3.58 ( 0.16 3.95 ( 0.05
7.69 2.68 17.76 23.82 12.99
Relative error (%) ) (simulation - experiment) × 100/experiment.
Figure 10. Simulated (top) and experimental (bottom) electrophoretic separations of the nucleotides. Simulation conditions are given in Figure 9. Experimental conditions: fused-silica capillary of 75.5 µm i.d. and 112.15 cm total length; phosphate buffer of pH 9.915 and ionic strength 1.25 × 10-2 M; constant current conditions of 12.5 µA with corresponding voltage 2.533 × 104 V; hydrodynamic injection at 2 cm height for 60 s; UV absorbance detection at 260 nm, 0.50 cm window at 43.4 cm length. Solutes: 1, AMP; 2, CMP; 3, GMP; 4, UMP, 1.0 × 10-4 M in phosphate buffer. The different molar absorptivities of the nucleotides lead to their unequal peak heights in the experimental electropherogram.
in Table 5, which summarizes the mean migration time and variance together with their precisions for each of the nucleotides. The simulated migration times appear to be in good agreement with the experimental values, having an average relative error of -1.65%. The simulated variances appear to be uniformly higher than the experimental values, having an average relative error of 12.99%. This error may arise from inaccuracy in the diffusion coefficient, which was not measured but was assumed to be 1.0 × 10-5 cm2 s-1 for all nucleotides. The error may also arise from inaccuracy in estimating the variance contributions from injection and detection. These contributions collectively comprise more than 75% of the total variance and are difficult to determine in the experimental system. In general, however, there is a high degree of qualitative and quantitative correlation between the stochastic simulation and the experimental results. 1074 Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
CONCLUSIONS The three-dimensional molecular simulation presented here provides several distinct advantages over previous stochastic simulations. First, the diffusion algorithm with variable step size, when coupled with small time increments (e.g., 10-4 s), yields discrete molecular steps that are very small and approach the appearance of continuous motion. This is necessary for statistically representative sampling of the electroosmotic flow velocity, which decreases rapidly in the vicinity of the wall. By incorporating the Rice-Whitehead flow profile in the electroosmotic flow algorithm, a potentially significant source of zone broadening is included that has been neglected in previous simulations. Finally, the electrophoretic migration algorithm considers the interconversion of solute species by acid-base and complexation equilibria and accounts directly for the equilibrium constant and mobility of each species. These algorithms for diffusion, electroosmotic flow, and electrophoretic migration have been validated by comparison with classical theory under a wide range of conditions representative of electrophoresis. Because of the modular nature of the simulation program, the algorithms may be simply and rapidly modified. The present algorithms may also be combined with those for laminar convection and surface interaction by partition and adsorption mechanisms, which were developed in previous studies for chromatographic separations.37,41 This combination of transport processes will enable the study of a wide variety of complex phenomena in electrophoretic separations. For example, any distortion or bias that may be introduced experimentally from either hydrodynamic or electrokinetic injection61,62 has been neglected in the present study. However, these injection techniques can be readily simulated by applying laminar flow or electrophoretic migration, respectively, to generate the initial distribution of molecules. In addition, the effects of multiple convection phenomena from laminar and electroosmotic flow can be examined and characterized with this simulation. In experimental electrophoresis systems, laminar flow can arise inadvertently if the fluid phase reservoirs are not maintained at precisely the same height and the same pressure.63 This will generally lead to detrimental broadening of the solute zone profiles and a corresponding loss of resolution. On the other hand, laminar flow has been applied intentionally in the opposite direction of the electroosmotic flow to enhance the separation.64,65 The surface interaction algorithms developed in the previous study can be used to study the effects (61) Grushka, E.; McCormick, R. M. J. Chromatogr. 1989, 471, 421-428. (62) Dose, E. V.; Guiochon, G. Anal. Chem. 1992, 64, 123-128. (63) Grushka, E. J. Chromatogr. 1991, 559, 81-93.
of adsorption at the fused-silica capillary wall36,66,67 or partition into a surface-bonded organic phase or polymer film.68-70 These surface interactions may have a deleterious effect on zone broadening and resolution or, if carefully controlled, may be used to advantage for electrochromatography separations. Finally, because the stochastic simulation follows the trajectories of individual molecules in their local environments, it can be used to study the effects of spatial and temporal transitions.41 Such transitions may arise from any of the simulation parameters (voltage, temperature, fluid phase pH or ionic strength, etc.), may be radial or axial within the electrophoretic system, and may be continuous or discontinuous in nature. There are a variety of electrophoretic methods in which boundaries and gradients in fluid composition are an inherent and essential feature of the separation, such as MBE, ITP, and IEF methods. In addition, (64) Gobie, W. A.; Ivory, C. F. J. Chromatogr. 1990, 516, 191-210. (65) Culbertson, C. T.; Jorgenson, J. W. Anal. Chem. 1994, 66, 955-962. (66) Green, J. S.; Jorgenson, J. W. J. Chromatogr. 1989, 478, 63-70. (67) Demana, T.; Chen, C. Y.; Morris, M. D. J. High Resolut. Chromatogr. 1990, 13, 587-589. (68) Cobb, K. A.; Dolnik, V.; Novotny, M. Anal. Chem. 1990, 62, 2478-2483. (69) Dougherty, A. M.; Woolley, C. L.; Williams, D. L.; Swaile, D. F.; Cole, R. O.; Sepaniak, M. J. J. Liq. Chromatogr. 1991, 14, 907-921. (70) Towns, J. K.; Regnier, F. E. Anal. Chem. 1992, 64, 2473-2478. (71) Mikkers, F. E. P.; Everaerts, F. M.; Verheggen, T. P. E. M. J. Chromatogr. 1979, 169, 1-10. (72) Mikkers, F. E. P.; Everaerts, F. M.; Verheggen, T. P. E. M. J. Chromatogr. 1979, 169, 11-20. (73) Chien, R. L.; Burgi, D. S. Anal. Chem. 1992, 64, 489A-496A.
the sample itself can serve to alter the local fluid composition in ZE, an effect that can be used advantageously to focus large sample volumes by stacking the adjacent zones.71-73 All of these situations can be addressed by means of a single, unified simulation program that incorporates the transport algorithms described above. Thus, the stochastic molecular simulation method can provide a powerful and versatile means to examine and characterize the complex behavior of solute molecules in electrophoretic separations. ACKNOWLEDGMENT The authors are grateful to Dr. Mark R. Schure (Rohm and Haas Co.) and Dr. Peiru Wu (Michigan State University) for insightful discussions. This research was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, under Contract No. DE-FG02-89ER14056. In addition, partial support for the IBM RS/6000 model 580 computer was provided by the Michigan State University Office of Computing and Technology and the Center for Fundamental Materials Research.
Received for review July 29, 1997. Accepted December 17, 1997. AC970802L
Analytical Chemistry, Vol. 70, No. 6, March 15, 1998
1075
