Next Article in Journal
Synthesis, Characterization and In Vitro Evaluation of Hybrid Monomeric Peptides Suited for Multimodal Imaging by PET/OI: Extending the Concept of Charge—Cell Binding Correlation
Previous Article in Journal
The Specificity and Broad Multitarget Properties of Ligands for the Free Fatty Acid Receptors FFA3/GPR41 and FFA2/GPR43 and the Related Hydroxycarboxylic Acid Receptor HCA2/GPR109A
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Pharmaceuticals
Volume 14
Issue 10
10.3390/ph14100988
pharmaceuticals-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Prerequisite Binding Modes Determine the Dynamics of Action of Covalent Agonists of Ion Channel TRPA1

by
Balázs Zoltán Zsidó
1,
Rita Börzsei
2,
Erika Pintér
1 and
Csaba Hetényi
1,*
1
Department of Pharmacology and Pharmacotherapy, Medical School, University of Pécs, Szigeti út 12, 7624 Pécs, Hungary
2
Department of Pharmacology, Faculty of Pharmacy, University of Pécs, Szigeti út 12, 7624 Pécs, Hungary
*
Author to whom correspondence should be addressed.
Pharmaceuticals 2021, 14(10), 988; https://0-doi-org.brum.beds.ac.uk/10.3390/ph14100988
Submission received: 23 August 2021 / Revised: 22 September 2021 / Accepted: 24 September 2021 / Published: 28 September 2021
(This article belongs to the Section Pharmacology)
Browse Figures
Review Reports Versions Notes

Abstract

:
Transient receptor potential ankyrin 1 (TRPA1) is a transmembrane protein channeling the influx of calcium ions. As a polymodal nocisensor, TRPA1 can be activated by thermal, mechanical stimuli and a wide range of chemically damaging molecules including small volatile environmental toxicants and endogenous algogenic lipids. After activation by such compounds, the ion channel opens up, its central pore widens allowing calcium influx into the cytosol inducing signal transduction pathways. Afterwards, the calcium influx desensitizes irritant evoked responses and results in an inactive state of the ion channel. Recent experimental determination of structures of apo and holo forms of TRPA1 opened the way towards the design of new agonists, which can activate the ion channel. The present study is aimed at the elucidation of binding dynamics of agonists using experimental structures of TRPA1-agonist complexes at the atomic level applying molecular docking and dynamics methods accounting for covalent and non-covalent interactions. Following a test of docking methods focused on the final, holo structures, prerequisite binding modes were detected involving the apo forms. It was shown how reversible interactions with prerequisite binding sites contribute to structural changes of TRPA1 leading to covalent bonding of agonists. The proposed dynamics of action allowed a mechanism-based forecast of new, druggable binding sites of potent agonists.

1. Introduction

Mammalian neurons of the pain pathway detect potentially dangerous environmental signals. In the peripheral nervous system, there are specialized nociceptive neurons, that recognize either noxious chemical signals, thermal or mechanical stimuli. Pathological processes, such as tissue damage and inflammation elicit the formation and subsequent release of a wide variety of mediators (arachidonic acid derivatives, free radicals, H2O2, H2S, etc.). These molecules depolarize the nerve terminals of nociceptors, which transmit the signals to the central nervous system. The transient receptor potential ankyrin 1 (TRPA1) is a Ca2+-permeable cation channel that was identified as the chemical nocisensor, expressed by primary afferent nerve fibers [1,2,3,4,5,6,7,8,9]. Activated TRPA1 promotes pain itching and induces local neurogenic inflammatory response via the release of neuropeptides, such as substance P, calcitonin gene-related peptide and neurokinins.
TRPA1 receptor is activated by the binding of electrophile ligands (Figure 1) to its N-terminus cytoplasmic binding site (Figure 2A), which is characterized by three nucleophilic cysteine residues (C621, C641 and C665) [6,7]. This binding event induces a local conformational change, that is translated to the whole of the receptor, and a 15° rotation of the transmembrane domain is observed [7], resulting in a pore widening, that facilitates Ca2+ influx, which first potentiates, then desensitizes agonist-induced responses [7], resulting in an inactive state of the TRPA1 ion channel [10]. On the local scale, a cytoplasmic A-loop near the transmembrane region of the receptor overlays the binding site cavity [7] and initially sterically hinders agonist binding. Upon the binding of the electrophile agent, a flip of the A-loop (residues 666–680, Figure 2B) was observed towards the cell membrane, leaving more binding space for the agonists. Thus, flipping of A-loop contributes to the widening of the pore of the ion channel, and the above-mentioned activation process [7] at the same time, and therefore, it is important in agonist design.
Known electrophile agonists of the TRPA1 receptor include dimethyl trisulfide (DMTS) and allyl isothiocyanate (AITC) [4,5,8]. However, given the size of these molecules, they would not bind selectively to only one cysteine amino acid residue in the body. Thus, the need for a selective site-specific electrophile agonist (JT010) was first met in 2015 by Takaya et al. [4], and a potent thiazol derivate agonist (EC50 = 0.65 nM, Figure 1) with a covalent alkyl-halide warhead was designed. Since then, the binding position of JT010 was experimentally found by cryo-electron microscopy [7]. Recent structural studies [2,3,6,7] provided additional details of the agonist binding mechanism and consequential receptor activation. Besides JT010, the binding of the other two site-specific covalent agonists BITC [6] and bodipy-iodoacetamide [7] (Figure 1) was investigated, which preferentially bind to the active site cysteine C621 of the TRPA1 receptor.
Covalently binding agonists are subjects of intense research [11], and they often form covalent bonds with nucleophilic cysteine residues. As cysteine is abundant in the human proteome, a careful design has to be performed to achieve binding site specificity [11,12] to avoid unwanted promiscuity of the agonists via non-selective covalent bonds with non-targeted cysteines. Thus, a common strategy of covalent agonist design adopts known agonists having selective non-covalent interactions [11] with the receptor. Both covalent and non-covalent interactions have been fully described [6,7] at the final, irreversible binding mode of the agonists in Figure 1. However, the binding routes leading to the final binding pocket have not been mapped at atomic resolution. The above-mentioned, agonist-induced structural changes of TRPA1 at the A-loop (Figure 2B) suggest that the agonist may form dynamic interactions with prerequisite sites on their route to the final, covalent positions.
The present study is focused on the elucidation of binding dynamics of covalent agonists using experimental structures of their complexes with TRPA1 (holo form) and compared to prerequisite interactions with the apo form. Covalent and non-covalent molecular docking techniques are tested at forecasting prerequisite and final, covalent binding modes of the agonists. The docked agonist-TRPA1 complexes are subjected to molecular dynamics calculations to complete the binding mechanism at the atomic level. We also aim at the mechanism-based forecasting of prerequisite binding sites which can become druggable targets of agonists in drug design projects.

2. Results and Discussion

2.1. Final Covalent Binding Modes

The atomic resolution structures of three covalent agonists JT010 [4], BITC [6] and bodipy-iodoacetamide [7] (Figure 1) bound to the TRPA1 ion channel are available in the Protein Databank (PDB). JT010 and bodipy-iodoacetamide have an alkyl-halide covalent warhead, the bonds made with chlorine and iodine atoms break up during covalent binding to C621 (Figure 3). BITC participates in an isothiocyanate bond with C621 (Figure 3), the N=C double bond diminishes, and the C atom forms a covalent bond with the S atom of the target C621, and the N atom of BITC gains a H atom. As the nucleophilic C621 binds all three electrophilic agonists covalently [6,7], our docking studies focused on the surrounding binding pocket.
As a first step of our investigation, the popular program package FITTED [13,14,15] was tested using the experimental PDB structures as references for comparison with the docking results. FITTED was first tested to reproduce the final binding modes of the covalently bound agonists. A standard evaluation protocol was applied to all covalent docking results. Firstly, the structural match of the calculated binding mode (bind position, orientation, and conformation of the ligand) to the crystallographic reference was calculated and the best match was expressed as a root mean squared deviation (RMSDbest, see Section 4 for the definition of all metrics used in the Tables) [16]. Secondly, it was tested if the docking method identified RMSDbest as an energetically favorable binding mode and ranked it at the top of the list of all binding modes. The second criterion reports on the applicability of the binding free energy (scoring) function of the docking method.
The covalent docking of the agonists to the holo form of TRPA1 was structurally successful as the RMSDbest values were comparable/below 2.5 Å (Table 1), a threshold accepted in the literature [17,18,19,20,21]. Only bodipy-iodoacetamide showed a slightly elevated RMSD (Table 1), which is due to the mobility of the -SH group around the Cβ of C621 during docking. The rotation of the S atom around the Cβ of C621 also turns bodipy-iodoacetamide from its experimental position around its longitudinal axis by ca. 180°. The Cβ-S-Cbodipy-iodoacetamide angle is also smaller than that observed in the experimental structure. The binding modes with RMSDbest were positioned to the first place on the ranking lists for all three agonists (Table 1). The calculated free energy of binding of JT010 is the most favorable, closely followed by that of BITC (Table 1), apparently, the alkyl-halide covalent bond formed by the alkyl-iodine pharmacophore of bodipy-iodoacetamide is less favorable compared to that formed by the alkyl-chlorine pharmacophore of JT010. The efficiency index (EINHA) of BITC is the best out of the three reference ligands.
The large conformational flexibility of target molecules is a challenge for fast docking programs and can be handled by the involvement of molecular dynamics approaches [22] which require longer calculation times. Conformational flexibility is important when induced fit occurs during agonist binding. Here, the A-loop is a flexible element that covers the binding site (Introduction, Figure 2B) in the apo TRPA1 conformation, and sterically prevents the agonist from reaching its destination (holo position) at the bottom of the pocket resulting in unacceptably large RMSD values and positive energies (Tables S1–S3). Therefore, the A-loop was removed from TRPA1 and covalent docking calculations were repeated using FITTED. Although the A-loop consists of the amino acids from 660–680, in the apo docking calculations, only the amino acids 665–677 were removed, as these are the ones that elicit the greatest movement during apo to holo transition. In the absence (Table 1, Table 2, Table 3) of the A-loop, all covalent bonds were formed with the apo TRPA1 (Table 1), with RMSD values larger than those observed in the case of the holo TRPA1 (Table 1). The corresponding ΔGbest values were lower than those of the apo TRPA1 by 6% on average. These findings emphasize the role of A-loop-agonist interactions in the final binding position. However, the removal of the A-loop did not influence the EI and ΔGbest order of the ligands.
The formation of the covalent bond between the agonist and TRPA1 (Figure 3) is a quantum mechanical phenomenon, which is hard to treat adequately by docking programs based on molecular mechanics scoring functions [23,24]. Despite the above challenges, the covalent docking methodology of FITTED performed well for the above test cases and supplied relevant structural and scoring (ranking) results. Encouraged by the above test results, we expect that FITTED will also help in mapping the prerequisite binding modes on route to the final binding pocket.

2.2. Prerequisite Binding Modes

As it was discussed in Section 1, the entrance to the binding cavity (outer prerequisite binding mode) and the formation of the final ligand-target covalent bond (inner prerequisite binding mode) is hindered by the position of the A-loop (Figure 2B) in the apo form of the TRPA1 target. During a successful binding process, the ligand initiates the flipping of the A-loop via intermolecular interactions with the loop. To develop such interactions, the ligand needs to occupy a prerequisite binding mode outside the final binding pocket. Two different programs, FITTED and AutoDock 4.2.6 (The Scripps Research Institute, La Jolla, CA, USA) [25], were involved in the mapping of possible prerequisite binding modes by non-covalent docking and the results are shown in Table 2 and Table 3, respectively. Both programs are based on physico-chemical principles. FITTED is a genetic algorithm-based docking method, that includes an ESFF [26] force field-based search engine, called CDiscoVer [27] to perform conjugate gradient minimizations [13]. AutoDock also uses a (Lamarckian) genetic algorithm and AMBER-based intermolecular force field terms for scoring [25].
For prerequisite binding modes RMSD was not calculated, the distance (d) between TRPA1 C621 S atom and the atom of the agonist that participates in the covalent binding was used as a measure of ligand position instead. The dbest value indicates the closest distance between the aforementioned atoms (black dots in Figure 3) achieved by subsequent docking calculations. The match of the docked binding mode was expressed as the percentage of matching amino acids (AAmatch) compared with the experimental binding pocket. Both programs ranked the binding mode of BITC with the dbest as the top 1st in all prerequisite docking calculations (Table 2 and Table 3). However, BITC is considerably smaller, than JT010 and bodipy-iodoacetamide. The head-to-tail docking orientation of larger agonists, like JT010 and bodipy-iodoacetamide might cause elevated dbest values.
In all cases and both scorings, the holo prerequisite docking calculations yielded better AAmatch and ΔGFD and ΔGAD, than the apo prerequisite docking calculations (Table 2 and Table 3). This was expected, as the holo conformation of the binding site is already prepared to accept the agonists. The FITTED prerequisite holo docking calculations yielded an AAmatch of 100% in the cases of all three agonists. The dbest values of the prerequisite docking calculations were under 4.0 Å in all cases, with the only exception of bodipy-iodoacetamide holo docking (Table 3). In the case of AutoDock, the dbest values of the holo prerequisite calculations were below 7.0 Å, and slightly above it in the apo prerequisite docking runs. The dbest ≥ 7.0 Å values are due to head-to-tail binding mode of the agonists (Table 2 and Table 3).
The comparison of the prerequisite binding modes produced by FITTED on both the holo and apo docking calculations of the three compounds resulted in C621 as the only common binding site amino acid for all three agonists. All holo prerequisite docking calculations found C621 with only the exception of bodipy-iodoacetamide prerequisite holo docking with AutoDock. These findings suggest that different prerequisite binding modes (Tables S5–S7) might result in good final covalently binding positions (Table 1). By investigating more agonists with two docking programs, one might expect to discover common amino acids that indicate a larger prerequisite binding area. If an agonist at least partially interacts with the prerequisite binding area it has a chance to find its way to the final binding pocket.
Amino acids C665, P666 and F669 of the A-loop are part of the binding pocket of bodipy-iodoacetamide, and also of most prerequisite binding modes of bodipy-iodoacetamide and JT010 (Tables S5–S7) found by both programs. C665 is also highlighted in the literature [7] as an important amino acid both in agonist binding and receptor activation. The absence of interaction of BITC with the above-mentioned amino acids might be due to the smaller size of BITC compared to the other two agonists (Figure 1). These results suggest a previously unexplored structural role of the amino acids P666 and F669 co-operating with C665 in agonist binding and consequent flipping of the A-loop, leading to conformational changes and receptor activation. These findings were also strengthened by virtual mutation and docking (Figure S1). However, BITC interacts with another part of the A-loop as prerequisite binding site found by apo docking with AutoDock. This BITC site was also sufficient to open the binding pocket in molecular dynamics (MD) simulations (Section 2.3).
Although not an accepted medicine, JT010 has a remarkable EC50 of 0.65 nM [4]. Thus, in novel drug design, JT010 can be regarded as a reference point, and therefore, it was further investigated from the mechanism viewpoint of this Section. During the transition from a non-covalent, prerequisite binding mode (d = 3.6 Å) to the final, covalent binding mode of JT010, its interactions with F669 and Y680 diminish (Table S4), and new interactions with K620, I623 and E625 are formed (cut-off distance of interaction of 3.5 Å for heavy atom–heavy atom distance). Interactions with the two binding site cysteines, C621 and C665 [7] are observed for both non-covalent and covalent binding modes of JT010. As it was highlighted in the previous section, F669 is part of the A-loop and has a possible role in the flipping of the loop during agonist binding to the TRPA1 receptor, based on the example of JT010. It can be hypothesized, that interaction with F669 is only important in the initial prerequisite binding of the agonist, later during covalent bond formation this interaction diminishes, and the agonist penetrates deeper into the binding site, interacting with amino acid residues that are in close proximity of C621. This observation is strengthened by MD simulations also (see Section 2.3). The covalent ΔGFD of JT010 almost doubles (and consequently its EI also), compared to that of the prerequisite binding mode.
Regardless of their potency, all three agonists can activate the TRPA1 ion channel, however, to prevent xenobiotic overload of the body it is advisable to administer the lowest possible dose of a drug, which is only effective if the agent is highly potent. If using JT010 as a reference for future studies, the following limit values can be concluded for the selection of potent agonists. A prerequisite EI value of at least 2 kcal/mol, and a prerequisite dbest ≤ 4.0 Å, and a prerequisite ΔGFD of at least -35 kcal/mol forecast a strong agonist. Notably, the ΔGFD value of JT010 is approximated by that of BITC, and the EI of BITC even surpasses that of JT010 (Table 1). bodipy-iodoacetamide somewhat lags behind JT010 and BITC. These findings are in good agreement with the literature, as the EC50 value of iodoacetamide (without the bodipy label) is 357 µM [28], which is substantially larger, than that of JT010. The EC50 of allyl-isothiocyanate (a similar compound to BITC) is 37 nM [7], which is also in the nanomolar range, as the EC50 of JT010.
The docking performance of FITTED slightly outperformed AutoDock as seen in Table 2 and Table 3. However, FITTED requires a probe previously placed within the binding site to select the binding site amino acids, which obviously helps the search. At the same time, AutoDock did not require such information, and an unrestricted search could be performed for the prerequisite binding mode within the docking box. Thus, we decided to use the prerequisite binding mode of BITC found by AutoDock apo calculation with A-loop (Table S3) for further MD simulations in the next, Section 3.

2.3. Ligand Migration Dynamics Connecting Prerequisite and Final Binding Modes

MD simulations (100 ns, unrestrained, with explicit waters and simulated annealing protocol as described in Section 4) were performed on both apo and holo forms of TRPA1 (Table 4) to further explore the binding dynamics of the agonist BITC of the best EI value (Table 1). As the results of the previous Section indicated that the prerequisite binding modes affect A-loop, we were particularly interested in the structural changes of the loop, and the communication between the distinct prerequisite and final binding modes. An MDapo simulation was used as a reference, to observe if there are any changes in the conformation of the A-loop in the absence of the agonist. Then, two MD simulations were started from two prerequisite binding modes of BITC on the apo TRPA1 found in the previous Section, one of them interacting with the loop (MDrank1). Finally, an MD simulation was started from the experimental binding position of BITC (MDholo), with the covalent bond cut and the geometry of the N=C=S bonds restored. In the MDholo simulation, an unbinding-binding event occurred and the A-loop remained stable throughout the simulation (Figure 4). The interaction with the original five amino acids of the TRPA1 pocket gradually diminished and appeared once again in a very short time interval of the first 0.7 ns (Figure 4). During the entire MDapo (Table 4) simulation, no significant changes were observed in the conformation of the A-loop, while, in the case of MDrank1, the loop moved upward (the red and blue arrows show the movement of A-loop and the teal arrows the movement of BITC on Figure 5). In the starting position of MDrank1, BITC interacted with the loop and was positioned beneath it (marked with 0 ns at Figure 5). After a very short time (0.3 ns) the ligand dissociates from the TRPA1 surface, dragging down the loop with itself. After 17.6 ns the loop moved upwards, approximating the open position of the binding site which is present in the holo structure (Figure 2 and Figure 5). Finally, (at 38 ns) BITC finds its way back into the binding pocket and resides there for 2 ns until dissociation. At the same time, in the case of MDrank3 (Table 4) in which the docked position of BITC did not interact with the A-loop, BITC dissociated after 1.6 ns and afterwards, no changes were observed in the structure of the loop. Thus, the above results showed how the binding of an agonist to the A-loop induces its motion towards opening the binding pocket and allowing the entrance of the agonists.
The interaction of the prerequisite binding mode with the P669 amino acid side chain (first mentioned in the previous Section) was observed in the starting frame of MDrank1 and was diminished both upon dissociation and the penetration of BITC towards its final binding pocket. BITC in the prerequisite binding mode forms mainly polar interactions with amino acids of the A-loop and the other loop (Figure 5), such as S613, P674, T675 and Q676. Towards the final binding mode, however, BITC interacts with hydrophobic amino acids such as V678, I679 and Y680. This latter observation is strengthened by the MDholo run (Figure 4), where interactions with I623 and Y662 dominate. The distances between the N atom of N615 of the A-loop and the O atom of Q676 of the opposite loop (Figure 5) appear to be good indicators for A-loop opening and closing.

3. Materials and Methods

3.1. Preparation of the Ligand Structures

Ligand conformations were obtained from their respective atomic coordinate structure files of TRPA1 receptor-ligand complexes. The ligands were modified to regain their original structure, before the covalent interaction. To JT010 [4] a chlorine was added. In the case of benzylisothiocyanate (BITC, [6]) the proper bond orders and hybridization states were restored, as R-N=C=S, the hybridization states of N and S were set to sp2 and that of C to sp. Finally, in the case of bodipy-iodoacetamide [7] an iodine was added. These modifications were carried out using the builder function of PyMol (Schrödinger, New York, NY, USA) [29]. All the ligands were energy minimized by a quantum chemistry program package, MOPAC [30,31] with PM7 parametrization [31]. Hydrogens and Gasteiger–Marsili [32] partial charges were added by OpenBabel [33]. In the case of FITTED program package [13,14,15,34], the built in preparation steps were used with default settings. For BITC the molecular mechanics force field parameters were obtained from the general AMBER force field (GAFF) [35]. The ligand was built in Maestro [36], then semi-empirical quantum mechanics optimization was performed with MOPAC [30,31] using PM7 parametrization [31], with the gradient norm set to 0.001. After energy minimization, a further force calculation step was included, the force constant matrices were positive definite. The RED-vIII.52 [37] software was used for restrained electrostatic potential (RESP) charge calculations, using RESP-A1B fitting (compatible with the AMBER99SB-ILDN force field) after ab inito geometry optimization by GAMESS [38]. Acpype [39] was used to assign atomtypes, bond and angle parameters for topology of ligand. The missing bond stretching, angle bending and torsional parameters were calculated by the antechamber [39,40] and parmchk utilities of AmberTools program package similarly as described in [41]. Torsional parameters for R-N=C=S moiety were manually added.

3.2. Target Preparation

The atomic coordinate file of the ligand free TRPA1 receptor was obtained from the Protein Data Bank (PDB, [42]), under the accession code 6V9W [7]. As the four chains of the target are symmetrical (homotetramer), only one chain was used to reduce computational costs. The amino acids of a chain do not interfere with the binding of the ligand to another chain. The missing atoms and residues [43] were rebuilt using SWISS MODEL [44], and energy minimized with GROMACS [45]. The convergence threshold of the steepest descent optimization was set to 103 kJ mol−1 nm−1, and that of the conjugate gradient optimization to 10 kJ mol−1 nm−1. AMBER99SB-ILDN force field [35] was used for the calculation, and a position restraint at a force constant of 103 kJ mol−1 nm−2 was applied on heavy atoms. The targets were further optimized by ProCESS tool of FITTED, with the original settings [13]. In the case of AutoDock (The Scripps Research Institute, La Jolla, CA, USA), the added H atoms and partial charges were kept from energy minimization.

3.3. Covalent Docking with FITTED

Covalent docking calculations were carried out using FITTED [13,14,15,34]. The covalent residue (C621) and adjacent basic residue (P622) were adjusted in the graphical user interface of the program. Root mean squared deviation (RMSD) values were calculated between the crystallographic and representative ligand conformations, if available. All other settings were used as the default of the program. In the PREPARE step of the program, the binding site interacting amino acids were identified by leaving the crystallographic ligand in the structure, which was then removed after this step. The non-covalent docking was performed similarly, with the exception, that the covalent mode of the program was switched off.

3.4. Prerequisite Docking with AutoDock 4.2

Prerequisite binding calculations were performed by AutoDock [25,46,47,48,49]. The number of grid points was set to 60 × 60 × 60 with a 0.375 Å grid spacing. Lamarckian genetic algorithm was used, flexibility on all active torsions was allowed on the ligands. Ten docking runs were performed for all ligands, and the resulting ligand conformations were ranked based on their calculated free energy of binding values. The binding mode with the most favorable calculated energy of binding was ranked in the lowest rank.

3.5. Molecular Dynamics Simulations

The apo TRPA1 and dry docked complexes of BITC were subject to a two step energy minimization, including steepest descent and conjugate gradient algorithms as described in “Target preparation”. After energy minimization, the apo and dry docked complexes were subject to 100-ns-long MD simulations. The simulation box was filled up with explicit TIP3P [50] water molecules, and to neutralize the systems, counter-ions (sodium or chloride) were added. The maximum step size of the steepest descent algorithm was 0.5 nm, and that of the conjugate gradient algorithm was 0.05 nm. The exit tolerance level of the steepest descent algorithm was set to and 103 that of the conjugate gradient algorithm to 10 kJ·mol−1·nm−1. Movement of the solute Cα atoms were restrained at a force constant of 103 kJmol−1nm−2, except for that of the A-loop. Calculations were performed with programs of the GROMACS [45] software package, using the AMBER99SB-ILDN [35] force field. After energy minimization, 100-ns-long NPSA MD simulation was carried out with a time step of 2 fs. Simulated annealing temperature scheme was applied as described in [22]. Simulated annealing temperature was rescaled and controlled for both solvent and solute. The temperature was gradually increased up to 323.15 K, then lowered back to 300 K in the first 20 ns, then the simulation was continued to 100 ns with constant temperature of 300 K. Pressure was coupled by the Parrinello–Rahman algorithm and a coupling time constant of 0.5 ps, compressibility of 4.5 × 10−5 bar−1 and reference pressure of 1 bar. Particle Mesh-Ewald summation was used for long range electrostatics. Van der Waals and Coulomb interactions had a cut-off at 11 Å. Coordinates were saved at regular time-intervals of 1 ps yielding 1 × 103 frames. Periodic boundary conditions were treated before analysis to center whole and recovered hydrated solute structures in the box. The original protein structure served as the basis of Cα fitting.

3.6. Scoring

AutoDock [25] estimates the binding free energy of the ligand (ΔGAD) with Equation (1) as a scoring function.
Δ G A D = W v d W ∑ i j ( A i j r i j 12 − B i j r i j 6 ) + W h b o u n d ∑ i j E ( t ) ( C i j r i j 12 − D i j r i j 10 ) + W e l e c ∑ i j q i q j ε ( r i j ) r i j + W s o l ∑ i j ( S i V j + S j V i ) e ( − r i j 2 2 ∂ 2 )
W terms are weighting constants calibrating to an experimentally determined set of free energies. Ligand atoms are represented by i, and protein atoms by j. A Lennard–Jones 12-6 dispersion/repulsion term, a directional 12–10 h-bonding term and a screened Coulombic electrostatic potential are included. A and B parameters are taken from the Amber force field. E(t) is a directional weight based on the angle, t, between the probe and the target atom. C and D parameters are assigned for well-depth calculations. The final term is a desolvation potential, V is the volume of the atoms surrounding a given atom and S is a solvation parameter for weighting [51]. δ is a distance weighting factor. The actual distance between the ith (ligand) and jth (target) atoms is marked with r.
The FITTED [13,14,15,34] scoring function estimates (ΔGFD) with the sum of various terms including the number of rotatable bonds, van der Waals and electrostatic interactions and directional H-bonding contributions as described in Equation (2).
Δ G F D = Δ G 0 + 0.14 N r o t + ∑ ( s c a l e   f a c t o r ) [ ( 0.26   U v d W i n h − p r o t + 0.035   U e l e c i n h − p r o t + 0.80 f h b ( Δ r , Δ α ) ) ]
Nrot is the number of rotatable bonds, Uvdw and Uelec are the van der Waals and electrostatic interactions based on the AMBER94 force field. The last term is the solvation contribution to the free energy of binding. Where fhb is the electrical field strength of hydrogen bonds, r is the length and α is the angle of hydrogen bonds.
E i j = ∑ i j N I N L [ A i j r i j 12 − B i j r i j 6 ] A i j = ε i j R i j 12 ; B i j = 2 ε i j R i j 6 ; R i j = R i + R j ; ε i j = ε i ε j
where εi and εj are the potential well depths in the equilibrium distance of atom pairs of identical types; εij is the potential well depth in equilibrium between the ith (ligand) and jth (target) atoms; Rij is the internuclear distance at equilibrium between ith (ligand) and jth (target) atoms; Ri and Rj are half equilibrium distances between ii and jj atom pairs of identical types, respectively; rij is the actual distance between the ith (ligand) and jth (target) atoms; NT is the number of target atoms; NL is the number of ligand atoms.

3.7. Ranking

