Computational MitoTarget-Scanning Based on Topological Vacancies

Subscriber access provided by ECU Libraries

Article

Computational MitoTarget-Scanning Based on Topological Vacancies of Single-Walled Carbon Nanotubes with Human Mitochondrial hVDAC1 Channel. Michael Gonzalez Durruthy, José M Monserrat, Patricia Viera de Oliveira, Solange Binotto Fagan, Adriano V. Werhli, Karina Machado, André Melo, Humberto González-Díaz, Riccardo Concu, and Maria Natália D.S. Dias Soeiro Cordeiro Chem. Res. Toxicol., Just Accepted Manuscript • Publication Date (Web): 14 Mar 2019 Downloaded from http://pubs.acs.org on March 14, 2019

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

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

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

Chemical Research in Toxicology

Computational MitoTarget-Scanning Based on Topological Vacancies of Single-Walled Carbon Nanotubes with Human Mitochondrial hVDAC1 Channel. Michael González-Durruthya, *, José M. Monserratb,c, Patricia Viera de Oliveirad, Solange Binotto Fagand, Adriano V. Werhlie, Karina Machadoe, André Meloa, Humberto González-Díazf, g, Riccardo Concua,*, M. Natália D.S. Cordeiroa,* aDepartment bInstitute

of Chemistry and Biochemistry, Faculty of Science, University of Porto, Porto, Portugal

of Biological Sciences (ICB), Universidade Federal do Rio Grande -FURG, 96270-900, Rio

Grande, RS, Brazil. cICB-FURG

Post-Graduate Program in Physiological Sciences, 96270-900, Rio Grande, RS, Brazil

dPost-Graduate

Program in Nanoscience-Franciscana University-UFN, RS, Brazil.

eCenter

of Computational Sciences (C3)- Federal University of Rio Grande - FURG, Cx. P. 474, CEP 96200-970, Rio Grande, RS, Brazil fDepartment of Organic Chemistry II, University of the Basque Country UPV/EHU, 48940, Leioa, Biscay, Spain. gIKERBASQUE, Basque Foundation for Science, 48011, Bilbao, Biscay, Spain. * To whom correspondence should be addressed: M. González-Durruthy, Email: [email protected], Fax: +351220402659 R. Concu, Email: [email protected], Fax: +351220402659 M. Natália D.S. Cordeiro, Email: [email protected], Fax: +351220402659

1

ACS Paragon Plus Environment

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

Page 2 of 29

'For TOC only'

2


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


Abstract We present an in-silico approach for modeling the non-covalent interactions between human mitochondrial voltage-dependent anion channel (hVDAC1) with a family of single walled carbon nanotubes (SWCNT) with a defined pattern of topological vacancies (from v = 1 to 16), obtained by removing atoms from the SWCNT surface. The general results showed more stable docking interaction complexes (SWCNThVDAC1), with more negative Gibbs free energy of binding-affinity values, and a strong dependence with the vacancy number (R2: 0.93) and vacancy formation energy (R2: 0.96). In addition, for most of the SWCNT-vacancy analyzed, the interatomic distances for the SWCNT-hVDAC1 interactions with the functional catalytic residues (i.e.: Pro7, Gln199, Gln182, Phe181, Val20, Asp19, Lys 15, Gly14, Asp12, Ala11, and Arg18) that form the hVDAC1 active site (i.e. the voltage-sensing N-terminal -helix segment) were very similar or lower than the interatomic distances of these residues for ATP-hVDAC1 interactions. Particularly, the hVDAC1 residues that can be phosphorylated like Tyr10, Tyr198 and Se16 were significantly perturbed by the interactions with SWCNT with a # of vacancies v ≥ 9. Besides, the SWCNTvacancy family members can affect the flexibility properties of hVDAC1-N-terminal-α-helix segment inducing different patterns of local-perturbations in the inter-residue communication. Finally, VacancyQuantitative-Structure Binding-Relationships (V-QSBR) were unveiled for setting up a robust model able to predict the strength of docking interactions between SWCNT with a specific topological vacancy and the hVDAC1. The developed V-QSBR model classified properly all the SWCNT with a different # of SWCNT-vacancies with exceptional sensitivity and specificity (both equal to 100%), indicating a strong potential to unequivocally predict the influence of SWCNT-vacancies in the mitochondrial channel interactions.

Keywords: nanotoxicology, single-walled carbon nanotubes, hVDAC1, molecular docking, QSAR models.

3



Page 4 of 29

Introduction The human mitochondrial voltage-dependent anion channel (hVDAC1) is the most abundant protein in the outer mitochondrial membrane of all cells, which in turn is the main communication route between this mitochondrial inter-membrane space and cytosol, generating the mitochondrial membrane potential, sucrose exchange, mitochondrial [Ca2+] dynamics and ADP/ATP influx/efflux. All these functions ensure the redox state balance and energy production between mitochondria-cytosol compartments.1-4 From a structural point of view, hVDAC1 presents an extensive secondary structure formed by 19 β-strands that are formed within residues 25–283 and a short α-helix located at the N-terminus containing residues 6–10. The N-terminal α-helix is folded horizontally inside the barrel wall approximately at the midpoint of the hydrophobic portion of the membrane with positively charged residues facing the barrel interior, whereas negatively charged interacting with the barrel wall. The height of the hVDAC1, including the loops, is approximately 30 Å and a circular conformer of hVDAC-1 has an open diameter of about 25 Å.5 Numerous reports have pointed out that hVDAC1 is involved in mitochondrial dysfunction leading to pathological conditions such as calcium overload in mitochondrial matrix that constitutes a common element condition in Alzheimer, Parkinson’s disease, epilepsy and many other chronic pathological processes.5-8 hVDAC1 is an essential component of the induced structure of mitochondrial permeability transition pore (MPTP), a multi-protein complex that is directly involved in apoptosis-based mitochondrial dysfunction.5-9 Since their discovery in the 20th century carbon nanomaterials have attracted a great deal of interest in the scientific community, stemming mainly from their wide variety of applications. Particularly, single walled carbon nanotubes (SWCNT), due to their flexible nature and versatility for chemical functionalization are being used in new nanomedicine applications, such as active ingredients and pharmaceutical excipients forward the design of several drug-delivery systems.10, 11 The main challenge for the safe use of these nanomaterials is the lack of data about their potential nanotoxicity,12 (mitochondrial channel toxicity-induced by SWCNT), a challenge that has prompt the use of modeling tools that can make an effective use of smaller datasets in the first instance. In this context, chemoinformatic methods may be also very useful to predict structure-property relationships for the carbon nanotubes or other nanomaterials.11 Concerning this, the so-called Quantitative Structure-Binding Relationships (QSBR) are used for models that employ ligand-protein docking interactions data and by analogy, for the particular case of nanomaterials, these are called Nano-QSAR (NQSBR) models. In general, these NQSBR models use as input physicochemical properties of nanomaterials (nanodescriptors) to predict their biological activity, nanotoxicity, and/or binding-affinity to specific targets.12-16 Despite the 4


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


