Multiscale molecular dynamics simulations predict arachidonic acid binding sites in human ASIC1a and ASIC3 transmembrane domains

Acid-sensing ion channels (ASICs) are trimeric proton-sensitive cation channels found in the central and peripheral nervous system. Identified as part of the family of degenerin/epithelial sodium channels, ASICs play an especially important role in inflammatory pain pathways (Cheng et al., 2018). Several ASIC subunit types exist, including the major types ASIC1, ASIC2, ASIC3, and ASIC4 and their variants. Different ASIC subtypes have distinct function and nervous system localization (Kweon and Suh, 2013) and may assemble as homo- or heterotrimers. All ASIC types share a similar structure (Fig. 1, A and B). Each subunit contains a large extracellular domain (ECD), two transmembrane helices (TM1 and TM2) forming the transmembrane domain (Jasti et al., 2007; Fig. 1, A and B), and a recently resolved re-entrant intracellular domain (ICD; Yoder and Gouaux, 2020). ASICs are activated at acidic pH, presumably through protonation of acidic residues in the ECD, leading to rearrangements and subsequent opening of the ion pore in the transmembrane domain (Jasti et al, 2007; Yoder et al, 2018). In this way, ASICs act as important sensors of tissue acidosis in the central and peripheral nervous systems.

The neuronal membrane environment plays an important role in the function of ion channels. For example, pentameric ligand-gated ion channels, such as the nicotinic acetylcholine receptors, are functionally sensitive to the lipid composition of their surrounding membrane (Baenziger et al., 2000; daCosta et al., 2013; Ananchenko et al., 2022). It has also been shown that several voltage-gated ion channels are functionally affected by polyunsaturated fatty acids (PUFAs), with the presence of PUFAs altering the voltage dependence and maximal current of these channels (Elinder and Liin, 2017; Tong et al., 2019). PUFAs such as docosahexaenoic acid (DHA) and arachidonic acid (AA; Fig. 1 C) are abundant in neuronal membranes in their esterified form as phospholipid tails (Salem et al., 1996). The concentration of PUFAs as free fatty acids has been estimated to be between 1 and 50 µM in plasma (De Caterina et al., 2000; Fraser et al., 2003); however, the exact concentration in the neuronal plasma membrane, and its physiological relevance, is unknown. It is suggested that the PUFA concentration, in particular DHA and AA, in neuronal membranes can be affected by factors such as diet or sex (Chen et al., 2020). Physiological conditions can also affect the free fatty acid abundance of AA, e.g., inflammation increases the AA concentration (Marra et al., 2016). The importance of PUFAs in neuronal function, in particular through their effects on ion channels, is an evolving field.

Recent studies have illustrated that ASICs can be modulated by fatty acids, most prominently by AA (Smith et al., 2007; Marra et al., 2016). Modulation by AA has been shown to occur in several ASIC subtypes (Smith et al., 2007). AA has been suggested to shift the pH dependence of the channel, leading to increased activation at higher pH, i.e., smaller proton concentration (Allen and Attwell, 2002; Smith et al., 2007; Klipp and Bankston, 2022). Evidence suggests that this occurs through a direct binding mechanism, rather than through changes in the bilayer properties (Smith et al., 2007). Interestingly, the ASIC3 subtype can be activated by AA at physiological pH (Marra et al., 2016), potentially due to a more extreme shift in pH₅₀ or a higher initial pH₅₀, whereas ASIC1a does not open without a drop in pH, despite the presence of AA (Marra et al., 2016). pH₅₀ is the pH that yields half maximal response. However, other studies suggest that AA acts in the same manner on both the ASIC1 and ASIC3 subtypes (Klipp and Bankston, 2022). Other PUFAs also have a modulatory effect on ASICs, the potency of which is suggested to correlate with tail length and saturation (Smith et al., 2007; Klipp and Bankston, 2022). AA itself is an important inflammatory marker in the nervous system, with increased levels of AA seen in exudates of patients with rheumatoid arthritis (Marra et al., 2016) and psoriasis (Hammarström et al., 1975).

