Молекулярная динамика

Моделирование самосборки липидного бислоя

Запуск молекулярной динамики

Скачать необходимые файлы

In [ ]:
> wget kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro # координаты одного липида  
> wget kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp # параметры для липидов
> wget kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp # дополнительная топология для липида DPPC
> wget kodomo.cmm.msu.ru/~golovin/bilayer/b.top # заготовка топологии системы
> wget kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp # параметры для минимизации энергии
> wget kodomo.cmm.msu.ru/~golovin/bilayer/ pr.mdp # параметры для "утряски" воды
> wget kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp # параметры для молекулярной динамики

Создать ячейку с 64 молекулами липида

In [ ]:
> genconf -f dppc.gro -o b_64.gro -nbox 4 4 4

Показать одинокий липид и ячейку в pymol

In [ ]:
> editconf -f dppc.gro -o dppc.pdb 
> editconf -f b_64.gro -o b_64.pdb
In [1]:
import __main__
__main__.pymol_argv = [ 'pymol', '-x' ]

import pymol
pymol.finish_launching()
from pymol import cmd
from IPython.display import Image
import time
In [4]:
cmd.load('dppc.pdb')
time.sleep(10)
cmd.png('prak11_1.png')
time.sleep(10)
Image(filename='prak11_1.png')
Out[4]:
In [5]:
cmd.reinitialize()
cmd.load('b_64.pdb')
time.sleep(10)
cmd.png('prak11_2.png')
time.sleep(10)
Image(filename='prak11_2.png')
Out[5]:

Установим правильное количество липидов в файле b.top

