Print

Simulate Water - Methanol Mixture (Amber-GAFF Gromacs 5.0)

Written by Super User. Posted in Uncategorised

recently, I got confronted with the task to simulate partial protein unfolding under denaturating conditions. Experimentally, it is known from Ubiquitin (for example pdb Code 1UBQ) that a 50% vol MeOH solution under acidic conditions gives rise to partial unfolding - the so-called A-state. The structure has been proposed e.g. by Bowers et al. from simulation data. A dertailed experimental structural characterization is still not available for this flexible system. The first task to obtain a computational model for a protein in such solution, however, is already not quite standard and involves the setup of a non-water solvation for the protein.

This little tutorial explains the setup for the Water - Methanol solvation.

For water I use the TIP3P parametrization. For MeOH I use GAFF. A 50 % Water Methanol (vol %) mixture corresponds to about 1:2.25 molecules of MeOH to H20.

Simple pdb files of a water molecule and a MeOH molecule are needed. Then I construct a box of 45*45*45 Angstr with the mixture using the GREAT software Packmol with the following input file.

# A mixture of water and meoh
tolerance 2.0
filetype pdb
output mixture.pdb
structure water.pdb
number 2250 
inside box 0. 0. 0. 45. 45. 45. 
end structure
structure meoh.pdb
number 1000
inside box 0. 0. 0. 45. 45. 45. 
end structure

The result is the following box

BOXwithsolvents

I created the topology by hand - using the GAFF parametrization with Acpype and adding simple modifications as shown below to include the TIP3P water model.

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes             0.5     0.8333

[ atomtypes ]
;name   bond_type     mass     charge   ptype   sigma         epsilon       Amb
 h1       h1          0.00000  0.00000   A     2.47135e-01   6.56888e-02 ; 1.39  0.0157
 c3       c3          0.00000  0.00000   A     3.39967e-01   4.57730e-01 ; 1.91  0.1094
 oh       oh          0.00000  0.00000   A     3.06647e-01   8.80314e-01 ; 1.72  0.2104
 ho       ho          0.00000  0.00000   A     0.00000e+00   0.00000e+00 ; 0.00  0.0000
 HW       HW          1.008    0.0000    A     0.00000e+00   0.00000e+00
 OW       OW   8      16.00    0.0000    A     3.15061e-01   6.36386e-01

[ moleculetype ]
;name            nrexcl
 MOH             3

[ atoms ]
;   nr  type  resi  res  atom  cgnr     charge      mass       ; qtot   bond_type
     1   h1     1   MOH   HC1    1     0.028700      1.00800 ; qtot 0.029
     2   c3     1   MOH    C1    2     0.116700     12.01000 ; qtot 0.145
     3   h1     1   MOH   HC2    3     0.028700      1.00800 ; qtot 0.174
     4   h1     1   MOH   HC3    4     0.028700      1.00800 ; qtot 0.203
     5   oh     1   MOH    O1    5    -0.598801     16.00000 ; qtot -0.396
     6   ho     1   MOH   HO1    6     0.396000      1.00800 ; qtot 0.000

[ bonds ]
;   ai     aj funct   r             k
     1      2   1    1.0930e-01    2.8108e+05 ;    HC1 - C1    
     2      3   1    1.0930e-01    2.8108e+05 ;     C1 - HC2   
     2      4   1    1.0930e-01    2.8108e+05 ;     C1 - HC3   
     2      5   1    1.4260e-01    2.6284e+05 ;     C1 - O1    
     5      6   1    9.7400e-02    3.0928e+05 ;     O1 - HO1   

[ pairs ]
;   ai     aj    funct
     1      6      1 ;    HC1 - HO1   
     3      6      1 ;    HC2 - HO1   
     4      6      1 ;    HC3 - HO1   

[ angles ]
;   ai     aj     ak    funct   theta         cth
     1      2      3      1    1.0955e+02    3.2786e+02 ;    HC1 - C1     - HC2   
     1      2      4      1    1.0955e+02    3.2786e+02 ;    HC1 - C1     - HC3   
     1      2      5      1    1.0988e+02    4.2652e+02 ;    HC1 - C1     - O1    
     2      5      6      1    1.0816e+02    3.9405e+02 ;     C1 - O1     - HO1   
     3      2      4      1    1.0955e+02    3.2786e+02 ;    HC2 - C1     - HC3   
     3      2      5      1    1.0988e+02    4.2652e+02 ;    HC2 - C1     - O1    
     4      2      5      1    1.0988e+02    4.2652e+02 ;    HC3 - C1     - O1    

[ dihedrals ] ; propers
; treated as RBs in GROMACS to use combine multiple AMBER torsions per quartet
;    i      j      k      l   func    C0         C1         C2         C3         C4         C5
     1      2      5      6      3    0.69733    2.09200    0.00000   -2.78933    0.00000    0.00000 ;    HC1-    C1-    O1-   HO1
     3      2      5      6      3    0.69733    2.09200    0.00000   -2.78933    0.00000    0.00000 ;    HC2-    C1-    O1-   HO1
     4      2      5      6      3    0.69733    2.09200    0.00000   -2.78933    0.00000    0.00000 ;    HC3-    C1-    O1-   HO1

[ system ]

methanol-water mixture

#include "./amber99.ff/tip3p.itp"

[ molecules ]
MOH    1000
SOL    2250

with the correct number of SOL and MOH molecules in the box, relaxation and careful equilibration is performed. First an NVT equilibration (monitor the pressure and adjust the box dimensions in order to obtain a reasonable order of magnitude), then an NPT simulation. Check the density of the solvent mixture! In my case I obtain 0.943 kg/m^3, which is quite reasonable (compare Mikhail et al. ;J. Chem. Engin. Data 1961). Now the interesting part: The equilibrated box of solvent can be used to solvate a system usig the GROMCS gmx solvate function with the -cs switch selecting the generated box of MeOH/H2O as template.

Notice, however, that the used forcefield may not be the best available for representing condensed phase properties of solvent mixtures.

SOLVATEDUBQ

Finally, a fully protonated Ubiquitin molecule in H2O/MeOH with Cl ions as countercharges.