Until very recently, the binding site, if any, for AA on ASICs has remained unknown. A recent study has suggested that R64_hASIC1 (R63_hASIC3, Fig. 1 B) could play an important role in binding of, and potentiation by, AA in both ASIC1 and ASIC3 (Klipp and Bankston, 2022). This functional study gives promising evidence of AA modulation of ASICs through structural effects; however, evidence of direct AA binding to this residue is yet to be shown. In addition, recent structures of chicken ASIC1 solved in styrene maleic anhydride nanodiscs reveal electron densities attributed to lipids (Yoder and Gouaux, 2020). However, the identity of these lipids and their functional relevance are unknown. Overall, the exact mechanism by which AA modulates or even activates ASICs remains unclear.

In this study, we identify a putative binding region for AA in the extracellular half of the TMD. This binding region is similar between the human ASIC3 (hASIC3) and hASIC1a subtypes, consistent with the general modulatory effect of AA on several ASIC isoforms. We observed subtle differences in the AA coordination pattern between the two isoforms, with the largest differences seen in residues which are conserved within their respective ASIC subtypes. Importantly, we suggest that AA can form stable, direct interactions with ASICs. This is, for example, shown by the observation from simulations that the hASIC3 residue R65_hASIC3, located near the transmembrane domain (TMD) fenestrations, forms a stable salt bridge with the carboxylic acid headgroup of AA. To expand upon recent studies such as Klipp and Bankston (2022), we speculate that modulation of ASIC by AA does not occur solely through effects on the structure of TMD fenestrations directly, but rather through a more complex mechanism involving changes in contacts and conformation between TM1 and TM2. We proposed that several coordinating interactions, rather than a single binding residue, may contribute to the overall ASIC-AA binding affinities and modulatory mechanism.

Homology models of hASIC1a and hASIC3 were constructed from crystal structures of chicken ASIC1 in both the open and resting states (PDB: 4NTX, 5WKU; Baconguis et al, 2014; Yoder and Gouaux, 2020) using Modeller (Eswar et al., 2006). The following sequences from the Uniprot database (Bateman et al., 2021) were used: Q1XA74 for chicken ASIC1, P78348 for hASIC1a (89.35% sequence identity with chicken ASIC1), and Q9UHC3 for hASIC3 (49.71% sequence identity with chicken ASIC1). A disordered extracellular loop (residues 291–304 in hASIC3 numbering) was removed from the hASIC3 structure due to poor modeling and to reduce the solvated box size. Rather, our models incorporate a direct link between residues 290 and 305 in hASIC3. This extended loop in hASIC3 relative to ASIC1 is found in the ECD and thus not believed to affect the potential interaction of PUFAs with the TMD. 100 structures were made for each homology model and the structure with the lowest molpdf score (Eswar et al., 2006) was chosen for simulations. The resulting all-atom homology models were converted to coarse-grain (CG) representation using the martinize.py script and the MARTINI2.2 CG forcefield (Marrink et al., 2007). The secondary, tertiary, and quaternary structures of the protein models were restrained using an elastic network with a force constant of 1,000 kJ/mol and a cut-off radius (R_c) value of 1.05 nm, constructed by the martinize.py script. These values were chosen to stabilize the specific functional states and to maintain the root mean square deviation of each CG structure within 4 Å. The insane.py script (Wassenaar et al., 2015) was used to embed the CG models in a bilayer in a 140 × 140 × 150 ų cubic box, solvated in MARTINI CG water with 0.15 mM NaCl, with additional neutralizing Na⁺ ions added to systems containing AA. Each system contained ∼25,000 CG beads. Insane.py randomly places lipids in the intra- and extracellular leaflets. Therefore, repeats of simulations were set up independently to randomize lipid placement.

Initial simulations, using CG structures of the open and resting states, were set up with a membrane containing an 8:1:1 ratio of palmitoyl-oleoyl-phosphatidylcholine (POPC), lysophosphatidylcholine (LPC) and AA. A 10 µM solution of AA has been used in electrophysiology recordings to invoke ASIC3 currents at physiological pH (Marra et al., 2016), a concentration which is not attainable in the membrane patch size used in our simulations. We initially set up systems containing six molecules each of AA and LPC in each leaflet (1% AA); however, we found this to be insufficient for sampling interactions over 30 µs. We therefore increased the number of modulatory lipids 10-fold to 60 molecules per leaflet (10% AA) to ensure adequate sampling of lipid–protein interactions. This AA concentration is well above that of a physiological membrane composition. However, it allowed us to probe a large number of contacts between AA and the ASIC protein. In this way, it provided us the ability to reproducibly measure lipid occupancy and interaction duration, without deforming the membrane. Additional control simulations were set up the same way but containing only POPC in the membrane. Headgroup density calculations were performed using the 8:1:1 POPC:LPC:AA membranes. No specific associations between LPC and ASIC models were seen in the initial simulations, so subsequent simulations were performed using a 9:1 ratio of POPC:AA.

Mutations were created in the atomistic homology models of hASIC1a and hASIC3 using the ChimeraX (Pettersen et al., 2021) rotamers tool, selecting rotamers to minimize clash. The resulting atomistic mutant structures were subsequently coarse-grained and embedded in a 90:10 POPC:AA membrane as described above. For all structures, all ionizable residues were kept in their standard state. ASICs are believed to be activated by protonation of a number of acidic residues in the ECD. In this case, however, where AA is proposed to be able to activate ASIC3 without a drop in pH, presumably by shifting pH₅₀, we did not want to bias our simulations toward the open state as we wanted to capture the “pure” AA effect. Therefore, we probed AA interactions with the channel in the absence of additional proton activation.

The GROMACS (van der Spoel et al., 2005; Abraham et al., 2015) 2019 simulation engine was used for all simulations. Energy minimization was performed for three rounds of 5,000 steps using steepest descent. Systems were equilibrated using the NPT ensemble (constant number of atoms and constant pressure and temperature) for a total of 4.75 ns with gradual lowering of protein position restraints from 1,000 to 50 kJmol⁻¹nm⁻² and lipid position restraints from 200 to 10 kJmol⁻¹nm⁻². Equilibration was performed with temperature held at 310 K using the v-rescale thermostat (Bussi et al., 2007) and pressure held at 1 bar with the Berendsen barostat (Berendsen et al., 1995). Production runs were simulated in triplicate for 30 µs each. Production runs were simulated with a timestep of 20 fs using an NPT ensemble with v-rescale set to 310 K and semi-isotropic Parinello-Rahman pressure coupling (Parrinello and Rahman, 1998) set to 1 bar.

Simulations with DHA and palmitic acid instead of AA were set up following the above protocol, using the open state CG structure, with a 9:1 ratio of POPC:FA where FA is the fatty acid of interest.

Frames from the last 100 ns of CG simulations were selected, for which AA was seen to interact with the protein at least two upper TMD sites. One frame was selected from each repeat of relevant systems. Frames were backmapped to atomistic resolution using the CG2AT tool (Vickery and Stansfeld, 2021). A higher resolution structure of the ASIC1 open state (PDB: 4NTW; Baconguis et al., 2014) was used to create an atomistic homology model as described in the homology modeling section above, and this model was subsequently used as the reference structure for alignment in the CG2AT protocol. One CG2AT system was created for each repeat of simulations containing the open and resting states of ASIC bound to at least two AA molecules at the extracellular TMD, solvated in TIP3P water. A total of 500 ns of simulation were performed for each starting system. AA parameters were obtained from the CHARMM-GUI lipid library (Lee et al., 2016).

Atomistic simulations were performed using the GROMACS 2019 simulation engine (van der Spoel et al., 2005; Abraham et al., 2015) and the CHARMM36 forcefield (Best et al., 2012; Huang and Mackerell, 2013). The resulting backmapped system was energy minimized for 5,000 steps or until a maximum force of 1,000 kJmol⁻¹nm⁻¹ on any atom was reached using the steepest descent algorithm. Equilibration was performed using the NPT ensemble for a total of 1 ns with temperature held at 310 K using the Berendsen thermostat, and pressure held at 1 bar using the Berendsen barostat. Production runs were performed for 500 ns in the NPT ensemble using the Nose-Hoover thermostat (Evans and Holian, 1998) with temperature set to 310 K and pressure held at 1 bar using the Parinello-Rahman barostat. The Verlet cut-off scheme (Verlet, 1967) was used throughout all steps with a force-switch modifier starting at 10 Å and a cutoff of 12 Å. The particle mesh Ewald (Darden et al., 1998) method was used for long-range electrostatics and a cutoff of 12 Å was used for short-range electrostatics. Covalent bonds including hydrogen atoms were constrained using the LINCS algorithm (Hess et al., 1997).

