On the Quest of Small Molecules that Can Mimic Psalmotoxin-1, a Powerful Peptidic Modulator of the Acid Sensing Channel ASIC1a

Acid-sensing ion channels (ASICs) are thought to play a key role in a number of pathologies, from neuronal injury to pain sensation, while no drug has yet been approved as their modulator. This work was devised to asses roughly, yet from rst principles, the relative energies of binding in the most important acidic pocket of cASIC1a, thereby paving the way to nd small molecules that can mimic Pctx1, the most powerful peptidic modulator of cASIC1a. To this end, MD simulations for the overall conformation, and QM-MM simulations specically for the location of hydrogen atoms, allowed disentangling the relative weight of the various non-bonded interactions between PcTx1 and cASIC1a. Main weight could be attributed to deeply buried salt bridges formed by the guanidinium end chains of residues Arg26 and Arg27 on PcTx1 and carboxylate end chain of distant residues Asp and Glu on cASIC1a. Rewardingly, in a preliminary attempt at exploiting these observations toward a small-molecule modulator, a Arg26-Arg27 stretch, excised from the PcTx1-cASIC1a complex and slightly simplied, on automated rigid docking on ligand-free receptor was observed to form most of the above salt bridges.


Introduction
Acid-sensing ion channels (ASICs) are low-pH activated Na + -permeable ion channels encoded, as far as it is presently known, by four genes to give homo-and hetero-trimeric ASIC1a, ASIC1b, ASIC2a, ASIC3, and ASIC4 [1]. They are widely expressed in central and peripheral nervous systems and are believed to be involved with many pathologies, like neuronal injury, pain sensation and seizure, while affecting synaptic function and plasticity [2]. Therefore, ASICs represent foreseeable drug targets, particularly as far as pain and ischemic stroke are concerned, although this goal has remained underdeveloped for a number of reasons. First of all, powerful exogenous modulators of ASIC channels had only been found with peptides of forty or more residues, which are prone to rapid degradation in blood plasma [3]. Second, the task was faced with a penury of structural data at atomic resolution, limited to the cASIC1a channel, achieved through X-ray diffraction at ca 3.0 Å resolution while leaving the architecture of the cytoplasmic terminal domains largely elusive [4].
Third, such structural limitations persisted for complexes with modulators, such as that of cASIC1a with the Tarantula peptide Psalmotoxin-1 (PcTx1) [5][6]. It is also to be noticed that the ASIC3 channel is particularly neglected, despite evidence that, by exhibiting high expression in sensory neurones while upon activation not being fully inactivated, it is a particularly interesting analgesic target [7].
These obstacles toward e cient small-molecule modulators of ASIC channels are now beginning to be removed. From one side, like in all areas of structural biology, cryo-electron microscopy is moving forward also with ASICs, having recently allowied to tackle for the rst time human ASICs. This led to unravel the structure of hASIC1a and its complex with the black mamba peptide mambalgin-1 at 3.56 and and 3.90 Å resolution [8]. This fast moving technique has also recently allowed to unravel the whole apo structure of the cASIC1a channel, both in the pH = 7.0 desensitized and the pH = 8.0 resting conformations, at resolutions of ~ 2.8 and 3.7 Å, clarifying details of the way the channel works [9]. Importantly, when present work was nearly completed, two powerful small-molecule modulators of cASIC1a and rASIC1a have been presented, which, however, involve an acidic pocket different from those known for the long peptides [10].
What is still completely lacking is a quantitative evaluation of the relative importance of the various non-bonded interactions between ASIC channels and their modulators, which could help modeling the long sought small molecules that can powerfully modulate the ASICs. This work was conceived to pave the way toward such small molecules by investigating in silico the complex between the ion channel cASIC1a and the spider peptidic toxin PcTx1. The choice of PcTx1 was suggested by it being the most powerful modulator of ASIC1a at the most important acidic pocket, while other acidic pockets proved to be of lesser importance [11] [12]. On the other hand, at the onset of this study, structural data at atomic resolution were only available for the PcTx1-cASIC1a complex [5].
As only protein residues are involved in the PcTx1-cASIC1a complex, for classical molecular dynamics (MD) simulations recourse was taken here to a long established xed-charge classical force eld rather than more elaborated force elds that are at an earlier stage of development. The same approach was taken for quantum mechanics-molecular mechanics (QM-MM) simulations, where the MD force eld was coupled to a widely established QM code, as detailed in the following.

Results And Discussion
This work has taken the challenge of evaluating in an aqueous medium the relative energies of the non-bonded interactions between the cASIC1a channel and its powerful modulator, the Tarantula peptide toxin Psalmotoxin-1 (PcTx1) [5][6]. Doing that by MD in aqueous solution allows viewing the cASIC1a-PcTx1 complex in a dynamic status, a task that the crystal data [5][6] cannot grant. In addition, to avoid any deterministic bias possibly introduced by the classical MD force eld, the geometry at the interaction loci was re-assessed from rst principles through QM-MM simulations with Density Functional Theory (DFT) with functional that are most suited for structural analysis, aided by account for dispersion, as detailed in the following.
To the above purpose, two models were built, as detailed in the Experimental Section, from the crystal structure at low pH (5.5) and high pH (7.25) of the cASIC1a-PcTx1 complex [5], i.e., primarily relying on the ligand-receptor interactions that resulted from the crystal study [5]. All missing atoms on both PcTx1 and the channel extracellular part were implemented while nearly all transmembrane domains were removed to alleviate the computational burden. While each subunit of the homotrimer hosts a PcTx1 molecule, only one was dealt with in this work. This is shown in Fig. 1 for the low pH complex, where PcTx1 (labeled as chain M, in red) inserts between chain A (yellow) and C (blue) of the homotrimer, digging deeply into the acidic pocket. To give a rough perspective of the spot, key residues are highlighted in Fig. 1 on both the ligand and the receptor. Residue R26 of PcTx1 gets close to D350 on chain C, while the adjacent R27 (not shown) does so with E220 on chain A of the receptor. Actually, it is only when the side chains are displayed in full for these residues that distances can be fully appreciated, as shown in the following. Highlighting residue F30 only serves to give a perception of the ligand sequence.
Each model was subjected to MD without any restraint on any atom, except for applying harmonic forces to where the transmembrane domains had been cut. For the QM-MM simulations all interacting residues between PcTx1 and cASIC1a were included, as resulting from the X-ray diffraction study of the crystal [5]. Water molecules that, from MD, turned out to be interacting with the selected peptidic residues were also added. In order to avoid duplication of statements, all interacting residues, as included in the QM-MM simulations, are indicated in the following Table 1 for the low pH and Table  2 for the high pH systems. Each complex was divided into the same four zones that were sketched in the original X-ray diffraction study for the low pH complex [5]. Each zone pivots on one or two PcTx1 residues, R26, R27, R28, and W7-W24, which interact with the closest residues on chains A and C of the homotrimer. Each zone was subjected separately to MD, followed by QM-MM simulations. The reason for the latter was to assess from rst principles the hydrogen bonds which are not treated in depth by the xed-charge force eld of classical MD.
The results for the low pH complex are compiled in Table 1. The rst column lists the PcTx1 residues that interact with those of cASIC1a in the second column, specifying, in the order, chain name, residue name and number, and atom name. Columns 3 to 6 tabulate the interatomic distances corresponding to the rst two columns, the third column listing crystal data [5]. In the fourth column data are reported for the pre-MD, minimized model in aqueous solution with restraints only applied to where the transmembrane domains had been cut. The fth column reports data for averaged distances from MD simulations. The sixth column lists averaged distances from QMMM simulations. The seventh column tabulates the interaction non-bond energies between the residues in columns 1-2, highlighting in boldface key values. Finally, the last column shows the HOMO-LUMO energy gap for each group of interactions. A similar compilation of data can be seen in Table 2 for the high pH complex.  H-bond and salt-bridge.
In parallel, molecular electronic potential maps (MEP) for each group of interactions, as in Table 1 and 2, are visualized in the next four illustrations. Comparisons are aided by placing side to side, for each group of interacting residues, MEPs for the low and the high pH systems. As these MEPs were built directly from the QM code, they require more attention than any a posteriori illustration from presentation graphics software. The extra attention is more than repaid by a direct view on what was going on with the computation from rst principles. What also helps, are chain, residue and atom labels set in a All that said, we are now in a position to analyze the results of this study. First about the very high values of the HOMO-LUMO energy gap, which speaks for high stability of the complex in solution, accompanied by some conformational changes, as discussed in the following. It also comes to the eyes the extremely small stabilizing non-bond energy for the aromatic 'embrace' between W7 and W24 of PcTx1 and F351 of cASIC1a for both the low (Table 1) and the high pH complex ( Table 2). Concerned is a residue, F351, that lies in a small region important for pH-dependent gating and PcTx1 modulation of ASIC1a [13].
What con dence is deserved by the non-bond energies in Table 1 and 2, as disentangled with the NAMDenergy plugin? The short answer is that the NAMDenergy plugin gives what the force eld used, CHARMM36, can give. This was applied to the ensemble in a periodic box under Particle Mesh Ewald (PME) conditions, thus taking also full electrostatic interactions into account. In addition, it is important noticing that the NAMDenergy plugin was applied to structures geometrically tempered with QM-MM simulations, i.e., by using rst principles, which, as said above, is particularly important for hydrogen bonds, which are not treated in depth by the xed-charge CHARMM36 force eld. With the aromatic 'embrace', visualized by the MEP maps of Fig. 2, the attractive forces are limited to dipole-dipole, dipole-induced dipole, and London (dispersion) interactions, treated in terms of the Lennard-Jones potential. No strong electrostatic forces come into play.  Fig. 3, left. Notably, there is no stabilizing interaction between residues R191 and D238, contrary to expectations from the crystal structure, (third column, same row) [5]. Under MD the side chains for these two residues were moved away from one another, the distance between R191/NH1 and D238/OD1 becoming as large as 6.8 Å. This remained nearly unaltered under QM-MM simulations, speaking for lack of stabilization. That is not the case of the high pH complex, where a salt bridge is retained between R191 and D238 ( Table 2 and Fig. 3, right). To this concern, it is to be noticed that a salt bridge without any hydrogen bond was established by QM-MM simulations between R191 and D238. Despite the lack of any hydrogen bond, a quite high non-bond energy was calculated for such interaction. [1] These observations indicate that PcTx1, in order to stick to the trimer, was mainly looking for interactions that generate a salt bridge.
With the group pivoting on residue R27 of PcTx1, a non-bonded interaction with E220 of chain A, generating a salt bridge, occurs with both the low and the high pH complexes. This is reported in Table 1 and 2 and visualized in Fig. 4. Otherwise, the two complexes differ markedly as to non-bonded interactions. The low pH complex is strongly stabilized by non-bonded interactions between R27 of PcTx1 and D408 of the receptor, while no such interaction occurs with the high pH complex, replaced by a hydrogen bond with G218.
With the last group of interacting residues, the one pivoting on R28 of PcTx1, it is seen from Table 1 and 2, and illustrated in Fig. 5, that, once the highly stabilizing salt bridges are set in place with R26 and R27, R28 can only interact with the trimer via hydrogen bonds. This explains why low non-bond energies were calculated in this case, especially low for the high pH complex.
From all these observations it was devised here to model a ligand that could form salt bridges with the two most important couples of residues as from Table 1, i.e., D238 and D350 on chain C for R26, and the distant E220 and D408 on chain A for R27. In a rst attempt, a R26-R27 stretch was excised from PcTx1 in the complex with cASIC1a for the low pH complex, as equilibrated with MD (Fig. 6, top). Un lled valencies were satis ed by changing both CHO and NH to H, while imposing the rigid conformation of the crystal that xes the two guanidinium ends at at distance of nearly 10 Å as required to interact with the widely separated carboxylates (Fig. 6, bottom). This modi ied-R26-R27 molecule was used for automated rigidbody docking on the complex freed from PcTx1. It was a pleasure observing that three of the expected salt bridges were formed, to D408 and E220 from one side and D238 from the other side (Fig. 7), mimicking the allosteric way of doing of PcTx1. 4 Admittedly, this is an arti cial rst approach to modeling a ligand that, by exploiting only a few points of binding to the receptor, docks into the most important acidic pocket of ASIC1a in the way that PcTx1 does. Nonetheless, this simple approach shows that modeling an apprppriately rigid small molecule that, by avoiding any great entropy loss, imitates the mode of binding of PcTx1 to ASIC1a should be possible.

Conclusions
A quantitative analysis carried out in this work of the distribution of non-bonded forces acting within the PcTx1 acidic pocket of cASIC1a in aqueous solution has shown a strict cooperativity of forces, like in a spider web. Nonetheless, MD and QM-MM simulations revealed that the non-bonded interactions that contribute most to the stability of the complex are salt bridges formed between end-chain guanidinium functionalities of PcTx1 and end-chain carboxylates of cASIC1a.
Ful lling the aim of furnishing indications as to model a small molecule that retains the ability of PcTx1 of allosterically modulating the receptor, the above results suggested to armor a small molecule of two guanidinium teeth that can grasp to carboxylate ends that are separated by as much as nearly 10 Å, in the acidic pocket. This was implemented by excising a R26-R27 stretch from PcTx1 and modifying it to saturate the un lled valencies. Rewardingly, on automated rigid-body docking of this stabilized R26-R27 molecule on the ligand free complex, most of the expected salt bridges were formed, as it occurs with the PcTx1-cASIC1a complex.
Treating in silico protein systems is a very challenging undertaking, if not else because their structural stability is strikingly dependent from the in vivo environment [15]. However, this simple in silico approach allowed observing the dynamics in the acidic pocket and the forces involved, complementing the static view from X-ray diffraction. This should prove of heuristic value in designing an effective modulator of ASICs.

