A Curvilinear-Path Umbrella Sampling Approach to Characterizing the Interactions Between Rapamycin and Three FKBP12 Variants

Abstract

Rapamycin is an immunosuppressant macrolide that exhibits anti-proliferative properties through inhibiting the mTOR kinase. In fact, the drug first associates with the FKBP12 enzyme before interacting with the FRB domain of its target. Despite the availability of structural and thermodynamic information on the interaction of FKBP12 with rapamycin, the energetic and mechanistic understanding of this process is still incomplete. We recently reported a multiple-walker umbrella sampling simulation approach to characterizing the protein–protein interaction energetics along curvilinear paths. In the present paper, we extend our investigations to a protein-small molecule duo, the FKBP12•rapamycin complex. We estimate the binding free energies of rapamycin with wild-type FKBP12 and two mutants in which a hydrogen bond has been removed, D37V and Y82F. Furthermore, the underlying mechanistic details are analyzed. The calculated standard free energies of binding agree well with the experimental data, and the roles of the hydrogen bonds are shown to be quite different for each of these two mutated residues. On one hand, removing the carboxylate group of D37 strongly destabilizes the association; on the other hand, the hydroxyl group of Y82 is nearly unnecessary for the stability of the complex because some nonconventional, cryptic, indirect interaction mechanisms seem to be at work.

Introduction

Protein–ligand interactions are central in modern drug-discovery, and their characterization by various approaches is crucial for better drug development. In this regard, computational investigation is one of the ways to acquire a deeper understanding of the interactions of interest. In particular, molecular dynamics (MD) simulations provide the physical connection between the structure and the function of biomolecules (Karplus and McCammon, 2002), especially at the atomic level; therefore, MD simulation-based techniques can often cast insights into such interactions, especially in the early stage drug-discovery (Mobley and Gilson, 2017). In addition to conformational dynamics of the interacting molecules, MD simulations are also employed to estimate thermodynamic properties such as the relative binding free energy and the binding affinity. Major ongoing challenges in this field are related to the reliability, the accuracy, and the rapidity of the estimation method/approach.

A reliable estimation of the free energy difference between thermodynamically well-defined end states of interacting molecules is one of the major goals of computational biophysics. Umbrella sampling along a chosen reaction coordinate, followed by potential of mean force calculations, is one of the ways widely used to estimate binding affinities (Torrie and Valleau, 1977; Kästner, 2011). However, umbrella sampling along predefined vectorial reaction coordinates followed by potential of mean force (PMF) profile constructions has some serious concerns in performing reliable energetic calculations (Doudou et al., 2009). Recently, we discussed them for the protein–protein systems and proposed a naïve multiple-walker approach, in which independent umbrella sampling simulations were conducted without predefined vectorial reaction coordinates. We observed similar large variations in the values of converged PMF profiles, resulting from different curvilinear paths. The variations were attributed as due to the different excessive dissipations in different paths taken, and therefore the lower-bound PMF was chosen, and by introducing a correction term derived from statistical mechanics, the standard free energy of binding was estimated for the protein–protein complex system (Joshi and Lin, 2019). The estimations were in good agreement with the experimental values for the barnase–barstar complex. Furthermore, the revealed mechanistic details from our simulations, e.g., the physical pathways/trajectories of dissociation/association, were quite consistent with two major physical pathways that were determined from several milliseconds-long adaptive molecular dynamics simulations reported by Plattner et al. (2017). Thus, the proposed approach is quite useful in maintaining a suitable balance between estimation of binding energetics and revealing underlying mechanistic details of the dissociation reaction within much less computational cost than the brute force approaches. Since the sampling is enhanced along a spontaneously evolved (i.e., non-predefined) curvilinear physical trajectories, the overall approach is referred to as curvilinear-path umbrella sampling (CPUS) approach. In the present paper, we extend our previous work to the interactions of a drug with its protein partner.

This implementation of CPUS will be validated on FKBP12, a model used in computational ligand binding studies for more than 20 years (Pearlman and Connelly, 1995; Lin et al., 2002; Lin et al., 2003; Sun et al., 2003; Swanson et al., 2004; Fujitani et al., 2005; Lee and Olson, 2006; Olivieri and Gardebien, 2011; Dickson and Lotz, 2016; Nerattini et al., 2016). This enzyme participates in protein folding by catalyzing the isomerization of amide bonds adjacent to proline residues (Hamilton and Steiner, 1998; Galat, 2013). Discovered because it could bind the immunosuppressant macrolide FK506, it was further shown that it could also complex with Rapamycin, a parent drug sharing the same pharmacophore and similar therapeutic indications (Van Duyne et al., 1993; Hamilton and Steiner, 1998; Galat, 2013). It later turned out that the inhibition of the FKBP12 enzymatic activity was irrelevant to account for the pharmaceutical properties of the two macrolides: In fact, this protein only potentializes the binding of each drug to its effective cellular target through the formation of a specific ternary complex, with calcineurin for FK506 and with mTOR for rapamycin (Choi et al., 1996; Banaszynski et al., 2005; Galat, 2013). Apart from its biological significance, the choice of FKBP12 as a model in the 1990s can be explained by its small size (107 residues) and its compact structure, both key factors at a time where computational power was limited. Moreover, several X-ray crystal structures were quickly released, in free form as well as in complex with FK506 and rapamycin (Van Duyne et al., 1991; Van Duyne et al., 1993; Wilson et al., 1995). Additionally, with time, a huge amount of biophysical data have been accumulated from thermodynamic and kinetic measurements on the complexation reaction (Bierer et al., 1990; Holt et al., 1993; Bossard et al., 1994; Connelly et al., 1994; Luengo et al., 1995; DeCenzo et al., 1996a; Schuler et al., 1997; Wagner et al., 1998; Graziani et al., 1999; Dickman et al., 2000; Banaszynski et al., 2005; Wear et al., 2007; Wear and Walkinshaw, 2007; Shor et al., 2008; Kozany et al., 2009; Wu et al., 2011; Tamura et al., 2013; Singh et al., 2015; Lu and Wang, 2017; Kostrz et al., 2019; Wang et al., 2019) to NMR investigations on the protein dynamics (Sapienza et al., 2011; Yang et al., 2015; Solomentsev et al., 2018).

As far as in silico studies are concerned, FKBP12 is extensively used as a model system to study protein–ligand interaction energetics. For instance, the relaxed complex scheme (Lin et al., 2002) was employed to estimate interaction energetics of FKBP12 with several different ligands (Lin et al., 2003). Another study using FKBP12 was performed to establish the groundwork for the end-point free energy methods, in which the theoretical framework was proposed to calculate the association free energy (Swanson et al., 2004). A study on estimation of absolute binding free energy calculations of FKBP12 and eight ligands was carried out, where the Bennett acceptance ration (BAR) method was employed in the direct calculations (Fujitani et al., 2005). The FKBP12•ligand system was used to estimate the binding free energies with two ligands, 4-hydroxy-2-butanone and FK506, respectively. Although the necessity of sampling along curvilinear paths was mentioned, the theoretical framework for the free energy estimation was developed only for linear/vectorial paths, and with some quadratic approximations for the variance along the principal axis (Swanson et al., 2004; Lee and Olson, 2006). The CPUS approach can be considered as an alternative approach without such approximations.

