pMD-Membrane: A Method for Ligand Binding Site Identification in Membrane-Bound Proteins

Probe-based or mixed solvent molecular dynamics simulation is a useful approach for the identification and characterization of druggable sites in drug targets. However, thus far the method has been applied only to soluble proteins. A major reason for this is the potential effect of the probe molecules on membrane structure. We have developed a technique to overcome this limitation that entails modification of force field parameters to reduce a few pairwise non-bonded interactions between selected atoms of the probe molecules and bilayer lipids. We used the resulting technique, termed pMD-membrane, to identify allosteric ligand binding sites on the G12D and G13D oncogenic mutants of the K-Ras protein bound to a negatively charged lipid bilayer. In addition, we show that differences in probe occupancy can be used to quantify changes in the accessibility of druggable sites due to conformational changes induced by membrane binding or mutation.

section_content: Identification of a suitable ligand-binding site on a drug target is a crucial first step in structure-based computer aided drug discovery [1]. This is not a trivial task if the desired target site is an allosteric one that is not readily observable in average experimental structures [2]. 

Recently, a number of techniques have been developed that allow for the identification of (allosteric) ligand binding sites in target proteins [3] [4] [5] [6]. Because ligand binding site identification usually requires sampling of the target's configurational space, considerable effort has also been made toward integrating molecular dynamics (MD) simulation into the site identification process (e. g. [6] ). In particular, MD-based computational solvent mapping [7] [8] [9] [10] [11] [12] is attracting wide attention as a convenient means of binding site identification in dynamic targets. Interest in this approach will likely increase with the expanding scope of MD simulations and because it recapitulates multi-solvent crystallographic [3] and fragment-based NMR screening experiments [5]. 

A typical MD-based computational solvent mapping entails carrying out MD simulations in the presence of small organic molecules in the solvent (e. g. [7, 8] ). The goal is to use the small organic molecules as probes to search for binding sites on an ensemble of MD-sampled target structures. The probability of contact (or interaction) between probe and protein atoms is then used to evaluate the druggability of sites. The method has been described in a number of recent reports under various names: probe-based MD [8], mixed-solvent MD [12], solvent competition [7], co-solvent MD [10] and ligand competitive saturation [9, 13]. We use the term probe-based MD (pMD) throughout this report. Surprisingly, thus far pMD has been applied only to soluble proteins despite the fact that some of the most important drug targets require membrane binding for their biological activity [14] [15] [16] [17] [18] [19]. A major goal of the current work is to extend the applicability of pMD to membrane-bound drug targets. This requires mitigating possible effects of the probe molecules on membrane structure and dynamics. For example, we previously found that small organic molecules such as ibuprofen, indomethacin and cholic acid partition into the hydrophobic core of DPC micelles [20] [21] [22]. Others found that similar small organic molecules partition into bilayers [23, 24]. Here we describe pMD-membrane, a method that avoids membrane partitioning of probe molecules and enables allosteric ligand binding site identification in proteins bound to a bilayer surface. 

Another challenge in current efforts of computational binding site identification is the difficulty in discriminating between closely related homologs or mutations that are associated with different disease phenotypes. Whether pMD can capture small changes in the properties of binding sites due to conformational changes induced by membrane-/substrate-binding or mutation has not been examined. We introduce analysis techniques to evaluate differential probe occupancy that inform on the changes in potential druggability of a site. 

We tested pMD-membrane and the new analysis tools on G12D and G13D mutants of K-Ras. We chose these K-Ras mutants as model systems for a number of reasons. First, K-Ras is a prototypical example of membrane-associated small GTPases for which there exist abundant experimental structure data [25]. Secondly, we recently found that the interaction of K-Ras with membrane involves at least two distinct conformations (Prakash and Gorfe, unpublished results). Third, K-Ras is a key regulator of numerous signaling pathways involved in cell division and proliferation [25] [26] [27], and therefore it is physiologically and therapeutically highly relevant. In fact, 15-25% of all cancer cases are associated with mutations in the homologous K-, N-and H-Ras proteins [28] ; K-Ras mutations represent 85% of these [29]. 

