{'fraction_non_alphanumeric': 0.040325784664282875, 'fraction_numerical': 0.05337041451463383, 'mean_word_length': 4.660607196401799, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A conserved residue in the P2X4 receptor has a nonconserved function in ATP recognition

Highly conserved amino acids are generally anticipated to have similar functions across a protein superfamily, including that of the P2X ion channels, which are gated by extracellular ATP. However, whether and how these functions are conserved becomes less clear when neighboring amino acids are not conserved. Here, we investigate one such case, focused on the highly conserved residue from P2X4, E118 (rat P2X4 numbering, rP2X4), a P2X subtype associated with human neuropathic pain. When we compared the crystal structures of P2X4 with those of other P2X subtypes, including P2X3, P2X7, and AmP2X, we observed a slightly altered side-chain orientation of E118. We used protein chimeras, double-mutant cycle analysis, and molecular modeling to reveal that E118 forms specific contacts with amino acids in the “beak” region, which facilitates ATP binding to rP2X4. These contacts are not present in other subtypes because of sequence variance in the beak region, resulting in decoupling of this conserved residue from ATP recognition and/or channel gating of P2X receptors. Our study provides an example of a conserved residue with a specific role in functional proteins enabled by adjacent nonconserved residues. The unique role established by the E118-beak region contact provides a blueprint for the development of subtype-specific inhibitors of P2X4.

Highly conserved amino acids are generally anticipated to have similar functions across a protein superfamily, including that of the P2X ion channels, which are gated by extracellular ATP. However, whether and how these functions are conserved becomes less clear when neighboring amino acids are not conserved. Here, we investigate one such case, focused on the highly conserved residue from P2X4, E118 (rat P2X4 numbering, rP2X4), a P2X subtype associated with human neuropathic pain. When we compared the crystal structures of P2X4 with those of other P2X subtypes, including P2X3, P2X7, and AmP2X, we observed a slightly altered side-chain orientation of E118. We used protein chimeras, double-mutant cycle analysis, and molecular modeling to reveal that E118 forms specific contacts with amino acids in the "beak" region, which facilitates ATP binding to rP2X4. These contacts are not present in other subtypes because of sequence variance in the beak region, resulting in decoupling of this conserved residue from ATP recognition and/or channel gating of P2X receptors. Our study provides an example of a conserved residue with a specific role in functional proteins enabled by adjacent nonconserved residues. The unique role established by the E118beak region contact provides a blueprint for the development of subtype-specific inhibitors of P2X4.
The side-chain orientation of a given residue is closely related to protein functions such as ligand recognition (1) and conformational transitions (2) of these proteins during allosteric processes. Fine interpretation of the relationship between residue function and side-chain orientation at atomic level will advance our understanding of protein function (3) and provide more accurate structural information for future drug design (4,5). In general, the relationship between 3D structural information and function is more conserved in a protein superfamily than that of primary amino acid sequence (6,7). And it has also been suggested that in a given functional protein superfamily, conserved amino acids have similar functions among different members (8,9). However, if there are nonconserved residues located around the highly conserved residues, how do these surrounding residues affect the function of those conserved residues? Do the conserved and nonconserved residues encode different structural information together? Although there has been a great deal of studies aimed at understanding how do primary sequences determine the 3D structure and function of proteins, studying the relationship between sequence, structure, and function (RSSF) remains a daunting task. However, in recent years, a number of different crystal or cryo-EM structures have been identified (10)(11)(12), including structures from different subtypes and species, or at different states (apo or open states), such as the P2X receptors (13)(14)(15)(16)(17)(18)(19)(20)(21), providing invaluable information used to study the RSSF of a particular protein superfamily.
In this study, we compared the conformation/side-chain orientation of some conserved residues in several recently structurally identified P2X subtypes and found that a highly conserved residue E121 (equal to E118 in rP2X4) in the head domain of zfP2X4 has a different side-chain orientation from the other subtypes. Through a combination of molecular dynamics (MD) simulations, chimera construction, mutagenesis, and electrophysiological studies, we suggest that the unique side-chain orientation of E121 provides additional contacts between this residue and amino acids in the beak region, thus establishing a special role in determining the apparent affinity of ATP in rP2X4 but not in other P2X subtypes.

Results
E121 is a highly conserved residue in zfP2X4 with a different side-chain orientation compared with other P2X subtypes To reassess the role of conserved residues in P2X channel function, we compared the side-chain orientations of all conserved residues in several P2X subtypes whose 3D structures have been determined by X-ray crystallography, including zfP2X4 (Protein Data Bank [PDB] IDs: 4DW0 and 4DW1), hP2X3 (PDB IDs: 5SVJ and 5SVK), AmP2X (PDB ID: 5F1C), ckP2X7 (PDB ID: 5XW6), and pdP2X7 (PDB ID: 5U1L). The cryo-EM structures of rP2X7 were not included into this comprehensive comparison since we could not exclude possible influence of two different methods (singleparticle cryo-EM and X-ray crystallography) in determining residue side-chain orientation. Through a comprehensive structural superimposition, we found that three highly conserved residues R298, K316, and E121 (zfP2X4 numbering) have different side-chain orientations in different P2X subtypes (Table 1). R298 and K316 produce direct contacts with the agonist ATP in the open structure of zfP2X4 (Fig. 1A), and it is therefore reasonable to infer that the different side-chain orientation of these two residues may contribute to distinct recognition patterns for free ATP 4− or MgATP 2− and/or distinct ATP affinities in various P2X subtypes (36). E121 is located far away from the ATP-binding site in the head domain of zfP2X4 (Fig. 1, A and B), and its role in the channel function of the P2X receptors remains unclear.
To test whether the conformation of E121 in zfP2X4 and equivalent residues in other P2X subtypes are stable in a physiological condition, a series of 0.5-μs MD simulations were performed on aforementioned crystal or cryo-EM structures. These simulations showed that during simulations of zfP2X4, E121 always maintained an extended conformation (χ Cα-Cβ-Cγ-Cδ = 180 or −180 ) both in the resting and open states, whereas the same residues in other structures remained in the bent conformation (χ Cα-Cβ-Cγ-Cδ = 60 -80 ) (Fig. 2). The χ Cα-Cβ-Cγ-Cδ of E121 also remained at ±180 and 60 to 80 when MD simulations were performed based on cryo-EM structures of rP2X7 at the resting and open states, respectively, indicating that the MD simulation did not change the initial χ Cα-Cβ-Cγ-Cδ of E121 of rP2X7 (Fig. 1H). Thus, the sidechain orientation of these equivalent residues is not a transient "snapshot" in the structure determination but is also relatively stable in the structure fluctuations.
Alanine substitutions of E121 in zfP2X4 and E118 in rP2X4 significantly alter the apparent affinity of ATP Mutagenesis performed in these residues showed that zfP2X4 E121A resulted in an approximately 15-fold increase Table 1 The side-chain dihedral angles of a few residues in various crystal structures of P2X receptors  Special role of E118 in the apparent affinity of P2X4 in the apparent affinity of ATP (concentration yielding half-maximal activation, EC 50 = 1426 ± 733 and 89.6 ± 18.1 μM for zfP2X4 WT and zfP2X4 E121A , respectively, p < 0.05;    Fig. 3G; p > 0.05). Notably, the ATP-induced maximum currents were not affected by these alanine substitutions (Fig. 3, J and K). These results suggest that the unique side-chain orientation of E121 in the zfP2X4 receptor has a special role in determining the apparent affinity for ATP. The underlying mechanism of the altered apparent ATP affinity in zfP2X4 E121A should be further investigated; however, previous studies have shown that zfP2X4 is poorly expressed in mammalian cells and insensitive to extracellular ATP (EC 50 > 1 mM) (13,14), making it difficult to study the mechanism by using comprehensive mutagenesis together with electrophysiological recordings. Therefore, the following chimera constructions, double-mutant cycle analysis, and interaction-function correlation analysis were carried out in rP2X4 rather than zfP2X4 receptors. In addition, to obtain more accurate information for rP2X4 studies, we built homology models of rP2X4 and performed MD simulations of rP2X4 WT and its mutants.
Similarly, alanine substitution on the equivalent residue E118A of rP2X4 significantly changed the apparent ATP affinity of rP2X4 (7.65 ± 1.56 and 45.9 ± 9.1 μM for rP2X4 WT and rP2X4 E118A , p < 0.05, respectively; Fig. 3B). As controls, the corresponding mutants in hP2X1 and hP2X2 had no effect on the apparent affinity of ATP (EC 50 = 0.86 ± 0.16 and 0.92 ± 0.17 μM for hP2X1 WT and hP2X1 E119A , p > 0.05; EC 50 = 13.16 ± 1.97 and 12.92 ± 1.67 μM for hP2X2 WT and hP2X2 E127A , p > 0.05; Fig. 3, H and I). Only negligible ATP current could be detected in cells expressing homomeric mammalian P2X5 and P2X6 receptors, and therefore, effects of the same mutants on these two P2X subtypes were not measured. Thus, the apparent affinity of ATP was only affected by the substitution of alanine in E118 or E121 in P2X4, but not in all other P2X receptor subtypes tested, although this glutamate is highly conserved.
Additional mutants in E118 were constructed, and their effects on the affinity of rP2X4 were measured by a sequential ATP application protocol. The ratio (R) of currents (I 1 mM /I 8 μM ) induced by 8 μM (EC 50 ) and 1 mM (saturated) ATP was used to simplify the effect of various mutants on the apparent affinity of rP2X4 (Fig. 4A). rP2X4 WT had an R value of 2.16 ± 0.22, with Special role of E118 in the apparent affinity of P2X4 higher R values indicating higher EC 50 values for the mutants (Fig. 4, A and B). Any small change in the polarity or size at position 118 significantly impaired the apparent ATP sensitivity of rP2X4, for example, by replacing the Glu by Asp (I 1 mM /I 8 μM = 8.20 ± 2.20), Lys (11.61 ± 7.12), Ile (15.16 ± 6.09), Asn (9.44 ± 2.54), and Gln (15.98 ± 7.60). These results imply that E118 exerts its role through a strict interaction with adjacent residues.
The sequence differences of the beak region result in higher and lower EC 50 values for zfP2X4 E121A and rP2X4 E118A compared with WT channels, respectively Although zfP2X4 E121A and rP2X4 E118A could significantly change the apparent affinity of ATP, zfP2X4 E121A decreased, whereas rP2X4 E118A increased the EC 50 values of ATP. This difference should be further addressed. Chimeras zfP2X4 r-beak (with residues 136-148 replacements, zfP2X4 numbering) and rP2X4 zf-beak (with residues 133-145 replacements, rP2X4 numbering) were built by swapping their beak region sequences (Fig. 5, A and B). The zfP2X4 r-beak chimera led 10-fold increases in the apparent affinity of ATP (EC 50 = 121.9 ± 11.5 μM and 1.43 ± 0.73 mM for zfP2X4 r-beak and WT, respectively; Fig. 5C). As control, the rP2X4 zf-beak chimera resulted in an approximately 40-fold decrease in the apparent affinity of ATP (EC 50 = 232.3 ± 50.3 and 7.65 ± 1.56 μM for rP2X4 zf-beak and WT, respectively; Fig. 5D). However, there was no significant difference in ATP-induced maximum currents between WT and rP2X4 mutants (Fig. 5E). These two chimeras indicated that the sequence of the beak region correlates well with the apparent affinity of ATP in P2X4 receptors.
The role of E118 in the apparent affinity is associated with S141 and S142, two nonconserved amino acids in the beak region of rP2X4 We further studied how E118, together with the adjacent beak region, determines the apparent ATP affinity of rP2X4. During MD simulations of apo P2X receptors, the head domain behaves with a tendency to spontaneously move downward and a closure of ATP-binding site jaw, at least for P2X2, P2X3, and P2X4 (38)(39)(40)(41)(42)(43). Here, rP2X4 E118A also exhibits similar head domain movement during MD simulations (Fig. 6A), indicating that E118A does not compromise the inherent dynamics of the head domain.
We then explored the possibility that polar contacts between E118 and adjacent residues may contribute the role of E118 in ATP recognition of P2X4 receptors. At the resting state, the polar residues within 5 Å of E118 are S124 and D138 (Fig. 6B). The alanine substitution of these two residues had no Special role of E118 in the apparent affinity of P2X4 significant effect on R value of rP2X4 (Fig. 6F). At the open state, polar residues within 5 Å are S124 and S141 (Fig. 6C). In 0.5-μs MD simulations, we could observe stable hydrogen bond (H-bond) contacts between E118 and S141/S142 (Fig. 6D). The fact that S141A/S142A, as well as E118A/ S141A/S142A, led to significantly changed R values of rP2X4 (I 1 mM /I 8 μM = 9.10 ± 5.13 and 7.93 ± 2.03 for S141A/S142A and E118A/S141A/S142A, respectively; Fig. 6F) suggested that the E118 … S141/S142 interaction observed at the open state is essential for the apparent affinity of ATP in rP2X4 receptors.
In addition, we introduced tryptophan into the V161, a residue located near the E118 and S141/S142. Because V161 is located in the rigid β-sheet (β7) of the head domain, V161W will certainly block the S141/S142 … E118 interactions during channel opening of rP2X4 (Fig. 6E). Indeed, the R value of V161W increased to 13.7 ± 2.4 (Fig. 6F), indicating a significant decrease in the apparent affinity of ATP in rP2X4.

H140 mediates the correlation of E118 with ATP recognition of rP2X4
H140 is located in the lower beak region and may be involved in ATP recognition and channel gating of rP2X4 ( Fig. 7A; right) (14,15). During 0.5-μs MD simulations of rP2X4 WT with bound ATP, a tightening of the cleft between the dorsal fin and the rostral/head region and a shortening of the distance between ATP and H140 can be observed (Fig. 7B). In contrast, A118 swinged upward with H140 together and away from the bound ATP during MD simulations of rP2X4 E118A (Fig. 7, C-E). This altered allostery will unwind the interaction between ATP and H140 and the subsequent coordinated movement of the beak and dorsal fin domains, ultimately affecting the ATP recognition and channel gating of P2X4.
Double-mutant cycle analysis was also applied to confirm the coupling between E118 and H140. Single mutations E118A and H140A and double-mutant E118A/H140A resulted in increased EC 50 values of ATP to 45.9 ± 9.1, 53.3 ± 9.1, and 55.1 ± 14.4 μM, respectively (Fig. 7F). The coupling free energy ΔΔG of E118 and H140 was calculated to be −1.02 kcal/mol (Fig. 7G), indicating that H140 and E118 are interdependent in contributing to ATP recognition and channel gating of P2X4.
In other P2X subtypes, there is no correlation between the highly conserved residue and ATP recognition Finally, we explored why this conserved residue has no effect on the apparent ATP affinity of other P2X subtypes (Fig. 8) by using MD simulations. In 0.5-μs MD simulations of hP2X3 WT or hP2X3 E109A at the open state, P128, the residue corresponding to H140 of rP2X4, remained away from ATP ranging from 8.4 to 14 Å (the distance between the C α of P128 and N9 atom of ATP; Fig. 8, A-C). Although the P128 moved closer to ATP in hP2X3 E109A , it is still impossible to establish a direct contact between the beak region and ATP at this distance (Fig. 8, A and C). Thus, sequence differences in the beak region and/or the structure of this region in hP2X3 might compromise the role of E109 in the apparent ATP affinity of hP2X3. Similarly, MD simulations of AmP2X and AmP2-X E129A reveal 8 Å distance between the beak region and ATP (Fig. 8, D-F). Therefore, it is possible that sequence differences of the beak regions in other P2X subtypes interrupt the coupling of this conserved residue, like E118 in rP2X4, to the final ATP recognition of the P2X receptor.

Discussion
In RSSF studies, highly conserved residues often thought to have a conserved role in protein architecture and function. In our study, E121 (zfP2X4 numbering) is a highly conserved residue, but its side-chain orientation differs from that of other P2X subtypes. Accordingly, the zfP2X4 E121A mutant significantly altered the apparent affinity of ATP rather than the corresponding mutants in other P2X subtypes. The specific role of this unique side-chain orientation-determined conserved residue in ATP recognition is due to the coupling established by this conserved residue and adjacent nonconserved residues to the ATP recognition of P2X4. Our findings suggest that a highly conserved residue in the primary amino acid sequence could "own" the nonconserved role in proteins, and that the conserved residue is able to establish a special function with adjacent nonconserved residues.
The ability of rP2X4 E118 in determining the ATP recognition and/or channel gating is possibly through the "E118 … S141/S142 … H 140 axis," especially for the H-bond contact between E118 and S141/S142 complexes. The following evidence can support this idea. First, several E118 mutants suggest that the strong polarity and proper side-chain size of this residue are important (Fig. 4), and stable H-bond contacts between E118 and S141/S142 could be detected during MD Figure 6. E118 interacts with S141/S142 in the beak region of rP2X4 receptors. A, downward motion of the head domain during 0.5-μs molecular dynamics (MD) simulation of rP2X4 E118A in the resting state. B, polar amino acids of the head domain near E118 (within 5 Å) in rP2X4 in the resting (B) and open (C) states. D, the H-bond interactions between E118 and S141/S142 during 0.5-μs MD simulation of open structure of rP2X4. The trajectory was sampled every 200 ps. E, superimposed homology models of V161 and its tryptophan substitution (W161) to show the position of this residue (at rigid β7) near to S141, S142, and E118. F, summary of R values of WT and rP2X4 mutants. n = 15, 7, 8, 7, 7, 5, 5, 6, and 5, from left to right, respectively. One-way ANOVA with Bonferroni post-test (F(8, 56) = 19.98, p < 0.0001). The data on WT and E118A are the same as those displayed for Figure 5E. G, concentrationresponse curves of ATP in rP2X4 WT , rP2X4 E118A , rP2X4 S141A/S142A , and rP2X4 E118A/S141A/S142A . The solid lines are fitted by Hill equation (mean ± SEM, n = 3-5, EC 50 = 51.5 ± 16.1 and 59.6 ± 13.2 μM, n H = 0.92 ± 0.14 and 1.07 ± 0.14 for S141A/S142A and S141A/S142A/E118A, respectively). The data of rP2X4 WT and rP2X4 E118A are the same as those displayed for Figure 4B and are shown here as controls. H, double-mutant cycle analysis between E118 and S141/S142. rP2X4, rat (Rattus norvegicus) P2X4.
Special role of E118 in the apparent affinity of P2X4 simulations (Fig. 6D). Second, the carboxyl group of E118 is oriented toward S141/S142 at the open state, which would facilitate the H-bond contacts among these residues (Fig. 6C). Third, the R values of S141A or S142A are slightly higher than that of rP2X4 WT, but the R value of double-mutant S141A/ S142A increases to a level equal to that of E118A, suggesting that both S141 and S142 can make contacts with E118 (Fig. 6F). This idea is further confirmed by the double-mutant analysis of E118 and S141/S142 (Fig. 6H). Fourth, during MD simulations of rP2X4 E118A , loss of E118 … S141/S142 contact was associated with upward movement of the beak region, which resulted in the displacement of H140 from the ATP  Special role of E118 in the apparent affinity of P2X4 (Fig. 7C). Double-mutant cycle analysis of E118 and H140 (ΔΔG = −1.02 kcal/mol) also demonstrated the functional coupling of these two residues (Fig. 7G). Thus, the established "E118 … S141/S142 … H140 axis" is the mechanism for the special role of E118 in determining the ATP recognition.
Our findings also support the idea that the nonconserved residue is able to significantly alter the conformation and function of highly conserved residues, especially when these conserved residues are located in a loop region. Based on the sequence alignment, we could find that although zfP2X4 E121 is highly conserved but the other residues in the beak region (a loop structure; Fig. 1C) are nonconserved throughout the P2X family. Accordingly, in the crystal structure of zfP2X4, the side chain of E121 points outside the head domain, showing an extended conformation with its carboxyl group facing to the lower beak region. However, in other P2X subtypes, the side chain of the equivalent residue bends inward to the head domain, making no contact with the lower beak region. Interestingly, differences in the beak region sequences of zfP2X4 and rP2X4, from the same subtype but different species, persist, which result in different changes in their EC 50 values under alanine substitutions in E118 and E121, respectively. However, the mere swap of the beak region failed to restore the change induced by E118A based on the rP2X4 zf-beak . It is possible that rP2X4 zf-beak adopts a different conformation that decouples the link between E118 and the ATP-binding site. What is the physiological function of this "designed" special role of E118/ E121 in the evolution of P2X4 receptors that should be further studied. We speculate that the unique role of E118/E121 in P2X4 receptors provides a different way for endogenous modulation in certain species. For example, previous studies have shown that some trace elements could regulate the function of P2X4, and that E118, H140, and their adjacent amino acids are likely to be involved in metal binding (44).
Finally, our study provides a prerequisite for the further development of new high-selective small molecule targeting P2X receptors. The upper lumen, formed by the head domains and the beak region, is an allosteric site, and several compounds/ metal ions have been shown to stimulate this site (45)(46)(47). Understanding the relationship between the side-chain orientation of E118 within this site and the ATP-recognition P2X receptors will facilitate the design and optimization of potent and highly selective new analogs based on the known compounds. Furthermore, there is no compound that has been shown to act directly on the P2X4 receptor at this site, and we can screen for new P2X4-selective lead compounds.

Cell culture and mutagenesis
Human embryonic kidney 293 (HEK293) cells were cultured in 5% CO 2 /95% air at 37 C in Dulbecco's modified Eagle's medium (Corning) supplemented with 10% fetal bovine serum (PAN), 1% Penicillin-Streptomycin (Gibco), and 1% Glutamax (Gibco). Calcium phosphate transfection was used for delivering plasmid DNA into HEK293 cells cultured on the coverslips placed in the 35-mm dishes. About 2.5 μg plasmid DNA was mixed directly with 100 μl 0.25 M CaCl 2 solution, and then, the mixed solution was slowly added into 2× Hepesbuffered saline solution with tip agitating to ensure the formation of fine calcium phosphate-DNA precipitate. Hepesbuffered saline solution (2×) contains (in millimolar) 140 NaCl, 1.5 Na 2 HPO 4 , 50 Hepes, and the pH was adjusted to 6.96. All mutants were constructed by KOD-Plus mutagenesis kit (Toyobo) and confirmed by DNA sequencing (Beijing Genomics Institute). rP2X4 and zfP2X4 plasmids were a gift from Drs Alan North and Lin-Hua Jiang, and AmP2X plasmid was a generous gift from Dr Motoyuki Hattori; the complementary DNAs of hP2X1 and hP2X7 were synthesized and subcloned into pcDNA3.1 by Shanghai Genechem; complementary DNAs for hP2X2, ckP2X7, and pdP2X7 were synthesized and subcloned to pcDNA3.0 by Beijing Genomics Institute; the plasmid for hP2X3 was purchased from Open Biosystems.

Electrophysiological recordings
Whole-cell patch clamp recordings were performed at room temperature (25 ± 2 C), using an Axopatch 200B amplifier (Molecular Devices). Current signals were sampled at 10 kHz and filtered at 2 kHz and analyzed by pClamp 10 (Molecular Devices). HEK293 cells were superfused with standard extracellular solution containing (in millimolar) 150 NaCl, 5 KCl, 10 glucose, 10 Hepes, 2 CaCl 2 , and 1 MgCl 2 (pH 7.4). The pipettes were pulled from glass capillaries through two-stage puller (Narishige PC-10) and filled in a general intracellular solution containing (in millimolar) 120 KCl, 30 NaCl, 0.5 CaCl 2 , 1 MgCl 2 , 10 Hepes, 5 EGTA (pH 7.2), and a resistance varied between 2 and 4 MΩ. During gap-free recordings, the membrane potential was held at −60 mV, and the ATP solution was applied with the Y-tube. As we previously described (19,48), nystatin-perforated patch clamp was carried out to reduce the dialysis of intracellular constituents and current rundown. Nystatin-perforated intracellular solution contains (in millimolar) 75 K 2 SO 4 , 55 KCl, 5 MgSO 4 , and 10 Hepes (pH 7.2). Nystatin was purchased from Sangon Biotech, and all other drugs were purchased from Sigma-Aldrich. Maximal currents or R ratios were obtained by whole-cell patching, whereas concentration-response curve recordings were performed using perforated patch clamp to avoid channel desensitization.

Homology modeling and MD simulations
The homology models of rP2X4 WT and rP2X4 E118A were constructed based on the crystal structures of zfP2X4 (PDB IDs: 4DW0 and 4DW1) by using Modeler as we described previously (42,49). The alignment of the target sequences was manually adjusted to match published alignments (13); the obtained model was checked and validated by ProCheck (50); the energy-minimized structures of P2X were used as the initial structures for MD simulations. A large 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine bilayer, available in System Builder of DESMOND (Schrödinger) (51), was used to generate a suitable membrane system in which the transmembrane domain of the P2X could be embedded. The P2X/ Special role of E118 in the apparent affinity of P2X4 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine system was dissolved in simple point charge water molecules. Counter ions were then added to compensate for the net negative charge of the system. NaCl (150 mM) was added into the simulation box that represents background salt at physiological condition. The default relaxation protocol of DES-MOND for each system was applied before the simulations were run: (1) 100 ps simulations in the canonical (constant moles, volume and temperature) ensemble with Brownian dynamics using a temperature of 10 K with solute heavy atoms restrained; (2) 12 ps simulations in the canonical (constant moles, volume and temperature) ensemble using Berendsen thermostat at 10 K with small-time step and solute heavy atoms restrained; (3) 12 ps simulations in the isothermal-isobaric (constant moles, pressure and temperature) ensemble using the Berendsen thermostat and barostat at 10 K and 1 atm with solute heavy atoms restrained; (4) 12 ps simulations in the isothermal-isobaric (constant moles, pressure and temperature) ensemble using the Berendsen thermostat and barostat at 300 K and 1 atm with solute heavy atoms restrained; (5) 24 ps simulations in the isothermalisobaric (constant moles, pressure and temperature) ensemble using the Berendsen thermostat and barostat at 300 K and 1 atm with no restraints. After equilibration, the MD simulations were performed for 500 ns. Long-range electrostatic interactions were computed using the smooth particle mesh Ewald method. The integration time step used was 1 fs, and the coordinate trajectories were saved every 200 ps. All simulations were performed by using DESMOND with a constant number of particles, pressure (1 bar), and temperature (300 K) and periodic boundary conditions by using Nose-Hoover chain thermostat. Proteins, ions, lipids, and the simple point charge waters were assigned to all-atom OPLS_2005 force field (52,53). Simulations were run on DELL T7910 graphic working station (with NVIDA Tesla K40C-GPU). Preparation, analysis, and visualization were performed on a 12-CPU CORE DELL T7500 graphic working station.

Data analysis
Bonferroni's multiple comparisons test (ANOVA) was used when there were more than two groups. *p < 0.05 or **p < 0.01 was considered significant. Concentration-response curves were fitted to the equation as I/I max = 1/(1 + (EC 50 / [ATP]) n ), where I is the normalized current at a given concentration of ATP, I max is the maximum normalized current, EC 50 is ATP dose for half-maximal effect, and n is the Hill coefficient. Statistical comparison between EC 50 values was performed using extra sum-of-squares F tests. The change in Gibbs free energy upon each mutant was described by the following equation: ΔG ¼ − RT ln EC 50 of WT EC 50 of mut , where the constants R and T are 1.987 cal/mol/K and 293 K, respectively (54). The interaction free energy in the double-mutant cycles was calculated from ΔΔG interaction ¼ ΔG mut1þ2 − ðΔG mut1 þ ΔG mut2 Þ. If the change in the free energy of the double mutant (ΔG mut1+2 ) equals the sum of the free energy changes of the corresponding single mutations (ΔG mut1 and ΔG mut2 ), which means that the value of ΔΔG interaction is close to zero, the two residues or clusters will not interact. In contrast, if the two residues interact, the change in the free energy of the double mutant will differ from the sum of the changes of the corresponding single mutations, and ΔΔG interaction should be larger than ±0.35 kcal/mol (37) that has been reported for noninteracting residues.

Data availability
All data are contained within the article.