Dynamic computer simulations of electrophoresis: 2010-2020.

The transport of components in liquid media under the influence of an applied electric field can be described with the continuity equation. It represents a nonlinear conservation law that is based upon the balance laws of continuous transport processes and can be solved in time and space numerically. This procedure is referred to as dynamic computer simulation. Since its inception four decades ago, the state of dynamic computer simulation software and its use has progressed significantly. Dynamic models are the most versatile tools to explore the fundamentals of electrokinetic separations and provide insights into the behavior of buffer systems and sample components of all electrophoretic separation methods, including moving boundary electrophoresis, CZE, CGE, ITP, IEF, EKC, ACE, and CEC. This article is a continuation of previous reviews (Electrophoresis 2009, 30, S16-S26 and Electrophoresis 2010, 31, 726-754) and summarizes the progress and achievements made during the 2010 to 2020 time period in which some of the existing dynamic simulators were extended and new simulation packages were developed. This review presents the basics and extensions of the three most used one-dimensional simulators, provides a survey of new one-dimensional simulators, outlines an overview of multi-dimensional models, and mentions models that were briefly reported in the literature. A comprehensive discussion of simulation applications and achievements of the 2010 to 2020 time period is also included.

conservation element and solution element

electromigration dispersion

R‐phycoerythrin

symmetric limited positive

Introduction

The development of simulation models for electrophoresis has been underway for more than 40 years. Shortly after computers became available, scientists at universities in the Czech Republic (Charles University in Prague, work of B. Gaš [1]), Switzerland (University of Bern in Bern, work of P. Ryser [2]), and the United States (University of Arizona in Tucson, work of G.T. Moore [3]) began to construct dynamic computer models for electrophoresis with the goal of exploring the basics of electrokinetic separations. These early models were restricted to strong electrolytes and can be regarded as the first dynamic electrophoretic simulation models. The first dynamic models predicting the behavior of strong and weak electrolyte systems were developed in the 1980s by Bier et al. [4, 5, 6], Radi and Schumacher [7], Roberts [8], and Schafer‐Nielsen [9]. Alternatively, computer models for the prediction of steady‐state distributions in IEF [10, 11, 12, 13] and ITP [14, 15] emerged in the same time period. Over the years, many dynamic simulation models of various degrees of complexity have been described in the literature. The history of dynamic computer simulation in electrophoresis together with a survey of the models published up to 2009 was given in a previous review [16]. The fundamentals of computer simulations of electrophoretic processes, including comprehensive simulation examples and applications, were the subject of a book [17], a review paper [18], and two book chapters [19, 20]. Furthermore, the use of simulations for the prediction of separations in IEF [21], instabilities in IEF [22], microfluidic ITP [23], and chiral separations [24] was discussed in the literature. Dynamic simulations were also mentioned in the context of theoretical principles of capillary electromigration methods [25].

Simulation efforts were and still are driven by a desire to understand the fundamentals of all electrophoretic separation methods, including moving boundary electrophoresis, CZE, CGE, ITP, IEF, EKC, ACE, and CEC. Underlying chemical and physical processes that are involved in the separation of compounds using these techniques can thereby be identified. The transport of components in liquid media under the influence of an applied electric field can be described with the continuity equation that represents a nonlinear conservation law that is based upon the balance laws of continuous transport processes [4, 5, 6, 17, 26, 27]. There is no simple solution to the continuity equation. However, it can be solved in time and space numerically that is referred to as dynamic computer simulation. Versatile dynamic computer models, such as GENTRANS [4, 5, 6, 28, 29, 30, 31], SIMUL5 [27], and SPRESSO [32, 33], are applicable to any initial buffer and sample configuration and follow the impact of the electric field on all components at every step from start to finish according to the underlying transport laws. They are able to predict the evolution of electrophoretic boundaries, sample components, system peaks, and all zones formed in electrolyte systems together with associated column properties such as pH, conductivity, electric field strength, ionic strength, and fluid flow.

In addition to the dynamic models that provide complete sets of concentration, conductivity, pH, and flow profiles as a function of time and location, there are also electrophoretic computer models in use that are based on the linear theory of electrophoresis. These models enable a rapid assessment of buffer systems and analyte separability in CZE and EKC and predict parameters of background electrolytes, system peaks, and analyte peaks (PeakMaster 5.4 [34, 35, 36, 37] and PeakMaster 6 [38], both available through the internet at https://echmet.natur.cuni.cz).

In the 2010 to 2020 time period, some of the existing dynamic simulators were extended and new simulation packages were developed. This review presents the basics and extensions of the three most used simulators GENTRANS, SIMUL5, and SPRESSO, provides a survey of new one‐dimensional simulators, gives an overview of multi‐dimensional models, and mentions models that were briefly reported in the literature. A comprehensive discussion of simulation applications is also included. A number of papers published in 2009 that were not mentioned in the review of 2010 [18] are also discussed.

Dynamic simulators for electrokinetic separations

Dynamic simulators are based upon algebraic acid‐base and/or complexation equations and continuity equations, which are partial differential equations in time and space that can only be solved numerically using computers. Such models calculate the transport of each component through the electrophoretic space as a result of electromigration, diffusion, imposed, and/or electrically driven bulk flow, solution‐based chemical reactions such as protolysis and, if incorporated, also interaction of solutes with electrolyte additives, the column walls, and the column matrix.

The inclusion of the diffusion current, i.e., the current carried by diffusion, has to be part of the transport equations of a dynamic electrophoretic simulator. This is illustrated with the simple cationic ITP example presented in Fig. 1. It represents a system with 10 mM potassium acetate as leader, 10 mM lithium acetate as terminator, and 10 mM sodium acetate as sample [20, 39, 40, 41]. Acetic acid acts as counter component. The simulation was performed with GENTRANS with a 1 cm column divided into 4000 segments (Δx = 2.5 μm) and a constant 2000 A/m2. The sample was applied between 5 and 10% of column length, the concentrations at the column ends were kept constant and no data smoothing was used. Upon application of power, the sodium zone becomes adjusted to a concentration of 8.48 mM as it gradually penetrates into the space originally occupied by the leader where it migrates between potassium and lithium. The same applies to lithium, the terminating constituent, which becomes adjusted to a concentration of 7.59 mM. Thereafter, the entire ITP zone structure is migrating at a constant velocity and without change of the zone structure toward the cathode. The data presented in Fig. 1 represent computer‐predicted profiles after 0.08 min (4.8 s) of the current application. These data illustrate that sharp boundary transitions with thicknesses of about 20 μm are formed that are not associated with numerical oscillations (Fig. 1A and D). If the current carried by diffusion is included in the calculations (solid black lines in Fig. 1), the pH change between the zones is predicted to be sigmoidal (Fig. 1C). Without the diffusion current in the transport equations, nonphysical spikes are predicted in the zone boundaries of the pH profiles [39, 40] that is illustrated with the red dotted line profiles in Fig. 1C. In all other panels, the overlaid profiles of data obtained with (solid black lines) and without (dotted red lines) the diffusion current also reveal small differences in the boundary shapes, as is best illustrated with the conductivity (insert in Fig. 1B) and concentration profiles (Fig. 1D) depicted for the transition between potassium and sodium. This system can be employed as a test to investigate whether the diffusion current is correctly included [41] and is recommended to developers of new models.

Figure 1

Example showing the importance of having the diffusion current included in the code. Simulated (A) concentration, (B) conductivity, and (C) pH distributions for cationic ITP of Na+ (sample) between K+ (leading component) and Li+ (terminating component). The concentration changes across the boundary between K+ and Na+ are presented in panel (D). In the bottom panel, data of the counter component (acetic acid) are not shown. The dotted line graphs in all panels depict the profiles obtained without inclusion of the diffusion current in the transport equations. The insert in the conductivity data graph depicts the transition between K+ and Na+. The data were produced with GENTRANS without using the smoothing option. Simulation conditions: 1 cm column divided into 4000 segments (Δx = 2.5 μm), 2000 A/m2 applied for 0.08 min. The mobilities of potassium, sodium, and lithium were taken as 7.91 × 10–8 m2/Vs, 5.19 × 10–8 m2/Vs, and 4.10 × 10–8 m2/Vs, respectively. The pKa and mobility of acetic acid were 4.76 and 4.12 × 10–8 m2/Vs, respectively. The cathode is to the right.

The execution of a dynamic simulation requires the specification of the initial distributions of all compounds, the physico‐chemical input data of the compounds in terms of mobilities and chemical equilibrium constants, the column geometry and its mesh, the input data for buffer flow or calculation of EOF, the applied constant voltage or constant current, and the time of power application together with the data storage interval. In performing a simulation, there are some points that the researcher must understand in order to maximize the utility of their simulation. Selection of the simulation parameters, such as the number of segments and the applied power level, is crucial for success. These aspects are well described in various publications [18, 20, 42] and are thus not discussed again in this review.

One‐dimensional models

A survey of frequently used and new one‐dimensional dynamic simulators is presented in Table 1. GENTRANS [4, 5, 6, 28, 29, 30, 31], SIMUL5 [27], and SPRESSO [32, 33] are by far the three most employed electrophoretic simulators. Although they differ in a number of specifications, they were shown to provide identical results when used with the same input conditions [42]. This is illustrated with the simulation data presented in Fig. 2. To compare the output of the three simulators, a uniform grid was employed to predict the shape of an ITP boundary at a high current density that is a demanding task. The data presented in Fig. 2 represent the anionically migrating ITP boundary between chloride and HEPES with Tris as a counter component, an example that was discussed in refs. [32, 42]. The anolyte comprised 100 mM hydrochloric acid and 200 mM Tris, whereas the catholyte was composed of 50 mM HEPES and 200 mM Tris. Simulations were performed in a 2 cm column divided into 16 000 segments of equal length (Δx = 1.25 μm) and via application of a constant current density of –5000 A/m2 for 6 s. GENTRANS was executed without using data smoothing and SPRESSO was run with the high‐resolution sixth‐order compact scheme. The initial boundary between anolyte and catholyte was at 20% of column length. The data presented reveal that the concentration profiles of the three components across this migrating steady‐state boundary of about 20 μm width are equal.

Table 1

Survey of selected dynamic one‐dimensional computer models for simulations of electrophoretic processes

Model	Numerical scheme	Features up to 2010	Features added 2010–2020 and new simulators
GENTRANS 1)	Second order centered (finite difference)	Basic version for monovalent components and simple ampholytes; options: various boundary conditions, data smoothing, detector profiles [4, 5, 6] Addition of proteins [28, 29, 43], imposed plug flow [44, 45], in situ calculated EOF [30, 31, 46, 47, 48], swapping of electrolyte in part of the column prior to continuation of a run [49, 50], flux corrected transport (flux limiter) for comparison of flux corrected transport with upwind and second order centered scheme [40, 51], high‐resolution version with addition of multivalent components and simple input of set of ampholytes for IEF [50, 52, 53, 54], micellar EKC interactions [55, 56]	Addition of Taylor‐Aris dispersion [57] Addition of complexation equilibria with monovalent components [58] Addition of chiral stationary phase in presence of EOF [59]
SIMUL5 2) SIMUL5complex 2)	Second order centered (finite difference)	SIMUL5: For monovalent and multivalent components, incl. ampholytes, and EOF as constant plug flow; data base of components; options: correction of mobilities for ionic strength and calculations of activity coefficients, input of set of ampholytes for IEF [27] Micellar EKC version 1,3) and special IEF version [60]	SIMUL5complex: Addition of fast complexation equilibria for monovalent and polyvalent components [61, 62]
SPRESSO 2)	6th order compact scheme combined with adaptive grid	Monovalent and multivalent components, incl. ampholytes; model for pressure driven flow with Taylor‐Aris dispersion; electroosmosis; data base of components; options: four numerical schemes and moving frame [32, 33]	Correction of mobilities for ionic strength [66] Quasi 1D model for variable cross‐section channels featuring a finite volume method based on a SLIP 6) scheme together with an adaptive grid algorithm [67] Coupling of nonlinear electromigration with multispecies nonequilibrium kinetic reactions in bulk solution and surfaces [68]
CESE models 1)	Space‐time CESE 4) CFL‐insensitive CESE 5)	Monovalent components [69] 1D reduced model for microchannels with uniform or variable cross‐sectional area, monovalent components and simple ampholytes [70]	Combination of adaptive grid redistribution and CESE scheme for ITP and IEF, mono‐ and multivalent components [71]
SIMUL6 2)	Second order centered (finite difference)		Model with multithreaded computation for monovalent and multivalent components, incl. ampholytes, data base of components, in situ detector profiles; options: variable cross‐section, electrolyte swapping, input of set of ampholytes for IEF, specific mobilities and pKa values for specified column segments [72]
SPYCE 2)	Fourier pseudo‐spectral method		The numerical scheme allows for fast high‐resolution numerical electrophoresis simulations on a periodic domain and is thus restricted to situations with equal compositions at anodic and cathodic column ends [75].

