Moleküler Dinamik Simülasyonları Bileşenleri

EpiSTEM Türkiye'deki yazıyı takiben biraz daha detaya girelim nedir bu simülasyonlar. Bileşenleri nelerdir?

Yazımda da değindiğim gibi Newton'un kanunlarından yola çıkılmıştır moleküler dinamik simülasyonları için. Bir cismin hareketi ona uygulanan kuvvet ile orantılıdır. Elinizle bir objeyi ne kadar uzağa fırlatabileceğiniz ona uyguladığınız kuvvet ile belirlenir. Moleküler dinamik simülasyonlarında da bir molekülün hareketi de ona uygulanan kuvvete bağlıdır. Bu kuvveti tanımlarken molekülün sahip olduğu geometrik yapıyı ve çevresini saran diğer molekülleri baz alıyoruz. İşte bu aşamada kuvvet alanı dediğimiz bir grup özellikler devreye giriyor. Nasıl insanların kol uzunlukları farklı ise, atomları birbirine bağlayan bağlarda farklı uzunluklara sahiptirler. Kimimizin kolu çabuk kırılırken, kimimiz sağlam. Bu bağların dayanıklılığı da aynı böyle değişkendir. Meraklısına aşağıdaki grafik bir kuvvet alanındaki tanımlı özellikleri özetliyor. Grafiğin uyarlandığı kaynağı okumanızı tavsiye ederim (Maalesef İngilizce).

Kuvvet alanı moleküler dinamik simülasyonlarının en temel ve en önemli bileşenleridir. Kuvvet alanına ek olarak bilmemiz gereken bir diğer bileşen ise koordinat. Yani atomların adresleri. Hepimizin bir ev ve iş adresi var. Koordinatlarıda atomların yerini belirleyen adresler olarak düşünebilirsiniz. Bu koordinatlar elbette 3 boyutlu. Basit moleküllerin koordinatlarını tahmin yolu ile çizmek mümkün olabilir. Fakat 5 tane yapı taşından (amino asitten) oluşan proteinin bile alabileceği çok sayıda şekil var. O nedenle deneysel yöntemlere başvuruyoruz bu koordinatları elde ederken.

İhtiyacımız olan üçüncü bileşen ise simülasyonu çalıştıracağımız program. Benim deneyimlerime göre GROMACS simülasyonlara yeni başlayanlar için en güvenilir program. Komut isteminden kullanmak birçok kişiyi ürkütse de yeni başlayan biri olarak yapacağınız hataların geri bildirimini almak çok önemli ve GROMACS bu konuda en iyilerden biri diyebilirim. Kullanabileceğiniz diğer bir alternatif NAMD, ki bunu VMD programından bir arayüz ile kullanabilirsiniz. Görsel açıdan öğrenmesi çok daha kolay bir program. Bir diğeri ise bir çok özelliği barındıran bir nevi kodlama dili diyebileceğimiz CHARMM; simülasyon yapabilirsiniz ama öğrenmesi kolay olmayan bir programdır kendisi (hata mesajlarının yoksunluğu da cabası). Ama GROMACS'a kıyasla sistemler üzerinde değişiklik yapması daha kolay bir program. Amber ismi verilen programına dair bir fikrim yok, çünkü hiç kullanmadım, ama ismi burada dursun ki sizin de haberiniz olsun. Lakin Amber'in aynı isimli kuvvet alanını kullandım ki bu rehberde onun üzerine. Alternatif kuvvet alanı arıyorsanız, CHARMM isimli (yukarıda bahsettiğim programı geliştirenler tarafından hazırlanan) kuvvet alanına bakabilirsiniz.

Programı belirledik, GROMACS. Kuvvet alanımız ise, Amber. Kaldı geriye koordinatlar... Bunun için PDB veri bankasına gidip PDB formatında bir protein yapısı indirdim: İnsülin. Linkten daha fazla bilgi edinebilirsiniz. Bu proteini seçmemin sebebi, muhtemelen bir çoğunuzun ismini duymuş olması. Şeker hastası var mı aramızda? Kendinize enjekte etmekte olduğunuz insülin ve insülin benzeri proteinler sayesinde hayat kaliteniz biraz daha artıyor, artık kendisi ile moleküler düzeyde de tanışmış oldunuz.