Experimental Section
Modeling: The receptor was modeled on the X-ray diffraction crystal structure of the complex between PcTx1 and cASIC1a at either low (5.5) or high pH (7.25), 2.8 and 3.0 Å resolution, Protein Data Bank, PDB code 4FZ0 and 4FZ1, respectively [5]. Both NAG (2-acetamido−2-deoxy-beta-d-glucopyranose) and GOL (propane−1,2,3-triol) residues were discarded and the following residues were eliminated from the intracellular domains, without capping:  [19], and solvating in a periodic box with TIP3 water with 12 Å padding, while neutralizing and adding 0.15M NaCl. For the low pH complex, this resulted in a box of dimensions x = 123, y = 125, z = 115 Å for a total of 167716 atoms, 51019 residues, of which 49520 waters, while for the high pH complex, box dimensions x = 126, y = 117, z = 111 Å, for a total 153479 atoms, 46249 residues, of which 44776 waters.
Methods : Minimization of the above ensembles was carried out with software NAMD, v 2.13 [20] by rst restraining all atoms except water and then only the residues where the transmembrane domains had been cut. Gradual heating to 310K was followed by constant number, volume and temperature (NVT) and then constant number, pressure and temperature (NPT) simulations with restraints like for minimization. The integrator timestep was kept as low as 1.0 fs to allow treating adequately vibrations involving H-atoms under Langevin stabilization with low damping and target 1.01325 bar, under PME and with 1-scaling. Rigid bonds were only used for water molecules. Constant RMS during NPT for the low pH ensemble was reached with 18 ns trajectory and was continued to 30 ns. For the high pH ensemble shorter trajectories of 10 and 15 ns were su cient. Averaged distances, as in Table 1 and 2, are from the avpos algorithm available with VMD software[18].
The QM-MM simulations were set up according to the same general methodology previously used with complexes between small molecules and ribosome [21], while subdividing the system into four zones for the interactions of PcTx1 key residues R26, R27, R28, and W7-W24 with the protein trimer, as delineated for the crystal structure [5]. The atoms included in the QM part were detailed above along the text and are speci ed in Tables 1 and 2

Declarations Acknowledgments
The author gratefully acknowledges the CINECA consortium for the award cASIC1-HP10C8KQFV under the ISCRA initiative, which allowed generous access to high-performance computing resources.
Ethical approval: The author declares to have followed all indications in the section "Compliancve with Ethycal Standards" and not to be infringing any of them.
Funding: Funding to this work was granted by the CINECA consortium, award cASIC1-HP10C8KQFV under the ISCRA initiative, which only granted access to to high-performance computing resources.
Con ict of interest/Competing interest: The author has no con ict of interest and no competing interest to declare.
Availability of data and material: Supplementary material is included in the submission.  Figure 1 Wall-eye stereo representation based on the crystal structure of the cASIC1a-PcTx1 complex at low pH (5.5) [5]. It is seen that peptide PcTx1 (in red) inserts into the homotrimer receptor between subunit A (in yellow) and C (in blue). Highlighted is residue R26, which interacts with both F30 on the same peptide chain and D350 on chain C of the receptor. Also highlighted on receptor chain A is residue E220, which interacts with R27 (not shown to avoid overlapping of labels) on PcTx1.

Figure 2
In parallel with data on Table 1, shown here are translucent MEP maps from QM-MM simulations for the group of interactions pivoting on residues W76 and W24 on PcTx1 in complex with cASIC1a at low pH (5.5). Red MEPs stand for high positive potential, degrading toward blue for negative potentials. It is seen that the aromatic 'embrace' formed by residue F351, on cASIC1a subunit C, between W7 and W24 for the low pH (left) and high pH complex (right) do not entail any major electrostatic interaction.

Figure 3
Left: In parallel with data on  Table 2, and for the same above system at high pH (7.25), translucent MEP maps are shown here, showing that the R191/NH1-D350/OD1 interaction at low pH is replaced at high pH by the R191/NH2-D238/OD2 interaction, like in the crystal. Notable involvement of electrostatics in these interactions is indicated by facing opposing colors Figure 4 Left: In parallel with data on Table 1 and 2, shown here are translucent MEP maps from QM-MM simulations for the group of interactions pivoting on R27 on PcTx1 in complex with cASIC1a at low pH (5.5). Interacting atoms are highlighted at 25% of the vdW size, in blue for nitrogen and red for oxygen. Red and blue MEPs stand for for high positive and high negative potential, respectively, with intermediate colors for intermediate potentials. It is seen that atom NH2 of R27 interacts with both atom OD1 of D408 and atom OE1 of E220 in the receptor, entailing great electrostatic attraction. Right: In parallel with data on Table 2, and for the same above system at high pH (7.25), it is seen that the receptor conformation has so much changed that the strongly negative end chain of D408 is out of reach by R27.

Figure 5
Left: In parallel with data on  Table 2 for the same above system, that at high pH (7.25) the O-N interaction is lost, replaced by interaction with a water molecule. All such interactions entail moderate electrostatics, as indicated by the colors.

Figure 6
Top: the R26_R27 stretch as excised from PcTx1 of the MD-equilibrated complex with cASIC1a at low pH (5.5). Bottom: The R26_R27 stretch simpli ed by changing lateral CHO and NH to H for automated docking.