Martini Coarse-Grained Models of Imidazolium-Based Ionic Liquids: From Nanostructural Organization to Liquid-Liquid Extraction

Ionic liquids (IL) are remarkable green solvents, which find applications in many areas of nano- and biotechnology including extraction and purification of value-added compounds or fine chemicals. These liquid salts possess versatile solvation properties that can be tuned by modifications in the cation or anion structure. So far, in contrast to the great success of theoretical and computational methodologies applied to other fields, only a few IL models have been able to bring insights towards the rational design of such solvents. In this work, we develop coarse-grained (CG) models for imidazolium-based ILs using a new version of the Martini force field. The model is able to reproduce the main structural properties of pure ILs, including spatial heterogeneity and global densities over a wide range of temperatures. More importantly, given the high intermolecular compatibility of the Martini force field, this new IL CG model opens the possibility of large-scale simulations of liquid-liquid extraction experiments. As examples, we show two applications, namely the extraction of aromatic molecules from a petroleum oil model and the extraction of omega-3 polyunsaturated fatty acids from a fish oil model. In semi-quantitative agreement with the experiments, we show how the extraction capacity and selectivity of the Abstract Ionic liquids (IL) are remarkable green solvents, which find applications in many areas of nano- and biotechnology including extraction and purification of value-added compounds or fine chemicals. These liquid salts possess versatile solvation properties that can be tuned by modifications in the cation or anion structure. So far, in contrast to the great success of theoretical and computational methodologies applied to other fields, only a few IL models have been able to bring insights towards the rational design of such solvents. In this work, we develop coarse-grained (CG) models for imidazolium-based ILs using a new version of the Martini force field. The model is able to reproduce the main structural properties of pure ILs, including spatial heterogeneity and global densities over a wide range of temperatures. More importantly, given the high intermolecular compatibility of the Martini force field, this new IL CG model opens the possibility of large-scale simulations of liquid-liquid extraction experiments. As examples, we show two applications, namely the extraction of aromatic molecules from a petroleum oil model and the extraction of omega-3 polyunsaturated fatty acids from a fish oil model. In semi-quantitative agreement with the experiments, we show how the extraction capacity and selectivity of the IL could be affected by the cation chain length or addition of co-solvents. C4: 1-butyl-3-methylimidazolium tetrafluoroborate, C8: 1-octyl-3-methylimidazolium tetraflouroborate, C12: 1-dodecyl-3-methylimidazolium tetrafluoroborate, CG :coarse-grained, DHA: docosahexaenoic acid, IL: ionic liquid, OLE, Oleic Acid, PAL: Palmitic Acid, RDF: radial distribution function.


