Using Gromacs to do computational electrophysiology simulations
With this add-on to Gromacs 4.5, ionic concentrations in two separate compartments of a simulation volume can be kept at constant values each. This is achieved by exchanging the ions from one compartment with water molecules from the other compartment if the concentrations deviate from the provided reference values. You will find more information on the method in this article [1].
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:
git clone git://git.gromacs.org/gromacs.git
cd gromacs
git checkout --track -b comp_el origin/comp_el
If you use autotools, you further need to run ./bootstrap in the Gromacs directory. Then, you can configure and install just as you would do with a regular Gromacs.
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.
When an exchange needs to be performed, the position of the ion is exchanged with the center of mass of a suitable water molecule from the other compartment. Out of all ions in a given compartment, the ion with the largest distance to any of the compartment boundaries is selected for the exchange. The same criterion is used for the water molecules. The swapions.xvg output file records the exchanges that have been performed (which on time-average are equal to the ion flux between the compartments), as well as ion permeations through the individual channels.
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, autoswapcoords = 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 frequencyswap_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 atomssplit_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 moleculesswap_group = NA+_CL-
Provide the index group with the exchangeable ions.
; Group name of solvent moleculessolvent_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 stepscoupl_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 reachedthreshold = 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.