Determining the Relative Binding Affinity of Ricin Toxin A (RTA

Subscriber access provided by University of Pennsylvania Libraries

Computational Chemistry

Determining the Relative Binding Affinity of Ricin Toxin A (RTA) Inhibitors by Using Molecular Docking and Non-Equilibrium Work Elton Chaves, Itácio Queiroz de Mello Padilha, Demetrius Antônio Machado Araújo, and Gerd Bruno Rocha J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00036 • Publication Date (Web): 11 May 2018 Downloaded from http://pubs.acs.org on May 13, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

Determining the Relative Binding Affinity of Ricin Toxin A (RTA) Inhibitors by Using Molecular Docking and Non-Equilibrium Work Elton J. F. Chaves,† Itácio Q. M. Padilha.,† Demétrius A. M. Araújo,†* Gerd B. Rocha,‡* †

Department of Biotechnology, Federal University of Paraíba, João Pessoa, Brazil ‡

Department of Chemistry, Federal University of Paraíba, João Pessoa, Brazil

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT

Ricin is a ribosome-inactivating protein (RIP-type2) consisting of two subunits, Ricin Toxin A (RTA) and Ricin Toxin B (RTB). Due to its cytotoxicity, ricin has worried world authorities for its potential use as a chemical weapon; therefore, its inhibition is of great biotechnological interest. RTA is the target for inhibitor synthesis and pterin derivatives are promising candidates to inhibit it. In this study, we used a combination of molecular docking approach and fast steered molecular dynamics in order to assess the correlation between non-equilibrium work, 〈〉, and IC50 of six RTA inhibitors. The results showed that molecular docking is a powerful tool to predict good bioactive pose of RTA inhibitors and the 〈〉 presented a strong correlation with IC50 (R² = 0.961). Such profile ranked the RTA inhibitors better than molecular docking approach. Therefore, the combination of docking and fast SMD simulation showed to be a promising tool to distinguish RTA active inhibitors from inactive ones, and could be used as post-docking filtering approach.


Page 2 of 43

Page 3 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TABLE OF CONTENTS GRAPHICS



INTRODUCTION Ricin is a cytotoxic protein produced in the seed of Ricinus communis.1,2 This protein has worried world authorities because of its potential use as a chemical weapon. Ricin is composed of two subunits, RTA and RTB with 267 and 262 aminoacids, respectively, connected by a single disulfide bond.3–5 RTA catalyzes a hydrolytic reaction that results in the depurination of adenine 4324 (A4324) located in the GAGA motif at the sarcin–ricin loop (SRL) region of eukaryotic 28S ribosomes. This deletion interrupts protein synthesis, compromising homeostasis of cells, which inevitably die.6 Structurally, the RTA active site is composed of five main residues: Glu177, Arg180, Trp211, Tyr80 and Tyr123. It is known that the natural substrate – adenine – interacts between Tyr80 and Tyr123, and forms hydrogen bonds with Val81, Gly121, Arg180, Glu177 and Glu208.7 RTB is a lectin which internalizes the RTA–RTB complex in the cytosol of target cells through recognition of terminal sugar residues, mainly galactose and N-6-acetylgalactosamine.8 Once internalized, the complex undergoes vesicular retrograde transport to the endoplasmic reticulum (ER).9 In the ER lumen, an isomerase reduces the disulfide bond between A and B subunits, and RTA is translocated to the cytoplasm, attacking ribosomes.10–13⁠ Since the natural substrate of the RTA is the ribosomal RNA (rRNA), some efforts to rational drug design defined rRNA as a template; e.g.: RNA-based inhibitors that bind strongly to RTA.14 Despite their high affinity, these inhibitors have not been used as antidotes due to their low permeability. Thus, it is believed that smaller organic molecules may play a more effective role in this sense.5⁠ In the initial steps of drug discovery, high throughput screening (HTS) and virtual screening (VS) are the most used approaches. HTS method has been used to evaluate the


Page 4 of 43

Page 5 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


protective ability of thousand compounds against ricin intoxication.15 However, the ricin intoxication is quite complex. It involves recognition and capture processes, retrograde vesicular transport to the ER, cytoplasm reallocation, and folding and ribosome inactivation.16 Cell-based screening approaches do not directly identify anti-RTA activities of compounds. Among the inhibitors that were identified by cell-based assays, only a small percentage have showed antiRTA activity in cell-free assays.17 The first small molecules with anti-RTA activity were identified by using VS methodologies such as molecular docking. After that, such molecules were validated by X-ray crystallography and cell-free assays, having the pterin derivatives as the most promising inhibitors.18–25 Molecular docking is a powerful tool to quickly screen large-scale chemical databases for hit identification. However, the protein-ligand binding affinity predictions by this method are uncertain and not efficient to ordering inhibitors according to their binding affinities, having as a contributing factor, the inability to consider receptor flexibility.26⁠ Therefore, there is a possibility of false positive results which, in order to be validated, require rigorous, timeconsuming and expensive experimental steps. In the structure-based methods, the consideration of protein-ligand interactions in a dynamic environment is a key factor in the success of VS protocol.27 Molecular docking in combination with Molecular Dynamics (MD) simulations and free energy calculations, such as Linear Interaction Energy (LIE),28 Molecular Mechanics/Poisson–Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics Generalized Born Surface Area (MM/GBSA),29⁠ are widely used and can predict binding energies with better accuracies. However, their computational cost are higher when compared to molecular docking.27,30



Recently, Steered Molecular Dynamics (SMD)31,32⁠ has been used as a promising tool for drug design.33⁠ This method applies a time-dependent external force on parts of the system, similar to Atomic Force Microscopy (AFM). This force is able to induce dissociation processes, which are not trivially achieved by conventional MD simulations. SMD was first designed for studying biophysical and biochemical processes such as protein unfolding34,35 and transportation of ions.36 However, it may be also used to discern active inhibitors from inactive ones. In this method, the scoring function is either the rupture force 〈 〉 or unbinding work 〈〉 generated by the external force applied to the ligand. Such idea was first reported by Colizzi and co-workers37⁠ and its potential impact was pointed out by Jorgensen.38 Since then, few groups have been used this method to evaluate the unbinding pathway and ligand binding affinities of some therapeutic targets.39–44 Ngo and co-workers40 evaluated the rupture force and pulling work profile to decouple eighteen inhibitors from the active site of HIV-1 protease using fast SMD simulations. Their results pointed out that the non-equilibrium work was better than the rupture force in ranking the relative binding free energy of HIV-1 protease inhibitors. The correlation with experimental binding free energies (∆Gexp) was of R = -0.95.40 Mai and co-workers41⁠ pulled inhibitors out of the binding site of the Influenza virus A/H5N1 neuraminidase using SMD simulations. Their results showed high correlation between rupture force 〈 〉 and experimental data.41 Vuong and co-workers43⁠ pulled the ligands out of binding site of three different target x-ray structures (α-thombin, neuraminidase and penicillopepsin) using SMD simulation. Their results showed strong correlations between nonequilibrium work 〈〉 and ∆Gexp with R < -0.8.43


Page 6 of 43

Page 7 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Singh and co-workers45 used a combination of ligand- and structure-based approaches, considering molecules from the DrugBank database, in order to identify inhibitors to Abl kinase, a very important biological target of the Chronic Myeloid Leukemia (CML). The ligand-based approach in combination with molecular docking identified four compounds. Then, the authors pulled the ligands out of the biding site using SMD simulations and the results for the rupture force, 〈 〉, showed that geftfinib could be a potential allosteric inhibitor of Abl kinase. In vitro validation using imatinib and gefitinib produced synergistic antiproliferative effects in K562 CML cell line.45 In our work, we performed a combination of molecular docking and fast SMD simulations for six RTA inhibitors (pterin derivatives) (Table 1) with binding mode and halfmaximum inhibitory concentration (IC50) determined. Using such strategy, we assessed the correlations of the rupture force 〈 〉 and non-equilibrium work 〈〉 with IC50 to discern active RTA inhibitors from inactive ones.

INSERT TABLE 1



THEORETICAL METHODS Selection of complexes for the test set. We queried the PDBbind-CN46,47⁠ ⁠ database using “ricin” as keyword for all protein–ligand complexes. The database returned 21 complexes; of these, we selected complexes that present a pterin group in the ligand structure and resolution ≤ 2.0 Å. We obtained nine complexes; however, for relative binding energy calculations, we used the following complexes: 4HUP, 4HUO, 4ESI, 4MX1, 3PX8 and 3PX9 (PDB codes) (Table 1). Preparation of the receptor. For RTA, we used the 1BR5 crystallographic structure (PDB ID). Native ligand (PDB code: NEO) and waters were removed, and the cleaned structure was submitted to the PDB2PQR48⁠ web service for residue protonation calculations at pH 7.0. Missing hydrogens and partial charges were added using the default option of the prepare_receptor4.py script, available in AutoDock Tools.49⁠ Preparation of the ligands. We used the following ligands: 19M, RS8, 0RB, 1MX, JP2 and JP3 (PDB code). Missing hydrogen were added using Open Babel v.2.3.2.50⁠ Partial charges and ligand flexibility were defined using the default option of the prepare_ligand4.py script, available in AutoDock Tools.49⁠ Redocking. Bioactive ligand pose predictions were obtained using AutoDock Vina v.1.1.1.51⁠ A box grid was positioned in RTA active site with volume of 20 × 14 × 20 Å3, considering a standard spacing of 0.375 Å. Then, the Vina Lamarckian genetic algorithm was applied 10 times for each ligand. Root mean square deviation (RMSD) measurement between predicted and native poses. We compared the lowest energy poses from redocking step to their respective


Page 8 of 43

Page 9 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


crystallographic poses using RMSD. For this analysis, it was necessary to perform structure alignment and, then, calculate RMSD. We used for that VMD v.1.9.252⁠ and FCONV v.1.24,53 respectively. MD simulations. Complexes with the lowest redocking energies were submitted to MD simulation. In this step, for RTA and ligand parameterizations, we used FF99SB54⁠ and GAFF55 force fields, respectively. Partial charges of the ligands were recalculated according to AM1BCC method by using ambertools package. Then, complexes were inserted into an octahedral box of 20.0 Å containing TIP3P waters. After that, we carried out MD simulations using NAMD v.2.956 having as configuration parameters: (i) periodic boundary conditions; (ii) restriction for vibration of the covalent bond involving hydrogen atoms, HOH angles and OH bond distance of TIP3P water molecules; (ii) time steps equal to 2 femtoseconds; (4) electrostatic interaction cutoff of 9.0Å for all steps of the simulation. The initial geometries of the solvated complexes were optimized and, then, used in heating and pressurizing stages. At the end of pressurization stage, the complexes were conducted to equilibration stage for 10 ns, with constant pressure and temperature of 1.0 atm and 310.0 K, respectively; their coordinates were stored every 1 ps. SMD simulations. The last frame of equilibrated RTA complexes were submitted to SMD simulations. Technically, each ligand was connected to a dummy atom by a spring of constant k. In our study, the spring was attached to a specific ligand atom (labeled as SMDAtom). Then, the ligand was pulled out from the active site following a pulling direction. Moving along the pulling direction with a constant velocity v, the dummy atom experiences the following external force: ext =k vt − ∆x ,

(1)



Page 10 of 43

where Fext is the external force, k is the spring constant, v the velocity, t the time and ∆xi the displacement of the atom i from its initial position (xi – x0). The rupture force and unbinding work are influenced by some parameters such as spring constant (k), pulling velocity (v) and pulling direction (or reaction coordinate).54 As such, we assessed the influence of different SMD settings and their correlations with IC50. Considering Eq. 1, we evaluated the following two scenarios: (1) k = 10 kcal/mol/Å2, v = 0.05, 0.25 and 1.25 Å/ps; (2) k = 2 kcal/mol/Å², v = 0.05, 0.25 and 1.25 Å/ps. In addition, the pulling direction was set using a vector that connects the FixedAtom and SMDAtom. The FixedAtom is the atom that has been kept fixed in a specific position in the protein and the SMDAtom is the atom that we defined previously. The work done by the force, on the isobaric–isothermal ensemble, can be measured using Eq. 2, and linked to the free energy using Jarzynski’s equality (JE)57⁠ (Eq. 3).

W=v ext dt,

(2)

β∆G = βW ,

(3)

where β = 1/kbT, kB is the Boltzmann constant, T is the absolute temperature, W the work given by Eq. 2, ∆G the difference in Gibbs free energy between two thermodynamic states of a bound and unbound system and the brackets … ! illustrate the mean work for infinite times of switch processes between the two thermodynamic states. Here, we did not perform the block-average analysis, thus, we got the free energy difference as ∆" = 〈〉# ,58 considering n = 25. In order to perform SMD simulation, we used NAMD v.2.956⁠ , taking as configuration parameters: (i) harmonic constraints with a constant force of 1.0 kcal/mol/Å2 on the Cα atoms of the RTA to


Page 11 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


limit abrupt movements; (ii) harmonic restraint cylinder of 10 Å diameter, to limit the movement of the ligand only to the axis of reaction coordinate. Except for these parameters, all SMD simulations had the same parameters as MD simulations. Scheme 1 shows a summary of the methodology used in our study.

INSERT SCHEME 1

RESULTS AND DISCUSSION Redocking. Figure 1 shows an illustration of best score poses superposed on their respective crystallographic ones for all studied RTA-ligand complexes. Four ligands had a RMSD ≤ 2.0 Å: 19M, RS8, 1MX and JP2, with RMSD ≈ 1.61, 1.15, 0.353 and 0.301 Å, respectively. 0RB and JP3 inhibitors presented RMSD > 2.0 Å. Inhibitors that had RMSD ≤ 2.0 Å presented total superposition of the pterin group and lateral chain (Figures 1A, 1B, 1D and 1E). In contrast to these, inhibitors that had RMSD > 2.0 Å presented superposition only of the pterin group (Figure 1C and 1F).

INSERT FIGURE 1

Some authors often report a cutoff value of 2.0 Å as the RMSD maximum value for correct pose prediction.59⁠ For the inhibitors 0RB and JP3, we did not obtain an RMSD value



within this criterion, probably due to the inefficiency of the force field used in AutoDock Vina. Nevertheless, the pterin moieties for these ligands were correctly predicted. The lowest energy complexes were submitted to 10 ns of MD simulations and, then, submitted to SMD simulations. After assessing the six complexes in different SMD settings, we noted that setting k = 2 kcal/mol/A² and v = 0.05 Å/ps, we achieved a good agreement with the experimental data. Therefore, in next two sections, we will present only the results considering such SMD setting. Information about results for the SMD force profiles considering different SMD settings are presented in supporting information (See Table S1, Table S2 and Table S3) and discussed in the third section. Mean force profile for decoupling inhibitors from RTA active site. We obtained equilibrated complexes after 6 ns of MD simulation (Supporting Information Figure S1). However, for SMD simulations we used the microstates corresponding to 10 ns (last frame of MD simulations). In addition, we used the Cα of I172 as FixedAtom and the ligand SMDAtoms are depicted in Figure S2. The pulling directions for each RTA-ligand complex are presented in Figure 2.

INSERT FIGURE 2

Figure 3 shows mean force profiles 〈〉 for all six RTA-ligand complexes. The individual force profiles were calculated according to Eq. 1 and their mean force profiles 〈〉 were measured according to the average of the individual forces. There are two parts in the 〈〉


Page 12 of 43

Page 13 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


profiles, before and after of 〈〉 , respectively (Figure 3A). In the part before 〈〉 , the force rapidly increases to 〈〉 , and then decreases to values close to 0.0 pN. Initially, RTA inhibitors are stable in their initial positions – in the RTA active site – but soon after, they are displaced due to the applied external force, which in turn is time-dependent and modulated according to the interactions that the inhibitor makes at the active site. In other words, a barrier to the inhibitor displacement is imposed by RTA–inhibitor interactions, which in turn is overcome by the force increase as a function of time (Figure 3B). That part ends when 〈〉 is reached, which corresponds to the maximum force required to break the interactions between inhibitor and RTA active site; thus, if there is no point of strong interactions – during displacement – the force decreases, which characterizes the second part – after 〈〉 – and the end of the decoupling process. Furthermore, from results presented in Figure 3B, we can note that the most potent inhibitors were pulled out of the RTA binding site slower than the weakest ones (e.g. RS8, 19M and 0RB).

INSERT FIGURE 3

Table 2 shows a summary of results for SMD forces. Regarding the maximum force of mean force profiles, 〈〉 , three complexes presented 〈〉 after 130 ps: RS8–RTA, 19M– RTA and 0RB-RTA with 〈〉 ≈ 884.57, 852.96, and 724.16 pN, respectively. For the same three complexes, we obtained 〈$〉 ≈ 170.82, 157.56, and 130.10 ps, respectively (Table 2). The complexes 1MX-RTA, JP2-RTA and JP3-RTA presented 〈〉 before 130 ps, with 〈〉



≈ 590.55, 499.71 and 309.27 pN, respectively, and 〈$〉 ≈ 129.18, 115.14 and 81.08 ps, respectively. In addition, 〈〉 presented correlation with IC50 of R² = 0.978 (Figure 4). Although 〈〉 shows good correlation with IC50, according to Ngo and coworkers,40⁠ this profile does not represent a relative binding affinity data. According to these authors, 〈〉 presents maximum force peaks at different times in individual simulations; thus, the mean is not calculated from the maximum force data, showing a value smaller than 〈 〉, contradicting rupture force theory.

INSERT TABLE 2

Rupture force and non-equilibrium work profiles rank RTA inhibitors better than molecular docking approach. In Table 2 we also present the rupture force 〈 〉 and nonequilibrium work 〈〉 results and their standard deviations. 〈 〉 was measured using the average of the maximum forces. RS8–RTA complex presented the largest value of 〈 〉, with 〈 〉 ≈ 1056.19 pN. In descending order of 〈 〉 we have the complexes: 19M–RTA, 0RB– RTA, 1MX–RTA, JP2–RTA and JP3–RTA, with 〈 〉 ≈ 979.70, 878.59, 741.15, 649.12 and 449.82 pN, respectively (Table 2). 〈 〉 presented correlation with IC50 of R² = 0.969 (Figure 4). Energy profiles for all RTA inhibitors calculated according to Eq. 3 are showed in Figure 5 (see the Figure S3 of the supporting information for more details about work fluctuation along the twenty-five SMD trajectories). The averaged non-equilibrium work profile presented


Page 14 of 43

Page 15 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


different energy levels for each complex, with values between ~ 31 and 95 kcal/mol. The largest value of 〈〉 was to the complex 19M−RTA, with 〈〉 ≈ 95.01 kcal/mol, followed by RS8–RTA, 0RB–RTA, 1MX–RTA, JP2–RTA and JP3–RTA, with 〈〉 ≈ 91.98, 73.51, 57.34, 48.52 and 31.50 kcal/mol, respectively (Table 2). 〈〉 presented correlation with IC50 of R² = 0.961 (Figure 4). The non-equilibrium work showed a strong association with the half-maximum inhibitory concentration of RTA inhibitors, meaning that a low IC50 corresponds to large non-equilibrium work required to unbinding the inhibitor from its active site. Furthermore, active RTA inhibitors showed relative binding affinities more than 73 kcal/mol. In opposite, inactive inhibitors showed relative binding affinities lower than 53 kcal/mol (Figure 5). This information may be useful for screen more potent RTA inhibitors in a hierarchical virtual screening protocol.

INSERT FIGURE 4 INSERT FIGURE 5

Table 2 also shows inhibitors ranking according all considered estimates of affinities, in descending order of their potencies. Regarding the performance of the estimates that correctly ranked the RTA inhibitors affinities, we can see that: (i) 〈〉 perfectly ranked all RTA inhibitors, including 19M and RS8, whose IC50 values differ only of 5 µM. This result indicates an excellent performance of 〈〉 and suitable selection of the SMD setting, mainly the pulling direction; (ii) 〈 〉 and 〈〉 failed only in ordering 19M and RS8, which are difficult to rank because their



IC50 values are too close. (iii) finally, the docking score function was able to correctly order only JP2 and JP3, which are the worst inhibitors of RTA among the studied ligands, having a poor performance with correlation of R² = 0.383 (Figure 4). In the present study, by using AutoDock Vina, we obtained good predictions of the ligand bioactive conformations (Figure 1). Nevertheless, as expected, the predicted relative binding affinities were not good enough to correctly order the inhibitors (Table 2). Anyway, using the correct starting bound pose of the ligands in undocking procedures, the SMD force profiles showed a good performance in ranking the six inhibitors. In addition, 〈〉 was better than 〈 〉 and 〈〉 in ranking 19M and RS8 inhibitors, because it considered the entire force profile along the time (Eq. 2). SMD force profiles using different spring constants and pulling velocities. When we carried out SMD simulations using spring constant (k) of 10 kcal/mol/A² and v = 0.05 Å/ps, the inhibitors were pulled out faster from the RTA active site (Supporting Information Figure S4). Using these parameters the averaged pulling work was reduced and the correlation with IC50 was R² = 0.954 (Supporting Information Table S3), a R2 slightly lower when compared to the SMD simulations using k = 2 kcal/mol/A² and v = 0.05 Å/ps (R² = 0.961). We noted that the correlations between different pulling velocities and 〈〉 were always larger than 0.9. Furthermore, we also noted that large pulling velocities cause an increase in 〈〉 and in their standard deviations (Figure 6). In addition, high pulling velocities also decrease the correlation between 〈〉 and IC50 (Supporting Information Table S3). Ngo and co-workers40⁠ also found similar correlations when pulled HIV-1 protease inhibitors out using different pulling velocities.


Page 16 of 43

Page 17 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INSERT FIGURE 6

According to Jarzynski’s equality, the free energy should be more accurately predicted when the inhibitor decoupling process occurs slowly, that is, in conditions close to equilibrium.40,43,57⁠ However, performing decoupling simulations at equilibrium demand high computational cost, which is not suitable during the drug discovery process.58 Thus, the relative binding affinities predictions using v = 0.05 Å/ps, despite its good agreement to IC50, are overestimated due to fast pulling speed. Anyway, the results seemingly showed a good linear dependence with the relative binding affinities calculated from equilibrium simulations (Figure 6). Thus, for drug discovery purposes, we suggest that this parameter set are sufficient to discern active RTA inhibitors from inactive ones. Binding mode analysis of pterin derivatives in the RTA active site and their correlations with SMD profiles. According to Ruiz-Carmona and co-workers60⁠ , the hydrogen bonds play a variable but substantial contribution to the binding free energy (∆Gbind). These authors pointed out also that one or more hydrogen bonds represent the minimal binding unit in protein-ligand complexes.60 In addition, breaking these strong non-covalent interactions through undocking procedures can provide important information about the binding ability of the compounds under study.60 In our work, we performed a visual inspection of the most relevant hydrogen bonds (we adopted the criterion of hydrogen bonds with distances shorter than ~ 2.0 Å) between the pterin



derivatives and RTA active site, using the starting bound pose. Such poses are depicted in Figure 7.

INSERT FIGURE 7 We were able to make a correlation between force and work profiles with RTA structure and the hydrogen bonds that are made with the pterin derivatives. From that analysis, we noted that the number of hydrogen bonds keep a strong relation with the force and work peaks obtained in the SMD simulations. The more potent the inhibitor is, the more hydrogen bonds it makes in the RTA active site, especially with G121 and V81. Hydrogen bonds with the backbone oxygen atom of G121 and V81 appear to be very important interactions for the pterin derivatives in the RTA active site. Because these HBs, the pterin moiety makes a strong π-π stacking interaction between Y80 and Y121, similar to adenine (natural substrate) on RTA active site. For the 1MX and JP3 inhibitors we have a different situation. These compounds present a protonation on the nitrogen atom of the pterin group (Figure S2), leading to a clash with the backbone nitrogen atom of V81 (Figures 7D and 7F). This clash causes a displacement of these inhibitors out of the cavity and a considerable decrease in the force and work profiles. 1MX and JP3 inhibitors had a mean decrease of 〈 〉 ≈ 237.31 and 469.69 pN and 〈〉 ≈ 27.68 and 42.29 kcal/mol, respectively, when compared to 19M, RS8 and 0RB. These observations suggest that such protonation (Figure S2) is an inefficient strategy to increase the anti-RTA activity of this class of compounds. In addition, some studies have been shown a significant correlation between 〈〉 and residence time.38,61 In this regard, we suggest


Page 18 of 43

Page 19 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


that these modifications cause a decrease in the residence time of this class of compounds in the RTA active site. The inhibitors that show such modification (1MX and JP3) presented fast displacement and low non-equilibrium work 〈〉. Molecular docking combined with fast SMD simulation. The combination of molecular docking with undocking simulations could be an effective strategy of post-docking filtering. In our study, we showed that the non-equilibrium work obtained with fast undocking simulations can be used as a score function to distinguish active compounds from inactive ones. However, for others protein-ligand complexes, the success of such strategy strongly depends of calibration and validation steps. As mentioned earlier, hydrogen bonds provide the minimal binding unit in protein-ligand complexes and contribute for the structural stability.60⁠ Thus, it is expected that an incorrect ligand conformation results in wrong prediction of relative binding affinity by undocking procedure. For that, redocking is commonly used to evaluate the effectiveness of docking algorithms in predicting ligand bioactive conformation in the binding pocket. Moreover, minimization and equilibration steps carried out by MD simulation are also commonly used for structure optimization (e.g.: bond geometries optimization and removing unfavorable contacts).61⁠ For undocking procedures, the main calibration steps include pulling direction, spring constant (k) and pulling velocity (v). From these, the pulling direction is the most important parameter to be tuned, because it directly alters the force profile, which strongly affects the work profile. According to Vuong and co-workers43⁠ , the optimal pulling direction should have a minimum steric hindrance along the movement of the ligand. A wrong pulling direction can



provide a pulling pathway with more hindrance and consequently a wrong description of the force profile. Thus, the choice of the best pulling direction is a critical point in SMD process, which affect the magnitude of the non-equilibrium work profile and, consequently, the rank of compounds. In this sense, there are some quantitative methods to determine the optimal direction for pulling ligand from binding pocket, e.g.: Random Accelerated Molecular Dynamics (RAMD),62 CAVER⁠ ,63 MOLE⁠ ,64 MolAxis and Minimal Steric Hindrance (MSH) method.43 Nevertheless, in our work, we used a simple analysis to determine the optimal direction for decoupling ligand; even so, by adopting such strategy, we got good correlations between SMD force profiles and half-maximum inhibitory concentration. A further point is the impact of pulling velocity and the CPU time consumption. For drug discovery purposes, discrimination between active and inactive inhibitors can be achieved using large pulling velocities, which in turn demands low computational cost and suitable to drug discovery process.44⁠ However, as we previously discussed, very high pulling velocities presented poor correlation between 〈〉 and IC50 (Supporting Information Table S3), thus, calibration steps are highly recommended. In addition, efforts to reduce the computational cost also can be achieved, by varying the number of SMD trajectories and, then, analyzing their impact in R2 correlation coefficient (Supporting Information Table S3). For example, when we used the second parameter set (k = 10 kcal/mol/Å²; v = 0.05 Å/ps) and only five SMD trajectories, we got a correlation of R² = 0.975 with 9 ns of non-equilibrium simulation (Supporting Information Table S3).


Page 20 of 43

Page 21 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CONCLUSION In this study, we used molecular docking and fast SMD simulations as a strategy to measure the relative binding affinities of six RTA inhibitors. We used the effectiveness of Autodock Vina to predict ligand bioactive poses into RTA active site and non-equilibrium work 〈〉 as scoring function in order to distinguish RTA active inhibitors from inactive ones. For 〈〉, we obtained distinct relative binding affinity profiles for the six RTA inhibitors and an excellent correlation with experimental data. In addition, using 〈〉 we were able to correctly rank the RTA inhibitors according to IC50. The rupture forces 〈〉 and 〈 〉 also showed reasonable correlation with experimental data; however, these profiles were not able to correctly rank all ligands. In addition, decreases in pulling velocity improved the correlation between 〈〉 and the experimental results. The choice of v must consider the balance between computational time and accuracy of results. Our results suggest that v = 0.05 Å/ps is suitable for estimating 〈〉 for ranking binding affinities, considering as spring constant k = 2 or 10 kcal/mol/Å2. Furthermore, SMD force profiles provided insights about the structure-activity relationship of the pterin derivatives. Finally, we believe that the computational strategy presented in our study as well as our results can contribute for the in silico assessment of new pterin derivatives for RTA.



Supporting Information

All theoretical results and additional information about SMDAtoms of the ligands are available in the Supporting Information.


Page 22 of 43

Page 23 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


FIGURES

Figure 1. Superposition of predicted (green) and native (purple) poses. PDB Codes: 19M (A); RS8 (B); 0RB (C); 1MX (D); JP2 (E); JP3 (F).



Figure 2. The pulling directions that we used for decoupling the six RTA inhibitors from their active sites. The red circle corresponds to Cα of I172, the FixedAtom. The blue arrow means the vector that connect FixedAtom to SMDAtom of the ligand (yellow). In green are depicted the residues around the ligand that represent possible steric hindrances to the reaction coordinate. Each plot corresponds to the following inhibitors: 19M-RTA (A); RS8-RTA (B); 0RB-RTA (C); 1MX-RTA (D); JP2-RTA (E); JP3-RTA (F).


Page 24 of 43

Page 25 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Mean force profiles 〈〉 of the six RTA inhibitors. (A) 〈〉 versus time. The maximum force corresponds to the mean force necessary to break protein-ligand interactions and decoupling inhibitors from their active site; (B) 〈〉 versus ligand displacement. Potent inhibitors presented small displacements in initial steps of SMD simulations. In contrast, weak inhibitors presented large displacements.



Figure 4. Correlations between IC50 and theoretical results for the six RTA inhibitors. Top-left (IC50 x 〈 〉), top-right (IC50 x 〈〉 ), bottom-left (IC50 x 〈〉) and bottom-right (IC50 x Edocking).


Page 26 of 43

Page 27 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Mean non-equilibrium work profiles 〈〉 versus time. The area under the curve corresponds to mean work needed to decouple the inhibitors from their active sites. These profiles showed different energy levels for each undocking procedure. We used them to rank the six RTA inhibitors according to IC50.



Figure 6. Correlations between mean non-equilibrium work 〈〉 and different pulling velocities (v). The plots correspond to the following ligands:19M (A); RS8 (B); 0RB (C); 1MX (D); JP2 (E); JP3 (F).


Page 28 of 43

Page 29 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7. Binding mode for all RTA-inhibitors (yellow). Hydrogen bonds are depicted by dotted black line. The HB bond lengths are also presented. Some atoms were omitted to make better the ligand pose visualization. Each plot corresponds to the following ligands: 19M (A); RS8 (B); 0RB (C); 1MX (D); JP2 (E); JP3 (F).



SCHEME

Scheme 1. Schematic diagram of the computational strategy used in our study to determine the relative binding affinities of RTA-inhibitors.


Page 30 of 43

Page 31 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TABLES Table 1. Chemical structures and IC50 values of the RTA inhibitors investigated in this study. Pterin group are highlighted in blue. Chemical structures were built using MarvinSketch from Chemaxon package. PDB ID Chemical Structure Ligand Complex

IC50 (µM)

Reference

19M

4HUP

15

24

RS8

4HUO

20

24

0RB

4ESI

70

22

1MX

4MX1

209

23

JP2

3PX8

230

25

JP3

3PX9

380

25



Page 32 of 43

Table 2. Calculated data of 〈〉 , 〈 〉 and 〈〉 and Edockinga for the ricin inhibitors used in our study. The proposed rankings are organized in descending order of RTA-inhibition potency. See table 1 for identification of the ligands.

PDB

IC50 (rank)

〈&〉'() (rank)

〈*〉'()

〈&'() 〉 (rank)

〈+〉 (rank)

Edocking (rank)

4HUO

15 (1)

852.96 (2)

157.56

979.70 ± 82.39 (2)

95.01 ± 9.10 (1)

-9.34 (3)

4HUP

20 (2)

884.57 (1)

170.82

1056.19 ± 129.92 (1)

91.98 ± 10.15 (2)

-9.40 (2)

4ESI

70 (3)

724.16 (3)

130.10

878.59 ± 89.13 (3)

73.51 ± 9.77 (3)

-9.20 (4)

4MX1

209 (4)

590.55 (4)

129.18

741.15 ± 93.04 (4)

57.34 ± 7.68 (4)

-10.09 (1)

3PX8

230 (5)

499.71 (5)

115.14

649.12 ± 57.32 (5)

48.52 ± 5.92 (5)

-7.90 (5)

3PX9

380 (6)

309.27 (6)

81.08

449.82 ± 61.08 (6)

31.50 ± 7.10 (6)

-7.85 (6)

a

〈〉 and 〈 〉 are measured in pN, while 〈〉 and Edocking are in kcal/mol. 〈$〉 is in ps.


Page 33 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


AUTHOR INFORMATION Corresponding Author ‡Gerd B. Rocha, ‡ Department of Chemistry, Federal University of Paraíba, CEP: 58051-900. João Pessoa, PB Brazil. E-mail: [email protected]

† Demétrius A. M. Araújo, † Department of Biotechnology, Federal University of Paraíba, CEP: 58051-900.João Pessoa, PB, Brazil. E-mail: [email protected]

Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.



Funding Sources The authors are grateful for financial support from the research project “Bioinformática Estrutural de Proteínas: Modelos, Algoritmos e Aplicações Biotecnológicas (Edital Biologia Computacional 51/2013, AUXPE1375/2014/CAPES)”. G.B.R. acknowledges support from the Brazilian National Council for Scientific and Technological Development (CNPq grant no. 305271/2013-0).

ACKNOWLEDGMENT The authors would also like to thank the Laboratory of Molecular Modeling and Chemical Reactions and Versatus HPC (High Performance Computing), for support in the computational resources used in this study.

ABBREVIATIONS RTA, Ricin Toxin A; RTB, Ricin Toxin B; SMD, Steered Molecular Dynamics; MD, Molecular Dynamics; SRL, Sarcin-Ricin Loop; ER, Endoplasmic Reticulum; RAMD, Random Molecular Dynamics; JE, Jarzinsky Equality; MSH, Minimal Steric Hindrance.


Page 34 of 43

Page 35 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


REFERENCES (1) Robertus, J. D.; Michael Piatak; Ferris, R.; Houston, L. L. Crystallization of Ricin A Chain Obtained from a Cloned Gene Expressed in E. Coli. J. Biol. Chem. 1986, 262 (1), 19–20. (2)

Bozza, W. P.; Tolleson, W. H.; Rivera Rosado, L. A.; Zhang, B. Ricin Detection: Tracking

Active Toxin. Biotechnol. Adv. 2015, 33 (1), 117–123. (3)

Montfort, W.; Villafrancas, J. E.; Monzingol, A. F.; Ernstll, S. R.; Katzinq, B.; Earl

Rutenberll, N.; XuonglJ, H.; HamlinII, R.; Robertusq, J. D. The Three-Dimensional Structure of Ricin at 2.8 A. J. Biol. Chem. 1987, 262 (11), 5396–5403. (4)

Rutenber, E.; Katzin, B. J.; Ernst, S.; Collins, E. J.; Mlsna, D.; Ready, M. P.; Robertus, J.

D. Crystallographic Refinement of Ricin to 2.5 A. Proteins 1991, 10 (3), 240–250. (5)

Robertus, J. D.; Monzingo, A. F. The Structure of Ribosome Inactivating Proteins. Mini

Rev. Med. Chem. 2004, 4 (5), 477–486. (6)

Endo, Y.; Tsurugi, K. RNA N-Glycosidase Activity of Ricin A-Chain. Mechanism of

Action of the Toxic Lectin Ricin on Eukaryotic Ribosomes. J. Biol. Chem. 1987, 262 (17), 8128– 8130. (7)

Monzingo, A. F.; Robertus, J. D. X-Ray Analysis of Substrate Analogs in the Ricin A-

Chain Active Site. J. Mol. Biol. 1992, 227 (4), 1136–1145. (8)

May, K. L.; Yan, Q.; Tumer, N. E. Targeting Ricin to the Ribosome. Toxicon 2013, 69,

143–151. (9)

Audi, J.; Belson, M.; Patel, M.; Schier, J.; Osterloh, J. Ricin Poisoning. Am. Med. Assoc.

2005, 294 (18), 2342–2351.



(10) Spooner, R. A.; Watson, P. D.; Marsden, C. J.; Smith, D. C.; Moore, K. A. H.; Cook, J. P.; Lord, J. M.; Roberts, L. M. Protein Disulphide-Isomerase Reduces Ricin to Its A and B Chains in the Endoplasmic Reticulum. Biochem. J 2004, 383, 285–293. (11) Argent, R. H.; Parrott, A. M.; Day, P. J.; Roberts, L. M.; Stockley, P. G.; Lord, J. M.; Radford, S. E. Ribosome-Mediated Folding of Partially Unfolded Ricin A-Chain. J. Biol. Chem. 2000, 275 (13), 9263–9269. (12) Simpson, J. C.; Roberts, L. M.; Römisch, K.; Davey, J.; Wolf, D. H.; Lord, J. M. Ricin A Chain Utilises the Endoplasmic Reticulum-Associated Protein Degradation Pathway to Enter the Cytosol of Yeast. FEBS Lett. 1999, 459 (1), 80–84. (13) Wesche, J.; Rapak, A.; Olsnes, S. Dependence of Ricin Toxicity on Translocation of the Toxin A-Chain from the Endoplasmic Reticulum to the Cytosol. J. Biol. Chem. 1999, 274 (48), 34443–34449. (14) Hesselberth, J. R.; Miller, D.; Robertus, J.; Ellington, A. D. In Vitro Selection of RNA Molecules That Inhibit the Activity of Ricin A-Chain. J. Biol. Chem. 2000, 275 (18), 4937–4942. (15) Wahome, P. G.; Bai, Y.; Neal, L. M.; Robertus, J. D.; Mantis, N. J. Identification of Small-Molecule Inhibitors of Ricin and Shiga Toxin Using a Cell-Based High-Throughput Screen. Toxicon 2010, 56 (3), 313–323. (16) Spooner, R.; Lord, J. Ricin Trafficking in Cells. Toxins (Basel). 2015, 7 (1), 49–65. (17) Jasheway, K.; Pruet, J.; Anslyn, E. V.; Robertus, J. D. Structure-Based Design of Ricin Inhibitors. Toxins (Basel). 2011, 3 (10), 1233–1248.


Page 36 of 43

Page 37 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(18) Yan, X.; Hollis, T.; Svinth, M.; Day, P.; Monzingo, a F.; Milne, G. W.; Robertus, J. D. Structure-Based Identification of a Ricin Inhibitor. J. Mol. Biol. 1997, 266 (5), 1043–1049. (19) Miller, D. J.; Ravikumar, K.; Shen, H.; Suh, J.-K.; Kerwin, S. M.; Robertus, J. D. Structure-Based Design and Characterization of Novel Platforms for Ricin and Shiga Toxin Inhibition. J. Med. Chem. 2002, 45 (1), 90–98. (20) Bai, Y.; Monzingo, A. F.; Robertus, J. D. The X-Ray Structure of Ricin A Chain with a Novel Inhibitor. Arch. Biochem. Biophys. 2009, 483 (1), 23–28. (21) Bai, Y.; Watt, B.; Wahome, P. G.; Mantis, N. J.; Robertus, J. D. Identification of New Classes of Ricin Toxin Inhibitors by Virtual Screening. Toxicon 2010, 56 (4), 526–534. (22) Pruet, J. M.; Saito, R.; Manzano, L. a; Jasheway, K. R.; Wiget, P. a; Kamat, I.; Anslyn, E. V; Robertus, J. D. Optimized 5-Membered Heterocycle-Linked Pterins for the Inhibition of Ricin Toxin A. ACS Med. Chem. Lett. 2012, 3 (7), 588–591. (23) Wiget, P. A.; Manzano, L. A.; Pruet, J. M.; Gao, G.; Saito, R.; Monzingo, A. F.; Jasheway, K. R.; Robertus, J. D.; Anslyn, E. V. Sulfur Incorporation Generally Improves Ricin Inhibition in Pterin-Appended Glycine-Phenylalanine Dipeptide Mimics. Bioorg. Med. Chem. Lett. 2013, 23 (24), 6799–6804. (24) Saito, R.; Pruet, J. M.; Manzano, L. A.; Jasheway, K.; Monzingo, A. F.; Wiget, P. A.; Kamat, I.; Anslyn, E. V; Robertus, J. D. Peptide-Conjugated Pterins as Inhibitors of Ricin Toxin A. J. Med. Chem. 2013, 56 (1), 320–329.



(25) Pruet, J. M.; Jasheway, K. R.; Manzano, L. A.; Bai, Y.; Anslyn, E. V; Robertus, J. D. 7Substituted Pterins Provide a New Direction for Ricin A Chain Inhibitors. Eur. J. Med. Chem. 2011, 46 (9), 3608–3615. (26) Ramírez, D. Computational Methods Applied to Rational Drug Design. Open Med. Chem. J. 2016, 10, 7–20. (27) Kumar, A.; Zhang, K. Y. J. Hierarchical Virtual Screening Approaches in Small Molecule Drug Discovery. Methods 2015, 71 (C), 26–37. (28) Aqvist, J. A New Method for Predicting Binding Affinity in Computer-Aided Drug Design. Protein Eng. 1994, 7 (3), 385–391. (29) Genheden, S.; Ryde, U. The MM/PBSA and MM/GBSA Methods to Estimate LigandBinding Affinities. Expert Opin. Drug Discov. 2015, 10 (5), 449–461. (30) Chipot, C. Frontiers in Free-Energy Calculations of Biological Systems. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4 (1), 71–89. (31) Isralewitz, B.; Gao, M.; Schulten, K. Steered Molecular Dynamics and Mechanical Functions of Proteins. Curr. Opin. Struct. Biol. 2001, 11 (2), 224–230. (32) Isralewitz, B.; Baudry, J.; Gullingsrud, J.; Kosztin, D.; Schulten, K. Steered Molecular Dynamics Investigations of Protein Function. J. Mol. Graph. Model. 2001, 19 (1), 13–25. (33) Suan Li, M.; Khanh Mai, B. Steered Molecular Dynamics-A Promising Tool for Drug Design. Curr. Bioinform. 2012, 7 (4), 342–351.


Page 38 of 43

Page 39 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(34) Lv, C.; Tan, C.; Qin, M.; Zou, D.; Cao, Y.; Wang, W. Low Folding Cooperativity of Hp35 Revealed by Single-Molecule Force Spectroscopy and Molecular Dynamics Simulation. Biophys. J. 2012, 102 (8), 1944–1951. (35) Zhang, Y.; Lou, J. The Ca2+ Influence on Calmodulin Unfolding Pathway: A Steered Molecular Dynamics Simulation Study. PLoS One 2012, 7 (11). (36) De Fabritiis, G.; Coveney, P. V.; Villà-Freixa, J. Energetics of K+ Permeability through Gramicidin A by Forward-Reverse Steered Molecular Dynamics. Proteins Struct. Funct. Genet. 2008, 73 (1), 185–194. (37) Colizzi, F.; Perozzo, R.; Scapozza, L.; Recanatini, M.; Cavalli, A. Single-Molecule Pulling Simulations Can Discern Active from Inactive Enzyme Inhibitors. J. Am. Chem. Soc. 2010, 132 (21), 7361–7371. (38) Jorgensen, W. L. Drug Discovery: Pulled from a Protein’s Embrace. Nature 2010, 466 (7302), 42–43. (39) Costantino, G.; Capelli, A. M.; Bruno, A.; Guadix, A. J. E. Unbinding Pathways from Glucocorticoid Receptor Shed Light on Reduced Sensitivity of Glucocorticoid Ligands to a Naturally Occurring, Clinically Relevant, Mutant Receptor. J. Med. Chem. 2013, 56 (17), 7003– 7014. (40) Ngo, S. T.; Hung, H. M.; Nguyen, M. T. Fast and Accurate Determination of the Relative Binding Affinities of Small Compounds to HIV-1 Protease Using Non-Equilibrium Work. J. Comput. Chem. 2016, 1–9.



(41) Mai, B. K.; Viet, M. H.; Li, M. S. Top Leads for Swine Influenza A/H1N1 Virus Revealed by Steered Molecular Dynamics Approach. J. Chem. Inf. Model. 2010, 50 (12), 2236–2247. (42) Mai, B. K.; Li, M. S. Neuraminidase Inhibitor R-125489 - A Promising Drug for Treating Influenza Virus: Steered Molecular Dynamics Approach. Biochem. Biophys. Res. Commun. 2011, 410 (3), 688–691. (43) Vuong, Q. Van; Nguyen, T. T.; Li, M. S. A New Method for Navigating Optimal Direction for Pulling Ligand from Binding Pocket: Application to Ranking Binding Affinity by Steered Molecular Dynamics. J. Chem. Inf. Model. 2015, 55 (12), 2731–2738. (44) Patel, J. S.; Berteotti, A.; Ronsisvalle, S.; Rocchia, W.; Cavalli, A. Steered Molecular Dynamics Simulations for Studying Protein-Ligand Interaction in Cyclin-Dependent Kinase 5. J. Chem. Inf. Model. 2014, 54 (2), 470–480. (45) Singh, V. K.; Chang, H. H.; Kuo, C. C.; Shiao, H. Y.; Hsieh, H. P.; Coumar, M. S. Drug Repurposing for Chronic Myeloid Leukemia: In Silico and in Vitro Investigation of DrugBank Database for Allosteric Bcr-Abl Inhibitors. J. Biomol. Struct. Dyn. 2016, 35 (8), 1833–1848. (46) Wang, R.; Fang, X.; Lu, Y.; Wang, S. The PDBbind Database: Collection of Binding Affinities for Protein-Ligand Complexes with Known Three-Dimensional Structures. J. Med. Chem. 2004, 47 (12), 2977–2980. (47) Wang, R.; Fang, X.; Lu, Y.; Yang, C.-Y.; Wang, S. The PDBbind Database: Methodologies and Updates. J. Med. Chem. 2005, 48 (12), 4111–4119.


Page 40 of 43

Page 41 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(48) Dolinsky, T. J.; Czodrowski, P.; Li, H.; Nielsen, J. E.; Jensen, J. H.; Klebe, G.; Baker, N. A. PDB2PQR: Expanding and Upgrading Automated Preparation of Biomolecular Structures for Molecular Simulations. Nucleic Acids Res. 2007, 35 (SUPPL.2), 522–525. (49) Morris, G. M. Huey, R. Lindstrom, W.; Sanner, M. F. Belew, R. K. Goodsell, D. S. Olson, A. J. AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility. J. Comput. Chem. 2009, 30 (16), 2785–2791. (50) O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An Open Chemical Toolbox. J. Cheminform. 2011, 3 (10), 1–14. (51) Steffen, C.; Thomas, K.; Huniar, U.; Hellweg, A.; Rubner, O.; Schroer, A. AutoDock Vina: Improving the Speed and Accuracy of Docking with a New Scoring Function, Efficient Optimization, and Multithreading. J. Comput. Chem. 2010, 31, 2967–2970. (52) Humphrey, W.; Dalke, A.; Schulten, K. {VMD} -- {V}isual {M}olecular {D}ynamics. J. Mol. Graph. 1996, 14, 33–38. (53) Neudert, G.; Klebe, G. Fconv: Format Conversion, Manipulation and Feature Computation of Molecular Data. Bioinformatics 2011, 27 (7), 1021–1022. (54) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins Struct. Funct. Bioinforma. 2006, 65 (3), 712–725. (55) Wang, J.; Wolf, R. M. .; Caldwell, J. W. .; Kollman, P. A. .; Case, D. A. Development and Testing of a General AMBER Force Field. J. Comput. Chem. 2004, 25, 1157–1174.



(56) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781–1802. (57) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free Energy Calculation from Steered Molecular Dynamics Simulations Using Jarzynski’s Equality. J. Chem. Phys. 2003, 119 (6), 3559. (58) Li, D.-C.; Ji, B.-H. Free Energy Calculation of Single Molecular Interaction Using Jarzynski’s Identity Method: The Case of HIV-1 Protease Inhibitor System. Acta Mech. Sin. 2012, 28 (3), 891–903. (59) Houston, D. R.; Walkinshaw, M. D. Consensus Docking: Improving the Reliability of Docking in a Virtual Screening Context. J. Chem. Inf. Model. 2013, 53 (2), 384–390. (60) Ruiz-Carmona, S.; Schmidtke, P.; Luque, F. J.; Baker, L.; Matassova, N.; Davis, B.; Roughley, S.; Murray, J.; Hubbard, R.; Barril, X. Dynamic Undocking and the Quasi-Bound State as Tools for Drug Discovery. Nat. Chem. 2017, 9 (3), 201–206. (61) Sliwoski, G.; Kothiwale, S.; Meiler, J.; Lowe, E. W. E. Computational Methods in Drug Discovery. Pharmacol. Rev. 2014, 66 (1), 334–395. (62) Niu, Y.; Li, S.; Pan, D.; Liu, H.; Yao, X. Computational Study on the Unbinding Pathways of B-RAF Inhibitors and Its Implication for the Difference of Residence Time: Insight from Random Acceleration and Steered Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2016, 18 (7), 5622–5629.


Page 42 of 43

Page 43 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(63) Chovancova, E.; Pavelka, A.; Benes, P.; Strnad, O.; Brezovsky, J.; Kozlikova, B.; Gora, A.; Sustr, V.; Klvana, M.; Medek, P.; Biedermannova, L.; Sochor, J.; Damborsky, J. CAVER 3.0: A Tool for the Analysis of Transport Pathways in Dynamic Protein Structures. PLoS Comput. Biol. 2012, 8 (10), 23–30. (64) Sehnal, D.MOLE 2.0 : Advanced Approach for Analysis of Biomacromolecular Channels. J. Cheminform. 2013, 5 (39), 1–13.


Determining the Relative Binding Affinity of Ricin Toxin A (RTA

Recommend Documents