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

Вгрузим же

In [ ]:
! wget https://kodomo.fbb.msu.ru/~golovin/bilayer/dppc.itp
! wget https://kodomo.fbb.msu.ru/~golovin/bilayer/lipid.itp
! wget https://kodomo.fbb.msu.ru/~golovin/bilayer/dppc.gro
! wget https://kodomo.fbb.msu.ru/~golovin/bilayer/b.top
! wget https://kodomo.fbb.msu.ru/~golovin/bilayer/em.mdp
! wget https://kodomo.fbb.msu.ru/~golovin/bilayer/pr.mdp
! wget https://kodomo.fbb.msu.ru/~golovin/bilayer/md.mdp

Перешли в рабочую директорию

In [9]:
! cd ../zhenia/
! source /home/preps/golovin/progs/bin/GMXRC.bash

На основе одного липида создали ячейку с 64 липидами

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

Создали pdb файлы

In [11]:
! editconf -f dppc.gro -o dppc.pdb
! editconf -f b_64.gro -o b_64.pdb
Read 50 atoms
Volume: 1.5477 nm^3, corresponds to roughly 600 electrons
No velocities found
Read 3200 atoms
Volume: 99.0529 nm^3, corresponds to roughly 44500 electrons
No velocities found
                         :-)  G  R  O  M  A  C  S  (-:

                       GRowing Old MAkes el Chrono Sweat

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                               :-)  editconf  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f ../Perevoshchikova/dppc.gro  Input        Structure file: gro g96 pdb
                                   tpr etc.
  -n      index.ndx  Input, Opt.  Index file
  -o ../Perevoshchikova/dppc.pdb  Output, Opt! Structure file: gro g96 pdb
                                   etc.
-mead      mead.pqr  Output, Opt. Coordinate file for MEAD
 -bf      bfact.dat  Input, Opt.  Generic data file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    0       Set the nicelevel
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-[no]ndef    bool   no      Choose output from default index groups
-bt          enum   triclinic  Box type for -box and -d: triclinic, cubic,
                            dodecahedron or octahedron
-box         vector 0 0 0   Box vector lengths (a,b,c)
-angles      vector 90 90 90  Angles between the box vectors (bc,ac,ab)
-d           real   0       Distance between the solute and the box
-[no]c       bool   no      Center molecule in box (implied by -box and -d)
-center      vector 0 0 0   Coordinates of geometrical center
-aligncenter vector 0 0 0   Center of rotation for alignment
-align       vector 0 0 0   Align to target vector
-translate   vector 0 0 0   Translation
-rotate      vector 0 0 0   Rotation around the X, Y and Z axes in degrees
-[no]princ   bool   no      Orient molecule(s) along their principal axes
-scale       vector 1 1 1   Scaling factor
-density     real   1000    Density (g/L) of the output box achieved by
                            scaling
-[no]pbc     bool   no      Remove the periodicity (make molecule whole again)
-resnr       int    -1       Renumber residues starting from resnr
-[no]grasp   bool   no      Store the charge of the atom in the B-factor
                            field and the radius of the atom in the occupancy
                            field
-rvdw        real   0.12    Default Van der Waals radius (in nm) if one can
                            not be found in the database or if no parameters
                            are present in the topology file
-[no]sig56   bool   no      Use rmin/2 (minimum in the Van der Waals
                            potential) rather than sigma/2 
-[no]vdwread bool   no      Read the Van der Waals radii from the file
                            vdwradii.dat rather than computing the radii
                            based on the force field
-[no]atom    bool   no      Force B-factor attachment per atom
-[no]legend  bool   no      Make B-factor legend
-label       string A       Add chain label for all residues
-[no]conect  bool   no      Add CONECT records to a .pdb file when written.
                            Can only be done when a topology is present


