Numerical Approximation of the Nonequilibrium Model of Gradient Elution Chromatography Considering Linear and Nonlinear Solvent Strength Models.

In both linear and nonlinear chromatography, the lumped kinetic model is a suitable model for predicting elution bands when appropriate equilibrium functions and mass transfer coefficients are accessible. This model also works well in the case of gradient elution chromatography if variations in the equilibrium functions due to changes in the mobile phase composition are known. The rational selection of an optimum gradient is explored in this study from three different perspectives using the lumped kinetic model. Elution profiles generated by using (a) linear solvent strength, (b) quadratic solvent strength, and (c) power law are investigated. The effectiveness and reliability of the suggested numerical approach, utilizing the flux-limiting finite volume method, are demonstrated through numerical simulations. The impacts of axial dispersion, nonlinearity coefficient, Henry's constant, mass transfer coefficient, and gradient parameters are studied on single and two-component elution profiles.

Introduction

In high-performance liquid chromatography (HPLC), the mobile phase composition mostly influences the separation of complicated mixtures, especially when several closely related analytes are to be measured.[1] It is well-known that partitioning dominates the retention of tiny nonpolar compounds in reversed-phase liquid chromatography (RP-LC) mode.[2] Desirable separation of a mixture (including wide-range retentive components) cannot be achieved in an appropriate time frame in isocratic mode. A gradient (in solvent strength) can be used to address such elution issues.[3] The retention factor has the greatest impact on the peak width in the gradient mode. A gradient separation produces narrow and almost consistent peak widths and, therefore, provides equivalent detection sensitivity. It is commonly employed in a wide-range of chromatographic separations. Gradient elution chromatography has a well-established theory,[4] and several recent publications and reviews provide a comprehensive description of the method from theory to application.[5,6]

In recent years, the gradient elution technique has been significantly applied for more complicated separations in analytical chromatography. It is more commonly applied in forward-phase liquid–solid chromatography and reversed-phase chromatography.[7−10] Optimal elution algorithms enable the separation of multicomponent mixtures including solutes with significantly different retention characteristics. Analysis times for complicated mixtures with a broad range of retention factors can be significantly reduced in gradient elution than in isocratic separation of the identical mixtures, as mobile phase strength can be raised either gradually (e.g., linear gradient) or abruptly (step gradient) during sample elution.

This work has two objectives. The first one is to adapt and assess three models of gradient elution considering mobile phase variation. The second one is to compare their ability to predict retention behavior in liquid chromatography. It is critical to build models (stepwise, linear, and nonlinear solvent strength) that describe the retention factors as functions of mobile-phase composition, having the aim to obtain required separations efficiently. Numerous efforts have been made to create retention-time models since the late 1970s. Various retention-time models in RP-LC have been created, each with its own set of expressions defining the relationship between retention and mobile phase composition, e.g., the linear solvent strength model (LSS),[11] the quadratic retention-time model (QSS),[12−14] and the Neue model.[15] The most popular model is the LSS, which Snyder and Dolan developed in the 1990s. It characterizes isocratic retention in RP-LC as[11]where kHr, is the reference (extrapolated) value of K for ϕ0 (specifically in pure water), and α is the solvent strength parameter that is only appropriate in a limited range of ϕ values.[15] To accommodate this issue and enhance retention-time modeling, multimodal retention mechanisms have been introduced, such as the quadratic model[6] and the Neue and the Kuss model.[16]

In nonlinear chromatography, a variety of kinetic models have been presented.[16−18] Various chromatographic models try to model band migration in numerous ways. The general rate model takes into account a more or less sophisticated set of kinetic equations and performs a comprehensive study of the many steps in the chromatographic process. The lumped rate model takes into account just one kinetic process, known as the rate-controlling step, or a few such procedures, and it includes the involvement of the kinetics of the other processes in the rate constant. The EDM is the simplest and basic, combining all kinetic effects into a single apparent dispersion value. It is generally known that under some fundamental assumptions simpler models may be obtained from the comprehensive general rate model.[19,20]

While gradient elution studies frequently require the use of tiny and dilute samples, in preparative gradient elution chromatography, the column is frequently overloaded. As a result, the adsorption isotherms are nonlinear and competitive, and interference effects become significant. Consequently, mass transfer resistances can be quite high, particularly for macromolecules.

Many researchers have already developed their simulation programs to grasp the specifics of chromatography modeling. These programs, however, seldom attain sufficient numerical performance since they are frequently created for particular purposes. A specific amount of software engineering, project management, maintenance, and support is necessary for publishing a code. Most open-source programs were created in academic institutions or government research facilities. For instance, CADET is one of the freely available multitask simulation frameworks for column liquid chromatography.[21]

When it comes to the simulation of chromatographic systems, there are certain unique challenges. Sharp fronts are likely to develop in a variety of situations, including when the column is extremely efficient[18] and whenever there are self-sharpening effects (Langmuir-type isotherm). When there are no analytical solutions to the model equations, which is the case in most of the scenarios, to ensure precision, stability, and speed when dealing with sharp fronts, appropriate numerical techniques must be explored. It has been demonstrated that traditional numerical techniques, such as simple finite difference (FD), are incapable of efficiently capturing the real sharp fronts.[22]

Thus, in this paper, a high-resolution flux-limiting finite volume scheme is proposed to solve the model equations.[23−25] This numerical approach is especially useful for convection-dominated problems in which sharp peaks or fronts are generated. It preserves the mass-conservation property of the current model equations and is capable of capturing sharp fronts and peaks in the solutions.[23−25] The conventional method for resolving discontinuities and sharp fronts has been the flux estimation technique. Because the finite volume (FV) approach incorporates the physical concept of flux, many flux estimate techniques may be easily applied to it.[23−25] In all techniques, the last step of the solution is the same: depending on the kind of discretization, they provide a set of algebraic equations or ordinary differential equations (ODEs).

The remaining portion of this article is organized in the steps outlined here. Section introduces the nonlinear lumped kinetic model (LKM) for three separate gradient elution strategies. In Section , the suggested finite volume method for solving the given model equations is formulated. Section contains discussions on numerical case studies that demonstrate the model’s and numerical scheme’s efficacy. Section presents conclusions based on the outcomes of the discussion.

Lumped Kinetic Model

Instead of examining the entire interparticle concentration profile, the lumped kinetic model uses a linear driving force in the solid phase and only examines one extra parameter to supplement the axial dispersion coefficient. It combines the effects of internal and external mass transport resistances (inside one mass transfer coefficient). Two kinetic parameters, the axial dispersion coefficient (D) and the mass-transfer rate coefficient (KL,i), are related to the overall mass-transfer rate in a column. The dynamical characteristics of the gradient elution chromatographic separation process are studied using the LKM.

The following mass balance equation is used in the mobile phase:For each component, i, c, q, and q* denote the solute concentration in the mobile phase, nonequilibrium solute concentration in the stationary phase, and the equilibrium solute concentration in the stationary phase, respectively. u is interstitial velocity; t and z symbolize the time coordinate and the axial coordinate; ϕ ≔ ϕ(t, z) is the solvent concentration; and F = (1 – ε/ε) is the phase ratio in which ε represents the external porosity.

In this model, it is assumed that the adsorption–desorption and the diffusion processes in the mobile phase are very rapid. Therefore, in order to complete the model, the accumulation in the solid phase is evaluated using the fundamental linear driving force approach shown below:Relationships between the liquid and solid phase equilibrium adsorbate concentrations are termed as isotherms. The equilibrium adsorption data may be simulated using the isotherms. Researchers use isotherms to explore adsorption information, such as adsorption mechanisms, maximum adsorption capacity, and the characteristics of adsorbents. In the literature, a variety of different adsorption isotherm models have been established. The multicomponent Langmuir equation is used to model the adsorption equilibrium in this study.[26,27] When the mobile or liquid phase composition varies during the gradient elution process as a result of variations in a modifier concentration, the local equilibrium can be described asKH, denotes Henry’s coefficient, and the degree of nonlinearity associated with the isotherm for the i-th component of the mixture is quantified by b. Most of the time, it is necessary to conduct experiments to determine the functional dependence of the two isotherm parameters on ϕ.

Non-LSS QSS Model

The LSS model has been used to separate tiny molecules as well as macromolecules such as proteins and peptides. Furthermore, the LSS model is only viable in a small range of ϕ values. Due to these constraints, the following two gradient models are employed in this study to simulate analyte retention as a function of solvent fraction and compare their results to previously published work using the LSS model.[25] Polynomial equations can sometimes be utilized to more correctly characterize retention behavior than linear models.

Non-LSS Power-Law Model

In this case, the model parameters are dependent on the solvent concentration asHere, D, kLr,, kH, and b symbolize the reference values of the axial dispersion, mass transfer, Henry’s, and nonlinearity coefficients. Furthermore, ϕ is the volume fraction of modifying nonretained solvent, and the specific solvent strength parameters are denoted by α1, α2, γ1, and γ2.

To determine the distribution of the strong solvent of the mobile phase across the column, the solvent is considered not to be retained. Consequently, the following model is the most accurate for estimating changes in the volume fraction (ϕ(t, z)) of the strong solvent in the mobile phase:[25]with initial and boundary conditions:For a linear gradient:Here, Φ and β are the implemented profile and slope of the gradient, while ϕ0, ts and ϕe, te are the initial and final volume fractions and time of the gradient, respectively.

To complete the lumped kinetic model, the initial and boundary conditions are expressed asThese inlet boundary conditions (BCs) are:where i = 1, 2, 3, ... Nc. The i-th component-injected concentration is denoted as c, and tinj is the injection time.

Numerical Solution

There are several known numerical methods for estimating chromatographic models.[28,29] A semidiscrete high-resolution flux-limiting finite volume technique is used in this section to solve the model equations.[30] This approach has been recently applied on two models of gradient elution chromatography.[24,25]

Domain Discretization

The first step for implementing this technique is to discretize the computing domain. The main purpose of this discretization is to produce a set of time-coordinated coupled ODEs.

Here, zj (mesh points) are covered by the cells for 1 ≤ h ≤ N in the intervals such thatandThe averaged initial data are defined aswhere .

Integrating eqs and 3 over Ω givesThe differential and dispersion terms in eq can be approximated asThe solutions of eqs –11 giveThere are many different approaches and computational schemes available to estimate concentrations (or fluxes) at cell interfaces. Here, just the first and second order approximations are shown.

First-Order Approximation

At the cell interfaces, concentration values are estimated in eq by a backward difference formula. The first-order approximation for the concentration can be represented as

Second-Order Approximation

The concentration values at the cell’s interface c are approximately calculated using the flux-limiting calculations given below.[29]andWe need η = 10–10. The flux limiting function, i.e., Ψ, is employed to maintain the numerical scheme’s local monotonicity in eq which is defined as[29]

Results and Discussion

Single- and two-component samples are used as examples to demonstrate the benefits of gradient elution over isocratic elution and the necessity of choosing the right gradient technique. For simplicity, it is assumed that the mass transfer coefficient KL, = kL; the axial dispersion coefficients D = D; the nonlinearity coefficient b = b; and solvent strength parameters α = γ (for LSS), α1 = γ1, and α2 = γ2 (for QSS). All other parameters are listed in Table . The values of parameters in Table are chosen from the ranges typically used in HPLC. The plots display ϕ (the modulator concentration) and c (the solute concentration) at the column outlet z = L against the time t.

Table 1

Parameters for Single-Component Elution

parameters	values
column length	L = 10 cm
interstitial velocity	u = 1.0 cm/min
porosity	ε = 0.4
reference axial dispersion coefficient	Dzr = 0.0002 cm2/min
reference Henry’s constant	kHr = 3.5
reference mass transfer coefficient	kLr = 10 min–1
reference nonlinearity coefficient	bbref = 2.0
gradient start time	ts = 5 min
gradient end time	te = 90 min
initial concentration	ϕo = 0.1
final concentration	ϕe = 0.9
solvent strength parameter	α = 10.0
solvent strength parameter (for the non-LSS model (QSS))	α1 = 8.0
solvent strength parameter (for the non-LSS model (QSS))	α2 = 10.0
order of the power-law model	n = 1

Single-Component Elution

Dual comparisons are done in these case studies, i.e., numerical profiles plotted for different values of parameters and a comparison of three different gradient models mentioned above. Figure shows the comparison of numerical solutions for different values of reference Henry’s constant kH using the LSS and non-LSS models of gradient elution chromatography. For kH = 2 the best results can be seen for all gradient models. By inspection it is found that the non-LSS QSS model predicts the elution profiles better; i.e., narrow and symmetric peaks are generated. The non-LSS power-law model gives highly asymmetric peaks (adsorption–desorption is slow), and the retention time decreases for small values of Henry’s constant.

Figure 1

Influence of reference Henry’s constant on nonlinear single-component elution profiles by taking non-LSS and LSS models of gradient elution.

The effects of reference nonlinearity coefficient bref are given in Figure . For bref = 0, the elution profiles are shaped like a Gaussian curve; i.e., symmetric peaks are obtained. By increasing bref, the well-known formation of sharp adsorption and dispersed desorption fronts having shorter retention periods are observed. Once again, the non-LSS QSS model better predicts the elution profiles.

Figure 2

Influence of the reference nonlinearity coefficient on nonlinear single-component elution profiles by taking non-LSS and LSS models of gradient elution.

Figure displays the effects of reference axial dispersion coefficient D considering LSS and non-LSS models of gradient elution. The mean retention period remains unaffected for LSS and non-LSS models. The concentration profile gets less broadened as the value of D decreases for the non-LSS power-law model, but for LSS and non-LSS QSS models decreasing D does not show prominent changes on elution. Overall, the non-LSS QSS model better predicts the elution profiles among others.

Figure 3

Influence of reference dispersion coefficient on nonlinear single-component elution profiles by taking non-LSS and LSS models of gradient elution.

Effects of the gradient start time on nonlinear single-component elution profiles can be seen in Figure . For ts = 15 min, a distortion in peak shape appears, i.e., the peak spilt for the LSS model. This distortion can be minimized by implementing corrective actions like reducing the sample size or using a diluted solution or by reverse flow of the mobile phase. If the gradient starts late, i.e., for ts = 30 min, ts = 60 min, and ts = 80 min, the LSS and non-LSS QSS models show better separation (elution peaks merge for both cases), whereas for the non-LSS power-law model an increase in gradient start time results in asymmetrical peak shapes, increased peak heights, and longer run durations, which make the detection more difficult.

Figure 4

Influence of gradient start time on nonlinear single-component elution profiles by taking non-LSS and LSS models of gradient elution.

The effects of gradient end time on nonlinear single-component elution profiles for LSS and non-LSS models are presented in Figure . When the gradient ends early, i.e., for te = 20 min, the LSS model shows narrower peaks as compared to the non-LSS QSS model, while the non-LSS pwer-law model gives a peak with a sharp front and pronounced Langmuir effect. Overall, when the gradient time is increased, the peak heights are lowered. This is characteristic of gradient elution with micromolecule samples and corresponds to elution of each peak at a lower volume fraction of solvent values as the gradient duration rises.

Figure 5

Influence of gradient end time on nonlinear single-component elution profiles by taking non-LSS and LSS models of gradient elution.

In Figure , comparisons of isocratic and gradient elution are shown. These comparisons are divided in two cases for convenience, isocratic conditions with ϕ = 0 and (i) gradient elution with ϕ0 = 0.1, ϕe = 0.9, ts = 5 min, and te = 80 min and (ii) gradient elution with ϕ0 = 0.1, ϕe = 0.5, ts = 5 min, and te = 120 min. Overlapping peaks are obtained for non-LSS and LSS models for isocratic cases ϕ = 0. For the non-LSS QSS model, gradient elution is the best choice in both cases. Overall, better separation can be seen for the gradient case as compared to isocratic cases. It is worth noting that the peak shapes differ in both cases due to the change in gradient techniques.

Figure 6

Comparison of gradient and isocratic elution for nonlinear single-component elution profiles by taking non-LSS and LSS models of gradient elution.

Effects of solvent strength parameters for LSS and non-LSS models and n of the non-LSS power-law model are shown in Figure . It can be seen that profiles become narrower with short retention time as solvent strength parameters increase (i.e., the sensitivity of solvent strength parameters to modulate concentration) for both LSS and QSS models and, hence, improved separation or purification. Also, for the non-LSS power-law model, increasing n results in more prominent peak tailings.

Figure 7

Influence of solvent strength parameters α of the LSS model, α1 and α2 of the QSS model, and n of the non-LSS model (power) on nonlinear single-component elution profiles.

Comparison and Error Analysis of the Numerical Schemes

This test problem quantitatively analyzes the performance of the Koren technique compared to other available flux-limiting finite volume schemes. Figure depicts that the Koren technique generates the lowest errors. The reference solution was obtained over a grid of 1000 mesh cells. The magnified plots show that the first-order scheme produces the most diffusive elution profile, whereas the Koren scheme generates the most resolved solution. These facts lead us to the conclusion that the Koren scheme is suitable for solving models of gradient elution chromatography. The same performance of the Koren has already been verified in our previous article analytically and numerically in the case of isocratic elution.[28]

Figure 8

Comparison and error analysis of numerical schemes for single-component elution profiles.

Two-Component Elution

Real separation problems typically involve more than two components in a feed. Moreover, as observed above, nonlinear gradient shapes have additional potential for enhancing the process performance compared to the linear gradient. Due to these facts, a theoretical study of two-component elution is also included here, i.e., Nc = 2. In this case, the retention factors of components differ significantly, which reflects real situations. For the connection between the component-specific adsorption equilibrium constants, mass transfer, axial dispersion coefficient, and modulator concentrations, the LSS and non-LSS models are used. Here we take L = 10 cm as the column length, porosity as ε = 0.4, and injected concentrations as c1,inj = c2,inj = 1.0 mol/L, which are injected for a duration of tinj = 2.0 min at u = 0.6 cm/min. All other parameters are listed in Table . The modulator concentration (ϕ) and the solute concentrations (c) are shown at the column outlet over time t. In addition, the figures show the impact of gradient parameters on the elution profiles considering the aforementioned three gradient models.

Table 2

Parameters for Two-Component Elution

parameters	values
reference axial dispersion coefficient	Dzr = 0.0002 cm2/min
reference mass transfer coefficient	kLr = 10 min–1
reference Henry’s constant for component I	kHr,1 = 1.0
reference Henry’s constant for component II	kHr,2 = 3.5
reference nonlinearity coefficient for component I	b1bref = 1
reference nonlinearity coefficient for component II	b2bref = 2
solvent strength parameter	α = 0.90
solvent strength parameter (for the non-LSS model 1 (QSS))	α1 = 0.40
solvent strength parameter (for non-LSS model 1 (QSS))	α2 = 0.90
initial concentration	ϕ0 = 0.10
final concentration	ϕe = 0.9
gradient start time	ts = 5 min
gradient end time	te = 80 min

Comparisons under Gradient Conditions

A comparison of different choices of gradient models is presented in Figure . For the situation ϕ = ϕ0 = 0.1, the QSS model achieves generally superior separation outcomes for both components than other models, with narrower elution peaks and improved retention time. Moreover, the non-LSS power-law model completely fails to produce desired results at low concentration (ϕ = ϕ0 = 0.1) of the mobile phase, and we did not achieve the baseline separations for component II. For the second case ϕ = ϕe, the results are almost similar with narrower peaks having shorter retention times for QSS and LSS models. While in the positive (β > 0) and negative (β < 0) gradient situations for two-component elution, the non-LSS QSS model appears to be efficient, while the LSS and non-LSS power-law models yield better results.

Figure 9

Comparison of the LSS and non-LSS models under gradient conditions on nonlinear two-component elution profiles.

Comparison of Numerical Schemes

This case study compares the Koren scheme’s performance to other flux-limiting finite volume schemes. Initially cinit = 0 mol/L; i.e., the column is equilibrated with the solvent. Afterward, pulses of c1,inj = c2,inj = 1 mol/L are injected for a duration of tinj = 2.0 min. The column L has a length of 10 cm, ε = 0.4, u = 1 cm/ min, and N = 300. In addition, a grid of 50 mesh cells was utilized. Figure shows the numerical outcomes at the column outlet. All results are obtained in MATLAB R2015a with an intel core (TM) i5-6200U CPU, RAM 8.00 GB, 2.30 GHz, and Window 10 pro graphic card. For all LSS and non-LSS models of gradient elution, the Koren scheme gives more resolved peaks for both components. As far as computational time is considered the Koren scheme has minimum computational time as compared to other flux-limiting finite volume schemes. According to the aforementioned observations, the Koren technique is preferable for solving these gradient elution models.

Figure 10

Comparison of different schemes for two components with a nonlinear isotherm.

Conclusion

The separation in overload columns utilizing three different techniques of gradient elution chromatography was simulated using the nonlinear lumped kinetic model. The influences of different model parameters were considered to be dependent on the modulator’s concentration. Comparisons of LSS and non-LSS models were demonstrated, and their possible applications in gradient elution chromatography were explored. Although most of our test problems were focused on single-component elution, real separation problems typically involve more than two components in a feed. Therefore, a theoretical study of two-component elution was also included in the test problems. It was observed that nonlinear gradient shapes have additional potential for enhancing the process performance compared to the linear gradient. In the current study, it was observed that gradient elution enhances the production rate by decreasing the retention and cycling times and by possibly increasing the loading factors. The results obtained verify that gradient elution has the potential to outperform isocratic operation in both analytical and preparative chromatography. It is particularly suitable for the separation of mixtures containing components to have strong retentions. In addition to productivity enhancement, the time required for column regeneration is relatively short in the gradient elution chromatography. In nonlinear multicomponent chromatography, the gradient elution conditions can be influenced by the dependence of separation factors between the target component and its neighbors on the solvent concentration. Thus, the considered gradient models are rigorous and allow the evaluation of the influence of gradient elution on concentration profiles for a wide variety of operating conditions which are generally difficult to analyze in experimental research.