.. _compel: Computational Electrophysiology ------------------------------- The Computational Electrophysiology (CompEL) protocol :ref:`147 ` allows the simulation of ion flux through membrane channels, driven by transmembrane potentials or ion concentration gradients. Just as in real cells, CompEL establishes transmembrane potentials by sustaining a small imbalance of charges :math:`\Delta q` across the membrane, which gives rise to a potential difference :math:`\Delta U` according to the membrane capacitance: .. math:: \Delta U = \Delta q / C_{membrane} :label: eqnmembcap The transmembrane electric field and concentration gradients are controlled by :ref:`mdp` options, which allow the user to set reference counts for the ions on either side of the membrane. If a difference between the actual and the reference numbers persists over a certain time span, specified by the user, a number of ion/water pairs are exchanged between the compartments until the reference numbers are restored. Alongside the calculation of channel conductance and ion selectivity, CompEL simulations also enable determination of the channel reversal potential, an important characteristic obtained in electrophysiology experiments. In a CompEL setup, the simulation system is divided into two compartments **A** and **B** with independent ion concentrations. This is best achieved by using double bilayer systems with a copy (or copies) of the channel/pore of interest in each bilayer (:numref:`Fig. %s ` A, B). If the channel axes point in the same direction, channel flux is observed simultaneously at positive and negative potentials in this way, which is for instance important for studying channel rectification. .. _fig-compelsetup: .. figure:: plots/compelsetup.* :width: 13.50000cm Typical double-membrane setup for CompEL simulations (A, B). Ion/water molecule exchanges will be performed as needed between the two light blue volumes around the dotted black lines (A). Plot (C) shows the potential difference :math:`\Delta U` resulting from the selected charge imbalance :math:`\Delta q_{ref}` between the compartments. The potential difference :math:`\Delta U` across the membrane is easily calculated with the :ref:`gmx potential ` utility. By this, the potential drop along :math:`z` or the pore axis is exactly known in each time interval of the simulation (:numref:`Fig. %s ` C). Type and number of ions :math:`n_i` of charge :math:`q_i`, traversing the channel in the simulation, are written to the swapions.xvg output file, from which the average channel conductance :math:`G` in each interval :math:`\Delta t` is determined by: .. math:: G = \frac{\sum_{i} n_{i}q_{i}}{\Delta t \, \Delta U} \, . :label: eqnchannelcond The ion selectivity is calculated as the number flux ratio of different species. Best results are obtained by averaging these values over several overlapping time intervals. The calculation of reversal potentials is best achieved using a small set of simulations in which a given transmembrane concentration gradient is complemented with small ion imbalances of varying magnitude. For example, if one compartment contains 1M salt and the other 0.1M, and given charge neutrality otherwise, a set of simulations with :math:`\Delta q = 0\,e`, :math:`\Delta q = 2\,e`, :math:`\Delta q = 4\,e` could be used. Fitting a straight line through the current-voltage relationship of all obtained :math:`I`-:math:`U` pairs near zero current will then yield :math:`U_{rev}`. Usage ^^^^^ The following :ref:`mdp` options control the CompEL protocol: :: swapcoords = Z ; Swap positions: no, X, Y, Z swap-frequency = 100 ; Swap attempt frequency Choose ``Z`` if your membrane is in the :math:`xy`-plane (:numref:`Fig. %s `). Ions will be exchanged between compartments depending on their :math:`z`-positions alone. ``swap-frequency`` determines how often a swap attempt will be made. This step requires that the positions of the split groups, the ions, and possibly the solvent molecules are communicated between the parallel processes, so if chosen too small it can decrease the simulation performance. The ``Position swapping`` entry in the cycle and time accounting table at the end of the ``md.log`` file summarizes the amount of runtime spent in the swap module. :: split-group0 = channel0 ; Defines compartment boundary split-group1 = channel1 ; Defines other compartment boundary massw-split0 = no ; use mass-weighted center? massw-split1 = no ``split-group0`` and ``split-group1`` are two index groups that define the boundaries between the two compartments, which are usually the centers of the channels. If ``massw-split0`` or ``massw-split1`` are set to ``yes``, the center of mass of each index group is used as boundary, here in :math:`z`-direction. Otherwise, the geometrical centers will be used (:math:`\times` in :numref:`Fig. %s ` A). If, such as here, a membrane channel is selected as split group, the center of the channel will define the dividing plane between the compartments (dashed horizontal lines). All index groups must be defined in the index file. If, to restore the requested ion counts, an ion from one compartment has to be exchanged with a water molecule from the other compartment, then those molecules are swapped which have the largest distance to the compartment-defining boundaries (dashed horizontal lines). Depending on the ion concentration, this effectively results in exchanges of molecules between the light blue volumes. If a channel is very asymmetric in :math:`z`-direction and would extend into one of the swap volumes, one can offset the swap exchange plane with the ``bulk-offset`` parameter. A value of 0.0 means no offset :math:`b`, values :math:`-1.0 < b < 0` move the swap exchange plane closer to the lower, values :math:`0 < b < 1.0` closer to the upper membrane. :numref:`Fig. %s ` A (left) depicts that for the **A** compartment. :: solvent-group = SOL ; Group containing the solvent molecules iontypes = 3 ; Number of different ion types to control iontype0-name = NA ; Group name of the ion type iontype0-in-A = 51 ; Reference count of ions of type 0 in A iontype0-in-B = 35 ; Reference count of ions of type 0 in B iontype1-name = K iontype1-in-A = 10 iontype1-in-B = 38 iontype2-name = CL iontype2-in-A = -1 iontype2-in-B = -1 The group name of solvent molecules acting as exchange partners for the ions has to be set with ``solvent-group``. The number of different ionic species under control of the CompEL protocol is given by the ``iontypes`` parameter, while ``iontype0-name`` gives the name of the index group containing the atoms of this ionic species. The reference number of ions of this type can be set with the ``iontype0-in-A`` and ``iontype0-in-B`` options for compartments **A** and **B**, respectively. Obviously, the sum of ``iontype0-in-A`` and ``iontype0-in-B`` needs to equal the number of ions in the group defined by ``iontype0-name``. A reference number of ``-1`` means: use the number of ions as found at the beginning of the simulation as the reference value. :: coupl-steps = 10 ; Average over these many swap steps threshold = 1 ; Do not swap if < threshold If ``coupl-steps`` is set to 1, then the momentary ion distribution determines whether ions are exchanged. ``coupl-steps > 1`` will use the time-average of ion distributions over the selected number of attempt steps instead. This can be useful, for example, when ions diffuse near compartment boundaries, which would lead to numerous unproductive ion exchanges. A ``threshold`` of 1 means that a swap is performed if the average ion count in a compartment differs by at least 1 from the requested values. Higher thresholds will lead to toleration of larger differences. Ions are exchanged until the requested number :math:`\pm` the threshold is reached. :: cyl0-r = 5.0 ; Split cylinder 0 radius (nm) cyl0-up = 0.75 ; Split cylinder 0 upper extension (nm) cyl0-down = 0.75 ; Split cylinder 0 lower extension (nm) cyl1-r = 5.0 ; same for other channel cyl1-up = 0.75 cyl1-down = 0.75 The cylinder options are used to define virtual geometric cylinders around the channel’s pore to track how many ions of which type have passed each channel. Ions will be counted as having traveled through a channel according to the definition of the channel’s cylinder radius, upper and lower extension, relative to the location of the respective split group. This will not affect the actual flux or exchange, but will provide you with the ion permeation numbers across each of the channels. Note that an ion can only be counted as passing through a particular channel if it is detected *within* the defined split cylinder in a swap step. If ``swap-frequency`` is chosen too high, a particular ion may be detected in compartment **A** in one swap step, and in compartment **B** in the following swap step, so it will be unclear through which of the channels it has passed. A double-layered system for CompEL simulations can be easily prepared by duplicating an existing membrane/channel MD system in the direction of the membrane normal (typically :math:`z`) with :ref:`gmx editconf` ``-translate 0 0 ``, where ``l_z`` is the box length in that direction. If you have already defined index groups for the channel for the single-layered system, :ref:`gmx make_ndx` ``-n index.ndx -twin`` will provide you with the groups for the double-layered system. To suppress large fluctuations of the membranes along the swap direction, it may be useful to apply a harmonic potential (acting only in the swap dimension) between each of the two channel and/or bilayer centers using umbrella pulling (see section :ref:`pull`). Multimeric channels ^^^^^^^^^^^^^^^^^^^ If a split group consists of more than one molecule, the correct PBC image of all molecules with respect to each other has to be chosen such that the channel center can be correctly determined. |Gromacs| assumes that the starting structure in the :ref:`tpr` file has the correct PBC representation. Set the following environment variable to check whether that is the case: - ``GMX_COMPELDUMP``: output the starting structure after it has been made whole to :ref:`pdb` file.