มาลองทำ Molecular Dynamics Simulation ของโปรตีนกันครับ 6

ต่อจากครั้งแล้ว http://www.sakngoi.com/?p=473
กลับมาดูระบบที่ผมพูดถึงในตอนแรก (HIV-RT + ยา) เนื่องจากว่าผมสงสัยว่าที่ปลายของสาย DNA บริเวณ active site ที่มีเจ้า inhibitor มาจับนี้จะมีการสั่นหรือเปลี่ยนแปลงอย่างไร ก็เลยออกแบบการทดลอง โดยลองเปลี่ยนบริเวณนั้นตามตารางด้านล่างนี้ครับ

ระบบ 1,2,5 และ 6 นั้น เป็นสถานการณ์ก่อนที่ inhibitors หรือ nucleotides จะเข้าจับกับสาย DNA ครับ ส่วนระบบ 3,4,7 และ 8 เป็นสถานการณ์หลังจากที่ inhibitors หรือ nucleotides จับกับสาย DNA แล้ว

อันนี้แสดงบริเวณที่สนใจครับ

ส่วนนี้แสดง pathway ของยาที่สนใจ (Adefovir/Tenofovir) ครับ

ระบบที่จะทำ MD ทั้ง 9 ระบบนี้ในแต่ละระบบจะมีอะตอมประมาณ 130,000  อะตอมครับ ส่วนการเตรียมโครงสร้างของแต่ล่ะระบบผมก็ทำคร่าวๆตามนี้ครับ โหลดโครงสร้างจาก PDB.org ที่ต้องการซึ่งก็คือ 1T05 เจ้าโครงสร้างนี้จะมีสาย DNA ที่มี tenofovir diphosphate (PMPApp) ติดมาด้วยแล้ว จากนั้นก็ใช้โปรแกรม LEAP ซึ่งเป็นหนึ่งในหลายโปรแกรมของ AmberTools (http://ambermd.org/#AmberTools) กับ CHIMERA (http://www.cgl.ucsf.edu/chimera/) มาช่วยในการเตรียมโครงสร้างไฟล์ (พวกtopology กับ initial coordinates)  ที่ต้องการของทั้ง 9 ระบบครับ (ถ้าไม่ขี้เกียจผมจะอธิบายให้ฟังครับเตรียมอย่างไรโดยละเอียดครับ :))

ภาพบริเวณ active site ของระบบที่ 1

ภาพบริเวณ active site ของระบบที่ 8

 

ขั้นตอนในการรัน MD ของระบบทั้งหมดนี้ก็คล้ายๆกันครับ เริ่มจากทุกระบบจะถูกทำ energy minimization โดยเทคนิค smooth Particel-Mesh Ewald (PME) ถูกใช้ด้วยสำหรับแรง electrostatic ส่วนระยะแรง van de Waals ถูกเซ็ตไว้ที่ 9 angstrom SHAKE ก็ถูกใช้ด้วยสำหรับพันธะของที่อีกด้านหนึ่งเป็น hydrogen  Langevin dynamics ก็ถูกประยุกต์ใช้ในการการควบคุมอุณหภูมิของระบบให้อยู่ที่ 310 K ความดันก็ถูกควบคุมให้เท่ากับ 1 atm ด้วยLangevin piston method และแต่ละระบบจะรันนาน 2.5 ns

ต่อครับ

ติดตั้ง driver ของ nvidia quadro บน Scientific Linux 6.1

พอดีว่าได้เครื่อง Dell Precision T1600 มาทำงานด้าน Molecular Simulation
เครื่องมีการ์ดจอที่สนับสนุน CUDA มาด้วยคือ Nvidia Quadro 2000
หลังจากติดตั้ง Scientifc Linux (SL6.1) เสร็จแล้ว ก็ทำการเพิ่ม Repository จาก elrepo.org เข้าไป
ตามนี้
rpm –import http://elrepo.org/RPM-GPG-KEY-elrepo.org

rpm -Uvh http://elrepo.org/elrepo-release-6-4.el6.elrepo.noarch.rpm

จากนั้นก็ติดตั้ง package ชื่อ kmod-nvidia

yum –disablerepo=\* –enablerepo=elrepo install kmod-nvidia 

แล้วก็ reboot เครื่องใหม่ เป็นอันเสร็จพิธี

ลอง render โครงสร้างของโปรตีนโดยใช้ VMD ปรากฏว่าเร็วมากอย่างเห็นได้ชัด

 

มาลองทำ Molecular Dynamics Simulation ของโปรตีนกันครับ 5

ต่อจากครั้งที่แล้ว (http://www.sakngoi.com/?p=450)

จากคุณสมบัติต่างๆทางเทอร์โมไดนามิคของระบบที่ต้องคำนึงถึงแล้ว ลองมาดูเรื่องของขั้นตอนของการทำ MD ว่าเราจะต้องมีอะไรคิดถึงหรือทำอะไรอีก  ในขั้นนี้ผมอยากจะให้มองว่าตัวโปรแกรมคอมพิวเตอร์สำหรับทำ MD อย่าง NAMD,GROMACS,CHARMM หรือ AMBER นั้นเป็นเหมือนกล่องดำที่จะคำนวณตำแหน่งของอะตอมต่างๆที่ประกอบเป็นโมเลกุลหรือระบบที่เราสนใจทุกๆช่วงเวลาที่เราสนใจ (รวมถึงเรื่องพลังงาน และอื่นๆอีก) เหมือนภาพถ่าย โดยที่เราจะป้อนค่าต่างๆ ของระบบที่ต้องการเข้าไป

การทำ MD เราอาจแบ่งออกได้เป็นขั้นตอนต่างๆตามลำดับดังนี้ครับ

Initial Input

ขั้นตอนนี้เป็นขั้นตอนของการเตรียมระบบ,ไฟล์ และ ค่าพารามีเตอร์ต่างๆสำหรับในการรันMD เช่น โครงสร้างโปรตีนที่เอามาจาก PDB นั้นมันจะไม่มีตำแหน่งของไฮโรเจนมาด้วยครับ เราก็ต้องเอามาใส่ไฮโดรเจนเองซึ่งการใส่ไฮโดรเจนเข้าไปก็มีโปรแกรมช่วยครับ (เดี๋ยวพูดถึงทีหลัง 🙂 ว่ามันใส่ยังไง ถ้าไม่ขี้เกียจก่อนนะครับ) , การใส่น้ำหรือ solvate โปรตีน หรือการในประจุของธาตุต่างๆเพื่อให้ระบบของเราที่มีน้ำ+โปรตีน มีประจุเป็นศูนย์ (neutralize)

 Energy minimization

ในขั้นตอนนี้เป็นการขยับอะตอมที่อาจจะถูกวางไว้ใกล้กันเกินไปจากขั้นตอนของการใส่อะตอมไฮโดนเจนหรือโมเลกุลของน้ำออกจากกันหน่อย เพราะถ้าหากเราเริ่มรันMDโดยที่มันมีอะตอมที่อยู่ใกล้กันเกินไปก็จะทำให้อะตอมอาจถูกผลักออกจากกันได้ ซึ่งอาจเป็นผลให้ระบบของเราพังได้ 🙂 ในขั้นตอนนี้ ถ้าหากเรา plot graph ระหว่างพลังงานของระบบกับจำนวน step จะเห็นว่าพลังงานมันค่อยๆลดลง

Heating

หลังจากที่เราทำ energy minimization แล้วเราก็ต้องทำการเพิ่มอุณหภูมิของระบบให้มีอุณหภูมิตามที่เราต้องการ โดยอาศัยเทคนิคตามที่ได้พูดถึงไปแล้ว

Equilibration

พอทำให้ระบบมีอุณหภูมิตามต้องการแล้วเราก็ต้องทำการควบคุมปริมาณอื่นๆ เช่น ความดัน หรือ ปริมาตรของระบบ ให้มีค่าตามที่เราต้องการด้วย ซึ่งในขั้นตอนนี้เราจะทำการรัน MD อยู่สักระยะให้ระบบมีค่าต่างๆ ตามที่เราต้องการ

MD Simulation

หลังจากที่ระบบพร้อมแล้วเราก็ทำการรัน MD ตามต้องการครับ

Data analysis

เสร็จแล้วก็ถึงเวลาวิเคราะห์ผลครับ

ต่อครับ http://www.sakngoi.com/?p=559

มาลองทำ Molecular Dynamics Simulation ของโปรตีนกันครับ 4

ต่อจากครั้งที่แล้ว http://www.sakngoi.com/?p=275
มาถึงเรื่องการควบคุมความดันของระบบครับ 🙂  ในโปรแกรม NAMD การควบคุมความดันให้คงที่มีหลายวิธีครับเช่น

๑) Berendsen pressure bath coupling
๒) Extended system method
๓) Langevin Piston Algorithm
๔) Nosé-Hoover constant pressure method

