Polyunsaturated fatty acids (PUFAs), but not saturated fatty acids, modulate ion channels such as the cardiac KCNQ1 channel, although the mechanism is not completely understood. Using both simulations and experiments, we find that PUFAs interact directly with the KCNQ1 channel via two different binding sites: one at the voltage sensor and one at the pore. These two amphiphilic binding pockets stabilize the negatively charged PUFA head group by electrostatic interactions with R218, R221, and K316, while the hydrophobic PUFA tail is selectively stabilized by cassettes of hydrophobic residues. The rigid saturated tail of stearic acid prevents close contacts with KCNQ1. By contrast, the mobile tail of PUFA linoleic acid can be accommodated in the crevice of the hydrophobic cassette, a defining feature of PUFA selectivity in KCNQ1. In addition, we identify Y268 as a critical PUFA anchor point underlying fatty acid selectivity. Combined, this study provides molecular models of direct interactions between PUFAs and KCNQ1 and identifies selectivity mechanisms. Long term, this understanding may open new avenues for drug development based on PUFA mechanisms.

## Introduction

Polyunsaturated fatty acids (PUFAs) affect many different ion channels, such as voltage-gated K+, Na+, and Ca2+ channels, as well as ryanodine receptors (Xiao et al., 2001; Hamilton et al., 2003; Oliver et al., 2004; Xiao et al., 2005; Ottosson et al., 2014; Farag et al., 2016; Tian et al., 2016; Liin et al., 2018; Bohannon et al., 2020). However, the molecular mechanisms behind the effects of PUFAs are not completely understood. The identity of PUFAs’ specific binding pockets remains uncharacterized for many ion channels; therefore, the molecular mechanisms of PUFA efficacy are still poorly understood. Atomistic understanding of PUFA binding sites and mechanisms of action on ion channels would further our ability to design PUFAs that modulate ion channel function and aid in the development of drugs based on PUFA mechanisms.

We previously showed that PUFAs act as activators of KCNQ1 (Liin et al., 2015), which is expressed, for instance, in cardiomyocytes (forming the KCNQ1/KCNE1 channel) and smooth muscle cells (Barrese et al., 2018). KCNQ1 forms a tetrameric, voltage-gated potassium channel with six transmembrane segments (S1–S6) in each subunit (Sun and MacKinnon, 2017). Recently, a cryo-EM structure of the Xenopus laevis KCNQ1 (xKCNQ1) channel was published (Sun and MacKinnon, 2017). Helices S1–S4 in each subunit form a voltage-sensing domain with positively charged arginine residues in S4 acting as the voltage sensor. S5–S6 from all four subunits together form the centrally located pore domain with the potassium selectivity filter and activation gate. Outward movement of the positively charged S4 in response to depolarization triggers opening of the activation gate and allows for the outward flux of potassium (Osteen et al., 2012; Barro-Soria et al., 2014; Zaydman et al., 2014; Hou et al., 2017; Westhoff et al., 2019).

Experiments combining electrophysiology and molecular biology provided evidence that PUFAs affect both the voltage sensor and the pore of KCNQ1 channels (Liin et al., 2016; Liin et al., 2018). PUFAs shift the voltage dependence of channel opening (detected as a shift in the voltage at which the conductance is half the maximal conductance [ΔV50] in Fig. 1 A). The outermost S4 arginine, R218 in xKCNQ1 (which corresponds to R228 in human KCNQ1; hKCNQ1), is important for this effect (Fig. 1 B; Liin et al., 2015; Liin et al., 2018). PUFAs also increase the maximum conductance of KCNQ1 channels (ΔGmax in Fig. 1 A). A conserved positively charged S6 residue, K316 in xKCNQ1 (which corresponds to K326 in hKCNQ1), is important for this increase (Fig. 1 B; Liin et al., 2018). The two effects of PUFAs on KCNQ1 channels are independent of each other (Liin et al., 2018). However, our previous studies have not allowed us to determine if there is one PUFA binding site on KCNQ1, from which one PUFA molecule exerts the two effects, or if there are two PUFA binding sites on KCNQ1: one close to the voltage sensor and one close to the pore. Moreover, molecular understanding of why PUFAs such as linoleic acid (LIN), but not monounsaturated or saturated fatty acids (FAs; such as the saturated stearic acid [STE]), activate KCNQ1 channels (Fig. 1, C and D; Liin et al., 2015) and how the PUFA tail interacts with KCNQ1 channels are still unclear. Here, we conducted computer simulations combined with electrophysiology experiments to address these open questions.

## Materials and methods

### Structure preparation

The 3-D structure used for the transmembrane domain of KCNQ1 channel (residues 94–348) was built using the recently published structure of xKCNQ1 channel solved by cryo-EM at 3.7 Å (Sun and MacKinnon, 2017). Throughout the Results section, xKCNQ1 numbering of residues is used if not stated otherwise. In this structure of KCNQ1, S4 is in an activated “up” position, whereas the pore is closed. Although PUFA interaction with the pore domain could potentially be different in closed and open channels, we anticipate any such differences to be minor for the following reasons. Previous work on KV channels show that the main structural differences between a closed and open pore is at the internal gate of the pore, which is in the intracellular half of the channel (Bavro et al., 2012; Jensen et al., 2012; Lenaeus et al., 2017; Sun and MacKinnon, 2017; Kuenze et al., 2019). In contrast, limited structural differences between closed and open pores are seen at the lipid-facing surface of the extracellular half of the pore (Bavro et al., 2012; Jensen et al., 2012; Lenaeus et al., 2017; Sun and MacKinnon, 2017; Kuenze et al., 2019), which is the part of the channel with which PUFA interacts. Moreover, previous experiments suggest that PUFA interaction with the pore domain increases the maximum conductance of KCNQ1 by affecting the selectivity filter (Liin et al., 2018), and the selectivity filter, which is located in the extracellular half of the channel, is not expected to change much between closed and open states (Cordero-Morales et al., 2006; Bavro et al., 2012; Jensen et al., 2012; Lenaeus et al., 2017).

CHARMM-GUI Martini Bilayer Maker (Jo et al., 2008) was used to prepare the system containing the protein embedded in a pure phosphatidylcholine bilayer solvated in 0.1 M KCl aqueous solution using CHARMM36 force field and TIP3P water model (Jorgensen et al., 1983; Noskov et al., 2004; Noskov and Roux, 2008). The system was initially subjected to steepest descent energy minimization and subsequently equilibrated for ∼90 ns before starting 1-µs production simulation using Gromacs with 2-fs time steps (Pronk et al., 2013). The LINCS algorithm (Hess et al., 1997) was applied for constraining bond lengths. Electrostatic interactions were calculated with the Particle-Mesh Ewald algorithm at every step (Essmann et al., 1995). A 1.2-nm cutoff was used both for electrostatics and van der Waals interactions, with neighbors list updated every 20 steps. The simulations were performed at constant pressure of 1.0 bar with Parrinello–Rahman pressure coupling (Parrinello and Rahman, 1981) and the semi-isotropic pressure scaling, time constant of 5.0 ps, and a system compressibility of 4.5e-5 bar-1. The temperature of the system was maintained at 300°K using the extended Nosé–Hoover thermostat (Nosé, 1984). The 1-µs fully equilibrated structure was subsequently used to seed coarse-grained (CG) simulations in multicomponent lipid membranes containing either the ω-6 PUFA LIN or the saturated FA STE. STE has an 18-carbon-long acyl tail, just like LIN. The key difference, however, is that the STE tail is fully saturated (compare structures in Fig. 2). In the CG model, STE differs from LIN in the flexibility between beads 2 and 3, counting from the tail end (Fig. 2). In contrast to LIN, STE, like other saturated FAs (Liin et al., 2015), neither shifts V50 nor increases Gmax of KCNQ1 (Fig. 1, C and D). Therefore, STE was used in the simulations as a negative control to LIN.

### CG MD simulations

The relaxed structure of the KCNQ1 channel was used as the starting structure for the CG simulations with and without FAs in the model bilayer. The channel was embedded into a lipid bilayer consisting of phosphatidylethanolamine:phosphatidylglycerol:cholesterol:FA (POPE:POPG:CHOL:FA) in a 3:1:1:1 ratio. KCNQ1 requires the presence of phosphatidylinositol 4,5-bisphosphate (PIP2) lipids for its function; therefore, PIP2 lipids were included in the inner leaflet in every bilayer simulation. The lipid composition for these systems was as follows: 307 POPE, 102 POPG, 88 CHOL, 16 PIP2, and 117 FA, of which 63 were distributed in the outer leaflet and 54 in the inner leaflet. Adding FA to the simulation systems did not alter the overall phospholipid distribution around the channel compared with the system free of FAs (FA-FREE; Fig. S1 B). Each system was simulated using standard Martini protocols with minor variations between systems to accommodate system-specific issues. The CG systems are circa 14.3 × 14.3 nm in the membrane plane (x and y), including 614 lipids and circa 10,300 water molecules. Simulations were performed using Gromacs (Pronk et al., 2013), with the Martini v2.2 force field parameters (de Jong et al., 2013) and standard simulation settings (de Jong et al., 2016). Briefly, all simulations were performed with a 10-fs time step, a temperature of 300°K set using a velocity-rescaling thermostat, and a time constant for coupling of 1 ps. A semi-isotropic pressure of 1 bar was maintained with the Parrinello–Rahman barostat, with a compressibility of 3·10−4 bar−1 and a relaxation time constant of 12 ps. Production runs of 50 µs were performed in the presence of position restraints applied to the backbone beads, with a force constant of 50 kJ mol−1 nm−2. The CG trajectories were clustered, and representative structures of the five-component lipid bilayer and an embedded channel were selected and back-mapped to the atomistic level using the CHARMM-GUI Martini to All-atom Converter (Jo et al., 2008). The clustering was limited to the FA head groups and performed with the built-in Gromacs tool gmx cluster, using the Jarvis–Patrick method (Jarvis and Patrick, 1973) with the distance cutoff set to 65 Å and the number of neighbors to three. The Jarvis–Patrick method adds a structure to a cluster when this structure and a structure in the cluster have each other as neighbors and they have a least a defined number of neighbors in common. The representative pose of every cluster, which is the structure with the smallest average RMSD from all other structures of the cluster, was inspected, and the frames in which FA head groups were within 10 Å of the identified high-contact and stable-contact residues at interaction sites I and II were chosen as representative structures. We selected four starting structures to seed our all-atom (AA) MD simulations on the ANTON2 platform: three structures were representative structures of the channel with LIN or STE bound to the different binding pockets mapped from the analysis of CG simulations and one was a frame with LIN fully dissociated from the channel. These structures were initially equilibrated for 500 ns with Gromacs, according to the protocol described above, with the exception of using the Berendsen barostat (Berendsen and Postma, 1984) for temperature and pressure control before starting AA runs on the ANTON2.

In addition to the above systems, two other systems were set up as described above and simulated for 50 µs: a five-component, membrane-only bilayer containing LIN and a five-component, membrane-only bilayer containing STE. The lipid composition for these systems was as follows: 572 POPE, 190 POPG, 164 CHOL, 218 FA, and 32 PIP2. The equilibrated 50-µs structure from each of these systems was selected as starting seeds for AA MD simulations, which were run for 1 µs with Gromacs.

### AA MD simulations

