Mechanistic basis of the dynamic response of TWIK1 ionic selectivity to pH

Highly selective for K+ at neutral pH, the TWIK1 channel becomes permeable to Na+ upon acidification. Using molecular dynamics simulations, we identify a network of residues involved in this unique property. Between the open and closed states previously observed by electron microscopy, molecular dynamics simulations show that the channel undergoes conformational changes between pH 7.5–6 involving residues His122, Glu235, Lys246 and Phe109. A complex network of interactions surrounding the selectivity filter at high pH transforms into a simple set of stronger interactions at low pH. In particular, His122 protonated by acidification moves away from Lys246 and engages in a salt bridge with Glu235. In addition, stacking interactions between Phe109 and His122, which stabilize the selectivity filter in its K+-selective state at high pH, disappear upon acidification. This leads to dissociation of the Phe109 aromatic side chain from this network, resulting in the Na+-permeable conformation of the channel.


SI-1. Comparison of the CHARMM PBEQ solver with the H++ automated system
As mentioned in the method section, our strategy to estimate the pKas of protein residues involves short MD trajectories (< 10 ns), from which protein conformation, positions of ions interacting with the protein and membrane characteristics (width and width of phospholipid headgroups) are extracted at times 0, 2, 4, 6, and 8 ns.These frames are then subjected to CHARMM's PBEQ module, which determines the total electrostatic potential (PB equation) within 1.5Å grids in a first step, and 0.5Å grids in a second step.For comparison, we subjected the same frames to the H++ automated system, which is faster and less labor-intensive [1].There is no perfect way of objectively compare these two approaches.Thus, we consider two independent features: • Residue 122: Histidine 122 has been identified by mutagenesis as a key TWIK1 proton sensor in the range covered by this study [2] (Fig. 1B).It can therefore be used as a reference.Work with CHARMM's PBEQ module yielded a pKa of 6.1 ± 0.7.Fig. S6 shows that the calculated pKa of His122 calculated by H++, using the default internal dielectric constant of 10, was 7.6 ± 2.7.Using internal dielectric constants of 20 resp. 4 resulted in His122 pKa values of 8.7 ± 1.0 resp.0. Thus, using H++ we would not have been able to validate His122 as a proton sensor.
• Lysine residues: CHARMM's PBEQ identified several lysine residues with a pKa below 6.5, which led us deprotonate them in our work.In our hands using H++, no lysine had a pKa below 10.8.Fig. S7 illustrates a case where Lys131 interacts with the Asp230 from the neighboring subunit.The PBEQ solver provides pKa values of 3.8 for the lysine and 7.3 for the aspartate, suggesting that a proton may be mediating the interaction between these two residues.Using H++, the pKa values are >12 for the lysine and 0.08 for the aspartate.Thus, the module selected for this work seems to better characterize salt bridges and similar interactions.Notably, in a previous computational work studying TWIK1 and in which the pKas were calculated with H++, no lysine residue were identified [3].These two features support the idea that the CHARMM PBEQ solver is more accurate, justifying its use despite the associated computational costs.

SI-2. Lipid membranes
Biological membranes have a wide variety of lipids and are richer in negative species on the cytoplasmic side.Our approach is based on a compromise between an oversimplified membrane (POPC only) and a "real" one that is difficult to manipulate.Similarly, we ensured that the cytoplasm-facing leaflet was negatively charged.This additional work was carried out only after the pKas had been calculated.To enable the CHARMM PBEQ solver to calculate the pKas of a protein inserted into the membrane, it is possible to specify the thickness of the membrane and the height of the protein relative to the bilayer (insertion level).To our knowledge, it is not possible to specify the lipid composition, but only an overall dielectric constant.In this context, membrane thickness is the main difference between the two membrane systems that is likely to affect the calculations.Thickness influences the overall dielectric constant and, without taking account adaptations in protein structure, a thick membrane will interact with certain residues, whereas these residues may interact with the solvent in the case of a thin membrane.We checked the thickness of the membrane.Fig. S8 shows that the membrane used for pKa calculations (left) and the membrane used for a few randomly selected trajectories (right) are identical.It also shows that membrane thickness neither increases nor decreases over time.The table shows the protonation state, pH7.4 or pH6.0 as described in the method section, whether TWIK1 or TWIK1.H122N was modelled, the membrane-embedded protein in the double bilayer system ('AB' or 'CD'), the replica, the length of the trajectory, and whether the trajectory was used in the RMSD resp.distance time series shown in Fig. S1.As four TWIK1 pH6.0 trajectories lasted only 200 ns, they were not included in the RMSD calculations.

Figure S5 .
Figure S5.Distances between Asp230 and Tyr217 at high and low pH.Distances between the centers of mass of the carboxylic acid functional group and the oxygen atom of the phenolic group.Means (horizontal lines), standard deviation (vertical dashed lines) and the results of an unpaired t-test assuming unequal variance and performed on 18 (WT) resp.8 (H122N) pairs of values are indicated.The results of statistical analyses are presented as follows: *, p< 0.05 (p = 0.018 (H122N)); ***, p<0.001 (p=1.00e - (WT)).

Figure S6 .
Figure S6.pKas calculated with H++.Mean and standard deviations of residues with pKa values close to or within the pH range in which the selectivity of TWIK1 varies are shown.The grey dotted lines indicate this pH range.Frames were extracted at t = 0, 2, 4, 6, and 8 ns of MD simulation.Thus, n = 10, is the number of identical residues in the dimer (2), multiplied by the number of frames extracted from the MD (5).

Figure S7 .
Figure S7.Proton sharing and pKa calculations using CHARMM and H++.Molecular representation showing a salt bridge between Lys131 of one subunit (backbone shown in cyan) and Asp230 of another subunit (in grey) during the short MD simulations used for pKa calculations.The distance between the carboxyl oxygen atom of Asp230 and a hydrogen atom of Lys131 is labelled.Only CHARMM pKa calculations suggest proton sharing in this case.

Figure S8 :
Figure S8: membrane thickness used for pKa calculations and MD trajectories.The distances between the phosphorous atoms of the two leaflets are indicated.The color codes, from light to dark green, reflect the values of the images taken at 0-10, resp.200 ns, where the darkest corresponds to the longest simulation time.n = 10, resp.18 for values extracted from the pKa calculations resp.MD simulations.