Previous efforts to inhibiting aberrant Ras function have failed [30, 31], but a number of allosteric Ras ligands have been discovered recently [32] [33] [34] [35] [36] [37] [38]. While these ligands are promising starting points, none have the necessary potency and selectivity to become a lead compound. Therefore, the search for Ras inhibitors continues. Desirable properties of a potential Ras inhibitor may include the following: (i) Ability to directly bind to membrane-associated cellular Ras. Inhibitor activity in solution is not sufficient because membrane binding is essential for the biological function of Ras, and there is evidence that signaling specificity among isoforms may involve distinct membrane localization and therefore differential accessibility to effectors and modulators [39] [40] [41] [42]. It is therefore important that changes in conformation and dynamics upon membrane binding are explicitly considered in binding site identification / drug discovery efforts against Ras proteins. (ii) Specificity toward a given Ras isoform. This is because, as we alluded to above, different Ras isoforms are associated with distinct cancer types [29] despite the fact that they share a catalytic domain that is nearly identical in sequence and average structure [17, 43]. For instance, K-Ras mutations are prevalent in lung, colorectal and pancreatic carcinomas [44] [45] [46] [47], N-Ras mutations in melanomas and hematologic malignancies [48] [49] [50], and H-Ras mutations in bladder and thyroid cancers [51, 52]. (iii) Specificity toward a mutation. This is because different Ras mutants, such as G12D and G13D, are associated with different cancer types [25, 29, 53]. 

The ability to identify unique ligand binding sites on each Ras isoform or mutant is the first step toward addressing the issues listed above. We demonstrate that pMD-membrane helps achieve this goal, and illustrate its robustness and sensitivity using two of the most common oncogenic mutants of K-Ras: G12D and G13D. We also show that differential membrane binding leads to altered propensity of probes for binding sites.
section_name: Introduction

section_content: The theoretical basis of pMD can be found in previous reports [7] [8] [9] [10], as well as its application to search for allosteric ligand binding sites on K-Ras in solution [11]. Here we focus on technical progresses that allow the use of pMD simulation to predict ligand-binding sites in membrane-bound lipid-modified proteins. We also describe post-simulation analysis techniques to correlate druggability of sites with changes in conformational dynamics induced by mutation or membrane reorientation.
section_name: Methods

section_content: Previous simulation [54] and experimental studies [55, 56] have shown that the full-length H-Ras protein interacts with membrane in a non-random fashion; it adopts two distinct orientations with respect to the membrane plane depending on the bound nucleotide [57]. In a separate work (Prakash and Gorfe, to be published), we examined the bilayer interaction of GTPbound G12D K-Ras based on a total of ~7. 5 μs all-atom MD simulations. The analysis yielded two predominant modes of membrane binding that differ in the membrane orientation of the catalytic domain (Fig 1 ). Because these orientation differences can potentially be exploited for the development of isoform-and mutation-specific small molecule inhibitors, we used the two distinct conformations as the starting structure for the current pMD-membrane analysis of G12D and G13D K-Ras, as described below.
section_name: Generation of initial configurations of membrane-bound mutant K-Ras

section_content: We performed two sets of simulations: reference set and target set. The reference set involved two 60 ns pMD-membrane simulations of G12D K-Ras starting from the two conformations shown in Fig 1. These simulations were performed without any modification to the force field parameters. The target set involved three 100 ns-long pMD-membrane simulations with the Lennard-Jones non-bonded interaction between selected atoms of the probe and lipid molecules modified as described in the following section. Two of the target simulations were on G12D K-Ras while the third was on G13D K-Ras. One of the G12D simulations was started from a conformation in which helices 3 and 4 directly interact with the bilayer (Fig 1A ) Since cellular K-Ras interacts with negatively charged inner leaflet of the plasma membrane via its polybasic and farnesylated C-terminus, we built a heterogeneous lipid bilayer made up of 320 POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) and 96 POPS (1-palmitoyl-2-oleoyl-sn-glycero-3-phoserine) lipids. The resulting symmetric bilayer was equilibrated through a 250 ns production simulation. We then embedded full-length G12D or G13D K-Ras4B in one leaflet of the pre-equilibrated bilayer. Membrane insertion was guided by previous reports [58] to determine the insertion depth of the farnesyl tail into the bilayer core. The resulting system was placed in a 114 x 112 x 110 Å 3 box containing 26299-27573 TIP3P water molecules and 1337-1423 isopropanol probe molecules. In each system, we maintained a 20:1 water to probe ratio and neutral charge by adding 96 sodium ions. Additional sodium and chloride ions were added to mimic physiological ionic strength. The total number of atoms varied between 153316 and 158178 depending on the initial conformation of the protein. 

Following system construction and 3000 steps of conjugate gradient energy minimization, we used simulated annealing to homogenize the probe and water molecules around the protein and the bilayer. The annealing process involved the application of a harmonic restraint with a force constant of 4 kcal/mol/Å 2 on the protein and lipid heavy atoms to prevent protein unfolding and bilayer instability, and incrementing the temperature every 5000 steps by 50 K until a temperature of 650 K was reached, followed by gradual cooling by 10 K every 5000 steps until a final temperature of 310 K was achieved. The resulting system was equilibrated for 1 ns while gradually decreasing the restraint force constant to zero. This was followed by a production run of 60 ns for the reference simulations and 100 ns for the target simulations. A nonbonded cutoff of 12 Å was used during both the equilibration and production phases of each simulation. Long-range electrostatic interactions were calculated by the Particle Mesh Ewald (PME) method [59] with SHAKE [60] restraints applied on bonds involving hydrogen atoms. Simulations were performed with a 2 fs timestep in the NPT ensemble (constant number of particle N, temperature T = 310 K, and pressure P = 1 bar). Nose-Hoover Langevin piston for pressure control was used to maintain constant pressure and Langevin thermostat to maintain constant temperature. Short-range non-bonded forces were computed every timestep and long-range electrostatic forces every other step. All simulations were performed with the NAMD2. 9 program [61] using the CHARMM27 force field for proteins [62] and CHARMM36 for lipids [63] ; isopropanol was parameterized as described in ref [11].
section_name: Molecular dynamics simulation of K-Ras mutants in mixed bilayer and mixed solvent

section_content: pMD-membrane simulation with the unmodified CHARMM force field led to partitioning of a fraction of the probe molecules into the bilayer (see Fig 2 and Results and Discussion). To prevent this diffusion of probe molecules into the bilayer interior, we modified the vdW interaction between the central carbon atom of the probe and the CTL2 atom type of the POPC and POPS lipids. After several tests (see for example S1 Fig), we arrived at the following protocol: the well depth of the Lennard-Jones potential was reduced to a very small value of 0. 01 kcal/ mol (see ref. 9), and the minimum inter-particle distance among the selected atoms was increased to 7 Å (S1 Fig). The former ensures that attraction between the apolar atoms of the lipids and probe molecules is almost completely eliminated, while the 7 Å distance yields a reasonable balance between allowing the probe to access the polar bilayer surface where the protein sits and preventing it from penetrating into the hydrophobic core (S1 Fig). This modification was made utilizing the NBFIX correction term in the CHARMM force field. This approach is similar in principle to the modifications MacKerell and colleagues made via dummy atoms to prevent aggregation of probe molecules [9].
section_name: Modification of selected non-bonded interaction terms

section_content: We calculated the extent of protein-probe interactions using both distance-based and gridbased measures of probe occupancy. 

(i) Distance-based probe occupancy. The distance-based probe occupancy measure quantifies the probability that a probe molecule is within a certain distance of a given protein heavy atom. Our implementation of this procedure entailed computing the probe occupancy R x i at protein atom i in simulation x as

