1. Введение

За последние 30 лет наиболее заметным достижением в современной компьютерной химии и биологии стало развитие вычислительных методов изучения динамического поведения молекулярных систем. По словам знаменитого физика Ричарда Фейнмана, для понимания принципов функционирования живой материи достаточно предположения, «что все состоит из атомов, и что все, что происходит в живом объекте, можно представить в терминах колебаний и смещений атомов» (курсив автора). Другими словами, знание динамического поведения молекулярной системы позволяет пролить свет на механизм её функционирования. С точки зрения физики, наиболее корректным описанием эволюции молекулярной системы во времени является квантово-механическое представление. Однако получить решение уравнения Шредингера из первых принципов (ab initio), даже учитывая современный уровень компьютерной техники, практически возможно только для атома гелия, но никак не более крупных систем. Существуют, правда, подходы, использующие различные приближения и так называемые полуэмирические методы квантовой химии, также весьма трудоемкие в случае больших систем атомов.

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

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

Сегодня использование МД в исследовании различных биомолекулярных систем является обычной практикой. Этот подход широко применяется для изучения структурно-динамичеких свойств белков, нуклеиновых кислот, липидных мембран. Однако всего 30 лет назад публикация в Nature кратковременного расчета МД бычьего панкреатического трипсинового ингибитора (БПТИ) [1] стала настоящим прорывом в структурной биологии и во многом определила развитие методов молекулярного моделирования. Впрочем, история началась задолго до этого эпохального исследования.

2. Молекулярная динамика биллиардных шаров

Рождение методов молекулярной механики, включающих МД и метод конформационного поиска Монте-Карло (МК), было абсолютно не связано с биологией и происходило в недрах американского ядерного проекта — национальных лабораториях Лос-Аламоса и Ливермора. Первая лаборатория основана в 1943 г. под руководством Роберта Оппенгеймера (Robert Oppenheimer) и стала известна всему миру в связи с проектом «Манхэттен» и последующей бомбардировкой Хиросимы и Нагасаки. Возникновение второй лаборатории в 1952 г. явилось следствием эскалации «холодной» войны между США и СССР. Именно в этих лабораториях ЭВМ стали впервые в мире интенсивно применяться в научных расчетах.

Одним из создателей методов молекулярной механики по праву можно считать Бэрни Алдера (Bernie Alder) (род. 1925). Работая в начале 50-х аспирантом у известного физика Джона Кирквуда (John Kirkwood) в Калифорнийском техническом институте (CALTECH), он занимался проблемой поведения системы твердых сфер, в которой, подобно биллиардным шарам, частицы взаимодействуют непосредственно при столкновении и двигаются между соударениями как свободные частицы. Возник вопрос — может ли в данной системе наблюдаться фазовый переход из жидкой фазы в твердую? Данная задача не поддавалась простому аналитическому решению, ввиду сложности и нелинейности получаемых интегральных уравнений. Тогда вместе со своим непосредственным руководителем Стэном Фрэнклом (Stan Frankel), который перешел в Калифорнийский институт из Лос-Аламоса, они решили использовать стохастический подход в решении данной задачи. Для этого ими был разработан алгоритм, который позволял исследовать поведение системы твердых сфер путем серии случайных «перетасовок» частиц в некоторой «виртуальной ячейке». Изначально в распоряжении Алдера и Фрэнкла был только механический компьютер IBM, который не позволял далеко продвинуться в расчетах. Однако поскольку Фрэнкл был широко известен в «компьютерных» кругах того времени, для него не составило труда провести расчеты в Англии на первом в мире ЭВМ FERRANTI. Это было летом 1950 г., ещё до того, как ЭВМ стали использовать в Лос-Аламосе. Тем не менее, Алдер так и не увидел фазового перехода в исследуемой им системе.