Introduction
Ionic liquids (IL) have been gaining importance during the last few years because of the many possible applications in different areas in chemistry like electrochemistry, chemical synthesis, catalysis, pharmaceutics and medicine, separation and extraction of fine chemicals and nanotechnology [1][2][3][4][5][6][7] . These unusual substances are a special class of liquids composed entirely of ions with melting points below 100 °C. Generally, an IL consists of one large organic cation and an organic or inorganic anion. Among the existing ILs, the most commonly used are the ones based on non-symmetrically substituted dialkylimidazolium cations and bulky anions 8,9 . ILs are also known as "designer solvents" because one can modify their physical properties by changing the anionic and cationic components 8 .
One of the most prominent applications of IL is in separation science, as ILs have characteristics like unusual selectivities, high extraction efficiencies, durability and resistance to thermal degradation which makes them ideal for performing extractions 10 . ILs have found a major application field in the extraction and separation of bioactive compounds 1011 . One of the most widely investigated IL for extractions is the cation 1-alkyl-3-methylimidazolium [CnImin] + in combination with [BF4] -, [PF6] -, Cl -, or Bras anions 10 .
Understanding the structure of ILs is one of the cornerstones for the development and improvement of applications of ILs 2,8,9,12 . In this direction, computational modelling, in particular through the use of molecular dynamics simulations, has proven to be a powerful technique to obtain insights into the nanostructure and dynamics of ILs [13][14][15][16][17][18] . However, the use of full atomistic simulations of IL is computationally intensive 13,[18][19][20] . Consequently, there is a need to develop coarse-grain (CG) models 11,14 that can address the problem by simulating ILs at a low computational cost but still capturing the details of the structure and properties of these polar solvents. A number of bottom-up CG approaches have been developed (e.g. Voth, 2005 21 , Voth, 2009 22 , Carbone, 2010 23 , Rotenberg, 2012 24 , Holm, 2017 25 and Aluru, 2018 26 ) that are capable of accurately reproducing reference data obtained from all-atom simulations. However, most of these CG force fields are specific for certain ILs or are not compatible with other molecules 16 , such that these models cannot be easily extended or adapted to describe complex mixtures with other compounds. Alternatively, top-down coarse-grained force fields such as Martini are more amenable to 'mixing-and-matching' several classes of compounds.
The Martini CG force field developed by Marrink and co-workers 27,28 is one of the most popular CG force fields for simulations of biomolecules 29,30 and nowadays is extending its applications to soft materials science [29][30][31][32] . This force field is constructed in a top-down approach by extensive calibration of the non-bonded interactions by comparison against thermodynamic data 29 . Because of the way it is constructed, Martini can be used in a broader range of applications without the necessity of reparametrize the model for every specific application. Recently the limits of the current version 2.2 of the Martini model have been reviewed 29,33 , which paved the way for developing a new version, Martini 3.0. Currently, only three classes of ILs models were developed for Martini 2.2, with the studies focused on pure IL, aqueous biphasic mixtures 34,35 and interactions with membranes 36 .
In this paper, we elaborate and test a new CG model for imidazolium-based IL in the framework of the Martini 3.0 force field, and we investigate applications of our model on the extraction of fine chemicals. The rest of the paper is structured as follows. The next section describes the concept of our CG model and methods used for the construction and simulation of pure IL and biphasic systems used for extractions. Thereafter, the results are presented, with the first part mainly devoted to description of the structure of IL in relation to experimental data, while the second part is focused on two extraction simulations: (i) extraction of aromatic compounds from a petroleum oil model system; (ii) extraction of omega-3 polyunsaturated fatty acids from a fish-oil model. The final section draws some conclusions and presents future perspectives for the simulation of these systems.

Simulation details
All simulations described were performed using GROMACS 5.1 37 and the open beta version of the Martini 3.0 CG force field 38 . We tested long-range Coulombic interactions by using both a particle mesh Ewald method (PME) 39 and a reaction field method 40 . The cutoff distance was defined as 1.1 nm for Coulombic and Lennard-Jones interactions as suggested by de Jong et al 41 . The Verlet neighbour search algorithm 42 was used in combination with the neighbour list updated every 20 steps. We use the LINCS algorithm 43 to constrain the bonds and the leapfrog integration algorithm for the solution of the equations of motion. The temperature was controlled using the vrescale thermostat 44 with a coupling time of 1 ps. For pressure coupling, the Parrinello-Rahman 45 barostat was used (compressibility 3x10 -4 bar -1 ). All initial configurations were built using GROMACS tools and the PACKMOL software 46 .