The basis of the structural clustering and ranking of the docked ligand conformations was their AutoDock 4.2 binding free energy values. In the respective Tables, the serial number of ranks are represented. To create one rank [41], the ligand structure with the lowest calculated free energy of binding, and its neighboring docked ligand structures within 2 Å [52] were selected. Then new ranks were opened for the remaining structures, and clustering was repeated with the same protocol. The low serial number of a rank indicates an energetically favorable binding conformation. The actual rank (N) selected from all the ranks (M) is given in the format N/M.
RMSD. In all cases, the structural match of the docked (D in Equation (4)) binding mode to the crystallographic reference (C) was expressed as a root mean squared deviation (RMSD) value according to Equation (1).
R M S D = 1 N ∑ n = 1 N | D n − C n | 2
In Equation (4), N is the number of ligand heavy atoms, C is the space vector of the nth heavy atom of the crystallographic reference ligand molecule, D is the space vector of the nth heavy atom of the calculated ligand conformation. RMSDbest is the RMSD value of the ligand binding mode with the lowest RMSD.
The distance (d) between the S atom of C621 amino acid and the ligand atom that participates in the covalent binding was also measured to check the presence of covalent bond and to estimate the degree of translation necessary to move the prerequisite binding mode into the covalent binding mode (Figure 3). The dbest value is the smallest distance observed.
The AAmatch (%) is the rate of identical AAs present in two different binding pockets interacting with the ligand in a 3.5 Å cut-off distance. It is calculated by the results of Tables S5–S7.
NHA Number of heavy atoms of the agonist counts all the atoms except for hydrogens.
EINHA Efficiency index, the calculated free energy of binding is divided by the NHA of the respective agonist. The dimension is kcal/mol.
dcovalent The length of the covalent bond in Å.
Rankbest The scoring function of the program collects the results into ranks based on their calculated free energy of binding. The lowest rank contains the best energy. The rank that contains the model with the best RMSD value is the Rankbest.

4. Conclusions

Agonist binding to TRPA1 is a dynamic process involving structural changes of the target, first at the smaller scale of the A-loop, necessary for binding site activation, then at the whole of the target, required for channel activation. The present study identified the prerequisite binding modes of three agonists and showed how the binding of a ligand to the prerequisite site can forecast its successful docking to the final binding pocket. The prerequisite binding sites proved to be milestones on the association/dissociation pathway of the agonists, important in mechanism-based design. The present study also showed how the prerequisite binding modes affect the opening of the A-loop region, a central scene of the agonist binding mechanism. The time step measured in nanoseconds necessary for binding site activation is currently hidden from experimental methods, and only the co-operation with in silico approaches can shed light on them. Thus, amino acids identified along the dynamic binding pathway will serve as new target sites for the design of reversible binding of future agonists, beyond the well-known target of the covalent binding pocket of TRPA1.

Supplementary Materials

The following are available online at www.mdpi.com/article/10.3390/xxx/s1. Figure S1: The close-up of binding of BITC to TRPA1 mutant structure, Table S1: Covalent docking calculations performed by FITTED on the apo target with Aloop, Table S2: Non-covalent docking calculations performed by FITTED on the apo target with A-loop, Table S3: Non-covalent docking calculations performed by AutoDock on the apo target with A-loop, Table S4: Interacting (≤3.5 Å) amino acid residues of 6PQO with the non-covalent top ranked binding mode of JT010 (FITTED), Table S5: The interacting amino acids (within 3.5 Å) of the experimental binding position, the covalently docked and prerequisite binding modes of bodipy-iodoacetamide, Table S6: The interacting amino acids (within 3.5 Å) of the experimental binding position, the covalently docked and prerequisite binding modes of BITC, Table S7: The interacting amino acids (within 3.5 Å) of the experimental binding position, the covalently docked and prerequisite binding modes of JT010.

Author Contributions

Conceptualization, B.Z.Z., C.H. and E.P.; methodology, B.Z.Z., C.H. and R.B.; formal analysis, C.H.; investigation, B.Z.Z., R.B. and C.H.; resources, C.H. and E.P.; writing—original draft preparation, B.Z.Z. and C.H.; writing—review and editing, C.H. and E.P.; visualization, B.Z.Z., C.H.; supervision, C.H.; project administration, C.H.; funding acquisition, C.H. and E.P. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by the Hungarian National Research, Development and Innovation Office (K123836, K134214), 2017-1.2.1 NKP-2017-00002 (NAP-2; Chronic Pain Research Group), EFOP-3.6.1.-16-2016-0004 and GINOP 2.3.2-15-2016-00050 “PEPSYS”. This work was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences. The work was supported by the ÚNKP-20-5 and ÚNKP-20-3-I. New National Excellence Program of the Ministry for Innovation and Technology. Supported by the PTE ÁOK-KA No:2019/KA-2019-31. The project was supported by the European Union, co-financed by the European Social Fund. Project name and code: Comprehensive Development for Implementing Smart Specialization Strategies at the University of Pécs, EFOP-3.6.1-16-2016-00004. This project is supported by Eötvös Loránd Research Network.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article and Supplementary Material.

Acknowledgments

