Evaluating Polarizable Biomembrane Simulations against Experiments

Owing to the increase of available computational capabilities and the potential for providing a more accurate description, polarizable molecular dynamics force fields are gaining popularity in modeling biomolecular systems. It is, however, crucial to evaluate how much precision is truly gained with increasing cost and complexity of the simulation. Here, we leverage the NMRlipids open collaboration and Databank to assess the performance of available polarizable lipid models—the CHARMM-Drude and the AMOEBA-based parameters—against high-fidelity experimental data and compare them to the top-performing nonpolarizable models. While some improvement in the description of ion binding to membranes is observed in the most recent CHARMM-Drude parameters, and the conformational dynamics of AMOEBA-based parameters are excellent, the best nonpolarizable models tend to outperform their polarizable counterparts for each property we explored. The identified shortcomings range from inaccuracies in describing the conformational space of lipids to excessively slow conformational dynamics. Our results provide valuable insights for the further refinement of polarizable lipid force fields and for selecting the best simulation parameters for specific applications.


Introduction
Classical molecular dynamics (MD) simulations are nowadays widely and almost routinely used to model wide range of biomolecular complexes. 1In conventional MD simulations, electrostatic interactions are described by assigning static point charges to the atoms and molecules at the beginning of the simulation as defined in the inputted model-the force field.Dynamic effects arising from electronic polarizability are traditionally not explicitly included and only considered in an averaged fashion within the force field parameterization process where parameters are obtained by fitting to macroscopic observables or to ab initio calculations.However, electronic polarization is perceived to be a key contribution to correctly describe many biomolecular systems.These include description of water, ion hydration and ion binding to molecules, cation-π and π-π interactions, 2 and general co-operativity in interactions. 3These lower-level interactions also translate to the behaviour of larger-scale biomolecular systems, such as ion channels, where ion-selectivity and ion currents may be affected by polarization, [4][5][6][7] and telomeric DNA 8 where the conformations adopted are mediated by ionic interactions.[11][12][13][14][15][16] In a bilayer membrane, specifically, the molecular (dielectric) environment varies dramatically while crossing the membrane from the water phase across the dipolar/charged lipid headgroup interface to the hydrophobic tail region.8][19][20][21][22][23] However, the quality of polarizable lipid models have not been evaluated on an equal footing with non-polarizable models.
The underlying strategy for including explicit polarizability differ among these models: 1) the classical Drude oscillator (CHARMM-Drude) models polarization by separating a core and a shell charge with a spring which stretches according to the environment giving the site a fluctuating dipole moment, 14 2) the induced point dipole/multipole, Atomic Multipole Optimized Energetics for BiomoleculAr (AMOEBA), uses polarizable point dipoles placed on chosen sites of the molecule, 27 and 3) the electronegativity equilization (FQ) employs atomic charges that are no longer constant but can change (redistribute within the molecule) according to the electronegativities of the molecule atoms and the electric fields from the molecular environment. 28All of these approaches result in an increasing computational cost, e.g, by introducing new types of interactions, more interactions sites, or by requiring a shorter time step in the simulations.As an alternative approach, electronic continuum correction (ECC) has been proposed as a computational efficient method to implicitly include polarizability by scaling the atom partial charges. 21,291][32][33][34][35][36] The ability to capture these membrane properties correctly is important on its own but also creates the basis for the description of more complex systems: for example, ion binding affinity regulates membrane surface charge and having wide variety of conformations available for lipid headgroups appears to be important for capturing realistic protein-lipid interactions. 33Therefore, such benchmark studies are urgently needed also for polarizable lipid force fields, considering the increased computational cost they come with, and most importantly, their pledge to capture a broader range of physical phenomena pertinent to enhancing our understanding of in particular polar membrane regions.
Here we assess the quality of the actively developed polarizable lipid force fields, CHARMM-Drude 17,24 and AMOEBA-based 18,25 parameters, using the resources and the framework from the NMRlipids open collaboration (http://nmrlipids.blogspot.fi/).These two force fields were selected for comparison because they are increasingly utilized in biomolecular simulations and have parameters available for several lipids for which corresponding experimental data is available in NMRlipids databank. 36We assess the structural quality of POPC, DOPC and POPE lipid bilayer simulations against experimental NMR and small angle X-ray scattering (SAXS) data using the quality metrics defined in the NMRlipids databank. 36Cation binding to membranes is evaluated against salt-induced changes NMR order parameters 31 and lipid conformational dynamics is benchmarked to data from NMR spin relaxation rate experiments. 34,37Furthermore, for each experimental benchmark we compare the results from polarizable models to those from the best performing non-polarizable simulations obtained from the NMRlipids Databank. 36Our results will act as a useful reference for selecting the best polarizable lipid models for wide range of applications and for guiding the development of polarizable force field parameters for lipids.

Using a polarizable force field for membrane simulations
While non-polarizable MD simulations of membranes can be nowadays routinely performed with several simulation packages and force fields, polarizable simulations bear many practical complications.Out of the currently popular MD simulation packages for membrane simulations, OpenMM 38 supports both AMOEBA and Drude force fields, NAMD 39 can only run the Drude force field, whereas GROMACS 40 has only limited support for the Drude polarizable force field via an unofficial Git-branch. 41TINKER is widely used with AMOEBA, but it does not support semi-isotropic pressure coupling, which is required in membrane simulations.Consequently, we have selected OpenMM for simulations in this work.
Another practical issue is availability of parameters for the molecules of interest.CHARMM-GUI can generate the topology and input parameters for the Drude force field, which greatly simplifies the employment of this model in particular. 42Standard protocols for AMOEBA lipid parameters are not available, but some parameters can be found online. 5,43stly, while conducting polarizable membrane simulations, one should consider the increased computational cost arising from the explicit treatment of electronic polarizability.
For the Drude-based models, this slowdown occurs both because the addition of the Drude particles increases the number of interaction pairs and the employed extended dual-Langevin thermostat requires a shorter 1 fs integration time step, compared to the 2 fs typically used for non-polarizable membrane simulations.The AMOEBA force field can use a multi timestep integration algorithm where the non-electronic interactions are iterated with a 2 fs time step whereas more computationally unstable polarization terms are iterated with a smaller time step.However, the multi-time step scheme only partly mitigates the computational cost: one typically experiences almost a ten-fold decrease in the simulation speed while using the polarizable force fields.
All simulations performed in this work are listed in Table 1 with links to the openly available trajectory data.Data not mentioned here was obtained by analysing pre-existing trajectories from the NMRlipids databank.

Simulations with CHARMM-Drude parameters
CHARMM-Drude2017 simulations were performed with OpenMM 7.5.0 38using parameters extracted with Membrane Builder [44][45][46][47] and Drude Prepper 42 from CHARMM-GUI. 48,49Before starting the simulations, membrane structures were equilibrated for 200 ns using the non-polarizable CHARMM36 force field, 50 and the last frames of these simulations were used to generate the starting structures for the polarizable force field simulations.Ion parameters are obtained from Ref. 51 and SWM4-NDP water model is employed in all Drude simulations. 52ARMM-Drude2023 force field parameters are currently not integrated into CHARMM-GUI.Therefore, the simulation setups with NaCl and CaCl 2 51 using SWM4-NDP water model 52 have been generated following the instructions in the original CHARMM-Drude2023 paper 24 using CHARMM program 53 and the last frame of 200 ns long CHARMM36 simulations (same for CHARMM-Drude2017).The salt-free CHARMM-Drude2023 simulations were obtained from Zenodo. 54,55dual Langevin thermostat was employed to keep the Drude particles and the rest of the system at 1.0 K and 303 K, respectively.A Drude-hardwall of 0.02 nm was utilized to keep the Drude particles close to their parent atoms.Semi-isotropic Monte Carlo barostat was used to couple pressure to 1 bar independently in membrane plane and normal directions.Bonds containing the hydrogens were constrained.For CHARMM-Drude2017, particle mesh Ewald was used to compute the Coulomb interactions and Van der Waals interactions were set to 0 between 1.0 nm and 1.2 nm using a switching function.For the CHARMM-Drude2023, Lennard-Jones particle-mesh Ewald (LJ-PME) method was used to compute the long-range dispersions. 56Simulation frames were saved in every 10 ps.

Simulations with AMOEBA-based parameters
For simulations with the AMOEBA-based force field, we used the OpenMM implementation of parameters developed by Chu et al. 25 available at GitHub. 5,43All AMOEBA simulations were run using OpenMM 7.1.1. 38The same initial structures as in CHARMM-Drude simulations were used.A multi-time step Langevin integrator was used to iterate the bonded and non-bonded interactions in every 0.5 fs and 2.0 fs, respectively.A non-bonded cutoff of 1.2 nm was applied while semi-isotropic Monte Carlo barostat was used to couple pressure to 1 bar independently in membrane plane and normal directions.Ion and water parameters were obtained from Ref. 57.Simulation frames were saved in every 10 ps.Further simulations details can be found in the input files of the respective simulations (Table 1).

Analysis of simulations
All simulations were first uploaded to the NMRlipids databank. 36Area per lipids, small angle X-ray scattering form factors, and order parameters (S CH ) are automatically calculated by the NMRlipids Databank 36 and were extracted from therein.Quality evaluation metrics were quantified as detailed in the NMRlipids databank 36 with the exception that POPC simulations at 303 K were paired with the experimental data measured at 300 K, while in the NMRlipids databank simulations are paired with experiments with the maximum temperature difference of two degrees.The order parameter qualities, P headgroup , P sn−1 , and P sn−2 , presents an average probability for order parameters of a given segment to locate within experimental values, taking the error bars of both values into account.Qualities of SAXS form factors, F F q , depict the distance of the first form factor minima between simulation and experiment.This avoids the effects from the simulation size dependency of the SAXS data. 36Note that in quantifying the bilayer electron densities for calculation of SAXS curves, the NMRlipids databank analysis algorithm places electrons as point charges at the atom centers without considering redistribution of charge density due to the polarizability.
Nevertheless, we expect this approximation not to have significant effect on the resulting SAXS form factors.
The density profiles were calculated using MDAnalysis 58,59 /NumPy 60 for CHARMM-Drude and AMOEBA simulations, and gmx density command from Gromacs package for the CHARMM36 and ECClipids simulations, for which the trajectories where extracted from the NMRlipids databank.In all calculations, the center of mass of the system is translated to origin coordinates.All data is normalized such that the area under the curve is equal to unity.The R 1 rates and effective correlation times, along with the accompanying error estimates, were quantified from the trajectories as elaborated in Ref. 34.In-house python script available at https://github.com/NMRLipids/NMRlipidsVIpolarizableFFs/tree/master/scripts/correlation_times was used for this purpose.3 Results and Discussion

Evaluation of lipid bilayer structural properties
To evaluate the structural properties of lipid bilayers in simulations with polarizable force fields, we simulated POPC and POPE lipid bilayers with CHARMM-Drude2017 17 and CHARMM-Drude2023 24 parameters, and POPE and DOPC bilayers with AMOEBA-based 18,25 parameters.These systems were selected based on simultaneous availability of force field parameters and experimental data. 36After the data production, simulation trajectories were added into the NMRlipids Databank and its quality metric were used to evaluate each simulation against experiments. 36The metric measures the quality of lipid bilayer simulations in two aspects: first, it evaluates the quality of conformational ensemble of individual lipids against C-H bond order parameters from NMR, and second, it examines consistency of membrane dimensions compared to small angle X-ray scattering data. 36,83The former is divided into three parts that separately describe the average quality of headgroup, glycerol backbone and acyl chain regions.While order parameters primarily reflect lipid conformations, acyl chain order parameters are also a good proxy for membrane packing since a smaller area per lipid tends to increase the magnitude of the order parameters. 36rect comparisons between simulated and experimental data are visualized in Fig 1 and the resulting quality metrics are shown in Table 2. Differences of headgroup conformations between force fields are shown in terms of dihedral distributions in Supplementary Information Fig. S2.Similarly to the non-polarizable CHARMM36 simulations, CHARMM-Drude2017 simulation predicts slightly too packed membrane with excessively negative order parameters with respect to experiments and simulations with the highest quality in the NMRlipids databank (OPLS3e). 36However, the quality of headgroup conformations in the CHARMM-Drude2017 simulations is worse than in its non-polarizable counterpart.This is likely because headgroup and glycerol backbone parameters were optimized to reproduce absolute average values of the experimentally determined order parameters, instead of taking   36 and area per lipids (APL) compared between the best simulations found from the NMRlipids databank (OPLS3e 84 for POPC and GROMOS-CKP 85-87 for POPE), simulations with CHARMM36 force field, and simulations with polarizable force fields.Order parameter quality measures, P headgroup , P sn−1 , and P sn−2 , present the average probability of order parameters within a given segment to agree with experiments (larger probabilities mean higher quality).Form factor quality measure, F F q , presents the difference in essential features between simulated and experimental form factors (smaller values indicate higher quality).Experimental estimations for areas per lipid are POPC: 64.3±1Å C-H bonds at a single carbon atom) and the order parameter sign. 35A better description of the headgroup and glycerol backbone is provided by CHARMM36-Drude2023, yet its quality is still below that offered by the non-polarizable CHARMM36 parameters (see also discussion about conformational dynamics below).Also the quality of membrane packing and acyl chain order are improved in CHARMM36-Drude2023 compared to the earlier version but again, it is outperformed by the best available simulations (OPLS3e in Table 2).
AMOEBA-based simulations capture the headgroup and glycerol backbone order parameters reasonably well with the exception of g 1 , where forking of order parameters is unacceptably large.However, experimentally observed low order parameters at double bond region in both DOPC acyl chains and sn-2 chain of POPE are not captured.These low order parameters signal an important mechanism through which acyl chain double bonds affect membrane properties 94 and are well reproduced in all state-of-the art MD simulation parameters. 83Furthermore, AMOEBA-based parameters substantially overestimate the area per lipid in POPE simulations, which is connected to too disordered acyl chains.Similar issues are evident also in the membrane data presented in the recent publication for the AMOEBA-based cholesterol model. 95The unsatisfactory description of the lipid tail region and the area per lipid is further reflected in the inability of AMOEBA-based parameters to capture the SAXS curves.Thus, we can conclude that the AMOEBA-based parameters used in our simulations did not reproduce essential membrane properties at the level of the state-of-the-art lipid parameters.

Evaluation of conformational dynamics of lipids
C-H bond order parameters are sensitive to the conformational ensemble but the correspondence is not unique.Instead, they essentially describe the average of the distributions and do not contain information of the dynamics of the conformational sampling.Therefore, a simulation that reproduces the order parameters has an ensemble that is only potentially correct, and even the correct ensemble may not be sampled at the experimentally observed dynamics.To further elucidate this aspect in the case of polarizable force fields, we compared their 13 C NMR spin relaxation rates, R 1 , and C-H bond effective correlation times, τ eff , with experiments 37 and best non-polarizable simulations from our previous study 34 in Fig. 2.
Here, we focus on PC headgroups and the glycerol backbone due to the availability of both experimental data and polarizable simulations.R 1 rates measured with typical magnetic field strengths are sensitive to rotational dynamics of C-H bonds with the timescales around ∼0.1-1 ns, while effective correlation times respond to a larger range of dynamical processes up to ∼1000 ns. 37The effective correlation time can be interpreted as an average measure for how fast molecular conformations sample the phase space that leads to the average order parameters.
The effective correlation times in CHARMM-Drude2017 and CHARMM-Drude2023 simulations are approximately two and one orders of magnitude slower, respectively, than the values extracted from experiments and best available simulations (Fig. 2).This indicates that not only are the simulations with this polarizable force field computationally costlier for equivalent lengths of trajectory, but one would also have to simulate longer to obtain converged results.The non-polarizable counterpart of the Drude models, CHARMM36, exhibits much more realistic, faster effective correlation times.Also the dynamics in the fast range (R 1 rates), are on average slightly more realistic in CHARMM36 simulations compared to both of its polarizable counterparts.Interestingly, CHARMM-Drude models have been reported to have slower water hydrogen bonding dynamics around amino acids compared to their non-polarizable counterpart, 96 which might align with an overall slower dynamics of the model in addition to enhanced water binding.
In contrast to the Drude-based models investigated here, the R 1 spin relaxation times and effective correlation times from DOPC simulations with AMOEBA-based force field are close to experimental data from POPC and on par with the best non-polarizable models (Slipids and CHARMM36).The small difference in acyl chain composition (DOPC vs. POPC) is not expected to affect headgroup dynamics due to effective decoupling between hydrophilic and hydrophobic membrane regions. 97,98

Cation binding to membranes in polarizable simulations
Given the abundance of cations in biological systems, accurately capturing their interactions with membranes in simulations has the uttermost importance.A wealth of experimental evidence shows that monovalent ions (except for lithium) exhibit very weak binding affinity of PC lipid bilayers, while multivalent ions like calcium bind more strongly. 31However, simulations with non-polarizable force fields without any additional corrections systematically overestimate cation binding to lipid bilayers. 31Implicit inclusion of electronic continuum correction (ECC) to both ion and lipid parameters can substantially improve the situation, 21,22,33 suggesting that electronic polarizability plays an important role in ion binding to membranes.
This might further lead one to expect that simulations with explicitly polarizable force fields will more accurately describe ion binding to membranes.To test this notion, we evaluated For simulated τ eff , the data point quantifies the average over C-H bonds.If τ eff could not be determined for all bonds due to slow convergence with respect to simulation time, only the range from the mean of the lower to the mean of the upper error estimates is shown.All simulations used here were salt free.
ion binding to membranes against the experimental NMR data using the "lipid electrometer", in which the amount of ion binding to the membrane is quantified by monitoring the change in lipid headgroup order parameters in response to increasing salt concentration. 31,99anges in these order parameters, as induced by increasing NaCl or CaCl 2 concentration, are shown in Fig. 3 for simulations and experiments, whereas Fig. 4 shows the density profiles of ions with respect to the bilayer normal.Results from AMBER-based ECClipids model are presented as a reference which gives a good agreement with experiments for cation binding. 21 also show data from the non-polarizable CHARMM36 where NBFIX correction for the ion models was specifically developed to address overbinding.The CHARMM-Drude2017 predicts similar calcium ion density profile to ECClipids-a model which is in good agreement with experiments.However, the sodium binding is equally strong in CHARMM-Drude2017 simulations in contrast with the experimental evidence. 31d the response of headgroup order parameters to the bound ions is not in qualitative agreement with experiments ( Fig. 3).In particular, increasing α carbon order parameters and forking upon ion binding are not observed in experiments.This is in contrast with the

Conclusions
Including electronic polarizability to MD simulation models of membranes is expected to improve the description of bilayer polar regions and their interactions with charged molecules, thereby making MD simulations of complex biomolecular systems more realistic.However, quality of polarizable membrane simulations have not been evaluated on equal footing with non-polarizable ones.Here, we used the quality evaluation metrics defined in the NMRlipids databank 36 together with additional analyzes on dynamics and ion binding to evaluate the performance of two available polarizable lipid models, CHARMM-Drude 17,24 and AMOEBAbased 18,25 parameters, against experimental data and results from the best-performing nonpolarizable force fields.Considering the complexity and additional computational cost of simulations with polarizable models, it is crucial to understand their accuracy with respect to experiments and choose the best models according to their respective strengths when planning simulations.
Our comparisons of lipid conformations and dynamics ahow there is room for improvement in current polarizable parameterizations to reach the level of the best currently available non-polarizable parameters.Although the most recent CHARMM-Drude23 has improved the description of the molecular conformations and dynamics, both tested CHARMM-Drude models predict slightly too ordered membrane core and vastly too slow headgroup dynamics.
The tested AMOEBA-based parameters have difficulties capturing membrane ordering in acyl chain region and other more general membrane properties interconnected to chain ordering.Yet, the PC and PE headgroup conformations and dynamics have similar or almost similar quality than the best non-polarizable force fields.
Sodium and calcium binding to membranes in simulations were evaluated using experimentally observed headgroup order parameter changes upon addition of NaCl or CaCl 2 .The binding in explicitly polarizable models was compared with the ECClipids model with implicitly included electronic polarizability that currently gives the most accurate response to ion binding against experiments, and to the non-polarizable CHARMM36, which when used with the NBFix ion model designed to provide more accurate binding also rivals its polarizable counterparts.That said, CHARMM-drude2023 model provides an improved description of the stronger calcium binding compared to sodium.This difference between sodium and calcium binding is also present in simulations with AMOEBA-based parameters.However, calcium binding depth, affinity, and the consequent structural response of the lipids do not exactly align with experiments or ECClipids results neither in CHARMM-Drude2023 nor AMOEBA simulations.The incorrect response to ion binding likely connects to the other discussed inaccuracies related to lipid conformational ensembles, dynamics, and membrane order.
In summary, the full potential and promise of explicitly polarizable lipid force fields to improve the description of bilayer membranes have not yet been fully realized.However, it seems likely that this is not an inherent flaw in polarizability but rather in the current parameterization and their incombapility with parameters describing other interactions, such as those for Lennard-Jones interactions.Non-polarizable force fields benefit from a longer history of development and intense scrutiny, it is not surprising that they currently outperform their polarizable counterparts in many respects, as we here show for the case of certain phospholipid bilayers.On the other hand, improvements from CHARMM-Drude2017 to CHARMM-Drude2023 demonstrate the potential for tuning these parameters to improve polarizable force fields.Such endeavors are expected to substantially ease with emerging automated methods for parameter development. 24,35

Figure 1 :
Figure 1: X-ray scattering form factors, and C-H bond order parameters for headgroup and glycerol backbone, and acyl chains (from top to bottom) compared between simulations and experiments using the NMRlipids databank.Experimental data in this figure are originally reported in Refs.33,36,91-93.For CHARMM-Drude2023 simulations we selected representative replicas among three available one (for all POPC replicas, see SI Fig. S1).Right most panel shows the modelled lipids and the used carbon naming scheme.10

Figure 2 :
Figure2: Effective correlation times (top) and R 1 rates (bottom) for the polarizable models and best performing non-polarizable force fields.Note, top panel axis is logarithmic to visualize the slow Drude dynamics.Experimental values are obtained from Ref.97.For simulated τ eff , the data point quantifies the average over C-H bonds.If τ eff could not be determined for all bonds due to slow convergence with respect to simulation time, only the range from the mean of the lower to the mean of the upper error estimates is shown.All simulations used here were salt free.

Figure 3 :
Figure 3: The change in the lipid head group order parameters β (top row) α (bottom row) upon increasing ion concentration with respect to the simulations without salt.CHARMM36 and ECClipids data is reproduced using the Zenodo repositories at Refs. 102-105 and Ref. 106, respectively.

1 ]Figure 4 :
Figure 4: Density profiles along the membrane normal for the CHARMM36, 31 CHARMM-Drude2017, CHARMM-Drude2023, AMOEBA, and ECClipids. 21Data with circle markers and dashed lines in first and second column denote density profiles calculated from the water oxygens.Solid and dashed lines in the third column denote Cl −1 density from NaCl and CaCl 2 , respectively.Note that for CaCl 2 350 mM (AMOEBA and CHARMM36) and 450 mM (Drude models and ECC) concentrations are shown, while 1000 mM NaCl concentration is shown for all force fields, except CHARMM36 simulations are at 950 mM NaCl.CHARMM36 data is reproduced using the Zenodo repositories Ref. 102-105.Results are from POPC simulations for all force fields other than AMOEBA (DOPC).

Table 1 :
Systems simulated for this work.Column N l gives the number of lipids, N w gives the number of water molecules, T(K) denotes the temperature in kelvins.Simulated time is listed in column t sim and time used for analysis in t analysis .Column "files" gives the reference to open access simulation data.The salt concentrations calculated as [salt]=N c ×[water] / N w , where [water] = 55.5 M w N ion T (K) t sim (ns) t analysis (ns) files [Ref.]

Table 2 :
NMRlipids databank quality metrics 90 88 DOPC: 67.5±1Å 2 , 89 and POPE: 56.7±3Å 2 .90 100,101 318107 with circle markers and dashed lines in first and second column denote density profiles calculated from the water oxygens.Solid and dashed lines in the third column denote Cl −1 density from NaCl and CaCl 2 , respectively.Note that for CaCl 2 350 mM (AMOEBA and CHARMM36) and 450 mM (Drude models and ECC) concentrations are shown, while 1000 mM NaCl concentration is shown for all force fields, except CHARMM36 simulations are at 950 mM NaCl.CHARMM36 data is reproduced using the Zenodo repositories Ref.102-105.Results are from POPC simulations for all force fields other than AMOEBA (DOPC).resultsfrom prbechmarking of non-polarizable simulations,31where experimentally observed decrease of order parameters to more negative values upon ion binding were observed in all simulations, even though the binding affinity was often inaccurately predicted.This qualitative discrepancy in CHARMM-Drude2017 simulations may result from incorrect lipid headgroup conformational ensemble (see previous sections), which leads to inaccuracies in the structural response of the ensemble to ion binding.Excessive sodium binding in the CHARMM-Drude model has been observed before for systems containing peptides or amino acids96,107and deep-eutectic solvents.108InsimulationswithCHARMM-Drude2023 parameters, sodium and calcium ion binding are in line with ECClipids simulations when comparing density profiles, although calcium binding affinity is slightly larger, ions penetrate deeper into the bilayer, and there is a stronger formation of ionic double layer and deeper penetration of chloride into the bilayer in CHARMM-Drude2023 than in ECClipids.Interestingly, these features are observed in all simulations with calcium using polarizable force fields.However, sodium or calcium ion binding seems to again induce forking of lipid headgroup order parameters in CHARMM-Drude2023 simulations in contrast to experiments.This might indicate inaccurate structural response to ion binding, but poor convergence of the simulations owing to the slow conformational dynamics detected above cannot be ruled out in this case or in the case of the older In AMOEBA-based simulations, sodium binds weakly and does not affect order parameters, consistent with experiments and ECClipids simulations.Calcium binding affinity is similar to ECClipids, but order parameters do not change upon binding contrasting the experiments.This can be explained by the binding position of calcium, which is outside phosphate density peak in AMOEBA simulations.In other simulations calcium penetrates to phosphate region or deeper and reorients the headgroup dipole, giving rise to changes in order parameters in line with the experimental data.31Inconclusion, incorporating explicit polarizability as implemented in CHARMM-Drude or AMOEBA-based parameters used here do not necessarily lead to improved description of cation binding to phospholipid membranes.These parameters do not correctly capture the response of lipid headgroup to cation binding to membrane, which most likely arises from inaccuracies in lipid parameters.Because such inaccuracies can also affect ion binding, it is difficult to isolate the explicit influence of polarizability per se on ion binding based on the results presented here.
the ion binding by scaling the Lennard-Jones parameters (NBFIX).That said, the overall structural response to ion binding in CHARMM36 appears more realistic.