Lipid occupancies and residence times at individual residues were calculated using the PyLipID tool (Song et al., 2022), with 0.475 nm and 1.1 nm lower and upper cutoffs, respectively. In-house tcl scripts and VMD (Humphrey et al., 1996) were used for calculating density maps of lipid occupancies and membrane curvature. All other analysis was performed using in-house MDAnalysis (Michaud-Agrawal et al., 2011; Gowers et al., 2016) scripts. Figures were created in VMD (Humphrey et al., 1996) and Chimera X (Pettersen et al., 2021).

For analysis, the charged “COO” bead was used to define the headgroup of fatty acids, and the “PO4” bead was used as the headgroup phosphate of POPC.

The supplementary data includes information on the optimization of cut-offs for use with the PyLipID software (Fig. S1). In addition, supplemental distance and interaction measurements for all chains, as shown in the main text Fig. 6, are illustrated in Fig. S2. Video 1 shows a representative CG simulation where AA interactions with the relevant sites are sampled. Video 2 illustrates a 200 ns section of a 500 ns simulation with one example of an AA molecule interacting with R65_hASIC3, R68_hASIC3, and Y58_hASIC3. Video 3 illustrates another 200 ns section of a 500 ns simulation with a different example of an AA molecule interacting with R65_hASIC3, R68_hASIC3, and Y58_hASIC3.

In order to computationally investigate lipid interactions with ASICs, we employed the MARTINI forcefield and performed CG MD simulations (Marrink et al., 2007). By restraining the channel to its original conformation and allowing lipids to sample the transmembrane surface of the protein, we were able to sample sites that are specific to the structurally resolved conformation of the protein. At the time when simulations were initiated, the only solved ASIC structures available were those of chicken ASIC1 without a resolved ICD. We therefore used these available structures to create homology models of hASIC1a and hASIC3. Subsequently, additional ASIC structures were solved of hASIC1a in an inhibited state (Wu et al., 2021) and chicken ASIC1 with additional ICD regions resolved in the form of re-entrant loops (Yoder and Gouaux, 2020). No structure of the ASIC protein is yet available in the open state with resolved re-entrant loops. For these reasons, in order to specifically sample AA interactions with the open and resting states of the channel, we continued to use homology models of hASIC1 and hASIC3 created from the chicken ASIC1 structures without the ICD, as outlined in the Materials and methods section. Furthermore, due to only the incomplete structures being available at the time, we chose to focus on AA interactions with the extracellular half of the TMD. In a recent study, it was demonstrated that fatty acids act on the channel in the extracellular membrane leaflet (Klipp and Bankston, 2022), making it likely that key interacting residues are located towards the extracellular end of the TMD. It was shown in the more recent structures of chicken ASIC1 in the desensitized state that the conformation of the TMD helices near the intracellular side appears unaffected by the presence of re-entrant ICD helices, although the side chain orientations in this region are altered (Yoder and Gouaux, 2020). Given the minimal effect of the re-entrant loop on the overall structure, we therefore assumed that lipid interactions in the upper TMD region are not affected by the missing ICD re-entrant helices in our models.

AA has been postulated to affect ASIC function through direct binding to the channel, rather than through membrane deformation (Smith et al., 2007). We wanted to test whether AA interacts specifically with ASICs in silico. We first investigated whether hASIC1a and hASIC3 preferentially interacted with AA in a membrane containing 10% LPC, 80% POPC, and 10% AA. As outlined above, we focus on the extracellular side of the TMD due to uncertainty of the structure of the intracellular component of the TMD. From a total of three 30 µs simulations per system, the average lipid occupancy was calculated using either all 60 AA carboxylic acid headgroups, or 60 randomly selected POPC headgroups for comparison. The AA carboxylic acid headgroups showed a higher average occupancy near the open state hASIC3 transmembrane domain in comparison to the POPC headgroups (Fig. 2). AA also specifically localized around hASIC1, albeit at a slightly lower apparent density than for hASIC3. In both cases, the AA occupancy was especially high around the outer region of TM1, the outer transmembrane helix of ASIC (Fig. 1). Interestingly, the general region of POPC localization correlated with electron densities attributed to lipids in recent structures (Yoder and Gouaux, 2020), showing possible lipid tail density at the lipid-exposed surface of TM2 (Fig. 2). On the contrary, AA density localized at TM1, toward the periphery of the channel (density from simulations only). This supports the idea that AA exhibits specific interactions with ASICs, which differs from that of phospholipids such as POPC.

First, using our models of hASIC3 in the open and resting state, we set out to identify transmembrane sites which participate in interactions with AA. We first compared AA occupancies between the open and resting states of hASIC3. Occupancy is measured as a percentage of total simulation time in which any AA molecule is within range of a lower distance cutoff (4.75 Å) from any bead of a residue. Distance cutoffs were previously optimized by sampling lipid interactions with a variety of distance cutoffs as described for the PyLipID method (Song et al., 2022). A lower cutoff of 4.75 Å was sufficiently short to select for specific interactions between AA and ASICs, while also allowing for sufficient sampling of lipid binding sites (Fig. S1). Surprisingly, AA occupancy did not appreciably differ between the two hASIC3 states (Fig. 3 A). An example of AA sampling the hASIC3 open state TMD is illustrated in Video 1. For further comparison, we measured the average AA interaction duration time for the TMD residues (Fig. 3 B), with an interaction considered to begin when an AA molecule is within range of the lower cutoff, and continuing until the molecule exits the upper cut-off range of 1.1 Å. Notably, the hASIC3 open state model exhibited longer interactions at R63, a residue lining the lateral fenestrations (Fig. 1 B) and previously suggested to be involved in AA interactions with ASIC (Klipp and Bankston, 2022), as well as in the proton sensitivity of the channel (Chen et al., 2021).

We then set out to elucidate the structural determinants of specific AA interactions with the hASIC1 and hASIC3 transmembrane domains in the open state. Such determinants were expected to account for specific interaction of AA with ASICs and furthermore illustrate a level of conservation of a potential AA binding site between hASIC1 and hASIC3. To determine which residues of the TM domain AA preferentially interacts with, if any, we calculated the occupancy of AA at each residue in TM1 and TM2 in the same way as for Fig. 3 A. The occupancy was calculated as an average for three homomeric chains in three 30 µs simulations in both hASIC1 and hASIC3, respectively, thereby predicting which residues are potential interacting partners for AA. Such residues could constitute a potential binding site. We compared occupancies at each residue in the upper region of TM1 (Fig. 4 A) and TM2 (Fig. 4 B) between hASIC1 and hASIC3 to identify conserved interactions but also predict potential interactions that may account for the subtle differences in overall occupancy as seen above (Fig. 2).

There were no clear differences in occupancies at residues in TM2 (Fig. 4 B) between hASIC1 and hASIC3. Although certain residues in TM2 have high occupancies, their large variability between chains and simulations, as shown by the standard deviations, make this data unclear. Likely, this is due to the reduction in lipid accessibility of TM2 compared to TM1, resulting in less frequent lipid interactions and therefore less sampling. However, there were residues of particular interest in TM1. Occupancy was high at both R68_hASIC3 in hASIC3 and F69_hASIC1 in hASIC1—the occupancy at these residues was above 10% and thus higher than expected given the AA concentration of the membrane. Occupancies at Y58_hASIC3, R65_hASIC3, and Y67_hASIC3 were noticeably higher in hASIC3 than for the equivalent residues in hASIC1 (Fig. 4, A and C). We also observed higher occupancy at V64_hASIC3; however, we rationalized that this higher occupancy at valine may be a product of its adjacency to R65_hASIC3. On the other hand, the AA occupancy at Y71_hASIC1 in hASIC1 was approximately twofold higher than at F70_hASIC3 in hASIC3. Interestingly, despite AA modulating, and therefore interacting with, both ASIC subtypes, Y58_hASIC3 and R65_hASIC3 are conserved in the ASIC3 subtype, but not ASIC1, while Y71_hASIC1 is conserved in the ASIC1 subtype, but not ASIC3 (Fig. 5 A). We next wondered whether the two identified hASIC3 residues, Y58_hASIC3 and R65_hASIC3, could constitute a single AA binding site on hASIC3. A representative binding pose from PyLipID analysis (Fig. 5 B) supports the presence of such a site. Briefly, binding sites are predicted based on community analysis of interaction networks. The most sampled binding poses are then clustered and scored based on density. This method is described in detail in Song et al. (2022). These data may suggest that while both hASIC1a and hASIC3 interact with AA at the outer leaflet region of TM1, the specific binding motif may vary somewhat between the two subtypes.