relative progress in the discovery of new nanodescriptors for modeling the toxicity caused by carbon nanotubes (i.e. mitochondrial nanotoxicity), the influence of topological vacancies has not been explored until now. Indeed, in most cases, a broad spectrum of topological vacancies are spontaneously generated in the side-walled SWCNT-surface derived from chemical oxidation processes (e.g. by OH, COOH), like vacancies, Stone–Wales and octagon–pentagon pair vacancies.17-21 Such vacancies significantly affect the geometric and electronic properties of SWCNTs and, at the same time, may increase their chemical reactivity and nanotoxicological potential.12 In fact, topological SWCNT-vacancies are often seen and may change redox-biochemical mechanisms at the subcellular level (mitochondria).10, 22 Previous theoretical studies have shown that the most stable vacancies correspond to monovacancy, divacancy, tetra-vacancy and hexa-vacancy in terms of chemical reactivity.23 Following this idea, we considered that the inhibition of the hVDAC1 channel by SWCNT with different # of topological vacancies could be an attractive strategy to tackle the channel nanotoxicity (i.e. hVDAC1 mitotoxicity-induced by carbon nanotubes with different reactivity). In view of that, the implementation of chemoinformatic tools such as molecular docking simulations (DS) can be an efficient in silico approach for the prediction/assessment of the potential nanotoxicity of SWCNTs based on relevant structural descriptors like the number of topological vacancies.24,25 Previous in vitro studies have demonstrated that SWCNT-vacancy can induce mitochondrial dysfunction (mitotoxicity) after interacting with mitochondrial membrane channels like: K+-channels, Na+/Ca2+-channels, mitochondrial ADP/ATP carrier and mitochondrial voltage-dependent anion channels.26-28 Due to this, the use of molecular DS can be a powerful platform for the elucidation of interaction mechanisms and the rational design of new carbon nanomaterials like SWCNT, as well as to efficiently identify hidden structural properties (as topological vacancies) linked to biophysical interactions with a low computational cost.29, 30 The main objectives of this study are to theoretically evaluate and predict the influence of a broad spectrum of SWCNT-topological vacancies to the mitochondrial channel nanotoxicity (hVDAC1-mitotoxicity). Firstly, molecular docking simulations based on density functional theory (DFT) computations are employed to predict the strength of SWCNT/hVDAC1-interactions from the Gibbs free energy of bindingaffinity (FEB) values. Then, the latter are used as a nanotoxicological endpoint forward-applying a new Vacancy-Quantitative Structure-Binding Relationships (V-QSBR) modeling approach. Lastly, the results of the present study shall open new opportunities towards the use of chemoinformatic tools based on ligand-mitotarget docking interactions for making regulatory decisions in nanotoxicology.

5



Page 6 of 29

Methods Docking simulations. To evaluate the interaction between the hVDAC1 and the SWCNT-vacancies family members (having from 1 to 16 topological vacancies), a theoretical study based on docking simulations was performed. The first step entailed preparing the hVDAC1 macromolecule structure files (receptor), obtained from the RCSB Protein Data Bank (PDB) X-ray structures,31 that is: the PDB ID: 2jk4 with resolution of 4.1 Å. The methodology includes the removal of the crystallographic water molecules and all the co-crystallized hVDAC1 ligand molecules if any. Further, based on built-in modules, hydrogen atoms are added according to the appropriate hybridization geometry, as well as partial charges and protonation states, followed by an assignment of bond orders and a setup of rotatable bonds for the hVDAC1 channel-PDB X-ray structure. In addition, before the molecular docking simulations, the hVDAC1 structure was optimized by using the AutoDock Tools 4 of the software AutoDock Vina.24, 32 In the second step, all the SWCNT-vacancy structures (ligands) containing topological vacancies (from v = 1 to 16) were carefully modeled and optimized starting from a zig-zag-SWCNT with Hamada index n = 8 and m= 0. The different SWCNT-vacancies were introduced by removing sp2-C atoms from the side-wall of the zig-zag-SWCNT pristine structure (v = 0). Details of the obtained SWCNT with vacancies are given in Figure S1 of the Supporting Information (SI). Next, in order to characterize the different electronic and structural properties of the several SWCNT_v systems understudy, periodic Density Functional Theory (DFT) calculations were carried out using the SIESTA software code.33 These DFT calculations were performed within the Perdew-Burke-Ernzerhof generalized gradient approach (GGA) exchange correlation potential34 to obtain the electron density,34 while the interactions between the nucleus and the valence electrons were described by the pseudopotential of Troullier and Martins.35 The SWCNT_v orbitals were accessed using a double-zeta polarization (DZP) basis set, which describe correctly the geometry and the electronic characteristics of the evaluated SWCNT_v-systems. The structures of the SWCNT_v systems were relaxed until the residual forces were less than 0.05 eV / Å for all involved atoms. In these calculations, a radius of 200 Ry to mesh-grid interactions was used for the representation of electronic space in real space, whereas the numerical integration in the reciprocal space was done considering a 112 Monkhorst-Pack grid of special kpoints.35,36 For each SWCNT_v optimized system, the electronic band structures were calculated, and we found that the valence band maximum (VBM) and the conduction band minimum (CBM) are both located at the Fermi level (see Figure S2 of the SI). The vacancy formation energy (Ef) for each SWCNT_v tested was calculated using the following equation: 6


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


𝐸𝑓(𝑆𝑊𝐶𝑁𝑇𝑣 = 1 𝑡𝑜16) = 𝐸𝑡(𝑆𝑊𝐶𝑁𝑇𝑣 = 1 𝑡𝑜 16) ― 𝐸𝑝(𝑆𝑊𝐶𝑁𝑇𝑣 = 0) ― 𝑛(sp2 ― 𝐶𝑎𝑡𝑜𝑚𝑠)µ(sp2 ― 𝐶𝑎𝑡𝑜𝑚𝑠) (1) where n is the number of removed sp2-carbon atoms (C-atoms) from the SWCNT-side wall, Et is the total energy of the SWCNT_v structure with # of vacancies (v) (from v = 1 to 16), Ep is the total energy for SWCNT-pristine (v = 0), and μ is the chemical potential of the C-atoms. To study the SWCNT_v-hVDAC1 docking interactions, we employed the flexible molecular docking implemented in the open source software Autodock Vina developed in 2010 by Trott and Olson.47 The cluster formed by ATP-VDAC-active binding site residues from the N-terminal α-helix domain (i.e. the residues Pro7, Gln199, Gln182, Phe181, Val20, Asp19, Lys15, Gly14, Asp12, and Ala11) and the phosphorylation binding site residues from hVDAC1 (i.e. the residues Tyr10*, Ser16*, and Tyr198*) were examined as flexible residues and the ligands (SWCNT_v) were kept rigid.32 Then, the free energy of binding (FEB) of the formed SWCNT_vhVDAC1 complexes was calculated based on the score function that approximates the chemical potential (ΔGbind). The optimization of these complexes was carried out by applying the Autodock Vina scoring function with default Amber force-field parameters. 37Conformational relaxation (flexible docking) favors a significant gain of enthalpy of SWCNT_v-hVDAC1 complexes non-associated with SWCNT_v intra-molecular deformation or vibrational decrease within the hVDAC1 active sites. As to the potential hVDAC1 active binding sites, these were previously predicted by DeepSite38 also considering the hVDAC1-pore calculated by excluding potential co-crystallized hVDACligands. This former step allowed us to identify and delimit the relevant hVDAC1 catalytic sites, potentially at the van der Waals surface that are likely to bind to a small ligand like the SWCNT-vacancy members. To this end, DeepSite considers all the molecular descriptor related to the receptor (in this case: hVDAC1), through 3D-deep convolutional neural networks validated using an extensive test set based on > 6500 proteins from the scPDB database.38, 39 The binding pocket predictions as well as the hVDACvolumetric map predictions40 were used to establish the Cartesian coordinates of the docking box simulation like hVDAC1 grid box size, with dimensions of X=44Å, Y= 16Å, Z= 18Å and the hVDAC1 grid box center X=24.1Å, Y= 10.4 Å, Z= - 4.0 Å. As to the SWCNT_v-hVDAC interactions, these were determined by considering the ATP biophysical environment (hVDAC1 active site) as reference control to evaluate the SWCNT_v and the affinity towards the hVDAC1-channel. Several runs starting from random conformations were performed, and the chosen iterations for each run adapted to the problem complexity. An exhaustiveness option set to 80 (average accuracy)32 in each docking calculation was used, since it was verified that setting higher values for this parameter (e.g. 100) highly increased the docking simulation time but gave the same FEB results. 7



