The ELBA Force Field for Coarse-Grain Modeling of Lipid Membranes

A new coarse-grain model for molecular dynamics simulation of lipid membranes is presented. Following a simple and conventional approach, lipid molecules are modeled by spherical sites, each representing a group of several atoms. In contrast to common coarse-grain methods, two original (interdependent) features are here adopted. First, the main electrostatics are modeled explicitly by charges and dipoles, which interact realistically through a relative dielectric constant of unity (). Second, water molecules are represented individually through a new parametrization of the simple Stockmayer potential for polar fluids; each water molecule is therefore described by a single spherical site embedded with a point dipole. The force field is shown to accurately reproduce the main physical properties of single-species phospholipid bilayers comprising dioleoylphosphatidylcholine (DOPC) and dioleoylphosphatidylethanolamine (DOPE) in the liquid crystal phase, as well as distearoylphosphatidylcholine (DSPC) in the liquid crystal and gel phases. Insights are presented into fundamental properties and phenomena that can be difficult or impossible to study with alternative computational or experimental methods. For example, we investigate the internal pressure distribution, dipole potential, lipid diffusion, and spontaneous self-assembly. Simulations lasting up to 1.5 microseconds were conducted for systems of different sizes (128, 512 and 1058 lipids); this also allowed us to identify size-dependent artifacts that are expected to affect membrane simulations in general. Future extensions and applications are discussed, particularly in relation to the methodology's inherent multiscale capabilities.

Shifted-force charge-charge potential

Shifted-force charge-dipole potential
Assuming i is the charge and j is the point-dipole: Since i is a point-particle, T ij = 0.

Figure 2.
Orientation-restraining potential. The "free" vector e and the "reference" vector v form the angle x.

Angle-bending
The potential U (A) associated with bond angle variation is [3]: with w the rigidity constant, A the actual angle and A 0 the reference angle (see Figure 1). Forces and torques can be found in the original reference [3, pp. 283-285].

Orientation-restraining potential
Consider two unit vectors v and e forming an angle x between them (Fig. 2). To restrain the "free" vector e to lie along the direction defined by the "reference" vector v we can use the potential U : with c a rigidity constant and cos x = v ·e. This orientation-restraining interaction is required to generate a torque on e to favor alignment along v. Such a torque T e is defined by [2, p. 333]: Applying the chain rule, we can write: Solving the derivatives, we obtain: and: The final expression is therefore: Note that T e always induces a rotation that tends to align e with the reference vector v. Also note that there is no force arising from the potential of equation 18. In fact, the interaction force is in general obtained from the gradient of the potential with respect to the vector r representing the distance from the origin. Since in this case the potential does not depend on such a vector, the interaction force is zero (f = −∇ r U = 0).

Dielectric constant of the ELBA water model
The ELBA force field does not include long-range interactions. Since these are necessary to calculate the dielectric constant, for this purpose we modified the standard ELBA potential U by the addition of a reaction field potential U RF ij , so that the new total potential energy U * becomes: where U includes the various components described in the main body of the paper.
Switched reaction field potential We consider a standard dipole-dipole reaction field potential [4,5] supplemented by a cubic switching function s ij : where r 3 α is a modified radius which takes into account the effects of a switching function, and the other symbols have the same meaning as in the equations for the dipole-dipole potential reported previously. In particular, r 3 α is defined as [6]: Considering our choice of switching function, equation 26 becomes: The dielectric constant of the of the continuum outside r c is set to infinity, that is, ǫ RF = ∞. The pair force is: The pair torques are:

Calculation methodology
The calculation of the dielectric constant was carried out using two standard methods. In particular, using the "fluctuation method", the dielectric constant can be obtained as [7]: where M is the total dipole moment of the system, V is the volume of the simulation region, and T is the temperature. The total dipole moment of the system is computed with: where µ i is the dipole moment of site i. Using the polarization method, the dielectric constant can be obtained as [8]: where P z is the system dipole moment along the direction of the applied field E ext z . We recall that the potential energy of a point-dipole µ in an external electric field E is: with θ the angle between the two vectors. This generates a torque on the point-dipole:

Results
To estimate the dielectric constant of the ELBA water model, we carried out simulations of water systems comprising 4000 sites. Unfortunately, we found that the addition of the RF potential (equation 25) to the ELBA model causes an unphysical phenomenon; in particular, the molecular dipoles align along a preferential orientation, incidentally causing the dielectric constant calculation to break down. Such artificial ordering effects could be removed by reducing the switching radius r s of the RF potential; the results obtained are collected in Table 1.
The dielectric constant calculations are carried out with both the fluctuation ("Fluct") and the polarization ("Pol") methods. For the polarization method, the external electric field E ext z (see equation 32) was set to 0.01 V.
It is clear that the results are highly sensitive to the cutoff treatment of the reaction field potential; such a finding was not completely unexpected, and in fact a similar phenomenon has been observed previously in the literature [9][10][11]. In particular, it can be seen from Table 1 that increasing values of the RF switching radius (from 0.6 r c to 0.8 r c ) result in increasing values of the dielectric constant. For switching radii set to 0.9 r c and 1.0 r c (corresponding to straight cutoff), the calculation of the dielectric constant breaks down, in that the reaction field causes the water dipole to align along a preferential direction. In general, any truncation scheme effectively changes the potential, and hence it is bound to alter the physical properties of the model. The ELBA model seems to be especially sensitive to the additional reaction field potential; while this behavior is not ideal, it only affects the calculation of the dielectric constant. For all other purposes, the ELBA force field does not require the addition of a reaction field.
An alternative method to estimate the dielectric constant would involve treating the long-range dipolar interactions through a lattice summation technique [7]. While outside the scope of this paper, future work will be devoted to this issue; this will require the implementation of a dipole-dipole lattice summation scheme in our software.

