Novel Inhibitors for a Novel Binding Site in ... - ACS Publications

Feb 23, 2016 - distinctly different from both Qo and Qi sites (hence designated as Non-. Q binding site, NQ). NQ site binding pocket extends up close ...
1 downloads 0 Views 4MB Size
Article pubs.acs.org/JPCB

Novel Inhibitors for a Novel Binding Site in Respiratory Complex III Muhammad A. Hagras and Alexei A. Stuchebrukhov* Department of Chemistry, University of California One Shields Avenue, Davis, California 95616 S Supporting Information *

ABSTRACT: A new binding site and potential novel inhibitors of the respiratory complex III are described. The site is located at the opposite side of the enzyme with respect to ubiquinol binding site (Qo site), and distinctly different from both Qo and Qi sites (hence designated as NonQ binding site, NQ). NQ site binding pocket extends up close to Phe90 residue, an internal switch (LH switch) that regulates electron transfer between heme bL and heme bH of the low potential redox chain. Docking studies and molecular dynamics simulations of different molecules to the NQ site revealed potential ligands which exhibit a novel inhibitory effect for bc1 complex by switching the LH switch to “off” conformation, thereby shutting down electron transfer in the low potential redox chain. Moreover, the novel inhibitors have lower binding affinity for both Qo and Qi sites, and hence do not interfere with binding of the natural ligands to those sites. The inhibitory activity of those novel ligands in bc1 complex is suggested to promote the production of reactive oxygen species (ROS) at the Qo site. Hence those ligands are potential candidates for designing new “mitocan” drugs.

1. INTRODUCTION Ideally, a smart anticancer drug should discriminate between cancer and normal cells. One of the recently discovered biochemical variations in cancer cells that distinguish them from the normal ones is the higher basal level of ROS, which make them more susceptible to ROS-induced apoptosis.1−4 But since cancer cells can adapt to such oxidative stress by upregulating antioxidant production,5 to make use of such a mechanism, the drug should induce rapid production and accumulation of ROS and trigger apoptosis in cancer cells before antioxidant up-regulation takes effect. The mitochondrial electron transport chain (METC) is one of the major sources of ROS in the cell,6,7 and respiratory complex III (also known as Ubiquinol:cytochrome c oxidoreductase complex or bc1 complex) is one of the two chief producers of ROS in the METC.8,9 Therefore, the bc1 complex is a natural target for any candidate anticancer drug whose function would be to increase ROS production and to bring it to an over-the-threshold level and thereby trigger apoptosis in cancer cells, while leaving the normal cells in under-thethreshold level of ROS.10 Moreover, the bc1 complex has a docking site for the water-soluble electron carrier cyctochrome c which plays a critical role in cell apoptosis.11 Thus, the bc1 complex has a dual role in this context: as a major site for ROS production which exerts a cytotoxic effect,12 and as a source of oxidative stress leading to cytochrome c-mediated cell apoptosis. Here we describe a novel binding site and its ligands that may inhibit electron transport and promote ROS production in bc1 complex. © XXXX American Chemical Society

Located in the inner-mitochondrial membrane, bc1 complex is the middle player in the electron transport proton pumping orchestra, which converts the free energy of redox reactions to electrochemical proton gradient.13,14 The mitochondrial bc1 complex, shown in Supporting Information, Figure S1, has an intertwined dimeric structure comprising 11 subunits in each monomer, but only three of them have catalytic function. The core subunits include the Rieske domain, which incorporates an iron−sulfur cluster [2Fe-2S]; trans-membrane cytochrome b domain, incorporating a low-potential heme group (heme bL) and high-potential heme group (heme bH); and cytochrome c1 domain, containing a heme c1 group and two separate binding sites, Qo (or QP) site where the hydrophobic electron carrier ubihydroquinol QH2 is oxidized, and Qi (or QN) site where the ubiquinone molecule Q is reduced.15−18 We previously discovered that residue Phe90 plays the role of an intramonomeric switch (designated as LH switch) which can regulate the rate of electron transfer between heme bL and heme bH.19 Phe90 can exist in two conformations, designated as OFFLH and ONLH. Being in the OFFLH conformation, the ET rate diminishes by about 3 orders of magnitude compared to the Phe90 ONLH conformation. As we describe below, there is a unique binding pocket that extends deep into the enzyme body reaching the region of Phe90 switch. We conducted a series of docking calculations and molecular dynamics Received: December 17, 2015 Revised: February 14, 2016

A

DOI: 10.1021/acs.jpcb.5b12347 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the different ligands were computed based on the MM-PBSA free-energy calculation method using the GROMACS g_mmpbsa tool34 for which different energy terms (electrostatic energy Eelec, van-der Waals energy Evdw, and nonelectrostatic solvation energy using the solvent-accessible surface area ESASA) were computed based on the force field parameters. The polar solvation energy (Epolar) was calculated using the Adaptive Possion-Boltzmann Solver (APBS) tool (www.poissonboltzmann.org).35 2.6. Electron Transfer Calculations. The tunneling calculations were done as described in refs 36 and 37. The initial guess molecular orbitals (MOs) were generated using the unrestricted Broken-Symmetry (BS) B3LYP method of the Gaussian 09 package,38 employing the ZINDO basis set for valence electrons and pseudopotentials for core electrons. The pruned systems were partitioned into donor, bridge, and acceptor(s) fragments of different charges and spin multiplicities, and consequently the tunneling electron was localized either on the donor or the acceptor site. The computed initialguess MOs were then utilized in a subsequent BS-ZINDO calculation to obtain the corresponding BS donor and acceptor diabatic ground states,

simulations at the new docking site, which lead to a set of potential ligands that fix the Phe90 switch to its OFF position, shutting down electron transfer between heme bL and heme bH. The inhibiting ligands should result in a major increase of ROS production by the enzyme, and thus can be utilized to design a smart mitocan drug.

2. METHODS 2.1. Protein−Membrane Setup. The orientation of different bc1 molecules with respect to the membrane was determined by using the Peripheral Proteins in Membranes (PPM) server (http://opm.phar.umich.edu/server.php)20 of the Orientations of Proteins in Membranes (OPM) database. The oriented proteins were inserted in the membrane using CHARMM-GUI Membrane Builder21 online tool using a membrane composition phosphatidylcholine (PC)/phosphatidylethanolamine (PE)/phosphatidylglycerol (PG) which substitutes for cardiolipin (CL) equal to 40:24:36.22 The protein− membrane system was placed in a water box with 0.15 M KCl neutralizing ions. The assembled system comprising 450 000 atoms was subjected to energy minimization using GROMACS23 program with CHARMM PARAM36 force field24,25 for 1500 steps with periodic boundary condition. 2.2. Binding Site Prediction. Cytochrome b of the native bc1 complex was evaluated by the AutoLigand26 program to search for binding pockets. A grid box of size 50 Å3 of spacing 1 Å and centered at Tyr95 hydroxyl oxygen atom was utilized. Grid box size was increased in 50 Å3 increments and EPV (energy per volume) vs volume was plotted to find the minimum volume describing the binding site. 2.3. Ligand Preparation. NCI diversity subset III containing 1597 molecules was retrieved from the ZINC database which is freely available from the Shoichet Laboratory at the University of California, San Francisco (http://zinc. docking.org/).27 The subset was processed further by AutoDockTools28 utility (prepareligand4.py) to add any missing polar hydrogens while preserving the original partial charges and to convert from Sybyl mol2 format to Autodock PDBQT format. 2.4. Docking and Virtual Screening. The virtual screening calculations were carried out using PyRx0.9.229 utilizing both Lamarckian Genetic Algorithm (LG) search method in Autodock 4.2 and Autodock Vina 4.2 advanced gradient optimization search method which was proven to be more accurate than Autodock.30 Default settings were used for Autodock Vina. For the Autodock LG search method, we chose the default settings with docking parameters of 100 for ga_run, 10 000 000 for ga_num_evals, and 1.5 for tstep. 2.5. Molecular Dynamics Simulation. Each of the five discovered ligands (excluding the isomer ZINC17147426) is inserted into the novel NQ site of the previously energyminimized bc1 structure (PDB ID: 1BE3) as discussed above. Ligand insertion was accomplished by structural superposition of the docked cytochrome b chain (with the ligand) with energy-minimized bc1 structure using the Chimera program31 and the default settings for the Needleman−Wunsch algorithm with BLOSUM-62. Subsequently, the whole system was energy minimized using the CHARMM PARAM3632 force field for 300 000 steps with periodic boundary condition and then equilibrated for 100 ps. Molecular dynamics simulation was run on the whole system for 50 ns. Parameterization of the different discovered ligands were calculated using SwissParam tool (www.swissparam.ch).33 Intermolecular interaction energies of

|ΨD⟩ = |φ1,Dσ , ..., φND, σ ⟩,

|ΨA⟩ = |φ1,Aσ , ..., φNA, σ ⟩

(1)

The Hartree−Fock approximation assumed above was validated in our recent study.37 To simplify the tunneling calculation, the reduction from a fully multielectronic picture to a one-electron approximation was accomplished by using the biorthogonalization scheme of corresponding orbitals, as described in refs 39 and 40. The donor and acceptor orbitals with the smallest overlap, and hence with the most significant change during the tunneling transition, represent the tunneling pair of the one-electron approximation; the remaining corresponding (core) orbitals experience only relatively small change in the tunneling transition, with almost unit overlap, and contribute to the tunneling matrix element via the electronic Franck−Condon factor.37 The transformed corresponding orbitals have the following form, |Ψ′D⟩ = |ξ1,Dα , ..., ξlD, α , ξ1,Dβ , ..., ξtD, β⟩, |Ψ′A⟩ = |ξ1,Aα , ..., ξl A, α , ξ1,Aβ , ..., ξtA, β⟩

(2)

⟨Ψ′D|Ψ′A⟩ = ⟨ξiD, σ |ξjA, σ ⟩ = δijsiσ

(3)

The product of the overlaps sσi of the core orbitals forms the electronic Franck−Condon factor ∏lisαi ∏tj ≠ tsβj , for which we assume that the tunneling electron has a β-spin, and the tunneling orbital index “t” is the last one of the β-spin orbitals. 2.7. Tunneling Current Calculation. Since the ZINDO canonical MOs are represented in a basis set which is orthogonalized (using Löwdin-orthogonalization41 or other schemes) and therefore delocalized over many atoms, the localized atomic picture of interatomic tunneling currents in a tunneling transition is lost in this basis set. To conduct the tunneling current calculations, the atomic basis set localization needs therefore to be restored. By introducing the localized atomic basis set and the Mulliken type coarse-graining of the tunneling current, the tunneling transition flux is expressed in terms of the interatomic currents, which (approximately37) has the following form:39,42 B

DOI: 10.1021/acs.jpcb.5b12347 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

∏ Siα ∏ Sjβ ∑ ∑ (Hvμ − E0Svμ)(θμv − θvμ)

Jab =

i

j≠t

heme bH to slow it down, or shut it off completely. Such interference would lead to a greater localization of the second electron (transferred from QH2 at Qo site) on heme bL and hence to a higher chance to be picked up by oxygen molecules leading to ROS production. Motivated by the discovery of ET switches, we identified a new binding pocket (designated as NonQ-site or NQ site for short) in the bc1 complex. As shown in Figure 1, the NQ site is

(4)

v∈a μ∈b

Here ν and μ are the atomic orbitals of atoms a and b; θνμ = AνDμ where Dμ and Aν are the expansion coefficients of the donor and acceptor tunneling orbitals, respectively; Hνμ and Sνμ are core Hamiltonian and overlap matrix, and E0 is a tunneling orbital energy defined by E0 =

∑ DλFλρDρ = ∑ AλFλρAρ λ,ρ

(5)

λ,ρ 39

where Fλρ is the reduced Fock matrix. The second equality in eq 7 corresponds to the resonance of the donor and acceptor energies at the transition state of electron transfer reaction; in practice, the resonance is achieved by applying a static electric field mimicking the action of the polar environment and solvation effects as described in refs 36 and 37. Equation 6 shows that the interatomic currents are primarily determined by the overlap of the tails of the two tunneling orbitals. The electronic Franck−Condon factor (∏isαi ∏j≠tsβj ) contributes as a uniform scaling factor for all interatomic currents. The total atomic current through a given atom a, Jatot ≡

1 2

∑ |Ja ,b | b

(6)

is proportional to the probability that the tunneling electron is passing through this atom; as such, it provides a convenient way (e.g., with color gradient coding, as in Figure S2) of identifying atoms that constitute the tunneling pathways. The tunneling matrix element that determines the rate of electron transfer (ET) is calculated as the total flux across the dividing surface between the donor and acceptor redox complexes (the flux theorem43), TDA = −ℏ ∑

∑ Jab

a∈S b∉S

Figure 1. bc1 homodimer complex (PDB 1NTZ) is shown embedded in the membrane where the chains of cytochrome b domain are either removed or transparent for clarity. Relative positions are shown of heme bL (light green color) and heme bH (dark green color) with the intervening internal switch Phe90 residue (yellow color). Different binding sites are visualized including the Qo binding pocket as computed with AutoLigand (dark blue color) occupied by a QH2 molecule (vdW spheres in magenta color), Q molecule occupying Qi site (vdW spheres in salmon color), and NQ site as computed by AutoLigand (vdW spheres in light blue color). Water molecules are displayed as red spheres. In the inset, zoomed-in visualization of the NQ site entrance is displayed.

(7)

Here, the summation of interatomic tunneling currents between atoms a on one side, denoted as S side, of the dividing surface and atoms b on the other side.44 The interatomic tunneling currents provide an internal assessment of the quality of calculation as a measure of conservation of the total tunneling flux between the donor and acceptor redox centers. An example is shown in Figure S2 for a model heme bL → heme bH system with Phe90 in the ON conformation. The flux is approximately conserved along the whole tunneling pathway especially in the middle part, between donor and acceptor, as it should;40 to account for small variations, we compute the tunneling matrix element as an average over the middle region. All the studied systems behave in a similar manner. The value of the tunneling flux was taken as a measure of electron transfer rate (see details in refs 36 and 37) in the analysis of different ligand-bound variants of the bc1 complex.

located at the opposite side of the enzyme with respect to the Qo site; like Qo, it is buried inside the inner- mitochondrial membrane, and its entrance is fully hydrated. In contrast to the Qo site, however, the NQ site penetrates deeply in the cytochrome b domain and reaches very closely the LH region. Hence the NQ site provides a suitable binding pocket for ligands that can influence the orientation of the Phe90 residue, and hence modulate the corresponding ET rate between heme bL and heme bH. We conducted a series of virtual screening computations using both Vina and LG search algorithms by docking 1597 ligands of the NCI Diversity Set III into the NQ site. To account for flexibility of the NQ-binding site and the different possible conformations of the Phe90 residue, we performed the docking calculation using both fixed and flexible residues (Phe90, Phe91, Ile92, Leu94, and Tyr95) as shown in Figure 2. Only 39 ligands show a preferential binding by at least 1 kcal/mol at the NQ site in the flexible-residues docking mode compared to the fixed mode in both docking search algorithms. In this way it was ensured that those picked ligands adjust the

3. DISCUSSION AND RESULTS We previously discovered two internal switches in the bc1 complex.19 One of them is Phe90 residue, which upon binding of aromatic ligands at the Qo site switches its own conformation from a tilted one (OFFLH) to a parallel one (ONLH) facilitating an electron transfer reaction between heme bL and heme bH. We built our drug design idea on this finding reasoning that by modifying the conformation of the LH switch, we can interfere with the trans-membrane ET reaction between heme bL and C

DOI: 10.1021/acs.jpcb.5b12347 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 2. Docking results of the NCI Diversity Set III at the NQ site using fixed residues (shown in black color) against flexible residues (shown in red color). The positive difference of the binding energies (fixed−flexible) is plotted at the top (shown in blue color). The white line represents the cutoff binding energy of the QH2 molecule at the NQ site.

NQ site residues and hence bind more firmly. Such modification of the NQ site residues is desirable because it affects the Phe90 internal switch orientation, which is our ultimate target. Subsequently we investigated the influence of those picked ligands on the orientation of Phe90 and chose those which orient Phe90 residue away from the intraheme b region. In addition, we performed docking calculation of the QH2 at the NQ site as above (indicated by white line in Figure 2) to exclude the ligands that have flexible-mode binding energy higher than that of the QH2 molecule at the NQ site (which is on average around −8.75 kcal/mol) in either Vina or LG search methods. We ended up with 30 ligands which have both preferential binding at the NQ site compared to QH2 molecule and induce the desired conformational change of Phe90 residue. In addition to the above criteria, we further restricted the ligand search such that the candidate ligands would not compete with the QH2 molecule for binding at the Qo site. The idea here is to minimize any possible inhibition of the enzyme by those ligands so as to retain heme bL reduction by QH2, while at the same time shutting down heme bL to heme bH ET and hence allowing for ROS production. Further docking analysis of the resultant 30 ligands into the active Qo site (as it exists in PDB 1NTZ) reduced the number of ligands down to six candidates. As shown in Figure 3 and listed in Table 1, those six candidate ligands exhibit binding energy lower than 8.8 kcal/mol and hence would bind preferentially at NQ site, but not at Qo site. The chemical structures of those ligands are shown in Figure 4. As displayed in Figure 5, the binding of those six candidate ligands at the NQ site induces a series of conformational changes. Tyr95 undergoes moderate-to-substantial conformational changes that induce (along with a strong aromatic interaction due to the docked ligand) a conformational change of Phe91 residue, which ultimately attracts the internal switch Phe90 residue from the intraheme b region to OFF position. In short, compared to the natural ONLH conformation (which is attained by binding of QH2 molecule at Qo site shown in red in Figure 5), upon binding of 6 candidate ligands at the NQ site,

Figure 3. Docking results at Qo site (1NTZ in PDB) of the selected 30 ligands using Vina (blue bars) and LG (light red bars) search algorithms. Magenta bars represent the overlap between Vina and LG corresponding bars.

Phe90 residue swings away by a torsion angle of 76.31° and clusters at a conformation which we designate as OFFNQ (torsion angle equals 9.77° compared to ONLH with a torsion angle of 86.08°), as indicated in Figure 5. As shown in Figure 6, in the OFFNQ conformation the electron tunneling flux is dramatically diminished, compared with the native ONLH conformation, which is an indication of the significantly reduced electron transfer rate between heme bL and heme bH. Subsequent to the above docking simulations, we conducted molecular dynamics (MD) simulations of the five discovered ligands (excluding the isomer ZINC17147426) bound to the ligand-free bc1 structure (PDB ID: 1BE3) in order to verify the binding affinities of those ligands to the NQ site and in addition to confirm the proposed induced OFFNQ conformational change of the internal switch Phe90 residue. As shown in Table 2 the binding energies of the different ligands vary significantly in contrast to the results obtained above in the docking simulations (Table 1). Ligands ZINC01691943 and D

DOI: 10.1021/acs.jpcb.5b12347 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 1. Vina and Lamarkian Genetic Docking Results for the Six Candidate Ligands and QH2 Molecule at NQ and Qo Sites binding energy at NQ site (kcal/mol) Vina Method

LG Method

binding energy at Qo site (kcal/mol)

ligand

fix

flex

fix

flex

Vina Method

LG Method

ZINC01661542 ZINC01691943 ZINC01748097 ZINC17465983 ZINC17147424 ZINC17147426 QH2

−8.6 −7.6 −8.7 −8.3 −7.5 −7.8 −7.9

−10.4 −8.7 −10.1 −9.5 −8.6 −8.8 −8.5

−7.2 −8.3 −8.3 −7.5 −8.6 −8.2 −8.2

−8.7 −9.4 −9.8 −12.4 −11.3 −10.9 −9.0

−8.0 −8.8 −8.3 −7.8 −8.1 −8.4 −8.5

−7.2 −8.3 −8.3 −7.5 −8.6 −8.2 −8.3

Figure 4. Distinct chemical structures of the candidate 5 ligands (ZINC17147424 and ZINC17147426 are isomers).

Figure 6. Calculated electron tunneling flux densities between heme bL and heme bH in two X-ray bc1 crystal structures (PDB 1NTZ in A and C and 1NTK in B). Intensity of red color corresponds to probability of a tunneling electron being on a given atom: (A) residue Phe90 exists in ON LH conformation; (B) residue Phe90 exists in OFF LH conformation; (C) residue Phe90 exits in OFFNQ conformation.

Figure 5. Docking of the selected ligands at the NQ site is displayed. Different conformations of the flexible residues including the internal switch Phe90 residue are shown. The arrows signify the change from native ONLH conformation to the new ones.

Analysis of the MD simulations trajectories for the different ligands further validates their binding mode and proposed influence on Phe90 conformation. As shown in Figure S3A, ligands ZINC01691943 and ZINC01748097 remain bound to the NQ site during the whole 50 ns simulation time and even get closer to the heme bL−heme bH system. Those two ligands show a high magnitude of stability as indicated by their minimal

ZINC01748097 seem to be the most promising ones showing the lowest stable binding energies while ligands ZINC01661542 and ZINC17465983 show comparably higher binding energies and finally ligand ZINC17147424 has the highest binding energy of all. E

DOI: 10.1021/acs.jpcb.5b12347 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 2. Binding Energy Averaged Components for the Five Ligands at the NQ Site during the MD Simulation Trajectory (in units of kcal/mol)a Eelec ZINC01661542 ZINC01691943 ZINC01748097 ZINC17147424 ZINC17465983 a

−36.47 −30.52 −35.57 −21.22 −29.44

(1.36) (0.56) (1.22) (1.66) (4.77)

Evdw

Epolar

0.65 (0.25) −56.19 (4.96) 2.44 (3.69) −8.70 (3.70) −3.38 (1.94)

10.36 (0.88) 17.91 (1.86) −9.16 (0.51) 14.22 (2.37) 13.34 (0.69)

ESASA −3.99 −3.38 −4.49 −2.81 −3.57

(0.11) (0.26) (0.23) (0.18) (0.16)

Ebind −29.45 −72.17 −46.77 −18.50 −23.05

(0.83) (6.24) (4.08) (3.82) (6.17)

The values in parentheses indicate the variation estimate of the corresponding quantity along the MD trajectory.

candidate ligands that satisfy all the criteria. Those ligands induce a cascade of conformational changes in the NQ-binding site that eventually, through aromatic−aromatic interactions, force Phe90 residue to OFFNQ conformation. Furthermore, compared to the binding of the QH2 molecule, those ligands bind more strongly at the NQ site but less strongly at the Qo site, thus providing for preferential binding at the NQ site, without interfering QH2 binding at Qo. As shown in Figure 4, the evident aromatic character and the bulkiness of the unique five ligands (two or more aromatic rings) are important factors to achieve the desired conformational changes of the Phe90 residue, to push it to OFFNQ conformation. Interestingly, docking the QH2 molecule in the flexible-mode at the NQ site (shown in Figure 5 in orange color) shows negligible conformational change of different residues at the NQ binding site, confirming the importance of aromatic character of the selected ligands. Out of the five candidate ligands, ligand ZINC01691943 is the most promising tightly bound one at the NQ site and also induces the desired conformational changes of the Phe90 residue. In addition, even suffering higher binding energy, ligand ZINC17465983 remains bound at the NQ site during the whole 50 ns simulation time and could be further modified to decrease the binding energy. It is interesting to notice that both ligands ZINC01691943 and ZINC17465983 share a common structure of a pair of two aromatic fused rings connected by a flexible rotatable region. It seems that this “motif” is required for binding at the newly discovered NQ site and in addition to induce the correct conformational changes of Phe90 through strong aromatic−aromatic interaction offered by the fused aromatic rings. Moving up in the number of the fused rings deteriorates the binding at the NQ site as indicated by the undocking of both ligands ZINC01661542 (three fused rings) and ZINC17147424 (four fused rings) which indicates that two fused rings is the most optimal for binding at the NQ site. On the other hand, moving down in the number of rings still leads to a tight binding to the NQ site as shown for ligand ZINC01748097 (which has only a pair of aromatic rings separated by a flexible region). Despite having a low stable binding energy, such a number of aromatic rings does not confer a strong aromatic−aromatic interaction to induce the desired conformational changes of Phe90 residue. In summary, a pair of two fused aromatic rings seems to be perfect for both binding at the NQ site and to induce the correct Phe90 conformational changes. As listed in Table 3 and displayed in Figure S2, for OFFNQ conformation, the donor/acceptor tunneling orbitals overlap for the heme bL → heme bH redox system is the lowest. Accordingly, the tunneling matrix element is lower by about 1.5 orders of magnitude, and the ET rate by about 3 orders of magnitude, in the OFFNQ conformation compared to that in the ONLH conformation. Such a dramatic reduction in ET rate

RMSD during the simulation time (Figure S3B). Nevertheless, only ligand ZINC01691943 induces the desired conformational changes of Phe90 as indicated in Figure S3C where the total computed distance of the displaced Phe90 is greater than that found in the Q-bound bc1 structure (PDB: 1NTZ, distance = 13.3 Å) along the whole simulation time. Out of the next two ligands ZINC01661542 and ZINC17465983, only ligand ZINC17465983 remains tightly bound and produces the desired conformational changes of the Phe90 residue. Finally, as expected for having the highest binding energy, ligand ZINC17147424 undocks from the NQ site, showing the highest RMSD at the NQ site which decreases as it leaves the site and also fails to retain the desired conformational changes of the Phe90 residue. Ultimately, ligand ZINC01691943 is the most promising one among the five discovered ligands while ligand ZINC17147424 is the least promising one. Visualization of MD snapshots of the two ligands along with the Phe90 residue (and compared to the other monomer Phe90 with a vacant Qo/Qi/NQ sites) would strongly assert our conclusions. As shown in Figure S4, ligand ZINC01691943 remains docked at the NQ site over all the simulation time and even binds tighter over time (as indicated by blue-to-red color gradient). In addition, it displaces the Phe90 residue away from the intraheme b region as compared to the other monomer Phe90 residue. On the other hand, ligand ZINC17147424 unbinds from the NQ site and leads to undesired conformational changes of Phe90 residue. In summary, we have identified two ligands (ZINC01691943 and ZINC17465983) out of the five discovered candidate ligands, that preferentially bind at the novel NQ site, induce the correct desirable conformational changes of the internal switch Phe90 residue, and hence significantly diminish ET between heme bL and heme bH of the low potential redox chain, while leaving the binding of QH2 at the Qo site intact, therefore allowing for the reduction of heme bL and possible production of ROS.

4. SUMMARY AND CONCLUSIONS An innovative inhibitory mechanism for bc1 enzyme is proposed in which the inhibitors do not compete with the natural ligands, but instead shut down ET between heme bL and heme bH of the low potential redox chain by modifying the conformation of the Phe90 residue. Because of the noncompetitive nature of those inhibitors at the Qo site, the normal function of bc1 is retained but with a slower ET rate between heme bL and heme bH. Such a slower ET rate is attained due to the modulating nature of the Phe90 residue. We have identified a novel NQ site that extends deep into the enzyme and reaches up to Phe90 residue, and hence provides the means to control its conformation. In various docking experiments that we conducted, we filtered the 1597 ligands in NCI Diversity Set III down to only five unique F

DOI: 10.1021/acs.jpcb.5b12347 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(7) Lenaz, G. The Mitochondrial Production of Reactive Oxygen Species: Mechanisms and Implications in Human Pathology. IUBMB Life 2001, 52, 159−164. (8) Turrens, J. F.; Boveris, A. Generation of Superoxide Anion by the Nadh Dehydrogenase of Bovine Heart Mitochondria. Biochem. J. 1980, 191, 421−427. (9) Sugioka, K.; Nakano, M.; Totsune-Nakano, H.; Minakami, H.; Tero-Kubota, S.; Ikegami, Y. Mechanism of O2- Generation in Reduction and Oxidation Cycle of Ubiquinones in a Model of Mitochondrial Electron Transport Systems. Biochim. Biophys. Acta, Bioenerg. 1988, 936, 377−85. (10) Trachootham, D.; Alexandre, J.; Huang, P. Targeting Cancer Cells by Ros-Mediated Mechanisms: A Radical Therapeutic Approach? Nat. Rev. Drug Discovery 2009, 8, 579−591. (11) Jiang, X.; Wang, X. Cytochrome C-Mediated Apoptosis. Annu. Rev. Biochem. 2004, 73, 87−106. (12) Fruehauf, J. P.; Meyskens, F. L., Jr. Reactive Oxygen Species: A Breath of Life or Death? Clin. Cancer Res. 2007, 13, 789−94. (13) Mitchell, P. Coupling of Phosphorylation to Electron and Hydrogen Transfer by a Chemi-Osmotic Type of Mechanism. Nature 1961, 191, 144−148. (14) Mitchell, P. Chemiosmotic Coupling in Oxidative and Photosynthetic Phosphorylation. Biol. Rev. Camb. Philos. Soc. 1966, 41, 445−502. (15) Xia, D.; Yu, C. A.; Kim, H.; Xia, J. Z.; Kachurin, A. M.; Zhang, L.; Yu, L.; Deisenhofer, J. Crystal Structure of the Cytochrome Bc1 Complex from Bovine Heart Mitochondria. Science 1997, 277, 60−6. (16) Schagger, H.; Link, T. A.; Engel, W. D.; von Jagow, G. Isolation of the Eleven Protein Subunits of the Bc1 Complex from Beef Heart. Methods Enzymol. 1986, 126, 224−37. (17) Zhang, Z.; Huang, L.; Shulmeister, V. M.; Chi, Y.-I.; Kim, K. K.; Hung, L.-W.; Crofts, A. R.; Berry, E. A.; Kim, S.-H. Electron Transfer by Domain Movement in Cytochrome Bc1. Nature 1998, 392, 677− 684. (18) Iwata, S.; Lee, J. W.; Okada, K.; Lee, J. K.; Iwata, M.; Rasmussen, B.; Link, T. A.; Ramaswamy, S.; Jap, B. K. Complete Structure of the 11-Subunit Bovine Mitochondrial Cytochrome Bc1 Complex. Science 1998, 281, 64−71. (19) Hagras, M. A.; Hayashi, T.; Stuchebrukhov, A. A. Quantum Calculations of Electron Tunneling in Respiratory Complex Iii. J. Phys. Chem. B 2015, 119, 14637−14651. (20) Lomize, M. A.; Pogozheva, I. D.; Joo, H.; Mosberg, H. I.; Lomize, A. L. Opm Database and Ppm Web Server: Resources for Positioning of Proteins in Membranes. Nucleic Acids Res. 2012, 40, D370−6. (21) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. Charmm-Gui: A Web-Based Graphical User Interface for Charmm. J. Comput. Chem. 2008, 29, 1859−1865. (22) Zinser, E.; Sperka-Gottlieb, C. D.; Fasch, E. V.; Kohlwein, S. D.; Paltauf, F.; Daum, G. Phospholipid Synthesis and Lipid Composition of Subcellular Membranes in the Unicellular Eukaryote Saccharomyces Cerevisiae. J. Bacteriol. 1991, 173, 2026−34. (23) Lindahl, E.; Hess, B.; van der Spoel, D. Gromacs 3.0: A Package for Molecular Simulation and Trajectory Analysis. J. Mol. Model. 2001, 7, 306−317. (24) Brooks, B. R.; Brooks, C. L., 3rd; Mackerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; Caflisch, A.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. Z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. Charmm: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−614. (25) Best, R. B.; Zhu, X.; Shim, J.; Lopes, P. E.; Mittal, J.; Feig, M.; Mackerell, A. D., Jr. Optimization of the Additive Charmm All-Atom Protein Force Field Targeting Improved Sampling of the Backbone Phi, Psi and Side-Chain Chi(1) and Chi(2) Dihedral Angles. J. Chem. Theory Comput. 2012, 8, 3257−3273.

Table 3. Tunneling Calculation Results for the New Phe90 Conformation against the Native ON/OFF Conformations Phe90 conformation ONLH OFFLH OFFNQ

TMO overlap −6

8.05 × 10 7.34 × 10−6 2.56 × 10−6

⟨TDA⟩ (cm−1) −1

2.68 × 10 7.86 × 10−3 7.50 × 10−3

ET rate (s−1) 2.41 × 10+5 2.06 × 10+2 1.88 × 10+2

will localize the electron at heme bL for a longer time and therefore will eventually promote an excessive ROS production, which can potentially trigger the apoptotic pathway in cancer cells.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b12347. Figure S1 shows the structure of the bc1 complex embedded in a membrane inside a water box (and showing the different redox centers and the Qo/Qi binding sites). Figure S2 shows the electron tunneling flux through a dividing surface in the heme bL → heme bH redox system at three different Phe90 conformations (ONLH, OFFLH, and OFFNQ). Figure S3 displays molecular dynamics simulation results of the five discovered ligands docked in pre-energy minimized bc1 complex. Figure S4 shows the visualization of 1000 snapshots o f t wo di ffe r e n t d o c k e d l i g a nd s (ZINC01691943 and ZINC17147424) at NQ-site in bc1 structure along with snapshots of the Phe90 residue (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors appreciate the useful comments of our colleague Benjamin Samudio, and his help with the docking calculations. Also the authors appreciate the guidance of Dr. Igor Leontyev in conducting the molecular dynamics simulations. This work has been supported by NIH Grant GM054052.



REFERENCES

(1) Szatrowski, T. P.; Nathan, C. F. Production of Large Amounts of Hydrogen Peroxide by Human Tumor Cells. Cancer Res. 1991, 51, 794−8. (2) Kawanishi, S.; Hiraku, Y.; Pinlaor, S.; Ma, N. Oxidative and Nitrative DNA Damage in Animals and Patients with Inflammatory Diseases in Relation to Inflammation-Related Carcinogenesis. Biol. Chem. 2006, 387, 365−72. (3) Toyokuni, S.; Okamoto, K.; Yodoi, J.; Hiai, H. Persistent Oxidative Stress in Cancer. FEBS Lett. 1995, 358, 1−3. (4) Trachootham, D.; Alexandre, J.; Huang, P. Targeting Cancer Cells by Ros-Mediated Mechanisms: A Radical Therapeutic Approach? Nat. Rev. Drug Discovery 2009, 8, 579−91. (5) Tiligada, E. Chemotherapy: Induction of Stress Responses. Endocr.-Relat. Cancer 2006, 13, S115−S124. (6) Adam-Vizi, V.; Chinopoulos, C. Bioenergetics and the Formation of Mitochondrial Reactive Oxygen Species. Trends Pharmacol. Sci. 2006, 27, 639−645. G

DOI: 10.1021/acs.jpcb.5b12347 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (26) Harris, R.; Olson, A. J.; Goodsell, D. S. Automated Prediction of Ligand-Binding Sites in Proteins. Proteins: Struct., Funct., Genet. 2008, 70, 1506−17. (27) Irwin, J. J.; Shoichet, B. K. Zinc-a Free Database of Commercially Available Compounds for Virtual Screening. J. Chem. Inf. Model. 2005, 45, 177−182. (28) 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, 2785−2791. (29) Dallakyan, S.; Olson, A. J. Small-Molecule Library Screening by Docking with Pyrx. Methods Mol. Biol. (N. Y., NY, U. S.) 2015, 1263, 243−50. (30) Trott, O.; Olson, A. J. Autodock Vina: Improving the Speed and Accuracy of Docking with a New Scoring Function, Efficient Optimization and Multithreading. J. Comput. Chem. 2010, 31, 455− 461. (31) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. Ucsf Chimera–a Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605−12. (32) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D., Jr.; Pastor, R. W. Update of the Charmm All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. J. Phys. Chem. B 2010, 114, 7830−43. (33) Zoete, V.; Cuendet, M. A.; Grosdidier, A.; Michielin, O. Swissparam: A Fast Force Field Generation Tool for Small Organic Molecules. J. Comput. Chem. 2011, 32, 2359−2368. (34) Kumari, R.; Kumar, R.; Lynn, A. G_Mmpbsa· a Gromacs Tool for High-Throughput Mm-Pbsa Calculations. J. Chem. Inf. Model. 2014, 54, 1951−1962. (35) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of Nanosystems: Application to Microtubules and the Ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (36) Hayashi, T.; Stuchebrukhov, A. A. Electron Tunneling in Respiratory Complex I. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 19157−19162. (37) Hagras, M. A.; Stuchebrukhov, A. A. Transition Flux Formula for the Electronic Coupling Matrix Element. J. Phys. Chem. B 2015, 119, 7712−7721. (38) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Gaussian, Inc.: Wallingford, CT, USA, 2009. (39) Stuchebrukhov, A. A. Tunneling Currents in Long-Distance Electron Transfer Reactions. V. Effective One Electron Approximation. J. Chem. Phys. 2003, 118, 7898−7906. (40) Stuchebrukhov, A. A. Long-Distance Electron Tunneling in Proteins. Theor. Chem. Acc. 2003, 110, 291−306. (41) Löwdin, P.-O. A Quantum Mechanical Calculation of the Cohesive Energy, the Interionic Distance, and the Elastic Constants of Some Ionic Crystals. Ark. Mater. Astr. Fys. A 1947, 35, 9. (42) Stuchebrukhov, A. A. Tunneling Currents in Long-Distance Electron Transfer Reactions. Iv. Many-Electron Formulation. Nonorthogonal Atomic Basis Sets and Mulliken Population Analysis. J. Chem. Phys. 1998, 108, 8510−8520.

(43) Stuchebrukhov, A. A. Toward Ab Initio Theory of LongDistance Electron Tunneling in Proteins: Tunneling Currents Approach. Adv. Chem. Phys. 2001, 118, 1−44. (44) Stuchebrukhov, A. A. Tunneling Currents in Electron Transfer Reactions in Proteins. J. Chem. Phys. 1996, 104, 8424−8432.

H

DOI: 10.1021/acs.jpcb.5b12347 J. Phys. Chem. B XXXX, XXX, XXX−XXX