From a structural point of view, the macrolide binding site (Van Duyne et al., 1991; Van Duyne et al., 1993) significantly overlaps with the FKBP12 catalytic cleft, in agreement with the observed catalytic inhibition (DeCenzo et al., 1996a; Ikura and Ito, 2007). Rapamycin, the only ligand we will study here, binds in the cavity located between the short α-helix and the five-stranded β-sheet that is wrapped around it. More specifically, the drug pipecolinyl ring is deeply buried inside the protein (Figure 1A) and is involved in hydrophobic interactions with the aromatic side-chains of residues Y26, F46, W59, and F99 (Figure 1B; Supplementary Figure S1)—see (Van Duyne et al., 1993; Sun et al., 2003) for lists of the atom-pair contacts. In addition to auxiliary hydrophobic interactions with some of the amino acids surrounding the cavity (i.e., F36, V55, I56, H87, Y82—Figure 1B), rapamycin is retained at the FKBP12 surface through five hydrogen bonds. Of the hydrogen bond-implied residues (Figures 1B,C), D37, I56, and Y82 are also positioned on the rim, whereas Q53 and E54 are more distant. Since the first three amino acids are the most conserved (Van Duyne et al., 1993; Hamilton and Steiner, 1998; Galat, 2013), they seemed the most attractive for performing an in silico mutagenesis program aiming at demonstrating the strengths of the CPUS approach. More precisely, our goal was to test if we could predict the influence of H-bond removal on the stability of the FKBP12•rapamycin complex. We finally selected to target D37 and Y82 because they represent significantly contrasting examples of the contribution of hydrogen bonding to ligand binding. Indeed, experimental measurements have demonstrated that the D37V substitution strongly affects the protein–drug interaction (DeCenzo et al., 1996a; Singh et al., 2015; Kostrz et al., 2019), whereas the Y82F substitution does not (Connelly et al., 1994; Pearlman and Connelly, 1995; Kostrz et al., 2019) (Figure 1D). Incidentally, beyond the determination of the binding free energies for the wild-type protein and the two mutants, we expect that the MD trajectories produced by CPUS will enable us to understand the intriguing role of the H-bond formed between Y82 and rapamycin, a bond that has an atypical geometry (Van Duyne et al., 1993), which displays a clear signature in NMR spectroscopy (Yang et al., 2015), but whose elimination has nearly no impact on the affinity (Bossard et al., 1994; DeCenzo et al., 1996b).

FIGURE 1

Materials and Methods

System Modeling

The FKBP12_WT•rapamycin binary complex was simulated using the PDB 1FKB crystal structure as a starting conformation (Van Duyne et al., 1991). For the complexes involving the Y82F and D37V mutants, we simply carried out in silico mutagenesis using the swapaa module of the Chimera molecular viewing platform (Pettersen et al., 2004) (Figure 1D). Thus, we here clearly assume that none of the two amino acid substitutions results in any significant structural change with respect to the wild-type reference. Such an assertion has already been harnessed in the FKBP12_Y82F case, crystallographic evidences being provided to support it (Pearlman and Connelly, 1995)—despite our efforts, we could not retrieve the original data from any of the usual repositories. As far as FKBP12_D37V is concerned, we performed ¹H-¹⁵N hetero-nuclear single quantum coherence NMR measurements and chemical shift perturbation analysis to back up our starting hypothesis (Supplementary Figure S2) (Williamson, 2018). Omitting the substituted amino acid, only three residues display changes larger than the mean values of the all changes plus two standard deviations: one is facing D37 in the adjacent β-strand and two are in the loop just downstream of the mutation site. As a consequence, we ruled out the possibility of any large-scale rearrangement.

The partial charges on the atoms of rapamycin were derived with the RESP scheme, calculated with Gaussian03 at the level of HF/6-31G basis set. The coordinates from the wild-type crystal structure and from the two mutant models were then processed using the antechamber and LEaP programs of the AMBER software suite (Case et al., 2005; Salomon-Ferrer et al., 2013). The well solvated system was built with 22926 TIP3P waters; 42 K⁺ and 43 Cl⁻ ions were added to neutralize the system and to mimic a 100 mM salt concentration. The box dimensions were 112.4 Å × 81.6 Å × 90.7 Å. The improved ff14SB force field along with the general amber force field (GAFF) (Wang et al., 2004) was used for all bonded and non-bonded parameters (Maier et al., 2015).

Molecular Dynamics Simulations

Production runs were conducted using the Graphical Processing Unit (GPU RTX 2080Ti) implementation of the AMBER pmemd program. All MD simulations were conducted under NPT conditions with SHAKE-enabled 2-fs time steps. The particle mesh Ewald (PME) algorithm of electrostatics was employed, and the non-bonded interaction cutoff was set to 10.0 Å. Prior to production run, all three systems were subjected to energy minimization, heated to 295 K, and then equilibrated for 100 ps. MD trajectories were analyzed using the cpptraj program of AmberTools-20 (Case et al., 2020) and several in-house shell and python scripts.

Multiple Walker Curvilinear-Path Umbrella Sampling Simulations

The successive steps necessary to implement multiple-walker umbrella sampling simulations along non-predefined curvilinear paths, as well as the corresponding theoretical framework and data treatment procedures, have been described in a previous article (Joshi and Lin, 2019). In brief, before each CPUS run, a short 10 ns-long unbiased MD simulation was conducted and the mean distance between the centers of geometry (CoG) of both FKBP12 and rapamycin was computed so as to provide a reference bound distance, . These CoG are defined by all Cα atoms for FKBP12 and all heavy atoms for rapamycin in the DISANG file, an input to the pmemd. The last frame of this unbiased MD simulation is used as the starting conformations for the CPUS simulation. In the umbrella sampling simulation, the reaction coordinate was chosen as the distance between the CoGs. The reaction path was divided into windows (i.e., the umbrella windows), and was restrained using a biasing potential with a spring constant of = 10.0 kcal/mol/Ų so as to keep the interacting molecules in the ^th distance window at distance D_ξ. The sampling in each umbrella window was enhanced sequentially so that the physical trajectory of dissociation could evolve spontaneously. To do so, the last conformation from each sampled window was chosen as the starting conformation for the next window for sampling enhancement. The distance sequence took the form and the distance per step, set to 0.1 Å. The PMF profile was constructed using well-established weighted histogram analysis methods (WHAMs) (Kumar et al., 1992; Grossfield, 2013), which removed the contribution of biased potential (Figures 2A–C). Finally, was chosen as the PMF value corresponding to the final distance of the constructed profile, . Incidentally, the PMF profile was completed for distances smaller than by conducting the umbrella sampling in the backward direction for 10 or 20 windows (depending on the cases).

