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