Predicting Melting Points of Biofriendly Choline-Based Ionic Liquids with Molecular Dynamics

Featured Application: The developed fully simulation-based melting point prediction method facilitates the design of novel green ionic liquids. Abstract: In this work, we introduce a simulation-based method for predicting the melting point of ionic liquids without prior knowledge of their crystal structure. We run molecular dynamics simulations of biofriendly, choline cation-based ionic liquids and apply the method to predict their melting point. The root-mean-square error of the predicted values is below 24 K. We advocate that such precision is sufﬁcient for designing ionic liquids with relatively low melting points. The workﬂow for simulations is available for everyone and can be adopted for any species from the wide chemical space of ionic liquids.


Introduction
The past decades have seen a steep increase of interest in ionic liquids (IL) such as solvents and electrolytes. ILs are commonly defined as salts that are liquid under 100 • C and are entirely composed of ions [1]. Due to their ionic character, they exhibit several desirable properties, for example, negligible vapor pressure, tunable polarity and high charge density. Moreover, there is a large number of different ions that can potentially yield an IL, which allows the fine-tuning of its properties for various applications [1]. Due to high thermal and (electro)chemical stability [2][3][4] some ILs are especially suitable candidates for energy applications, where higher cell voltage directly translates into increased energy density [5]. Furthermore, ILs are widely utilized as catalysts, lubricants, solvents in organic synthesis and countless other applications [6]. However, common ILs containing imidazolium cations and fluorinated anions are far from being biofriendly, which prevents their use in medicine and biotechnology [7]. For example, there is an increasing interest for usage of electroactive polymer actuators containing IL electrolytes in various medical devices [8][9][10].
On the one hand, it is true that the stability and low volatility of ILs decrease the risk of pollution by evaporation. On the other hand, even ILs of the latest generation are still toxic and often non-biodegradable. For example, fluorinated and cyano-based anions can emit hydrofluoric acid and hydrogen cyanide, respectively [7]. Toxicity depends on both cation and anion and thus less toxic ions have been considered. In case of anion, increase in the hydrophilicity as well as exclusion of fluoride-and cyano-substituents generally results in practically harmless and more biodegradable ILs [7]. In case of cation, toxicity is strongly influenced by alkyl chain length-shorter alkyl chains increase hydrophilicity and result in more benign ILs [7]. Petkovic et al. showed that choline cation is considerably less toxic than other traditional quaternary nitrogen based cations [11]. The primary drawback of introducing a more biofriendly ion to an IL composition is the resulting higher melting point; in many cases, the IL turns solid at ambient temperatures [11].
The experimental data on melting points of the biofriendly choline ILs is scarce and sometimes inconsistent due to variable water content and the presence of anion conjugated acid in the salt [12,13]. For example, for choline acetate [11,14,15] and choline citrate [16,17] significantly different melting points have been reported. Herewith, any experimental high-throughput characterization of promising ILs is unacceptably time-consuming-evaporation of water in vacuo, extraction or crystallization of ILs require much time and effort. A computational approach to the assessment of ILs properties looks most appealing. Previously, the melting points of ILs were computationally predicted with molecular dynamics (MD) simulations and quantitative structure-activity relationship (QSAR) models [18]. An overview of these methods is given in the discussion section. Their comparison reveals that the most accurate is the so-called pseudo supercritical path (PSCP) method, which has a significant limitation-requirement for the knowledge of the IL crystal structure [19,20].
In this study, we have developed simulation-based melting point prediction method for ILs with unknown crystal structure. Using this method, we have computationally investigated fourteen biofriendly, choline cation-based ILs and choline bis[(trifluoromethyl)sulfonyl]imide (TFSI). We have validated the predicted results against experimental melting points, which were available in the literature and/or measured in the current work.

Models
The choice of anions was made based on low toxicity and biodegradability of the corresponding ILs. We paired choline cation with fourteen biofriendly carboxylic acid, benzoic acid and sulfonimide-derived anions displayed in Figure 1. The harmful TFSI anion was included for comparison purpose and because choline TFSI is relatively well studied in the literature. To the best of our knowledge, the crystal structure is known only for choline acetate, choline acesulfamate, choline saccharinate and choline TFSI among the studied ILs [21][22][23][24]. Thus, for each IL, we generated a solid phase using Coulombic potential wells inspired by the PSCP method [19,20]. The solid phases were then heated and their melting points were determined by the change in the self-diffusion coefficient of the ions. The developed method is summarized in Figure 2 and explained in detail in the following sections. An example workflow is also available online in Reference [25].

Simulations Details
The MD simulations were run using the GROMACS code (GROMACS 5.1.4, Uppsala University, Uppsala, Sweden, 2016) with a time step 0.5 fs. The cut-off length of 1.3 nm for Lennard-Jones and electrostatic interactions was used. The OPLS-AA (All Atom Optimized Potentials for Liquid Simulations) force field with refitted partial charges was employed as discussed in Section 2.6. The force field parameters along with the refitted partial charges are available online in Reference [25]. The number of IL ions in systems was varied so that the packed solid phase lattice, which is discussed in the next section, would extend more than 2.8 nm in all dimensions to accommodate at least 128 or at most 350 ion pairs. Two different initial densities were tested for each IL, which were assumed from the IL or corresponding anion acid densities. Sulfonimides were packed at 1.4 kg/dm 3 and 1.6 kg/dm 3 ; carboxylic and benzoic acids were packed at 1.0 kg/dm 3 , 1.2 kg/dm 3 or 1.4 kg/dm 3 . This density determined only the initial package density as in the later simulations an NPT ensemble was used.

Solid Phase Structures
The starting structures were generated randomly with the Packmol (Packmol 1.1.2.272, University of Campinas, São Paulo, Brazil, 2009) algorithm [26]. First, these structures were equilibriated in a short MD run with the steep integrator. Subsequently, the solid phases were obtained in 200 ps NVT simulations, where the temperature was changed from 1000 K to 1 K in 100 ps and then kept at 1 K for 100 ps.
200 ps, or even 200 µs, is too short timescale for the formation of a realistic solid phase-a process that is known to be problematic and time-consuming for the ILs in experiments. The necessary macroscopic timescales are practically unreachable in simulations. Thus, Coulombic potentials wells were added to the system to facilitate the generation of more ordered solid phase structures. Similar approach was used in the PSCP method by Eike et al. to obtain a single known solid phase during a cycle of simulations for the purposes of melting point prediction [19,20].
In contrast to the PSCP method, in this work, as the placement of ions in the real solid phase was unknown for our ILs, eight different potential well lattices were constructed for each IL. Different lattices were generated by the combination of CsCl and NaCl ionic lattices with 3:2:2, 3:3:2, 3:3:3 and 4:3:2 lattice vector ratios. Thus, all subsequent simulations were conducted with eight differently packed systems for each IL. Even though cubic or orthorombic initial lattices were used, the ions were not explicitly restricted to be at the centers of potential wells even if it was the most energetically favored. Also ions were allowed to relax prior to annealing simulations, so that the IL lattice could adapt any of the Bravais lattices.
A potential well is a virtual site which interacts with the rest of the system only through a linear Coulombic potential. The linear potential scaling from distance is important to avoid the forces to approach infinity at close distances since they have no Lennard-Jones repulsion. The presence of the potential wells enable the use of elevated temperatures without the decomposition of an ordered lattice. The elevated temperatures are beneficial for allowing ions to overcome rotational energy barriers in short amounts of time, after which the system is cooled down.
After the initial structuring, six subsequent 500 ps-long NVT simulations were run again with potential wells. In simulations one, three and five only the anions were allowed to move while the positions of the cations were restrained. In simulations two, four and six the cations moved while the anions were restrained. The temperature control was similar to the initial potential well simulation: the temperature was decreased from 1000 K to 1 K within 400 ps and was kept there for the final 100 ps. However, weaker potential wells were used in each consecutive simulation so that the yielded solid phase structure could improve iteratively as the perturbation from the potential wells was weaker. An iterative approach and using position restrictions is necessary for obtaining a more organized solid phase, because the position of every ion in the system depends on all other ions.

Annealing Simulations
The potential wells were removed from the system and the obtained solid phases were equilibriated at 175 K for 2.5 ns with an NPT ensemble. Then the temperature was increased to 475 K over 30 ns (10 K · ns −1 ). The NPT simulations were run at 1 bar with the compressibility of 4.5 × 10 −5 bar −1 . The isotropic Berendsen pressure coupling algorithm and the velocity rescaling temperature coupling algorithm were used. The cut-off for van der Waals and electrostatic interactions was 1.3 nm. The particle-mesh Ewald model was used to account for the long-range electrostatic interactions with 0.10 nm Fourier spacing. The resulting trajectories were used to predict the melting points of the systems.

Melting Point Prediction
Throughout the annealing simulations, the mean square displacement of particles was calculated. From the mean square displacement the average diffusion coefficients were estimated for every consecutive 10 ps time-frame. This approach does not provide accurate diffusion coefficients, as they are usually simulated over longer time-frames at constant temperature. However, the approach allows to monitor the relative change of the diffusion coefficient during annealing. At first approximation, the diffusion coefficient, D, follows the Arrhenius' equation: where E A denotes the diffusion activation energy, R is universal gas constant and T is temperature. Assuming that for a given phase the diffusion mechanism does not change with the temperature, we take that the diffusion coefficient depends linearly on the inverse temperature. As the liquid and solid phases have markedly different pre-exponent factors (D 0 ) and activation energies, melting can be identified during an annealing simulation by a significant change in slope of the D-T dependence.
Melting points were identified for all systems as the point between the solid and liquid phase D-T dependencies. The solid and liquid phase D-T dependencies were linearized according to a few criteria: (1) the regions were distributed so that the sum of the correlation coefficients of both phases (R 2 ) was maximum; (2) systems with no distinct change in slope were excluded as there may have been only one phase throughout the simulation; (3) systems with moderate and low correlation of the solid and liquid phase fits (R 2 solid < 0.7 and R 2 liquid < 0.9) were excluded from the analysis. Subsequently, the melting point of the highest melting solid phase of each IL was reported as the predicted melting point of the corresponding IL.

Force Fields and Partial Charges
The OPLS-AA force field was employed. The force field parameters for choline and TFSI ions were obtained from References [27,28]. The rest of the force fields were developed using the generic OPLS-AA atom types. The Ligandbook repository was used as the starting point for most anions [29]. The force fields that were used in the final method are available in Reference [25].
Ion charges in MD are often taken as ±1e. As charge distribution is expected to be a dominant factor for ILs, where Coulombic forces dominate, such approximation would introduce a large error. The first estimation of the atomic partial charges was made with single point DFT (density functional theory) calculations of single ions in vacuum. CHELPG (charges from electrostatic potentials using a grid based method) charges were chosen as they are developed to accurately reproduce the electrostatic potential [30]. The calculations were done on PBE/Def2-SVP level [31,32] with the Orca (Orca 4.0.0.2, Max Planck Institute for Chemical Energy Conversion, Mülheim an der Ruhr, Germany, 2017) quantum chemical program suite [33]. As the GHELPG analysis does not depend much on the basis set, one of the simplest option (Def2-SVP) was chosen. For the same reason PBE density functional, showing a moderate performance on ILs [34,35], was used.
The CHELPG charges were utilized in MD simulations to construct solid phases using potential wells as described in the previous sections. Potential wells were removed and after a 2.5 ns period of relaxation, snapshots of the systems were taken. These snapshots served as starting points for semi-empirical GFN2-xTB (GFN2-xTB 5.9, Universität Bonn, Bonn, Germany, 2018) calculations [36,37]. Instead of CHELPG charges, Hirshfeld partial charges were extracted from the single point GFN2-xTB calculations as the CHELPG method is not available in the GFN2-xTB code. New atomic partial charges were calculated by averaging the charges of the corresponding atoms within all differently packed systems of the same IL. The ions that were further than 15 Å from the geometrical centre of the box were excluded, because the GFN2-xTB calculations were run without periodic boundary conditions and thus such ions might have been affected by the boundary effects.
Finally, the obtained averaged charges were inserted into the force fields and the workflow was restarted using the new charges. This approach yields more realistic atomic charge distribution for condensed phases than gas-phase DFT calculations of single ions.

Experimental
The ILs experimentally studied in the current work were synthesized by neutralization of choline bicarbonate solution with the respective acids [15]. Structures of the obtained ILs were confirmed by 1 H nuclear magnetic resonance (NMR) spectroscopy (Bruker Avance II  The experimental measurements of normal melting points (p = 1 bar) were done on a NETZSCH DSC 204 F1 Phoenix ® (NETZSCH Group, Selb, Germany) differential scanning calorimeter (DSC). Samples were placed in aluminum pans and heated under nitrogen with a temperature screening rate of 5 K · min −1 . The reported melting points correspond to the endothermic peak maximum temperatures. The uncertainty of measurement was ±1 K.

Experimental Results
Several simulated ILs were synthesized and characterized experimentally in the current work to better validate the simulation results. We prioritized ILs with missing values (e.g., choline 2-methylbutanoate, choline isovalerate) and those with conflicting values in the literature (choline acetate [11,14,15], choline citrate [16,17]). The glass transition temperatures (T G ) and melting points (T M ), obtained with differential scanning calorimetry, are compiled in Table 1. The water content in the synthesized ILs is also provided in Table 1.

Predicted Melting Points
For each IL, the solid phases were packed and the systems determining the melting points were identified to the criteria defined in the methods section. Following this, solid and liquid phase linear fits were made. For example, Figure 3 displays the D-T dependence of choline ibuprofenate solid phase during annealing along with the fitted lines-red for solid phase and green for liquid phase. Similar plots for all of the studied ILs are provided in the Supplementary Materials.  It is evident from Figure 3 that although the liquid and solid phases' fit is quite strong, the phase transition from solid to liquid introduces a slightly curved region. A more severe example of non-linearity among the simulated systems is choline saccharinate (Figure 4). This phenomenon occurs most likely due to the presence of more than one competing diffusion mechanism. In this case, the D-T dependence cannot be properly linearized using the Arrhenius' equation (Equation (1)), as different mechanisms may have different activation energies E A and pre-exponent factors D 0 .  Regardless of whether this situation arises in pure phases, it is unavoidable at phase transfer conditions where different phases and/or their transition states coexist. Given an infinitely long simulation, the melting process goes to completion at a discrete temperature provided that the system is crystalline. However, at simulation size and time scales, melting occurs over a range of temperatures. The diffusion mechanism is different and complex during the time when the system has begun to melt, but its structure is still dependent on the initial solid phase. This effect is manifested as non-linearity in the D-T dependence and it lowers the accuracy of the melting point prediction.
Despite the described non-linearity, the prediction of melting points is possible through an assumption, that melting point is the temperature at which the fitted solid phase dependence ends. The system at the end of the solid fit may be a complex transitional phase instead of a proper solid or liquid phase, but the end-point can still be mathematically determined. To increase robustness, we excluded from analysis systems with correlation coefficients below 0.7 and 0.9 for the solid and liquid phase fits, respectively. This way we avoided systems with large non-linearity and oscillations at phase transition. For choline saccharinate, a 0.7 minimum liquid phase fit was allowed, as none of its systems produced a liquid phase fit of above 0.9. Note that simulations with different partial charges or temperature rates have weaker correlations requiring lower R 2 thresholds.
Melting points were found as the intersection point of solid and liquid phase fits. The predicted melting points are presented in Table 2. The root-mean-square error of the predicted melting points is 18-24 K, depending on the choice of the experimental reference. It is seen that the experimental results can significantly vary among themselves, which is probably caused by different heating rates used in differential scanning calorimetry and more importantly the presence of impurities (incl. water), on which the data often remains unpublished. Regardless of the choice of reference value, the mean error of 2-10 K indicates a positive skew in the predicted values, for example, the prediction method tends to slightly overestimate the melting points.

Annealing Temperature Rate Effect
To be assured that the hysteresis effect does not affect our determined melting "onset" temperatures, we carried out annealing simulations with different temperature rates for those solid phases that yielded the predicted melting points. If the delayed phase transition, that we observe as non-linearity in D-T dependencies, is accompanied by superheating, then the predicted melting point should strongly depend on the simulation temperature increase regime.
As can be seen in Figure 5, the annealing temperature rate indeed influences the evaluated melting point values. For most ILs, these values seem to converge, but for Choline Acesulfamate and Choline 2-Methylbutanoate the resulting dependence was unclear. Within subgroups shown in Figure 5A-D, except for these two ILs, the order of predicted T M agrees with the experimental one. In case of superheating due to hysteresis, faster temperature rates should systematically produce higher melting point estimates which is not the case for our systems. Hence, it is possible that there is a hysteresis effect, but it does not severely affect our obtained results.
Interestingly, in most cases faster temperature rates tended to yield underestimated predictions rather overestimated predictions, while the latter would be the expected if superheating effect occurred. We can speculate that rapid heating either causes local overheating or the velocity rescaling temperature coupling algorithm misbehaves at such extreme conditions. In any case it is concluded, that the temperature rate of 10 K · ns −1 used in this work is sufficient to provide accurate results.
A linear regression line was added to Figure 5 to indicate trend in predicted melting point when the rate of melting is slowed. The most straightforward albeit expensive way to minimize the effect of hysteresis is to conduct longer simulations. Another way is to conduct simulations with different and not necessarily slower duration and extrapolate the results. We added a regression line to simulations with different time rates in Figure 5 to enable the estimation of extrapolated melting point by their intercept values at 0 K · ns −1 .
The root-mean-squared error of the extrapolated melting points prediction was 19-23 K, which is almost identical to the error of the regular predictions (18-24 K). It appears that most regression lines in Figure 5 tend to steer close to the experimental values than predictions at 10 K · ns −1 , but there is a severe outlier-choline 2-methylbutanoate. With the exclusion of choline 2-methylbutanoate the root-mean-squared error of the extrapolated predictions is 12-18K (17-22 K for regular predictions). Individual melting points predicted from the extrapolated results are found in Table S1 in the Supplementary Materials.

The Choice of Partial Charges
Throughout the development of the melting point prediction method, the partial charge distribution was established as a key factor. It is often approximated in the MD studies, that the ionic charges are ±1 e. Following this example, the CHELPG atomic charges were calculated for single ions in vacuum with DFT on PBE/Def2-SVP level. Using the obtained atomic charges, the determined melting point values were severely overestimated ( Table 3). The melting point estimations improved for some cases when all CHELPG charges were scaled by a factor of 0.8 to account for the possibility of partial charge transfer [41].
It is natural that lower charges lead to decrease in melting point, because the Coulombic interactions become weaker. However, arbitrarily reducing the partial charges is not justifiable and still led to overestimated melting point predictions. Partial charge transfer manifests itself in DFT calculations of large sets of ion pairs [30,42,43]. Yet, as previous studies have shown, larger associates are require to obtain more realistic ionic charges and partial charge distributions within ions [44,45]. Thus, we calculated partial charges in IL clusters using the GFN2-xTB method and Hirshfeld partitioning scheme. These charges, when incorporated into the force fields, yielding more accurate melting point predictions in most of the cases ( Table 3). The cluster geometries for the GFN2-xTB calculations were extracted from the workflow simulation trajectories obtained with CHELPG charges.

Discussion
The prediction of IL melting point was performed before with varying success via quantitative structure-property relationship (QSPR) or MD-based [18] approach. In QSPR methods, several simple descriptors are correlated to the melting point. Numerous QSPR [46][47][48][49][50][51][52] and other statistical models [53,54] were assembled, with perhaps the most impressive being a work by Torrecilla et al., who report a mean prediction error of only 1.30% [46]. However, the QSPR and QSAR-based models have an intrinsic shortcoming-they are reliable as far as the characterized substances are similar to those of the training set that were used to build the model. Therefore, they are not best suited for screening for novel compounds or for studying more complex systems, for example, mixtures of ILs.
There are several MD-based methods for melting point determination such as the free-energy based PSCP method [19,20,55], hysteresis method [56], voids method [57,58] and solid-liquid interface method [59]. They are in principle more reliable and universally applicable than the QSPR methods [18]. In practice, they have a severe limitation-requirement for the crystal structure data. The main motivation for the melting point prediction of ILs is to identify new low melting ILs prior to their synthesis. If the prediction method requires crystal structure data, then it makes synthesis and crystallization of the IL inevitable. At that point, a little effort put in measuring the melting point makes the post-prediction useless.
Our method does not require the knowledge of IL crystal structure. Instead, the solid phases are packed with the help of Coulombic potential wells. Differently from the PSCP method [19,20,55], where the wells are placed in known lattice vertices, we simulate different lattices and choose the highest melting point.
Our analysis of D-T dependencies can be considered direct as opposed to constructing phase transition free energy cycles as is done in PSCP. The authors of the PSCP method, Zhang and Maginn claimed that direct prediction of melting points of ILs is not accurate due to the superheating phenomenon which leads to severe overestimation of the melting points [18]. Contrary to that claim, in our results (Table 2), there is a relatively small skew (2-10 K mean average error) towards higher predicted values. In addition, we did not observe clear superheating hysteresis when changing the temperature rate of annealing simulations ( Figure 5). Interestingly, faster temperature rates did yield inaccurate results, but they tended to underestimate rather than overestimate the melting points. This might be happen due to local overheating or problems with the temperature coupling algorithm at more rapid temperature rates. Conducting annealing simulations with different temperature rates also yields a slightly more accurate prediction for most cases. This was shown with the regression lines in Figure 5 via extrapolation towards temperature rate of 0 K · ns −1 , which would correspond to an infinitely slow annealing.
Another factor, which did lead to markedly overestimated results for our method, was using inaccurate partial charges. In the work by Zhang and Maginn [18], the IL force field charges used to test different methods were obtained similarly to our initial DFT CHELPG charges, where the ionic charges are forced to be ±1 e. In current work, we observed large errors when using these partial charges ( Table 3). The situation improved slightly with scaling the charges uniformly down with the factor of 0.8-a typical value obtained in DFT calculations [43,60]. Arbitrary down-scaling is better than using charges of unity, but an accurate charge distribution is also required as demonstrated in the previous works [61][62][63]. The error of our predicted melting points greatly decreased when the charges were extracted from calculations with clusters instead of single ions. In this work, we have used GFN-xTB2 semi-empirical calculations to estimate charges in condensed phases, but perhaps even more accurate charges can be extracted from ab initio MD if required (see e.g., Reference [55]).
Another viable alternative is to use polarizable force fields to circumvent the need for such static partial charges altogether. Previously, this was shown to have a significant effect on interactions between ions in aqueous solutions [64][65][66][67][68] and the dynamics of the ILs [69][70][71][72]. Running simulations with polarizable force field was not feasible in current work, because polarizable force fields were not readily available for all studied ILs. Nevertheless, recent progress in the development of transferable polarizable force fields reassures that in the nearest future the presented method for predicting melting points could be applied to a larger set of ionic liquids without the need to refine charges.

Conclusions
This work introduces an MD-based method for estimating the melting points of ILs. While the method is suitable for any ionic compound, we have focused on the melting point of biofriendly choline ILs. Application of these liquids requires them to be liquid at ambient temperature, yet most of them tend to be solid. We have tested the developed method on 15 ILs comprising of choline cation and various biofriendly anions.
Compared to previously published methods, the accuracy of our method is excellent with a root-mean-square error of only 18-24 K. Moreover, compared to other approaches, our method has a relevant advantage-it does not require the knowledge of the crystal structure nor does it depend on the choice of a training set as QSPR-based methods do. Hence, our method is useful for the discovery of ionic compounds for example, for computational screening and characterization of novel biofriendly ILs prior to their synthesis. Instead of replicating an experimental crystal, the solid phases in our method are constructed with the help of Coulombic potential wells. The potential well positions are systematically varied for each IL to create optimal structures which are then melted in annealing simulations. Melting point is then determined via analysis of the D − T dependence of the annealing simulation.
Direct melting point determination methods such as ours have previously received criticism that they suffer from overheating due to hysteresis resulting in overestimated melting points. On varying the temperature rate of melting simulations, we observed no clear overheating tendency. However, we did observe significant overestimation of predicted melting points when inaccurate charges partial charge distribution was used, for example, when charges were extracted from the DFT calculations of single ions. We obtained significantly better melting point predictions only after using partial charges extracted from semi-empirical calculations of condensed phase. This indicates that both the effective ionic charges and the charge distributions within ions are detrimental for accurate simulation of ILs.
The publicly available automated workflow in Reference [25], which implements the developed method, can be readily extended to other ILs and their mixtures. As similar calculations are readily parallelizable, the only practical limitation is the availability of the force fields. We believe that that disadvantage is minor, as it is being eliminated with the development of novel transferable force field parameters.

Funding:
The EU supported this research through the European Regional Development Fund (TK141 "Advanced materials and high-technology devices for energy recuperation systems"), Estonian Institutional Research Grant IUT20-13, Estonian Personal Research Project PUT1107, and the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 793377 (BIOACT). Computational results were obtained in part using the High Performance Computing Center of the University of Tartu and in part using the EPSRC funded ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk, EPSRC grant no. EP/K000586/1).