Tag: GROMACS

มาลองทำ 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

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

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

Molecular dynamics (MD) ก็เป็นเทคนิคที่ฮิตฮอตสำคัญอันหนึ่งที่นักวิทยาศาสตร์ใช้ในการศึกษาเกี่ยวกับการเคลื่อนไหวหรือสั่นไหวของโมเลกุลต่างๆ โดยเค้าจะมองว่าพันธะที่คอยเชื่อมต่ออะตอมเข้าด้วยกันเป็นโมเลกุลนี้เป็นเหมือนสปริงที่เชื่อมลูกบอลเข้าไว้ด้วยกัน  ซึ่งพลังงานศักย์ทั้งหมดของระบบก็จะมาจากผลรวมทั้งหมดของ interaction ระหว่างลูกบอลที่เชื่อมกันผ่านสปริง(bond) และระหว่างสปริงด้วยกัน(non-bonded) ถ้าใครที่เรียนมาก็คงจะจำได้ว่าแรงของปริงนี้เป็นแรงอนุรักษ์ซึ่งเราจะเขียนได้ว่าแรงทั้งหมดที่กระทำกับลูกบอล(อะตอม) i นี้เขียนได้ว่า
F_i=-\nabla E

โดย E ก็คือพลังงานศักย์ ซึ่งเราเราสามารถเขียนได้ว่า E_{total}=E_{bond}+E_{non-bonded} โดยเทอมที่ไม่เป็นพันธะสามารถเขียนได้ว่าเป็นผลรวมศักย์จากแรงวันเดอร์วาลส์กับแรงจากประจุไฟฟ้า E_{non-bond}=E_{vdW}+E_{elec} โดย E_{elec}=\sum_{i=1}^{N} \sum_{j=i+1}^{N}\frac{q_iq_j}{4\pi\varepsilon_0 r_{ij}} และ E_{vdW}= \sum_{i=1}^{N} \sum_{j=i+1}^{N}4\varepsilon\left(\frac{\sigma_{ij}^{12}}{r_{ij}^{12}}-\frac{\sigma_{ij}^{6}}{r_{ij}^{6}}\right)

ส่วนพลังงานศักย์ของแรงจากพันธะก็สามารถเขียนได้ว่าเป็นรวมจากรูปแบบของพันธะที่เป็นไปได้ทั้งหมด (ดูรูป เดี๋ยวมาใส่ให้ครับ)

E_{bond}=E_{stretch}+E_{angle}+E_{out-of-plane}+E_{str-str}+...
หรือ
E_{bond}=\frac{1}{2}\sum_\text{angles}k_a\left(\theta-\theta_0\right)^2+\frac{1}{2}\sum_\text{bonds}k_b\left(l-l_0\right)^2+

\sum_\text{bonds}k_d\left(1 + \cos\left(n\phi-\phi_0\right)\right)+\sum_{\text{improper}}k_i\left(\omega-\omega_0\right)^2
สมการที่เขียนมานี่ก็เป็นเพียงตัวอย่างเท่านั้นครับซึ่งมันจะแตกต่างกันไป ขึ้นกับว่าเราเป็นสาวกของสำนักไหน 🙂 แต่ล่ะสำนักก็จะมีรูปแบบของฟังก์ชั่นและค่า constants ต่างๆที่ใช้แตกต่างกัน (เราเรียก พวกฟังก์ชั่นนี้ว่า force field)
force field ที่เป็นที่นิยมกันก็ได้แก่

