Charged Small Molecule Binding to Membranes in MD Simulations Evaluated against NMR Experiments.

Interactions of charged molecules with biomembranes regulate many of their biological activities, but their binding affinities to lipid bilayers are difficult to measure experimentally and model theoretically. Classical molecular dynamics (MD) simulations have the potential to capture the complex interactions determining how charged biomolecules interact with membranes, but systematic overbinding of sodium and calcium cations in standard MD simulations raises the question of how accurately force fields capture the interactions between lipid membranes and charged biomolecules. Here, we evaluate the binding of positively charged small molecules, etidocaine, and tetraphenylphosphonium to a phosphatidylcholine (POPC) lipid bilayer using the changes in lipid head-group order parameters. We observed that these molecules behave oppositely to calcium and sodium ions when binding to membranes: (i) their binding affinities are not overestimated by standard force field parameters, (ii) implicit inclusion of electronic polarizability increases their binding affinity, and (iii) they penetrate into the hydrophobic membrane core. Our results can be explained by distinct binding mechanisms of charged small molecules with hydrophobic moieties and monoatomic ions. The binding of the former is driven by hydrophobic effects, while the latter has direct electrostatic interactions with lipids. In addition to elucidating how different kinds of charged biomolecules bind to membranes, we deliver tools for further development of MD simulation parameters and methodology.

Introduction

Binding affinities of charged molecules, such as drugs, amino acids, ions, and pollutants, on cellular membranes, regulate many of their biological functions.[1−6] For example, signaling domains in proteins often contain charged residues and interact with membranes in an ion-dependent manner,[1,4−6] translocation of drugs through membranes depends on their charge state,[7] and bioaccumulation of charged pollutants can be related to their membrane affinity.[3,8] Such processes are particularly poorly understood for charged molecules because their binding to membranes is significantly more difficult to study than for neutral molecules. Therefore, better understanding and predictive models for membrane binding of charged molecules would benefit applications in a wide range of fields, such as designing drugs with better translocation properties, and understanding cell signaling and bioaccumulation of potentially toxic molecules.[3,6,7,9,10]

For neutral molecules, the water–oil (often octanol) partition coefficients correlate well with the binding data on model membranes and their binding affinities can be captured equally well by theoretical models with atomistic and continuum-level descriptions.[9,11−13] However, measurements of water–oil partition coefficients are problematic for charged molecules due to the charge neutralization in both phases.[14,15] Furthermore, experimental data on lipid binding affinity is more scarce, and a complex electrostatic environment at membrane interfaces, created by dipoles of lipids and water molecules, complicates theoretical analyses.[14,15] Incorporation of the membrane dipole potential in continuum models improves the results,[14] but atomistic resolution description is needed to fully capture the complexity of electrostatic and other interactions between lipids, water, and charged molecules at membrane interfaces.

However, the correct description of interactions between charged water-soluble molecules and membranes has been challenging for atomistic resolution models employed in molecular dynamics (MD) simulations. Canonical force fields tend to overestimate the binding affinities of sodium and calcium to membranes,[16,17] yet this can be improved by electronic continuum correction (ECC)[16,18−20] or additional repulsive potential (NBFIX).[17,21,22] Also lipid–protein interactions depend on force field parameters,[23−26] particularly for charged residues[27,28] and their interactions with lipid head groups.[29,30] Furthermore, charged residues in strongly membrane-bound peptides appear as potential sources of discrepancies in comparisons with NMR data[23,31] and often lead to the largest deviations from experimental hydrophobicity scales.[27,28,32] These findings raise the question of whether the interactions between charged molecules and membranes can be correctly captured in MD simulations without including electronic polarizability or other corrections. Therefore, it is currently not clear how accurately atomistic MD simulations can predict interactions between membranes and charged molecules, such as drugs, amino acids, or pollutants.

Sodium and calcium ion binding affinities to membranes with various compositions have been quantitatively evaluated in simulations using the changes in lipid head-group C–H bond order parameters,[16,20,33] which depend on the accumulation of charges on the membrane.[34] Here, we apply the same approach to quantify the binding affinities of etidocaine and tetraphenylphosphonium (TPP) ions to a phosphatidylcholine (POPC) lipid bilayer, serving as a model cellular membrane. Etidocaine is a clinically used local anesthetic and serves here as a model for charged drug binding to a membrane. TPP is a hydrophobic ion historically used to establish the concept of membrane dipole potential[35] and serves as a model for charged aromatic compounds which are common among potentially bioaccumulating ions.[3] Our results compare the relative binding affinities of these molecules to membranes with respect to sodium and calcium and test the quality of MD simulations regarding these affinities. Furthermore, we also investigate the binding mechanisms of these small molecules to membranes and the effect of implicit inclusion of electronic polarizability using ECC.

Methods

Force Field Parameters

For simulations with the standard force fields, CHARMM36 parameters[36] from CHARMM-GUI[37−39] were used for lipids. Parameters for etidocaine were generated with two standard approaches: SwissParam[40] from https://www.swissparam.ch/, denoted here as CHARMM36-SwissParam, and Cgenff (CHARMM general force field)[41−43] from https://cgenff.umaryland.edu/, denoted here as CHARMM36-ParamChem. For TPP, the standard approaches cannot be used since there are no parameters for the phosphonium atom. Therefore, we used CHARMM atom types together with the charges calculated using Gaussian using the density functional theory B3LYP functional and CHELPG scheme (charges from electrostatic potentials using a grid-based method),[44] denoted here as CHARMM36-Qmcharges. Partial charges from our Gaussian calculation slightly deviated from the ones reported in the literature,[45] see Figure S1 in the Supporting Information. Binding affinities to membranes were tested with both charges, but the results were similar, and only the results with our parameters from Gaussian (CHARMM36-Qmcharges) are shown. In another approach, denoted here as CHARMM36-ProteinFF, charges for TPP were obtained from the phenyl ring of phenylalanine in the CHARMM36 protein force field. All the systems were run with the TIP3P CHARMM36 water model.[46]

For simulations with ECC, we used our recently parameterized model compatible with CHARMM36 for POPC (PN-model in ref (47)). Following the simplest ECC approach,[48,49] all partial charges in CHARMM36-ParamChem parameters of etidocaine (Table S1 in the Supporting Information) and CHARMM36-Qmcharges of TPP (Figure S1 in the Supporting Information) were scaled by 0.75. Partial charges used in all etidocaine and TPP simulations are listed in Figure S1 and Table S1 in the Supporting Information.

Simulation Details

Lipid bilayers, consisting of 200 POPC molecules, were generated in CHARMM-GUI.[37−39] Structures of etidocaine and TPP were built in Chimera or downloaded from https://zinc.docking.org/. Lipid bilayers were then surrounded by the desired number of etidocaine or TPP, which were placed in a 3D grid with equidistant positions to fill the whole space of the selected box size. The systems were then solvated with water molecules and neutralized with chloride counterions using GROMACS tools. The simulations were started from two different kinds of starting structures. The small molecules were either placed in the water phase or a snapshot from an already existing simulation was used. The latter approach was used for systems where all molecules were bound to the membrane, as described below. The numbers of water molecules were tuned to produce the desired concentrations. Equilibration of the systems was monitored by analyzing the number of bound particles as a function of time. The simulation times of individual trajectories ranged between approximately 100 ns and 5 μs, see Tables S3–S13. For selected systems, we prepared also starting structures with all molecules bound in a membrane and confirmed that the number of bound molecules converged to the same value as from a starting configuration with all molecules in the water phase, see Figure S5. The exact number of molecules and other details are listed in Tables S3–S13.

Systems were run using GROMACS software (versions 2018.6, 2018.8, 2019.3, 2021.1, 2021.4, and 2021.5).[50,51] The temperature was coupled to 298 K corresponding to temperatures in experiments using a Nosé–Hoover thermostat.[52,53] Particle mesh Ewald (PME) was used to calculate electrostatic interactions at distances longer than 1.2 nm.[54,55] Lennard-Jones interactions were cut off at 1.2 nm. The pressure of 1 bar was maintained using a semi-isotropic Parrinello–Rahman barostat.[56]

Analysis of Simulations

C–H bond order parameters were calculated from equation SCH = ⟨3 cos θ – 1⟩/2, where θ is the angle between the C–H bond vector and the average is taken over the ensemble. Codes available in the NMRlipids project were used for order parameter calculations.[33] To estimate the number of bound molecules to the bilayer, we first calculated how strongly the number of bound molecules to the membrane depends on the selected criteria for the distance from lipids and fraction of bound atoms. Then, the values at the point of weakest dependence (minimum of the gradient) were selected; see the Supporting Information for details. The density distributions along the membrane normal were calculated using a histogram method. Atom coordinates were centered around the center of mass of the POPC lipid molecules for every time frame, and a histogram of these centered positions was calculated with the bin width of 1/3 Å. The script utilizing the MDAnalysis module[57,58] in Python is available at https://github.com/nencini/charged_molecules_binding. The potential of mean force (PMF) profiles were calculated from the density profiles ρ(z) using the inverse Boltzmann formula PMF = −kbT ln ρ(z). When calculating the number of contacts to specific regions (such as the phosphate group), two atoms were considered to be in contact when the distance between them was smaller than 0.325 nm. Permeation of particles through the membrane bilayer was analyzed by the program from Camilo et al.[59] Instead of water oxygen used in the original work, the N14 atom was used for the analysis of permeation of etidocaine and the P atom was used for TPP molecules.

Results and Discussion

Evaluating the Charged Small Molecule Binding to Membranes in MD Simulations Using Lipid Head-Group Order Parameters

Partition coefficients have been typically used to compare binding affinities of small molecules between MD simulations and experiments.[12,13,60] Such values are available also for some charged small molecules[15] but not for simple ions such as sodium or calcium. On the other hand, the affinities of ions to membranes have been reported using binding constants based on models assuming a certain binding stoichiometry for lipid–ion interactions, but these values strongly depend on the model used to interpret the experimental data.[61] For example, the reported binding constants of calcium to neutral POPC membranes vary between 7 and 441 M–1.[61] Furthermore, such binding constants may not even reflect the relative binding affinities of different molecules to membranes. This is exemplified in Figure where the number of bound TPP (reported binding constant of 21 M–1) and calcium (reported binding constant of 14 M–1) ions per lipid molecule are similar, while the number of bound etidocaine molecules (reported binding constant of 11 M–1) is larger at similar bulk concentrations.

Figure 1

(a) Number of bound molecules per POPC lipid as a function of equilibrium concentration in the supernatant and binding coefficients for different charged molecules from the literature. The number of bound particles is determined by the difference in supernatant concentration before and after the addition of lipids. Concentrations are determined by UV spectroscopy (SMS,[62] TPP,[63] etidocaine,[64] and dibucaine[64]) or atomic absorption spectroscopy (calcium ions[65]). The models used to determine binding coefficients and further details are shown in Table S2 in the Supporting Information. (b) Structures and charge distributions of the small molecules. For the chemical structure of etidocaine and TPP, see Figures S1 and S2 in the Supporting Information. (c) Structure of the POPC lipid molecule with the head-group α and β carbons labeled with blue and phosphate oxygens prone to calcium binding labeled with red.

The unambiguity in reported binding constants can be circumvented by directly comparing how lipid head-group C–H bond order parameters change upon the addition of ions between simulations and experiments.[16] These order parameters are proportional to the amount of charged molecules bound to the bilayer because binding of positive charges decreases their values (due to the orientation of the lipid head-group dipole more parallel to the membrane normal) and vice versa for negative charges.[16,34] As we are interested in the relative binding affinities between charged small molecules and biologically relevant simple ions, such as calcium and sodium, we use this approach here. From the available experimental data for head-group order parameter changes upon binding of charged molecules,[66] we chose to evaluate the binding affinities of etidocaine, TPP, dibucaine, and cyclic somatostatin (SMS) binding to a POPC lipid bilayer in simulations.

The approach has been successfully applied for sodium and calcium binding to various membranes,[16,18−20,33] but care must be taken to confirm that the simulation setup is sufficiently close to experiments, particularly in terms of molecular concentrations in the system and equilibration of the binding. While changing the hydration level with constant ion concentration did not affect the conclusions in previous sodium and calcium simulations,[18] the situation is more complex for small molecules with stronger binding affinities. The concentration of these molecules substantially decreases in bulk water upon binding to the membrane. In simulations with the hydration levels lower than in experiments, all solute molecules may bind to the membrane, thereby leading to the artificially underestimated bulk concentration. An ideal solution would be to use the same hydration level as in experiments, but this often leads to large boxes with substantial computational costs. Therefore, we estimated the minimal level of hydration, after which the results are not affected by the further addition of water while keeping the ion concentration in water constant.

To quantify the effect of hydration on the binding, we simulated systems with a fixed etidocaine concentration in simulation boxes with the z-dimension ranging from 5 to 81 nm. The lower limit corresponds to a typical size for membrane simulations, and the upper limit corresponds to the water–lipid ratio used in the experimental setup for etidocaine.[64] For parameters predicting the weakest binding affinity, CHARMM36-SwissParam, the results are similar with box sizes above approximately 12 nm in the z-direction (Figures S4 and S5 in the Supporting Information). However, in simulations with CHARMM36-ParamChem and CHARMM36-ParamChem-ECC predicting stronger binding affinities, the dependence on the box size is observed also above z-dimensions of 12 nm. Further simulations with these parameters were, therefore, run using the large box with the z-dimension of approximately 81 nm, while boxes with the z-dimension of approximately 12 nm were used for other systems with weaker binding affinities to reduce the computational cost. Because binding of charged molecules to zwitterionic PC membranes substantially increases the lamellar repeat distance in a lipid bilayer stack,[67] we consider such a large interbilayer space reasonable for the studies of strongly membrane-bound charged molecules.

The combination of the required large simulation boxes with the slow equilibration times makes simulations unfeasible for some strongly bound molecules. For example, bulk concentrations of the SMS peptide are in the sub-millimolar range in experiments[62] which would require approximately 200,000 water molecules per peptide in simulations. On the other hand, binding or unbinding events were not observed for dibucaine during the 1360 ns simulation. Nevertheless, for etidocaine and TPP, the binding dynamic was sufficiently fast to enable to equilibrate the simulations within feasible time (1–5 μs) and simulation box z-dimension (12–80 nm) scales (Figures S4 and S5 in the Supporting Information). Therefore, we proceeded to evaluate the binding affinities of these molecules in MD simulations to membranes against the available lipid head-group order parameter data.

Comparison of Etidocaine and TPP Binding to POPC Membranes between Simulations and Experiments

The binding behavior of etidocaine and TPP molecules to POPC membranes predicted by different force field parameters is illustrated in Figure a using the density profiles along the membrane normal and PMF curves derived from therein. For etidocaine, CHARMM36-SwissParam predicts the weakest binding affinity, followed by the CHARMM36-ParamChem and CHARMM36-Paramchem ECC. The weakest affinity for TPP ions is predicted by CHARMM36-QMcharges, followed by CHARMM36-ProteinFF and CHARMM36-QMcharges ECC.

Figure 2

(a) Number densities and PMF profiles of etidocaine and TPP along the membrane normal from simulations with the small molecule concentration of 140 mM at 298 K. To emphasize the membrane region, the x-axis is limited between −6 and +6 nm, for the full density profile of etidocaine, see Figure S8 in the Supporting Information. Etidocaine results are on the left and TPP on the right throughout the figure. (b) Changes in the α and the β order parameters as a function of the bound etidocaine and TPP per lipid in the POPC membrane. (c) Changes in the POPC head-group α-carbon order parameters as a function of small molecule concentration with respect to water from simulations and experiments for etidocaine[64] and TPP.[63] Error bars from simulations are not visible for most order parameters as they are smaller than the point size.

To evaluate the simulation predictions against experiments, we calculated the C–H bond order parameters of α and β segments in the POPC head group, which are proportional to the amount of bound charge in a membrane.[34] As shown in Figure b, a linear decrease in these order parameters is observed as a function of the number of bound etidocaine or TPP molecules to a POPC membrane. This trend is observed in both MD simulations and NMR experiments,[63,66] but the slope of the decrease depends on force field parameters used in a simulation and deviates from experiments in all cases except for the α carbon in ECC simulations. This indicates that the interactions of these molecules with a POPC lipid membrane are not accurately captured by any of these simulations.

Nevertheless, we consider that the responses of the α-carbon order parameters are sufficiently close to experiments to be used for the validation of etidocaine and TPP binding affinities in simulations, yet the observed inaccuracy should be taken into account when interpreting the results. The decrease in the α-carbon order parameter as a function of etidocaine concentration in water in Figure c correlates with the binding affinity in Figure a, being weakest for the CHARMM36-SwissParam, followed by the CHARMM36-ParamChem and CHARMM36-ParamChem ECC. Comparison with the experimental data suggests that CHARMM36-ParamChem simulations predict the etidocaine binding affinity closest to experiments, while CHARMM36-SwissParam predicts too weak and CHARMM36-ParamChem ECC too strong binding. A similar comparison for TPP suggests that CHARMM36-ProteinFF prediction is closest to experiments, while CHARMM36-QMcharges predicts too weak binding, which is then overestimated after applying the ECC. However, because QM-derived partial charges are presumably more realistic, the better results with CHARMM36-ProteinFF parameters probably originate from the cancellation of errors.

In conclusion, none of the parameters generated with the standard approaches for etidocaine or TPP predicted the overbinding of these molecules to POPC membranes, while some parameters predicted too weak binding affinity. Therefore, our results suggest that the previously observed overbinding of sodium and calcium ions to membranes in canonical MD simulation force fields[16] is not a general feature for all positively charged biomolecules.

Effect of Electronic Polarizability and Binding Mechanism of Etidocaine and TPP to Membranes

The most likely source for the discrepancies in MD simulations of charged small-molecule binding on membranes is the missing electronic polarizability. While polarizable force fields are available, such as CHARMM36-Drude,[68] they cannot yet capture the lipid head-group conformational ensembles with sufficient accuracy to enable the usage of head-group order parameters for the evaluation of small-molecule binding affinities.[69] On the other hand, the implicit inclusion of electronic polarizability by scaling the partial charges of atoms in an approach known as ECC[48] has been shown to correct the overestimated calcium binding to membranes containing charged and zwitterionic lipids.[18−20] Therefore, we decided to study the effect of electronic polarizability by applying the ECC to CHARMM36-ParamChem parameters of etidocaine and to CHARMM36-QMcharges of TPP. These were used with the recently introduced parameters for POPC where ECC was applied to the CHARMM36 lipid force field.[47]

Applying ECC increases the binding affinities of both the etidocaine and the TPP molecules to POPC bilayers, leading to an intensified decrease in order parameters in Figure a. On the other hand, ECC makes the β-carbon order parameter less sensitive, and α-carbon slightly more sensitive, to both etidocaine and TPP. This brings the responses to the number of bound molecules in Figure b closer to experiments, although the decrease in β-carbon upon the addition of etidocaine is now slightly underestimated. In conclusion, ECC potentially improves the interactions of etidocaine and TPP with the POPC head group. However, it introduces too strong binding affinity, emphasizing the need for further optimization of force field parameters to accurately capture the binding details of TPP, etidocaine, and other charged small molecules to lipid membranes.

Notably, the effect of ECC on the binding affinities of etidocaine and TPP is opposite to its effect on the calcium-binding affinity which decreased upon applying ECC.[18−20] This can be explained by the different binding mechanisms of calcium and small molecules, such as etidocaine and TPP. The main driving forces for calcium binding to a membrane with PC lipids are the direct electrostatic interactions with phosphate oxygens,[18] as also seen in Table , where the majority of calcium–lipid interactions occur with phosphate oxygens. For small molecules with charges surrounded by hydrophobic moieties, such as etidocaine and TPP, the hydrophobic effect is the most likely driving force for their membrane binding. In the case of calcium, the inclusion of ECC reduces the electrostatic attraction with phosphate oxygen, which is the most likely reason for the reduced binding affinity. Indeed, the reduction in interactions with phosphate oxygen due to ECC is observed for calcium in Table  but not for etidocaine or TPP. On the other hand, the scaling of total charge by ECC makes etidocaine and TPP less soluble in water, thereby increasing their binding affinity to membranes. In conclusion, the binding behavior of monoatomic ions, such as calcium, to membranes differs from small molecules, such as etidocaine and TPP, due to their different binding mechanisms.

Table 1

Number of Contacts with Any Lipid Atom per Bound Molecule and Percentage of the Contacts with Phosphate Oxygensa

molecule	contacts/bound molecule	P-contacts/all contacts [%]
Ca2+ CHARMM36 no NBFIX	3.9	90.6
Ca2+ CHARMM36 ECC	2.2	64.6
Ca2+ Lipid14	7.6	67.3
Ca2+ Lipid14 ECC	3.5	64.9
ETI CHARMM36-ParamChem	17.4	10.0
ETI CHARMM36-ParamChem ECC	18.5	11.1
TPP CHARMM36-QMcharges	17.0	7.9
TPP CHARMM36-QMcharges ECC	18.74	11.2

CHARMM36-based simulations with 450 mM CaCl2[47] available at refs (70) and (71) and Amber-based Lipid14 simulations[18] with 467 mM CaCl2 available at refs (72) and (73) were used. Results from simulations with 70 mM Etidocaine and 140 mM TPP are shown.

Comparison of Sodium, Calcium, Etidocaine, and TPP Binding to Membranes

Experimental methods can provide accurate information on binding affinities of charged molecules to membranes, but it is not straightforward to define a concentration and model-independent binding coefficient that would correctly describe the binding affinity of charged biomolecules to membranes with a single number, as demonstrated in Figure . While the relative binding affinities can be judged by examining the whole experimental binding isotherm, detailed information such as the depth of binding or the potential mechanism of penetration remains inaccessible from the experimental data alone. On the other hand, simulations that reproduce the experimental data can be used to interpret such properties from the experimental data.

Because a single force field that would correctly reproduce the binding of all charged molecules to membranes is not yet available, we selected the best available models from this and our previous study[47] to compare the binding of sodium, calcium, etidocaine, and TPP to a POPC membrane. The density profiles, PMFs, and α-carbon order parameter decrease compared with experiments for these simulations are shown in Figure . In these simulations, the sodium exhibits weaker binding affinity than other molecules as expected for simple monovalent ions.[16] Calcium and TPP bind with similar affinities, while etidocaine exhibits the strongest binding. These relative binding affinities are in line with the experimentally determined number of bound particles in Figure and measured order parameters but not with the reported binding coefficients. Despite the similar binding affinities, TPP penetrates deeper in the membrane than calcium. This can be explained by the different binding mechanisms of these molecules. Calcium ions bind to phosphate oxygens in the POPC head group, whereas the binding of TPP is driven by the hydrophobic effect. Similarly to TPP, also etidocaine penetrates deeper into a membrane. The PMF profiles in Figure show a lower energy barrier at the membrane center for etidocaine and TPP than for water, sodium, and calcium. However, we did not observe any permeation events of etidocaine or TPP through the membrane. This is in contrast to water molecules, for which dozens of events were observed in each simulation. This can be explained by the larger size of etidocaine and TPP molecules increasing their probability to locate at the membrane center without permeating through the membrane. While we consider these observations as reasonable interpretations of the experimental data with the best available MD simulation models, we emphasize that the detailed interactions in these simulations may not be exactly correct since the β-carbon order parameter response to studied charged molecules was not in agreement with the experiment in Figure .

Figure 3

(a) Molecule number density profiles and PMF profile from simulations giving the most realistic binding affinities of etidocaine (CHARMM36-ParamChem), TPP (CHARMM36-ProteinFF), calcium, and sodium (simulations with ECC correction[47]) according to the head-group order parameters responses. Results from CHARMM36 with ECC applied[47] are shown for calcium and sodium,[70,74] but simulation with ECC applied to Amber parameters would give essentially the same conclusions.[18] (b) Experimental NMR α-carbon order parameter responses to etidocaine, TPP, calcium, and sodium ions are shown as full lines together with selected most realistic simulations shown as dashed lines.

In conclusion, our results indicate that the binding affinities of charged drugs, etidocaine, dibucaine, and SMS, are significantly stronger to zwitterionic membranes than the binding affinities of sodium and calcium. Also, their binding mechanisms differ from that of monoatomic ions. Therefore, the binding of these drugs to membranes would be most likely not interfered by physiological ions. On the other hand, the binding of TPP could be interfered by calcium as their affinities are similar, yet their interrelated binding may be complicated as their binding mechanisms are different. MD simulations have the potential to model such complex behavior, but force field parameters correctly capturing both TPP and calcium-binding are not yet available. Furthermore, improvements in force fields would be needed to correctly model the binding of etidocaine to POPC head groups as the experimental response of β-carbon order parameters to the number of bound molecules is not captured by any simulation models in Figure b.

Conclusions

While concentration and model-independent binding coefficients that would enable the comparison between the binding affinity of simple ions and charged small molecules to membranes are not available, we managed to evaluate the binding affinity of etidocaine and TPP to POPC lipid bilayers in MD simulations against experiments using the changes in lipid head-group-order parameters as done previously for sodium and calcium.[16] However, the required size and length of simulations increase with increasing binding affinity, thereby setting practical limitations for the molecules that can be studied with this approach. For etidocaine and TPP, the simulation box sizes of 12–81 nm in membrane normal direction and time scales of 1–5 μs were sufficient, but the evaluation of SMS or dibucaine binding with higher affinities was not feasible within this work.

According to the evaluation based on lipid head-group order parameter changes, etidocaine and TPP binding affinities to a POPC membrane were slightly underestimated or close to experiments in simulations with standard CHARMM36-based force field parameters. This is in contrast with previous studies for sodium and calcium ions, where binding affinities were typically overestimated by canonical force fields.[16] Therefore, our results suggest that force fields do not generally overestimate the binding of all positively charged molecules to membranes. Furthermore, the implicit inclusion of electronic polarizability using ECC increased the binding affinities of etidocaine and TPP to a POPC membrane, whereas calcium is known to behave oppositely.[18−20] These observations can be explained by the different binding mechanisms of calcium and the small molecules with hydrophobic moieties. Calcium binds to the lipid phosphate oxygens via direct electrostatic interactions. The binding of small molecules is, on the other hand, driven by hydrophobic interactions. While ECC reduces the direct attraction between lipids and ions, it also reduces the solubility to water, which is more important for small molecules. Different binding mechanisms also explain the deeper penetration of etidocaine and TPP into the membrane core than calcium ions.

The relatively weak binding of metal cations (also other than sodium and calcium[75]) with a distinct mechanism to PC membranes suggests that they probably do not interfere with the binding of charged drugs with higher affinities, such as etidocaine, dibucaine, and SMS. The situation may be more complex for molecules with similar binding affinities but different mechanisms, such as calcium and TPP. While MD simulations are a promising tool to model such complicated systems, the accuracy of current force fields is not sufficient for such applications as none of the available force fields correctly captures both the calcium and charged small-molecule binding to membranes.

The evaluation of charged small-molecule binding affinity to membranes using changes in the lipid head-group-order parameters offers a tool to evaluate and develop force fields that would correctly capture interactions between charged biomolecules and membranes. This approach is complementary to the comparisons of partition coefficients[13,15] as it gives information also on the detailed interactions between lipids and small molecules in addition to the binding affinity. While straightforward application of ECC to standard force fields has substantially improved the ion binding behavior,[18,19] the small molecules seem to require further optimization. Nevertheless, we believe that our results pave the way toward force fields that would correctly capture lipid membrane interactions with charged biomolecules and amino acids. Such force fields have potential applications in a wide range of fields, from drug design to molecular biology and toxicology.