Introduction

Ionic liquids (IL) are a unique class of high-performance chemical compounds with many applications in electrochemistry (Armand et al. 2009), synthesis (Itoh 2017), catalysis (Pârvulescu and Hardacre 2007; van Rantwijk and Sheldon 2007) as well as solvation (Hallett and Welton 2011) and extraction processes (Ventura et al. 2017). Their tunable properties via variation or modification of either the cation or the anion as well as their shared properties such as low vapor pressure, low flammability, and high thermal and electrochemical stability make them interesting solvents for various applications and started a hype in the beginning of the twenty-first century.

Kunz and Häckl (2016) as well as Wasserscheid pointed out that the expectations followed a Gartner cycle as visible in Fig. 1. The initial success of ionic liquids in (bio-) catalysis (van Rantwijk and Sheldon 2007) and solvation (Welton 1999; Gutowski et al. 2003), e.g., dissolution of cellulose (Swatloski et al. 2002), served as a technology trigger. The vast number of \(10^{18}\) possible combinations of cations and anions (Holbrey and Seddon 1999) seemed to promise the design of an optimal solvent for a particular application. As a result, the term “designer solvent” was coined (Freemantle 1998) and soon joined by the rash promise “ILs are green solvents” (Earle and Seddon 2009) leading to an exponential number of paper submissions to the journal Green Chemistry. Also, the availability of ILs increased due to beginning industrial production (Plechkova and Seddon 2008).

Fig. 1
figure 1

Gartner hype cycle on ionic liquids adapted from Kunz and Häckl (2016), Wasserscheid (2017)

Full size image

However, soon, first doubts were cast on the environmental benefit of ILs (Thuy Pham et al. 2010) and problems on the hydrolysis of some ILs were reported (Steudte et al. 2012). Additionally, many ionic liquids did not have a unique selling point or their improved performance did not justify their increased cost (Plechkova and Seddon 2008; Kunz and Häckl 2016). As a result, the hype cooled down in the past decade reaching the so-called Trough of Disillusionment of a Gartner cycle.

Nevertheless, scientific research did not abate, visible by the roughly exponential increase in publications indicated by the blue dashed line in Fig. 1. Here, we normalized the publication growth obtained from a search in the “Web of knowledge” to the current number of publications up to 2017. This way, we can compare the development with publication efforts concerning simulations of ionic liquids shown as an orange dashed line in Fig. 1 which agrees almost quantitatively with the overall publication output on ILs.

Quite generally, hype cycles can be decomposed into two processes (Sasaki 2015): a hype stage and an implementation stage. The latter can be fitted by an S-shaped Gompertz function

$$ f(t) = a \cdot e^{-e^{-k (t-t_{0})}}. $$
(1)

Taking the publication efforts depicted in Fig. 1 as a measure of implementation, the “Slope of Enlightenment” phase can be extrapolated. In addition to the amplitude a and the stretching factor k, the inflection point \(t_{0}\) marks the year in which the implementations start to level off. In the case of the ILs, the year 2020 is obtained for both the total publications and the publications concerning the simulation of ionic liquids. The “Plateau of Productivity” will be reached in twenty to thirty years based on the current data.

In order to get there, many scientific questions have yet to be answered. In particular, the knowledge gained by simulations on biomolecules in ionic liquids and their mixtures is lagging behind the general ionic liquid trend as visible by the green dotted line in Fig. 1. The current review tries to summarize the current efforts of biomolecular solvation in ionic liquids and their mixtures from a molecular dynamics (MD) perspective.

Solvation starts at the first solvation layer which is the transition region between the bulk solvent and the solute as depicted in Fig. 2. Here, coordination numbers of particular solvent molecules and hydrogen bonding are interesting features to be computed from MD trajectories. Moreover, radial distribution functions (solid black line in Fig. 2) detecting the accumulation or depletion of solvent species at the biomolecular surface and their interpretation in terms of Kirkwood-Buff theory (Lesch et al. 2015; Diddens et al. 2017; Smiatek 2017) are reported in “The solvation layer”.

Fig. 2
figure 2

Computational analysis of the solvation interaction between biomolecules, their solvation layers, and the bulk solvent

Full size image

Due to the surrounding hydrophilic or hydrophobic solvent, the structure and stability of the biomolecular solute may adapt to the environment. Depending on the size and shape of a solute, different mechanisms contribute to solvent-induced changes. Hence, we distinguish between large flexible biomolecules and small rigid drugs, employing different analysis routines as sketched in Fig. 2:

  1. 1.

    Large proteins, polypeptides, or polymeric biomolecules possess a secondary structure which can deform during the solvation process (indicated by the gray arrows in Fig. 2). In a hydrophilic environment, for example, hydrophilic subunits are supposed to be at the surface of the biomolecule and in contact with the solvent. Hydrophobic subunits are covered in the core of the biomolecule. In addition to refolding, the size of a biomolecule may change when moving from a hydrophilic to a hydrophobic environment or vice versa. Computational size and shape observables of this biomolecular class are discussed in “Large flexible proteins”.

  2. 2.

    In contrast, small rigid drugs cannot increase their size or change their shape significantly. Here, other mechanisms like solute aggregation can be monitored by molecular dynamics. The inclusion of the solute in the solvent network can also be computed by free energy calculations. Furthermore, functional groups and their respective influence on the solvation are of significant importance. This will be discussed in “Small rigid drugs.”

The solvation layer

The tunable solvation properties (Anderson et al. 2002; Gutowski et al. 2003) of ILs gave reasons for the hype on ILs at the start of this century. They were considered as green replacements for organic solvents in general (Cull et al. 2000) and for biological processes in particular (Yang and Pan 2005; Park and Kazlauskas 2003). For example, Rogers and co-workers reported on the dissolution of cellulose in 2002 (Swatloski et al. 2002). In 2005, BASF licensed the exclusive use of various intellectual property rights from this group. In 2007, the Rogers group continued their research on dissolution of lignocellulosic materials (Fort et al. 2007). Mikkola et al. (2007), Anugwom et al. (2012), Raut et al. (2015), Kilpeläinen et al. (2007), King et al. (2011), and Parviainen et al. (2013) worked on cellulose processing. Excellent reviews summarize the current development (Pinkert et al. 2009; Brandt et al. 2013; Mäki-Arvela et al. 2010; Zhang et al. 2017) and show that we are still in the “Slope of Enlightenment” phase during the implementation stage.