Given that hASIC3 had slightly higher AA occupancy overall compared to hASIC1a, we backmapped selected simulation frames from our hASIC3 CG simulations to atomistic resolution for further simulations of 500 ns to determine whether the interactions observed in coarse-grained simulations remained stable. We observed further time-resolved AA interactions with R65_hASIC3 and Y58_hASIC3 (Videos 2 and 3), as well as notable salt bridge interactions with R68_hASCI3. Defining an interaction as a distance <4.5 Å between a non-hydrogen atom of AA and a non-hydrogen atom of a protein residue, we observed that an interaction is present >50% of the total trajectory length between R65_hASIC3 and an AA headgroup in some individual chains of ASIC. Taking the chains together, this interaction is present at 100% of the trajectory in at least one of three chains (Fig. 6 A and Fig. S2). Taking the initially bound AA molecules, we measured the distance between the carboxylic acid carbon atom and the C-ζ atom of the arginine that it was interacting with. Stable interactions (distance <4.5 Å) were observed in all chains, with stable interactions lasting as long as 300 ns (Fig. 6 B and Fig. S2). This interaction appears to be a salt bridge between the carboxylic acid of AA and the guanidinium group of R65_hASIC3 (Fig. 6 C). The AA carboxylic acid headgroup can also form occasional simultaneous salt bridges with R65_hASIC3 and R68_hASIC3 (Videos 2 and 3). These observations from our atomistic simulations support that ASICs can stably bind AA.

To investigate whether the two identified residues with high AA occupancy (Y58_hASIC3, R65_hASIC3), and the additional coordinating residue R68_hASIC3, influence the overall interaction pattern of AA at the hASIC1/3 TM1, we mutated the three sites in hASIC3 to their equivalent residues in hASIC1 and vice versa. We saw that simulating hASIC3 with the mutations Y58C, R65Q, and R68F in a 90:10 POPC:AA membrane yielded a TM1 interaction pattern similar to wild-type hASIC1 (Fig. 7 A). hASIC1 with the reverse mutations (C59Y, Q66R, F69R) yielded a similar pattern to wild-type hASIC3 (Fig. 7 B). Significantly, these three mutations changed the overall occupancy at all residues in the upper portion of TM1, and not only at these three positions (Fig. 7). We speculate therefore that these three sites could be responsible for the subtle differences in overall occupancy between hASIC1a and hASIC3 observed in Fig. 2 and give rise to different binding motifs on the two channels, despite the similarity in modulation.

Early studies showed that a larger number of double bonds in the fatty acid tail increased ASIC potentiation by fatty acids, whereas a smaller number of double bonds decreased potentiation (Smith et al., 2007). Similarly, it has been shown that at least three double bonds are required for potentiation by fatty acids with 20- or 22-carbon tails (Klipp and Bankston, 2022). It was therefore suggested that the number of double bonds determines the degree of potentiation (Klipp and Bankston, 2022). Furthermore, cis double bonds, specifically, increase the flexibility of fatty acid tails and can potentiate, e.g., voltage-gated ion channels (Elinder and Liin, 2017). Fatty acids with trans double bonds, on the other hand, are ineffective.

Thus, we probed whether changing the acyl tail length or saturation affected the occupancy of fatty acids at the ASIC transmembrane domain. We speculated that increasing PUFA tail unsaturation increased tail flexibility, and thus allowed for fatty acid positioning at the slightly concave transmembrane domain of ASICs. Indeed, DHA, which has six cis double bonds and a 22-carbon tail, showed an overall slightly higher occupancy along the upper TM1 residues with a conserved binding pattern, compared to AA. In contrast, palmitic acid, which has 16 carbon atoms in its tail and no double bonds, has significantly decreased occupancy for all residues (Fig. 8). Our data therefore support that for AA interactions with ASICs, a flexible tail with multiple cis double bonds is indeed required. This is consistent with recent functional data (Klipp and Bankston, 2022) showing that multiple cis double bonds as well as a negatively charged headgroup are important in interactions with the ASIC channel.

