Using Gromacs to do computational electrophysiology simulations

Downloading the code

The code can be downloaded from the official Gromacs GIT repository. Download instructions can be found at www.gromacs.org. Here is a quick version of what you need to do to get the most recent version of the Gromacs 4.5 Computational Electrophysiology patch:

Usage instructions

You need to:

• define two compartments by giving two index groups (the so-called split groups),
• define a group of ions as well as a group of solvent molecules,
• define the reference number of anions and cations for each compartment.

For the latter functionality, a 'split cylinder' must be defined around the split group's center by providing a cylinder radius and an upper and lower extension. Ions that pass this cylinder are counted as having passed this channel. Note that when starting a new simulation it may take some time until these counters show ion passages since only ions can be counted that really have changed their compartment, i.e. A->channel(0/1)->B or B->channel(0/1)->A. Ions that reside inside a channel at the start of a simulation can not be counted as having passed the channel, since we can not know from which compartment they originated. This is however only an issue at the very start of a simulation, for restarts we keep the history of every ion in the checkpoint file.

Preprocessing

You can control ion/water position exchanges with the following input file options:

; Swap coordinates: no, X, Y, Z, auto
swapcoords               = Z

Choose Z if your membrane is in the X-Y-plane. Ions will be swapped depending on their Z-positions alone. Note that the 'auto' option for heavily tilted membranes or vesicles is not yet implemented.

; Swap attempt frequency
swap_frequency           = 100

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.

; Two index groups that contain the compartment-partitioning atoms
split_group0             = channelA
split_group1             = channelB
; Use center of mass of split groups (yes/no), otherwise geometrical center is used
massw_split0             = no
massw_split1             = no

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.

; Group name of ions that can be exchanged with solvent molecules
swap_group               = NA+_CL-

Provide the index group with the exchangeable ions.

; Group name of solvent molecules
solvent_group            = SOL

Provide the group of exchangeable solvent molecules.

; Split cylinder: radius, upper and lower extension [nm] (this will define the channels)
cyl0_r                   = 5.0
cyl0_up                  = 0.75
cyl0_down                = 0.75
cyl1_r                   = 5.0
cyl1_up                  = 0.75
cyl1_down                = 0.75

; Average the number of ions per compartment over these many swap attempt steps
coupl_steps              = 10

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.

; Requested number of anions and cations for each of the two compartments
; -1 means fix the numbers as found in time step 0
anions0                  = -1
cations0                 = -1
anions1                  = -1
cations1                 = -1

Set the reference number of anions and cations for each of the compartments. They have to add up to the total number of ions present in the swap group.

; Start to swap ions if threshold difference to requested count is reached
threshold                = 1

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.

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.

References

1. Kutzner C, Grubmüller H, de Groot BL, Zachariae U. Computational Electrophysiology: The Molecular Dynamics of Ion Channel Permeation and Selectivity in Atomistic Detail. Biophys J. 101: 809-817 (2011) [pdf] [New&Notable]