5 M . Dyna Mix Studies of Solvation , Solubility and Permeability

During the last four decades Molecular Dynamics (MD) simulations have developed to a powerful discipline and finally very close to the early vision from early 80’s that it would mature and become a computer laboratory to study molecular systems in conditions similar to that valid in experimental studies using instruments giving information about molecular structure, interactions and dynamics in condensed phases and at interfaces between different phases. Today MD simulations are more or less routinely used by many scientists originally educated and trained towards experimental work which later have found simulations (along with Quantum Chemistry and other Computational methods) as a powerful complement to their experimental studies to obtain molecular insight and thereby interpretation of their results.


Introduction
During the last four decades Molecular Dynamics (MD) simulations have developed to a powerful discipline and finally very close to the early vision from early 80's that it would mature and become a computer laboratory to study molecular systems in conditions similar to that valid in experimental studies using instruments giving information about molecular structure, interactions and dynamics in condensed phases and at interfaces between different phases. Today MD simulations are more or less routinely used by many scientists originally educated and trained towards experimental work which later have found simulations (along with Quantum Chemistry and other Computational methods) as a powerful complement to their experimental studies to obtain molecular insight and thereby interpretation of their results.
In this chapter we wish to introduce a powerful methodology to obtain detailed and accurate information about solvation and solubility of different categories of solute molecules and ions in water (and other solvents and phases including mixed solvents) and also about permeability and transport of solutes in different non-aqueous phases. Among the most challenging problems today are still the computations of the free energy and many to it related problems. The methodology used in our studies for free energy calculations is our Expanded Ensemble scheme recently implemented in a general MD simulation package called "M.DynaMix".
Although several very powerful MD simulation programs exist today we wish to highlight some of the more unique features of M.DynaMix which is both general purpose package but also designed towards applications which can be directly connected to spectroscopic experiments with solvation state NMR as example. As MD is the only method to obtain dynamical information in realistic (but still limited) time scales we wish to give examples of motional characteristics of molecules in both bulk solvents and in complex biomolecular systems (containing for example DNA, proteins and lipid membranes) in specific noncovalent interactions of varying strength. In addition we will demonstrate the capabilities of Even if .mmol files can be created with a simple text editor using available data on molecule structure and force field parameters, their preparation can be significantly simplified by the makemol utility. This utility creates a .mmol file from two sources: a "simple molecular structure" (.smol) file and "force field" (.ff) file containing force field parameters. A .smol file contains list of atoms with their initial coordinates, partial charges and force field atom types, and list of bonds. An example of .smol file for a methanol molecule is given in Figure 2. A "force field" file contains force field parameters for a number of atom types (which can be common for different molecules). The makemol utility, using atom types specified in the .smol file, substitutes corresponding parameters from the force field file and creates the resulting .mmol file. It also generates automatically the list of covalent and torsion angles from the list of bonds (present in the .smol file) and substitutes parameters for them. The distribution of M.DynaMix contains examples of Amber94 and CHARMM27 force fields. It also contains utility char2mdx to transform general Charmm parameter file to M.DynaMix force field file. There exists a number of open source molecular editors (for example, kalziumhttp://edu.kde.org/kalzium -included into Fedora Linux distribution, or Dundee ProDrg server -http://davapc1.bioch.dundee.ac.uk/prodrg/), that allow to draw a molecule with a mouse and generate a preliminary optimized molecular structure, as well as generate the list of bonds. Next, one needs to assign the atom types for each atom (the notations for atom types are different in different force fields), and specify partial atomic charges, which are usually not included into general force fields, but need to be calculated by a quantum-chemical program for each specific molecule. After the .smol file is created, one can use Makemol utility to create the structure-parameter .mmol file for the given type of molecules.

Simulations
Molecular dynamics simulation parameters are described in the main input file. From version 5.0, the format of the main input file was changed significantly compared with the previous publication of the program (Lyubartsev & Laaksonen, 2000). In the new format, all simulation parameters and features are specified in the form "Keyword value(s)", where keyword describes one specific parameter or property of the simulation (for example, time step, temperature, restart control, etc). For missing keywords, default values are used. A complete description of available keywords and their parameters is given in the manual at the M.DynaMix distribution site (http://www.mmk.su.se/mdynamix).
Finally, in order to start simulation, one needs to specify the start configuration of the simulated system. It is always possible to start from an already prepared (by various means) configuration written in a file in the xmol (xyz) format. This option may be preferable if one wish to simulate a specifically organized system, e.g. lipid bilayer or protein in a folded form. It is also possible to generate the start structure automatically. In this case, molecules are placed with random orientation and their centers of mass on FCC lattice. Clearly, in such case one can expect numerous overlappings of atoms belonging to different molecules. In order to deal with the problem, the following procedure can be used: 1. One starts from a low enough density, and run a short run (of order of 1000 MD steps) with the options "Cut large forces", "Velocity scaling" for temperature control and with pressure barostate switched off. The "Cut large forces" option limits forces acting on atoms by a specified cutoff value, while keeping the direction of forces. As a result, bad atom contacts and overlappings disappear quickly. 2. One proceeds for another few thousands steps specifying expected density of the system and using option "Change volume slowly". The system gradually change volume and increase the density to expected, providing in the end preliminary equilibrated system 3. After the second stage it is possible to switch from "velocity scaling" control of temperature to "Nosé thermostat". It is advisable to continue for some time to run simulation at constant volume, and when the pressure stabilizes, turn on the barostat if constant-pressure simulations are desirable.

Post-simulation analysis
The main result of a simulation is a set of trajectory files which contain configurations of the system (Cartesian coordinates and optionally velocities of all the atoms in the system) as a function of the simulation time. The trajectory analysis suite Tranal includes a number of utilities to compute a great variety of structural, thermodynamical and dynamical properties of the system, among which are radial and spatial distribution functions, various time correlation functions, diffusion, residence times, dielectric permittivity, order parameters and lateral diffusion for bilayer-like systems, etc. The suite includes the base module which can read trajectories of various formats (including those of other simulation packages), and specialized modules to compute specific properties.

Expanded ensemble mode
This special feature allows computations of solvation free energies by the expanded ensemble (EE) methodology (Lyubartsev et al., 1992(Lyubartsev et al., , 1998 EE method was implemented together with Wang-Landau algorithm for automatic optimization of weighting factors in expanded ensemble, which allows obtaining accurate free energies in a single MD run. A short description of the algorithm is given below. The EE method implies a gradual insertion/deletion of the studied solute molecule into/from the solvent. The insertion parameter is introduced which describes the degree of insertion of the solute, so that = 1 describes the solute molecule fully interacting with the solvent while = 0 represents the solute completely decoupled from the solvent (that is, case = 0 describes pure solvent and the solute molecule in a gas phase of the same volume). The insertion parameter can accept a number of (fixed) values { i }, i = 0,...,M in the range between 0 and 1. During the simulation, which can be carried out by either Monte Carlo (Lyubartsev et al., 1992) or molecular dynamics algorithm (Lyubartsev et al.,1994), attempts are made to change the insertion parameter to another (normally, neighboring) value, with acceptance probability: where V Ss ( ) is the interaction energy of the solute particle with the solvent corresponding to the given insertion parameter and η i are so-called balancing (weighting) factors introduced with the purpose to make distribution over subensembles close enough to the uniform one (Lyubartsev et al., 1992). During the simulation, probabilities of different subensembles ρ i are defined, and free energy difference between states with fully inserted and fully deleted solute (excess solvation free energy) is computed by: In previous works on computation of the solvation free energies (Lyubartsev et al., 1994, a linear scaling of solute-solvent interaction with was used, namely V Ss ( ) = V Ss (1). Such a scheme has a shortage that the repulsive core of the Lennard-Jones potential decreases very slowly with decrease of , which leads to necessity to consider very small values. Now a new scheme is implemented in which interaction V Ss ( ) between the solute and solvent atoms (which is supposed to be a sum of the Lennard-Jones and electrostatic terms) is scaled according to the following: where V Ss LJ and V Ss el are Lennard-Jones and electrostatic parts of the solute-solvent interaction respectively. The rationale between 4 scaling of the Lennard-Jones interaction is that the effective core radius of the LJ potential scales then as ( 4 ) 1/12 = 1/3 , and thus the effective volume of the core (approximately proportional to the free energy of cavity formation) scales linearly with . Scaling 2 for the electrostatic interactions is chosen in order to switch off them faster than the repulsive LJ interaction, to minimize problems with hydrogen atoms of water molecules which do not have LJ potential. With scaling of the solute-solvent interaction described by (3), it is possible to choose i points uniformly distributed in the range [0:1] so that a reasonable acceptance probability of the transitions is maintained through the whole range of subensembles.

91
Another modification of the EE procedure implemented in the new version of M.DynaMix is the automatic choice of the balancing factors η i by the Wang-Landau algorithm (Wang & Landau, 2001). In the application to the expanded ensemble technique the Wang-Landau algorithm can be formulated as follows. We start simulations with zero balancing factors. After visiting a subensemble i, a small increment Δη is added to the corresponding value of the balancing factor: η i (t+dt) = η i (t)+Δη, which decreases the probability to go to already visited states and favors to attaining a uniform distribution. After a certain number of steps (a sweep), when the system visited all subensembles and passed several times (default 2) between the end points, the value of the increment Δη is decreased twice. After a certain number of sweeps (usually 10-12), the value of the increment is becoming very low, and simultaneously the profile of balancing factors is tuned in a way providing uniform walking in the space of subensembles. After that the equilibration stage ends, and production run with fixed balancing factors is made yielding the solvation free energies.

Examples of M.DynaMix applications
In this chapter we concentrate only on some selected types of MD simulations which are just very few of many M.DynaMix has been used in the past.

Solubility and permeability
Solubility (or lack of it), in particular in water, but also in other solvents and solvent mixtures is among the most important problems in the development of new drug molecules today, in particular, for poorly water-soluble, orally bioavailable drugs, with a controllable release rate. Computer modelling and simulations (together with combinatorial chemistry, high-throughput screening and genomics/proteomics) have long been key tools in modern drug design but has not been utilized much in solubility studies. After drug molecules are dissolved in the beginning of the administration they start their journey in living organisms to find the target. During this process they are adsorbed and distributed and doing so they need penetrate numerous obstacles, including the blood-brain barrier. To study the molecular details solubility and permeability during this transport, MD simulations together with free energy calculations are highly useful (Lyubartsev et al., 1998c. Solubility of drug molecules and their permeability across lipid membranes have been studied using M.DynaMix by calculating logP partition coefficients from the solvation free energies (Lyubartsev et al. 1999;Åberg et al. 2004).

Solvation in mixed solvents
Many substances turn out to be readily soluble in solvent mixtures (mixed solvents) but not in the pure components and vice versa (Reichardt, 1988). Selective solvation of ions from salts in mixed solvents is another, early discovered, property. Both IR measurements and computer simulations have confirmed the existence of micro-heterogeneities and conglomerate water structures at higher concentrations of the more organic co-solvent. Many amphiphilic and biomolecular systems are better dissolved in mixtures, as a result of preferential solvation. This Section will highlight the advantage of using MD simulation to study preferential and selective solvation to obtain insight into the complex molecular processes behind them and to give a picture and rules in the delicate balance of interactions between the different molecular species ). We will first use a small disaccharide as a solute in discussing the preferential solvation. However, before discussing the solvation in mixed solvents we highlight a number of key features in hydration and solvation of the disaccharide in several pure solvents.
An important question in studying carbohydrates is the role of the solvent, most often water, and in particular its influence on the structure. Several aspects, such as conformational states, rotation barriers and intra/intermolecular hydrogen bonds, flexibility/rigidity, internal dynamics, hydrophilicity/phobicity, as well as electrostatic interactions are important in studying and describing the full three-dimensional structure of this category of molecules. Even small differences in the molecular geometry or conformation can easily affect the hydration around these inherently flexible molecules. Radial distribution functions (RDF) as well angular distribution functions (ADF) are the most commonly used tools to study solvent structure and conformational populations in computer simulations. Carbohydrates typically have a large number of hydroxyl groups ("half water molecules") which easily engage themselves in the hydrogen bond network of water around them leading to very distinct and sometimes rather spectacular solvation structures.
We have studied a small prototypic carbohydrate molecule -D-Manp-(1→3)--D-Glcp-OMe containing both a mannose and a glucose ring connected with a glycosidic bridge, by carrying out molecular dynamics simulations in water, methanol, dimethylsulfoxide (DMSO) and also in the mixture of water and DMSO as the solvents (Vishnyakov et al., 1999(Vishnyakov et al., , 2000a(Vishnyakov et al., , 2000b. The mixed solvent allows experimentally low temperature studies at a molar ratio of 3:1 of water and DMSO making the solution more viscous. We use the same concentration in our MD simulations. We will simply call the molecule as "disaccharide" in the text below. Conformations of our disaccharide are well described using the two torsional angles across its glycosidic linkage, namely φ and ψ (Figure 3). Our simulations exhibit clearly two main minima (potential wells) on the adiabatic energy surface, one close to g-(-60°, we call it "well A") and the other one close to g+ (+60°, we call it "well B") with respect to the φ torsional angle. All simulations are carried out using the M.DynaMix package (Lyubartsev & Laaksonen, 2000). The solvation structure is presented using spatial distribution functions (SDF) (Kusalik et al., 1999) visualized with the gOpenMol software (Laaksonen 1992;Bergman et al. 1997) http://www.csc.fi/english/pages/ g0penMol/tutorials. In all SDFs below the iso-density threshold is three times the bulk density in the simulations except when stated otherwise.
In water solution the disaccharide shows a clear preference for the well A. The hydration is dominated by the first-shell structuring in the solvent. The large number of hydrophilic centers in carbohydrates imposes a strong anisotropic structuring on the surrounding solvent. The hydroxyl groups were found to be extensively hydrogen-bonded to surrounding water molecules (Figure 4). For an optimal hydration both the solute and the solvent should have a compatible topology for the hydration requirements of the involved functional groups. The disaccharide molecule is found to be fairly rigid when dissolved in water showing no rotation around the glycosidic linkage. In methanol solution (Figure 5), the solvation structure, as seen in the SDFs, reflects the ability of methanol to both donate and accept hydrogen bonds. Besides, having in principle the same possibility to form hydrogen bond combinations as water, methanols can even hydrogen bond to the non-hydroxyl oxygens. The methyl groups in methanol restrict the formation of similar solute-solvent hydrogen bonding pattern as in water. However, the hydrogen bonding can be easily seen in the SDFs; e.g. the donation of hydroxyl protons in methanol around O2m (upper left), O6m (middle right), and O4g (lower left). Note, m denotes the mannose ring and g denotes the glucose ring, while the ordering of the oxygens follows the standard rules. The solvation structure of methanol around the disaccharide is found less extended than in water, but rather more like an average between that of water and that of DMSO presented in the next paragraph.
In DMSO solution the most probable hydrogen-bonding sites for the DMSO oxygens around the disaccharide are depicted in Figure 6. Distinct features showing hydrogen bonding to the hydroxyl groups are present. Compared with water, where the corresponding SDFs show distinct belt-like shapes around the entire saccharide, DMSO behaves differently. While four possible combinations of hydrogen bonds between the hydroxyl groups and water are possible, DMSO can only be engaged in one type of H bond as it cannot donate hydrogen. However, distinct regions are observed for HO4g (middle left in Figure 4), HO2m (upper right), and HO2g (lower right). Furthermore, the bigger size and much heavier mass slow down the mobility of DMSO and contribute to a more distinct type structure.  In Water-DMSO mixture (Figure 7) we can see that the saccharide is selectively solvated. By combining the information obtained from using the traditional RDFs and the threedimensional SDFs, we can see the different parts of disaccharide preferentially solvated by either water or DMSO. We can also see that in some parts there is a clear competition between the two components to different hydroxyl groups. For example, DMSO acts as a hydrogen bond acceptor and solvates predominantly the regions close to HO2g and HO2m (see Figure). Water, acting as a hydrogen bond donor, hydrates selectively the regions around O6m and O2g. The SDFs provide a very detailed picture of solvation. Different SDF representations and techniques to analyze solvation graphically are given by Bergman and co-workers. An interesting detail is that when DMSO is acting as a hydrogen bond acceptor to the HO2g hydroxyl group, it can also simultaneously act as an acceptor a specific water molecule hydrogen bonding to O6m. This leads to a hydrogen bonded complex which involves two hydroxyl groups of the disaccharide and one water molecule and DMSO molecule.
We have studied preferential solvation of phenol in equimolar water/acetonitrile and water/ethanol mixtures (Dahlberg & Laaksonen, 2006). In this work we introduce a special "difference spatial distribution function" (ΔSDF) to illustrate the excess densities over the cosolvents in preferential solvation locations around the phenol molecule in these two solvent mixtures. This work was inspired by intermolecular 1 H NOESY experiments carried out by Bagno (Bagno, 2002) nicely showing preferential solvation around organic molecules.

www.intechopen.com
Simulations and experiments agree well with each other but reveal also some differences. In fact, this is just another example where NMR and MD simulation techniques can be successfully combined (Odelius & Laaksonen, 1999).

Mixed solvents
The Water-DMSO system as such is a commonly used mixed solvent. Other examples of such solvents are the binary mixtures of water, alcohols and acetonitrile, including both aqueous and nonaqueous mixtures of them. Mixed solvents are important for example in chemical industry. We have studied the DMSO-water binary system which is known to exhibit strongly non-ideal behavior as a high positive heat of mixing and negative excess mixing volume. It was suggested earlier that two water molecules first hydrogen-bond simultaneously to the sulfoxide oxygen with a third water molecule found hydrogen bonding and bridging the two water molecules, in this way forming a 3:1 water-DMSO complex. It should be straightforward to study if this is the case using computer simulations. For this purpose we introduced a "multi-particle SDF" in which the local coordinate system was defined based on the three oxygens; the sulfoxide oxygen and the two oxygens in water molecules involved in hydrogen bonding to the sulfoxide oxygen . This was possible because of the long enough residence times for these particular hydrogen bonds. It is easy to see from the MP-SDF (Figure 8) where the second shell water molecules to the sulfoxide group are found. It becomes clear from the simulations that no third water molecule hydrogen bonding to the first shell water molecules was found.

Hydration of Ni 2+
MD simulations of a Ni 2+ ion in water have been carried out to study the hydration structure around nickel (Egorov et al., 2006). The analysis is extended to the second hydration shell around the divalent ion. The nickel aqua-complex has been treated without any constraints in order to analyze the structure and dynamics, as well as molecular mechanisms of water substitutions and exchange in the first and second shells around the cation. The simulations show that the structure of [Ni(H 2 O) 6 ] 2+ complex is very stable (Figure 9). The main molecular mechanisms contributing to reorientational motion of the whole complex are the "pushes and kicks" given by the second hydration shell water oxygens. The rotations of water molecules themselves in the second shell are only changing the hydrogen positions and do not affect the reorientation of the complex. In spite of frequent exchange of water molecules between the bulk and the second hydration shell the overall dynamical behavior in the first and second shell is very similar.

Hydration and coordination of counterions around DNA
Nucleic acids are highly negatively charged polyelectrolytes, displaying a considerable sensitivity to their ionic surroundings while undergoing different structural transitions and interacting with other charged species in their surroundings. Since early studies of DNA structure it has been shown that a variation in the counter ion nature or concentration, and in the water content, lead to different conformational properties; these properties depend also on the base pair composition and sequence (Saenger, 1984;Leslie et al., 1980).
Computer simulations and MD in particular are a very powerful technique to obtain detailed information about processes behind hydration and counterion coordination, including the dynamical behaviour of cations around the charged surface of DNA. However the highly charged nature of DNA has lead to a later application of MD simulation compared to that of proteins. In fact, one of the main problems in the simulation of these polyelectrolytes was the treatment of the long range electrostatic interactions, which has been solved with the use of the Ewald summation methods (Allen & Tildesley, 1987;Laaksonen et al., 1989). Proper treatment of these interactions, together with the developments of Force Fields designed for explicit solvent models (Cornell et al., 1995;MacKerell et al., 1995;MacKerell, 2004) and advances in computer power have lead to large use of molecular dynamic simulations in the study of nucleic acids. Here we focus on the modeling, at different levels of sophistication, of nucleic acids in physiological environments with different species of counterion, either mono and multivalent ( (Lyubartsev & Laaksonen, 1998b;van Dam et al., 1998;Mocci et al, 2004;Bunta et al., 2007;Korolev et al., 2001Korolev et al., , 2002Korolev et al., , 2004aKorolev et al., , 2004b and we give some examples of how the MD trajectories can be analyzed and the results checked against experimental data using M.DynaMix software.

Interactions with monovalent counterions
Until almost the end of the last century it was believed that monovalent counterions did not show any "sign of localization or dehydration due to interactions with double-stranded (ds) nucleic acids" (Anderson et al., 1995) and that the interactions were not affected by the base pair sequence. Since that time, many investigations have lead to a revision of this view, and it is now believed that monovalent cations binds to DNA partially losing the solvation shell, often in a sequence specific way which also depends on the nature of the cation. As examples of this concept we discuss some of our studies focused on DNA interactions with alkali cations. In order to understand the differences in the behaviour of cations of different size, we performed MD simulations on DNA solutions containing either the physiologically relevant DNA counter ion Na + , and the smallest and the largest among the alkali ions: Li + and Cs + (Lyubartsev & Laaksonen, 1998b;van Dam et al., 1998). The simulations were performed mimicking an infinite array of in parallel ordered DNA duplexes, by placing a DNA decamer [d(ATGCAGTCAG) 2 ] along the Z direction of the simulation box and applying periodical boundary conditions. Important information on ion-DNA interactions can be obtained inspecting the radial and spatial distribution function (RDF and SDF). The RDF gives the probability of finding a counterion type at a distance r from selected DNA atoms, or from DNA axis, relative to the probability expected from a random distribution at the same density. The SDF is the three dimensional distribution of a given atom species in a local coordinate system fixed on some (portion of a) molecule. In the SDF calculated using Tranal utility of M.DynaMix the value of the function in each point corresponds to the probability of finding there an atom compared to the bulk. Representation of a four dimensional function in two or three dimensions can be done in different ways Kusalik et al. 1999); one of the most used to visualize solvation structure and ion binding modes is the iso-intensity representation, where a surface links the points where the SDF function has a chosen value. RDF calculated between the ions and selected atoms either in the minor or major groove, or in the backbone, are able to show at a glance important differences in the behaviour of the counter ions. For example in Figure 10 it is shown the RDF between counterions and P atoms.
The smallest ion has a remarkable higher probability to be directly bound to the phosphates compared to the other counterions. On the other hand the RDFs calculated between the counterions and selected atoms in the minor and major grooves, clearly indicate that the smallest ion do not penetrate in the grooves, while the largest has a remarkable preference www.intechopen.com for the minor groove. Analysis of the SDFs of the water and of the counterions around the binding site is very helpful for understanding these differences in the interactions and to understand the conformational preference of DNA in presence of different counter-ions. In the case of Lithium the SDFs of the water atoms and of the cation around the phosphates groups explained the higher preferences of Li + for this site compared to other alkaline ions.
In fact the average organization of the water molecules around the phosphates' oxygen atoms contains a hole whose dimensions is perfect for the small Li + ions to fit in, keeping a tetrahedral hydration shell, as shown by the SDF in Figure 11. The larger alkali ions cannot bind in the same position without largely perturbing this solvation structure. The observed difference in the binding modes offered also an important key to the explanation of why Li + counterions promote the transitions of DNA from the B conformation to the C-form. In absence of experimental data that could confirm in a direct way the binding preferences observed in the simulations, comparison with experimental data on the diffusion of studied ions in presence of DNA gave an important validation on the used potential model.
In the field of DNA -monovalent counterions, in the last 15 years a large attention has been devoted to the dependence of the binding mode of alkaline counterions on the sequence of DNA base pairs. Among the most studied DNA sequences are the so called A-tracts: uninterrupted sequences of at least 4 A·T base pairs without 5'TpA3' steps. The interest on these sequences is motivated by their peculiar structural features: when inserted in phase with the helix turn A-tracts induce the macroscopic bending of DNA, they are characterized by a narrow minor groove which is known to possess a very ordered hydration structure, called the hydration spine. It has been proposed that sequence specific interactions of monovalent counterions in the minor groove of A-tracts were responsible of these structural features (McFail-Isom et al., 1999). In fact an electrostatic collapse of the A-tract around ions substituting one or more water molecules of the hydration spine in the minor groove of Atract could induce either the bending, through an electrostatic collapse of the structure around the ion, and the narrowing of minor groove. Unfortunately experimental investigations on the sequence specificity of the binding to macromolecule of ions characterized by high mobility, like Na + or K + , give results whose analysis are not straightforward and leave space to multiple interpretations.
Interesting experimental results had been obtained by NMR dispersion experiments (Denisov & Halle, 2000), measuring the NMR relaxation times of Na-23 in A-T rich sequences solutions at several magnetic fields; therefore we modified the M.DynaMix code to calculate the NMR relaxation time of this ion. It is important to note that even if the observed relaxation times are in the second time scale, the motion responsible of such decay are in the subpicosecond to the nanosecond time scale, thus perfectly accessible to MD simulations of DNA oligomers in aqueous solution (Mocci et al, 2004;Odelius & Laaksonen, 1999). The coupling of these techniques was used for studying the interactions between Na + and a DNA oligomer containing either A-tracts and non A-tract regions, the sequence [d(CTTTTAAAAG) 2 ]. In agreement with previous simulations on A-T rich oligomers (Mocci & Saba, 2003) we observed direct binding with long residence time only in the minor groove of one of the A-tracts. More specifically, a sodium ion, partially losing its hydration waters, intruded in the minor groove substituting one water molecule of the hydration spine (see Figure 12), and remained therein stuck for almost the entire simulation. Fig. 12. SDFs of Na+ (yellow), water's oxygen (violet), and hydrogen (cyan) around DNA calculated during the binding to adenine N3 and thymine O2 in the minor groove (ca. 7 ns). The intensity levels are 10 for H and O of water, 1000 for Na + . Reproduced with permission from J. Phys. Chem. B 108, 16295-16302 (2004).

Copyright 2004 American Chemical Society
Combination of the MD results with the quadrupolar relaxation experimental data offered a very good way to verify whether the high occupancy of that binding site, observed in the simulation, was reliable. In fact, the combined approach revealed that the occupancy of the binding sites in the minor groove of uninterrupted adenine sequences should be actually pretty low, indicating that longer samplings would be required for proper evaluation of the occupancy of those binding site.
The MD-NMR combination allowed also obtaining information at the microscopic level on the NMR relaxation processes, showing how the polyion affects the NMR parameters of Na-23. Importantly it indicated that the speed-up of the Na-23 relaxation rate, compared to simple electrolyte solution, is mainly due to the contribution of ions directly bound to the oligomer's surface, while for a long while it has been thought that no sign of dehydration due to the interaction with DNA was exhibited by monovalent counterions (Anderson et al., 1995).
The MD studies discussed above revealed that the preferential binding modes of each cation at the DNA surface are dependent on the hydration structure of DNA and of the cations; such binding modes could have not been investigated without explicit inclusion of the solvent in the simulations. Unfortunately explicit modelling of the solvent is the most computationally demanding part of the calculations and to extend the MD simulations to spatial dimension beyond those of oligomeric DNA fragments it is necessary to use simplified models for the solvent. A method to keep a realistic picture of the DNA-cations specific interactions while simplifying the description of the solvent is based on the use of effective potentials for the interactions between solutes (Lyubartsev & Laaksonen, 1995. This method is based on reconstruction of effective, solvent-mediated interaction potentials between DNA atoms or group of atoms and the cations from MD simulation with the explicit solvent. More in detail the construction of effective potentials start from the RDF between the solute molecules calculated from a fully atomistic MD simulation with explicit solvent; by means of a reverse Monte Carlo procedure (Lyubartsev & Laaksonen, 1995) the RDF functions are used to construct a set of effective interaction potentials that, when used to run a MD simulation without the solvent should reproduce the same solute-solute RDF. This approach allows obtaining potentials able to maintain the ion specific information also when using an implicit model of the solvent.

Interactions with multivalent counterions
The difficulties in experimentally determining the favourite binding modes with DNA are not limited to monovalent counterions; in facts, also the study of interactions with multivalent molecular cations as polyamines (PA) are affected by the same problem. Polyamines, which are positively charged organic compound having two or more primary amino groups, are present in the living cells and interact with cellular polyelectrolytes like DNA. Despite the fact that PA, especially spermine4+, are largely used as DNA crystallization agents, only a very limited number of crystallographic studies report their presence in the crystallographic cell. Simulations have been extremely useful in understanding why these components are "invisible" in X-ray studies. MD simulations performed mimicking an infinite array of parallel BDNAs (Korolev et al., 2001(Korolev et al., , 2002(Korolev et al., , 2004a(Korolev et al., , 2004b in presence of either the natural polyamines spermine 4+ (H 3 N + -(CH 2 ) 3 -NH 2 + -(CH 2 ) 4 -NH 2 + -(CH 2 ) 3 -NH 3 + ), spermidine 3+ (H 3 N + -(CH 2 ) 3 -NH 2 + -(CH 2 ) 4 -NH 3 + ), putrescine 2+ (NH 3 + -(CH 2 ) 4 -NH 3 + ) or the synthetic polyamine diaminopropane 2+ (NH 3 + -(CH 2 ) 3 -NH 3 + ) revealed that these highly charged cations interact strongly with different groups on DNA. However all of the polyamines adopt disparate binding modes that make their detection pretty difficult with x-ray diffraction. The simulations also showed that PA and Na + have a different impact on DNA hydration and structure: while the Sodium cation attracts and organizes the water molecules around DNA, the polyamine pushes water away from the minor groove and induce a significant narrowing of the same. Differences in the binding preferences are observed also among the PA: while a small fraction of divalent polyamines can be found in the major groove, the other two PA are nearly absent. Furthermore differences in the binding modes where observed also between the synthetic and naturally occurring divalent PA, giving some hints on why nature selected only one of the two.

Molecular dynamics of lipid bilayers
In this section we describe application of M.DynaMix package for simulation of lipid bilayers. Lipid bilayers represent a framework of biological membranes, which are very complex heterogeneous systems consisting of many different types of lipids, sterols, proteins, carbohydrates and various membrane associated molecules which are involved in a variety of cellular processes. Biomembranes surround cells: a membrane separates the interior of a cell from the outside environment. Being selectively permeable, membranes participate in control of the movement of various compounds (substances) into and out of cells. Permeability of drug molecules through membrane is one of the key properties defining efficacy of the drug, since drug molecules have to penetrate numerous membrane barriers in order to reach their targets. Studies of drug solubility in both aqueous and lipophilic environment are thus important for understanding drug transport in living organisms.

Fig. 13. Dipalmitoylphosphatidylcholine (DMPC) lipid
Lipid molecules constituting biomembranes differ with respect to the type of hydrophilic head-group and occur with a wide variety of hydrophobic hydrocarbon chains. Usually the most abundant phospholipid in animal and plants is phosphatidylcholine which is the key building block of membrane bilayers. An example of such lipid, dipalmitoylphosphatidylcholine (DMPC) is shown in Figure 13. It is not surprising that lipid membrane bilayers composed of various phosphatidylcholine lipids have been extensively studied by molecular dynamics as soon as computer hardware allowed to do such simulations. However, obtaining of reliable information on physical-chemical properties of lipid bilayers was, in many studies of the last decades, seriously limited by two factors: 1) necessity to simulate a large (of the order of few tens thousand) number of atoms during long (hundreds of nanosecond) simulation time and 2) not satisfactory reliability of the commonly available force fields to describe behaviour of lipid bilayers consistently, with respect to variation of lipid chemical structure, composition, thermodynamic conditions (Lyubartsev & Rabinovich, 2011). While solution of the first problem is determined by advances in the development of computer hardware, the issue of proper parametrization of the force field requires extensive computational work including numerous test simulation of the bilayer systems and calculations of different experimentally measured properties.
In the last several years we have used M.DynaMix package in order to improve CHARMM force field (Feller & MacKerell, 2000) which is one of the frequently used in biomolecular simulations. While been fully atomistic, this force field have potential advantages in comparison with the Gromos force field (Berger et al., 1997), which is based on the united atom model. However, recent detailed investigations have shown that the CHARMM force field has not-negligible disagreement with experiment in description of number important properties of lipid bilayers. For example, the CHARMM force field favour to rigid gel-like structures of bilayers composed of saturated lipids, and in order to keep bilayer in a natural liquid crystalline phase one need to apply surface tension (Lyubartsev & Rabinovich, 2011: and references therein). In paper (Högberg et al., 2008) a solution was suggested how to improve the CHARMM force field in order to simulate DMPC lipid bilayer in constantpressure tensionless simulations. The correction included two changes: 1) scaling of the socalled 1-4 electrostatic interactions (of atoms separated by exact 3 covalent bonds) was introduced, which was tuned to reproduce experimentally measured ration of trans-and gauche conformations in hydrocarbon chains, and 2) partial atom charges in the lipid headgroup were recalculated from high-quality ab-initio calculations. It was demonstrated in paper (Högberg et al., 2008) that these modifications of the CHARMM force field allowed to obtain, in 100 ns constant-pressure simulations, excellent agreement with experimentally measured properties of DMPC bilayer such as area per lipid at zero tension, electron density, structure factor and order parameters. An important feature of these simulations was also that long-range correction to the Lennard-Jones potential was included into pressure calculations. This correction, having physical origin in the dispersion (van-der -Waals) forces, can have an order of 100 -200 bar in the typical range of cut-off distances, and is known to be important in correct determination of the average bilayer area, as well as affecting phase behaviour of the bilayer (Lyubartsev & Rabinovich, 2011).
In continuation of paper (Högberg et al., 2008), the modified CHARMM force field was used to simulate bilayers consisting of two similar lipids: DSPC (which differ from DMPC lipid displayed in Figure 13 that it contains 18 carbon atoms in each tail) and DOPC lipid (which exactly as DSPC contains 18 carbon atoms in each tail but have a double bond between 9-th and 10-th carbons in the second tail). Despite very similar chemical structure, bilayers composed of these lipids show different behaviour, and have noticeably different temperatures of the gel -liquid crystalline transition: 24 o C for DMPC, 53 o C for DSPC, and 5 o C for DOPC, that is why in the physiological range of temperatures DMPC and DOPC bilayers exist in a liquid crystalline phase while DSPC is in a gel phase. Simulations of the three bilayers were carried out at 30 o C, for 128 lipids in the presence of 3840 water molecules (which correspond to fully hydrated state of phosphatydylcholine lipids with 30 water molecules per lipid). In the beginning of simulations, the lipids were arranged in two leaflets, each leaflet being generated by translation of coordinates of one lipid molecule in X and Y direction. The lipid spacing was chosen to correspond to a typical value of area per lipid for liquid-crystalline phase, 64 Å 2 . The necessary amount of water molecules was added outside the lipid bilayer, and the system was put into the periodic boundary conditions. The systems were simulated 1 ns under constant volume and then 1 ns under constant pressure and isotropic cell fluctuations. The obtained configurations were considered as starting points for longer simulations with independent cell fluctuations in Z and XY directions. All the systems were simulated after that for 100 ns, with the first 20 ns considered as equilibration. Simulations showed the picture which corresponds well to the behaviour expected from experimental observations: while in DMPC and DOPC bilayers the lipids formed quickly a liquid-crystalline phase, with the area per lipid 59.5 and 62.7 Å 2 respectively, the DSPC lipids became clearly ordered in a tilted structure, with much lower area per lipid (51 Å 2 ) which is a typical value for a gel phase. Final snapshots of DOPC and DSPC bilayers are shown in Figure 14. The structures of two bilayers are strikingly different, taking in mind the fact that the only difference between the two kinds of lipids is presence of one double bond in the middle of the second tails of a DOPC lipid. Nevertheless, this behaviour is what one can expect from experimental observations.
As our exploratory simulations show that we can get reliable, consistent with experiment, behaviour of different bilayer systems, we can use the developed methodologies to address to more challenging problems, related to permeability of different substances across membranes, simulations of membrane proteins, ion channel, effects of other membrane associated molecules (cholesterol, polypeptides, anaesthetics) on membrane properties , and relating observation of these studies with the features essential for biological functioning, thus implementing the idea of a "computer laboratory" for biomembrane research.