Digging into Lipid Membrane Permeation for Cardiac Ion Channel Blocker d-Sotalol with All-Atom Simulations

Interactions of drug molecules with lipid membranes play crucial role in their accessibility of cellular targets and can be an important predictor of their therapeutic and safety profiles. Very little is known about spatial localization of various drugs in the lipid bilayers, their active form (ionization state) or translocation rates and therefore potency to bind to different sites in membrane proteins. All-atom molecular simulations may help to map drug partitioning kinetics and thermodynamics, thus providing in-depth assessment of drug lipophilicity. As a proof of principle, we evaluated extensively lipid membrane partitioning of d-sotalol, well-known blocker of a cardiac potassium channel Kv11.1 encoded by the hERG gene, with reported substantial proclivity for arrhythmogenesis. We developed the positively charged (cationic) and neutral d-sotalol models, compatible with the biomolecular CHARMM force field, and subjected them to all-atom molecular dynamics (MD) simulations of drug partitioning through hydrated lipid membranes, aiming to elucidate thermodynamics and kinetics of their translocation and thus putative propensities for hydrophobic and aqueous hERG access. We found that only a neutral form of d-sotalol accumulates in the membrane interior and can move across the bilayer within millisecond time scale, and can be relevant to a lipophilic channel access. The computed water-membrane partitioning coefficient for this form is in good agreement with experiment. There is a large energetic barrier for a cationic form of the drug, dominant in water, to cross the membrane, resulting in slow membrane translocation kinetics. However, this form of the drug can be important for an aqueous access pathway through the intracellular gate of hERG. This route will likely occur after a neutral form of a drug crosses the membrane and subsequently re-protonates. Our study serves to demonstrate a first step toward a framework for multi-scale in silico safety pharmacology, and identifies some of the challenges that lie therein.


Sotalol solvation and protonation energetics based on implicit solvent QM calculations
We used quantum mechanical (QM) calculations in vacuum as well as in widely used implicit solvent model SMD (Marenich et al., 2009) with results shown in Supporting Tables S1 and S2. According to MP2/6-31G(d) QM calculations in vacuum SotN form is favored over SotZ by around 47 kcal/mol. This preference drops to ~30 kcal/mol in hexadecane, a membrane-mimetic alkane and to just ~1.3 kcal/mol in water due to substantially more favorable solvation free energy for SotZ. However, even in water the neutral species is slightly favored over the charged one (Table S1). We also tested the robustness of such assessment using HF, MP2 and a popular hybrid density functional method B3LYP levels of theory and two basis sets: standard 6-31G(d) and substantially more complex aug-cc-pVDZ, as shown in Supporting Table 2. In most cases, SotN is still predicted to be favored in water with the only exception of B3LYP/ aug-cc-pVDZ calculations predicting nearly equal stabilities of SotZ and SotN (Table  S2). These calculations do not guarantee accurate assessment since they do not take into account specific solute-solvent interactions explicitly, but nevertheless provide us with a valuable approximation in the absence of experimental data.
Comparison of aqueous and hexadecane solvation free energies obtained via QM calculations helped us predict energetic penalties for different sotalol ionization forms as they transfer from water to a non-polar hydrophobic environment such as that of membrane interior ( Table S1). All the ionized drug species, i.e. SotC, SotZ and SotA, are expected to have very large water-hydrocarbon translocation penalties, exceeding 30 kcal/mol, whereas that for SotN is predicted to be ~8 kcal/mol. Interestingly, SotZ is predicted to have a larger energetic penalty than either SotA or SotC despite the fact that the latter two forms carry a net charge. The high energetic cost for SotZ is commensurate with its huge dipole moment, which is almost twice the magnitude of the dipoles for cationic or anionic forms. All dipole moments were computed with the same drug orientation and position for all of them.

