docs/reference-manual/special/comp-electrophys.rst

   1 .. _compel:
   2
   3 Computational Electrophysiology
   4 -------------------------------
   5
   6 The Computational Electrophysiology (CompEL) protocol
   7 :ref:`147 <refKutzner2011b>` allows the simulation of ion flux through membrane channels,
   8 driven by transmembrane potentials or ion concentration gradients. Just
   9 as in real cells, CompEL establishes transmembrane potentials by
  10 sustaining a small imbalance of charges :math:`\Delta q` across the
  11 membrane, which gives rise to a potential difference :math:`\Delta U`
  12 according to the membrane capacitance:
  13
  14 .. math:: \Delta U = \Delta q / C_{membrane}
  15           :label: eqnmembcap
  16
  17 The transmembrane electric field and concentration gradients are
  18 controlled by :ref:`mdp` options, which allow the user to set
  19 reference counts for the ions on either side of the membrane. If a
  20 difference between the actual and the reference numbers persists over a
  21 certain time span, specified by the user, a number of ion/water pairs
  22 are exchanged between the compartments until the reference numbers are
  23 restored. Alongside the calculation of channel conductance and ion
  24 selectivity, CompEL simulations also enable determination of the channel
  25 reversal potential, an important characteristic obtained in
  26 electrophysiology experiments.
  27
  28 In a CompEL setup, the simulation system is divided into two
  29 compartments **A** and **B** with independent ion concentrations. This
  30 is best achieved by using double bilayer systems with a copy (or copies)
  31 of the channel/pore of interest in each bilayer
  32 (:numref:`Fig. %s <fig-compelsetup>` A, B). If the channel axes
  33 point in the same direction, channel flux is observed simultaneously at
  34 positive and negative potentials in this way, which is for instance
  35 important for studying channel rectification.
  36
  37 .. _fig-compelsetup:
  38
  39 .. figure:: plots/compelsetup.*
  40    :width: 13.50000cm
  41
  42    Typical double-membrane setup for CompEL simulations (A, B).
  43    Ion/water molecule exchanges will be performed as needed between the
  44    two light blue volumes around the dotted black lines (A). Plot (C)
  45    shows the potential difference :math:`\Delta U` resulting from the
  46    selected charge imbalance :math:`\Delta q_{ref}` between the
  47    compartments.
  48
  49 The potential difference :math:`\Delta U` across the membrane is easily
  50 calculated with the :ref:`gmx potential <gmx potential>` utility. By this, the potential drop
  51 along :math:`z` or the pore axis is exactly known in each time interval
  52 of the simulation (:numref:`Fig. %s <fig-compelsetup>` C). Type and number of ions
  53 :math:`n_i` of charge :math:`q_i`, traversing the channel in the
  54 simulation, are written to the swapions.xvg output file, from which the
  55 average channel conductance :math:`G` in each interval :math:`\Delta t`
  56 is determined by:
  57
  58 .. math:: G = \frac{\sum_{i} n_{i}q_{i}}{\Delta t \, \Delta U} \, .
  59           :label: eqnchannelcond
  60
  61 The ion selectivity is calculated as the number flux ratio of different
  62 species. Best results are obtained by averaging these values over
  63 several overlapping time intervals.
  64
  65 The calculation of reversal potentials is best achieved using a small
  66 set of simulations in which a given transmembrane concentration gradient
  67 is complemented with small ion imbalances of varying magnitude. For
  68 example, if one compartment contains 1M salt and the other 0.1M, and
  69 given charge neutrality otherwise, a set of simulations with
  70 :math:`\Delta q = 0\,e`, :math:`\Delta q = 2\,e`,
  71 :math:`\Delta q = 4\,e` could be used. Fitting a straight line through
  72 the current-voltage relationship of all obtained :math:`I`-:math:`U`
  73 pairs near zero current will then yield :math:`U_{rev}`.
  74
  75 Usage
  76 ^^^^^
  77
  78 The following :ref:`mdp` options control the CompEL protocol:
  79
  80 ::
  81
  82     swapcoords     = Z        ; Swap positions: no, X, Y, Z
  83     swap-frequency = 100      ; Swap attempt frequency
  84
  85 Choose ``Z`` if your membrane is in the :math:`xy`-plane
  86 (:numref:`Fig. %s <fig-compelsetup>`). Ions will be exchanged
  87 between compartments depending on their :math:`z`-positions alone.
  88 ``swap-frequency`` determines how often a swap attempt will
  89 be made. This step requires that the positions of the split groups, the
  90 ions, and possibly the solvent molecules are communicated between the
  91 parallel processes, so if chosen too small it can decrease the
  92 simulation performance. The ``Position swapping`` entry in
  93 the cycle and time accounting table at the end of the
  94 ``md.log`` file summarizes the amount of runtime spent in
  95 the swap module.
  96
  97 ::
  98
  99     split-group0   = channel0 ; Defines compartment boundary
 100     split-group1   = channel1 ; Defines other compartment boundary
 101     massw-split0   = no       ; use mass-weighted center?
 102     massw-split1   = no
 103
 104 ``split-group0`` and ``split-group1`` are two
 105 index groups that define the boundaries between the two compartments,
 106 which are usually the centers of the channels. If
 107 ``massw-split0`` or ``massw-split1`` are set to
 108 ``yes``, the center of mass of each index group is used as
 109 boundary, here in :math:`z`-direction. Otherwise, the geometrical
 110 centers will be used (:math:`\times` in
 111 :numref:`Fig. %s <fig-compelsetup>` A). If, such as here, a membrane
 112 channel is selected as split group, the center of the channel will
 113 define the dividing plane between the compartments (dashed horizontal
 114 lines). All index groups must be defined in the index file.
 115
 116 If, to restore the requested ion counts, an ion from one compartment has
 117 to be exchanged with a water molecule from the other compartment, then
 118 those molecules are swapped which have the largest distance to the
 119 compartment-defining boundaries (dashed horizontal lines). Depending on
 120 the ion concentration, this effectively results in exchanges of
 121 molecules between the light blue volumes. If a channel is very
 122 asymmetric in :math:`z`-direction and would extend into one of the swap
 123 volumes, one can offset the swap exchange plane with the
 124 ``bulk-offset`` parameter. A value of 0.0 means no offset
 125 :math:`b`, values :math:`-1.0 < b < 0` move the swap exchange plane
 126 closer to the lower, values :math:`0 < b < 1.0` closer to the upper
 127 membrane. :numref:`Fig. %s <fig-compelsetup>` A (left) depicts that
 128 for the **A** compartment.
 129
 130 ::
 131
 132     solvent-group  = SOL      ; Group containing the solvent molecules
 133     iontypes       = 3        ; Number of different ion types to control
 134     iontype0-name  = NA       ; Group name of the ion type
 135     iontype0-in-A  = 51       ; Reference count of ions of type 0 in A
 136     iontype0-in-B  = 35       ; Reference count of ions of type 0 in B
 137     iontype1-name  = K
 138     iontype1-in-A  = 10
 139     iontype1-in-B  = 38
 140     iontype2-name  = CL
 141     iontype2-in-A  = -1
 142     iontype2-in-B  = -1
 143
 144 The group name of solvent molecules acting as exchange partners for the
 145 ions has to be set with ``solvent-group``. The number of
 146 different ionic species under control of the CompEL protocol is given by
 147 the ``iontypes`` parameter, while
 148 ``iontype0-name`` gives the name of the index group
 149 containing the atoms of this ionic species. The reference number of ions
 150 of this type can be set with the ``iontype0-in-A`` and
 151 ``iontype0-in-B`` options for compartments **A** and **B**,
 152 respectively. Obviously, the sum of ``iontype0-in-A`` and
 153 ``iontype0-in-B`` needs to equal the number of ions in the
 154 group defined by ``iontype0-name``. A reference number of
 155 ``-1`` means: use the number of ions as found at the
 156 beginning of the simulation as the reference value.
 157
 158 ::
 159
 160     coupl-steps    = 10       ; Average over these many swap steps
 161     threshold      = 1        ; Do not swap if < threshold
 162
 163 If ``coupl-steps`` is set to 1, then the momentary ion
 164 distribution determines whether ions are exchanged.
 165 ``coupl-steps > 1`` will use the time-average of ion
 166 distributions over the selected number of attempt steps instead. This
 167 can be useful, for example, when ions diffuse near compartment
 168 boundaries, which would lead to numerous unproductive ion exchanges. A
 169 ``threshold`` of 1 means that a swap is performed if the
 170 average ion count in a compartment differs by at least 1 from the
 171 requested values. Higher thresholds will lead to toleration of larger
 172 differences. Ions are exchanged until the requested number :math:`\pm`
 173 the threshold is reached.
 174
 175 ::
 176
 177     cyl0-r         = 5.0      ; Split cylinder 0 radius (nm)
 178     cyl0-up        = 0.75     ; Split cylinder 0 upper extension (nm)
 179     cyl0-down      = 0.75     ; Split cylinder 0 lower extension (nm)
 180     cyl1-r         = 5.0      ; same for other channel
 181     cyl1-up        = 0.75
 182     cyl1-down      = 0.75
 183
 184 The cylinder options are used to define virtual geometric cylinders
 185 around the channel’s pore to track how many ions of which type have
 186 passed each channel. Ions will be counted as having traveled through a
 187 channel according to the definition of the channel’s cylinder radius,
 188 upper and lower extension, relative to the location of the respective
 189 split group. This will not affect the actual flux or exchange, but will
 190 provide you with the ion permeation numbers across each of the channels.
 191 Note that an ion can only be counted as passing through a particular
 192 channel if it is detected *within* the defined split cylinder in a swap
 193 step. If ``swap-frequency`` is chosen too high, a particular
 194 ion may be detected in compartment **A** in one swap step, and in
 195 compartment **B** in the following swap step, so it will be unclear
 196 through which of the channels it has passed.
 197
 198 A double-layered system for CompEL simulations can be easily prepared by
 199 duplicating an existing membrane/channel MD system in the direction of
 200 the membrane normal (typically :math:`z`) with
 201 :ref:`gmx editconf` ``-translate 0 0 <l_z>``, where ``l_z`` is the box
 202 length in that direction. If you have already defined index groups for
 203 the channel for the single-layered system, :ref:`gmx make_ndx`
 204 ``-n index.ndx -twin`` will provide you with the groups for the
 205 double-layered system.
 206
 207 To suppress large fluctuations of the membranes along the swap
 208 direction, it may be useful to apply a harmonic potential (acting only
 209 in the swap dimension) between each of the two channel and/or bilayer
 210 centers using umbrella pulling (see section :ref:`pull`).
 211
 212 Multimeric channels
 213 ^^^^^^^^^^^^^^^^^^^
 214
 215 If a split group consists of more than one molecule, the correct PBC
 216 image of all molecules with respect to each other has to be chosen such
 217 that the channel center can be correctly determined. |Gromacs| assumes
 218 that the starting structure in the :ref:`tpr` file has the
 219 correct PBC representation. Set the following environment variable to
 220 check whether that is the case:
 221
 222 -  ``GMX_COMPELDUMP``: output the starting structure after
 223    it has been made whole to :ref:`pdb` file.