Evaluation of Aqueous and Nonaqueous Binary Solvent Mixtures as

Article pubs.acs.org/JPCC

Evaluation of Aqueous and Nonaqueous Binary Solvent Mixtures as Mobile Phase Alternatives to Water−Acetonitrile Mixtures for Hydrophilic Interaction Liquid Chromatography by Molecular Dynamics Simulations Sergey M. Melnikov,†,‡ Alexandra Höltzel,† Andreas Seidel-Morgenstern,‡ and Ulrich Tallarek*,† †

Department of Chemistry, Philipps-Universität Marburg, Hans-Meerwein-Strasse, 35032 Marburg, Germany Max-Planck-Institut für Dynamik komplexer technischer Systeme, Sandtorstrasse 1, 39106 Magdeburg, Germany

‡

ABSTRACT: The retentive principle of hydrophilic interaction liquid chromatography (HILIC) with a water (W)−acetonitrile (ACN) mobile phase (MP) and a hydrophilic stationary phase is the formation of a W-rich layer with a rigid and diffuse part in the immediate and extended surface region, respectively. Through molecular dynamics simulations of the adsorption of W−acetone (Ace), W−methanol (MeOH), and MeOH−ACN mixtures to a hydrophilic silica surface, we evaluate their MP potential for aqueous and nonaqueous HILIC. We analyze solvent and hydrogen-bond density profiles, solvent−surface coordination, as well as local solvent orientation, mobility, and composition. Our data show that W−Ace mixtures closely mimic the behavior of W− ACN mixtures, whereas W−MeOH mixtures fail as HILIC MP because the similar affinity of the silica surface for W and MeOH prevents preferential adsorption of W. MeOH−ACN mixtures form a rigid MeOH layer that reverses the surface polarity and prohibits formation of a diffuse MeOH layer in the extended surface region. Generally, a rigid layer at a hydrophilic surface is formed by binary mixtures whose solvents have sufficiently different hydrogen-bonding abilities, and a diffuse layer is formed when the rigid layer maintains hydrophilic properties.

I. INTRODUCTION Reversed-phase and hydrophilic interaction liquid chromatography (RPLC and HILIC, respectively) are the two most important liquid chromatography modes for chemical separations today. In both modes, partitioning between two phases of different polarity and different mobility (mobile versus stationary) plays a central role for the retention and separation of analytes. In RPLC, analytes of low-to-moderate polarity partition from an aqueous−organic mobile phase (MP), where water (W) is the weak solvent, into a layer of organic solvent molecules enriched on top of long hydrocarbon chains that are tethered to the solid support structure.1 In HILIC, analytes of moderate-to-high polarity partition from an aqueous−organic MP, where water (W) is the strong solvent, into an extended W-rich layer adsorbed at the hydrophilic surface of the stationary phase (SP).2 Analyte retention is considerably more difficult to predict in HILIC than in RPLC, because several mechanisms operate simultaneously. The relative contribution of each mechanism is generally unknown and subject to small changes in the experimental parameters. In the simplest case, HILIC retention results (i) from partitioning of analytes from the organic solvent-rich MP into the W-rich layer and (ii) from weak electrostatic interactions between analytes and the SP surface.3,4 If analyte and SP surface carry permanent and opposite charges, strong electrostatic interactions will likely dominate retention. Partitioning, weak and strong electrostatic interactions involve © 2014 American Chemical Society

