Mixed-Solvent Molecular Dynamics Simulation-Based Discovery of a Putative Allosteric Site on Regulator of G Protein Signaling 4

Structure Preparation and System Setup

The ensemble-averaged, solution NMR structure of rat RGS4 (PDB: 1EZT) was employed for the current study,¹⁷ and its structure was prepared with Molecular Operating Environment (MOE, version 2019.0101) for molecular dynamics simulations. The N- and C-termini were capped with acetate and N-methyl amide, respectively. Hydrogens were added, and the protonation states of all residues were set to pH 7.0. In silico mutations (K99A/K100A) for RGS4 were performed using the Protein Builder module from MOE.

From AmberTools18, tleap was utilized to prepare the relevant input files using the prepared structure with the AMBER ff14SB force field.¹⁸ The protein was surrounded with a layer of chemical probes (pyrimidine, acetonitrile, or isopropyl alcohol) before layering the rest of the system with TIP3P water. The concentration of probe solvent was 5% v/v.¹⁶ The resulting systems were neutralized with counterions and had a dimension of approximately 78 Å × 75 Å × 71 Å. An example is shown for RGS4 layered with pyrimidine in Figure S1.

Mixed-Solvent Molecular Dynamics (MixMD) Simulations

Molecular dynamics simulations were run using the GPU implementation of PMEMD from AMBER18.¹⁹ A first round of minimization was performed on the systems with a 10 kcal/mol·Å² harmonic restraint on the protein, which consisted of 2,500 steps of steepest descent followed by 2,500 steps of conjugate gradient. The second round of minimization was the same as the first apart from removing the harmonic restraint. Subsequently, the systems were gradually heated in the NVT ensemble from 0 K to 100 K for 12.5 ps, then from 100 K to 310.15 K in the NPT ensemble for 125 ps at 1 bar; a 1 fs timestep was used, along with a 10 kcal/mol·Å² harmonic restraint on the protein. After reaching the target temperature, the systems were equilibrated in the NPT ensemble at 310.15 K and 1 bar with a 2 fs time step. Beginning with a 5 kcal/mol·Å² harmonic restraint on the protein, the restraint was reduced by 1 kcal/mol·Å² every 500 ps for a total of 2.5 ns, then by 0.1 kcal/mol·Å² every 500 ps for a total of 5 ns. Production simulations had no constraints and were run in the NPT ensemble at 310.15 K and 1 bar for 50 ns with a 2 fs time step. All bond lengths concerning hydrogen were constrained with the SHAKE algorithm. Non-bonded interactions were set to cut off at 9.0 Å, while Particle mesh Ewald summation was used for long-range electrostatics. The Langevin thermostat and Monte Carlo barostat were used where applicable. Periodic boundary conditions were applied to all heating, equilibration, and production runs. For both wild-type and mutant RGS4, ten independent production runs were carried out for each of the three probes (pyrimidine, acetonitrile, and isopropyl alcohol), resulting in a total of 3 μs simulation time.

Probe Occupancy Calculations

Using the ptraj program from AmberTools18, all of the resulting trajectories for each probe were centered, imaged, and aligned, then the locations of the probe atoms from the last 25 ns of each independent simulation were combined and binned onto a grid with 0.5 Å spacing. To compare across conditions, Z-scores were then calculated using the following equation to normalize the occupancies:

where x_i is the occupancy at grid point i, μ is the mean occupancy of all grid points, and σ is the standard deviation of occupancy at all grid points. The normalized occupancy maps for each probe could then be visualized like electron densities with PyMOL (version 2.3.5).²⁰

Assessment of Predicted Allosteric Sites

MixMD Probeview²¹ was used to identify potential allosteric sites from the occupancy data using the averaged protein structure from each probe in PyMOL. The default parameters for DBSCAN clustering were set at occupancy cutoff= 0.1, ε= 3, and minimum number of points= 10.

To simulate the effect of having a bound ligand in each site, an octahedron composed of six carbon dummy atoms was manually placed in regions where all three probes had overlapping occupancies, adapted from that performed in Panjkovich and Daura (2012).²² The edges of the octahedron (i.e. carbon-carbon distance) were 1.5 Å in length. With these pseudo-ligands positioned on the protein, further analyses involving its effects on protein dynamics could then be performed.

Custom R scripts were written using the Bio3D package (version 2.4–1) to assess the impact of the pseudo-ligand on overall protein motion.^23,24 Normal mode analysis (NMA) calculations were performed using an all-atom elastic network model (ENM) with the aaenm force field. The lowest-frequency modes, which correspond with large-amplitude conformational changes, were used in the study. Theoretical temperature factors were calculated as follows for each non-hydrogen atom:

where 〈μ²〉 is the mean squared displacement. Wilcoxon signed-rank tests were used for an atom-wise comparison of predicted temperature factors between the bound (holo) and unbound states (apo), generating two-tailed p-values. P-values < 0.001 were considered significant. Additionally, coarse-grained ENM using the calpha force field was used to generate the dynamic cross-correlation matrix.

Probe Kinetics

Fragments have a low binding affinity and so it is challenging to assess the binding affinities of the probes from simulation studies. However, direct evaluation of on-off rates and binding kinetics could be done if there is constant exchange of probes with bulk solvent. This resembles an equilibrium scenario of enough binding and unbinding phenomena. Caflisch and coworkers have reported an elegant route for assessment of probe binding kinetics based on on-off rates from MD simulations.²⁵ In a recent publication, Pan et al extracted on-off rates for drug like fragments with millimolar affinity using long scale simulation data.²⁶ They performed microsecond long simulations and observed spontaneous binding and unbinding repeatedly occurring throughout the production runs. In this work, the large number of binding and unbinding events to the allosteric site detected, allows for extraction on on-off rates from the MixMD simulations.

The K_D values (mM) were estimated using the equation:

Conformational Clustering of Site 3

POVME²⁷ was employed to explore the possible binding site conformations that Site 3 could adopt. One hundred evenly-spaced frames were extracted from the last 25 ns of each trajectory and aligned to the starting structure. An inclusion sphere with a radius of 12 Å was centered on Site 3 with a grid spacing of 1 Å. Convex hull exclusion was applied with the first frame as reference to remove extraneous points outside the potential allosteric site (i.e. out in solution). Tanimoto coefficients were calculated on overlapping grid points between each frame, and hierarchical clustering with a Kelly penalty was performed.

Molecular Docking

An SDF file consisting of 10,240 compounds from the ‘Discovery Diversity Set’ was obtained from Enamine, which was then imported into a MOE database. All compounds were subjected to the ‘Wash’ module; the dominant protonation state at pH 7.4 was set for each compound, followed by the generation and minimization of the 3D structure.

GOLD²⁸ was used for molecular docking experiments. Representative structures from the top 5 most populated conformational clusters of Site 3 from RGS4 were used for ensemble docking. A point was specified in the center of Site 3, and all residues within 12 Å were considered part of the binding site. Ring conformations were explored using the ‘flip ring corners’ and ‘match template conformations’ options, while pyramidal nitrogens and amide bonds were allowed to flip. For each compound, 10 independent genetic algorithm runs were used with the ‘ensemble’ search efficiency, and conformations were assessed with the ChemPLP scoring function. The top-scoring compound was ultimately saved for analysis, in addition to the RGS4 conformer with which it was associated.

Prediction of Allosteric Sites on RGS4

MixMD simulations were performed using the solution NMR structure of rat RGS4 (PDB: 1EZT).¹⁷ Utilizing the three small organic probes, pyrimidine, acetonitrile, and isopropyl alcohol, binding hotspots were identified through visual analysis of probe occupancies from the simulations (Figure 1A). In order to identify and rank the predicted binding sites, we applied MixMD Probeview²¹ to the probe occupancies, and subsequently the top-five predicted binding sites were chosen for further analysis (Figure 1B). As compared to using a more stringent level of significance (35σ), a more complete manifestation of the sites was observed at a lower level (20σ). In particular, sites 1, 2, and 3 possessed very strong probe density, while sites 4 and 5 did not map as well (Figure 1). This suggests that the former 3 sites have a greater propensity for being physiologically relevant. A small pocket was mapped to 35σ by all three probe hotspots at Site 2 (Figure S3). This site is known from a previous study of the structure of G_iα1 in complex with RGS4 (PDB: 1AGR) where G_iα1 interacted with RGS4 primarily through three switch regions.²⁹ The Thr-182 from switch I (residues 176–184), which is completely buried by a pocket formed by conserved residues from RGS4,²⁹ overlaps with the predicted pocket from the 35σ contour map (Figure S2). By lowering the contour to 20σ, we revealed additional sites where hotspots from multiple probes overlapped, albeit with increased noise (Figure 1A, bottom). The binding site designated as Site 2 was consistent with the pocket that accommodated Thr-182 from switch I of G_αi1.

In addition to Site 2, we also examined the remaining sites in detail. Site 1 was star shaped and contained two distinct hotspots; one was strongly mapped by all three probes, while the other was mapped by pyrimidine and acetonitrile (Figure S3). Similarly, Site 3 was located on a cleft opposite that of Site 2 and also contained two distinct hotspots; one was strongly mapped by all three probes, while the other was mapped only by acetonitrile (Figure 2A). Though hotspots within Sites 4 and 5 were present at 20σ, their probe occupancies were much more spurious at 35σ, unlike the first three sites (Figure 1, Figure S3).

To supplement our MixMD-based predictions, we submitted the RGS4 structure (PDB: 1EZT) to both the FTSite and PARS servers, both of which operate on a static protein structure.^30,31 The PARS server uses the method from Panjkovich and Daura (2012), which predicts allosteric sites based on NMA.^22,31 Conversely, FTSite docks chemical probes onto a static protein structure, and sites are identified by consensus clusters of probes.³² Three sites were found using the FTSite algorithm (Figure S4A); the top two sites could be considered one continuous site and overlapped with Site 3 (Figure S4A, pink and green), while the last site overlapped with Site 2 (Figure S4A, blue). Independently, the PARS server identified binding sites that also corresponded with the same sites (Figure S4B, orange). Combination of the three approaches resulted in a consensus of Sites 2 and 3.

Effect of Octahedron Pseudo-ligands on Protein Flexibility with NMA

To examine the impact of a simulated allosteric modulator on protein flexibility, a single carbon octahedron pseudo-ligand was manually placed into locations in which the probe occupancies had the greatest agreement (Figure 2Ba; Figure S3, top). This method was previously shown by Panjkovich and Daura (2012) to aid in the positive identification of known allosteric sites in a retrospective study without a priori knowledge.²² Subsequently, NMA was performed to compare differences in flexibility between the apo and holo structures. Using the first 10, 20, and 30 non-redundant normal modes, no difference was observed for Sites 1, 4, and 5 (Table S1). However, the placement of an octahedron at the high-confidence hotspot at Site 3 produced a statistically significant change in flexibility for RGS4 within all relevant ranges of normal modes (Table S1, Site 3 - a). Known to interact with G_iα1, Site 2 produced an appreciable difference in flexibility using the first 10 non-redundant normal modes (p-value= 0.076), though it was not statistically significant. The mechanism by which G_iα1 interacts with RGS4 to enhance its GTPase activity requires interaction with all three of its switch regions,²⁹ so it is expected that the sole perturbation of Site 2 on switch region 1 will not change the dynamics of RGS4 substantially, though it did appear to have a marginal impact on protein flexibility. Both Sites 2 and 3 were consistent with predictions from the PARS server and FTSite (Figure S3). Taken together, we hypothesized that Site 3 potentially harbors allosteric activity.

Identification of Site 3 of RGS4 as Potential Allosteric Site