Page 8 of 29

The Gibbs free energy of binding output results (or FEB values) is therefore defined by the ΔGbind affinity values for all docked poses of the formed SWCNT_v and hVDAC1 complexes, including the internal steric energy of a given ligand (SWCNT_v with v = 1-16) expressed as the sum of individual molecular mechanics force field terms, such as van der Waals interactions, hydrogen bonding, electrostatic interactions, and intramolecular ligand interactions.32 Additional details about the thermodynamic docking force field parameters applied in this study can be consulted in Autodock Vina scoring function.32 Docking results were found as energetically unfavorable when the FEB-value for the complexes was greater or equal than 0 kcal/mol (worst crystallographic pose), thus showing either an extremely low or complete absence of affinity (repulsive pattern of interactions). According to this criterion, SWCNT_v conformations with the lowest Gibbs docking free energy of binding (FEB negatives value) were obtained. The best root-mean-square deviation (RMSD) defined below (eq. 2) was stored as a criterion of correct docking pose accuracy for atomic positions (pi-SWCNT_v, Pj-hVDAC1) below 2Å.41, 42

𝑅𝑀𝑆𝐷(𝑝(𝑖),𝑝(𝑗)) =

∑ (𝑎𝑡𝑜𝑚(𝑖) ― 𝑎𝑡𝑜𝑚(𝑗))2 𝑛 𝑛

(2)

The next step consists of analyzing the results obtained from the docking approach with respect to the final Gibbs free energy of binding for the SWCNT_v-hVDAC1 complexes of each simulation. The minimum interatomic distances were calculated between hVDAC1-atoms from the critical amino acids of the hVDAC1-receptor and SWCNT_v-atoms at the best crystallographic binding position for ligands (SWCNT_v-family) The only considered interatomic distances were related to ATP-hVDAC1 active binding site residues from the N-terminal α-helix domain (i.e.: residues Pro7, Gln199, Gln182, Phe181, Val20, Asp19, Lys15, Gly14, Asp12, Ala11) involved in the ATP transport to compare with the ATPhVDAC1 interatomic distances as reference control of potential ATP efflux inhibition. The ATP ligand model (C10 H16 N5 O13 P3) in .sdf format was retrieved from PubChem (CID: 5957). In addition, we calculated the interatomic distances of SWCNT_v related to hVDAC1-phosphorylation binding site residues (Tyr10*, Ser16*, and Tyr198*). For this instance, the Euclidean distances (𝑑𝑖𝑗: SWCNT_vhVDAC1 interatomic distances) were calculated for all atoms in the SWCNT-vacancies to all atoms in the hVDAC1 channel. To evaluate the FEB SWCNT_v(i)-hVDAC1(j) interactions, a cutoff value of 7 Å was used, i.e. all atoms with distances dij ≤ 7Å were considered as interacting-atoms (see Table S1 of the SI).41, 42Considering

the interatomic distance (minimum Euclidean distance, dij) between i and j interacting atoms,

the Cartesians coordinates (hVDAC1_x(j), hVDAC1_y(j), hVDAC1_z(j)) were taken from all the j-atoms of the critical amino acids of the hVDAC1-receptor binding site as an input, while the coordinates 8


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


(SWCNT_v x(i), SWCNT_v y(i), SWCNT_v z(i)) of all the i-atoms of SWCNT_v systems evaluated were taken as output (see Eq. 3) 𝑑𝑖𝑗 = (𝑥(𝑖) ― 𝑥(𝑗))2 + (𝑦(𝑖) ― 𝑦(𝑗))2 + (𝑧(𝑖) ― 𝑧(𝑗))2

(3)

Statistical analysis. In order to explore the relationships between FEB-values of the formed docking complexes (SWCNT_v-hVDAC1) with the structural parameters (# of SWCNT-vacancies and the SWCNT-vacancy formation energy (Ef)), different statistical regression models (linear, second order polynomial and plateau followed by one phase decay function) were implemented to seek the best fit curves. V-QSBR model. To develop the Vacancy-Quantitative Structure-Binding Relationships (V-QSBR model), 3D-molecular SWCNT_v descriptors (more than 500) were calculated using the DRAGON 7.0 software,43 and belonging to the following main classes: 3D matrix-based descriptors, 3D autocorrelations, 3D atom pairs, and CATS 3D. The more relevant SWCNT_v descriptors included in the V-QSBR model were found to belong to the 3D matrix-based class. These nanodescriptors are commonly known as topographic indices and are defined from the graph representation of molecules but using geometric Euclidean distances between atoms instead of topological distances. For more information regarding the 3D matrix-based descriptors definition, please refer to the Dragon user’s manual.43 A binary V-QSBR classification model was setup using the Linear Discriminant Analysis (LDA) algorithm implemented in the Statistica 12 software,44 and by applying a forward stepwise regression procedure to select the best 3D-molecular SWCNT_v descriptors to be included in the final model. The original dataset of 17 SWCNT_v (SWCNT-pristine_v = 0 and SWCNT_v = 1-16) was split into two binary classes of SWCNT_v-hVDAC1 nano-interactions: strong and weak. For this instance, a threshold value of 7.0 kcal/mol was used for the free energy binding obtained from the molecular docking of the ATP-ligand (native hVDAC1-substrate) with the hVDAC1-active binding site. In so doing, if the FEB values obtained were lower than the cut-off, the SWCNT_v were assigned to the class of strong hVDAC1 interactions, otherwise they were considered as having a weak hVDAC1 interaction. To assess the goodness-of-fit of our V-QSBR classification model, a deep validation process was carried out by means of the crossvalidation procedure, as well as by assessing both the specificity and sensitivity (percentage of weak and strong SWCNT_v interactions correctly classified, respectively), the receiver-operating characteristic (ROC) curve,45 Wilk’s lambda () and the Matthews correlation coefficient (MCC).46 In order to perform the cross-validation procedure, we have randomly chosen 10 out of 13 cases for the training series and 3 9



Page 10 of 29

out of 13 cases for the validation series in the strong interaction class. Additionally, we randomly assigned 3 out of 4 cases for the training series and 1 out of 4 cases for the validation series in the weak class. This means that the V-QSBR model was built using 9 training cases and validated using 4 validation cases. Then, we carried out a cross-validation test using the specific module included in the LDA of the STATISTICA software.44 Notice that cross-validation is a classical procedure for assessing the performance of a statistical analysis, based on performing the analysis on the training set, and then validating the results with the validation subset. More details on the V-QSBR model dataset can be seen in Table S2 of SI, where it is specified if a given SWCNT_v belongs to the training or the validation series.

Results and discussion Mechanistic interpretation of molecular docking. An important task to ensure accuracy of our theoretical data is to check the prediction of feasible hVDAC1 binding sites. As referred to before, in the present study the prediction of binding active sites of hVDAC1 was performed using a machine learning tool based on deep convolutional neural networks (DeepSite-CNNs cheminformatics’ tool) with maximum score = 1.0.38, 39 The resulting hVDAC1 binding sites are shown in Figure 1.

10


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


Figure 1. A) DeepSite prediction of the topology of hydrophobic cavities for the hVDAC1 active binding sites (volumetric orange regions); B) hVDAC1 binding sites depicted like van der Waals surfaces (volumetric orange regions); C) Representation of the X-ray crystallographic molecular structure of the voltage-sensing N-terminal-α-helix segment; and D) Structural analysis of the voltage-sensing N-terminalα-helix segment, considering relevant properties like structural class (helix, turn and sheet) of the residue sequence (from 1 to 35) and conformational properties based on the Ramachandran plot scoring, the latter showing that all voltage-sensing N-terminal ɑ-helix segment residues are conformationally favored according to the torsion dihedral angles defined by the different (Psi) and (Phi) classes around the peptidebond and within the Ramachandran colored contour (gray color label). Further details on the solvent accessibility for different residues of the voltage-sensing N-terminal-α-helix segment are also shown (blue shades). Let us now evaluate the relationships between topological vacancies of the SWCNT_v-family members and their influence in the binding affinity to hVDAC1, based on the Gibbs free energy of binding obtained from molecular docking simulations. In fact, it is crucial to understand how these vacancies can modulate the biochemical interactions since it is well known that, during the synthesis of SWCNT, topological vacancies can be generated spontaneously,47 particularly at the subcellular level.22 In this regard, relevant information was obtained due to a solid correlation between the Gibbs free energy of binding and the chosen SWCNT_v nanodescriptors (i.e. the # of SWCNT-vacancies (v = 0-16) and the vacancy formation energy Ef. The pure pristine SWCNT_v = 0) has a stable structure with sp2-sp3 bonding hybridization with π and π*-bonding perpendicular pz orbitals, which can be quite large specially for small diameter SWCNTs23 as those here studied (zig-zag SWCNT(8,0) with vacancies). However, the appearance of vacancies induces some valence electrons and bonding with other sp2-C atoms. Hence, the variation of Ef with the # of SWCNT-vacancies can be explained by the stabilization of the SWCNT_v conﬁgurations due to a fully relaxation of the sp2-C atoms around the vacancy and related geometrical distortions of the SWCNT_v systems (see Table 1). The increment of Ef for the highest # of vacancies essentially results from the existence of more sp3-C-atoms with dangling bonds, which in turn are expected to produce a larger increase in the reactivity of these SWCNTs_v.19, 23 Besides, the presence of different # of topological vacancies in the SWCNT-family studied also affects the values of the spin polarization (∆S), which vary from 0 to 9.94 μB (Table 1).

11



Page 12 of 29

Table 1. DFT-characterization parameters of the SWCNT-family with topological vacancies, namely: the vacancy energy of formation (Ef), spin polarization (∆S) and band gap energy of CBM and VBM (Gap). # of SWCNT-vacancies

Ef (eV)

∆S (µB)

Gap (eV)

SWCNT_v = 0

0

0

0.66

SWCNT_v = 1 SWCNT_v = 2 SWCNT_v = 3 SWCNT_v = 4 SWCNT_v = 5 SWCNT_v = 6 SWCNT_v = 7 SWCNT_v = 8 SWCNT_v = 9 SWCNT_v = 10 SWCNT_v = 11 SWCNT_v = 12 SWCNT_v = 13 SWCNT_v = 14 SWCNT_v = 15 SWCNT_v = 16

6.00 4.69 10.83 12.50 14.34 15.60 14.39 14.56 14.25 16.56 17.63 19.35 21.59 22.40 24.67 25.14

0 0 1.60 1.79 2.31 5.86 5.02 3.69 2.53 3.79 5.15 5.75 4.43 7.89 6.39 9.94

0.84 0.51 0 0.25 0.24 0.36 0 0 0 0 0 0.15 0 0.10 0 0

Altogether, the interaction results suggest that more stable docking complexes are formed as the # SWCNT-vacancies and vacancy formation energy increase, showing an R2 of 0.93 and 0.96 for FEB and # SWCNT-vacancies and FEB and vacancy formation energy relationships, respectively with p < 0.05 (see Figure 2).

12


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


Figure 2. Relationships found between the hVDAC1-binding affinity (FEB) and the nanodescriptors used for the SWCNT_v-family member studied. A) FEB values vs. # of SWCNT-vacancies; the dotted red line shows the free energy of binding of the ATP native substrate of hVDAC1 (= 7.0 kcal/mol), used as control reference for the docking simulations. B) FEB values vs. the vacancy formation energy Ef. The best fit regression curves were obtained by statistical regression analysis, based on second order polynomial (A) and linear (B) functions (p < 0.05). The SWCNT_v with v = 0, 9 and 16 are depicted to show examples of systems with low, medium and strong docking affinity interactions. Furthermore, it was observed that the topological vacancies induce appreciable changes in the structural and electronic properties in all SWCNT_v systems evaluated due to the fact that, when the atoms of the 13



Page 14 of 29

structure are removed, these spontaneously reorganize forming more new stable structures.23 In addition, the gap between the conduction band minimum and the valence band maximum (CBM-VBM values) shows that those go from semi metallicSWCNT-pristine with v = 0 (0.66 eV)  to metallic  SWCNT with v =16 (0 eV)  behavior. Regarding this, the presence of topological vacancies in the side-wall of the SWCNTs tested could mimic the recognized phenomenon of “edge effects” that only appears in semiconducting zig-zag topologies of SWCNTs which present open extremes from the cycloparaphenylene aromatic system; the latter playing an important role in the mitochondrial channel nanotoxicity as showed in previous studies.23, 48 It should be noticed here that the band gap energy values for the evaluated SWCNT_v family members have a non-periodic behavior with the vacancy size increase (see Table 1), mainly due to the structural reorganization in the SWCNT_v side-wall when the carbon atoms are removed. In fact, from a structural point of view, the number of SWCNT vacancies alters the Fermi energy depending on their relative positions, and that drastically modifies the SWCNT_v electronic properties leading to a non-uniform pattern with the increasing of topological vacancies size. In fact, the vacancy formation energy induces a flat-band magnetism in the Fermi level and the π-electrons contribute to magnetic polarization as well as to the increase in the intrinsic vacancy reactivity due to more stable docking interactions (more negative FEB values).19, 23, 48 Therefore, the docking results showed interesting aspects about the influence of SWCNT-topological vacancies (# SWCNT-vacancies) to potentially modulate the VDAC channel activity. The defined SWCNT vacancy-moieties seem to be new relevant nanodescriptors to explain the interactions occurring within the hVDAC1-hydrophobic binding active site, which is formed by cationic amino acid residues (histidine, lysine, and arginine) involved in each of the phases of mitochondrial ATP-transportation.1-3, 49 In addition, considering the Euclidean distances between the centroids of atoms from SWCNTs_v and hVDAC1, we suggest that topological vacancies affect the permeability and bioenergetic function of the hVDAC1-channel. Indeed, the taken cut-off distance value for interatomic distances (dij ≤ 7Å) implies a high probability for i,j-interactions and a low probability of obtaining false contacts.39 An illustration of the critical docking interactions with the # of vacancies and Ef parameters for SWCNTs with v = 0, 9 and 16 are displayed in Figure 3. As can be seen, the topological vacancies are shown once more to be valuable SWCNT nanodescriptors, since they predict the mitochondrial channel interactions (hVDAC1) with high significant correlation, judging from the obtained determination coefficient (R2 = 0.98; p < 0.05) of a second order polynomial function based on both # of vacancies and Ef parameters. Here, it is also important to remind that the FEB docking results were obtained considering the diameter and zig-zag chirality (Hamada index: n = 8; m = 0) as invariant properties for all tested SWCNT_v systems. 14


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


Figure 3. Representative snapshots of the obtained docking complexes following an increment in the critical # of SWCNT-vacancies (v = 0, 9, 16): A) SWCNT_v = 0 with hVDAC1, B) SWCNT_v = 9 with hVDAC1 and C) SWCNT_v = 16 with hVDAC1. In all cases, it is shown the voltage-sensing N-terminal α-helix of hVDAC1 involved in the ATP-transport. D) 2D-contour plot showing the relationships between the FEB-values of all the formed docking complexes and the chosen parameters (# of SWCNT-vacancies and vacancy formation energy Ef). The function FEB(Z) was obtained by a statistical regression analysis based on a second order polynomial function relating the FEB values with those parameters (R2 = 0.98). Dark red corresponds to low interaction affinity (≈ 0 kcal/mol) and orange to high affinity (≈ -20 kcal/mol). As can be seen in Figure 3, the hVDAC1 affinity for the SWCNT with v = 16 is higher (more negative FEB value) compared to its SWCNT congeners with v = 0 and v = 9. The high proximity of the location of the SWCNT vacancies to the voltage-sensing N-terminal ɑ-helix segment could explain the docking-mechanisms based on both the vacancy properties and the interatomic distances to the key hVDAC1-functional residues (Pro7, Gln199, Gln182, Phe181, Val20, Asp19, Lys15, Gly14, Asp12, Ala11) as well as with the hVDAC1-phosphorylation binding site amino acid residues 15



Page 16 of 29

(Tyr10*, Ser16*, and Tyr198*).9 Notice that, for most of the SWCNT_v systems, the interatomic distances for the SWCNT_v-hVDAC1 interactions with the amino acid residues cited above were generally very close but, in some cases, lower than the interatomic distance of these residues for the ATP-hVDAC1 interaction. The latter case was verified at least for the SWCNT-vacancy systems with the highest number of topological vacancies (from SWCNT_v with v = 13 to 16; see Figure 4).

Figure 4. Radar plots showing the relative interatomic distances (in Å) of the SWCNT_v with different # of topological vacancies (v) after docking simulations. For comparison purposes, the interatomic distances from the ATP-VDAC-N-terminal α-helix domain active binding site residues are included (red line). The hVDAC1-phosphorylation binding site residues are marked with asterisks (*). Overall, the presence of vacancies in the SWCNTs could theoretically lead to interactions with any amino acid present in the N-terminal α-helix domain, which in turn are ideally positioned to modulate metabolites flux and ions conductance (Ca2+) through the hVDAC1 channel when the interatomic distance is less than 7Å.40,41 So, topological vacancies of SWCNTs are expected to induce channel mitotoxicity through the aberrant ATP-transport efflux just as it occurs in various diseases, such as cancer-induced mitochondrial 16


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


dysfunction where the hVDAC1 channel plays an important role in the redox and bioenergetics balance in several pathophysiological circumstances.8, 9, 50 According to the theoretical results, it should be reasonable to hypothesize that the increment of the # of SWCNT-vacancies with sp3-C-atom with dangling bonds (with negative charges) can establish more stable complexes with positive residues of hVDAC1 and, in this way, could disrupt the association of the ATP4anion (phosphate tail) with the positively charged ε-amino group (N+ primary amines) of Arg18 or Lys15 residues present at the hVDAC1-bottom cavity (N-terminal ɑ-helix segment). We suggest that this interaction critically depends on the # of SWCNT-vacancies, explaining the low hVDAC1 affinity values observed for the SWCNT_v = 0. Here, it is important to note that the SWCNT-vacancy interatomic distance values relative to the hVDAC1phosphorylation binding site residues (Tyr10*, Ser16*, and Tyr198*) shows a similar trend in terms of magnitude to the interatomic distance values relative to ATP-hVDAC1 active binding site (with similar or lower distance values), based on the applied Euclidean distance criterion (dij < 7Å). Particularly, a low interatomic distance was observed for the regulatory Tyr10*, which is also located in the N-terminal αhelix segment.9 These results were verified for all the SWCNTs with the highest number of topological vacancies (SWCNT_v ≥ 9). Phosphorylation is an important mechanism for controlling the physiological function of proteins involved in cellular signaling pathways. hVDAC1 post-translational modifications are associated with pathophysiological conditions such as metabolic disorders, cancer, and cardiovascular diseases.8,9 So, we have identified an interesting relationship linking low interatomic distances from hVDAC1-phosphorylation binding site residues of the N-terminal α-helix segment with an increment on the # of SWCNT-topological vacancies, which in turn has important implications from the point of view of the structure-binding affinity relationships for single walled-carbon nanotubes. In addition, an unusual critical number of topological vacancies (SWCNT_v = 9) from which the interatomic SWCNT_v(i)hVDAC1(j)-phosphorylation residues distance significantly decreases with the increase of SWCNT topological vacancies was identified after statistical regression analysis based on a plateau followed by one phase decay model (see Figure 5).

17



Page 18 of 29

Figure 5. Correlation found between the # of SWCNT-vacancies and the interatomic distances (Å) of critical hVDAC1-phosphorylation binding site residues like: A) Tyr10* (red line; R2 = 0.93), B) Ser16*(purple line; R2 = 0.93), C) Tyr198* (green line; R2 = 0.95). In all cases, data were fitted to a plateau followed by one phase decay model. Note that, for the three hVDAC1-phosphorylation binding site residues examined, the interatomic distance significantly decreases (p < 0.05) from a cutoff-vacancy number (X0 or v ≈ 9). Herein the Arg18 (black) from the hVDAC1-active binding site was used as control to show the absence significant correlation between the # of SWCNT-vacancies for a non-regulatory 18


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


hVDAC1-phosphorylation binding site residue with interatomic distance of 3.5Å. On the right side, a Pymol cartoon representation is shown of the formed complex with each hVDAC1-phosphorylation binding site residues and the corresponding interatomic distance for the cutoff-vacancy number (X0 or v = 9) for Tyr10* (red), Ser16* (purple) and Tyr198* (green). Also, note that the Arg18-hVDAC1-residue (black) of the N-terminal α-helix domain interacts with the vacancy center of mass of the SWCNT_v (light green) with vacancies SWCNT_v = 9 (light blue). Therefore, we suggest that the sp3-C-atom with dangling bonds from SWCNT_v systems could drastically affect the phosphorylation mechanisms of serine and tyrosine of the N-terminal α-helix domain of the hVDAC1-channel with several implications in the mitochondrial bioenergetics (ATP efflux).42 Further, the SWCNT_v systems could affect the intrinsic flexibility of the hVDAC1 residues from the voltagesensing N-terminal α-helix segment that plays an important physiological role in the mitochondrial ATPtransport. It is well known that vibration modes between residue pairs (i,j-residues flexibility) are linked to the distance residue fluctuations and deformation energy, which in turn are a function of the changes in the residue flexibility induced by the ligands like SWCNT_v family members (SWCNT_v with v = 9; see Figure 6).

19



Page 20 of 29

Figure 6. Distance matrix-cross correlation map (Cij) showing the correlations between residue fluctuations from hVDAC1/Voltage-N-terminal -helix segment, plotted as a function of residue indices i and j. A) Residues with unperturbed-flexibility of hVDAC1/Voltage-N-terminal ɑ-helix segment (Pro7, Val20, Asp19, Lys15, Gly14, Asp12, Ala11) in the absence of SWCNT_v systems. B) An example of the aforementioned residues from hVDAC1/Voltage-N-terminal -helix segment with perturbed-flexibility by docking interactions with SWCNT_v = 9. The pairs subject to fully correlated motions or strong correlations (0.2 ≤ Cij ≤ 1) appear in red color indicating the same direction for the i and j residuesflexibility, while the flexibility based on anti-correlated motion (i.e.: Cij ≤ 0) are colored blue (oppositedirection for the fluctuation-motions of residues i and j) and the moderately correlated and uncorrelated (Cij ≈ 0) regions are colored in light red and blue, respectively. This behavior is representative for all performed docking simulations. 20


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


Inspection of figure 6 lead us to suggest that the molecular docking mechanism of SWCNT_v family members (v = 1-16) are most likely based on local-perturbations (interactions) with effector residues from the hVDAC1/voltage-N-terminal -helix segment. In fact, the transitions between the functional states (hVDAC1open/hVDAC1 close) are strongly dependent on the hVDAC1 flexibility properties (see Figure 6) and require the biochemical and biophysical regulation based on allosteric inter-communication residues.43 Following this idea, we performed a theoretical modeling on local perturbation-induced by SWCNT-vacancy family members only considering the hVDAC1 binding residues that are simultaneously affected by the same anisotropic-fluctuations according to low flexibility normal modes, allowing allosteric communication across large distances in its voltage-N-terminal -helix segment (or ATPcatalytic binding site). This fact explains the different FEB results found for the several SWCNT_v-family members here studied, by considering minimal differences in the strength of residues communication (hVDAC1 local perturbations-based interaction) when the SWCNT_v ligands are present and inducing highly anisotropic behavior in the residues of the hVDAC1 voltage-N-terminal -helix segment based on non-covalent docking interactions (see Figure S3 of the SI). V-QSBR Model to predict the strength of hVDAC1-interactions To evaluate the relationships between the influence of # of SWCNT-vacancies that are relevant to predict the strength of hVDAC1-interactions, a Vacancy-Quantitative-Structure Binding-Relationship (V-QSBR) approach was applied for the prediction of the docking free energy of binding (or FEB) values between SWCNT- vacancy and hVDAC1. The output of the V-QSBR model is the scoring function f(FEB). The main aim of the herein reported V-QSBR model was to find a set of vacancy descriptors related to the topological SWCNT-vacancies that were able to clearly connect the structure of the SWCNT_v systems with their nanotoxicity-based binding interactions. In this context, the variable selection was the key point of this theoretical procedure. The developed V-QSBR model is able to correct classify all the 17 SWCNT_v systems (from v = 0 to v = 16), with a sensitivity and specificity of 100%. In addition, 13 out of 13 were predicted as weak while 4 out of 4 as strong. The resulting best fit equation found using linear discriminant analysis (LDA) is shown below: 𝐹𝐸𝐵 = 181.3 ∗ 𝐶ℎ𝑖𝐴_𝑅𝐺 ― 36.7

(𝟒)

As can be seen, the V-QSBR model is based only on one 3D-molecular descriptor of the SWCNT_v systems. This nanodescriptor is the 3D-average Randić-like index, a geometrical descriptor included into 21



Page 22 of 29

the block of 3D matrix-based descriptors that are commonly known as topographic indices.51 The ChiA_RG being calculated from the graph representation of SWCNTs_v but using geometric Euclidean distances between their atoms instead of topological distances. The reciprocal squared geometrical distance matrix (RG) collects reciprocal squared distances between all the between the i-th and j-th atom pairs of the system under study,51 such as the SWCNT_v systems with different topological vacancies (0 ≤ v ≤ 16). Thus, considering that the 3D-molecular structure (SWCNT_v) and the type of nanodescriptor (ChiA_RG) found as the most relevant is straightforward to suggest how the # of vacancies of the SWCNT_v family studied can directly influence their potential hVDAC channel nanotoxicity based on the binding interactions. The ChiA_RG nanodescriptor allows decoding the corresponding structure for each SWCNT_v system evaluated, since it is independent of the # of SWCNT_v atoms, i.e. the SWCNT_v size, and it is unique regarding the 3D arrangement of the SWCNT_v atoms; the latter being an invariant property according to the SWCNT translational and rotational invariance of the entire SWCNT_v structure. It is important to note that the calculated ChiA_RG show a low degeneracy in the SWCNT_v family studied according to the 3D-nanodescriptor properties while the number of removed atoms of the SWCNT (# of vacancies) shows a high degeneracy according to constitutional descriptor properties. The degeneracy refers to the ability of a descriptor to avoid equal values for different molecules in a dataset. Then, as a control test, Figure 7 shows a plot of the relationship between FEB-values and both ChiA_RG and # of SWCNT-vacancies.

22


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


Figure 7. 2D-contour plot based on the relationships between FEB-values in all the formed docking complexes and both the constitutional vacancy descriptors (# of SWCNT-vacancies) and the 3D-average Randić-like index ChiA_RG ones for each modeled SWCNT-vacancy system. The function (Z) or free energy of binding (FEB-values) was obtained by a statistical regression analysis based on a second order polynomial function of such descriptors (R2 = 0.93). Green corresponds to high strength of docking interaction affinity, and red to low strength interaction affinity. It should be mentioned here that, the ChiA_RG 3D-nanodescriptor contains relevant information about interatomic distances in the entire SWCNT-system and structural properties like bond distances in the presence of vacancies, rings, planar and non-planar systems,43 as well as atom types such as sp3-C-atoms with dangling bonds, which are expected to occur with the vacancies’ increase. It is thus expected to reproduce well the FEB behavior previously obtained for the # of vacancy and Ef parameters (Figure 3). In fact, the derived V-QSBR model (eq. 4) is able to specifically discriminate between the SWCNT_v with topological vacancies that strongly bind to hVDAC1 (FEB < 7.0 kcal/mol) from those having a weak hVDAC1-binding affinity (FEB ≥ 7.0 kcal/mol), and that only by including the 3D-average Randićlike index of the SWCNT_v structures. In addition, it is important to note that, the FEB value of 7 kcal/mol refers to the Gibbs free binding of the ATP-native substrate of hVDAC1 that was used as reference control 23



Page 24 of 29

for comparison purposes with the docking simulation experiments. Moreover, regarding the validation of the derived V-QSBR model, by performing a classical cross-validation test, it was found that this model correctly classifies 4 out 4 cases. In addition, as can be seen in Figure 8, the areas under the ROC curves for both training and validation sets is one, meaning that our model is clearly not a random classifier (area = 0.5). Furthermore, the calculated MCC45 of the V-QSBR model was, as expected, equal to one, along with the sensitivity and specificity that were also equal to 100%. Finally, the main statistical validation parameters of the predictive model including the 2 parameter and other relevant statistical parameters are shown in Table S3 of the SI.

Figure 8. Results for the ROC curves analysis of the V-QSBR model. Lastly, we have also checked whether the developed V-QSBR model comply with the Organization for Economic Co-operation and Development (OECD) guidelines for QSARs.52 According to these OECD guidelines, the following principles must be obeyed for a successful QSAR modeling, and further allowing the derived models to be used in regulatory purposes. Herein, one can easily see that our V-QSBR model complies with such principles, since it has: 1) a definite dataset with a well-defined end point  the prediction of the channel mitotoxicity of SWCNT_v; 2) an unambiguous algorithm for its development  24


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


a linear equation with only one variable; 3) a well-defined domain of applicability and appropriate measures of goodness-of-fit, robustness and predictivity; and finally 4) a proper mechanistic interpretation.

Conclusions The binding-affinity properties of single-walled carbon nanotubes with different number of topological vacancies (v: from 0 to 16) with human mitochondrial voltage-dependent anion channel 1 were explored using molecular docking simulations. The results showed that the free energy of binding (FEB) was significantly more favorable with the increase in the # of SWCNT-vacancies and vacancy formation energy (Ef). Besides, we observed a high correlation between the vacancy-parameters (# of SWCNTvacancies, Ef) and the FEB-values obtained for the docking complexes of hVDAC1 with SWCNT_v systems (v = 0-16), describing a singular behavior like plateau followed by one phase decay. These theoretical results thus provided detailed information about the influence of topological vacancies in the interaction mechanism that governs the mitochondrial channels nanotoxicity of carbon nanotubes which has largely been ignored until the present. Furthermore, an interesting relationship was found for the SWCNT-vacancy systems involving low interatomic distances from hVDAC1-binding sites and the regulatory phosphorylation residues of the N-terminal α-helix segment (Tyr10*, Ser16*, and Tyr198*) following an increment in the # of SWCNT-vacancies, a result that has important implications from the bioenergetic point of view (mitochondrial ATP-efflux) and permeability function of hVDAC1. Besides, the SWCNT_v family members can affect the flexibility properties of the VDAC N-terminal-α-helix segment by inducing different patterns of local perturbations in the inter-residue communication. Lastly, a new Vacancy-Quantitative-Structure Binding-Relationship (V-QSBR) modeling analysis allowed us to unequivocally classify the strength of docking interactions (i.e.: weak or strong) of the formed SWCNT_v/hVDAC1 complexes with exceptional predictive classification performance (= 100%), based only on the 3D-average Randić-like index (ChiA_RG) of a given SWCNT_v. Moreover, the V-QSBR model derived fully comply with the OECD principles for predictive QSARs. Finally, these in silico results open a way for the rational design of novel carbon nanomaterials and also to disclose nanotoxicological problems considering the influence of topological vacancies in single-walled carbon nanotubes, as well as regarding the use of human proteins like the human mitochondrial hVDAC1 channel. The latter in turn having potential relevance towards identifying optimal benefit-risk relationships with application in mitotarget Precision Medicine and expand the spectrum of structural attributes of carbon nanomaterials to make regulatory decisions in nanotoxicology. 25



Page 26 of 29

Notes: The authors declare no competing financial interests. Funding This work received financial support from Fundação para a Ciência e a Tecnologia (FCT/MEC) through national funds and co-financed by the European Union (FEDER funds) under the Partnership Agreement PT2020, through projects UID/ QUI/50006/2013, POCI/01/0145/FEDER/007265, NORTE-01-0145FEDER-000011 (LAQV@REQUIMTE), and the Interreg SUDOE NanoDesk (SOE1/P1/E0215; UP). RC acknowledges FCT and the European Social Fund for financial support (Grant SFRH/BPD/80605/2011). J.M. Monserrat also acknowledges funds from INCT (Process Number: 421701/2017-0). H.G.D. acknowledges financial support from MINECO (CTQ2016-74881-P) - 2017 - 2019 and Basque Government (Eusko Jaurlaritza) Consolidation Grant (IT1045-16) - 2016 - 2021. P.V.O and S.B.F acknowledge the support from National Processing Center of High Computing Performance/Brazil/SP (CENAPAD-SP).

Supporting Information: Figure S1, Optimized structures of the SWCNT-vacancy systems (v = 1-16); Figure S2, Energy bands and charge plots for the SWCNT-vacancy systems (v = 0-16); Table S1, DFT and docking interaction results obtained for the SWCNT-vacancy systems (v = 0-16); Table S2, V-QSBR model raw dataset; Figure S3. Anisotropic behavior and local perturbation response scanning (LPRS map) induced by SWCNT-vacancy systems (v = 0-16); and Table S3. Relevant statistical parameters obtained for the V-QSBR-LDA model. This material is available free of charge via the internet at http://pubs.acs.org.

26


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


References (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)

(13) (14) (15) (16) (17) (18) (19)

Okada, S. F., O'Neal, W. K., Huang, P., Nicholas, R. A., Ostrowski, L. E., Craigen, W. J., Lazarowski, E. R., and Boucher, R. C. (2004) Voltage-dependent anion channel-1 (VDAC-1) contributes to ATP release and cell volume regulation in murine cells. J Gen Physiol 124, 513-526. Bayrhuber, M., Meins, T., Habeck, M., Becker, S., Giller, K., Villinger, S., Vonrhein, C., Griesinger, C., Zweckstetter, M., and Zeth, K. (2008) Structure of the human voltage-dependent anion channel. Proc Natl Acad Sci U S A 105, 15370-15375. Ujwal, R., Cascio, D., Colletier, J. P., Faham, S., Zhang, J., Toro, L., Ping, P., and Abramson, J. (2008) The crystal structure of mouse VDAC1 at 2.3 A resolution reveals mechanistic insights into metabolite gating. Proc Natl Acad Sci U S A 105, 17742-17747. Choudhary, O. P., Paz, A., Adelman, J. L., Colletier, J.-P., Abramson, J., and Grabe, M. (2014) Structureguided simulations illuminate the mechanism of ATP transport through VDAC1. Nature Structural &Amp; Molecular Biology 21, 626. Colombini, M. (2009) The published 3D structure of the VDAC channel: native or not? Trends Biochem Sci 34, 382-389. Summers, W. A., and Court, D. A. (2010) Origami in outer membrane mimetics: correlating the first detailed images of refolded VDAC with over 20 years of biochemical data. Biochem Cell Biol 88, 425-438. Bernardi, P., and Di Lisa, F. (2015) The mitochondrial permeability transition pore: Molecular nature and role as a target in cardioprotection. Journal of Molecular and Cellular Cardiology 78, 100-106. Pi, Y., Goldenthal, M. J., and Marin-Garcia, J. (2007) Mitochondrial channelopathies in aging. J Mol Med (Berl) 85, 937-951. Martel, C., Wang, Z., and Brenner, C. (2014) VDAC phosphorylation, a lipid sensor influencing the cell fate. Mitochondrion 19 Pt A, 69-77. Foldvari, M., and Bagonluri, M. (2008) Carbon nanotubes as functional excipients for nanomedicines: II. Drug delivery and biocompatibility issues. Nanomedicine: Nanotechnology, Biology and Medicine 4, 183200. Toropova, A. P., Toropov, A. A., Veselinović, A. M., Veselinović, J. B., Leszczynska, D., and Leszczynski, J. (2017) Chapter 8 - Quasi-SMILES as a Novel Tool for Prediction of Nanomaterials′ Endpoints, In Multi-Scale Approaches in Drug Discovery (Speck-Planche, A., Ed.) pp 191-221, Elsevier. Gonzalez-Durruthy, M., Alberici, L. C., Curti, C., Naal, Z., Atique-Sawazaki, D. T., Vazquez-Naya, J. M., Gonzalez-Diaz, H., and Munteanu, C. R. (2017) Experimental-Computational Study of Carbon Nanotube Effects on Mitochondrial Respiration: In Silico Nano-QSPR Machine Learning Models Based on New Raman Spectra Transform with Markov-Shannon Entropy Invariants. J Chem Inf Model 57, 1029-1044. Toropova, A. P., and Toropov, A. A. (2017) Nano-QSAR in cell biology: Model of cell viability as a mathematical function of available eclectic data. Journal of Theoretical Biology 416, 113-118. Singh, K. P., and Gupta, S. (2014) Nano-QSAR modeling for predicting biological activity of diverse nanomaterials. RSC Advances 4, 13215-13230. Fourches, D., Pu, D., Tassa, C., Weissleder, R., Shaw, S. Y., Mumper, R. J., and Tropsha, A. (2010) Quantitative nanostructure-activity relationship modeling. ACS Nano 4, 5703-5712. Concu, R., Kleandrova, V. V., Speck-Planche, A., and Cordeiro, M. (2017) Probing the toxicity of nanoparticles: a unified in silico machine learning model based on perturbation theory. Nanotoxicology 11, 891-906. Fagan, S. B., da Silva, L. B., and Mota, R. (2003) Ab initio Study of Radial Deformation Plus Vacancy on Carbon Nanotubes: Energetics and Electronic Properties. Nano Letters 3, 289-291. Ma, Y., Ma, J., Lv, Y., Liao, J., Ji, Y., and Bai, H. (2016) Effect of mono vacancy defect on the charge carrier mobility of carbon nanotubes: A case study on (10, 0) tube from first-principles. Superlattices and Microstructures 99, 140-144. Yuchen, M., Lehtinen, P. O., Foster, A. S., and Nieminen, R. M. (2004) Magnetic properties of vacancies in graphene and single-walled carbon nanotubes. New Journal of Physics 6, 68. 27



(20) (21) (22)