the polar functional groups of analytes and SPs. Analytes (often) and SPs (occasionally) may also contain apolar moieties, so that hydrophobic interactions between analyte and SP cannot be excluded. Nevertheless, hydrophobic interactions are not a dominant force in aqueous (AQ)HILIC, whose defining characteristic is the retentive W-rich layer. In chromatographic practice, the choice of the organic solvent in the MP has turned out to be a much more critical parameter in HILIC than RPLC. Whereas several organic solvents are suitable MP components for RPLC, HILIC has essentially been restricted to ACN as the organic MP component. Attempts to replace the toxic and costly ACN with other organic solvents, pursued especially during the ACN shortage (2008−2009), have had limited success so far.5−7 Alcohols, methanol (MeOH) in particular, rarely gave sufficient retention,8−12 except with a zwitterionic SP that retains analytes through strong electrostatic interactions rather than partitioning.13 Tetrahydrofuran and dioxane worked in some cases10,14 but were ineffective in others.8,11 So far, acetone (Ace) remains the only solvent to consistently achieve sufficient HILIC retention.9,12,15,16 Received: November 6, 2014 Revised: December 11, 2014 Published: December 11, 2014 512

dx.doi.org/10.1021/jp511111z | J. Phys. Chem. C 2015, 119, 512−523

The Journal of Physical Chemistry C

Article

Figure 1. Adsorption of AQ and NA binary solvent mixtures to a hydrophilic silica surface. The MD simulation snapshots show the equilibrated state for a nominal ratio of 10/90 (v/v) W/ACN (A), W/Ace (B), W/MeOH (C), or MeOH/ACN (D) in the bulk reservoirs. Atoms are color-coded as follows: Si − yellow, O − red, H − white, N − blue, and C and CH3 − cyan.

by H atoms. With very few applications published,15−19 empirical data on retention in NA-HILIC are sparse. From a theoretical viewpoint, a retentive principle is difficult to envision for NA-HILIC considering that a W-rich layer cannot form under the experimental conditions. The theory of the W-rich layer has been developed by a series of molecular dynamics (MD) simulations.2,20,21 The first MD studies investigated the physicochemical properties of the W-rich layer that forms under HILIC conditions with a W− ACN MP inside the porous silica particles constituting the column bed.2,20 This setting was modeled by a silica nanopore equilibrated with W−ACN mixtures at varying solvent ratios. Variation of the W/ACN ratio proved to be a vital element of the study, as the HILIC effect in chromatographic practice only sets in at high ACN fraction in the MP. The studies revealed a

A complementary approach to varying the weak solvent of the MP was introduced by the group of Lindner as nonaqueous (NA)-HILIC.17 Bicker et al.17 used MeOH, ethanol (EtOH), or 1,2-ethanediol (Et(OH) 2 ) in place of W to separate nucleobases and nucleosides. The observed effect on retention depended on the type of SP, the type of alcohol, and the analyte. Retention of most analytes on diol SPs increased in the order W < Et(OH)2 < MeOH < EtOH as the strong solvent. A similar trend did not emerge with a bare-silica SP, where Et(OH)2 was less retentive than W, MeOH retained some but not all analytes stronger than W, and EtOH produced severely distorted peak shapes. Soukup and Jandera18 obtained increased retention of phenolic acids with a MeOH−ACN compared to a W−ACN MP but on a silica hydride SP where up to 95% of the silica surface OH groups have been replaced 513

dx.doi.org/10.1021/jp511111z | J. Phys. Chem. C 2015, 119, 512−523

The Journal of Physical Chemistry C

Article

relevant only at high W fraction in the MP (cf. Figure 11 in ref 26). Therefore, a surface of single silanol groups suffices to approximate the hydrophilic properties of a bare-silica SP under HILIC conditions. Apart from a significant reduction of computation time, the use of a planar surface instead of a curved pore model as in our earlier work2,20,21 eliminates the effects of a nanometer-sized confinement and the surface curvature. Confinement effects in hydrophilic silica pores (cylindrical or slit-shaped) filled with the neat solvents ACN, Ace, and MeOH27−33 as well as their AQ mixtures2,20,21,34−37 have been studied previously. Surface curvature has been addressed directly by the group of Siepmann, in the context of RPLC, for a (111) β-cristobalite SiO2 surface modified with C18 hydrocarbon chains.38 The hydrocarbon phase was found to extend deeper into the liquid phase when the C18 chains were grafted on the curved surface of a cylindrical pore with a 6 nm diameter than when grafted on a planar surface (cf. Figure 11 in ref 38). Because of the increased width of the bonded phase, the influence of the C18-modified silica on solvent density could be observed over a longer distance than with a planar surface. For an unmodified silica surface under HILIC conditions, the absence of curvature shortens the distance over which solvent density profiles recover bulk behavior, that is, the extension of the W-rich layer (both rigid and diffuse part). A possible influence of the surface curvature on the solvent distribution in its immediate vicinity cannot be excluded at this point but is not expected to be a decisive factor in the formation of a hydrophilicity gradient. The simple, planar model of a hydrophilic silica surface allows us to focus on the fundamental differences in the adsorptive behavior between W−ACN and W−Ace mixtures on the one hand and W−MeOH mixtures on the other hand to find a physicochemical explanation for their experimentally observed success or failure as MP for AQHILIC and to provide molecular-level details of the conditions in NA-HILIC with MeOH−ACN mixtures as MP.

strong resemblance between the W-rich layer and the electrical double layer formed by an electrolyte solution in contact with a charged surface. In analogy to the electrical double layer, the W-rich layer was divided into a rigid part in the immediate surface region (up to a distance of ca. 0.4 nm from the surface Si atoms of a cylindrical pore with 9 nm diameter) and a diffuse part in the extended surface region (from the rigid layer up to a distance of 1.5−2 nm from the surface Si atoms of the 9 nm pore). The rigid part consists predominantly of nearly immobile W molecules attached through hydrogen bonds (HBs) to the silica surface; the diffuse part begins with W molecules attached through hydrogen bonds (HBs) to the rigid layer, before solvent composition, structure, and mobility smoothly transition into bulk liquid properties. The terms “rigid” and “diffuse” refer not only to the translational mobility of the solvent molecules but also to the perceived influence of the hydrophilic surface and the nominal W/ACN ratio on the solvent composition in each layer. The composition in the rigid W layer is governed by the silica surface and is nearly independent from the nominal W/ACN ratio, whereas the composition in the diffuse W layer is influenced by surface properties and the nominal W/ACN ratio. The results of previous MD simulations involving a hydrophilic silica nanopore and W−ACN mixtures2,21 show that a hydrophilicity gradient runs from the bulk liquid (the MP) to the silica surface (the SP surface) and that this gradient becomes steeper with increasing ACN volume fraction in the MP. The hydrophilicity gradient is quantified as the ratio between the local W mole fraction (in rigid and diffuse W layer) and the W mole fraction in the bulk. Because analyte retention shows a similar dependence on the ACN volume fraction in the MP as the local-to-bulk ratio of the W mole fraction, the steepness of the hydrophilicity gradient is supposed to indicate the magnitude of the system’s retentive power. On the basis of the results of previous MD simulations, we have formed the hypothesis that the suitability of an AQ mixture as a HILIC MP depends on its ability to form a hydrophilicity gradient. Transferred to NA-HILIC, the hypothesis requires a gradient in the local-to-bulk ratio of the MeOH (instead of the W) mole fraction. In this work, we test our hypothesis by MD simulations of the adsorption of W− ACN, W−Ace, W−MeOH, and MeOH−ACN mixtures to a hydrophilic silica surface. To mimic the conditions of typical HILIC separations the simulations cover ratios between 30/70 and 2/98 (v/v) W/organic solvent or MeOH/ACN. Our simulation scheme recovers the conditions inside a particlepacked HILIC column equilibrated with the MP and ready for sample injection. The simulation box consists of a β-cristobalite SiO2 slab exposing the (111) face to solvent reservoirs at both sides (Figure 1). The (111) surface of β-cristobalite SiO2 is the simplest and most frequently used model for simulating the adsorption of neat organic solvents and their aqueous mixtures to a hydrophilic silica surface.22−24 The surface hydroxylation (4.55 OH groups/nm2), which closely resembles that of chromatographic silica substrates,25 results from only one type of functional group, namely single, isolated silanol groups. We recently investigated the contribution of different functional groups on the amorphous silica surface to overall surface hydrophilicity.26 We found that single and geminal silanol groups produce highly similar W excess adsorption isotherms and are largely responsible for overall surface hydrophilicity, whereas the specific properties of siloxane groups become

II. SIMULATIONS Simulation Box and Force-Field Parameters. Our fully periodic simulation box with a cross-section of 3.542 × 3.52 nm2 and a length of 5.93 nm contained a three-layer slab of βcristobalite SiO2 (0.93 nm thick) exposing a surface of single, isolated silanol groups to a solvent reservoir (2.5 nm long) at each side. The length of the solvent reservoirs was sufficient to capture the transition of the liquid properties from surfaceinfluenced to bulk behavior. The single silanols surface (with a hydroxylation of 4.5 OH groups/nm2 or 7.6 μmol OH groups/ m2) was created following an approach of Coasne et al.39 After cutting the slab from the initial crystal structure parallel to the (111) face, all undercoordinated Si atoms and unconnected O atoms were removed, and dangling bonds of O atoms were completed with H atoms. For the silica surface, we used the force-field parameters of Gulmen and Thompson,40 which provide van der Waals and electrostatic parameters for all silica surface atoms (Si, O, and H). The motion of all silica atoms was frozen during simulations except for the free rotation of H atoms. Solvent molecules were treated as rigid, three-site (W, MeOH, ACN) or four-site molecules (Ace). For W we used the force-field parameters of the simple point charge/extended (SPC/E) model,41 for ACN and MeOH the united-atom version of the transferable potentials for phase equilibria (UATraPPE).42,43 Why these particular force fields were chosen is explained in detail in previous works.2,20,21,44 For Ace, we used the Kirkwood−Buff derived force field (KBFF),45 which when 514

dx.doi.org/10.1021/jp511111z | J. Phys. Chem. C 2015, 119, 512−523

The Journal of Physical Chemistry C

Article

Determination of Solvent Translational Mobility. The translational mobility of solvent molecules was determined for the bulk region (Z > 2.0 nm) and for each peak in the solvent density profiles (for a space interval of ±0.01 nm around the center coordinate). The translational movement of solvent molecules perpendicular to the surface (the residence time τ) was calculated as the average time the solvent molecules remained within a given space interval. A shift of ±0.10 nm around the initial Z-position of a molecule was tolerated to account for perpetual molecular motion. The translational molecular movement of solvent molecules parallel to the surface (the parallel diffusion coefficient D∥) was calculated from the mean square displacements ⟨r2(t)⟩ that the O atom of W, the central C atoms of ACN and Ace, and the center of mass of MeOH molecules made parallel to the surface during the observation time t.49 The linear slope of the curve during the interval t = 4−10 ps was used to determine D∥ via the Einstein equation:

used in combination with the SPC/E model for W was shown to reproduce the thermodynamics and aggregation behavior of W−Ace mixtures in bulk solution reasonably well.46 Long-range electrostatic interactions were treated with the particle-mesh Ewald algorithm (with a real space cutoff of 1.2 nm). Nonbonded interactions were modeled with a 12-6 LennardJones potential and truncated at 1.2 nm. Lennard-Jones parameters for unlike interactions were calculated using Lorentz−Berthelot combination rules. MD Simulations. MD simulations were carried out with Gromacs 4.5.2.47 A Nosé−Hoover thermostat with a coupling constant of 0.1 ps was used for temperature control. The equations of motion were integrated with a time step of 1.0 fs. Each simulation run was performed for a canonical NVT ensemble, but the number of solvent molecules (N) in the simulation box was manually adjusted between successive simulation runs. At the start of a simulation run, the respective number of solvent molecules for each species of a binary mixture was determined from the nominal solvent ratio. Solvent molecules of each species were then distributed randomly in the reservoirs on the basis of nonoverlapping van der Waals radii. After energy minimization, initial velocities were randomly assigned according to a Maxwell−Boltzmann distribution before a simulation run was carried out at T = 300 K. After equilibration of the system, which took between 2 and 3 ns, the current solvent ratio in the bulk region of the reservoirs (at a distance of Z > 2 nm from the surface, whereby Z = 0 nm is located at the position of the Si atoms of the silanol groups) was checked and manually corrected toward the target value (tolerating a deviation of ±1%) by removing or adding the necessary number of solvent molecules from each species. Four to five of these preliminary simulations runs (each 12 ns long) followed by manual adjustment of the solvent ratio were necessary, before a productive simulation run (24 ns long) could be performed. For data analysis, configurations in the second half of the trajectory were saved in 1 ps intervals. Ratios of 30/70, 20/80, 15/85, 10/90, 5/95, and 2/98 (v/v) W/ organic solvent or MeOH/ACN were simulated for each investigated mixture. Data analysis was performed for mixtures of 30/70, 10/90, and 2/98 (v/v) W/organic solvent or MeOH/ ACN; the number of molecules for each solvent species in the productive simulations runs is given in Table 1. The intermediary solvent ratios appear as additional data points in Figure 6.

D = lim t →∞

W/ACN

W/Ace

W/MeOH

III. RESULTS AND DISCUSSION Solvent Density Profiles. Figure 2 shows solvent density profiles for the investigated binary mixtures at nominal ratios of 30/70, 10/90, and 2/98 (v/v) W/organic solvent or MeOH/ ACN. All solvent species in the study contain a HB-accepting heteroatom X (X = N, O), which is capable of direct interaction with the surface OH groups. Therefore, solvent density was profiled through the number density of the X atom. We begin by comparing the established W−ACN system (Figure 2A) with W−Ace (Figure 2B) and W−MeOH mixtures (Figure 2C). At first glance, the W profiles look quite similar for all AQ mixtures. Three peaks can be identified: a surface peak (1A) with an accompanying shoulder peak (1B) at Z = 0.25 and 0.35 nm, respectively, a smaller peak (2) at Z = 0.51 nm, and finally a weakly defined peak (3) at Z = 0.73 nm. At all solvent ratios, the maximum W density is contained in the surface peak (1A), as expected for AQ mixtures at a hydrophilic surface. At second glance, an important difference between the W profiles of the W−MeOH system and those of the W−ACN and the W−Ace systems can be detected, concerning the dependence of the surface W density from the nominal solvent ratio as well as the absolute peak densities at and near the surface. In W−ACN and W−Ace systems, the density of the surface W peak (1A) responds only mildly to a sinking W fraction in the bulk liquid, whereas the W−MeOH system shows a strong correlation between surface W density and bulk W fraction. Between 30/70 and 2/98 (v/v) W/organic solvent, the surface W density decreases by less than 20% in the W−ACN and W−Ace systems but by more than 90% in the W−MeOH system. The comparison of the W surface peak densities also shows that W− ACN and W−Ace mixtures accumulate more W at the hydrophilic silica surface than W−MeOH mixtures. Even at 2/98 (v/v) W/ACN or W/Ace, the W surface peak density

MeOH/ACN

NW

NACN

NW

NAce

NW

NMeOH

NMeOH

NACN

30/70 10/90 2/98

680 295 135

433 562 604

737 345 141

312 405 453

621 200 40

602 757 826

325 165 60

403 530 606

(1)

Following standard Gromacs protocol, an error estimate was calculated as the difference between values of the diffusion coefficient obtained from the slope of the curve over the two halves of the fit interval (i.e., during time intervals of 4−7 ps and 7−10 ps). Only solvent molecules that remained in the tolerated space interval (±0.1 nm) during the whole observation time were counted.

Table 1. Number of Solvent Molecules (N) in Productive Simulation Runs (v/v)

⟨r 2(t )⟩ 4t

Analysis of Hydrogen Bonding. We applied three geometric criteria to determine whether two solvent molecules (of the same or another species) or a solvent molecule and a surface OH group were positioned within hydrogen-bonding distance to each other: (i) distance between donor O and acceptor X atom (X = N, O)