where N frames is the number of frames in the trajectory, N probes is the number of probe molecules in the system, d ijk is the minimum distance between protein atom i and any heavy atom of probe molecule j at frame k. n (an integer) and α (a small positive real number) are adjustable parameters carefully chosen to include all probes within a desired cutoff value and then smoothly switch off the selection within a small distance range. In this study we used n = 20 and α = 1×10 -17 so that all probes within ~6 Å are fully included and those between about 6 and 8 Å are included in a weighted fashion as shown in S2 Fig. We used an arbitrary scaling factor of 250 to facilitate visualization. To compare changes in the probe-binding potential of a site in two different simulations, we calculated difference probe occupancy (ΔR i ) defined as the probe occupancy of atom i in trajectory x (R x i ) and the probe occupancy of the same atom in trajectory y (R y i ):

(ii) Grid-based probe occupancy. The Grid-based probe occupancy calculation technique has been extensively discussed in previous reports by others [7, 8], and has been used in our previous work on K-Ras in solution [11]. In the current work, we calculated grid occupancy in two slightly different ways. In the first approach, only probe molecules that are within a cutoff distance from the protein were included for grid occupancy calculation. We used a cutoff value of 4 Å between any heavy atom of the probe molecule and any heavy atom of the protein. After aligning the trajectory frames based on backbone atoms excluding the flexible switch regions and the termini, we calculated time-averaged number densities per 1 Å 3 grid using the Volmap and Volutil plugins of VMD [64]. The grid number densities were converted to grid free energies as described previously [11]. 

In the second approach, we first removed translational and rotational motions by aligning the trajectory to the initial frame of each production run using backbone atoms of the catalytic domain (residues 1-166). Then, a spatial concentration profile of the probe molecules was evaluated on a hexahedral mesh using the following formula

where C 0 is the standard state concentration (1M = (1 mole/liter = 1/1660. 3 molecule/Å 3 ), δ is the Kronecker delta function, bc is the integer floor function, c ijk is probe concentration in M at node (i,j,k), N heavy = 4 is number of heavy atoms in the probe molecule, dx, dy and dz are the hexahedral grid spacing in the x, y and z direction, respectively. We used a uniform grid spacing of 0. 37 Å in the spatial concentration profile calculations. 

The probe concentration value at each node was smoothed out by assigning an average value based on its neighboring node values,

This averaging procedure eliminates the noise in the constructed spatial concentration profile and assists in making a more vivid iso-surface visualization.
section_name: Probe occupancy calculations

section_content: To assess the relative polarity of putative ligand binding sites, we calculated the ensemble-averaged orientation of probe molecules with respect to the surface of the protein. Specifically, we calculated <cos(θ r )>, where θ r is defined as the angle between the vector radiating from the central carbon atom to the O atom of an isopropanol probe molecule and the vector connecting the center of the catalytic domain of the protein (residues 1-166) and a given atom on the protein surface. Only protein atoms whose probe occupancy was above an empirically determined threshold were considered for this analysis (in this work we used R i ! 0. 24).
section_name: Analysis of probe orientation

section_content: Convergence of protein-probe interactions was evaluated by monitoring the time evolution of the atom-averaged probe occupancy R ave R ave ¼

where the summation is over all protein atoms M with R > 0. 05. The 0. 05 cutoff ensures that approximately all of the atoms with non-zero R-value are included in the statistics. (Note that the profile of R ave would be unaffected even if atoms with zero R-values were also included. ) S3 Fig (left) shows that R ave has equilibrated within the first ~20 ns in all three of our target simulations. Moreover, comparison of the mean (<R ave >) in 20 ns blocks yielded very small differences among all blocks except the first one. For instance we obtained <R ave >±S. D of 0. 20 ±0. 05, 0. 21±0. 02, 0. 22±0. 04, 0. 24±0. 02 and 0. 24±0. 03 for the 1-20, 20-40, 40-60, 60-80 and 80-100 ns blocks of the G12D simulation in membrane binding mode 1; similar results were obtained for G12D mode 2 and G13D simulations. This suggests that the simulations were well equilibrated in terms of probe binding, and that any portion of the last 80 ns data can be used to compute average occupancies. For a better statistics, however, we used all frames in the last 80 ns of the trajectories unless stated otherwise. Statistical uncertainty in R i was estimated based on an analysis of block standard error (BSE) following the procedure described by Grossfield and Zuckerman [65]. In this analysis, BSE was calculated at different sizes of time blocks (b n ). The convergence of BSE versus b n was used to evaluate the quality of our sampling and to determine the maximum value of the BSE, which serves as a measure of our sampling error. In S3 Fig (right), we show several BSE vs. b n plots for a few atoms chosen for illustration of the diverse convergence rates and error values. For a more rigorous analysis of errors in the relative probe accessibility of each atom in trajectories x and y, i. e., the uncertainty in ΔR i, we first calculated BSE in R x i and R y i for multiple block sizes b n. We then used nonlinear fitting of the dependence of BSE on b n to the following functional form: a 0 (1-a 1 exp(-a 2 b n )), where a 0, a 1 and a 2 are three fitting parameters with a 0 being an estimate of the asymptotic value of BSE. We then computed the uncertainty in ΔR i by combining the BSEs in R x i and R y i (Eq 6).
section_name: Convergence and estimation of sampling error

section_content: As noted previously, probe-based MD has been successfully applied to a number of soluble drug targets [8] [9] [10] 12], including the soluble catalytic domain of K-Ras [11]. Yet in cells K-Ras and a large number of other targets are membrane-bound. Therefore, we first describe the challenge of using probes in membrane simulations and how to overcome the challenge. We will then discuss the application of pMD-membrane on G12D and G13D K-Ras proteins in bilayer.
section_name: Results/Discussion

section_content: To our knowledge pMD has not been previously applied to membrane proteins. Because molecular probes are completely or mostly nonpolar and small, we reasoned that one of the challenges might be the possibility that they partition into and modulate the structure of the bilayer. To test this hypothesis, we ran two 60 ns pMD-membrane simulations without the Lennard-Jones non-bonded interaction modifications described above. We found that a significant fraction of our test probe, isopropanol, quickly partitioned to the POPC/POPS bilayer While the impact of the isopropanol-lipid interaction on the bilayer structure might appear relatively small, it can have substantial impact on the dynamics of the bound protein. Moreover, we anticipate larger effect for bigger and more lipophilic probes such as benzene or cyclohexane. Therefore, it is desirable to avoid probe partitioning into the bilayer in order to ensure that pMD-membrane will have broad application. To achieve this, we modified the Lennard-Jones potential between selected atoms of the probe and lipids (see Methods). As shown in Fig 3, this modification led to a smooth decline in the number of probe molecules that approach the polar head group, and exclusion of probe molecules from the hydrophobic core. This is reflected in the average area per lipid and bilayer thickness, both of which remained unaffected (59. 7 ± 0. 4 Å 2 vs. 60. 8 ± 0. 7 Å 2 and 41. 2 ± 0. 2 Å and 40. 5 ± 0. 3 Å in the presence and absence of probes, respectively). The very small polarization of the few probe molecules near the bilayer surface (Fig 3C ) can be eliminated if needed by using a larger non-bonded inter-particle distance for the modified atom pairs. Our choice of parameters was meant to ensure that probe molecules can approach the protein from the side of the membrane surface as well as from bulk. 

In summary, comparison of Figs 2 and 3 makes it clear that a simple modification of some of the vdW terms on selected atoms of the probe and bilayer extends the applicability of pMD to membrane proteins, a major focus of many drug discovery campaigns (e. g. [19] ). In the subsequent sections we demonstrate the application of the method on membrane-bound K-Ras, a highly sought after anti-cancer drug target [30, 66].
section_name: Preventing partitioning of probe molecules into bilayer

section_content: Protein motion can be affected by the composition of the surrounding solvent [67, 68]. As can be surmised from Fig 1, the dynamics of the catalytic domain of K-Ras G12D is different before and after it formed direct contact with the bilayer (i. e., when fully in water and after part of the surface is restrained by interaction with lipids). Therefore, we checked if (i) pMD-membrane qualitatively reproduces binding sites on Ras that have been previously characterized by other solvent mapping techniques [69, 70], and (ii) these sites/sub-sites are modulated by conformational change induced by membrane binding. 

In the previous pMD analysis of G12D in solution [11], we identified five druggable sites and three sub-sites (see Fig 3 in ref. [11] ). A detailed comparison with known ligand binding pockets indicated that three of the predicted druggable sites overlap well with pockets p1, p2 and p3 while two of the predicted sub-sites were found to be parts of p1 and p4, respectively (see Fig 4 of ref. [11] ). Allosteric pockets p1 to p4 have been previously described in detail: p1 is the binding site of ligands reported by Maurer et al [32], Sun et al [33] and Shima et al [34] ; p3 is near the C-terminus of the catalytic domain where Cu 2+ -cylen binds [35, 36] ; p4 is the proposed target site of Andrographolide derivatives [38] and Zn(II)-bis(2 picolyl)amines [36] ; there is no known non-covalent binder that targets p2 but covalent ligands that target this region have been reported [37]. Direct comparison of the calculated probe occupancies in membrane-bound K-Ras with experimental results is not possible because one cannot turn on Previously characterized pockets p1 (near β1-β3/h2), p2 (between h2/h3), p3 (near h5 and loop7) and p4 (behind s1) are labeled. The last 10 ns of the current pMDmembrane runs and last 10 ns of the longest run of G12D K-Ras in solution from ref [11] were used to calculate the probe densities and grid free energies. The protein from solution was first aligned to each of mode 1 and mode 2 membrane-bound forms prior to grid-based calculations. The protein is in yellow and a portion of the inner leaflet of the bilayer is shown as iceblue lines. See ref [11] for details. doi:10. 1371/journal. pcbi. 1004469. g004 in experiment non-physical repulsive interactions to prevent bilayer partitioning of isopropanol. Nonetheless, the following analysis provides a strong evidence that pMD-membrane is able to identify true drug binding sites. 

In order to compare the current pMD-membrane with the previous pMD-solution, we performed grid free energy analysis on the last 10 ns data from the two current simulations of G12D K-Ras and the previous pMD-solution run of G12D K-Ras (we chose the longest, 100 ns run). The results are displayed in Fig 4, where grid densities that yielded grid free energies of -1 kcal/mol or lower are shown in blue and green for the simulations in membrane and solution, respectively. Pockets p1 to p4 are labeled where possible. One can see that there is a remarkable overall agreement between the membrane and solution simulations (see overlaps between the blue and green iso-surfaces). We take this as validation of pMD-membrane, because pockets p1-p4 are all confirmed ligand binding sites for which there exist crystallographic or solution NMR structures of K-Ras in complex with ligands [32] [33] [34] [35] [36]. However, there are also clear differences. The most significant differences include the following. (i) In membrane binding mode 1 (Fig 4A ), pocket p1 is completely invisible. Instead, a new putative site appears very close to the P-loop but distinct from the nucleotide-binding site. (ii) In membrane binding mode 2 (Fig 4B ), p2 is absent but no new site is discovered. Taken together, these results demonstrate that pMD-membrane captures known druggable sites, and that protein-membrane interaction modulates binding site accessibility. Whether the observed differences in some of the sites will translate into differential ligand binding in the soluble and membrane-bound K-Ras is yet to be determined. Nonetheless, this observation highlights the importance of incorporating the effect of membrane in Ras drug discovery efforts.
section_name: Comparison of druggable sites on K-Ras predicted by pMD in membrane and in solution

section_content: For the analysis in the previous section, we included only probe molecules that lie within 4 Å of protein heavy atoms. In principle, grid occupancy can be calculated over the entire system (see Methods, Eqs 3 and 4), so that all high-density grid points around the protein can be considered. Then, the probe density can be visualized at the desired concentration cutoff. An example of this is shown in S4 Fig. There is a clear overlap between the high-density iso-surfaces and known pockets p1 to p4, as well as the sites highlighted in Fig 4A. However, there are also a large number of other high-density regions that, though unlikely to bind drug-like molecules, clutter the picture. 