1)Available in the laboratories of the program developers; 2)available for free via internet; 3)B. Gaš, personnel communication; 4)CESE: conservation element and solution element; 5)CFL = v Δt/Δx (product of local velocity and time interval/grid size ratio); 6)SLIP: symmetric limited positive

Figure 2

Equality of predicted ITP boundary shape by three simulators. The data depict the boundary structure between chloride and HEPES with Tris as counter component simulated with GENTRANS, SIMUL5, and SPRESSO. The anolyte comprised 100 mM hydrochloric acid and 200 mM Tris, whereas the catholyte was composed of 50 mM HEPES and 200 mM Tris. Simulations were performed in a 2 cm column divided into 16 000 segments of equal length (Δx = 1.25 μm), at a constant current density of –5000 A/m2, and for 6 s of electrophoresis. The initial boundary was at 20% of column length. The pKa and mobility value used for the simulation were 7.50 and 2.35 × 10 8 m2/Vs, respectively, for HEPES, and 8.05 and 2.95 × 10 8 m2/Vs for Tris. The mobility of chloride was 7.90 × 10 8 m2/Vs. The anode is to the right.

GENTRANS is the oldest of the three models [4, 5, 6, 28]. It was adapted to be executed on PCs, extended with new features [29, 30, 31, 43, 44, 45, 46, 47, 48, 49, 50, 51], made applicable for high‐resolution simulations [50, 52, 53, 54], and extended for the simulation of micellar EKC [55, 56]. Current versions run under Windows (XP, 7 [64‐bit; 32‐bit until 2018 only], and 10) and are available from the corresponding author upon request. GENTRANS does not feature a shell with a user‐friendly possibility to input or change parameters and to graphically display the progress of an ongoing simulation. GENTRANS comprises seven separate modules, namely those for simple biprotic ampholytes, monovalent weak acids, monovalent weak bases, monovalent strong acids, monovalent strong bases, peptides (used for peptides and other multivalent components), and proteins. This is an attractive approach because the operation of the monovalent codes is faster compared to that for multivalent components that reduces computation time for monovalent components and simple biprotic ampholytes and thereby proved to be beneficial for the simulation of IEF in presence of a large number of carrier ampholytes [50, 52, 53]. Simulations with GENTRANS can be executed using several options for the column boundaries: (1) with open column ends which allow mass transport into and out of the separation space; (2) with fixed concentrations at column ends; (3) with column ends that are impermeable to any buffer and sample compounds; and (4) with mixed conditions with no transport or free transport of each component at both left and right boundaries. GENTRANS offers the possibility of using data smoothing (removal of negative concentrations), which can speed up a simulation and avoid numerical oscillations [42]. In GENTRANS, electrophoretic mobilities are considered to be independent of the ionic strength and temperature. Furthermore, GENTRANS features algorithms that estimate the magnitude and impact of electroosmosis in CE via the use of inputs that are dependent on the column wall material and the ionic strength of the electrolyte [30, 31, 46, 47, 48] and allows swapping the content in a part of the column with a new electrolyte prior to continuation of a run [49, 50].

Extensions of GENTRANS that were implemented during the 2010–2020 time period include Taylor–Aris dispersion to account for dispersion due to the parabolic flow profile associated with pressure‐driven flow [57] and chemical interactions between components required for the prediction of chiral separations [58]. With the former feature, effective diffusivity of analyte and system zones as functions of the capillary diameter and the amount of flow in comparison to molecular diffusion alone can be studied for configurations with concomitant action of imposed hydrodynamic flow and electroosmosis. Data obtained under realistic experimental conditions, for example, a 50 μm id fused‐silica capillary of 90 cm total length, revealed that inclusion of flow profile‐based Taylor‐Aris diffusivity provides simulation patterns of analyte and system peaks that compare well with those monitored experimentally with UV and conductivity detection [57].

The GENTRANS EKC code used to simulate separations, transient trapping, and sweeping in micellar EKC [54, 55] was extended by Breadmore et al. for chiral separations via consideration of complexation constants and specific mobilities of 1:1 analyte‐selector complexes [58]. The model handles interactions between monovalent weak and strong acids and bases with a single monovalent weak or strong acid or base additive, including a neutral CD, under real experimental conditions. It is a tool to investigate the dynamics of chiral separations and to provide insight into the buffer systems used in chiral EKC and ITP together with features of analyte stacking and destacking. Furthermore, an option that permits the simulation of chiral CEC in which the selector zone remains immobile even under conditions of EOF was incorporated. Together with complex mobilities set to zero to provide an electrophoretically immobilized selector, this approach can be employed to study the CEC behavior of analytes within the stationary phase and at the transitions between the stationary phase and free solution [59].

SIMUL5 of Hruška et al. appeared in 2006 [27]. It is based on one equation that handles protolysis of all components [26], comprised a newly designed model, and replaced previous SIMUL models that were of limited scope (see ref. [16] for overview). SIMUL5 runs under WINDOWS and can be downloaded from https://echmet.natur.cuni.cz. It features the use of a reduced calculation space with moving borders that bracket the separation space with changing concentrations (part of the column where considerable changes from the initial values are expected), an approach that leads to faster simulations. SIMUL5 has a comfortable windows environment for data input, data evaluation, visual control of the ongoing simulation, and visualization of a completed simulation in a movie format. SIMUL5 offers the option of considering mobilities as a function of ionic strength and permits the implementation of a constant buffer flow with a flat flow profile to mimic EOF. It features boundary conditions for open column ends, which allow mass transport into and out of the separation space. Non‐released versions of SIMUL5 include closed boundary conditions at column ends in a column divided into an array of compartments with identical or different cross sections in which the solutions are mixed and permitted simulation of isoelectric trapping separations and desalting that take place in recirculating multicompartmental electrolyzers [60], and the possibility of simulating micellar EKC systems (Table 1).

During the 2010 to 2020 time period, SIMUL5 was extended with algorithms that describe 1:1 chemical equilibria between solutes and a buffer additive with fast interactions that can be considered instantaneous in comparison to the time scale of peak movement [61, 62]. This version of SIMUL5, referred to as SIMUL5complex, handles equilibria with neutral, singly charged, and multiply charged complexation reagents, is more versatile than the EKC possibility of GENTRANS, and was extensively applied to describe chiral interactions, enantiomer separations, and electromigration dispersion effects caused by complexation [24, 36, 37, 63, 64, 65].

SPRESSO is a MATLAB based, fast open‐source nonlinear electrophoresis solver that appeared in 2009 and is available for free at http://microfluidics.stanford.edu/spresso [32, 33]. It is based on the same general component equation as SIMUL5, features an adaptive grid system to reduce the computational time interval in configurations with a few components but not in presence of a large number of constituents that are unevenly distributed throughout the column, and has visual control of the ongoing simulation. Thus, the use of SPRESSO with its adaptive grid can be beneficial for rapid predictions of configurations with sharp sample or buffer boundaries that are confined to a small section of the electrophoretic column. SPRESSO is written in the secondary programming language of MATLAB (MathWorks). The compiled version of SPRESSO can be run on Windows‐based PCs via installation of the MATLAB Component Runtime (MCR) library. Data evaluation, however, requires MATLAB.

In the first version of SPRESSO, selections for a simulation included four different numerical schemes for spatial discretization, e.g., a high‐resolution sixth‐order compact scheme, a model for pressure‐driven flow and Taylor‐Aris dispersion, an approach to estimate electroosmosis, as well as a moving frame [32, 33]. The second version featured the option of considering mobilities as a function of ionic strength and permitted the input of a convergence tolerance value [66]. Thereafter, SPRESSO was extended to simulate electrokinetic processes in channels with nonuniform cross‐sectional areas and was thereby converted into a quasi 1D model [67]. In this third version, the quasi 1D governing equations are solved with a dissipative finite volume method based on a symmetric limited positive (SLIP) scheme together with an adaptive grid algorithm. This approach is claimed to ensure both unconditional stability and high accuracy and allows to perform fast simulations at high electric field strengths with a small number of grid points. Subsequently, the nonlinear electromigration part of SPRESSO was coupled to multispecies nonequilibrium kinetic reactions [68]. This 1D simulation tool is applicable to both homogeneous reactions and surface‐based reactions. It provides an aid to assist the development, analysis, and optimization of electrophoresis‐based biosensors and assays involving nonequilibrium chemical reactions, such as immunoassays.

The space–time conservation element and solution element (CESE) model of Yu et al. is claimed to be more accurate than conventional numerical schemes, suppresses the numerical oscillations or peaks observed in the results obtained using traditional second‐order finite difference schemes, and was applied to high‐resolution simulations of CZE and ITP [69]. This simulator was first extended to handle columns with uniform or variable cross‐sectional areas (1D reduced model for microchannels) and applied to high‐resolution simulation of IEF [70]. The last development encompasses the combination of an adaptive mesh redistribution with CESE that was successfully applied to high‐resolution ITP and IEF simulations [71]. No further developments or applications of the CESE method were found in the literature.

SIMUL6, the successor of SIMUL5, was recently released as free software that can be downloaded from https://echmet.natur.cuni.cz and https://www.simul6.app. A completely new source code was written. It includes faster procedures for the numerical integration of partial differential equations and uses multithreaded computation. This made the simulator up to 15 times faster compared to SIMUL5 [72]. SIMUL6 runs on Windows‐based 64‐bit computers, Linux 64‐bit, and macOS 64‐bit, and it has a user‐friendly interface to input all parameters. Simulations in progress can be followed graphically with three in situ data frames for (i) the distributions of all components, pH, conductivity, and the electric field strength along the column, (ii) the current and voltage as a function of time, and (iii) the signal of a detector placed at a selected position along the column. The detector graph encompasses signals for components, conductivity, and pH. SIMUL6 has a number of very useful features. The separation column can be divided into a number of individual sections and the diameter of each section can be set individually. In all segments, transversal diffusion is assumed to be infinitely fast. Thist represents a quasi 1D approach that was used previously in other models [60, 67, 70]. Having larger diameters at the column ends permits the simulation of the impact of electrode compartments. Furthermore, mobilities and pKa values of the components can be specified for different column segments that, e.g., allows the study of mass transport across a boundary between free solution and a zone containing a neutral gel in which migration is retarded [73]. Swapping of a part of the electrolyte with another electrolyte after some time during the run is a procedure that is useful to study electrophoretic mobilization in IEF [50] or to characterize ITP condensation followed by CZE separation [49]. These optional features are of interest in various situations of chip electrophoresis. The current version of SIMUL6 (version 6.1.2) neither includes electroosmotic and/or hydrodynamic flow along the column nor complexation reactions.