In [ ]:
import nglview #NGLview paketini ortama yükle (https://github.com/arose/nglview)
view = nglview.show_pdbid("2ins")  #"2ins" kodlu yapıyı PDB veri bankasından indir ve NGLview uygulamasına gönder
view #görüntüle
In [24]:
%%html --isolated
<iframe src="2INS.html" allowTransparency="true" scrolling="no" width=100% height=400px> </iframe>'
'

Yukarıda görmüş olduğunuz yapı insülinin aktif formu. Proteinlerin aktif formları birbirinin aynı veya farklı bir çok zincirden oluşabilir. Hız kazanmak açısından bu rehberde biz sadece aşağıda göreceğiniz tek bir zincir ile devam edeceğiz. Bu arada üstteki resim hareketli. Başka bir html sayfayı yükleyerek size bu interaktif görseli sundum, ama eğer bu rehberi indirip kendiniz çalıştıracaksanız html ile başlayan hücreleri silip python komutundan devam edebilirsiniz (Aaa söylemedim mi bu rehber Python programlama dili ile çalışıyor).

Bu rehberi hazırlarken kullandığım programlar ile ilgili bilgileri yazının en sonunda biraz anlattım. Benim hazırladığım dosyalara ve bu programlara nereden ulaşabileceğiniz de orada.

1. Simülasyon ortamının hazırlanması

a. Protein yapısının temizlenmesi

Her ne kadar sizi burada angarya dediğim bu kısım ile uğraştırmayacak olsam da bilmeniz gerek. Protein veri bankasından indirip de öyle hemen simülasyon başlatabileceğiniz yapı sayısı azdır. Bunun en birincil nedeni deneysel yöntemlerin her atoma konum atamakta eksik kalmasından kaynaklı. Günümüzde detayına inmediğim bu yöntemler (X-ray, NMR... vs.) daha yüksek kalitede protein resimleri çizebiliyor olsa da proteinin her bölgesi görüntülenemeyebiliyor. İkincil neden protein ile beraber başka moleküllerin de resminin çıkması ve bizim bunları simülasyonda istemiyor olmamız. Bu moleküllerden bazıları deneyin başarılı olması için gerekli (mesela deterjanlar), bazıları da proteinlere bağlanan ilaç gibi diğer bileşenlerdir. Genelde deterjan gibi bileşenleri simülasyonda istemiyoruz, çünkü proteinin doğal ortamında aslında onlar yok. İlaç benzeri molekülleri ise kuvvet alanında tanımları olmadığı için çıkarmanız gerekebilir (veya tanımlayarak kuvvet alanına eklememiz gerekir). Diğer nedenler ki bunlar aslında ikincil nedenden daha fazla sıklıkta görülüyor, proteinin yapı taşlarının (amino asitlerin) numaralandırılmasının farklı olması, bazı bağların eksik olması, amino asitlerin eksik olması ama eksik olarak rapor edilmemiş olması, vs. vs. Hayatınızı kolaylaştırmanın tek yolu bir protein ailesi ile çalışmaya başlamadan önce numaralandırma ve doğal amino asit dizisi kontrollerini yapıp aklı selim bir rota belirlemeniz. Ne çok laf ettin diyeceksiniz, inanın bu yapıları hazırlamak simülasyonların çalışmasından daha uzun sürüyor. Bu aşamada bir hata yapmış iseniz... hoooop... bütün veriler çöpe. Ben size aşağıda, buyurun, hazır temizlenmişini sunuyorum.

In [ ]:
import nglview
view = nglview.show_file("2ins_zincirAB.pdb")
view #görüntüle
In [25]:
%%html
<iframe src="2INS_zincir.html" allowTransparency="true" scrolling="no" width=100% height=400px> </iframe>

b. Protein yapısının kuvvet alanı ile hazırlanması

Fark ettiniz mi hala iki zincir var insülinde. Aslında tek bir peptit zincirinden oluşur insülin, yani A ve B zincirleri tek bir gen tarafından üretilmiştir, ancak insülün aktif halini almadan önce bazı değişimlerden geçer ve sonucunda sanki birbirinden bağımsızmış gibi iki zincir oluşur. Biz şimdi bu iki zinciri kuvvet alanındaki tanımları ile "topoloji" denilen bir dosyaya yeniden yazacağız. PDB gibi dosyalar proteindeki atomların uzaydaki yerini belirlerken, topoloji bu atomların birbirleri ve çevresi ile ilişkisi hakkında nicel bilgiler içerir. Karbon atomu, oksijenden 0.15 nanometre uzakta (evet geçiniz kilometreyi santimetreyi biz nano boyuttayız artık, yani $10^{-9}$m) ve oksijen dışında 1'i hidrojen olmak üzere 3 atoma daha bağlı. Karbon-hidrojen-oksijen açısı 120 derece. ... Bu şekilde bir sürü bilgi var, bunların el ile hazırlanması mümkün olsa da, aklımızı peynir ekmek ile yemedik bu devirde. Sağ olsunlar simülasyon programlarını yazanlar bize bu dosyaları da hazırlayacak bir bileşen eklemişler. Aşağıdaki komut ile bu dosyayı çok kolay hazırlayabiliriz. Bu komutları bir komut istemi (terminalde) doğrudan çalıştırıyorum. Farklı terminaller var, benim burada kullandığım terminam bash.