AMBER (Assisted Model Building with Energy Refinement, http://ambermd.org/)

OPLS, OPLS-AA (Optimized Potentials for Liquid Simulations,http://zarbi.chem.yale.edu/)

CHARMM (Chemistry at HARvard Macromolecular Mechanic, http://mackerell.umaryland.edu/MacKerell_Lab.html)

GROMOS (GROningen Molecular Simulation, http://www.gromos.net/)

โดยปกติเค้าจะไม่เอาแต่ล่ะ force field มาใช้ปนกันครับ

มาต่อที่เรื่องของสมการนิวตันที่เราได้ครับ F_i=m_i\frac{d^2r_i}{dt^2}=-\nabla E  ซึ่งเป็นสมการที่แสดงความสัมพันธ์ของพลังงานศักย์กับระยะทางหรือตำแหน่งของอะตอมที่เราสนใจ เราสามารถแก้สมการนี้โดยอาศัยวิธีที่เรียกว่า Verlet หรือ Leap-fr0g ใครที่เรียนเรื่อง numerical analysis มาคงเข้าใจนะครับ เพราะมันก็เป็นวิธีการอินทีเกรตทำธรรมดา แตกกันก็เพียงว่าจะเอากี่เทอม เริ่มที่ตรงไหนก่อนเท่านั้นเอง

ใครที่สนใจเรื่องการเขียนโปรแกรมเกี่ยวกับ MD นี้ลองดูที่นี่ครับ http://www.sakngoi.com/?p=64

ต่อตอนที่ 3

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

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

โปรแกรมที่ผมใช้ในการคำนวณคือ NAMD ส่วนที่เอาไว้ใช้ดูโครงสร้างสวยๆของโปรตีนคือ VMD

ทั้งสองโปรแกรมนี้ฟรีครับ

โปรตีนที่ผมพูดถึงนี้ก็คือ HIV 1-RT (reverse transcriptase) ครับ เราจะมาลองดูครับว่าเจ้าสาย DNA หรือโครงสร้างของโปรตีนแถวๆ Active site ของมันจะมีการเปลี่ยนแปลงไปยังไงบ้างถ้าเราใส่ตัว Inhibitor หรือยาอย่าง adefovir หรือ tenofovir เข้าไป โครงสร้างของโปรตีนที่ผมใช้คือ 1T05 ซึ่งสามารถdownload ได้จากเวบ www.pdb.org

รูปด้านล่างนี้แสดงลำดับของนิวคลีโอไทด์กับตำแหน่งของกรดอะมิโนจากโปรตีนที่ผมว่าจะเอาไว้ดูว่ามันเจ้าสายDNA กับ HIV-RT มันมีอะไรกันหรือเปล่า จริงๆแล้วมันก็คือสะพานเกลือครับ 🙂 เรียกตรงๆเลย ผมสนใจว่าเจ้าฟอสเฟตจากสาย DNA กับพวกขั้วบวกจากกรดอะมิโนอย่างไลซีนกับอาร์จินินมันจะแนบแน่นขนาดไหนตลอดเวลาหรือเปล่า ซึ่งผมก็หาอยู่ตั้งนานก็เจอตัวที่ใกล้ที่กันมากที่สุดก็ตามรูปเลยครับ

ส่วนเจ้าสาย DNA เราก็จะดูว่า furanose ring มันมีการผับตัวหรือบิดยังไง (sugar puckering) ในตำแหน่งต่างๆ

แต่ก่อนจะไปดูว่าเราจะทำ Molecular dynamics ได้อย่างไรเราก็ควรจะเข้าใจก่อนสักนิดหนึ่งว่า Molecular dynamics คืออะไรซึ่งเดี๋ยวผมจะเล่าให้ฟังคร่าวๆนะครับส่วนใครที่สนรายละเอียดทางด้านทฤษฎีของ Molecular dynamics ผมแนะนำให้อ่านหนังสือตามนี้เลยครับ

The Art of Molecular Dynamics Simulation by D. C. Rapaport

Understanding Molecular Simulation, Second Edition: From Algorithms to Applications (Computational Science) by Daan Frenkel and Berend Smit

Molecular Modeling and Simulation by Tamar Schlick

Molecular Modelling: Principles and Applications (2nd Edition) by Andrew R. Leach

Computer Simulation of Liquids by M. P. Allen and D. J. Tildesley

ทำไมถึงต้องเป็นห้าเล่มนี้ ง่ายๆเลยครับเพราะผมมีแค่ห้าเล่มนี้เอง 🙂

ต่อตอน 2

การทำงานของโปรแกรมสำหรับ Molecular dynamics simulation

ใครที่สนใจว่าโปรแกรมที่เราใช้ศึกษาด้าน molecular dynamics simulation เช่น AMBER, GROMACS, CHARMM หรือ NAMD ทำงานอย่างไร ผมขอแนะนำให้ลองศึกษาจาก source codes ของ  Mindy ครับ http://www.ks.uiuc.edu/Development/MDTools/mindy มันเป็นเหมือนเวอร์ชั่นย่อๆ ของ NAMD ครับ ซึ่งมันทำงานหลักๆอย่างที่ NAMD ทำเลยครับ จะไม่มีก็ในส่วนของการคำนวณแบบขนาน

%d bloggers like this: