Molecular Insights on the Interactions of Nitrosamines from Cigarette Smoking with CYP2A13 using Molecular Docking and Molecular Dynamics Simulation

Cigarette smoking continues to be one of the leading causes of death and disease worldwide. In the United States alone, approximately 433,000 adults die each year from smoking related diseases, causing medical expenses of up to $100 billion [1]. In Brazil, tobacco consumption is still very high, especially among young students [2]. However, worldwide statistics show that the number of smokers has fallen, mainly between the years 1980 and 2012 [3,4]. One of the consequences of smoking is the various lung diseases that its chemical components can cause, for example chronic bronchitis, chronic obstructive pulmonary disease, obstruction of small airways, emphysema and pulmonary hypertension [1].

Another serious illness that affects most smokers is cancer [5], especially lung cancer. This disease may be associated with high frequency of cigarette smoking as well as its prolonged use during life [6,7]. Increased incidences of lung cancer have been reported since the 1930s, and by 2012 it was estimated to have been responsible for the death of approximately 1,590,000 people. It remains one of leading causes for cancer deaths globally. Among men, the deaths occur mainly in Central and Eastern Europe and among women in North America, Northern Europe, Eastern Asia, Western Europe, Oceania and the Caribbean [8,9]. Overall, the mortality rate due to lung cancer is the highest in countries where people start smoking earlier, such as in North America and Europe [10].

Figure 1 shows the chemical compounds N'-Nitrosoanabasine (NAB) and N'-Nitrosoanatabine (NAT) formed from anabasine and anatabine, respectively [11]. One of the likely causes of lung cancer is the nitrosamines derived from tobacco smoke [12].

 

Figure 1: 2D structural formulas of nitrosamines NAB and NAT.

The molecules present in tobacco and tobacco smoke begin to act in the body only after enzymatic activation that can be achieved by cytochrome P450 2A13 (P450 or CYP) and is considered a key reaction pathway for initiation of bioactivation of carcinogenesis, and consequently of the progressive changes in the biological profile of the cell [13,14].

A useful tool to demonstrate interactions between nitrosamine linkers and CYP2A catalytic sites is the molecular docking study, which comprises the mechanism of activation and/or oxidation of substrates that may aid in the development of new drugs. In this sense, the present study aims to demonstrate the mechanism of interaction for the activation of NAB and NAT molecules by cytochrome P450 2A13 [15,16].

In this paper, we have described the molecular docking results of NAB and NAT ligands interaction with cytochrome P450 2A13. We also performed Molecular Dynamics (MD) simulation to evaluate the intermolecular interactions during the period of the simulation. Finally, we studied the spontaneity of the formation of the systems and their energetic components by means of free energy calculations using the Molecular Mechanics/Generalized-Born Surface Area (MM/GBSA) methods.

The chemical structures of the NAB and NAT were designed using the GaussView software and then optimized with the DFT method at the B3LYP [17,18] using the Gaussian 09 software [19].

Molecular docking studies were performed using the software Molegro Virtual Docker 5.5 (MVD 5.5) [20]. The crystallographic structure of cytochrome P450 2A13, used as a receptor, can be located in the Protein Data Bank (PDB) with ID: 4EJG [16]. The MolDock scoring function (GRID) was used with a grid resolution of 0.30 Å and a radius of 6 Å encompassing the entire bond cavity.

Every inhibitor atom had its atomic charge calculated using the Restrained Electrostatic Potential (RESP) from the HF/6-31G* basis set [21-23]. The protonation state of each protein residue was determined using the PDB2PQR server [24].

The protein inhibitor system was constructed for the simulation using the tleap modulus of AMBER 16 [25,26]. The General Amber Force Field (GAFF) [27] and ff14SB [28] force field was used for all the solvent-solvent simulations with the TIP3P [29] model water molecules. The system was solvated in a truncated octahedron periodic box. The cut-off radius was 12 Å for all solvent atoms relative to the solute in the x, y and z axes and a suitable number of counter-ions were added to the system to neutralize their charge.

