Оценка профилей потенциала средней силы методом Umbrella sampling для трансмембранного переноса молекулы воды

Математические заметки СВФУ Июль—сентябрь, 2014. Том 21, №3

УДК 51.76

ОЦЕНКА ПРОФИЛЕЙ ПОТЕНЦИАЛА СРЕДНЕЙ СИЛЫ МЕТОДОМ UMBRELLA SAMPLING ДЛЯ ТРАНСМЕМБРАННОГО ПЕРЕНОСА МОЛЕКУЛЫ ВОДЫ

М. Ю. Антонов, Т. В. Науменкова, А. В. Попинако, И. Н. Николаев, К. В. Шайтан

Аннотация. Моделирование пассивного трансмембранного транспорта малых молекул методами молекулярного моделирования сопряжено с рядом сложностей по той причине, что за доступное время моделирования практически невозможно наблюдать процесс спонтанной диффузии и измерять усредненные макроскопические характеристики процесса. В связи с этим интерес вызывают подходы и методы, используя которые, можно за достижимое время получать результат, позволяющий проводить сравнительное изучение или сравнение с натурными экспериментальными данными. В данной работе с помощью методов управляемой молекулярной динамики и метода Umbrella sampling изучен процесс трансмембранного транспорта молекул воды через липидный бислой. Использован бислой 1,2-димиристоил^п-глицеро-3-фосфатидилхолина (ДМФХ) в силовом поле OPLS-AA. Профили средней силы построены из данных расчетов с длинами траекторий 6—10 нс. Произведено сравнение получаемых профилей как при различном пространственном положении пробной (проникающей через бислой) молекулы воды, так и при разном шаге выбора стартовых конфигураций, а также при различных начальных значениях скоростей атомов в системе. Общее время моделирования составило более 10 мс. Показана значительная зависимость результатов расчетов от начальных данных и протокола вычислительного эксперимента. В рамках подхода произведены оценка кинетических параметров проникновения и оценка профиля потенциала средней силы, а также изучена зависимость вычисляемых значений от начальных данных эксперимента. Оценена погрешность вычислений на основе анализа данных ряда численных экспериментов.

Ключевые слова: молекулярное моделирование, метод зонтичной выборки, биомембраны.

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 12—04—31934\13), Научно-образовательного фонда поддержки молодых ученых Республики Саха (Якутия) (грант 2014—01—0008) и Российского научного фонда (код проекта 14-24-00172).

© 2014 Антонов М. Ю., Науменкова Т. В., Попинако А. В., Николаев И. Н., Шайтан К. В.

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

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

Достаточно широко распространенной методикой является изучение трансмембранного транспорта с использованием методов управляемой молекулярной динамики (УМД). В рамках данного подхода к системе прикладывается дополнительный потенциал, позволяющий стимулировать систему по выбранным степеням свободы. В частности, метод позволяет получать значения «эффективной» микровязкости бислоя, а также оценку коэффициента диффузии, которые можно использовать для сравнительного анализа. Недостатком этого подхода является наблюдаемая зависимость результатов от протокола вычислительного эксперимента, а также отсутствие оценки коэффициента разделения растворитель-мембрана [4,5].

Другим подходом является использование метода зонтичной выборки (Umbrella Sampling, US). Использование этого метода позволяет проводить статистический анализ высокоэнергетических областей конфигурационного пространства, которые практически не заполнены в равновесных МД экспериментах. В данной работе с помощью методов управляемой молекулярной динамики и метода US изучен процесс трансмембранного транспорта молекул воды через липидный бислой. Использован бислой 1,2-димиристоил-8п-глицеро-3-фосфатидилхолина (ДМФХ) в силовом поле OPLS-AA. Произведено сравнение получаемых профилей как при различном пространственном положении пробной (проникающей через бислой) молекулы воды, так и при разном шаге выбора стартовых конфигураций, а также при различных начальных значениях скоростей атомов в системе.

1. Материалы и методы

Суть метода молекулярной динамики сводится к представлению изучаемой системы как множества материальных точек, взаимодействие которых описывается классическими законами механики. Каждая точка при этом представляет собой атом или группу атомов изучаемой системы. На рис. 1 показан вид изучаемой системы, а также вид молекулы ДМФХ и ее химическое строение.

Взаимодействие в системе описывается в виде суммы частичных потенци-

Рис. 1. Вид изучаемой системы вид и строение молекулы ДМФХ (Ь).

алов:

и

+ £ 3 (ф<>з->к>,) + £ (гг>з) + 3 (Г4,3-

г,з,к,1 г=з

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

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

величина константы жесткости выбирается таким образом, чтобы воспроизводить геометрию и колебательные частоты модельных соединений.

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

3 (Фг,з,м) = £ 4Ып(1 + осз(пф - ф0)).

В методе МД электроны не учитываются явно, каждый атом обладает парциальным зарядом, а кулоновская энергия взаимодействия между атомами опи сывается стандартным выражением:

№з

Цс (Гг,з ) =

4пеоГгз'

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

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

ЧШ -ШУ <«

В целом общий вид функции потенциальной энергии практически не изменяется в разных вариантах силовых полей, однако параметризация сил может различаться. В данной работе использовалось силовое поле ОРЬ8-ЛЛ [6, 7].

Таким образом, движение атомов в потенциальном поле описывается системой обыкновенных дифференциальных уравнений второго порядка:

= — (6) \дхдуг дгг)'

Данная система уравнений решается численно. Могут применяться различные численные методы решения уравнений, различающиеся точностью и вычислительной сложностью. Компромиссом между точностью процедуры и скоростью вычислений является широко используемый в МД метод Верле [1] и его варианты, при этом значения координат атомов на новом временном слое вычисляются по следующим формулам:

п{1, + Аг) = 2п(г) - п{1, - М) +

тг

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

Исследование систем в ансамблях с постоянной энергией, как правило, не представляет значительного интереса, поскольку большинство экспериментальных данных получено в условиях постоянства температуры и/или давления. Поэтому для поддержания условий постоянства температуры и давления в МД экспериментах используются специальные алгоритмы: термостаты и баростаты. В данной работе использовался термостат стохастической динамики. Его суть заключается в том, что к атомам системы прилагается дополнительная случайная (стохастическая) сила, описывающая взаимодействие системы с внешней средой заданной температуры [8]. Для контроля условий постоянства давления использовался баростат Берендсена [9]. В рамках данного подхода объем системы может изменяться: к основному уравнению движения добавляются уравнения движения границ периодической ячейки.

В экспериментальных работах при исследовании пассивного трансмембранного транспорта малых молекул удается измерить коэффициент проницаемости мембраны Р, а также коэффициент распределения вода-гидрофобная среда Кр [2, 3,10,11]. Эти величины связаны соотношением

6

где D — коэффициент диффузии в мембране, Kp — коэффициент межфазного разделения вещества, Дх — толщина мембраны. В МД удается различными способами прямо измерить коэффициент диффузии D [4], однако измерить прямыми способами коэффициент межфазного распределения Kp не удается.

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

G

где p(£) — плотность вероятности нахождения системы в состоянии c заданным значением При моделировании трансмембранного транспорта в качестве как правило, выбирается координата z центра массы малой молекулы, где ось Oz перпендикулярна латеральной плоскости мембраны (рис. 2a).

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

Выходом здесь является метод зонтичной выборки, позволяющий за достижимое время моделирования оценить значение коэффициента распределения вода-мембрана. В рамках данного подхода к системе прилагается дополнительный потенциал, способный удержать исследуемую молекулу внутри потенциально невыгодной области [12-14] (рис. 2b). Модификация исходной функции потенциальной энергии заключается в добавлении внешнего потенциала — как правило, используется гармонический потенциал, приложенный к координате реакции:

U

Поскольку вид функции W (r) известен, значение свободной энергии G для невозмущенной системы может быть получено по формуле

G

где p(£) — плотность вероятности нахождения системы в состоянии c заданным значением £ в возмущенной системе.

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

Исходный вид профиля свободной энергии может быть восстановлен из множества полученных распределений вокруг каждой начальной точки одним из методов. В данной работе использовался метод анализа взвешенных гистограмм (WHAM, [15]).

Рис. 2. Ось координаты реакции (а), вид гармонических потенциалов для удержания молекулы воды внутри мембраны (области высоких значений свободной энергии) (Ь).

2. Протокол молекулярной динамики

В работе использовалась модель воды Т1Р4Р/2005, которая достаточно хорошо описывает жидкое и кристаллическое состояние, диффузные и вязкостные свойства, а также эффекты поверхностного натяжения [16-18]. Парциальные заряды атомов в липидах оценивались с использованием неограниченного метода Хартри — Фока, базиса 6-3Ю* и метода Малликена. Радиус обрезания невалентных взаимодействий составлял 1.8 нм. Вычисления проводились в NPT-ансамбле, для поддержания условий постоянства температуры и давления использовался термостат стохастической динамики и баростат Берендсена [9]. Для учета эффектов поверхностного натяжения бислоя величина давления в латеральной плоскости выбиралась отрицательной, достаточной для поддержания величины удельной площади на липид в известных из эксперимента пределах [19]. Давление баростата вдоль нормали мембраны составляло 1 бар, перпендикулярно нормали мембраны — 50 бар. В качестве модельной мембраны эукариотической клетки использовался бислой из 80 молекул ДМФХ. Для построения бислоя использовались по 5 конформаций молекул каждого типа. Молекулы липидов помещались в узлы гексагональной решетки и поворачивались на случайный угол вокруг своей оси. Площадь поверхности на одну молекулу липида составляла 0.6 нм2. Для сборки мембранных систем использовалось оригинальное программное обеспечение Б1Ьауег.

Вычислительные эксперименты проводились в программном пакете Сгошася 4.6.7 [20]. Перед проведением процедуры ИЯ производилась релаксация модельной мембраны в течение 20 нс в NPT ансамбле при температуре 300К. Начальные скорости атомов определялись с помощью генератора случайных

чисел по распределению Максвелла. Шаг интегрирования составлял 1 фс. Для процедуры зонтичной выборки стартовые конфигурации для молекулы воды выбирались с шагом 0.05 и 0.1 нм вдоль оси Ог (координаты реакции). Использовался гармонический потенциал с силовой константой 1000 кДж/(моль-нм2). Длины траекторий для набора статистики составляли 6-10 нс. Проводился анализ зависимости вычисляемых профилей энергии от начальных данных, координат, длин траекторий, а также зависимость от начальных скоростей атомов в системе.

Обработка и анализ траекторий осуществлялись с помощью встроенных модулей Сгошаев [20], для визуализации использовался пакет УМБ [21].

3. Результаты

В табл. 1 приведены основные структурные характеристики бислоев после релаксации в течение 20 нс. Усреднение значений проводилось по последним 2 нс траекторий. Для сравнения приведены данные, полученные экспериментально (Э) [22].

Таблица 1. Основные структурные характеристики моделируемых бислоев

Параметр Вислой ДМФХ

Э мд

Толщина бислоя, Е 33.8 37.1±0.2

Площадь поверхности на липид, Е2 65.4 55.2±0.2

Параметр порядка 0.184 0.218

Как видно, использованный протокол МД позволяет моделировать мембранные системы, структурные характеристики которых близки к экспериментально наблюдаемым. Полученные мембраны и разработанный протокол молекулярной динамики использовались для дальнейших вычислений.

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

Была исследована зависимость вычисляемых результатов от положений начальных координат молекул воды в толще липидного бислоя. Для этого было выбрано 4 возможных траектории прохождения молекулы воды через липид-ный бислой, и начальные точки были выбраны на этих траекториях с шагом 0.1 нм, длины траекторий составляли 6-10 нс. Было проведено 9 расчетов с разными начальными значениями скоростей атомов в системе (рис. 3). Рассчитанные высоты потенциального барьера, приведенные в табл. 2, составили для разных расчетов от 20 до 36 кДж/моль, среднее рассчитанное значение — 27 кДж/моль, среднеквадратическое отклонение — 5.5 кДж/моль. Известные экспериментальные оценки коэффициента межфазного распределения

Кр вода-гидрофобный растворитель для воды лежат в пределах от 4.2•Ю-5 (вода-гексадекан [2]) до 1.4•Ю-3 (вода-оливковое масло [2]). Полученная нами средняя величина высоты барьера 27 кДж/моль по 9 измерениям соответству-

т-^ша-Ь/тетЬгапе п 1Л—4

ет расчетному значению Кр ' =1.6 • 104 и не противоречит экспе-

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

Таблица 2. Рассчитанная высота барьера в экспериментах при шаге 0.1 нм и К^агт = 1000 кДж/

Номер вычислительного эксперимента 1 2 3 4 5 6 7 8 9

Величина барьера, кДж/моль 23.5 19.6 25.0 22.8 32.0 27.7 35.5 25.0 36.0

На рис. 4 приведен вид гистограмм, характеризующих статистическое покрытие конформационного пространства вдоль координаты реакции в вычислительном эксперименте. Можно отметить, что при выборе шага 0.1 нм и величине

Рис. 4. Вид гистограммы для одного из экспериментов с шагом 0.1 нм.

Заметны области низкой заселенности

силовой константы К^агт гармонического потенциала 1000 кДж/(моль-нм2) в гистограмме присутствуют области с низкой заселенностью.

Для исследования влияния качества выборки конформационного пространства проведены дополнительные вычислительные эксперименты с шагом 0.05 нм, величина силовой константы К^агт не менялась. Было проведено 4 измерения с различными начальными скоростями атомов. Как видно из табл. 3, уменьшение величины шага до 0.05 нм не приводит к значительному изменению величины среднеквадратического отклонения (рассчитанная величина 4.64 кДж/моль). Исследование гистограмм, характеризующих статистическое покрытие конформационного пространства вдоль координаты реакции (не показано), не указывает на наличие областей, слабо представленных в итоговой выборке.

Таблица 3. Рассчитанная высота барьера в экспериментах при шаге 0.05 нм и К^агт = 1000 кДж/

Номер вычислительного эксперимента 1 2 3 4

Величина барьера, кДж/моль 32.0 19.5 28.0 24.0

4. Выводы

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

Показано, что метод Umbrella Sampling позволяет проводить исследование трансмембранного транспорта молекул воды через бислой ДМФХ и произво-

дить оценку профилей средней силы. При этом вычисленная средняя величина энергетического барьера порядка 27 кДж/моль лежит в пределах допустимых значений, известных из натурных экспериментов по измерению Кр^а1 межфазного распределения вода-гидрофобный растворитель.

Сэмплирование с шагом 0.1 нм и длинами траекторий 6-10 нс может приводить к появлению областей конформационного пространства с низкой заселенностью вдоль координаты реакции, при шаге 0.05 нм такого не наблюдалось. Однако изменение шага незначительно повлияло на измеряемое значение среднеквадратического отклонения, что может указывать в данном случае на стохастичность процесса трансмембранного транспорта.

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

Вычислительные расчеты выполнены с использованием ресурсов суперкомпьютерного комплекса МГУ им. М. В. Ломоносова [23] и суперкомпьютерного комплекса СВФУ «Ариан Кузьмин».

ЛИТЕРАТУРА

1. Verlet L. Computer experiments on classical fluids. I. Thermodynamical properties of Len-nard-Jones molecules // Phys. Rev. 1967. V. 159, N 1. P. 98-103.

2. Walter A., Gutknecht J. Permeability of small nonelectrolytes through lipid bilayer membranes // J. Membr. Biol. 1986. V. 90, N 3. P. 207-217.

3. Diamond J. M., Katz Y. Interpretation of nonelectrolyte partition coefficients between dimyri-stoyl lecithin and water // J. Membr. Biol. 1974. V. 17, N 2. P. 121-154.

4. Shaitan K. V., Antonov M. Y., Tourleigh Y. V., Levtsova O. V., Tereshkina K. B., Nikolaev I. N., Kirpichnikov M. P. Comparative study of molecular dynamics, diffusion, and permeability for ligands in biomembranes of different lipid composition // Biochem. Suppl. Ser. A Membr. Cell Biol. 2008. N 2. P. 73-81.

5. Антонов М. Ю., Науменкова Т. В., Левцова О. В., Шайтан К. В. Исследование динамики и трансмембранной диффузии методами компьютерного моделирования // Суперкомпьютерные технологии математического моделирования: тр. II Междунар. конф. 2014. С. 16-29.

6. Jorgensen W. L., Chandrasekhar J., Madura J. D., Impey R. W., Klein M. L. Comparison of simple potential functions for simulating liquid water //J. Chem. Phys. 1983. V. 79, N 2. P. 926.

7. Jorgensen W. L., Maxwell D. S., Tirado-Rives J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids //J. Amer. Chem. Soc. 1996. V. 118, N 45. P. 11225-11236.

8. Holm C., Kremer K. Advanced computer simulation. V. 173. Berlin; Heidelberg: Springer-Verl., 2005.

9. Berendsen H. J. C., Postma J. P. M., van Gunsteren W. F., DiNola A., Haak J. R. Molecular dynamics with coupling to an external bath // J. Chem. Phys. 1984. V. 81, N 8. P. 3684.

10. Olbrich K., Rawicz W., Needham D., Evans E. Water permeability and mechanical strength of polyunsaturated lipid bilayers // Biophys. J. 2000. V. 79, N 1. P. 321-327.

11. Bean R. C., Shepherd W. C., Chan H. Permeability of lipid bilayer membranes to organic solutes // J. Gen. Physiol. 1968. V. 52, N 3. P. 495-508.

12. Buch I., Sadiq S. K., de Fabritiis G. Optimized potential of mean force calculations for standard binding free energies // J. Chem. Theory Comput. 2011. V. 7, N 6. P. 1765-1772.

13. Kastner J. Umbrella sampling // Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011. V. 1, N 6. P. 932-942.

14. Torrie G. M., Valleau J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling // J. Comput. Phys. 1977. V. 23, N 2. P. 187-199.

15. Kumar S., Rosenberg J. M., Bouzida D., Swendsen R. H., Kollman P. A. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method // J. Comput. Chem. 1992. V. 13, N 8. P. 1011-1021.

16. Vega C., de Miguel E. Surface tension of the most popular models of water by using the test-area simulation method // J. Chem. Phys. 2007. V. 126, N 15. P. 154707.

17. Abascal J. L. F., Vega C. A general purpose model for the condensed phases of water: TIP4P/2005 // J. Chem. Phys. 2005. V. 123, N 23. P. 234-505.

18. Conde M. M., Gonzalez M. A., Abascal J. L. F., Vega C. Determining the phase diagram of water from direct coexistence simulations: the phase diagram of the TIP4P/2005 model revisited // J. Chem. Phys. 2013. V. 139, N 15. P. 154-505.

19. Chiu S. W., Clark M., Balaji V., Subramaniam S., Scott H. L., Jakobsson E. Incorporation of surface tension into molecular dynamics simulation of an interface: a fluid phase lipid bilayer membrane // Biophys. J. 1995. V. 69, N 4. P. 1230-1245.

20. Hess B., Kutzner C, van der Spoel D., Lindahl E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation //J. Chem. Theory Comput. 2008. V. 4, N 3. P. 435-447.

21. Humphrey W., Dalke A., Schulten K. VMD: visual molecular dynamics // J. Mol. Graph. 1996. V. 14, N 1. P. 33-38.

22. Petrache H. I., Dodd S. W., Brown M. F. Area per lipid and acyl length distributions in fluid phosphatidylcholines determined by (2)H NMR spectroscopy // Biophys. J. 2000. V. 79, N 6. P. 3172-3192.

23. Sadovnichy V., Tikhonravov A., Voevodin V., Opanasenko V. "Lomonosov": Supercomputing at Moscow State University // Contemporary high performance computing: From petascale toward exascale. 2013. P. 283-307.

Статья поступила 13 ноября 2014 г.

Антонов Михаил Юрьевич, Николаев Иван Никитич Северо-Восточный федеральный университет, ул. Белинского, 58, Якутск 677000, Республика Саха (Якутия) mikhail@s-vfu.ru, n_ivan_n@mail .ru

Науменкова Татьяна Витальевна, Шайтан Константин Вольдемарович Московский гос. университет им. М. В. Ломоносова Москва, 119991

tnaumenkova@gmail.com, shaytan49@yandex.ru

Попинако Анна Владимировна

Институт биохимии им. А. Н. Баха Российской академии наук, Ленинский пр., 33, Москва 119071 popinakoav@gmail.com