FIGURE 2

For each of the FKBP12 variants, 15 umbrella sampling simulations were carried out independently by assigning different starting velocities, i.e., setting ig = −1 in the pmemd input file. Sampling in each distance window was enhanced for 1.0 ns, and CPUS MD simulations were 3.15 µs-long in total.

PMF Correction to Standard Binding Free Energy

To determine the standard free energy of binding, , each pair composed of a CPUS MD trajectory and of the associated PMF profile was further processed according to Eq. 1 (Joshi and Lin, 2019),with being the separation distance, being the standard state volume (1 663 Å³), being the PMF value at the separation distance , and . The evaluation of the second term on the right-hand side requires to know , i.e., the separation distance at which the interaction between the protein and its ligand vanished. To determine this cutoff, an interface interaction analysis was conducted using the linear interaction energy (LIE) module of the AMBER cpptraj program (Roe and Cheatham, 2013). The first cancellation of the van der Waals component provided . Then, using home-grown python and UNIX shell scripts, the bound integral included in the second term of Eq. 1 was computed and was finally determined (Table 1; Supplementary Table S1).

Curvilinear Path Tracing

The physical paths of dissociation were traced from the trajectories issued from the umbrella sampling simulations. First, for each run, the conformations were extracted at the 10-ps interval using the cpptraj module of AmberTools-20 and aligned with respect to FKBP12 only, the reference conformation being the one obtained after equilibration of the system. Doing so, the traversing of rapamycin from a bound to an unbound state could be plotted in the reference frame of the protein. Next, the aligned conformations were sorted with respect to the umbrella windows, which were sampled for 1 ns and thus contained a total of 100 conformations. For each window, the CoG of all 100 rapamycin molecules were computed and the geometric center of this ensemble (i.e., the window-center) was computed. The aligned conformation for which rapamycin was the nearest to this geometric center was chosen as the representative conformation for that umbrella window. All such umbrella window representative conformations were determined and their CoG, represented by colored spheres, used for tracing the physical path of dissociation. Finally, for each run, a black curve was drawn manually as a guide for the eye evidencing the separation process. Paths were only traced up to = 10 Å, i.e., for the first 100 umbrella windows, so as to provide a clear vision of the initial steps of the dissociation reaction.

Results

Binding Free Energy Computations

All 10 ns-long unbiased MD simulations prior to CPUS showed root-mean-square deviations (RMSDs) within 1.0 Å and root-mean-square fluctuations (RMSFs) in the 0.4–3.6 Å range, which indicates the absence of any large conformational transition (Supplementary Figure S6). Thus, the FKBP12 variants in complexes with rapamycin are stable while well equilibrated. If this result is not surprising for the FKBP12_WT•rapamycin complex, it also validates our simple modeling of FKBP12_Y82F•rapamycin and FKBP12_D37V•rapamycin.

Next, for each variant, 15 independent CPUS MD simulations were conducted and the corresponding PMF profiles were constructed. All curves flatten before reaching a CoG separation distance of 30.0 Å, which is a clear signature of complete dissociations (Figures 2A–C). The binding free energies could thus be identified as the last obtained PMF values. Additionally, in all cases, we could determine a cutoff separation distance at which the van der Waals interactions had decreased to zero (Supplementary Figure S7): roughly ranges between 13 and 19 Å (Supplementary Table S1). With these cutoffs in hand we could further correct the according to Eq. 1, which yielded the standard free energies (Supplementary Table S1). In application of the variational principle for each of the three complexes, we retained the lower-bound PMF profile for comparison with experimental data (Table 1, Figure 2D, Supplementary Figure S8, and Supplementary Table S2).

We first compared our simulations results with the data contained in the sole article we found that experimentally evaluated, with the same technique, the binding of rapamycin to all of the three FKBP12 variants hereby considered (Table 1) (DeCenzo et al., 1996b). With our CPUS approach, the differences of binding free energies between the mutant proteins and the wild-type one are = 0.76 and 3.20 kcal/mol for Y82F and D37V, respectively. With the enzymatic inhibition assay reported in the literature, we have = 0.55 ± 0.33 and 2.92 ± 0.25 kcal/mol. Thus, both datasets are consistent if we consider uncertainties: We could predict in silico that Y82F is a nearly neutral mutation whereas D37V significantly weakens the interaction with rapamycin. Interestingly, the thermodynamic integration technique applied to the FKBP12•FK506 complex could also account for the very small destabilizing effect of Y82F, in similar agreement with experimental measurements (Pearlman and Connelly, 1995).

In a second step, we aimed at testing if our CPUS strategy can also provide reliable estimates of binding free energies (and not only of their differences). Therefore, we collected dissociation equilibrium constants for the FKBP12•rapamycin complex from a little less than 20 original research articles, spanning more than a 30 years-long period and describing measurements ranging from enzymatic inhibition assays to single-molecule ones, from isothermal calorimetry to surface plasmon resonance (Supplementary Table S2). The experimental data were first temperature corrected using published binding enthalpies (Supplementary Figure S8) (Connelly et al., 1993; Connelly et al., 1994) and then gathered in Figure 2D. The results that could be anticipated from Table 1 were confirmed: The in silico stability seems to be ≈1 kcal/mol higher than the one observed at the bench–although the reported values for FKBP12_WT are spread over more than 3 kcal/mol, the less negative binding free energies are possibly artifactual. Since the same discrepancy between simulations and measurements is observed for all three variants, the chances are high that it is a systematic effect. It could be due to the force-field inaccuracy and/or to the inability to fully account for the experimental system.

Atom-Pair Distance Analysis

Even though both Y82 and D37 residues form a hydrogen bond with rapamycin, their role in ligand binding seems to differ: Removing the hydroxyl from the tyrosine has nearly no impact on the affinity (Connelly et al., 1994; Pearlman and Connelly, 1995; Kostrz et al., 2019) whereas exchanging the ethanoate side chain of the aspartate for an isopropyl one is significantly destabilizing (DeCenzo et al., 1996a; Singh et al., 2015; Kostrz et al., 2019). Since CPUS simulations could reproduce these energetic measurements, we next tried to see if the obtained trajectories could shed some light on the molecular mechanisms at work. More precisely, we performed a detailed analysis of both hydrogen bonded atom-pairs in the FKBP12_WT•rapamycin complex, during the starting 10 ns-long unbiased simulations and during the following 190 ns-long CPUS ones (Figures 3, 4; Supplementary Figures S9, S10). Furthermore, the resulting patterns were compared with the ones obtained on the two mutant complexes, for atoms located at equivalent positions on the side chains.

FIGURE 3

FIGURE 4

During the first 50 ns simulating the wild-type complex (10 ns unbiased and 40 ns of CPUS), the length of the H-bond between Y82 and rapamycin remains unchanged, with ground level values around 3 Å for all 15 simulations. This fact is illustrated by the single position of the main peaks in the atom-pair distance histograms of Figure 3A and Supplementary Figure S9A (to illustrate our description different markers have been posted on the time-trace and on the histogram associated with the lower-bound PMF profile). Then, step-increases are observed and the atom-pair distances stabilize on plateaus located between 5 and 10 Å, depending on the simulation run. It thus yields secondary peaks with narrowly dispersed positions. As pulling progresses, these rapid transitions are sometimes followed by others, reflecting temporary interactions that last from ten to tens of ns and atom-pair distances that correspond to plateaus with higher average values. Such an evolution is evidenced by additional peaks in the histograms, peaks that are broader and that now lie in the 10 to 25 Å range. Incidentally, temporary back motions to a previous position are also possible, as exemplified by the R01 time-trace in Figure 3A. As one can see in Figure 3B and Supplementary Figure S9B, a similar behavior is observed for the carbon atom at the tip of the F82 phenyl ring in the FKBP12_Y82F•rapamycin complex; the only differences are that the first step-increases occur roughly 5–10 ns earlier and that the starting atom-pair distance is 1–2 Å shifted (because we are one more covalent bond apart from the oxygen of the ligand). Such resembling stepwise dissociations for both the wild-type and the mutant proteins indicate that the aromaticity of the side-chain may play a role here. To probe this assumption, we plotted the atom-pair distance for all six carbons of the ring (Figures 3C,D and Supplementary Figure S9) and indeed evidenced a concerted motion making us think that a persistent, long-range pseudo-bond could be relevant to anion-π interactions between the side-chain benzyl ring and the oxygen atoms of the rapamycin (Schottel et al., 2008; Chifotides and Dunbar, 2013; Giese et al., 2016; Lucas et al., 2016).

The same detailed approach was subsequently applied to study the role of D37 in rapamycin binding. In the wild-type complex, the length of the H-bond between this residue and the ligand remains nearly unchanged during the first 50–70 ns of simulation: It just fluctuates between 2.7 and 4.0 Å (Figure 4A and Supplementary Figures S10A). The origin of this fluctuation is clearly revealed when considering the other oxygen of the aspartate side-chain (Figure 4B and Supplementary Figures S10A): The two heteroatoms alternatively interact with rapamycin upon rotation of the carboxylate group, which results in two sharp peaks in the atom-pair distance histograms. Once the H-bond is broken, the oxygens of D37 and the one of the macrolide are torn in a continuous manner, showing very little step-like distance increments. Moreover, if a plateau is observed, it manifests itself in the histogram through a low and shallow peak, i.e., a somehow short and weak interaction. Comparatively, in the FKBP12_D37V•rapamycin complex the methyl groups of valine do not display any switching (Figure 4C and Supplementary Figures S10B). Furthermore, all peaks in the histograms are broader than in the wild-type case, the one corresponding to the complex at equilibrium as well as the ones corresponding to transient interactions occurring upon dissociation. These two observations are in line with the inability of the alkyl moieties to form strong, directional intermolecular bonds. Finally, the most obvious signature of the weakening of the association between rapamycin and FKBP12 after the D37V substitution is that the atom-pair distances start to increase earlier, nearly as soon as the bias is applied.

Dissociation Path Tracing

The umbrella sampling simulation trajectories were processed to trace the physical paths of dissociation for the FKBP12_WT•rapamycin, FKBP12_Y82F•rapamycin, and FKBP12_D37V•rapamycin complexes (Supplementary Figures S11–S13, respectively). Rapamycin is always seen to traverse to the solvent in a curvilinear way, exiting through either the top left or the top right of the binding cavity. For the paths corresponding to the lower-bound PMF profiles, the macrolide initially slides on the surface of the protein smoothly and then wind a lot. In contrast, the paths associated with the highest PMF profiles are relatively less curvy (Supplementary Figure S14). The highly curved nature of lower-bound paths suggests that rapamycin does not get stuck on FKBP12 and tends to free from its surface towards the solvent quite early, or at least it does not get trapped into local energy minima. On the other hand, the highest PMF paths for FKBP12_WT and FKBP12_Y82F display accumulations of spheres that could be the signature of some interactions outside of the binding site.

Discussion and Conclusions

In this work, we present a curvilinear-path umbrella sampling simulation approach to estimating the free energy of binding for the FKBP12•rapamycin complex. The strategy consists of several crucial steps that includes sampling enhancement along multiple independent curvilinear paths, construction of PMFs, determination of the lower-bound profile and of the corresponding term, and then correction of this term to obtain the standard free energy of binding. The corresponding theoretical framework was recently developed and successfully implemented to estimate protein–protein interaction energetics (Joshi and Lin, 2019). Here, we extended the approach to a protein–ligand system, and our results are in good agreement with reported experimental data. Although overall the CPUS approach is quite effective to estimate reliably interaction energetics, there are a few things one needs to keep in mind before computing the correction term. The umbrella implementation needs to decide the distance range, i.e., starting reference bound distance and the end separation distance. Upon constructing PMF, the profile may depict a signature of complete dissociation in terms of converging of the curve to some value. However, in the computation of the correction term, we suggested interaction energetics (VDW component) criterion to judge the complete dissociation (see the Methods and Materials section). The simulations should be run to that separation distance where the VDW component converges to zero. In the case of protein–small ligand dissociation reaction, this distance may vary a bit depending upon the path that the ligand could take. Therefore, while setting the end separation distance, it is advised to set it to sufficiently large values, considering that the dissociating compound may slide/rebind on the protein surface after an initial dissociation.

Although Lee and Olson mentioned the necessity of sampling along a curvilinear physical trajectory, they estimated the absolute binding affinities using a vectorial path (Lee and Olson, 2006). Furthermore, their energetic results were in good agreement with the experiments but quadratic approximations were required. Doudou et al. (2009) discussed rigorously the problem associated with energetics obtained by sampling along a predefined vectorial reaction coordinate and pointed out related inconsistencies in the free energy estimation along different linear paths/directions. The CPUS approach is a naïve approach that tries to address these issues by conducting successive umbrella sampling simulations. The approach is applicable especially when a precise input, such as the detailed shape of the binding pocket or a predefined vector of confinement, is unavailable. It relies on the minimal use of restraint potentials (namely, only one biasing potential in the umbrella sampling simulation); hence, it fairly bypasses some major challenges one can face for appropriate de-biasing. However, the cost has to be paid in terms of conducting a sufficiently large number of umbrella sampling simulations to find the lower-bound of the PMF profiles. In addition, due to the complex nature of biomolecular systems it is rather difficult, if not impossible, to decide a priori the number of multiple-walker umbrella sampling simulations to be conducted. This difficulty is actually also shared by any variational-based approaches, e.g., quantum Monte Carlo (Ceperley, 1978; Ceperley and Alder, 1980). It would be always a good sign to find some robust lower-bound from the multiple PMF profiles. Thus, the current work delivers an important message that the PMF profile from any single umbrella sampling simulation (irrespective of whether predefined vector-based or non-predefined curved-path-based sampling enhancement) may likely be misleading. With ongoing advances in GPU architecture, the simulation time is expected to keep decreasing, enabling more simulations to be performed in parallel. Indeed, the implicit/continuum solvent model would rather be a more natural choice for rapid estimation and reducing the computation cost. However, some serious concerns are yet to be suitably resolved, such as too large numerical ranges of estimated energies and strong dependence of the energetics on the employed continuum solvation (such as the choice of protein dielectric constant, the definition of protein boundary, etc.) (Genheden and Ryde, 2015). Therefore, explicit solvent modes are highly preferred.