Due the high cost of ILs compared to the conventional alternatives, not many large-scale applications were introduced. In order to reduce the costs, ILs were mixed with other cheap solvents, in particular water. Since only few proteins are soluble in pure ionic liquids (Shao 2013; Kragl et al. 2002), homogeneous and heterogeneous aqueous ionic liquid mixtures were used (Kragl et al. 2002). The first ionic liquid-based aqueous biphasic system consisted of 1-butyl-3-imidazolium chloride mixture in aqueous K3PO4 and was reported by Rogers and co-workers in 2003 (Gutowski et al. 2003) and may be considered as a “Technology Trigger” (see Fig. 1). Recently, biphasic systems regained interest, in particular by the group of Coutinho for extraction (Ventura et al. 2009, 2017; Pei et al. 2009; Dreyer and Kragl2008; Pereira et al. 2010; Lee et al. 2017; Cláudio et al. 2010) and for purification of biomolecules (Pereira et al. 2013). An excellent review was given by Freire et al. (2012).

Although of great importance for the rational choice of an ionic liquid for any given application (Weingärtner et al. 2012), predicting the effect of an ionic liquid on the stability of a specific protein still remains difficult (Tung and Pfaendtner 2016) due to the complex and little understood molecular interactions between ionic liquid and protein (Burney et al. 2015) demonstrating that we are still in the “Slope of Enlightenment” phase and have not reached the “Plateau of Productivity” in Fig. 1. A popular classification of the effect of ionic liquids on proteins already mentioned is the extension of the well-known Hofmeister series (Zhang and Cremer 2006) to aqueous ionic liquid solutions (Constantinescu et al. 2007). The Hofmeister series ranks ions by their specific ion effect on protein properties such as stability or solubility (Zhang and Cremer 2006; Schröder 2017). Ions are classified either as kosmotropes (structure-makers) or chaotropes (structure-breakers) by their ability to influence water structure (Tung and Pfaendtner 2016; Constantinescu et al. 2007; Zhang and Cremer 2006). Generally, kosmotropic anions stabilize proteins, while chaotropic anions lead to destabilization (Tung and Pfaendtner 2016). Still, exceptions to the Hofmeister series have been observed (Zhang and Cremer 2009), not only due to the fact that specific protein properties such as charge and surface characteristics (Schröder 2017; Constantinescu et al. 2007; Tung and Pfaendtner 2016) are not considered. Moreover, the ability of ions to influence bulk water structure in a kosmo- or chaotropic manner has been questioned (Omta et al. 2003).

Less publications exist on the computational research of aqueous ionic liquid mixtures (Bhargava and Klein 2009; Varela et al. 2015; Chang et al. 2010) without additional solutes as these studies focus on more fundamental properties of IL mixtures rather than a particular application. The effect of water on the viscosity of an IL/water mixture was studied by Kelkar and Maginn (2007). Furthermore, the alkyl chain length of the IL cations changed the collective structure of the aqueous IL mixture (Bhargava and Klein 2009). We started our computational analysis on the collective network between water and various ionic liquids in 2007 (Schröder et al. 2007, 2008, 2009, 2014) with a focus on dielectric properties. Smiatek and co-workers published interesting papers on Kirkwood-Buff analysis of protein solvation in IL/water mixtures (Diddens et al. 2017; Smiatek 2017). Gupta and co-workers reported MD simulations on the interaction of cellulose with ionic liquids (Gupta et al. 2013). As a complete review of computational analysis of interactions of biomolecules with IL mixtures is out of the scope of this review, we will focus on the basic techniques to encourage other authors to contribute to the understanding of biomolecular solvation in ionic liquid mixtures by MD simulations . Hereafter, we will give a compact overview of solvation layer properties obtainable by MD simulations as well as their interpretation.

Spatial distribution of the solvent

The spatial distribution of molecular species j around a reference site i is often characterized by the radial distribution function \(g_{ij}(r)\),

$$ g_{ij}(r) = \frac{ \sum\limits_{j} \delta(r-r_{ij})}{4 \pi r^{2} dr \cdot} \cdot \frac{1}{\rho_{j}} $$
(2)

which is the ratio between the local density of j in a spherical shell \(4 \pi r^{2} dr\) and the global density \(\rho _{j}\) of species j. The distance \(r_{ij}\) is usually defined between the center of masses of the molecules i and j or between the coordinates of particular atoms belonging to these molecules.

Coordination numbers

To calculate the coordination number \(n_{ij}\) of species around the biomolecule i, spherical integration of the radial distribution function is commonly used:

$$ n_{ij} = \int\limits_{0}^{R_{1}} 4 \pi r^{2} g_{ij}(r) dr $$
(3)

Micaêlo and Soares found that the coordination number of water around cutinase in aqueous 1-butyl-3-methylimidazolium nitrate [C4mim][NO3] and 1-butyl-3-methylimidazolium hexafluorophosphate [C4mim][PF6] mixtures is reduced at the protein surface (Micaêlo and Soares 2008). Similar removal of water was observed for hydrated Candida rugosa lipase 1 in [C4mim][PF6] and [C4mim][NO3] (Burney and Pfaendtner 2013).

A problem that may arise when comparing coordination numbers of different IL ions is the fact that the coordination number depends on ion density. For the analysis of cation and anion coordination numbers around CAL-B, Klähn et al. (2011a) normalized the coordination number by the average number of the respective species within 10 Å of the protein to facilitate comparison. They observed the highest ion densities around charged residues. Interestingly, although density of cations and anions was highest around oppositely charged residues, densities around residues with the same charge was also higher than for polar and non-charged residues (Schröder 2017; Klähn et al. 2011a), which is an effect of the strong interionic interactions in ionic liquids (Klähn et al. 2011b). Furthermore, a diffusion of cations into the active site of the protein occurred, possibly hindering substrate access. A similar occupation of the active site by cations was observed in simulations of xylanase II from trichoderma longibrachiatum in aqueous [C2mim][EtSO4] and [C2mim][OAc] (Jaeger and Pfaendtner 2013).

Steinhauser and co-workers noted that Eq. 3 not only requires a meaningful definition of a solvation shell radius \(R_{1}\) but also assumes sphericity of the solutes (Haberler and Steinhauser 2011; Haberler et al. 2011; Zeindlhofer et al. 2018). Although this approach is feasible for many small solutes, it may give highly inaccurate results in case of anisotropic solutes (Zeindlhofer et al. 2017, 2018) like large proteins (Neumayr et al. 2010) due to non-symmetric excluded volume effects. Furthermore, for heterogeneous solvents with different molecular sizes, more than one shell radius would be necessary to describe a solvation shell appropriately (Haberler et al. 2011).

Voronoi-based coordination numbers

