- Email: [email protected]
Theoretical studies on the mechanism of sugammadex for the reversal of aminosteroid-induced neuromuscular blockade
Journal of Molecular Liquids 265 (2018) 450–456
Contents lists available at ScienceDirect
Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq
Theoretical studies on the mechanism of sugammadex for the reversal of aminosteroid-induced neuromuscular blockade Linwei Li a, Yanan Zhou a, Zhengjun Wang a, Chengjun Wu a, Zhen Li a, Changshan Sun b,⁎, Tiemin Sun a,⁎ a b
Key Laboratory of Structure-Based Drug Design and Discovery, Ministry of Education, Shenyang Pharmaceutical University, Shenyang 110016, PR China Pharmacy Department, Shenyang Pharmaceutical University, Shenyang 110016, PR China
a r t i c l e
i n f o
Article history: Received 3 May 2018 Received in revised form 1 June 2018 Accepted 8 June 2018 Available online 15 June 2018 Keywords: Sugammadex Neuromuscular blockers Mechanism DFT MD
a b s t r a c t The mechanism of action of sugammadex in reversing neuromuscular blockers (NMBs) was theoretically studied by extensive molecular modeling (molecular docking, molecular dynamic and molecular mechanics-quantum mechanics). The NMBs-sugammadex inclusion complex was established by an enthalpy-driven process, which stabilized by hydrophobic interaction, hydrogen bonding and electrostatic interaction between NMBs and sugammadex. Furthermore, the binding energy of NMBs-sugammadex indicated that vecuronium (ΔE = −25.29 kJ mol−1) had lower binding affinity to sugammadex than rocuronium (ΔE = −34.61 kJ mol−1), thus the former one would be more difficult to be reversed which was consist with the phenotype of clinical trials. Molecular electrostatic potential, natural bond orbital and frontier molecular orbitals were also utilized to give explanation for the discrepancy of interaction strength between rocuronium-sugammadex and vecuroniumsugammadex inclusion complexes. In summary, the molecular modeling study made it possible to explain the different reversal efficiency of NMBs by sugammadex, which could be used to predict the drug-receptor interactions during the design phase of novel compounds development and also provide great opportunity to design novel γ-cyclodextrin-based muscle relaxant antagonists. © 2018 Published by Elsevier B.V.
1. Introduction Neuromuscular blockers (NMBs) are routinely used for facilitating surgical procedures and tracheal intubation during anesthesia. They usually bind to nicotinic acetylcholine receptor (nAChR) at the neuromuscular junction, and then block neurotransmission induced by acetylcholine (ACh) [1–7]. One such typical and widely used species is the non-depolarizing aminosteroid-based NMBs, such as pancuronium (proved in 1973) [8], vecuronium (1982) [9], and rocuronium (1994) [10], which behaviors fast onset timers and short duration of action. Although aminosteroid-based NMBs have few adverse effects during anesthesia, the residual duration of muscle relaxants beyond the end of the surgery, also referred to as postoperative residual neuromuscular blockers (RNMBs), is a well-known problem [11–16]. The postoperative residual neuromuscular blockers associate with respiratory insufficiency, pulmonary complications, upper airway obstruction and hypoxia, which can increase the incidence of tracheal re-intubation and mortality [17–23]. Therefore, it is necessary to exploit rapid and complete agents for controlling the duration of the residual neuromuscular blockers at the end of operation, which be used to accelerate the recovery of patients' muscle function. ⁎ Corresponding authors. E-mail addresses: [email protected] (C. Sun), [email protected] (T. Sun).
https://doi.org/10.1016/j.molliq.2018.06.033 0167-7322/© 2018 Published by Elsevier B.V.
Chemical encapsulation (Fig. 1) is a recent mechanism to remove the RNMBs drugs; and there have been investigations into exogenous macrocyclic host molecules, such as cyclodextrins (CDs) and anionic cyclophanes, as innovative and potential reversal agents [19, 24–27]. Inclusion complexes can be formed by aminosteroidal NMBs (guest molecule) and chemically modified γ-cyclodextrin (host molecule) with high binding affinity, resulting in a decrease of NMBs in plasma, then rapidly reverse the effect of residual neuromuscular blockers. Furthermore, γ-cyclodextrin derivatives (γ-CDs) do not bind to the muscarinic receptors, which lead its reversal effect more reliable than acetylcholinesterase inhibitors. Therefore, a series of modified γ-CDs have been synthesized with the aim of constructing host molecules capable of forming host-guest complexes with neuromuscular blockers [17, 25, 28]. Among them, sugammadex (Bridion®, octakis-(6-deoxy-6-Smercaptopropionyl)-γ-CD sodium salt) is demonstrated to be the first potent and selective relaxant binding agent, exerts its effect by forming very tight water-soluble inclusion complex at a 1: 1 ratio with aminosteroidal NMBs, approved by FDA at 2015 [28–32]. (See Figs. 2–4.) It was reported that the binding constant for vecuroniumsugammadex complex is somewhat smaller (Ki = 5.72 × 106 M−1) than rocuronium (Ki = 1.79 × 107 M−1) inclusion complex, while still weaker binding to pancuronium (Ki = 2.62 × 106 M−1) is detected [32]. This, in turn, lead to the suggestion that rocuronium is removed
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
451
Fig. 1. The mechanism of action of γ-cyclodextrin derivatives on the nAChR.
more effectively by sugammadex than the other aminosteroids, such as vecuronium and pancuronium. Further, in order to achieve equivalent reversal potency for low binding affinity NMBs, excessive sugammadex would be employed in clinical, and therefore, it could increase the possibility of bleeding risk induced by sugammadex [33]. Thus, it is better to exploit new potent γ-CD-based muscle relaxant antagonists, especially, screen a specific γ-CD derivative for a given aminosteroidal NMB. Different interaction modes and strengths between host and various guest molecules are attributing to distinct reversal efficiency. Molecular recognition mechanism is the prerequisite foundation to learn about the binding mode of inclusion complex, which would obviously help in designing new synthetic γ-CDs. Considering these concepts, we attempted to give explanations of the exact nature of reversal divergence in reversal process based on reversal mechanism at molecular level. Up to data, among NMBs-sugammadex complexes, only the geometric structure of rocuronium-sugammadex inclusion complex has been established by X-ray diffraction (XRD) [34] and NMR technique [35]. Additionally, the growing condition of single crystals and the limitation of XRD make XRD method infeasible [36, 37]; while NMR can't supply the exact distances of bonds in encapsulation complex. Therefore, the development of novel γ-CD-based muscle relaxant antagonists is still retard; and the evaluation of reversal potency is mainly based on pharmacological test. Luckily, molecular modeling methods have been recently proposed as powerful tools to unravel the interaction between
host and guest molecules and assist in elucidating the reversal process, and then evaluate the reversal efficacy. The methods involve molecular mechanics (MM), molecular docking [38–41], molecular dynamics (MD) [42–46], quantum mechanics (QM) [47–50] and molecular mechanics-quantum mechanics (QM-MM) calculations [51–57]. Molecular docking could rapidly provide the binding model of guest molecule with host molecule at low molecular mechanics level [58, 59]. Molecular dynamics simulation is employed which provides the equilibrium and transport properties of a classical many-body system in liquid circumstance. The major advantage of molecular dynamics is the simulated conformations that are less dependent on the initial structure than those obtained by molecular mechanic [58–61]. Thus, the MD results are very similar to real experiments. A multi-layered hybrid approach, our own N-layered integrated molecular orbital and molecular mechanics method (ONIOM), is capturing for macromolecular system since it can simultaneously treat different parts of a system with multiple approaches of varying accuracy and computational cost (QM and MM). In summary, the combination of the mentioned theoretical methods (molecular docking, MD and QM-MM) offers successful and rational data in solving structural, energetic and dynamic problems. Herein, with the aid of computational simulation, we aim to study the reversal mechanism of sugammadex with NMBs and explain the reversal efficacy among different NMBs at molecular level. Rocuronium (Rocu) and vecuronium (Vecu) were selected as the examples to
Fig. 2. The chemical structure and schematic diagram of sugammadex.
452
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
Fig. 3. The chemical structures of rocuronium and vecuronium.
perform our investigation. Molecular docking, MD, QM-MM methods were simultaneously utilized to establish the energetically most favorable encapsulation complexes of NMBs with sugammadex. Subsequently, QNIOM, natural bond orbital (NBO), molecular electrostatic potential (MEP) and frontier molecular orbitals (FMOs) were conducted to evaluate binding efficacy of inclusion complexes from electrostatic interaction, hydrogen bonding and chemical stability aspects. Finally, thermodynamic parameters were calculated to describe encapsulation process of NMBs-sugammadex complexes along with the involved driving force. This work presents a deep insight into the antagonistic mechanism by sugammadex, which can be adopted to design novel γ-CDbased muscle relaxant antagonists for a given steroidal neuromuscular blocker, as well as can be used to rapidly predict potency before pharmacological test.
2. Computational methods The initial structures of rocuronium and vecuronium depicted by ChemBioDraw Ultra 13.0 were imported to the Spartan 14.0 for conformational searching by using the MM+ force field method. Then selected a series of low-energy conformers that were in equilibrium with each other, and their relative energies were within a region of 4 kcal mol−1. To confirm the more stable structures of guest molecules, the above conformers were further calculated at the density functional theory (DFT) level with a Becke-3-Lee-Yang-Parr (B3LYP) functional and 6–31 + G(d, p) basis set by using Gaussian 09 [62] suite of programs. Then redistributed energy values according to the Boltzmann distribution, permitting decline to the most stable structures of isolated rocuronium and vecuronium for further simulations. The γcyclodextrin structure was obtained from the crystallographic parameters provided by the Cambridge crystallographic Data center, then based on the initial structure, sugammadex was constructed and then fully
optimized at the parameterized model 3 (PM3) level of theory without imposing any symmetrical. Molecular docking studies were undertaken using AutoDock-Vina program [63]. rocuronium and vecuronium were allowed to dock anywhere in a 35 Å × 35 Å × 35 Å grid box with 1 Å spacing centered the center of mass of sugammadex, with the Lamarckian genetic algorithm (LGA). The initial torsions and positions of rocuronium and vecuronium were generated randomly. The maximum numbers of energy evaluation and generation were set to 250,000 and 27,000, respectively. A total of 20 separate docking runs were performed with the initial population of 20 individuals, and the lowest energy structure from each run was retained. Final docked conformations were clustered by using a tolerance of 1.0 Å root mean square deviation (RMSD). AMBERTOOL 12 [64] software package was used for the molecular dynamic (MD) simulations using the general amber force field (GAFF) parameter set [65]. Topology and parameter files were generated with the LEaP program on structure of the complexes obtained by the aforementioned docking procedures. The Sander module of AMBERTOOL 12 package was used to perform MD simulations. Each system was solvated in a truncated octahedral box of TIP3P water molecules with a minimum solute-water distance of 10 Å. All bonds involving hydrogen atoms were constrained using the SHAKE flag [66]. Long-range electrostatic interactions were calculated using the particle-mesh Ewald (PME) algorithm [67], with a cutoff of 12 Å for Lennard-Jones interactions. Periodic boundary conditions were applied to avoid edge effects. Prior to MD production, the steepest-descent method was used for the 500 iterations and the conjugated gradient method for the subsequent 2500 iteration in each step for the solvent and the entire model system, respectively. Afterward, the entire system was gradually heated up to 300 K gradually over 50 ps using the constant volume and normal temperature (NVT) [68] ensemble with the solutes restrained by a weak harmonic potential. During the heating, time constant for heat bath coupling for the solute was set as 1.0. Subsequently, the constant normal
Fig. 4. Docked geometries of energetically most favorable inclusion complexes of rocuronium (A) and vecuronium (B) with sugammadex. The polar parts of NMBs protrude outside the cavity from the opening of the secondary OH side, and all four steroidal rings are embodied within the hydrophobic cavity of sugammadex extended by side chains.
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
pressure and normal temperature (NPT) ensemble was carried out for 500 ps, keeping a temperature of 300 K by Berendsen thermostat, and maintaining pressure constant at 1.0 atm. Finally, 5 ns MD simulations were conducted at 1 atm and 300 K under the NPT ensemble with a time integration step of 2 fs. Dynamic information (co-ordinates, velocities and forces) of inclusion complexes with integration step of 2 fs was preserved at every 1 ps in a trajectory file (.mdcrd) during the entire MD simulation. Analysis of MD trajectories generated was performed by Ptraj module in AMBERTOOL 12. Subsequently, to improve the accuracy of the MD results, the abovementioned complexes were optimized by using a two layered hybrid ONIOM method (ONIOM2). In ONIOM2 calculations, rocuronium and vecuronium were defined as the layer of high-level, whereas sugammadex was defined as the layer of low-level, which were calculated by B3LYP/6–31 + G(d, p) level of theory and PM3 method, respectively. The complexation energy (ΔE) between host and guest molecules was calculated according Eq. (1): Where Εcomplex, Εguest and Εhost stand for the total energy of the inclusion complex, the free energy of guest and host molecules, respectively. ΔΕ ¼ Εcomplex − Εhost þ Εguest
ð1Þ
The natural bond orbital (NBO) [69] calculation for the most stable complexes optimized by ONIOM2 method was performed to elucidate the intermolecular interaction between the host and guest molecules via the determination of the stabilization energy E(2), which was correlated with the delocalization tendency of electrons on the basis of perturbation method, and this method can be used to analyze the presence of hydrogen bonding. The harmonic frequency of the above stable complexes was calculated by using ONIOM2 method, which would contribute to analyze enthalpy (H), entropy (S) and Gibbs free energy (G). Furthermore, to understand the information about the region from where the host and guest molecules could have electrostatic
453
interaction, the molecular electrostatic potential (MEP) [70] was also calculated for the isolated guest and host molecules. Frontier molecular orbitals (FMOs) [71] were also applied to confirm the chemical stability of NMBs-sugammadex complexes. 3. Results and discussion 3.1. Molecular docking results Molecular docking was performed to establish the inclusion models of NMBs-sugammadex complexes initially. The predominant molecular docking configurations of rocuronium and vecuronium with sugammadex are shown in Fig. 1. Rocuronium and vecuronium possesses similar binding poses with sugammadex. The hydrophobic steroidal skeleton of NMBs penetrates into the hydrophobic cavity of sugammadex and the polar moieties are found to protrude outside the cyclodextrin rim; especially, the quaternary ammonium salt group closes to the side chain of sugammadex that is guarded by the carboxyl groups in a position that favors electrostatic interaction formation. The binding free energy of rocuronium and vecuronium inclusion complexes with sugammadex obtained from AutoDock-Vina is −11.2 and − 9.8 kcal mol−1, respectively, which confirms that rocuroniumsugammadex complex is more stable than that of vecuroniumsugammadex. 3.2. Molecular dynamics analysis Molecular docking is usually limited in the description of host-guest interactions, solvent effects, charge transfer, electrostatic effects, among others [58, 60]. In order to achieve accurate geometric structure of NMBs-sugammadex complexes in solution, we further ran MD simulation based on docked data fully accounting for the flexibility of the host-guest interactions. To monitor the system stability, the time
Fig. 5. RMSD plots for the simulation of Rocu-sugammadex (A) and Vecu-sugammadex (B).
454
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
Table 1 Hydrogen bond occupancy and distance (and standard deviation) calculated during the last 3 ns of the MD trajectories for Rocu-sugammadex. Donor
Acceptor
Occupancy%
Diatance/Å (SD)
OH (sugammadex)
OH (rocuronium)
19.6
2.56 (0.12)
dependence of the root mean square deviation (RMSD) of the MD simulation for all atomic position of NMBs-sugammadex inclusion complexes relative to those of the initial structures are illustrated in Fig. 5. The MD trajectory of Rocu-sugammadex inclusion complex rises sharply to 3.5 Å within 1 ns then averages at around 3.0 Å, while the RMSD of Vecu-sugammadex complex fluctuates to an average at 2.7 Å after 1.5 ns. Therefore, from the RMSD plots, it seems that all the simulated systems have reached equilibrium at 2 ns, and so the MDs trajectories from the last 3 ns were extracted for further analysis. We further investigated the intermolecular hydrogen bonding along the trajectories of the MD simulation for both complexes (See Table 1). It is imperative to mention here that the length of intermolecular hydrogen bonds in inclusion complexes is an indicator of the stability of these complexes. For Rocu-sugammadex complex, a significant hydrogen bond formed involves the secondary hydroxyl moiety of sugammadex and the oxygen atom of morpholino in Rocu acting as acceptor atom. These results suggest that the polar group of the rocuronium molecule is located at the entrance of CD and is participating in hydrogen bonding interaction. Surprisingly, no reliable hydrogen bonding is found in Vecusugammadex inclusion complex during the equilibrium stage of simulation indicating that a less stable inclusion is formed for vecuronium compared to rocuronium. This data is in line with the docked results of NMBs-sugammadex complexes in the gas phase. 3.3. Quantum mechanical simulation 3.3.1. Geometric structures of NMBs-sugammadex complexes QM-MM calculations analyze intermolecular interactions of hostguest system more accurately, allowing the energy values and intermolecular interactions between the host and guest molecule to be determined [53–56]. The minimum state conformers of NMBssugammadex complexes obtained by MD were further calculated by a multi-layered hybrid ONIOM method to acquire exact geometric structures. In the inclusion process, rocuronium and vecuronium prefer to enter from the wide side of sugammadex to assemble stable inclusion complexes. As can be seen from Fig. 6, the rocuronium-sugammadex and vcuronium-sugammadex inclusion complexes show similar binding modes comparing to molecular docking and dynamic results. The
Table 2 The E, ΔE and ΔΔE calculated by ONIOM2 methods for rocuronium-sugammadex and vecuronium-sugammadex complexes. Sugammadex complex
EONIOM2a (kcal mol−1)
ΔEONIOM2 (kJ mol−1)
ΔΔEONIOM2 (kJ mol−1)
Rocuronium Vecuronium
−899,170.33 −899,091.57
−34.61 −25.29
9.32
a
The total energy of inclusion complex.
most significant difference between docking and QM-MM results mainly locates on sugammadex, with CD ring somewhat puckering and carboxyl-containing alkyl side chains unordered extending. The mentioned discrepancy could be attributed to the neglection of the flexibility of the host as well as solvent and electrostatic effects (this interaction will be discussed in below) in molecular docking algorithm. Moreover, the refined simulated binding model of rocuroniumsugammadex complex demonstrates highly accordance with the reported X-ray crystal geometric structure. Thus, the three-step molecular modeling procedure (molecular docking + MD + QM-MM) is suitable for structural elucidating of NMBs-sugammadex encapsulation complexes. The performed NBO analyses of hydrogen bond system suggests that rocuronium-sugammadex inclusion complex is stabilized by a strong intermolecular hydrogen bond [LP(O36):BD × (O190 − H278)] between host and guest, verus vecuronium-sugammadex complex. The interaction energies are one of the important parameters to quantify the binding strength between host and guest molecule, which are summarized in Table 2. The total energies of the inclusion complexes are lower than the sum energies of free host and guest molecules, indicates the two complexes are both easily to be formed and thereby sugammadex can rapidly reverse their effect of neuromuscular blockade. Furthermore, the value of binding energies for rocuronium and vecuronium with sugammadex obtained by ONIOM2 method are respectively to −34.61 and −25.29 kJ mol−1. The binding energy difference (ΔΔE) between rocuronium-sugammadex and vecuroniumsugammadex complexes is −9.32 kJ mol−1, resulting the complexes formed by rocuronium is more stable than vecuronium, thus the latter neuromuscular blocker has lower binding affinity to sugammadex which is identical with the clinical trials. 3.3.2. Molecular electrostatic potential To verify whether the host and guest molecules can have electrostatic interactions, the molecular electrostatic potential (MEP) was calculated for the abovementioned isolated rocuronium, vecuronium and
Fig. 6. The detailed structures of inclusion complexes formed by rocuronium (A) and vecuronium (B) with sugammadex obtained by ONIOM2 method. Possible intermolecular hydrogen bonds are shown as a dotted line. The sugammadex ring is somewhat puckered, with the unordered carboxyl-containing alkyl side chains forming electrostatic interactions with NMBs.
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
455
Fig. 7. Molecular electrostatic potential maps of sugammadex (A) rocuronium (B) and vecuronium (C).
sugammadex. The importance of MEP lies in the fact that it simultaneously displays molecular size, shape and electrostatic potential value and is very useful in research of molecular structure with its physiochemical property relationship. In the 3D plots of MEP, it can be seen that the quaternary nitrogen atom is surrounded by a greater surface of positive charge, indicating a possible site for nucleophilic attack. The maximum values of the positive region on the above nitrogen atoms are determined as +92.65 and +81.13 kcal mol−1 for rocuronium and vecuronium, respectively. Furthermore, a maximum negative region localized on the carboxyls of sugammadex has an average value of −103.71 kcal mol−1, is the preferred site for electrophilic attack. Considering these calculated results, the MEP map shows that the electrostatic interaction between sugammadex and rocuronium/vecuronium could be occurred from the region of negative carboxyl group and positive quaternary ammonium moiety. In summary, hydrophobic interaction, hydrogen bonding, electrostatic interaction and host-guest size-fit concept are the essential interactions of NMBs-sugammadex inclusion complexes considering in drug design stage. Moreover, due to the electrostatic potential value of nitrogen atom in rocuronium (+92.65 kcal mol−1) is larger than vecuronium (+81.13 kcal mol−1), the electrostatic interaction between rocuronium and sugammadex is predicted to be stronger than the latter. (See Fig. 7.) 3.3.3. Frontier molecular orbitals analysis Frontier molecular orbitals parameters including the energy of highest occupied molecular orbital (HOMO), the lowest unoccupied molecular orbital (LUMO) and energy difference between HOMO and LUMO, can be used to determine some of the structural and physical properties, the chemical reactivity, the kinetic stability of molecules [71]. Therefore, so as to investigate the chemical stability of the above inclusion complexes, the energies of HOMO (EHOMO), LUMO (ELUMO) and their orbital energy gap (Egap) were calculated by using the above ONIOM2 method, and the results are listed in Table 3. The value of the energy separation between the HOMO and LUMO is −5.31 and −4.96 eV for rocuronium-sugammadex and vecuroniumsugammadex complexes, respectively. The larger HOMO-LUMO gap
Table 3 ELUMO, EHOMO and Egap of the NMBs-sugammadex complexes.
automatically means higher excitation energies of the excited states, more stability and a larger chemical hardness. Thus vecuronium has lower binding affinity to sugammadex than rocuronium, resulting in rocuronium is more effectively reversed by sugammadex.
3.3.4. Thermodynamic analysis To investigate the thermodynamics of the binding process, the statistical thermodynamic calculation was carried out at 1 atm and 298.15 K by ONIOM2 approach. The thermodynamic parameters, such as Enthalpy (H), Entropy (S) and Gibbs free energy (G) are the most important parameters in thermodynamic system, due to the changes of them can be used to confirm the exothermic or endothermic reaction, and predict the equilibrium and spontaneity. The thermodynamic quantities of the complexation process for sugammadex and the two NMBs, including Enthalpy change, Entropy change and Gibbs free energy change are summarized in Table 4. The complex formation of NMBs (rocuronium/vecuronium) with sugammadex is an exothermic and spontaneous reaction judged from the negative enthalpy change and Gibbs free energy change. The negative ΔH and positive ΔS values indicate that the formation of the above inclusion complexes is a both enthalpy and entropy favored process. The statistical thermodynamic parameters for the above calculated complexes reveal a determinacy of forming inclusion complex between host and guest molecules. Suitable less polar molecule NMBs inserts into the cavity of CD accompanying by formation of strong hydrophobic, hydrogen bonding and electrostatic interaction, and then remarkably lower the energy of the whole system. Furthermore, the Gibbs free energy change of the inclusion process is −9.47 and −7.36 kcal mol−1 for rocuronium and vecuronium, respectively, indicating that rocuronium is more favorable to dock into sugammadex than vecuronium. These results are accordance with those obtained by the theoretical procedures used in primary part.
Table 4 Statistical thermodynamic parameters of rocuronium and vecuronium complexed with sugammadex. Sugammadex complex
ΔHa (kcal mol−1)
ΔSb (cal mol−1 K−1)
ΔGc (kcal mol−1)
−6.73 −5.11
9.20 7.59
−9.47 −7.36
Sugammadex complex
ELUMO (eV)
EHOMO (eV)
Egapa (eV)
Rocuronium Vecuronium
Rocuronium Vecuronium
−0.89 −0.71
−6.32 −5.67
−5.31 −4.96
a
a
Egap = ELUMO − EHOMO.
b c
ΔH = Hcomplex − Hhost − Hguest. ΔS = Scomplex − Shost − Sguest. ΔG = ΔH − TΔS.
456
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
4. Conclusions An extensive molecular modeling was conducted here to elaborate the mechanism of action of sugammadex in reversing NMBs. The establishment of NMBs-sugammadex inclusion complexes was a spontaneous and enthalpy-driven process, which resulting in the hydrophobic steroidal skeleton of NMBs preferentially dipped to sugammadex, and polar moieties pointed to the cyclodextrin rims. The nature of bonding between host and guest molecules was investigated using MEP and NBO, indicating that theses complexes were stabilized by hydrogen bonding and electrostatic interaction. Furthermore, the reversal potency of rocuronium and vecuronium with sugammadex was also evaluated, which showed a good agreement with clinical trials, suggesting that rocuronium is more effectively reversed by sugammadex than vecuronium. This work highlights the antagonistic mechanism by sugammadex and essential interaction sites of NMBs-sugammadex inclusion complexes. The investigation data provides foundations to design novel γ-CD-based muscle relaxant antagonists, especially for a designated NMB; and the calculation model could be used to rapidly predict drugs' efficacy before pharmacological test. Acknowledgements The authors of this manuscript acknowledge the Innovative Research Team of the Ministry of Education and Program for Liaoning Innovative Research Team in University for providing financial assistance. The theoretical calculations were conducted on the ScGrid and Deepcomp 7000 the Supercomputing Center, Computer Network Information Center of Chinese Academy of Sciences. References [1] E.E. Boros, E.C. Bigham, G.E. Boswell, M.R. Jr, S.S. Patel, J.J. Savarese, J.A. Ray, J.B. Thompson, M.A. Hashim, J.C. Wisowaty, J. Med. Chem. 42 (1999) 206. [2] L.I. Chuan-Xiang, S.L. Yao, N. Hui, Y.Q. Zhang, L. Bin, Chin. Pharmacol. Bull. 24 (2005) 605. [3] D. Gusmão, N. Engl. J. Med. 363 (2010) 2562. [4] J.B. Hall, W. Schweickert, J.P. Kress, Crit. Care Med. 37 (2009) 416. [5] L. Hawryluck, J. Med. Ethics 28 (2002) 170. [6] J.M. Hunter, N. Engl. J. Med. 332 (1995) 1691. [7] B. Yegneswaran, R. Murugan, Crit. Care 15 (2011) 311. [8] E.S. Kilchevsky, J.T. Wung, Pediatrics 76 (1985) 653. [9] D.R. Bevan, C.E. Smith, F. Donati, Anesthesiology 69 (1988) 272. [10] T. Magorian, K.B. Flannery, R.D. Miller, Anesthesiology 79 (1993) 913. [11] A. Bom, J.K. Clark, R. Palin, Curr. Opin. Drug Discov. Devel. 5 (2002) 793. [12] M. Eikermann, H. Groeben, J. Hüsing, J. Peters, Anesthesiology 98 (2003) 1333. [13] M. Grosse-Sundrup, J.P. Henneman, W.S. Sandberg, B.T. Bateman, J.V. Uribe, N.T. Nguyen, J.M. Ehrenfeld, E.A. Martinez, T. Kurth, BMJ 345 (2017) 599. [14] G.S. Murphy, J.W. Szokol, M.J. Avram, S.B. Greenberg, T. Shear, J.S. Vender, J. Gray, E. Landry, Anesth. Analg. 117 (2013) 133. [15] G.S. Murphy, J.W. Szokol, J.H. Marymont, S.B. Greenberg, M.J. Avram, J.S. Vender, M. Nisman, Anesthesiology 109 (2008) 389. [16] I.F. Sorgenfrei, J. Viby-Mogensen, F.A. Swiatek, Ugeskr. Laeger 167 (2005) 3878. [17] J.M. Adam, D.J. Bennett, A. Bom, J.K. Clark, H. Feilden, E.J. Hutchinson, R. Palin, A. Prosser, D.C. Rees, G.M. Rosair, J. Med. Chem. 45 (2002) 1806. [18] A.S. Akha, J.S. Jahr, A. Li, K. Kiai, Anesthesiol. Clin. 28 (2010) 691. [19] A.H.A. Bom, A.W. Muir, D. Rees, US, 2007. [20] A. Gamaleldinmona, H. Macartneydonal, Can. J. Chem. 92 (2014) 243. [21] R.M.G. Hogg, R.K. Mirakhur, J. Anaesthesiol. Clin. Pharmacol. 25 (2009) 403. [22] D.H. Macartney, Future Med. Chem. 5 (2013) 2075.
[23] H.J. Shu, C.M. Zeng, C. Wang, D.F. Covey, C.F. Zorumski, S. Mennerick, Br. J. Pharmacol. 150 (2007) 164. [24] S. Ganapati, (2017). [25] G.J. Tarver, S.J.A. Grove, K. Buchanan, A. Bom, A. Cooke, S.J. Rutherford, M.Q. Zhang, Bioorg. Med. Chem. 10 (2002) 1819. [26] M. Zhang, R. Palin, D.J. Bennett, CA, 2005. [27] M. Zhang, G. Tarver, US, 2004. [28] M. Naguib, Anesth. Analg. 104 (2007) 575. [29] A.L. Kovac, J. Clin. Anesth. 21 (2009) 444. [30] S.J. Schaller, H. Fink, Core Evid. 8 (2013) 57. [31] L.P.H. Yang, S.J. Keam, Drugs 69 (2009) 919. [32] M.A. Zwiers, M.V.D. Heuve, J. Smeets, S. Rutherford, Clin. Drug Investig. 31 (2011) 101. [33] J. Raft, P. Guerci, V. Harter, T. Fuchs-Buder, C. Meistelman, Korean J. Anesthesiol. 68 (2015) 17. [34] A. Bom, M. Bradley, K. Cameron, J.K. Clark, E.J. Van, H. Feilden, E.J. Maclean, A.W. Muir, R. Palin, D.C. Rees, Angew. Chem. 41 (2002) 266. [35] K.S. Cameron, F. Dan, L. Fielding, Magn. Reson. Chem. 40 (2002). [36] Z. Zhou, L. Li, N. Yan, L. Du, C. Sun, T. Sun, Molecules 20 (2015) 20926. [37] L. Li, C. Wu, Z. Li, W. Xie, X. Guo, Z. Zhou, T. Sun, J. Mol. Struct. 1164 (2018) 493–500. [38] S. Basak, S. Mondal, S. Dey, P. Bhattacharya, A. Saha, P.V. Deep, A. Abbas, S.N. Gopal, J. Pharm. Pharmacol. 69 (2017) 1304–1317. [39] Y. Chu, Y. Li, Y. Tong, X. Li, Comput. Appl. Chem. 7 (2014) 17. [40] S. Vaishnav, G.K. Kulkarni, P. Perumal, Int. J. Pharm. Phytopharmacol. Res. (2013) 2. [41] J.Q. Zhang, K.M. Jiang, K. An, S.H. Ren, X.G. Xie, Y. Jin, J. Lin, Carbohydr. Res. 418 (2015) 20. [42] Y. Hasegawa, Y. Inoue, K. Deguchi, S. Ohki, M. Tansho, T. Shimizu, K. Yazawa, J. Phys. Chem. B 116 (2012) 1758. [43] S.S. Tallury, M.B. Smyth, E. Cakmak, M.A. Pasquinelli, J. Phys. Chem. B 116 (2012) 2023. [44] K. Boonyarattanakalin, P. Wolschann, P. Toochinda, L. Lawtrakul, Eur. J. Pharm. Sci. 47 (2012) 752. [45] M. Jana, S. Bandyopadhyay, J. Phys. Chem. B 117 (2013) 9280. [46] N.L. Pacioni, A.B. Pierini, A.V. Veglia, Spectrochim. Acta A Mol. Biomol. Spectrosc. 103 (2013) 319. [47] A. Zaboub, F. Madi, R. Merdes, M. Mohamedi, L. Nouar, J. Mol. Liq. 216 (2016) 716. [48] W. Y, W. X, L. N, J. Mol. Model. 20 (2014) 2431. [49] M. Dehestani, S. Pourestarabadi, L. Zeidabadinejad, Russ. J. Phys. Chem. A 90 (2016) 1192. [50] L. Yao, Y. Mori, K. Takano, Bull. Chem. Soc. Jpn. 87 (2014) 258. [51] A.L. Ivansyah, M.A. Martoprawiro, Buchari, 2017, p. 012070. [52] S. Benghodbane, D. Khatmi, J. Incl. Phenom. Mol. Recognit. Chem. 77 (2013) 231. [53] R. Djemil, D. Khatmi, J. Comput. Theor. Nanosci. 9 (2012) 1571. [54] B. Abdelaziz, N. Leila, H. Sakina, D. Imene, M. Fatiha, K. Djameleddine, J. Incl. Phenom. Mol. Recognit. Chem. 77 (2013) 455. [55] M. Fatiha, L. Leila, N. Leila, K.D. Eddine, J. Taiwan Inst. Chem. Eng. 43 (2012) 868. [56] L. Li, C. Wu, Y. Ma, S. Zhou, Z. Li, T. Sun, Analyst 142 (2017). [57] X. Guo, Z. Wang, L. Zuo, Z. Zhou, X. Guo, T. Sun, Analyst 139 (2014) 6511. [58] S.K. Tripathi, R. Muttineni, S.K. Singh, J. Theor. Biol. 334 (2013) 87. [59] S.Y. Lu, Y.J. Jiang, J. Lv, T.X. Wu, Q.S. Yu, W.L. Zhu, J. Mol. Graph. Model. 28 (2010) 766. [60] D.C. Rapaport, Comput. Sci. Eng. 1 (2002) 70. [61] M. Gholampour, N. Edraki, M. Khoshneviszadeh, Iranian Pharmaceutical Science Congress and Symposium of Biopharmaceutics and Pharmacokinetics, 2015. [62] J.F.M.G W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, (2009). [63] O. Trott, A.J. Olson, J. Comput. Chem. 31 (2010) 455. [64] D.A. Case, R. Anandakrishnan, K. Ayers, R. Betz, Y.J. Bomble, R.A. Brown, S. Brozell, Q. Cai, D.S. Cerutti, T.E.C. Iii, (2012). [65] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1998) 926. [66] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [67] H.J.C. Berendsen, J.P.M. Postma, W.F.V. Gunsteren, A. Dinola, J.R. Haak, J. Chem. Phys. 81 (1998) 3684. [68] U. Essmann, L. Perera, M. Berkowitz, T. Darden, H. Lee, L. Pedersen, S.H. Lee, (1995). [69] F. Weinhold, (2008). [70] J.S. Murray, Z. Peralta-Inga, P. Politzer, K. Ekanayake, P. Lebreton, Int. J. Quantum Chem. 83 (2001) 245. [71] K. Fukui, The Role of Frontier Orbitals in Chemical Reactions, 1982.
Contents lists available at ScienceDirect
Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq
Theoretical studies on the mechanism of sugammadex for the reversal of aminosteroid-induced neuromuscular blockade Linwei Li a, Yanan Zhou a, Zhengjun Wang a, Chengjun Wu a, Zhen Li a, Changshan Sun b,⁎, Tiemin Sun a,⁎ a b
Key Laboratory of Structure-Based Drug Design and Discovery, Ministry of Education, Shenyang Pharmaceutical University, Shenyang 110016, PR China Pharmacy Department, Shenyang Pharmaceutical University, Shenyang 110016, PR China
a r t i c l e
i n f o
Article history: Received 3 May 2018 Received in revised form 1 June 2018 Accepted 8 June 2018 Available online 15 June 2018 Keywords: Sugammadex Neuromuscular blockers Mechanism DFT MD
a b s t r a c t The mechanism of action of sugammadex in reversing neuromuscular blockers (NMBs) was theoretically studied by extensive molecular modeling (molecular docking, molecular dynamic and molecular mechanics-quantum mechanics). The NMBs-sugammadex inclusion complex was established by an enthalpy-driven process, which stabilized by hydrophobic interaction, hydrogen bonding and electrostatic interaction between NMBs and sugammadex. Furthermore, the binding energy of NMBs-sugammadex indicated that vecuronium (ΔE = −25.29 kJ mol−1) had lower binding affinity to sugammadex than rocuronium (ΔE = −34.61 kJ mol−1), thus the former one would be more difficult to be reversed which was consist with the phenotype of clinical trials. Molecular electrostatic potential, natural bond orbital and frontier molecular orbitals were also utilized to give explanation for the discrepancy of interaction strength between rocuronium-sugammadex and vecuroniumsugammadex inclusion complexes. In summary, the molecular modeling study made it possible to explain the different reversal efficiency of NMBs by sugammadex, which could be used to predict the drug-receptor interactions during the design phase of novel compounds development and also provide great opportunity to design novel γ-cyclodextrin-based muscle relaxant antagonists. © 2018 Published by Elsevier B.V.
1. Introduction Neuromuscular blockers (NMBs) are routinely used for facilitating surgical procedures and tracheal intubation during anesthesia. They usually bind to nicotinic acetylcholine receptor (nAChR) at the neuromuscular junction, and then block neurotransmission induced by acetylcholine (ACh) [1–7]. One such typical and widely used species is the non-depolarizing aminosteroid-based NMBs, such as pancuronium (proved in 1973) [8], vecuronium (1982) [9], and rocuronium (1994) [10], which behaviors fast onset timers and short duration of action. Although aminosteroid-based NMBs have few adverse effects during anesthesia, the residual duration of muscle relaxants beyond the end of the surgery, also referred to as postoperative residual neuromuscular blockers (RNMBs), is a well-known problem [11–16]. The postoperative residual neuromuscular blockers associate with respiratory insufficiency, pulmonary complications, upper airway obstruction and hypoxia, which can increase the incidence of tracheal re-intubation and mortality [17–23]. Therefore, it is necessary to exploit rapid and complete agents for controlling the duration of the residual neuromuscular blockers at the end of operation, which be used to accelerate the recovery of patients' muscle function. ⁎ Corresponding authors. E-mail addresses: [email protected] (C. Sun), [email protected] (T. Sun).
https://doi.org/10.1016/j.molliq.2018.06.033 0167-7322/© 2018 Published by Elsevier B.V.
Chemical encapsulation (Fig. 1) is a recent mechanism to remove the RNMBs drugs; and there have been investigations into exogenous macrocyclic host molecules, such as cyclodextrins (CDs) and anionic cyclophanes, as innovative and potential reversal agents [19, 24–27]. Inclusion complexes can be formed by aminosteroidal NMBs (guest molecule) and chemically modified γ-cyclodextrin (host molecule) with high binding affinity, resulting in a decrease of NMBs in plasma, then rapidly reverse the effect of residual neuromuscular blockers. Furthermore, γ-cyclodextrin derivatives (γ-CDs) do not bind to the muscarinic receptors, which lead its reversal effect more reliable than acetylcholinesterase inhibitors. Therefore, a series of modified γ-CDs have been synthesized with the aim of constructing host molecules capable of forming host-guest complexes with neuromuscular blockers [17, 25, 28]. Among them, sugammadex (Bridion®, octakis-(6-deoxy-6-Smercaptopropionyl)-γ-CD sodium salt) is demonstrated to be the first potent and selective relaxant binding agent, exerts its effect by forming very tight water-soluble inclusion complex at a 1: 1 ratio with aminosteroidal NMBs, approved by FDA at 2015 [28–32]. (See Figs. 2–4.) It was reported that the binding constant for vecuroniumsugammadex complex is somewhat smaller (Ki = 5.72 × 106 M−1) than rocuronium (Ki = 1.79 × 107 M−1) inclusion complex, while still weaker binding to pancuronium (Ki = 2.62 × 106 M−1) is detected [32]. This, in turn, lead to the suggestion that rocuronium is removed
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
451
Fig. 1. The mechanism of action of γ-cyclodextrin derivatives on the nAChR.
more effectively by sugammadex than the other aminosteroids, such as vecuronium and pancuronium. Further, in order to achieve equivalent reversal potency for low binding affinity NMBs, excessive sugammadex would be employed in clinical, and therefore, it could increase the possibility of bleeding risk induced by sugammadex [33]. Thus, it is better to exploit new potent γ-CD-based muscle relaxant antagonists, especially, screen a specific γ-CD derivative for a given aminosteroidal NMB. Different interaction modes and strengths between host and various guest molecules are attributing to distinct reversal efficiency. Molecular recognition mechanism is the prerequisite foundation to learn about the binding mode of inclusion complex, which would obviously help in designing new synthetic γ-CDs. Considering these concepts, we attempted to give explanations of the exact nature of reversal divergence in reversal process based on reversal mechanism at molecular level. Up to data, among NMBs-sugammadex complexes, only the geometric structure of rocuronium-sugammadex inclusion complex has been established by X-ray diffraction (XRD) [34] and NMR technique [35]. Additionally, the growing condition of single crystals and the limitation of XRD make XRD method infeasible [36, 37]; while NMR can't supply the exact distances of bonds in encapsulation complex. Therefore, the development of novel γ-CD-based muscle relaxant antagonists is still retard; and the evaluation of reversal potency is mainly based on pharmacological test. Luckily, molecular modeling methods have been recently proposed as powerful tools to unravel the interaction between
host and guest molecules and assist in elucidating the reversal process, and then evaluate the reversal efficacy. The methods involve molecular mechanics (MM), molecular docking [38–41], molecular dynamics (MD) [42–46], quantum mechanics (QM) [47–50] and molecular mechanics-quantum mechanics (QM-MM) calculations [51–57]. Molecular docking could rapidly provide the binding model of guest molecule with host molecule at low molecular mechanics level [58, 59]. Molecular dynamics simulation is employed which provides the equilibrium and transport properties of a classical many-body system in liquid circumstance. The major advantage of molecular dynamics is the simulated conformations that are less dependent on the initial structure than those obtained by molecular mechanic [58–61]. Thus, the MD results are very similar to real experiments. A multi-layered hybrid approach, our own N-layered integrated molecular orbital and molecular mechanics method (ONIOM), is capturing for macromolecular system since it can simultaneously treat different parts of a system with multiple approaches of varying accuracy and computational cost (QM and MM). In summary, the combination of the mentioned theoretical methods (molecular docking, MD and QM-MM) offers successful and rational data in solving structural, energetic and dynamic problems. Herein, with the aid of computational simulation, we aim to study the reversal mechanism of sugammadex with NMBs and explain the reversal efficacy among different NMBs at molecular level. Rocuronium (Rocu) and vecuronium (Vecu) were selected as the examples to
Fig. 2. The chemical structure and schematic diagram of sugammadex.
452
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
Fig. 3. The chemical structures of rocuronium and vecuronium.
perform our investigation. Molecular docking, MD, QM-MM methods were simultaneously utilized to establish the energetically most favorable encapsulation complexes of NMBs with sugammadex. Subsequently, QNIOM, natural bond orbital (NBO), molecular electrostatic potential (MEP) and frontier molecular orbitals (FMOs) were conducted to evaluate binding efficacy of inclusion complexes from electrostatic interaction, hydrogen bonding and chemical stability aspects. Finally, thermodynamic parameters were calculated to describe encapsulation process of NMBs-sugammadex complexes along with the involved driving force. This work presents a deep insight into the antagonistic mechanism by sugammadex, which can be adopted to design novel γ-CDbased muscle relaxant antagonists for a given steroidal neuromuscular blocker, as well as can be used to rapidly predict potency before pharmacological test.
2. Computational methods The initial structures of rocuronium and vecuronium depicted by ChemBioDraw Ultra 13.0 were imported to the Spartan 14.0 for conformational searching by using the MM+ force field method. Then selected a series of low-energy conformers that were in equilibrium with each other, and their relative energies were within a region of 4 kcal mol−1. To confirm the more stable structures of guest molecules, the above conformers were further calculated at the density functional theory (DFT) level with a Becke-3-Lee-Yang-Parr (B3LYP) functional and 6–31 + G(d, p) basis set by using Gaussian 09 [62] suite of programs. Then redistributed energy values according to the Boltzmann distribution, permitting decline to the most stable structures of isolated rocuronium and vecuronium for further simulations. The γcyclodextrin structure was obtained from the crystallographic parameters provided by the Cambridge crystallographic Data center, then based on the initial structure, sugammadex was constructed and then fully
optimized at the parameterized model 3 (PM3) level of theory without imposing any symmetrical. Molecular docking studies were undertaken using AutoDock-Vina program [63]. rocuronium and vecuronium were allowed to dock anywhere in a 35 Å × 35 Å × 35 Å grid box with 1 Å spacing centered the center of mass of sugammadex, with the Lamarckian genetic algorithm (LGA). The initial torsions and positions of rocuronium and vecuronium were generated randomly. The maximum numbers of energy evaluation and generation were set to 250,000 and 27,000, respectively. A total of 20 separate docking runs were performed with the initial population of 20 individuals, and the lowest energy structure from each run was retained. Final docked conformations were clustered by using a tolerance of 1.0 Å root mean square deviation (RMSD). AMBERTOOL 12 [64] software package was used for the molecular dynamic (MD) simulations using the general amber force field (GAFF) parameter set [65]. Topology and parameter files were generated with the LEaP program on structure of the complexes obtained by the aforementioned docking procedures. The Sander module of AMBERTOOL 12 package was used to perform MD simulations. Each system was solvated in a truncated octahedral box of TIP3P water molecules with a minimum solute-water distance of 10 Å. All bonds involving hydrogen atoms were constrained using the SHAKE flag [66]. Long-range electrostatic interactions were calculated using the particle-mesh Ewald (PME) algorithm [67], with a cutoff of 12 Å for Lennard-Jones interactions. Periodic boundary conditions were applied to avoid edge effects. Prior to MD production, the steepest-descent method was used for the 500 iterations and the conjugated gradient method for the subsequent 2500 iteration in each step for the solvent and the entire model system, respectively. Afterward, the entire system was gradually heated up to 300 K gradually over 50 ps using the constant volume and normal temperature (NVT) [68] ensemble with the solutes restrained by a weak harmonic potential. During the heating, time constant for heat bath coupling for the solute was set as 1.0. Subsequently, the constant normal
Fig. 4. Docked geometries of energetically most favorable inclusion complexes of rocuronium (A) and vecuronium (B) with sugammadex. The polar parts of NMBs protrude outside the cavity from the opening of the secondary OH side, and all four steroidal rings are embodied within the hydrophobic cavity of sugammadex extended by side chains.
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
pressure and normal temperature (NPT) ensemble was carried out for 500 ps, keeping a temperature of 300 K by Berendsen thermostat, and maintaining pressure constant at 1.0 atm. Finally, 5 ns MD simulations were conducted at 1 atm and 300 K under the NPT ensemble with a time integration step of 2 fs. Dynamic information (co-ordinates, velocities and forces) of inclusion complexes with integration step of 2 fs was preserved at every 1 ps in a trajectory file (.mdcrd) during the entire MD simulation. Analysis of MD trajectories generated was performed by Ptraj module in AMBERTOOL 12. Subsequently, to improve the accuracy of the MD results, the abovementioned complexes were optimized by using a two layered hybrid ONIOM method (ONIOM2). In ONIOM2 calculations, rocuronium and vecuronium were defined as the layer of high-level, whereas sugammadex was defined as the layer of low-level, which were calculated by B3LYP/6–31 + G(d, p) level of theory and PM3 method, respectively. The complexation energy (ΔE) between host and guest molecules was calculated according Eq. (1): Where Εcomplex, Εguest and Εhost stand for the total energy of the inclusion complex, the free energy of guest and host molecules, respectively. ΔΕ ¼ Εcomplex − Εhost þ Εguest
ð1Þ
The natural bond orbital (NBO) [69] calculation for the most stable complexes optimized by ONIOM2 method was performed to elucidate the intermolecular interaction between the host and guest molecules via the determination of the stabilization energy E(2), which was correlated with the delocalization tendency of electrons on the basis of perturbation method, and this method can be used to analyze the presence of hydrogen bonding. The harmonic frequency of the above stable complexes was calculated by using ONIOM2 method, which would contribute to analyze enthalpy (H), entropy (S) and Gibbs free energy (G). Furthermore, to understand the information about the region from where the host and guest molecules could have electrostatic
453
interaction, the molecular electrostatic potential (MEP) [70] was also calculated for the isolated guest and host molecules. Frontier molecular orbitals (FMOs) [71] were also applied to confirm the chemical stability of NMBs-sugammadex complexes. 3. Results and discussion 3.1. Molecular docking results Molecular docking was performed to establish the inclusion models of NMBs-sugammadex complexes initially. The predominant molecular docking configurations of rocuronium and vecuronium with sugammadex are shown in Fig. 1. Rocuronium and vecuronium possesses similar binding poses with sugammadex. The hydrophobic steroidal skeleton of NMBs penetrates into the hydrophobic cavity of sugammadex and the polar moieties are found to protrude outside the cyclodextrin rim; especially, the quaternary ammonium salt group closes to the side chain of sugammadex that is guarded by the carboxyl groups in a position that favors electrostatic interaction formation. The binding free energy of rocuronium and vecuronium inclusion complexes with sugammadex obtained from AutoDock-Vina is −11.2 and − 9.8 kcal mol−1, respectively, which confirms that rocuroniumsugammadex complex is more stable than that of vecuroniumsugammadex. 3.2. Molecular dynamics analysis Molecular docking is usually limited in the description of host-guest interactions, solvent effects, charge transfer, electrostatic effects, among others [58, 60]. In order to achieve accurate geometric structure of NMBs-sugammadex complexes in solution, we further ran MD simulation based on docked data fully accounting for the flexibility of the host-guest interactions. To monitor the system stability, the time
Fig. 5. RMSD plots for the simulation of Rocu-sugammadex (A) and Vecu-sugammadex (B).
454
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
Table 1 Hydrogen bond occupancy and distance (and standard deviation) calculated during the last 3 ns of the MD trajectories for Rocu-sugammadex. Donor
Acceptor
Occupancy%
Diatance/Å (SD)
OH (sugammadex)
OH (rocuronium)
19.6
2.56 (0.12)
dependence of the root mean square deviation (RMSD) of the MD simulation for all atomic position of NMBs-sugammadex inclusion complexes relative to those of the initial structures are illustrated in Fig. 5. The MD trajectory of Rocu-sugammadex inclusion complex rises sharply to 3.5 Å within 1 ns then averages at around 3.0 Å, while the RMSD of Vecu-sugammadex complex fluctuates to an average at 2.7 Å after 1.5 ns. Therefore, from the RMSD plots, it seems that all the simulated systems have reached equilibrium at 2 ns, and so the MDs trajectories from the last 3 ns were extracted for further analysis. We further investigated the intermolecular hydrogen bonding along the trajectories of the MD simulation for both complexes (See Table 1). It is imperative to mention here that the length of intermolecular hydrogen bonds in inclusion complexes is an indicator of the stability of these complexes. For Rocu-sugammadex complex, a significant hydrogen bond formed involves the secondary hydroxyl moiety of sugammadex and the oxygen atom of morpholino in Rocu acting as acceptor atom. These results suggest that the polar group of the rocuronium molecule is located at the entrance of CD and is participating in hydrogen bonding interaction. Surprisingly, no reliable hydrogen bonding is found in Vecusugammadex inclusion complex during the equilibrium stage of simulation indicating that a less stable inclusion is formed for vecuronium compared to rocuronium. This data is in line with the docked results of NMBs-sugammadex complexes in the gas phase. 3.3. Quantum mechanical simulation 3.3.1. Geometric structures of NMBs-sugammadex complexes QM-MM calculations analyze intermolecular interactions of hostguest system more accurately, allowing the energy values and intermolecular interactions between the host and guest molecule to be determined [53–56]. The minimum state conformers of NMBssugammadex complexes obtained by MD were further calculated by a multi-layered hybrid ONIOM method to acquire exact geometric structures. In the inclusion process, rocuronium and vecuronium prefer to enter from the wide side of sugammadex to assemble stable inclusion complexes. As can be seen from Fig. 6, the rocuronium-sugammadex and vcuronium-sugammadex inclusion complexes show similar binding modes comparing to molecular docking and dynamic results. The
Table 2 The E, ΔE and ΔΔE calculated by ONIOM2 methods for rocuronium-sugammadex and vecuronium-sugammadex complexes. Sugammadex complex
EONIOM2a (kcal mol−1)
ΔEONIOM2 (kJ mol−1)
ΔΔEONIOM2 (kJ mol−1)
Rocuronium Vecuronium
−899,170.33 −899,091.57
−34.61 −25.29
9.32
a
The total energy of inclusion complex.
most significant difference between docking and QM-MM results mainly locates on sugammadex, with CD ring somewhat puckering and carboxyl-containing alkyl side chains unordered extending. The mentioned discrepancy could be attributed to the neglection of the flexibility of the host as well as solvent and electrostatic effects (this interaction will be discussed in below) in molecular docking algorithm. Moreover, the refined simulated binding model of rocuroniumsugammadex complex demonstrates highly accordance with the reported X-ray crystal geometric structure. Thus, the three-step molecular modeling procedure (molecular docking + MD + QM-MM) is suitable for structural elucidating of NMBs-sugammadex encapsulation complexes. The performed NBO analyses of hydrogen bond system suggests that rocuronium-sugammadex inclusion complex is stabilized by a strong intermolecular hydrogen bond [LP(O36):BD × (O190 − H278)] between host and guest, verus vecuronium-sugammadex complex. The interaction energies are one of the important parameters to quantify the binding strength between host and guest molecule, which are summarized in Table 2. The total energies of the inclusion complexes are lower than the sum energies of free host and guest molecules, indicates the two complexes are both easily to be formed and thereby sugammadex can rapidly reverse their effect of neuromuscular blockade. Furthermore, the value of binding energies for rocuronium and vecuronium with sugammadex obtained by ONIOM2 method are respectively to −34.61 and −25.29 kJ mol−1. The binding energy difference (ΔΔE) between rocuronium-sugammadex and vecuroniumsugammadex complexes is −9.32 kJ mol−1, resulting the complexes formed by rocuronium is more stable than vecuronium, thus the latter neuromuscular blocker has lower binding affinity to sugammadex which is identical with the clinical trials. 3.3.2. Molecular electrostatic potential To verify whether the host and guest molecules can have electrostatic interactions, the molecular electrostatic potential (MEP) was calculated for the abovementioned isolated rocuronium, vecuronium and
Fig. 6. The detailed structures of inclusion complexes formed by rocuronium (A) and vecuronium (B) with sugammadex obtained by ONIOM2 method. Possible intermolecular hydrogen bonds are shown as a dotted line. The sugammadex ring is somewhat puckered, with the unordered carboxyl-containing alkyl side chains forming electrostatic interactions with NMBs.
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
455
Fig. 7. Molecular electrostatic potential maps of sugammadex (A) rocuronium (B) and vecuronium (C).
sugammadex. The importance of MEP lies in the fact that it simultaneously displays molecular size, shape and electrostatic potential value and is very useful in research of molecular structure with its physiochemical property relationship. In the 3D plots of MEP, it can be seen that the quaternary nitrogen atom is surrounded by a greater surface of positive charge, indicating a possible site for nucleophilic attack. The maximum values of the positive region on the above nitrogen atoms are determined as +92.65 and +81.13 kcal mol−1 for rocuronium and vecuronium, respectively. Furthermore, a maximum negative region localized on the carboxyls of sugammadex has an average value of −103.71 kcal mol−1, is the preferred site for electrophilic attack. Considering these calculated results, the MEP map shows that the electrostatic interaction between sugammadex and rocuronium/vecuronium could be occurred from the region of negative carboxyl group and positive quaternary ammonium moiety. In summary, hydrophobic interaction, hydrogen bonding, electrostatic interaction and host-guest size-fit concept are the essential interactions of NMBs-sugammadex inclusion complexes considering in drug design stage. Moreover, due to the electrostatic potential value of nitrogen atom in rocuronium (+92.65 kcal mol−1) is larger than vecuronium (+81.13 kcal mol−1), the electrostatic interaction between rocuronium and sugammadex is predicted to be stronger than the latter. (See Fig. 7.) 3.3.3. Frontier molecular orbitals analysis Frontier molecular orbitals parameters including the energy of highest occupied molecular orbital (HOMO), the lowest unoccupied molecular orbital (LUMO) and energy difference between HOMO and LUMO, can be used to determine some of the structural and physical properties, the chemical reactivity, the kinetic stability of molecules [71]. Therefore, so as to investigate the chemical stability of the above inclusion complexes, the energies of HOMO (EHOMO), LUMO (ELUMO) and their orbital energy gap (Egap) were calculated by using the above ONIOM2 method, and the results are listed in Table 3. The value of the energy separation between the HOMO and LUMO is −5.31 and −4.96 eV for rocuronium-sugammadex and vecuroniumsugammadex complexes, respectively. The larger HOMO-LUMO gap
Table 3 ELUMO, EHOMO and Egap of the NMBs-sugammadex complexes.
automatically means higher excitation energies of the excited states, more stability and a larger chemical hardness. Thus vecuronium has lower binding affinity to sugammadex than rocuronium, resulting in rocuronium is more effectively reversed by sugammadex.
3.3.4. Thermodynamic analysis To investigate the thermodynamics of the binding process, the statistical thermodynamic calculation was carried out at 1 atm and 298.15 K by ONIOM2 approach. The thermodynamic parameters, such as Enthalpy (H), Entropy (S) and Gibbs free energy (G) are the most important parameters in thermodynamic system, due to the changes of them can be used to confirm the exothermic or endothermic reaction, and predict the equilibrium and spontaneity. The thermodynamic quantities of the complexation process for sugammadex and the two NMBs, including Enthalpy change, Entropy change and Gibbs free energy change are summarized in Table 4. The complex formation of NMBs (rocuronium/vecuronium) with sugammadex is an exothermic and spontaneous reaction judged from the negative enthalpy change and Gibbs free energy change. The negative ΔH and positive ΔS values indicate that the formation of the above inclusion complexes is a both enthalpy and entropy favored process. The statistical thermodynamic parameters for the above calculated complexes reveal a determinacy of forming inclusion complex between host and guest molecules. Suitable less polar molecule NMBs inserts into the cavity of CD accompanying by formation of strong hydrophobic, hydrogen bonding and electrostatic interaction, and then remarkably lower the energy of the whole system. Furthermore, the Gibbs free energy change of the inclusion process is −9.47 and −7.36 kcal mol−1 for rocuronium and vecuronium, respectively, indicating that rocuronium is more favorable to dock into sugammadex than vecuronium. These results are accordance with those obtained by the theoretical procedures used in primary part.
Table 4 Statistical thermodynamic parameters of rocuronium and vecuronium complexed with sugammadex. Sugammadex complex
ΔHa (kcal mol−1)
ΔSb (cal mol−1 K−1)
ΔGc (kcal mol−1)
−6.73 −5.11
9.20 7.59
−9.47 −7.36
Sugammadex complex
ELUMO (eV)
EHOMO (eV)
Egapa (eV)
Rocuronium Vecuronium
Rocuronium Vecuronium
−0.89 −0.71
−6.32 −5.67
−5.31 −4.96
a
a
Egap = ELUMO − EHOMO.
b c
ΔH = Hcomplex − Hhost − Hguest. ΔS = Scomplex − Shost − Sguest. ΔG = ΔH − TΔS.
456
L. Li et al. / Journal of Molecular Liquids 265 (2018) 450–456
4. Conclusions An extensive molecular modeling was conducted here to elaborate the mechanism of action of sugammadex in reversing NMBs. The establishment of NMBs-sugammadex inclusion complexes was a spontaneous and enthalpy-driven process, which resulting in the hydrophobic steroidal skeleton of NMBs preferentially dipped to sugammadex, and polar moieties pointed to the cyclodextrin rims. The nature of bonding between host and guest molecules was investigated using MEP and NBO, indicating that theses complexes were stabilized by hydrogen bonding and electrostatic interaction. Furthermore, the reversal potency of rocuronium and vecuronium with sugammadex was also evaluated, which showed a good agreement with clinical trials, suggesting that rocuronium is more effectively reversed by sugammadex than vecuronium. This work highlights the antagonistic mechanism by sugammadex and essential interaction sites of NMBs-sugammadex inclusion complexes. The investigation data provides foundations to design novel γ-CD-based muscle relaxant antagonists, especially for a designated NMB; and the calculation model could be used to rapidly predict drugs' efficacy before pharmacological test. Acknowledgements The authors of this manuscript acknowledge the Innovative Research Team of the Ministry of Education and Program for Liaoning Innovative Research Team in University for providing financial assistance. The theoretical calculations were conducted on the ScGrid and Deepcomp 7000 the Supercomputing Center, Computer Network Information Center of Chinese Academy of Sciences. References [1] E.E. Boros, E.C. Bigham, G.E. Boswell, M.R. Jr, S.S. Patel, J.J. Savarese, J.A. Ray, J.B. Thompson, M.A. Hashim, J.C. Wisowaty, J. Med. Chem. 42 (1999) 206. [2] L.I. Chuan-Xiang, S.L. Yao, N. Hui, Y.Q. Zhang, L. Bin, Chin. Pharmacol. Bull. 24 (2005) 605. [3] D. Gusmão, N. Engl. J. Med. 363 (2010) 2562. [4] J.B. Hall, W. Schweickert, J.P. Kress, Crit. Care Med. 37 (2009) 416. [5] L. Hawryluck, J. Med. Ethics 28 (2002) 170. [6] J.M. Hunter, N. Engl. J. Med. 332 (1995) 1691. [7] B. Yegneswaran, R. Murugan, Crit. Care 15 (2011) 311. [8] E.S. Kilchevsky, J.T. Wung, Pediatrics 76 (1985) 653. [9] D.R. Bevan, C.E. Smith, F. Donati, Anesthesiology 69 (1988) 272. [10] T. Magorian, K.B. Flannery, R.D. Miller, Anesthesiology 79 (1993) 913. [11] A. Bom, J.K. Clark, R. Palin, Curr. Opin. Drug Discov. Devel. 5 (2002) 793. [12] M. Eikermann, H. Groeben, J. Hüsing, J. Peters, Anesthesiology 98 (2003) 1333. [13] M. Grosse-Sundrup, J.P. Henneman, W.S. Sandberg, B.T. Bateman, J.V. Uribe, N.T. Nguyen, J.M. Ehrenfeld, E.A. Martinez, T. Kurth, BMJ 345 (2017) 599. [14] G.S. Murphy, J.W. Szokol, M.J. Avram, S.B. Greenberg, T. Shear, J.S. Vender, J. Gray, E. Landry, Anesth. Analg. 117 (2013) 133. [15] G.S. Murphy, J.W. Szokol, J.H. Marymont, S.B. Greenberg, M.J. Avram, J.S. Vender, M. Nisman, Anesthesiology 109 (2008) 389. [16] I.F. Sorgenfrei, J. Viby-Mogensen, F.A. Swiatek, Ugeskr. Laeger 167 (2005) 3878. [17] J.M. Adam, D.J. Bennett, A. Bom, J.K. Clark, H. Feilden, E.J. Hutchinson, R. Palin, A. Prosser, D.C. Rees, G.M. Rosair, J. Med. Chem. 45 (2002) 1806. [18] A.S. Akha, J.S. Jahr, A. Li, K. Kiai, Anesthesiol. Clin. 28 (2010) 691. [19] A.H.A. Bom, A.W. Muir, D. Rees, US, 2007. [20] A. Gamaleldinmona, H. Macartneydonal, Can. J. Chem. 92 (2014) 243. [21] R.M.G. Hogg, R.K. Mirakhur, J. Anaesthesiol. Clin. Pharmacol. 25 (2009) 403. [22] D.H. Macartney, Future Med. Chem. 5 (2013) 2075.
[23] H.J. Shu, C.M. Zeng, C. Wang, D.F. Covey, C.F. Zorumski, S. Mennerick, Br. J. Pharmacol. 150 (2007) 164. [24] S. Ganapati, (2017). [25] G.J. Tarver, S.J.A. Grove, K. Buchanan, A. Bom, A. Cooke, S.J. Rutherford, M.Q. Zhang, Bioorg. Med. Chem. 10 (2002) 1819. [26] M. Zhang, R. Palin, D.J. Bennett, CA, 2005. [27] M. Zhang, G. Tarver, US, 2004. [28] M. Naguib, Anesth. Analg. 104 (2007) 575. [29] A.L. Kovac, J. Clin. Anesth. 21 (2009) 444. [30] S.J. Schaller, H. Fink, Core Evid. 8 (2013) 57. [31] L.P.H. Yang, S.J. Keam, Drugs 69 (2009) 919. [32] M.A. Zwiers, M.V.D. Heuve, J. Smeets, S. Rutherford, Clin. Drug Investig. 31 (2011) 101. [33] J. Raft, P. Guerci, V. Harter, T. Fuchs-Buder, C. Meistelman, Korean J. Anesthesiol. 68 (2015) 17. [34] A. Bom, M. Bradley, K. Cameron, J.K. Clark, E.J. Van, H. Feilden, E.J. Maclean, A.W. Muir, R. Palin, D.C. Rees, Angew. Chem. 41 (2002) 266. [35] K.S. Cameron, F. Dan, L. Fielding, Magn. Reson. Chem. 40 (2002). [36] Z. Zhou, L. Li, N. Yan, L. Du, C. Sun, T. Sun, Molecules 20 (2015) 20926. [37] L. Li, C. Wu, Z. Li, W. Xie, X. Guo, Z. Zhou, T. Sun, J. Mol. Struct. 1164 (2018) 493–500. [38] S. Basak, S. Mondal, S. Dey, P. Bhattacharya, A. Saha, P.V. Deep, A. Abbas, S.N. Gopal, J. Pharm. Pharmacol. 69 (2017) 1304–1317. [39] Y. Chu, Y. Li, Y. Tong, X. Li, Comput. Appl. Chem. 7 (2014) 17. [40] S. Vaishnav, G.K. Kulkarni, P. Perumal, Int. J. Pharm. Phytopharmacol. Res. (2013) 2. [41] J.Q. Zhang, K.M. Jiang, K. An, S.H. Ren, X.G. Xie, Y. Jin, J. Lin, Carbohydr. Res. 418 (2015) 20. [42] Y. Hasegawa, Y. Inoue, K. Deguchi, S. Ohki, M. Tansho, T. Shimizu, K. Yazawa, J. Phys. Chem. B 116 (2012) 1758. [43] S.S. Tallury, M.B. Smyth, E. Cakmak, M.A. Pasquinelli, J. Phys. Chem. B 116 (2012) 2023. [44] K. Boonyarattanakalin, P. Wolschann, P. Toochinda, L. Lawtrakul, Eur. J. Pharm. Sci. 47 (2012) 752. [45] M. Jana, S. Bandyopadhyay, J. Phys. Chem. B 117 (2013) 9280. [46] N.L. Pacioni, A.B. Pierini, A.V. Veglia, Spectrochim. Acta A Mol. Biomol. Spectrosc. 103 (2013) 319. [47] A. Zaboub, F. Madi, R. Merdes, M. Mohamedi, L. Nouar, J. Mol. Liq. 216 (2016) 716. [48] W. Y, W. X, L. N, J. Mol. Model. 20 (2014) 2431. [49] M. Dehestani, S. Pourestarabadi, L. Zeidabadinejad, Russ. J. Phys. Chem. A 90 (2016) 1192. [50] L. Yao, Y. Mori, K. Takano, Bull. Chem. Soc. Jpn. 87 (2014) 258. [51] A.L. Ivansyah, M.A. Martoprawiro, Buchari, 2017, p. 012070. [52] S. Benghodbane, D. Khatmi, J. Incl. Phenom. Mol. Recognit. Chem. 77 (2013) 231. [53] R. Djemil, D. Khatmi, J. Comput. Theor. Nanosci. 9 (2012) 1571. [54] B. Abdelaziz, N. Leila, H. Sakina, D. Imene, M. Fatiha, K. Djameleddine, J. Incl. Phenom. Mol. Recognit. Chem. 77 (2013) 455. [55] M. Fatiha, L. Leila, N. Leila, K.D. Eddine, J. Taiwan Inst. Chem. Eng. 43 (2012) 868. [56] L. Li, C. Wu, Y. Ma, S. Zhou, Z. Li, T. Sun, Analyst 142 (2017). [57] X. Guo, Z. Wang, L. Zuo, Z. Zhou, X. Guo, T. Sun, Analyst 139 (2014) 6511. [58] S.K. Tripathi, R. Muttineni, S.K. Singh, J. Theor. Biol. 334 (2013) 87. [59] S.Y. Lu, Y.J. Jiang, J. Lv, T.X. Wu, Q.S. Yu, W.L. Zhu, J. Mol. Graph. Model. 28 (2010) 766. [60] D.C. Rapaport, Comput. Sci. Eng. 1 (2002) 70. [61] M. Gholampour, N. Edraki, M. Khoshneviszadeh, Iranian Pharmaceutical Science Congress and Symposium of Biopharmaceutics and Pharmacokinetics, 2015. [62] J.F.M.G W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, (2009). [63] O. Trott, A.J. Olson, J. Comput. Chem. 31 (2010) 455. [64] D.A. Case, R. Anandakrishnan, K. Ayers, R. Betz, Y.J. Bomble, R.A. Brown, S. Brozell, Q. Cai, D.S. Cerutti, T.E.C. Iii, (2012). [65] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1998) 926. [66] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23 (1977) 327. [67] H.J.C. Berendsen, J.P.M. Postma, W.F.V. Gunsteren, A. Dinola, J.R. Haak, J. Chem. Phys. 81 (1998) 3684. [68] U. Essmann, L. Perera, M. Berkowitz, T. Darden, H. Lee, L. Pedersen, S.H. Lee, (1995). [69] F. Weinhold, (2008). [70] J.S. Murray, Z. Peralta-Inga, P. Politzer, K. Ekanayake, P. Lebreton, Int. J. Quantum Chem. 83 (2001) 245. [71] K. Fukui, The Role of Frontier Orbitals in Chemical Reactions, 1982.