[gmx-users] Net charge in implicit solvent simulations

[saber naderi]

Dear gmx-users,

We currently performed simulations to study the effect of repulsive
electrostatic interactions on the stability of a silk-like protein. The
amino acid sequence of our protein is the following: Lys [ (Gly-Ala)_3 Gly
(Gln or Lys)]_3 where (Gln or Lys) means that we alternatively have Lys or
Gln in that position starting with Gln. We start our simulations with three
of these proteins on top of each other (each having a beta-sheet
structure).

We performed two set of Replica exchange simulations with implicit solvent
(see below for an mdp file). In one the lysine residues are positively
charged and in the other they are neutral (LYN). Surprisingly, we see that
the charged protein is more stable than the neutral one, i.e., repulsive
electrostatic interactions are helping the charged protein not to unfold!!!

I suppose this is not a normal behavior. What could go wrong in such
simulations that results in this?

Thanks for your time,

Regards,
Saber


mdp file for one of the replicas:

;  RUN CONTROL PARAMETERS
integrator               = sd
tinit                    = 0       ; Starting time
dt                       = 0.002   ; 2 femtosecond time step for integration
nsteps                   = 10000000

; OUTPUT CONTROL OPTIONS
nstxout                  = 50000 ; Writing full precision coordinates every
nanosecond
nstvout                  = 50000 ; Writing velocities every nanosecond
nstfout                  = 0     ; Not writing forces
nstlog                   = 500  ; Writing to the log file every step
nstenergy                = 500  ; Writing out energy information every step
nstxtcout                = 500  ; Writing coordinates every 5 ps
energygrps               = System

; NEIGHBORSEARCHING PARAMETERS
nstlist                  = 0
ns-type                  = simple
pbc                      = no
rlist                    = 0

; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype = cut-off
rvdw = 0
rcoulomb = 0


; Temperature coupling
tc-grps                  = System
ld_seed             = -1
tau_t                    = 0.1
ref_t                    = 298.00,
; implicit solvent

implicit_solvent    = GBSA
gb_algorithm    = OBC
rgbradii                = 0
gb_epsilon_solvent = 80

; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel                  = yes
gen_temp                = 300
comm_mode                = Angular

; OPTIONS FOR BONDS
constraints              = all-bonds
constraint-algorithm     = Lincs
unconstrained-start      = yes
lincs-order              = 4
lincs-iter               = 1
lincs-warnangle          = 30