Error analysis
The numerical properties reported in the "Results" section of the main body of the paper are reported in the form average ± standard error, where the average was calculated from the entire simulation, and the standard errors were obtained by dividing each trajectory into two consecutive blocks. Considering a property A, evaluated in a number M of measurements, the average A was obtained as: where A µ represents a single measurement. The following properties were evaluated at every step: A L , V L , µ HG , θ HG . In this case, the number of measurements M is obviously the same as the number of simulation steps. The following properties, and associated parameters, were instead evaluated every 100 steps (corresponding to 1.5 ps using a 15 fs timestep): electron density, lateral pressure, electrostatic potential. In this case, M = total number of steps/100. For every property, two subaverages A 1 and A 2 were calculated over the two consecutive halves of the trajectory of each simulation. The standard deviation s d was then computed as: where N is the number of subaverages. The standard error s e was eventually computed with: We are aware that, for the error estimation to be reliable, the subaverages must be statistically independent; this can be proved rigorously using a block averaging (BA) procedure, as described by Flyvbjierg and Petersen [12]. The BA method should in principle identify the longest correlation time present in the data, thus providing us with a reasonable value for the minimum block length that guarantees statistical independence of the block averages. Unfortunately, while this procedure is relatively straightforward for simple systems where only a few quantities of interest are calculated, it becomes problematic in cases like ours. The various membrane properties calculated are expected to be characterized by different correlation times, also in relation to the specific system size, temperature and presumably hydration; hence the optimal block length would have to be estimated independently for each property in each simulation. In fact, the BA method is not strictly guaranteed to work for large and heterogeneous systems such as lipid bilayers, as the possible existence of very long correlation times might prevent the method from converging over typical simulation timescales. In general, for all these reasons, we should note that BA is hardly ever applied in its rigorous form to the analysis of molecular dynamics data from complex biomolecular simulations. Overall, we believe that our sampling times are long enough to provide reliable averages and error estimations, as they are 1-2 orders of magnitude larger than what is normally achieved by standard atomistic models.

DOPE electron density profiles
The electron density profiles calculated from our simulations of DOPE bilayers are reported in Figure 3.  A slight size dependence can be noticed; the origin of this artifact is discussed in the main body of the paper.  A qualitative comparison between typical conformations of single lipids in the fluid and gel phase is displayed in Figure 5. The gel-phase conformation is characterized by parallel and "stretched out" tails, in contrast with the "disordered" fluid-phase conformation.

Electrostatic potential profiles
The electrostatic potential profiles calculated from our simulations of DOPC, DOPE and DSPC are reported respectively in Figure 6, 7 and 8.

Diffusion coefficients
The diffusion coefficients as a function of the measurement time for runs A, D G and H are reported in Figure 9.  It can be seen that the lateral diffusion coefficient for gel DSPC is 2 − 3 orders of magnitude lower than those obtained for fluid bilayers. Table 2 collects the ratios between water permeability coefficients through PC and PE bilayers from our calculations and from the measurements reported in the literature.  ,D), 2.8 (runs B,E), 3.6 (runs C,F) Exp [13] 16  Table 3 collects the ratios between water permeability coefficients through gel and fluid DSPC bilayers from our calculations and from experimental measurements reported in the literature. Pre-assembled Self-assembled Figure 10. DOPC electron density profiles. Comparison between electron density profiles obtained from self-assembled and pre-assembled systems. Pre-assembled Self-assembled Figure 11. DOPC lateral pressure profiles. Comparison between lateral pressure profiles obtained from self-assembled and pre-assembled systems. Pre-assembled Self-assembled Figure 12. DOPC electrostatic potential profiles. Comparison between electrostatic potential profiles obtained from self-assembled and pre-assembled systems.

Water permeation
5 Sensitivity to timestep size

Energy conservation
The total and potential energies over 1000 steps for different timestep sizes are plotted in Figure 13.  It can be seen that the energy is conserved up to a timestep size of 15 fs.  Figure 14. DOPC electron density profiles. Comparison between electron density profiles obtained from self-assembled and pre-assembled systems.

Benchmarks
This section presents an assessment of the computational efficiency of the ELBA coarse-grain model. The model was run with the serial program Brahms [15] and with the parallel program Lammps [16]. The benchmark system employed comprises 1152 lipids and 17328 water sites. To compare the ELBA performances with standard atomistic models, we prepared two roughly equivalent systems using a unitedatom (UA) and an all-atom (AA) force field [17,18]. These atomistic systems were run with the popular Gromacs software [19]. The various simulation parameters (timestep, cutoff radii, long-range schemes) were set following standard protocols [18]. All data reported in this section were obtained on the Iridis supercomputer [20] using identical processors (2.27 Ghz Nehalem).
For each system, we quantified the efficiency by calculating the "intrinsic speed" of the molecular dynamics simulation as: MD speed = # molecular dynamics steps second (38) where "second" is intended in terms of calculation (wall-clock) time. Moreover, we calculated the overall sampling efficiency as: Sampling speed = # nanoseconds sampled day (39) where "day" is intended in terms of calculation (wall-clock) time.
The results obtained for single-processor simulations are collected in table 4. Further tests were carried out to investigate the "parallel performances" of the same calculations when run on multiple processors using standard domain-decomposition message-passing methods (this was not done for the "ELBA-Brahms" simulation because Brahms does not run in parallel). These tests were run on the Iridis supercomputer [20], which is equipped with 8-processor nodes (each processor being a 2.27 Ghz Nehalem). Figure 18 reports the parallel efficiency as a function of the number of processors used (where "parallel efficiency" is the ratio of ideal to actual run time 1 ). Figure 19 reports a comparison in terms of sampling speed (as defined in equation 39).