A related issue in the CPUS approach is to evaluate if the molecular behaviors associated with the lower-bound PMF profile significantly differ from the ones associated with the higher profiles. In our previous work, we traced the protein–protein dissociation paths and observed that there was a clear difference in the direction of traverse between the lower-bound case and the rest (Joshi and Lin, 2019). This distinction between the paths was made possible because the two prominent paths had been previously evidenced using an extensive brute force adaptive MD simulation approach (Plattner et al., 2017). In the case of the FKBP12•rapamycin system, we observed that the rapamycin initially moves up and is then released in the solvent by taking either a left or a right path (Supplementary Figures S11–S14). However, we could not find any supportive correlation between the plateau value of the PMF profile and the path taken. Similarly, we looked for correlations between the PMF final value and the presence of sphere clusters along the way to dissociation, the latter signature being interpreted as the temporary immobilization of rapamycin in a local energy minimum present at the surface of FKBP12. Although few of such potential wells could be identified (see, for instance, the red circles in Supplementary Figure S14), these finding are at the moment only qualitative. As a consequence, the physical paths of dissociation do not enable us to differentiate one or a subset of simulations from the other ones. In fact, it would be interesting to see whether the paths traced from the CPUS approach are consistent with the ones revealed by other pathway sampling methods, such as the string method (Weinan et al., 2002), the minimum free energy path (Maragliano et al., 2006), the adaptive and unconstrained enhanced sampling (Miao and McCammon, 2016), the nudged elastic band (NEB) method (Jónsson et al., 1998), or stochastic difference equation (in length version) SDEL (Arora and Schlick, 2005).

In this work, we also performed a detailed analysis on the atom-pair distance distribution patterns (Figures 3, 4 and Supplementary Figures S9, S10). For all three variants the data that yielded the lower-bound PMF profile and the ones issued from the 14 other runs were rather similar in terms of global behavior, i.e., notwithstanding the variability one can expect when looking at single molecules. Since the ligand conformation may also play a crucial role in the dissociation mechanism, we plotted the temporal evolution of the RMSD for the heavy atoms (C, N, O) of rapamycin (Supplementary Figures S15–S17). Although some peaks or steps are present, the overall RMSD profiles are quite stable for all CPUS simulations and the dissociation process does not seem to be coupled with important conformational transitions in the macrolide. The two hydrogen-bonded atom-pairs investigated here, i.e., Y(F)82-O3 and residue D(V)37-O6, have significantly different patterns. Moreover, the comparison of distance distribution patters between atom-pairs of the side-chain benzyl ring of FKBP12 variants (Y82 and F82) and the O6 of rapamycin pointed out that they might be involved in the anion-π type of interactions. These atomistic-level factors could contribute to the excess dissipation along the curvilinear trajectory of dissociation, and further detailed investigations are needed to have a comprehensive understanding of the relationship between estimated energetics and atom-pair interactions.

Thus, the CPUS approach paves the way to reveal binding energetics along with mechanistic details. The approach is highly generalized and can be implemented to almost any protein–ligand systems, along with protein–protein systems, as well. The approach offers a strong platform to perform several types of in-depth analysis towards revealing the underlying mechanistic details. The drug discovery process often deals with known drug targets and with new or modified drug molecules. In the case of rapamycin-FRB interactions, one such study recently identified DL001 compound that reduces the side effects in vivo (Schreiber et al., 2019). The CPUS approach-based characterization could be helpful for further investigation and for the design of such effective drugs.

References

• 1

  AroraK.SchlickT. (2005). Conformational Transition Pathway of Polymerase β/DNA upon Binding Correct Incoming Substrate. J. Phys. Chem. B109, 5358–5367. 10.1021/jp0446377

  • CrossRef
  • Google Scholar

• 2

  BanaszynskiL. A.LiuC. W.WandlessT. J. (2005). Characterization of the FKBP·Rapamycin·FRB Ternary Complex. J. Am. Chem. Soc.127, 4715–4721. 10.1021/ja043277y

  • CrossRef
  • Google Scholar

• 3

  BiererB. E.MattilaP. S.StandaertR. F.HerzenbergL. A.BurakoffS. J.CrabtreeG.et al (1990). Two Distinct Signal Transmission Pathways in T Lymphocytes Are Inhibited by Complexes Formed between an Immunophilin and Either FK506 or Rapamycin. Proc. Natl. Acad. Sci. U.S.A.87, 9231–9235. 10.1073/pnas.87.23.9231

  • CrossRef
  • Google Scholar

• 4

  BossardM. J.BergsmaD. J.BrandtM.LiviG. P.EngW. K.JohnsonR. K.et al (1994). Catalytic and Ligand Binding Properties of the FK506 Binding Protein FKBP12: Effects of the Single Amino Acid Substitution of Tyr82 to Leu. Biochem. J.297, 365–372. 10.1042/bj2970365

  • CrossRef
  • Google Scholar

• 5

  CaseD. A.BelfonK.Ben-ShalomI. Y.BrozellS. R.CeruttiD. S.CheathamT. E.Iiiet al (2020). AMBER 2020. San Francisco: University of California.

  • Google Scholar

• 6

  CaseD. A.CheathamT. E.IiiDardenT.GohlkeH.LuoR.MerzK. M.Jret al (2005). The Amber Biomolecular Simulation Programs. J. Comput. Chem.26, 1668–1688. 10.1002/jcc.20290

  • CrossRef
  • Google Scholar

• 7

  CeperleyD. (1978). Ground State of the Fermion One-Component Plasma: A Monte Carlo Study in Two and Three Dimensions. Phys. Rev. B18, 3126–3138. 10.1103/physrevb.18.3126

  • CrossRef
  • Google Scholar

• 8

  CeperleyD. M.AlderB. J. (1980). Ground State of the Electron Gas by a Stochastic Method. Phys. Rev. Lett.45, 566–569. 10.1103/physrevlett.45.566

  • CrossRef
  • Google Scholar

• 9

  ChifotidesH. T.DunbarK. R. (2013). Anion−π Interactions in Supramolecular Architectures. Acc. Chem. Res.46, 894–906. 10.1021/ar300251k

  • CrossRef
  • Google Scholar