We acknowledge the support from the Governmental Information Technology Development Agency, Hungary.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. De Logu, F.; Nassini, R.; Materazzi, S.; Carvalho, G.M.; Nosi, D.; Rossi, D.D.; Marone, I.M.; Ferreira, J.; Li Puma, S.; Benemei, S.; et al. Schwann cell TRPA1 mediates neuroinflammation that sustains macrophage-dependent neuropathic pain in mice. Nat. Commun. 2017, 8, 1–16. [Google Scholar] [CrossRef] [PubMed]
  2. Paulsen, C.E.; Armache, J.P.; Gao, Y.; Cheng, Y.; Julius, D. Structure of the TRPA1 ion channel suggests regulatory mechanisms. Nature 2015, 520, 511–517. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Liu, C.; Reese, R.; Vu, S.; Rougé, L.; Shields, S.D.; Kakiuchi-Kiyota, S.; Chen, H.; Johnson, K.; Shi, Y.P.; Chernov-Rogan, T.; et al. A Non-covalent Ligand Reveals Biased Agonism of the TRPA1 Ion Channel. Neuron 2020, 109, 273–284. [Google Scholar] [CrossRef] [PubMed]
  4. Takaya, J.; Mio, K.; Shiraishi, T.; Kurokawa, T.; Otsuka, S.; Mori, Y.; Uesugi, M. A Potent and Site-Selective Agonist of TRPA1. J. Am. Chem. Soc. 2015, 137, 15859–15864. [Google Scholar] [CrossRef] [Green Version]
  5. Pozsgai, G.; Bátai, I.Z.; Pintér, E. Effects of sulfide and polysulfides transmitted by direct or signal transduction-mediated activation of TRPA1 channels. Br. J. Pharmacol. 2019, 176, 628–645. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Suo, Y.; Wang, Z.; Zubcevic, L.; Hsu, A.L.; He, Q.; Borgnia, M.J.; Ji, R.R.; Lee, S.Y. Structural Insights into Electrophile Irritant Sensing by the Human TRPA1 Channel. Neuron 2020, 105, 882–894. [Google Scholar] [CrossRef]
  7. Zhao, J.; Lin King, J.V.; Paulsen, C.E.; Cheng, Y.; Julius, D. Irritant-evoked activation and calcium modulation of the TRPA1 receptor. Nature 2020, 585, 141–145. [Google Scholar] [CrossRef]
  8. Chernov-Rogan, T.; Gianti, E.; Liu, C.; Villemure, E.; Cridland, A.P.; Hu, X.; Ballini, E.; Lange, W.; Deisemann, H.; Li, T.; et al. TRPA1 modulation by piperidine carboxamides suggests an evolutionarily conserved binding site and gating mechanism. Proc. Natl. Acad. Sci. USA 2019, 116, 26008–26019. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Tseng, W.C.; Pryde, D.C.; Yoger, K.E.; Padilla, K.M.; Antonio, B.M.; Han, S.; Shanmugasundaram, V.; Gerlach, A.C. TRPA1 ankyrin repeat six interacts with a small molecule inhibitor chemotype. Proc. Natl. Acad. Sci. USA 2018, 115, 12301–12306. [Google Scholar] [CrossRef] [Green Version]
  10. Wang, Y.Y.; Chang, R.B.; Waters, H.N.; McKemy, D.D.; Liman, E.R. The nociceptor ion channel TRPA1 is potentiated and inactivated by permeating calcium ions. J. Biol. Chem. 2008, 283, 32691–32703. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  11. Ábrányi-Balogh, P.; Petri, L.; Imre, T.; Szijj, P.; Scarpino, A.; Hrast, M.; Mitrović, A.; Fonovič, U.P.; Németh, K.; Barreteau, H.; et al. A road map for prioritizing warheads for cysteine targeting covalent inhibitors. Eur. J. Med. Chem. 2018, 160, 94–107. [Google Scholar] [CrossRef] [PubMed]
  12. Petri, L.; Egyed, A.; Bajusz, D.; Imre, T.; Hetényi, A.; Martinek, T.; Ábrányi-Balogh, P.; Keserű, G.M. An electrophilic warhead library for mapping the reactivity and accessibility of tractable cysteines in protein kinases. Eur. J. Med. Chem. 2020, 207, 1–9. [Google Scholar] [CrossRef] [PubMed]
  13. Corbeil, C.R.; Englebienne, P.; Moitessier, N. Docking ligands into flexible and solvated macromolecules. 1. Development and validation of FITTED 1.0. J. Chem. Inf. Model. 2007, 47, 435–449. [Google Scholar] [CrossRef]
  14. Pottel, J.; Therrien, E.; Gleason, J.L.; Moitessier, N. Docking ligands into flexible and solvated macromolecules. 6. Development and application to the docking of HDACs and other zinc metalloenzymes inhibitors. J. Chem. Inf. Model. 2014, 54, 254–265. [Google Scholar] [CrossRef] [PubMed]
  15. Therrien, E.; Englebienne, P.; Arrowsmith, A.G.; Mendoza-Sanchez, R.; Corbeil, C.R.; Weill, N.; Campagna-Slater, V.; Moitessier, N. Integrating medicinal chemistry, organic/combinatorial chemistry, and computational chemistry for the discovery of selective estrogen receptor modulatorswith FORECASTER, a novel platform for drug discovery. J. Chem. Inf. Model. 2012, 52, 210–224. [Google Scholar] [CrossRef]
  16. Bálint, M.; Horváth, I.; Mészáros, N.; Hetényi, C. Towards Unraveling the Histone Code by Fragment Blind Docking. Int. J. Mol. Sci. 2019, 20, 422. [Google Scholar] [CrossRef] [Green Version]
  17. Kevener, H.E.; Zhao, W.; Ball, D.M.; Babaoglu, K.; Qi, J.; White, S.W.; Lee, R.E. Validation of Molecular Docking Programs for Virtual Screening against Dihydropteroate Synthase. J. Chem. Inf. Model. 2009, 49, 444–460. [Google Scholar] [CrossRef] [PubMed]
  18. Castro-Alvarez, A.; Costa, A.M.; Vilarrasa, J. The Performance of several docking programs at reproducing protein-macrolide-like crystal structures. Molecules 2017, 22, 136. [Google Scholar] [CrossRef] [Green Version]
  19. Mena-Ulecia, K.; Tiznado, W.; Caballero, J. Study of the differential activity of thrombin inhibitors using docking, QSAR, molecular dynamics, and MM-GBSA. PLoS ONE 2015, 10, 1–21. [Google Scholar]
  20. Ramírez, D.; Caballero, J. Is It Reliable to Take the Molecular Docking Top Scoring Position as the Best Solution without Considering Available Structural Data? Molecules 2018, 23, 1038. [Google Scholar] [CrossRef] [Green Version]
  21. Gohlke, H.; Hendlich, M.; Klebe, G. Knowledge-based scoring function to predict protein-ligand interactions. J. Mol. Biol. 2000, 295, 337–356. [Google Scholar] [CrossRef]
  22. Bálint, M.; Jeszenoi, N.; Horváth, I.; Van Der Spoel, D.; Hetényi, C. Systematic exploration of multiple drug binding sites. J. Cheminform. 2017, 9, 65–79. [Google Scholar] [CrossRef]
  23. Sotriffer, C. Docking of Covalent Ligands: Challenges and Approaches. Mol. Inform. 2018, 37, 1–12. [Google Scholar] [CrossRef]
  24. Kumalo, H.M.; Bhakat, S.; Soliman, M.E.S. Theory and applications of covalent docking in drug discovery: Merits and pitfalls. Molecules 2015, 20, 1984–2000. [Google Scholar] [CrossRef] [Green Version]
  25. Morris, G.M.; Goodsell, D.S.; Halliday, R.S.; Huey, R.; Hart, W.E.; Belew, R.K.; Olson, A.J. Automated docking using a Lamarckian Genetic Algorithm and empirical binding free energy function. J. Comput. Chem. 1998, 19, 1639–1662. [Google Scholar] [CrossRef] [Green Version]
  26. Shi, S.; Yan, L.; Yang, Y.; Fisher-Shaulsky, J.; Thacher, T. An extensible and systematic force field, ESFF, for molecular modeling of organic, inorganic, and organometallic systems. J. Comput. Chem. 2003, 24, 1059–1076. [Google Scholar] [CrossRef] [PubMed]
  27. CDiscoVer, 98.0; Accelrys, Inc.: San Diego, CA, USA, 2001.
  28. Macpherson, L.J.; Dubin, A.E.; Evans, M.J.; Marr, F.; Schultz, P.G.; Cravatt, B.F.; Patapoutian, A. Noxious compounds activate TRPA1 ion channels through covalent modification of cysteines. Nature 2007, 445, 541–545. [Google Scholar] [CrossRef] [PubMed]
  29. Warren, L.D. The PyMOL Molecular Graphics System; Version 2.0; Schrödinger, LLC.: New York, NY, USA, 2002. [Google Scholar]
  30. Stewart, J.J.P. Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements. J. Mol. Model. 2007, 13, 1173–1213. [Google Scholar] [CrossRef] [Green Version]
  31. Stewart, J.J.P. Optimization of parameters for semiempirical methods VI: More modifications to the NDDO approximations and re-optimization of parameters. J. Mol. Model. 2013, 19, 1–32. [Google Scholar] [CrossRef] [Green Version]
  32. Gasteiger, J.; Marsili, M. Iterative partial equalization of orbital electronegativity-a rapid access to atomic charges. Tetrahedron 1980, 36, 3219–3228. [Google Scholar] [CrossRef]
  33. O’Boyle, N.M.; Banck, M.; James, C.A.; Morley, C.; Vandermeersch, T.; Hutchison, G.R. Open Babel. J. Cheminform. 2011, 3, 1–14. [Google Scholar]
  34. Moitessier, N.; Therrien, E.; Hanessian, S. A method for induced-fit docking, scoring, and ranking of flexible ligands. Application to peptidic and pseudopeptidic β-secretase (BACE 1) inhibitors. J. Med. Chem. 2006, 49, 5885–5894. [Google Scholar] [CrossRef]
  35. Wang, J.; Wolf, R.M.; Caldwell, J.W.; Kollman, P.A.; Case, D.A. Development and testing of a general Amber force field. J. Comput. Chem. 2004, 25, 1157–1174. [Google Scholar] [CrossRef]
  36. Schrödinger Release 2020-4: Maestro; Schrödinger, LLC.: New York, NY, USA, 2021.
  37. Dupradeau, F.Y.; Pigache, A.; Zaffran, T.; Savineau, C.; Lelong, R.; Grivel, N.; Lelong, D.; Rosanski, W.; Cieplak, P. The R.E.D. tools: Advances in RESP and ESP charge derivation and force field library building. Phys. Chem. Chem. Phys. 2010, 12, 7821–7839. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Schmidt, M.W.; Baldridge, K.K.; Boatz, J.A.; Elbert, S.T.; Gordon, M.S.; Jensen, J.H.; Koseki, S.; Matsunaga, N.; Nguyen, K.A.; Su, S.; et al. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347–1363. [Google Scholar] [CrossRef]
  39. Sousa Da Silva, A.W.; Vranken, W.F. ACPYPE—AnteChamber PYthon Parser interfacE. BMC Res. Notes 2012, 5, 1–8. [Google Scholar] [CrossRef] [Green Version]
  40. Wang, J.; Wang, W.; Kollman, P.A.; Case, D.A. Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 2006, 25, 247–260. [Google Scholar] [CrossRef]
  41. Zsidó, B.Z.; Börzsei, R.; Szél, V.; Hetényi, C. Determination of Ligand Binding Modes in Hydrated Viral Ion Channels to Foster Drug Design and Repositioning. J. Chem. Inf. Model. 2021, 8, 4011–4022. [Google Scholar] [CrossRef]
  42. Berman, H.M.; Battistuz, T.; Bhat, T.N.; Bluhm, W.F.; Bourne, P.E.; Burkhardt, K.; Feng, Z.; Gilliland, G.L.; Iype, L.; Jain, S.; et al. The protein data bank. Acta Crystallogr. Sect. D Biol. Crystallogr. 2002, 58, 899–907. [Google Scholar] [CrossRef]
  43. Zsidó, B.Z.; Hetényi, C. Molecular structure, binding affinity, and biological activity in the epigenome. Int. J. Mol. Sci. 2020, 21, 4143. [Google Scholar] [CrossRef]
  44. Waterhouse, A.; Bertoni, M.; Bienert, S.; Studer, G.; Tauriello, G.; Gumienny, R.; Heer, F.T.; De Beer, T.A.P.; Rempfer, C.; Bordoli, L.; et al. SWISS-MODEL: Homology modelling of protein structures and complexes. Nucleic Acids Res. 2018, 46, W296–W303. [Google Scholar] [CrossRef] [Green Version]
  45. Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A.E.; Berendsen, H.J.C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701–1718. [Google Scholar] [CrossRef]
  46. Fliszár-Nyúl, E.; Faisal, Z.; Mohos, V.; Derdák, D.; Lemli, B.; Kálai, T.; Sár, C.; Zsidó, B.Z.; Hetényi, C.; Horváth, Á.I.; et al. Interaction of SZV 1287, a novel oxime analgesic drug candidate, and its metabolites with serum albumin. J. Mol. Liq. 2021, 333, 1–10. [Google Scholar] [CrossRef]
  47. Zsidó, B.Z.; Balog, M.; Erős, N.; Poór, M.; Mohos, V.; Fliszár-Nyúl, E.; Hetényi, C.; Nagane, M.; Hideg, K.; Kálai, T.; et al. Synthesis of spin-labelled bergamottin: A potent CYP3A4 inhibitor with antiproliferative activity. Int. J. Mol. Sci. 2020, 21, 508. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Mohos, V.; Fliszár-Nyúl, E.; Ungvári, O.; Bakos, É.; Kuffa, K.; Bencsik, T.; Zsidó, B.Z.; Hetényi, C.; Telbisz, Á.; Özvegy-Laczka, C.; et al. Effects of chrysin and its major conjugated metabolites chrysin-7-sulfate and chrysin-7-glucuronide on cytochrome P450 enzymes and on OATP, P-gp, BCRP, and MRP2 transporters. Drug Metab. Dispos. 2020, 48, 1064–1073. [Google Scholar] [CrossRef] [PubMed]
  49. Mohos, V.; Fliszár-Nyúl, E.; Lemli, B.; Zsidó, B.Z.; Hetényi, C.; Mladěnka, P.; Horký, P.; Pour, M.; Poór, M. Testing the pharmacokinetic interactions of 24 colonic flavonoid metabolites with human serum albumin and cytochrome P450 enzymes. Biomolecules 2020, 10, 409. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  50. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926–935. [Google Scholar] [CrossRef]
  51. Huey, R.; Morris, G.M.; Olson, A.J.; Goodsell, D.S. A semiempirical free energy force field with charge-based desolvation. J. Comput. Chem. 2007, 28, 1145–1152. [Google Scholar] [CrossRef] [PubMed]
  52. Hetényi, C.; Van Der Spoel, D. Blind docking of drug-sized compounds to proteins with up to a thousand residues. FEBS Lett. 2006, 580, 1447–1450. [Google Scholar] [CrossRef] [Green Version]
