Understanding the large solubility of lidocaine in 1-n-butyl-3-methylimidazolium based ionic liquids using molecular simulation Ryan T. Ley and Andrew S. Paluch
Citation: J. Chem. Phys. 144, 084501 (2016); doi: 10.1063/1.4942025 View online: http://dx.doi.org/10.1063/1.4942025 View Table of Contents: http://aip.scitation.org/toc/jcp/144/8 Published by the American Institute of Physics
THE JOURNAL OF CHEMICAL PHYSICS 144, 084501 (2016)
Understanding the large solubility of lidocaine in 1-n-butyl-3-methylimidazolium based ionic liquids using molecular simulation Ryan T. Ley and Andrew S. Paluch Department of Chemical, Paper and Biomedical Engineering, Miami University, Oxford, Ohio 45056, USA
(Received 3 December 2015; accepted 3 February 2016; published online 23 February 2016) Room temperature ionic liquids have been proposed as replacement solvents in a wide range of industrial separation processes. Here, we focus on the use of ionic liquids as solvents for the pharmaceutical compound lidocaine. We show that the solubility of lidocaine in seven common 1-n-butyl-3-methylimidazolium based ionic liquids is greatly enhanced relative to water. The predicted solubility is greatest in [BMIM]+[CH3CO2]−, which we find results from favorable hydrogen bonding between the lidocaine amine hydrogen and the [CH3CO2]− oxygen, favorable electrostatic interactions between the lidocaine amide oxygen with the [BMIM]+ aromatic ring hydrogens, while lidocaine does not interfere with the association of [BMIM]+ with [CH3CO2]−. Additionally, by removing functional groups from the lidocaine scaffold while maintaining the important amide group, we found that as the van der Waals volume increases, solubility in [BMIM]+[CH3CO2]− relative to water increases. C 2016 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4942025]
I. INTRODUCTION
Room temperature ionic liquids (ILs) have received tremendous attention for use as replacement solvents in a wide range of industrial processes.1,2 The cation and anion may be changed and modified to tune their chemical and physical properties. Additionally, due to their ionic nature, ILs have negligible vapor pressures, minimizing solvent losses. Recently, the application of ILs for biological processes has emerged.3,4 One application has been to use ILs as solvents in crystallization processes.5–8 Weber et al.8 found the equilibrium solubility of acetaminophen in the IL 1-ethyl3-methylimidazolium acetate to be very large due to hydrogen bonding between acetaminophen (hydrogen bond donor) and acetate (hydrogen bond acceptor). They found that the solubility of acetaminophen could therefore be tuned with the addition of a hydrogen bond donating co-solvent (an “antisolvent”) that would compete with acetaminophen to hydrogen bond with the acetate anion. In this fashion efficient crystallization processes could be developed. In addition, the development of biologically active ILs has been proposed.3,9 In addition to advantages of using an IL over a conventional liquid, biologically active ILs are advantageous as one can obtain a high concentration of the active ingredient(s). Lidocaine is an important local anesthetic typically applied topically. A great deal of research is devoted to developing advanced controlled topical delivery methods.10,11 It has been shown that similar results may be obtained using a lidocaine cation based IL.3,9 These lidocaine based ILs form a gel that may readily be applied topically with a high lidocaine concentration. Following Refs. 8 and 9, here we study the solubility of (neutral) lidocaine in seven common ILs with the 0021-9606/2016/144(8)/084501/10/$30.00
1-n-butyl-3-methylimidazolium cation (see Figs. 1 and 2). Using molecular simulation free energy calculations, we predict that in all seven ILs, the solubility is greatly increased relative to water, with enhancements ranging from 103 to 105. Like acetaminophen, lidocaine contains an amide group capable of donating hydrogen bonds. However, as compared to Ref. 8, upon detailed investigation using molecular structural analysis, we find that while the ability of lidocaine to hydrogen bond with the anion is important, so too are the interactions between the ionic liquid cation and anion. For the case of lidocaine in [BMIM]+[CH3CO2]−, we find that the favorable solvation results from favorable hydrogen bonding between the lidocaine amine hydrogen and the [CH3CO2]− oxygen, favorable electrostatic interactions between the lidocaine amide oxygen with the [BMIM]+ aromatic ring hydrogens, while lidocaine does not interfere with the association of [BMIM]+ with [CH3CO2]−. Essentially, lidocaine is able to interact favorable with two separate associated ion pairs. The result is interesting for lidocaine given its large size and chemical complexity. Additionally, maintaining the amide group of lidocaine responsible for the favorable interactions with [BMIM]+ and [CH3CO2]−, we considered three additional solutes wherein functional groups were removed from the lidocaine scaffold. This had the effect of decreasing the solute van der Waals volume. We found that as the van der Waals volume increased, the solubility in [BMIM]+[CH3CO2]− relative to water increased. The application of ILs as solvents in controlled release applications, such as IL supported membranes or ILs separated from skin by a semipermeable membrane are extremely promising. In addition, the results demonstrate the ability to use ILs for crystallization and the extraction of lidocaine from aqueous environments.
144, 084501-1
© 2016 AIP Publishing LLC
084501-2
R. T. Ley and A. S. Paluch
J. Chem. Phys. 144, 084501 (2016)
FIG. 1. The chemical structure of the studied solutes consisting of lidocaine (CAS 137-58-6) and three solutes resulting from removing functional groups from the lidocaine scaffold.
achieved in a simulation by having a single solute molecule in a system with a large number of solvent molecules and offers the advantage of isolating solute-solvent intermolecular interactions. The present study focuses on the pharmaceutical solute lidocaine (see Fig. 1) in the seven ILs composed of the 1-n-butyl-3-methyl-imidazolium cation ([BMIM]+) with one of the seven anions shown in Fig. 2. As we will show, we find that the solubility of lidocaine in the IL [BMIM]+[CH3CO2]− is 105 times greater than in water, which is larger than for any of the other studied ILs. The mechanism responsible for this favorable solvation of lidocaine in [BMIM]+[CH3CO2]− is elucidated by performing structural analyses and additional free energy calculations. To facilitate the structural analysis, the concentration of lidocaine in [BMIM]+[CH3CO2]− was increased to a mole fraction of approximately 0.1 and 0.3. Next, solubilities in [BMIM]+[CH3CO2]− relative to water were computed for three additional solutes wherein functional groups were systematically removed from the lidocaine scaffold. The chemical structure of the additional solutes is provided in Fig. 1.
III. COMPUTATIONAL DETAILS A. Force fields
FIG. 2. The chemical structure of the seven anions and the 1-n-butyl-3methylimidazolium cation ([BMIM]+) which constitute the seven studied ionic liquids.
II. METHODOLOGY
The solubility of a solid solute (component 1) in an IL (solvent i) relative to water may be described by the classical equations of phase-equilibria thermodynamics. Assuming the solute concentration is sufficiently small so as to be considered infinitely dilute, the relative solubility may be computed using molecular simulation free energy calculations as12,13 c1,i = β µres,∞ ln 1,water (T, p, N1 = 1, Nwater) c1,water − β µres,∞ 1,i (T, p, N1 = 1, Ni ) ,
(1)
where c1,i and c1,water are the molar or mass concentrations of the solute in solvent i and water at equilibrium, respectively, β µres,∞ and β µres,∞ 1,i 1,water are the dimensionless residual chemical potential of the solute infinitely dilute in solvent i and water, respectively, where β −1 = k BT and k B is the Boltzmann constant, T is the temperature, p is the pressure, and N1, Ni , and Nwater are the number of solute, solvent i, and water molecules in the simulation. The infinite dilution limit is
Interactions were modeled using a “class I” potential energy function where all non-bonded intermolecular interactions (Unb) were accounted for using a combined Lennard-Jones (LJ) plus fixed point charge model of the form14,15 ( σ ) 12 ( σ ) 6 1 qi q j ij ij − + , (2) Unb r i j = 4ε i j r i j 4πε 0 r i j r i j where r i j is the separation distance between sites i and j, ε i j is the well-depth of the LJ potential, σi j is the distance at which the LJ potential is zero, and qi and q j are the partial charges of sites i and j, respectively. The united atom force field of Zhong et al.16 was used to model all of the ILs and water was modeled with TIP4P,17 which has been shown to work well using TraPPE force field models.18,19 The force field for lidocaine and the studied solutes wherein functional groups were removed from the lidocaine scaffold were all constructed based on the Explicit Hydrogen Transferable Potentials for Phase Equilibria (TraPPE-EH) force field.20–22,86 The LJ parameters for the phenyl group were taken directly from TraPPE-EH.20 The LJ parameters for N and H in the secondary amide were taken to be the same as for N and H in a secondary amine, and parameters for the primary and tertiary amine were all taken from TraPPE-EH.23,24 The carbonyl group LJ parameters were taken from the TraPPE-UA amide model,23,25 and all of the alkyl (CH3 and CH2) groups were modeled as a single alkyl pseudoatom with TraPPE-UA LJ parameters.26,27 Partial atomic charges were obtained for the TraPPEEH based models in a similar fashion as the original TraPPE-EH parameterization.20–22 Slight differences result due to the capabilities of the electronic structure software package available to us at the time of this study and the
084501-3
R. T. Ley and A. S. Paluch
use of united-atom alkyl sites. Charges were parameterized following a three step procedure. First, the gas-phase structure for each solute was optimized at the M06-2X/cc-pVTZ level of theory/basis set.28,29 Second, a single point energy calculation was performed on the gas-phase optimized structure at the M06-2X/6-31G(d) level of theory/basis set28,29 in a self-consistent reaction field (SCRF) using the SMD universal solvation model for 1-octanol.30–32 All of the electronic structure calculations were performed using Gaussian 09.33 Third, partial atomic charges were then obtained from the electrostatic potential (obtained in step 2) using the restrained electrostatic potential (RESP)34,35 method in ANTECHAMBER (part of the AMBER 12 simulation suite).36,37 During the fitting procedure, the partial charges of alkyl hydrogens were fixed to zero; this had the effect of using a single charge site at the center of the alkyl group appropriate for use with a united atom model. All of the intramolecular parameters for the solutes were taken from the General AMBER Force Field (GAFF).38,39 Parameters were generated using ANTECHAMBER and converted from AMBER to GROMACS format using ACPYPE.40,41 Since alkyl groups were modeled as a single alkyl pseudoatom, intramolecular parameters involving alkyl hydrogens were removed. Throughout the present study, solute and IL bonds involving hydrogens were held fixed. All of the GROMACS force field files used in the present study are provided in the supplementary material of this manuscript.42 B. Molecular dynamics 1. Ionic liquids
As a result of their sluggish dynamics, extensive equilibrations were performed for the IL systems. Simulations were performed for neat (pure) ILs, ILs plus a single lidocaine molecule, which corresponds to infinite dilution, [BMIM]+[CH3CO2]− with a lidocaine mole fraction of approximately 0.1 and 0.3, and [BMIM]+[CH3CO2]− plus a single solute molecule which results from systematically removing functional groups from lidocaine. All of these simulations were carried out following the same procedure. First, Packmol was used to generate initial structures.43,44 The number of IL pairs for the neat IL simulations were chosen to give a cubic box length of approximately 4.5 nm at 298.15 K based on the simulation results of Ref. 16. The number of IL pairs remained the same for systems containing solutes and the box size was adjusted accordingly. Subsequently, between 3000 and 6000 steepest descent minimization steps were performed to remove any bad contacts that might have resulted from the packing. The next two steps involved dynamics with the equations of motion integrated using the Verlet leap-frog algorithm.14,15,45,46 The system was first equilibrated in an NVT ensemble at 738.15 K for 1 ns using stochastic velocity rescaling.45,47–49 Next, the system was brought from 738.15 K to 298.15 K and 1 bar using simulated annealing using stochastic velocity rescaling and the Berendsen barostat.45,46,50
J. Chem. Phys. 144, 084501 (2016)
The system was initially carried out for 1 ns at 738.15 K and 1 bar. Over the next 8.8 ns, the temperature was decreased linearly at a rate of 50 K/ns to reach the final temperature of 298.15 K and was then held at 298.15 K for 1 ns. Finally, a 13 ns NpT simulation was performed at 298.15 K and 1 bar with the equations of motion integrated using the velocity Verlet algorithm.14,15,45,46 The simulation used the Nosé-Hoover chain thermostat (with a chain length of 10) and the Martyna-Tuckerman-TobiasKlein (MTTK) barostat.15,45,51,52 To ensure accurate statistics when performing the structural analyses, the simulations for neat [BMIM]+[CH3CO2]− and [BMIM]+[CH3CO2]− plus approximately 0.1 and 0.3 mol fraction lidocaine were carried out for an additional 101 ns. For all three steps involving dynamics, the time constant for the thermostat was 1 ps and the time constant for the barostat was 4 ps. Free energy calculations in the NpT ensemble were conducted for systems with a single solute molecule in the IL, which will be elaborated upon momentarily. 2. Water
The initial structure for water was taken from a previously equilibrated system at ambient conditions. A single solute molecule (lidocaine or a solute resulting from removing functional groups from the lidocaine scaffold) was solvated in these systems in a cubic box approximately 4.5 nm in length using the GROMACS tool “genbox.” The system was then equilibrated for 2.5 ns using a Nosé-Hoover chain thermostat (with a chain length of 10) and the MTTK barostat with a time constant of 1 and 4 ps for the thermostat and barostat, respectively, using the velocity Verlet algorithm. Free energy calculations in the NpT ensemble were then conducted which will be elaborated upon momentarily. All of the molecular dynamics simulations in this study were performed using GROMACS 4.6.3.53–55 Atom bond lengths involving a hydrogen were constrained using P-LINCS45,56,57 with the leap-frog Verlet algorithm and using SHAKE45,58,59 with the velocity Verlet algorithm. The Verlet neighbor list was used45 and LJ interactions were cutoff at 1.4 nm. Long-range analytic dispersion corrections were applied to the energy and pressure to accommodate the truncation.14,15,45,46 Lorentz-Berthelot mixing rules were used for unlike LJ sites.14 The electrostatic terms were evaluated with the smooth particle-mesh-Ewald (SPME) method with tin-foil boundary conditions45,46,60 with real space interactions truncated at 1.4 nm. The SPME B-spline was order 4, the Fourier spacing was 0.12 nm, and the relative tolerance between long and short-range energies was 10−4. 3. Free energy calculations
The free energy calculations were performed at 298.15 K and 1 bar using a similar protocol as our recent work.61 The infinite dilution residual chemical potential, µres,∞ 1,i , of the solute in the ILs and water was calculated using a multi-stage free energy perturbation method62–66 with the multi-state Bennett’s acceptance ratio (MBAR) method.67–70,87
084501-4
R. T. Ley and A. S. Paluch
J. Chem. Phys. 144, 084501 (2016)
A “soft-core” potential was used to decouple the solutesolvent intermolecular LJ interactions. Stage (m) dependent elec decoupling parameters, λ LJ controlled the LJ m and λ m and electrostatic intermolecular interactions, respectively.
The decoupling parameters varied from 0 to 1. When elec λ LJ m = λ m = 1, the solute is fully coupled to the system. When elec λ LJ = λ m m = 0, the solute is decoupled from the system. The “soft-core” potential had the form71–73,87
σi6j σi12j sc , − ULJ r i j ; m = 4λ LJ ε m ij 2 LJ 6 6 1 − λ LJ α σ 6 + r 6 1 − λ α σ + r LJ m ij ij LJ i j m ij
where αLJ is a constant, which had a value of 1/2. The advantage of using a “soft-core” potential to decouple the LJ interactions is that while it yields the correct limiting value of the potential (when λ LJ m = 0 and 1), it additionally allows nearly decoupled molecules to overlap with a finite energy (and hence finite probability). The electrostatic term in the intermolecular potential was decoupled linearly as 1 qi q j . Uelec r i j ; m = λ elec m 4πε 0 r i j
(4)
At each stage m, an independent MD simulation was performed. The simulation time for each stage m was 18.5 ns, where the first 0.5 ns was discarded from analysis as equilibration. The initial structure for these simulations was taken as the final structure from the NpT simulation of a single solute in solution (either IL or water). The change in the Hamiltonian with the current configuration between stage m and the other stages is computed every 0.20 ps. This is saved for subsequent post-simulation analysis with MBAR70 to determine µres,∞ 1,i . This analysis was performed using the Python implementation of MBAR (PyMBAR) and the GROMACS analysis script distributed with it.74–76 A total of 15 different stages were used for the free energy calculations where m = 0 corresponds to a non-interacting (ideal gas) state and m = 14 is a fully interacting system. From m = 1 to 10, the LJ interactions were increased from λ LJ m = 0.1 to 1.0 in 10 equal increments of 0.1. Electrostatic interactions were increased in a square root fashion following 77 λ elec m = {0.50, 0.71, 0.87, 1.00} from m = 11 to 14. The simulation parameters for the free energy calculations were the same as the last step of equilibration (i.e., using velocity Verlet with a Nosé-Hoover chain thermostat and MTTK barostat) with one exception: the efficient Verlet cut-off scheme is not available for use with free energy calculations in GROMACS 4.6.3. Instead, the LJ interactions were switched off smoothly between 1.15 and 1.2 nm, and the real space part of the electrostatic interactions was switched off smoothly between 1.18 and 1.2 nm. Long range corrections were applied as described previously. C. Structural analysis
As will be shown in Section IV, lidocaine was found to have the greatest affinity for the IL [BMIM]+[CH3CO2]−. To understand the underlying mechanism, detailed structural analysis was performed. Specifically, radial distribution
(3)
function (RDF) and spatial distribution function (SDF) calculations were performed for both the neat IL and the IL with approximately 0.1 and 0.3 mol fraction lidocaine. In both cases, the number of IL pairs was constant and 273, and the number of lidocaine molecules was 30 in the 0.1 mole fraction case and 90 in the 0.3 mol fraction case. To ensure accurate statistics, the RDFs and SDFs were computed from long 101 ns trajectories. The RDFs were computed using the GROMACS tool “g_rdf” over the last 100 ns of the trajectory. The SDFs were computed using the TRAVIS software.78–80,88,89 IV. RESULTS AND DISCUSSION A. Free energy calculations
First, we predict the solubility of lidocaine in the seven ILs relative to water. The ILs all consisted of a [BMIM]+ cation paired with one of seven common anions (see Fig. 2). In the order of least to most soluble, the anions are ranked as [BF4]− < [CH3SO4]− < [CF3SO3]− = [PF6]− < [TF2N]− < [CF3CO2]− < [CH3CO2]−. The corresponding values of ln ci /cwater are, respectively, {7.64, 8.55, 9.95, 9.94, 10.24, 10.63, 11.84} where the subscript corresponds to the uncertainty in the last digit. In all cases, the solubility of lidocaine in the ILs is larger than in water, with enhancements ranging from 103 to 105. This result is extremely encouraging. Even for the case of the hydrophobic ILs, [BMIM]+[BF4]−, [BMIM]+[PF6]−, and [BMIM]+[TF2N]−, the solubility relative to water is very large. These three ILs would be extremely promising for the extraction of lidocaine from aqueous environments. The largest solubility observed was for lidocaine in [BMIM]+[CH3CO2]−. For this case, the solubility was approximately 105 orders of magnitude greater than in water.81,90 Recently, Weber et al.8 demonstrated how mixtures of [EMIM]+[CH3CO2]− with a hydrogen bond donating co-solvent can be used to develop advanced separation processes. In that study, they found the solubility of acetaminophen in a range of ILs increased with increasing IL hydrogen bond basicity.82 The hydrogen bond basicity corresponds to the ILs ability to accept hydrogen bonds. Similar to acetaminophen, lidocaine also contains an amide group which can act as a hydrogen bond donor (N—H), so we might expect similar results here. If hydrogen bonding between lidocaine and the IL was the
084501-5
R. T. Ley and A. S. Paluch
J. Chem. Phys. 144, 084501 (2016)
dominate (or determining) interaction present, we would predict the trend [PF6]− < [TF2N]− < [BF4]− < [CF3SO3]− < [CH3SO4]− < [CF3CO2]− < [CH3CO2]−.82 However, this ranking differs from that observed. Therefore, while solubility is the greatest in [BMIM]+[CH3CO2]−, this suggests that it does not solely result from hydrogen bonding between lidocaine and the IL (anion). There must be additional molecular level driving forces responsible for the favorable solubility of lidocaine. To elucidate the underlying driving forces, we will next perform detailed structural analyses for neat [BMIM]+[CH3CO2]− and for systems of [BMIM]+[CH3CO2]− with approximately 0.1 and 0.3 mol fractions of lidocaine. B. Structural analysis
The long 101 ns simulations of neat [BMIM]+[CH3CO2]− and [BMIM]+[CH3CO2]− with 0.1 and 0.3 mol fractions of lidocaine were used to compute RDFs and SDFs to investigate the underlying solvation mechanism. 1. Ionic liquid
First we consider the ionic liquid structure. RDFs associated with the interaction of [BMIM]+ with [CH3CO2]− are shown in Fig. 4 and are computed both with respect to the center-of-ring (COR) of [BMIM]+ and the center-ofmass (COM) of [CH3CO2]−, and between the aromatic ring hydrogens of [BMIM]+ and the [CH3CO2]− oxygens. We find that [BMIM]+ and [CH3CO2]− are well ordered at short distances. The dominant interaction is between the aromatic hydrogens of [BMIM]+ with the [CH3CO2]− oxygens. We observe a high intensity peak in the RDF between 0.2 and 0.21 nm which decreases sharply. This short, well defined interaction is characteristic of a weak hydrogen bond (due to electrostatic interactions).83 In the supplementary material (Fig. S2),42 we provide RDFs for each aromatic hydrogen of [BMIM]+ with the [CH3CO2]− oxygens. These RDFs suggest that the strongest interaction is with aromatic ring hydrogen H1 (see Fig. 3), although the difference with the other hydrogens is rather small. After the first peak, g (r) remains close to one. Additionally, in Fig. 4, we find that the RDFs between [BMIM]+ and [CH3CO2]− change little with increasing lidocaine concentration. The peak locations all remain the same, with the intensity of the first peak increasing slightly with increasing concentration. SDFs for [CH3CO2]− around [BMIM]+ in the neat IL and with a lidocaine concentration of 0.1 and 0.3 mol fractions are provided in Fig. S3 of the supplementary material.42 The SDFs appear to change little with increasing lidocaine concentration, with [CH3CO2]− populating the area around the [BMIM]+ aromatic ring hydrogens, in agreement with the RDFs. Figure 5 shows the coordination number (CN) associated with the interaction of [BMIM]+ with [CH3CO2]− computed with respect to the COR of [BMIM]+ and the COM of [CH3CO2]−. We find that as the lidocaine concentration increases, the coordination number remains equal to 1 at
FIG. 3. Ball-and-stick representation of the united atom model of [BMIM]+ with the aromatic ring hydrogens labeled to facilitate discussion of the results.
a distance of approximately 0.47 nm. This distance is equal to the location of the first peak in the corresponding RDF. We therefore find that in the neat IL and with a lidocaine concentration of 0.1 and 0.3 mol fractions, we always have a [BMIM]+ associated with a [CH3CO2]− via a weak hydrogen bond.
FIG. 4. Radial distribution function (RDF) for the interaction between [BMIM]+ and [CH3CO2]− where the distance is computed with respect to the center-of-ring (COR) of [BMIM]+ and the center-of-mass (COM) of [CH3CO2]− (top pane), and between the aromatic hydrogens of [BMIM]+ with the [CH3CO2]− oxygens (bottom pane). A RDF is shown for the neat IL (x 1 = 0) and the IL plus 0.1 and 0.3 mol fractions of lidocaine.
084501-6
R. T. Ley and A. S. Paluch
FIG. 5. Coordination number (CN) for the interaction between [BMIM]+ and [CH3CO2]− where the distance is computed with respect to the center-of-ring (COR) of [BMIM]+ and the center-of-mass (COM) of [CH3CO2]−. A CN is shown for the neat IL (x 1 = 0) and the IL plus 0.1 and 0.3 mol fractions of lidocaine.
From the [BMIM]+ COR–[CH3CO2]− COM RDF in Fig. 4, the first local minimum occurs at approximately 0.71 nm, which corresponds to the radius of the first solvation shell. From Fig. 5, we have a coordination number of approximately 5.15, 4.75, and 4.08 at 0.71 nm for the neat IL (x 1 = 0) and with a lidocaine mole fraction of 0.1 (x 1 = 0.1) and 0.3 (x 1 = 0.3), respectively. As compared to the neat IL, this corresponds to a percent difference of 7.8 and 20.8% with a lidocaine mole fraction of 0.1 and 0.3, respectively. As the distance increases, the percent difference continues to increase toward 10% and 30%, the corresponding lidocaine mole fraction. Before looking at lidocaine, we consider [BMIM]+ –[BMIM]+ and [CH3CO2]−–[CH3CO2]− self-interactions in Fig. 6. Using the COR interaction distance, we do not observe ordering of [BMIM]+ with itself. We observe ordering for [CH3CO2]−. However, the first peak in the RDF reaches a height of just g (r) ≈ 1.5, suggesting a relatively weak interaction. 2. Lidocaine
Next, we consider lidocaine in [BMIM]+[CH3CO2]−. We first consider lidocaine–lidocaine self-interactions. The RDF computed with respect to the COR of lidocaine is shown in Fig. 7. We find that self-interactions of lidocaine are relatively weak. This is confirmed in Fig. S4 of the supplementary material42 where we observe a coordination number of 1 at approximately 0.89 and 0.76 nm for
J. Chem. Phys. 144, 084501 (2016)
FIG. 6. Radial distribution function (RDF) for the (self) interaction between [BMIM]+ and itself where the distance is computed with respect to the center-of-ring of [BMIM]+ (top pane) and for the (self) interaction between [CH3CO2]− and itself where the distance is computed with respect to the center-of-mass of [CH3CO2]−. An RDF is shown for the neat IL (x 1 = 0) and the IL plus 0.1 and 0.3 mol fractions of lidocaine.
a lidocaine mole fraction of 0.1 and 0.3, respectively. Lidocaine therefore is not associating with itself suggesting that the studied mole fractions are less than the equilibrium solubility of lidocaine. In Fig. 7, we additionally plot the RDF for lidocaine–[BMIM]+ interactions, with the distance computed with respect to the COR. We find that these interactions are weak and not ordered, suggesting that π–π interactions between the lidocaine and [BMIM]+ ring are not present. We found that the aromatic ring hydrogens of [BMIM]+ formed weak hydrogen bonds with [CH3CO2]− and that at all of the studied lidocaine concentrations we had an associated [BMIM]+[CH3CO2]− pair. Lidocaine is capable of both donating (N—H) and accepting (==O) hydrogen bonds, and we have found that it is extremely soluble in [BMIM]+[CH3CO2]−. Following this, we next consider hydrogen bonding interactions of lidocaine with [BMIM]+ and [CH3CO2]−. The CN and RDF for the lidocaine amine H (N—H) and the [CH3CO2]− oxygen, along with the lidocaine amide O (==O) and the [BMIM]+ aromatic ring hydrogens are shown in Fig. 8. Considering first the interaction between lidocaine and [CH3CO2]−, we observe a high intensity peak in the RDF occurring at r ≈ 0.19 nm with g(r) ≈ 11 and 13 for 0.1 and 0.3 mol fractions of lidocaine, respectively. Immediately after this peak is a sharp decrease to nearly 0. This short, strong, and well defined interaction is characteristic of a moderate hydrogen bond (due mostly to electrostatic interactions).83 For both lidocaine concentrations studied, the peak locations in the
084501-7
R. T. Ley and A. S. Paluch
FIG. 7. Radial distribution function (RDF) for lidocaine · · · · lidocaine and lidocaine · · · · [BMIM]+ intermolecular interactions where the distance is computed with respect to the center-of-ring of lidocaine and [BMIM]+.
RDF do not change. From the CN at both concentrations, we observe a coordination number of 1 at r ≈ 0.29 nm. Therefore, we find that at both lidocaine concentrations, lidocaine is associated with a single [CH3CO2]− via a moderate hydrogen bond. In the RDF, we observe a second peak that is less pronounced yet still has a well defined depletion well at each concentration. This second peak is the result of the second (equivalent) oxygen in the [CH3CO2]− associated with lidocaine and is reflected in the “S” shape of the coordination number between a value of 1 and 2. Next, we consider the possibility of lidocaine forming hydrogen bonds with [BMIM]+. Looking at the RDF for the interaction of the lidocaine amide oxygen (==O) with the aromatic ring hydrogens of [BMIM]+, we observe a peak at r ≈ 0.21 nm with a height of g(r) ≈ 2.8, followed by a well defined depletion well. In Fig. S5 of the supplementary material,42 we show RDFs for each aromatic ring hydrogen of [BMIM]+ and find that the interaction with aromatic ring hydrogen H1 (see Fig. 3) is strongest, with a first peak height of g(r) ≈ 3.5. While weaker than the interaction between the aromatic ring hydrogens and the [CH3CO2]− oxygen, we nonetheless have a well defined electrostatic interaction. From the CN, we observe a coordination number of 1 at approximately 0.35 and 0.38 nm for a lidocaine concentration of 0.1 and 0.3 mole fractions, respectively. The distance at which the coordination number is 2 and 3 is likewise offset and correspond to the two other (intramolecular) aromatic ring hydrogens.
J. Chem. Phys. 144, 084501 (2016)
FIG. 8. CNs (top pane) and RDFs (bottom pane) for the lidocaine amine hydrogen and the [CH3CO2]− oxygen interaction (N—H · · · · O) and for the interaction between the lidocaine amide O and the [BMIM]+ aromatic ring hydrogens (O · · · · C—H). Solid and dashed lines correspond to a lidocaine mole fraction of 0.1 and 0.3, respectively.
3. Solvation mechanism
We now have sufficient information to piece together the underlying molecular details or solvation mechanism. We found that [BMIM]+ associated with a [CH3CO2]− via a weak hydrogen bond. Likewise, lidocaine associated with a single [CH3CO2]− via a moderate hydrogen bond and with a single [BMIM]+ via an electrostatic interaction between the lidocaine amide oxygen and the aromatic ring hydrogens of [BMIM]+. Therefore, it must be the case that lidocaine associates with the [CH3CO2]− from one associated ion pair and with the [BMIM]+ of another (different) ion pair. In this light, let us re-examine the coordination number plot of Fig. 8. Lidocaine associates with a single [BMIM]+. The coordination number is one at the distance corresponding to the electrostatic interaction of the lidocaine amide oxygen and the aromatic ring hydrogens of [BMIM]+. The coordination number is equal to two and three at well defined distances corresponding to the second and third (intramolecular) hydrogen of the associated [BMIM]+ (i.e., on the same aromatic ring). This same [BMIM]+ should be associated with a single [CH3CO2]−. While lidocaine associates with a single [CH3CO2]−, we expect the next closest [CH3CO2]− oxygen to come from the [CH3CO2]− associated with [BMIM]+. From Fig. 8, the coordination number for the lidocaine amine hydrogen with the [CH3CO2]− oxygen is 3 (corresponding to nearest interaction with the second [CH3CO2]−) at a distance of approximately 0.64
084501-8
R. T. Ley and A. S. Paluch
J. Chem. Phys. 144, 084501 (2016)
C. Alchemical transformations
FIG. 9. Spatial distribution function (SDF) for the [CH3CO2]− oxygen (yellow) and the [BMIM]+ aromatic ring hydrogen H1 (red) around lidocaine and a snapshot from the simulations with 0.1 mol fractions of lidocaine.
and 0.66 nm for 0.1 and 0.3 mol fractions of lidocaine, respectively. This is approximately 0.15 and 0.10 nm further than the distance at which the lidocaine amide oxygen interacts with the second and third aromatic ring hydrogen of [BMIM]+ (i.e., where the coordination number is 2 and 3 for the lidocaine oxygen–[BMIM]+ hydrogen interaction). In Fig. 9 we provide the SDF for the [CH3CO2]− oxygen and aromatic ring hydrogen H1 (see Fig. 3) of [BMIM]+ around lidocaine. We observe that the [CH3CO2]− oxygen is highly populated around the lidocaine amine hydrogen, while simultaneously H1 of [BMIM]+ is populated adjacent to the [CH3CO2]− oxygen. Note that the [CH3CO2]− oxygen population density area is relatively large because it includes both [CH3CO2]− oxygens. Additionally, the H1 hydrogen of [BMIM]+ is highly populated around the lidocaine amide oxygen, while simultaneously the [CH3CO2]− oxygen is populated adjacent to the H1 hydrogen of [BMIM]+. In Fig. 9 we additionally provide a snapshot from the simulation trajectory to further illustrate the underlying mechanism. We clearly see that the [CH3CO2]− oxygen interacts with the lidocaine amine hydrogen, while the [BMIM]+ H1 hydrogen interacts with the lidocaine amide oxygen. At the same time, the oxygen of a different [CH3CO2]− interacts with a different aromatic ring hydrogen of the [BMIM]+ associated with lidocaine.
We have found that the favorable solvation of lidocaine in [BMIM]+[CH3CO2]− results from favorable hydrogen bonding between the lidocaine amine hydrogen and the [CH3CO2]− oxygen, favorable electrostatic interactions between the lidocaine amide oxygen with the [BMIM]+ aromatic ring hydrogens, while lidocaine does not interfere with the association of [BMIM]+ with [CH3CO2]−. We next seek to understand the effect of the solute size. We did this by computing the relative solubility of three solutes resulting from removing functional groups from the lidocaine scaffold (see Fig. 1) in [BMIM]+[CH3CO2]− relative to water. In all cases, the amide group is maintained so that we expect to maintain the favorable interactions between the solute with [BMIM]+ and [CH3CO2]−. The removal of functional groups decreases the solute size. Note that we have only decreased the solute size in the present study to ensure adequate solute conformational sampling. In each case, the solubility of the modified solute in the IL relative to water decreased. The solubility of the solutes in the IL relative to water from lowest to highest follows the trend: mod-lido C < mod-lido B < mod-lido A < lidocaine. We find that the relative solubility increases as the van der Waals volume increases. We therefore expect the performance of [BMIM]+[CH3CO2]− to continue to increase as the solute size increases. Of course, there must be an upper bound on the size of the solute at which point it begins to perturb the IL structure. A detailed study of solute size is reserved for a future study. Additionally, the free energy calculations of lidocaine in water were repeated with the TIP4P/200584 and TIP4P-Ew85 models to confirm that our choice of water model did not influence the findings of the present study. These and the results of all of the free energy calculations for all of the studied solutes are tabulated in the supplementary material of this manuscript.42
V. SUMMARY AND CONCLUSIONS
Using molecular simulation free energy calculations for a single lidocaine molecule in solution, the solubility of lidocaine in seven common room temperature ionic liquids was predicted to be 103 to 105 orders of magnitude greater than in water. While the infinite dilution approximation would break down at such large concentrations and the actual solubility would be smaller than predicted, the results none-the-less suggest that the solvation of lidocaine in the studied ionic liquids is highly favorable. Even for the case of the hydrophobic ionic liquids, the solubility is greatly enhanced relative to water. If the main reason for the favorable solubility was hydrogen bonding between lidocaine and the ionic liquid, we would expect the solubility to increase with increasing ionic liquid hydrogen bond basicity. We found this not to be the case for lidocaine, suggesting that hydrogen bonding alone does not explain the underlying mechanism. For the most favorable case, [BMIM]+[CH3CO2]−, long trajectories were computed with a lidocaine mole fraction of 0.1 and 0.3 to study solute-solvent and solvent-solvent interactions.
084501-9
R. T. Ley and A. S. Paluch
We found that the favorable solvation of lidocaine in [BMIM]+[CH3CO2]− results from favorable hydrogen bonding between the lidocaine amine hydrogen and the [CH3CO2]− oxygen, favorable electrostatic interactions between the lidocaine amide oxygen with the [BMIM]+ aromatic ring hydrogens, while lidocaine does not interfere with the association of [BMIM]+ with [CH3CO2]−. Essentially, lidocaine is able to interact favorable with two separate associated ion pairs. Therefore, while the ability of the anion to accept hydrogen bonds (which is related to the ionic liquid hydrogen bond basicity) is important, so too are the interactions between the ionic liquid cation and anion. Maintaining the amide group of lidocaine responsible for the favorable interactions with [BMIM]+ and [CH3CO2]−, three additional solutes were studied wherein functional groups were removed from the lidocaine scaffold. This had the effect of decreasing the solute van der Waals volume. We found that as the van der Waals volume decreased, the solubility in [BMIM]+[CH3CO2]− relative to water decreased. This suggests that as the solute size increases, IL performance should increase. At some point we would expect this trend to breakdown, but that is subject to another study. The underlying mechanism is very interesting and promising for future investigations. The favorable solvation of lidocaine in ionic liquids, including hydrophobic ionic liquids, has promising applications for the development of improved separation processes for pharmaceutical compounds. Additionally, ionic liquid use is promising for controlled release and delivery applications. ACKNOWLEDGMENTS
R.T.L. acknowledges summer financial support from the Miami University College of Engineering and Computing. A.S.P. gratefully acknowledges funding from the Miami University Committee on Faculty Research and start-up support from the Miami University College of Engineering and Computing. Computing support was provided by the Ohio Supercomputer Center and Miami University’s Research Computing Support group. 1C.
A. Angell, Y. Ansari, and Z. Zhao, Faraday Discuss. 154, 9 (2012). F. Brennecke and E. J. Maginn, AIChE J. 47, 2384 (2001). 3W. L. Hough, M. Smiglak, H. Rodríguez, R. P. Swatloski, S. K. Spear, D. T. Daly, J. Pernak, J. E. Grisel, R. D. Carliss, M. D. Soutullo et al., New J. Chem. 31, 1429 (2007). 4I. M. Marrucho, L. C. Branco, and L. P. N. Rebelo, Annu. Rev. Chem. Biomol. Eng. 5, 527 (2014). 5K. B. Smith, R. H. Bridson, and G. A. Leeke, J. Chem. Eng. Data 56, 2039 (2011). 6A. Forte, C. I. Melo, R. Bogel-Lukasik, and E. Bogel-Lukasik, Fluid Phase Equilib. 318, 89 (2012). 7A. D. dos Santos, A. R. C. Morais, C. I. Melo, R. Bogel-Lukasik, and E. Bogel-Lukasik, Fluid Phase Equilib. 356, 18 (2013). 8C. C. Weber, A. J. Kunov-Kruse, R. D. Rogers, and A. S. Myerson, Chem. Commun. 51, 4294 (2015). 9H. J. Park and M. R. Prausnitz, AIChE J. 61, 2732 (2015). 10Y. C. Kau, C. C. Liao, Y. C. Chen, and S. J. Liu, Materials 7, 6660 (2014). 11D. E. Ickowicz, L. Golovanevski, A. J. Domb, and C. F. Weiniger, Int. J. Pharm. 468, 152 (2014). 12A. S. Paluch, S. Parameswaran, S. Liu, A. Kolavennu, and D. L. Mobley, J. Chem. Phys. 142, 044508 (2015). 13A. S. Paluch, C. A. Vitter, J. K. Shah, and E. J. Maginn, J. Chem. Phys. 137, 184504 (2012). 2J.
J. Chem. Phys. 144, 084501 (2016) 14A.
R. Leach, Molecular Modelling: Principles and Applications, 2nd ed. (Pearson Education Limited, Harlow, England, 2001). 15D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. (Academic Press, San Diego, CA, 2002). 16X. Zhong, Z. Liu, and D. Cao, J. Phys. Chem. B 115, 10027 (2011). 17W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). 18B. Chen and J. I. Siepmann, J. Phys. Chem. B 110, 3555 (2006). 19J. L. Rafferty, L. Sun, J. I. Siepmann, and M. R. Schure, Fluid Phase Equilib. 290, 25 (2010). 20N. Rai and J. I. Siepmann, J. Phys. Chem. B 111, 10790 (2007). 21N. Rai, D. Bhatt, J. I. Siepmann, and L. E. Fried, J. Chem. Phys. 129, 194510 (2008). 22N. Rai and J. I. Siepmann, J. Phys. Chem. B 117, 273 (2013). 23C. D. Wick, J. M. Stubbs, N. Rai, and J. I. Siepmann, J. Phys. Chem. B 109, 18974 (2005). 24Note that TraPPE-EH does parameterize aniline in Ref. 21. We chose not to use the N LJ parameters from aniline because it is a primary amine, and it was previously shown in Ref. 23 that there is an appreciable change in LJ parameters in going from a primary to secondary amine. Also, in Ref. 23, the primary amide LJ parameters are the same as the primary amine (i.e., only the charges change). Therefore, we took the secondary amide LJ parameters for N to be the same as for a secondary amine. The H has LJ parameters of 0 in primary and secondary amines and primary amides. The LJ parameters for H in a secondary amide were therefore taken to also be 0 in the present study. 25J. M. Stubbs, J. J. Potoff, and J. I. Siepmann, J. Phys. Chem. B 108, 17596 (2004). 26M. G. Martin and J. I. Siepmann, J. Phys. Chem. B 102, 2569 (1998). 27The alkyl groups (CH and CH ) were modeled as a single united-atom 3 2 pseudoatom as a result of the parameterization of the TraPPE-EH force field for n-alkanes which places the LJ site for a hydrogen atom at the center of the corresponding bond,86 and the complication of implementing such a model in a molecular dynamics framework. 28Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008). 29C. J. Cramer, Essentials of Computational Chemistry (John Wiley & Sons Ltd., Chichester, West Sussex, England, 2002). 30A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. Phys. Chem. B 113, 6378 (2009). 31A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. Phys. Chem. B 113, 4538 (2009). 32R. F. Ribeiro, A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. Comput.Aided Mol. Des. 24, 317 (2010). 33M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson et al., 09, Revision A.02, Gaussian, Inc., Wallingford, CT, 2009. 34C. I. Bayly, P. Cieplak, W. D. Cornell, and P. A. Kollman, J. Phys. Chem. 97, 10269 (1993). 35P. Cieplak, W. D. Cornell, C. Bayly, and P. A. Kollman, J. Comput. Chem. 16, 1357 (1995). 36D. A. Case, T. Cheatham, T. Darden, H. Gohlke, R. Luo, K. M. Merz, A. Onufriev, C. Simmerling, B. Wang, and R. Woods, J. Comput. Chem. 26, 1668 (2005). 37D. A. Case, T. A. Darden, T. E. Cheatham III, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, R. C. Walker, W. Zhang, K. M. Merz et al., Amber 12 (University of California San Francisco, San Francisco, CA, 2012), p. 12. 38J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, J. Comput. Chem. 25, 1157 (2004). 39J. Wang, W. Wang, P. A. Kollman, and D. A. Case, J. Mol. Graphics Modell. 25, 247 (2006). 40A. W. Sousa da Silva and W. F. Vranken, BMC Res. Notes 5, 367 (2012). 41A. W. Sousa da Silva and W. F. Vranken, acpype: AnteChamber PYthon Parser interfacE, accessed May 1, 2014, http://www.gromacs.org/. 42See supplementary material at http://dx.doi.org/10.1063/1.4942025 for GROMACS force field files (type “itp” and “xyz”) for the studied ionic liquids and solutes are provided in gromacs_itp.txt and gromacs_xyz.txt; additional SDFs, RDFs, and coordination numbers are provided to facilitate discussion along with tabulated results of the free energy calculations are provided in supp_info.pdf. 43L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, J. Comput. Chem. 30, 2157 (2009). 44See http://www.ime.unicamp.br/∼martinez/packmol/ for Packmol: Packing Optimization for Molecular Dynamics Simulations, accessed August 1, 2013.
084501-10 45D.
R. T. Ley and A. S. Paluch
van der Spoel, E. Lindahl, B. Hess, and the GROMACS development team, GROMACS User Manual, version 4.6.3, 2013, ftp://ftp.gromacs.org/ pub/manual/manual-4.6.3.pdf. 46H. J. C. Berendsen, Simulating the Physical World: Hierarchial Modeling from Quantum Mechanics to Fluid Dynamics (Cambridge University Press, New York, NY, 2007). 47G. Bussi, D. Donadio, and M. Parrinello, J. Chem. Phys. 126, 014101 (2007). 48G. Bussi and M. Parrinello, Comput. Phys. Commun. 179, 26 (2008). 49G. Bussi, T. Zykova-Timan, and M. Parrinello, J. Chem. Phys. 130, 074101 (2009). 50H. J. C. Berendsen, J. P. M. Postma, A. DiNola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984). 51G. J. Martyna, D. J. Tobias, and M. L. Klein, J. Chem. Phys. 101, 4177 (1994). 52G. J. Martyna, M. L. Klein, and M. E. Tuckerman, J. Chem. Phys. 97, 2635 (1992). 53B. Hess, C. Kutzner, D. van der Spoel, and E. Lindal, J. Chem. Theory Comput. 4, 435 (2008). 54S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel et al., Bioinformatics 29, 845 (2013). 55See http://www.gromacs.org/ for GROMACS: Fast, flexible, free, accessed August 1, 2013. 56B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, J. Comput. Chem. 18, 1463 (1997). 57B. Hess, J. Chem. Theory Comput. 4, 116 (2008). 58J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comput. Phys. 23, 327 (1977). 59H. C. Andersen, J. Comput. Phys. 52, 24 (1983). 60M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 (1998). 61G. B. Fuerst, R. T. Ley, and A. S. Paluch, Ind. Eng. Chem. Res. 54, 9027 (2015). 62K. S. Shing and S. T. Chung, J. Phys. Chem. 91, 1674 (1987). 63D. A. Kofke and P. T. Cummings, Mol. Phys. 92, 973 (1997). 64M. R. Shirts, J. W. Pitera, W. C. Swope, and V. S. Pande, J. Chem. Phys. 119, 5740 (2003). 65D. A. Kofke and P. T. Cummings, Fluid Phase Equilib. 150-151, 41 (1998). 66Free Energy Calculations: Theory and Applications in Chemistry and Biology, Springer Series in Chemical Physics Vol. 86, edited by C. Chipot and A. Pohorille (Springer, New York, NY, 2007). 67C. H. Bennett, J. Comput. Phys. 22, 245 (1976). 68M. R. Shirts, E. Bair, G. Hooker, and V. S. Pande, Phys. Rev. Lett. 91, 140601 (2003). 69N. Lu, J. K. Singh, and D. A. Kofke, J. Chem. Phys. 118, 2977 (2003). 70M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129, 124105 (2008). 71T. C. Beutler, A. E. Mark, R. C. van Schaik, P. R. Gerber, and W. F. van Gunsteren, Chem. Phys. Lett. 222, 529 (1994).
J. Chem. Phys. 144, 084501 (2016) 72M.
R. Shirts and V. S. Pande, J. Chem. Phys. 122, 134508 (2005). Steinbrecher, D. L. Mobley, and D. A. Case, J. Chem. Phys. 127, 214108 (2007). 74See https://github.com/choderalab/pymbar for PyMBAR: Python implementation of the Multistate Bennett Acceptance Ratio (MBAR), accessed May 1, 2014. 75J. D. Chodera, W. C. Swope, J. W. Pitera, C. Seok, and K. A. Dill, J. Chem. Theory Comput. 3, 26 (2007). 76P. V. Klimovich, M. R. Shirts, and D. L. Mobley, J. Comput.-Aided Mol. Des. 29, 397 (2015). 77In summary, from m = 0 to 14, we have λ LJ = {0.0, 0.1, 0.2, 0.3, 0.4, m 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0} and λ elec m = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.50, 0.71, 0.87, 1.00}. A detailed discussion of best practices when coupling/decoupling intermolecular interactions during free energy calculations is provided in the recent review by Shirts and Mobley.87 78M. Brehm and B. Kirchner, J. Chem. Inf. Model. 51, 2007 (2011). 79B. Kirchner, TRAVIS–Trajectory Analyzer and Visualizer, accessed November 13, 2015, http://www.travis-analyzer.de/. 80The conformation of the central molecule was generated by TRAVIS from the initial system configuration. The SDF was visualized using VMD88,89 with the central molecule represented as a CPK model and the second species or interaction site was modeled as an “isosurface.” 81Note, the experimental mole fraction solubility of lidocaine in water at 298.15 K is 3 × 10−4.90 Lidocaine in the ILs at the large mole fraction solubilities predicted here would be beyond the infinite dilution limit, and we would expect the actual solubility to be less than predicted. Regardless, it stands that there are extremely favorable interactions present in the ILs leading to extremely favorable solubilities over water. 82A. F. M. Cláudio, L. Swift, J. P. Hallett, T. Welton, J. A. P. Coutinho, and M. G. Freire, Phys. Chem. Chem. Phys. 16, 6593 (2014). 83G. A. Jeffrey, An Introduction to Hydrogen Bonding (Oxford University Press, New York, NY, 1997). 84J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005). 85H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys. 120, 9665 (2004). 86B. Chen and J. I. Siepmann, J. Phys. Chem. B 103, 5370 (1999). 87M. R. Shirts and D. L. Mobley, “An introduction to best practices in free energy calculations,” in Biomolecular Simulations: Methods and Protocols (Humana Press, Totowa, NJ, 2013), pp. 271–311, ISBN: 978-1-62703017-5, http://dx.doi.org/10.1007/978-1-62703-017-5_11. 88W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics 14, 33 (1996). 89See http://www.ks.uiuc.edu/Research/vmd/ for VMD: Visual Molecular Dynamics, accessed November 25, 2015. 90Solubility and Related Properties of Large Complex Chemicals Part 1: Organic Solutes Ranging from C4 to C40, edited by J. Marrero and J. Abildskov (DECHEMA, Frankfurt am Main, Germany, 2003). 73T.
Citation: J. Chem. Phys. 144, 084501 (2016); doi: 10.1063/1.4942025 View online: http://dx.doi.org/10.1063/1.4942025 View Table of Contents: http://aip.scitation.org/toc/jcp/144/8 Published by the American Institute of Physics
THE JOURNAL OF CHEMICAL PHYSICS 144, 084501 (2016)
Understanding the large solubility of lidocaine in 1-n-butyl-3-methylimidazolium based ionic liquids using molecular simulation Ryan T. Ley and Andrew S. Paluch Department of Chemical, Paper and Biomedical Engineering, Miami University, Oxford, Ohio 45056, USA
(Received 3 December 2015; accepted 3 February 2016; published online 23 February 2016) Room temperature ionic liquids have been proposed as replacement solvents in a wide range of industrial separation processes. Here, we focus on the use of ionic liquids as solvents for the pharmaceutical compound lidocaine. We show that the solubility of lidocaine in seven common 1-n-butyl-3-methylimidazolium based ionic liquids is greatly enhanced relative to water. The predicted solubility is greatest in [BMIM]+[CH3CO2]−, which we find results from favorable hydrogen bonding between the lidocaine amine hydrogen and the [CH3CO2]− oxygen, favorable electrostatic interactions between the lidocaine amide oxygen with the [BMIM]+ aromatic ring hydrogens, while lidocaine does not interfere with the association of [BMIM]+ with [CH3CO2]−. Additionally, by removing functional groups from the lidocaine scaffold while maintaining the important amide group, we found that as the van der Waals volume increases, solubility in [BMIM]+[CH3CO2]− relative to water increases. C 2016 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4942025]
I. INTRODUCTION
Room temperature ionic liquids (ILs) have received tremendous attention for use as replacement solvents in a wide range of industrial processes.1,2 The cation and anion may be changed and modified to tune their chemical and physical properties. Additionally, due to their ionic nature, ILs have negligible vapor pressures, minimizing solvent losses. Recently, the application of ILs for biological processes has emerged.3,4 One application has been to use ILs as solvents in crystallization processes.5–8 Weber et al.8 found the equilibrium solubility of acetaminophen in the IL 1-ethyl3-methylimidazolium acetate to be very large due to hydrogen bonding between acetaminophen (hydrogen bond donor) and acetate (hydrogen bond acceptor). They found that the solubility of acetaminophen could therefore be tuned with the addition of a hydrogen bond donating co-solvent (an “antisolvent”) that would compete with acetaminophen to hydrogen bond with the acetate anion. In this fashion efficient crystallization processes could be developed. In addition, the development of biologically active ILs has been proposed.3,9 In addition to advantages of using an IL over a conventional liquid, biologically active ILs are advantageous as one can obtain a high concentration of the active ingredient(s). Lidocaine is an important local anesthetic typically applied topically. A great deal of research is devoted to developing advanced controlled topical delivery methods.10,11 It has been shown that similar results may be obtained using a lidocaine cation based IL.3,9 These lidocaine based ILs form a gel that may readily be applied topically with a high lidocaine concentration. Following Refs. 8 and 9, here we study the solubility of (neutral) lidocaine in seven common ILs with the 0021-9606/2016/144(8)/084501/10/$30.00
1-n-butyl-3-methylimidazolium cation (see Figs. 1 and 2). Using molecular simulation free energy calculations, we predict that in all seven ILs, the solubility is greatly increased relative to water, with enhancements ranging from 103 to 105. Like acetaminophen, lidocaine contains an amide group capable of donating hydrogen bonds. However, as compared to Ref. 8, upon detailed investigation using molecular structural analysis, we find that while the ability of lidocaine to hydrogen bond with the anion is important, so too are the interactions between the ionic liquid cation and anion. For the case of lidocaine in [BMIM]+[CH3CO2]−, we find that the favorable solvation results from favorable hydrogen bonding between the lidocaine amine hydrogen and the [CH3CO2]− oxygen, favorable electrostatic interactions between the lidocaine amide oxygen with the [BMIM]+ aromatic ring hydrogens, while lidocaine does not interfere with the association of [BMIM]+ with [CH3CO2]−. Essentially, lidocaine is able to interact favorable with two separate associated ion pairs. The result is interesting for lidocaine given its large size and chemical complexity. Additionally, maintaining the amide group of lidocaine responsible for the favorable interactions with [BMIM]+ and [CH3CO2]−, we considered three additional solutes wherein functional groups were removed from the lidocaine scaffold. This had the effect of decreasing the solute van der Waals volume. We found that as the van der Waals volume increased, the solubility in [BMIM]+[CH3CO2]− relative to water increased. The application of ILs as solvents in controlled release applications, such as IL supported membranes or ILs separated from skin by a semipermeable membrane are extremely promising. In addition, the results demonstrate the ability to use ILs for crystallization and the extraction of lidocaine from aqueous environments.
144, 084501-1
© 2016 AIP Publishing LLC
084501-2
R. T. Ley and A. S. Paluch
J. Chem. Phys. 144, 084501 (2016)
FIG. 1. The chemical structure of the studied solutes consisting of lidocaine (CAS 137-58-6) and three solutes resulting from removing functional groups from the lidocaine scaffold.
achieved in a simulation by having a single solute molecule in a system with a large number of solvent molecules and offers the advantage of isolating solute-solvent intermolecular interactions. The present study focuses on the pharmaceutical solute lidocaine (see Fig. 1) in the seven ILs composed of the 1-n-butyl-3-methyl-imidazolium cation ([BMIM]+) with one of the seven anions shown in Fig. 2. As we will show, we find that the solubility of lidocaine in the IL [BMIM]+[CH3CO2]− is 105 times greater than in water, which is larger than for any of the other studied ILs. The mechanism responsible for this favorable solvation of lidocaine in [BMIM]+[CH3CO2]− is elucidated by performing structural analyses and additional free energy calculations. To facilitate the structural analysis, the concentration of lidocaine in [BMIM]+[CH3CO2]− was increased to a mole fraction of approximately 0.1 and 0.3. Next, solubilities in [BMIM]+[CH3CO2]− relative to water were computed for three additional solutes wherein functional groups were systematically removed from the lidocaine scaffold. The chemical structure of the additional solutes is provided in Fig. 1.
III. COMPUTATIONAL DETAILS A. Force fields
FIG. 2. The chemical structure of the seven anions and the 1-n-butyl-3methylimidazolium cation ([BMIM]+) which constitute the seven studied ionic liquids.
II. METHODOLOGY
The solubility of a solid solute (component 1) in an IL (solvent i) relative to water may be described by the classical equations of phase-equilibria thermodynamics. Assuming the solute concentration is sufficiently small so as to be considered infinitely dilute, the relative solubility may be computed using molecular simulation free energy calculations as12,13 c1,i = β µres,∞ ln 1,water (T, p, N1 = 1, Nwater) c1,water − β µres,∞ 1,i (T, p, N1 = 1, Ni ) ,
(1)
where c1,i and c1,water are the molar or mass concentrations of the solute in solvent i and water at equilibrium, respectively, β µres,∞ and β µres,∞ 1,i 1,water are the dimensionless residual chemical potential of the solute infinitely dilute in solvent i and water, respectively, where β −1 = k BT and k B is the Boltzmann constant, T is the temperature, p is the pressure, and N1, Ni , and Nwater are the number of solute, solvent i, and water molecules in the simulation. The infinite dilution limit is
Interactions were modeled using a “class I” potential energy function where all non-bonded intermolecular interactions (Unb) were accounted for using a combined Lennard-Jones (LJ) plus fixed point charge model of the form14,15 ( σ ) 12 ( σ ) 6 1 qi q j ij ij − + , (2) Unb r i j = 4ε i j r i j 4πε 0 r i j r i j where r i j is the separation distance between sites i and j, ε i j is the well-depth of the LJ potential, σi j is the distance at which the LJ potential is zero, and qi and q j are the partial charges of sites i and j, respectively. The united atom force field of Zhong et al.16 was used to model all of the ILs and water was modeled with TIP4P,17 which has been shown to work well using TraPPE force field models.18,19 The force field for lidocaine and the studied solutes wherein functional groups were removed from the lidocaine scaffold were all constructed based on the Explicit Hydrogen Transferable Potentials for Phase Equilibria (TraPPE-EH) force field.20–22,86 The LJ parameters for the phenyl group were taken directly from TraPPE-EH.20 The LJ parameters for N and H in the secondary amide were taken to be the same as for N and H in a secondary amine, and parameters for the primary and tertiary amine were all taken from TraPPE-EH.23,24 The carbonyl group LJ parameters were taken from the TraPPE-UA amide model,23,25 and all of the alkyl (CH3 and CH2) groups were modeled as a single alkyl pseudoatom with TraPPE-UA LJ parameters.26,27 Partial atomic charges were obtained for the TraPPEEH based models in a similar fashion as the original TraPPE-EH parameterization.20–22 Slight differences result due to the capabilities of the electronic structure software package available to us at the time of this study and the
084501-3
R. T. Ley and A. S. Paluch
use of united-atom alkyl sites. Charges were parameterized following a three step procedure. First, the gas-phase structure for each solute was optimized at the M06-2X/cc-pVTZ level of theory/basis set.28,29 Second, a single point energy calculation was performed on the gas-phase optimized structure at the M06-2X/6-31G(d) level of theory/basis set28,29 in a self-consistent reaction field (SCRF) using the SMD universal solvation model for 1-octanol.30–32 All of the electronic structure calculations were performed using Gaussian 09.33 Third, partial atomic charges were then obtained from the electrostatic potential (obtained in step 2) using the restrained electrostatic potential (RESP)34,35 method in ANTECHAMBER (part of the AMBER 12 simulation suite).36,37 During the fitting procedure, the partial charges of alkyl hydrogens were fixed to zero; this had the effect of using a single charge site at the center of the alkyl group appropriate for use with a united atom model. All of the intramolecular parameters for the solutes were taken from the General AMBER Force Field (GAFF).38,39 Parameters were generated using ANTECHAMBER and converted from AMBER to GROMACS format using ACPYPE.40,41 Since alkyl groups were modeled as a single alkyl pseudoatom, intramolecular parameters involving alkyl hydrogens were removed. Throughout the present study, solute and IL bonds involving hydrogens were held fixed. All of the GROMACS force field files used in the present study are provided in the supplementary material of this manuscript.42 B. Molecular dynamics 1. Ionic liquids
As a result of their sluggish dynamics, extensive equilibrations were performed for the IL systems. Simulations were performed for neat (pure) ILs, ILs plus a single lidocaine molecule, which corresponds to infinite dilution, [BMIM]+[CH3CO2]− with a lidocaine mole fraction of approximately 0.1 and 0.3, and [BMIM]+[CH3CO2]− plus a single solute molecule which results from systematically removing functional groups from lidocaine. All of these simulations were carried out following the same procedure. First, Packmol was used to generate initial structures.43,44 The number of IL pairs for the neat IL simulations were chosen to give a cubic box length of approximately 4.5 nm at 298.15 K based on the simulation results of Ref. 16. The number of IL pairs remained the same for systems containing solutes and the box size was adjusted accordingly. Subsequently, between 3000 and 6000 steepest descent minimization steps were performed to remove any bad contacts that might have resulted from the packing. The next two steps involved dynamics with the equations of motion integrated using the Verlet leap-frog algorithm.14,15,45,46 The system was first equilibrated in an NVT ensemble at 738.15 K for 1 ns using stochastic velocity rescaling.45,47–49 Next, the system was brought from 738.15 K to 298.15 K and 1 bar using simulated annealing using stochastic velocity rescaling and the Berendsen barostat.45,46,50
J. Chem. Phys. 144, 084501 (2016)
The system was initially carried out for 1 ns at 738.15 K and 1 bar. Over the next 8.8 ns, the temperature was decreased linearly at a rate of 50 K/ns to reach the final temperature of 298.15 K and was then held at 298.15 K for 1 ns. Finally, a 13 ns NpT simulation was performed at 298.15 K and 1 bar with the equations of motion integrated using the velocity Verlet algorithm.14,15,45,46 The simulation used the Nosé-Hoover chain thermostat (with a chain length of 10) and the Martyna-Tuckerman-TobiasKlein (MTTK) barostat.15,45,51,52 To ensure accurate statistics when performing the structural analyses, the simulations for neat [BMIM]+[CH3CO2]− and [BMIM]+[CH3CO2]− plus approximately 0.1 and 0.3 mol fraction lidocaine were carried out for an additional 101 ns. For all three steps involving dynamics, the time constant for the thermostat was 1 ps and the time constant for the barostat was 4 ps. Free energy calculations in the NpT ensemble were conducted for systems with a single solute molecule in the IL, which will be elaborated upon momentarily. 2. Water
The initial structure for water was taken from a previously equilibrated system at ambient conditions. A single solute molecule (lidocaine or a solute resulting from removing functional groups from the lidocaine scaffold) was solvated in these systems in a cubic box approximately 4.5 nm in length using the GROMACS tool “genbox.” The system was then equilibrated for 2.5 ns using a Nosé-Hoover chain thermostat (with a chain length of 10) and the MTTK barostat with a time constant of 1 and 4 ps for the thermostat and barostat, respectively, using the velocity Verlet algorithm. Free energy calculations in the NpT ensemble were then conducted which will be elaborated upon momentarily. All of the molecular dynamics simulations in this study were performed using GROMACS 4.6.3.53–55 Atom bond lengths involving a hydrogen were constrained using P-LINCS45,56,57 with the leap-frog Verlet algorithm and using SHAKE45,58,59 with the velocity Verlet algorithm. The Verlet neighbor list was used45 and LJ interactions were cutoff at 1.4 nm. Long-range analytic dispersion corrections were applied to the energy and pressure to accommodate the truncation.14,15,45,46 Lorentz-Berthelot mixing rules were used for unlike LJ sites.14 The electrostatic terms were evaluated with the smooth particle-mesh-Ewald (SPME) method with tin-foil boundary conditions45,46,60 with real space interactions truncated at 1.4 nm. The SPME B-spline was order 4, the Fourier spacing was 0.12 nm, and the relative tolerance between long and short-range energies was 10−4. 3. Free energy calculations
The free energy calculations were performed at 298.15 K and 1 bar using a similar protocol as our recent work.61 The infinite dilution residual chemical potential, µres,∞ 1,i , of the solute in the ILs and water was calculated using a multi-stage free energy perturbation method62–66 with the multi-state Bennett’s acceptance ratio (MBAR) method.67–70,87
084501-4
R. T. Ley and A. S. Paluch
J. Chem. Phys. 144, 084501 (2016)
A “soft-core” potential was used to decouple the solutesolvent intermolecular LJ interactions. Stage (m) dependent elec decoupling parameters, λ LJ controlled the LJ m and λ m and electrostatic intermolecular interactions, respectively.
The decoupling parameters varied from 0 to 1. When elec λ LJ m = λ m = 1, the solute is fully coupled to the system. When elec λ LJ = λ m m = 0, the solute is decoupled from the system. The “soft-core” potential had the form71–73,87
σi6j σi12j sc , − ULJ r i j ; m = 4λ LJ ε m ij 2 LJ 6 6 1 − λ LJ α σ 6 + r 6 1 − λ α σ + r LJ m ij ij LJ i j m ij
where αLJ is a constant, which had a value of 1/2. The advantage of using a “soft-core” potential to decouple the LJ interactions is that while it yields the correct limiting value of the potential (when λ LJ m = 0 and 1), it additionally allows nearly decoupled molecules to overlap with a finite energy (and hence finite probability). The electrostatic term in the intermolecular potential was decoupled linearly as 1 qi q j . Uelec r i j ; m = λ elec m 4πε 0 r i j
(4)
At each stage m, an independent MD simulation was performed. The simulation time for each stage m was 18.5 ns, where the first 0.5 ns was discarded from analysis as equilibration. The initial structure for these simulations was taken as the final structure from the NpT simulation of a single solute in solution (either IL or water). The change in the Hamiltonian with the current configuration between stage m and the other stages is computed every 0.20 ps. This is saved for subsequent post-simulation analysis with MBAR70 to determine µres,∞ 1,i . This analysis was performed using the Python implementation of MBAR (PyMBAR) and the GROMACS analysis script distributed with it.74–76 A total of 15 different stages were used for the free energy calculations where m = 0 corresponds to a non-interacting (ideal gas) state and m = 14 is a fully interacting system. From m = 1 to 10, the LJ interactions were increased from λ LJ m = 0.1 to 1.0 in 10 equal increments of 0.1. Electrostatic interactions were increased in a square root fashion following 77 λ elec m = {0.50, 0.71, 0.87, 1.00} from m = 11 to 14. The simulation parameters for the free energy calculations were the same as the last step of equilibration (i.e., using velocity Verlet with a Nosé-Hoover chain thermostat and MTTK barostat) with one exception: the efficient Verlet cut-off scheme is not available for use with free energy calculations in GROMACS 4.6.3. Instead, the LJ interactions were switched off smoothly between 1.15 and 1.2 nm, and the real space part of the electrostatic interactions was switched off smoothly between 1.18 and 1.2 nm. Long range corrections were applied as described previously. C. Structural analysis
As will be shown in Section IV, lidocaine was found to have the greatest affinity for the IL [BMIM]+[CH3CO2]−. To understand the underlying mechanism, detailed structural analysis was performed. Specifically, radial distribution
(3)
function (RDF) and spatial distribution function (SDF) calculations were performed for both the neat IL and the IL with approximately 0.1 and 0.3 mol fraction lidocaine. In both cases, the number of IL pairs was constant and 273, and the number of lidocaine molecules was 30 in the 0.1 mole fraction case and 90 in the 0.3 mol fraction case. To ensure accurate statistics, the RDFs and SDFs were computed from long 101 ns trajectories. The RDFs were computed using the GROMACS tool “g_rdf” over the last 100 ns of the trajectory. The SDFs were computed using the TRAVIS software.78–80,88,89 IV. RESULTS AND DISCUSSION A. Free energy calculations
First, we predict the solubility of lidocaine in the seven ILs relative to water. The ILs all consisted of a [BMIM]+ cation paired with one of seven common anions (see Fig. 2). In the order of least to most soluble, the anions are ranked as [BF4]− < [CH3SO4]− < [CF3SO3]− = [PF6]− < [TF2N]− < [CF3CO2]− < [CH3CO2]−. The corresponding values of ln ci /cwater are, respectively, {7.64, 8.55, 9.95, 9.94, 10.24, 10.63, 11.84} where the subscript corresponds to the uncertainty in the last digit. In all cases, the solubility of lidocaine in the ILs is larger than in water, with enhancements ranging from 103 to 105. This result is extremely encouraging. Even for the case of the hydrophobic ILs, [BMIM]+[BF4]−, [BMIM]+[PF6]−, and [BMIM]+[TF2N]−, the solubility relative to water is very large. These three ILs would be extremely promising for the extraction of lidocaine from aqueous environments. The largest solubility observed was for lidocaine in [BMIM]+[CH3CO2]−. For this case, the solubility was approximately 105 orders of magnitude greater than in water.81,90 Recently, Weber et al.8 demonstrated how mixtures of [EMIM]+[CH3CO2]− with a hydrogen bond donating co-solvent can be used to develop advanced separation processes. In that study, they found the solubility of acetaminophen in a range of ILs increased with increasing IL hydrogen bond basicity.82 The hydrogen bond basicity corresponds to the ILs ability to accept hydrogen bonds. Similar to acetaminophen, lidocaine also contains an amide group which can act as a hydrogen bond donor (N—H), so we might expect similar results here. If hydrogen bonding between lidocaine and the IL was the
084501-5
R. T. Ley and A. S. Paluch
J. Chem. Phys. 144, 084501 (2016)
dominate (or determining) interaction present, we would predict the trend [PF6]− < [TF2N]− < [BF4]− < [CF3SO3]− < [CH3SO4]− < [CF3CO2]− < [CH3CO2]−.82 However, this ranking differs from that observed. Therefore, while solubility is the greatest in [BMIM]+[CH3CO2]−, this suggests that it does not solely result from hydrogen bonding between lidocaine and the IL (anion). There must be additional molecular level driving forces responsible for the favorable solubility of lidocaine. To elucidate the underlying driving forces, we will next perform detailed structural analyses for neat [BMIM]+[CH3CO2]− and for systems of [BMIM]+[CH3CO2]− with approximately 0.1 and 0.3 mol fractions of lidocaine. B. Structural analysis
The long 101 ns simulations of neat [BMIM]+[CH3CO2]− and [BMIM]+[CH3CO2]− with 0.1 and 0.3 mol fractions of lidocaine were used to compute RDFs and SDFs to investigate the underlying solvation mechanism. 1. Ionic liquid
First we consider the ionic liquid structure. RDFs associated with the interaction of [BMIM]+ with [CH3CO2]− are shown in Fig. 4 and are computed both with respect to the center-of-ring (COR) of [BMIM]+ and the center-ofmass (COM) of [CH3CO2]−, and between the aromatic ring hydrogens of [BMIM]+ and the [CH3CO2]− oxygens. We find that [BMIM]+ and [CH3CO2]− are well ordered at short distances. The dominant interaction is between the aromatic hydrogens of [BMIM]+ with the [CH3CO2]− oxygens. We observe a high intensity peak in the RDF between 0.2 and 0.21 nm which decreases sharply. This short, well defined interaction is characteristic of a weak hydrogen bond (due to electrostatic interactions).83 In the supplementary material (Fig. S2),42 we provide RDFs for each aromatic hydrogen of [BMIM]+ with the [CH3CO2]− oxygens. These RDFs suggest that the strongest interaction is with aromatic ring hydrogen H1 (see Fig. 3), although the difference with the other hydrogens is rather small. After the first peak, g (r) remains close to one. Additionally, in Fig. 4, we find that the RDFs between [BMIM]+ and [CH3CO2]− change little with increasing lidocaine concentration. The peak locations all remain the same, with the intensity of the first peak increasing slightly with increasing concentration. SDFs for [CH3CO2]− around [BMIM]+ in the neat IL and with a lidocaine concentration of 0.1 and 0.3 mol fractions are provided in Fig. S3 of the supplementary material.42 The SDFs appear to change little with increasing lidocaine concentration, with [CH3CO2]− populating the area around the [BMIM]+ aromatic ring hydrogens, in agreement with the RDFs. Figure 5 shows the coordination number (CN) associated with the interaction of [BMIM]+ with [CH3CO2]− computed with respect to the COR of [BMIM]+ and the COM of [CH3CO2]−. We find that as the lidocaine concentration increases, the coordination number remains equal to 1 at
FIG. 3. Ball-and-stick representation of the united atom model of [BMIM]+ with the aromatic ring hydrogens labeled to facilitate discussion of the results.
a distance of approximately 0.47 nm. This distance is equal to the location of the first peak in the corresponding RDF. We therefore find that in the neat IL and with a lidocaine concentration of 0.1 and 0.3 mol fractions, we always have a [BMIM]+ associated with a [CH3CO2]− via a weak hydrogen bond.
FIG. 4. Radial distribution function (RDF) for the interaction between [BMIM]+ and [CH3CO2]− where the distance is computed with respect to the center-of-ring (COR) of [BMIM]+ and the center-of-mass (COM) of [CH3CO2]− (top pane), and between the aromatic hydrogens of [BMIM]+ with the [CH3CO2]− oxygens (bottom pane). A RDF is shown for the neat IL (x 1 = 0) and the IL plus 0.1 and 0.3 mol fractions of lidocaine.
084501-6
R. T. Ley and A. S. Paluch
FIG. 5. Coordination number (CN) for the interaction between [BMIM]+ and [CH3CO2]− where the distance is computed with respect to the center-of-ring (COR) of [BMIM]+ and the center-of-mass (COM) of [CH3CO2]−. A CN is shown for the neat IL (x 1 = 0) and the IL plus 0.1 and 0.3 mol fractions of lidocaine.
From the [BMIM]+ COR–[CH3CO2]− COM RDF in Fig. 4, the first local minimum occurs at approximately 0.71 nm, which corresponds to the radius of the first solvation shell. From Fig. 5, we have a coordination number of approximately 5.15, 4.75, and 4.08 at 0.71 nm for the neat IL (x 1 = 0) and with a lidocaine mole fraction of 0.1 (x 1 = 0.1) and 0.3 (x 1 = 0.3), respectively. As compared to the neat IL, this corresponds to a percent difference of 7.8 and 20.8% with a lidocaine mole fraction of 0.1 and 0.3, respectively. As the distance increases, the percent difference continues to increase toward 10% and 30%, the corresponding lidocaine mole fraction. Before looking at lidocaine, we consider [BMIM]+ –[BMIM]+ and [CH3CO2]−–[CH3CO2]− self-interactions in Fig. 6. Using the COR interaction distance, we do not observe ordering of [BMIM]+ with itself. We observe ordering for [CH3CO2]−. However, the first peak in the RDF reaches a height of just g (r) ≈ 1.5, suggesting a relatively weak interaction. 2. Lidocaine
Next, we consider lidocaine in [BMIM]+[CH3CO2]−. We first consider lidocaine–lidocaine self-interactions. The RDF computed with respect to the COR of lidocaine is shown in Fig. 7. We find that self-interactions of lidocaine are relatively weak. This is confirmed in Fig. S4 of the supplementary material42 where we observe a coordination number of 1 at approximately 0.89 and 0.76 nm for
J. Chem. Phys. 144, 084501 (2016)
FIG. 6. Radial distribution function (RDF) for the (self) interaction between [BMIM]+ and itself where the distance is computed with respect to the center-of-ring of [BMIM]+ (top pane) and for the (self) interaction between [CH3CO2]− and itself where the distance is computed with respect to the center-of-mass of [CH3CO2]−. An RDF is shown for the neat IL (x 1 = 0) and the IL plus 0.1 and 0.3 mol fractions of lidocaine.
a lidocaine mole fraction of 0.1 and 0.3, respectively. Lidocaine therefore is not associating with itself suggesting that the studied mole fractions are less than the equilibrium solubility of lidocaine. In Fig. 7, we additionally plot the RDF for lidocaine–[BMIM]+ interactions, with the distance computed with respect to the COR. We find that these interactions are weak and not ordered, suggesting that π–π interactions between the lidocaine and [BMIM]+ ring are not present. We found that the aromatic ring hydrogens of [BMIM]+ formed weak hydrogen bonds with [CH3CO2]− and that at all of the studied lidocaine concentrations we had an associated [BMIM]+[CH3CO2]− pair. Lidocaine is capable of both donating (N—H) and accepting (==O) hydrogen bonds, and we have found that it is extremely soluble in [BMIM]+[CH3CO2]−. Following this, we next consider hydrogen bonding interactions of lidocaine with [BMIM]+ and [CH3CO2]−. The CN and RDF for the lidocaine amine H (N—H) and the [CH3CO2]− oxygen, along with the lidocaine amide O (==O) and the [BMIM]+ aromatic ring hydrogens are shown in Fig. 8. Considering first the interaction between lidocaine and [CH3CO2]−, we observe a high intensity peak in the RDF occurring at r ≈ 0.19 nm with g(r) ≈ 11 and 13 for 0.1 and 0.3 mol fractions of lidocaine, respectively. Immediately after this peak is a sharp decrease to nearly 0. This short, strong, and well defined interaction is characteristic of a moderate hydrogen bond (due mostly to electrostatic interactions).83 For both lidocaine concentrations studied, the peak locations in the
084501-7
R. T. Ley and A. S. Paluch
FIG. 7. Radial distribution function (RDF) for lidocaine · · · · lidocaine and lidocaine · · · · [BMIM]+ intermolecular interactions where the distance is computed with respect to the center-of-ring of lidocaine and [BMIM]+.
RDF do not change. From the CN at both concentrations, we observe a coordination number of 1 at r ≈ 0.29 nm. Therefore, we find that at both lidocaine concentrations, lidocaine is associated with a single [CH3CO2]− via a moderate hydrogen bond. In the RDF, we observe a second peak that is less pronounced yet still has a well defined depletion well at each concentration. This second peak is the result of the second (equivalent) oxygen in the [CH3CO2]− associated with lidocaine and is reflected in the “S” shape of the coordination number between a value of 1 and 2. Next, we consider the possibility of lidocaine forming hydrogen bonds with [BMIM]+. Looking at the RDF for the interaction of the lidocaine amide oxygen (==O) with the aromatic ring hydrogens of [BMIM]+, we observe a peak at r ≈ 0.21 nm with a height of g(r) ≈ 2.8, followed by a well defined depletion well. In Fig. S5 of the supplementary material,42 we show RDFs for each aromatic ring hydrogen of [BMIM]+ and find that the interaction with aromatic ring hydrogen H1 (see Fig. 3) is strongest, with a first peak height of g(r) ≈ 3.5. While weaker than the interaction between the aromatic ring hydrogens and the [CH3CO2]− oxygen, we nonetheless have a well defined electrostatic interaction. From the CN, we observe a coordination number of 1 at approximately 0.35 and 0.38 nm for a lidocaine concentration of 0.1 and 0.3 mole fractions, respectively. The distance at which the coordination number is 2 and 3 is likewise offset and correspond to the two other (intramolecular) aromatic ring hydrogens.
J. Chem. Phys. 144, 084501 (2016)
FIG. 8. CNs (top pane) and RDFs (bottom pane) for the lidocaine amine hydrogen and the [CH3CO2]− oxygen interaction (N—H · · · · O) and for the interaction between the lidocaine amide O and the [BMIM]+ aromatic ring hydrogens (O · · · · C—H). Solid and dashed lines correspond to a lidocaine mole fraction of 0.1 and 0.3, respectively.
3. Solvation mechanism
We now have sufficient information to piece together the underlying molecular details or solvation mechanism. We found that [BMIM]+ associated with a [CH3CO2]− via a weak hydrogen bond. Likewise, lidocaine associated with a single [CH3CO2]− via a moderate hydrogen bond and with a single [BMIM]+ via an electrostatic interaction between the lidocaine amide oxygen and the aromatic ring hydrogens of [BMIM]+. Therefore, it must be the case that lidocaine associates with the [CH3CO2]− from one associated ion pair and with the [BMIM]+ of another (different) ion pair. In this light, let us re-examine the coordination number plot of Fig. 8. Lidocaine associates with a single [BMIM]+. The coordination number is one at the distance corresponding to the electrostatic interaction of the lidocaine amide oxygen and the aromatic ring hydrogens of [BMIM]+. The coordination number is equal to two and three at well defined distances corresponding to the second and third (intramolecular) hydrogen of the associated [BMIM]+ (i.e., on the same aromatic ring). This same [BMIM]+ should be associated with a single [CH3CO2]−. While lidocaine associates with a single [CH3CO2]−, we expect the next closest [CH3CO2]− oxygen to come from the [CH3CO2]− associated with [BMIM]+. From Fig. 8, the coordination number for the lidocaine amine hydrogen with the [CH3CO2]− oxygen is 3 (corresponding to nearest interaction with the second [CH3CO2]−) at a distance of approximately 0.64
084501-8
R. T. Ley and A. S. Paluch
J. Chem. Phys. 144, 084501 (2016)
C. Alchemical transformations
FIG. 9. Spatial distribution function (SDF) for the [CH3CO2]− oxygen (yellow) and the [BMIM]+ aromatic ring hydrogen H1 (red) around lidocaine and a snapshot from the simulations with 0.1 mol fractions of lidocaine.
and 0.66 nm for 0.1 and 0.3 mol fractions of lidocaine, respectively. This is approximately 0.15 and 0.10 nm further than the distance at which the lidocaine amide oxygen interacts with the second and third aromatic ring hydrogen of [BMIM]+ (i.e., where the coordination number is 2 and 3 for the lidocaine oxygen–[BMIM]+ hydrogen interaction). In Fig. 9 we provide the SDF for the [CH3CO2]− oxygen and aromatic ring hydrogen H1 (see Fig. 3) of [BMIM]+ around lidocaine. We observe that the [CH3CO2]− oxygen is highly populated around the lidocaine amine hydrogen, while simultaneously H1 of [BMIM]+ is populated adjacent to the [CH3CO2]− oxygen. Note that the [CH3CO2]− oxygen population density area is relatively large because it includes both [CH3CO2]− oxygens. Additionally, the H1 hydrogen of [BMIM]+ is highly populated around the lidocaine amide oxygen, while simultaneously the [CH3CO2]− oxygen is populated adjacent to the H1 hydrogen of [BMIM]+. In Fig. 9 we additionally provide a snapshot from the simulation trajectory to further illustrate the underlying mechanism. We clearly see that the [CH3CO2]− oxygen interacts with the lidocaine amine hydrogen, while the [BMIM]+ H1 hydrogen interacts with the lidocaine amide oxygen. At the same time, the oxygen of a different [CH3CO2]− interacts with a different aromatic ring hydrogen of the [BMIM]+ associated with lidocaine.
We have found that the favorable solvation of lidocaine in [BMIM]+[CH3CO2]− results from favorable hydrogen bonding between the lidocaine amine hydrogen and the [CH3CO2]− oxygen, favorable electrostatic interactions between the lidocaine amide oxygen with the [BMIM]+ aromatic ring hydrogens, while lidocaine does not interfere with the association of [BMIM]+ with [CH3CO2]−. We next seek to understand the effect of the solute size. We did this by computing the relative solubility of three solutes resulting from removing functional groups from the lidocaine scaffold (see Fig. 1) in [BMIM]+[CH3CO2]− relative to water. In all cases, the amide group is maintained so that we expect to maintain the favorable interactions between the solute with [BMIM]+ and [CH3CO2]−. The removal of functional groups decreases the solute size. Note that we have only decreased the solute size in the present study to ensure adequate solute conformational sampling. In each case, the solubility of the modified solute in the IL relative to water decreased. The solubility of the solutes in the IL relative to water from lowest to highest follows the trend: mod-lido C < mod-lido B < mod-lido A < lidocaine. We find that the relative solubility increases as the van der Waals volume increases. We therefore expect the performance of [BMIM]+[CH3CO2]− to continue to increase as the solute size increases. Of course, there must be an upper bound on the size of the solute at which point it begins to perturb the IL structure. A detailed study of solute size is reserved for a future study. Additionally, the free energy calculations of lidocaine in water were repeated with the TIP4P/200584 and TIP4P-Ew85 models to confirm that our choice of water model did not influence the findings of the present study. These and the results of all of the free energy calculations for all of the studied solutes are tabulated in the supplementary material of this manuscript.42
V. SUMMARY AND CONCLUSIONS
Using molecular simulation free energy calculations for a single lidocaine molecule in solution, the solubility of lidocaine in seven common room temperature ionic liquids was predicted to be 103 to 105 orders of magnitude greater than in water. While the infinite dilution approximation would break down at such large concentrations and the actual solubility would be smaller than predicted, the results none-the-less suggest that the solvation of lidocaine in the studied ionic liquids is highly favorable. Even for the case of the hydrophobic ionic liquids, the solubility is greatly enhanced relative to water. If the main reason for the favorable solubility was hydrogen bonding between lidocaine and the ionic liquid, we would expect the solubility to increase with increasing ionic liquid hydrogen bond basicity. We found this not to be the case for lidocaine, suggesting that hydrogen bonding alone does not explain the underlying mechanism. For the most favorable case, [BMIM]+[CH3CO2]−, long trajectories were computed with a lidocaine mole fraction of 0.1 and 0.3 to study solute-solvent and solvent-solvent interactions.
084501-9
R. T. Ley and A. S. Paluch
We found that the favorable solvation of lidocaine in [BMIM]+[CH3CO2]− results from favorable hydrogen bonding between the lidocaine amine hydrogen and the [CH3CO2]− oxygen, favorable electrostatic interactions between the lidocaine amide oxygen with the [BMIM]+ aromatic ring hydrogens, while lidocaine does not interfere with the association of [BMIM]+ with [CH3CO2]−. Essentially, lidocaine is able to interact favorable with two separate associated ion pairs. Therefore, while the ability of the anion to accept hydrogen bonds (which is related to the ionic liquid hydrogen bond basicity) is important, so too are the interactions between the ionic liquid cation and anion. Maintaining the amide group of lidocaine responsible for the favorable interactions with [BMIM]+ and [CH3CO2]−, three additional solutes were studied wherein functional groups were removed from the lidocaine scaffold. This had the effect of decreasing the solute van der Waals volume. We found that as the van der Waals volume decreased, the solubility in [BMIM]+[CH3CO2]− relative to water decreased. This suggests that as the solute size increases, IL performance should increase. At some point we would expect this trend to breakdown, but that is subject to another study. The underlying mechanism is very interesting and promising for future investigations. The favorable solvation of lidocaine in ionic liquids, including hydrophobic ionic liquids, has promising applications for the development of improved separation processes for pharmaceutical compounds. Additionally, ionic liquid use is promising for controlled release and delivery applications. ACKNOWLEDGMENTS
R.T.L. acknowledges summer financial support from the Miami University College of Engineering and Computing. A.S.P. gratefully acknowledges funding from the Miami University Committee on Faculty Research and start-up support from the Miami University College of Engineering and Computing. Computing support was provided by the Ohio Supercomputer Center and Miami University’s Research Computing Support group. 1C.
A. Angell, Y. Ansari, and Z. Zhao, Faraday Discuss. 154, 9 (2012). F. Brennecke and E. J. Maginn, AIChE J. 47, 2384 (2001). 3W. L. Hough, M. Smiglak, H. Rodríguez, R. P. Swatloski, S. K. Spear, D. T. Daly, J. Pernak, J. E. Grisel, R. D. Carliss, M. D. Soutullo et al., New J. Chem. 31, 1429 (2007). 4I. M. Marrucho, L. C. Branco, and L. P. N. Rebelo, Annu. Rev. Chem. Biomol. Eng. 5, 527 (2014). 5K. B. Smith, R. H. Bridson, and G. A. Leeke, J. Chem. Eng. Data 56, 2039 (2011). 6A. Forte, C. I. Melo, R. Bogel-Lukasik, and E. Bogel-Lukasik, Fluid Phase Equilib. 318, 89 (2012). 7A. D. dos Santos, A. R. C. Morais, C. I. Melo, R. Bogel-Lukasik, and E. Bogel-Lukasik, Fluid Phase Equilib. 356, 18 (2013). 8C. C. Weber, A. J. Kunov-Kruse, R. D. Rogers, and A. S. Myerson, Chem. Commun. 51, 4294 (2015). 9H. J. Park and M. R. Prausnitz, AIChE J. 61, 2732 (2015). 10Y. C. Kau, C. C. Liao, Y. C. Chen, and S. J. Liu, Materials 7, 6660 (2014). 11D. E. Ickowicz, L. Golovanevski, A. J. Domb, and C. F. Weiniger, Int. J. Pharm. 468, 152 (2014). 12A. S. Paluch, S. Parameswaran, S. Liu, A. Kolavennu, and D. L. Mobley, J. Chem. Phys. 142, 044508 (2015). 13A. S. Paluch, C. A. Vitter, J. K. Shah, and E. J. Maginn, J. Chem. Phys. 137, 184504 (2012). 2J.
J. Chem. Phys. 144, 084501 (2016) 14A.
R. Leach, Molecular Modelling: Principles and Applications, 2nd ed. (Pearson Education Limited, Harlow, England, 2001). 15D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. (Academic Press, San Diego, CA, 2002). 16X. Zhong, Z. Liu, and D. Cao, J. Phys. Chem. B 115, 10027 (2011). 17W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). 18B. Chen and J. I. Siepmann, J. Phys. Chem. B 110, 3555 (2006). 19J. L. Rafferty, L. Sun, J. I. Siepmann, and M. R. Schure, Fluid Phase Equilib. 290, 25 (2010). 20N. Rai and J. I. Siepmann, J. Phys. Chem. B 111, 10790 (2007). 21N. Rai, D. Bhatt, J. I. Siepmann, and L. E. Fried, J. Chem. Phys. 129, 194510 (2008). 22N. Rai and J. I. Siepmann, J. Phys. Chem. B 117, 273 (2013). 23C. D. Wick, J. M. Stubbs, N. Rai, and J. I. Siepmann, J. Phys. Chem. B 109, 18974 (2005). 24Note that TraPPE-EH does parameterize aniline in Ref. 21. We chose not to use the N LJ parameters from aniline because it is a primary amine, and it was previously shown in Ref. 23 that there is an appreciable change in LJ parameters in going from a primary to secondary amine. Also, in Ref. 23, the primary amide LJ parameters are the same as the primary amine (i.e., only the charges change). Therefore, we took the secondary amide LJ parameters for N to be the same as for a secondary amine. The H has LJ parameters of 0 in primary and secondary amines and primary amides. The LJ parameters for H in a secondary amide were therefore taken to also be 0 in the present study. 25J. M. Stubbs, J. J. Potoff, and J. I. Siepmann, J. Phys. Chem. B 108, 17596 (2004). 26M. G. Martin and J. I. Siepmann, J. Phys. Chem. B 102, 2569 (1998). 27The alkyl groups (CH and CH ) were modeled as a single united-atom 3 2 pseudoatom as a result of the parameterization of the TraPPE-EH force field for n-alkanes which places the LJ site for a hydrogen atom at the center of the corresponding bond,86 and the complication of implementing such a model in a molecular dynamics framework. 28Y. Zhao and D. G. Truhlar, Theor. Chem. Acc. 120, 215 (2008). 29C. J. Cramer, Essentials of Computational Chemistry (John Wiley & Sons Ltd., Chichester, West Sussex, England, 2002). 30A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. Phys. Chem. B 113, 6378 (2009). 31A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. Phys. Chem. B 113, 4538 (2009). 32R. F. Ribeiro, A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. Comput.Aided Mol. Des. 24, 317 (2010). 33M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson et al., 09, Revision A.02, Gaussian, Inc., Wallingford, CT, 2009. 34C. I. Bayly, P. Cieplak, W. D. Cornell, and P. A. Kollman, J. Phys. Chem. 97, 10269 (1993). 35P. Cieplak, W. D. Cornell, C. Bayly, and P. A. Kollman, J. Comput. Chem. 16, 1357 (1995). 36D. A. Case, T. Cheatham, T. Darden, H. Gohlke, R. Luo, K. M. Merz, A. Onufriev, C. Simmerling, B. Wang, and R. Woods, J. Comput. Chem. 26, 1668 (2005). 37D. A. Case, T. A. Darden, T. E. Cheatham III, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, R. C. Walker, W. Zhang, K. M. Merz et al., Amber 12 (University of California San Francisco, San Francisco, CA, 2012), p. 12. 38J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, J. Comput. Chem. 25, 1157 (2004). 39J. Wang, W. Wang, P. A. Kollman, and D. A. Case, J. Mol. Graphics Modell. 25, 247 (2006). 40A. W. Sousa da Silva and W. F. Vranken, BMC Res. Notes 5, 367 (2012). 41A. W. Sousa da Silva and W. F. Vranken, acpype: AnteChamber PYthon Parser interfacE, accessed May 1, 2014, http://www.gromacs.org/. 42See supplementary material at http://dx.doi.org/10.1063/1.4942025 for GROMACS force field files (type “itp” and “xyz”) for the studied ionic liquids and solutes are provided in gromacs_itp.txt and gromacs_xyz.txt; additional SDFs, RDFs, and coordination numbers are provided to facilitate discussion along with tabulated results of the free energy calculations are provided in supp_info.pdf. 43L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, J. Comput. Chem. 30, 2157 (2009). 44See http://www.ime.unicamp.br/∼martinez/packmol/ for Packmol: Packing Optimization for Molecular Dynamics Simulations, accessed August 1, 2013.
084501-10 45D.
R. T. Ley and A. S. Paluch
van der Spoel, E. Lindahl, B. Hess, and the GROMACS development team, GROMACS User Manual, version 4.6.3, 2013, ftp://ftp.gromacs.org/ pub/manual/manual-4.6.3.pdf. 46H. J. C. Berendsen, Simulating the Physical World: Hierarchial Modeling from Quantum Mechanics to Fluid Dynamics (Cambridge University Press, New York, NY, 2007). 47G. Bussi, D. Donadio, and M. Parrinello, J. Chem. Phys. 126, 014101 (2007). 48G. Bussi and M. Parrinello, Comput. Phys. Commun. 179, 26 (2008). 49G. Bussi, T. Zykova-Timan, and M. Parrinello, J. Chem. Phys. 130, 074101 (2009). 50H. J. C. Berendsen, J. P. M. Postma, A. DiNola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984). 51G. J. Martyna, D. J. Tobias, and M. L. Klein, J. Chem. Phys. 101, 4177 (1994). 52G. J. Martyna, M. L. Klein, and M. E. Tuckerman, J. Chem. Phys. 97, 2635 (1992). 53B. Hess, C. Kutzner, D. van der Spoel, and E. Lindal, J. Chem. Theory Comput. 4, 435 (2008). 54S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel et al., Bioinformatics 29, 845 (2013). 55See http://www.gromacs.org/ for GROMACS: Fast, flexible, free, accessed August 1, 2013. 56B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, J. Comput. Chem. 18, 1463 (1997). 57B. Hess, J. Chem. Theory Comput. 4, 116 (2008). 58J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comput. Phys. 23, 327 (1977). 59H. C. Andersen, J. Comput. Phys. 52, 24 (1983). 60M. Deserno and C. Holm, J. Chem. Phys. 109, 7678 (1998). 61G. B. Fuerst, R. T. Ley, and A. S. Paluch, Ind. Eng. Chem. Res. 54, 9027 (2015). 62K. S. Shing and S. T. Chung, J. Phys. Chem. 91, 1674 (1987). 63D. A. Kofke and P. T. Cummings, Mol. Phys. 92, 973 (1997). 64M. R. Shirts, J. W. Pitera, W. C. Swope, and V. S. Pande, J. Chem. Phys. 119, 5740 (2003). 65D. A. Kofke and P. T. Cummings, Fluid Phase Equilib. 150-151, 41 (1998). 66Free Energy Calculations: Theory and Applications in Chemistry and Biology, Springer Series in Chemical Physics Vol. 86, edited by C. Chipot and A. Pohorille (Springer, New York, NY, 2007). 67C. H. Bennett, J. Comput. Phys. 22, 245 (1976). 68M. R. Shirts, E. Bair, G. Hooker, and V. S. Pande, Phys. Rev. Lett. 91, 140601 (2003). 69N. Lu, J. K. Singh, and D. A. Kofke, J. Chem. Phys. 118, 2977 (2003). 70M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129, 124105 (2008). 71T. C. Beutler, A. E. Mark, R. C. van Schaik, P. R. Gerber, and W. F. van Gunsteren, Chem. Phys. Lett. 222, 529 (1994).
J. Chem. Phys. 144, 084501 (2016) 72M.
R. Shirts and V. S. Pande, J. Chem. Phys. 122, 134508 (2005). Steinbrecher, D. L. Mobley, and D. A. Case, J. Chem. Phys. 127, 214108 (2007). 74See https://github.com/choderalab/pymbar for PyMBAR: Python implementation of the Multistate Bennett Acceptance Ratio (MBAR), accessed May 1, 2014. 75J. D. Chodera, W. C. Swope, J. W. Pitera, C. Seok, and K. A. Dill, J. Chem. Theory Comput. 3, 26 (2007). 76P. V. Klimovich, M. R. Shirts, and D. L. Mobley, J. Comput.-Aided Mol. Des. 29, 397 (2015). 77In summary, from m = 0 to 14, we have λ LJ = {0.0, 0.1, 0.2, 0.3, 0.4, m 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0} and λ elec m = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.50, 0.71, 0.87, 1.00}. A detailed discussion of best practices when coupling/decoupling intermolecular interactions during free energy calculations is provided in the recent review by Shirts and Mobley.87 78M. Brehm and B. Kirchner, J. Chem. Inf. Model. 51, 2007 (2011). 79B. Kirchner, TRAVIS–Trajectory Analyzer and Visualizer, accessed November 13, 2015, http://www.travis-analyzer.de/. 80The conformation of the central molecule was generated by TRAVIS from the initial system configuration. The SDF was visualized using VMD88,89 with the central molecule represented as a CPK model and the second species or interaction site was modeled as an “isosurface.” 81Note, the experimental mole fraction solubility of lidocaine in water at 298.15 K is 3 × 10−4.90 Lidocaine in the ILs at the large mole fraction solubilities predicted here would be beyond the infinite dilution limit, and we would expect the actual solubility to be less than predicted. Regardless, it stands that there are extremely favorable interactions present in the ILs leading to extremely favorable solubilities over water. 82A. F. M. Cláudio, L. Swift, J. P. Hallett, T. Welton, J. A. P. Coutinho, and M. G. Freire, Phys. Chem. Chem. Phys. 16, 6593 (2014). 83G. A. Jeffrey, An Introduction to Hydrogen Bonding (Oxford University Press, New York, NY, 1997). 84J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005). 85H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys. 120, 9665 (2004). 86B. Chen and J. I. Siepmann, J. Phys. Chem. B 103, 5370 (1999). 87M. R. Shirts and D. L. Mobley, “An introduction to best practices in free energy calculations,” in Biomolecular Simulations: Methods and Protocols (Humana Press, Totowa, NJ, 2013), pp. 271–311, ISBN: 978-1-62703017-5, http://dx.doi.org/10.1007/978-1-62703-017-5_11. 88W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics 14, 33 (1996). 89See http://www.ks.uiuc.edu/Research/vmd/ for VMD: Visual Molecular Dynamics, accessed November 25, 2015. 90Solubility and Related Properties of Large Complex Chemicals Part 1: Organic Solutes Ranging from C4 to C40, edited by J. Marrero and J. Abildskov (DECHEMA, Frankfurt am Main, Germany, 2003). 73T.
