[gmx-users] Computational electrophysiology issues with membrane leakage
Kutzner, Carsten
ckutzne at gwdg.de
Wed Apr 12 07:55:25 CEST 2017
Hi,
the output file gives you the numbers of the atoms that seem to move from
compartment A to B without passing a channel. I would look at their
trajectories and inspect which path they actually take.
The channel trimer also needs to be whole at the beginning of the
simulation, otherwise the channel counting won't work. You can
use the GMX_COMPELDUMP environment variable to check this.
http://manual.gromacs.org/documentation/5.1.1/user-guide/environment-variables.html
Best,
Carsten
> On 12. Apr 2017, at 05:06, Carlos Navarro <carlos.navarro87 at gmail.com> wrote:
>
> 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
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
More information about the gromacs.org_gmx-users
mailing list