• 10

  ChoiJ.ChenJ.SchreiberS. L.ClardyJ. (1996). Structure of the FKBP12-Rapamycin Complex Interacting with Binding Domain of Human FRAP. Science273, 239–242. 10.1126/science.273.5272.239

  • CrossRef
  • Google Scholar

• 11

  ConnellyP. R.AldapeR. A.BruzzeseF. J.ChambersS. P.FitzgibbonM. J.FlemingM. A.et al (1994). Enthalpy of Hydrogen Bond Formation in Aprotein-Ligand Binding Reaction. Proc. Natl. Acad. Sci. U.S.A.91, 1964–1968. 10.1073/pnas.91.5.1964

  • CrossRef
  • Google Scholar

• 12

  ConnellyP. R.ThomsonJ. A.FitzgibbonM. J.BruzzeseF. J. (1993). Probing Hydration Contributions to the Thermodynamics of Ligand Binding by Proteins. Enthalpy and Heat Capacity Changes of Tacrolimus and Rapamycin Binding to FK506 Binding Protein in Deuterium Oxide and Water. Biochemistry32, 5583–5590. 10.1021/bi00072a013

  • CrossRef
  • Google Scholar

• 13

  DecenzoM. T.ParkS. T.JarrettB. P.AldapeR. A.FuterO.MurckoM. A.et al (1996a). FK506-binding Protein Mutational Analysis: Defining the Active-Site Residue Contributions to Catalysis and the Stability of Ligand Complexes. Protein Eng. Des. Sel.9, 173–180. 10.1093/protein/9.2.173

  • CrossRef
  • Google Scholar

• 14

  DecenzoM. T.ParkS. T.JarrettB. P.AldapeR. A.FuterO.MurckoM. A.et al (1996b). FK506-binding Protein Mutational Analysis: Defining the Active-Site Residue Contributions to Catalysis and the Stability of Ligand Complexes. Protein Eng. Des. Sel.9, 173–180. 10.1093/protein/9.2.173

  • CrossRef
  • Google Scholar

• 15

  DickmanD. A.DingH.LiQ.NiliusA. M.BalliD. J.BallaronS. J.et al (2000). Antifungal Rapamycin Analogues with Reduced Immunosuppressive Activity. Bioorg. Med. Chem. Lett.10, 1405–1408. 10.1016/s0960-894x(00)00184-0

  • CrossRef
  • Google Scholar

• 16

  DicksonA.LotzS. D. (2016). Ligand Release Pathways Obtained with WExplore: Residence Times and Mechanisms. J. Phys. Chem. B120, 5377–5385. 10.1021/acs.jpcb.6b04012

  • CrossRef
  • Google Scholar

• 17

  DoudouS.BurtonN. A.HenchmanR. H. (2009). Standard Free Energy of Binding from a One-Dimensional Potential of Mean Force. J. Chem. Theory Comput.5, 909–918. 10.1021/ct8002354

  • CrossRef
  • Google Scholar

• 18

  FujitaniH.TanidaY.ItoM.JayachandranG.SnowC. D.ShirtsM. R.et al (2005). Direct Calculation of the Binding Free Energies of FKBP Ligands. J. Chem. Phys.123, 084108. 10.1063/1.1999637

  • CrossRef
  • Google Scholar

• 19

  GalatA. (2013). Functional Diversity and Pharmacological Profiles of the FKBPs and Their Complexes with Small Natural Ligands. Cell. Mol. Life Sci.70, 3243–3275. 10.1007/s00018-012-1206-z

  • CrossRef
  • Google Scholar

• 20

  GenhedenS.RydeU. (2015). The MM/PBSA and MM/GBSA Methods to Estimate Ligand-Binding Affinities. Expert Opin. drug Discov.10, 449–461. 10.1517/17460441.2015.1032936

  • CrossRef
  • Google Scholar

• 21

  GieseM.AlbrechtM.RissanenK. (2016). Experimental Investigation of Anion-π Interactions - Applications and Biochemical Relevance. Chem. Commun.52, 1778–1795. 10.1039/c5cc09072e

  • CrossRef
  • Google Scholar

• 22

  GrazianiF.AldegheriL.TerstappenG. C. (1999). High Throughput Scintillation Proximity Assay for the Identification of FKBP-12 Ligands. SLAS Discov.4, 3–7. 10.1177/108705719900400102

  • CrossRef
  • Google Scholar

• 23

  GrossfieldA. (2013). WHAM: the Weighted Histogram Analysis Method, Version 2.0. 9. Available at membrane. urmc. rochester. edu/content/wham. Accessed November 15, 2013.

  • Google Scholar

• 24

  HamiltonG. S.SteinerJ. P. (1998). Immunophilins: Beyond Immunosuppression. J. Med. Chem.41, 5119–5143. 10.1021/jm980307x

  • CrossRef
  • Google Scholar

• 25

  HoltD. A.LuengoJ. I.YamashitaD. S.OhH. J.KonialianA. L.YenH. K.et al (1993). Design, Synthesis, and Kinetic Evaluation of High-Affinity FKBP Ligands and the X-Ray Crystal Structures of Their Complexes with FKBP12. J. Am. Chem. Soc.115, 9925–9938. 10.1021/ja00075a008

  • CrossRef
  • Google Scholar

• 26

  IkuraT.ItoN. (2007). Requirements for Peptidyl-Prolyl Isomerization Activity: A Comprehensive Mutational Analysis of the Substrate-Binding Cavity of FK506-Binding Protein 12. Protein Sci.16, 2618–2625. 10.1110/ps.073203707

  • CrossRef
  • Google Scholar

• 27

  JónssonH.MillsG.JacobsenK. W. (1998). Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions.

  • Google Scholar

• 28

  JoshiD. C.LinJ. H. (2019). Delineating Protein-Protein Curvilinear Dissociation Pathways and Energetics with Naïve Multiple‐Walker Umbrella Sampling Simulations. J. Comput. Chem.40, 1652–1663. 10.1002/jcc.25821

  • CrossRef
  • Google Scholar

• 29

  KarplusM.MccammonJ. A. (2002). Molecular Dynamics Simulations of Biomolecules. Nat. Struct. Biol.9, 646–652. 10.1038/nsb0902-646

  • CrossRef
  • Google Scholar

• 30

  KästnerJ. (2011). Umbrella Sampling. Wiley Interdiscip. Rev. Comput. Mol. Sci.1, 932–942. 10.1002/wcms.66

  • CrossRef
  • Google Scholar

• 31

  KostrzD.Wayment-SteeleH. K.WangJ. L.FollenfantM.PandeV. S.StrickT. R.et al (2019). A Modular DNA Scaffold to Study Protein-Protein Interactions at Single-Molecule Resolution. Nat. Nanotechnol.14, 988–993. 10.1038/s41565-019-0542-7

  • CrossRef
  • Google Scholar

