A criterion to identify the equilibration time in lipid bilayer simulations

With the aim of establishing a criterion for identifying when a lipid bilayer has reached steady state using the molecular dynamics simulation technique, lipid bilayers of different composition in their liquid crystalline phase were simulated in aqueous solution in presence of CaCl$_2$ as electrolyte, at different concentration levels. In this regard, we used two different lipid bilayer systems: one composed by 288 DPPC (DiPalmitoylPhosphatidylCholine) and another constituted by 288 DPPS (DiPalmitoylPhosphatidylSerine). In this sense, for both type of lipid bilayers, we have studied the temporal evolution of some lipids properties, such as the surface area per lipid, the deuterium order parameter, the lipid hydration and the lipid-calcium coordination. From their analysis, it became evident how each property has a different time to achieve equilibrium. The following order was found, from faster property to slower property: coordination of ions $\approx$ deuterium order parameter $>$ area per lipid $\approx$ hydration. Consequently, when the hydration of lipids or the mean area per lipid are stable, we can ensure that the lipid membrane has reached the steady state.


I. Introduction
Over the last few decades, different computational techniques have emerged in different fields of science, some of them being extensively implemented and used by a great number of scientists around the globe. Among others, the Molecular Dynamics (MD) simulation is a very popular computational technique, which is widely used to obtain insight with atomic detail of steady and dynamic properties in the fields of biology, physics and chemistry. In this regard, a critical aspect that must be identified in all the MD simulations is related to the required equilibration time to achieve a steady state.
This point is crucial in order to avoid simulation artifacts that could lead to wrong conclusions. Currently, with the increment of the computing power accessible to different investigation groups, much longer simulation trajectories are being carried out to obtain reliable information about the systems, with the purpose of approaching the time scale of the experimental phenomena. However, even when this fact is objectively desirable without further objections, nowadays, much longer equilibration times are arbitrarily being required by certain reviewers during the revision process. From our viewpoint, this should be thoroughly revised due to the following two main reasons: first, because it results in a limiting factor in the use of this technique by other research groups which cannot access to very expensive computing centers (assuming that authors provide enough evidence of the equilibration of the system). Second, to avoid wasting expensive computing time in the study of certain properties which do not require such long equilibration times, once the steady state of the system has been properly identified.
As mentioned above, the Molecular Dynamics (MD) simulations have emerged during the last decades as a powerful tool to obtain insight with atomic detail of the structure and dynamics of lipid bilayers [34][35][36]. Several MD simulations of membranes under the influence of different salt concentrations have been carried out. One of the main obstacles related to these studies has been the time scale associated to the binding process of ions to the lipid bilayer. Considering the literature, a vast dispersion of equilibration times associated to the binding of ions to the membrane has been reported, where values ranging from 5 to 100 ns have been suggested for monovalent and divalent cations [25, 27-29, 32, 37, 38]. In this regard, we carried out four independent simulations of a lipid bilayer formed by 288 DPPC in aqueous solutions, for different concentrations of CaCl 2 to provide an overview of their equilibration times. Among other properties, the surface area per lipid, the deuterium order parameters, lipid hydration and lipid-calcium coordination were studied.
Finally, in order to generalize our results, a bilayer formed by 288 DPPS in its liquid crystalline phase, in presence of CaCl 2 at 0.25N, was simulated as well .

II. Methodology
Different molecular Dynamics (MD) simulations of lipid bilayer formed by 288 DPPC were carried out in aqueous solutions for different concentrations of CaCl 2 , from 0, up to 0.50 N. Furthermore, with the aim of generalizing our results, a bilayer of 288 DPPS in presence of CaCl 2 at 0.25 N was simulated as well. Note that the concentration of CaCl 2 in terms of normality is defined as: where n equivalent grams = gr(solute) equivalent weight and equivalent weight = Molecular weight n , being n the charge of the ions in the solution.
In Table 1, the number of molecules that constitute each system, applying Eq. (1), is summarized.
To build up the original system, a single DPPC lipid molecule, or DPPS lipid (Fig. 1), was placed with its molecular axis perpendicular to the membrane surface (xy plane). Next, each DPPC, or DPPS, was randomly rotated and copied 144 times on each leaflet of the bilayer. Finally, the gaps existing in the computational box (above and below the phospholipid bilayer) were filled using an equilibrated box containing 216 water molecules of the extended simple point charge (SPC/E) [39] water model.
Thus, the starting point of the first system of Table 1 was formed by 288 DPPC in absence of 040005-2 CaCl 2 . Once this first system was generated, the whole system was subjected to the steepest descent minimization process to remove any excess of strain associated with overlaps between neighboring atoms of the system. Thereby, the following DPPC systems in presence of CaCl 2 were generated as follows: to obtain a [CaCl 2 ]=0.06 N, 15 water molecules were randomly substituted by 5 Ca 2+ and 10 Cl − . An analogous procedure was applied to the rest of the systems, where 36, 69 and 138 water molecules were substituted by 12, 23 and 46 Ca 2+ and 24, 46 and 92 Cl − , to obtain a [CaCl 2 ] concentration of 0.13 N, 0.25 N and 0.50 N, respectively. Finally, the DPPS bilayer was generated following the same procedure described above for the DPPC, starting from a single DPPS molecule and once the lipid bilayer in presence of water passed the minimization process, 324 water molecules were substituted by 204 Ca 2+ and 120 Cl − (note that 144 of the 204 calcium ions were added to balance the negative charge associated with the DPPS). GROMACS 3.3.3 package [40,41] was used in the simulations, and the properties showed in this work were obtained using our own code. The force field proposed by Egberts et al. [2] was used for the lipids, and a time step of 2 fs was used as integration time in all the simulations. A cut-off of 1.0 nm was used for calculating the Lennard-Jones interactions. The electrostatic interaction was evaluated using the particle mesh Ewald method [42,43]. The real space interaction was evaluated using a 0.9 nm cut-off, and the reciprocal space interaction using a 0.12 nm grid with a fourth-order spline interpolation. A semi-isotropic coupling pressure was used for the coupling pressure bath, with a reference pressure of 1 atm which allowed the fluctuation of each axis of the computer box independently. For the DPPC bilayer, each component of the system (i.e., lipids, ions and water) was coupled to an external temperature coupling bath at 330 K, which is well above the transition temperature of 314 K [44,45]. For DPPS bilayer, each component of the system was coupled to an external temperature coupling bath at 350 K, which is above the transition temperature [46,47]. All the MD simulations were carried out using periodic boundary conditions. The total trajectory length of each simulated system was of 80 ns of MD simulation, where the coordinates of the system were recorded every 5 ps for their appropriate analysis.
Finally, in order to study the effect of the temperature, only the case corresponding to 0.25 N CaCl 2 was investigated at two additional temperatures, 340 K and 350 K.

III. Results and discussion
i. Effect of the CaCl 2 concentration a. Surface area per lipid Surface area per lipid A is a property of lipid bilayers which has been accurately measured from experiments [48]. The calculation of mean area per lipid can be determined from the MD simulation as: where x and y represent the box sizes in the direction x and y (perpendicular to the membrane surface) over the simulation, and N is the number of lipids contained in one leaflet, in our case N = 144.
Focusing on the study of the time evolution of the area per lipid, Figure 2  surface area per lipid for different concentrations of CaCl 2 and type of lipid. In general, for the 5 bilayers formed by DPPC or DPPS, the area per lipid achieved a steady state after 10 ns of simulation, being this equilibration time almost independent of the concentration of CaCl 2 and type of lipid which composed the membrane. In absence of salt, an average area per lipid of A = 0.663 ± 0.008 nm 2 was calculated from the last 70 ns of the simulated trajectory, discarding the first 10 ns corresponding to the equilibration time. This value agrees with experimental data, where values in a range from 0.55 to 0.72 nm 2 have been measured [10,11,[48][49][50][51]. Table 2 shows the mean surface area per lipid (again, after discarding the equilibration time of 10 ns) with their corresponding error bar. From the simulation results, a shrinking in the surface area per lipid with the increment of the ionic strength of the solution is observed. This shrinking is expected and attributed to the complexation of lipid molecules by calcium, such as it has been pointed out in previous studies [28,29,52].

b. Deuterium order parameter
The deuterium order parameter, S CD , is measured from 2 H-NMR experiments. This parameter provides relevant information related to the disorder of the hydrocarbon region in the interior of the lipid bilayers by measuring the orientation of the hydrogen dipole of the methylene groups with respect to the perpendicular axis to the lipid bilayer. Due to the fact that hydrogens of the lipid methylene groups (CH 2 ) have not been taken into account (in an explicit way) in our simulations, the order parameter −S CD on the i + 1 methylene group was defined as the normal unitary vector to the vector defined from the i to the i + 2 CH 2 group and contained in the plane formed by the methylene groups i, i + 1 and i + 2. Thus, the deuterium order parameter −S CD on the i − th of the CH 2 group can be estimated by Molecular Dynamics simulations as follows: where θ is the angle formed between the unitary vector defined above and the z axis. The expression in brackets . . . denotes an average over all  the lipids and time. Hence, note that the −S CD can adopt any value between -0.5 (corresponding to a parallel orientation to the lipid/water interface) and 1 (oriented along the normal axis to the lipid bilayer). Figure 3 shows the running −S CD for different carbons of the DPPC and DPPS tails and salt concentrations. Only the carbons which correspond to the initial (hydrocarbons 2 and 6), the middle (hydrocarbon 10) and final (hydrocarbons 13 and 15) methylene groups of the lipid tails were depicted in this figure. Each point of the figure represents the average values of −S CD over 5 ns of subtrajectory length, and the lines represent the mean values calculated from the last 70 ns of the trajectories simulated. From this figure, it is observed how in all the cases, the required equilibration time is less than 10 ns of simulation time, independently of the salt concentration and the type of lipid. Finally, it is noted that Figure 3 exhibits an increase in the deuterium order parameters with the salt concentration, consistent with the shrinking of the area per lipid described above.

c. Lipid hydration
To analyze the lipid hydration, the radial distribution function g(r) of water around one of the oxygens of the phosphate group (atom number 10 in Fig. 1 for DPPC and DPPS) was calculated. The radial distribution function g (r) is defined as follows: g(r) = N (r) 4πr 2 ρδr (4) where N (r) is the number of atoms in a spherical shell at distance r and thickness δr from a reference atom. ρ is the density number taken as the ratio of atoms to the volume of the total computing box. From numerical integration of the first peak of the radial distribution function, the hydration numbers can be estimated for different atoms of the DPPC or DPPS. Figure 4 depicts the hydration number of phosphate oxygen (atom 10 in Fig. 1 for DPPC and DPPS) in presence of CaCl 2 , where each point represents the average of 5 ns subtrajectory length. These results show how this property reached a steady state for the cases (A), (B) and (E), after 10 ns of simulation. However, for the cases (C) and (D), 5 ns of extra simulation trajectory were required to reach a steady state.  shows the hydration numbers for the last 70 ns of the trajectory length, corresponding to four concentrations of CaCl 2 and both types of lipids, DPPC and DPPS. In this regard, from Fig. 4, the significant lipid dehydration with the increment of the ionic strength of the solution is evident, in good accordance with previous results [52].

d. Phospholipid-calcium coordination
Some authors have reported how the lipid coordination by divalent cations widely varies . Thus, on the one hand, some authors [25] have reported that this is a very slow process, which requires about 85 ns of simulation time, but, on the other hand, other authors [26] have suggested that this process results much more rapid, taking less than 1 ns. In this sense, the coordination of DPPC-Ca 2+ was studied by monitoring the oxygen-calcium coordination of the carbonyl oxygens (atoms 16 and 35 in Fig. 1) and phosphate oxygens (atoms 9 and 10 of DPPC in Fig. 1), as a function of time. The left column in Fig. 5 represents the oxygen coordination number, while the right one depicts the percentage of calcium ions involved in the coordination process with respect to the total number of calcium ions present in the aqueous solution. Figure 5 shows how the DPPC coordination by calcium is a quick process, taking less than 5 ns of simulation time to achieve a steady state. The kinetic of this process appears to be related to the ratio between calcium/lipid. Thus, after the first 5 ns of simulation time, the Ca-lipid coordination presents some fluctuation along the rest of the simulated trajectory. However, in Fig. 5 (A) and (B) (for the cases of lower concentration), it is observed how the percentage of coordination fluctuates between a 60% and a 100%. We consider that this broad fluctuation is related to the limited sample size of our simulations that introduces a certain noise in our results.

ii. Effect of temperature
This section focuses on the study of the role played by temperature on the equilibration process. In this regard, only the system corresponding to a  concentration of 0.25 N in CaCl 2 was studied, for a range of temperatures from 330 K to 350 K (all of them above the transition temperature of 314 K [44,45] for the DPPC). Figure 6 shows the running area along the trajectory. In this case, it was noticed how the systems achieve a steady state after a trajectory length of roughly 10 ns, where Table 3 shows the mean area per lipid calculated from the last 70 ns of simulation time. Figure 7 shows the deuterium order parameter of the methylene groups along the lipid tails, calculated from Eq. (3). Figure 7, on the one hand, clearly shows that for the three temperatures the systems have reached the steady state before the first 10 ns of simulation. On the other hand, it shows an increase in the disorder of the lipid tails with temperature, which is closely related with the increase of the area per lipid, such as it was pointed out above. Figure 8 depicts the results of the hydration numbers of DPPC for the three temperatures studied, where the equilibrated state was achieved after 10 ns of simulation time. Table 3 provides the calculated hydration numbers in the equilibrium, showing how the lipid hydration remained invariable with the rising of the temperature. Concerning the lipid-calcium coordination, Fig. 9 represents the lipid-calcium coordination number, and the right column represents the calcium that participates in the coordination expressed in percentage respect the total of calcium ions in solution. From simulation, it becomes evident how calcium ions required less than 5 ns to achieve an equilibrated state for the three temperatures studied. In summary, for all the properties studied in this section, a slight decrease in the equilibration time with the increasing temperature was appreciated.

IV. Conclusions
The present work deals with the simulation time required to achieve the steady state for a lipid bilayer system in presence of CaCl 2 . In this regard, we studied two different systems: one with DPPC and another one with DPPS bilayer; both systems in  presence of CaCl 2 (at different level concentration). The salt free case was also studied, as control. The analysis of various lipid properties studied here indicates that some properties reach the steady state more quickly than others. In this sense, we found that the area per lipid and the hydration number are slower than the deuterium order parameter and the coordination of cations. Consequently, to ensure that a system composed by a lipid bilayer has reached a steady state, the criterion that we propose is to show that the area per lipid or the hydration number have reached the equilibrium. From our results, two important aspects should be remarked: 1. The equilibration time is strongly dependent on the starting conformation of the system. Wrong starting conformations will require much longer equilibration times, even of one order of magnitude higher than the requested from a more refined starting conformation.
2. Temperature is a critical parameter for reducing the equilibration time in our simulations, due to the fact that higher temperatures increase the kinetic processes, i.e., the sampling of the configurational space of the system.