AA production runs were performed for the following KCNQ1 systems: LIN dissociated, LIN site I, LIN site II, and STE site I. The ANTON2 software version 1.27.0 from D. E. Shaw Research was used for production runs on the purpose-built ANTON2 supercomputer (Shaw, D. E., et al. 2014. Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. https://doi.org/10.1109/SC.2014.9). The production runs were performed for 5 to 15 µs each using the CHARMM36M force field to assess structural dynamics of KCNQ1 in various membranes. The latest CHARMM36 lipid parameters were used to describe multicomponent membranes (Klauda et al., 2012). CHARMM-NBFIX LJ parameters for K+ and Cl (Noskov and Roux, 2008; Luo and Roux, 2010) were used to simulate counter-ion dynamics with standard TIP3P water model (Jorgensen et al., 1983). The production runs were executed in a semi-isotropic (NPaT) ensemble at temperature 315°K maintained by the Nosé–Hoover thermostat (Martyna et al., 1994). Nonbonded and long-range electrostatic interactions were evaluated every 2 and 6 fs, respectively, using the RESPA multiple–time-step algorithm (Tuckerman et al., 1992). Long-range electrostatics was calculated using the k-Gaussian Ewald method implemented to enhance performance on ANTON2 platform (Shaw, D. E., et al. 2014. Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis) with a 64 × 64 × 64 grid. SHAKE was used to constrain all bonds involving hydrogen atoms. The multi-integrator (multigrator) algorithm (Lippert et al., 2013) developed in-house by D. E. Shaw Research was used for temperature and semi-isotropic pressure coupling (Shaw, D. E., et al. 2014. Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis). The time step for production runs was set to 2 fs, and trajectories were saved every 500 ps. To assess the structural integrity of the KCNQ1 channel, the RMSD of the protein backbone was calculated against the initial structure for all four systems. Overall, the channels maintain stable structures, settling around an RMSD of 2.0–4.0 Å throughout the simulation time (Fig. S2).

### AA MD simulations of the R218Q/R221Q (“pseudo S4 down”) system

To construct a pseudo S4 down state of S4, we neutralized in silico the two top charges of S4 by introducing the R218Q/R221Q mutations in the model. Representative poses of LIN and STE constant binders at site I were first superimposed across each channel subunit, after which Pymol (Schrodinger, 2010) was used to introduce the R218Q/R221Q double mutation in every channel subunit. CHARMM-GUI Bilayer Maker (Jo et al., 2008) was subsequently used to prepare the system containing the protein embedded in a multicomponent lipid bilayer, according to the protocol described above. The systems, referred to as R218Q/R221Q systems, were initially equilibrated for 20 ns with position restraint acting on the bound LIN and STE molecules to allow them to relax into their placed positions, using a force constant of 1,000 kJ mol−1 nm−2. Production runs of 500 ns were performed for three replicated systems for each FA using Gromacs (Pronk et al., 2013). As control, we set up identical systems for the WT channel, referred to as “S4 up” systems, of which three replicated runs were simulated for 500 ns each. We pushed the simulations to 5 µs for two of the R218Q/R221Q systems, one for the LIN system and one for STE, in order to look at the occupancy of these FAs across different layers around binding site I. Occupancies across each 0.5-Å layer along the bilayer normal were determined by counting the total number of contacts between each FA head group and any atom on residues R218, R221, F222, Y268, F269, and L272 constituting site I.

### Analysis of FA enrichment regions and interaction sites in CG simulations

The cutoff distance between two carbon atoms for a hydrophobic interaction is 4.0 Å, and the carbon-hydrogen bond length is at ∼1.09 Å. To capture the long-range hydrophobic interactions between LIN or STE molecules and the channel, the distance cutoff between any two atoms in our systems was set to 6 Å in CG simulations, which would in effect capture any long-range molecular interactions. The cutoff distance of 6 Å is used widely in contact analysis, as it captures all important interactions between two particles and is dependent on the distance of the first minimum of the Lennard–Jones potential of the Martini force field (Marrink et al., 2007; de Jong et al., 2012; Barbera et al., 2018). The number of LIN or STE molecules within 6 Å of the channel was calculated by counting the number of negatively charged head group beads within the cutoff distance to the protein surface during the entire 50 µs of the simulation, as previous experiments have demonstrated that it is the negatively charged PUFA head group that exerts an electrostatic effect on the KCNQ1 channel (Liin et al., 2015; Larsson et al., 2018; Liin et al., 2018). Enrichment regions for LIN, STE, or other membrane components (POPE, POPG, and CHOL) were calculated by their occupancies over the last 15 µs of each CG simulation using the VolMap plugin in VMD 1.9.1 with a grid resolution of 1 Å (Humphrey et al., 1996). The 2-D occupancy maps were based on the head group beads, and the contour maps were based on entire LIN and STE molecules.

To identify the FA interaction sites in KCNQ1, we performed lipid–protein interaction analysis. An interaction is defined as when an entire (i.e., both head group and tail) lipid molecule (LIN or STE) is located within 6 Å of the channel. PUFA interaction sites were identified by two separate measures: the total interaction time and the longest lifetime. The total interaction time calculates the total time of interaction observed between each residue and any FA (could be more than one FA at a time) over the last 45 µs of the simulation, regardless of the duration of each interaction. The longest lifetime, on the other hand, is a measure of the longest-lasting contact between each residue and any FA without any interruptions in the interaction. The total interaction time and longest lifetime were averaged across the four channel subunits to account for the fourfold symmetry of the channel. Those observations that fall above Q3 (third quartile) + 1.5 interquartile range in the total interaction time and longest lifetime measurements are referred to as high-contact and stable-contact residues, respectively. Distance measurements between the FAs and the channel were calculated with the built-in Gromacs tool gmx mindist, and the total interaction time and longest lifetime measurements were done by in-house R scripts. The radial distribution functions were calculated between the LIN and STE head group beads and the channel using the built-in Gromacs tool gmx RDF.

### Analysis of FA enrichment regions and interaction sites in AA simulations

Distances between FAs and the channel for the AA simulations were calculated with the MDTraj python library. In-house R scripts were used to calculate the total interaction time made by each FA molecule over the course of the simulations. The FA molecules that had total interaction times that were in the upper quartile of total interaction time were selected for the lipid-protein analysis and are referred to as the “UQ binders.” The measurements of total interaction time and longest lifetime of each residue in the AA simulations were carried out in a manner like the CG simulation data as stated above. The FAs that stayed bound to the channel during the entire 5-µs simulation were identified by calculating the longest lifetime for each FA and are referred to as the “constant binders.” Clustering was performed with the built-in Gromacs tool gmx cluster using the GROMOS method (Daura et al., 1999), with the distance cutoff set to 5 Å. By setting a distance cutoff, the GROMOS method takes the structure with the largest number of neighbors as the centroid of the first cluster and eliminates it from the pool of structures together with all its neighbors. This procedure is repeated for all remaining structures until all structures have been designated to a cluster. The center of each cluster, which is the structure with the smallest average RMSD from all other structures of the cluster, was taken as the representative pose. Hydrogen bonds were considered to form when a hydrogen-bond donor, that is, a hydrogen atom bound covalently to an electronegative atom, was within a cutoff distance of 3.5 Å of an electronegative hydrogen-bond acceptor. As a definition of salt bridge formation, we considered a cutoff of 4 Å between N-O atom pairs.

### Xenopus oocyte experiments

hKCNQ1 (also called KV7.1; GenBank accession no. NM_000218) was in expression plasmid pXOOM. Mutations were introduced using site-directed mutagenesis (QuikChange II XL with 10 XL Gold cells; Agilent), and constructs were sequenced at the core facility at Linköping University to ensure correct sequence. Complementary RNA was prepared using T7 mMessage mMachine transcription kit (Ambion/Invitrogen). The RNA concentration was quantified using spectrophotometry (NanoDrop 2000c; Thermo Fisher Scientific). Xenopus oocytes were surgically isolated at Linköping University. Animal experiments were approved by the local ethics committee. Isolated Xenopus oocytes were injected with 50 nl of KCNQ1 RNA (50 ng RNA), and oocytes were incubated at 16°C for 2–5 d before performing two-electrode voltage clamp experiments. The two-electrode voltage clamp recordings were performed using a Dagan CA-1B Amplifier. Currents were filtered at 500 Hz and sampled at 5 kHz. The holding voltage was generally set to −80 mV. Activation curves were generally generated in steps between −80 and +60 mV in increments of 10 mV (3-s duration). The tail voltage was set to −30 mV. The control solution contained 88 mM NaCl, 1 mM KCl, 15 mM HEPES, 0.4 mM CaCl2, and 0.8 mM MgCl2. pH was set to 7.4 using NaOH. All compounds were bought from Sigma-Aldrich if not stated otherwise. Stock solutions of 100 mM LIN (9-cis, 12-cis LIN) or docosahexaenoic acid (DHA; cis-4,7,10,13,16,19–DHA), 50 mM STE, or 25 mM N-arachidonoyl taurine (N-AT; Cayman Chemical) were prepared in 99.5% ethanol. The final test solution was prepared shortly before experiments. Control or LIN or STE solution was applied using a Minipuls 3 peristaltic pump (Gilson). LIN or STE solution was applied until the effect on current amplitude reached steady state or for a minimum of 5 min. The chamber was cleaned in between each oocyte using ethanol-supplemented control solution.

### Electrophysiological analysis

Electrophysiological analysis was performed in GraphPad Prism 8 (GraphPad Software Inc.). To quantify the voltage dependence for channel opening, tail currents were measured shortly after stepping to the tail voltage and plotted against the preceding activation voltage. A Boltzmann function was fitted to the data to generate the conductance versus voltage (G(V)) curve:
$G(V)=Gmin+(Gmax−Gmin)/{1+exp[(V50−V)s]},$
(1)
where Gmin is the minimal conductance, Gmax the maximal conductance, V50 the midpoint (i.e., the voltage at which the conductance is half the maximal conductance determined from the fit), and s the slope of the curve. The difference in V50 induced by LIN in each oocyte (i.e., ΔV50) was calculated to quantify the shift in the voltage dependence for channel opening. The difference in Gmax induced by LIN in each oocyte (i.e., ΔGmax) was calculated to quantify the change in the maximal conductance. Biophysical properties of tested constructs are summarized in Table S1.
To plot the concentration dependence of the LIN-induced effect on V50 or Gmax as a function of the LIN concentration, the following concentration-response curve was fitted to the data:
$ΔEffect=ΔEffectmax/[1+([LIN]50[LIN])]H,$
(2)
where ΔEffectmax is the maximal shift in V50 or change in Gmax, [LIN]50 the LIN concentration needed to cause 50% of the maximal effect, and H the Hill coefficient. Limited solubility of LIN at concentrations above 70 µM prevented us from quantifying the LIN effect at higher concentrations. Therefore, the Hill coefficient of the concentration-response curves was constrained to either 1 or −1 to make the fits more robust. Average values are expressed as mean ± SEM. Statistics was calculated using either one-way ANOVA followed by Dunnett’s multiple comparisons test (to compare with WT) or one-sample t test (to compare with a hypothetical value). P < 0.05 was considered statistically significant.

### Online supplemental material

Fig. S1 shows the lipid distribution around KCNQ1 in different systems. Fig. S2 shows channel stability in AA simulations. Fig. S3 presents contact analysis for inner leaflet LIN and STE in CG simulations. Fig. S4 presents contact analysis for AA LIN simulations. Fig. S5 shows that the pattern of LIN interaction with KCNQ1 in AA simulations is preserved in longer and initial seed unbiased simulations. Fig. S6 presents contact analysis for LIN constant binders in AA simulations.Fig. S7 shows binder residues at sites Iα, Iβ, and II. Fig. S8 illustrates the representative pose of a LIN constant binder at sites Iα and Iβ. Fig. S9 shows that the mutation of nearby residues at LIN interaction sites displays lack of effect. Fig. S10 presents numerical values corresponding to Fig. 5, C and D. Fig. S11 presents numerical values for differences in LIN and STE occupancy in the R218Q/R221Q system. Fig. S12 shows the experimental effect of LIN on channel opening and closing kinetics. Fig. S13 shows the impact of Y268 for experimental DHA and N-AT effects. Video 1 shows the top view of the 5-µs AA simulation of LIN in the site I system. Video 2 shows the top view of the 5-µs AA simulation of LIN in the site II system. Table S1 lists the biophysical properties of KCNQ1 mutants.

## Results

We studied the interaction of the ω6 PUFA LIN or saturated FA STE with the KCNQ1 channel by combining long-time-scale CG MD and microseconds-long AA MD simulations in a multicomponent membrane bilayer with a lipid composition of POPE, POPG, CHOL, and PIP2 (inner leaflet only; Fig. S2, A and B) and PUFAs added to both leaflets. The levels of the different lipids were chosen to be comparable to the concentrations in simplified models of the cardiomyocyte membrane (Marrink et al., 2019). 50-µs-long CG simulations allowed us to directly map putative PUFA–channel interaction sites by evaluating total interaction time and the longest lifetime of each residue. To further characterize these sites, we proceeded with 5–15-µs-long AA MD simulations and identified the so called constant binders, which are molecules that stay bound to the putative interaction sites during the entire simulation, and binder residues, which are residues that interact in a direct and continuous manner with the constant binder. Based on our previous experimental findings that residues in the extracellular part of the channel are important for PUFA activation of the KCNQ1 channel (Liin et al., 2015; Larsson et al., 2018; Liin et al., 2018), we focused our analysis on PUFA–channel interactions in the outer leaflet of the lipid bilayer. As a control, we further tested putative PUFA binding sites in the inner leaflet (see LIN–KCNQ1 interactions).

### LIN, but not STE, is highly enriched in two distinct regions of the KCNQ1 channel

50-µs CG MD simulations of KCNQ1 in a multicomponent lipid bilayer, with or without LIN, were performed to map the lipid distribution around the channel. The number of LIN molecules within 6 Å of the KCNQ1 channel remained nearly constant throughout the 50-µs-long simulation, with an average of 6.6 LIN molecules within 6 Å around the entire channel (Fig. 2 A). 2-D occupancy maps of the LIN head group in different layers along the z axis of the membrane (bilayer normal) revealed two distinct high-occupancy regions in the different layers across all four voltage-sensing domains (Fig. 2 B). A high-occupancy layer marked by R218 in the KCNQ1 model contributes to binding region I formed by residues from S3, S4, and S5 helices (Fig. 2 B, left). A layer marked by K316 in the KCNQ1 model defines binding region II, composed of residues from S1 and S6 helices (Fig. 2 B, right). These highly LIN-enriched binding regions are clearly visible in the contour map of occupancy of entire LIN molecules, highlighting LIN enrichment at both regions I and II across at least three voltage-sensing domains (Fig. 2 C). The two sites in proximity of R218 and K316, respectively, are consistent with previous experimental findings of two independent PUFA effects on KCNQ1 reliant on charge conservation in these positions (Liin et al., 2018). CG MD simulations identified important quantitative differences in STE enrichment at the KCNQ1 channel surface: the number of STE compared with LIN molecules near the channel differed by a factor of two to three (Fig. 2 and Fig. S1 C), with no enrichment of STE at region II and limited enrichment at region I (Fig. 2, E and F).

### CG MD simulations identify two main interaction sites for LIN on KCNQ1

To further characterize the identified regions of PUFA–channel interaction, we measured the total interaction time of each residue with every LIN molecule (interaction = a LIN within 6 Å of the residue) on the external leaflet of the membrane bilayer (Fig. 3 A), an approach that allowed us to identify two distinct sites defining LIN–channel interaction (Fig. 3, B and C). High-contact residues located on the upper parts of transmembrane segments S4 and S5 constitute interaction site I, which forms a pocket between S4 and S5. The high-contact residues at site I include R218 and F222 in S4 and F269 and L272 in S5. High-contact residues located in the upper parts of transmembrane segments S6 and S1 constitute interaction site II, which is formed by a cavity in between these two helices. The high-contact residues at site II are I314 in S6 and F129 and Y138 in S1. Another measure of PUFA–channel interaction is the residues–LIN interaction with the longest lifetime for each residue (Fig. 3 D). The longest lifetime measurements distinguished long-lived, stable contacts from short-lived, but frequent, low-affinity contacts and highlight the regions of higher affinity binding on the channel surface (Fig. 3, E and F). The stable-contact residues in site I include R218, G219, R221, and F222 in S4 and Y268, F269, and L272 in S5. At interaction site II, the stable-contact residues are W313, I314, and T317 in S6 and I128 in S1.

Thus, the combined approach of measuring total interaction time with the longest lifetime of interacting residues of KCNQ1 enabled us to identify two plausible PUFA–channel interaction sites within the regions of PUFA enrichment. These two regions are in excellent spatial agreement with regions I and II identified by the occupancy maps. In contrast, for STE the total interaction time as well as the longest lifetime were reduced by at least a factor of 2 compared with that of LIN, also with differences in which residues constitute high-contact residues (Fig. 3, G–L).

### LIN–KCNQ1 interactions in inner leaflet do not contribute to functional effects

Although our analysis was focused on interactions in the outer leaflet, the CG MD simulations were performed with either LIN or STE also in the inner leaflet. In analogy with interactions in the outer leaflet, the total interaction time as well as the longest lifetime were longer for LIN compared with STE in the inner leaflet (Fig. S3, A–D). Clustering of LIN molecules interacting with the inner part of KCNQ1 identified R239 and R249 as putative interaction residues with the LIN head (Fig. S3, E and F). However, two-electrode voltage clamp experiments on hR249Q (xR239) and hR259Q (xR249) mutants did not support any functional roles of potential LIN interactions with these residues for LIN effects: both mutants responded to LIN with induced V50 and Gmax effects comparable to those of WT (Fig. S3 G and figure legend). These data suggest that although LIN shows more prominent accumulation than STE near KCNQ1 in the inner leaflet, LIN interactions with KCNQ1 residues in the inner leaflet do not contribute significantly to the functional effects on KCNQ1. This finding highlights the importance of functional testing of putative binding sites identified in simulations.

### AA simulations resolve atomistic details for the interaction of the LIN head and tail with KCNQ1 at site I

We next examined LIN and STE interactions with KCNQ1 at sites I and II on the atomistic level. We first focused the analysis on LIN interaction with KCNQ1 in two AA simulations of systems, where CG MD frames in which LIN molecules that were in the vicinity of each of the two identified interaction sites were converted to atomistic structures and were used to seed 5-µs AA MD simulations (referred to as LIN site I and site II systems in Materials and methods). Long-lived interactions between LIN molecules and KCNQ1 over the course of the AA MD simulations were monitored by determining the longest interaction time and longest lifetime of LIN molecules in the upper quartile. This analysis of upper-quartile LIN molecules reinforces the previously identified sites I and II. Fig. 4 A displays typical density of LIN at site I, with the LIN head group concentrated at the top of S4, whereas the density associated with the mobile tail is distributed with a higher intensity near S5 (see also Fig. S4 A and Video 1). We refer to this location of LIN molecules on the S5 side of S4 as site Iα. At site I, we also observed LIN molecules on the S3 side of S4, referred to as site Iβ, with the density of the LIN head group at the top of S4 and tail density scattered between S3 and S4 (Fig. S4 C). The LIN density at site II is concentrated in the cavity formed by S1, S5, and S6 in a similar manner as in the CG simulations for both the head and tail groups of the molecule (Fig. 4 B, Fig. S4 B, and Video 2).

To ensure sufficient sampling in the AA MD simulations, we extended the AA simulations for the LIN site I system for an additional 10 µs, which showed a similar pattern of LIN interaction with KCNQ1 (Fig. S5 A). To rule out any potential initial seed system bias in AA MD simulations, one CG frame with no LIN molecules in close proximity of site I and II was selected as a seed system for a 5-µs-long AA MD simulation. In this system, started from a LIN-dissociated state, LIN diffuses and stably binds to each of the identified sites I and II (Fig. S5, B and C).

We next turned our attention to three LIN constant binders, one at each distinct interaction site (sites Iα, Iβ, and II), to identify channel residues important for LIN interaction. Through distance measurements of each residue identified as interacting with the constant binders, we were able to discern between proximal residues (i.e., those residues that are merely within 6 Å of the binder but do not interact directly with LIN) and binder residues (i.e., residues that are in direct contact with the constant binder LIN molecule and form important interactions of either an electrostatic or hydrophobic nature; Fig. S6 and Fig. S7). The promiscuous LIN tail required clustering of all the LIN binder conformations to allow for selection of representative poses of the LIN constant binder. The representative pose of the largest cluster of LIN binder conformations at site Iα (encompassing 70% of all LIN conformations) shows a LIN constant binder that sits snugly in between S4 and S5: its negatively charged head group forms direct interaction with R218, while its curled-up tail is lined up along residues situated on S4 and S5 deep into the membrane (Fig. 4, C and D). The binder residues encompass the basic residues R218 and R221 on S4 and the aromatic residue Y268 on S5, which through a network of salt bridges and hydrogen bonds stabilize the LIN head group in between S4 and S5. Meanwhile, the mobility of the tail allows for hydrophobic interactions with residues F222, L226, and L229 on S4 and the S5 residues F265, F269, and L272 (Fig. 4, D and E; and Fig. S7 A). The representative poses from additional clusters at site Iα and Iβ display an overall comparable pattern of interaction between the LIN and KCNQ1 (see Fig. S6, Fig. S7, and Fig. S8 for further details). Based on the similarity of interactions between the LIN head group and the top arginines of S4 in sites Iα and Iβ, we consider these sites as variants of site I with anticipated similar effects on KCNQ1 voltage dependence.

### AA simulations resolve atomistic details for the interaction of the LIN head and tail with KCNQ1 at site II

At site II, the properties of the binder residues lining this pocket are, as for site I, of either charged, aromatic, or hydrophobic nature (Fig. S7 C). At site II, the LIN constant binder conformations form one large cluster encompassing 98% of all conformations. The representative pose displays a LIN constant binder that sits comfortably in the cavity formed by S5, the pore helix, and S6 (Fig. 4 F). While the head group interacts mainly with K316 on S6 via a salt bridge, the tail forms a number of hydrophobic interactions with residues F260 and L263 on S5; A290, L293, W294, and V297 on the pore helix; and S320 and V324 on S6 (Fig. 4, F and G; and Fig. S7 C). Here, we see how the LIN constant binder maintains largely the same conformation and binds quite stably in this narrow interaction site (Fig. 4 H).

Altogether, the detailed probing of LIN–KCNQ1 interactions in the AA simulations has enabled a finer characterization of the LIN interaction sites. A common feature found for all sites is the bi-component mechanism of specific interactions: the stable binding relies on LIN head group electrostatic interactions with a positively charged residue located at the top of the channel close to water–membrane interface, while the LIN tail interacts with a network of hydrophobic residues embedded deeper into the membrane. In contrast, we distinctly see in AA simulations how the STE fails to interact with the channel at site II (Fig. S6 B), and is unable to maintain close interactions with the S5 residues Y268 and L272 at site I (Fig. S7). Overall, the stiffness of the STE tail prevents this saturated FA from interacting as efficiently as the more mobile unsaturated tail of LIN with the network of hydrophobic residues at site I or from occupying the more restricted cavity at site II (Fig. S6 and Fig. S7). The involvement of aromatic residues in the stabilization of the unsaturated LIN tail is in line with the recent observation of a specific interaction between unsaturated tails and aromatic residues of a sensor protein (Ballweg et al., 2020). Using an array of biophysical and modeling techniques, Ballweg et al. (2020) showed that the presence of unsaturated alkyl tails in lipids stabilizes lipid–aromatic group interactions via, presumably, specific CH-π stacking interactions. At the same time, the presence of saturated lipid tails is a destabilizing factor for protein–lipid interactions.

### Electrophysiology experiments verify LIN tail–KCNQ1 interactions predicted from simulations

Previous electrophysiology experiments have shown that the basic residues on S4 and S6 identified as binder residues for the PUFA head in the AA simulations (i.e., R218 and R221 on S4 and K316 on S6) are important for PUFA effects on V50 and Gmax, respectively (Liin et al., 2018). However, the importance of the predicted binder residues for the PUFA tail from the AA simulations remains experimentally unstudied. We therefore used two-electrode voltage clamp electrophysiology on KCNQ1 to experimentally test the effect of mutating specific binder residues.

For site I, we focused on F222, F269, and L272, which interacted most consistently with the LIN tail (Fig. S7 A). We reasoned that reducing the size of the side chain of these hydrophobic residues through substitution with the hydrophobic but smaller alanine may impair the ability of the LIN tail to interact with these residues to enable the shift in V50. Whereas alanine mutation of hF232 (xF222) did not affect the ability of LIN to shift V50, alanine mutation of hF279 (xF269) or hL282 (xL272) impaired the ability of LIN to shift V50, which was detected as an apparent lower LIN affinity for the V50 effect seen as a shift in the concentration-response curve toward higher LIN concentrations (Fig. 4 I). At 70 µM of LIN, the hF279A mutant shifted V50 by only −7.7 mV compared with −14.3 mV WT (P < 0.05 with one-way ANOVA followed by Dunnett’s multiple comparisons test). In contrast, neither of the site I mutants for which Gmax could be determined significantly altered the ability of LIN to increase Gmax (Fig. S9 A).

For site II, we focused on L263, A290, L293, V297, and S320, which interacted consistently with the LIN tail (Fig. S7 C). Because of the narrow nature of site II, we reasoned that increasing the size of the side chain of these residues through substitution with hydrophobic but bulkier phenylalanine or tryptophan may impair the ability of the LIN tail to interact with this site, a property of PUFA essential for an increase in Gmax. In general, mutation of this site was poorly tolerated by the channel, but the effect of mutation of hL273 (xL263) and hA300 (xA290) could be experimentally tested. Tryptophan mutation of hL273 (xL263) and hA300 (xA290) impaired the ability of LIN to increase Gmax, which was detected as a smaller maximal effect induced by LIN (Fig. 4 J). At 70 µM of LIN, the hL273W mutant increased Gmax by <20% compared with 49% WT (P < 0.05 with one-way ANOVA followed by Dunnett’s multiple comparisons test), and hA300W did not significantly increase Gmax (P > 0.05 with one-sample t test compared with a hypothetical value of 0). In contrast, neither of the site II mutants significantly altered the ability of LIN to shift V50 (Fig. S9 B).

As a negative control, we introduced similar mutations in the vicinity of the proposed binding residues to test the potential impact of mutating nearby residues predicted to not directly engage with LIN. The hL271A (xL261, near site I) and hI274W (xI264, near site II) mutations did not impair the ability of LIN to shift V50 (Fig. S9, C–E) or increase Gmax compared with WT (Fig. S9 E, legend). These results support the conclusion that mutation of the predicted hydrophobic binder residues specifically alter LIN–channel interactions. Altogether, these experimental results are in concordance with the prediction from the CG and AA simulations that several hydrophobic residues constitute important binder residues for the PUFA tail in both sites I and II.

### Ability of FA head group to hydrogen bond with S5 tyrosine Y268 is important for LIN/STE selectivity in site I

The predicted interactions between Y268 and the LIN head group through hydrogen bonding (illustrated in Fig. 5 A) suggests that this tyrosine is a novel potential PUFA binding partner. The limited ability of STE to retain close distance to Y268 indicates that the tyrosine may play a role in FA selectivity (Fig. S7 A). We experimentally mutated Y268 to either a phenylalanine (to remove just the OH group) or an alanine to probe the functional impact of this residue. Phenylalanine or alanine mutation of hY278 (xY268) clearly impaired the ability of LIN to shift V50 (Fig. 5 B), emphasizing the importance of hydrogen-bonding interactions of the amphiphilic side chain of Y268. At 70 µM of LIN, the hY278F and hY278A mutants shifted V50 by only −5.0 and −3.8 mV compared with −14.3 mV for WT (P < 0.05 with one-way ANOVA followed by Dunnett’s multiple comparisons test to compare with WT). The small remaining effect of LIN on hY278F and hY278A mutants (−5 mV or smaller) is more like the lack of effect of STE at this concentration (0.9 mV) and can be attributed to the residual compensation of the hydrogen-bonding interactions by aromatic or hydrophobic stabilization (Lees-Miller et al., 2009). Thus, experiments support an important role of the S5 tyrosine for LIN effects at site I and suggest that the different ability of LIN and STE to hydrogen bond with Y268 might underlie selectivity at site I.

### Large difference in ability of LIN and STE to interact when outermost S4 arginines are neutralized

Our experimental data show that LIN induces a negative shift in the voltage dependence of activation, suggesting that LIN binds stronger in the activated S4 up state of the voltage sensor than in a resting “S4 down” state of the voltage sensor. In addition, STE does not induce a significant shift in the voltage dependence of activation. To better understand the difference in effect of LIN and STE, we estimated LIN and STE binding to different S4 states. An experimentally determined structure of the S4 down state of KCNQ1 is not yet available. Therefore, we constructed a pseudo S4 down state by neutralizing in silico the two top charges of S4 by introducing the R218Q/R221Q mutations in the model. If S4 moves as in the sliding helix model of S4 (Broomand and Elinder, 2008), where each S4 charge takes the place of its next S4 charge, then, from an electrostatic point of view, removing the two top S4 charges is similar to moving S4 down two “clicks” as a sliding helix (Fig. S10). While this is not a real S4 down state, it provides a simple model to test the importance of R218, R221, and Y268 for PUFA interaction with KCNQ1. We positioned either LIN or STE in each of the channel’s site I (i.e., four sites per channel) according to representative binding poses identified in the AA simulations and performed similar simulations for the R218Q/R221Q system as for the WT (S4 up) system. Three replicates of each system were simulated for 500 ns each (providing simulation data for a total of 12 LIN or STE molecules, respectively, in site I). In the WT system, there is a 2–3-fold difference in ability of LIN and STE to hydrogen bond with Y268 (Fig. 5, C and D; and Fig. S10). Neutralizing R218 and R221 reduces LIN’s ability to hydrogen bond with Y268, which is observed as reduced total interaction time and longest lifetime of roughly a factor of 6–12 in the R218Q/R221Q system compared with the WT system (Fig. 5, C and D; and Fig. S10). This likely reflects decreased stabilization of LIN in site I upon removal of the electrostatic coordination of LIN head group by S4 arginines. In the R218Q/R221Q system, there is a dramatic difference in the ability of LIN compared with STE to hydrogen bond with Y268: there are 13- and 37-fold differences in total interaction time and longest lifetime, respectively, for LIN compared with STE (Fig. 5, C and D; and Fig. S10). This suggests that once the electrostatic coordination of FA head group by S4 arginines is lost (as in the predicted down state of S4), STE largely fails to hydrogen bond with Y268.

We extended the simulation of one LIN and one STE R218Q/R221Q system to 5 µs to analyze occupancies of LIN and STE head groups at site I across the membrane normal with 0.5-Å slabs averaging for each consecutive 1-µs interval. These simulations revealed a 20–300-fold difference in occupancy of LIN compared with STE within 3 Å of site I (Fig. S11 A). The differences in interaction of LIN and STE in the R218Q/R221Q systems are further illustrated by occupancies of head groups around R218 and Y268. LIN displays an 8–30-fold difference compared with STE in occupancy within 3 Å of R218 (Fig. S11 B) and 60–300-fold difference in head group occupancy within 3 Å of Y268 (Fig. S11 C). Altogether, these data suggest that LIN binds better to KCNQ1 in the S4 up state than in the S4 down state and that LIN binds better than STE in both states with the largest difference observed for the S4 down state. These data also suggest that Y268 is a critical anchor point for the FA head group and an important determinant of FA selectivity at site I.

## Discussion

Our results here clearly suggest that there are two functionally important PUFA binding sites present in KCNQ1, one close to the voltage sensor and one close to the pore. In both regions, the negatively charged head group of the PUFA has long-lived and frequent contacts with a positively charged KCNQ1 residue, R218 in the voltage sensor and K316 in the pore domain. We discovered a novel head group interaction—a hydrogen bond between Y268 and the negative PUFA head group—with KCNQ1 at the voltage sensor site. We also discovered several hydrophobic interactions between the PUFA tail and KCNQ1 residues at both sites. The hydrophobic tail of the PUFA is more mobile and samples contacts with many hydrophobic KCNQ1 residues during one interaction event of the negatively charged PUFA head group with the positively charged KCNQ1 residue.

We propose that it is mainly an electrostatic interaction between the negative PUFA head group and the top positive S4 charges that causes the shift in voltage dependence of activation. The electrostatic force between PUFA and R218 in the down state will attract S4 from the down state to the up state, thereby speeding up activation of the channel. Similarly, the electrostatic force between PUFA and R218 in the up state will try to retain S4 from leaving the up state and go to the down state, thereby slowing deactivation of the channel. In experiments, we observed LIN-induced acceleration of KCNQ1 opening kinetics and slowing down of closing kinetics (Fig. S12). These two effects, speeding up of activation and slowing down of deactivation, will result in a shift in the voltage dependence of activation. We propose that Y268 functions as an anchor point for the PUFA head group, thereby allowing a more immobile PUFA head group to attract the mobile S4 charges toward it. For these effects to work, the PUFAs do not have to be bound to the channel sites at all times. The measurable effects on macroscopic parameters are due to time averaging of the PUFA-induced changes in the probabilities of transitions between conformational states. Even if the PUFAs are not constantly bound to the channel, they will still affect the macroscopic parameters by increasing or decreasing the probability of transitions during the times they are bound. Better binding of PUFA to the S4 up state than to the S4 down state is consistent with PUFA stabilization of the activated S4.

What could be a possible molecular mechanism explaining PUFA impact on the permeation properties of KCNQ1? We previously proposed that the PUFA effect on permeation is due to conformational changes in the pore due to electrostatic interaction between the PUFA head group and K316 in S6 (Liin et al., 2018). These conformational changes are unlikely to be accurately sampled in an equilibrium MD simulation. In a previous simulation study (Liin et al., 2018), we instead used a strong force on K316 in the direction of the putative PUFA localization in the lipid bilayer to speed up these conformational changes. We found that this force induced conformational changes in the pore that would explain an increased K+ conduction observed experimentally. PUFA binding to site II identified in this study places the PUFA head group in a position that is expected to cause electrostatic pulling on K316 in a direction aligning with the previous pulling experiments (Liin et al., 2018). The role of each site for distinct PUFA effects is further supported by our finding that mutation of residues predicted to constitute site I predominantly affect the ability of LIN to shift V50, whereas mutation of residues predicted to constitute site II predominantly affect the ability of LIN to increase Gmax. It should, however, be noted that the effect on Gmax may be influenced also by other factors, such as the intrinsic open probability of a mutant (Liin et al., 2018).

In comparison to the PUFA LIN, we find that the saturated STE interacts much less with KCNQ1. This is most likely due to the more rigid nature of the saturated alkane chain of STE. The saturated tail of STE does not explore as many different conformations in a short time period as the polyunsaturated tail of LIN and therefore does not explore as many stabilizing interactions with KCNQ1 residues as shown here for LIN with KCNQ1 residues. For site II, there are clearly fewer interactions for STE compared with LIN in all simulations. For site I, the apparent difference in protein–lipid interactions between STE and LIN is illustrated by the simulations of the R218Q/R221Q system (pseudo S4 down state). STE interactions with site I are 10 and 40 times less likely and long lasting, respectively, than for LIN. The impairment of STE interactions causes large differences in LIN/STE occupancy at site I in the R218Q/R221Q system. The observed reduction in STE–KCNQ1 interactions in the R218Q/R221Q system may explain the lack of effect of STE on the activation time course in KCNQ1. The mobility enabled by the level of unsaturation in the alkyl tail has previously been shown through NMR studies, where the DHA has been observed to exhibit rapid conformational transitions with short correlation times as well as exceptionally low-order parameters (Eldho et al., 2003; Soubias and Gawrisch, 2007). An enhanced mobility of alkyl tails in membrane environment has also been demonstrated by ab initio calculations in combination with AA MD simulations of DHA. The presence of double bonds in the lipid tails leads to low torsional energy barriers for the rotatable bonds in the DHA chain and explains the differences in mobility between polyunsaturated and saturated lipid chains (Feller, 2008). The differences in LIN versus STE tail interactions with a hydrophobic-aromatic cassette in sites present in KCNQ1 is in general agreement with the experimental and modeling results of Ballweg et al. (2020) and Palermo et al. (2015). The previous studies and our results emphasize the critical importance of chain unsaturation, which impacts mobility of the tail and hence an entropic component of binding.

What could be expected for other types of FAs? Monounsaturated FAs have a tail with mobility properties more similar to those of saturated tails than PUFAs have (Palermo et al., 2015; Leng et al., 2018). We therefore expect that monounsaturated FAs would bind poorly to sites I and II, similarly to STE as shown here. This is in agreement with previous experimental data showing a lack of activating effect of the monounsaturated FA oleic acid on KCNQ1 (Liin et al., 2015). In contrast, we expect that other PUFAs such as DHA or eicosapentaenoic acid or PUFA analogues such as N-AT would bind well overall to these sites, similarly to LIN in this study, due to the more mobile tail. However, there might be some differences in how PUFA interacts with KCNQ1, depending on the spatial location of the double bonds in the PUFA tails, which may play a role in specific interactions with aromatic moieties in proteins (Palermo et al., 2015). This is in agreement with previous experimental data showing activating effects on KCNQ1 of a large set of PUFAs (including DHA, eicosapentaenoic acid, and N-AT) but with some variability depending on the location of double bonds (Bohannon et al., 2019). When assessing the importance of Y268 for the V50 effect of the PUFA DHA and PUFA analogue N-AT on KCNQ1, we find it to be critical for the DHA effect (Fig. S13; DHA has no effect on hY278F) but not for the N-AT effect (Fig. S13; N-AT has preserved effect on hY278F). This suggests that for PUFAs with carboxyl head groups, Y268 is an important anchor point, while PUFA analogues with other head groups might interact with the channel in a slightly different way. Future studies are needed to resolve the atomistic details of how double bond location in the PUFA tail affects PUFA interactions with the aromatic cassette and how PUFA analogues with different head groups interact with KCNQ1.

Many different mechanisms have been proposed to explain the effects of PUFA on different ion channels. Indirect effects, such as altering the membrane fluidity or affecting second messenger pathways, have been proposed to explain the relative broad selectivity of PUFAs (Boland and Drzewiecki, 2008; Elinder and Liin, 2017). Moreover, different effects of polyunsaturated and saturated FAs on mechanical properties of the membrane may underlie ion channel effects, but also differences in their abilities to make specific interactions with aromatic residues (Lundbæk et al., 2010; Romero et al., 2019; Ballweg et al., 2020). Although we cannot rule out that PUFAs have indirect effects on KCNQ1 function, our results here are consistent with direct PUFA effects on KCNQ1 by binding to two separate sites on KCNQ1. Moreover, our previous experimental data showing that altering the charge of the PUFA head group without altering the PUFA tail or mutating specific residues in KCNQ1 changes the effect of PUFA on KCNQ1 channels (Liin et al., 2015; Larsson et al., 2018) are more consistent with a direct electrostatic interaction between the head group and KCNQ1 residues. Direct PUFA effects on ion channels were proposed earlier. For instance, electrostatic interactions are important for PUFA effects on the Shaker K+ channel (Börjesson et al., 2008; Ottosson et al., 2014; Yazdi et al., 2016), whereas an aromatic interaction is important for PUFA effects on the BK channel (Hoshi et al., 2013; Tian et al., 2016).

However, for many ion channels, the mechanistic basis for direct interaction remains unknown. We anticipate that the approach of combined CG and AA simulations used in this study, which allows an unbiased exploration of interaction across longer timescales as well as capturing finer atomistic details of interaction, will aid in resolving lipid interaction with the diverse ion channels, whose function has been shown experimentally to be modulated by FAs. Altogether, the characterization of the interactions between PUFAs and KCNQ1 performed in this study enables mechanistic understanding of PUFA-mediated activation of KCNQ1. These insights may be used in future development of drugs based on PUFA mechanisms. For instance, KCNQ1 together with its β subunit KCNE1 form voltage-gated potassium channels that generate the important repolarizing IKs current in the heart (Barhanin et al., 1996; Sanguinetti et al., 1996; Nerbonne and Kass, 2005). Because loss-of-function mutations in KCNQ1 or KCNE1 predispose the heart to arrhythmia, development of KCNQ1/KCNE1 channel activators using PUFA binding sites and similar underlying electrostatic mechanisms may be an anti-arrhythmic strategy. However, many important questions remain to be addressed, such as how KCNE1 impacts PUFA binding to KCNQ1 and how PUFAs target other cardiac ion channels.

## Acknowledgments

Jeanne M. Nerbonne served as editor.

Computational resources were provided by the Swedish National Infrastructure for Computing (2017/1-290; 2018/3-316). Anton2 computer time was provided by the Pittsburgh Supercomputing Center through a grant from the National Institutes of Health (R01GM116961). The Anton2 machine at the Pittsburgh Supercomputing Center was generously made available by D.E. Shaw Research. The computational support for this work was partially provided by Compute Canada through a resource allocation award to S.Yu. Noskov and D.P. Tieleman. This work was supported by the Swedish Society for Medical Research (to S.I. Liin), the Swedish Research Council (2017-02040, to S.I. Liin), the European Research Council under the European Union’s Horizon 2020 research and innovation program (grant agreement no. 850622), and the National Institutes of Health (R01HL131461, to H.P. Larsson). S. Yazdi was supported by the Center for Systems Neurobiology at Linköping University. The work in Calgary was supported by the Canadian Institutes of Health Research Project Grant FRN-CIHR: 156236 (D.P. Tieleman and S.Yu. Noskov) and Discovery grants from the Natural Scientific and Engineering Research Council of Canada (to D.P. Tieleman and S.Yu. Noskov). Further support came from the Canada Research Chairs program (D.P. Tieleman). W. Miranda would like to acknowledge support from a Vanier Canada Graduate scholarship, a Killam scholarship, and an Alberta Innovates Health Solutions studentship.

A patent application (#62/032,739), including a description of the interaction of charged lipophilic compounds with the KCNQ1 channel, has been submitted by the University of Miami with H.P. Larsson and S.I. Liin identified as inventors. The remaining authors declare no competing financial interests.

Author contributions: S. Yazdi, W. Miranda, V. Corradi, D.P. Tieleman, S.Yu. Noskov, H.P. Larsson, and S.I. Liin contributed to study design and analysis approaches. S. Yazdi performed computational simulations. J. Nikesjö and S.I. Liin performed electrophysiological experiments. All authors commented on the manuscript and approved the final version.

## References

Ballweg
,
S.
,
E.
Sezgin
,
M.
Doktorova
,
R.
Covino
,
J.
Reinhard
,
D.
Wunnicke
,
I.
Hänelt
,
I.
Levental
,
G.
Hummer
, and
R.
Ernst
.
2020
.
Regulation of lipid saturation without sensing membrane fluidity
.
Nat. Commun.
11
:
756
.
Barbera
,
N.
,
M.A.A.
Ayee
,
B.S.
Akpa
, and
I.
Levitan
.
2018
.
Molecular dynamics simulations of Kir2.2 interactions with an ensemble of cholesterol molecules
.
Biophys. J.
115
:
1264
1280
.
Barhanin
,
J.
,
F.
Lesage
,
E.
Guillemare
,
M.
Fink
,
M.
Lazdunski
, and
G.
Romey
.
1996
.
K(V)LQT1 and lsK (minK) proteins associate to form the I(Ks) cardiac potassium current
.
Nature.
384
:
78
80
.
Barrese
,
V.
,
J.B.
Stott
, and
I.A.
Greenwood
.
2018
.
KCNQ-encoded potassium channels as therapeutic targets
.
Annu. Rev. Pharmacol. Toxicol.
58
:
625
648
.
Barro-Soria
,
R.
,
S.
Rebolledo
,
S.I.
Liin
,
M.E.
Perez
,
K.J.
Sampson
,
R.S.
Kass
, and
H.P.
.
2014
.
KCNE1 divides the voltage sensor movement in KCNQ1/KCNE1 channels into two steps
.
Nat. Commun.
5
:
3750
.
Bavro
,
V.N.
,
R.
De Zorzi
,
M.R.
Schmidt
,
J.R.
Muniz
,
L.
Zubcevic
,
M.S.
Sansom
,
C.
Vénien-Bryan
, and
S.J.
Tucker
.
2012
.
Structure of a KirBac potassium channel with an open bundle crossing indicates a mechanism of channel gating
.
Nat. Struct. Mol. Biol.
19
:
158
163
.
Berendsen
,
H.J.C.
, and
J.P.M.
Postma
.
1984
.
Molecular dynamics with coupling to an external bath
.
J. Chem. Phys.
81
:
3684
3690
.
Bohannon
,
B.M.
,
M.E.
Perez
,
S.I.
Liin
, and
H.P.
.
2019
.
ω-6 and ω-9 polyunsaturated fatty acids with double bonds near the carboxyl head have the highest affinity and largest effects on the cardiac IKs potassium channel
.
Acta Physiol. (Oxf.).
225
:e13186.
Bohannon
,
B.M.
,
A.
de la Cruz
,
X.
Wu
,
J.J.
Jowais
,
M.E.
Perez
,
D.M.
Dykxhoorn
,
S.I.
Liin
, and
H.P.
.
2020
.
Polyunsaturated fatty acid analogues differentially affect cardiac NaV, CaV, and KV channels through unique mechanisms
.
eLife.
9
:e51453.
Boland
,
L.M.
, and
M.M.
Drzewiecki
.
2008
.
Polyunsaturated fatty acid modulation of voltage-gated ion channels
.
Cell Biochem. Biophys.
52
:
59
84
.
Börjesson
,
S.I.
,
S.
Hammarström
, and
F.
Elinder
.
2008
.
Lipoelectric modification of ion channel voltage gating by polyunsaturated fatty acids
.
Biophys. J.
95
:
2242
2253
.
Broomand
,
A.
, and
F.
Elinder
.
2008
.
Large-scale movement within the voltage-sensor paddle of a potassium channel-support for a helical-screw motion
.
Neuron.
59
:
770
777
.
Cordero-Morales
,
J.F.
,
L.G.
Cuello
,
Y.
Zhao
,
V.
Jogini
,
D.M.
Cortes
,
B.
Roux
, and
E.
Perozo
.
2006
.
Molecular determinants of gating at the potassium-channel selectivity filter
.
Nat. Struct. Mol. Biol.
13
:
311
318
.
Daura
,
X.
,
K.
,
B.
Jaun
,
D.
Seebach
,
W.F.
van Gunsteren
, and
A.E.
Mark
.
1999
.
Peptide folding: When simulation meets experiment
.
Angew. Chem. Int. Ed.
38
:
236
240
.
de Jong
,
D.H.
,
S.
Baoukina
,
H.I.
Ingólfsson
, and
S.J.
Marrink
.
2016
.
Martini straight: Boosting performance using a shorter cutoff and GPUs
.
Comput. Phys. Commun.
199
:
1
7
.
de Jong
,
D.H.
,
X.
Periole
, and
S.J.
Marrink
.
2012
.
Dimerization of amino acid side chains: Lessons from the comparison of different force fields
.
J. Chem. Theory Comput.
8
:
1003
1014
.
de Jong
,
D.H.
,
G.
Singh
,
W.F.
Bennett
,
C.
Arnarez
,
T.A.
Wassenaar
,
L.V.
Schäfer
,
X.
Periole
,
D.P.
Tieleman
, and
S.J.
Marrink
.
2013
.
Improved parameters for the Martini coarse-grained protein force field
.
J. Chem. Theory Comput.
9
:
687
697
.
Eldho
,
N.V.
,
S.E.
Feller
,
S.
Tristram-Nagle
,
I.V.
Polozov
, and
K.
Gawrisch
.
2003
.
Polyunsaturated docosahexaenoic vs docosapentaenoic acid-differences in lipid matrix properties from the loss of one double bond
.
J. Am. Chem. Soc.
125
:
6409
6421
.
Elinder
,
F.
, and
S.I.
Liin
.
2017
.
Actions and mechanisms of polyunsaturated fatty acids on voltage-gated ion channels
.
Front. Physiol.
8
:
43
.
Essmann
,
U.
,
L.
Perera
,
M.L.
Berkowitz
,
T.
Darden
,
H.
Lee
, and
L.G.
Pedersen
.
1995
.
A smooth particle mesh Ewald method
.
J. Chem. Phys.
103
:
8577
8593
.
Farag
,
N.E.
,
D.
Jeong
,
T.
Claydon
,
J.
Warwicker
, and
M.R.
Boyett
.
2016
.
Polyunsaturated fatty acids inhibit Kv1.4 by interacting with positively charged extracellular pore residues
.
Am. J. Physiol. Cell Physiol.
311
:
C255
C268
.
Feller
,
S.E.
2008
.
Acyl chain conformations in phospholipid bilayers: a comparative study of docosahexaenoic acid and saturated fatty acids
.
Chem. Phys. Lipids.
153
:
76
80
.
Hamilton
,
K.L.
,
C.A.
Syme
, and
D.C.
Devor
.
2003
.
Molecular localization of the inhibitory arachidonic acid binding site to the pore of hIK1
.
J. Biol. Chem.
278
:
16690
16697
.
Hess
,
B.
,
H.
Bekker
,
H.J.C.
Berendsen
, and
J.G.E.M.
Fraaije
.
1997
.
LINCS: A linear constraint solver for molecular simulations
.
J. Comput. Chem.
18
:
1463
1472
.
Hoshi
,
T.
,
R.
Xu
,
S.
Hou
,
S.H.
Heinemann
, and
Y.
Tian
.
2013
.
A point mutation in the human Slo1 channel that impairs its sensitivity to omega-3 docosahexaenoic acid
.
J. Gen. Physiol.
142
:
507
522
.
Hou
,
P.
,
J.
Eldstrom
,
J.
Shi
,
L.
Zhong
,
K.
McFarland
,
Y.
Gao
,
D.
Fedida
, and
J.
Cui
.
2017
.
Inactivation of KCNQ1 potassium channels reveals dynamic coupling between voltage sensing and pore opening
.
Nat. Commun.
8
:
1730
.
Humphrey
,
W.
,
A.
Dalke
, and
K.
Schulten
.
1996
.
VMD: visual molecular dynamics
.
J. Mol. Graph.
14
:
33
38: 27–28
.
Jarvis
,
R.A.
, and
E.A.
Patrick
.
1973
.
Clustering using a similarity measure based on shared near neighbors
.
IEEE Trans. Comput.
C-22
:
1025
1034
.
Jensen
,
M.O.
,
V.
Jogini
,
D.W.
Borhani
,
A.E.
Leffler
,
R.O.
Dror
, and
D.E.
Shaw
.
2012
.
Mechanism of voltage gating in potassium channels
.
Science.
336
:
229
233
.
Jo
,
S.
,
T.
Kim
,
V.G.
Iyer
, and
W.
Im
.
2008
.
CHARMM-GUI: a web-based graphical user interface for CHARMM
.
J. Comput. Chem.
29
:
1859
1865
.
Jorgensen
,
W.L.
,
J.
Chandrasekhar
,
J.D.
,
R.W.
Impey
, and
M.L.
Klein
.
1983
.
Comparison of simple potential functions for simulating liquid water
.
J. Chem. Phys.
79
:
926
935
.
Klauda
,
J.B.
,
V.
Monje
,
T.
Kim
, and
W.
Im
.
2012
.
Improving the CHARMM force field for polyunsaturated fatty acid chains
.
J. Phys. Chem. B.
116
:
9424
9431
.
Kuenze
,
G.
,
A.M.
Duran
,
H.
Woods
,
K.R.
Brewer
,
E.F.
McDonald
,
C.G.
Vanoye
,
A.L.
George
Jr
.,
C.R.
Sanders
, and
J.
Meiler
.
2019
.
Upgraded molecular models of the human KCNQ1 potassium channel
.
PLoS One.
14
:e0220415.
,
J.E.
,
H.P.
, and
S.I.
Liin
.
2018
.
KCNE1 tunes the sensitivity of KV7.1 to polyunsaturated fatty acids by moving turret residues close to the binding site
.
eLife.
7
:e37257.
Lees-Miller
,
J.P.
,
J.O.
Subbotina
,
J.
Guo
,
V.
Yarov-Yarovoy
,
S.Y.
Noskov
, and
H.J.
Duff
.
2009
.
Interactions of H562 in the S5 helix with T618 and S621 in the pore helix are important determinants of hERG1 potassium channel structure and function
.
Biophys. J.
96
:
3600
3610
.
Lenaeus
,
M.J.
,
T.M.
Gamal El-Din
,
C.
Ing
,
K.
,
R.
Pomès
,
N.
Zheng
, and
W.A.
Catterall
.
2017
.
Structures of closed and open states of a voltage-gated sodium channel
.
114
:
E3051
E3060
.
Leng
,
X.
,
J.J.
Kinnun
,
A.T.
Cavazos
,
S.W.
Canner
,
S.R.
Shaikh
,
S.E.
Feller
, and
S.R.
Wassall
.
2018
.
All n-3 PUFA are not the same: MD simulations reveal differences in membrane organization for EPA, DHA and DPA
.
Biochim. Biophys. Acta Biomembr.
1860
:
1125
1134
.
Liin
,
S.I.
,
J.E.
,
R.
Barro-Soria
,
B.H.
Bentzen
, and
H.P.
.
2016
.
Fatty acid analogue N-arachidonoyl taurine restores function of IKs channels with diverse long QT mutations
.
eLife.
5
:e20272.
Liin
,
S.I.
,
M.
Silverå Ejneby
,
R.
Barro-Soria
,
M.A.
Skarsfeldt
,
J.E.
,
F.
Starck Härlin
,
T.
Parkkari
,
B.H.
Bentzen
,
N.
Schmitt
,
H.P.
, and
F.
Elinder
.
2015
.
Polyunsaturated fatty acid analogs act antiarrhythmically on the cardiac IKs channel
.
112
:
5714
5719
.
Liin
,
S.I.
,
S.
Yazdi
,
R.
Ramentol
,
R.
Barro-Soria
, and
H.P.
.
2018
.
Mechanisms underlying the dual effect of polyunsaturated fatty acid analogs on Kv7.1
.
Cell Rep.
24
:
2908
2918
.
Lippert
,
R.A.
,
C.
Predescu
,
D.J.
Ierardi
,
K.M.
Mackenzie
,
M.P.
Eastwood
,
R.O.
Dror
, and
D.E.
Shaw
.
2013
.
Accurate and efficient integration for molecular dynamics simulations at constant temperature and pressure
.
J. Chem. Phys.
139
:164106.
Lundbæk
,
J.A.
,
S.A.
Collingwood
,
H.I.
Ingólfsson
,
R.
Kapoor
, and
O.S.
Andersen
.
2010
.
Lipid bilayer regulation of membrane protein function: gramicidin channels as molecular force probes
.
J. R. Soc. Interface.
7
:
373
395
.
Marrink
,
S.J.
,
V.
,
P.C.T.
Souza
,
H.I.
Ingólfsson
,
D.P.
Tieleman
, and
M.S.P.
Sansom
.
2019
.
Computational modeling of realistic cell membranes
.
Chem. Rev.
119
:
6184
6226
.
Luo
,
Y.
, and
B.
Roux
.
2010
.
Simulation of osmotic pressure in concentrated aqueous salt solutions
.
J. Phys. Chem. Lett
.
1
:
183
189
.
Marrink
,
S.J.
,
H.J.
,
S.
Yefimov
,
D.P.
Tieleman
, and
A.H.
de Vries
.
2007
.
The MARTINI force field: coarse grained model for biomolecular simulations
.
J. Phys. Chem. B.
111
:
7812
7824
.
Martyna
,
G.J.
,
D.J.
Tobias
, and
M.L.
Klein
.
1994
.
Constant-pressure molecular-dynamics algorithms
.
J. Chem. Phys.
101
:
4177
4189
.
Nerbonne
,
J.M.
, and
R.S.
Kass
.
2005
.
Molecular physiology of cardiac repolarization
.
Physiol. Rev.
85
:
1205
1253
.
Nosé
,
S.
1984
.
A molecular dynamics method for simulations in the canonical ensemble
.
Mol. Phys.
52
:
255
268
.
Noskov
,
S.Y.
,
S.
Berneche
, and
B.
Roux
.
2004
.
Control of ion selectivity in potassium channels by electrostatic and dynamic properties of carbonyl ligands
.
Nature
.
431
:
830
834
.
Noskov
,
S.Y.
, and
B.
Roux
.
2008
.
Control of ion selectivity in LeuT: two Na+ binding sites with two different mechanisms
.
J. Mol. Biol.
377
:
804
818
.
Oliver
,
D.
,
C.C.
Lien
,
M.
Soom
,
T.
Baukrowitz
,
P.
Jonas
, and
B.
Fakler
.
2004
.
Functional conversion between A-type and delayed rectifier K+ channels by membrane lipids
.
Science.
304
:
265
270
.
Osteen
,
J.D.
,
R.
Barro-Soria
,
S.
Robey
,
K.J.
Sampson
,
R.S.
Kass
, and
H.P.
.
2012
.
Allosteric gating mechanism underlies the flexible gating of KCNQ1 potassium channels
.
109
:
7103
7108
.
Ottosson
,
N.E.
,
S.I.
Liin
, and
F.
Elinder
.
2014
.
Drug-induced ion channel opening tuned by the voltage sensor charge profile
.
J. Gen. Physiol.
143
:
173
182
.
Palermo
,
G.
,
I.
Bauer
,
P.
Campomanes
,
A.
Cavalli
,
A.
Armirotti
,
S.
Girotto
,
U.
Rothlisberger
, and
M.
De Vivo
.
2015
.
Keys to lipid selection in fatty acid amide hydrolase catalysis: Structural flexibility, gating residues and multiple binding pockets
.
PLOS Comput. Biol.
11
:e1004231.
Parrinello
,
M.
, and
A.
Rahman
.
1981
.
Polymorphic transitions in single-crystals - a new molecular-dynamics method
.
J. Appl. Phys.
52
:
7182
7190
.
Pronk
,
S.
,
S.
Páll
,
R.
Schulz
,
P.
,
P.
Bjelkmar
,
R.
Apostolov
,
M.R.
Shirts
,
J.C.
Smith
,
P.M.
Kasson
,
D.
van der Spoel
, et al
.
2013
.
GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit
.
Bioinformatics.
29
:
845
854
.
Romero
,
L.O.
,
A.E.
Massey
,
A.D.
Mata-Daboin
,
F.J.
Sierra-Valdez
,
S.C.
Chauhan
,
J.F.
Cordero-Morales
, and
V.
Vásquez
.
2019
.
Dietary fatty acids fine-tune Piezo1 mechanical response
.
Nat. Commun.
10
:
1200
.
Sanguinetti
,
M.C.
,
M.E.
Curran
,
A.
Zou
,
J.
Shen
,
P.S.
Spector
,
D.L.
Atkinson
, and
M.T.
Keating
.
1996
.
Coassembly of K(V)LQT1 and minK (IsK) proteins to form cardiac I(Ks) potassium channel
.
Nature.
384
:
80
83
.
Schrodinger
,
L.L.C.
2010
. The PyMOL Molecular Graphics System, Version 1.3r1. https://pymol.org/2/support.html? (Accessed March 23, 2019).
Soubias
,
O.
, and
K.
Gawrisch
.
2007
.
Docosahexaenoyl chains isomerize on the sub-nanosecond time scale
.
J. Am. Chem. Soc.
129
:
6678
6679
.
Sun
,
J.
, and
R.
MacKinnon
.
2017
.
Cryo-EM structure of a KCNQ1/CaM complex reveals insights into congenital long QT syndrome
.
Cell.
169
:
1042
1050.e9
.
Tian
,
Y.
,
M.
Aursnes
,
T.V.
Hansen
,
J.E.
Tungen
,
J.D.
Galpin
,
L.
Leisle
,
C.A.
Ahern
,
R.
Xu
,
S.H.
Heinemann
, and
T.
Hoshi
.
2016
.
Atomic determinants of BK channel activation by polyunsaturated fatty acids
.
113
:
13905
13910
.
Tuckerman
,
M.
,
B.J.
Berne
, and
G.J.
Martyna
.
1992
.
Reversible multiple time scale molecular dynamics
.
J. Chem. Phys.
97
:
1990
2001
.
Westhoff
,
M.
,
J.
Eldstrom
,
C.I.
Murray
,
E.
Thompson
, and
D.
Fedida
.
2019
.
I Ks ion-channel pore conductance can result from individual voltage sensor movements
.
116
:
7879
7888
.
Xiao
,
Y.F.
,
Q.
Ke
,
S.Y.
Wang
,
K.
Auktor
,
Y.
Yang
,
G.K.
Wang
,
J.P.
Morgan
, and
A.
Leaf
.
2001
.
Single point mutations affect fatty acid block of human myocardial sodium channel α subunit Na+ channels
.
98
:
3606
3611
.
Xiao
,
Y.F.
,
D.C.
Sigg
, and
A.
Leaf
.
2005
.
The antiarrhythmic effect of n-3 polyunsaturated fatty acids: modulation of cardiac ion channels as a potential mechanism
.
J. Membr. Biol.
206
:
141
154
.
Yazdi
,
S.
,
M.
Stein
,
F.
Elinder
,
M.
, and
E.
Lindahl
.
2016
.
The molecular basis of polyunsaturated fatty acid interactions with the shaker voltage-gated potassium channel
.
PLOS Comput. Biol.
12
:e1004704.
Zaydman
,
M.A.
,
M.A.
Kasimova
,
K.
McFarland
,
Z.
Beller
,
P.
Hou
,
H.E.
Kinser
,
H.
Liang
,
G.
Zhang
,
J.
Shi
,
M.
Tarek
, and
J.
Cui
.
2014
.
Domain-domain interactions determine the gating, permeation, pharmacology, and subunit modulation of the IKs ion channel
.
eLife.
3
:e03606.