«Молекулярная динамика. Расчёт ньютоновских траекторий движения»

Здесь также приводятся примеры молекулярного кино, созданного с помощью этих методов. Молекулярное кино, это отображение молекулярной конформационной подвижности в реальном времени. Приведённые далее примеры иллюстрируют "жизнь молекул" в диапазонах от 0,5 до 10 пикосекунд:

Расчёт ньютоновских траекторий движения.

В методе молекулярной динамики рассчитываются классические (ньютоновские) траектории движения атомов макромолекулы в силовом поле эмпирического атом-атомного потенциала, т. е. моделируется детальная микроскопическая картина внутренней тепловой подвижности макромолекулы в субнаносекундных интервалах времен. Основу метода составляет численное решение классических уравнений Ньютона для системы взаимодействующих частиц:

где - радиус-вектор i-го атома, - его масса, суммарная сила, действующая на i-ый атом со стороны остальных частиц:

Здесь: -потенциальная энергия, зависящая от взаимного расположения всех атомов; n - число атомов.

Задав координаты и скорости всех частиц в начальный момент времени, числено решают уравнения движения, вычисляя на каждом шаге все силы и новые координаты и скорости частиц. Температура определяется как средняя кинетическая энергия, приходящаяся на одну степень свободы системы:

Здесь N - полное число степеней свободы молекулы, - постоянная Больцмана. В случае изолированной системы N=3n-6, поскольку сохраняется ее полный импульс и момент импульса. Кроме того, в этом случае сохраняется полная энергия системы, а температура получается усреднением ее мгновенных значений T(t) по некоторому интервалу времени.

Потенциальная энергия молекулы задается в виде:

+ + + + + +

где слагаемые отвечают следующим типам взаимодействий:

- химическим связям; - валентным углам; - торсионным углам; - плоским группам; - ван-дер-ваальсовым контактам; - электростатике; - водородным связям

Указанные слагаемые имеют различный функциональный вид. Валентные длины поддерживаются за счет потенциала:

где суммирование проводится по всем химическим связям, - обозначение для равновесных валентных длин, r - текущие длины связей, - соответствующие силовые константы Уравнение, описывающее потенциал валентных связей следует из закона Гука

Валентные углы задаются потенциалами

где - равновесные значения углов, - их текущие значения, - силовые константы. Уравнение, описывающее потенциал валентных связей также следует из закона Гука.

Энергия торсионных взаимодействий и потенциалов, отвечающих плоским группам, записываются в одинаковом виде:

где n - кратность торсионного барьера, - сдвиг фазы, константы определяют высоты потенциальных барьеров двугранных углов.

Ван-дер-ваальсовые взаимодействия атомов, разделенных тремя и более валентными связями описываются с помощью потенциалов Леннард-Джонса:

Параметры потенциала A и B зависят от типов атомов i и j, участвующих во взаимодействии; - расстояние между этими атомами.

Очевидно, что вид потенциала зависит от свойств атомов, участвующих в образовании Ван-дер-Ваальсовых связей.

Электростатические взаимодействия задаются кулоновским потенциалом

где , - парциальные заряды на атомах, - диэлектрическая проницаемость среды.

Водородные связи возникают и исчезают в процессе движения атомов между теми из них, которые имеют донорно-акцепторный статус. Функциональный вид потенциала водородной связи аналогичен потенциалу ван-дер-ваальсовым взаимодействий:

Существуют различные наборы параметров для потенциалов взаимодействий. Их значения определяются из учета различных типов экспериментальных данных (спектральные, калориметрические, кристаллографические) и квантовомеханических расчетов.

Последние исследования в этой области показали, что колебания валентных углов, по-видимому, связаны с коллективными движениями в молекуле и, из-за плотной упаковки атомов внутри белка небольшие флуктуации валентных углов ( ) играют существенную роль в движениях, включающих другие степени свободы. При фиксации валентных углов амплитуда флуктуаций торсионных углов уменьшается в 2 раза, а конформационные переходы по торсионным углам из одного минимума энергии в другой исчезают совсем.

В некоторых случаях степени свободы, соответствующие изменениям значений валентных углов, учитываются неявно. Этот учет валентных углов незначительно увеличивает время счета, но значительно увеличивает конформационную подвижность, благодаря чему полученные результаты в большей степени соответствуют экспериментальным данным.

Иногда используют алгоритмы, в которых переменные, соответствующие медленным степеням свободы, постоянны на протяжении некоторого числа шагов. Однако, при таких расчетах происходит довольно быстрое накопление ошибки. Этого недостатка лишены методы MTS (multiple-time-step methods). В этих методах для вычисления сил, соответствующих быстрым и медленным степеням свободы, используются разные временные интервалы. Время счета при этом сокращается в 4 - 5 раз.

Как правило, в методе молекулярной динамике для ускорения расчетов ван-дер-ваальсовые, водородные и электростатические взаимодействия рассчитываются только между атомами, находящимися на расстоянии меньшем, чем радиус обрезания (10-15 ). Недавно был предложен метод PPPC (particle-particle and particle-cell) расчета кулоновских взаимодействий. В методе РРРС каждый атом взаимодействует с ближними атомами путем обычных кулоновских взаимодействий, а с далеко отстоящими ячейками через общий заряд и дипольный момент этой ячейки. Размеры ячеек возрастают как функции расстояния от атома:

Учёт влияния среды в молекулярной динамике.

Проблема учёта влияния кружающей среды на конформационную подвижность изучаемых молекул в рамках метода молекулярной динамики решается двумя способами: явным введением в уравнения движения дополнительного потенциала и явным учётом окружения с добавлением в систему молекул окружающей среды. Добавление в систему дополнительных молекул среды существенно увеличивает объём расчётов (пропорционально квадрату числа атомов) и поэтому мало применимо к расчёту молекулярной динамики биологических макромолекул. Однако существует модификация этого подхода, связанная с вводом периодических граничных условий, что позволяет уменьшить объём системы. Дополнительный потенциал в уравнениях движения не может быть универсальным для разных систем и его тип сильно привязан к конкретной задаче.

Периодические граничные условия.

Периодические граничные условия позволяют рассматривать сравнительно небольшой "кубик" пространства, в котором расположена изучаемая молекула. Молекулы, расположенные внутри кубика со временем претерпевают конформационные движения и перемещаются в пространстве, причём могут пересечь границы кубика. Суть метода заключается в том, что пространство разбивается на одинаковые кубики, причём предполагается, что содержимое кубиков одинаково и границы кубиков соприкасаются. При пересечении молекулой границы одного кубика, она попадает в другой, но это значит, что в первый кубик с противоположной стороны попадает такая же молекула. При этом моделируется динамика лишь одного такого кубика. Естественно, что размер кубика должен быть достаточно большим для исключения возможности краевых эффектов.

Термостаты.

Часто взаимодействие с тепловым резервуаром моделируется дополнительной силой трения . Коэффициент  выбирается таким образом, чтобы сила обеспечила изменение энергии системы по закону

Здесь E – энергия изолированной системы (при отсутствии взаимодействия с резервуаром сохраняется),  E – характерное время взаимодействия с резервуаром, – кинетическая энергия системы, – константа, равная средней кинетической энергии, соответствующей температуре резервуара T0. Уравнения движения метода имеют вид

 =1,…,N

Расчёт траекторий движения в молекулярной динамике по этим уравнениям носит название WCEB (weak coupling to an external bath). Однако, среди специалистов такой метод задания теплового резервуара более известен как метод термостата Берендсена. Этот метод широко применяется для моделирования молекулярной динамики молекул с большим числом степеней свободы, в частности полипептидов и белков.

В броуновской динамике сила, осуществляющая взаимодействие системы с тепловым резервуаром, состоит из двух частей: систематической силы трения FT и шума FC.

i=1,…,3N

 i – коэффициент трения, соответствующий координате Xi. Сила FC(t) –  -коррелированный по времени, гауссовский случайный процесс. Первый и второй моменты этого процесса равны

Интенсивность шума Dij называют тензором диффузии. Уравнения движения метода броуновской динамики называются уравнениями Ланжевена, а метод расчёта молекулярной динамики по этим уравнениям - методом Ланжевеновской динамики

i=1,…,3N

Добавление случайной силы превращает все динамические переменные в случайные величины. Меняется сам способ описания состояния системы: бессмысленно говорить о нахождении системы в точке фазового пространства. Теперь под состоянием системы понимается плотность распределения P(X,V,t) на фазовом пространстве: P(X,V,t) dX dV равно вероятности нахождения системы в момент времени t в малой окрестности точки (X,V). С вышеуказанным уравнением связано дифференциальное уравнение второго порядка для функции P(X,V,t) – уравнение Фоккера-Планка.

Это уравнение описывает эволюцию состояния нашей системы. При выполнении этих условий, любое начальное состояние P0(X,V) будет стремиться к единственному стационарному состоянию P (X,V). Естественно потребовать, чтобы это стационарное состояние совпадало с равновесным ансамблем Гиббса, плотность распределения которого равна

где H(X,V) – гамильтониан системы. Достаточное условие того, чтобы распределение  было стационарным решением уравнения Фоккера-Планка, суть соотношение

Следовательно, интенсивность шума в броуновской динамике не произвольна, а должна выбираться в соответствии с этим уравнением.

В методе Андерсона (МА) взаимодействие системы с тепловым резервуаром моделируется следующим образом. В определённые моменты времени tk движения изолированной системы происходит замена её скоростей V на новые скорости U. Скорости U суть случайные величины, распределённые в соответствии с импульсной частью PM(X,V) равновесного распределения Гиббса. Если на систему не наложено геометрических связей, то PM(X,V) не зависит от X и совпадает с распределением Максвелла.

Если система состоит из одной частицы массы m, то предлагаемая замена скоростей эквивалентна соударению с виртуальным атомом резервуара, который имеет ту же массу m и скорость U, выбранную случайно из распределения Максвелла. Но если частица обладает внутренними степенями свободы, то предлагаемая замена не эквивалентна одному соударению.

Метод Нозе. Во всех предыдущих методах, моделирующих динамику системы в тепловом резервуаре, тепловой резервуар явно не описывался. Обычно вводится явно дополнительная степень свободы s, которая описывает резервуар. Лагранжиан расширенной системы имеет вид

Здесь Q – некоторая константа. Из уравнений Лагранжа получаются уравнения движения как для атомов системы, так и для переменной s. Хотя физический смысл переменной s не ясен, было показано, что такое искусственное расширение системы формально позволяет оценивать средние по Гиббсу величины от функций, определённых на фазовом пространстве системы (без резервуара), путём усреднения их вдоль траекторий расширенной системы.

Часто изучаемые процессы происходят в локализованной области. Например, в ферментативных реакциях, при связывании лигандов для транспортировки. В этих случаях биологическая активность связана с динамикой в окрестности активного центра или места связывания. Для моделирования таких процессов наиболее широко используется метод МД со стохастическими граничными условиями (СГУ). В этом методе изучаемая система разделяется на две области: «реакционную» зону (РЗ) и резервуарную область (РО). Реакционная зона содержит ту часть системы, которая участвует в интересуемом нас процессе. Резервуарная область исключается из вычислений и её влияние на атомы в РЗ учитывается эффективно через среднюю силу и стохастическую силу . Причём стохастическая сила действует только на атомы в окрестности границы РЗ, которая называется буферной зоной.

end;

dxn[0]:=1/(100-vmax)*random(10-5);

dyn[0]:=1/(100-vmax)*random(10-5);

dzn[0]:=1/(100-vmax)*random(10-5);

assign(fx,'c:\x.txt');

assign(fy,'c:\y.txt');

assign(fz,'c:\z.txt');

rewrite(fx);

rewrite(fy);

rewrite(fz);

for i:=0 to 100 do

begin

x0:=x0+dxn[0];

y0:=y0+dyn[0];

z0:=z0+dzn[0];

for j:=1 to k-1 do

begin

n[j].x:=n[j].x+1/(100-vmax)*random(10-5);;

n[j].y:=n[j].y+1/(100-vmax)*random(10-5);;

n[j].z:=n[j].z+1/(100-vmax)*random(10-5);;

if (trunc(x0)=trunc(n[j].x))

and (trunc(y0)=trunc(n[j].y))

and (trunc(z0)=trunc(n[j].z)) then

begin

dxn[0]:=1/(100-vmax)*random(10-5);

dyn[0]:=1/(100-vmax)*random(10-5);

dzn[0]:=1/(100-vmax)*random(10-5);

end;

end;

writeln(i);

writeln(fx,inttostr(trunc(x0)));

writeln(fy,inttostr(trunc(y0)));

writeln(fz,inttostr(trunc(z0)));

end;