In Fig 5, we show overlays of probe-occupancies derived from distance-based and gridbased analyses (see Methods). As expected, the two techniques yielded very similar results (notice the overlap between the red and white contours representing high probe densities from distance-based and grid-based analysis, respectively). White isosurfaces circumscribed by the red contours are likely pocket-like, suggesting that a combined use of distance-based and gridbased occupancy analyses would help localize relevant sites somewhat. However, it is still difficult to unambiguously isolate potentially druggable pocket-like sites. This can be regarded as a limitation of probe-based analyses in cases where there is no prior knowledge of druggable sites. This suggests that it is prudent to complement pMD with analysis of geometric/chemical features such as curvature, volume and polarity. There are a number of useful tools to perform structure-based pocket analysis, such as SiteMap [71], MDpocket [72] and AutoGrow [73]. The use of more than one type of probes or mixtures thereof may also be helpful. 

For the purposes of the current work, where we are interested in the relative druggabilty of different mutants/orientations of membrane-bound K-Ras, an inherently less cluttered differential probe occupancy (or density) is most relevant. Both grid-based and distance-based occupancy measures can be used for such an analysis, but we found the latter to be more convenient. The following sections will therefore focus on changes in distance-based probe occupancies.
section_name: Visualizing probe densities--challenges and proposed solutions

section_content: As noted earlier, the current work was motivated in part by the observation that bilayer interactions of H-Ras G12V [54] and K-Ras G12D (Prakash and Gorfe, to be published) involve at least two distinct modes. The prime difference between the modes is the orientation of the catalytic domain with respect to the membrane plane so that, in the case of K-Ras G12D, either helix 3/4 or helix 2 directly contact the bilayer (Fig 1 ). We wanted to see if these two membrane-binding modes differ in ligand binding potential when assessed by pMD-membrane. To this end, we calculated the difference in atomic probe occupancy between the simulations started from the conformation in Fig 1B (mode 2 ) and from the conformation in Fig 1A (mode 1 ): ΔR conf = R 2 -R 1. Thus, negative ΔR conf at a given atom means that the atom is more accessible to probes in membrane binding mode 1 than 2 (positive ΔR conf is the opposite). 

The data in Fig 6A shows that the two membrane binding modes substantially differ in probe occupancy, particularly at helices h2, h3, h4, the hvr and to a lesser extent between h1 and β2. Coloring the 3D structure by ΔR further shows that the differences are confined to four surface patches (Fig 6B ). Three of these patches correspond to previously described pockets, including p1, p3 and part of p2. Pocket p3 is more accessible in mode 1 than in mode 2 (negative ΔR conf ) while p1 is more accessible in mode 2 than in mode 1. There are some changes in the accessibility of p4 as well. Apart from these pockets, variations in probe accessibility include surface sites that may not be druggable, such as the sharp positive ΔR conf peak on helix 4 arising from the fact that it is engaged with the bilayer in mode 1 but not in mode 2. 

Overall, this analysis demonstrates that sites' accessibility to probe molecules is a function of membrane orientation of the protein, and that pMD-membrane is capable of capturing those differences. We propose that, at least for Ras proteins, ligand accessibility will likewise depend on the details of membrane binding.
section_name: Ligand binding sites are modulated by membrane binding

