[gmx-users] Computational electrophysiology issues with membrane leakage
Carlos Navarro
carlos.navarro87 at gmail.com
Wed Apr 12 05:06:50 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
<http://airmail.calendar/2017-04-11%2016:05:21%20GMT-3> 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