From visual inspection of the hotspots at Site 3, we observed additional acetonitrile occupancy at 20σ (Figure 2Bb), which could be connected to the high-confidence hotspot based on the occupancies using MixMD Probeview to form a contiguous site (Figure 2A). Interestingly, this site encompasses one (Lys-99) of the two residues (Lys-99/Lys-100) that are known to interact with the endogenous modulators, PIP₃ and calmodulin.¹⁴ Hence, we performed MixMD on a mutant RGS4 with both residues mutated to alanine. This resulted in Site 3 largely remaining the same, though a slight increase in probe occupancy was observed where Lys-99 was replaced with alanine (Figure S5); this suggests that the lysine residues are likely not solely responsible for forming an allosteric site. Since the high-confidence hotspot at Site 3 was shown earlier to affect protein flexibility with the presence of the octahedron (Figure 2Ba), we decided on placing an additional octahedron on the acetonitrile hotspot in Site 3 to simulate a ligand interacting at both subsites (Figure 2Bb). This produced a more significant effect on protein flexibility than placing an octahedron at the high-confidence hotspot on Site 3 alone and decreased the overall predicted temperature factors in the holo state as compared with the apo state (Figure 2C; Table S1, Site 3 – a+b). Additionally, it is evident that RGS4 does not exhibit drastic global movements in the apo state based on NMA, but the octahedrons appear to reduce RMSF in most residues in the holo state (Figure 3). Consequently, we sought to elucidate the effect of the pseudo-ligands on perturbation among each of the individual modes.

As the low-frequency modes represent the collective dynamics of a protein,³³ we examined the impact of the octahedrons on each of the first ten non-redundant modes. Normal modes 8 and 12 resulted in the most statistically significant changes in protein flexibility, though modes 7, 10, 13, and 14 had a noticeable impact (Table 1). NMA trajectories for modes 8 and 12 are shown in Figure 4. The motion of RGS4 could be summarized by designating dynamic units comprised by two halves of the protein: H1 contains both N- and C-termini (Figure 4, left half), while H2 consists of the remainder of the protein (Figure 4, right half). Mode 8 describes a twisting motion between H1 and H2 of the protein demarcated by the cleft formed by Site 3. Conversely, mode 12 demonstrates a similar twisting motion as mode 8 with H1, while H2 simultaneously moves inward towards the cleft. From the NMA trajectories, we observed what appeared to be correlated motions between H1 and H2 (Figure 4). The movement between the two halves from both modes 8 and 12 are such that the presence of the octahedrons could potentially impinge overall movement. Therefore, we next examined the effect on regions of RGS4 with correlated motion.

We generated a dynamic cross-correlation map from NMA to assess whether the presence of the octahedron pseudo-ligands in Site 3 would affect coupled atomic fluctuation from non-contiguous residues (Figure 5A). Overall, correlated motion was decreased in the holo state, which was consistent with our observation of decreased predicted temperature factors. More specifically, coupling was markedly reduced in three regions: 1.) between α2 and α4-α5, 2.) between α3 and α5, and 3.) between α4 and α6 (Figure 5A and 5B). Intriguingly, these regions correspond with the motion described with normal modes 8 and 12, which were negatively impacted by the presence of the octahedron pseudo-ligands (Figure 4). Moreover, α2, α4, and α7 form Site 3, of which the former two helices overlapped with the regions of reduced coupling. Overall, this indicates a potential role of Site 3 in an allosteric mechanism.

Dissociation Constants and Analyses on On-Off Rates Kinetics of Site 3