CG models
For the imidazolium IL type [CnImin] + [BF4] -, henceforth called Cn, the CG models described in Figure 1 were constructed. The substituted imidazolium rings were initially modelled with three beads, using the bead type TC6 for the nitrogen-containing aromatic fragments while the rest of the carbon-based aromatic ring is represented by a TC5 bead. Tiny bead sizes were used to better mimic the overall shape and packing of the molecule, which is important for the stacking of the imidazolium cation with other aromatic compounds. From this template, that represents the C1 cation, the different tail lengths were modelled. For C2, the size of one of the TC6 beads was increased to its small "S" version (SC6), which includes an ethyl group attached to one of the nitrogen atoms. For C4, an SC3 bead was added to represent the additional n-propyl fragment while the subsequent tail lengths in C8 and C12 were obtained adding a regular size C1 bead for each additional n-butyl fragment (which is analogous to the strategy used to build different lipid tail lengths in Martini 47 ). A representation of the coarse-grained models of imidazolium cations is shown in Figure 1.
In contrast to most Martini models, partial (instead of integer) charges are used to describe the Cn cations. The values of the charges in the imidazolium ring (+0.5 on each nitrogen-based bead) were computed according to a variant of the Dipole Preserving Charge 48 method described in the Supplementary Material. Based on the open-beta release of Martini 3.0, the Lennard-Jones interactions of beads in the imidazolium ring are slightly increased with themselves (in one interaction level) and with the rest of the neutral beads (two interaction levels) compared to neutral versions of the same bead type. At the moment, this is implemented by introducing bead-types carrying the label "h", which will be updated in the final release of Martini 3. The bond lengths of the imidazolium ring were tuned to get the best compromise in terms of bulk density and molecular surface area. The parameters for the tails are inspired by the lipid CG model 47 . All structural parameters are described in the supplementary material (Table S1).
The benzene molecule used is constructed with three TC4 beads with a constrained bond between them of 0.291 nm. For octane and the fatty acids of the fish oil, we use the same bonded parameters and mapping from the previous version of Martini 27,28 , adapting some bead types to the new ones of the beta version of Martini 3.0. For the neutral carboxylic heads of the fatty acids, we use a P2 bead. All the parameters are available on our web portal cgmartini.nl.

System setup
For simulations of pure ionic liquids, a cubic box of approx. 6.4 x 6.4 x 6.4 nm 3 was used. ILs of C2, C4, C8, and C12 ( Figure 1) were each simulated over a temperature range of 293-393K. The number of molecules for each system is described in Table S2 in the supplementary material. The energy of the system was minimized for 500 steps using the steepest descent algorithm. The system was then equilibrated in two steps: 500 ps at 298 K and 100 ns at the target temperature. The length of the production simulation was 400 ns. For the simulations used to study the temperature-driven phase transition in C12, a bigger and rectangular box (12.7 x 12.7 x 6.4 nm 3 ) was also used, avoiding possible artefacts induced by periodic boundary conditions. For these simulations, the initial configuration was built in a lamellar form, resembling a smectic A phase. For comparison of structural properties, atomistic simulations of pure C2, C4, and C8 ILs were performed according to the procedure described in the supplementary material.
Biphasic systems containing two liquid phases, namely an oil and an IL, were simulated at 300 K and 1 bar in the NPzT ensemble, with X and Y dimensions fixed and pressure applied in z-direction (normal to the liquid/liquid interface). The energy of the system was first minimized using 5000 steps. Equilibration of the system was performed using a time step of 5 fs during 1 ns. Production was generated using a time step of 20 ps during 4.4 μs. The ILs considered are C2, C4, C8, and C12 ( Fig. 1). For the extraction simulation of aromatic compounds from a model of petroleum oil, the system was built as a mixture of octane and benzene at a weight ratio of 90:10 49 . Subsequently, this model system was mixed with IL, adding the same number of cation and anion pairs as the number of octane molecules. For the construction of the system with fish oil, we follow an approximate composition of a model fish oil from Chinook Salmon according to the data reported by Stansby et. al. 50 The composition of the fish oil is 48% of Oleic Acid (OLE), 48% of Palmitic Acid (PAL) and 4% of Docosahexaenoic Acid (DHA) in weight percentage. The biphasic systems were constructed in a proportion of 2:1 IL/fish oil in volume percentage. To test the stability of the biphasic system through the addition of co-solvents, we add octane in the same volume percentage of IL as suggested in experimental tests 51 . The number of particles and the size of the boxes used for all systems modelated are described in the supplementary material (Table S2).

Analysis for pure ionic liquids
Radial Distribution Functions were computed using GROMACS tool gmx rdf 37 . Snapshots of the simulations were generated using the VMD (Visual Molecular Dynamics) suite 52 . To characterize long range order in the C12 system, we computed the order parameter of the alkyl tail defined as: where θ is the angle between the molecular axis (the vector formed by the connection of the beads of the alkyl tail) and the z-axis which was used as a reference. The order parameter was computed using a modified version of the script found in cgmartini.nl 53 .
Heat capacity was computed as the derivative of the energy with respect of the temperature as: Where, is the change in the total energy and is the change in the temperature.

Oil/IL transfer free energy and selectivity calculations
For the biphasic systems which conventional MD simulations showed enough sampling of the all the relevant components in both phases, the oil/IL transfer free energy where R denotes the gas constant, T the temperature and Dx distribution coefficient of the x component. Errors in &'()* were calculated using the block averaging. Dx is defined in equation 4 as the ratio of the average partial density ( ,-) of the compound x in the zone where the normalized partial density of IL is close to one and the same quantity in the zone of the oil phase ( ./0 ). Partial densities were computed using GROMACS tool gmx density 37 .
For systems where more sampling was required, the &'()* was obtained from umbrella sampling (US) simulations 54 , extracting the free energies with the weighted histogram analysis method (WHAM) 55 , as implemented in the GROMACS tool gmx wham 37 . In this case, errors were estimated using bootstrapping. The reaction coordinate was the distance between the centers of one molecule of component x and the IL phase; 56 windows spaced 0.1 nm apart were used, with a simple harmonic umbrella potential with a force constant of 2000 kJ·mol -1 ·nm -2 . The sampling time for each window was 100 to 200 ns, depending on the convergence. The same biphasic systems (in terms of box size and composition) used for the conventional MD simulations were also used for the US.
To quantify the selectivity of the extraction, we use the following expression 51,56 : Where + is the distribution coefficient for the desired component to extract and / are the distribution coefficients for the other components on the mixture.

Pure ionic liquids
Experimentally, it has been documented that the density of imidazolium-based ILs decreases as the temperature rises 57,58 . Density is also affected by the length of the alkyl tail of the imidazolium cation, with long alkyl tail cations showing lower density than the counterparts with shorter alkyl tails 57,58 . The results of our MD simulations with the pure IL indicate that the density profiles obtained of the ILs studied (Figure 2A and Figures S1-S2) follow these experimental trends. In all cases, we observe that as the temperature decreases, the agreement of the densities obtained with our model and the experimental values improves. There is not a notable difference in the density if PME or reaction-field are used for the description of the electrostatic interactions (see Figures S3-S5 and Tables S3 and S4). This was also observed for all the other systems simulated in this work. So, from this point on, we will only describe the results using PME.
One crucial aspect that is well described by our model is the formation of nanodomains in the ILs. We can see in Figure 2B that as we increase the tail of the alkyl fragment, the organization of the liquid becomes more heterogeneous. For C12, we observe the formation of a smectic phase. The formation of domains as a function of the length of the alkyl tail of the imidazolium cation was further quantified using the radial distribution function (RDF). In the RDF profiles between the centre of mass of the alkyl tails and itself ( Figure 2C), we observe that the size of the domain increases from C4 to C8 and has a maximum for C12. Despite the changes in long range organization, the short-range structure of the IL charged groups is conserved for different cations, as evidenced by the RDFs between the charged groups ( Figures S6-S8). The RDFs between the imidazolium ring and the anion are also in close agreement to those obtained from our reference atomistic simulations ( Figure 2D), as well as to those obtained with another CG model by Wang et al. 22 . However, it is also evident that the aggregation of nonpolar tails affects the local structure to some extent, in particular a closer approach of the anions to the imidazolium ring as the tail length increases ( Figure 2D). These features clearly imply that our model can describe the spatial heterogeneity of the ILs. It is known from the literature that the range of temperatures simulated in this work includes a phase transition for C12 59,60 . To study this phase transition in more detail, additional simulations of the C12 system were performed with a bigger simulation box (see Methods). On the analysis of the RDFs of the cation C12 IL at different temperatures, we observed a change in the longrange structure of the IL (Figures 3A and S9). In particular, we observe a change of the intensity of the second peak (around 1.05 nm) of the RDFs between the alkyl tails, which is lower at a high temperature (see insert, Figure 3A). A plot of the intensity of the 2nd peak of the RDF with respect to temperature ( Figure 3B) clearly shows a transition occurring around 330 K. The phase transition was further quantified by measuring the order parameter of the alkyl tail with respect to the system's z-axis ( Figure 3C). From this data, it is clear that the phase transition of C12 is accompanied by a change in the orientation of the alkyl tails of the C12 cation, passing from a smectic (order parameter around 0.3) to a completely isotropic orientation (order parameter close to zero). Finally, in order to obtain a more accurate estimate of the temperature at which this phase transition occurs, we analyze the derivative of the energy with respect to the temperature, which corresponds to the heat capacity of the system ( Figure 3D). From this figure, it is immediately clear that the phase transition from the smectic phase to the isotropic phase occurs at 325 K, a value that is in very good agreement with the experimental reported value for C12 of 324.75 K 60 . Snapshots from the simulation box at different temperatures ( Figure 3E) also provide visual evidence for the phase transition. Next to the structure and phase behavior, we analyzed another commonly studied propriety of ILs, namely the diffusion coefficient. This dynamical property gives us an insight into the mobility of cations and anions in the ionic liquid. The existing coarse-grain models used to study ILs 26,61 , usually overestimate the value of the diffusion coefficients by several orders of magnitude. This is a natural consequence of the reduction of the degrees of freedom of the system. In the cases studied in this work the obtained diffusion coefficients are of the order 10 -5 cm 2 s -1 ( Figure S10-S11), between 2 to 3 orders of magnitude bigger than the values obtained from atomistic simulations 20,62 . As expected, the diffusivity of the ions increases with temperature, and decreases with the length of the alkyl tails. For short tail ILs (C2-C4), the cations and anions diffuse at roughly the same speed (Figures S10-S11) reproducing the experimental trend 19 and results obtained at the atomistic level 19 and reflecting the behavior that the anion prefers to stay close to the positively charged imidazolium ring. For the longer tailed ILs, however, the anion is diffusing relatively fast compared to the cation, presumably because the larger mass of the long tail cation. In addition, the anion becomes increasingly trapped inside the nanodomains (or lamellae in case of C12) as can be observed in Figure 2B.

Petroleum oil extraction
ILs have been proposed as an alternative of conventional solvents for the extraction of aromatic compounds from mixtures of them with aliphatic species 63,64 . This proposal is based on the fact that the solubility of aliphatic compounds in ILs can be modulated by making changes on the structure of the cation or anion; meanwhile, the solubility of aromatic substances in ILs is always considerable 65 . Additionally to these advantages, it has been proven that the use of IL on the extraction of aromatic compounds is economically feasible 64 . Selectivity and the distribution ratio of the aromatic molecules between the IL and the aliphatic solvent are the main criteria for the economic sustainability of the extraction process 65 of aromatic compounds. These parameters depend on many factors such as the type of hydrocarbon compounds in the mixture, the type of IL used for the extraction, the length of the alkyl tail on the cation, the temperature of the process, among other considerations 66 . For instance, an increase in the length of the alkyl tail on the cation results in an increased solubility of aliphatic molecules, but a reduced selectivity 65 . Furthermore, ILs can be used for the extraction of aromatic molecules from mixtures in which those compounds are in low concentration with a good selectivity 66 .
In order to test our model in the simulation of the extraction of molecules from complex mixtures, we build a simple system analogous to petroleum oil to analyze the extraction of aromatic molecules, following a methodology which was previously used experimentally 49 . Our petroleum oil model consists of octane with 10% (weight) of benzene. The oil phase was put into contact with a specific IL phase (either C2, C4, C8, or C12), creating a biphasic system as shown in Figure  4 and Figures S12-S13. From this system, we evaluate the efficiency of the extraction of benzene molecules to the IL and the contamination by the migration of octane to the IL. We did not consider the effect of increased temperature as it is known that an increase in the temperature of the mixture reduces the selectivity on the extraction of aromatic molecules 63 . To analyze the extraction efficiency of the various ILs considered, we computed the normalized partial density profiles of the IL/oil mixtures, as shown in Figures 4A, 4B and 4C for C2, C4 and C8, respectively. Additional results for C12 are shown in Figure S12. Except for C12, all IL/oil systems remain clearly biphasic. In case of C12, the solubility of octane in the IL is considerable and the phase separation becomes less clear (See Figures S12 and S13). In terms of benzene extraction, comparison of the normalized densities ( Figure 4) as well as the transfer free energy of benzene between the oil phase and the IL ( Table 1), indicates that the longer tail ILs extract the largest fraction of benzene. However, also the fraction of octane that is extracted increases with the size of the alkyl tail (Figure 4), making the long tail ILs less selective for the extraction of benzene. The selectivity is quantified by computing the ratio of the distribution coefficients of benzene over octane in the IL (Table 1). From this data it appears that C4 would be the optimal choice for extraction of benzene, with both a high extraction capacity and high selectivity for benzene. In summary, with the growth of the hydrocarbon tail, our CG simulations show that more benzene can migrate to the bulk of the IL. At the same time, the selectivity starts to be reduced, and part of the octane migrates to the IL phase. This is exactly the behaviour experimentally shown for C4, C6 and C10 IL/octane mixtures containing 10% of benzene 49 . Consequently, the extraction efficiency and selectivity of benzene can be accurately captured by our new Martini IL model. Except for C12, we also noticed that the molecules of benzene have the highest concentration at the oil/IL interface. In our view, this behaviour suggests that the full separation of benzene from octane using imidazolium-based IL may need experimental protocols using multiple extraction steps.

Fish oil extraction
Encouraged by the results obtained for the extraction of benzene from the model petroleum oil, we study a more complex system: a model fish oil, which consists of a mixture of saturated and polyunsaturated fatty acids. The saturated ones were represented by palmitic (PAL) and oleic (OLE) acids, while the docosahexaenoic acid (DHA) was used as a model of typical omega−3 fatty acids present in salmon oil. DHA has many recognized human health benefits for cardiovascular diseases, diabetes, cancer, depression and various mental illnesses, age-related cognitive decline, periodontal disease, and rheumatoid arthritis 67,68 . Unfortunately, DHA cannot be produced by the human body because of an absence of desaturase enzymes 69 and consequently, DHA needs to be ingested from the food. Therefore, the extraction of DHA from different sources has won interest.
Extraction of polyunsaturated compounds with the aid of imidazolium-based ionic liquids has been previously studied experimentally 51,70,71 . These studies tested the extraction of omega-3 fatty acids (PUFAs) or ester derivatives from fish oil in the presence of silver salts 51,70,71 . However, silver salts are toxic and corrosive making its use in food application not suitable 51 . As an alternative for the use of silver salts, the addition of solvents has been proposed or the use of sorbents 71 . Some of the experimental findings from these studies are i) the extraction of omega-3 compounds is governed by the π-π stacking between the aromatic ring of the imidazolium cation and the double bonds of the omega-3 fatty acids, ii) the increase in the size of the alkyl length with its concomitant increase in hydrophobicity of the compound aids to the extraction of the PUFAs, and iii) the use of alkanes as co-solvents is recommended to enhance the extraction because ILs are highly miscible in polar solvents but immiscible in non-polar solvents.
To investigate to what extent, we can reproduce these experimental findings, we simulated the extraction efficiency of DHA from fish oil by the imidazolium-based ILs, both in the presence of a non-polar co-solvent (octane) and without it. It is observed in Figures 5 and S14 that when the length of the alkyl tail grows, also the miscibility with the OLE and PAL increases. The final configurations of the systems C2 and C4 are shown in Figure 5. Cations C8 and C12 are mixed with the fish oil, losing the biphasic nature of the system ( Figure S14). In case of C4, mixing is also observed ( Figure 5C) but only in the absence of the co-solvent octane. Apparently, octane can stabilize the fish oil with respect to the IL ( Figure 5D), probably by making it more hydrophobic. To quantify the level of DHA extraction as well as selectivity, we computed the transfer free energy of the fish oil components between the two phases in the cases of C2 and C4 (see Methods for more details). The results, presented in Table 2, show that the free energy of transfer of DHA from fish oil to IL is favorable in all cases. The most efficient extraction (i.e. lowest free energy) for mixtures without octane is observed for C2. Upon addition of octane to the mixture, the transfer of DHA becomes even more favorable, with the largest effect observed for C4. For the other components of our fish oil model, OLE and PAL, transfer to the IL is unfavorable. The differences between OLE and PAL are small in mixtures without octane (around 4 kJ/mol), which is to be expected given their comparable chemical nature. With octane added, transfer of these compounds to the IL becomes even less favorable. As a consequence, the addition of the cosolvent decreases the miscibility of the short tail ILs with fish oil, implying that less OLE and PAL can migrate to the IL phase. The consequence of this behavior on the extraction selectivity is shown in Table 2. The selectivity is computed as the ratio between the distribution coefficients of DHA and OLE+PAL in the IL (see Methods). The selectivity for the extraction has its maximum for C2. We can understand the trend in selectivity considering the interactions between the components of the mixture. In experimental studies 51,70 , it was found that π interactions between the double bonds of DHA and the imidazole aromatic ring are one of the leading forces for the extraction of this compound. These interactions are parameterized explicitly in the current Martini model. An opposing force is also present, as the nonpolar character of DHA gives rise to a repulsive interaction with the polar IL. Progressive addition of methyl groups to the IL could reduce that repulsive interaction and increase extraction of DHA, but this structural change on the imidazolium cation also results in more mixing with the other non-polar components of fish oil. Then, the selectivity of extraction is modulated by a compromise between the π interactions of DHA and the IL and hydrophobic interactions between the IL and the other molecules in the mixture. In the presence of octane, the balance is shifted and C4 becomes the IL with the highest extraction capacity (i.e. lowest free energy, see Table 2), with also an increase of selectivity of 2 order of magnitude Taken together, our data show that we can capture the experimental trends on the extraction of PUFAs from fish oil with our IL Martini model. In particular, we show the migration of the PUFA DHA from the oil to the IL phase, driven by favorable π interactions between the unsaturated DHA bonds and the imidazolium rings. We also showed a favorable effect on the extraction efficiency and selectivity by addition of co-solvents. Consequently, the proper selection of co-solvents and chain length of the cation is a critical point for the extraction of a desired compound. However, the overall efficiency is also governed by the stability of the biphasic system. As the hydrophobicity of the IL increases (i.e., longer tailed cations), or the hydrophobicity of the oil phase is reduced (in the absence of co-solvent), the biphasic system becomes less stable.