(23) (24) (25) (26) (27) (28) (29)

(30)

(31) (32) (33) (34) (35) (36) (37) (38) (39) (40)

Page 28 of 29

Ruiz-Tijerina, D. A., and da Silva, L. G. G. V. D. (2016) Symmetry-protected coherent transport for diluted vacancies and adatoms in graphene. Physical Review B 94, 085425. Wang, C., Zhou, G., Liu, H., Wu, J., Qiu, Y., Gu, B.-L., and Duan, W. (2006) Chemical Functionalization of Carbon Nanotubes by Carboxyl Groups on Stone-Wales Defects: A Density Functional Theory Study. The Journal of Physical Chemistry B 110, 10266-10271. González-Durruthy, M., Castro, M., Nunes, S. M., Ventura-Lima, J., Alberici, L. C., Naal, Z., Atique-Sawazaki, D. T., Curti, C., Ruas, C. P., Gelesky, M. A., Roy, K., González-Díaz, H., and Monserrat, J. M. (2017) QSPR/QSAR-based Perturbation Theory approach and mechanistic electrochemical assays on carbon nanotubes with optimal properties against mitochondrial Fenton reaction experimentally induced by Fe2+-overload. Carbon 115, 312-330. Padilha, J. E., Amorim, R. G., Rocha, A. R., da Silva, A. J. R., and Fazzio, A. (2011) Energetics and stability of vacancies in carbon nanotubes. Solid State Communications 151, 482-486. Forli, S., Huey, R., Pique, M. E., Sanner, M. F., Goodsell, D. S., and Olson, A. J. (2016) Computational proteinligand docking and virtual drug screening with the AutoDock suite. Nat Protoc 11, 905-919. Shoichet, B. K. (2004) Virtual screening of chemical libraries. Nature 432, 862-865. Wang, X., Guo, J., Chen, T., Nie, H., Wang, H., Zang, J., Cui, X., and Jia, G. (2012) Multi-walled carbon nanotubes induce apoptosis via mitochondrial pathway and scavenger receptor. Toxicol In Vitro 26, 799806. Jia, G., Wang, H., Yan, L., Wang, X., Pei, R., Yan, T., Zhao, Y., and Guo, X. (2005) Cytotoxicity of carbon nanomaterials: single-wall nanotube, multi-wall nanotube, and fullerene. Environ Sci Technol 39, 13781383. Park, K. H., Chhowalla, M., Iqbal, Z., and Sesti, F. (2003) Single-walled carbon nanotubes are a new class of ion channel blockers. J Biol Chem 278, 50212-50216. González-Durruthy, M., Werhli, A. V., Cornetet, L., Machado, K. S., González-Díaz, H., Wasiliesky, W., Ruas, C. P., Gelesky, M. A., and Monserrat, J. M. (2016) Predicting the binding properties of single walled carbon nanotubes (SWCNT) with an ADP/ATP mitochondrial carrier using molecular docking, chemoinformatics, and nano-QSBR perturbation theory. RSC Advances 6, 58680-58693. González-Durruthy, M., Werhli, A. V., Seus, V., Machado, K. S., Pazos, A., Munteanu, C. R., González-Díaz, H., and Monserrat, J. M. (2017) Decrypting Strong and Weak Single-Walled Carbon Nanotubes Interactions with Mitochondrial Voltage-Dependent Anion Channels Using Molecular Docking and Perturbation Theory. Scientific Reports 7, 13271. Berman, H. M., Westbrook, J., Feng, Z., Gilliland, G., Bhat, T. N., Weissig, H., Shindyalov, I. N., and Bourne, P. E. (2000) The Protein Data Bank. Nucleic Acids Res 28, 235-242. Trott, O., and Olson, A. J. (2010) AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 31, 455-461. José, M. S., Emilio, A., Julian, D. G., Alberto, G., Javier, J., Pablo, O., and Daniel, S.-P. (2002) The SIESTA method for ab initio order- N materials simulation. Journal of Physics: Condensed Matter 14, 2745. Perdew, J. P., Burke, K., and Ernzerhof, M. (1996) Generalized Gradient Approximation Made Simple. Physical Review Letters 77, 3865-3868. Troullier, N., and Martins, J. L. (1991) Efficient pseudopotentials for plane-wave calculations. Physical Review B 43, 1993-2006. Monkhorst, H. J., and Pack, J. D. (1976) Special points for Brillouin-zone integrations. Physical Review B 13, 5188-5192. Ponder, J. W., and Case, D. A. (2003) Force fields for protein simulations. Adv Protein Chem 66, 27-85. Jiménez, J., Doerr, S., Martínez-Rosell, G., Rose, A. S., and De Fabritiis, G. (2017) DeepSite: protein-binding site predictor using 3D-convolutional neural networks. Bioinformatics 33, 3036-3042. Meslamani, J., Rognan, D., and Kellenberger, E. (2011) sc-PDB: a database for identifying variations and multiplicity of ‘druggable’ binding sites in proteins. Bioinformatics 27, 1324-1326. Feinstein, W. P., and Brylinski, M. (2015) Calculating an optimal box size for ligand docking and virtual screening against experimental and predicted binding pockets. J Cheminform 7, 18. 28


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

(41) (42) (43) (44) (45) (46) (47) (48) (49) (50) (51) (52)


da Silveira, C. H., Pires, D. E., Minardi, R. C., Ribeiro, C., Veloso, C. J., Lopes, J. C., Meira, W., Jr., Neshich, G., Ramos, C. H., Habesch, R., and Santoro, M. M. (2009) Protein cutoff scanning: A comparative analysis of cutoff dependent and cutoff free methods for prospecting contacts in proteins. Proteins 74, 727-743. Das, S., Wong, R., Rajapakse, N., Murphy, E., and Steenbergen, C. (2008) Glycogen synthase kinase 3 inhibition slows mitochondrial adenine nucleotide transport and regulates voltage-dependent anion channel phosphorylation. Circ Res 103, 983-991. Mitternacht, S., and Berezovsky, I. N. (2011) Coherent Conformational Degrees of Freedom as a Structural Basis for Allosteric Communication. PLOS Computational Biology 7, e1002301. Lewicki, T. H. (2006) Statistics: Methods and Applications: A Comprehensive Reference for Science, Industry, and Data Mining. Statsoft. Hanley, J. A., and McNeil, B. J. (1982) The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology 143, 29-36. Matthews, B. W. (1975) Comparison of the predicted and observed secondary structure of T4 phage lysozyme. Biochim Biophys Acta 405, 442-451. S Dresselhaus, M., Jorio, A., Filho, A., and Saito, R. (2010) Defect characterization in graphene and carbon nanotubes using Raman spectroscopy. Vol. 368. Xu, Z., Liang, Z., and Ding, F. (2017) Isomerization of sp2-hybridized carbon nanomaterials: structural transformation and topological defects of fullerene, carbon nanotube, and graphene. Wiley Interdisciplinary Reviews: Computational Molecular Science 7, e1283. Choudhary, O. P., Paz, A., Adelman, J. L., Colletier, J. P., Abramson, J., and Grabe, M. (2014) Structureguided simulations illuminate the mechanism of ATP transport through VDAC1. Nat Struct Mol Biol 21, 626-632. Rasola, A., and Bernardi, P. (2007) The mitochondrial permeability transition pore and its involvement in cell death and in disease pathogenesis. Apoptosis 12, 815-833. Roberto Todeschini, V. C. (2010) Molecular Descriptors for Chemoinformatics. Methods and Principles in Medicinal Chemistry I, 967. Patlewicz, G., Worth, A. P., and Ball, N. (2016) Validation of Computational Methods. Adv Exp Med Biol 856, 165-187.

29


Computational MitoTarget-Scanning Based on Topological Vacancies

Recommend Documents