Energy minimization was performed in the AMBER 16 sander module. At each stage, 5000 cycles were performed using the steepest descent method and conjugate gradient algorithm. In the first, the protein and inhibitor were fixed by applying a harmonic force constant of 200 kcalmol-1Å-2, the water molecules being free to move. In the following two, the restraint strength was reduced, so that in the final stage the entire system would be able to move freely.

After minimization, the systems were gradually heated from 0 to 300k in five steps, for a total time of 600ps. In four of these, a constant force harmonic of 25kcalmol-1Å-2 was used to restrict heavy atoms. In the final step, the restraining force was reduced to zero so that the whole solute was free to move. These stages were performed using the Langevin thermostat [30] with a collision frequency of 3.0 ps-1.

Following this, in order to obtain system equilibrium, a simulation of 2ns was performed at a temperature of 300k, with constant pressure without restriction. The calculation of the electrostatic interactions was performed using the Particle Mesh Ewald method [31] and the SHAKE [32] algorithm was utilized to restrict the bond lengths involving the hydrogen atoms. Afterward, molecular dynamics simulations of production up to 100ns were performed for all systems, unrestrained in the NVT ensemble.

The binding free energies were calculated using the MM/GBSA [33,34] method implemented in the AMBER 16. For each system, 500 snapshots of the last 5ns were used on the MD trajectory. For each snapshot, free energy is calculated for the protein, ligand, and complex using a single trajectory approach. The binding free energy was computed as the difference:
ΔG_bind = G_complex - G_receptor - G_ligand (1)
ΔG_bind = ΔH - TΔS = E_MM + ΔG_solv – TΔS (2)
ΔE_MM = ΔE_internal + ΔE_{electrostatic} + ΔE_vdW (3)
ΔG_solv = ΔG_GB + ΔG_nonpol (4)

In which ΔG_bind is the free energy of the protein-ligand binding, resulting from the sum of the molecular mechanics energy (ΔE_MM), the desolvation free energy (ΔG_solv), and the entropy (-TΔS). The energy of molecular gas phase mechanics (ΔE_MM) can be described by the sum of the internal energy contributions (ΔE_internal), the sum of the connection, angle and dihedral energies), electrostatic contributions (ΔEeletrostatic) and van der Waals terms (ΔE_vdW). Desolvation free energy (ΔG_solv) is the sum of the polar (ΔG_GB) and non-polar (ΔG_nonpol) contributions. The polar desolvation term was calculated using the implicit Generalized Born (GB) approaches.

For studying the interactions between NAB and NAT molecules with CYP450 2A13, the molecular docking simulation at the cytochrome binding site was performed. This site has a surface area of 166.40 Å2, volume 67.58 Å3 and position x: 22.55, y: -11.51 and z: 37.43 for accommodation of the molecules (Figure 2).

 

Figure 2: The binding pocket of the molecules is represented by green color.

This cytochrome was chosen because it is expressed in the respiratory tract [35] and is responsible for the metabolism of nitrosamines [36].

The molecular docking simulation of these nitrosamines with cytochrome was performed because NAB [37] and NAT [38] are carcinogenic agents present in tobacco smoke. In several studies they have been reported as being responsible for causing lung cancer [39-41]. Therefore, elucidation of the binding mechanism of these molecules is important for the study of their metabolism. The molecular docking results are summarized in table 1.

 

Molecules	MolDock Scorea	Rerank Scoreb	Interaction (kcal/mol)c	Contribution of the Cofactord	Internale	LE1f	LE3g
NAB	-86.74	-72.07	-97.96	-17.28	11.22	-6.19	-5.14
NAT	-82.86	-71.82	-99.8	-11.03	16.93	-5.91	-5.12

Table 1: Molecular docking results.

Note: ^aThe scoring function used by MolDock is derived from the PLP score and has a new term for hydrogen bonds and new charge methods [20]^bRerank Score is a linear combination of E-inter (stereo, van der Waals, hydrogen bonding and electrostatic) between the binder and the protein, and E-intra (torsion, sp2-sp2, hydrogen bonding, van der Waals and electrostatic) of the binder corrected by pre-defined coefficients [20]^cEnergy of the total interaction between the pose and the protein^dEnergy of interaction between the pose and the prosthetic group^eEnergy of the pose^fMolDock Score divided by the count of the heavy atoms^gRerank score divided by the count of the heavy atoms