Pharmaceuticals 14 00988 g001 550
Figure 1. The Lewis structures of three TRPA1 agonists, JT010, benzyl-isothiocyanate (BITC) and bodipy-iodoacetamide, and the PDB IDs of their complexes with the TRPA1 receptor. The original molecular structures were restored prior to the covalent bond formation. A chlorine was added to JT010, an iodine to bodipy-iodoacetamide and the geometry of the N=C=S bond of BITC was restored. The atoms that participate in the formation of the covalent bond are marked by black dots.
Figure 1. The Lewis structures of three TRPA1 agonists, JT010, benzyl-isothiocyanate (BITC) and bodipy-iodoacetamide, and the PDB IDs of their complexes with the TRPA1 receptor. The original molecular structures were restored prior to the covalent bond formation. A chlorine was added to JT010, an iodine to bodipy-iodoacetamide and the geometry of the N=C=S bond of BITC was restored. The atoms that participate in the formation of the covalent bond are marked by black dots.
Pharmaceuticals 14 00988 g001
Pharmaceuticals 14 00988 g002 550
Figure 2. (A) The TRPA1 ion channel shown as grey cartoon representation, the agonist binding site is shown by the binding of bodipy-iodoacetamide (PDB: 6v9v) as teal spheres. The A-loop is highlighted with blue; (B) the movement of the A-loop during ligand binding. The figure was prepared with the superposition of 6v9w on 6pqo. The blue loop is A-loop in the holo, open form, and the red is A-loop in the apo, closed form. The rest of the binding site is shown with grey cartoon. bodipy-iodoacetamide is shown in teal sticks, and C621 in all atom representation grey sticks. (C) The close-up of the binding of bodipy-iodoacetamide to the agonist binding site of TRPA1, interacting amino acids are shown as grey thick lines in all atom representation, bodipy-iodoacetamide is shown as teal all atom representation sticks.
Figure 2. (A) The TRPA1 ion channel shown as grey cartoon representation, the agonist binding site is shown by the binding of bodipy-iodoacetamide (PDB: 6v9v) as teal spheres. The A-loop is highlighted with blue; (B) the movement of the A-loop during ligand binding. The figure was prepared with the superposition of 6v9w on 6pqo. The blue loop is A-loop in the holo, open form, and the red is A-loop in the apo, closed form. The rest of the binding site is shown with grey cartoon. bodipy-iodoacetamide is shown in teal sticks, and C621 in all atom representation grey sticks. (C) The close-up of the binding of bodipy-iodoacetamide to the agonist binding site of TRPA1, interacting amino acids are shown as grey thick lines in all atom representation, bodipy-iodoacetamide is shown as teal all atom representation sticks.
Pharmaceuticals 14 00988 g002
Pharmaceuticals 14 00988 g003 550
Figure 3. The reaction schemes of JT010, BITC, bodipy-iodoacetamide. Atoms, that participate in the formation of the covalent bond are highlighted by a black dot (●). The distance of these two atoms is referred to as d.
Figure 3. The reaction schemes of JT010, BITC, bodipy-iodoacetamide. Atoms, that participate in the formation of the covalent bond are highlighted by a black dot (●). The distance of these two atoms is referred to as d.
Pharmaceuticals 14 00988 g003
Pharmaceuticals 14 00988 g004 550
Figure 4. The MDholo simulation starting from the experimental binding position, with the covalent bond cut. Interaction energy distribution of the interacting amino acids of the pocket (inner prerequisite site) is shown during the MD simulation. The structural figures are snapshots of the binding position of BITC at the stated time frame of the MD simulation. The protein is shown in grey cartoon and BITC with teal sticks. The A-loop is marked with blue in the open conformation. C621 amino acid is also shown as all atom sticks representation. The teal arrows indicate the movement of BITC. Lennard–Jones interaction energies calculated between BITC and the TRPA1 target amino acids are shown per residue on the diagram.
Figure 4. The MDholo simulation starting from the experimental binding position, with the covalent bond cut. Interaction energy distribution of the interacting amino acids of the pocket (inner prerequisite site) is shown during the MD simulation. The structural figures are snapshots of the binding position of BITC at the stated time frame of the MD simulation. The protein is shown in grey cartoon and BITC with teal sticks. The A-loop is marked with blue in the open conformation. C621 amino acid is also shown as all atom sticks representation. The teal arrows indicate the movement of BITC. Lennard–Jones interaction energies calculated between BITC and the TRPA1 target amino acids are shown per residue on the diagram.
Pharmaceuticals 14 00988 g004
Pharmaceuticals 14 00988 g005 550
Figure 5. The MDrank1 simulation starting from the top 1st docked prerequisite binding mode of BITC on the apo TRPA1. Interaction energy distribution of the interacting amino acids of the outer prerequisite site are shown during the MD simulation. The loop motion is quantified by the distances between the N atom of N615 of the A-loop and O atom of Q676 of the opposite loop (black dotted line on the t = 38 ns structural plot, dLOOPS) shown as joint black boxes. The structural figures are snapshots of the binding position of BITC at the stated time frame of the MD simulation. The protein is shown in grey cartoon and BITC with teal sticks. The A-loop is marked with blue and red in open and closed conformations, respectively. C621 amino acid is also shown as all atom sticks representation. The red and blue arrows show the movement of the A-loop and the teal arrows the movement of BITC. Lennard–Jones interaction energies calculated between BITC and the TRPA1 target amino acids are shown per residue on the diagram.
Figure 5. The MDrank1 simulation starting from the top 1st docked prerequisite binding mode of BITC on the apo TRPA1. Interaction energy distribution of the interacting amino acids of the outer prerequisite site are shown during the MD simulation. The loop motion is quantified by the distances between the N atom of N615 of the A-loop and O atom of Q676 of the opposite loop (black dotted line on the t = 38 ns structural plot, dLOOPS) shown as joint black boxes. The structural figures are snapshots of the binding position of BITC at the stated time frame of the MD simulation. The protein is shown in grey cartoon and BITC with teal sticks. The A-loop is marked with blue and red in open and closed conformations, respectively. C621 amino acid is also shown as all atom sticks representation. The red and blue arrows show the movement of the A-loop and the teal arrows the movement of BITC. Lennard–Jones interaction energies calculated between BITC and the TRPA1 target amino acids are shown per residue on the diagram.
Pharmaceuticals 14 00988 g005
Table
Table 1. Covalent docking calculations performed by FITTED [13,14,15].
Table 1. Covalent docking calculations performed by FITTED [13,14,15].
Ligand NameJT010BITCBodipy-Iodoacetamide
HOLO target
AAmatch (%)100%100%100%
RMSDbest (Å)2.282.053.87
Rankbest1/31/31/3
ΔGFD (kcal/mol)−84.1−77.7−44.3
NHA c231022
EINHA d (kcal/mol)3.667.772.01
dcovalent (Å)1.8 (1.8) a1.8 (1.8) a1.8 (1.8) a
APO target b
AAmatch (%)100%60%66.6%
RMSDbest (Å)6.824.756.55
Rankbest1/51/51/5
ΔGFD (kcal/mol)−77.4−73.8−43.1
NHA c231022
EINHA d (kcal/mol)3.367.381.96
dcovalent (Å)1.81.81.8
a The experimental covalent bond lengths are shown in brackets; b Without A-loop; c Number of heavy atoms; d Efficiency index.
Table
Table 2. Non-covalent docking calculations performed by FITTED.
Table 2. Non-covalent docking calculations performed by FITTED.
Ligand NameJT010BITCBodipy-Iodoacetamide
HOLO target
ΔGFD (kcal/mol)−46.1−32.4−13.7
Rankbest10/101/108/10
AAmatch (%)100%100%100%
dbest (Å)3.64.08.7
APO target a
ΔGFD (kcal/mol)−33.4−26.70.5
Rankbest3/51/54/5
AAmatch (%)100%60%33.3%
dbest (Å)3.53.93.3
a Without A-loop.
Table
Table 3. Non-covalent docking calculations performed by AutoDock 4.2 [25].
Table 3. Non-covalent docking calculations performed by AutoDock 4.2 [25].
Ligand NameJT010BITCBodipy-Iodoacetamide
HOLO target
ΔGAD (kcal/mol)−6.8−3.8−5.9
Rankbest1/51/14/4
AAmatch (%)100%80%66%
dbest (Å)3.66.54.0
APO target a
ΔGAD (kcal/mol)−5.16−3.74−5.26
Rankbest1/31/23/5
AAmatch (%)50%40%0%
dbest (Å)7.57.27.3
a Without A-loop.
Table
Table 4. The details of the MD simulations performed to unravel the binding mechanism of BITC.
Table 4. The details of the MD simulations performed to unravel the binding mechanism of BITC.
Simulation NameTRPA1LigandChange in A-LoopMovement of the Agonist
MDapoApo protein-No change in A-loop conformation-
MDholo,PSAHolo proteinExperimentalNo change in A-loop conformationUnbinding–binding
MDrank1Apo proteinRank 1 docked ligand binding modeA-loop flipping to the active conformationDissociation–association
MDrank3Apo proteinRank 3 docked ligand binding modeNo change in A-loop conformationDissociation
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zsidó, B.Z.; Börzsei, R.; Pintér, E.; Hetényi, C. Prerequisite Binding Modes Determine the Dynamics of Action of Covalent Agonists of Ion Channel TRPA1. Pharmaceuticals 2021, 14, 988. https://0-doi-org.brum.beds.ac.uk/10.3390/ph14100988

AMA Style

Zsidó BZ, Börzsei R, Pintér E, Hetényi C. Prerequisite Binding Modes Determine the Dynamics of Action of Covalent Agonists of Ion Channel TRPA1. Pharmaceuticals. 2021; 14(10):988. https://0-doi-org.brum.beds.ac.uk/10.3390/ph14100988

Chicago/Turabian Style

Zsidó, Balázs Zoltán, Rita Börzsei, Erika Pintér, and Csaba Hetényi. 2021. "Prerequisite Binding Modes Determine the Dynamics of Action of Covalent Agonists of Ion Channel TRPA1" Pharmaceuticals 14, no. 10: 988. https://0-doi-org.brum.beds.ac.uk/10.3390/ph14100988

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop