Molecular dynamics simulation of ion channel

We constructed models of the channel based on the crystal structure with code 3UKM, which was prepared at a pH >7.4 and is believed to represent an open conformation. These models were protonated to simulate either pH6.0 or pH7.4 and were embedded in a double bilayer system, following a previously established method^{11 (link)}. It’s important to note that this system was different from the one used for the pKas calculations. In each extracellular-facing leaflet of the bilayer, there were 105 molecules of 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphocholine (POPC) and 10 of 1-palmitoyl-2-oleoylphosphatidylethanolamine (POPE). To partially recreate the negative charge inside biological membranes, the cytoplasmic-facing leaflet contained 10 negatively charged dimyristoylphosphatidylglycerol (DMPG) molecules and 105 POPC molecules only. We show in SI-2 section that this membrane modification is unlikely to significantly affect the pKas of the membrane inserted protein. The numbers of lipids were chosen to ensure that both leaflets had a comparable surface area, considering the structure of the channel. After solubilisation in 150 mmol $\cdot$ L-1 KCl and neutralisation, a typical simulation box contained $~$ 243’000 atoms. For the wild-type (WT) sequence, we performed a total of 18 trajectories at each pH, resulting in a cumulative simulation time of 5.8 μs at pH 7.4 and 4.8 μs at pH6.0. A similar procedure was followed for the H122N mutant, except that in this case, we conducted four simulations, each lasting 300 ns, at each pH value. Therefore, the total trajectory time dedicated to the H122N mutant was 2.4 μs. The all-atom MD simulations were performed with the CHARMM36 force field^{31 (link)} using the GROMACS package version 2021.5^{24 (link)}. The TIP3P water model was used^{28 (link)}. Bond and angle lengths involving hydrogen atoms were constrained using the LINCS algorithm allowing an integration time step of 2 fs. Short-range electrostatics were cut off at 1.2 nm. Van der Waals interactions were calculated explicitly up to 10 Å, beyond which a switch function was used to smoothly switch off the interactions to reach zero at 12 Å. Long-range electrostatic interactions were calculated by the PME algorithm^{32 (link)}. The protein, lipids, and water/ions were coupled separately to a temperature bath at 310 K with the Nose-Hoover method with a time constant of 1.0 ps^{33 (link),34 (link)}. The system pressure was kept constant by semi-isotropic Parrinello-Rahman coupling to a reference value of 1 bar as implemented in the GROMACS suite. GROMACS and in-house Python scripts were used to analyse the data. Time series analyses reveal that the Cα RMSD approached values close to 3 Å over the course of the simulations, but without achieving prefect convergence. More precisely, the interactions studied in this article converged significantly and rapidly, enabling us to incorporate all values after 50 ns of stabilisation in our calculations. Some distances involving contact losses required more time to converge (Fig. S1). To determine the orientation of residue 122 with respect to the main axis of the channel, the structures were aligned along their main axis. The vertical red arrow in Fig. 4C shows this orientation. A second vector links the Cα atom of residue 122 with the center of mass of the selected side chain pairs, which are ND1 and NE2 in the case of histidine and OD1 and ND2 in the case of asparagine. The angle between these two vectors determines the orientation of the residue.