SIMUL6 comprises a convenient input routine for a large number of carrier ampholytes that produce the pH gradient in IEF. As an example, focusing data of 14 amphoteric analytes together with 182 carrier ampholytes and two spacing components between 300 mM H3PO4 (anolyte, between 0 and 9 mm of column length) and 200 mM LiOH (catholyte, between 63 and 72 mm of column length) are presented in Fig. 3. The input data of all components used are those provided as example in SIMUL6 [72] and stemmed from ref. [74]. Focusing occurred during 500 s at 600 V in a 50 μm id column of 72 mm length that was divided into 6000 segments of equal length (12 μm mesh). Thereafter, the cathodic electrode solution was replaced with the anolyte and electrophoretic (chemical) mobilization was induced by a continuation of power application (600 V) for 700 s. These data reveal that (i) the 14 markers are separated according to their pI and (ii) electrophoretic mobilization can change the detection sequence of amphoteric analytes, i.e., the analytes do not necessarily pass the point of detection in the order of decreasing pI values when mobilization is induced by exchanging the catholyte with the anolyte.

Figure 3

IEF pattern simulated by SIMUL6 with 14 analytes in a pH 3.0–10.3 gradient formed by 182 carrier ampholytes. Fourteen amphoteric analytes (0.01 mM each), 182 carrier ampholytes (0.1 mM each), iminodiacetic acid (IDA, 10.5 mM), and arginine (ARG, 8.79 mM) were supplied as homogeneous mixture between 300 mM H3PO4 (anolyte, between 0 and 9 mm of column length) and 200 mM LiOH (catholyte, between 63 and 72 mm of column length). A 50 μm id column of 72 mm length that was divided into 6000 segments of equal length (12 μm mesh) was employed. A constant voltage of 600 V was applied for 500 s. For all carrier ampholytes, ΔpKa was 1 and the pKa step was 0.0406 for each following ampholyte. The cationic and anionic mobilities of all ampholytes were set to 20 × 10–9 m2/Vs. The input data of the analytes used are those provided in ref. [74]. The data presented include the distributions of (A) all components, (B) pH and electric field strength, (C) the foci of the 14 analytes between IDA and ARG, and (D) the electrolytes and spacing components. Electrophoretic mobilization toward the cathode was induced by replacing the cathodic electrode solution with the anolyte followed by a continuation of power application at 600 V for 700 s. The detector plot for the 14 analytes and a detector placed at 57.6 mm is presented as insert in panel D. The cathode is to the right. Key: 1, serotonin (pI 10.58); 2, tyramine (10.17); 3, metanephrine (9.72); 4, epinephrine (9.32); 5, norepinephrine (9.21); 6, labetalol (8.49); 7, 3‐methylhistidine (7.71); 8, glycyl‐histidine (7.55); 9, leucoberbelin blue I dye (5.27); 10, 4‐(4‐aminophenyl) butyric acid (4.86); 11, dansylated γ‐aminobutyric acid (4.20); 12, dansylated glutamic acid (3.49); 13, dansylated aspartic acid (3.34); 14, dansylated iminodiacetic acid (2.99).

The SPYCE (Pseudo‐spectral Python Code for Electrophoresis) simulator is a new model that uses the Fourier pseudo‐spectral method for fast high‐resolution numerical simulations of a periodic situation in which the value of any spatially varying physical quantity at the beginning of the column is equal to that at the end [75]. It can thus be employed for the study of CZE, transient ITP within a sample zone, field‐amplified sample stacking, and oscillating systems, but not for ITP, moving boundary electrophoresis and IEF. The Fourier pseudo‐spectral method yields accurate and stable solutions on coarser computational grids compared with other nondissipative spatial discretization schemes. SPYCE can be obtained for free as open‐source code at http://web.iitd.ac.in/~bahga/SPYCE.html. Furthermore, a dynamic mathematical simulator of CE of unreported sophistication was developed in Russia by Sursyakova et al. and applied to the study of system peak formation and physical properties of migrating analytes [76, 77]. Bahga et al. described a diffusion‐free ITP model for channels with axially varying cross‐sections that predicts properties of homogeneous zones but not those of electrophoretic boundaries [78]. Finally, Zhang et al. developed a stump‐like mathematical model for computer simulation of CZE [79].

Multi‐dimensional models

A survey of frequently used and new multi‐dimensional dynamic simulators is presented in Table 2. The finite volume‐based 2D model of Shim et al. was developed from scratch at the Washington State University [80, 81, 82, 83]. It represents a solver that predicts IEF separations in closed microchannels (absence of electrode solutions, for details about these boundary conditions, see also ref. [22]) and was, e.g., used to study the behavior of protein bands in a horseshoe microchannel [84]. These simulations require an enormous number of calculations and, with a large number of components, can thus last up to a couple of months of calculation time. The 2D simulator was extended for faster IEF simulations with the parallel use of multiple CPUs [85], with an efficient algorithm of a parallel scheme [86], and was used for the validation of a steady‐state protein focusing model [87]. Other developments included the integration of Joule heating to study its effect on the IEF of proteins in a microchannel [88] and a mathematical and numerical model to study 2D free flow IEF [89]. With a proper change of the boundary conditions at the column ends, the same approach could be employed for studying the ITP process in straight channels [90] and in those featuring changes of the cross‐sectional area [91]. Furthermore, the 2D model was applied to determine a quasi‐steady state in complex microgeometries operated under constant voltage [92] and was extended to include temperature effects in the constant voltage mode of ITP in a microchannel [93].

Table 2

Survey of multi‐dimensional computer models for simulations of electrophoretic processes

Model	Numerical scheme	Features up to 2010	Features added 2010–2020 and new simulators
2D model	Finite‐volume	Model for monovalent components, simple ampholytes and model proteins [80, 81, 82] Mobility as function of ionic strength [83] IEF protein bands in horseshoe microchannel [84] Model for ITP in straight microchannels [90]	Parallel implementation with multiple CPUs for IEF simulations [85] and parallel scheme efficient algorithm for 2D IEF [86] Effect of Joule heating on IEF [88] 2D model for free flow IEF [89] 2D model of ITP in channels of changing cross‐sectional area [91]; quasi‐steady state of ITP in complex microgeometries [92] and effect of Joule heating in ITP [93]
CFD‐ACE+ with user models	Finite‐element	Effects of electrode configuration and setting on electrokinetic analyte injection in CZE [94, 95]	Analyte behavior in a curved channel [96] Electrokinetic supercharging [97], behavior of DNA in CGE after electrokinetic injection [98], conditions at the capillary tip during field‐amplified sample injection [99] pH gradient formation and cathodic drift in microchip [100]
COMSOL Multiphysics with electrophoresis interface	Finite‐element		1D simulations of electrophoretic separations [41]
COMSOL Multiphysics with user implemented models	Finite‐element	ITP of proteins in a networked microfluidic chip [101]	Models for studying various aspects of ITP and other CE modes in microchannels [102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113] Nanochannel electrophoresis model allowing deviation from electroneutrality [115]
PETSc‐FEM	Finite‐element	3D simulation model for electrokinetic flow and transport in microfluidic chips [116]	3D model for high‐performance simulation of electrophoretic techniques in microfluidic chips [117]
OpenFOAM with user implemented models	Finite‐volume		Open‐source toolbox for 3D electrophoretic separations [118] Model of electrokinetic transport in microfluidic paper‐ based systems [119, 120]

Besides the dedicated 2D simulator for electrophoresis that was developed at the Washington State University, multiphysics simulation packages, such as CFD‐ACE+, COMSOL Multiphysics, and OpenFOAM, were employed for the simulation of a variety of nonlinear electrophoresis techniques (Table 2). They enable both one‐ and multi‐dimensional high‐resolution dynamic simulations of electrophoresis via the implementation of individual models into the framework of the multiphysics packages. Furthermore, in order to increase the accessibility of multidimensional electrophoresis simulations, COMSOL Multiphysics was extended by an electrophoretic transport interface that was first released in COMSOL Multiphysics version 5.3. The performance of the interface in COMSOL Multiphysics was investigated by comparing 1D results with those generated with GENTRANS and SIMUL5. The profiles were found to be essentially identical that confirms that the algorithms incorporated into COMSOL Multiphysics are able to correctly describe the dynamics of electrophoretic processes [41]. The assembly and solver routines used by COMSOL Multiphysics are multithreaded. Using fast multicore PCs, the time intervals required for the simulation of the investigated examples were found to be comparable to those with GENTRANS. The COMSOL Multiphysics electrophoretic transport interface can be expected to be useful as a general model for the investigation of electrophoretic phenomena in 1D, 2D, and 3D geometries.

The CFD‐ACE+ software (version 2006, CFDRC, Huntsville, AL, USA) with a 2D finite element based electrophoretic model was previously used to characterize the effects of electrode configuration and setting on electrokinetic analyte injection in high‐resolution CZE [94, 95]. During the 2010 to 2020 period, this model was applied to explore the migration behavior of the analytes with or without ITP stacking in a curved channel [96], to study electrokinetic supercharging with a system‐induced terminator and an optimized capillary versus electrode configuration for parts‐per‐trillion detection of rare‐earth elements in CZE [97], to investigate DNA aggregation and cleavage in CGE induced by a high electric field during electrokinetic sample injection [98], and to describe the conditions at the capillary tip during field‐amplified sample injection with a mobility‐boost effect [99]. Furthermore, the same simulator was employed to investigate the pH gradient formation between anolyte and catholyte reservoirs together with its cathodic drift in microchip IEF with imaged UV detection [100].

COMSOL Multiphysics is a commercial multiphysics finite‐element simulator from Sweden (COMSOL AB, Stockholm, Sweden). This package was employed by various researchers working in the field of electrophoretic separations who applied their own models. At the Washington State University, COMSOL Multiphysics was employed to investigate various aspects of ITP in microchannels, including the behavior of proteins in comparison to CZE after they pass through a T‐junction (2D simulations, [101]), a 10 000‐fold concentration increase of proteins in a cascade microchip (3D, [102]), separation of lanthanides in presence of hydroxyisobutyric acid and acetate as complexing agents (1D, [103]), separation of eight lanthanides on a serpentine PMMA microchip (2D, [104]), and separation and concentration of fluorescently tagged cardiac troponin I from two proteins with similar isoelectric properties in a PMMA straight‐channel microfluidic chip (1D, [105]). Another effort involving COMSOL Multiphysics simulations led to the assessment of the effect of counterflow in ITP performed in an open capillary (2D, [106]) and in a monolithic column [107]. In the latter work, the fluid‐flow profile in a porous medium was simulated using the Brinkman Equation built into COMSOL Multiphysics. Counterflow ITP in the monolithic column showed undistorted analyte zones with significantly reduced dispersion compared to the severe dispersion observed in an open capillary. The data presented in Fig. 4 represent those obtained for free solution anionic ITP of fluorescein with high molecular diffusivity (4.25 × 10–10 m2/s) and R‐phycoerythrin (RPE) with a low molecular diffusivity (0.157 × 10–10 m2/s) in the absence of counterflow flow (Fig. 4A) and in presence of a counterflow that compensates anionic migration such that analyte zones become immobile (Fig. 4B and C). Figure 4D depicts a schematic illustration of the electromigration of an analyte counteracted by a parabolic counterflow. The electromigration velocity of the analyte is uniform whereas the flow induced by pressure or gravity is of parabolic shape. The composite migration velocity of the analyte is the superposition of these two velocities. The leading electrolyte is to the left, and the terminating electrolyte to the right. The applied current density is 226.4 A/m2 and the average flow velocity is 1.45 × 10−4 m/s. In absence of an applied counterflow, fluorescein and RPE are migrating as narrow peaks with fluorescein ahead of RPE (Fig. 4A). With a parabolic counterflow profile and molecular diffusivity (images (i) and solid line graphs in panels B and C) somewhat broader and distorted zones are predicted with the flow effect being significantly larger for the case of RPE with the much lower diffusivity. With a plug counterflow profile and Taylor‐Aris effective diffusivity (images (ii) and broken line graphs in panels B and C), peak distortion for both cases is predicted to be larger compared to the predictions for the parabolic counterflow. This indicates that Taylor‐Aris diffusivity is somewhat overestimating peak broadening [106].

Figure 4

2D simulation of ITP in the absence and presence of counterflow using COMSOL Multiphysics. (A) Anionic ITP of fluorescein (analyte with a high diffusivity) and R‐phycoerythrin (RPE, analyte with a low diffusivity) without application of counterflow. (B) Anionic ITP of fluorescein under stationary counterflow with (i) a parabolic counterflow profile and molecular diffusivity and (ii) a flat counterflow profile and Taylor‐Aris effective diffusivity. (C) same as (B) but for RPE. (D) Schematic illustration of electromigration of an analyte counteracted by a parabolic counterflow. The leading electrolyte was composed of 10 mM HCl and ethanolamine (pH 9.5) and the terminating electrolyte comprised 20 mM barium hydroxide. A 75 μm id coated fused‐silica capillary with minimized EOF served as separation space. The applied current density was 226.4 A/m2 and the average flow velocity was 1.45 × 10−4 m/s. Other input data are given in [106]. The anode is to the left. From ref. [106]; copyright Wiley‐VCH GmbH; reproduced with permission.

In other research groups, COMSOL Multiphysics simulations were used to investigate diffusion‐dependent focusing regimes in peak mode counterflow ITP [108], sample dispersion in ITP [109], analyte stacking in a continuous sample flow interface under conditions approaching quantitative electrokinetic injection from the entire sample volume [110], mass transport in a micro flow‐through vial of a junction‐at‐the‐tip in CE‐MS interface [111], pressure‐assisted capillary electrophoresis frontal analysis for faster binding constant determination [112], the scaling behavior in on‐chip field‐amplified sample stacking [113], and electrolysis phenomena occurring at the interface between electrode and electrolyte that have an impact on the electrode environment [114]. Very recently, a 2D mathematical model of electromigration that considers the deviation from electroneutrality in the diffuse layer of the double layer when the liquid phase is composed of a solution of weak electrolytes of any valence and complexity was developed, integrated into COMSOL Multiphysics and applied to the prediction of electromigration in nanochannels [115]. Its outcome was demonstrated by the numerical simulation of the double layer composed of a charged silica surface and an adjacent liquid solution composed of weak multivalent electrolytes. The validity of this model is not limited to the diffuse part of the double layer but is valid for electromigration of electrolytes in general.

Kler et al. designed 3D high‐performance simulation models to study electrokinetic flow and transport [116] and electrophoretic separation techniques [117] in microfluidic chips that were performed with a parallel multiphysics code referred to as PETs‐FEM within a Python programming environment. These models employed parallel computing and were found to compute faster than COMSOL Multiphysics. As this kind of approach is not widely available to the scientific community, an open‐source toolbox for electromigration separations using the platform OpenFOAM was developed [118]. The OpenFOAM (Open Field Operation and Manipulation) platform of OpenCFD Ltd (https://www.openfoam.com) is a free open‐source project for the solution of multiphysics problems based on the finite volume method. It offers native 3D support, EOF calculation, automatic parallel and supercomputing support, and is licensed under the GNU general public license. Recently, a complete mathematical model for electromigration in paper‐based analytical devices was constructed [119] and applied under the same framework [120].

Flow dynamics in ITP, i.e., the impact of the mismatch of EOF between leading and trailing electrolytes on boundary distortion in ITP with concomitant EOF, was assessed in the research group of Hardt et al. by means of 2D numerical modeling. A finite‐element model implemented on COMSOL Multiphysics provides EOF‐driven flow patterns and ion concentration profiles [121]. The same program was employed to validate an analytical approximation for the flow field in the vicinity of an ITP transition between electrolytes of different mobilities. The resulting convective ion transport inherently reduces the resolution of ITP separations [122]. Sample dispersion in ITP with Poiseuille counterflow was studied with a 2D finite‐volume model and used to validate a 1D model for the area‐averaged concentrations that is based on a Taylor–Aris‐type effective axial diffusivity [123].

In other efforts, Monte Carlo simulations were used to describe electrodynamic dispersion of sample streams in free‐flow zone electrophoresis [124] and broadening of analyte streams due to a transverse pressure gradient in free‐flow IEF [125]. An OpenFOAM fluid dynamic simulation model was developed for the description of the local interaction of hydrodynamics and Joule heating in annular CEC [126]. IEF in a multielectrode OFFGEL setup was simulated in a 2D domain with COMSOL by solving the time‐dependent Nernst–Planck differential equation [127]. Previous modeling of IEF in an OFFGEL multicompartment cell was implemented on the finite element commercial software Flux‐Expert (Astek Rhône‐Alpes, Grenoble, France) [128]. A simple, diffusion‐based mathematical model for dynamic computer simulation of free‐flow zone electrophoresis was developed and implemented in the Delphi XE2 environment. The model was used to simulate operational parameters (e.g., electric field, flow rate, and pH) for the prediction of amino acid [129] and protein [130] separations. Models to study the electrokinetic transport (mainly EOF) at intersections of microfluidic chips, such as that of Yang et al. implemented into the CFD‐ACE+ solver and used to study geometry and voltage parameters to avoid sample leakage [131], are not elaborated here.

Applications

The applications discussed in this chapter refer to simulations performed with the one‐dimensional simulators GENTRANS, SIMUL5, and SPRESSO. Work with multi‐dimensional applications is referred to in Section 2.2. together with the respective solvers.

Zone electrophoresis, isotachophoresis, moving boundaries, and analyte stacking

Dynamic simulation with SIMUL5 [132] and GENTRANS [133] was employed to predict the developments of peaks during the CZE separation of analytes that were monitored with setups comprising arrays of 16 and 8 contactless conductivity detectors, respectively, along the capillary. In the second case, the GENTRANS code featuring Taylor–Aris diffusivity to account for dispersion due to the parabolic flow profile associated with pressure‐driven laminar flow [57] was employed and found to predict realistic conductivity detector profiles for the migration and separation of cationic and anionic analyte and system zones in presence of an imposed constant buffer flow and a small EOF. For configurations with discontinuous buffer systems, including ITP, experimental data obtained with the array detector revealed that the EOF is not constant. Comparison of simulation and experimental data of ITP systems provided the insight that the EOF can be estimated with an ionic strength‐dependent model similar to that previously used to describe EOF in fused‐silica capillaries dynamically double coated with Polybrene and poly(vinylsulfonate) [133]. SIMUL5 was employed to provide a comparison of peak shapes predicted by simulation with those of a new nonlinear model of electromigration for the analysis of two co‐migrating fully dissociated analytes [134]. For the systems studied, an almost perfect agreement was obtained. When the velocity of the separating analyte depends on the concentration of the co‐analyte, the consequence is a mutual influence and a distortion of both analyte zones.

Unwanted effects of carbonate in CZE of anions at high pH were assessed with SIMUL5 and validated experimentally [135]. Carbon dioxide absorbed from air into BGEs and samples is thereby shown to form additional zones and/or boundaries that may induce strong and pronounced temporary changes in the migration of analytes. The presence of millimolar amounts of carbonate in an alkaline BGE (i) changes the pH of the BGE which results in changes of effective mobilities, (ii) produces system peaks, dips, or more complex disturbances in the detection signal, and (iii) interacts with the sample components that can substantially modify the course of the separation process and its result. Detailed investigation of possible effects of carbonate for a given sample and BGE by computer simulation is very useful to prevent problems associated with the analysis of specific analytes. The computer prediction of the formation and migration of system peaks was also the subject in other CZE investigations [24, 57, 58, 76, 136, 137]. In one particular study, online preconcentration of weak electrolytes at the pH boundary induced by a migrating system zone was investigated using SIMUL5 which provided a detailed understanding of the process [136]. In another interesting work, system peaks were deliberately generated via the addition of Li and Cs into the BGE for use as internal standards to improve quantification in microchip electrophoresis with conductivity detection [137].

As discussed in our previous review, dynamic simulation can effectively be used to investigate the electrokinetic processes occurring during the different modes of analyte stacking [18]. This application was continued in the 2010–2020 time period. Wang et al. [138] used nicotine as the test compound to compare predictions made with SIMUL5 with real experiments for two types of dynamic pH junctions. It could be demonstrated that focusing of at least 95% of the injected target molecule was achievable for both cases. In another project, simulations of pH barrage junction focusing were employed for its optimization and application to weakly alkaline or zwitterionic analytes [139]. A simulation study performed with the same simulator dealt with the concentration of weak acids in a pH boundary formed between sodium borate pH 9.5 and sodium phosphate pH 2.5 electrolytes [140]. It revealed that the preconcentration is due to dissociation changes of the analytes and transient ITP. Furthermore, computer simulations with SIMUL5 were employed to characterize the field‐enhanced sample injection process and associated formation of moving boundaries in a process that was coupled to sweeping [141]. In the context of the CZE analysis of imatinib in plasma samples, computer simulations with SIMUL5 were used to confirm the experimental results and to understand the electrophoretic behavior of imatinib in the presence of NaCl [142].

In another laboratory, SIMUL5 simulations were employed to examine several electrolyte configurations: investigation of the role of counter‐ions in the BGE for analysis of cationic weak bases and amino acids in neutral aqueous solutions by CZE with electrokinetic injection [143]; to study the role of the counter‐ions in the BGE in transient ITP‐CZE [99]; to stack phenol from seawater samples using an improved dynamic pH junction for the determination of submicromolar levels of this analyte [144]; and to characterize transient ITP for the concentration of aniline and pyridine from sewage samples while having water as anolyte [145]. During the latter process, a zone comprising H+ as system‐induced terminating compound is produced and analytes become stacked between this zone and the BGE. The simulation data presented in Fig. 5 were obtained with SIMUL5 and illustrate how micromolar concentrations of phenol in presence of 550 mM NaCl can be stacked and analyzed in a BGE comprising 160 mM boric acid and 126 mM NaOH (pH 9.8). A zone of sodium propionate injected after the seawater sample (Fig. 5C and D) provides improved sensitivity compared to the use of a conventional dynamic pH junction (Fig. 5A and B) [144].

Figure 5

Investigation for the improvement of the dynamic pH junction for analysis of phenol in seawater using SIMUL5. A sample composed of 550 mM NaCl and 1 μM phenol was applied between a BGE comprising 160 mM boric acid and 126 mM NaOH (pH = 9.8) for analysis with (A,B) a conventional dynamic pH junction and (C,D) an improved configuration with a mixture of 600 mM propionic acid and 600 mM NaOH as saturated fatty acid (SFA) solution on the cathodic side of the applied sample. Panels (A,C) depict the initial distributions of all components together with pH and potential gradient, whereas panels (B,D) show the distributions after (B) 200 s and (C) 230 s of application of a constant 50 V. The cathode is to the left. A 50 μm id capillary of 50 mm length divided by a mesh of 10 μm without EOF served as separation space. Other input data are given in [144]. From ref. [144]; copyright Wiley‐VCH GmbH; reproduced with permission.

SIMUL5 served as a tool to assess the removal of sample background buffering ions and myoglobin enrichment via a pH junction created by discontinuous buffers [146], and to investigate the stacking process of analytes (peptides and proteins) in a discontinuous buffer system comprising a pH 9.75 ammonium buffer as catholyte and a pH 4.25 acetate buffer as anolyte [147]. Stacking is shown to occur in a sharp, cationically migrating boundary referred to as the neutralization reaction boundary. The simulated results closely resembled the experimental data, and together, they effectively revealed the characteristics of the discontinuous buffers. SIMUL5 was employed to optimize the ITP concentration conditions of nucleic acid fragments via analysis of various operating factors, including electric current, sample amount, and sample matrix [148]. Under the elucidated conditions, DNA focusing was studied in a commercial ITP instrument with preparative fraction collection. SIMUL5 was applied to investigate the effects of local conductivity differences between analyte plugs, reagent plugs, and the BGE on EMMA analyses, i.e., via analysis of the ionic boundaries that are formed upon power application [149]. Furthermore, GENTRANS that encompasses Taylor–Aris diffusivity was used in the diffusion mode to predict transverse diffusional reactant mixing occurring during hydrodynamic plug injection of configurations featuring four and seven plugs. The same simulator in the electrophoretic mode was applied to study electrophoretic reactant mixing caused by voltage application in absence of buffer flow. Simulations also visualized buffer changes that occur upon power application between incubation buffer and background electrolyte that have an influence on the reaction mixture [150].

SIMUL5 was employed to simulate analyte peak splitting of salt‐containing samples that may occur in CZE with sample self‐stacking [151]. In this comprehensive work, theoretical considerations based on velocity diagrams and computer simulations reveal that these effects originate in the transient phase of separation where electromigration dispersion profiles and sharp boundaries are formed and evolve. Multiple transient sharp boundaries (including system boundaries) may exist that are simultaneously capable of stacking an analyte resulting in permanent or transient multiple peaks. This is illustrated with the simulation data presented in Fig. 6. A 50 mm long sample with 0.01 mM benzoic acid in 59.80 mM HCl and 59.88 mM NaOH (pH 9.80) was applied into a BGE comprising 50 mM boric acid and 32 mM NaOH (pH 9.48). Upon power application, three transient, sharp migrating anionic boundaries are formed and benzoic acid becomes stacked in two of these boundaries and therefore becomes split into two peaks. Eventually, two of these sharp boundaries disappear and the two benzoic acid peaks become destacked and migrate zone electrophoretically in the BGE behind the chloride zone (1138 s time point of Fig. 6). The computer predicted peak splitting of benzoic acid in this system could be validated experimentally. The process of multiple peak formation is complex, depends on the amount of sample, the composition of the sample, and the composition of the BGE. Thus, the detected analyte pattern may vary from sample to sample and may depend on detector location. The elucidation of the peak‐splitting mechanism by simulation allows both the identification of its presence in a given BGE and sample, and to find ways to avoid it [151].

Figure 6

Peak splitting in CZE with sample self‐stacking. Computer simulated concentration profiles of chloride, borate, and benzoate at indicated time points (0 to 1138 s) with a sample composed of 59.80 mM HCl, 59.88 mM NaOH, and 0.01 mM benzoic acid (pH 9.80) applied into a BGE comprising 50 mM boric acid and 32 mM NaOH (pH 9.48). The initial sample length was 50 mm. The vertical axes show the concentration values for chloride and borate at the left and for benzoate at the right. The simulations were performed with SIMUL5 in a capillary of 200 mm length divided into 6000 segments and application of a constant 500 V. The input data of all components are given in ref. [151]. The anode is to the right. From ref. [151]; copyright Wiley‐VCH GmbH; reproduced with permission.

The fundamentals of electrokinetic injection of a weak acid across a short water plug into a phosphate buffer at low pH were studied by computer simulation with GENTRANS and validated experimentally [152, 153]. The current during electrokinetic injection, the formation of the analyte zone, changes occurring within and around the water plug, and mass transport of all compounds in the electric field were investigated. The simulation provided insight into these changes, including the nature of the migrating boundaries and the stacking of the weak base. The data confirm the role of the water plug to prevent contamination of the sample by components of the background electrolyte and suggest that mixing caused by electrohydrodynamic instabilities increases the water plug conductivity. The sample conductivity must be controlled by the addition of an acid to create a buffering environment and to prevent the generation of reversed flow that removes the water plug [152]. Simulations of the behavior of the electrophoretic system after electrokinetic injection of cationic compounds across a short water plug revealed that a phosphoric acid zone at the plug‐buffer interface is formed and becomes converted into a migrating phosphate buffer plug that corresponds to the cationically migrating system zone of the phosphate buffer system. Its mobility is higher than that of the analytes such that they migrate behind the system zone in a phosphate buffer comparable to the applied BGE [153]. The involvement of GENTRANS simulations encompassed also studies of electroosmotic flow‐balanced isotachophoretic stacking with continuous electrokinetic injection for the concentration of anions in high conductivity samples [154], pressure‐assisted electrokinetic supercharging for the enhancement of non‐steroidal anti‐inflammatory drugs [155], acid‐induced transient isotachophoretic stacking of basic drugs in co‐electroosmotic flow CZE [156], and stacking of monosaccharides after hydrolysis of a single wood fiber in an open microchannel [157].

In the 2010–2020 time period, SPRESSO was extensively used for the simulation of ITP systems. SPRESSO was extended with the Onsager‐Fuoss model for actual ionic mobility and the extended Debye–Hückel theory for correction of ionic activity [66]. For model ITP and CZE separations performed at low and high ionic strength, this version of SPRESSO is reported to predict data that compare well with those monitored experimentally. The quasi 1D version for variable cross‐section channels was applied to simulate focusing and separation of two analytes in plateau mode ITP and peak mode ITP prior, during, and after traveling through a converging channel with a fivefold cross‐sectional area reduction [67]. This is shown with the cationic ITP simulation data presented in Fig. 7. Pyridine and aniline are separated between a leading electrolyte composed of 18 mM NaOH and 20 mM acetic acid and an anolyte comprising 40 mM β‐alanine and 20 mM acetic acid. Figure 7A depicts simulated patterns with sample amounts that are sufficient to form plateau zones. Initially, the two analytes are injected as overlapping Gaussian peaks at the interface of leading and terminating electrolytes. Upon application of a high electric field strength (1207.7 A/m2), the two analytes separate quickly in the large cross‐section region and form zones with a plateau concentration (t = 27 s). Pyridine has a higher mobility and forms its zone in front of aniline. As the analyte zones migrate from the large to the small cross‐section region (t = 46 s), their zone lengths increase inversely with a decrease in cross‐sectional area. When analyte zones fully migrate into the narrow cross‐section region (t = 53 s), their zone lengths reach new steady‐state values that are five times longer compared to those in the large cross‐section region. It is important to note that the plateau concentrations of the analytes do not depend on channel cross section. They are given by the area‐independent Kohlrausch adjustment. The thickness of the steady‐state migrating boundaries, however, becomes smaller with the increase of current density from 1207.7 to 6038.6 A/m2. This is, however, not visible in the graphs with the used x‐axis scale. Figure 7B depicts isotachophoretic focusing of the same two analytes in the peak mode and along the same converging channel. Compared to the case of Fig. 7A, the amounts of analytes are 20‐fold smaller and the applied current 10‐fold lower. Data for three locations are presented. Due to their low initial concentrations, the two analytes focus in the peak mode in the large cross‐section region (t = 241 s). Analyte peaks largely overlap and migrate essentially within the boundary formed by the leading and the terminating ions. As the analyte zones migrate from the large to the small cross‐section region the analyte peak concentrations increase due to increased separation and zone boundaries become sharper that is associated with the higher local electric field strength (t = 452 s). In the small cross‐section region, their zone concentrations nearly reach the plateau values (t = 517 s). During the transition from the large to the small cross‐section regions, the electric field strength changes from 120.8 to 603.9 A/m2, and a higher detection sensitivity is reached. The quasi‐1D approach of SPRESSO is computationally efficient and highly suited to simulate systems comprising very sharp boundaries at high electric field strengths and in channels with a variable cross‐sectional area [67].

Figure 7

Simulation data obtained with the SPRESSO quasi‐1D model that accounts for the dispersive effect of a nonuniform electric field in channels with a variable cross‐sectional area. (A) Plateau mode cationic ITP and (B) peak mode cationic ITP separation of pyridine (S1) and aniline (S2) between sodium (LE) and β‐alanine (TE) while traveling through a converging channel with fivefold cross‐sectional area reduction. The leading electrolyte was composed of 18 mM NaOH and 20 mM acetic acid, and the terminating electrolyte comprised 40 mM β‐alanine and 20 mM acetic acid. Simulations used a 60 mm long computational domain with 450 grid points and a constant applied current of (A) 1 μA and (B) 0.1 μA. The amounts of sampled analytes for the peak‐mode case were 20‐fold lower compared to that of the plateau mode. The cathode is to the right. From ref. [67]; copyright Wiley‐VCH GmbH; reproduced with permission.

The coupling of transient ITP and CZE for increasing sensitivity and resolution was comprehensively reviewed by Bahga and Santiago [158]. The methods are compared and discussed with simulation examples generated with SPRESSO. The same simulator was employed for validation of analytical solutions derived for analyte distributions in peak mode ITP [159], the plateau concentrations predicted by a diffusion free ITP model for channels with axially varying cross‐section [78], and indirect detection in ITP assays on a miniaturized system with LIF detection [160]. In other efforts, SPRESSO was used to simulate focusing and separation of ionic species using bidirectional ITP [161]. In the described process, anionic sample species are focused isotachophoretically prior to the interaction with a cationic ITP shock that changes the local pH and other properties such that analytes become defocused and begin to separate under zone electrophoretic conditions. This principle was applied to the high‐resolution separation of a 1 kbp DNA ladder. Bidirectional ITP was also applied to the formation of a sudden decrease in the concentration of the leading electrolyte ahead of the focused anions and thus a corresponding decrease in analyte zone concentrations and increase in the length of the analyte zone with plateau concentration. The latter aspect comes along with an increased detection sensitivity as the zone length is proportional to the amount of analyte. Simulations were used to describe the species transport of the involved processes in such a system with a concentration cascade of leading electrolytes [162].

Gebauer and co‐workers reported a number of interesting contributions in which SIMUL5 was used to verify the behavior of electrolyte systems that are compatible for highly sensitive analysis of cations [163, 164, 165] and anions [166, 167, 168, 169] by capillary ITP and moving boundary systems coupled to ESI‐MS detection. They represent innovative contributions in which the properties of the stacking ITP boundary can be tuned by the composition of both the leading and terminating zone, as well as the proper selection of a suitable spacing component. These approaches permit the monitoring of sub‐ppb analyte concentrations in the peak mode stacking format in real‐world samples, such as in drinking and river water. Simulation data of three moving boundary systems comprising ITP stacking areas at the front and rear boundaries of a spacer are presented in Fig. 8. These examples demonstrate how the spacer technique in moving‐boundary ITP allows playing with selectivity [168]. A mixture of three model analytes with pKa values in the medium acidic region, benzoic acid, salicylic acid, and p‐hydroxybenzoic acid, together with a spacer served as analytes. All three examples represent moving boundary ITP systems in which one or both electrolytes contained also the other component. Formic acid and propionic acid were anionic co‐ionic constituents and NH4 + served as a counter component. Data of a system with the leading zone containing a mixture of formic acid and propionic acid, the terminating zone propionic acid only, and lactic acid as a spacer, are depicted in Fig. 8A. In this configuration, benzoic acid becomes stacked at the rear boundary of the spacer zone, salicylic acid at the front boundary of the spacer zone, whereas p‐hydroxybenzoic acid does not form an ITP zone and migrates as a broad zone in the terminating electrolyte. Figure 8B shows the simulation result for another moving‐boundary ITP system where both electrolytes comprise mixtures of formic acid and propionic acid, and acetic acid is used as a spacer. In this system, p‐hydroxybenzoic acid and benzoic acid are focused in the rear and front boundaries of the spacer zone, respectively, and salicylic acid is not stacked and migrates zone electrophoretically within the leading electrolyte. In the case where the leading zone contained formic acid only, the terminating zone both anions, and lactic acid served as a spacer (Fig. 8C), the simulation predicts the stacking of salicylic acid in the front boundary and the two other analytes in the rear boundary. These simulations could be validated experimentally by CE‐MS [168]. They illustrate that proper tuning of the electrolytes and the selection of the spacer compounds are crucial for ITP stacking of a target analyte. Such systems enhance the application range of ITP‐MS.

Figure 8

Computer simulation of separations in moving‐boundary ITP systems with spacers. The concentration profiles depict the situations after 3 s of electrophoresis from the starting situation with a 0.25 mm sample zone containing 0.01 mM benzoic acid (Ben), 0.01 mM salicylic acid (Sal), 0.02 mM p‐hydroxybenzoic acid (PHBen) and (A) 6 mM lactic acid (Lac), (B) 4 mM acetic acid (Ac), and (C) 4 mM Lac. The compositions of the leading (L) and terminating (T) zones were: (A) L: 8 mM formic acid (For), 2 mM propionic acid (Pro), and 6 mM NH4 + (pH 4.16), T: 10 mM Pro (Kohlrausch adjusted values: 8.24 mM Pro, 3.87 mM NH4 +, pH 4.82); (B) L: 9.5 mM For, 0.5 mM Pro, and 8 mM NH4 + (pH 4.45), T: 6 mM Pro and 4 mM For (adjusted values: 5.38 mM Pro, 3.23 mM For, 6.35 mM NH4 +, pH 5.07); (C) L: 10 mM For and 8 mM NH4 + (pH 4.37), T: 7 mM Pro and 3 mM For (adjusted values: 6.02 mM Pro, 2.28 mM For, 5.97 mM NH4 +, pH 5.11). The concentration profiles of the analytes, spacers, and system anions are shown by solid, dashed, and dotted lines, respectively. The simulations were performed with SIMUL5 using 2000 segments along a 5 mm separation space and 200 V. Input parameters of all components are given in ref. [168]. The anode is to the right. From ref. [168]; copyright Wiley‐VCH GmbH; reproduced with permission.

A new capillary electrophoretic separation and focusing principle in which weak nonamphoteric ionogenic species are focused and transported into or through the detector was explored by Gebauer and co‐workers [170, 171, 172, 173]. In that work, SIMUL5 was used to simulate the entire process and validate the theoretical prediction based on calculated velocity diagrams. The prerequisite condition for the application of this principle is the existence of an inverse electromigration dispersion (EMD) profile, i.e., a profile in which pH is decreasing toward the anode or cathode for focusing of anionic or cationic weak analytes, respectively. The conditions under which a weak acid is focused on a profile of this type are depicted schematically in Fig. 9. It shows the distributions of pH, conductivity κ, and velocity v along the separation column at time t after the application of current. The starting configuration is assumed to be a sharp boundary between the leading solution (zone L) and the trailing solution (zone T). The compositions of both zones are chosen such that they are equal to the compositions of the front and rear edges of the EMD profile that evolves between them after the application of current. Under a constant current density, the front and rear edges of the profile move with the constant velocities v L and v T, respectively, where v L > v T. All points along the gradient have velocities in between those values. In Fig. 9, this is represented by the v j line that illustrates the velocity profile for the depicted time point. The line marked v A shows the dependence of the migration velocity for a weak acid A that becomes focused at the location of the intersection of the two velocity graphs. The orientation of the pH gradient with higher pH at the trailing edge and lower pH at the fronting edge is the major prerequisite for the focusing of the analyte. The EMD zone with the focused analyte is migrating and dispersing such that all property profiles are flattening with time. Unlike in ITP, no migrating steady‐state distributions are produced. The peak height of the focused analyte zone decreases as it continues to migrate toward a detector. It can be assumed that its concentration profile at any time is the result of the balance between the flux of electromigration and diffusion.

Figure 9

Focusing of an anionic analyte within an inverse electromigration dispersion profile. Schematic representation of pH, conductivity κ, velocity v, and analyte concentration cA along the migration coordinate for an anionic electromigration dispersion profile after time t evolving from a sharp initial boundary between the leading zone L and the terminating zone T. The anionic analyte focuses on the profile at point xA0 where its electrophoretic velocity v A is equal to v j of the gradient. EMD refers to electromigration dispersion. The anode is to the right. From ref. [173]; copyright Wiley‐VCH GmbH; reproduced with permission.

The simulation data presented in Fig. 10 depict the concentration profiles of 12 weak acids after 10 s of power application that were applied in small quantities between a leading electrolyte composed of 15.2 mM maleic acid and 21.7 mM 2,6‐lutidine (pH 5.81) whereas the trailing electrolyte contained 2 mM maleic acid and 27.6 mM 2,6‐lutidine (pH 7.53) [171]. The data illustrate that (i) weak acids with pKa values between 9 and 5.5 can be quickly focused in the produced EMD profile, (ii) the focusing power becomes smaller with decreasing pKa value, and (iii) migration proceeded toward the anode on the right. Furthermore, the effect of carbonate on the focusing properties was simulated in the same way and shown to have an appreciable impact for weak acids with pKa values < 6.5 [171]. This configuration was found to be robust and was successfully applied to the analysis of nanograms per milliliter concentrations of sulfonamides in waters using CE–MS [171]. The fundamental resolution equation and pressure‐assisted performance enhancement were reported in ref. [172] and a complete theory of analyte focusing on an inverse EMD profile is given in ref. [173]. It is important to note that this separation principle is based upon a mechanism that is different from both CZE and IEF where separation is based on the difference in mobilities and pI values, respectively.

Figure 10

Focusing of weak acids within an inverse electromigration dispersion profile predicted by SIMUL5. Computer simulated concentration profiles of 12 analytes consisting of weak acids with an ionic mobility of 30 × 10 9 m2/Vs and with pKa values between 4.5 and 10 (peak labels) after 10 s of electrophoresis time with an applied voltage of 2500 V across a separation column of 25 mm length. The leading electrolyte was composed of 15.2 mM maleic acid and 21.7 mM 2,6‐lutidine (pH 5.81) whereas the trailing electrolyte contained 2 mM maleic acid and 27.6 mM 2,6‐lutidine (pH 7.53). The sample contained 2 × 10 8 M of each analyte and leading electrolyte, was placed between the two electrolytes at the cathodic column side and had a length of 1.25 mm. The distributions of 2,6‐lutidine (Lut, lower panel), conductivity and pH (both upper panel) are presented for the same time point. The anode is to the right. From ref. [171]; copyright Elsevier; reproduced with permission.

Isoelectric focusing

Although IEF simulations can be performed with many different simulators [27, 41, 60, 71, 72, 74, 85, 86, 87, 88, 100, 117, 118], GENTRANS was by far the most used solver to predict pH gradient formation, stabilization, destabilization, and migration, as well as focusing and mobilization of analytes in IEF [16, 17, 18, 19, 20, 21, 22]. The value of using dynamic simulation for the investigation of pH gradient drifts and instabilities is the subject of a recent comprehensive review with a large number of simulation examples produced by GENTRANS [22]. It discusses the electrokinetic processes that lead to pH gradient instabilities in carrier ampholyte‐based IEF. In addition to electroosmosis, there are four of electrophoretic nature, namely (i) the stabilizing phase with the plateau phenomenon, (ii) the gradual isotachophoretic loss of carrier ampholytes at the two column ends in presence of electrode solutions, (iii) the inequality of the mobilities of positively and negatively charged species of ampholytes (the anionic form of an ampholyte is more hydrated than the cationic form and thus possesses slightly smaller mobility [74, 174]), and (iv) the continuous penetration of carbonate from the catholyte into the focusing column. The impact of these factors on cathodic and anodic drifts was analyzed by simulation of carrier ampholyte‐based focusing in closed and open columns. Focusing under realistic conditions within a 5 cm long capillary in which three amphoteric low molecular mass dyes were focused in a pH 3–10 gradient formed by 140 carrier ampholytes was investigated and compared to experimental results. In open columns, electroosmosis displaces the entire gradient toward the cathode or anode whereas the electrophoretic processes act bidirectionally with a transition around pH 4 (drifts for pI > 4 and pI < 4 typically toward the cathode and anode, respectively). The data illustrate that focused zones of carrier ampholytes have an electrophoretic flux and that dynamic simulation can be effectively used to assess the magnitude of each of the electrokinetic destabilizing factors and the resulting drift for a combination of these effects. The resulting electrokinetic transport indicates that a true steady‐state is never attained in carrier ampholyte‐based IEF, that is, in IEF without an immobilized pH gradient [22].

In other efforts, GENTRANS was used to study three additional aspects of IEF. The first one comprises various sampling strategies for capillary IEF with concomitant electroosmotic zone mobilization [175]. The separation and focusing of analytes in a pH 3–11 gradient formed by 101 biprotic carrier ampholytes was investigated with the application of the analytes (i) mixed with the carrier ampholytes (as is customarily done), (ii) as a short zone within the initial carrier ampholyte zone, (iii) sandwiched between zones of carrier ampholytes, and (iv) introduced before or after the initial carrier ampholyte zone. This is illustrated with the data presented in Figs. 11 and 12 that were obtained with the ionic strength‐dependent EOF model for a fused‐silica capillary and focusing between 10 mM phosphoric acid as anolyte and 20 mM NaOH as catholyte. Analyte separation dynamics for sampling in presence of carrier ampholytes are depicted in Fig. 11 and those for the sample being applied before, after, and between carrier ampholytes are presented in Fig. 12. It is important to note that current density and EOF toward the cathode are not constant in these situations [175]. With sampling as a short zone within or adjacent to the carrier ampholytes, separation and focusing of analytes proceed as a cationic, anionic, or mixed process, and separation of the analytes is predicted to be much faster than the separation of the carrier components. Thus, after the initial separation, analytes continue to separate and eventually reach their focusing locations. This is different from the double‐peak approach to equilibrium that takes place when analytes and carrier ampholytes are applied as a homogenous mixture (Fig. 11A). Simulation data reveal that sample application between two zones of carrier ampholytes results in the formation of a pH gradient disturbance. As a consequence, the properties of this region are sample matrix dependent, the pH gradient is flatter, and the region is likely to represent a conductance gap (hot spot). Simulation data suggest that samples placed at the anodic side or at the anodic end of the initial carrier ampholyte zone are the favorable configurations for capillary IEF focusing with electroosmotic zone mobilization [175].

Figure 11

Dynamics of analytes in IEF sampled in presence of carrier ampholytes. Computer predicted analyte dynamics after 0, 0.2, 0.6, 1.0, and 1.4 min of power application for (A) analytes mixed with carrier components and placed between 3% and 23% of column length, (B) a 2% sample zone at the anodic end of the carrier ampholyte zone, (C) a sample zone between 9% and 11% of column length, and (D) a 2% sample zone at the cathodic end of the carrier ampholyte zone. GENTRANS with the ionic strength dependent EOF model for a fused‐silica capillary was used for simulation. Phosphoric acid (10 mM) and 20 mM NaOH served as anolyte and catholyte, respectively. A 10 cm focusing space divided into 4000 segments of equal length and a constant voltage of 1000 V were employed. A total of 101 hypothetical biprotic carrier ampholytes were used to establish a pH gradient between anode and cathode. Their pI values uniformly span the range 3.0–11.0 (ΔpI = 0.08). For each ampholyte, ΔpK was 2.5, the ionic mobility was 2.5 × 10−8 m2/Vs, and the initial concentration was 0.2 mM. Physicochemical input data for the seven amphoteric sample components and the other constituents are given in [175]. The data are presented as the sum of analyte concentrations. In the graphs at the bottom of each panel, the dotted lines demarcate the initial positions of the carrier ampholyte zones (sum of carrier component concentrations divided by 10). S and CA refer to sample and carrier ampholytes, respectively. Analyte peaks are labelled with their pI values. Data are depicted with y‐scale offsets of 4 mM. The cathode is to the right. From ref. [175]; copyright Wiley‐VCH GmbH; reproduced with permission.

Figure 12

Dynamics of analytes in IEF sampled outside carrier ampholytes. Computer predicted analyte dynamics after 0, 0.2, 0.6, 1.0, and 1.4 min of power application for a sample zone (A) at the anodic side of the carrier ampholyte zone, (B) between two zones of carrier ampholytes, and (C) at the cathodic side of the initial carrier zone. Other conditions, data presentation and key are the same as for Fig. 11. The cathode is to the right. From ref. [175]; copyright Wiley‐VCH GmbH; reproduced with permission.

In a second project, simulations with the sample being applied between zones of carrier ampholytes or on the anodic side of the carrier ampholytes were used to characterize the behavior of sample components with pI values outside a pH 6–8 gradient formed by 101 hypothetical biprotic carrier ampholytes [176]. Application of power leads to a situation in which the pH gradient is bracketed by two isotachophoretic zone structures comprising selected sample and carrier components as isotachophoretic zones. The isotachophoretic structures electrophoretically migrate in opposite direction. When electroosmosis or an imposed flow is present, the overall pattern is transported toward the capillary end for detection of the entire zone pattern. Sample components whose pI values are outside the established pH gradient are demonstrated to form isotachophoretic zones behind the leading cation of the catholyte (components with pI values larger than the upper edge of the pH gradient) and the leading anion of the anolyte (components with pI values smaller than the lower edge of the pH gradient). Amphoteric compounds with appropriate pI values or nonamphoteric components can act as isotachophoretic spacer compounds between sample compounds or between the leader and the sample with the highest mobility [176].

The third project dealt with the modeling of the formation and prevention of a pure water zone in IEF with narrow pH range carrier ampholytes [177]. Characteristics of gradients covering two pH units that end or begin around neutrality were investigated. Data obtained revealed that a zone of water is formed in focusing with carrier ampholytes when the applied pH range does not cover the neutral region, ends at pH 7.00 or begins at pH 7.00. The presence of additional amphoteric components that cover the neutrality region prevent water zone formation under current flow. Furthermore, no water zone evolves when atmospheric carbon dioxide dissolved in the catholyte causes the migration of carbonic acid (in the form of carbonate and/or hydrogen carbonate ions) from the catholyte through the focusing structure [177].

Simulation with GENTRANS was also used to predict the dynamics of pH gradient formation and protein separation in simple buffers that were applied in a binary system micropreparative IEF approach featuring a 12 μL suspended drop between two palladium electrodes. During the electrophoretic process at low voltages (1.5–5 V), fluid was allowed to evaporate until splitting into two fractions [178]. Computer simulation of protein and peptide preconcentration in carrier‐free systems and IEF in microchannels using simple ampholytes was investigated with GENTRANS [179]. In the configurations studied, the sample was uniformly distributed throughout the channel with driving electrodes used as column ends. Experimental results from carrier‐free systems are compared to simulation results, and the effects of atmospheric carbon dioxide and impurities in the sample solution are examined via simulation. Simulation data provided insight into the dynamics of the transport of all components under the applied electric field and revealed the formation of a pure water zone in the channel center. IEF with simple well defined amino acids as carrier components was assessed for concentration and fractionation of peptides and component distributions in the channel were assessed using MALDI–MS and nano‐ESI–MS [179]. Finally, important information about the focusing dynamics and location of the foci of Alzheimer's disease‐related amyloid‐beta peptides by IEF in 75 nL microchannels using simple, well‐defined buffers were obtained via computer simulation [180]. A stepwise pH gradient was tailored for focusing the C‐terminal peptides with a pI of 5.3 in the boundary formed between cycloserine and aspartyl‐histidine. Detection was performed by direct sampling of a nanoliter volume containing the focused peptides from the microchannel, followed by deposition of this volume onto a chip with micropillar MALDI targets. In addition to purification, IEF preconcentration was found to provide at least a 10‐fold increase of the MALDI–MS signal.

The dynamic simulation was used for the validation of analytical approximations designed to predict steady‐state protein distributions in IEF [87, 181], to investigate pH gradient formation and cathodic drift in microchip IEF [100], to study the effects of electromigration and electroosmosis in IEF within a silica nanofluidic channel [182], and to gain qualitative insight into the behavior of different chemical mobilization schemes in microchannel IEF [183]. In the latter work, a non‐released developer version of SIMUL5 which features the options to fix concentrations at specific boundaries and to mimic the electrode reservoirs with varying capillary cross‐sections was employed. Focusing and mobilization as a one‐step process [184] was not further explored in the 2010–2020 time period.

Electrokinetic chromatography, chiral separations, and affinity electrophoresis

In this section, simulations of electrophoretic systems comprising chemical reactions other than protolysis that appeared in the 2010–2020 time period are discussed. They include MEKC, complexation involved in chiral separations, affinity electrophoresis, and CE‐based online reaction methods.

A few reports dealing with dynamic MEKC simulations can be found in the literature. Work with the MEKC version of GENTRANS [55] provided insight into the mechanism of transient trapping in MEKC [56]. Transient trapping is a mechanism of online sample concentration and separation. It involves the injection of a short length of micellar solution in front of the sample, making it similar to sweeping in partial‐filling MEKC. Simulations revealed that the mechanism for concentration in transient trapping is indeed similar to sweeping since the analytes interact and accumulate in the micelles that penetrate the sample zone. The mechanism for separation is, however, quite unique since the concentrated analytes are trapped for a few seconds on the sample/micelle boundary before they are released as the micelle concentration becomes reduced. This induces electromigration dispersion and the separation of the analytes down a micelle gradient. Compared to sweeping MEKC, transient trapping occurs faster (1/10 of the time) and within a shorter capillary length (1/4 of the capillary length), which results in two to three times increase in sensitivity [56]. Furthermore, SIMUL5 that does not include SDS complexation reactions was employed to study boundary formation of buffer components, SDS (implemented as anion), and sample matrix components in MEKC systems in which anionic micelles are focused using sample induced transient ITP [185], sweeping under conditions with an inhomogeneous electric field and low surfactant concentration [186], and systems in which stacking induces migration time shifts of analytes [187].

The incorporation of complexation equilibria with monovalent components into GENTRANS [58] and with monovalent and polyvalent components into SIMUL5 [61, 62] provided the possibility of studying the impact of 1:1 chemical equilibria between solutes and a buffer additive with fast interactions such that they can be considered instantaneous in comparison to the time scale of peak movement. Complexation constants and specific mobilities of the formed analyte‐selector complexes are required as additional inputs. Both of the 1D dynamic simulators provide equal results when used with identical inputs. They allow the prediction of the dynamics of analyte migration and separation, the elucidation of the origin and dynamics of system peaks, and the interference of analyte and system peak migration [24, 58]. SIMUL5complex was used with uncharged and charged selectors [62], to investigate electromigration dispersion effects caused by complexation [36, 37, 62, 63], to validate the occurrence and shape of analyte and system peaks in situations with complex‐forming equilibria predicted by a generalized model of the linear theory of electromigration [24, 38, 188], and to study the impact of complex mobilities, complexation constants, pH, and selector concentration on migration order and separation of drug enantiomers [63] and profens [64].

Simulations performed with SIMUL5complex revealed the existence of unexpected and previously unexplained electromigration dispersion effects that are caused by the complexation process itself. This dispersion may occur with interactions between a charged analyte and a neutral selector for the case of a low concentration of the chiral selector and a high complexation constant [37, 62, 63]. Three examples are given in Fig. 13. The presented data depict the migration of R‐flurbiprofen in presence of β‐CD (Fig. 13A; complexation constant of 4037 L/mol), heptakis(2,6‐di‐O‐methyl)‐β‐CD (Fig. 13B; 4800 L/mol), and heptakis(2,3,6‐tri‐O‐methyl)‐β‐CD (Fig. 13C, 552 L/mol) in a pH 8.13 buffer composed of 50 mM Tris and 50 mM tricine. The concentration of the selectors was varied between 0 and 100 mM. Excellent agreement between simulation (lower graphs) and experimental (upper graphs) data was obtained [37]. SIMUL5complex was also employed to elucidate the impact of the complexation of buffer constituents with neutral agents on common buffer properties [189] and analyte and system peaks [190, 191], to provide insight into the sweeping process of a drug in presence of a neutral CD [191], to study affinity capillary electrophoresis and vacancy affinity capillary electrophoresis methods that are used for the determination of complexation constants for two cases, a fully charged analyte with a neutral selector and an uncharged analyte with a charged selector [192], and to validate the applicability of a new theoretical formula derived from partial‐filling affinity capillary electrophoresis for the determination of apparent stability constants of analyte–ligand complexes [193].

Figure 13

Complexation can be a source of electromigration dispersion. Comparison of the experimental (upper panels) and simulated (lower panels; SIMUL5complex) electropherograms for R‐flurbiprofen as the analyte with (A) β‐CD, (B) heptakis(2,6‐di‐O‐methyl)‐β‐CD, and (C) heptakis(2,3,6‐tri‐O‐methyl)‐β‐CD as the complexing agent. The complexation constants were 4037, 4800, and 552 L/mol, respectively, and the mobilities of the complexes were 8.82 × 10 9, 7.54 × 10 9, and 6.50 × 10 9 m2/Vs, respectively. The BGE was composed of 50 mM Tris and 50 mM tricine having an experimental pH of 8.13 and an ionic strength of 25.76 mM. The samples contained 0.3 mM R‐flurbiprofen. The peaks are labeled with the concentration of the selector used in the BGE. From ref. [37]; copyright Elsevier; reproduced with permission.

Systems comprising sulfated β‐CD, a multiply negatively charged selector, were simulated with SIMUL5. The executed projects included the assessment of the fundamentals of field‐amplified electrokinetic injection of weak bases for enantioselective CE [194], the characterization of isomer mixtures of this selector [65], the investigation of the enantioselective separation of weak bases in an online microanalysis configuration comprising this selector [195], and the exploration of inverse cationic ITP for separation of methadone enantiomers with this chiral selector [196]. Good agreement between simulation and experimental data was obtained in all these studies. For a weak base, simulation properly predicted the inversion of analyte migration direction in presence of a multiply negatively charged selector when the selector concentration is increased and the selector concentration interval at which the enantiomers of a weak base migrate in opposite directions [24, 65].

GENTRANS was applied to characterize the migration and separation of enantiomers in systems with neutral CDs as selectors and at relevant power levels that are used in CE experiments [58]. The separation of the enantiomers of two weak bases, methadone and 2‐ethylidene‐1,5‐dimethyl‐3,3‐diphenylpyrrolidine, and codeine (achiral compound without complexation) in presence of heptakis(2,6‐di‐O‐methyl)‐β‐CD as a selector and without selector are depicted as lower and upper graphs, respectively, in Fig. 14. Complexation constants of the enantiomers (between 344 and 474 L/mol) and all other input data are given in ref. [58]. The simulation was run at a constant current density of 32.37 kA/m2 that corresponds to a current of 63.6 μA in a 50 μm id capillary. Predicted concentration data and those adjusted for differences in absorbance at a detection wavelength of 195 nm are presented in Fig. 14A and B, respectively. Simulation data are in agreement with those monitored experimentally as is illustrated with the data of Fig. 14C. GENTRANS was also used to assess the migration order change of ketoconazole enantiomers at low pH in presence of increasing amounts of (2‐hydroxypropyl)‐β‐CD [24, 197], and to investigate electrophoretic aspects of enantiomer migration and separation of methadone in capillary CEC at low pH with an immobilized neutral selector [59]. The latter work was accomplished by using zero mobilities for the formed complexes and otherwise the same input parameters as used for free solution simulations. This mimics conditions in absence of unspecific interactions between analytes and the chiral stationary phase. Simulation data revealed that separations are quicker, electrophoretic displacement rates are reduced, and sensitivity is enhanced in CEC with on‐column detection in comparison to free solution conditions. In the same work, simulation was used to study electrophoretic analyte behavior at the interface between the sample and the CEC column with the immobilized chiral selector (analyte stacking) and at the rear end when analytes leave the environment with complexation (analyte destacking). Furthermore, simulations provided insight into an approach to counteract analyte dilution at the column end via use of a BGE with higher conductivity, and the impact of EOF on analyte migration, separation, and detection for configurations with the selector zone being displaced or remaining immobilized under buffer flow [59].

Figure 14

Simulation of chiral CZE separations at a relevant current density. Computer predicted detector responses at 6 Hz showing the peaks (A) in concentration units and (B) in concentration units adjusted to differences in absorption. Panel C depicts the experimental data obtained in a 50 μm id capillary. These data were generated for the separation of codein, methadone, and 2‐ethylidene‐1,5‐dimethyl‐3,3‐diphenylpyrrolidine in a capillary of 50 cm length having an absorption detector placed at 45 cm (90% of column length). The BGE comprised 90 mM NaOH, 132 mM phosphoric acid (calculated pH: 2.39), and 1.8 mM heptakis(2,6‐di‐O‐methyl)‐β‐CD, whereas the sample was composed of racemic methadone chloride (28.90 μM), racemic 2‐ethylidene‐1,5‐dimethyl‐3,3‐diphenylpyrrolidine iodide (24.68 μM) and codeine (33.4 μM) in 10‐fold diluted BGE without additive. The sample initially occupied 1% of column length and was placed between 3 and 4% of column length. Simulations were performed with GENTRANS using a 4 μm mesh, a constant current density of 32.37 kA/m2, and a constant EOF of 160 μm/s as is described in ref. [58]. The upper graphs, depicted with a y‐axis offset, are corresponding data obtained in absence of the chiral selector. For simulation without selector, the EOF was assumed to be 180 μm/s. MET, EDDP, and COD refer to methadone, 2‐ethylidene‐1,5‐dimethyl‐3,3‐diphenylpyrrolidine, and codeine, respectively. From ref. [58]; copyright Wiley‐VCH GmbH; reproduced with permission.

GENTRANS was employed to characterize isotachophoretic enantiomer separation and zone stability of weak bases in presence of a neutral CD as chiral selector [58, 198, 199]. The systems studied comprised acidic cationic electrolyte systems with sodium and H30+ as leading and terminating components, respectively, and acetic acid as a counter component. One contribution focused on the investigation of zone formation and stability in free solution with methadone enantiomers as analytes [198]. The simulation data presented in Fig. 15 depict the cationic separation of methadone enantiomers with (2‐hydroxypropyl)‐β‐CD complexation in a system comprising a pH 4.60 leader composed of 10.0 mM NaOH, 24.6 mM acetic acid, and 5 mM (2‐hydroxypropyl)‐β‐CD. Ten millimolar acetic acid was used as an anolyte. Computer‐predicted methadone, acetic acid, sodium (dashed line), and (2‐hydroxypropyl)‐β‐CD (dotted line) profiles along the 10 cm column after 0, 5.0, and 10.0 min are presented in Fig. 15A. Concentrations (lower graphs) and pH/conductivity (upper graphs) of the migrating zone structures after 2.5, 5.0, 7.5, and 10.0 min of power application are presented in Fig. 15B, C, D, and E, respectively. These data depict the separation dynamics of the methadone enantiomers via the formation of a mixed zone which becomes shorter as a function of time and vanishes at the completion of the separation. They also reveal that all zones within the migrating isotachophoretic zone structure have higher concentrations of the selector compared to that applied in the leader, a deviation that is caused by the migration of the charged complexes.

Figure 15

Separation of enantiomers of a weak base by chiral ITP. Computer predicted separation of methadone enantiomers with (2‐hydroxypropyl)‐β‐CD complexation in a pH 4.60 leading electrolyte composed of 10.0 mM NaOH, 24.6 mM acetic acid, and 5 mM (2‐hydroxypropyl)‐β‐CD. The anolyte contained 10 mM acetic acid. (A) Profiles of methadone, acetic acid, sodium (dashed line), and (2‐hydroxypropyl)‐β‐CD (dotted line) along the column after 0, 5.0, and 10.0 min (displayed with a y‐scale offset of 50 mM). Arrows mark the position of the sample. Concentrations (lower graphs) and pH/conductivity (upper graphs) of the migrating zone structures after 2.5, 5.0, 7.5, and 10.0 min of power application are presented in (B), (C), (D), and (E), respectively. The inset in (A) depicts 0.5 min data showing the formation of the mixed zone at the initial sample/leader interface. Simulations were performed in a 10 cm column divided into 20 000 segments (5 μm mesh) with the sample being placed between 5 and 6% of column length, a constant current density of 250 A/m2 and without any EOF. The cathode is to the right. S, R, M, CD, HAc, L, and T refer to S‐methadone, R‐methadone, mixed analyte zone, (2‐hydroxypropyl)‐β‐CD, acetic acid, leader, and adjusted terminator, respectively. From ref. [198]; copyright Wiley‐VCH GmbH; reproduced with permission.

Another study employed norpseudoephedrine stereoisomers as analytes and dealt with zone formation, enantiomer separation, and migration for cases with the neutral selector added to the leader, immobilized to the capillary wall or support, or partially present in the separation column [199]. For the cases of a free and an immobilized selector, the conducted study focused on the electrophoretic transport of the analytes from the sampling compartment into the separation medium with the selector, the formation of a transient mixed zone, the separation dynamics of the stereoisomers, and the zone changes occurring during the transition from the chiral environment into a selector free leader. Dynamic computer simulation also provided a mean to investigate the dependence of the leader pH, the ionic mobility of the weak base, the mobility of the complex, the complexation constant, and selector immobilization on steady‐state plateau zone properties in a hitherto unexplored way.

In the 2010–2020 time period, other software packages were used to characterize interactive systems under CE conditions. One is SimulChir that represents a special, unreleased version of SIMUL that includes the full dynamics of interconverting enantiomers [200]. This software is based upon solving a complete set of continuity equations for all constituents of the separation system together with complexation and acid–base equilibria. It was used to simulate the dynamics of the interconversion of enantiomers in chiral separation systems and to determine the rate constants of interconversion of oxazepam enantiomers separated with carboxymethyl‐β‐CD as chiral selector. Furthermore, a simplified version of SimulChir was developed for use with multiple chiral selectors in which the electrophoretic velocity of each enantiomer is regarded constant. This solver was applied to study the interconversion of lorazepam enantiomers in presence of highly sulfated β‐CD that comprises multiple charged chiral selectors [201]. The two SimulChir versions were applied to assess accuracy and sensitivity of the determination of rate constants of interconversion in achiral and chiral environments featuring a single, well‐defined chiral selector and a mixture of selectors [202]. An overview of simulation models that describe interconversion is given by Trapp [203], a survey that also includes the direct calculation method based on approximation functions and a unified equation. DCXplorer is a software tool for the determination of interconversion barriers that utilizes the analytical solution of the unified equation and can be applied to dynamic chromatography and electrophoresis [204]. This software was recently used to investigate the enantiomerization barriers of the phthalimidone derivatives EM12 and lenalidomide by dynamic EKC [205].

Finally, the affinity electrophoresis model of Fang and Chen [206, 207, 208] that describes affinity interactions in CE under simplified electromigration conditions (assumption of constant electric field strength throughout the column) was employed to study the migration of interacting drug enantiomers in CE [209] and to test a mobility‐based correction method for accurate determination of binding constants by CE frontal analysis [210]. Furthermore, a simulator that predicts the characteristics of a moving chelation boundary was developed and applied to the sweeping of metal ions, such as Cu2+, with ethylenediaminetetraacetic acid as a complexing agent [211].

Concluding remarks

Dynamic simulators are able to predict the movement of ions in solution under the influence of an electric field and are as such the most versatile tools to explore the fundamentals of electrokinetic separations. The state of dynamic computer simulation software and its use has progressed significantly over the 2010–2020 decade. In the considered time period, the three most used simulators and the first 2D model were extended and new one‐ and multi‐dimensional models were developed (Tables 1 and 2). Efforts were geared towards the design of computation schemes that enable faster simulations (e.g., SIMUL6 and SPYCE), simpler access to multi‐dimensional simulations, and the availability of new features, including the incorporation of complexation for simulation of chiral separations and the deviation from electroneutrality for simulation of electrophoresis in nanochannels. While new solvers were mainly applied to benchmark simulations, GENTRANS, SIMUL5, and SPRESSO were applied to a large number of relevant investigations that provided insights into the behavior of analytes and buffer systems in moving boundary electrophoresis, CZE, CGE, ITP, IEF, EKC, ACE, and CEC.

Simulations led to the exploration of basic phenomena in CZE (most notably the occurrence and use of system peaks, impact of flow with a parabolic flow profile, and analyte stacking and destacking), in ITP (the behavior of analytes in plateau and peak mode ITP before, during and after traveling through a converging channel, and the effect of counterflow on zone shape), and in IEF (pH gradient drifts and instabilities, sampling strategies, and occurrence of water zones). Other topics include the characterization and optimization of new separation and focusing systems (focusing of weak electrolytes in an inverse EMD profile and in specifically tuned moving boundary systems comprising ITP stacking areas for CE–MS analysis of trace amounts of selected analytes), the discovery of analyte dispersion based on complexation, the dynamics of chiral separations in CE, ITP, and CEC with neutral and charged selectors, and the detailed description of EMMA and online systems, including reagent mixing, product separation, and formation of the complex system of moving boundaries that evolve upon current application.

With the availability of user‐friendly solvers, dynamic computer simulations can be employed by almost anyone with a basic knowledge of electrophoresis. Furthermore, based on the recent developments and achievements, it is anticipated that this area will continue to grow over the next years, to lead to new discoveries and to provide important insights that will allow the optimization of electrophoretic systems on any scale. In addition, dynamic simulation can be used as attractive tool for educational purposes.

The authors have declared no conflict of interest.

The authors have declared no conflict of interest.