[latexpage]
วิธีที่ ๑
ไอเดียก็คล้ายกับวิธี weak coupling with a heat bath สำหรับการควบคุมอุณหภูมิครับ (แต่เป็น pressure bath แทน) โดยที่อัตราการเปลี่ยนของความดันคือ

\begin{eqnarray*}
\frac{dP(t)}{dt} &=& \frac{1}{\tau_P}\left(P_{bath}-P(t)\right)\nonumber
\end{eqnarray}
โดย $\tau_P$ คือ coupling constant, $P_{bath}$ คือความดันของ bath โดยที่ปริมาตรของ simulation box จะถูกคูณด้วย $\lambda$
\begin{eqnarray*}
\lambda&=&1-\kappa\frac{\delta t}{\tau_P}\left(P-P_{bath}\right)\nonumber
\end{eqnarray}

ใน NAMD ก็ต้องเซ็ต BerendsenPressure เป็น on ก่อน ส่วน BerendsenPressureTarget คือความดันที่ต้องการและ BerendsenPressureRelaxationTime ก็คือ Relaxation time

วิธีที่๒-๔ ไปหาอ่านเอานะครับ 😛
ใน NAMD จะมีวิธีควบคุมความดันที่เป็นการรวมเอาวิธีที่๓กับ๔เข้าด้วยกัน ซึ่งถ้าต้องการใช้วิธีที่ว่านี้ก็ต้องเซ็ต LangevinPiston เป็น on โดย LangevinPistonTarget กับ LangevinPistonTemp ก็เป็นความดันกับอุณหภูมิที่ต้องการตามลำดับครับ

ในการทำ MD นั้นนอกเหนือจากการควบคุมอุณหภูมิกับความดันแล้วยังต้องคำนึงถึงเรื่องอื่นๆอีกเช่น
1) ระบบที่สนใจจริงๆในธรรมชาติมันมีหลายอนุภาคหรือดูเหมือนไม่มีขอบหรือที่สิ้นสุด เราก็สามารถที่จะลดขนาดของการคำนวณระบบใหญ่ๆนี้ได้โดยใช้เทคนิค periodic boundary conditions โดยมันก็เป็นการ copy ระบบที่เราย่อมาแล้วนำไปวางต่อทุกทิศทาง เช่นถ้าระบบของเราเป็นกล่อง(cell) เทคนิคนี้ก็เหมือนเอาตัวcopyกล่องนี้ไปวางรอบทั้งหกด้าน ใน NAMD การเซ็ตตำแหน่งหรือขนาดของกล่องก็คือ cellOrigin กับ cellBasisVector1(2,3)

2) electrostatic potential ลองไปดูเพิ่มเติมเรื่อง Ewald method ซึ่งเป็นวิธีที่ช่วยคิด potentialนี้ครับ ใน NAMD ถ้าจะใช้ Ewald method ก็คือ PME (Particle-Mesh Ewald) มันก็จะมีkeywords ที่เกี่ยวข้องเช่น PMETolerance, PMEGridSizeX(Y,Z) และ PMEInterpOrder

3) บางที่เราอาจจะลดเวลาที่ใช้ในการคำนวณได้โดยการใช้เทคนิค SHAKE เพื่อfixขนาดของพันธะที่อีกข้างเป็นอะตอมเบาๆเช่นไฮโดรเจนไม่ต้องสั่นมาก ถ้าจะใช้เทคนิคนี้ใน NAMD ก็ลองดูที่ rigidBonds, rigidTolerance, rigidInterations

4) ในกรณีที่ระบบของเรามีน้ำ ก็อาจจะต้องดูด้วยว่าจะใช้โมเดลของน้ำอันไหน (http://en.wikipedia.org/wiki/Water_model) ในกรณีของเรานี้ใช้ TIP3P ครับ

ต่อ ครับ 😛  http://www.sakngoi.com/?p=473

มาลองทำ Molecular Dynamics Simulation ของโปรตีนกันครับ 3

ต่อจากตอนที่แล้ว http://www.sakngoi.com/?p=238

พอเห็นภาพคร่าวๆของสมการการเคลื่อนที่ที่เราใช้แล้วนะครับ ตอนนี้เราจะมาลองดูเรื่องคุณสมบัติทางเทอโมไดนามิค (อุณหภูมิ, ความดัน) ของระบบที่เรากำลังสนใจครับ เพราะในการคำนวณจริงๆ เราต้องคิดถึงว่าระบบของเรานั้นอยู่ในสภาวะแบบไหน อุณหภูมิคงที่หรือเปล่า ความดันล่ะคงที่ด้วยหรือเปล่า จำนวนของอะตอมคงที่หรือมีการหายไปด้วยอีกหรือเปล่า เพื่อที่ว่าจะได้ดูสมจริงมากที่สุดเมื่อเอาไปเปรียบเทียบกับของจริงในธรรมชาติ

ในโปรแกรม NAMD ที่ผมใช้นี้ การควบคุมในอุณหภูมิคงของระบบ ก็มีอยู่สามวิธีครับ คือ
1) rescaling of velocities
2) weak coupling with a heat bath
3) Langevin dynamics

ในวิธีที่ 1
ก็คือค่อยๆเปลี่ยนอุณหภูมิไปเรื่อยๆที่ล่ะนิด ด้วยการคูณอุณหภูมิของระบบด้วยค่าคงที่ (จริงๆแล้วคูณกับความเร็ว)ครับ
อุณภูมิมันสัมพันธ์กับพลังงานจลน์ของระบบ ดูสมการด้านล่างครับ

\Delta T = \frac{1}{3}\sum_{i=1}^{N}\frac{m_i (\lambda v_i)^2}{Nk_b}

-\frac{1}{3}\sum_{i=1}^{N}\frac{m_i\lambda  v_i^2}{Nk_b}

\Delta T = \left( \lambda^2 - 1\right) T(t)

\lambda  = \sqrt{T_{req}/T(t)}

ใน NAMD จะมีพารามีเตอร์คือ rescaleTemp – อุณหภูมิที่ต้องการ กับ rescaleFreq -ความถี่ที่จะคูณเจ้าค่าคงที่นี้ประมาณว่าจะคูณมันกี่step

วิธีที่ 2

วิธีนี่ก็คล้ายกับวิธีแรก แต่จะมองเหมือนว่าระบบของเราถูกแช่อยู่ในอ่างของความร้อนที่มีอุณหภูมิที่เราต้องการ ระบบของเรากับอ่างนี้ก็จะมีการแลกเปลี่ยนความร้อนกัน โดยที่อัตราการแลกเปลี่ยนความร้อนหรืออุณหภูมิสามารถเขียนได้ว่า

\frac{dT(t)}{dt}=\frac{1}{\tau}\left(T_{bath}-T(t)\right)

หรือในแต่ล่ะstep อุณหภูมิจะต่างกัน

\Delta T  = \frac{\delta t}{\tau}\left(T_{bath}-T(t)\right)

\lambda^2 = 1+ \frac{\delta t}{\tau} \left( \frac{T_{bath}}{T(t)}-1 \right)

วิธีนี้ NAMD มีพารามีเตอร์  tCoupleTemp- ค่าของอุณหภูมิที่ต้องการ (T_{bath}) แต่ต้องมีการเรียก tCouple = on ก่อนนะครับถึงจะใช้ได้

วิธีที่ 3

การควบคุมอุณหภูมิด้วยวิธีที่เรียกว่า Langevin dynamics เป็นอีกวิธีที่นิยมใช้กัน
วิธีนี้จะไม่ได้แก้ไขความเร็วของอะตอมหรือโมเลกุลโดยตรงเหมือนสองวิธีแรก แต่จะมีการเพิ่มแรงนอกเหนือจากแรงที่คำนวณจาก
force field ที่กระทำกับอะตอมหรือโมเลกุลที่สนใจแบบสุ่มเข้าไป (R_i(t))
แล้วลบออกด้วยแรงเสียดทานที่มีค่าสัมประสิทธิ์เสียดทาน หรือ damp constant \eta แรงที่เพิ่มเข้านี้จะเป็นเหมือนการเพิ่มพลังงานหรืออุณหภูมิให้ระบบในขณะที่แรงเสียดทานจะเหมือนการเอาพลังงานที่เพิ่มเข้าไปนี้ออก

m_i\frac{dr_i^2}{dt} = F_{intra}+R_i(t)-\eta m_i\frac{dr_i}{dt}

โดยที่แรงสุ่ม (R_i(t)) จะมาจาก Gaussian distribution ครับ ซึ่งจะมีคุณสมบัติตามนี้

R_i(t)= 0

R_i(t)R_i(t) = 2k_b \eta T_o \deta(t)

โดย T_0 คืออุณหภูมิที่ต้องการครับ 🙂

ถ้าใช้ใน NAMD จะต้องตั้งค่า langevin เป็น on และ
langevinTemp เป็นอุณหภูมิที่ต้องการ ส่วนค่าของ \eta ก็คือ langevinDamping ครับ ลองเล่นดู

ต่อตอนที่ 4