[gmx-users] Computational electrophysiology issues with leakage

Carlos Navarro carlos.navarro87 at gmail.com
Wed Apr 12 04:31:15 CEST 2017


Dear all,
I’m trying to set up the computational electrophysiology protocol. To do
this, I’m trying to replicate the following article
http://www.sciencedirect.com/science/article/pii/S0006349511007065
but instead of using g_membed I used a CG approximation (using MARTINI ff +
insane script to set up the system + back mapping to recover the resolution
of the system).
which at the end generate different numbers of atoms with respect to the
systems of the paper (but with the same salt concentration).
After the CG simulations and backpamming I run a simulation of 50ns to
equilibrate the system.
Followed this, I set up the double bilayer by using the makeSandwish script
(https://www.mpibpc.mpg.de/grubmueller/compel)
To test the protocol, I set up an 12e charge imbalance, but I keep seeing
‘leakage’ in the swapions.xvg file practically at the beginning of the
simulation

*# This file was created Tue Apr 11 16:05:21 2017*
*# Created by:*
*# GROMACS:      gmx mdrun, VERSION 5.0.1*
*# Executable:   /usr/local/gromacs5.0.1/bin/gmx*
*# Library dir:  /usr/local/gromacs5.0.1/share/gromacs/top*
*# Command line:*
*#   mdrun -v*
*# mdrun is part of G R O M A C S:*
*#*
*# Good gRace! Old Maple Actually Chews Slate*
*#*
*@    title "Ion counts"*
*@    xaxis  label "Time (ps)"*
*@    yaxis  label "counts"*
*@TYPE xy*
*# ion group contains 4628 atoms with 1 atom in each molecule.*
*# split0 group contains 3069 atoms.*
*# split1 group contains 3069 atoms.*
*# solvent group contains 211728 atoms with 3 atoms in each molecule.*
*#*
*# Initial positions of split groups:*
*# split0 group Z-center 4.945151 nm*
*# split1 group Z-center 14.793151 nm*
*#*
*# split0 cylinder radius 6.000000 nm, up 1.000000 nm, down 1.000000 nm*
*# split1 cylinder radius 6.000000 nm, up 1.000000 nm, down 1.000000 nm*
*#*
*# Coupling constant (number of swap attempt steps to average over): 10
 (translates to 0.800000 ps).*
*# Threshold is 1.000000*
*#*
*# Remarks about which atoms passed which channel use global atoms numbers
starting at one.*
*# Requested charge imbalance is Q(A)-Q(B) = 12z.*
*@ view 0.15, 0.15, 0.75, 0.85*
*@ legend on*
*@ legend box on*
*@ legend loctype view*
*@ legend 0.78, 0.8*
*@ legend length 2*
*@ s0 legend "A anions"*
*@ s1 legend "A av. mismatch to 1166-"*
*@ s2 legend "A netto anion influx"*
*@ s3 legend "A cations"*
*@ s4 legend "A av. mismatch to 1142+"*
*@ s5 legend "A netto cation influx"*
*@ s6 legend "B anions"*
*@ s7 legend "B av. mismatch to 1178-"*
*@ s8 legend "B netto anion influx"*
*@ s9 legend "B cations"*
*@ s10 legend "B av. mismatch to 1142+"*
*@ s11 legend "B netto cation influx"*
*@ s12 legend "Z-center of mass of split group 0"*
*@ s13 legend "Z-center of mass of split group 1"*
*@ s14 legend "A->ch0->B anion permeations"*
*@ s15 legend "A->ch0->B cation permeations"*
*@ s16 legend "A->ch1->B anion permeations"*
*@ s17 legend "A->ch1->B cation permeations"*
*@ s18 legend "leakage"*
*# Instantaneous ion counts and time-averaged differences to requested
numbers*
*#   time[ps]   A_an   diff   t_in  A_cat   diff   t_in   B_an   diff
t_in  B_cat   diff   t_in    Z-Split0    Z-Split1  A-ch0-B_an A-ch0-B_cat
 A-ch1-B_an A-ch1-B_cat ion_leakage*
* 8.00000e-02   1171    5.9      0   1143    0.1      0   1173   -5.9
 0   1141   -0.1      0   4.953e+00   1.482e+01           0           0
      0           0           0*
*# Solv. molecules in comp.A: 35291   comp.B: 35285*
* 8.00000e-02   1166    0.9     -5   1143    0.1      0   1178   -0.9
 5   1141   -0.1      0   4.953e+00   1.482e+01           0           0
      0           0           0  # after swap*
* # Warning: step 80, ion 184990 (-) moved from Domain_B to Domain_A
(probably through the membrane)*
* # Warning: step 80, ion 185463 (-) moved from Domain_A to Domain_B
(probably through the membrane)*
* 1.60000e-01   1167    0.9     -5   1143    0.2      0   1177   -0.9
 5   1141   -0.2      0   4.960e+00   1.484e+01           0           0
      0           0           2*
* # Warning: step 120, ion 183757 (+) moved from Domain_B to Domain_A
(probably through the membrane)*
* # Warning: step 120, ion 185463 (-) moved from Domain_B to Domain_A
(probably through the membrane)*
* # Warning: step 120, ion 369231 (+) moved from Domain_A to Domain_B
(probably through the membrane)*
* 2.40000e-01   1169    1.1     -5   1142    0.2      0   1175   -1.1
 5   1142   -0.2      0   4.966e+00   1.486e+01           0           0
      0           0           5*

I have even tried with a charge imbalance of 2e but with the same problem,
so I suspect the issue is related to the cyl0-r and cyl1-r variables.
Here, and even using values of 7 nm of radius I kept seeing leakage, so I
don’t know what can be happening.

Here are the CompEL settings:

; Ion/water position swapping for computational electrophysiology setups

; Swap positions along direction: no, X, Y, Z

swapcoords = Z

; Swap attempt frequency

swap-frequency           = 40

; Two index groups that contain the compartment-partitioning atoms

split-group0             = Backbone_&_TM-helices

split-group1             = Backbone_copy_&_TM-helices_copy

; Use center of mass of split groups (yes/no), otherwise center of geometry
is used

massw-split0             = no

massw-split1             = no

; Group name of ions that can be exchanged with solvent molecules

swap-group               = Ion_Ion_copy

; Group name of solvent molecules

solvent-group            = SOL_SOL_copy

; Split cylinder: radius, upper and lower extension (nm) (this will define
the channels)

; Note that the split cylinder settings do not have an influence on the
swapping protocol,

; however, if correctly defined, the ion permeation events are counted per
channel

cyl0-r                   = 7.0

cyl0-up                  = 2.5

cyl0-down                = 2.5

cyl1-r                   = 7.0

cyl1-up                  = 2.5

cyl1-down                = 2.5

; Average the number of ions per compartment over these many swap attempt
steps

coupl-steps              = 10

; Requested number of anions and cations for each of the two compartments

; -1 means fix the numbers as found in time step 0

anionsA                  = 1166

cationsA                 = -1

anionsB                  = 1178

cationsB                 = -1

; Start to swap ions if threshold difference to requested count is reached

threshold                = 1


And here is the .ndx file:


  0 System              : 185556 atoms

  1 Protein             : 15336 atoms

  2 Protein-H           :  7779 atoms

  3 C-alpha             :  1023 atoms

  4 Backbone            :  3069 atoms

  5 MainChain           :  4089 atoms

  6 MainChain+Cb        :  4986 atoms

  7 MainChain+H         :  5106 atoms

  8 SideChain           : 10230 atoms

  9 SideChain-H         :  3690 atoms

 10 Prot-Masses         : 15336 atoms

 11 non-Protein         : 170220 atoms

 12 Other               : 62042 atoms

 13 POPC                : 62042 atoms

 14 NA                  :  1142 atoms

 15 CL                  :  1172 atoms

 16 Water               : 105864 atoms

 17 SOL                 : 105864 atoms

 18 non-Water           : 79692 atoms

 19 Ion                 :  2314 atoms

 20 POPC                : 62042 atoms

 21 NA                  :  1142 atoms

 22 CL                  :  1172 atoms

 23 Water_and_ions      : 108178 atoms

 24 TM-helices          : 15336 atoms

 25 System_copy         : 185556 atoms

 26 Protein_copy        : 15336 atoms

 27 Protein-H_copy      :  7779 atoms

 28 C-alpha_copy        :  1023 atoms

 29 Backbone_copy       :  3069 atoms

 30 MainChain_copy      :  4089 atoms

 31 MainChain+Cb_copy   :  4986 atoms

 32 MainChain+H_copy    :  5106 atoms

 33 SideChain_copy      : 10230 atoms

 34 SideChain-H_copy    :  3690 atoms

 35 Prot-Masses_copy    : 15336 atoms

 36 non-Protein_copy    : 170220 atoms

 37 Other_copy          : 62042 atoms

 38 POPC_copy           : 62042 atoms

 39 NA_copy             :  1142 atoms

 40 CL_copy             :  1172 atoms

 41 Water_copy          : 105864 atoms

 42 SOL_copy            : 105864 atoms

 43 non-Water_copy      : 79692 atoms

 44 Ion_copy            :  2314 atoms

 45 POPC_copy           : 62042 atoms

 46 NA_copy             :  1142 atoms

 47 CL_copy             :  1172 atoms

 48 Water_and_ions_copy : 108178 atoms

 49 TM-helices_copy     : 15336 atoms

 50 SOL_SOL_copy        : 211728 atoms

 51 Ion_Ion_copy        :  4628 atoms

 52 TM-helices_TM-helices_copy: 30672 atoms

 53 Protein_Protein_copy: 30672 atoms

 54 non-Protein_non-Protein_copy: 340440 atoms

 55 System_System_copy  : 371112 atoms

 56 Backbone_&_TM-helices:  3069 atoms

 57 Backbone_copy_&_TM-helices_copy:  3069 atoms


Here are the .gro and .mdp files:

.gro:   https://www.dropbox.com/s/e4j74zff7mn93bg/membrane.double.gro?dl=0

.mdp: https://www.dropbox.com/s/rsoiv5nkao2ut96/membrane.double.mdp?dl=0

If you need more info just let me know.

Hope someone can help me to find out where the error is

Best regards,

Carlos


-- 
Carlos Navarro Retamal
Bioinformatic Engineering. PhD
Postdoctoral Researcher in Center for Bioinformatics and Molecular
Simulations
Universidad de Talca
Av. Lircay S/N, Talca, Chile
T: (+56) 712201 <//T: (+56) 712201> 798
E: carlos.navarro87 at gmail.com o cnavarro at utalca.cl


More information about the gromacs.org_gmx-users mailing list