In [26]:
%%bash
export GMX_MAXBACKUP=-1
pdb2gmx -f 2ins_zincirAB.pdb -o kuvvetli_protein.pdb -p topoloji.top -ff amber99sb-ildn -water tip3p -merge all -quiet
Using the Amber99sb-ildn force field in directory amber99sb-ildn.ff

Reading 2ins_zincirAB.pdb...
Read 393 atoms
Analyzing pdb file
Splitting chemical chains based on TER records or chain id changing.

Merged chains into joint molecule definitions at 1 places.

There are 1 chains and 0 blocks of water and 30 residues with 393 atoms

  chain  #res #atoms
  1 'A'    50    393  

Reading residue database... (amber99sb-ildn)
Processing chain 1 'A' (393 atoms, 50 residues)
Identified residue GLY1 as a starting terminus.
Identified residue ASN21 as a ending terminus.
Identified residue VAL2 as a starting terminus.
Identified residue ALA30 as a ending terminus.
Checking for duplicate atoms....
Now there are 389 atoms. Deleted 4 duplicates.
Generating any missing hydrogen atoms and/or adding termini.
Now there are 50 residues with 755 atoms
		--------- PLEASE NOTE ------------
You have successfully generated a topology from: 2ins_zincirAB.pdb.
The Amber99sb-ildn force field and the tip3p water model are used.
		--------- ETON ESAELP ------------

