Max Planck Institute for Multidisciplinary Sciences
Studying ion channel permeation and selectivity with »Computational Electrophysiology«
Figure 1: Simulation Setup. (A) Using periodic boundaries, a double membrane partitions the system into two compartments α and β. (B) Simulation system of PorB in a lipid membrane at 200 mM NaCl concentration. (C) Time average of the electrostatic potential U along the z-axis for wild type PorB at 1 M NaCl, resulting from charge imbalances ∆q of up to 12 elementary charges.
Figure 1: Simulation Setup. (A) Using periodic boundaries, a double membrane partitions the system into two compartments α and β. (B) Simulation system of PorB in a lipid membrane at 200 mM NaCl concentration. (C) Time average of the electrostatic potential U along the z-axis for wild type PorB at 1 M NaCl, resulting from charge imbalances ∆q of up to 12 elementary charges.
Ion channels play an essential role in cellular control and signalling processes. Ion selectivity and permeation rates are key quantities measured in experiments. We show that these quantities can be directly measured in molecular dynamics (MD) simulations by establishing a preserved ion gradient across the membrane, which leads to sustained ion flux through the channel. To preserve a chosen ionic gradient, we exchange ion/water pairs between both sides of the membrane.
We demonstrate the validity of this method using the bacterial channel PorB from pathogenic Neisseria meningitidis. During Neisserial infection, PorB inserts into the mitochondrial membrane of target cells, where it leads to apoptosis by dissipating the membrane potential. Our approach accurately predicts PorB ion conductance and selectivity and elucidates conduction mechanisms in atomistic detail.
Figure 2: Ion flux through the PorB channels as recorded in the 1 M NaCl PorB simulation.
Every 0.2 ps, the number of anions and cations is determined in the compartments α and β (Fig. 1) and, if needed, the chosen reference numbers of ions are restored by exchanging ion/water pairs, thus providing a fixed charge imbalance ∆q between α and β.
Calculation of ion selectivity and conductance
The transmembrane potential difference ∆U(z) is calculated by integrating the charge distribution twice (Fig. 1C). The conductance is G = I/∆U, the selectivity Ianion/Ication.
Figure 3: The ionic current can be computed using the slopes of the flux curves (Fig. 2). (A) Net ion flux at the end of the simulation. (B) Ionic current as a function of ∆U for multiple 20-ns-long slices of trajectories shown in Fig. 2.
Figure 3: The ionic current can be computed using the slopes of the flux curves (Fig. 2). (A) Net ion flux at the end of the simulation. (B) Ionic current as a function of ∆U for multiple 20-ns-long slices of trajectories shown in Fig. 2.
To obtain an error estimate, and since potential and flux are time-dependent, multiple ∆U(t) and I(t) are determined from 20 ns time windows (10 ns overlap). For each window, ∆U(t) is determined as shown in Fig. 1C, while I(t) is derived by fitting a straight line to the ion flux in that time window. This way, each data point from Fig. 3A results in a cluster of data points in Fig. 3B. Combining the results from the three different ∆q values and then fitting linear functions to those data gives the straight lines shown in Fig. 3B. The selectivity ratio (Cl−/Na+) is provided by the slope ratios.
The conductance of wild type PorB was calculated to be 0.8 ± 0.1 nS, with an anion selectivity of 3.0 ± 0.6 (experiment (Olesky et al., 2006): 1.0 nS, selectivity 2.8, both in 200 mM NaCl)—an improvement over potential of mean force calculations (Fig. 4B).
Figure 4: PorB ion pathway. (A) Overlay of ion positions from 500 snapshots of a 100 ns simulation at ∆q=8 e. (B) Potential of mean force for Cl- and Na+ ions for wild type PorB. (C) as A, but for PorB mutant G103K, which shows a disrupted cation pathway.
Figure 4: PorB ion pathway. (A) Overlay of ion positions from 500 snapshots of a 100 ns simulation at ∆q=8 e. (B) Potential of mean force for Cl- and Na+ ions for wild type PorB. (C) as A, but for PorB mutant G103K, which shows a disrupted cation pathway.
The ion transfer mechanism inside PorB has been suggested to be unusual. Individual, hardly interacting pathways for the diffusion of anions and cations were postulated based on its crystal structure (Tanabe et al., 2010). As shown by Fig. 4A, we find that there is indeed virtually no overlap of the pathway that cations and anions take along the pore. Anions pass along a cluster of basic residues on one side of the eyelet, while cations move along acidic residues lining the opposite wall.
Summary
Computational electrophysiology allows a direct estimation of ion channel conductance and selectivity. Simultaneously, it provides mechanistic insight into channel function. As the method is applicable to channels of small and large conductance and size, we expect it to be useful to study the molecular mechanism of ion passage, e.g. for the improvement of drugs to overcome resistance mutations in bacterial porins, and to design drugs for a range of ion channel targets.
The "Computational Electrophysiology" protocol is integrated in the official 5.0 release of GROMACS, so you can use it with any 5.0 or later version.
Check out the "Computational Electrophysiology" section in the GROMACS .pdf manual for more information.
Kutzner, C.; Köpfer, D.; Machtens, J. P.; de Groot, B. L.; Song, C.; Zachariae, U.: Insights into the function of ion channels by computational electrophysiology simulations. Biochimica et Biophysica Acta - Biomembranes 1858 (7, B), pp. 1741 - 1752 (2016)
Kutzner, C.; Grubmüller, H.; de Groot, B. L.; Zachariae, U.: Computational Electrophysiology: The molecular dynamics of ion channel Permeation and selectivity in atomistic detail. Biophysical Journal 101 (4), pp. 809 - 817 (2011)
Köpfer, D.; Song, C.; Gruene, T.; Sheldrick, G. M.; Zachariae, U.; de Groot, B. L.: Ion permeation in K+ channels occurs by direct Coulomb knock-on. Science 345 (6207), pp. 352 - 355 (2014)
Song, C.; Weichbrodt, C.; Salnikov, E.; Dynowski, M.; Forsberg, B.; Bechinger, B.; Steinem, C.; de Groot, B. L.; Zachariae, U.; Zeth, K.: Crystal structure and functional mechanism of a human antimicrobial membrane channel. Proceedings of the National Academy of Sciences of the United States of America 110 (12), pp. 4586 - 4591 (2013)
Exemplary GROMACS .mdp file entries controlling the deterministic protocol
You can control ion/water position exchanges with the following input file options:
; Swap coordinates: no, X, Y, Z. Choose Z if your membrane is in the X-Y-plane. Ions will be swapped depending on their Z-positions alone. swapcoords = Z
; Swap attempt frequency. This determines how often a swap attempt will be made. Since therefore the positions of the ions, solvent, and swap groups are communicated around, this is a time-consuming operation. Do not try to swap every step, this will slow down the simulation a lot. swap-frequency = 100
; Two index groups that contain the compartment-partitioning atoms split-group0 = channel0 split-group1 = channel1
; Use center of mass of split groups (yes/no), otherwise geometrical center is used. Choose two index groups that divide the MD system into two compartments, here in Z-direction. If massw-split is activated, this group's center of mass will be used as the dividing point, if not, the geometrical center is used. If you choose a membrane channel as split group, then the center of the channel will define the compartment-dividing plane. massw-split0 = no massw-split1 = no
; Group name of solvent molecules solvent-group = SOL
; Average the number of ions per compartment over these many swap attempt steps. If coupl-steps is set to 1, then the instantaneous ion distribution will determine whether ions are exchanged. coupl-steps > 1 will use the time-averaged ion distribution instead. This is useful when ions are diffusing around near compartment boundaries (in the channel for example) which would lead to numerous in- and outswaps for coupl-steps=1. coupl-steps = 10
; Number of ion types to be controlled iontypes = 3
; Names of the ion types that can be exchanged with solvent molecules ; -1 means fix the numbers as found in time step 0. These numbers have to add up to the total number of ions present in the swap group. iontype0-name = NA iontype0-in-A = 51 ; requested number of NA ions in compartment A iontype0-in-B = 35 ; requested number of NA ions in compartment B iontype1-name = K iontype1-in-A = 10 ; requested number of K ions in compartment A iontype1-in-B = 38 ; requested number of K ions in compartment B iontype2-name = CL iontype2-in-A = -1 ; -1 means: use the number of ions iontype2-in-B = -1 ; as found at time step 0
; Split cylinder: radius, upper and lower extension (nm) (for counting on-the-fly diagnostics, has no influence on whether or not ions are swapped) cyl0-r = 5.0 cyl0-up = 0.75 cyl0-down = 0.75 cyl1-r = 5.0 cyl1-up = 0.75 cyl1-down = 0.75
; Start to swap ions if threshold difference to requested count is reached. A threshold of 1 means that a swap is performed if the average ion count in a compartment differs by 1 ore more from the requested values. Higher thresholds mean that larger differences are accepted. Ions are also only swapped until the requested number +/- the threshold is reached.
threshold = 1
Output
Apart from the usual mdrun output files, there is a file called swapions.xvg which can be directly plotted with xmgrace and which contains information about the ion counts per compartment, the differences to the reference values, and the number of swapped-in and swapped-out ions.
Build a double-layered system from a single-layered one
helps you to set up your index groups for the double layered system if you already have index groups for a single layered system.
Remarks about which atoms passed which channel is printed to the swapions.xvg output file (using global atoms numbers starting at one). A warning is also output if an ion is causing leakage because it is in the membrane (again, with atom number so you can check in a PDB viewer).
grompp outputs the initial, whole channel configuration to .pdb file. If in this file your channels do not have the correct (whole) PBC representation, ion permeations can not be counted properly! You will need to correct your channels to have the proper ("unbroken") PBC representation.
