Molecular Dynamics Simulations and MM/PBSA Analysis of Annocatacin B in ND1 Subunit of Human Mitochondrial Respiratory Complex I

ND1 subunit possesses the majority of the inhibitor binding domain of the human 1 MRC-I. This is an attractive target for the search for new inhibitors that seek mitochondrial 2 dysfunction. It is known, from in vitro experiments, some metabolites from Annona muricata called 3 acetogenins have important biological activities such as anticancer, antiparasitic, and insecticide. 4 Previous studies propose an inhibitory activity of bovine MRC-I by bis-THF acetogenins such 5 as annocatacin B, however, there are few studies on its inhibitory effect on human MRC-I. In 6 this work, we evaluate the molecular and energetic affinity of the annocatacin B molecule with 7 the human ND1 subunit in order to elucidate its potential capacity to be a good inhibitor of this 8 subunit. For this purpose, QM optimizations, MD simulations and MM/PBSA analysis were 9 performed. As a control to compare our outcomes, the molecule rotenone, which is a known 10 MRC-I inhibitor, was chosen. Our results show that annocatacin B has a greater affinity for the 11 ND1 structure, its size and folding were probably the main characteristics that contributed to 12 stabilize the molecular complex. Furthermore, the MM/PBSA calculations showed a 35% stronger 13 BFE compared to the rotenone complex. Detailed analysis of the BFE shows that the aliphatic 14 chains of annocatacin B play a key role in molecular coupling by distributing favorable interactions 15 throughout the major part of the ND1 structure. These results are consistent with experimental 16 studies that mention that acetogenins may be good inhibitors of MRC-I. 17


Introduction
It has been almost 100 years since Warburg presented the first connection between 21 the mitochondria and tumors appearance [1]. The mitochondria fulfill an energetic role 22 in cells, specifically in cancer cells; this role is essential for developing tumors through 23 glycolysis [2,3]. On that basis, several mechanisms associated with tumor generation, 24 such as loss of enzymatic function, mitochondrial genome mutation, reprogramming of 25 mitochondrial metabolism, have been studied [4,5]. Some hypotheses and studies show 26 that, to a greater or lesser extent, neoplastic cells have many phenotypes related to their 27 energy production, from high aerobic glycolysis, through a partially active oxidative 28 phosphorylation, to a highly productive one [6,7]. Although this issue is controversial [8]. For instance, the mitochondrial respiratory complex I (MRC-I) is directly involved 30 in the appearance of colorectal cancer [9], prostate cancer [10], endometrial cancer [11], 31 breast cancer [12], and melanoma [13]. Thus, this complex protein has become a thera- 32 peutic target to develop anticancer drugs. Besides, the MRC-I catalyze the formation of 33 reactive oxygen species (ROS). 34 MRC-I, also named ubiquinone oxidoreductase , has a molecular mass of approxi- 35 mately 1 MDa; its structural conformation is composed by fourteen central subunits. ND1 36 subunit is one of those and has most of the inhibitor binding domain in the ubiquinone 37 oxidoreductase. The main inhibitor of the MRC-I is the rotenone molecule. [14]. It is 38 an isoflavone compound and has been found in many Fabaceae plants. Furthermore, 39 it was used as a pesticide and piscicide [15] due to its high toxicity [16,17]. Its effect 40 on cancer cell lines has been evaluated in vitro, showing inhibition of proliferation and 41 induction of apoptosis [18,19]. Nevertheless, its toxicity in cells complicates its use as an 42 anticancer drug, mainly because it is highly neurotoxic due to its lipophilic nature and 43 also the fact that it does not need an extra metabolism to be active or transporter to enter 44 neurons [14,20,21] Consequently, the challenge is to find new inhibitors that could be 45 less toxic than rotenone. 46 Murai et al. analyzed rotenone and a synthetic acetogenin as an inhibitor of the 47 bovine heart MRC-I [22]. They revealed that acetogenins are involved in the binding 48 domain of several inhibitors as rotenone does.In fact, acetogenins with two adjacent 49 tetrahydrofurans (THF) rings were reported to show higher antitumor activity and 50 toxicity than those that had only one THF [23], and have been found in the family of 51 Annonaceae, i.e., soursop (Annonamuricata) [24]. 52 In traditional medicine, soursop also has important uses, including anticonvulsant, 53 antiarthritic, antiparasitic, hepatoprotective, etc. Many of these beneficial attributes have 54 been ascribed to acetogenins [24,25]. One of the most studied properties in soursop is its 55 potential anticarcinogenic effect due, in a way, to its powerful cytotoxic features [25,26]. 56 It has been possible to isolate more than 100 acetogenins from different parts of the 57 Annonaceae plants [24,25,27]. The effect of acetogenins as inhibitors of the MRC-I has 58 been suggested and demonstrated for more than 20 years [28]. 59 Acetogenins have showed important behaviors when evaluating their potential 60 cytotoxic activity against cancer cells; some of these molecules already have proven 61 anti-cancer properties such as bullatacin, motrilin, assimin, trilobacin, annonacin, gi- 62 anttronenin, and squamocin. However, we still do not have enough information about 63 most of the acetogenins [29]. The main characteristic of acetogenins' molecular structure 64 is their linear 32 to 34 carbon chains containing oxygen-containing functional groups.

76
We analyzed two molecules as ligands to the ND1 complex, rotenone (PubChem 77 ID 6758) and annocatacin B (PubChem ID 10483312) (Figure 1a). The structures of both 78 molecules were built using GaussView v.6 software package [31], and optimized by 79 DFT calculations using Gaussian 16 software package [32] (Figure 1b). The optimization 80 process were performed using the CAM-B3LYP exchange-correlation functional [33], and 81 the TZVP basis set [34]. The vibrational frequencies were calculated to ensure that the 82 geometries were those of the minimum energy. In order to investigate the electrostatic 83 effect of the ligands on ND1 complex, atomic charges were calculated using the Hirshfeld 84 population analysis [35][36][37], and molecular Electrostatic Potential (ESP) surfaces were 85 used to visualize the polar and non polar regions of these ligands. To obtain the MD 86 parameters and topologies of the ligands, we used the TPPMKOP server [38], which uses 87 the parameters of the OPLS-AA force field to generates them [39,40]. These topologies 88 were reparametrized using the optimized structures and atomic charges obtained in 89 previous quantum calculations. 90 On the other hand, the phospholipid bilayer membrane was built with 512 dipalmit- of ND1 protein in the lipid membrane [41].

95
The three-dimensional crystallographic structure of the human MRC-I was con-  [48].

135
Based on the g_mmpbsa program [49]. Calculations of free energies and energy contri-136 butions by residue were carried out in order to localize the main residue interactions 137 and to assess the effect of each residue on the ND1 -ligand complexes. The last 200 ns 138 of the MD trajectories were analyzed at a 1 ns time interval to estimate the binding free 139 energy (∆G bind ) which was calculated by the following equation: where G complex is the total free energy of the ND1-ligand complexes; G ND1 and 141 G lig , are the free energies of isolated ND1 structure and rotenone or annocatacin B in 142 solvent. ∆E MM , represents the molecular mechanics energy contibutions; ∆G sol is the 143 free energy solvation requiered to tranfer a solute from vacuum into the solvent. The

144
T∆S term refers the entropic contribution and was not included in this calculation due 145 to the computational costs [49][50][51]. Therefore, individual E MM , and G sol terms were 146 calculated as follow:  [43]. The graphs were plotted using XMGrace software [53]. 2D

173
The human MRC-I belong to a highly organized supercomplex, named respirasome.

174
The complexes I, III, and IV arise more stability at that supercomplex and have the special 175 task of channeling electrons effectively through the electron transport chain [59]. Never-176 theless, Guo et al. proposed an even larger system called megacomplex that includes 177 complex II at the previous respirasome [42]. They suggested that a quinone/quinol 178 (oxidized/reduced forms of the same molecule) pool maximize the oxide-reduction reac-179 tions. Recent studies suggest there are around 100 Å between complex I and complex 180 III when actively translocating electrons, proposing with this that there is no need for a 181 mediating protein to help the electron channeling through these complexes [42,60]. respectively [61]. The ND1 subunit has a predominantly structural conformation of 188 alpha-helices that provides the hydrophobic environment expected of a membrane 189 protein and owns the quinone binding domain which is in its core (Table 1).  very hydrophobic, being annocatacin B more lipophilic than rotenone, due mostly to its 220 alkyl chain.

221
The pharmacokinetic properties of absorption, distribution, metabolism, excretion, 222 and toxicity (ADMET) are in Table 2. The absorption is similar in both compounds; 223 however, rotenone is not a P-glycoprotein substrate giving a slim advantage to the other 224 molecule. The distribution property is slightly higher for rotenone, which implies that its  As we said earlier, ND1 subunit is located in the transmembrane region of human 233 MRC-I [42]. In order to understand the ligand effect on its structure, we carried out  For H-bonds calculations, we considered determining those formed between ND1 subunit itself (intra); the ND1 subunit and solvent molecules (inter/solv); and the ND1 subunit and lipid bilayer membrane (inter/mem).
Values between parenthesis were calculated on global minimum energy structures obtained in the FEL analysis.
All values were obtained from the last 300 ns of the MD simulations.

318
Despite the high stability of the active site, a fluctuations analysis was performed to 319 understand the ligand effect in these region. For this purpose, in addition to RMSF cal-320 culations, we analyzed the B-factor, also called thermal factor or Debye-Walle factor [69] 321 and we mapped the values on the active site surfaces.      ns of the MD trajectories. In addition, an energy decomposition analysis per residue 372 was performed to highlight the main residues that contribute to the stability of the 373 complexes. As shown in Table 5, the binding free energy (BFE) of the two complexes 374 was energetically favorable, however, the interaction energy of the ND1 -annocatacin 375 B complex was more spontaneous (-333.18 ± 2.14 kJ/mol) than the rotenone complex     interactions between the protein-ligand complexes, Table 7 shows the residues with the 403 highest positive values. In the rotenone complex, these residues were located previous 404 to the active site, namely S115, W118, and S119, being the W118 residue with the highest contributions was the oxygen of the hydroxyl group located at the pyranoloid ring.

424
Residues with the lowest contribution to the BFE were located on the methoxy groups, 425 being S115, W118, and S119 residues with high positive values. The non-polar solvation 426 contribution was the predominant energy term in these residues. In the case of the annocatacin B complex, the aliphatic chains of this ligand involve 428 the majority of the interactions with its amino acid environment. This includes the 429 interactions with V113 and I116, that were the main contribution and no-contribution 430 residues to the binding energy, respectively (Figure 10c, Figure 10d)