• 32

  KozanyC.MärzA.KressC.HauschF. (2009). Fluorescent Probes to Characterise FK506-Binding Proteins. Chembiochem10, 1402–1410. 10.1002/cbic.200800806

  • CrossRef
  • Google Scholar

• 33

  KumarS.RosenbergJ. M.BouzidaD.SwendsenR. H.KollmanP. A. (1992). THE Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem.13, 1011–1021. 10.1002/jcc.540130812

  • CrossRef
  • Google Scholar

• 34

  LeeM. S.OlsonM. A. (2006). Calculation of Absolute Protein-Ligand Binding Affinity Using Path and Endpoint Approaches. Biophysical J.90, 864–877. 10.1529/biophysj.105.071589

  • CrossRef
  • Google Scholar

• 35

  LinJ.-H.PerrymanA. L.SchamesJ. R.MccammonJ. A. (2002). Computational Drug Design Accommodating Receptor Flexibility: the Relaxed Complex Scheme. J. Am. Chem. Soc.124, 5632–5633. 10.1021/ja0260162

  • CrossRef
  • Google Scholar

• 36

  LinJ.-H.PerrymanA. L.SchamesJ. R.MccammonJ. A. (2003). The Relaxed Complex Method: Accommodating Receptor Flexibility for Drug Design with an Improved Scoring Scheme. Biopolymers68, 47–62. 10.1002/bip.10218

  • CrossRef
  • Google Scholar

• 37

  LuC.WangZ.-X. (2017). Quantitative Analysis of Ligand Induced Heterodimerization of Two Distinct Receptors. Anal. Chem.89, 6926–6930. 10.1021/acs.analchem.7b01274

  • CrossRef
  • Google Scholar

• 38

  LucasX.BauzáA.FronteraA.QuiñoneroD. (2016). A Thorough Anion-π Interaction Study in Biomolecules: on the Importance of Cooperativity Effects. Chem. Sci.7, 1038–1050. 10.1039/c5sc01386k

  • CrossRef
  • Google Scholar

• 39

  LuengoJ. I.YamashitaD. S.DunningtonD.BeckA. K.RozamusL. W.YenH.-K.et al (1995). Structure-activity Studies of Rapamycin Analogs: Evidence that the C-7 Methoxy Group Is Part of the Effector Domain and Positioned at the FKBP12-FRAP Interface. Chem. Biol.2, 471–481. 10.1016/1074-5521(95)90264-3

  • CrossRef
  • Google Scholar

• 40

  MaierJ. A.MartinezC.KasavajhalaK.WickstromL.HauserK. E.SimmerlingC. (2015). ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput.11, 3696–3713. 10.1021/acs.jctc.5b00255

  • CrossRef
  • Google Scholar

• 41

  MaraglianoL.FischerA.Vanden-EijndenE.CiccottiG. (2006). String Method in Collective Variables: Minimum Free Energy Paths and Isocommittor Surfaces. J. Chem. Phys.125, 024106. 10.1063/1.2212942

  • CrossRef
  • Google Scholar

• 42

  MiaoY.MccammonJ. A. (2016). Unconstrained Enhanced Sampling for Free Energy Calculations of Biomolecules: a Review. Mol. Simul.42, 1046–1055. 10.1080/08927022.2015.1121541

  • CrossRef
  • Google Scholar

• 43

  MobleyD. L.GilsonM. K. (2017). Predicting Binding Free Energies: Frontiers and Benchmarks. Annu. Rev. Biophys.46, 531–558. 10.1146/annurev-biophys-070816-033654

  • CrossRef
  • Google Scholar

• 44

  NerattiniF.ChelliR.ProcacciP. (2016). II. Dissociation Free Energies in Drug-Receptor Systems via Nonequilibrium Alchemical Simulations: Application to the FK506-Related Immunophilin Ligands. Phys. Chem. Chem. Phys.18, 15005–15018. 10.1039/c5cp05521k

  • CrossRef
  • Google Scholar

• 45

  OlivieriL.GardebienF. (2011). Molecular Dynamics Simulations of a Binding Intermediate between FKBP12 and a High-Affinity Ligand. J. Chem. Theory Comput.7, 725–741. 10.1021/ct100394d

  • CrossRef
  • Google Scholar

• 46

  PearlmanD. A.ConnellyP. R. (1995). Determination of the Differential Effects of Hydrogen Bonding and Water Release on the Binding of FK506 to Native and Tyr82→Phe82 FKBP-12 Proteins Using Free Energy Simulations. J. Mol. Biol.248, 696–717. 10.1006/jmbi.1995.0252

  • CrossRef
  • Google Scholar

• 47

  PettersenE. F.GoddardT. D.HuangC. C.CouchG. S.GreenblattD. M.MengE. C.et al (2004). UCSF Chimera?A Visualization System for Exploratory Research and Analysis. J. Comput. Chem.25, 1605–1612. 10.1002/jcc.20084

  • CrossRef
  • Google Scholar

• 48

  PlattnerN.DoerrS.De FabritiisG.NoéF. (2017). Complete Protein-Protein Association Kinetics in Atomic Detail Revealed by Molecular Dynamics Simulations and Markov Modelling. Nat. Chem.9, 1005–1011. 10.1038/nchem.2785

  • CrossRef
  • Google Scholar

• 49

  RoeD. R.CheathamT. E.Iii (2013). PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput.9, 3084–3095. 10.1021/ct400341p

  • CrossRef
  • Google Scholar

• 50

  Salomon‐FerrerR.CaseD. A.WalkerR. C. (2013). An Overview of the Amber Biomolecular Simulation Package. Wiley Interdiscip. Rev. Comput. Mol. Sci.3, 198–210. 10.1002/wcms.1121

  • CrossRef
  • Google Scholar

• 51

  SapienzaP. J.MauldinR. V.LeeA. L. (2011). Multi-Timescale Dynamics Study of FKBP12 along the Rapamycin-mTOR Binding Coordinate. J. Mol. Biol.405, 378–394. 10.1016/j.jmb.2010.10.037

  • CrossRef
  • Google Scholar

• 52

  SchottelB. L.ChifotidesH. T.DunbarK. R. (2008). Anion-π Interactions. Chem. Soc. Rev.37, 68–83. 10.1039/b614208g

  • CrossRef
  • Google Scholar

• 53

  SchreiberK. H.Arriola ApeloS. I.YuD.BrinkmanJ. A.VelardeM. C.SyedF. A.et al (2019). A Novel Rapamycin Analog Is Highly Selective for mTORC1 In Vivo. Nat. Commun.10, 3194. 10.1038/s41467-019-11174-0

  • CrossRef
  • Google Scholar

• 54

  SchulerW.SedraniR.CottensS.H?berlinB.SchulzM.SchuurmanH.-J.et al (1997). Sdz Rad, a New Rapamycin Derivative. Transplantation64, 36–42. 10.1097/00007890-199707150-00008

  • CrossRef
  • Google Scholar