section_content: In an unpublished study, we observed that the active sites of G12D and G13D K-Ras significantly differ. While the active site of G12D K-Ras is similar to wild type, switch 1 is open and some functionally critical residues, such as Tyr32, have re-oriented in G13D K-Ras. At the functional level, G12D and G13D K-Ras differ in intrinsic GTPase activity [74] and oncogenicity [25]. Therefore, we ran a pMD-membrane simulation on G13D starting from the conformation shown in Fig 1A (mode 1 ) and calculated ΔR seq = R G13D -R G12D. We found major differences in probe accessibility of the two mutants (Fig 7A ), indicating that the two simulations started from the same initial configuration have drifted apart, leading to different probe binding propensities. The differences are largely confined to helices 2, 3, and 4 (Fig 7A ), representing two surfaces on the 3D structure (Fig 7B ). Part of the surface of helix 2 where p1 is located is more accessible in G13D than G12D whereas the region between helices 3 and 4 is significantly more probe-accessible in G12D than G13D. This is despite the fact that these regions are far away from the site of the mutation. We propose that these observations highlight potential differences in the druggability of the two mutants and thereby the possibility of isoform-specific drug leads. We find the region between helices 3 and 4 especially interesting as it might represent a potentially unique new ligand-binding site.
section_name: Impact of mutation on druggability of sites

section_content: We have seen that isopropanol has preferred orientations at the glycerol and head group regions of the POPC/POPS bilayer (Figs 2 and 3 ). This was because the OH functional group of the probe prefers to interact with lipid oxygen atoms which the CH3 groups tend to avoid. For the same reason, interaction of the probes with protein atoms is likely to be polarized so that the OH group points away from hydrophobic surface cavities but points toward polar cavities. Therefore, we wondered if the local orientation of the probe contains information about the polarity/hydrophobicy of individual sites. To check this, we calculated the average orientation of the probe taking into account every protein atom that is in contact with a probe molecule. In this analysis, positive <cos(θ r )> indicates that the hydroxyl oxygen points away from the protein (see Methods). We found that the ensemble averaged cos(θ r ) is positive for the vast majority of the highly probe-accessible surface protein atoms (Fig 8A and 8B ), suggesting that the probe-binding sites are mostly hydrophobic and therefore potentially druggable. Negative <cos(θ r )> was found only at a couple of surface sites that are unlikely to bind to ligands.
section_name: Analysis of probe orientation

section_content: We introduced a technique referred to as pMD-membrane as a novel approach for the analysis of ligand binding potential of surface cavities in membrane proteins. This represents a major expansion of the scope of probe-based molecular dynamics simulation approaches. The goal of pMD-membrane, and probe-based simulations in general, is to map fragment positions in pockets that may then be used by medicinal chemists to design specific binders. Extensive analysis involving multiple types of co-solvents and multiple drug targets found a strong correlation between probe occupancies and binding affinities of true binders in most cases (refs 7-9). Similarly, we have shown that, following modification of selected vdW interaction terms in the force field, pMD-membrane was able to identify allosteric ligand binding sites (including known binding sites) on the surface-bound K-Ras without any significant effect on the structure and dynamics of the bilayer or the protein. We have also demonstrated that pMD-membrane can capture the impact of conformational changes induced by membrane binding or mutation on the probe accessibility of putative druggable sites. This is important because the ultimate goal of any site identification scheme is to differentiate cryptic binding sites based on changes in size, location or chemical feature. Such changes can result from small differences in protein motion in water versus membrane environments, as well as from mutations, substrate binding, posttranslational modification etc. Our method thus extends the scope of probe-based molecular dynamics simulation in two majors ways: as a novel means by which to find (allosteric) ligand binding sites in membrane proteins and as a tool with which to probe differential ligand accessibility in closely related targets. The method can be easily expanded to any type of probe or mixture of probes through similar modifications of non-bonded terms. Therefore, pMD-membrane and the analysis tools described in this study are applicable to a wide variety of membrane proteins, whether trans-membrane or surface-bound.
section_name: Conclusion

section_name: Supporting Information

section_name: Acknowledgments

section_content: We thank the Texas Advanced Computing Center (TACC) for computational resources.
section_name: Acknowledgments

section_content: This work was funded by a Postdoctoral training fellowship from the Keck Center Computational Cancer Biology Training Program of the Gulf Coast Consortia from Cancer Prevention and Research Institute of Texas (Grant No. RP140113 ); Cancer Prevention and Research Institute of Texas (Grant No. DP150093 ); and National Institutes of Health General Medical Sciences (Grant No. R01GM100078 ). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

section_content: All relevant data are within the paper and its Supporting Information files.