Conclusions and further considerations
We presented a new Martini CG model for imidazolium-based ionic liquids. We showed that the model can describe the macroscopic density of the ILs studied in good agreement with the experimental values. The model can also accurately describe their nanostructure, comparable to atomistic models. Additionally, our model is able to capture collective behavior of the IL, like the temperature induced phase transition of C12, a process that cannot be easily simulated using fullatomistic simulations.
A significant advantage of our new model in comparison with other CG models is that it can be used in combination with other Martini 3 CG molecules, which allows simulations of extraction experiments. To exemplify such features, two extraction studies were performed: benzene extraction from a petroleum oil model and omega-3 extraction from a fish oil model. In both cases, we find that the length of the alkyl tail of the IL was fundamental for the extraction and the selectivity. Concerning the fish oil, we also tested the effect of co-solvents in our extraction predictions. We show that the addition of octane can be a useful resource to increase the IL extraction capacity, selectivity for omega-3 fatty acids extraction, and the stability of the phase separation in the IL/oil biphasic mixtures.
Our results largely agree with previous experiments reported, indicating that the new Martini 3 CG model could be useful for the computational design of extraction experiments using ionic liquids. Addition of cosolvents, polarity modulation of tails attached to the charged groups or mixture of different IL cations could be straightforwardly employed, enabling the virtual screening of thousands of IL combinations. In the future, further development of the model should allow its usage in aqueous biphasic extraction of other biomolecules of interest like amino acids, peptides or proteins.

Acknowledgements
L.I. Vazquez-Salazar acknowledges financial support by the EACEA in the form of an Erasmus+ scholarship (Project 159680-EPP-1-2015-ES-EPPKA1-EPQR). All the authors thank the Center for Information Technology of the University of Groningen for providing access to the Peregrine highperformance computing cluster.
Groningen Biomolecular Sciences and Biotechnology Institute and Zernike Institute for Advanced Materials, University of Groningen, Nijenborgh 7, 9747 AG Groningen, The Netherlands Content: Tables   Table S1. Structural parameters of the CG models for the ionic liquids. Table S2. Size of simulation box and number of molecules for the systems studied. Table S3. IL density as a function of temperature using PME. Table S4. IL density as a function of temperature using Reaction Field. Figure S1. Density Vs Temperature for PME. Figure S2. Number of carbons vs Density for PME.         Figure S12. Partial density profile of the system petroleum oil-C12. Figure S13. Snapshots of the biphasic system petroleum oil-IL. Figure S14. Snapshots of the biphasic system fish oil-IL. Figure S15. Partial density profile of the system fish oil-C2 without octane. Figure S16. Partial density profile of the system fish oil-C2 with octane.

Supplementary Methods
3.1 Calculation of partial charges with quantum mechanics models.
3.2 Details of full atom simulation of IL.

Supplementary Tables
Ring constraints (nm)   Table S2. Size of simulation box and number of molecules for the systems simulated on this work. In the case of the fish oil between brackets are the values of the system for the case when octane is added to the simulation box. Note: Oleic Acid (OLE), Palmitic Acid (PAL), Docosahexaenoic Acid (DHA).  Figure S1. Profile of Density Vs Temperature for PME Figure S2. Profile of number of carbons vs Density for PME         . Partial density profile of the system petroleum oil-IL for the cation C12. Figure S13. Snapshots of the final state of the simulation boxes for the system IL-octanebenzene. Colour code: In green, the cation composed for the imidazolium ring and the alkyl tail. In pink, the beads correspond to the anion. In violet, the octane molecules and in red, the molecules of benzene. Figure S14. Snapshots of the final simulation box of the ILs and the fish oil without (A) and with (B) octane. Colour code: In green, the cation composed for the imidazolium ring and the alkyl tail. In pink, the beads correspond to the anion. For the fatty acids, DHA is shown in blue, OLE and PAL in yellow and orange, respectively. Octane is depicted in cyan. Figure S15. Partial density profile of the system fish oil-IL without octane for the cation C2 Figure S16. Partial density profile of the system fish oil-IL with octane for the cation C2