• 55

  ShorB.ZhangW.-G.Toral-BarzaL.LucasJ.AbrahamR. T.GibbonsJ. J.et al (2008). A New Pharmacologic Action of CCI-779 Involves FKBP12-independent Inhibition of mTOR Kinase Activity and Profound Repression of Global Protein Synthesis. Cancer Res.68, 2934–2943. 10.1158/0008-5472.can-07-6487

  • CrossRef
  • Google Scholar

• 56

  SinghV.NandA.SaritaS. (2015). Universal Screening Platform Using Three-Dimensional Small Molecule Microarray Based on Surface Plasmon Resonance Imaging. RSC Adv.5, 87259–87265. 10.1039/c5ra15637h

  • CrossRef
  • Google Scholar

• 57

  SolomentsevG.DiehlC.AkkeM. (2018). Conformational Entropy of Fk506 Binding to Fkbp12 Determined by Nuclear Magnetic Resonance Relaxation and Molecular Dynamics Simulations. Biochemistry57, 1451–1461. 10.1021/acs.biochem.7b01256

  • CrossRef
  • Google Scholar

• 58

  SunF.LiP.DingY.WangL.BartlamM.ShuC.et al (2003). Design and Structure-Based Study of New Potential FKBP12 Inhibitors. Biophysical J.85, 3194–3201. 10.1016/s0006-3495(03)74737-7

  • CrossRef
  • Google Scholar

• 59

  SwansonJ. M. J.HenchmanR. H.MccammonJ. A. (2004). Revisiting Free Energy Calculations: a Theoretical Connection to MM/PBSA and Direct Calculation of the Association Free Energy. Biophysical J.86, 67–74. 10.1016/s0006-3495(04)74084-9

  • CrossRef
  • Google Scholar

• 60

  TamuraT.KioiY.MikiT.TsukijiS.HamachiI. (2013). Fluorophore Labeling of Native FKBP12 by Ligand-Directed Tosyl Chemistry Allows Detection of its Molecular Interactions In Vitro and in Living Cells. J. Am. Chem. Soc.135, 6782–6785. 10.1021/ja401956b

  • CrossRef
  • Google Scholar

• 61

  TorrieG. M.ValleauJ. P. (1977). Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys.23, 187–199. 10.1016/0021-9991(77)90121-8

  • CrossRef
  • Google Scholar

• 62

  Van DuyneG. D.StandaertR. F.KarplusP. A.SchreiberS. L.ClardyJ. (1993). Atomic Structures of the Human Immunophilin FKBP-12 Complexes with FK506 and Rapamycin. J. Mol. Biol.229, 105–124. 10.1006/jmbi.1993.1012

  • CrossRef
  • Google Scholar

• 63

  Van DuyneG. D.StandaertR. F.SchreiberS. L.ClardyJ. (1991). Atomic Structure of the Rapamycin Human Immunophilin FKBP-12 Complex. J. Am. Chem. Soc.113, 7433–7434. 10.1021/ja00019a057

  • CrossRef
  • Google Scholar

• 64

  WagnerR.RhoadesT. A.OrY. S.LaneB. C.HsiehG.MollisonK. W.et al (1998). 32-Ascomycinyloxyacetic Acid Derived Immunosuppressants. Independence of Immunophilin Binding and Immunosuppressive Potency. J. Med. Chem.41, 1764–1776. 10.1021/jm960066y

  • CrossRef
  • Google Scholar

• 65

  WallaceA. C.LaskowskiR. A.ThorntonJ. M. (1995). LIGPLOT: a Program to Generate Schematic Diagrams of Protein-Ligand Interactions. Protein Eng. Des. Sel.8, 127–134. 10.1093/protein/8.2.127

  • CrossRef
  • Google Scholar

• 66

  WangJ.WolfR. M.CaldwellJ. W.KollmanP. A.CaseD. A. (2004). Development and Testing of a General Amber Force Field. J. Comput. Chem.25, 1157–1174. 10.1002/jcc.20035

  • CrossRef
  • Google Scholar

• 67

  WangY.BarnettS. F. H.LeS.GuoZ.ZhongX.KanchanawongP.et al (2019). Label-free Single-Molecule Quantification of Rapamycin-Induced FKBP-FRB Dimerization for Direct Control of Cellular Mechanotransduction. Nano Lett.19, 7514–7525. 10.1021/acs.nanolett.9b03364

  • CrossRef
  • Google Scholar

• 68

  WearM. A.PattersonA.WalkinshawM. D. (2007). A Kinetically Trapped Intermediate of FK506 Binding Protein Forms In Vitro: Chaperone Machinery Dominates Protein Folding In Vivo. Protein Expr. Purif.51, 80–95. 10.1016/j.pep.2006.06.019

  • CrossRef
  • Google Scholar

• 69

  WearM. A.WalkinshawM. D. (2007). Determination of the Rate Constants for the FK506 Binding Protein/rapamycin Interaction Using Surface Plasmon Resonance: an Alternative Sensor Surface for Ni2+-Nitrilotriacetic Acid Immobilization of His-Tagged Proteins. Anal. Biochem.371, 250–252. 10.1016/j.ab.2007.06.034

  • CrossRef
  • Google Scholar

• 70

  WeinanE.RenW.Vanden-EijndenE. (2002). String Method for the Study of Rare Events. Phys. Rev. B66, 052301. 10.1103/physrevb.66.052301

  • CrossRef
  • Google Scholar

• 71

  WilliamsonM. P. (2018). “Chemical Shift Perturbation,” in Chemical Shift perturbationModern Magnetic Resonance. Editor WebbG.A. (Cham: Springer International Publishing), 995–1012. 10.1007/978-3-319-28388-3_76

  • CrossRef
  • Google Scholar

• 72

  WilsonK. P.YamashitaM. M.SintchakM. D.RotsteinS. H.MurckoM. A.BogerJ.et al (1995). Comparative X-Ray Structures of the Major Binding Protein for the Immunosuppressant FK506 (Tacrolimus) in Unliganded Form and in Complex with FK506 and Rapamycin. Acta Cryst. D.51, 511–521. 10.1107/s0907444994014514

  • CrossRef
  • Google Scholar

• 73

  WuX.WangL.HanY.ReganN.LiP.-K.VillalonaM. A.et al (2011). Creating Diverse Target-Binding Surfaces on FKBP12: Synthesis and Evaluation of a Rapamycin Analogue Library. ACS Comb. Sci.13, 486–495. 10.1021/co200057n

  • CrossRef
  • Google Scholar

• 74

  YangC.-J.TakedaM.TerauchiT.JeeJ.KainoshoM. (2015). Differential Large-Amplitude Breathing Motions in the Interface of FKBP12-Drug Complexes. Biochemistry54, 6983–6995. 10.1021/acs.biochem.5b00820

  • CrossRef
  • Google Scholar