gcq#36: "All Work and No Play Makes Jack a Dull Boy" (The Shining)

                         :-)  G  R  O  M  A  C  S  (-:

                       GRowing Old MAkes el Chrono Sweat

                            :-)  VERSION 4.5.5  (-:

        Written by Emile Apol, Rossen Apostolov, Herman J.C. Berendsen,
      Aldert van Buuren, Pär Bjelkmar, Rudi van Drunen, Anton Feenstra, 
        Gerrit Groenhof, Peter Kasson, Per Larsson, Pieter Meulenhoff, 
           Teemu Murtola, Szilard Pall, Sander Pronk, Roland Schulz, 
                Michael Shirts, Alfons Sijbers, Peter Tieleman,

               Berk Hess, David van der Spoel, and Erik Lindahl.

       Copyright (c) 1991-2000, University of Groningen, The Netherlands.
            Copyright (c) 2001-2010, The GROMACS development team at
        Uppsala University & The Royal Institute of Technology, Sweden.
            check out http://www.gromacs.org for more information.

         This program is free software; you can redistribute it and/or
          modify it under the terms of the GNU General Public License
         as published by the Free Software Foundation; either version 2
             of the License, or (at your option) any later version.

                               :-)  editconf  (-:

Option     Filename  Type         Description
------------------------------------------------------------
  -f ../Perevoshchikova/b_64.gro  Input        Structure file: gro g96 pdb
                                   tpr etc.
  -n      index.ndx  Input, Opt.  Index file
  -o ../Perevoshchikova/b_64.pdb  Output, Opt! Structure file: gro g96 pdb
                                   etc.
-mead      mead.pqr  Output, Opt. Coordinate file for MEAD
 -bf      bfact.dat  Input, Opt.  Generic data file

Option       Type   Value   Description
------------------------------------------------------
-[no]h       bool   no      Print help info and quit
-[no]version bool   no      Print version info and quit
-nice        int    0       Set the nicelevel
-[no]w       bool   no      View output .xvg, .xpm, .eps and .pdb files
-[no]ndef    bool   no      Choose output from default index groups
-bt          enum   triclinic  Box type for -box and -d: triclinic, cubic,
                            dodecahedron or octahedron
-box         vector 0 0 0   Box vector lengths (a,b,c)
-angles      vector 90 90 90  Angles between the box vectors (bc,ac,ab)
-d           real   0       Distance between the solute and the box
-[no]c       bool   no      Center molecule in box (implied by -box and -d)
-center      vector 0 0 0   Coordinates of geometrical center
-aligncenter vector 0 0 0   Center of rotation for alignment
-align       vector 0 0 0   Align to target vector
-translate   vector 0 0 0   Translation
-rotate      vector 0 0 0   Rotation around the X, Y and Z axes in degrees
-[no]princ   bool   no      Orient molecule(s) along their principal axes
-scale       vector 1 1 1   Scaling factor
-density     real   1000    Density (g/L) of the output box achieved by
                            scaling
-[no]pbc     bool   no      Remove the periodicity (make molecule whole again)
-resnr       int    -1       Renumber residues starting from resnr
-[no]grasp   bool   no      Store the charge of the atom in the B-factor
                            field and the radius of the atom in the occupancy
                            field
-rvdw        real   0.12    Default Van der Waals radius (in nm) if one can
                            not be found in the database or if no parameters
                            are present in the topology file
-[no]sig56   bool   no      Use rmin/2 (minimum in the Van der Waals
                            potential) rather than sigma/2 
-[no]vdwread bool   no      Read the Van der Waals radii from the file
                            vdwradii.dat rather than computing the radii
                            based on the force field
-[no]atom    bool   no      Force B-factor attachment per atom
-[no]legend  bool   no      Make B-factor legend
-label       string A       Add chain label for all residues
-[no]conect  bool   no      Add CONECT records to a .pdb file when written.
                            Can only be done when a topology is present


gcq#36: "All Work and No Play Makes Jack a Dull Boy" (The Shining)

Теперь посмотрим на то, что у нас получилось

In [3]:
import __main__
from IPython.display import Image
__main__.pymol_argv = [ 'pymol', '-c' ] # открыть pymol  без использования внешнего GUI
import pymol
pymol.finish_launching()
from pymol import cmd,stored
import numpy as np
import time
import pandas as pd
import matplotlib.pyplot as plt
In [11]:
cmd.load('dppc.pdb')
cmd.bg_color('white')
cmd.do('''set stick_radius, 0.5
show stick, all
orient dppc
zoom dppc, complete = 1
''')
cmd.do('png lipid.png, width=1080, height=1000, ray=1')
In [12]:
Image(filename='lipid.png')
Out[12]:

Видим, что у меня истекла лицензия, а файл dppc содержит структуру фосфатидилхолина у которого одна из жирных кислот содержит циклопентановый элемент.

In [15]:
cmd.delete("all")
cmd.bg_color('white')
cmd.load('b_64.pdb')
cmd.do('''as cartoon, all
set stick_radius, 0.5
show stick, all
zoom b_64, complete = 1
turn x, 20
turn y, 20
turn z, 20'''
cmd.do('png lipids.png, width=1080, height=1000, ray=1')
In [16]:
Image(filename='lipids.png')
Out[16]:

Получилась ячейка!

В текстовом редакторе в файле b.top установили правильное количество липидов в системе - 64. В файле с параметрами для молекулярной динамики установили шаг по времени 0.004 ps.

Сделали небольшой отступ в ячейке от липидов, что бы добавить примерно 2500 молекул воды.

In [ ]:
! gmx editconf -f b_64.gro -o b_ec -d 0.5 

Проведём оптимизацию геометрии системы, что бы удалить "плохие" контакты молекул

In [ ]:
! gmx grompp -f em -c b_ec -p b -o b_em -maxwarn 2
! gmx mdrun -deffnm b_em -v

Потенциальная энергия в ходе 68 шагов уменьшилась с 4.63573e+05 Кдж/моль до -3.11687e+04 Кдж/моль.

Добавим в ячейку молекулы воды

In [ ]:
! gmx solvate -cp b_em -p b -cs spc216 -o b_s

Проведём "утряску" воды:

In [ ]:
%%bash
gmx grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
gmx mdrun -deffnm b_pr -v
#если водяной взрыв, то
gmx grompp -f em -c b_s -p b -o b_empr -maxwarn 1
gmx mdrun -deffnm b_empr -v

gmx grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
gmx mdrun -deffnm b_pr -v

Теперь переведем файл с добавленной водой b_s.gro и файл с оптимизированной водой b_pr.gro в pdb формат и посмотрим, как это выглядит

In [ ]:
! gmx editconf -f b_s.gro -o b_s.pdb
! gmx editconf -f b_pr.gro -o b_pr.pdb

Посмотрим на бислой до "утряски" воды

In [18]:
cmd.do("reinitialize")
cmd.delete("all")
cmd.bg_color('white')
cmd.load('b_s.pdb')
cmd.do('''as cartoon, all
set stick_radius, 0.5
show stick, all
show nonbonded
zoom b_s, complete = 1
turn x, 20
turn y, 20
turn z, 20''')
cmd.do('png before.png, width=1080, height=1000, ray=1')
In [19]:
Image(filename='before.png')
Out[19]:

... и после "утряски" воды

In [20]:
cmd.do("reinitialize")
cmd.delete("all")
cmd.bg_color('white')
cmd.load('b_pr.pdb')
cmd.do('''as cartoon, all
set stick_radius, 0.5
show stick, all
show nonbonded
zoom b_pr, complete = 1
turn x, 20
turn y, 20
turn z, 20''')
cmd.do('png after.png, width=1080, height=1000, ray=1')
In [21]:
Image(filename='after.png')
Out[21]:

В структуре после оптимизации молекулы воды расположены более равномерно и, соответственно, более оптимально для образования водородных связей. В структуре до оптимизации молекулы воды расположены чуть более "клочками".