Несмотря на это, созданный ими алгоритм получил в последствие широкое распространение, и известен миру под названием «метод Монте-Карло» (МК). Правда, автором МК считается Николас Метрополис (рис. 1), который параллельно с Алдером работал вместе с Эдвардом Теллером и Маршалом Розенблатом над сходной проблемой в Лос-Аламосе. По поводу истории с авторством Алдер в одном из своих интервью рассказывает, что они вместе с Фрэнклом разработали этот метод независимо от Метрополиса и даже раньше него провели расчеты на ЭВМ. Однако их начальник Джон Кирквуд скептически отнесся к таким исследованиям, что не позволило опубликовать полученные результаты. Метрополис и соавторы опубликовали работу в 1953 г. [2], и вся последующая известность досталась именно им. Это достаточно типичное явление в науке, когда одно и то же открытие делается сразу несколькими независимыми исследователями, но приоритет остается за первым опубликовавшим. Стоит отметить, что Метрополис указал в сносках к своей знаменитой статье, что также над этой проблемой независимо работают Алдер и Фрэнкл (рис. 1.)

После основания Ливерморской лаборатории в 1952 г. Алдер перешел в эту организацию, работая также по совместительству в университете Беркли (Berkeley University of California). Ему пришлось сменить область научных интересов и даже заняться экспериментальной работой на станции по изучению взрывов. Тем не менее, проблема фазового перехода в системе тверды сфер и неудачи с использованием МК в ее исследовании не давали Алдеру покоя. В 1955 году он вновь занялся этой задачей. Совместно с Томом Вайнрайтом (Tom Wainwright) он использовал динамический подход в исследовании фазовой диаграммы системы твердых сфер (рис. 2). Расчеты МД даже для простейшей системы в то время были значительно более трудоемкими, чем стохастический поиск методом МК.

Поначалу многие коллеги отнеслись скептически к их затее. Группой Алдера была разработана первая программа для расчета МД — STEP. Были потрачены тысячи часов расчетов с использованием STEP для накопления необходимой статистики. Вычисления проводились на компьютерах UNIVAC и на IBM 704 (рис. 3). Поскольку компьютерное время было крайне ограниченно, Алдеру приходилось даже обращаться к своим знакомым, имеющим доступ к ЭВМ, с просьбой использовать все свободное машинное время для расчетов STEP.

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

3. От биллиардных шаров к реальным атомам

Дальнейшее развитие метода МД шло по нескольким направлениям. Одно из них было связано с изучением поведения не абстрактных сфер, а реалистических атомарных систем. Пионерской работой в данной области стали расчеты жидкого аргона, проведенные Анесуром Рахманом (Aneesur Rahman, рис. 5) [4]. В системе, состоящей из 864 атомов, взаимодействие частиц определялось не только их столкновениями, но также дополнительным парным потенциалом Ленарда-Джонса, соответствующим ван-дер-ваальсовый энергии взаимодействия атомов. Поскольку Рахман работал в американской Национальной Аргоннской лаборатории, занимающейся стратегическими вопросами ядерной энергетики, в его доступе были мощные вычислительные ресурсы. Стоит отметить, что в отличие от лабораторий в Лос-Аламосе и Ливерморе, Аргоннская лаборатория никогда не занималась разработкой оружия и исследовала «мирный атом», в частности, проектировала ядерные реакторы электростанций. Расчеты жидкого аргона были проведены с использованием новейшего компьютера CDC 3600, уже к тому моменту имевшего компилятор для языка программирования Фортран-66. Общая длительность МД-траектории (последовательного набора координат и скоростей системы), полученного Рахманом составила ~1 пс (10⁻¹² с).

Исследования ван-дер-ваальсовых жидкостей с помощью метода МД нашли продолжения в работах французского физика Лу Верле (Loup Verlet, рис. 6) из Университета Пари-Сюд (L’Université Paris-Sud). В 1967 г. Верле рассчитал фазовую диаграмму для жидкого аргона, а также на основании компьютерного эксперимента определил для такой жидкости ряд характерных корреляционных функций, которые давали хорошее согласие с разработанной ранее теорией жидкого состояния [5]. Кроме того, заслугой Верле является разработка эффективного метода численного интегрирования уравнений движения (алгоритм Верле), согласно которому находят координаты каждого атома на основании значения координат в два предшествующих момента времени. Это позволяет не проводить решение дифференциальных уравнений для вычисления скоростей на каждом шаге МД, а получать их на основании разностного уравнения, также предложенного Верле. Такой алгоритм помогает существенно экономить компьютерное время при расчете МД.

4. Метод эмпирических силовых полей. Моделирование биомолекул

Существенный прорыв в развитии вычислительных подходов в изучении реальных молекул стал возможен после того, как израильским физиком Шнеиром Лифсоном (Shneior Lifson) из Вейсманского института в середине 60-х была предложена концепция самосогласованного силового поля (CFF — consistent force field).

Идея заключалась в использовании набора простых термов потенциальной энергии для расчета различных свойств молекул. Таким образом, потенциальная энергия молекулы описывается набором уравнений (рис. 8), параметрами которых служат константы (длины и энергии межатомных связей, значения валентных и двугранных углов и др.), полученные экспериментально или с помощью методов квантовой химии. Нобелевский лауреат, английский биохимик Джон Кендрю (John Kendrew [6]), узнав об идее Лифсона, увидел большие возможности использования CFF в изучении молекул белков и нуклеиновых кислот. В 1967 г. он направил одного из своих аспирантов работать в группу Лифсона. Этим аспирантом был Майкл Левит (Michael Levitt, рис. 9). Он начал работу по созданию компьютерной программы, в которой была бы реализована идея Лифсона. Вейсманский институт в то время обладал мощными вычислительными ресурсами. В частности, в 1963 г. была введена в пользование ЭВМ GOLEM I, разработанная и собранная местными инженерами. Особенность этого компьютера, названного по мотивам еврейского фольклора, заключалась в возможности обработки слов размером 75 бит (мировой рекорд того времени).

За два месяца Левиту удалось написать программу, также названную CFF, которая позволяла рассчитывать энергию и силы в любой молекулярной системе. Следуя наставлениям Кендрю, Левит провел первые в мире расчеты потенциальной энергии белков. Для этого были использованы уже полученные к тому времени с помощью рентгеноструктурного анализа пространственные структуры миоглобина и лизоцима. Вернувшись в Кэмбридж в 1968 году, Левит продолжил свою работу в Лаборатории молекулярной биологии, в которой тогда одновременно работало три нобелевских лауреата — Френсис Крик (Francis Crick), Фрэд Сэнгер (Fred Sanger) и Кендрю. Дальнейшая работа Левита была связана с построением модели пространственной структуры транспортной РНК (тРНК). Существование этой молекулы было предсказано Криком за десять лет до того, как начали проводиться активные исследования ее структуры и функций. тРНК содержит ~2000 атомов, поэтому, когда Левит у себя дома построил CPK-модель из конструктора [6], вес «молекулы» составил почти 50 килограмм. Транспортировать такое сооружение, по словам Левита, было нелегко — ему даже пришлось воспользоваться окном, чтобы спустить «молекулу» из своего домашнего кабинета на первый этаж. Позднее в лаборатории Левит заново построил тРНК уже в виде проволочной модели [6], которую использовал для определения взаимного расположения всех атомов в молекуле. Оптимизированная с помощью программы CFF, модель тРНК была вскоре опубликована [7]. Однако структура нуклеиновых кислот в то время казалась достаточно простой и очевидной, и поэтому более интригующей тематикой были принципы организации полипептидной цепи. В дальнейшем Левит активно работал над проблемой сворачивания белковых молекул и механизмов ферментативных реакций, используя различные методы конформационного анализа, основанные на концепции силового поля. Его пионерские работы предваряли множество последующих исследований в области структурной биологии, активно проводимые уже на протяжении четверти века.

Программа CFF, написанная Левитом ещё в годы аспирантуры, оказалась «впутанной» в историю развития методов МД. Как он писал позднее, CFF «путешествовала по миру самостоятельно», появляясь в разных лабораториях [8]. В 1971 году в лаборатории Мартина Карплюса (Martin Karplus, рис. 10) в Гарварде, Брюс Джелин (Bruce Gelin), только что вернувшийся со службы в американской армии, начал работу по переписыванию кода CFF. Впоследствии Левит, увидевший новую версию своей программы, был поражен элегантностью написанного Джелиным кода: «Этот вариант сильно отличался: пока я осваивал программирование по стандартному руководству IBM FORTRAN II, Брюс Джелин сильно преуспел в этом вопросе. Я до сих пор помню свое восхищение, когда увидел эту версию своей программы — многие названия переменных остались от меня, однако все в целом выглядело намного элегантнее!» [8]. Программа, написанная Джелином, позволила провести первый расчет молекулярной динамики белка, который был опубликован Эндрю Мак-Кэммоном (Andrew McCammon, рис. 10) и Карплюсом в 1977 г. [1]. Работа стала настоящим прорывом, свидетельствуя о наступлении новой эры в компьютерных исследованиях молекулярных систем. С этого момента использование силовых полей для расчетов МД биомолекул стало стандартной практикой.

Как ни странно, свои вычисления Мак-Кэммон и Карплюс провели не в Америке, а во французском международном вычислительном центре CECAM (Centre Européen de Calcul Atomique et Moléculaire), организованном в конце 60-х при университете Пари-Сюд, где в то время активно занимались изучением динамики жидкостей (в частности, Верле — см. выше). Работа в этом центре была организована по весьма оригинальной схеме — на несколько месяцев собиралась рабочая группа ученых из разных стран для решения какой-нибудь конкретной задачи.

Код, написанный Брюсом Джелином, впоследствии лег в основу нескольких МД-программ следующего поколения, которые активно используются и по сей день: CHARMM (Chemistry at HARvard Molecular Mechanics) — разработка группы Карплюса, AMBER — группа Питера Коллманна (Peter Kollmann) из Калифорниского университета (UCSF) и Discover — компания-производитель научного программного обеспечения Biosym (сейчас — Accelrys), организованная Арнольдом Хаглером (Arnold Hagler).

5. Динамика белковой глобулы

В чем же все-таки заключалась революционность работы Мак-Кэммона и Карплюса? Вряд ли современного читателя впечатлят полученные ими результаты и используемые методы расчетов. Объектом исследования авторы выбрали небольшой белок БПТИ, состоящий из 58 аминокислотных остатков (рис. 11) и содержащий три дисульфидные связи, которые стабилизируют молекулу. В качестве стартовой конформации авторами была взята кристаллографическая структура БПТИ, полученная несколькими годами ранее.

Выбранная ими схема расчета МД в вакууме, в целом напоминает современные методики. Шаг интегрирования динамических уравнений был выбран равным 1 фс (10⁻¹⁵ с) — заведомо меньше самых быстрых движений в белке (колебания протонов). Поскольку исходная конформация БПТИ не находилась в энергетическом минимуме, чтобы понизить потенциальную энергию, часть ее была переведена в кинетическую. В результате белок за несколько сотен шагов МД «нагрелся» от 0 до 285 K. Последующие 9 пс (9000 шагов) основного расчета температура в системе колебалась в районе 295 K (22° C — комнатная температура). Стоит отметить, что «температура» в МД является достаточно абстрактным понятием и в соответствии с молекулярно-кинетической теорией определяется кинетической энергией атомов (их скоростями). В современных расчетах используют специальные алгоритмы — «термостаты», которые позволяют усреднять температуру в системе около заданного значения, избегая появления «нагретых» и «замороженных» областей. Их принцип основан на введении в уравнения МД дополнительной знакопеременной силы трения, зависящей от скорости атома. Эта сила может либо разгонять «замороженные» атомы, либо тормозить «нагретые». Один из первых подобных алгоритмов, позволяющих проводить расчеты в реалистичных условиях с постоянной температурой несколькими годами позже был предложен Германом Берендсеном (Herman J.C. Berendsen, рис. 12) [9].

Несмотря на то, что расчеты БПТИ Мак-Кэммона и Карплюса были весьма кратковременными и проводились в нереалистичных для белковой глобулы условиях (отсутствие растворителя, непостоянная температура и пр.), ряд выводов, сделанных авторами на основании свой работы, имел большое значение для развития теоретических представлений о принципах организации молекулы белка. В частности, анализ динамического поведения атомов БПТИ в процессе МД позволил им установить, что внутренняя организация белковой глобулы по своим динамическим свойствам ближе к жидкости, чем к твердому телу. При этом ими было постулировано, что внутримолекулярные движения в белке имеют диффузионный характер, или как более точно было сформулировано позднее — характер ограниченной диффузии. Кроме того, сам факт, что с помощью компьютерных расчетов можно наблюдать за структурно-динамическими свойствами биомолекул на атомном уровне, подвиг многих ученых на активные исследования в данной области. Количество работ по МД белков и, чуть позже — ДНК, стремительно выросло за несколько следующих лет. Уже через два года Мак-Кэммон и Карплюс опубликовали статью, в которой длина траектории МД для БПТИ была на порядок больше (100 пс) [10]. В те годы ими также были исследованы вопросы об учёте растворителя в расчетах [11], характере связывания кислорода с молекулами гемоглобина и миоглобина [12].

6. Заключение

В настоящее время использование методов МД в исследовании различных биомолекулярных систем стало обычным явлением. Ежегодно выходят тысячи публикаций, посвященных данной тематике. Длительность расчетов часто достигает несколько микросекунд (10⁻⁶ с), а размеры моделируемой системы — нескольких миллионов атомов. При этом объектом исследования являются уже не отдельные молекулы, а сложные молекулярные комплексы, например, белки, встроенные в многокомпонентные мембраны в ячейках с молекулами воды и различными ионами [13], и даже такие сложные молекулярные машины, как рибосома [14] (рис. 13). Безусловно, наблюдаемый прогресс в области моделирования биологических систем, по сравнению с работами тридцатилетней (и тем более — полувековой) давности, стал возможен благодаря стремительной эволюции компьютерной техники. Однако стоит отметить, что основные концепции и практически весь математический аппарат МД были разработаны в те далекие годы, когда исследователи даже не представляли фантастических вычислительных возможностей будущего.

Литература

1. J. Andrew McCammon, Bruce R. Gelin, Martin Karplus. (1977). Dynamics of folded proteins. Nature. 267, 585-590;
2. Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, Edward Teller. (1953). Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics. 21, 1087-1092;
3. B. J. Alder, T. E. Wainwright. (1957). Phase Transition for a Hard Sphere System. The Journal of Chemical Physics. 27, 1208-1209;
4. A. Rahman. (1964). Correlations in the Motion of Atoms in Liquid Argon. Phys. Rev.. 136, A405-A411;
5. Loup Verlet. (1967). Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev.. 159, 98-103;
6. На заре молекулярной графики;
7. MICHAEL LEVITT. (1969). Detailed Molecular Model for Transfer Ribonucleic Acid. Nature. 224, 759-763;
8. Michael Levitt. (2001). . Nat. Struct Biol.. 8, 392-393;
9. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak. (1984). Molecular dynamics with coupling to an external bath. The Journal of Chemical Physics. 81, 3684-3690;
10. M. KARPLUS, J. A. MCCAMMON. (1979). Protein structural fluctuations during a period of 100 ps. Nature. 277, 578-578;
11. Peter J. Rossky, Martin Karplus. (1979). Solvation. A molecular dynamics study of a dipeptide in water. J. Am. Chem. Soc.. 101, 1913-1937;
12. D.A. Case, M. Karplus. (1979). Dynamics of ligand binding to heme proteins. Journal of Molecular Biology. 132, 343-368;
13. David L. Bostick, Charles L. Brooks III. (2007). Deprotonation by Dehydration: The Origin of Ammonium Sensing in the AmtB Channel. PLoS Comput Biol. 3, e22;
14. K.Y. Sanbonmatsu, C.-S. Tung. (2007). High performance computing in biology: Multimillion atom simulations of nanoscale systems. Journal of Structural Biology. 157, 470-480;
15. Michael G. (1997). An interview with Bernie Alder. Lawrence Livermore national laboratory;
16. David A. Case. (2000). Perspective on "Dynamics of folded proteins". Theoretical Chemistry Accounts: Theory, Computation, and Modeling (Theoretica Chimica Acta). 103, 332-334.