Accuracy of P21 point group (symmetry) application for CHARMM MD simulations of drug -lipid membrane systems
Two US simulations for each drug form (SotC and SotN) were performed -using NAMD and CHARMM programs with the same empirical force field models for all the system components and practically the same simulation parameters. A major difference is that in CHARMM simulations we used P21 symmetry, which allows sampling of lipid redistribution between bilayer leaflets (Dolan et al., 2002) to maintain the same cross-sectional area of the membrane as the drug moves across it. In the absence of P21 symmetry an area difference may arise between the top and bottom leaflets over the course of drug translocation. A strain contribution that stems from the area difference can increase free energy barrier, and P21 symmetry might help alleviate this problem. Indeed, we observed some variation in the average number of lipids in the top and bottom leaflets, which is small when a drug is in bulk water and can reach up to 4 lipids, when it is in the membrane interior (see top panels in Figure S3). We also tested if P21 symmetry can lead to additional membrane distortions and thus measured average area per lipid and membrane thickness (middle and bottom panels in Figure S3), both of which are in good agreement with experimental values at 303 K, 64.3±1.3 Å 2 and 39.1±0.78 Å, respectively (Kucerka et al., 2011). Deuterium order parameter profiles, |SCD|, shown in Figure S4, are also in good agreement with experimental estimates (Ferreira et al., 2013).

Drug reorientation sampling
We found that hindered drug tumbling discussed in the main text can lead to substantial asymmetries in the PMF profiles, especially for SotC. In our original NAMD and CHARMM US MD simulations we obtained asymmetries as high as 4 and 12 kcal/mol ( Figure S9), which are clearly unacceptable for any simulation analysis. We found that such asymmetries are mostly caused by SotC being stuck in a particular orientation with respect to membrane normal and perturbing the membrane by forming an interfacial connection with water and lipid head groups. Since the cationic ammonium group is a main attractor of polar moieties in SotC, the interfacial connection forms with the side of the membrane the functional group is pointing toward. Due to a random initial orientation of drug molecules and their slow reorientation in the membrane interior, drug molecules may form interfacial connections with an opposite bilayer leaflet. For instance, a drug with COM at z = 2 Å can coordinate water molecules and lipid head groups from the lower leaflet. Such interfacial connection will likely lead to high strain energy due to substantial membrane perturbations required. Thus, we need to balance this out with connections to an adjacent bilayer leaflet for the same drug position. This can be done in several ways as was described in our previous study on arginine side chain membrane partitioning (Li et al., 2008). Here for CHARMM simulations, we created additional US windows for |z| £3 Å by starting with the neighboring windows with interfacial connection to an opposite bilayer leaflet. This led to a reduction in asymmetry from ~12 to ~3 kcal/mol ( Figure S9). For SotC NAMD simulations a different approach was used where a second set of US windows for |z| £20 Å (roughly corresponding to POPC membrane thickness) was generated with a different initial drug orientation, namely rotated around z-axis by 180°. Using two sets of simulations for central US windows also helped to substantially reduce asymmetry from ~4 kcal/mol to nearly 0 kcal/mol in this case (see Figure S9). Interestingly, such approach worked differently for SotN NAMD simulation leading to an increase in asymmetry in bulk solvent by ~1 kcal/mol but a similar decrease in the interfacial binding region ( Figure S9). Free energy profiles shown in Figure 7A in the main text included this additional sampling except for SotN CHARMM run, where such simulations were not performed. We also tested slightly different approach, were we flipped initial drug orientation around the xy plane by 180° for central windows (-20≤z≤20 Å for SotC and -4≤z≤4 Å for SotN). PMF convergence plots with those simulations included are shown in Figure S10. There is ~2 kcal/mol asymmetry for SotC and practically no asymmetry for SotN simulations. Final symmetrized profiles from those runs (from 4 to 15 ns) are nearly identical to ones shown in main text Figure (1991). Determination of partition coefficients between dimyristoylphosphatidylcholine and water using differential scanning calorimetry. Analytica chimica acta 251, 79-81.       Non-symmetrized profiles are shown. They were computed using weighted histogram analysis method (WHAM) either using all the simulation data or discarding first 4 ns for each umbrella sampling window based on solvation number profiles in Figure S5. In some cases, additional simulations with initial drug orientation rotated around the z axis (indicated as +flip) were performed to improve drug re-orientation sampling.

SotN (NAMD)
CisC ( Both non-symmetrized (left) and symmetrized (right) profiles are shown. They were computed using weighted histogram analysis method (WHAM) discarding first 4 ns for each umbrella sampling window based on solvation number profiles in Figure S5. Additional simulations with initial drug orientation flipped around the xy plane were included for central windows (-20≤z≤20 Å for SotC and -4≤z≤4 Å for SotN) were performed to test drug re-orientation sampling. Final symmetrized profiles from those runs (from 4 to 15 ns) are nearly identical to ones shown in main text Figure  7A.