CC BY
14
УДК 541.124
РЕАЛИЗАЦИЯ МЕТОДА МОЛЕКУЛЯРНОЙ ДИНАМИКИ С ЖЕСТКИМИ ФРАГМЕНТАМИ ДЛЯ МОДЕЛИРОВАНИЯ ХИМИЧЕСКИХ РЕАКЦИЙ В БИОМОЛЕКУЛЯРНЫХ СИСТЕМАХ
A.A. Московский, И.А. Калиман, A.B. Акимов, С.С. Конюхов, Б.Л. Григоренко, A.B. Немухин
(кафедра физической химии; e-mail: moskov@lcc.chem. msu. ru)
Описана новая реализация метода молекулярной динамики, направленного на моделирование свойств биомолекулярных систем, в которых могут происходить химические реакции. Для вычислений энергий и сил вдоль траекторий используется приближение квантовой и молекулярной механики на основе теории потенциалов эффективных фрагментов. В связи с особенностями метода эффективных фрагментов поведение молекулярно-механической подсистемы описывается при помощи динамики твердых тел. Методика применена для моделирования протонного транспорта вдоль молекул воды внутри канала грамицидина.
Метод молекулярной динамики (МД) представляет собой компьютерный эксперимент по моделированию поведения молекулярной системы во времени с помощью численного интегрирования уравнений движения [1, 2]. Подавляющее большинство МД-расчетов выполняется в рамках классической механики с использованием эмпирических потенциалов взаимодействия атомов, так называемых силовых полей. В приложении к биомолекулярным системам [3] метод классической МД позволяет решать большое число практически важных задач, прежде всего моделировать конформационные превращения белков и нуклеиновых кислот, но не позволяет рассматривать химические превращения, происходящие в системе. Одно из перспективных направлений развития МД, в рамках которого можно моделировать разрывы и образования химических связей в биомолекулярной системе, ориентируется на использование гибридного метода квантовой и молекулярной механики (КМ/ММ) [4-7] для вычислений энергий и сил, действующих на атомы. В этом подходе относительно небольшая часть всей системы описывается в рамках квантовомеханических приближений, тогда как для описания оставшейся части используются классические силовые поля.
Практическая реализация нового комбинированного метода квантовой механики и молекулярной динамики (КМ/МД), в котором силы, рассчитываемые в квантовой подсистеме по алгоритмам квантовой химии, потребляются компьютерными программами
молекулярной динамики и вместе с классическими силами позволяют вычислять траектории всех частиц в системе, весьма сложна, и к настоящему времени известны лишь немногочисленные попытки решения подобной задачи [8].
В данной работе описана реализация программы КМ/МД, использующая для расчетов энергии и градиентов энергии вариант метода квантовой и молекулярной механики (КМ/ММ) на основе теории потенциалов эффективных фрагментов [6, 7, 9]. В этом подходе молекулярные группы, относящиеся к моле-кулярно-механической (ММ) подсистеме, представляются эффективными фрагментами, которые вносят электростатический, поляризационный и обменный вклады в гамильтониан квантовой подсистемы. Взаимодействие между молекулярными группами (конфор-мационно-подвижными эффективными фрагментами), которые относят к ММ-подсистеме, описывается классическими силовыми полями. Соответственно для расчетов траекторий частиц в ММ-подсистеме необходима молекулярная динамика, описывающая перемещения твердых тел. В нашей предшествующей статье [10] обсуждались первые попытки применения молекулярной динамики с жесткими фрагментами для поиска конформаций минимальной энергии с использованием алгоритма обмена репликами. В настоящем сообщении приводятся примеры расчетов и анализа траекторий с новой реализацией программы КМ/МД.
Для описания движения твердых тел с применением явных алгоритмов интегрирования вводят две си-
стемы отсчета: внешнюю неподвижную систему отсчета Oxyz и внутреннюю жестко связанную с телом Oxy' z (рис. 1), базисные векторы которой коллине-арны главным осям инерции тела, а начало системы O совпадает с центром масс фрагмента [11]. Движение отдельного фрагмента, рассматриваемое как совокупность перемещения центра масс q = (qx, qy, qz) относительно Oxyz и вращения, описывается с помощью векторов полного импульса p = т\цм, определяемого относительно неподвижной системы отсчета и момента импульса п = 1ю, задаваемого относительно подвижной системы отсчета, где m, уцм - масса тела и скорость центра масс; I, ю - соответственно тензор момента инерции и угловая скорость тела. Взаимодействие фрагмента с остальной частью системы описывается суммарной силой
f = ^f. , действующей на тело, и суммарным мо-
i
ментом сил т = ^т., вычисленным относительно
i
центра масс.
Уравнения динамики, описывающие движение центра масс dp/dt = f, dq/dt = p/m, могут быть проинтегрированы либо в кватернионном приближении [12\, либо по методу матриц поворота [13]. Нами выбран последний вариант, в котором ориентация фрагмента задается с помощью матрицы Q, элементами которой являются направляющие косинусы между внешней системой отсчета Oxyz и внутренней Oxyz:
Q
cos(x л x') cos(x л y') cos(x л z ) COS(y л x') COs(y л y') C0s(y л z') C0s(z л x') C0s(z л y') C0s(z л z ).
где aAb - угол между векторами a и b. В этом случае уравнения динамики жесткого фрагмента принимают следующий вид [13]:
dp/dt = f,
dn/dt = т + пх( I"1 п),
dq/dt = p/m,
dQ/dt = Qskew( I"1 п),
где skew(a)
0
a„ 0
a.
3 —a
a 2 —a
0
Рис. 1. Выбор системы отсчета в динамике твердого тела
Нами был реализован численный алгоритм, позволяющий рассчитывать МД-траектории в молекулярных системах, разделенных на квантовую часть и молекулярно-механическую часть, представляющую собой конформационно-подвижные цепи эффективных фрагментов, перемещающихся как твердые тела, движение которых описывается вышеприведенными уравнениями. Значения энергии и силы в квантовой подсистеме вычисляются при помощи пакета программ РС ОАМЕ88 [14].
В качестве первого приложения была рассмотрена модельная система, состоящая из дипептида И-ацетил-Ь-аланин-К'-метиламида и четырех молекул воды [5] (рис. 2, а). К квантовой части отнесены молекулы воды, и расчеты энергий и сил проводились в приближении Хартри-Фока с базисом 6-ЗШ.
Молекула дипептида отнесена к молекулярно-ме-ханической подсистеме, подразделенной в свою очередь на 6 эффективных фрагментов (по два фрагмента СНЗ, СОИН и СН) и описываемой параметрами силового поля ОРЬ8АА [15]. Интегрирование уравнений КМ/МД проводилось для условий ИУЕ ансамбля с шагом интегрирования 0,5 фс с длинами траекторий до З пс. Было продемонстрировано удовлетворительное сохранение полной энергии на всем пути траекторий.
В другом примере мы рассматривали перемещение протона по ориентированной цепи молекул воды внутри двойной спирали пептида грамицидина, построенной на основе структуры ЦИО из Банка данных бел-
^ \
Рис. 2. Модельная система, состоящая из дипептида N-аце-тил-Ь-алаиии-Ы'-метиламида и четырех молекул воды (а); подразделение дипептида на эффективные фрагменты (б)
ковых структур (рис. 3). В процессе движения протона по цепи молекул воды происходят разрывы и образования химических связей (внутри молекулы воды) и водородных связей (между молекулами воды). В расчетах траекторий молекулы воды и дополнительный протон были отнесены к квантовой подсистеме, описываемой в приближении Хартри-Фока с базисом 6-31G. Молекулярные группы грамицидина были разделены на 144 эффективных фрагмента, взаимодействие между которыми описывалось параметрами силового поля AMBER [16]. В каждой точке молекулярно-динамической траектории каждый фрагмент из ММ-подсистемы вносит вклад в квантовый гамильтониан квантовой подсистемы. Интегрирование уравнений КМ/МД проводили для условий NVE-ансамбля с шагом интегрирования 0,5 фс с
длинами траекторий до 3 пс. Все расчеты проводили на обычных персональных компьютерах: Pentium IV 3,03 Ггц и Athlon XP 1800+. Наибольшие временные затраты приходятся на квантовохимические вычисления. Результаты данной работы доказывают возможность практической реализации комбинированного метода квантовой и молекулярной динамики. В дальнейшем планируется совершенствовать методику применения молекулярной динамики. Поскольку основные временные затраты приходятся на расчет в квантовомеханической части системы, это делает задачу КМ/МД моделирования достаточно нетипичной, так как в большинстве случаев время, затраченное на вычисление энергии и силы для части исследуемой системы, пропорционально размеру подсистемы. Эти особенности задачи могут быть использованы для сокращения временных затрат на проведение расчетов. В первую очередь могли бы быть использованы следующие приемы: 1) использование метода интегрирования уравнений движения с переменным шагом по времени (например, Бюрли-ха-Штоера [11]), 2) аппроксимация потенциальной поверхности в тех областях фазового пространства,
Рис. 3. Димер грамицидина с ориентированной цепью молекул воды внутри канала
молекулярно-динамической траекторией потенциального барьера [18].
10. MoskovskyA.A., Vanovschi V.V., Konyukhov S.S., Nemukhin A.V. // Intern. J. Quant. Chem. 2006. 106. P. 2208
11. ГолдстейнГ. Классическая механика. M., 1915. С. 109.
12. Miller T.F., Eleftheriou M., Pattnaik P. et. al. // J. Chem. Phys. 2002.116. P. 8649.
13. Dullweber A., Leimkuhler B., McLachlan R. // J. Chem. Phys. 1991.107. P. 5840.
14. Granovsky A. URL http://lcc.chem.msu.ru/gran/gamess/ index.html .
15. Kaminski G., Friesner R.A., Tirado-Rives J., Jorgensen W.L. // J . Phys. Chem. B 2001. 105. P. 6414.
16. Cornell W.D., Cieplak P., Bayly C.I. et al. //J. Amer. Chem. Soc. 1995. 117. P. 5119.
11. Press W.H., FlanneryB.P., TeukolskyS.A., Vetterling W.T.
Numerical Recipes. Cambridge, 1989. 18. Lange O.F., Schafer L.V., Grubmuller H.// J. Comput. Chem. 2006.27. P. 1693.
Поступила в редакцию 04.12.06
REALIZATION OF MOLECULAR DYNAMICS APPROACH WITH RIGID FRAGMENTS FOR SIMULATION OF CHEMICAL REACTIONS IN BIOMOLECULAR SYSTEMS
A.A. Mosckovskii, I.A. Kaliman, A.V. Akimov, S.S. Konyukhov, B.L. Grigorenko, A.V. Nemukhin
(Division of Physical Chemistry)
We describe a new implementation of the molecular dynamics method targeting simulations of biomolecular systems in which chemical reactions may occur. The quantum mechanical -molecular mechanical method based on the effective fragment potential theory is used for c a lculations of energies a nd forces along trajectories. In this approa ch, the molecular mechanical subsystem is described by the dynamics of rigid bodies, due to peculiarities of effective fragments theory. The method has been applied for modeling proton transfer along water molecules inside the gramicidin channel.
где траектория уже проходила, З) деформация потенциальной поверхности для ускоренного прохождения
СПИСОК ЛИТЕРАТУРЫ
1. Tildesley D.J. // Сотр. Sim. Chem. Phys. 1993. P. 23.
2. Tuckerman M.E. // J. Phys. Chem. 2000. 104. P. 159.
3. Karplus M., McCammon A.J. // Nature S tract. Biol. 2002. 9.
P. 646.
4. WarshelA., Levitt M. // J. Molec. Biol. 1976. 103. P. 2 27.
5. Nemukhin A.V., Grigorenko B.L., Bochenkova A.V. et. al. // J.
Molec. Struct. (Theochem) 2002. 581. P. 167.
6. Grigorenko B.L., Nemukhin A.V., TopolI.A., Burt S.K. //J.
Phys. Chem. A 2002. 106. P. 10663.
7. Nemukhin A.V., Grigorenko B.L., TopolI.A., Burt S.K. // J.
Comp. Chem. 2003. 24. P. 1410.
8. Parker C., Ventura O., Burt S., Cachau R. //Molec. Phys. 2003.
101. P. 2659.
9. Gordon M.S., Freitag M.A., Bandyopadhyay P. et. al. //
J. Phys. Chem. A 2001. 105. P. 293.
