Parent dataset and code from: Atomistic Mechanisms of the regulation of small conductance Ca 2+ -activated K + channel (SK2) by PIP2

Small conductance Ca 2+ -activated K + channels (SK, K Ca 2) are gated solely by intracellular microdomain Ca 2+. The channel has emerged as a therapeutic target for cardiac arrhythmias. Calmodulin (CaM) interacts with the CaM binding domain (CaMBD) of the SK channels, serving as the obligatory Ca 2+ sensor to gate the channels. In heterologous expression systems, phosphatidylinositol 4,5-bisphosphate (PIP2) coordinates with CaM in regulating SK channels. However, the roles and mechanisms of PIP2 in regulating SK channels in cardiomyocytes remain unknown. Here, optogenetics, magnetic nanoparticles, combined with Rosetta structural modeling and molecular dynamics (MD) simulations revealed the atomistic mechanisms of how PIP2 works in concert with Ca 2+ -CaM in the SK channel activation. Our computational study affords evidence for the critical role of the amino acid residue R395 in the S6 transmembrane domain, which is localized in propinquity to the intracellular hydrophobic gate. This residue forms a salt bridge with residue E398 in the S6 transmembrane domain from the adjacent subunit. Both R395 and E398 are conserved in all known isoforms of SK channels. Our findings suggest that the binding of PIP2 to R395 residue disrupts the R395:E398 salt bridge, increasing the flexibility of the transmembrane segment S6 and the activation of the channel. Importantly, our findings serve as a new platform for testing structural-based drug designs for therapeutic inhibitors and activators of the SK channel family. The study is timely since inhibitors of SK channels are currently in clinical trials to treat atrial arrhythmias.

Methods
See the PNAS publication for a full description and references

Computational modeling of hSK2 channel to generate starting structures for MD simulations:

The generation of structural models for the hSK2 channel was achieved via three stages using Rosetta molecular modeling (Online Supplemental Fig. S1) as described below.

Electron density map refinement:

In the first stage, we refined and converted three states of the hSK4-CaM cryo-EM structures (PDB IDs: 6CNM, 6CNN, and 6CNO) into Rosetta-optimized energy with Rosetta 2021 software. Cryo-EM refinement was performed with side chains only being optimized. This conversion facilitated the energy terms of these structures to align with Rosetta's scoring function (Score.gd2), leading to improved homology modeling. We used a modified version of the Rosetta demo included in the Rosetta software
(https://new.rosettacommons.org/demos/latest/public/electron_density_structure_refinement/structure_refinement) (7, 8).

Homology modeling of hSK2 channels in closed, intermediate, and open states:

In the second stage, we employed the modified Rosetta Comparative Modeling (RosettaCM) protocol to generate hSK2 homology models (8-15). This involved conducting sequence alignments between hSK4 and hSK2 channels, and subsequently formatting the aligned sequences into a Grishin format suitable for RosettaCM. Regions in hSK2 that were not present in hSK4 structures were modeled using the loop modeling protocol, (https://www.rosettacommons.org/docs/latest/application_documentation/structure_prediction/loop_modeling/KIC_with_fragments) (16) primarily the S3-S4 linker. The template protocol can be found at this link: (https://new.rosettacommons.org/demos/latest/Home). The top model from the cryo-EM refinement of the hSK4 step was then used as the template for homology modeling for hSK2. The first attempt was to only model hSK2, then dock CaM onto the channel. However, the absence of CaM in the models triggered large movements in the CaM binding domain (CaMBD) in the C-terminal domain of hSK2 (main-text Fig. 2A) due to two factors: A) CaMBD is perpendicular to the S1-S6 transmembrane segments of hSK2 and B) the linker region that connects the CaMBD and S1-S6 is very flexible. In contrast, the inclusion of CaM led to convergence and improved agreement in the top models for the CaMBD (main-text Fig. 2B-D), resulting in a model that is much closer to the template. In our attempt to create reliable homology models, several features were included in the protocol, namely implicit lipid membrane environment to accurately model membrane-spanning protein segments, enforcing symmetry to preserve a four-fold homotetrameric hSK2-CaM complex symmetry, explicit inclusion of metal ions to preserve the Ca 2+ binding loops, treatment of multiple chains that include hSK2 and CaM, and loop modeling for the flexible regions of the protein missing from the cryo-EM structures. However, an error occurred between the symmetry function and the metal binding feature, specifically, the Ca 2+ ions assumed the exact coordinates of the first C α atom of the first residue of CaM. This necessitated the homology modeling to be performed without the symmetry function. In addition, we meticulously monitored for possible displacement of the backbone carbonyl oxygen atoms in the selectivity filter (SF) of hSK2 during the homology modeling since deformation of the SF may lead to a non-conducting channel in the MD simulations. Indeed, repulsion of the oxygen atoms results in the SF deformation at amino acid residues I359 with widening and shortening of the SF. This necessitated the inclusion of harmonic restraints that were determined empirically, and a weighted value of 100.0 kcal/mol/Å 2 was used on the C α atoms of the backbone of each amino acid residue in the SF. Additionally, the deformation was minimized with the explicit inclusion of K + ions in the SF. Interestingly, we did observe that the open state of hSK2 required the largest restraints to maintain the SF structure. 50,000 models were created for each conformational state of hSK2-CaM, and a top model was selected with a standard clustering selection process.

Clustering and top model selection:

Standard Rosetta clustering was performed for each cryo-EM refinement and homology modeling step with minor differences in the filtering process, based on the models that were used as inputs. The cryo-EM refinement output structures were filtered first by sorting using the “r15” term (REF2015 terms: SCORE – elec_dens_fast) and keeping the lowest 50% of r15 ranked decoys. The resulting list of decoys was then sorted by elec_dens_fast term, and the lowest 20% were kept. This final selection was then clustered with a radius determined empirically to provide a distribution, where the most decoys (~30%) were in the first cluster, and subsequent clusters were with progressively fewer decoys. An attempt to obtain upwards of 90% of all models in the first 20 clusters was made. The centers of the top 10 clusters were then evaluated for retention of similar critical structural features described above. The top model from these 10 was used as a template for the homology modeling. The results of the homology modeling were sorted by the total score term with the lowest 10% used for clustering. Clustering of the homology models was performed in the same manner as the cryo-EM refinement with the top model used as the starting structure for MD simulations.

Molecular dynamics (MD) simulations:

The final stage involved molecular dynamics (MD) simulations. The hSK2 models derived from homology modeling were initially visualized without a membrane and then embedded into 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipid bilayer and solvated by 0.15 M KCl aqueous solution using CHARMM-GUI (17-19) (Fig. 3A). Online Supplemental Table S1 provides a summary of 18 5-µs-long MD simulations on Anton 2 supercomputer (20) of hSK2-CaM complex in POPC membrane with or without mono-protonated state of phosphatidylinositol -(4,5)-bisphosphate (PIP2) with protonation on P4 oxygen atom (SAPI24) at 2.5, 5, and 10% in the lower leaflet of the lipid bilayer. The concentrations of PIP2 were chosen based on recent estimations of PIP2 in the lower leaflet of the lipid bilayer that can be as high as 2-5% (21). In addition, we performed a total of 27 simulations at 1 µs using either NAMD 3.0 alpha (22) or AMBER18 (23) on the high-performance computing (HPC) EXPANSE platform (San Diego Supercomputer Center at the University of California, San Diego) with computational time granted through Extreme Science and Engineering Discovery Environment, XSEDE (now Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support, ACCESS). MD simulations were run in the NPT ensemble at 310 K and 1 atm pressure using tetragonal periodic boundary conditions using a standard set of non-bonded cutoffs and other options as in our previous studies (24, 25). All-atom biomolecular CHARMM36m protein (26), C36 lipid (27, 28), and TIP3 water (29) were used. Each MD simulation system was equilibrated for 2.27 ns with suggested gradually diminishing positional and dihedral restraints provided by CHARMM-GUI scripts. Due to the relatively large size of the hSK2-CaM protein complex in our MD simulations and to ensure its conformational stability (as described above), a follow-up extended equilibration MD simulation was performed for 100 ns with gradually reduced restraints on protein backbone atoms of the hSK2-CaM complex and its components as shown in main-text Fig. 3B. The extended equilibration protocol was performed on the EXPANSE platform with AMBER18 (Fig. 3B). The protein backbone restraints were maintained using the force constant of 1.0 kcal·mol -1 ·Å -2 during the initial environment equilibration until the start of the extended equilibration stage, when the restraints were successively reduced in a 5 or 10 ns stepwise fashion, starting from the periphery of the protein and moving towards the center of the protein over a period of 100 ns (main-text Fig. 3B). The extended equilibration protocol was determined empirically to maintain the stability of essential structural features such as the pore domain or SF. Each successive step of either 5 or 10 ns was determined based on root-mean-square deviation (RMSD) profiles reaching a plateau indicating system equilibration. Post 100 ns of this extended equilibration, the production MD simulation runs were performed for 1000 ns using AMBER18 on EXPANSE or for 5000 ns on Anton 2. Small restraints were used on the C α backbone atoms of the selectivity filter residues SIGYGD throughout the production stage (Fig. 3B). The 27 MD simulations created were divided as follows: 3 distinct hSK2 conformational states (closed, intermediate, and open) embedded into POPC:PIP2 complex membranes with PIP2 residing in the lower leaflet only (100:0, 95:5, 90:10). All simulations were run in triplicate to ensure consistency. The 18 most stable simulations of the 27 runs using AMBER18 were identified and transferred to Anton 2 by using molecular coordinates and velocities from the final frame of the 100 ns protein equilibration stage as a starting point. The MD simulations on Anton 2 are unbiased 5-µs-long runs and are described in Online Supplemental Table S1. The MD simulation trajectories were analyzed via root-mean-square deviation (RMSD) calculations (main-text Fig. 3C).

Pore volume representation and initial state characterization:

The distinct initial channel states were characterized by central channel pore radius mapping using HOLE analysis script (30) (main-text Fig. 4A-C). We thus monitored channel conformational states by measuring the minimum hydrophobic gate diameter for each frame (main-text Fig. 4D-E). This was performed by measuring the distance between side-chain C g1 (CG1) atoms of the V390 residues on the pore-lining S6 helices from opposing subunits.

Ion conduction:

Ion conduction was measured by monitoring the Z-axis coordinates of K + ions through the SF (main-text Fig. 4F-G). We then quantified the number of ions that completed a connected path through the SF to avoid counting ions that would “jump” from the intracellular aqueous compartment to the extracellular one using periodic boundary conditions (PBC) image recentering.

PIP2 movement and specific residues of SK2-CaM complex involved in binding:

This was determined by tracking its XY coordinates and plotting 50 points per graph for each PIP2 molecule. The identification of all possible unbiased PIP2 binding sites was achieved using modified scripts from MD analysis GitHub (https://github.com/MDAnalysis/mdanalysis), (31) and amino acid residues of interest were located by identifying salt bridges formation between the hSK2-CaM complex and PIP2 head group. 

Visual, structural, and computational analysis:

Completed using Visual Molecular Dynamics (VMD) (32), UCSF Chimera (33), and ChimeraX (34) (University of California San Francisco). 

Clustering Analysis:

After the completion of the MD simulations, we conducted a clustering analysis of the simulation data. The primary objective of this analysis was to categorize conformations based on their similarity and to pinpoint the most frequently occurring conformational states of hSK2. Additionally, this analysis aimed to investigate the coordination of hSK2 and PIP2 clusters and the formation of the R395:PIP2 salt bridge. To achieve these goals, we employed the TTClust program, a specialized tool specifically designed for trajectory clustering (35) (https://github.com/tubiana/TTClust). We first aligned the PIP2 molecule with respect to the hSK2-CaM. The trajectory of the aligned PIP2 molecule was saved as a binary DCD trajectory file. Subsequently, we executed the TTClust program using the command: ttclust -f HETA-c.dcd -t HETA-c.pdb -sa "none" -sr "all". The -sa parameter was set to "none" to indicate no further alignment was required as it was already performed in the first step. The -sr parameter was set to "all" to select all atoms in the analysis. We then analyzed the clustering results to identify clusters with PIP2:R395 salt bridges, focusing on the major residues involved. We visualized and compared these clusters against the hSK2-CaM structure to observe their distribution and identify their location. Furthermore, we calculated the time point at which a salt bridge was detected via a 3.6 Å cutoff. 

Solvent-Accessible Surface Area (SASA):

Analyses Solvent-accessible surface area (SASA) analyses were performed to quantify the solvent exposure of each amino acid residue in the transient, transfer, and activation PIP2 binding sites from the closed, intermediate, and open states. There were no significant differences in the SASA over time or averaged SASA for amino acid residues within the same state for MD simulations with the hSK2-CaM embedded in a POPC or a POPC/PIP2 membrane if a PIP2 molecule did not directly interact with that amino acid residue. Therefore, we calculated the SASA for each amino acid residue on the four subunits for each frame based on binary classification: bound vs. unbound by PIP2 (the distance between any anionic oxygen on the PIP2 head group and any cationic nitrogen on the side chain of basic amino acids was <4 Å was defined as bound and >4 Å was defined as unbound). For the open hSK2 state, we further classified the MD simulations into 2 additional groups with an intact hydrophobic gate and a collapsed hydrophobic gate, which became nonconductive after the initial K + ions in the channel pore were depleted. These simulations were excluded from the analyses since there was major asymmetry of the pore.