- Email: [email protected]
Numerical determination of competitive adsorption isotherm of mandelic acid enantiomers on cellulose-based chiral stationary phase
Journal of Chromatography A, 1202 (2008) 34–39
Contents lists available at ScienceDirect
Journal of Chromatography A journal homepage: www.elsevier.com/locate/chroma
Numerical determination of competitive adsorption isotherm of mandelic acid enantiomers on cellulose-based chiral stationary phase Yan Zhang, Sohrab Rohani, Ajay K. Ray ∗ Department of Chemical and Biochemical Engineering, The University of Western Ontario, London, Ontario N6A 5B9, Canada
a r t i c l e
i n f o
Article history: Received 10 December 2007 Received in revised form 6 June 2008 Accepted 16 June 2008 Available online 24 June 2008 Keywords: Inverse method Competitive adsorption isotherm Least-square fitting Genetic algorithm Orthogonal collocation on finite element
a b s t r a c t The use of inverse method for the determination of competitive adsorption isotherm of mandelic acid enantiomers on cellulose tris(3,5-diethylphenyl carbamate) stationary phase is proposed in this work. Non-dominated sorting genetic algorithm with jumping genes (NSGA-II-JG) was applied to acquire the isotherm parameters by minimizing the sum of square deviations of the model predictions from the measured elution profiles. Three different competitive isotherm models, i.e., Langmuir, biLangmuir and ´ Toth, combined with transport-dispersive chromatographic model were used in predicting the elution profiles. Orthogonal collocation on finite element (OCFE) method was applied to obtain the calculated ´ elution profiles. Results indicate that biLangmuir isotherm and Toth isotherm give remarkably similar equilibrium isotherms within the investigated liquid concentration range. Band profiles calculated from both isotherm models are in good agreement with the experimental data. The validity of the determined parameters was verified by comparing the model predictions with experimental elution profiles at various experimental conditions. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Preparative high performance liquid chromatography is currently the most widely used method for separation and purification of optical enantiomers [1,2]. Preparative chromatographic enantioseparations are performed essentially under nonlinear conditions and separations are governed by the equilibrium isotherm of the solutes in the feedstock [3]. The study of phase equilibriums and more specifically the determination of adsorption isotherm of the optical enantiomers are extremely important, either from a scientific or an industrial point of view. Frontal analysis [4,5] and perturbation peak [6–8] method are frequently used for determination of binary or multi-component isotherm. In frontal analysis, successive or abrupt step changes in the feed are introduced and the breakthrough curves are monitored from the beginning until final compositions are reached. For a system with N components, (N − 1) intermediate plateau concentrations and N retention times of shock fronts from each breakthrough curve have to be analyzed so that the equilibrium data can be calculated by the integral mass balance equations. Adsorption parameters are then determined by fitting these equi-
∗ Corresponding author. Tel.: +1 519 661 2111x81279; fax: +1 519 661 3498. E-mail address: [email protected] (A.K. Ray). 0021-9673/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.chroma.2008.06.026
librium data to the isotherm models using linear or nonlinear regression. Perturbation peak method determines the isotherm by measuring the retention times of small-size perturbations injected into the column equilibrated with sample solutions at different concentrations. A least-square fitting of the retention times of calculated perturbation peaks to those of measured ones is used to obtain the isotherm parameters. Although both methods can help to acquire very accurate isotherm parameters, they share a drawback of consuming a large amount of expensive pure enantiomers [5,9]. Inverse method [10,11] is another choice for isotherm determination of binary systems. This method allows the rapid estimation of the best values of the isotherm parameters by minimizing the differences between experimental overloaded profiles of binary mixtures and the profiles calculated by solving the mass balance equation for the system under examination. Due to the low consumption of enantiopure chemicals, inverse method has been applied in many preparative enantioseparations, particularly simulated moving bed separations, to determine nonlinear competitive equilibrium isotherms [12–15]. Numerical determination of the competitive isotherm parameters of recemic mandelic acid on cellulose tris(3,5-diethylphenyl carbamate) stationary phase (ChiralCel OD) is presented in this work. Cellulose tris(3,5-diethylphenyl carbamate) is the most powerful cellulose derivative chiral stationary phases (CSPs). It has been
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
successfully used to separate more than 20 racemates. However, equilibrium studies aimed at acquiring the adsorption isotherm on this kind of CSP are seldom reported [16]. The competitive adsorption isotherm of (R)- and (S)-mandelic acids, valuable chiral building blocks and very important auxiliaries and resolving agents in the fine-chemical and pharmaceutical industry, is studied in this work. The current study also aims to provide the prerequisite isotherm parameters for the simulation and optimization of a hybrid SMB and fractional crystallization process. 2. Theoretical 2.1. Isotherm models Several models have been suggested to describe the distribution of dilute sample solutes between the stationary and the mobile phase, and which one applies to a particular binary separation is hard to be determined from the experimental elution profiles. To evaluate the possible adsorption isotherm, three competitive ´ isotherm models [17], i.e., Langmuir, biLangmuir and Toth, were selected in this study to predict the equilibrium behaviors of the mandelic acid on cellulose tris(3,5-diethylphenyl carbamate) stationary phase. 2.1.1. The competitive Langmuir isotherm This model assumes monolayer coverage of solute molecules over a homogeneous adsorbent surface [18]. It is the most common model to describe nonlinear adsorption, although it might be too simple to account for the adsorption on chiral stationary phase. The competitive Langmuir isotherm for each component, i, is given as q∗i =
qs bi ci 1 + b1 c1 + b2 c2
(1)
where q∗i is adsorbed amount of each component in the stationary phase at equilibrium of the mobile phase concentration ci ; qs is the monolayer capacity of the CSP and bi is the equilibrium constant of adsorption; subscript 1 refers to the less retained component ((S)-mandelic acid) while subscript 2 the more retained component ((R)-mandelic acid). 2.1.2. The competitive biLangmuir isotherm This model assumes that the surface of chiral stationary phase is paved with two different types of binding sites, either one being homogeneous. The equilibrium isotherm is the addition of two independent Langmuir isotherms: q∗i =
q1s b1,i ci 1+
2
b c i=1 1,i i
+
q2s b2,i ci 1+
2
b c i=1 2,i i
(2)
where q1s and q2s are the saturation capacities of the two sites; b1,i and b2,i are the equilibrium constants of the two sites for component i. 2.1.3. The competitive T´ oth isotherm ´ isotherm [19] was developed to improve the fitting of the Toth Langmuir isotherm to experimental data based on the premise that there exists site heterogeneity on the sorbent and that most sites have sorption energy lower than the maximum adsorption energy. It differs from the Langmuir isotherm only by an exponent, (0 < < 1). q∗i =
qs bi ci 1/
[1 + (b1 c1 + b2 c2 ) ]
(3)
35
2.2. Column model Among the various mathematical models which are used to describe the migration of molecules through a chromatographic column, general rate model (GRM) is considered to be the most complete and realistic model of chromatography [3]. GRM takes into account all the phenomena that may have an influence on the band profile. However, requirement for the independent determination of too many parameters and long computation time limits its usage. The simplified lumped kinetic models, particularly the transport-dispersive (TD) model, are frequently used in liquid chromatography [20–22]. Mass balance equations of the TD model are given as ∂ci ∂2 c ∂q ∂c + ϕ i + u i = DL 2i ∂t ∂t ∂z ∂z
(4)
∂qi = km,i (q∗i − qi ) ∂t
(5)
where ci and qi are the mobile and stationary phase concentrations of each component i at time and space coordinates t and z; q* is the solid phase concentration in equilibrium with the liquid phase concentrations by isotherm models; ϕ is the phase ratio related to the total column porosity εt , by ϕ = (1 − εt )/εt , u is the interstitial mobile phase velocity; DL is the axial dispersion coefficient and km,i is the lumped rate coefficient. Initial and boundary conditions (Danckwerts conditions) are as follows: ci (0, z) = 0
qi (0, z) = 0 0 < t < tF t > tF
(7)
= u(ci (t, 0+ ) − ci (t, 0− ))
(8)
ci (t, 0− ) =
cF,i 0
(6)
∂c DL i ∂z
∂ci ∂z
z=0
=0
(9)
z=L
where cF,i is the sample feed concentration and tF is the injection time, L is the column length. The assumption that the sample is introduced into the column as a rectangular pulse of length tF may not be valid in most practical applications. However, as the injection time is very small compared to the retention time, such a simplification still seems to be applicable [13,23]. The effect of the simplified injection profile on the derived isotherm parameters will be discussed in Section 4.4. 2.3. OCFE for calculation of the elution profiles Orthogonal collocation on finite element (OCFE) was used in this work to obtain the simulated band profiles by solving Eqs. (4)–(9) simultaneously. In the OCFE method [24], the spatial domain was divided into “NE ” elements, then “NP ” internal collocation points are specified in each element. The internal collocation points are chosen as the roots of an NP th degree orthogonal polynomial. The first and second order differentials with respect to space are then approximated using two Matrices A and B which are obtained by solving the Gaussian–Jacobi quadratures. In order to apply OCFE method, dimensionless parameters must be used in spatial domain. Eqs. (4) and (5) thus transform into the following ODEs at internal collocation point j (j = 2 to NP + 1) in element l (l = 1 to NE ).
36
∂c ∂t
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
k=NP +2 c=cjl
= −ϕkm (q∗jl − qjl ) −
k=1
Ajk ckl +
Pe
k=NP +2
Table 1 Axial dispersion coefficient (DL ) and lumped rate coefficient (km ) used in calculation
Bjk ckl ,
k=1
k = 1, 2, . . . , NP + 2
Q (ml/min)
DL (×10−5 cm2 /s)
km (s−1 )
(10)
0.7 0.9
4.17 5.5
1.4 1.52
(11)
(R)-mandelic acid (98%) and (S)-mandelic acid (99%) were obtained from Alfa Aesar (Heysham, Lancashire). Solutions of solid compounds for the determination of the isotherm data were prepared by weighing the required amounts of solutes in the mobile phase.
∂q ϕ ∂t
q=qjl
= km (q∗jl − qjl )
where ␥ = u/LE , LE = L/NE ; Pe = uLE /DL . The continuity of the first derivative between the various elements gives us (NE − 1) algebraic equations as follows:
k=NP +2
k=NP +2
ANP +2,k ck,l−1 =
k=1
A1,k ck,l ,
l = 2, 3, . . . ., NE
(12)
k=1
Boundary conditions at the inlet and outlet of the column (Eqs. (8) and (9)) provide the last two algebraic equations,
k=NP +2
A1,k ck,1 = Pe (c1,1 − c(t, 0− ))
(13)
ANP +2,k ck,NE = 0
(14)
k=1 k=NP +2
k=1
The resulting system of stiff algebraic-differential equations (a total number of NE (2NP + 1) + 1 equations for each component) was solved using routine DASPG in the IMSL based on Petzold–Gear BDF method. In this study, 50 elements for NE and two internal collocation points for NP were selected in all the calculations. 2.4. Inverse method Inverse method involves in tuning the coefficients of the preselected adsorption isotherm model so that the difference between one or several experimental band profiles and the profiles calculated using a model of chromatography is minimum. An error function F(p), defined as the sum of square deviations of the calculated concentration profiles from the experimentally measured curves, is used as the objective function to obtain the best-fit values of isotherm parameters. F(p) = min
p m
exp
mod [ck,j − ck,j ]
2
(15)
k=1 j=1 exp
mod are the measured and model predicted concenwhere ck,j and ck,j trations of the racemic mandelic acid at point j of the elution curves of the kth injection. F(p) value can be used as a kind of degree of overlap between the experimental and simulated band profiles. By minimizing the F(p) value, isotherm parameters which lead to the highest degree of overlap between the experimental and simulated band profiles are determined. In this study, a new adaptation of the state-of-the-art optimization method, non-dominated sorting genetic algorithm with jumping genes (NSGA-II-JG) [25] was used to obtain the best-fit values of the isotherm parameters. Details of the NSGA can be found elsewhere [26].
3. Experimental 3.1. Materials The mobile phase was a mixture of hexane and isopropyl alcohol (85:15, v/v), with 0.1% of acetic acid used as a modifier. Hexane and IPA (both HPLC grade) as well as acetic acid (ACS grade) were purchased from EMD Chemicals (Gibbstown, NJ, USA). Enantiopure
3.2. Apparatus and methods All the overloaded band profiles were obtained using Waters Breeze HPLC system (Waters Corporation, Milford, MA, USA). This system consists of a binary HPLC pump, an in-line degasser, an autosampler with a 100-l sample loop, a dual wavelength absorbance detector and a data workstation. UV absorbance was measured at a wavelength of 254 nm. The detector was calibrated before experiments and the responses from the detectors were linearly proportional to the concentrations. A ChiralCel OD-H column (250 mm × 4.6 mm, 5 m), purchased from Chiral Technologies (West Chester, PA, USA), was used for all the measurement. This column was packed with cellulose tris(3,5diethylphenyl carbamate) coated on 5-m silica-gel substrate. The column hold-up volume, Vm , was evaluated from the record of the detector set to 254 nm after injection of pure IPA into mixed mobile phase. The total column porosity (εt = 0.74) corresponding to the hold-up volume was determined from εt = Vm /V, where V is the geometric volume of the column. The extra-column volume between the autosampler to the inlet of column (Vex ≈ 180.0 l) was measured by removing the column from the system. All the retention data were corrected for this contribution. The overloaded band profiles of racemic mixture were used to derive the competitive isotherm parameters. Racemic feed solutions were obtained by dissolving the same amount of (R)and (S)-mandelic acids into the mobile phase. Experiments were performed at 25 ± 0.1 ◦ C with the elution flow rate of 0.7 and 0.9 ml/min. The injection volumes were 30.0 and 45.0 l and the total sample concentrations were 20.0, 30.0 and 40.0 g/l. 4. Results and discussion 4.1. Determination of the kinetic parameters Accurate determination of the two kinetic parameters, DL and km,i is of great importance for obtaining the correct isotherm parameters. In this study axial dispersion coefficient, DL, was determined from theoretical correlation [18]; lumped rate coefficient, km,i was determined along with the isotherm parameters by fitting the model predicted band profiles to the experimental data. For simplicity, same km value was used in the calculation for both components. DL and km at various flow rates are listed in Table 1. 4.2. Best-fit isotherm parameters Four overloaded elution band profiles of the racemic mixture with different injection volumes and feed concentrations were used to derive the isotherm parameters of different models. Each band profile used for fitting contains 560 points in the nonzero concentration range. Initially, some simulations were carried out to determine the possible ranges of each decision variable which were prerequisite in the optimization. Systematic optimization was then performed to determine the isotherm parameters by fitting the
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
37
Table 2 Competitive isotherm parameters by inverse method Langmuir qs (g/l) b1 (× 10−2 l/g) b2 (× 10−2 l/g) F(p)
62.46 5.6 7.5 48.9
biLangmuir q1s (g/l) b11 (× 10−2 l/g) b12 (× 10−2 l/g) q2s (g/l) b21 (l/g) b22 (l/g) F(p)
89.8 2.92 3.82 1.65 0.94 1.32 36.65
´ Toth qs (g/l) b1 (× 10−2 l/g) b2 (× 10−2 l/g) F(p)
88.56 4.4 5.84 0.702 37.39
model predicted band profiles to the four chromatograms simultaneously. Results of the isotherm parameters as well as the F(p) values for the three isotherm models are summarized in Table 2. Comparison of the three isotherm models is illustrated in Fig. 1. ´ isotherm Results from Table 2 indicate that biLangmuir and Toth ´ modfit better than Langmuir isotherm while biLangmuir and Toth els give rather similar representations of the experimental data. The second conclusion coincides with the observation from Fig. 1, ´ in which biLangmuir and Toth models generate almost the same equilibrium isotherm with the given mobile phase concentration. Such observation confirms the heterogeneous nature of this kind of chiral stationary phase. Fig. 2 demonstrates the results of fitting biLangmuir model to three of the measured chromatograms. It can be seen that the fit is less satisfactory at small sample size, while bands with large sample size fit better. The main reason of the error made is that the influence of high concentration data on the values derived for the numerical coefficients is higher than that of the low concentration points. However, this kind of error is somewhat acceptable, since inverse method provides more reliable results at higher column loading.
Fig. 1. Comparison of the three competitive isotherm models for (R)- and (S)mandelic acids.
Fig. 2. Best-fit of competitive biLangmuir model.
4.3. Validation of the acquired parameters In order to verify the validity of the derived isotherm parameters, experimental elution band profile with higher injection volume (Vinj = 60.0 l, cF,t = 40.0 g/l) was compared with the profile calculated from the derived isotherm parameters listed in Table 2. Good agreement between the calculated and measured band profiles are also acquired this time. The comparison of the calculated (with biLangmuir isotherm) and measured profiles is shown in Fig. 3. Experiments using various solute compositions were conducted as well to see whether the acquired isotherm parameters could well predict the elution band profiles at different operating conditions. Fig. 4 illustrates the comparison of the calculated and measured elution profiles at various solute compositions. Results from Fig. 4 indicate that isotherm parameters acquired from the inverse method can give rather good predictions of the elution band profiles at a wide range of solute compositions. Both the simulated concentration profiles (using the isotherm parameters of competitive biLangmuir model) and the measured elution profiles at higher elution flow rate (Q = 0.9 ml/min) are plotted in Fig. 5. Remarkably good agreements are observed for all the four chromatograms with different solute compositions. Except for the slight deviation of the tail of the second component, numerical solutions from the chromatography model with the acquired isotherm parameters can provide accurate
38
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
Fig. 6. Measured inlet concentration profile. Fig. 3. Comparison of model predictions (biLangmuir) to experimental band profiles at higher injection volume.
between the measured and predicted concentration profiles are observed in Figs. 2–5. Errors from the following four aspects may lead to such discrepancies. First, the simplified boundary condition, which assumes an ideal rectangular pulse injection, may result in the deviations. Measurement errors could be the second cause. It is known that the experimentally measured concentrations were derived from the original chromatograms (absorbance versus time) through a calibration. Both the drift of baseline and the noise of the signals will bring about some errors in the measured concentration profiles. The third reason could be the selected isotherm models. Because of the highly complex structure of the selector, chromatographic mechanism involved in the chiral discrimination has not yet been fully clarified. The selected isotherm models may lead to some theoretical errors in predicting the elution curves. Finally, method used to solve the PDEs may also generate some numerical errors in the simulated band profiles. Compared with the results obtained using finite difference method, calculated band profiles from OCFE ´ model. are not so smooth when smaller is used for Toth Fig. 4. Comparison of the model predicted concentration profiles (biLangmuir) with measured elution profiles at different solute compositions.
predictions on the elution times and peak shapes of the two components. The achievement of the good agreement between the calculated and measured elution profiles under various experimental conditions indicates that isotherm parameters obtained from the inverse method are quite robust and reliable. However, small deviations
Fig. 5. Comparison of the model predicted concentration profiles (biLangmuir) with the measured elution profiles at higher elution flow rate and different solute compositions.
4.4. Impact of the boundary condition The inlet profile constitutes the boundary condition of the differential mass balance equation, thus it should have some influence on the simulated band profiles [13]. Usually, the shape of the inlet profile was assumed to be a rectangle for the sake of simplicity. However, significant deviations from this assumption are always observed. Fig. 6 illustrates a measured inlet profile of mandelic acid, determined by injecting the sample with column removed from the system and replaced by a zerovolume connector. Obvious asymmetry of the injection profile can be seen from Fig. 6. How this asymmetry affects the elution
Fig. 7. Calculated elution profiles using real and simplified boundary conditions.
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
profiles and the final derived isotherm parameters remains a problem. In order to quantify the effect of the simplified boundary conditions, optimization runs were also performed with experimentally measured injection profile as boundary conditions. This time the comparison is confined to the two chromatograms with injection amounts of 0.9 and 1.8 mg. Results of the isotherm parameters using the true inlet profile are very close to the values listed in Table 2. Deviations of the final isotherm parameters using the true inlet profile from those using the simplified rectangular injection profiles are less than 2%. This demonstrates that the simplified boundary condition do not have a significant impact on the final results of the isotherm parameters and such a simplification can be used in calculation. This observation coincides with that reported recently by [15]. From our simulation results, isotherm parameters acquired with simplified boundary conditions are still valid in predicting the elution profiles when real injection profile is applied. Comparison of the calculated elution profiles using both real and simplified boundary conditions with isotherm parameters listed in Table 2 is illustrated in Fig. 7. 5. Conclusions Inverse method was used in this study to determine the competitive isotherm of (R)- and (S)-mandelic acids on a kind of cellulose chiral stationary phase using NSGA-II-JG as the optimization method. By fitting four chromatograms simultaneously with different injection volumes and different feed concentrations, the best-fit isotherm parameters for three different competitive ´ isotherm models, i.e., Langmuir, biLangmuir and Toth, were deter´ models both give mined. Results indicate that Bilangmuir and Toth good representations of the experimental data. Accuracy and reliability of the derived isotherm parameters were verified by comparing the model predicted band profiles with the experimental elution profiles under various experimental conditions, such as different injection volume, different solute compositions and different elution flow rates. Results indicate that the band profiles calculated using the derived isotherm parameters can provide rather accurate predictions for all experimental data. Isotherm parameters derived from the inverse method are reliable.
39
The inverse method offers a rather quick method of adsorption isotherm determination of racemic mandelic acid with minimal sample and solvent consumption. However, one important disadvantage of inverse method is that the acquired data depend on the column efficiency, which cannot be accessed conveniently under nonlinear conditions. Therefore, for more accurate isotherm, frontal analysis is recommended. It should also be stressed that the isotherm parameters acquired in this study are satisfactory only for dilute solutions of analytes. References [1] G. Subramanian, Chiral Separation Techniques—A Practical Approach, Wiley–VCH, Weinheim, 2001. [2] S. Ahuja, in: S. Ahuja (Ed.), Chiral Separations—Applications and Technology, American Chemical Society, Washington, DC, 1997, p. 1. [3] G. Guiochon, A. Felinger, D.G. Shirazi, A.M. Katti, Fundamentals of Preparative and Nonlinear Chromatography, Elsevier Academic Press, San Diego, CA, 2006. [4] F. Gritti, G. Guiochon, J. Chromatogr. A 1028 (2004) 197. [5] A. Seidel-Morgenstern, J. Chromatogr. A 1037 (2004) 255. ¨ [6] C. Blumel, P. Hugo, A. Seidel-Morgenstern, J. Chromatogr. A 865 (1999) 51. ´ T. Fornstedt, Anal. Chem. 76 (2004) 4856. [7] J. Lindholm, P. Forssen, ´ T. Fornstedt, Anal. Chem. 76 (2004) 5472. [8] J. Lindholm, P. Forssen, [9] A. Cavazzini, A. Felinger, G. Guiochon, J. Chromatogr. A 1012 (2003) 139. ´ [10] F. James, M. Sepulveda, Inv. Prob. 10 (1994) 1299. [11] A. Felinger, A. Cavazzini, G. Guiochon, J. Chromatogr. A 986 (2003) 207. [12] D. Antos, W. Piatkowski, K. Kaczmarski, J. Chromatogr. A 874 (2000) 1. [13] A. Felinger, D. Zhou, G. Guiochon, J. Chromatogr. A 1005 (2003) 35. [14] N. Marchetti, F. Dondi, A. Felinger, R. Guerrini, S. Salvadori, A. Cavazzini, J. Chromatogr. A 1079 (2005) 162. ´ T. Fornstedt, J. Chromatogr. A 1099 (2005) 167. [15] R. Arnell, P. Forssen, ¨ [16] C. Heuer, E. Kusters, T. Plattner, A. Seidel-Morgenstern, J. Chromatogr. A 827 (1998) 175. [17] G. Guiochon, B. Lin, Modeling for Preparative Chromatography, Academic Press, San Diego, CA, 2003. [18] D.M. Ruthven, Principle of Adsorption and Adsorption Process, Wiley, New York, 1984. [19] J. Toth, Acta Chem. Hung. 32 (1962) 31. [20] D.C.S. Azevedo, L.S. Paris, A.E. Rodrigues, J. Chromatogr. A 865 (1999) 187. [21] A. Felinger, G. Guochon, Chromatographia 60 (2004) S175. [22] K. Kaczmarski, G. Guiochon, Anal. Chem. 79 (2007) 4648. ˜ [23] I. Quinones, C.M. Grill, L. Miller, G. Guiochon, J. Chromatogr. A 867 (2000) 1. [24] B.A. Finlayson, Nonlinear Analysis in Chemical Engineering, McGraw-Hill, New York, NY, 1980. [25] R.B. Kasat, S.K. Gupta, Comp. Chem. Eng. 27 (2003) (1785). [26] V. Bhaskar, S.K. Gupta, A.K. Ray, Rev. Chem. Eng. 16 (2000) 1.
Contents lists available at ScienceDirect
Journal of Chromatography A journal homepage: www.elsevier.com/locate/chroma
Numerical determination of competitive adsorption isotherm of mandelic acid enantiomers on cellulose-based chiral stationary phase Yan Zhang, Sohrab Rohani, Ajay K. Ray ∗ Department of Chemical and Biochemical Engineering, The University of Western Ontario, London, Ontario N6A 5B9, Canada
a r t i c l e
i n f o
Article history: Received 10 December 2007 Received in revised form 6 June 2008 Accepted 16 June 2008 Available online 24 June 2008 Keywords: Inverse method Competitive adsorption isotherm Least-square fitting Genetic algorithm Orthogonal collocation on finite element
a b s t r a c t The use of inverse method for the determination of competitive adsorption isotherm of mandelic acid enantiomers on cellulose tris(3,5-diethylphenyl carbamate) stationary phase is proposed in this work. Non-dominated sorting genetic algorithm with jumping genes (NSGA-II-JG) was applied to acquire the isotherm parameters by minimizing the sum of square deviations of the model predictions from the measured elution profiles. Three different competitive isotherm models, i.e., Langmuir, biLangmuir and ´ Toth, combined with transport-dispersive chromatographic model were used in predicting the elution profiles. Orthogonal collocation on finite element (OCFE) method was applied to obtain the calculated ´ elution profiles. Results indicate that biLangmuir isotherm and Toth isotherm give remarkably similar equilibrium isotherms within the investigated liquid concentration range. Band profiles calculated from both isotherm models are in good agreement with the experimental data. The validity of the determined parameters was verified by comparing the model predictions with experimental elution profiles at various experimental conditions. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Preparative high performance liquid chromatography is currently the most widely used method for separation and purification of optical enantiomers [1,2]. Preparative chromatographic enantioseparations are performed essentially under nonlinear conditions and separations are governed by the equilibrium isotherm of the solutes in the feedstock [3]. The study of phase equilibriums and more specifically the determination of adsorption isotherm of the optical enantiomers are extremely important, either from a scientific or an industrial point of view. Frontal analysis [4,5] and perturbation peak [6–8] method are frequently used for determination of binary or multi-component isotherm. In frontal analysis, successive or abrupt step changes in the feed are introduced and the breakthrough curves are monitored from the beginning until final compositions are reached. For a system with N components, (N − 1) intermediate plateau concentrations and N retention times of shock fronts from each breakthrough curve have to be analyzed so that the equilibrium data can be calculated by the integral mass balance equations. Adsorption parameters are then determined by fitting these equi-
∗ Corresponding author. Tel.: +1 519 661 2111x81279; fax: +1 519 661 3498. E-mail address: [email protected] (A.K. Ray). 0021-9673/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.chroma.2008.06.026
librium data to the isotherm models using linear or nonlinear regression. Perturbation peak method determines the isotherm by measuring the retention times of small-size perturbations injected into the column equilibrated with sample solutions at different concentrations. A least-square fitting of the retention times of calculated perturbation peaks to those of measured ones is used to obtain the isotherm parameters. Although both methods can help to acquire very accurate isotherm parameters, they share a drawback of consuming a large amount of expensive pure enantiomers [5,9]. Inverse method [10,11] is another choice for isotherm determination of binary systems. This method allows the rapid estimation of the best values of the isotherm parameters by minimizing the differences between experimental overloaded profiles of binary mixtures and the profiles calculated by solving the mass balance equation for the system under examination. Due to the low consumption of enantiopure chemicals, inverse method has been applied in many preparative enantioseparations, particularly simulated moving bed separations, to determine nonlinear competitive equilibrium isotherms [12–15]. Numerical determination of the competitive isotherm parameters of recemic mandelic acid on cellulose tris(3,5-diethylphenyl carbamate) stationary phase (ChiralCel OD) is presented in this work. Cellulose tris(3,5-diethylphenyl carbamate) is the most powerful cellulose derivative chiral stationary phases (CSPs). It has been
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
successfully used to separate more than 20 racemates. However, equilibrium studies aimed at acquiring the adsorption isotherm on this kind of CSP are seldom reported [16]. The competitive adsorption isotherm of (R)- and (S)-mandelic acids, valuable chiral building blocks and very important auxiliaries and resolving agents in the fine-chemical and pharmaceutical industry, is studied in this work. The current study also aims to provide the prerequisite isotherm parameters for the simulation and optimization of a hybrid SMB and fractional crystallization process. 2. Theoretical 2.1. Isotherm models Several models have been suggested to describe the distribution of dilute sample solutes between the stationary and the mobile phase, and which one applies to a particular binary separation is hard to be determined from the experimental elution profiles. To evaluate the possible adsorption isotherm, three competitive ´ isotherm models [17], i.e., Langmuir, biLangmuir and Toth, were selected in this study to predict the equilibrium behaviors of the mandelic acid on cellulose tris(3,5-diethylphenyl carbamate) stationary phase. 2.1.1. The competitive Langmuir isotherm This model assumes monolayer coverage of solute molecules over a homogeneous adsorbent surface [18]. It is the most common model to describe nonlinear adsorption, although it might be too simple to account for the adsorption on chiral stationary phase. The competitive Langmuir isotherm for each component, i, is given as q∗i =
qs bi ci 1 + b1 c1 + b2 c2
(1)
where q∗i is adsorbed amount of each component in the stationary phase at equilibrium of the mobile phase concentration ci ; qs is the monolayer capacity of the CSP and bi is the equilibrium constant of adsorption; subscript 1 refers to the less retained component ((S)-mandelic acid) while subscript 2 the more retained component ((R)-mandelic acid). 2.1.2. The competitive biLangmuir isotherm This model assumes that the surface of chiral stationary phase is paved with two different types of binding sites, either one being homogeneous. The equilibrium isotherm is the addition of two independent Langmuir isotherms: q∗i =
q1s b1,i ci 1+
2
b c i=1 1,i i
+
q2s b2,i ci 1+
2
b c i=1 2,i i
(2)
where q1s and q2s are the saturation capacities of the two sites; b1,i and b2,i are the equilibrium constants of the two sites for component i. 2.1.3. The competitive T´ oth isotherm ´ isotherm [19] was developed to improve the fitting of the Toth Langmuir isotherm to experimental data based on the premise that there exists site heterogeneity on the sorbent and that most sites have sorption energy lower than the maximum adsorption energy. It differs from the Langmuir isotherm only by an exponent, (0 < < 1). q∗i =
qs bi ci 1/
[1 + (b1 c1 + b2 c2 ) ]
(3)
35
2.2. Column model Among the various mathematical models which are used to describe the migration of molecules through a chromatographic column, general rate model (GRM) is considered to be the most complete and realistic model of chromatography [3]. GRM takes into account all the phenomena that may have an influence on the band profile. However, requirement for the independent determination of too many parameters and long computation time limits its usage. The simplified lumped kinetic models, particularly the transport-dispersive (TD) model, are frequently used in liquid chromatography [20–22]. Mass balance equations of the TD model are given as ∂ci ∂2 c ∂q ∂c + ϕ i + u i = DL 2i ∂t ∂t ∂z ∂z
(4)
∂qi = km,i (q∗i − qi ) ∂t
(5)
where ci and qi are the mobile and stationary phase concentrations of each component i at time and space coordinates t and z; q* is the solid phase concentration in equilibrium with the liquid phase concentrations by isotherm models; ϕ is the phase ratio related to the total column porosity εt , by ϕ = (1 − εt )/εt , u is the interstitial mobile phase velocity; DL is the axial dispersion coefficient and km,i is the lumped rate coefficient. Initial and boundary conditions (Danckwerts conditions) are as follows: ci (0, z) = 0
qi (0, z) = 0 0 < t < tF t > tF
(7)
= u(ci (t, 0+ ) − ci (t, 0− ))
(8)
ci (t, 0− ) =
cF,i 0
(6)
∂c DL i ∂z
∂ci ∂z
z=0
=0
(9)
z=L
where cF,i is the sample feed concentration and tF is the injection time, L is the column length. The assumption that the sample is introduced into the column as a rectangular pulse of length tF may not be valid in most practical applications. However, as the injection time is very small compared to the retention time, such a simplification still seems to be applicable [13,23]. The effect of the simplified injection profile on the derived isotherm parameters will be discussed in Section 4.4. 2.3. OCFE for calculation of the elution profiles Orthogonal collocation on finite element (OCFE) was used in this work to obtain the simulated band profiles by solving Eqs. (4)–(9) simultaneously. In the OCFE method [24], the spatial domain was divided into “NE ” elements, then “NP ” internal collocation points are specified in each element. The internal collocation points are chosen as the roots of an NP th degree orthogonal polynomial. The first and second order differentials with respect to space are then approximated using two Matrices A and B which are obtained by solving the Gaussian–Jacobi quadratures. In order to apply OCFE method, dimensionless parameters must be used in spatial domain. Eqs. (4) and (5) thus transform into the following ODEs at internal collocation point j (j = 2 to NP + 1) in element l (l = 1 to NE ).
36
∂c ∂t
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
k=NP +2 c=cjl
= −ϕkm (q∗jl − qjl ) −
k=1
Ajk ckl +
Pe
k=NP +2
Table 1 Axial dispersion coefficient (DL ) and lumped rate coefficient (km ) used in calculation
Bjk ckl ,
k=1
k = 1, 2, . . . , NP + 2
Q (ml/min)
DL (×10−5 cm2 /s)
km (s−1 )
(10)
0.7 0.9
4.17 5.5
1.4 1.52
(11)
(R)-mandelic acid (98%) and (S)-mandelic acid (99%) were obtained from Alfa Aesar (Heysham, Lancashire). Solutions of solid compounds for the determination of the isotherm data were prepared by weighing the required amounts of solutes in the mobile phase.
∂q ϕ ∂t
q=qjl
= km (q∗jl − qjl )
where ␥ = u/LE , LE = L/NE ; Pe = uLE /DL . The continuity of the first derivative between the various elements gives us (NE − 1) algebraic equations as follows:
k=NP +2
k=NP +2
ANP +2,k ck,l−1 =
k=1
A1,k ck,l ,
l = 2, 3, . . . ., NE
(12)
k=1
Boundary conditions at the inlet and outlet of the column (Eqs. (8) and (9)) provide the last two algebraic equations,
k=NP +2
A1,k ck,1 = Pe (c1,1 − c(t, 0− ))
(13)
ANP +2,k ck,NE = 0
(14)
k=1 k=NP +2
k=1
The resulting system of stiff algebraic-differential equations (a total number of NE (2NP + 1) + 1 equations for each component) was solved using routine DASPG in the IMSL based on Petzold–Gear BDF method. In this study, 50 elements for NE and two internal collocation points for NP were selected in all the calculations. 2.4. Inverse method Inverse method involves in tuning the coefficients of the preselected adsorption isotherm model so that the difference between one or several experimental band profiles and the profiles calculated using a model of chromatography is minimum. An error function F(p), defined as the sum of square deviations of the calculated concentration profiles from the experimentally measured curves, is used as the objective function to obtain the best-fit values of isotherm parameters. F(p) = min
p m
exp
mod [ck,j − ck,j ]
2
(15)
k=1 j=1 exp
mod are the measured and model predicted concenwhere ck,j and ck,j trations of the racemic mandelic acid at point j of the elution curves of the kth injection. F(p) value can be used as a kind of degree of overlap between the experimental and simulated band profiles. By minimizing the F(p) value, isotherm parameters which lead to the highest degree of overlap between the experimental and simulated band profiles are determined. In this study, a new adaptation of the state-of-the-art optimization method, non-dominated sorting genetic algorithm with jumping genes (NSGA-II-JG) [25] was used to obtain the best-fit values of the isotherm parameters. Details of the NSGA can be found elsewhere [26].
3. Experimental 3.1. Materials The mobile phase was a mixture of hexane and isopropyl alcohol (85:15, v/v), with 0.1% of acetic acid used as a modifier. Hexane and IPA (both HPLC grade) as well as acetic acid (ACS grade) were purchased from EMD Chemicals (Gibbstown, NJ, USA). Enantiopure
3.2. Apparatus and methods All the overloaded band profiles were obtained using Waters Breeze HPLC system (Waters Corporation, Milford, MA, USA). This system consists of a binary HPLC pump, an in-line degasser, an autosampler with a 100-l sample loop, a dual wavelength absorbance detector and a data workstation. UV absorbance was measured at a wavelength of 254 nm. The detector was calibrated before experiments and the responses from the detectors were linearly proportional to the concentrations. A ChiralCel OD-H column (250 mm × 4.6 mm, 5 m), purchased from Chiral Technologies (West Chester, PA, USA), was used for all the measurement. This column was packed with cellulose tris(3,5diethylphenyl carbamate) coated on 5-m silica-gel substrate. The column hold-up volume, Vm , was evaluated from the record of the detector set to 254 nm after injection of pure IPA into mixed mobile phase. The total column porosity (εt = 0.74) corresponding to the hold-up volume was determined from εt = Vm /V, where V is the geometric volume of the column. The extra-column volume between the autosampler to the inlet of column (Vex ≈ 180.0 l) was measured by removing the column from the system. All the retention data were corrected for this contribution. The overloaded band profiles of racemic mixture were used to derive the competitive isotherm parameters. Racemic feed solutions were obtained by dissolving the same amount of (R)and (S)-mandelic acids into the mobile phase. Experiments were performed at 25 ± 0.1 ◦ C with the elution flow rate of 0.7 and 0.9 ml/min. The injection volumes were 30.0 and 45.0 l and the total sample concentrations were 20.0, 30.0 and 40.0 g/l. 4. Results and discussion 4.1. Determination of the kinetic parameters Accurate determination of the two kinetic parameters, DL and km,i is of great importance for obtaining the correct isotherm parameters. In this study axial dispersion coefficient, DL, was determined from theoretical correlation [18]; lumped rate coefficient, km,i was determined along with the isotherm parameters by fitting the model predicted band profiles to the experimental data. For simplicity, same km value was used in the calculation for both components. DL and km at various flow rates are listed in Table 1. 4.2. Best-fit isotherm parameters Four overloaded elution band profiles of the racemic mixture with different injection volumes and feed concentrations were used to derive the isotherm parameters of different models. Each band profile used for fitting contains 560 points in the nonzero concentration range. Initially, some simulations were carried out to determine the possible ranges of each decision variable which were prerequisite in the optimization. Systematic optimization was then performed to determine the isotherm parameters by fitting the
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
37
Table 2 Competitive isotherm parameters by inverse method Langmuir qs (g/l) b1 (× 10−2 l/g) b2 (× 10−2 l/g) F(p)
62.46 5.6 7.5 48.9
biLangmuir q1s (g/l) b11 (× 10−2 l/g) b12 (× 10−2 l/g) q2s (g/l) b21 (l/g) b22 (l/g) F(p)
89.8 2.92 3.82 1.65 0.94 1.32 36.65
´ Toth qs (g/l) b1 (× 10−2 l/g) b2 (× 10−2 l/g) F(p)
88.56 4.4 5.84 0.702 37.39
model predicted band profiles to the four chromatograms simultaneously. Results of the isotherm parameters as well as the F(p) values for the three isotherm models are summarized in Table 2. Comparison of the three isotherm models is illustrated in Fig. 1. ´ isotherm Results from Table 2 indicate that biLangmuir and Toth ´ modfit better than Langmuir isotherm while biLangmuir and Toth els give rather similar representations of the experimental data. The second conclusion coincides with the observation from Fig. 1, ´ in which biLangmuir and Toth models generate almost the same equilibrium isotherm with the given mobile phase concentration. Such observation confirms the heterogeneous nature of this kind of chiral stationary phase. Fig. 2 demonstrates the results of fitting biLangmuir model to three of the measured chromatograms. It can be seen that the fit is less satisfactory at small sample size, while bands with large sample size fit better. The main reason of the error made is that the influence of high concentration data on the values derived for the numerical coefficients is higher than that of the low concentration points. However, this kind of error is somewhat acceptable, since inverse method provides more reliable results at higher column loading.
Fig. 1. Comparison of the three competitive isotherm models for (R)- and (S)mandelic acids.
Fig. 2. Best-fit of competitive biLangmuir model.
4.3. Validation of the acquired parameters In order to verify the validity of the derived isotherm parameters, experimental elution band profile with higher injection volume (Vinj = 60.0 l, cF,t = 40.0 g/l) was compared with the profile calculated from the derived isotherm parameters listed in Table 2. Good agreement between the calculated and measured band profiles are also acquired this time. The comparison of the calculated (with biLangmuir isotherm) and measured profiles is shown in Fig. 3. Experiments using various solute compositions were conducted as well to see whether the acquired isotherm parameters could well predict the elution band profiles at different operating conditions. Fig. 4 illustrates the comparison of the calculated and measured elution profiles at various solute compositions. Results from Fig. 4 indicate that isotherm parameters acquired from the inverse method can give rather good predictions of the elution band profiles at a wide range of solute compositions. Both the simulated concentration profiles (using the isotherm parameters of competitive biLangmuir model) and the measured elution profiles at higher elution flow rate (Q = 0.9 ml/min) are plotted in Fig. 5. Remarkably good agreements are observed for all the four chromatograms with different solute compositions. Except for the slight deviation of the tail of the second component, numerical solutions from the chromatography model with the acquired isotherm parameters can provide accurate
38
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
Fig. 6. Measured inlet concentration profile. Fig. 3. Comparison of model predictions (biLangmuir) to experimental band profiles at higher injection volume.
between the measured and predicted concentration profiles are observed in Figs. 2–5. Errors from the following four aspects may lead to such discrepancies. First, the simplified boundary condition, which assumes an ideal rectangular pulse injection, may result in the deviations. Measurement errors could be the second cause. It is known that the experimentally measured concentrations were derived from the original chromatograms (absorbance versus time) through a calibration. Both the drift of baseline and the noise of the signals will bring about some errors in the measured concentration profiles. The third reason could be the selected isotherm models. Because of the highly complex structure of the selector, chromatographic mechanism involved in the chiral discrimination has not yet been fully clarified. The selected isotherm models may lead to some theoretical errors in predicting the elution curves. Finally, method used to solve the PDEs may also generate some numerical errors in the simulated band profiles. Compared with the results obtained using finite difference method, calculated band profiles from OCFE ´ model. are not so smooth when smaller is used for Toth Fig. 4. Comparison of the model predicted concentration profiles (biLangmuir) with measured elution profiles at different solute compositions.
predictions on the elution times and peak shapes of the two components. The achievement of the good agreement between the calculated and measured elution profiles under various experimental conditions indicates that isotherm parameters obtained from the inverse method are quite robust and reliable. However, small deviations
Fig. 5. Comparison of the model predicted concentration profiles (biLangmuir) with the measured elution profiles at higher elution flow rate and different solute compositions.
4.4. Impact of the boundary condition The inlet profile constitutes the boundary condition of the differential mass balance equation, thus it should have some influence on the simulated band profiles [13]. Usually, the shape of the inlet profile was assumed to be a rectangle for the sake of simplicity. However, significant deviations from this assumption are always observed. Fig. 6 illustrates a measured inlet profile of mandelic acid, determined by injecting the sample with column removed from the system and replaced by a zerovolume connector. Obvious asymmetry of the injection profile can be seen from Fig. 6. How this asymmetry affects the elution
Fig. 7. Calculated elution profiles using real and simplified boundary conditions.
Y. Zhang et al. / J. Chromatogr. A 1202 (2008) 34–39
profiles and the final derived isotherm parameters remains a problem. In order to quantify the effect of the simplified boundary conditions, optimization runs were also performed with experimentally measured injection profile as boundary conditions. This time the comparison is confined to the two chromatograms with injection amounts of 0.9 and 1.8 mg. Results of the isotherm parameters using the true inlet profile are very close to the values listed in Table 2. Deviations of the final isotherm parameters using the true inlet profile from those using the simplified rectangular injection profiles are less than 2%. This demonstrates that the simplified boundary condition do not have a significant impact on the final results of the isotherm parameters and such a simplification can be used in calculation. This observation coincides with that reported recently by [15]. From our simulation results, isotherm parameters acquired with simplified boundary conditions are still valid in predicting the elution profiles when real injection profile is applied. Comparison of the calculated elution profiles using both real and simplified boundary conditions with isotherm parameters listed in Table 2 is illustrated in Fig. 7. 5. Conclusions Inverse method was used in this study to determine the competitive isotherm of (R)- and (S)-mandelic acids on a kind of cellulose chiral stationary phase using NSGA-II-JG as the optimization method. By fitting four chromatograms simultaneously with different injection volumes and different feed concentrations, the best-fit isotherm parameters for three different competitive ´ isotherm models, i.e., Langmuir, biLangmuir and Toth, were deter´ models both give mined. Results indicate that Bilangmuir and Toth good representations of the experimental data. Accuracy and reliability of the derived isotherm parameters were verified by comparing the model predicted band profiles with the experimental elution profiles under various experimental conditions, such as different injection volume, different solute compositions and different elution flow rates. Results indicate that the band profiles calculated using the derived isotherm parameters can provide rather accurate predictions for all experimental data. Isotherm parameters derived from the inverse method are reliable.
39
The inverse method offers a rather quick method of adsorption isotherm determination of racemic mandelic acid with minimal sample and solvent consumption. However, one important disadvantage of inverse method is that the acquired data depend on the column efficiency, which cannot be accessed conveniently under nonlinear conditions. Therefore, for more accurate isotherm, frontal analysis is recommended. It should also be stressed that the isotherm parameters acquired in this study are satisfactory only for dilute solutions of analytes. References [1] G. Subramanian, Chiral Separation Techniques—A Practical Approach, Wiley–VCH, Weinheim, 2001. [2] S. Ahuja, in: S. Ahuja (Ed.), Chiral Separations—Applications and Technology, American Chemical Society, Washington, DC, 1997, p. 1. [3] G. Guiochon, A. Felinger, D.G. Shirazi, A.M. Katti, Fundamentals of Preparative and Nonlinear Chromatography, Elsevier Academic Press, San Diego, CA, 2006. [4] F. Gritti, G. Guiochon, J. Chromatogr. A 1028 (2004) 197. [5] A. Seidel-Morgenstern, J. Chromatogr. A 1037 (2004) 255. ¨ [6] C. Blumel, P. Hugo, A. Seidel-Morgenstern, J. Chromatogr. A 865 (1999) 51. ´ T. Fornstedt, Anal. Chem. 76 (2004) 4856. [7] J. Lindholm, P. Forssen, ´ T. Fornstedt, Anal. Chem. 76 (2004) 5472. [8] J. Lindholm, P. Forssen, [9] A. Cavazzini, A. Felinger, G. Guiochon, J. Chromatogr. A 1012 (2003) 139. ´ [10] F. James, M. Sepulveda, Inv. Prob. 10 (1994) 1299. [11] A. Felinger, A. Cavazzini, G. Guiochon, J. Chromatogr. A 986 (2003) 207. [12] D. Antos, W. Piatkowski, K. Kaczmarski, J. Chromatogr. A 874 (2000) 1. [13] A. Felinger, D. Zhou, G. Guiochon, J. Chromatogr. A 1005 (2003) 35. [14] N. Marchetti, F. Dondi, A. Felinger, R. Guerrini, S. Salvadori, A. Cavazzini, J. Chromatogr. A 1079 (2005) 162. ´ T. Fornstedt, J. Chromatogr. A 1099 (2005) 167. [15] R. Arnell, P. Forssen, ¨ [16] C. Heuer, E. Kusters, T. Plattner, A. Seidel-Morgenstern, J. Chromatogr. A 827 (1998) 175. [17] G. Guiochon, B. Lin, Modeling for Preparative Chromatography, Academic Press, San Diego, CA, 2003. [18] D.M. Ruthven, Principle of Adsorption and Adsorption Process, Wiley, New York, 1984. [19] J. Toth, Acta Chem. Hung. 32 (1962) 31. [20] D.C.S. Azevedo, L.S. Paris, A.E. Rodrigues, J. Chromatogr. A 865 (1999) 187. [21] A. Felinger, G. Guochon, Chromatographia 60 (2004) S175. [22] K. Kaczmarski, G. Guiochon, Anal. Chem. 79 (2007) 4648. ˜ [23] I. Quinones, C.M. Grill, L. Miller, G. Guiochon, J. Chromatogr. A 867 (2000) 1. [24] B.A. Finlayson, Nonlinear Analysis in Chemical Engineering, McGraw-Hill, New York, NY, 1980. [25] R.B. Kasat, S.K. Gupta, Comp. Chem. Eng. 27 (2003) (1785). [26] V. Bhaskar, S.K. Gupta, A.K. Ray, Rev. Chem. Eng. 16 (2000) 1.