gcq#355: "I like to wait, then I feel like I do something" (Carl Caleman)

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

               Gromacs Runs One Microsecond At Cannonball Speeds

                            :-)  VERSION 4.6.5  (-:

        Contributions from Mark Abraham, Emile Apol, Rossen Apostolov, 
           Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar,  
     Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans, 
        Peter Kasson, Carsten Kutzner, 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-2012,2013, 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 Lesser General Public License
        as published by the Free Software Foundation; either version 2.1
             of the License, or (at your option) any later version.

                               :-)  pdb2gmx  (-:

Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/aminoacids.r2b
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/dna.r2b
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/rna.r2b

WARNING: there were 0 atoms with zero occupancy and 8 atoms with
         occupancy unequal to one (out of 393 atoms). Check your pdb file.

Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/atomtypes.atp
Atomtype 1
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/aminoacids.rtp
Residue 93
Sorting it all out...
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/dna.rtp
Residue 109
Sorting it all out...
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/rna.rtp
Residue 125
Sorting it all out...
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/aminoacids.hdb
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/dna.hdb
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/rna.hdb
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/aminoacids.n.tdb
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/aminoacids.c.tdb
Analysing hydrogen-bonding network for automated assignment of histidine
 protonation. 75 donors and 72 acceptors were found.
There are 103 hydrogen bonds
Will use HISE for residue 5
Will use HISE for residue 10
8 out of 8 lines of specbond.dat converted successfully
Special Atom Distance matrix:
                    CYS6    CYS7   CYS11   CYS20    HIS5    CYS7   HIS10
                    SG43    SG49    SG73   SG151  NE2194   SG208  NE2228
    CYS7    SG49   0.867
   CYS11    SG73   0.214   0.942
   CYS20   SG151   1.196   1.944   1.240
    HIS5  NE2194   0.906   0.445   0.880   2.049
    CYS7   SG208   0.878   0.188   0.918   1.937   0.364
   HIS10  NE2228   1.307   0.946   1.312   1.964   1.035   0.800
   CYS19   SG298   1.116   1.813   1.169   0.197   1.935   1.800   1.786
                
                
Linking CYS-6 SG-43 and CYS-11 SG-73...
Linking CYS-7 SG-49 and CYS-7 SG-208...
Linking CYS-20 SG-151 and CYS-19 SG-298...
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/aminoacids.arn
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/dna.arn
Opening force field file /Users/cattibrie_fr/anaconda/anaconda/share/gromacs/top/amber99sb-ildn.ff/rna.arn
Making bonds...
Number of bonds was 767, now 765
Generating angles, dihedrals and pairs...
Before cleaning: 1993 pairs
Before cleaning: 2051 dihedrals
Keeping all generated dihedrals
Making cmap torsions...There are 2051 dihedrals,  162 impropers, 1372 angles
          1975 pairs,      765 bonds and     0 virtual sites
Total mass 5584.378 a.m.u.
Total charge -2.000 e
Writing topology

Writing coordinate file...

Bu komutta -f ile girdiyi, -o ile çıktıyı tanımladık (bu atomları kuvvet alanına göre isimlendirilmiş bir dosya). -p topoloji dosyasını tanımladı, -ff kuvvet alanını, -water su modelini seçti (anlaşılması en zor molekül sudur, birden fazla modeli vardır. Onur'un yazısından daha fazla bilgi edinebilirisniz su ile ilgili. ). -merge all ise iki zincirin topoloji bilgisini tek bir dosyaya (topoloji.top) yazdı. Üstteki seyir defterini (log) incelerseniz, insülinin toplam yükünü (charge) öğrenebilirsiniz.

c. Simülasyon kutusu ve çözelti

Simülasyonları bir deney tüpü olarak düşünebilirsiniz. Sadece bir sürü protein ile aynı anda deney yapmak yerine biz her bir proteini tek tek tüplere yerleştiriyoruz. Bu tüplerin sınırını biz çiziyoruz. Gönül isterdi ki kocaman hazırlayalım bu simülasyonları ama çok pahalı bir hesaplama olduğu için mümkün olduğunca ufak hazırlıyoruz. Aşağıdaki komutlar önce kutunun sınırlarını çizip, sonra da o kutuyu su ile dolduracak. Bir önceki aşamada seçtiğimiz su modeli de burada devreye giriyor. Kutunun büyüklüğü proteinin etrafında yaklaşık 1 nanometre boşluk kalacak şekilde 5.4 nanometre olarak seçildi. Su ekledikten sonra ise sistemin toplam yükünü sıfırlıyor (1) ve belirli miktarda bir çözelti hazırlıyoruz (2).

1) Atomlar arasındaki etkileşimi hesaplamada kullanılan yöntem sistemin toplam yükünün sıfır olmasını gerektiriyor. Bu yük dengeleme işi sisteme iyonlar eklenerek yapılıyor. İnsülinin yükünü öğrendi iseniz sisteme artı yüklü sodyum ($Na^{+1}$) mu yoksa eksi yüklü klorür ($Cl^{-1}$) mü eklememiz gerekiyor tahmin edebilirsiniz, çünkü kutuya eklenen suyun toplam yükü sıfır. Yani sistem insülinin yükünü taşıyor.

2) İnsülinin içinde bulunduğu doğal ortamı sistemi mümkün olduğunca basit tutarak hazırlıyoruz. Bunun en temel yolu iste sisteme daha fazla iyon yüklemek. Yani sadece net yükü sıfırlamakla kalmıyoruz, üzerine daha fazla sayıda sodyum ve klorür ekliyoruz.

-p topoloji.top bu dosyanın devamlı güncellenmesini sağlıyor, böylece sistemde ekleme/çıkarmalar sonucunda ne kaldı biliyoruz. Bu işlem önemli çünkü topoloji dosyasının PDB dosyası ile uyuşması gerek.

In [27]:
%%bash
export GMX_MAXBACKUP=-1
editconf -f kuvvetli_protein.pdb -box 5.4 5.4 5.4 -bt cubic -o kutuda_protein.pdb -quiet
genbox -cp kutuda_protein.pdb -p topoloji.top -cs -o sulu_protein.pdb -quiet
Read 755 atoms
No velocities found
    system size :  2.478  2.554  3.217 (nm)
    diameter    :  3.356               (nm)
    center      : -0.996  1.065  0.534 (nm)
    box vectors :  0.000  0.000  0.000 (nm)
    box angles  :   0.00   0.00   0.00 (degrees)
    box volume  :   0.00               (nm^3)
    shift       :  3.696  1.635  2.166 (nm)
