[gmx-users] Computational electrophysiology (compEL) setup issues

[Francesco Petrizzelli]

Dear gmx users,

I'm trying to simulate change in the ion flux upon mutations using Computational electrophysiology (compEL).
I have used packmol-memgen in order to generate the double bilayer using a salt concentration of 1 M (3216 Cl- and 3168 K+).
Then, I have setted up the simulation following the compEL protocol. However, after 50 ns of simulation, I still have two big issues:

1) In the swapions logfile I've got a huge number of warnings regarding ions that moved from Domain_B to Domain_A probably throught the membrane (checking the movie with VMD it seems they move across the boundaries). I thought to tackle this issue by increasing the water layer.

2) In the same way, watching  the movie I have notice that ions enter the channel (the pore radius is about 3 A) but they do not cross from one domain to the other. I think I may have misunderstood the protocol, because in numerous article they talk about  setting an ionic imbalance between compartments. However, I have not performed it in any step. Do I have to manually set an ionic imbalance (and, in this case there are a tool to perform it?) or it should be set by the input compEL file? I have used the following inputs:

; Ion/water position swapping for computational electrophysiology setups
; Swap positions along direction: no, X, Y, Z
swapcoords = Z
adress                   = no
; Swap attempt frequency. 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 opera
tion. Do not try to swap every step, this will slow down the simulation a lot.
swap-frequency = 100
; Two index groups that contain the compartment-partitioning atoms
split-group0 = channel1
split-group1 = channel2
; Use center of mass of split groups (yes/no), otherwise geometrical center is used. Choose two index groups that divide the MD system into two compartments, here in Z-direction. If massw-split is activat
ed, 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.
massw-split0 = no
massw-split1 = no
; Group name of solvent molecules
solvent-group = Water
; Average the number of ions per compartment over these many swap attempt steps. 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 outswap
s for coupl-steps=1.
coupl-steps = 10
; Number of ion types to be controlled
iontypes = 2
; Names of the ion types that can be exchanged with solvent molecules
; -1 means fix the numbers as found in time step 0. These numbers have to add up to the total number of ions present in the swap group.
iontype0-name = Cl-_Cl-
iontype0-in-A = -1   ; requested number of Cl ions in compartment A
iontype0-in-B = -1   ; requested number of Cl ions in compartment B
iontype1-name = K+_K+
iontype1-in-A = -1   ; requested number of K+ ions in compartment A
iontype1-in-B = -1   ; requested number of K+ ions in compartment B
; Offset compartment-defining layers
bulk-offsetA = 0.0
bulk-offsetB = 0.0
; Split cylinder: radius, upper and lower extension (nm) (for counting on-the-fly diagnostics, has no influence on whether or not ions are swapped)
cyl0-r    = 0.7
cyl0-up   = 3
cyl0-down = 3
cyl1-r    = 0.7
cyl1-up   = 3
cyl1-down = 3
; Start to swap ions if threshold difference to requested count is reached. 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 request
ed values. Higher thresholds mean that larger differences are accepted. Ions are also only swapped until the requested number +/- the threshold is reached.

threshold = 1


; User defined thingies
user1-grps               =
user2-grps               =
userint1                 = 0
userint2                 = 0
userint3                 = 0
userint4                 = 0
userreal1                = 0
userreal2                = 0
userreal3                = 0
userreal4                = 0
; Electric fields
; Format for electric-field-x, etc. is: four real variables:
; amplitude (V/nm), frequency omega (1/ps), time for the pulse peak (ps),
; and sigma (ps) width of the pulse. Omega = 0 means static field,
; sigma = 0 means no pulse, leaving the field to be a cosine function.
electric-field-x         = 0 0 0 0
electric-field-y         = 0 0 0 0
electric-field-z         = 0 0 0 0

Thank you very much for your help and I'm sorry if some questions sound silly, but I am a novel gmx user.

Regards.