We noticed that R65 is located slightly below the polar headgroup region of the membrane bilayer, with AA therefore interacting with this residue below the overall headgroup region. We wondered whether binding of AA at the ASIC TMD could create a local deviation in the membrane curvature immediately adjacent to the lateral fenestrations located in the upper region of the TMD (Fig. 1), which have been proposed to conduct ions through the pore (Gründer and Chen, 2010). This could subsequently increase the ion flow through the fenestrations and present a possible mechanism for ASIC modulation by AA and alter pH dependence of the channel by influencing I_max overall and thereby shifting the pH₅₀ activation curve. To probe this, we analyzed the local membrane curvature surrounding the open state of hASIC3 in a 10% AA membrane as well as that in a pure POPC membrane. We did not observe a noticeable difference in local membrane curvature between the two membrane types (Fig. 9). We therefore proposed that modulation of ASICs by AA occurs through direct structural changes due to AA biding, which are likely to occur at a timescale beyond the current capabilities of our atomistic simulations.

In this molecular dynamics simulation study, we have described a potential AA binding pattern for ASICs. Both hASIC3 and hASIC1a exhibit accumulation of AA and high-occupancy AA interactions at the outer leaflet region of TM1, albeit with differences in the pattern of coordinating residues. AA interacts with hASIC3 primarily at three residues; Y58_hASIC3, R65_hASIC3, and R68_hASIC3, with Y58_hASIC3 and R65_hASIC3 conserved among the ASIC3 subtype. On the other hand, the residue Y71_hASIC1, which is conserved in the ASIC1 subtype, exhibits high occupancy interactions in hASIC1. If these interactions play a role in the modulatory effects of AA, these residues could therefore cause potential differences in the extent of lipid modulation of these two ASIC subtypes, affecting perhaps affinity or efficacy. Furthermore, by comparing AA occupancies to those of other fatty acids, we show that fatty acid tail length and number of double bonds may be important to their affinity to ASICs. This finding agrees with experimental data illustrating functional differences depending on the number of double bonds or the tail length (Smith et al., 2007; Elinder and Liin, 2017; Klipp and Bankston, 2022).

In both ASIC subtypes, our simulations suggest that the charged fatty acid headgroup participates in electrostatic/polar interactions with residues in the upper region of TM1, while the tail interacts with hydrophobic residues several turns lower on the same helix. We proposed that in hASIC3, the charged AA headgroup forms salt bridges with R65_hASIC3 and R68_hASIC3, while in hASIC1, AA could form an ion–dipole interaction with Y71_hASIC1. Our atomistic simulations support that AA participates in stable interactions with R65_hASIC3.

Since R65 and R68 are on the same face of TM1 in hASIC3, we speculate that these two residues could contribute to a region of positive electrostatic potential in this area. Compared to hASIC1, which does not have positively charged residues in this region, this could create a larger accumulation of AA in the upper TM1 region, as was seen in our simulations (Fig. 2). We proposed that this slightly higher degree of attraction, and therefore greater degree of AA accumulation, could affect the affinity or efficacy of hASIC3 modulation by AA.

Surprisingly, we did not see high AA occupancy at the arginine at position 63_hASIC3, which is conserved among both ASIC3 and ASIC1 types. We did, however, despite this low occupancy, see relatively longer-lived interactions at this residue in the open state of our hASIC3 model compared to the resting state. This could suggest that, although MD simulations are not able to capture extensive interactions between AA and R63_hASIC3, this residue could participate in AA interactions and constitute part of the modulatory mechanism. In a recent study, mutation of this arginine residue to glutamine in rat ASIC3 and ASIC1a eliminated potentiation by the fatty acid DHA (Klipp and Bankston, 2022). We proposed two reasons for this discrepancy between our CG molecular dynamics (CG-MD) simulations and the electrophysiology results. Firstly, R63_hASIC3, compared to R65_hASIC3 and R68_hASIC3, is oriented toward the TMD helix bundle rather than toward the membrane environment. Since we restrained the tertiary and quaternary structure of ASIC to the conformation of the solved crystal structure, the position of the TMD helices cannot readily change to accommodate lipid binding between tightly packed helices. However, without restraints, protein secondary and tertiary structure does not remain stable in CG simulations. Therefore, in our simulations, R63_hASIC3 is less solvent exposed and does not interact readily with AA, potentially due to the inherent properties of CG-MD. Furthermore, we did not observe interactions between AA and R63_hASIC3 in our atomistic simulations, in which the protein was unrestrained, although this lack of sampling could additionally be caused by the relatively shorter timescales of atomistic simulations. This may also suggest that this residue remains poorly solvent exposed even when restraints are removed, at least in our hands.

Moreover, the structure of the resting state of chicken ASIC1 shows a possible salt bridge between R65_cASIC1 and E413_cASIC1, which is broken in the open state structure (Fig. 10, A and B). The corresponding interaction between R63_hASIC3 and E430_hASIC3 remains stable in our atomistic simulations of the resting state of hASIC3 (Fig. 10 C) but only sometimes re-forms in simulations of the open state (Fig. 10 D). An R63Q_hASIC3 mutation will eliminate this salt bridge, and in doing so, potentially destabilize the resting state of the channel. This is supported by the fact that, in the absence of potentiating fatty acids, this mutation still shifted the pH₅₀ of both rat ASIC1a and ASIC3 to slightly more alkaline pH (Klipp and Bankston, 2022). In rat ASIC3, this shift in pH₅₀ caused by the R63Q_hASIC3 mutation was similar to that observed for WT ASIC3 when potentiated by DHA (Klipp and Bankston, 2022). The equivalent residue has also been previously shown to be involved in proton-mediated gating of hASIC1 (Chen et al., 2021). It is therefore unclear whether this mutation is directly involved in AA binding, or whether it is involved in structural changes relating to the activation and gating of the ASIC channel. To build on these experimental findings, as well as our simulation data, we suggest that R63_hASIC3 could be involved in ASIC modulation by AA, but the mechanism may be more complex and involve a conformational bias toward the open state of the channel.

We finally predict that the effect of AA on ASICs does not arise from changes in local membrane curvature, and thus further support that AA binding may have direct effects on the conformation of the ion channel. However, such conformational changes are likely to occur at timescales inaccessible to atomistic simulations. Further functional studies are needed to confirm the presence of the proposed regulatory sites, as well as to investigate whether these sites cause small differences in AA modulation between hASIC1a and hASIC3. A combination of mutagenesis and electrophysiological recordings in membranes containing AA could determine whether these identified residues play a functional role in ASIC activation, desensitization, ion conductance and PUFA modulation, after which more detailed simulation studies could shed further light on the molecular mechanism of fatty acid modulation of ASICs. In combination with functional studies, additional in silico mutants, such as the additional mutations (e.g., Y67 in both ASIC subtypes) studied in the work of Klipp and Bankston, could give a more detailed overview of AA interactions in this region. This could further be elaborated on by performing free energy calculations to study the identified ASIC-AA interactions. Given the complex nature of fatty acid effects on ASIC function, and more broadly, lipid effects on ion channel function, it would furthermore be valuable to simulate these channels in complex bilayers containing other lipid species found in neuronal membranes to better understand the interplay between the protein and its cell environment.

Christopher J. Lingle served as editor.

WestGrid and Compute Canada are acknowledged for providing the computing facilities. We thank David MacLean for helpful discussions.

We thank Natural Sciences and Engineering Research Council and the Canada Research Chairs Program for providing funding for this research (Natural Sciences and Engineering Research Council Discovery Grant RGPIN 2019-06864 and the Canada Research Chairs Program 950-232154 to M. Musgaard).

The authors declare no competing financial interests.

Author contributions: Both authors designed the study. A. Ananchenko performed the simulations and analyzed the data; both authors interpreted the results. M. Musgaard obtained the funding. A. Ananchenko wrote the manuscript with input from M. Musgaard. Both authors approved the final version of the manuscript.