Anal. Chem. 2000, 72, 4353-4362
Monte Carlo Model of Nonlinear Chromatography Francesco Dondi,* Pietro Munari, Maurizio Remelli, and Alberto Cavazzini
Department of Chemistry, University of Ferrara, via L. Borsari 46, I-44100 Ferrara, Italy
A stochastic approach to the nonlinear chromatography theory, based on the Monte Carlo simulation method, is presented. A computer program, acting as a “virtual chromatograph” and performing a discrete event simulation, is described. Such a program allows one to choose the column type, operating conditions, sample composition, injection method, mobile-phase dispersion model, and stationary-phase sorption-desorption kinetics. Nonlinearity is accounted for by continuously monitoring and updating both the column and the solute status and by moving individual molecules step by step along the column according to specific random modes. The program has been validated through a series of statistical tests and comparing the results with the well-known achievements of the classical stochastic theory. A first application is presented, referred to a real case benzene elution on a gas solid capillary column, where the Langmuir adsorption isotherm is assumed. The effect of both the sorption modes and the site capacity are investigated. Possible applications to investigate open problems in several fields of separation science are emphasized. In addition, several specific points such as the down-scaling of a real case and the correspondence of specific adsorption dynamics with the equilibrium Langmuir isotherm are described. Chromatography is universally recognized as the most powerful fractionation technique for molecular species in the wide range from 2 to 106 Da. Theoretical modeling of the many physicochemical aspects of chromatography has always been a vital input for its development.1-5 The most significant results have been obtained in describing the chromatographic process under the linear conditions usually followed in the “analytical” applications, where the amount injected is very low with respect to the amount of the stationary phase. Several physicomathematical approaches are employed in analyzing the above-mentioned simplified chromatographic problem, called the linear model. They belong to three categories: mass balance approach; nonequilibrium ap* Corresponding author: (e-mail)
[email protected]; (fax) +39 0532 240709. (1) Giddings, J. C. Dynamics of Chromatography; M. Dekker: New York, 1965. (2) Giddings, J. C. Unified Separation Science; J. Wiley and Sons, Inc.: New York, 1991. (3) Guiochon, G.; Golshan Shirazi, S.; Katti, A. M. Fundamentals of Preparative and Nonlinear Chromatography; Academic Press: Boston, 1994. (4) Dondi, F.; Blo, G.; Remelli, M.; Reschiglian, P. In Theoretical Advancement in Chromatography and Related Separation Techniques; Dondi, F., Guiochon, G., Eds.; NATO ASI Series C 383; Kluwer Academic: Dordrech, 1992; pp 173-210. (5) Shure, M. R. In Advances in Chromatography; Brown, P. R., Grushka, E., Eds.; Marcel Dekker: New York, 1998; Vol. 39, pp 139-200. 10.1021/ac0003347 CCC: $19.00 Published on Web 08/12/2000
© 2000 American Chemical Society
proach; the random walk and stochastic approach.1-3 Recently the stochastic approach developed through the powerful characteristic function mathematical method has proved able to handle column heterogeneity effects even in combination with mobile-phase dispersion.6,7 When the injected amount increases, the partition isotherm is no longer linear and peak shape significantly deviates from Gaussian symmetry. In general, the theory of nonlinear chromatography is more complex because of the interdependence of the stationary-phase and mobile-phase concentrations. Significant theoretical advances have been attained in both theoretical and practical understanding of nonlinear chromatography,3 but several unsolved problems still remain, primarily the problems of multicomponent mixture elution and the combined effects of nonlinearity coupled with stationary-phase heterogeneity. The abovementioned stochastic approach based on the characteristic function4 does not appear able to handle the combined nonlinearity and heterogeneity effects since it is strictly based on the assumption that the individual molecules behave independently. The Monte Carlo (MC) methods,5,8 although based on stochastic description of chromatography, are in principle able to handle the above-mentioned nonlinearity effects. In fact, the method consists of simulating the behavior of a particle set progressing through the column as objects performing phase changes, between a stationary phasesconsidered as an array of sorption sitessand a mobile phase, at random time intervals. There is no a priori limit in considering short columns with different energy sites over which chemically different solute molecules move. If the simulation model is designed to take into consideration particle-particle interaction in the stationary phase, MC computation should be able to investigate the many aspects of both linear and nonlinear chromatography still unexplored. To date, only MC stochastic transport simulation of a single particle has been reported: the whole simulation must be repeated nr times in order to build up a peak corresponding to the elution of nr particles.5,8 This corresponds to a linear transport condition since no particle-particle interaction can be considered. Under these conditions, MC simulation only represents an easy form of the classical stochastic model of linear chromatography, freed of mathematical complexity. A MC model of nonlinear chromatography should instead be able to master the simultaneous stochastic transport of nm molecules, accounting for their reciprocal interac(6) Cavazzini, A.; Remelli, M.; Dondi, F.; Felinger, A. Anal. Chem. 1999, 71, 3453-3462 (7) Felinger, A.; Cavazzini, A.; Remelli, M.; Dondi, F. Anal. Chem. 1999, 71, 4472-4479. (8) Phillips, J. B.; Wright, N. A.; Burke, M. F. Sep.Sci. Technol. 1981, 16, 861884.
Analytical Chemistry, Vol. 72, No. 18, September 15, 2000 4353
tions. Since the analytical stochastic theory of nonlinear chromatography is still only a dream, the MC version is thus an appealing “surrogate”. In the present paper, a MC model of nonlinear chromatography based on discrete events simulation9 is presented, where nonlinear isotherm effects coming from competition of particles with respect to the sorption site are accounted for. A model of “virtual” MC chromatograph, where not only can the column type, operating conditions, sample type, size and composition, and injection mode be changed and chromatographic peak collected but the progress of the separation process can also be monitored. Moreover, the problem of thermodynamic correspondence of the proposed model with the classical Langmuir adsorption mode is discussed and modeling of a real case is performed. It must be pointed out that the computational approach here presented belongs to the theoretical field of pure chromatography. In fact, among the possible theoretical fields of chromatography, one can distinguish according to Giddings three basic levels of theoretical investigation in chromatography.10 A first level deals with specific structural concepts of physics and chemistry such as the phase distribution and equilibria, flow hydrodynamics, and kinetics of phase-transfer processes. In this level, the chromatography is not yet recognized. The second level deals with the pure chromatographic theory considering the dynamics of zone evolution. The third level contains the theory applied to laboratory and practical problem solving (e.g., programmed elution, optimization). MC methods are able to contribute to all these theoretical and applied aspects. The MC model developed here however belongs to level two; i.e., it deals with pure chromatographic theory.10,11 For this reason it is reported here as a “model” of chromatography. The main effort was to validate the MC computational program, and this is attained by comparing the “experimental” MC peak shapes to the exact peak shapes and peak parameters provided by the stochastic theory of linear chromatography10 and by other more sophisticated investigations. These results are, in fact, the prerequisite for future systematic application to the many open questions in separation science, not only in chromatography but also in closely related fields such as electromigration and field flow fractionation techniques. Thus, even though nonlinear effects are not quantitatively investigated, nor are other effectsssuch as multisite adsorption, longitudinal diffusion, flow pattern effects, nonimpulsive injection profile and multicomponent systems, structured stationary phases, etc.sexplicitly considered, the present model could account for them with only minor improvements. The aim of the present study is to contribute to the growing area of computational separation science,5 the most promising for further advances in separation methods and technologies. THEORY The General Stochastic Model of Linear and Nonlinear Chromatography. The basic chromatography model for MC simulation is the well-known stochastic model introduced by Giddings and Eyring,12,13 further developed by other authors,14-17 (9) Iazeolla, G. Introduzione alla simulazione discreta; Boringhieri: Torino, 1978. (10) Giddings, J. C. In Gas Chromatography; Goldup, A., Ed.; The Institute of Petroleum: London, 1965; p 3. (11) Dondi, F.; Cavazzini, A.; Remelli, M. In Advances in Chromatography; Brown, P. R., Grushka, E., Eds.; Marcel Dekker: New York, 1998; Vol. 38, pp 5174. (12) Giddings, J. C.; Eyring, H. J. Phys. Chem. 1955, 59, 416-421.
4354
Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
and recently reviewed.11 The chromatographic migration process is built up by two random processes: (a) entry into the stationary phase from the mobile phase; (b) desorption from the stationary phase. Process a controls process b, and the total time spent in the stationary phase is the sum of a random number of random intervals: it is a composite stochastic process.17 The case where different molecules do not interact with each other corresponds to linear chromatography, and the corresponding stochastic description will be called the “linear composite stochastic process”. For the sake of simplicity it is assumed that both the mobilephase velocity vm and the total time spent in the mobile phase tM are constant:
vm ) L/tM
(1)
where L is the column length. The general model of chromatography as a linear composite stochastic process has been mathematically handled with the “characteristic function” method under the sole hypothesis of randomness of the sorption site distribution along the column.4,10 The simplest case is the linear composite process of Poisson type in which both τm and τs, respectively, the “fly time” (i.e., the time between two subsequent sorption processes) and the desorption time from a single sorption site, are exponentially distributed random variables according to
f(τm) ) λm exp[-τmλm]
(2)
f(τs) ) λ∞s exp[-τsλ∞s ]
(3)
where λm and λ∞s are the entry and the desorption frequencies, respectively; the superscript ∞ denotes the infinite dilution condition. The mean fly time is
τjm ) 1/λm
(4)
whereas the mean desorption time is
jτ∞s ) 1/λ∞s
(5)
This case was first investigated by Giddings and Eyring12 and the resulting peak shape function is
f(t) ) exp[-λ∞s t - mm]xλ∞s mmtI1 (x4λ∞s mmt)
(6)
where I1 is the modified Bessel function of first kind and mm is the mean number of entries:
mm ) tMλm
(7)
The mean time spent in the stationary phase is (13) Giddings, J. C. J. Chem. Phys. 1957, 26, 169-173. (14) McQuarrie, D. A. J. Chem. Phys. 1963, 38, 437-445. (15) Dondi, F.; Remelli, M. J. Phys. Chem. 1986, 90, 1885-1891. (16) van Kampen, N. G. Physica 1979, 96A, 435-453. (17) van Kampen, N. G. Stochastic Processes in Physics and Chemistry; NorthHolland Physics Publishing: Amsterdam, 1981.
tS ) mm/λ∞s
(8)
MC simulation of the linear model of chromatography as a composite Poissonian stochastic process would simply consist of generating a chain of exponentially distributed fly times followed by exponentially distributed desorption times according to eqs 2 and 3.5,8 During fly time, the molecule advances with a fly velocity equal to vm and it enters the stationary phase having performed a random step equal to vm × τm. The single-molecule elution ends when it reaches the end of the column. It is possible to easily reproduce a process of entry into the stationary phase different from a Poisson entry or to account for mobile-phase diffusion by generating proper fly times and displacement steps.5 The sum of these random variables represents the residence time of one molecule in the chromatographic column. If each molecule behaves independentlyswhich is the basic assumption of the linear chromatographysa great enough number of single-molecule elutions can subsequently be generated one at a time: in this case, the column residence time histogram corresponds to the chromatogram. The case where molecules can interact in the sorption site is more complex and it must be handled differently: in this case, each sorption site must be continuously monitored since the molecule’s sorption behavior on that must account for the possible presence of other molecules. In practice, the chromatographic column must be considered as an array of labeled sorption sites and all the molecules, individually labeled, must be simultaneously displaced. The problem can be handled by using the same approach developed for simulating traffic problems or service stations, known as the “discrete event simulation technique”,9 which consists of processing a list of “events” ordered on the time axis (see Figure 1) and continuously updating “system status”. Appendix 1 reports the flow diagram for the computer program of MC model of nonlinear chromatography. Only two kinds of “events” are possible: (A) mobile phase w stationary phase transition (adsorption); (B) stationary phase w mobile phase transition (desorption). An “event” is thus specified by the following: (1) transition type (A or B); (2) molecule undergoing to transition; (3) site involved. Processing the event involves the following steps: (1) taking into consideration the next event; (2) generation of the corresponding random time (adsorption or fly time); (3) scheduling such event in the future at the proper time; (4) updating the “system status”. The “system status” comprises the following information: (1) molecule location (site/coordinate l along the column); (2) number of molecules in the different sites. It can be seen that this technique makes it possible to take nonlinearity into account. For example, if the event is the “mobile w stationary phase transition”, an adsorption time must be generated. Because the “system status” is known, we know at which sorption site the molecule is falling in adsorption and the coverage factor (θ) for that site at that time: thus, it is possible to generate a proper sorption time, a function of the local coverage. In this way, the nonlinearity is accounted for. For example, in the case of a sorption site having only monomolecular site capacity (CAP ) 1), the site coverage factor can be either 1 or 0 depending
Figure 1. Schematic representation of the “list of events” on the axis time. Event ep (for instance, an adsorption) takes place at the time tp (present) and it is processed. The corresponding random adsorption time τs is generated and a new desorption event is scheduled in the future at the proper time tp + τs. Since other n events (ep+1, ep+2, ..., deriving from previous adsorption-desorption processes) are already scheduled in the time interval [tp, tp + τs], such new desorption will become the ep+n event and its corresponding occurrence time (tp + τs) referred as tp+n.
on whether the site is free or already occupied. When θ ) 0 (i.e., if the site is empty), the random variable “sorption time” (τs) will be distributed according to eq 3; otherwise, the molecule will not stop in the site if it is occupied (i.e., τs is equal to zero for θ ) 1). Consequently, only when θ ) 0 is the particle behavior the same as in the Giddings-Eyring model (see eq 6). In the case of CAP > 1 (e.g., a pore where several molecules can be simultaneously trapped at a given time), the coverage factor spans between 0 and 1:
θ ) nads/CAP
(9)
where nads is the number of particles in the site. The mean time spent on the site is now a function of θ (denoted as τjθs ):
jτθs ) 1/λθs ) Φ(θ)τj∞s
(10)
and eq 3 can be rewritten as
f(τs) ) λθs exp[-τsλθs ]
CAP > 1
for
(11)
The function Φ(θ), here called “discipline factor”, modulates the mean time spent on the site according to the hypothesized model. For example, in the case when the Langmuir case is considered, the discipline factor can be expressed as (see Appendix 2)
Φ(θ) ) (1 - θ)
for for
θ 1 (parallel type) (13)
Chromatographic Quantities and Stochastic Quantities. The chromatographic peak is essentially a time frequency function, characterized by the mean (the retention time tR) and standard Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
4355
deviation σt. In addition, skewness18
Table 1. Standard Run Conditions 3
S ) κ3/σt
(14)
and excess
E ) κ4/σt4
(15)
are used to characterize peak tailing and flatness vs the standard Gaussian peak. κ3 and κ4 in eqs 14 and 15 are respectively the third and the fourth cumulants. The two basic chromatographic quantities are capacity factor k′ and plate height H:
k′ ) tS/tM
(16)
H ) L(σt/tR)2
(17)
tR ) tS + tM
(18)
with
k′ ) λm/λ∞s ) τ∞s /τm
Injected Sample sample type no. of injected molecules, nm Site capacity, CAP loading factor, L mean time spent in the site, τj∞s ) 1/λ∞s mean fly time, τjm ) 1/λm mm, mean number of sorption steps
L ) 2000 µm 3.5 µm × 0.5 µm 4000 carbon black benzene 1000 100 0.25% 10 µs 100 µs (jump over 40 sites) 100
Chromatographic Quantities mobile-phase velocity, vm 0.2 µm µs-1 retention, tS; tM; k′; R 1000 µs; 10000 µs; 0.1; 0.909 efficiency, N; H ) 2R(1 - R)vmτjm 6050; 0.33 µm Scale Factors with Respect to the Real Case; Benzene Elution on Capillary Column, 169 °C19 column length 1 × 10-4 slice/column sectiona 3 × 10-4 site capacity 1 × 10-4 a The slice/column section scale factor is the ratio between a considered slice area and the column cross-sectional area.
However, both k′ and H are defined and have meaning only under linearity conditions, i.e., under conditions of infinite dilution, and thus, both tR and σt must be evaluated under conditions of linearity. Combining eqs 7, 8, and 16 gives us
(19)
which is the sought after relationship between mean, single-step dynamic quantities (τjm and τj∞s ) and the result of the chromatographic experiment, k′. Appendix 2 also reports the relationship of the Henry constant in the case of gas solid chromatography with the above-reported dynamical quantities. EXPERIMENTAL AND COMPUTATION SECTION In this paper, the expressions “experimental data” and “simulation results” have the same meaning and they will be employed without distinction. Column Model and Standard Run Conditions. The complexity of a chromatographic column has been reduced here to a linear array of sorption sites. Assigning array length, number of sorption sites, the sorption time, and so on, means specifying column type, solute type, and temperature. The model should obviously correspond to a realistic case in order to obtain meaningful results from the MC experiment. Since in a real case the numbers of particles and sites are enormous, a proper scale reduction was applied, as extensively discussed in Appendix 3. The standard case considered is described in Table 1: it refers to a capillary column covered with an internal layer of carbon black.19 Carbon black particles are supposed to create a hole (the sorption site) where the sample molecules undergo the sorption process having the assigned mean time. (18) Crame´r, H. Mathematical Methods of Statistics, 13th ed.; Princeton University Press: Princeton, 1974. (19) Vidal-Madjar, C.; Bekassy, S.; Gonnord, M. F.; Arpino, P.; Guiochon, G. Anal. Chem. 1977, 49, 768-772
4356
Column Type column length single-site dimension no. of sites stationary-phase type
Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
For both the entry and desorption processes, exponential distributions (eqs 2 and 3) have been considered because the peak shape equation is known (eq 6) and thus the MC chromatography simulator can be validated for this model (linear or analytical conditions). In this case, the capacity factor is given by eq 16 and the other parameters are6,15
vm k′ H)2 2 ∞ (1 + k′) λs
(20)
σt ) x2mm/λ∞s
(21)
S ) 3/x2mm
(22)
E ) 6/mm
(23)
Since the mobile-phase standard velocity value is 0.2 µm µs-1 and the mean fly time is 100 µs, the mean fly path (τjm × vm) is 20 µm. The width of the site is 0.25 µm, and a distance of 0.25 µm spaces out the sites each other. A fly path length of 20 µm in the mobile phase means that an average of 1 out of 40 sites is visited by the molecule. Since the column contains 4000 sites (see Table 1), a molecule will undergo an average of 200 phase transition events () 2 × 4000/40). In the standard case of Table 1, 1000 molecules are handled in a single run and this is repeated at least 100 times in order to collect meaningful statistics, each standard run comprising ∼2 × 107 events and 100 000 molecules. It must be observed that 1000 molecules are enough to investigate the nonlinearity effect, which depends on the column loading. The percent loading factor L% defined as
L (%) )
number of injected molecules × 100 (24) number of sites × site capacity
expresses the column overloading. In the standard run, L ) 0.25%. This value ought to be enough to induce a detectable displacement of the peak mean with respect to the condition of infinite dilution, under Langmuir conditions. A rough estimation of this “isotherm effect” is given by the following equation:3
tmax,θ/tS,∞ ) (1 - x(L%/100))2
Table 2. Repeated Runs, nr ) 100, Standard Case of Table 1a
a
parameter
true value
Monte Carlo values
tS σ2 S E
1000 20000 0.212 0.060
997.6 ( 4.6b 19700 ( 9400 0.206 0.026
Site capacity, 100 (L ) 0.25%). b Confidence interval, 90%.
(25)
where tmax,θ and tS,∞ are, respectively, the time corresponding to the peak maximum under the condition of finite coverage and infinite dilution, in a very efficient column with zero band broadening. For L ) 0.25%, tmax,θ/tS,∞ ) 0.9. However, since the stochastic band broadening is also present and the peak maximum is behind the peak mean, eq 25 overestimates the peak-mean shift. Programming. The computer programsthe chromatography simulatorswas written in FORTRAN 77 and run both on mainframes (Cray X-MP 48, Casalecchio, Bologna, Italy; IBM 3090 University of Tennessee, TN) and RISC work stations, by using either the CERN or ISML libraries for computing mathematical functions. The RAND20 subroutine was employed to generate equally distributed (0-1) random numbers from which the exponential τm and τs times (eqs 2 and 3) were obtained by taking the logarithm. Separate RAND calls were organized for each individual molecule and for the different event types: introduction of the molecule into the column, fly time, and desorption time. In this way, the previously generated random time acted as a seed for the subsequent event of the same type involving the same molecule: the random series were thus disjoined, a condition needed for optimal statistical goodness. The starting series of seeds was moreover checked according to the uniformity χ2 test.18 These conditions were safe enough since in any case the maximum sequence length (107) was still lower than the starting sequence generation period (227). Event list management, with optimum access and ordering times, was a two-level list with access pointers, as suggested by Franta and Maly.21 In practice, the secondary list made it possible to rapidly scan a given time field, partitioned into a number of steps, to search/locate the proper event. Indexes were employed to establish the relationships among molecule type, event type, and their succession. The performance of this computer program was evaluated in terms of computational time and its dependence on the various parameters: the mean number of sorption-desorption steps (mm), the number of injected molecules (nm), and the number of (20) System/360 Scientific Subroutine Package (360A-CM-03X) Version III, IBM, 1966. (21) Franta, W. R.; Maly, K. Commun. ACM 20 1977, 8, 596-602. (22) van Deemter, J. J.; Zuiderweg, F. J.; Kninkenberg, A. Chem. Eng. Sci. 1956, 5, 271-289. (23) Feller, W. An Introduction to Probability Theory and Its Applications; J.Wiley and Sons: New York, 1966; Vols. I and II. (24) de Boer, J. H. The Dynamical Character of Adsorption; Clarendon Press: Oxford, 1968. (25) Dondi, F.; Gonnord, M.-F.; Guiochon, G. J. Colloid Interface Sci. 1977, 62, 303-315. (26) Fuller, E. N.; Schettler, P. D.; Giddings, J. C. Ind. Eng. Chem. 1966, 58, 19-27.
repeated runs (nr). The main improvement (4:1) was obtained by introducing a two-level indexed list instead of a simple linear list into the event list management, which originally occupied 80% of the total execution time. The CPU time resulted approximately proportional to the product (mm × nm × nr). RESULTS AND DISCUSSION Validation. MC simulator validation was performed under linear conditions and by comparing the experimental output data (histogram of the residence times and parameters) to the expected values (eq 6 for the histogram and eqs 14-17 and 20-23 for the parameters, respectively). The linearity conditions were assured by injecting a number of sample molecules not greater than the single-site capacity as a Dirac pulse (i.e., of zero width) and by using the “parallel” discipline factor of eq 13. Table 2 reports an example of the peak parameter validation. One can see that all the experimental parameterssincluding S and Esare close to those expected. Moreover, the confidence intervals of both tS and σ2 include the expected values. The 90% experimental percentile range of 4.6 of the retention time (see Table 2) is, in addition, close to this expected value:
1 2zR/2 n x mnr
∆tR ) σt
(26)
In this expression, the first factor is the standard deviation of the expected distribution (see eq 21), the second accounts for the number of the repeated observationssnumber of injected molecules (nm) times that of the repeated runs (nr)sand the third is the confidence interval of the standardized normal quantity z at R significance. Note that ∆tR and ∆tS are equal since a constant mobile-phase velocity (and thus a constant tM) is assumed. Figure 2 shows the excellent agreement between a series of experimental tR values, reported, and their experimental confidence intervals (full lines), as a function of 1/(nmnr)1/2, and the corresponding expected values. Validation of the global residence time histogram, i.e., of the chromatographic peak, is presented in Table 3. By using the χ2 test,18 100 repeated independent experimental histograms were compared (H0 hypothesis) both to the exact solution of the Giddings-Eyring function (eq 6, Bessel-type function) and to a Gaussian function of the same expected mean and standard deviation. The results of the H0 hypothesis test were classified at various degrees of acceptance (“accepted”, “ g suspicious”, “suspicious”, or “rejected”) according to the χ2 quantile range into which the χ2exp experimental value fell. One can see that the test is accepted Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
4357
Table 4. Study of the Site Capacity (Number of Runs, 100) site capacity loading factor L (%) stationary-phase retention time, tS,θ σt2 S E τjs expected tmax,θ/tS,∞ experimental tS,θ/tS,∞ H0 Bessel test H0 Gauss test a
100a 0.25 999.5 (1.1 19970 0.205 0.050 10.00 0.90 0.999 Y N
10 2.5 994.4 (1.2 19971 0.201 0.053 9.99 0.71 0.994 N N
1 25 519.2 (0.9 40749 1.076 0.790 9.99 0.25 0.520 N N
expected values linear case 1000 ()tS,∞) 20000 0.212 0.060 10
Standard case of Table 1.
Figure 2. Experimental retention time validation (see text). Table 3. Validation by χ2 Testa sample dimension nm ) 1000 nr ) 100 nm ) 100 000, nr ) 1
H0 hypothesisb
Bessel
Gauss
accepted gsuspicious suspicious rejected accepted
77/100 11/100 8/100 4/100 1/1
69/100 13/100 11/100 7/100 0/1
a Standard case of Table 1; L ) 0.25%; k′) 0.1; N ) 6050. b Accepted: P10 e χ2exp e P90. gSuspicious: P5 e χ2exp e P10 or P90 e χ2exp e P95. Suspicious: P1 e χ2exp e P5 or P95 e χ2exp e P99. Rejected: χ2exp e P1 or χ2exp g P99. Where Pn is n/100 percentile of the χ2 distribution at r degrees of freedom; χ2exp is the experimental χ2 value at r + 1 categories; and r + 1 ) 18.
in the majority of cases for both reference functions. However, close examination of the data shows that the results of the H0 hypothesis test distribution in the Bessel case fit the expected probability values (0.8, 0.1, 0.08, 0.02) a little bit more closely than does the Gaussian case comparison (see data of column 4 in Table 3). Moreover, a goodness-of-fit χ2 test performed by comparing each H0 test distribution to all the others18 gives positive results only for the Bessel case. A further validation is obtained by pooling all the 100 repeated runs into a single residence time distribution: only the H0 hypothesis with respect to the Bessel distribution proves accepted. One can thus conclude that the chromatography simulator is validated. Onset of Nonlinearity Effects. The onset of nonlinearity was investigated with reference to both site capacity and the amount injected and the discipline factor type. Table 4 reports the effect of increasing the loading factor, assuming a Langmuir-type discipline factor. The increase of L was obtained by keeping the number of injected molecules constant (1000) and by changing the site capacity. One can see that the Bessel test is no longer passed when L goes from 0.25 to g 2.5%. The observed experimental shift in retention time is also significant starting from that condition. The observed shifts are always lower than those expected from eq 25, as explained; however, the differences tend to be less significant for higher loading factors where the band dispersion effect stemming from the stochastic behavior loses relevance to the pure isotherm effect. The increase in S (from 0.205 to 1.076) looks particularly significant: it accounts for the strong tailing determined by the Langmuir-type isotherm effect. A more ex4358 Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
Figure 3. Experimental peak profiles at different column loading factors. The inset shows an expanded view of the lower curves, together with the proper Gaussian curve. Langmuir adsorption model. Standard run conditions, Table 1.
tended and quantitative evaluation of the isotherm effects would require a detailed, extensive work, since the MC results would have to be compared with those attainable by using other models and also taking into account the approximation degrees introduced in the respective model solutions. In Figure 3 the onset of nonlinearity is investigated by increasing the number of injected molecules (from 50 to 1000) and by keeping constant all other conditions including site capacity (CAP ) 1) and number of sites (15 000). Under such conditions, the peak shape deviates from the symmetrical Gaussian form and the peak maximum time shortens, with the tail becoming more pronounced and a vertical front appearing (the “shock”), as expected. Moreover, when the injected quantity is increased, the peaks build up the envelope typical of the Langmuir case.3 One can see that even at an L of 0.3% the peak deviates from that reported with L ) 0%, taken as reference. The latter was obtained by using a site capacity (CAP ) 2000) high enough to prevent any overloading effect; as expected, it is located at the right position of k′ ) 0.1. In Figure 4, a parallel-type discipline factor was instead applied (eq 13) and the loading factor was increased by lowering the CAP, while keeping all the other quantities constant. The same effect of peak distortion is evident, but it converges more quickly to a limit peak shape when the loading factor is decreased: when the site capacity is 8 (molecules), the resulting peak is almost Gaussian.
Figure 4. Experimental peak profiles at different site capacities. Parallel adsorption model (see eq 13). Standard run conditions, Table 1.
Figure 6. Convergence to the Gaussian shape as a function of column length (L) and mobile-phase velocity (vm). Circles: vm ) 0.2 µm µs-1; L ) 1000, 2000, 4000%. Squares: vm ) 0.3 µm µs-1; ) 2000, 3000%.
Figure 5. Convergence to the Gaussian shape: effect of the column length (µm). Mobile-phase velocity: vm ) 0.2 µm µs-1.
Exploitation of the MC Model. The MC model can be employed to investigate a large number of aspects of the chromatographic practice. Two examples are reported here. The effect of column length on the peak shape is shown in Figure 5, where the difference ∆Y between the “experimental” peak shape and the Gaussian peak shape (having the same mean and standard deviation) is reported in a standardized plot. One can see that the differences have the same shape as the first Gaussian-curve derivative and their level lowers when column length is increased.18 This effect is the content of the central limit theorem, which establishes that the Gaussian shape is reached during the process of the addition of random variables: increasing the column length increases the number of the sorptiondesorption steps and the peak shape comes closer and closer to the Gaussian shape. It must be observed that all the simulations were done under condition of loading factor lower than 0.1% in order to exclude the nonlinearity effect. Accordingly, when the column was lowered, the number of molecules was lowered to the same degree but the number of replicas was increased to keep statistical precision constant. The rate of convergence to the Gaussian law is described in a quantitative manner in Figure 6 reporting the total squared deviations F ) (∑(∆Y)2)1/2 versus L/νm, in a bilogarithmic plot. In this plot, (∑(∆Y)2)1/2 is the “distance” between the actual peak shape and the limit shape. One can see that this distance lowers when L/νm is increased, thus representing the functional conver-
Figure 7. Comparison between “experimental” van Deemter plots obtained through the Monte Carlo simulation (symbols) and the corresponding straight lines computed for the Giddings-Eyring model (eq 20), at different k′ values.
gence to Gaussian law.4,15,18 The L/νm variable is proportional to the mean number of sorption-desorption processes performed by the particle. The slope value of -0.448 in Figure 6 is close to the theoretical value of -0.5. This value is related to the content of the central limit theorem, which establishes that in a process of addition of n random variables the rate of convergence to the Gaussian law is proportional to 1/n1/2. Alternatively, the bilogarithmic plot transformation of this relationship turns out in a linear plot of slope equal to -0.5 (as shown in Figure 6).4,15 Figure 7 reports a van Deemter plot22 under low-dilution conditions, at different values of k′ within the vm range 0.1-0.5 µm µs-1. This plot accounts only for the C term of the van Deemter equation, since the stochastic model here considered only includes the sorption-desorption kinetics effect.3,6,7 One can see that, over an extended range of both k′ and vm values, all the points fall along the expected theoretical line of eq 20. CONCLUSIONS The MC model presented here has been validated by a series of statistical tests and by some examples based on the classical stochastic theory (e.g., the van Deemter plot) and on some Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
4359
advancements of probability theory applied to chromatography. Moreover, it has been shown how some problems connected with the classical central limit theorem, not yet fully elucidated even by the most advanced theory, can be dealt with using the present MC method. For example, the theorem establishes the proportionality between the rate of convergence to the Gaussian law and 1/n1/2 (and thus 1/(L/vm)1/2). The corresponding intercept value, establishing the actual rate of convergence, is instead not known but could be evaluated by the present approach (see Figure 6). Finally, it is worth underlining that, since all the chromatographic parameters depend on the sorption time distribution, as established by the Berry-Esseen theorem,15,23 once the pertinent sorption-desorption kinetics is assumed, the resulting chromatographic behavior can be fully investigated by the MC simulation approach. In the present study, exponential time distributions of both the entry and desorption processes were employed. In Appendix 3,6 an explanation is given of how the time constants were chosen. More advanced and complex sorption-desorption kinetics and the effect of column heterogeneity will be exploited in forthcoming papers. ACKNOWLEDGMENT The authors gratefully acknowledge the University of Tennessee, for computation resources, NATO (CRG Linkage Grant 971480) and Italian MURST, for financial support. APPENDIX 1 A flow diagram for the computer program of MC model of nonlinear chromatography is given in Chart 1. APPENDIX 2 Equivalence between the MC Model of Nonlinear Chromatography under the “Langmuir” Discipline Factor Condition and the Langmuir Adsorption Equilibrium. In the present treatment, reference to the gas adsorption on flat, homogeneous surfaces is made, the adsorbed and the gas phases being the stationary and the mobile phases, respectively. The basic equations reported by deBoer24
js ) n τj∞s
(linear conditions)
(II.1)
The nonlinearity comes from the fact that not all the n striking molecules are actually adsorbed on the surface, but only those impinging on the bare surface, which are n(1 - θ). Another equivalent point of view to interpret the nonlinearity effect is to consider that the effective mean sorption timesi.e., the time common to all striking moleculessis dependent on coverage according to
jτθs ) (1 - θ)τj∞s
(II.4)
jτθs ) Φ(θ)τj∞s
(II.5)
or, in general
which corresponds to eq 10 in the text. If τjm is the mean fly time, A the surface area, and N the number of moles, the n quantity in eqs II.1 and II.2 can be expressed as
n ) NNA/τjmA
where NA is the Avogadro number. If the gas is perfect, from eqs II.1 and II.6 one obtains
js )
pVmNA τj∞s ART τjm
(linear conditions)
js ) (1 - θ)n τj∞s
(nonlinear Langmuir conditions) (II.2)
relate the surface equilibrium concentration, js (number of molecules/surface unit), to two dynamical quantities n and τj∞s , respectively, the mean number of molecules striking the unit surface per unit of time and the mean time spent by a molecule remaining in the adsorption state, for the two cases of linear (eq II.1) and nonlinear Langmuir-type (eq II.2) adsorption. The superscript ∞ underlines the fact that τj∞s refers to the infinite dilution conditions and it is assumed to be independent of js, the effects of the lateral interaction among adsorbed molecules being neglected. θ is the surface coverage:
θ ) js /sj0
(II.7)
where Vm is the mobile (gas) phase volume, p partial pressure of the considered solute, R the gas constant and T the absolute temperature. By putting
K∞H )
Vm NA τj∞s A RT τjm
(II.8)
eq II.7 becomes
js ) K∞Hp
and
(II.6)
(II.9)
which is the well-known Henry law for gas adsorption. If the phase ratio
β ) Vm/A
(II.10)
is introduced, and eq 19 of the text is employed, one obtains the well-known equation relating the Henry constant to the capacity factor: ∞
k′ )
RT KH NA β
(II.11)
(II.3) The Langmuir equation is derived from eqs II.1, II.2, and II.9:
where js0 is the maximum surface concentration. Equation II.1 is the limit of eq II.2 for js , js0 or when θ f 0. 4360
Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
js ) s0K∞Hp/(sj0 + K∞Hp)
(II.12)
Chart 1. Flow Chart
The conventional thermodynamic description of Langmuir adsorption (eq II.2) is thus related to the mean times τj∞s and τjm spent by a single molecule respectively in the stationary and mobile phases (eq II.8). Moreover, eq II.5 tells us that in order to use the MC method to simulate any type of nonlinearity, the
proper Φ(θ) factor must be set as a multiplying factor for the mean sorption time τj∞s . APPENDIX 3 The Chromatographic Case Model. In building the model, reference is made to the real case of benzene elution on a capillary Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
4361
column (0.5 mm i.d. × 19 m), internally coated with a thin layer of graphitized thermal carbon black19 with an estimated thickness of 3.2 µm. This example was chosen since, in our opinion, it appeared the simplest to represent as a stochastic model. The adsorbing material is a conglomerate of 0.25-µm grains of graphite, a homogeneous, nonporous adsorbent. Its adsorption behavior has already been experimentally investigated both under linear and nonlinear conditions25 and theoretically interpreted by statistical thermodynamics. The single sorption desorption process can be represented as first-order kinetics with the mean adsorption time computed according to the Frenkel equation.24 For the benzene molecule, adsorption energy (Q) is about 6-7 kcal/mol and thus the mean sorption time is ∼10 ns. A simple reduction of a longitudinal cross section of this column was thus imagined as a sequence of full and void squared prisms 0.25 µm × 0.25 µm × 3.2 µm, the void prisms being the “sites” in which the solute molecule adsorbs. Inside the adsorption site, the mobile phase is “stagnant” and the molecule performs a series of elementary sorption processes, each of which can be assumed to follow the Frenkel equation. The global behavior in the adsorption site is thus complex. It was assumed that the mean number of single sorption processes is approximately the ratio between the time spent by a molecule to diffuse toward the bottom of the site, and the time required to diffuse laterally:
nsteps ≈ tbottom/tlateral
(III.1)
This quantity can easily be estimated using the Einstein equation valid for a two-dimensional diffusion:
d2 ) 4DMt
(III.2)
where d2 measures the mean distance covered by the molecule in time t and DM is the diffusion coefficient in the mobile phase. By taking into account the dimensions of the hole (3.2 × 0.25 µm)
nsteps ≈ (3.2/0.125)2 ≈ 1000
(III.3)
The effective mean sorption time in such sites should then be 10 ns × 1000 ) 10 µs. The global process of adsorption in the site ought to be the convolution of 1000 individual exponential processessa Γ distribution function23swhich, by virtue of the central limit theorem, is practically a Gaussian of mean equal to 10 µs and standard deviation of 0.31 µs. However, the number of sorption-desorption steps is not always the same and thus the resulting distribution is tailed. It was assumed that the sorption times at the site are exponentially distributed (see eq 3 in the text) with 1/λ∞s ) 0.1 µs-1. This assumption accounts for the fact that the real condition
4362
Analytical Chemistry, Vol. 72, No. 18, September 15, 2000
is geometrically much more complex than a simple reduction to the assumed sequence of void and full prisms and the resulting sorption process in the site is much more dispersed than the simple Gaussian function. Moreover, this hypothesis allows reference to the stochastic model of Giddings-Eyring, whose peak shape is known, thus enabling us to use eq 6 (see text) to validate the MC result. A mean mobile-phase fly time τjm value of 100 µs was assumed. This value is comparable to that derived from the diffusion Einstein equation (625 µs), assuming that the benzene molecule can reach the surface of the sorption layer by diffusion from the center of the column of radius 0.25 mm with a diffusion coefficient DM ) 0.5 cm2 s-1.26 The distribution function f(τm)was assumed of exponential type (see eq 2 of the text), again to make use of the Giddings-Eyring expression (eq 6) for peak shape. The model is obviously approximated. For example, it neglects the fact that the already desorbed molecule has a higher probability of being adsorbed than does a molecule which has remained in the mobile phase longer (i.e., that adsorption is a diffusion-controlled process). These aspects, together with a better description of the time spent in the sorption site lie beyond the aim of this initial approach. It is worth noting that, despite these approximations, the k′ value of 0.1 at v ) 0.2 µm µs-1 is not significantly different from those experimentally measured. Moreover the H value of 0.33 µm (under the standard conditions in Table 1) is roughlysalthough lower, possibly because of the neglected band-broadening factorssthat found in the reference work.19 The real case explained above is impossible to manage, both in terms of computation time and amount of memory required, because of the enormous number of both sorption sites (=1012) and individual injected molecules (=1015). Moreover, each sorption site has a very great capacity, expressed as the ratio of the area (= 3 × 108 Å2) with respect to the benzene surface area (=25 Å2), amounting thus to 107. Under these conditions, L was however 0.01%, the right condition for “analytical” separations where nonlinearity effects are to be avoided. In effect, the maximum retention time distortion is evaluated as no greater than 2% (eq 25 in the text). To make the problem tractable in terms of computer resources, preserving a certain meaning with respect to the reality and the content of the investigated problems, several scaling factor were applied to column length (10-4), column radius (3 × 10-4, i.e., roughly the ratio between dimension of the pore, 0.5 µm, and column circumference, 1500 µm) and site capacity (10-5-10-7). The number of the injected molecules were accordingly reduced by comparable factors. In this way, realistic values of the L factor of practical interest (0.1-10%) were exploited. Received for review March 21, 2000. Accepted June 28, 2000. AC0003347