An alternative method of spatial decomposition is Voronoi decomposition (Schröder et al. 2009). It is parameter-free and also allows for straightforward decomposition of multiple solvation shells. In Voronoi tessellation, each point in space representing an atomic coordinate is assigned an irregular polyhedron that contains all space closer to its associated reference point than to any other point in the system. If two points share a polyhedron face, they are defined as neighbors. Solvent molecules that share a face with the biomolecule are considered first-shell members, molecules that share a face with first-shell members are considered second shell, and so on. This analysis can be performed for each time step during a trajectory to yield mean residence times of particular species at the surface of the biomolecule (Haberler and Steinhauser 2011; Haberler et al. 2011). This shell assignment can also be applied to the radial distribution function (Neumayr et al. 2010; Zeindlhofer et al. 2018) as depicted in Fig. 3. The decomposition of the displayed \(g_{ij}(r)\) into spherical solvation shells discussed in the last section is not possible as no clear first minimum of the radial distribution can be detected. Furthermore, the different shells overlap to a large extent. Consequently, the spherical integration performed in Eq. 3 will not yield a correct coordination number. However, integrating the shell contributions in Fig. 3 results in the exact coordination number obtained from counting neighbors by shared faces of the solute and solvent polyhedra.

Fig. 3
figure 3

Radial distribution function gij(r) of water oxygens (j) around the center of mass of calbindin i and the corresponding decomposition into Voronoi shells. Reproduced from Neumayr et al. (2010) with permission of AIP Publishing

Full size image

Employing Voronoi decomposition, Haberler and Steinhauser confirmed a surplus of C\(_{2}\textit {mim}^{+}\) in the first solvation shell of various proteins despite the positive net charge of the protein (Haberler and Steinhauser 2011; Haberler et al. 2011). This is very likely due to expulsion of the cations from the strong anion-water network, as also observed in other studies (Schröder et al. 2007, 2009; Haberler and Steinhauser2011). Interestingly, the mean residence time of the anions at the protein surface is longer compared to the cations (Haberler and Steinhauser 2011) which can be explained by stronger interactions with the oppositely charged amino acid (Schröder 2017; Jaeger and Pfaendtner 2013; Klähn et al. 2011a; Micaêlo and Soares 2008).

Kirkwood-Buff theory

Another approach utilizing radial distribution functions is the so-called Kirkwood-Buff theory, which relates structural properties of the solutions

$$ G_{ij}(R) = {\int\limits_{0}^{R}} \left( g_{ij}(r)-1 \right) 4 \pi r^{2} dr, $$
(4)

to thermodynamic quantities (Kirkwood and Buff 1951). The Kirkwood-Buff integral is often interpreted as excess volume around a solute species when compared to an ideal solution. It is also possible to compute Kirkwood-Buff integrals of distinct solvation shells based on Voronoi tessellation (Zeindlhofer et al. 2018).

As a standard notation, \(i,j = 1 {\ldots } 3\) denote water, the biomolecular solute, and the co-solvent (in our case the ionic liquid), respectively. Distance-dependent binding coefficients ν23(R) between the solute \(i = 2\) and the ionic liquid \(j = 3\) can be obtained as follows:

$$\begin{array}{@{}rcl@{}} \nu_{23}(R) &=& \frac{1}{n\pm}\left( \rho_{3}\left( G_{23}(R) - G_{21}(R) \right)-\frac{q^{\text{solute}}}{e}\right) \end{array} $$
(5)
$$\begin{array}{@{}rcl@{}} &=&\frac{1}{n\pm} \left( n_{23}(R)-\frac{\rho_{3}}{\rho_{1}}n_{21}(R)-\frac{q^{\text{solute}}}{e}\right) \end{array} $$
(6)

using coordination numbers defined in Eq.3 as well as the Donnan equilibrium condition (Pierce et al. 2008; Smith 2004; Smiatek 2017). Here \({q^{\text {solute}}}\) corresponds to the charge of the biomolecule in units of the elementary charge \(e = 1.6022 \times 10^{-19}\) A s and \(n{\pm }\) equals the number of indistinguishable ion species. The binding coefficient determines the transfer free energy (Smiatek 2017).

$$ F_{23}(R) = - k_{B} T \cdot \nu_{23}(R) $$
(7)

which estimates the excess energy at a given temperature T that is needed to transfer the ionic liquid from infinite distance to close proximity around the biomolecule.

In a MD study of a small \(\beta \)-hairpin peptide in aqueous [C2mim][OAc], where both native and denaturated conformations of the peptide were simulated, Kirkwood-Buff preferential binding coefficients revealed that the unfolded conformation attracts more anions while no difference in cation affinity is observed (Lesch et al. 2015). This agrees with the fact that the anionic influence on the denaturation process is stronger (Klähn et al. 2011a). The unfolding effect of [C2mim][OAc] on the peptide was also confirmed by experiment (Patel et al. 2014). In a follow-up publication by the same group, a similar approach was used to study the behavior of the same \(\beta \)-hairpin peptide in aqueous [C2mim][BF4] and [C2mim]Cl. The binding strength order of the anions was identified as OAc−≫ BF\(_{4}^{-} >\) Cl−, and again it was found that OAc− strongly binds to the denatured state, which is not observable for BF\(_{4}^{-}\) and Cl−. The authors suggested that the larger acetate anion exhibits a conformation-dependent binding behavior due to specific interactions with the protein (Diddens et al. 2017). Furthermore, they found only negligible dependence of the binding behavior of C\(_{2}\textit {mim}^{+}\) on the anion.

The interactions of caffeine, gallic acid, protocatechuic acid, and quercetin with aqueous mixtures of the ionic liquid [C2mim][OAc] were studied by a shell-resolved approach using Voronoi tessellation. As shown in the respective publications (Zeindlhofer et al. 2017, 2018), the investigated biomolecules are flat, aromatic systems of high anisotropy, which makes Voronoi tessellation necessary to distinguish between several solvation shells. Since Voronoi tessellation allows straightforward calculation of solvent shell volumes, a shell-wise calculation of Kirkwood-Buff interaction parameters analogous to Kirkwood-Buff integrals (Zeindlhofer et al. 2018) can be computed. However, the comparison of Kirkwood-Buff interaction parameters between solutes of different size is complicated, since they depend on the solute size. In contrast, the shell-wise potential of mean force (Zeindlhofer et al. 2017) measures solvent affinities of solutes differing in size via the ratio of shell to bulk concentration of a species.

Solvation dynamics

The role of ionic liquids in solvation dynamics (Mandal and Samanta 2005; Samanta 2010; Roy and Maroncelli 2012; Terranova and Corcelli 2013; Liang et al. 2014; Heid and Schröder 2018) as well as the solvation dynamics of biomolecules (Sajadi et al. 2014; Heid and Schröder 2016; 2017) has attracted significant interest over the past decade. The time-dependent Stokes shift

$$ S(t) = \frac{\nu(t)-\nu(\infty)}{\nu(0)-\nu(\infty)} $$
(8)

monitors the transient behavior of the solvent relaxation after an electronic excitation of a dissolved chromophore. THz absorption bands of biomolecular hydration layers are generally swamped by absorption from bulk solvent. However, linking the chromophore to the biomolecule, the solvation properties of the immediate solvent molecules can be experimentally studied (Sajadi et al. 2014) as shown for the chromophore oxyquinolinium betaine attached to trehalose in Fig. 4.

Fig. 4
figure 4

The chromophore oxyquinolinium betaine can be directly linked to the biomolecule of interest. Upon laser excitation, the local electric field (red lines) is changed and the response of the solvent molecules in the direct vicinity of the trehalose can be studied (Sajadi et al. 2014). Copyright 2014, American Chemical Society

Full size image

In MD simulations, this experiment can be modeled by equilibrium (Roy and Maroncelli 2012; Schmollngruber et al. 2013; Terranova and Corcelli 2013) and non-equilibrium simulations (Heid and Schröder 2016, 2017, 2018). However, the latter approach seems to be more appropriate (Heid and Schröder 2018). The change of the Coulomb interaction energy \({\Delta } U^{elec}(t)\) of the chromophore atoms j with the solvent atoms i

$$ {\Delta} U^{elec}(t) = \frac{1}{4 \pi \epsilon_{0}} \sum\limits_{j} \sum\limits_{i} \frac{{\Delta} q_{j} \cdot q_{i}}{r_{ij}(t)} $$
(9)

relaxes during the non-equilibrium simulation after switching the partial charges \(q_{j}\) of the solute from ground to excited state values:

$$ S(t)= \frac{{\Delta} U^{elec}(t)-{\Delta} U^{elec}(\infty)}{{\Delta} U^{elec}(0)-{\Delta} U^{elec}(\infty)} $$
(10)

In contrast to neutral liquids like water, ionic liquids can also respond to the change in the local electric field via translation (Terranova and Corcelli 2013; Schmollngruber et al. 2013), in particular the anions. The response of induced solvent dipoles is restricted to the first solvation shell (Schmollngruber et al. 2013). Using polarizable non-equilibrium molecular dynamics simulations, not only the normalized Stokes shift relaxation function but also the pure shift \({\Delta } U^{elec}(t)-{\Delta } U^{elec}(\infty )\) can be reproduced (Heid and Schröder 2018).

Large flexible proteins

One of the “Technology Triggers” in Fig. 1 for protein stability was published in 2000 by Summers et al. who showed that ethylammonium nitrate can successfully refold denatured egg-white lysozyme (Summers and Flowers 2000). Several groups observed that ionic liquids can significantly stabilize certain proteins and enzymes (Lozano et al. 2001; Fujita et al. 2007; Weingärtner et al. 2012), especially lipases (Kaar et al. 2003; Ulbert et al. 2005; Lai et al. 2011), even at high temperatures (Ulbert et al. 2005). In some cases, enzyme activity, enantioselectivity, and regioselectivity in ionic liquids are increased compared to conventional organic solvents (Schöfer et al. 2001) as discussed in more detail in several reviews (Kragl et al. 2002; Park and Kazlauskas 2003; Yang and Pan 2005; van Rantwijk and Sheldon 2007; Moniruzzaman et al. 2010; Lai et al. 2011; Naushad et al. 2012; Gao et al. 2015; Chen et al. 2010; Zhao 2010; Weingärtner et al. 2012; Patel et al. 2014; Benedetto and Ballone 2016; Smiatek 2017; Schröder 2017; Zhao 2016).

An observation that many groups agree with is that the presence of water in the ionic liquid is crucial for protein stability (Klähn et al. 2011a; Lozano et al. 2001). It was also observed that the effect of the anion on protein structure is more pronounced than the effect of the cation in most cases (Weingärtner et al. 2012; Constantinescu et al. 2007; Kaar et al. 2003). Furthermore, the destabilizing effect of ionic liquids in the aqueous mixture is a function of the ion concentration (Senske et al. 2016) as shown in Fig. 5. Stabilizing effects (ΔTm > 0) were only observed for dihydrogenphosphate salts as well as sodium and potassium chloride above 2 M. Here, “enlightenment” from MD simulations on the mechanism is highly welcome.

Fig. 5
figure 5

Salt-induced shifts of the melting temperature ΔTm of RNase in aqueous ionic liquid mixtures with various ion concentrations (Senske et al. 2016)—Published by the PCCP Owner Societies. (a,b) chlorides (c) bromides (d) sodium salts (e) C2mim based salts (f) stabilizing oxo-anions

Full size image

In the quest of rationalizing ionic liquid protein interactions, MD simulations provide a useful tool for investigations at a molecular level and can complement experimental findings (Haberler and Steinhauser 2011). Although the first MD simulations of proteins were reported in the late 1970s (McCammon et al. 1977), the first MD simulation of a protein in an aqueous ionic liquid was published only in 2008 (Micaêlo and Soares 2008) revealing the time delay of computational analysis in Fig. 1 concerning interaction of proteins with ionic liquids. In Micaêlo and Soares (2008), investigated the behavior of the serine protease cutinase in [C4mim][PF6] and [C4mim][NO3] at various water contents. For [C4mim][PF6], 5–10% of water was necessary for optimal enzyme stabilization (Micaêlo and Soares 2008) indicating the importance of even small amounts of water.

A subsequent publication by the same group further investigated the influence of the protein surface characteristics on stability in ILs by conducting simulations of Candida rugosa lipase and Bos taurus\(\alpha \)-chrymotopsine with modified surface charges (Burney et al. 2015). By replacing surface lysines with glutamate, artificial mutations with positive-to-negative surface charge ratios were examined in simulation. On a 50-ns timescale, no stability difference for the different enzymes was found, but the authors explicitly state that such a timescale is likely too short (Burney et al. 2015). Due to the slow dynamics brought about by the high viscosity of ionic liquids, the need for large timescales when conducting simulations in ionic liquids is one of the great obstacles, and unfolding times are often too large to be accessible by simulation (Tung and Pfaendtner 2016). The problem of insufficient sampling in pure ionic liquids is also mentioned by Burney and Pfaendtner in their study on the behavior of Candida rugosa lipase in [C4mim][NO3] and [C4mim][PF6]. Despite sampling times of more than 100 ns, no visible effects of the ionic liquid on the protein could be observed, despite experimental evidence that the enzyme is heavily denatured in [C4mim][NO3] (Burney et al. 2015). Tung and Pfaendtner address this problem of large timescales by introducing an enhanced sampling method called “infrequent metadynamics,” an algorithm that allows estimation of unfolding times on a greatly accelerated timescale (Tung and Pfaendtner 2016). For assessing protein stability, several structural quantities are accessible via MD simulations:

Radius of gyration

As unfolding increases protein volume, the loss of secondary structure leads to an increase of the radius of gyration \(\mathrm {r_{gyr}}(t)\)

$$ r_{gyr}(t) = \sqrt{\frac{\sum\limits_{i} m_{i} \left( \mathbf{r}_{i}(t) - \mathbf{r}_{\text{CM}}(t) \right)^{2}}{M}} $$
(11)

calculated with respect to the current center of mass \(\mathbf {r}_{\text {CM}}(t)\) of the protein, where \(m_{i}\) is the mass of the i th atom and M the total mass of the protein. Klähn et al. (2011a) reported that \(\mathrm {r_{gyr}}(t)\) of Candida antarctica lipase B (CAL-B) in aqueous IL mixtures (hydration level 20% w/w) increased by 0.8, 1.1, and 2.8 Å for [C4mim][PF6], [C4mim][BF4], and [C4mim][NO3], respectively. They suggested that the strength of interactions between anions and enzyme is mainly governed by ion density and ion surface charge, favoring interactions with smaller ions which can form stronger hydrogen bonds with the protein. This finding goes along with the denaturation capabilities of classical ILs in Fig. 5e.

However, Chaban and co-workers (Chevrot et al. 2015) showed that the radius of gyration does not increase in ionic liquids with amino acid-based anions during the simulation period of 10 ns which may not be long enough to monitor the unfolding process due to the high viscosity of these ionic liquids. Interestingly, Shao (2013) reported that rgyr of the B domain of protein A from Staphylococcus aureus actually decreased with increasing [C4mim]Cl content of the aqueous solution. Shao performed three independent simulations of 200 ns length to get meaningful values. A decrease of the protein radius was also observed by Bhattacharyya and co-workers (Ghosh et al. 2015) for lysozyme in [C5mim]Br. However, they computed the radius of hydration from the Stokes-Einstein relation:

$$ D = \frac{k_{B} T}{6 \pi \eta r_{H}} $$
(12)

using experimental diffusion coefficients from fluorescence correlation spectroscopy and the experimental viscosity \(\eta \). The radius of gyration and the radius of hydration are linked by a factor of \(\sqrt {5/3}\) (Schröder et al. 2015).

Root-mean-squared displacements

The time-dependent root-mean-square deviation RMSD (Kufareva and Abagyan 2012)

$$ RMSD =\left \langle \sqrt{\frac{1}{N} {\sum\limits_{i}^{N}} \left( \mathbf{r}_{i}(t) - \mathbf{r}_{i}^{\text{ref}}\right)^{2}} \right\rangle_{t} $$
(13)

is defined as a function of the deviation of the current atomic coordinates \(\mathbf {r}_{i}(t)\) from the coordinates \(\mathbf {r}_{i}^{\text {ref}}\) of a reference structure with N atoms, usually at the beginning of the simulation or the native configuration.

In the first MD simulation of a protein in an ionic liquid by Micaêlo and Soares mentioned above, the authors analyzed the RMSD of cutinase in [C4mim][PF6] and [C4mim][NO3] at different hydration levels and temperatures (Micaêlo and Soares 2008). In accordance with experimental findings reporting that many enzymes retain activity in ILs with PF\(_{6}^{-}\) anions and lose activity in ILs with NO\(_{3}^{-}\) anions (Lozano et al. 2001; Kaar et al. 2003; Persson and Bornscheuer 2003; Lou et al. 2006), they observed generally larger RMSD values in [C4mim][NO3]. In contrast, Burney et al. did not find a significant change in the RMSD of Candida rugosa lipase 1 in [C4mim][PF6] and [C4mim][NO3] (Burney and Pfaendtner 2013). Pfaendtner and co-workers (Burney et al. 2015) demonstrated that the RMSD is not only a function of the ionic liquid but also depends on the nature of the protein. For the enzyme Candida rugosa lipase 1, the stability decreases in the following order of the solvent: water \(>\) [C4mim]Cl \(>\) [C2mim][EtSO4]. However, no RMSD trend is visible for Bos taurus \(\alpha \)-chymotrypsin in the very same solvents (Burney et al. 2015).

The stronger destabilizing effect of nitrate compared to hexafluorophosphate was confirmed by Klähn et al. (2011a) for CAL-B in various aqueous ionic liquid mixtures by means of RMSD-values and agreed well with the concurrent observation of an increased radius of gyration \(r_{gyr}(t)\) mentioned in the last section. Although to less extent than the anions, longer alkyl chains lead to higher values of these quantities, indicating destabilization (Klähn et al. 2011b). Replacing butyl with a methoxyethyl group in C\(_{4}\textit {mim}^{+}\) also resulted in substantial protein destabilization, possibly by the increased hydrogen-bonding ability of the methoxyethyl group.

Diddens et al. (2017) analyzed the RMSD of a \(\beta \)-hairpin peptide by metadynamics simulations in aqueous mixtures of [C2mim][OAc], [C2mim][BF4], and [C2mim]Cl. The \(\beta \)-hairpin peptide unfolded in aqueous mixture of [C2mim][OAc] and [C2mim][BF4] as visible in the insets of Fig. 6. However, the global energetic minima of the two unfolded configurations in Fig. 6b and c are completely different. The authors argued that the strong denaturating acetate attacks the intramolecular hydrogen bonds of the peptide. As a result, the peptide structure unfolds completely. In aqueous [C2mim][BF4], the global minimum is still located near that of the pure water solution; however, the RMSD is quite large. Since the anions are excluded from the surface of the peptide, this effect was attributed to the impact of the cations.

Fig. 6
figure 6

Free energy landscape with global energetic minimum conformation of a β-hairpin. a Native state of the peptide in water. b Denaturated state in aqueous [C2mim]. c Denaturated state in [C2mim]. d Denaturated state in [C2mim]Cl. Re is the end-to-end distance of the peptide. (Diddens et al. 2017) - Published by the PCCP Owner Societies

Full size image

In contrast to the aqueous mixture of [C2mim][OAc] and [C2mim][BF4], the \(\beta \)-hairpin peptide did not unfold in aqueous [C2mim]Cl. The minimum structure of the peptide resembles the structure found in pure water (see Fig. 6a and d). This finding is in agreement with the experiment in Fig. 5b, where [C2mim]Cl was characterized to be only a weak denaturant.

In MD simulations of Shao (2013), the RMSD values of the B domain of protein A from Staphylococcus aureus decreased with increasing content of [C4mim]Cl which was interpreted as a stabilizing effect that may also be an effect of the increased viscosity.

Root-mean-squared fluctuations

The root-mean-squared fluctuation

