[gmx-users] Strange bilayer behavior in protein-multimembrane models

[Justin A. Lemkul]

Thomas Schmidt wrote:
> Dear all,
> 
> by doing MD simulations of a protein complex embedded in 2 membranes
> (inner and outer membrane of bacteria, POPE), we observe a bilayer
> splitting in one of the membranes. This has the effect that the bilayer
> forms "bubbles" with vacuum inside.
> 

We have seen this before with POPE membranes, which appear to be particularly 
problematic.  We solved it by using a better equilibration protocol that 
involved simulated annealing (from a very low temperature up to 310 K).  For us, 
it seemed to happen independently of just about anything we tried.

-Justin

> It might have something to do with the PME handling of GroMACS 4.0.3
> (GROMOS96-53a6 ff).
> - there is no splitting using a 1 nm smaller boxsize in exactly the same
>     model
> - if we switch off the PME mode and use only Cut-off's we don't have
>     any splitting effect
> - the behavior of the bilayer splitting depends on the number of used
>     cluster/cpu nodes
> 
> Changing the "fourierspacing" parameter to create different PME grids
> has no effect to avoid bilayer splitting.
> 
> Changing: thermostat | barostat (semiisotropic):
> - berendsen | berendsen:         "splitting"
> - v-rescale | parrinello-rahman: "splitting"/"keep together" (50:50)
> 
> 
> Here's the mdp file of our last try:
> title           = Yep, sometimes I will cause a bilayer splitting;
> 
> ; The following lines tell the program the standard locations where to
> find certain files
> cpp             = cpp; Preprocessor
> include         = -I../top; Directories to include in the topology
> format
> define          = -DPOSRES; Apply position restraints.
> 
> 
> 
> ; RUN CONTROL PARAMETERS
> integrator               = MD
> ; Start time and timestep in ps
> tinit                    = 0
> dt                       = 0.002
> nsteps                   = 10000000
> ; For exact run continuation or redoing part of a run
> init_step                = 0
> ; mode for center of mass motion removal
> comm-mode                = Linear
> ; number of steps for center of mass motion removal
> nstcomm                  = 1
> ; group(s) for center of mass motion removal
> comm-grps                = Protein POPE_center POPE_inner POPE_outer
> SOL_NA+_CL-
> 
> 
> 
> ; ENERGY MINIMIZATION OPTIONS
> ; Force tolerance and initial step-size
> emtol                    = 500
> emstep                   = 0.01
> ; Max number of iterations in relax_shells
> niter                    = 0
> ; Step size (1/ps^2) for minimization of flexible constraints
> fcstep                   = 0
> 
> 
> 
> ; OUTPUT CONTROL OPTIONS
> ; Output frequency for coords (x), velocities (v) and forces (f)
> nstxout                  = 500000
> nstvout                  = 500000
> nstfout                  = 0
> ; Checkpointing helps you continue after crashes
> nstcheckpoint            = 50000
> ; Output frequency for energies to log file and energy file
> nstlog                   = 5000
> nstenergy                = 5000
> ; Output frequency and precision for xtc file
> nstxtcout                = 50000
> xtc-precision            = 1000
> ; This selects the subset of atoms for the xtc file. You can
> ; select multiple groups. By default all atoms will be written.
> xtc-grps                 = Protein POPE_center POPE_inner POPE_outer
> SOL_NA+_CL-
> ; Selection of energy groups
> energygrps               = Protein POPE_center POPE_inner POPE_outer
> SOL_NA+_CL-
> 
> 
> 
> ; NEIGHBORSEARCHING PARAMETERS
> ; nblist update frequency
> nstlist                  = 10
> ; ns algorithm (simple or grid)
> ns_type                  = grid
> ; Periodic boundary conditions: xyz (default), no (vacuum)
> ; or full (infinite systems only)
> pbc                      = xyz
> ; nblist cut-off
> rlist                    = 1.15
> domain-decomposition     = no
> 
> 
> 
> ; OPTIONS FOR ELECTROSTATICS AND VDW
> ; Method for doing electrostatics
> coulombtype              = PME
> rcoulomb-switch          = 0
> rcoulomb                 = 1.15
> ; Dielectric constant (DC) for cut-off or DC of reaction field
> epsilon-r                = 1
> ; Method for doing Van der Waals
> vdw-type                 = Cut-off
> ; cut-off lengths
> rvdw-switch              = 0
> rvdw                     = 1.4
> ; Apply long range dispersion corrections for Energy and Pressure
> DispCorr                 = No
> ; Extension of the potential lookup tables beyond the cut-off
> table-extension          = 1
> ; Spacing for the PME/PPPM FFT grid
> fourierspacing           = 0.12
> ; FFT grid size, when a value is 0 fourierspacing will be used
> fourier_nx               = 0
> fourier_ny               = 0
> fourier_nz               = 0
> ; EWALD/PME/PPPM parameters
> pme_order                = 4
> ewald_rtol               = 1e-05
> ewald_geometry           = 3d
> epsilon_surface          = 0
> optimize_fft             = yes
> 
> 
> 
> ; OPTIONS FOR WEAK COUPLING ALGORITHMS
> ; Temperature coupling
> ; Tcoupl                   = berendsen
> Tcoupl                   = V-rescale
> ; Groups to couple separately
> tc-grps                  = Protein POPE_center POPE_inner POPE_outer
> SOL_NA+_CL-
> ; Time constant (ps) and reference temperature (K)
> tau-t                    = 0.1 0.1 0.1 0.1 0.1
> ref-t                    = 310 310 310 310 310
> ; Pressure coupling
> ; Pcoupl                   = berendsen
> Pcoupl                   = Parrinello-Rahman
> Pcoupltype               = semiisotropic
> ; Time constant (ps), compressibility (1/bar) and reference P (bar)
> tau-p                    = 4 4
> compressibility          = 4.5e-5 4.5e-5
> ref-p                    = 1.0 1.0
> ; Random seed for Andersen thermostat
> andersen_seed            = 815131
> 
> 
> 
> ; GENERATE VELOCITIES FOR STARTUP RUN
> gen_vel                  = yes
> gen-temp                 = 310
> gen-seed                 = 24071998
> 
> 
> 
> ; OPTIONS FOR BONDS
> constraints              = all-bonds
> ; Type of constraint algorithm
> constraint-algorithm     = Lincs
> ; Do not constrain the start configuration
> unconstrained-start      = yes
> ; Use successive overrelaxation to reduce the number of shake iterations
> Shake-SOR                = no
> ; Relative tolerance of shake
> shake-tol                = 1e-04
> ; Highest order in the expansion of the constraint coupling matrix
> lincs-order              = 4
> ; Number of iterations in the final step of LINCS. 1 is fine for
> ; normal simulations, but use 2 to conserve energy in NVE runs.
> ; For energy minimization with constraints it should be 4 to 8.
> lincs-iter               = 4
> ; Lincs will write a warning to the stderr if in one step a bond
> ; rotates over more degrees than
> lincs-warnangle          = 90
> ; Convert harmonic bonds to morse potentials
> morse                    = no
> 
> 
> Does anyone have any idea?
> Many thanks for your time and any help.
> 
> Cheers,
> Thomas
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================