new center      :  2.700  2.700  2.700 (nm)
new box vectors :  5.400  5.400  5.400 (nm)
new box angles  :  90.00  90.00  90.00 (degrees)
new box volume  : 157.46               (nm^3)

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.

Neighborsearching with a cut-off of 0.48
Table routines are used for coulomb: FALSE
Table routines are used for vdw:     FALSE
Cut-off's:   NS: 0.48   Coulomb: 0.48   LJ: 0.48
System total charge: 0.000
Potential shift: LJ r^-12: 0.000 r^-6 0.000, Coulomb 0.000

Grid: 12 x 12 x 12 cells
Adding line for 4924 solvent molecules to topology file (topoloji.top)
                         :-)  G  R  O  M  A  C  S  (-:

               GRoups of Organic Molecules in ACtion for Science

                            :-)  VERSION 4.6.5  (-:

        Contributions from Mark Abraham, Emile Apol, Rossen Apostolov, 
           Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar,  
     Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans, 
        Peter Kasson, Carsten Kutzner, 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-2012,2013, 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 Lesser General Public License
        as published by the Free Software Foundation; either version 2.1
             of the License, or (at your option) any later version.

                               :-)  editconf  (-:


gcq#26: "Art For Arts Sake, Money For Gods Sake" (10 CC)

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

               GRoups of Organic Molecules in ACtion for Science

                            :-)  VERSION 4.6.5  (-:

        Contributions from Mark Abraham, Emile Apol, Rossen Apostolov, 
           Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar,  
     Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans, 
        Peter Kasson, Carsten Kutzner, 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-2012,2013, 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 Lesser General Public License
        as published by the Free Software Foundation; either version 2.1
             of the License, or (at your option) any later version.

                                :-)  genbox  (-:

Reading solute configuration
Gromacs Runs One Microsecond At Cannonball Speeds
Containing 755 atoms in 50 residues
Initialising van der waals distances...
Reading solvent configuration
"216H2O,WATJP01,SPC216,SPC-MODEL,300K,BOX(M)=1.86206NM,WFVG,MAR. 1984"
solvent configuration contains 648 atoms in 216 residues

Initialising van der waals distances...
Will generate new solvent configuration of 3x3x3 boxes
Generating configuration
Sorting configuration
Found 1 molecule type:
    SOL (   3 atoms):  5832 residues
Calculating Overlap...
box_margin = 0.315
Removed 0 atoms that were outside the box
Successfully made neighbourlist
nri = 27112, nrj = 507986
Checking Protein-Solvent overlap: tested 16068 pairs, removed 786 atoms.
Checking Solvent-Solvent overlap: tested 199501 pairs, removed 1938 atoms.
Added 4924 molecules
Generated solvent containing 14772 atoms in 4924 residues
Writing generated configuration to sulu_protein.pdb
Gromacs Runs One Microsecond At Cannonball Speeds

Output configuration contains 15527 atoms in 4974 residues
Volume                 :     157.464 (nm^3)
Density                :      998.24 (g/l)
Number of SOL molecules:   4924   

Processing topology

gcq#26: "Art For Arts Sake, Money For Gods Sake" (10 CC)

İyonları eklemeden önce sisteme bakalım mı? Üstteki seyir defterinde kaç tane su molekülü eklenmiş bulabilirsiniz. Aşağıda da sistemin görsel halini inceleyebilirsiniz...

In [ ]:
import nglview #bu komut üstte çalıştığı için tekrar çalıştırmaya gerek yok ama üstteki veri kayboldu ise tekrar yükledim.
view = nglview.show_file("sulu_protein.pdb")
view.add_ball_and_stick("sol",opacity="0.5")
view #görüntüle
In [28]:
%%html --isolated
<iframe src="Simulasyon_sistemi_yuklu.html" allowTransparency="true" scrolling="no" width=100% height=600px> </iframe>

İyonları eklerken GROMACS'ın bir gerekliliği var, o da çalıştırma dosyası yaratmak. Bunun detaylarına girmiyorum, her simülasyon programının farklı bir yolu var bunun için. Ama aynı işlevi görüyorlar: rastgele su moleküllerini iyonlar ile değiştirmek. Yani bu aşamada sistemden su çıkararak iyonlara yer açıyoruz. Eklediğim iyon çözeltisinin miktarını aşağıdaki komut isteminden tahmin edebilirsiniz (birim mol/litre).

In [29]:
%%bash
export GMX_MAXBACKUP=-1
touch 4ion.mdp
grompp -f 4ion.mdp -c sulu_protein.pdb -p topoloji.top -o 4ion.tpr -quiet
genion -s 4ion.tpr -conc 0.15 -neutral -o cok_iyonize_protein.pdb -p topoloji.top -quiet << EOF
13
EOF
Analysing residue names:
There are:    50    Protein residues
There are:  4924      Water residues
Analysing Protein...
Largest charge group radii for Van der Waals: 0.039, 0.039 nm
Largest charge group radii for Coulomb:       0.084, 0.084 nm
This run will generate roughly 1 Mb of data
Will try to add 16 NA ions and 14 CL ions.
Select a continuous group of solvent molecules
Selected 13: 'SOL'

Processing topology
Replacing 30 solute molecules in topology file (topoloji.top)  by 16 NA and 14 CL ions.
                         :-)  G  R  O  M  A  C  S  (-:

                   GROningen MAchine for Chemical Simulation

                            :-)  VERSION 4.6.5  (-:

        Contributions from Mark Abraham, Emile Apol, Rossen Apostolov, 
           Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar,  
     Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans, 
        Peter Kasson, Carsten Kutzner, 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-2012,2013, 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 Lesser General Public License
        as published by the Free Software Foundation; either version 2.1
             of the License, or (at your option) any later version.

                                :-)  grompp  (-:


NOTE 1 [file 4ion.mdp]:
  You are using a cut-off for VdW interactions with NVE, for good energy
  conservation use vdwtype = Shift (possibly with DispCorr)


NOTE 2 [file 4ion.mdp]:
  You are using a cut-off for electrostatics with NVE, for good energy
  conservation use coulombtype = PME-Switch or Reaction-Field-zero

Generated 2211 of the 2211 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 2211 of the 2211 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type 'Protein_chain_A'
Excluding 2 bonded neighbours molecule type 'SOL'

NOTE 3 [file topoloji.top, line 7198]:
  System has non-zero total charge: -2.000000
  Total charge should normally be an integer. See
  http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
  for discussion on how close it should be to an integer.
  



NOTE 4 [file topoloji.top, line 7198]:
  The bond in molecule-type Protein_chain_A between atoms 114 OG and 115 HG
  has an estimated oscillational period of 9.0e-03 ps, which is less than
  10 times the time step of 1.0e-03 ps.
  Maybe you forgot to change the constraints mdp option.

Number of degrees of freedom in T-Coupling group rest is 31806.00

NOTE 5 [file 4ion.mdp]:
  You are using a plain Coulomb cut-off, which might produce artifacts.
  You might want to consider using PME electrostatics.



There were 5 notes

gcq#77: "Every Sperm is Sacred" (Monty Python)

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

                   GROningen MAchine for Chemical Simulation

                            :-)  VERSION 4.6.5  (-:

        Contributions from Mark Abraham, Emile Apol, Rossen Apostolov, 
           Herman J.C. Berendsen, Aldert van Buuren, Pär Bjelkmar,  
     Rudi van Drunen, Anton Feenstra, Gerrit Groenhof, Christoph Junghans, 
        Peter Kasson, Carsten Kutzner, 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-2012,2013, 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 Lesser General Public License
        as published by the Free Software Foundation; either version 2.1
             of the License, or (at your option) any later version.

                                :-)  genion  (-:

Reading file 4ion.tpr, VERSION 4.6.5 (single precision)
Reading file 4ion.tpr, VERSION 4.6.5 (single precision)
Group     0 (         System) has 15527 elements
Group     1 (        Protein) has   755 elements
Group     2 (      Protein-H) has   389 elements
Group     3 (        C-alpha) has    50 elements
Group     4 (       Backbone) has   150 elements
Group     5 (      MainChain) has   202 elements
Group     6 (   MainChain+Cb) has   248 elements
Group     7 (    MainChain+H) has   255 elements
Group     8 (      SideChain) has   500 elements
Group     9 (    SideChain-H) has   187 elements
Group    10 (    Prot-Masses) has   755 elements
Group    11 (    non-Protein) has 14772 elements
Group    12 (          Water) has 14772 elements
Group    13 (            SOL) has 14772 elements
Group    14 (      non-Water) has   755 elements
Select a group: Number of (3-atomic) solvent molecules: 4924
Replacing solvent molecule 577 (atom 2486) with NA
Replacing solvent molecule 3730 (atom 11945) with NA
Replacing solvent molecule 2403 (atom 7964) with NA
Replacing solvent molecule 4166 (atom 13253) with NA
Replacing solvent molecule 1639 (atom 5672) with NA
Replacing solvent molecule 788 (atom 3119) with NA
Replacing solvent molecule 1901 (atom 6458) with NA
Replacing solvent molecule 4363 (atom 13844) with NA
Replacing solvent molecule 1939 (atom 6572) with NA
Replacing solvent molecule 3228 (atom 10439) with NA
Replacing solvent molecule 2633 (atom 8654) with NA
Replacing solvent molecule 3497 (atom 11246) with NA
Replacing solvent molecule 1998 (atom 6749) with NA
Replacing solvent molecule 1929 (atom 6542) with NA
Replacing solvent molecule 1469 (atom 5162) with NA
Replacing solvent molecule 402 (atom 1961) with NA
Replacing solvent molecule 1383 (atom 4904) with CL
Replacing solvent molecule 2502 (atom 8261) with CL
Replacing solvent molecule 2020 (atom 6815) with CL
Replacing solvent molecule 699 (atom 2852) with CL
Replacing solvent molecule 942 (atom 3581) with CL
Replacing solvent molecule 977 (atom 3686) with CL
Replacing solvent molecule 1002 (atom 3761) with CL
Replacing solvent molecule 2806 (atom 9173) with CL
Replacing solvent molecule 3036 (atom 9863) with CL
Replacing solvent molecule 3960 (atom 12635) with CL
Replacing solvent molecule 557 (atom 2426) with CL
Replacing solvent molecule 257 (atom 1526) with CL
Replacing solvent molecule 2664 (atom 8747) with CL
Replacing solvent molecule 3460 (atom 11135) with CL


gcq#77: "Every Sperm is Sacred" (Monty Python)

In [ ]:
import nglview
view = nglview.show_file("cok_iyonize_protein.pdb")
view.add_surface("sol and not within 3 of protein",color="cyan",opacity=0.1)
view.add_ball_and_stick("NA",radiusSize=1) #mor renkli
view.add_ball_and_stick("CL",radiusSize=1) # yeşil renkli. 
view #görüntüle
In [30]:
%%html --isolated
<iframe src="Simulasyon_sistemi.html" allowTransparency="true" scrolling="no" width=100% height=600px> </iframe>

2.Simülasyon

Bu kısım için maalesef size "bende hazır yapılmışı var" diyeceğim. Ama size tabii ki bu simülasyonu nasıl elde ettiğimi anlatacağım. Ben bu rehberi genişletene kadar sizi İngilizce de olsa başka bir kaynağa yönlendirmek istiyorum. Justin Lemkul bu rehberleri GROMACS için hazırlamış ve düzenli olarak güncelliyor. Farklı bir sürü örnek bulabilirsiniz simülasyonlar ile neler yapılabileceğine dair. Dönelim bizim insüline...

Sağlıklı bir veri elde etmemiz için öncelikle sistemi biraz rahatlatmamız lazım. Biz proteini hazırlarken, suyu ve iyonları eklerken hiç sormadık acaba keyfiniz yerinde mi, komşu atomlarınız sizi çok rahatsız ediyor mu diye. Bu aşamada sistemdeki her bir bileşenin hem kendi içindeki hem de çevresindeki kavgalarını çözeceğiz. Bir diğer deyişle, enerjisini azaltacağız. Rahatlatma birbirini takip etmesi gereken iki ayrı basamak ile sağlanıyor.

  1. Enerji İndirgenmesi: Protein yapıları protein veri bankasına yüklenirken bizim kullandığımız kuvvet alanları ile hazırlanmadıkları için sahip oldukları enerji yüksek olabiliyor. Bu enerji, yazının başındaki grafikteki terimler ile hesaplanıyor. Sistemin enerjisi yüksek olduğunda hareketi de aşırı olacaktır. Biz zaten proteini hareket ettirmek istiyorduk diyeceksiniz ama bu aşırı hareket fiziki olmayan bir durum. Mesela bizim tanımladığımız hidrojen atomları (deneysel yöntemler hidrojen atomlarını göremezler) diğer atomlara gereğinden fazla yakın yerleştirilmiş olabilir. Bu nedenle bu aşamada her molekül bulunduğu yerde titreştiriliyorlar (aynen telefonun titreşimi gibi bu; oldukları yerde hafifçe hareket ediyorlar) ki bizim kuvvet alanına en uygun şekli alabilsinler.
  2. Çözelti Dengelenmesi: Biz bu su ve iyonları eklerken onları mümkün olduğunca üst üste gelmeyecek şekilde yerleştirdik. Bunu yapabilmenin en iyi yolu ise bir ağ dizgesi ya da ızgara(grid) şeklinde dizmektir. O nedenledir ki üstteki resimde bu moleküller sıra sıra dizili gibiler, özellikle su molekülleri proteini bir kılıf gibi sarmıyorlar. Bu fiziki bir durum değil. Çözeltinin protein ile daha sağlıklı bir etkileşime girmesi ve bütün boşlukları doldurması için kısa bir simülasyon yapılır. Bu kısa simülasyonda proteinin henüz hareket etmesini istemiyoruz. O nedenle ekstra bir güç uyguluyoruz. Bu simülasyon sadece cihazları kalibre etmek gibi bir şey o nedenle bu aşamada veri toplanması yapılmıyor.

Veri toplanmasını her molekülün serbest dolaşım hakkına sahip olduğu üretim simülasyonlarından yapıyoruz. Bu simülasyonda proteine uygulanan ekstra gücü kaldırıyoruz. İşte size aşağıda bu simülasyonu gösteriyorum. Pek bir şey olmadığını göreceksiniz, malum ancak 1-2 nanosaniye uzunlukta.Gülün gülün, bende gülüyorum, o bir iki nanosaniyeyi alabilmek için eskiden kaç yıllar geçerdi. Bir deney tüpünde var 1 milyon protein bense gitmiş tek bir tane proteini birkaç nanosaniye kıpraştıracağım diye debeleniyorum.

Öyle bir delilik hali işte bilim. Ama yine de güzel değil mi, insülün ile tanışıp onu sulu bir ortama yerleştirdiniz. İleride simülasyon kısmını detaylandıracağım bir rehberde de koşullar (misal sıcaklık) ile oynayarak neler yapabileciğimizi öğreniriz. Bu alana ilk başladığımda ben böyle film izler gibi izliyordum videoları, keşif yapıyordum. Sonra öğrendim ki meğer grafikler alınıyormuş, her şey matematik her şey istatistikmiş. Şimdi o üstteki resimlere ancak makale yazacağım zaman bakıyorum ben.

Neyse uzun bir rehber oldu ama umarım eğlenmişsinizdir. Bilim ve sanatla kalın! Hoşçakalın. (Aşağıda başka bilgilerde var, ben kaçtım ama siz kaçmayın devam edin.)

In [31]:
%%html --isolated
<iframe src="Simulasyon.html" allowTransparency="true" scrolling="no" width=100% height=550px> </iframe>

Gerekli programlar ve dosyalar

Programlama öğrenmek istiyorsanız Python öğrenmek güzel bir adım. İsterseniz Python'ı tek başına indirip kurabilir isterseniz aşağıdaki program üzerinden kurabilirsiniz.

Kendi bilgisayarınızda benim yaptığım gibi bu rehberi çalıştırmak, başka rehberler oluşturmak için Jupyter Notebook kurabilirsiniz. Bu rehberdeki tüm bileşenler Python ve Jupiter Notebook da dahil Anaconda'da mevcut. Tavsiyem Anaconda kurmanız. Komut istemi, kurma etme vs. de yeterli olmama rağmen bütün bu programlar ve gereklilikleri ile tek tek ilgilenmek çok vakit kaybettirici bir iş. O nedenle ben de dahil bir çok kişi Anaconda gibi gibi uygulamalar kullanıyor. Anacondayı kurduktan sonra ek olarak yüklemeniz gereken paketler ise şöyle:

conda install -c bioconda gromacs

conda install -c bioconda nglview

Biraz yorucu olabilir bu işe başlamak, korkutabilir. Bir çok işlemi komut isteminden yapmak daha kolay ama ilk adımı atıp başlaması zor olabilir. Hiç çekinmeyin! Tavsiyem şimdiden elinizi bulaştırıp bu işe komut isteminden girişmeniz.

Benim hazırladığım bu rehberi ve tüm dosyaları github.com/ozyo/MDRehber adresinden indirebilirsiniz. Rehber henüz emekleme aşamadında ama vaktim oldukça güncelleyeceğim. Yorumlarınız veya önerileriniz varsa bana ozgeyoluk_at_proteinart.net adresinden veya EpiSTEM Türkiye üzerinden ulaşabilirsiniz.