Tag: พลวัตเชิงโมเลกุล

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

ต่อครับ

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

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

%d bloggers like this: