MỤC LỤC
MỞ ĐẦU 1
CHƯƠNG 1 - TỔNG QUAN 2
1.1. Đối tượng nghiên cứu .2
1.1.1. Acetylcholinesterase 2
1.1.2. Đặc điểm xúc tác 7
1.2. Phương pháp nghiên cứu .9
1.2.1. Protein docking 9
1.2.2. Phương pháp phiếm hàm mật độ .13
1.2.3. Cơ học phân tử .20
1.2.4. Kết hợp phương pháp cơ học lượng tử-cơ học phân tử .24
CHƯƠNG 2. NGUỒN DỮ LIỆU VÀ CÔNG CỤ TÍNH TOÁN 27
2.1. Nguồn dữ liệu 27
2.2. AutoDock 4.2 và AutoDockTools 1.5.4 .30
2.3. AutoDock Vina 1.1.1 33
2.4. Gaussian 03W và GaussView 3.0 33
CHƯƠNG 3. KẾT QUẢ VÀ THẢO LUẬN .36
3.1. Protein docking .36
3.2. Áp dụng phương pháp QM/MM đối với hệ phản ứng 43
3.2.1. Cấu trúc enzyme 43
3.2.2. Cơ chất trong hốc phản ứng ở trạng thái chưa liên kết .46
3.2.3. Cấu trúc phức enzyme-cơ chất .50
3.2.4. Cấu trúc sản phẩm .53
KẾT LUẬN 56
TÀI LIỆU THAM KHẢO 58
PHỤ LỤC .60
63 trang |
Chia sẻ: banmai | Lượt xem: 1897 | Lượt tải: 0
Bạn đang xem trước 20 trang tài liệu Luận văn Góp phần nghiên cứu cơ chế phản ứng ester bằng phương pháp tính lượng tử, để xem tài liệu hoàn chỉnh bạn click vào nút DOWNLOAD ở trên
0ρNeE (1.10)
Ta có thể phân tách phần năng lượng do tương tác hạt nhân -electron (có
dạng phụ thuộc vào hệ cụ thể) và phần còn lại có dạng chung độc lập với N, vị trí và
điện tích hạt nhân.
[ ]00 ρE = ( ) drVr Ne∫ 0ρ + [ ]0ρT + [ ]0ρeeE (1.11)
15
Định nghĩa phiếm hàm Hohenberg-Kohn
[ ]0ρHKF = [ ]0ρT + [ ]0ρeeE (1.12)
ta viết lại
[ ]00 ρE = ( ) rdVr Ne
∫ 0ρ + [ ]0ρHKF (1.13)
Nếu biết chính xác phiếm hàm [ ]0ρHKF chúng ta có thể có lời giải chính xác
cho phương trình Schrodinger. Vì phiếm hàm này có một dạng chung độc lập với
kích thước hệ nên rất thuận tiện khi áp dụng cho các hệ lớn. Nhưng chúng ta vẫn
chưa biết được dạng chính xác của cả [ ]0ρT và [ ]0ρeeE .
[ ]ρeeE = [ ]ρJ + [ ]ρK = 2
1 ( ) ( )
∫∫ 21
12
rdrd
r
rr ρρ + [ ]ρK (1.14)
Việc tìm dạng của các phiếm hàm [ ]ρT và [ ]ρK là thách thức lớn trong lý
thuyết phiếm hàm mật độ.
Định đề Hohenberg-Kohn thứ hai: nguyên lý biến phân
Chúng ta đã khẳng định rằng mật độ trạng thái cơ bản về nguyên tắc là đủ để
xác định đầy đủ tính chất của hệ. Nhưng, làm thế nào để chắc chắn rằng một mật độ
nào đó thực sự là mật độ ở trạng thái cơ bản.
Định đề hai thực chất có dạng của nguyên lý biến phân
0E ≤ [ ]ρ~E = [ ]ρ~T + [ ]ρ~NeE + [ ]ρ~eeE (1.15)
Có thể phát biểu như sau: với bất kì một hàm thử ( )rρ~ nào thỏa mãn các điều kiện
biên như ( )rρ~ ≥0, ( )∫ rdr
ρ~ =N tương ứng với một thế ngoài extV
~ nào đó, thì năng
lượng nhận được cũng không thể nhỏ hơn năng lượng trạng thái cơ bản 0E . Dấu
bằng chỉ nhận được nếu và chỉ nếu mật độ trong công thức đúng là mật độ trạng thái
cơ bản.
16
1.2.2.2. Phương trình Kohn-Sham
Hai định đề trên là cơ sở của phương pháp phiếm hàm mật độ nhưng chưa
chỉ ra được cách áp dụng vào hệ cụ thể vì chưa đưa ra được một dạng phiếm hàm
phù hợp liên hệ giữa năng lượng và mật độ electron. Năm 1965, Kohn và Sham đã
đề xuất một cách thức để xác lập phiếm hàm đã nói ở trên, trước hết là để tính động
năng với độ chính xác tương đối. Để xác định phần động năng này, Kohn, Sham
đưa vào khái niệm hệ quy chiếu không tương tác được xây dựng từ một tập hợp các
orbital là các hàm một electron. Phần sai số cùng với tương tác giữa các electron
[ ]ρK khá nhỏ sẽ được xác định bằng một phiếm hàm xấp xỉ.
Orbital và hệ quy chiếu không tương tác
Với mô hình hệ khí đồng nhất không tương tác, Thomas, Fermi đã xây dựng
trực tiếp các phiếm hàm động năng và tương tác electron-electron nhưng kết quả áp
dụng lại không phù hợp với thực tế, không mô tả được liên kết hóa học. Kohn và
Sham đã tìm một cách tiếp cận khác, đó là dựa vào các hàm sóng và liên hệ với
cách thức tiếp cận của Hatree-Fock.
Giả sử một hệ electron không tương tác, ta có thể viết toán tử Hamilton ở
dạng:
SHˆ = 2
1
− ∑∇
N
i
i
2 + ( )∑
N
i
iS rV
(1.16)
Và liên hệ với phương pháp Hartree-Fock thì định thức Slater là hàm sóng
chính xác. Đối với hệ tương tác thì thế SV là thế hiệu dụng địa phương tương tự như
thế hiệu dụng trong phương trình Hartree-Fock.
Ta đưa vào hàm sóng dưới dạng định thức Slater
SΦ =
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( )NNN
N
N
N
N
ϕϕϕ
ϕϕϕ
ϕϕϕ
...
............
2...22
1...11
!
1
21
21
21
(1.17)
17
Các orbital iϕ thỏa mãn phương trình
iii
KSf ϕεϕ =ˆ (1.18)
với KSfˆ là toán tử Kohn-Sham một electron
KSfˆ =
2
1
− 2∇ + SV (1.19)
Để áp dụng cho hệ thực là hệ tương tác ta phải tìm được thế hiệu dụng thích
hợp thỏa mãn điều kiện tổng các bình phương modun hàm sóng phải bằng mật độ
trạng thái cơ bản của hệ tương tác.
Phương trình Kohn-Sham
Kohn-Sham đề nghị dùng biểu thức dưới đây để nhận được động năng chính
xác của hệ không tương tác có cùng mật độ như hệ thực có tương tác
ST = 2
1
− ∑ ∇
N
i
ii ϕϕ
2 (1.20)
Tất nhiên, động năng của hệ không tương tác không thể bằng động năng của
hệ thực có tương tác dù chúng có chung một mật độ. Bao gồm phần sai khác này,
Kohn-Sham đưa vào số hạng năng lượng tương quan-trao đổi xcE .
( )[ ]rF ρ = ( )[ ]rTS
ρ + ( )[ ]rJ ρ + ( )[ ]rExc
ρ (1.21)
nghĩa là
( )[ ]rExc
ρ =( ( )[ ]rT ρ - ( )[ ]rTS
ρ )+( ( )[ ]rEee
ρ - ( )[ ]rJ ρ )= ( )[ ]rTer
ρ + ( )[ ]rK ρ (1.22)
Từ biểu thức trên ta thấy số hạng năng lượng tương quan trao đổi trong
phương pháp Kohn-Sham không đồng nhất với phần tương quan-trao đổi trong
phương pháp Hartree-Fock mà nó còn bao gồm cả một phần động năng không được
xác định chính xác.
18
Vấn đề đặt ra là làm thế nào để xác định duy nhất các orbital iϕ trong hệ
không tương tác hay nói cách khác là làm thế nào để định nghĩa được SV để có thể
nhận được định thức Slater tương ứng với mật độ điện tích đúng như hệ thực.
Ta viết lại biểu thức năng lượng có sự phụ vào orbital
( )[ ]rE ρ = ( )[ ]rTS
ρ + ( )[ ]rJ ρ + ( )[ ]rExc
ρ + ( )[ ]rENe
ρ
= ( )[ ]rTS
ρ +
2
1 ( ) ( )
∫∫ 21
12
rdrd
r
rr ρρ + ( )[ ]rExc
ρ + ( ) rdVr Ne
∫ 0ρ
=
2
1
− ∑ ∇
N
i
ii ϕϕ
2 +
2
1 ( ) ( )∑∑∫∫
N
i
N
j
ji rdrdrr
r 21
2
2
12
2
1
1 ϕϕ + ( )[ ]rExc
ρ -
( )∑∫∑
N
i
M
A
i
A
A rdr
r
Z 2ϕ (1.23)
Trong biểu thức trên, số hạng duy nhất không có dạng phụ thuộc rõ ràng là
xcE .
Tiếp theo áp dụng nguyên lý biến phân với điều kiện ràng buộc
ji ϕϕ = ijδ (1.24)
dẫn đến phương trình
( ) ( ) i
M
A A
A
xc r
ZrVrd
r
r
ϕ
ρ
−++∇− ∫ ∑12
12
22
2
1
= ( ) ieff rV ϕ
+∇−
2
2
1 = iε iϕ (1.25)
Xem xét các thành phần trong thế hiệu dụng Kohn-Sham
( )rVS
= ( )rVeff
= ( )∫ 2
12
2 rd
r
r ρ + ( )rVxc
-∑
M
A A
A
r
Z (1.26)
xcV = δρ
δ xcE (1.27)
Vì chúng ta chưa biết dạng phụ thuộc của xcE nên cũng chưa biết được dạng
của xcV . Nếu biết dạng chính xác xcE ( xcV ), phương trình Kohn-Sham sẽ cho ra trị
riêng chính xác nhưng cho đến nay, phiếm hàm này mới chỉ được đưa ra một cách
19
gần đúng. Và sự phát triển của lý thuyết phiếm hàm mật độ tập trung vào việc tìm ra
dạng tốt hơn của phiếm hàm tương quan-trao đổi.
Khi đã có dạng của phiếm hàm tương quan-trao đổi thì việc giải phương
trình Kohn-Sham được thực hiện bằng cách giải lặp tương tự như phương trình
Hartree-Fock.
Phiếm hàm tương quan-trao đổi phải được xác định chung cho tất cả các hệ.
Nhiều dạng phiếm hàm tương quan-trao đổi đã được đưa ra, việc xây dựng các dạng
gần đúng này dựa vào so sánh với thực nghiệm hoặc trên cơ sở so sánh với một
phương pháp hàm sóng mức cao. Thông thường phiếm hàm được tách thành 2 phần
riêng rẽ, phần tương quan và phần trao đổi.
[ ] [ ] [ ] ( ) ( )[ ] ( ) ( )[ ]drrrdrrrEEE cxcxxc ∫∫ +=+= ρερρερρρρ (1.28)
εx và εc được đưa vào biểu thức với ý nghĩa là mật độ năng lượng.
Thế tương quan-trao đổi tương ứng được xác định là đạo hàm của năng
lượng theo mật độ:
( ) [ ]( ) ( )[ ] ( )
( )
ρ
ε
ρρε
ρ
ρ
∂
∂
+=
∂
∂
=
rrr
r
ErV xcxcxcxc
(1.29)
Tương quan giữa các electron có spin song song khác với tương quan giữa
các electron có spin đối song, năng lượng trao đổi theo định nghĩa chỉ liên quan đến
các electron có cùng spin.
[ ] [ ] [ ]ββαα ρρρ xxx EEE += (1.30)
[ ] [ ] [ ] [ ]βααββββααα ρρρρρ ,cccc EEEE ++= (1.31)
Mật độ tổng là tổng của phần đóng góp của các electron α và β: ρ = ρα + ρβ.
Tuy nhiên, các phiếm hàm thường được viết theo độ phân cực hóa spin ζ và bán
kính của thể tích hiệu dụng chứa một electron rS.
20
βα
βα
ρρ
ρρζ
+
−
=
và
3
4
3
πρ
=Sr
(1.32)
Có nhiều dạng phiếm hàm tương quan -trao đổi đã được thiết lập, một số
dạng thường dùng trong hệ hóa học như BPW91, BLYP, B3LYP, B3PW91...
1.2.2.3. Phiếm hàm B3LYP
Phiếm hàm B3LYP là dạng phiếm hàm thường dùng trong khi nghiên cứu
các hệ hóa học theo phương pháp DFT. Đây là dạng kết hợp tuyến tính của phiếm
hàm trao đổi Hartree-Fock và các phiếm hàm tương quan, trao đổi dạng khác. Các
tham số xác định trọng số của mỗi phiếm hàm thành phần được xác định bằng cách
khớp với thực nghiệm hay các dữ kiện nhiệt hóa học được tính toán một cách chính
xác.
( ) ( ) ( )LDAcGGAccLDAxGGAxxLDAxHFxLDAxcLYPBxc EEaEEaEEaEE −+−+−+= 03 (1.33)
với 0a = 0.20, xa = 0.72, ca = 0.81. Các tham số này được xác định bằng cách
khớp các giá trị dự đoán với một tập hợp các giá trị năng lượng nguyên tử hóa, thế
ion hóa, ái lực proton, các giá trị năng lượng tổng của nguyên tử thể hiện trong các
nghiên cứu của Becke, và Lee-Yang-Parr [4, 9].
1.2.3. Cơ học phân tử
Trong các phương pháp trường lực, năng lượng điện tử với một cấu hình hạt
nhân cho trước được tính bằng cách viết lại Ee dưới dạng một hàm tham số của các
tọa độ hạt nhân. Các tham số đưa vào hàm này được lấy phù hợp với thực nghiệm
hay từ các phương pháp tính toán mức cao hơn. Phân tử được mô hình như hệ các
quả cầu và lò xo, gồm các nguyên tử được giữ với nhau bằng các liên kết. Các
nguyên tử được xử lí như trong cơ học cổ điển theo định luật II Newton.
Năng lượng trường lực được phân tách thành các số hạng mô tả năng lượng
cần thiết để làm biến dạng phân tử theo những kiểu riêng khác nhau.
crosselvdwtorsbendstrFF EEEEEEE +++++= (1.34)
21
Trong đó,
Estr : năng lượng cần để làm thay đổi độ dài liên kết giữa 2 nguyên tử
Ebond : năng lượng cần để làm thay đổi góc liên kết
Etors : năng lượng xoắn, cần để quay quanh một liên kết
Evdw , Eel : năng lượng tương tác nguyên tử-nguyên tử không liên kết
Ecross : mô tả ảnh hưởng qua lại giữa ba số hạng đầu tiên.
Để tìm cấu hình bền của phân tử tương ứng với cực tiểu trên bề mặt thế
năng, ta tiến hành cực tiểu hóa EFF theo các tọa độ hạt nhân.
1.2.3.1. Năng lượng thay đổi độ dài liên kết Estr
Ta khai triển năng lượng làm thay đổi độ dài liên kết giữa 2 nguyên tử A và
B theo chuỗi Taylor đến bậc 2 (đây là dạng đơn giản nhất).
( ) ( ) ( ) ( )202
2
00 2
10 RR
dR
EdRR
dR
dEERREstr −+−+=−
(1.35)
R0 là độ dài liên kết tự nhiên hay cân bằng giữa 2 nguyên tử A và B, điểm không
được xác định tại R0. Khi khai triển gần giá trị cân bằng, số hạng bậc nhất bằng
không, Estr có dạng đơn giản
( ) ( ) ( )2200 RkRRkRREstr ∆=−=− (1.36)
với k là hằng số lực của liên kết giữa A và B. Biểu thức có dạng mô tả dao động tử
điều hòa và với mỗi liên kết cần 2 tham số là k và R0. Để chính xác hơn cần thêm
các số hạng bậc cao hơn vào biểu thức khai triển Taylor.
( ) ( ) ( ) ( ) ...443322 +∆+∆+∆=∆ RkRkRkREstr (1.37)
Khi đó số tham số cần đưa vào cũng tăng lên. Ngoài ra còn có thể đưa vào các dạng
hàm khác để hiệu chỉnh.
1.2.3.2. Năng lượng làm thay đổi góc liên kết Ebend
22
Tương tự như Estr, ta cũng có thể viết Ebend ở dạng khai triển Taylor.
( )
2
00 )( θθθθ −=− kEbend (1.38)
θ0 là góc liên kết tự nhiên giữa 3 nguyên tử A -B-C. Cũng giống như E str, đây là
dạng đơn giản nhất của Ebend, để chính xác hơn có thể đưa thêm vào các số hạng
khai triển bậc cao hơn và hiệu chỉnh cho phù hợp với dữ kiện thực nghiệm hoặc các
phương pháp tính toán lượng tử mức cao.
1.2.3.3. Năng lượng xoắn Etors
Trong chuỗi 4 nguyên tử liên kết A-B-C-D, xét năng lượng làm quay quanh
liên kết B – C. Để đảm bảo tính tuần hoàn của phép quay, ta sử dụng khai triển
Fourier cho Etors.
( ) ( )∑
=
=
1
cos
n
ntors nVE ωω
(1.39)
Phụ thuộc vào tính đối xứng mà một vài hằng số Vn có thể bằng 0. Đối với
hệ phân tử hữu cơ biểu thức thông dụng cho năng lượng xoắn được viết ở dạng:
( ) ( )[ ] ( )[ ] ( )[ ]ωωωω 3cos12cos1cos1 321221121 ++−++= VVVEtors (1.40)
1.2.3.4. Năng lượng Van der Waals Evdw
Năng lượng Van der Waals mô tả tương tác đẩy hay hút giữa các nguyên tử
không liên kết trực tiếp với nhau và không tính đến phần tĩnh điện. Nếu khoảng
cách giữa các nguyên tử lớn thì Evdw bằng 0, nếu khoảng cách nhỏ thì chúng đẩy
nhau. Tương tác Van der Waals bao hàm tương tác khuếch tán, tương tác cảm ứng,
tương tác lưỡng cực-lưỡng cực, tứ cực-lưỡng cực, ... Evdw rất dương ở khoảng cách
nhỏ, có cực tiểu hơi âm tại khoảng cách tương ứng khi hai nguyên tử chỉ vừa chạm
nhau, và tiến tới 0 tại khoảng cách vô cùng. Một dạng hàm thỏa mãn tính chất này
là
( ) ( ) 6R
CRERE repulvdw −=
(1.41)
23
Ta không thiết lập được biểu thức chính xác cho Erepul. Số hạng này phải tiến
tới 0 khi R tiến tới vô cùng và phải tiệm cận 0 nhanh hơn số hạng thứ hai. Một trong
những dạng thường sử dụng trong tính toán là biểu thức Lennard-Jones:
( ) 62121 R
C
R
CRELJ −=
(1.42)
Dạng hàm này cũng được hiệu chỉnh phù hợp với các phương pháp tính mức
cao để cho kết quả tin cậy hơn.
1.2.3.5. Năng lượng tĩnh điện Eel
Do sự phân bố các electron trên phân tử mà hình thành nên những phần tích
điện dương và âm. Ta có thể mô tả tương tác này như tương tác giữa các điểm tích
điện bằng cách phân bổ điện tích cho mỗi nguyên tử hoặc xem liên kết như một
lưỡng cực. Hai mô hình này tương đương về mặt vật lí nhưng cho kết quả không
hoàn toàn giống nhau trong tính toán số.
Với tương tác giữa các điện tích điểm
( )
R
QQRE
BA
el ε
=
(1.43)
(ε là hằng số điện môi)
Điện tích của nguyên tử cũng được phân bổ cho phù hợp với các phương pháp tính
mức cao hoặc dữ kiện thực nghiệm giống như các tương tác khác trong trường lực.
Với mô hình liên kết lưỡng cực, biểu thức viết cho năng lượng tương tác
giữa 2 lưỡng cực có dạng sau:
( ) ( )BA
BA
el R
RE ααχ
ε
µµ coscos3cos3 −= (1.44)
Không có cơ sở chặt chẽ nào cho việc chọn giá trị hằng số điện môi ε, thông
thường giá trị này được lấy trong khoảng từ 1 đến 4.
1.2.3.6. Các số hạng chéo Ecross
24
Thực tế không có sự tách biệt hoàn toàn giữa các tương tác đã nêu ở trên khi
cho phân tử biến dạng để xây dựng bề mặt thế năng. Để mô tả ảnh hưởng qua lại
của các tương tác này người ta đưa thêm vào số hạng Ecross và thường viết ở dạng
tích của các khai triển Taylor.
Như vậy, đối với các phương pháp trường lực, các phép tính đều được viết ở
dạng cơ học cổ điển, do đó cho kết quả tính toán nhanh chóng, vấn đề cốt yếu của
các phương pháp này là xác định các tham số để đưa vào biểu thức tính. Khi tính
cho hệ các phân tử lớn, không thể xác định tham số cho từng nguyên tử, từng liên
kết cụ thể và cũng không thể xác định lại các tham số khi nghiên cứu các hệ phân tử
khác nhau. Vì vậy cần xây dựng bộ tham số có tính chất khái quát và rút gọn. Trong
mỗi bộ tham số được xây dựng, người ta xác định tham số cho các dạng nguyên tử
theo số hiệu nguyên tử và tính chất liên kết mà nó tham gia. Những bộ tham số này
đều phải phù hợp tương đối với thực nghiệm hoặc các phương pháp tính toán lượng
tử mức cao. Không có bộ tham số nào tuyệt đối tốt hơn các bộ tham số khác và
không thể đưa tham số từ trường lực này vào trường lực khác. Một số trường lực
phổ biến như UFF, Dreiding, Amber, CHARM ...
1.2.4. Kết hợp phương pháp cơ học lượng tử-cơ học phân tử
Một trong những khó khăn chủ yếu của hóa học tính toán khi nghiên cứu các
hệ lớn là cân bằng giữa độ chính xác của kết quả và thời gian tính toán. Các phương
pháp tính toán lượng tử tuy cho kết quả chính xác nhưng lại không thích hợp với
những hệ như vậy do khối lượng tính toán quá lớn và việc thực hiện tính lượng tử
cho những hệ hàng nghìn nguyên tử là điều hoàn toàn không khả thi.
Phương pháp kết hợp lợi dụng đặc điểm là trong hầu hết các phản ứng với
xúc tác enzyme, quá trình phá vỡ và hình thành liên kết chỉ xảy ra trên tâm hoạt
động có sự tham gia của một số ít các nguyên tử trong phân tử protein, ảnh hưởng
của phần còn lại trên protein thường chỉ là về mặt không gian và tương tác tĩnh
điện. Trong phương pháp kết hợp, mỗi vùng được xử lí bằng một phương pháp tính
khác nhau. Phần hoạt động hóa học được xử lí bằng phương pháp tính toán lượng tử
25
mô tả chính xác sự phá vỡ, hình thành liên kết hóa học, phần còn lại có thể được xử
lí bằng các phương pháp đỡ tốn kém thời gian hơn. Nhờ đó vừa mô tả được quá
trình hóa học vừa tiết kiệm thời gian.
Trước khi phương pháp kết hợp được áp dụng rộng rãi, mô hình tâm hoạt
động là một giải pháp nghiên cứu hoạt tính xúc tác của enzyme. Khi đó, chỉ một
phần phân tử có tâm hoạt động được sử dụng trong mô hình nghiên cứu và việc tính
toán bằng các phương pháp lượng tử không mấy khó khă n. Tuy nhiên, nhiều
enzyme có hoạt tính xúc tác cao mà các xúc tác khác không có được và đặc biệt là
có tính đặc thù, đặc điểm này không thể giải thích bằng một phần nhỏ trong phân
tử. Mô hình tâm hoạt động không chỉ ra được ảnh hưởng của toàn bộ phân tử
enzyme và trong nhiều trường hợp không thể hiện được vai trò xúc tác của enzyme.
Vì vậy, hiện tại phương pháp kết hợp đang là phương pháp hiệu quả nghiên cứu hệ
xúc tác enzyme.
Về nguyên tắc, có nhiều kiểu kết hợp các phương pháp tính với nhau nhưng
phổ biến là kết hợp giữa phương pháp lượng tử và cơ học phân tử - QM/MM.
Phương pháp kết hợp lượng tử/cơ học phân tử được tiên phong bởi Warshel
và Levitt vào năm 1976. Trong đó, Warshel đưa ra biểu thức nă ng lượng khi dùng
kết hợp phương pháp như sau
EEMMQME −/ = QME ,ν + MME + MMQME − (1.45)
với QME ,ν là năng lượng của vùng QM trong trường ν tạo ra bởi điện tích riêng
phần của vùng MM
MME là năng lượng của vùng MM chứa tất cả các số hạng MM liên kết và
không liên kết liên quan đến các tâm nằm gọn trong vùng MM).
MMQME − thể hiện tương tác giữa hai vùng và gồm hai thành phần: một là nếu có
liên kết cộng hóa trị giữa vùng QM và MM thì nó sẽ chứa các số hạng MM liên kết
qua biên (liên quan đến cả các tâm trong vùng QM và MM); thứ hai, nó gồm tất cả
26
các số hạng MM cho tương tác Van der Waals liên quan đến một tâm QM và một
tâm MM. MMQME − không chứa tương tác tĩnh điện giữa vùng QM và MM.
Kollman đưa ra một biểu thức dạng khác như sau
MEMMQME −/ = QME + MME + MMQMQE −, (1.46)
Ở đây, QME không còn tính đến ảnh hưởng từ vùng MM nữa. Thay vào đó, tương
tác tĩnh điện giữa các vùng được tính vào MMQMQE −, bằng cách áp điện tích riêng
phần cho các nguyên tử trong vùng QM, và dùng các biểu thức tính thông thường
cho tương tác giữa các điện tích điểm theo trường lực MM.
Trong phần mềm Gausian, phương pháp kết hợp được thực hiện với kĩ thuật
ONIOM. Phân tử được chia thành các vùng ở mức cao và mức thấp, mỗi vùng áp
dụng một phương pháp tính, trong phương pháp QM/MM thì vùng cao áp dụng một
phương pháp cơ học lượng tử, vùng thấp áp dụng phương pháp cơ học phân tử.
Năng lượng được tính như sau
LowModelLowalHighModelONIOM EEEE ,,Re, −+= (1.47)
ở đây Real chỉ toàn bộ hệ thực, Model chỉ vùng QM, High chỉ phương pháp áp
dụng ở mức cao, Low chỉ phương pháp áp dụng cho mức thấp. Khi phân vùng, một
số liên kết có thể bị cắt, do đó cần phải đưa một nguyên tử “ảo” vào để thay thế
phần bị cắt, những nguyên tử này ghép vào với phần mức cao để tạo Model.
Do áp dụng các phương pháp khác nhau cho mỗi vùng nên có sự gián đoạn
qua phần phân cắt. Vì thế để đảm bảo tính chính xác và liên tục trên bề mặt thế
năng, vị trí phân cắt cần phải xa tâm phản ứng hóa học và không cắt qua các bộ
phận cứng nhắc trong phân tử.
27
CHƯƠNG 2. NGUỒN DỮ LIỆU VÀ CÔNG CỤ TÍNH TOÁN
2.1. Nguồn dữ liệu
Protein Data Bank (PDB) là kho lưu trữ dữ liệu cấu trúc 3-D của các phân tử
sinh học lớn như là protein và axit nucleic. Các file dữ liệu cấu trúc được đưa lên
bởi các nhà sinh học và hóa sinh từ khắp thế giới, có thể truy cập và tải về miễn phí
qua các trang web thành viên PDBe, PDBj, RCSB. Dữ liệu đưa lên Protein Data
Bank được kiểm tra lại bằng phần mềm PDB Validation Suite.
Các phương pháp thường được sử dụng để xác định cấu trúc trên PDB là xác
định cấu trúc tinh thể dùng tia X, phương pháp phổ cộng hưởng từ NMR, và
phương pháp hiển vi điện tử nhiệt độ thấp. Các cấu trúc sử dụng trong luận văn đều
thu được bằng phương pháp nhiễu xạ tia X. Mỗi cấu trúc đều ghi chú rõ ràng độ
phân giải của dữ liệu.
Độ phân giải là thước đo chất lượng của dữ liệu được tập hợp. Nếu tất cả các
protein ở những điểm tương đương trong các tinh thể định hướng giống nhau thì ta
sẽ thu được tinh thể hoàn hảo, khi đó tất cả các protein sẽ phân tán tia X cùng một
kiểu như nhau và nhiễu xạ đồ thu được sẽ thể hiện được thông tin chi tiết về tinh
thể, vị trí của các nguyên tử có thể xác định được rõ ràng. Nhưng nếu không có
được tinh thể hoàn hảo, do tính mềm dẻo của từng phần trong protein và do các
phân tử protein lớn, khi kết tinh không định hướng như nhau thì nhiễu xạ đồ sẽ thể
hiện thông tin cấu trúc kém chi tiết hơn. Nói cách khác, độ phân giải là thước đo
mức độ thể hiện chi tiết của nhiễu xạ đồ và do đó là thước đo mức độ chi tiết khi
tính mật độ electron. Với độ phân giải cao thì có thể nhìn thấy ngay vị trí của mọi
nguyên tử từ bản đồ mật độ electron, còn với độ phân giải thấp thì chỉ thấy được
khung của chuỗi protein. Từ nhiễu xạ đồ có thể lập được bản đồ mật độ electron và
dự đoán được vị trí của các nguyên tử. Từ cấu trúc dự đoán tính ngược lại mật độ
electron để khớp với bản đồ mật độ electron từ nhiễu xạ đồ. Quy trình này được lặp
đi lặp lại cho đến khi có độ phù hợp mong muốn. Độ phân giải xác định giới hạn
của dữ liệu nhiễu xạ.
28
Nếu độ phân giải lớn hơn 4.0 Ǻ thì không thể xác định được tọa độ riêng rẽ
của các nguyên tử.
Độ phân giải trong khoảng 3.0 – 4.0 Ǻ: cấu trúc bộ khung có thể đúng nhưng
phần mạch nhánh có cấu dạng không đáng tin cậy.
Độ phân giải trong khoảng 2.5 – 3.0 Ǻ: bộ khung có thể xác định đúng, với
các nhánh dài và mảnh của một số aminoaxit như Lys, Glu, Gln, ... và các nhánh
nhỏ của Ser, Val, Thr, ... có cấu dạng không đáng tin cậy.
Độ phân giải trong khoảng 2.0 – 2.5 Ǻ: số nhánh có cấu dạng sai ít hơn đáng
kể. Có thể xác định được các phân tử nước và các phối tử nhỏ.
Độ phân giải trong khoảng 1.5 – 2.0 Ǻ: chỉ còn lỗi nhỏ về cấu dạng.
Độ phân giải trong khoảng 0.5 – 1.5 Ǻ: tọa độ của các nguyên tử được xác
định với độ tin cậy cao.
Hình 2.1 minh họa dữ kiện nhiễu xạ ở các mức độ phân giải khác nhau.
29
Hình 2. 1 Ảnh hưởng của độ phân giải đến khả năng xác định chính xác cấu trúc từ
nhiễu xạ đồ
Ở hình 2.1, các đường màu xanh và màu vàng bao quanh vùng có mật độ
electron lớn. Với độ phân giải 1.0 Å, bản đồ mật độ thể hiện ngay vị trí các nguyên
tử. Với độ phân giải 2.0 Å có thể xác định được các cấu trúc vòng, dự đoán được
các đơn vị aminoaxit. Với độ phân giải 2.7 Å có thể dự đoán được cấu trúc vòng,
mạch nhánh khó xác định. Còn với độ phân giải 3.0 Å, mật độ có dạng hình ống,
khó xác định cấu trúc chính xác.
30
Cấu trúc 2ACE xác định bằng phương pháp nhiễu xạ tia X với độ phân giải
2.5 Ǻ, bản đồ mật độ electron và cấu trúc dự đoán được ghi lại trong file 2ACE.pdb
thể hiện trong hình 2.2.
Hình 2. 2 Bản đồ mật độ electron và cấu trúc dự đoán của 2ACE
Từ hình trên ta thấy trên cơ chất acetylcholine có một số nguyên tử không
xác định được chính xác vị trí từ nhiễu xạ đồ. Cấu trúc trên được dự đoán và tinh
chỉnh cho khớp với bản đồ mật độ electron và không thể tránh khỏi sai số.
Dữ liệu cấu trúc của protein còn được lấy từ PDBsum, UniProt. UniProt
cung cấp thông tin chi tiết về cấu trúc sơ cấp, thứ cấp và ghi chú về chức năng sinh
học, hoạt tính xúc tác, các tâm hoạt động.
Đối chiếu với thứ tự các aminoaxit, cấu trúc 2ACE thiếu các aminoaxit từ
485 đến 489. Đoạn còn thiếu và một số sai sót ở các nhánh cũng được chỉnh sửa
dùng phần mềm Alcelrys MS Modeling 4.0, Acelrys Discovery Studio Visualizer
2.5 và GaussView 3.0 dựa vào cấu trúc thứ cấp và tham khảo cấu trúc 1CFJ.
2.2. AutoDock 4.2 và AutoDockTools 1.5.4
31
AutoDock 4.2 là phiên bản mới nhất trong chuỗi phần mềm AutoDock, sản
phẩm của The Scripps Research Institute. Đây là phần mềm mã nguồn mở được sử
dụng trong luận văn với mục đích khảo sát docking đối với acetylcholine lên phân
tử enzyme acetylcholinesterase. AutoDock 4.2 được dùng kèm với AutoDockTools
1.5.4 để hỗ trợ giao diện đồ họa.
Các bước tiến hành docking:
• Bước 1: Chuẩn bị cấu trúc enzyme và cơ chất dưới dạng file .pdbqt
Cấu trúc acetylcholinesterase từ file 2ACE được bổ sung, chỉnh sửa, loại bỏ
cơ chất, các phân tử nước, thêm H (dùng phần mềm Alcelrys Discovery
Studio Visualizer 2.5), dạng nguyên tử (dùng để xác định tham số trường
lực) và điện tích sẽ được tự động thêm vào trong AutoDockTools và được
ghi lại dưới dạng file .pdbqt.
Cấu trúc của acetylcholine được tối ưu hóa sơ bộ dùng Gaussian 03W, sau
đó dạng của các nguyên tử trong phân tử cũng được tự động gán cho và ghi
lại dưới dạng file .pdbqt trong AutoDockTools.
Xác định phần phân tử enzyme có khả năng chuyển động linh hoạt, ghi riêng
rẽ cấu trúc cứng và phần cấu trúc có thể chuyển động dưới dạng file .pdbqt.
• Bước 2: Tính trên AutoGrid
Xác định không gian khảo sát, khoảng cách trong lưới điểm. Trong phân tử
docking có bao nhiêu dạng phân tử thì AutoGrid sẽ cho ra kết quả là bấy
nhiêu bản đồ ghi thế của các dạng nguyên tử đó trong không gian mạng lưới
khảo sát dưới tác dụng của phần cấu trúc cứng. AutoGrid cũng tạo ra 2 bản
đồ thế là thế tĩnh điện và thế khử solvat hóa cho mọi trường hợp. Các thế này
sẽ được dùng khi đánh giá năng lượng trong quá trình docking.
Trong trường hợp khảo sát phân tử acetylcholine gắn kết lên
acetylcholinesterase với nhánh của Ser(200) có thể chuyển động, AutoGrid
32
tạo ra 4 bản đồ cho 4 loại nguyên tử C, OA, N, HD và 2 bản đồ thế tĩnh điện
và thế khử solvat hóa.
Hình 2.3 minh họa các bản đồ thế trên được ghép chung và bản đồ thế cho C.
• Bước 3: Docking
Thuật toán tìm kiếm được sử dụng là kết hợp thuật giải di truyền với tối ưu
cục bộ.
Các thông số cho thuật giải di truyền
Kích thước quần thể: 150 (autodock cho phép trong khoảng 50 đến
200)
Đánh giá năng lượng tối đa: 2500000 lần
Số thế hệ khảo sát tối đa: 27000
Tỉ lệ đột biến: 0.02
Tỉ lệ lai ghép: 0.8
Hình 2. 3 Bản đồ thế tạo ra bởi AutoGrid ghép chồng và bản đồ thế riêng cho C
33
Lai ghép và đột biến xảy ra tại 2 điểm trên nhiễm sắc thể.
Số cấu dạng đầu ra được xác định trong từng trường hợp cụ thể.
Các thông số còn lại sử dụng mặc định của phần mềm.
2.3. AutoDock Vina 1.1.1
AutoDock Vina là phần mềm mã nguồn mở từ The Scripps Research
Institute được dùng trong luận văn để khảo sát docking.
AutoDock Vina cũng dùng các file cấu trúc như AutoDock, có thể dùng kết
hợp AutoDockTools để hỗ trợ nhưng không cần tính với AutoGrid, việc tính toán
với các lưới điểm được thực hiện tự động. AutoDock Vina có ưu điểm là tính toán
nhanh, kết quả được trình bày rõ ràng, tự động sắp thứ tự, số lượng cấu dạng đầu ra
được tự động xác định theo khoảng năng lượng tránh lỗi rơi vào vòng lặp không xác
định như trong AutoDock, thích hợp cho việc khảo sát sơ bộ. Các thông số cho
thuật toán tìm kiếm được đặt sẵn từ kết quả khớp với ngân hàng cơ sở dữ liệu PDB
khi xây dựng phần mềm.
2.4. Gaussian 03W và GaussView 3.0
Gaussian là phần mềm tính toán hóa học được sử dụng rộng rãi. Phương
pháp QM/MM được thực hiện trên phần mềm này. Gaussian dùng kết hợp với
GaussView để hỗ trợ giao diện đồ họa, thuận tiện cho việc chỉnh sửa cấu trúc, gán
các tham số khi chuẩn bị file đầu vào.
Các bước tiến hành:
• Bước 1: Chuẩn bị cấu trúc
File cấu trúc protein ở dạng .pdb được bổ sung, chỉnh sửa dùng Accelrys MS
Modeling 4.0 và Accelrys Discovery Studio 2.5.
• Bước 2: Bổ sung các tham số trường lực
34
Bổ sung các dạng nguyên tử và các tham số trường lực theo bộ tham số
AMBER parm96 trong gói Amber10. Amber là bộ tham số được dùng phổ
biến khi nghiên cứu các phân tử protein và axit nucleic.
Vì trường lực Amber không có đủ các tham số cho nhiều phân tử hữu cơ
(thường gặp đối với các cơ chất khi nghiên cứu hệ xúc tác enzyme) nên tham
số cho phân tử cơ chất được bổ sung theo bộ tham số GAFF dùng
AmberTools-1.4. GAFF là bộ tham số phù hợp với trường lực Amber và nó
bao gồm hầu hết tham số cho các phân tử hữu cơ chứa C, N, O, H, S, P và
các nguyên tử halogen.
• Bước 3: Phân mức tính
Để tính toán đạt hiệu quả, việc phân mức cần tuân theo một số quy tắc
o Quy tắc 1:
Phần xảy ra phản ứng hóa học cần nằm trọn trong vùng được tính toán
bằng phương pháp chính xác (mức cao). Trong trường hợp nghiên
cứu, tương tác hóa học trực tiếp xảy ra giữa acetylcholine và nhóm
OH trên Ser(200). Như vậy việc phân vùng sẽ cắt qua 2 liên kết trên
phân tử acetylcholinesterase.
o Quy tắc 2:
Phương pháp dùng ở mức thấp phải nhanh nhưng vẫn đảm bảo mô tả
được các hiệu ứng gây ra bởi vùng này. Phần còn lại của phân tử
enzyme chủ yếu gây ra hiệu ứng về mặt không gian và tĩnh điện đối
với cơ chất nên phương pháp MM với các thành phần trường lực đã
nói ở chương 1 có thể đáp ứng được quy tắc này.
o Quy tắc 3:
Vì khi cắt qua liên kết, để đảm bảo tính chất của hệ hóa học, liên kết
bị cắt trong phần QM phải được thay thế bằng liên kết với một
nguyên tử khác. Và vì chỉ một nguyên tử duy nhất được dùng để thay
35
thế cho mỗi vị trí bị cắt, liên kết bị cắt nên lựa chọn là liên kết đơn và
nguyên tử thay thế là nguyên tử tạo liên kết đơn duy nhất.
o Quy tắc 4:
Nguyên tử thay thế không tham gia vào quá trình hóa học và cách xa
tâm phản ứng hóa học (tức là phần QM có kích thước lớn) để đảm bảo
phần sai lệch ảnh hưởng không đáng kể.
Cụ thể, quá trình phá vỡ và hình thành liên kết nên xảy ra cách vùng
MM từ 3 liên kết trở lên để đảm bảo tính liên tục của phần MM. Vì
các số hạng liên kết trong phần MM kéo dài trên 4 tâm nên phần có
liên kết mới hình thành phải cách phần MM 3 liên kết trở lên.
Ngoài ra, phần phân cắt cũng phải giống nhau trong hệ chất phản ứng
và sản phẩm để khi so sánh giữa hệ phản ứng và hệ sản phẩm, sai số
có thể được loại trừ.
o Quy tắc 5:
Không cắt qua các cấu trúc vòng, đặc biệt là vòng nhỏ vì các cấu trúc
vòng thường cứng nhắc, đối với vòng bé, sức căng vòng lớn, thay thế
liên kết trong vòng bằng một liên kết đơn làm mất tính chất này.
o Quy tắc 6:
Trong trường hợp nghiên cứu phải cắt qua 2 liên kết, do đó cần thay
thế bằng 2 nguyên tử, sai số sẽ không đáng kể nếu sai khác giữa
nguyên tử thay thế và nguyên tử bị thay thế bù trừ cho nhau ở 2 phía
bị cắt.
Áp dụng các quy tắc trên vào hệ nghiên cứu, phần QM được chọn là
Ser(200) và phân tử acetylcholine, hai nguyên tử tiếp giáp với vùng QM được thay
thế bằng 2 nguyên tử H.
36
CHƯƠNG 3. KẾT QUẢ VÀ THẢO LUẬN
3.1. Protein docking
Tiến hành Docking với AutoDock Vina, ban đầu toàn bộ phân tử
acetylcholinesterase đưa vào dưới dạng cấu trúc cứng. Không gian khảo sát xác
định như trong hình 3.1.
Hình 3. 1 Không gian khảo sát docking với AutoDock Vina
Không gian khảo sát gồm toàn bộ phân tử acetylcholinesterase. Các thông số
như sau:
center_x = 5.116
center_y = 64.613
center_z = 55.588
size_x = 60
size_y = 60
37
size_z = 60
Thực hiện docking phân tử acetylcholine, kết quả cho ra 9 cấu dạng.
mode
affinity
(kcal/mol)
dist from best mode
(rmsd l.b. Ǻ)
1 -4.9 0
2 -4.6 2.174
3 -4.5 2.801
4 -4.4 5.672
5 -4.3 3.69
6 -4.3 4.82
7 -4.3 2.635
8 -4.2 2.807
9 -4.2 5.052
Trong cả 9 trường hợp acetylcholine đều nằm trong cùng một hốc của phân
tử acetylcholinesterase (hình 3.2).
Hình 3. 2 Các cấu dạng gắn kết có ái lực âm nhất của acetylcholine lên
acetylcholinesterase
Đây chính là hốc phản ứng có đơn vị Ser(200) ở đáy. Kết quả này phù hợp
với dữ kiện thực nghiệm về vị trí phản ứng.
38
Trường hợp tiếp theo, tiến hành docking cũng với AutoDock Vina, không
gian khảo sát lấy cùng cỡ, số cấu dạng tối đa là 100, khoảng cách mức năng lượng
gắn kết giữa cấu hình có năng lượng thấp nhất với cấu hình có năng lượng cao nhất
tối đa là 10kcal/mol (ở trong trường hợp trên, số cấu dạng tối đa là 9, khoảng cách
mức năng lượng gắn kết tối đa là 3kcal/mol). Mạch nhánh của đơn vị aminoaxit
Ser(200) được để tự do. Kết quả cho ra 20 cấu dạng.
mode
affinity
(kcal/mol)
dist from best mode
(rmsd l.b. Ǻ)
1 -5.5 0
2 -5.1 2.007
3 -4.8 5.288
4 -4.7 3.353
5 -4.6 4.457
6 -4.5 4.454
7 -4.5 7.18
8 -4.5 8.49
9 -4.4 3.972
10 -4.3 5.391
11 -4.3 2.425
12 -4.2 12.932
13 -4.1 8.245
14 -4.1 8.16
15 -4 8.948
16 -4 5.483
17 -4 12.809
18 -3.8 14.354
19 -3.8 16.418
20 -3.8 24.441
Trong đó ở 11 dạng đầu, phân tử acetylcholine cũng nằm trong hốc có
Ser(200), 3 trạng thái nằm ở cửa hốc phản ứng; ở các trạng thái còn lại
acetylcholine neo đậu tại các vị trí bên ngoài của phân tử acetylcholinesterase
nhưng tại các vị trí này ngoài yếu tố tĩnh điện và liên kết hydro, không có điều kiện
để phản ứng hóa học xảy ra.
Trong những cấu dạng đầu, acetylcholine định hướng ngay trên đơn vị
aminoaxit Ser(200). Ta nhận xét thấy khi Ser(200) được cho chuyển động thì ái lực
của các cấu dạng đầu âm hơn so với trường hợp để cả phân tử acetylcholinesterase
cứng nhắc.
39
Hình 3.3, 3.4 là 2 cấu dạng có ái lực âm nhất.
Hình 3. 3 Cấu dạng có ái lực âm nhất
40
Hình 3. 4 Cấu dạng có ái lực âm thứ hai
Ngoài yếu tố tĩnh điện và liên kết hydro, phân tử acetylcholine có xu thế định
hướng nhóm amin bậc 4 về phía vòng Tryptophan.
Tiến hành docking dùng AutoDock với kích thước không gian khảo sát nhỏ
hơn, bao hốc phản ứng, mạch nhánh của Ser(200) có thể quay tự do.
41
Hình 3. 5 Không gian khảo sát với AutoDock
Tiến hành docking với 2 trường hợp, tìm 10 cấu dạng và 20 cấu dạng. Kết
quả trong cả 2 trường hợp các, phân tử acetylcholine đều nằm trong hốc phản ứng,
định hướng trên Ser(200) và hướng nhóm amin bậc 4 về vòng Tryptophan.
Trường hợp tìm 10 cấu dạng
Mode
Affinity
(kcal/mol)
dist from best mode
(rmsd l.b. Ǻ)
1 -4.57 0
2 -4.56 0.44
3 -4.55 0.82
4 -4.52 0.77
5 -4.36 0.46
6 -4.26 1.19
7 -4.09 0.84
8 -4.08 1.05
9 -3.88 1.34
10 -3.88 1.27
Trường hợp tìm 20 cấu dạng
42
Mode
Affinity
(kcal/mol)
dist from best mode
(rmsd l.b. Ǻ)
1 -5.22 0
2 -5.22 0.69
3 -5.16 0.81
4 -4.92 0.6
5 -4.83 0.86
6 -4.8 0.87
7 -4.7 0.48
8 -4.69 0.84
9 -4.59 0.9
10 -4.57 0.9
11 -4.54 0.8
12 -4.39 0.94
13 -4.25 0.97
14 -4.24 1.32
15 -4.19 0.96
16 -3.98 1.42
17 -3.96 1.32
18 -3.83 1.32
19 -3.77 1.48
20 -3.61 1.4
Dưới đây là 2 cấu dạng có năng lượng gắn kết âm nhất trong mỗi trường hợp
(hình 3.6, 3.7).
Trường hợp 10 cấu dạng:
Hình 3. 6 Hai cấu dạng bền nhất trong trường hợp 10 cấu dạng
43
Trường hợp 20 cấu dạng:
Kết quả trong cả 2 trường hợp tương đối giống nhau cả về năng lượng và các
cấu dạng. Mặc dù phương pháp tìm kiếm là ngẫu nhiên nhưng kết quả có tính ổn
định và đáng tin cậy.
Như vậy từ kết quả docking trên cả 2 phần mềm AutoDock và AutoDock
Vina với 2 phương pháp khác nhau ta có thể bước đầu nhận định acetylcholine có
khả năng gắn kết cao trong hốc chứa đơn vị aminoaxit Ser(200) và có xu hướng
quay nhóm amin bậc 4 về phía vòng Tryptophan. Thực nghiệm đã xác nhận kết quả
này. Đồng thời, dựa vào kết quả khảo sát trên AutoDock Vina trong không gian lớn
chứa toàn bộ phân tử có thể đưa ra nhận định là quá trình dịch chuyển của
acetylcholine từ bên ngoài vào trong hốc phản ứng không chỉ tuân theo định luật
khuếch tán thông thường mà có thể qua các trạng thái neo đậu là những trạng thái
bền cục bộ.
3.2. Áp dụng phương pháp QM/MM đối với hệ phản ứng
3.2.1. Cấu trúc enzyme
Hình 3. 7 Hai cấu dạng bền nhất trong trường hợp 20 cấu dạng
44
Cấu trúc của acetylcholinesterase xây dựng từ dữ kiện nhiễu xạ tia X ở dạng
file pdb (2ACE.pdb và 1CFJ.pdb) được chỉnh sửa, bổ sung cho đầy đủ các
aminoaxit.
Phân vùng tính toán: đơn vị aminoaxit Ser(200) được tính ở mức cao theo
phương pháp DFT, phiếm hàm B3LYP, phần còn lại trong phân tử
acetylcholinesterase được tính ở mức thấp, phương pháp cơ học phân tử, trường lực
AMBER. Các tham số trường lực được bổ sung theo bộ tham số parm96 (đi kèm
trong antechamber).
Khảo sát với một số bộ hàm cơ sở được kết quả tóm tắt như trong bảng sau:
Bảng 1. Kết quả và thời gian tính tối ưu cấu trúc enzyme với một số bộ hàm cơ sở
Bộ hàm cơ sở RMSD(Ǻ) E Thời gian tính (Job cpu time)
6-311++g(d,p) -358.75734276 a.u.
6-311++g(d) 0.045 -358.73338665 a.u. 7 hours 54 minutes 20.0 seconds.
6-311++g 0.026 -358.65200654 a.u. 9 hours 31 minutes 44.0 seconds.
6-311+g(d,p) 0.003 -358.75706059 a.u. 7 hours 4 minutes 39.0 seconds.
6-311+g(d) 0.006 -358.74010777 a.u. 9 hours 23 minutes 54.0 seconds.
6-311+g 0.026 -358.65148321 a.u. 8 hours 13 minutes 59.0 seconds.
6-311g(d,p) 0.014 -358.73760160 a.u. 6 hours 51 minutes 32.0 seconds.
6-311g(d) 0.011 -358.72085986 a.u. 4 hours 58 minutes 34.0 seconds.
6-311g 0.020 -358.63378965 a.u. 6 hours 2 minutes 25.0 seconds.
(Giá trị RMSD chỉ tính cho phần mức cao).
N
d
RMSD
N
i
ii∑
== 1
2
trong đó dii là khoảng cách giữa 2 nguyên tử tương ứng
ở 2 cấu trúc. RMSD được tính trên phầm mềm PyMol, 2 cấu trúc sẽ được sắp đặt để
độ lệch giữa chúng là nhỏ nhất. Hình học được so sánh với cấu trúc 1QID.
Từ bảng thống kê trên, để cân bằng giữa thời gian tính toán cũng như độ
chính xác về năng lượng và hình học ta có thể chọn bộ hàm cơ sở 6-311g(d) để
khảo sát.
45
Cơ chất acetylcholine cũng được tối ưu hóa theo phương pháp DFT, phiếm
hàm B3LYP, bộ hàm cơ sở 6-311g(d).
Hình 3.8 mô tả bề mặt của enzyme và acetylcholine đã tối ưu.
Hình 3. 8 Bề mặt của acetylcholine và hốc phản ứng trên acetylcholinesterase (hai cấu
trúc đã được tối ưu)
Ta thấy kích thước nhóm amin bậc bốn khá lớn so với kích thước cửa hốc
phản ứng và cơ chất không thể đi qua để vào hốc phản ứng nếu cấu trúc của protein
cứng nhắc do sự cản trở của Tyr(121) và Phe(330) (xem hình 3.12). Cùng với nhận
xét đã nêu khi tiến hành docking ta có thể khẳng định là enzyme
acetylcholinesterase có cấu trúc thay đổi linh hoạt trong quá trình phản ứng và
acetylcholine di chuyển vào trong hốc phản ứng không chỉ tuân theo định luật
khuếch tán thông thường mà qua các trạng thái neo đậu bền cục bộ. Acetylcholine
phải có ái lực tương đối lớn với một số điểm trên phân tử enzyme [5,6].
46
3.2.2. Cơ chất trong hốc phản ứng ở trạng thái chưa liên kết
Trạng thái của acetylcholine trong hốc phản ứng cũng được tối ưu dựa vào
định hướng trong quá trình docking. Tham số trường lực của các nguyên tử trên cơ
chất được bổ sung theo bộ tham số GAFF như đã nói ở chương 2.
Hình 3. 9 Trạng thái acetylcholine trong hốc phản ứng đã được tối ưu
Năng lượng tính được
ONIOM Total Energy = -840.19781353 a.u.
Sự biến đổi cấu trúc của cơ chất trong hốc phản ứng so với trạng thái tự do
(đã tối ưu) thể hiện trong hình 3.10.
47
Hình 3. 10 Sự biến đổi cấu trúc acetylcholine trong hốc phản ứng so với
trạng thái tự do
với giá trị độ lệch tiêu chuẩn
RMSD = 0.232 Ǻ
Cấu trúc của acetylcholine đã bị biến đổi dưới tác dụng của enzyme và ở
đây mới chỉ có ảnh hưởng của hiệu ứng không gian và tĩnh điện. Sự thay đổi cấu
trúc thể hiện chi tiết trong bảng 2.
(Cách đánh số các nguyên tử trong acetylcholine như hình 3.11).
48
Hình 3. 11 Thứ tự các nguyên tử trong acetylcholine
Từ các số liệu trong bảng 2 ta thấy độ dài của các liên kết khi acetylcholine
vào trong hốc phản ứng nói chung lớn hơn so với trạng thái tự do, sự biến đổi góc
liên kết thể hiện không rõ ràng. Một số góc xoắn có sự thay đổi đáng kể như góc
O(6)-C(5)-C(1)-H(4), C(11)-C(8)-O(7)-C(5). Độ dài liên kết C(5)-O(7) tăng lên đến
1.398 Ǻ, gần với độ dài liên kết C(8)-O(7) (1.404 Ǻ-độ dài của liên kết đơn ở trạng
thái tự do), tính chất đồng phẳng của các nguyên tử C(1), C(5), O(6) và O(7) giảm.
Những thay đổi trên tạo thuận lợi về hình học cho sự tiếp cận của C(5) lên O trong
nhóm OH trên Ser(200).
49
Bảng 2. Độ dài liên kết, góc liên kết và góc nhị diện của acetylcholine trong hốc phản ứng và trạng thái tự do.
STT Kí hiệu nguyên tử NA NB NC
Acetylcholine trong hốc phản ứng Acetylcholine tự do
Độ dài liên kết Góc liên kết Góc nhị diện Độ dài liên kết Góc liên kết Góc nhị diện
1 C
2 H 1 1.099 1.084
3 H 1 2 1.094 106.400 1.084 107.797
4 H 1 2 3 1.085 110.372 120.407 1.079 110.505 120.380
5 C 1 4 3 1.503 111.039 122.063 1.496 109.477 119.976
6 O 5 1 4 1.204 127.795 31.252 1.184 126.668 -3.525
7 O 5 1 6 1.398 109.599 176.044 1.342 112.014 179.862
8 C 7 5 1 1.421 119.154 -168.039 1.404 117.199 174.293
9 H 8 7 5 1.096 107.812 -134.332 1.080 105.487 -159.565
10 H 8 7 5 1.089 113.010 -14.261 1.079 109.597 -42.168
11 C 8 7 5 1.537 104.517 109.076 1.535 109.012 80.969
12 H 11 8 7 1.089 111.648 -64.094 1.078 109.661 -40.187
13 H 11 8 7 1.093 108.620 53.697 1.082 109.972 79.340
14 N 11 8 7 1.523 118.241 172.883 1.508 115.235 -159.868
15 C 14 11 8 1.504 106.808 -177.237 1.496 107.740 -172.925
16 H 15 14 11 1.089 109.895 -175.250 1.080 108.927 -177.908
17 H 15 14 11 1.091 109.457 -55.186 1.080 108.971 -58.027
18 H 15 14 11 1.091 109.242 63.652 1.080 109.140 62.099
19 C 14 11 8 1.506 112.351 63.361 1.499 110.964 68.688
20 H 19 14 11 1.091 108.722 59.319 1.079 108.916 55.755
21 H 19 14 11 1.090 108.007 177.632 1.080 108.329 175.842
22 H 19 14 11 1.090 109.145 -62.725 1.077 110.001 -64.198
23 C 14 11 8 1.505 110.549 -60.257 1.493 111.145 -53.851
24 H 23 14 11 1.092 108.126 -49.486 1.081 109.056 -54.885
25 H 23 14 11 1.089 109.628 70.606 1.079 109.652 65.755
26 H 23 14 11 1.090 108.674 -169.438 1.080 108.816 -174.632
50
Cấu trúc hốc phản ứng khi có acetylcholine cũng biến đổi so với khi chưa có
acetylcholine (hình 3.12).
Hình 3. 12 Sự biến đổi cấu trúc hốc phản ứng khi có (nét đậm) và
không có cơ chất (nét mảnh)
Các đơn vị aminoaxit trong hốc phản ứng có cấu trúc biến đổi nhi ều nhất, có xu
hướng giãn ra, các mạch nhánh được định hướng lại.
3.2.3. Cấu trúc phức enzyme-cơ chất
Phức giữa acetylcholine với enzyme cũng được tối ưu cho kết quả như hình
3.13.
51
Hình 3. 13 Cấu trúc phức enzyme-cơ chất
Tách riêng phần mức cao
Hình 3. 14 Cấu trúc phức giữa acetylcholine trên nhóm OH của Ser(200)
52
So sánh với cấu hình phức trong file pdb (trích từ 2ACE.pdb) được
RMSD = 0.851 Ǻ
Độ lệch này trong giới hạn chấp nhận được khi xem xét độ phân giải của dữ
liệu từ file pdb. Nếu chỉ tính phần Ser(200) thì RMSD = 0.117 Ǻ, sai lệch chủ yếu
là ở phần acetylcholine. Nếu chỉ xét phần cấu trúc khung của protein trong hốc phản
ứng thì độ lệch chỉ nhỏ hơn RMSD = 0.006 Å. (Chú ý rằng phần cơ chất không
được định vị chính xác từ dữ kiện nhiễu xạ tia X).
Năng lượng tính được cho cấu hình phức đã tối ưu
ONIOM Total Energy = -840.18639692 a.u.
Ở đây có sự biến đổi hoàn toàn góc liên kết và độ dài liên kết trong trạng thái
chưa tạo phức và trạng thái tạo phức đối với acetylcholine tại nhóm chức este. Các
góc liên kết quanh C(5) đều có giá trị gần với góc tứ diện đều (109o28’), cấu trúc
phẳng của nhóm C=O không còn nữa.
Hiệu năng lượng giữa trạng thái chưa liên kết và trạng thái phức:
∆E1 ≈ 29.98 (kJ/mol)
Tuy nhiên cấu trúc hốc phản ứng biến đổi không nhiều, chỉ có sự thay đổi ở
nhánh của Ser(200) là rõ ràng (hình 3.15). Như vậy ở giai đoạn này, cấu trúc chỉ
biến đổi mạnh tại phần có tương tác hóa học.
53
Hình 3. 15 Cấu trúc hốc phản ứng trong trạng thái tạo phức (nét đậm) và
chưa tạo phức (nét mảnh)
3.2.4. Cấu trúc sản phẩm
Cấu trúc tối ưu khi một phân tử choline tách ra:
54
Hình 3. 16 Cấu trúc phân tử choline tách ra và nhóm OH trên Ser(200) bị acetyl hóa
Năng lượng tính được cho cấu trúc này
ONIOM Total Energy = -840.27146076 a.u.
Lúc này cũng có sự biến đổi mạnh bản chất góc liên kết và độ dài liên kết
giữa trạng thái phức và trạng thái sản phẩm tại phần có tương tác hóa học. Các
nguyên tử C(1), C(5), O(7) và O(Ser) nằm trên mặt phẳng.
Hiệu năng lượng giữa trạng thái sản phẩm và trạng thái phức:
∆E2 ≈ -223.37 (kJ/mol)
Từ giá trị này ta thấy giai đoạn này tỏa nhiệt mạnh, thuận lợi về mặt nhiệt
động học. Phân tử choline được tạo ra tương đối bền trong hốc phản ứng.
Cấu trúc hốc phản ứng thay đổi không nhiều so với trạng thái phức, chỉ có sự
thay đổi ở mạch nhánh của một số aminoaxit phía cửa hốc (hình 3.17).
55
Hình 3. 17 Hốc phản ứng trong trạng thái sản phẩm (nét mảnh) và
trạng thái phức (nét đậm)
Trong cả 3 trạng thái với sự có mặt của cơ chất trong hốc phản ứng, tuy cấu
trúc cơ chất thay đổi trong khi cấu trúc hốc phản ứng không biến đổi nhiều nhưng
nhóm amin bậc 4 luôn hướng vào vòng Trp(84) (hình 3.9, 3.13, 3,16). Kết quả này
phù hợp với thực nghiệm [6].
Như vậy, dùng phương pháp ONIOM bước đầu đã khảo sát được cấu trúc
của trạng thái acetylcholine trong hốc phản ứng khi chưa tạo phức, trạng thái phức
giữa acetylcholine và enzyme, cấu trúc sản phẩm choline được tạo ra trong hốc
phản ứng và hiệu năng lượng của các trạng thái trên. Ta có thể dùng phương pháp
này để nghiên cứu chi tiết hơn đường phản ứng.
56
KẾT LUẬN
Từ các kết quả và những nhận xét trong chương 3 ta có thể đi đến một vài
kết luận sau:
• Sử dụng các công cụ tính toán AutoDock, AutoDock Vina cùng với
Gaussian, bước đầu đã thu được những kết quả phù hợp với thực nghiệm về
phản ứng thủy phân acetylcholine nhờ enzyme acetylcholinesterase.
• Bằng phương pháp QM/MM đã khảo sát được trạng thái phức giữa enzyme
acetylcholinesterase với acetylcholine, vai trò của enzyme trong phả n ứng
thủy phân acetylcholine, trong đó có vai trò hóa học của tâm xúc tác và vai
trò không gian, tĩnh điện của phần còn lại trong phân tử enzyme.
• Dự đoán được sự biến đổi cấu trúc enzyme cũng như cơ chất trong quá trình
phản ứng. Những kết quả này đã được xác nhận từ thực nghiệm.
Những kết quả thu được đã mở ra hướng nghiên cứu chi tiết cơ chế phản ứng
xúc tác của enzyme acetylcholinesterase, từ đó có thể dự đoán được hoạt tính ức
chế hay hoạt hóa của một số chất khác.
Định hướng phát triển đề tài:
• Nghiên cứu đầy đủ cơ chế phản ứng gồm các bước
o Quá trình di chuyển của acetylcholine vào bên trong hốc phản ứng.
o Quá trình di chuyển từ trạng thái chưa liên kết trong hốc phản ứng đến
trạng thái tạo phức.
o Quá trình từ trạng thái phức đến khi tạo thành sản phẩm.
Các bước trên sẽ được lập đường phản ứng chi tiết, tìm trạng thái chuyển
tiếp, năng lượng hoạt hóa của mỗi giai đoạn.
• Nghiên cứu phản ứng xúc tác của acetylcholinesterase thuộc các loài khác
nhau một mặt làm sáng tỏ thêm ảnh hưởng của toàn bộ phân tử enzyme, một
57
mặt chú trọng nghiên cứu enzyme acetylcholinesterase ở người vì đây là loại
enzyme rất quan trọng trong cơ thể.
• Từ kết quả về năng lượng thu được ở trên ta nhận thấy phương pháp
QM/MM chưa thể giải thích được thỏa đáng hoạt tính xúc tác đặc biệt cao
của acetylcholinesterase. Do giai đoạn tạo phức của phản ứng được dự đoán
là giai đoạn quan trọng, quyết định tốc độ phản ứng lại có hiệu năng lượng
dương nên có thể tiên đoán ảnh hưởng cung cấp nhiệt góp phần vào hoạt tính
xúc tác, làm tăng tốc độ phản ứn g. Để giải thích hiệu ứng này không thể chỉ
xem xét phương trình Schrodinger ở trạng thái dừng, mà c ó thể phải xét cả
phương trình Schrodinger phụ thuộc thời gian. Hiện tại chưa có công cụ nào
hỗ trợ hướng trên nên cần nghiên cứu thêm.
• Nghiên cứu ảnh hưởng của một số chất đến hoạt tính xúc tác của
acetylcholinesterase nhằm tìm ra các chất có lợi tác động đến hoạt động thần
kinh để khắc phục những hư tổn trong truyền dẫn thần kinh.
58
TÀI LIỆU THAM KHẢO
Tiếng Việt
1. Lâm Ngọc Thiềm (2008), Cơ sở hóa học lượng tử, NXB Khoa học và Kỹ
thuật, Hà Nội.
2. Đặng Ứng Vận (2003), Động lực học các phản ứng hóa học, NXB Giáo dục,
Hà Nội.
3. Đặng Ứng Vận (2007), Giáo trình Hóa tin cơ sở, NXB Đại học Quốc gia Hà
Nội.
Tiếng Anh
4. Axel D. Becke (1992), “Density-functional thermochemistry. III. The role of
exact exchange”, J. Chem. Phys., 98(7), pp. 5648-5652.
5. Yves Bourne, Palmer Taylor, Zoran Radic, Pascale Marchot (2003),
“Structural insights into ligand interactions at the acetylcholinesterase
peripheral anionic site”, The EMBO Journal, 22(1), pp. 1-12.
6. Jacques-Philippe Colletier et al. (2006), “Structural insights into substrate
traffic and inhibition in acetylcholinesterase”, The EMBO Journal, 25(12),
pp. 2746-2756.
7. P. Hohenberg, W. Kohn (1964), “Inhomogeneous Electron Gas”, Physical
Review, Vol.136(3B), pp. B864-B871.
8. Wolfram Koch, Max C. Holthausen, A Chemist’s Guide to Density
Functional Theory, Wiley-VCH, Germany.
9. Chengteh Lee, Weitao Yang, Robert G. Parr (1988), “Development of the
Colle-Salvetti correlation-energy formula into a functional of the electron
density”, Physical Review, 37(2), pp. 785-789.
59
10. Chérif F. Matta (2010), Quantum Biochemistry, Wiley-VCH, Weinheim,
Germany.
11. John P. Perdew, Kieron Burke, Matthias Ernzerhof (1996), “Generalized
Gradient Approximation Made Simple”, Physical Review Letters, 77(18), pp.
3865-3868.
12. Joel L. Sussman et al. (1991), “Atomic Structure of Acetylcholinesterase
from Torpedo californica: A Prototypic Acetylcholine-binding protein”,
Science, 253, pp. 872-879.
13. O. Trott, A. J. Olson (2010), “AutoDock Vina: improving the speed and
accuracy of docking with a new scoring function, efficient optimization and
multithreading”, Journal of Computational Chemistry, 31, pp. 455-461.
14. Darrin M. York, Tai-Sung Lee (2009), Multi-scale quantum models for
Biocatalysis: Modern Techniques and Applications, Springer.
Websites
15.
16.
17.
Structures/
60
PHỤ LỤC
BẢNG KÍ HIỆU AMINOAXIT TRONG CHUỖI POLYPEPTIT
Alanine Ala A
Arginine Arg R
Asparagine Asn N
Aspartic acid Asp D
Cysteine Cys C
Glutamic acid Glu E
Glutamine Gln Q
Glycine Gly G
Histidine His H
Isoleucine Ile I
Leucine Leu L
Lysine Lys K
Methionine Met M
Phenylalanine Phe F
Proline Pro P
Serine Ser S
Threonine Thr T
Tryptophan Trp W
Tyrosine Tyr Y
Valine Val V
Các file đính kèm theo tài liệu này:
- a1.PDF