$$ RMSF(i) = \sqrt{ \frac{1}{t} {\sum\limits_{0}^{t}} \left( \mathbf{r}_{i}(t) - \langle \mathbf{r}_{i}\rangle_{t} \right)^{2}} $$
(14)

measures the fluctuation of atom i around its mean position \(\langle \mathbf {r}_{i}\rangle _{t}\) (Pitera 2014). Experimentally, atomic flexibility is described with the crystallographic B-factor (Bornot et al. 2011). Amino acids with high RMSF values have an increased mobility which may also be an indicator for interactions with the ionic liquid (Attri et al. 2015; Burney and Pfaendtner 2013; Jaeger and Pfaendtner 2013) which breaks or at least disturbs intramolecular hydrogen bonding of the biomolecule.

For the analysis of simulations of three different cellulases in aqueous mixtures of the ionic liquid [C2mim][OAc], Jaeger et al. (2015) used the RMSD and an RMSF (averaged by residue) of the Cα atoms of the protein. They found that the three investigated cellulases (a cellulase from Thrichderma viride, a cellulase 5A from Thermotoga maritima, and an endoglucanase from Pyrococcus horikoshii) exhibited remarkably different stabilities in the aqueous ionic liquid [C2mim][OAc], highlighting that specific features of the protein, such as surface characteristics, and not only the ionic liquid itself determine stability.

Electrostatic interactions

Non-bonded energies in MD simulation are most generally calculated by the electrostatic energy \(U^{\text {elec}}(t)\)

$$ U^{\text{elec}}(t) = \frac{1}{4 \pi \epsilon_{0}} \sum\limits_{i} \sum\limits_{j} \frac{q_{i} \cdot q_{j}}{r_{ij}(t)} $$
(15)

between the partial charges \(q_{i}\) and \(q_{j}\) of the atoms i and j at a distance of \(r_{ij}(t)\) as well as van-der-Waals energies \(U^{\text {vdW}}(t)\)

$$ U^{\text{vdW}}(t) = \sum\limits_{i} \sum\limits_{j} 4 \epsilon_{ij} \left( \left( \frac{\sigma_{ij}}{r_{ij}(t)} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}(t)}\right)^{6} \right) $$
(16)

based on the interatomic Lennard-Jones parameter \(\sigma _{ij}\) and \(\epsilon _{ij}\).

Usually, the first atomic index i concerns all atoms in the biomolecule P under investigation. The second atomic index j may also comprise all atoms (separated by at least three bonds) in the very same molecule to obtain intramolecular energies \(U^{\text {elec}}_{\text {PP}}(t)\) and \(U^{\text {vdW}}_{\text {PP}}(t)\). Steinhauser and co-workers (Haberler and Steinhauser 2011) computed \(\langle U^{\text {elec}}_{\text {PP}} \rangle \) and \(\langle U^{\text {vdW}}_{\text {PP}}\rangle \) for a zinc finger protein in aqueous solutions of [C2mim][CF3SO3], and found a rise in \(\langle U^{\text {elec}}_{\text {PP}} \rangle (x_{\text {IL}})\) with increasing ionic mole fraction \(x_{\text {IL}}\), whereas \(\langle U^{\text {vdW}}_{\text {PP}} \rangle (x_{\text {IL}})\) decreased with increasing mole fraction. Moreover, the sum of electrostatic and van-der-Waals energies exhibited a maximum at an IL mole fraction of \(x_{\text {IL}}= 0.073\). The authors termed this mole fraction a “magic point” since other secondary structure descriptors also showed extrema at this point. This behavior can be attributed to a change from dipolar screening by water to charge screening by the ionic liquid ions when the IL mole fraction increases, which is consistent with a rise of \(\langle U^{\text {elec}}_{\text {PP}} \rangle (x_{\text {IL}})\).

Klähn et al. analyzed the electrostatic and van-der-Waals-interaction energy of hydrated Candida antarctica lipase B (CAL-B) with a set of ionic liquids composed of imidazolium- and guadinium-based cations and nitrate, tetrafluoroborate or hexfluorophosphate anions (Klähn et al. 2011a, 2011b). Here, the second atomic index j of Eq. 15 contains all solvent atoms of the particular solvent species. They found that interaction energies are dominated by Coulombic interactions \(\langle U^{\text {elec}}_{\text {PA}} \rangle \) between the anions (A) and the protein, while smaller cation (C) contributions stem from van-der-Waals-interactions \(\langle U^{\text {vdW}}_{\text {PC}} \rangle \). The charge on cations is more delocalized while the charge density on the smaller anions is higher, allowing them to form stronger hydrogen bonds with the protein. In their study on xylanase II in [C2mim][EtSO4] and [C2mim][OAc], Jaeger and Pfaendtner also attributed the low affinity of cations for negatively charged residues to the high charge delocalization (Jaeger and Pfaendtner 2013). Klähn et al. reported the highest values of 〈UvdW〉 + 〈Uelec〉 (strongest hydrogen-bond interactions) for the NO\(_{3}^{-}\) anion, which is consistent with experimental and simulation results that NO\(_{3}^{-}\) destabilizes CAL-B (Klähn et al. 2011a; Micaêlo and Soares 2008). Weakest interactions were observed for the PF\(_{6}^{-}\) anion, which stabilizes CAL-B. From experiment, it is also known that [C4mim][NO3] can dissolve solid CAL-B, while [C4mim][PF6] cannot (Klähn et al. 2011a; Lau et al. 2004).

Solvation free energy

The general procedure of solvation free energy calculation is briefly described in the next section on small molecules. As already mentioned before, the simulations of proteins in ionic liquids are hampered by the increased viscosity in ionic liquids, requiring long simulation times for sufficient sampling (Tung and Pfaendtner 2016). Under these circumstances, it is impossible to acquire sufficient sampling for the calculation of solvation free energies (Klähn et al. 2011a). However, some other methods for obtaining thermodynamic information about proteins in ionic liquids have been employed. Klähn et al. (2011a) calculated the solvation free enthalpy \(H_{\text {solv}}\) of CAL-B based on a linear response formalism (Klähn et al. 2011a),

$$ H_{\text{solv}} = \frac{1}{N^{\text{solute}}} \left( U_{\text{solv}} - \left( U^{\text{solute}}_{\text{gas}} + U^{\text{solvent}} \right) + P {\Delta} V \right), $$
(17)

where \(U_{\text {solv}}\) is the potential of the system containing \(N^{\text {solute}}\) solvated solutes in \(N^{\text {solvent}}\) solvent molecules, and \(U^{\text {solute}}_{\text {gas}}\) and \(U_{\text {solvent}}\) are the potential energies of the N solutes in gas phase and the pure solvent, respectively. The pressure-volume term is small and may be neglected, and when \(N^{\text {solute}}= 1\) is assumed, Eq. 17 may be written as

$$\begin{array}{@{}rcl@{}} H_{\text{solv}} &=& U^{\text{solute}}_{\text{solv}} + U^{\text{solvent}}_{\text{solv}} + U^{\mathrm{solute-solvent}} \end{array} $$
$$\begin{array}{@{}rcl@{}} &&- \left( U^{\text{solute}}_{\text{gas}} + U^{\text{solvent}} \right) \end{array} $$
(18)
$$\begin{array}{@{}rcl@{}} &=& {\Delta} U_{\text{cage}} + {\Delta} U^{\text{solute}} + U^{\mathrm{solute-solvent}}. \end{array} $$
(19)

Here, \(U^{\text {solute}}_{\text {solv}}\) is the potential energy of only the solute in solvent, \(U^{\text {solvent}}_{\text {solv}}\) is the potential energy of only the solvent in the system, and \(U^{\mathrm {solute-solvent}}\) is the solute-solvent interaction energy. \({\Delta } U_{\text {cage}} = U^{\text {solute}}_{\text {solv}} - U^{\text {solvent}}\) is termed the cage formation energy and measures the energy required to form a solute cavity in the solvent. ΔUsolute = U solvsolvent − U gassolute is the change in internal energy of the solute upon insertion into solvent (Klähn et al. 2011a). It was found that solvation enthalpies in ILs increased noticeably in ionic liquids compared to water, which is in accordance with experimental observations that CAL-B is less soluble in ionic liquids compared to water. It was also observed that these high solvation enthalpies were mainly due to large cage formation energies. These large cage formation energies are due to the strong interactions between ionic liquid ions compared to the weaker interactions in water. An exception was the ionic liquid [C4mim][PF6], in which the cage formation energy was even lower than in water. Furthermore, upon transfer from water to IL, the internal energy of the protein increased, except in [C4mim][PF6], corresponding to the experimentally observed stability increase in [C4mim][PF6]. Another approach to estimate the free energy landscape of a protein is the use of metadynamics (Barducci et al. 2008), which was already mentioned in the context of a small \(\beta \)-hairpin protein whose unfolding free energy landscape in aqueous [C2mim][OAc], [C2mim][BF4], and [C2mim][Cl] was calculated via metadynamics (Lesch et al. 2015; Diddens et al. 2017). As mentioned earlier, the strong denaturation of the peptide by [C2mim][OAc] apparent from preferential binding parameters was also visible in the free energy landscapes.

Small rigid drugs

Ionic liquids are of special interest for the extraction of compounds from biomass. They can dissolve a wide range of biomatrices, even materials such as cellulose (Pinkert et al. 2009; Wang et al. 2012; Zhu et al. 2006) or chitin (Qin et al. 2010), and can enhance solubility of hydrophobic compounds in aqueous solution by acting as hydrotropes (Cláudio et al. 2015) or surfactants.

From the computational point of view, these molecules are much more rigid compared to a large flexible protein. Consequently, observables concerning size and shape, e.g., the radius of gyration as well as root-mean-squared displacements and fluctuations, are of less use for the small biomolecular drugs dissolved in the aqueous mixtures of ionic liquids. Hence, particular interactions like \(\pi \)-π-stacking or hydrogen bonding attract more attention. Also, free energy calculations, which are not feasible for complete proteins, may help to understand the biomolecular solvation of smaller solutes. However, the current knowledge on the interaction of small biomolecules with ILs is not as far progressed compared to protein solvation.

Hydrotropic theory

The hydrotropic effect describes the increased water solubility of hydrophobic organic compounds by the addition of a co-solvent (Shimizu and Matubayasi2014; Eastoe et al. 2011; Cláudio et al. 2015; Zeindlhofer et al. 2017, 2018) which usually consists of small, amphiphilic molecules (Booth et al. 2012). In principle, the following mechanisms are discussed:

  1. 1.

    The solute forms a stable complex with few hydrotropic molecules via hydrogen bonding or \(\pi \)-π-stacking. The complex has an increased polarity, hence increased solubility in water.

  2. 2.

    The hydrotrope assembles in micelle-like structures in which the solute can be accommodated.

  3. 3.

    The water-immiscible solute induces a particular structure of the hydrotrope in the aqueous solution. This new structure promotes the solubility of the solute.

Imidazolium-based ILs seem to be promising candidates for hydrotropy. Their charged, hydrophilic ring is accompanied by more or less hydrophobic alkyl side chains (Tan et al. 2009; Cláudio et al. 2015). We recently reported (Zeindlhofer et al. 2018) on the solvation of coffee ingredients in [C2mim][OAc] mixtures. The cations seem to repel water in general and also anions from the surface of the solute at higher concentration. Even in the second solvation shell determined by Voronoi tessellation, the concentration of water is below bulk density which argues against mechanism 1. Nevertheless, water still has access to the solute surface, even at higher concentrations, and micellar cationic structures are not observed at the surface. Also, the last hydrotropic hypothesis 3 could not be proven. The inclusion of acetate in the water network is usually not a problem and does not change the water structure in the proximity of the solute very much. However, a Voronoi shell-wise computation of the dielectric permittivities reveals a subtle transition from the solute permittivity to the permittivity of the bulk solution (Zeindlhofer et al. 2018).

Furthermore, hydrogen bonding between the water and the solute is not changed by the cations (as suggested by hypothesis 3) but in competition with the anions. Interestingly, the solvation efficiency drops from tosylate to chloride for gallic acid (Cláudio et al. 2015): Cl− <[MeSO\(_{4}^{-}\)] \(<\) [CF\(_{3}\textit {SO}_{3}^{-}\)] \(<\) [N(CN)\(_{2}^{-}\)] \(<\) [Tos−]. A similar ranking was found in vanillin. Although the cations are the most prominent species at the solute surface, particular anion interactions are very important for the solvation behavior of the ionic liquid, i.e., hydrogen bonding. Cationic \(\pi \)-π-stacking is of minor importance.

Solute-solvent interactions and solute aggregation

The extraction of shikimic acid from star anise pods using [C2mim]-based ionic liquids with a set of different anions (PF\(_{6}^{-}\), NTf\(_{2}^{-}\), BF\(_{4}^{-}\), Cl−, OTf−, OAc−) was reported in an experimental study (Zirbs et al. 2013). Accompanying MD simulations were performed to identify the relevant molecular interactions in the extraction process: No correlation of the experimental extraction yield with the coordination number or the contact surface obtained from simulations could be established (Zirbs et al. 2013). Coordination numbers around shikimic acid are higher for cations, but anions seem to interact stronger since they are able to form hydrogen bonds, and their residence time in the first solvation shell around shikimic acid is longer than for cations. The number of anion hydrogen bonds correlates with the extraction yield, highlighting the importance of the anions in the extraction process.

In agreement with electrostatic analysis of the IL ion with proteins (Klähn et al. 2011a; Klähn et al. 2011b), \(\langle U^{\text {elec}}_{\text {SA}} \rangle \) of the anions (A) with the solute (S) is much larger compared to \(\langle U^{\text {elec}}_{\text {SC}} \rangle \) of the cations (C) (Zeindlhofer et al. 2017; Tomé et al. 2012). Furthermore, the van-der-Waals interaction \(\langle U^{\text {vdW}}_{\text {SC}} \rangle \) between cations and the solute exceeds the corresponding Coulomb interaction \(\langle U^{\text {elec}}_{\text {SC}} \rangle \) for C\(_{2}\textit {mim}^{+}\) at the surface of various coffee ingredients (Zeindlhofer et al. 2017). Cations are mainly found parallel to the aromatic solute rings in a stacked orientation, while anions preferred positions around the ring that allowed hydrogen bonding to ring substituents as visible from \(\langle U^{\text {elec}}_{\text {SA}} \rangle \). The same finding was observed for five different amino acids (glycine, alanine, valine, isoleucine, glutamic acid) in aqueous solutions of [C4mim][NTf2] (Tomé et al. 2012).

Another experimental study complemented by MD simulations concerned the extraction of vanillin in aqueous solutions of [C4mim][N(CN)2] and [C4mim][SCN] (Cláudio et al. 2015). In MD simulations, the ionic liquid ions form several filamentous ionic strands and interact with water in an anion-water hydrogen bond network. The hydrophobic vanillin molecules form a large droplet in pure water as visible in Fig. 7c, which is broken up into smaller aggregates in the aqueous IL solution by formation of larger vanillin-cation clusters and smaller vanillin-anion clusters in Fig. 7d and e.

Fig. 7
figure 7

Distribution functions P(na) of aggregate sizes in pure aqueous ionic liquid mixtures (a, b) and with vanillin (d, e). c depicts a vanillin-water mixture. Reproduced from Cláudio et al. (2015) with permission of The Royal Society of Chemistry

Full size image

Solvation free energy

Enthalpic and entropic effects of the dissolution of a solute in a solvent can be best characterized by the computation of solvation free energies. In contrast to large proteins, the amount of sampling needed for these calculations is feasible for small compounds. A common approach is thermodynamic integration. During several simulations, a coupling parameter λ is stepwise increased from zero to one corresponding to a switching off of solute-solvent interactions \(U^{\mathrm {solute-solvent}}\)

$$ U(\lambda) = (1-\lambda) \; U^{\mathrm{solute-solvent}} + \lambda \cdot U^{\text{solvent}} $$
(20)

Usolvent corresponds to the total energy of the solvent only. The free energy difference between the initial and final \(\lambda \)-state is given by

$$ {\Delta} G_{\text{solv}} = \int \left\langle \frac{\partial U(\lambda)}{\partial \lambda} \right\rangle_{\lambda} d\lambda $$
(21)

To study molecular details of extraction of amino acids from aqueous solution using ILs (Absalan et al. 2010; Tomé et al. 2012), Seduraman et al. (2012) calculated the solvation free energy of L-tryptophan in pure and aqueous [C4mim][PF6], [C8mim][PF6], and [C8mim][BF4]. In contrast to experimental studies, they found only negligible differences in the free energies in the pure ILs. However, including residual water in the IL simulations, which is higher in [C4mim][PF6], the experimental trend is reproduced illustrating the importance of water for the extraction process. The solvation in BF\(_{4}^{-}\)-based ILs is more favorable, probably due to the fact that the water density is 7.5 times higher compared to the hexafluorophosphate systems. Furthermore, electrostatic interactions of tryptophan with BF\(_{4}^{-}\) are stronger than with PF\(_{6}^{-}\). This anion dependency was observed for other amino acids as well (Absalan et al. 2010).

Conclusion

After the initial hype in the beginning of the twenty-first century, expectations on ionic liquids decreased because of issues with the price, viscosity, the stability in water, the toxicity, and biodegradability. Often, the unique selling point of ionic liquids was missing in applications to justify the higher cost compared to conventional solvents. Moreover, exorbitant promises during the hype prompted several authors to question the current status of knowledge and the usefulness of IL in general (Kunz and Häckl 2016). Consequently, to prevent making mistakes again, which were made during the hype, good practice routines should go without saying. In the current review, we briefly summarized meaningful computational methods to investigate the solvation of large and small biomolecules in aqueous ionic liquid mixtures.

The increasing number of publications concerning ionic liquids demonstrates their unbowed interest. However, in order to make real progress during the implementation of ionic liquids and their computational analysis, several ideas should be taken into account:

  1. 1.

    The vast majority of computational works deal with aprotic imidazolium-based ionic liquids. Maybe, protic imidazoliums as well as other cations like choline may be more biodegradable and cheaper. Therefore, juxtaposition of computational results of these “new” ionic liquids to the aprotic imidazolium-based ionic liquid are interesting.

  2. 2.

    The same is true for the anions. Tetrafluoroborate and hexafluorophosphate can easily be coarse-grained and much experimental data exists, but these anions are of less interest for the experimenter. Therefore, other anions such as dihydrogenphosphate (dhp), methyl- and ethylsulfate, hydrophobic bis-(trifluoromethylsulfonyl)-imide NTf\(_{2}^{-}\), or hydrophilic dicyanamide N(CN)\(_{2}^{-}\) should also be included in the portfolio of the theoreticians. Even deep eutectic solvents are interesting but currently still mostly neglected in MD simulations on solvation of biomolecules.

  3. 3.

    Many solvation effects depend on the concentration of the co-solvent and should be given much more room in current theoretical investigations as these simulations can be simply parallelized.

  4. 4.

    Particular effects of ionic liquids should be tested on more than one protein as the observation made may be a peculiarity of the investigated protein.

  5. 5.

    Polarizable force fields to include electronic effects may also be an option.

  6. 6.

    The computational setup is of utmost importance for the reliability of the MD results. Hence, simulation system sizes should be large enough to prevent artifacts from self-interactions. For example, at least five complete solvation shells around the biomolecule should be considered. The simulation period should be at least 50 ns in low viscous systems, i.e., mixtures with high water content, and even much longer with increasing ionic liquid concentration. Of course, the simulation period has also to be prolonged when analyzing folding/unfolding processes. Here, smaller proteins may be more useful than larger proteins to cover the complete process or at least major parts of it.

Applying these suggestions to future simulations on biomolecular solvation in ionic liquids should help to reach the “Plateau of Productivity” within reasonable time.