In [ ]:
> cat b.top 
; topology for 1 alm molecule, 128 popc lipids, water and 1 counter ion
; alm.itp can be made in a straightforward manner with pdb2gmx, starting
; with a pdb file of alamethicin.
; make sure lipid.itp, popc.itp and alm.itp are in a location where
; grompp can find them (GMXLIB, current directory, or directory given in
; the .mdp file with the include option.
#include "ffgmx.itp"
#include "lipid.itp"
#include "dppc.itp"

#ifdef FLEX_SPC
#include "flexspc.itp"
#else
#include "spc.itp"
#endif

[ system ]
; name
bilayer in water

[ molecules ]
; name  number
DPPC    64

Добавим свободное пространство вокруг ячейки и оптимизируем геометрию

In [ ]:
> editconf -f b_64.gro -o b_ec -d 0.5 
> grompp -f em -c b_ec -p b -o b_em -maxwarn 2
> mdrun -deffnm b_em -v
                        
Getting Loaded...
Reading file b_em.tpr, VERSION 4.5.5 (single precision)
Starting 8 threads
Loaded with Money

Making 3D domain decomposition 2 x 2 x 2

Steepest Descents:
   Tolerance (Fmax)   =  1.00000e+00
   Number of steps    =      1000000
Step=    0, Dmax= 2.0e-02 nm, Epot=  4.74007e+05 Fmax= 4.37970e+05, atom= 1842
Step=    1, Dmax= 2.0e-02 nm, Epot=  7.19293e+04 Fmax= 3.57933e+04, atom= 2588
Step=    2, Dmax= 2.4e-02 nm, Epot=  3.56795e+04 Fmax= 9.52568e+03, atom= 2993
Step=    3, Dmax= 2.9e-02 nm, Epot=  3.46459e+04 Fmax= 1.50508e+04, atom= 2992
Step=    4, Dmax= 3.5e-02 nm, Epot=  2.99879e+04 Fmax= 1.37059e+04, atom= 2592
Step=    6, Dmax= 2.1e-02 nm, Epot=  1.30782e+04 Fmax= 4.02791e+03, atom= 1342
Step=    9, Dmax= 6.2e-03 nm, Epot=  1.04672e+04 Fmax= 1.87751e+03, atom= 2392
Step=   10, Dmax= 7.5e-03 nm, Epot=  1.02199e+04 Fmax= 3.76147e+03, atom= 1792
Step=   11, Dmax= 9.0e-03 nm, Epot=  1.01147e+04 Fmax= 4.34365e+03, atom= 2392
Step=   13, Dmax= 5.4e-03 nm, Epot=  7.19423e+03 Fmax= 1.91912e+03, atom= 2765
Step=   15, Dmax= 3.2e-03 nm, Epot=  6.81202e+03 Fmax= 2.03533e+03, atom= 2765
Step=   16, Dmax= 3.9e-03 nm, Epot=  6.37647e+03 Fmax= 3.00818e+03, atom= 2765
Step=   17, Dmax= 4.6e-03 nm, Epot=  6.07244e+03 Fmax= 2.75324e+03, atom= 2765
Step=   19, Dmax= 2.8e-03 nm, Epot=  5.60144e+03 Fmax= 8.43337e+02, atom= 2765
Step=   20, Dmax= 3.3e-03 nm, Epot=  5.49136e+03 Fmax= 3.08327e+03, atom= 2765
Step=   21, Dmax= 4.0e-03 nm, Epot=  5.09456e+03 Fmax= 2.09083e+03, atom= 2765
Step=   23, Dmax= 2.4e-03 nm, Epot=  4.83141e+03 Fmax= 8.78173e+02, atom= 2765
Step=   25, Dmax= 1.4e-03 nm, Epot=  4.78063e+03 Fmax= 1.07042e+03, atom= 2765
Step=   26, Dmax= 1.7e-03 nm, Epot=  4.51750e+03 Fmax= 1.00741e+03, atom= 2765
Step=   27, Dmax= 2.1e-03 nm, Epot=  4.46707e+03 Fmax= 1.77657e+03, atom= 2765
Step=   28, Dmax= 2.5e-03 nm, Epot=  4.35788e+03 Fmax= 1.28547e+03, atom= 915
Step=   30, Dmax= 1.5e-03 nm, Epot=  4.20725e+03 Fmax= 6.94208e+02, atom= 915
Step=   31, Dmax= 1.8e-03 nm, Epot=  3.97280e+03 Fmax= 1.39504e+03, atom= 915
Step=   32, Dmax= 2.2e-03 nm, Epot=  3.97165e+03 Fmax= 1.44422e+03, atom= 915
Step=   33, Dmax= 2.6e-03 nm, Epot=  3.85030e+03 Fmax= 1.70904e+03, atom= 915
Step=   35, Dmax= 1.6e-03 nm, Epot=  3.71519e+03 Fmax= 3.29197e+02, atom= 472
Step=   36, Dmax= 1.9e-03 nm, Epot=  3.32876e+03 Fmax= 1.57250e+03, atom= 915
Step=   39, Dmax= 5.6e-04 nm, Epot=  3.23180e+03 Fmax= 8.45033e+02, atom= 2765
Step=   42, Dmax= 1.7e-04 nm, Epot=  3.20617e+03 Fmax= 6.21543e+02, atom= 9153
Step=   50, Dmax= 1.6e-06 nm, Epot=  3.20591e+03 Fmax= 6.19379e+02, atom= 9155
Step=   51, Dmax= 1.9e-06 nm, Epot=  3.20911e+03 Fmax= 6.16887e+02, atom= 2765
Stepsize too small, or no change in energy.
Converged to machine precision,
but not to the requested precision Fmax < 1

Double precision normally gives you higher accuracy.

writing lowest energy coordinates.

Steepest Descents converged to machine precision in 52 steps,
but did not reach the requested Fmax < 1.
Potential Energy  =  3.2059077e+03
Maximum force     =  6.1937860e+02 on atom 915
Norm of force     =  1.7839124e+02

gcq#105: "Take Your Medications and Preparations and Ram It Up Your Snout" (F. Zappa)

В результате оптимизации свободная энергия системы уменьшилась на 2 порядка

Добавим в систему молекулы воды и оптимизируем (утряска воды).

In [ ]:
> genbox -cp b_em -p b -cs spc216 -o b_s
> grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
> mdrun -deffnm b_pr -v

Getting Loaded...
Reading file b_pr.tpr, VERSION 4.5.5 (single precision)
Starting 8 threads
Loaded with Money

Making 2D domain decomposition 4 x 1 x 2


starting mdrun 'bilayer in water'
1000 steps,      0.2 ps.
step 900, remaining runtime:     0 s          imb F 23%
NOTE: Turning on dynamic load balancing


Writing final coordinates.

step 1000, remaining runtime:     0 s
 Average load imbalance: 18.6 %
 Part of the total run time spent waiting due to load imbalance: 5.8 %
 Steps where the load balancing was limited by -rdd, -rcon and/or -dds: X 0 % Z 0 %

NOTE: 5.8 % performance was lost due to load imbalance
      in the domain decomposition.


        Parallel run - timing based on wallclock.

               NODE (s)   Real (s)      (%)
       Time:      6.685      6.685    100.0
               (Mnbf/s)   (GFlops)   (ns/day)  (hour/ns)
Performance:    212.299     15.611      2.587      9.276

Посмотрим, как выглядела система до и после утряски

In [ ]:
> editconf -f b_pr.gro -o b_pr.pdb 
> editconf -f .gro -o b_s.pdb
In [7]:
cmd.reinitialize()
cmd.load('b_s.pdb')
time.sleep(10)
cmd.png('prak11_3.png')
time.sleep(10)
Image(filename='prak11_3.png')
Out[7]:
In [8]:
cmd.reinitialize()
time.sleep(5)
cmd.load('b_pr.pdb')
time.sleep(10)
cmd.png('prak11_4.png')
time.sleep(10)
Image(filename='prak11_4.png')
Out[8]:

Разницы особой не видно

На суперкомпьютере запустим тестовое моделирование молекулярной динамики

In [ ]:
> cp /home/users/golovin/progs/share/gromacs/top/residuetypes.dat .
> grompp -f md -c b_pr -p b -o b_md -maxwarn 1
> sbatch -n 4 -e error.log -o output.log -t 5 -p test impi /opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi mdrun -testverlet -deffnm b_md -v

Submitted batch job 1811968

Основное моделирование запускалось автоматически

In [ ]:
> sbatch -N1 --ntasks-per-node=2 -e error-gpu.log -o output.log -t 350 -p gpu impi /opt/ccoe/gromacs-5.0.4/build/bin/gmx_mpi mdrun -testverlet -deffnm  b_md -v

Анализ результатов молекулярной динамики

Параметры системы:

  • силовое поле: amber99sb
  • заряд системы: 0
  • Параметры ячейки:
In [ ]:
> cat b_pr.pdb | grep CRYST
CRYST1   62.600   44.430   57.780  90.00  90.00  90.00 P 1           1
  • параметры оптимизации геометрии
In [ ]:
> cat em.mdp
title               =  d_helix     ; it is title!
cpp                 =  /lib/cpp    ; my preprocessor
define              =  -DFLEX_SPC  ;
constraints         =  none        ; no constraints
integrator          =  steep       ; steepest descent algoritm for energy minimization
dt                  =  0.001       ; ps ! ;time step for integration
nstxout             =  500
nstenergy           =  500
nsteps              =  1000000        ; maximum number of steps to integrate

;      Neighbor searchering
nstlist             =  1          ; frequency to update the neighbor list
ns_type             =  grid        ;
rlist               =  0.35        ; cut-off distanse for the short-range neighbor list
rcoulomb            =  0.35        ; distanse for the Coulomb cut-off
rvdw                =  0.35        ; distanse for the LJ or Bcukingham cut-off
;
;       Energy minimizing stuff
;
emtol               =  1     ; minimization is converget when the max. force is
                                   ; smaller then this value
emstep              =  0.02        ; initial step size
  • параметры утряски воды:
In [ ]:
> cat pr.mdp
title               =  d_helix     ; it is title!
cpp                 =  /lib/cpp    ; my preprocessor
define              =  -DFLEX_SPC  ;
constraints         =  none        ; no constraints
integrator          =  steep       ; steepest descent algoritm for energy minimization
dt                  =  0.001       ; ps ! ;time step for integration
nstxout             =  500
nstenergy           =  500
nsteps              =  1000000        ; maximum number of steps to integrate

;      Neighbor searchering
nstlist             =  1          ; frequency to update the neighbor list
ns_type             =  grid        ;
rlist               =  0.35        ; cut-off distanse for the short-range neighbor list
rcoulomb            =  0.35        ; distanse for the Coulomb cut-off
rvdw                =  0.35        ; distanse for the LJ or Bcukingham cut-off
;
;       Energy minimizing stuff
;
emtol               =  1     ; minimization is converget when the max. force is
                                   ; smaller then this value
emstep              =  0.02        ; initial step size
  • параметры основного моделирования
In [ ]:
> cat md.mdp
title               = MD
cpp                 =  /lib/cpp
constraints         =  all-bonds
integrator          =  md
dt                  =  0.005   ; ps !
nsteps              =  10000000     ; total 1.0 ps.
nstcomm             =  1
nstxout             =  50000
nstvout             =  50000
nstxtcout           =  5000
;xtc-grps            = System
nstfout             =  0
nstlog              =  1000
nstenergy           =  10000
nstlist             =  10
ns_type             =  grid
coulombtype              = pme 
rlist               =  1.0
rcoulomb            =  1.0
rvdw                =  1.0
; Berendsen temperature coupling is on in two groups
Tcoupl              =  v-rescale 
tau_t               =  0.1   0.1  
tc-grps             =  DPPC    SOL
ref_t               =  320   320  
; Pressure coupling is not on
Pcoupl                   = Berendsen
Pcoupltype               = semiisotropic
nstpcouple               = -1
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p                    = 5
compressibility          = 4.5e-5  4.5e-5 4.5e-5 0.0 0.0     0.0
ref_p                    = 1.0     1.0    1.0    1.0 1.0    1.0
; Generate velocites is on at 300 K.
gen_vel             =  yes
gen_temp            =  300.0
gen_seed            =  173529

Создадим pdb файл с движением молекул

In [ ]:
> trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 5 -pbc mol
...
Select group for output
Group     0 (         System) has 11054 elements
Group     1 (          Other) has  3200 elements
Group     2 (           DPPC) has  3200 elements
Group     3 (          Water) has  7854 elements
Group     4 (            SOL) has  7854 elements
Group     5 (      non-Water) has  3200 elements
Select a group: DPPC
Selected: 'DPPC'
...

Рассчитаем скорость движения молекул в файле

In [ ]:
> cat b_pbc_1.pdb | egrep MODEL\ *50 -B3
TITLE     bilayer in water t= 24500.00000
REMARK    THIS IS A SIMULATION BOX
CRYST1   53.695   38.109   78.431  90.00  90.00  90.00 P 1           1
MODEL       50

50 модель отражает временную точку 24500 пикосекунд, то есть мы имеем 500 пкс на кадр.

Cохраним видео в формате mpg, конвертируем в mp4 и покажем.

Как можно заметить, ось Z является нормалью к поверхности бислоя.

Получим координаты ячеек в каждый момент времени

In [ ]:
> g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
...
Reading file b_md.tpr, VERSION 4.5.4 (single precision)
Reading file b_md.tpr, VERSION 4.5.4 (single precision)
Group     0 (         System) has 11054 elements
Group     1 (          Other) has  3200 elements
Group     2 (           DPPC) has  3200 elements
Group     3 (          Water) has  7854 elements
Group     4 (            SOL) has  7854 elements
Group     5 (      non-Water) has  3200 elements
Select a group: DPPC
Selected 2: 'DPPC'
Last frame       1000 time 25000.000

gcq#310: "Shoot them in the back now" (The Ramones)

Построим график зависимости площади бислоя от времени

In [2]:
import pandas as pd
import matplotlib.pyplot as plt
In [13]:
s = pd.read_csv('geraseva1/box_1.xvg', sep='\t', skiprows=24, header=None )
In [30]:
time=s[1]
square=s[2]*s[3]/32

plt.plot(time, square)
plt.xlim(0,25000)
plt.ylim(0.6,0.9)
Out[30]:
(0.6, 0.9)

Соответственно, сборка бислоя происходит через примерно 15000 пкс, а площадь бислоя в расчете на 1 липид составляет 0.65 кв. ангстрем.

Определим изменение гидрофобной и гидрофильной поверхности в ходе самосборки и построим график.

In [ ]:
> g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg
...
Reading frame       0 time    0.000   Select a group for calculation of surface and a group for output:
Group     0 (         System) has 11054 elements
Group     1 (          Other) has  3200 elements
Group     2 (           DPPC) has  3200 elements
Group     3 (          Water) has  7854 elements
Group     4 (            SOL) has  7854 elements
Group     5 (      non-Water) has  3200 elements
Select a group: DPPC
Selected 2: 'DPPC'

WARNING: Masses and atomic (Van der Waals) radii will be guessed
         based on residue and atom names, since they could not be
         definitively assigned from the information in your input
         files. These guessed numbers might deviate from the mass
         and radius of the atom type. Please check the output
         files if necessary.

WARNING: could not find a Van der Waals radius for 64 atoms
1920 out of 3200 atoms were classified as hydrophobic
Last frame       1000 time 25000.000


gcq#238: "Uh-oh" (Tinky Winky)

> cat sas_b.xvg | egrep -v '#|@' | awk '{print $1","$2","$3","$4}' > sas_b1.xvg
In [66]:
s = pd.read_csv('geraseva1/sas_b1.xvg', header=None)
time=s[0]
hphil=s[2]
hphob=s[1]
all=s[3]

plt.plot(time, hphil,c='blue')
plt.plot(time, hphob,c='red')
plt.plot(time, all, c='black')
plt.xlim(0,25000)
Out[66]:
(0, 25000)

В процессе самосборки общая поверхность уменьшается. Это касается конкретно гидрофобной и гидрофильной поверхностей. Однако уменьшение площади гидрофильной поверхности происходит медленнее. Также можно заметить увеличение гидрофильной площади в самом начале, за счет чего резко уменьшается гидрофобная площадь.