CC BY
120
□
УДК 577.322 И.Н. Николаев, КВ. Шайтан
МОЛЕКУЛЯРНОЕ МОДЕЛИРОВАНИЕ СЛОЖНЫХ БИОСТРУКТУР, СОДЕРЖАЩИХ МЕМБРАННЫЙ КОМПОНЕНТ
Приводится изучение динамики мембранных структур с помощью метода молекулярной динамики. Проводится классификация динамических свойств аминокислотных остатков в системе вода/вакуум.
Ключевые слова: молекулярное моделирование, молекулярная динамика, мембраны, пептиды, сложные биоструктуры, аминокислотные остатки.
Роль молекулярного моделирования для решения задач в различных областях науки и при совершенствовании технологии производства трудно переоценить. Теоретические методы молекулярной динамики (МД), являющиеся одними из методов молекулярного моделирования, широко применяются при изучении фундаментальных проблем естествознания, а также в прикладных задачах молекулярной биоинженерии, биотехнологии, нанотехнологии, материаловедения и др. Успехи современных молекулярных технологий вызвали значительный интерес к МД биомем-бранных структур. Однако до сих пор остаются проблемы с определением различных макроскопических параметров, позволяющих установить целостную картину поведения пептидов в мембранном окружении. Кроме того, часто при изучении таких структур не уделяется должное внимание методическим аспектам моделирования. Однако проработка протоколов молекулярной динамики и тщательный подход к выбору параметров моделирования, несомненно, необходим.
В данной статье проведено изучение динамики пептид-мембранных структур на основе проработанного и определённого на предыдущих этапах изучения протокола МД, позволяющего получать значения параметров мембран, наиболее близкие к экспериментальным. Накопление траекторий проводилось как с помощью обычного метода МД, так и метода управляемой динамики.
Одной из проблем при изучении столь масштабных объектов, как пептиды и белки в мембранном окружении, ограниченные водной фазой, является ещё недостаточное развитие вычислительных мощностей. Поэтому помимо
НИКОЛАЕВ Иван Никитич - к.ф.-м.н., доцент, проректор по учебной работе ЯГУ имени М.К. Аммосова
ШАЙТАН Константин Вольдемарович - д.ф.-м.н., профессор, зав. кафедрой биоинженерии МГУ имени М.В. Ломоносова
явно заданных мембран в полноатомном окружении, в работе предлагается модель двухфазной системы, состоящей из воды и неполярного растворителя. Модель позволяет получить детальную картину поведения аминокислотных остатков на границе раздела фаз.
Кроме того, в работе продолжено накопление траекторий молекулярной динамики для мембран. Проведена корректировка ряда параметров мембран, таких как коэффициенты латеральной диффузии липидов, коэффициентов самодиффузии молекул воды, флуктуационные характеристики объёма системы, величины латерального давления.
Резуньтаты вычислительных экспериментов позволили оценить полноту и эффективность методики молекулярно-динамических расчётов и моделирования сложных липидных структур при решении задач молекулярного дизайна в сравнении с современным научно-техническим уровнем, а также продемонстрировать возможность создания конкурентоспособной продукции и услуг в области молекулярного дизайна.
1. Модель виртуальной среды
Известно, что структура и функции белка главным образом зависят от того окружения, в которое он помещён. При растворении белка в воде немаловажную роль играет так называемый гидрофобный эффект [1]. Существование трансмембранных белков связано с их специфическими взаимодействиями: как с полярным водным окружением, так и с неполярной средой гидрофобных хвостов билипид-ного слоя мембраны. Способность аминокислотных остатков и различных их последовательностей по-разному взаимодействовать как друг с другом, так и с тем или иным типом окружения, по-видимому, лежит в основе формирования специфической структуры белка и его функциональной активности. Поэтому выявление характера взаи-
модействия аминокислотных остатков с разными типами сред является необходимым для понимания тех факторов, под действием которых существуют белковые молекулы.
Наиболее разумным вариантом является изучение взаимодействия аминокислотных остатков (АКО) с двумя типами окружений с полярной фазой - водой и неполярной фазой. В качестве неполярной фазы иногда используется октанол или н-гексан.
С макроскопической точки зрения поведение веществ по отношению к воде и неполярному окружению можно выразить через свободную энергию переноса из фазы в фазу и свободную энергию адсорбции на границу раздела двух фаз. Экспериментальные работы по определению этих характеристик были неоднократно выполнены, например [2, 3, 4]. В работе [5] предложена классификация веществ одновременно по свободной энергии переноса и адсорбции.
Следует отметить, что таким образом для каждого вещества можно ввести так называемый эмпирический параметр гидрофобности. Он используется, в частности, в фармокинетике, занимающейся предсказанием поведения и распределения лекарственных веществ в организме. Один из таких подходов называется QSAR (Quantitative Structure-Activity Relationship). Гидрофобность является одним из важнейших параметров для реакций, катализируемых ферментами, так как сайты связывания последних очень восприимчивы к гидрофобным участкам молекул субстрата [6].
Для самих аминокислот можно создать шкалу гидрофобности. С помощью шкалы гидрофобности можно построить профили гидрофобности белка, которые используются для предсказания структуры и трансмембранных частей белка [7, 8].
2. Постановка задачи
К сожалению, экспериментальные данные, в основном, позволяют нам лишь судить о макроскопических параметрах взаимодействия аминокислотных остатков с различными средами и границами раздела фаз, в то время как поведение аминокислотных остатков вблизи границы раздела фаз на микроскопическом уровне остаётся недостаточно изученным. Поэтому нами было предложено применить методы МД к исследованию поведения аминокислотных остатков вблизи границы раздела фаз воды и неполярной фазы.
Безусловный интерес представляет рассмотрение динамики АКО на границе билипидной мембраны, однако даже рассмотрение простейшей системы с разделом фаз может дать нам информацию о взаимодействии как молекулы в целом, так и её частей с разными типами фаз.
Следует также отметить, что полноатомное молекулярное моделирование является достаточно трудоёмкой задачей. В настоящее время при существующих вычислительных ресурсах зачастую не удаётся достичь требуемых времён моделирования больших систем. Поэтому одной из
приоритетных задач является построение огрублённой модели двухфазной системы, которая должна хорошо отражать качественные особенности поведения молекул (АКО) в двухфазной среде.
Поставленная в этой работе задача состоит в определении, качественном и количественном описании поведения АКО вблизи границы раздела фаз вода/неполярный растворитель. Она подразделяется на следующие этапы:
1. Создание системы с границей раздела фаз вода/гексан и её моделирование методом МД.
2. МД основных типов аминокислотных остатков вблизи границы раздела фаз вода/гексан.
3. Разработка моделей гидрофобного потенциала. Разработка метода по определению свободной энергии сольватации АКО и его апробация.
4. Систематизация по лученных характеристик для различных групп АКО.
3. Изучение поведения АКО вблизи границы раздела фаз методами равновесной молекулярной динамики
3.1. Модель аминокислотных остатков
Как известно, в водном растворе со стандартным pH аминокислотные остатки, взятые по отдельности, будут представлять из себя ионы с зарядами на обоих концах. N-конец будет существовать в форме NH3+, а C-конец в форме COO-. В такой форме молекулы будут представлять из себя сольватированные диполи и вероятно проявлять сильные гидрофильные свойства. В данной же задаче для нас больший интерес представляло моделирование АКО в тех формах, в которых они встраиваются в белковую последовательность. Поэтому в нашей модели остатки представляли из себя молекулы с формально ненасыщенными связями на N и C концах. Нами использованы параметры силового поля и взаимодействий между молекулами, которые применяются при моделировании пептидных последовательностей целиком.
3.2. Модели гидрофобной мембраны
Необходимым условием моделирования является выбор подходящей модели гидрофобной мембраны. Наиболее приемлемым веществом для моделирования гидрофобной среды явился н-гексан. С одной стороны, в некоторых работах он используется как гидрофобная среда для экспериментального определения свободной энергии переноса, с другой, - молекула гексана обладает небольшим размером, что существенно для молекулярно динамических расчётов [5].
Одновременно с системой/вода гексан нами была предложена система вода/вакуум + столкновительная среда. Несмотря на явное отсутствие моделируемой гидрофобной среды, эта модель дала неплохие результаты по динамике молекул на качественном уровне.
3.3. Модель с виртуальной гидрофобной средой
Данная модель представляет собой следующую систему (рис. 1). Раздел фаз моделировался в виде виртуальной плоскости, которая отделяла водный слой от гидрофобного слоя. Гидрофобный слой в этой модели реальными ча-
стицами не моделировался. Непроницаемая для молекул воды виртуальная плоскость задавалась как отталкивающий потенциал мягких сфер для всех атомов воды. Таким образом, молекулы воды не могут пройти через эту стенку, тогда как исследуемые молекулы АКО могут свободно проходить через эту стенку.
Рис. 1. Система для моделирования вода/вакуум
Практически система была реализована следующим образом. Базовая ячейка размером 80*46*46 А была разделена двумя виртуальными плоскостями на три части: два водных слоя по 20 Ап гидрофобный слой шириной в сере -дине 36 А. Следует учитывать, что моделирование производилось с использованием периодических граничных условий , то есть для каждой молекулы, вылетевшей из ячей -ки, её образ влетал с другой стороны ячейки. Таким обра-зом, система эффективно представляла из себя один водный слой толщиной 40 А, граничащий с гидрофобным слоем. Отметим, что из-за выбранного нами способа моделирования раздела фаз в виде непроницаемых для воды стенок периодические граничные условия для воды выполнялись лишь по двум направлениям, в то время как для исследуемых молекул АКО - по всем направлениям. Исследуемые АКО в большинстве случаев изначально помеща-лись на интерфейс.
Моделирование производилось в NVT-aнcaмблe, то есть количество частиц, объём базовой ячейки и темпера -тура поддерживались постоянными.
Для поддержания постоянной температуры в 300 К использовался столкновительный термостат (метод столкно-вительной динамики).
Суть процедуры термостатирования заключалась в том, что через некоторые промежутки времени случайным образом выбирались атомы и им приписывались новые скорости. Новая скорость вычислялась как скорость, полученная атомом в результате упругого столкновения с вир -туальной частицей термостата, скорость которой выбиралась случайным образом в соответствии с распределением Максвелла по скоростям для заданной температуры. Температура, масса виртуальных частиц, а также частота соударений задавались изначально.
Таким образом, столкновительный термостат частично восполнял отсутствие гидрофобной среды.
3.4. Модель вода/н-гексан
Более реалистичной моделью является прямое моделирование гидрофобной среды с помощью молекул гек-сана. Для этого нами была создана система (рис. 2), состоящая из двух слоёв (ламелей): воды и гексана. Периодические граничные условия наложены по всем осям. По оси, перпендикулярной к интерфейсу, одномерное баростати-рованиев 1 атм. (баростат Берендсена). Размеры системы по XY, фиксированные 40 А, по Z переменный около 41 А.
Слой гексана (~20 А) состоит из 147 молекул, изначально упорядоченных в три слоя по 7*7 молекул в слое, ориентированных перпендикулярно к интерфейсу. Слой воды (~20А) состоит из 1070 молекул. Одна молекула АКО размещается на одном интерфейсе.
і-'* V . г ІД -V
!;■ ■■ t,'-
Л ~ №s ,Vv" f‘
■':A ,
Рис. 2. Система вода/гексан
Существенное отличие этой системы от модели вода/ вакуум + столкновительный термостат состоит в отсутствии отталкивающей стенки, поэтому следовало бы ожидать некоторого перемешивания двух фаз. Однако при построении данной системы мы исходили из того, что для определения структурных особенностей поверхностей жидкостей с начала 1980-х годов используется методика исследования отражённого рентгеновского излучения [9], а в работе [Mitrinovic D.M. et all] приведены данные по толщине границы раздела фаз для системы вода/н-гексан. Профиль плотности одной из фаз при этом обычно аппроксимируется интегралом ошибок по следующей формуле:
р(х) =
р0
"->/2:
п
е 2°2 ds
3.5. Сравнение профилей плотности моделей
Рассмотрим получившиеся профили плотности для двух типов моделей.
Если рассмотреть сначала профиль плотности двухфазной системы вода/гексан (рис. 3), то данный профиль получен путём усреднения статистики за время, равное 2 не. Эти значения плотности в середине фаз с достаточной точностью совпадают с экспериментальными (660 кг/мл3 для гексана и 997 кг/мЛ3 для воды). Из профилей видно, что толщина на границе раздела фаз составляет порядка 5 ангстрем, что согласуется с экспериментальными данными.
Вблизи границы раздела фаз наблюдаются некоторые неоднородности плотности как воды, так и гексана. Это, по-видимому, связано с поверхностными конформацион-ными эффектами упорядочения. Так, измерения нелинейной восприимчивости с использованием генерации второй гармоники допускают, что раздел фаз вода/алкан может быть высоко упорядочен [11]. Для длинных цепей полученные результаты работы исключают возможность перпендикулярного к интерфейсу упорядочения н-алка-нов [10]. Результаты компьютерного моделирования допускают возможность упаковки молекул параллельно интерфейсу [12]. Однако в случае вода/гексан этот вопрос ещё требует дополнительного изучения.
Следует отметить, что данный профиль плотности построен для системы, в которой на правом интерфейсе присутствовал аминокислотный остаток - глицин, фактически не обладающий радикалом и представляющий собой лишь кусок основной цепи полипептида. При этом заметна явная асимметрия профиля плотности воды вблизи двух интерфейсов. Вблизи границы раздела, где АКО отсутствова-ли, более выражен первый максимум, тогда как вблизи другой границы раздела он сглажен. Если мы теперь рассмотрим профиль плотности гексана, то он в большей степени сохраняет свою симметрию. Таким образом, можно
где а - параметр толщины границы раздела. Для системы вода/н-гексан он равен 3,5 А [10].
При этом она соответствует толщине, которую предсказывает капиллярно волновая теория. Вследствие капиллярных волн происходит эффективное уширение границы раздела.
Таким образом, исходя из приведенных экспериментальных данных следует ожидать, что в нашей модели не будет происходить перемешивания фаз и граница раздела будет поддерживаться достаточно чётко.
Рис. 3. Профили усредненной плотности каждой из фаз (воды и гексана) вдоль оси, перпендикулярной границе раздела фаз
2
S
предположить, что присутствие АКО на границе раздела в большей степени разрушает приповерхностную структуру воды, чем гексана. Похожее поведение наблюдалось и для других АКО, адсорбирующихся на интерфейсе.
Рассмотрим профиль плотности нашей модельной системы вода/вакуум (рис. 4). Как видно из этого рисунка, для такой системы профиль плотности по форме и ширине практически совпадает с профилем плотности для воды предыдущей системы. И хотя первый максимум профиля плотности воды вблизи границы раздела фаз в данном случае выражен не так сильно, но вполне возможно, что сходства профилей плотности воды будет достаточно, чтобы уловить основные закономерности поведения АКО и в системе вода/вакуум (рис. 4).
ковыми цепями и с основными или кислотными цепями: ОЬУ, ОЬН ЬУ8, АКО
• Глицин (ОЬУ), который хотя и оставался вблизи границы раздела, но были заметны флуктуации, когда он уходил в водную фазу.
Таким образом, в системе вода/гексан для нас представляло интерес моделирование АКО из первой и третьей групп, так как разумно предположить, что для остатков, десорбирующихся с границы раздела в воду в системе вода/ вакуум, аналогичное поведение было бы свойственно и для системы вода/гексан.
В системе вода/гексан нами была рассмотрена динамика трёх АКО с различными по величине гидрофобными боковыми группами: ОЬУ, АЬА, РНЕ.
го
<
В
(Л
М
4-1
■Н
щ
а
ш
р
І "
7
- / 1 -
- 1 / -
/ /
/ 1 і , . —X , : і
10 :12 14 10 16 20 22 24 28 2S 30
Z-axis,А
Рис. 4. Профиль плотности воды вдоль оси, перпендикулярной границе раздела фаз, для системы вода/ вакуум
Для дальнейшего сравнения систем и рассмотрения поведения АКО вблизи интерфейса рассмотрим его динамику в этих двух системах.
Динамика аминокислотных остатков
Мы исследовали поведение различных типов остатков в системе вода/вакуум. И как следует из сделанных выше замечаний, мы можем использовать результаты данной модели для грубой оценки поведения АКО, помещённых на границе раздела.
По динамическому поведению в системе вода/вакуум исследованные АКО были разделены нами на три группы:
• АКО, которые во временном масштабе компьютерного моделирования остаются на границе раздела. К ним относятся АКО с неполярными, алифатическими, ароматическими боковыми цепями: ALA, VAL, LEU, PHE, MET, TRP.
• АКО, которые десорбируются с границы раздела и уходят в водную фазу - с полярными незаряженными бо-
Динамика аминокислотных остатков с разными радикалами
Рассмотрим полученные результаты по динамике АКО.
Глицин относится к самой маленькой аминокислоте, у которой боковая цепь представляет собой один атом водорода. На профиле плотности вероятности для центра масс АКО (рис. 5) глицин занимает различное положение вдоль оси, перпендикулярной к границе раздела. Линией обозначен центр переходного слоя к границе раздела. За центр нами взята середина профиля плотности воды. Мы видим, что максимум распределения практически точно находится на середине границы раздела.
Рис. 5. Гистограмма плотности вероятности для центра масс молекулы глицина вдоль оси, перпендикулярной к границе раздела
Рассмотрим динамику более детально, для этого на рис. 6 приведены гистограммы отдельно для групп СО и молекулы глицина. Видно, что максимум для СО
группы сдвинут относительно максимума группы примерно на 1 ангстрем в сторону водной фазы.
Рис. 6. Гистограмма плотности вероятности для центров масс групп СО и МИ молекулы глицина вдоль оси, перпендикулярной к границе раздела
системе вода/вакуум иногда отдалялся от границы раздела в водную фазу.
Рассмотрим динамику двух других АКО: аланина и фенилаланина. Для этих остатков также построены гистограммы положения центра масс, а также гистограммы положений центров масс радикала и остова АКО. Остовом АКО мы называем группу -NH-CH-CO-, которая присутствует во всех АКО. Для нас представлял интерес проследить, как будет изменяться динамика АКО при увеличении гидрофобного радикала. В табл. приведены данные о размерах и типах радикалов аланина и фенилаланина.
Как видно из материалов рис. 8, где приведены гистограммы распределения центров масс молекул для всех трёх аминокислотных остатков (GLY, ALA, PHE), увеличение неполярной боковой группы АКО приводит к смещению гистограммы и её максимума в сторону гексана.
Для более детального анализа происходящих изменений необходимо рассмотреть гистограммы для центров
Для более подробного изучения ориентационной динамики молекулы нами была построена гистограмма для вектора, соединяющего центр масс СО группы и центр масс МИ группы (рис. 7). Из этой гистограммы видно, что при динамике молекула располагается преимущественно так, что её остов ориентирован по отношении к поверхности раздела под небольшим углом (порядка 20 градусов). При этом МИ группа несколько больше смещена в сторону гексана, а СО - в сторону воды.
Одним из отличий динамики глицина в системе вода/ гексан явилось то, что он так же, как и другие остатки с гидрофобными радикалами, адсорбировался на интерфейсе и не отходил от него во время динамики. Тогда как в
Таблица
Характеристики некоторых аминокислотных остатков
Название Формула боковой цепи Средняя длина боковой цепи, А Средняя длина всех АКО, А
ALA СНз 1.6 4.4
PHE CH2-C6H5 6 7.2
Рис. 7. Ориентационная гистограмма для вектора, соединяющего центр масс СО группы и центр масс МИ группы молекулы глицина
Рис. 8. Гистограммы плотности вероятности для центра масс молекул АКО вдоль оси, перпендикулярной к границе раздела. Слева направо следуют кривые для PHE, ALA, GLY
масс остова и боковой группы (рис. 7, 8), а также ориентационные гистограммы для вектора, соединяющего центры масс остова и радикала (рис. 9, 10). На основе анализа материалов этих рисунков можно сделать следующие выводы.
Рис. 9. Гистограммы для положения центра масс радикала и остова у молекулы ALA
Рис. 10. Гистограммы для положения центра масс радикала и остова у молекулы PHE
У аланина, имеющего небольшой гидрофобный радикал, гистограмма для центра масс боковой цепи смещена больше в сторону гексана, чем гистограмма для центра масс остова. При этом заметно, что в сравнении с поведением глицина, положение остова у аланина несколько смещается в сторону гексана. Это позволяет утверждать, что гидрофобный радикал стремится находиться в гексане и эффективно несколько вытягивает за собой остов.
Для фенилаланина, имеющего достаточно крупную ароматическую боковую группу, наблюдается аналогичное поведение, однако наиболее вероятное положение остова совпадает со случаем глицина. Вероятно, это можно объяснить стерическими эффектами. Несмотря на то, что радикалу фенилаланина выгодно находиться в гексане, его большой размер не позволяет ему уходить глубоко в гек-сан. Дополнительным аргументом может являться рассмотрение ориентационных гистограмм (рис. 9, 10), где
боковая цепь молекулы фенилаланина направлена иод некоторым углом в сторону фазы гексана. Таким образом, ориентация остатков фенилаланина и аланина не является произвольной, в среднем их боковые группы обращены в сторону гексана иод некоторым углом. При этом ориентационная подвижность молекулы фенилаланина значительно меньше, чем у аланина, что объясняется различными размерами боковых групп (рис. 9, 10).
Более того, если рассмотреть рис. 7 и 8 в тех областях, где молекулы отдаляются от границы раздела в воду, мож-но заметить, что у аланина средние положения центров масс остова и радикала совпадают, тогда как у фенилаланина нет. В отдельных случаях, когда молекула аланина несколько отходит от границы раздела в водную фазу, ориентационная подвижность ее сильно возрастает, что не выражено у молекулы фенилаланина (рис. 11, 12).
90
Рис. 11. Ориентационная гистограмма для PHE
Рис. 12. Ориентационная гистограмма для ALA
Влияние взаимодействия с гексаном
Рассмотрим вопрос о вкладе в динамику АКО и взаимодействие их с гексаном. Для этого рассмотрим гистограммы для положения центра масс остова и радикала аланина в системе вода/вакуум (рис. 13).
Из анализа этого рисунка видно, что появление в моделировании гексана вносит достаточное изменение в положение максимумов гистограмм, хотя и незначительно изменяет их форму. Отсюда можно сделать вывод, что взаи-модействие молекулы АКО с фазой гексана приводит к смещению центра масс молекулы в сторону гексана в среднем на 2 А. При этом расстояние между максимумами гистограмм для центра масс остова и радикала уменьшается. Это можно объяснить тем, что гексан эффективно поджимает молекулу и ориентирует ее в большей степени вдоль интерфейса.
Distance, А
Рис. 13. Гистограммы для положения центра масс радикала и остова у молекулы ALA в системе вода/вакуум
Литература
1. Финкельштейн А.В., Птицын О.Б. Физика белка: Курс лекций. М.:Книжный дом университет, 2002. 455 с.
2. Bull H.B., Breese K. (1974) Arch Biochem Biophys 161:665
3. Wolfenden R., Anderson I., Cullis P.M., Southgate C.C.B. (1981) Biochemistry 20:849.
4. Leodidis E.B., Hutton T.A. (1990) J Phys Chem 94:6411.
5. Ivan M., Okhapkin I.M., Makhaeva E.E., Khokhlov A.R. Colloid Polym Sci (2005).
6. Stryer, Biochemistry L. 4th ed. W.H. Freeman & Company: New York, 1995; Part I, chapter 1.
7. Kyte J. and Doolittle R.F. A Simple Method for Displaying the Hydropathic Character of a Protein . Journal of Molecular Biology 157(6): 105-142, 1982.
8. Hopp T.P., and Woods K.R. Prediction of protein antigenic determinants from amino acid sequences. Proc. Nat. Acad. Sci. USA 78(6): 3824-3828, June 1981.
9. Als-Nielsen, Christensen F and Pershan P.S. Phys. Rev. Lett. 48, 1107 (1982).
10. Mitrinovic D.M., Tikhonov A.M, Li M., Huang Z. and Schlossman M.L. Noncapillary-Wave Structure at the Water-Alkane Interface, Phys. Rev. Lett. 85, 585 (2000).
11. Conboy J.C., Daschbach J.L. and Richmond G.L. Appl. Phys. A 59, 623 (1994).
12. Rivera J.L, McCabe C., Cummings P.T. Molecular simulations of liquid-liquid interfacial properties: Water-n-alkane and water-methanol-n-alkane systems, Phys. Rev. E 67, 011603 (2003).
I.N. Nikolaev, K. V Shaitan
Molecular modeling of complex bio-structures containing a membrane component
The article presents a study of dynamics of membrane structures by method of molecular dynamics. A classification of dynamics, properties of amino acid residues in the system water/vacuum has been suggested.
Key-words: molecular modeling, molecular dynamics, membranes, peptides, complex bio-structures, amino-acid residues.