To assess the probe binding kinetics to the allosteric pocket, we used the method Caflisch proposed for analyzing DMSO binding-unbinding rates to the FKBP protein.²⁵ Our water miscible probes are all less than 6 heavy atoms, and they diffuse well in the 20–25 ns simulation time. The production runs for each probe type were analyzed to track the closest distance between the center of the hotspot and the probe in question (Figure 6). The analyses were done using two different distance thresholds; probes were considered bound when they were within 3Å of the hotspot and unbound when they were over 6Å. For verifying robustness of the calculations, a 7Å distance cut-off was also employed for defining unbound states. Apart from our main goal of estimating K_D values from the simulations, usage of two different distance cut-offs showed that the affinities computed were not sensitive to the distance metric used to define the bound and unbound states. Using these thresholds, the trajectories were partitioned into bound and unbound states. The time lengths of each event were tallied, and a distribution of dwell times was estimated. The cumulative distributions of the binding and unbinding events were plotted and fitted to exponential decay plots and on-off times extracted. The fitting equation and plots are in shown in Figures S6 – S9, while the values are reported in Table 2. We also used 7Å as the distance threshold to probe unbinding and the values obtained were similar (Table 2). This reinforces Caflsich’s hypothesis²⁵ that the binding affinities are not extremely sensitive to the cut-off distances employed in our calculations. The values for on–off rates for using 3Å and 7 Å cutoffs are tabulated in Table 2. We do find that the concentration of the cosolvent in the simulation directly impacts the calculated affinities.

Characterization of Site 3 of RGS4

The full predicted allosteric site (Site 3) is shown in Figure 7 and is comprised of the following residues: Glu-64, Ile-67, Asn-68, Trp-92, Glu-96, Lys-99, Asp-150, Lys-154, Asn-158, and Glu-161. Three of the α-helices from RGS4 (α2, α4, and α7) contributed residues to form this cleft. Site 3 consists largely of polar residues, while also harboring 5 charged residues (Glu-64, Glu-96, Lys-99, Lys-154, and Glu-161). Also, the bottom of the cleft is lined with aromatic and nonpolar residues. The cleft itself is saddle shaped and appears as though it would best accommodate a ‘linear’ small molecule (Figure 1B), as will be examined in the next sections.

To assess the conformational diversity of Site 3, we performed clustering on the cleft from the last 25 ns of all trajectories for each probe with POVME.²⁷ Overall, the cleft appeared to maintain an average volume of ~1,400 ų (Figure S10A). Moreover, much conformational diversity was observed among the 25 clusters, though the top 5 most populated clusters harbored the majority of structures (Figure S10B). Cluster representatives from each of the top 5 clusters are displayed in surface representation in Figure S10C.

Apart from the conformational diversity of Site 3, we were also interested in identifying potential compounds which seek out this site. As a result, we performed molecular docking with GOLD²⁸ against the ‘Discovery Diversity Set’ from Enamine using a representative structure from each of the top 5 conformational clusters. The docking scores followed a normal distribution centered around a docking score of ~55 (Figure S11). Moreover, manual inspection of the top 100 ranked compounds revealed that 65 of them interacted with RGS4 where the 2 pseudoligand-bound hotspots were located in Site 3 (Figure 2B, Figure 8A). Perhaps unsurprisingly, most of them could be described as ‘linear’ small molecules. , The top-scoring ‘hit’, Z346633068, had a docking score of 82.48, significantly higher than the next highest compound, Z2495891309, which had a docking score of 75.34 (Figure 8B, Figure S11). While the docking score for Z346633068 was very high, we observed a cis-peptide bond with its pose (Figure 8C); though this is generally an unfavorable conformation, the remaining interactions made with RGS4 likely offset the energetic penalty, given the high docking score. Much of the predicted interaction between No. and RGS4 appears to be governed by hydrophobic interactions with the nonpolar/aromatic residues in the cleft. Furthermore, the phenyl-2-pyrrolidinone portion of the compound overlaid with all three probe densities of Site 3a, while the -tetrahydroisoquinolone moiety interacted with Site 3b (Figure 8C, Figure 2B). However, an important insight from this pose is the potential for the formation of a salt bridge between the tertiary amine of the tetrahydroisoquinolone moiety of Z346633068, which is charged at pH 7.4, and Glu-64 of RGS4 (Figure 8C). In fact, this same salt bridge is seen for 17 other compounds out of the 65 potential hits, suggesting that Glu-64 could play a key role in ligand binding to Site 3.