The MolDock score found for NAB (-86.74) and NAT (-82.86) at the P450 binding site was favorable for the ligand receptor interaction, as it favored the molecular stability of the system. In addition, the affinity values for the binding of the molecules to the cytochrome were favorable, demonstrating that these compounds were capable of binding to cytochrome to perform productive interactions. For NAB and NAT, these values were -97,963 kcal/mol and -99,803 kcal/mol respectively. The energy contribution of the HEME prosthetic group was also calculated and its values were found favorable in both the systems. The energy contribution for the molecules NAB and NAT were -17.281 kcal/mol and -11.036 kcal/mol respectively.

With the above obtained results it became possible to verify the interactions between the molecules and the residues of the protein and the HEME prosthetic group.

These interactions are described in the following tables 2 and 3 and figure 3a and 3b.

 

Molecule	Residue	Group Heme	Interaction	Distance (Å)
P30	Ala301-CB	-	Hidrophobic: Alkyl	4.18
P31	Ala301-CB	-	Hidrophobic: Pi-Alkyl	5.32
-	Ala301-CB	B	Hidrophobic: Pi-Alkyl	4.05
-	Ala301-CB	A	Hidrophobic: Pi-Alkyl	4.69
P30	-	A	Hidrophobic: Pi-Alkyl	4.42
P30	-	C	Hidrophobic: Pi-Alkyl	5.36

Table 2: Interactions performed by NAB with CYP2A13 and with the HEME group.

 

Molecule	Residue	Group Heme	Interaction	Distance (Å)
O1	Ala301-CA	-	Hydrogen Bond	3.47
P20	Ala301-CB	-	Hydrophobic: Alkyl	3.78
P21	Ala301-CB	-	Hydrophobic: Pi-Alkyl	5.21
-	Ala301-CB	A	Hydrophobic: Pi-Alkyl	4.69
-	Ala301-CB	B	Hydrophobic: Pi-Alkyl	4.05

Table 3: Interactions performed by the NAT with CYP2A13 and with the HEME group.

 

Figure 3: Molecular interactions. (a) Molecular interactions between NAB and CYP2A13 and HEME group. (b) Molecular interactions between NAT and CYP2A13 and HEME group.

The value of binding affinity energy for the NAB compound was found to be -97,963 kcal/mol. The Pyridine ring (P31) interacted with the CB atom of alanine 301, at a distance of 5.32 Å through a hydrophobic pi-alkyl interaction. The same alanine atom interacted with the Piperidine ring (P30) at a distance of 4.18 Å by means of an alkyl interaction. The piperidine ring interacted with the pyrrole rings A and C of the HEME prosthetic group via hydrophobic pi-alkyl interactions at the distances of 4.42 and 5.36 Å, respectively. The alanine residue also interacted with the HEME group via hydrophobic pi-alkyl interactions with its pyrrole rings A and B at the distances of 4.69 and 4.05 Å, respectively.

The NAT molecule, having an interaction energy value of -99.803 kcal/mol, showed the best interaction with the active site of the protein in the docking simulation. This molecule showed three interactions with the alanine residue. The atom O1 of the molecule interacted with the CA atom of alanine 301 by means of hydrogen bonding at a distance of 3.47 Å. In this interaction the NAT molecule is the acceptor and the alanine molecule is the donor of the hydrogen atom. The CB atom interacted with the P21 ring by means of a pi-alkyl interaction at a distance of 5.21 Å while it interacted with the Pyridine ring (P20) through an alkyl interaction at a distance of 3.78 Å. The pyrrole rings A and B of the HEME group interacted with the alanine residue by means of hydrophobic pi-alkyl type interaction at the distances of 4.69 and 4.05 Å, respectively.

In molecular docking simulations, with the Mocker Virtual Docker, the protein is considered rigid and the binder flexible [20]. Thus, the molecular dynamics simulation was adopted to refine the results obtained for the cytochrome interaction with the NAB and NAT ligands. In addition, the stability of the systems was also analyzed. For this analysis the Root Mean Square Deviation (RMSD) value was used in relation to the starting structure (minimized, heated and balanced) for molecular dynamics simulation. The RMSD results for the systems over a total time of 100 nanoseconds are shown in figure 4.

 

Figure 4: RMSD of the CYP2A13-ligands systems. The RMSD of the protein backbone is represented in black color while the RMSD of the binder is represented in various other colors.

The mean values obtained for the RMSD of the alpha carbon of the protein backbone were as follows, CYP-NAB: 1.03 Å and CYP-NAT: 1.14 Å. These values demonstrated that the protein structure during the molecular dynamics simulation did not deviate significantly from its initial conformation.

Throughout the simulation of the NAT its structural conformation underwent changes alternating its state between moments of molecular balance and small conformational variations. The NAB molecule also maintained its equilibrium at the CYP site, undergoing small conformational changes during the simulation. Throughout the simulation time the NAB and NAT compounds remained bound to the CYP binding pocket.

To obtain the values of the binding free energy components of all the ligands interacting with CYP, we used the MM/GBSA methods [33-35]. The results obtained are summarized in table 4.mation underwent changes alternating its state between moments of molecular balance and small conformational variations. The NAB molecule also maintained its equilibrium at the CYP site, undergoing small conformational changes during the simulation. Throughout the simulation time the NAB and NAT compounds remained bound to the CYP binding pocket.

 

Ligand	ΔEvdw	ΔEele	ΔGele,sol	ΔGnpol,sol	ΔEnon-polar	ΔEpolar	ΔGbind
NAB	-34.97 ± 0.04	-9.50 ± 0.05	10.31 ± 0.02	-4.02 ± 0.01	-38.99	1.07	-38.19 ± 0.05
NAT	-31.78 ± 0.04	-7.92 ± 0.05	12.62 ± 0.03	-4.01 ± 0.01	-35.79	4.70	-31.08 ± 0.05

Table 4: Binding free energy and its components.

According to the values obtained with the free energy calculations, the binders showed the following bond affinity values: CYP-NAB: -38.19 and CYP-NAT: -31.08.

In addition to the absolute values of free energy, its components were also calculated so that the mechanism of interaction of the molecules with the protein is better understood. In all the systems the van de Waals interactions (ΔE_vdW) and the electrostatic interactions of the gas phase (ΔE_ele) favored the bond free energy (ΔG_bind) of the complexes. The values for nonpolar solvation energies (ΔG_npol,sol) and the non-polar (ΔE_non-polar) contributions were similar to each other and also favored the ligand receptor interaction.

Our molecular docking results demonstrated that the major interactions that occurred between the NAB and NAT molecules and the enzyme and the HEME group were hydrophobic in nature. During the entire period of simulation, the ligands continued to interact with the CYP, undergoing some conformational changes. Our MM/GBSA results demonstrated that the formation of the complexes was spontaneous and the values obtained for ΔG_bind were as follows: CYP-NAB: -38.19 kcal/mol and NAT: -31.08 kcal/mol. For the complexes the van de Waals (ΔE_vdW) interactions and the electrostatic interactions of the gas phase (ΔEele) favored the bonding energy (ΔG_bind). The values for nonpolar solvation energies (ΔGnpol,sol) and the non-polar (ΔEnon-polar) contributions were similar to each other and also favored the ligand-receptor interaction. With our results it was possible to understand the mechanism of interaction of NAB and NAT that after being metabolized by CYP2A13 will form DNA adducts capable of inducing lung cancer.

Jorddy N. Cruz and Antonio M. J. C. Neto appreciate the support of the Federal University of Pará and the National Council for Scientific and Technological Development (CNPq). Oliveira M. S (Process Number: 1662230), Costa W. A. (Process Number: 1427204) and Bezerra F. W. F. (Process Number: 1718822) thank CAPES for the doctorate scholarship.

Conflict of Interest: The authors declare that there are no conflicts of interest regarding the publication of this paper.

Ethical Approval: The present work was approved by the Ethics Council of the Federal University of Pará

J.N. Cruz is the principle author. The study includes no identifiable information.

All authors contributed equally to the development of this article.