Моделирование молекулярной подвижности липидов и проницаемости бислойных мембран
На правах рукописи
Зленко Дмитрий Владимирович МОДЕЛИРОВАНИЕ МОЛЕКУЛЯРНОЙ ПОДВИЖНОСТИ ЛИПИДОВ И ПРОНИЦАЕМОСТИ БИСЛОЙНЫХ МЕМБРАН 03.00.02 — биофизика
Автореферат диссертации на соискание ученой степени кандидата биологических наук
Москва — 2008
Работа выполнена на кафедре биофизики биологического факультета Мос ковского государственного университета имени М.В. Ломоносова.
Научные руководители: кандидат физико-математических наук, доцент Красильников П.М.
доктор биологических наук, профессор, член-корреспондент РАН Рубин А.Б.
Официальные оппоненты: доктор химических наук, профессор, член-корреспондент РАН Чизмаджев Юрий Александрович.
доктор физико-математических наук, профессор Ефремов Роман Гербертович.
Ведущая организация: Кафедра медицинской и биологической фи зики Московской Медицинской Академии им. Сеченова.
Защита состоится 12 марта 2009 года в 14 ч. 00 мин. на заседании Дис сертационного совета Д 501.001.96 при кафедре биофизики биологического факультета Московского государственного университета имени М.В. Ломо носова по адресу: 119991, Москва, Воробъевы горы, МГУ, биологический факультет, кафедра биофизики, аудитория “Новая”.
С диссертацией можно ознакомиться в библиотеке биологического факуль тета Московского государственного университета имени М.В. Ломоносова.
Отзывы на автореферат, заверенные гербовой печатью организации, просьба направлять по указанному адресу в двух экземплярах не позднее, чем за две недели до защиты.
Автореферат разослан 16 января 2009 года.
Ученый секретарь Диссертационного совета доктор биологических наук, профессор Кренделева Т.Е.
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Актуальность темы. Молекулярное моделирование липидных бис лойных мембран представляет значительный интерес в связи с тем, что эти объекты достаточно трудны для экспериментальных исследований. Фунда ментальные закономерности микроскопической динамики молекул липида остаются все еще не вполне ясными, несмотря на прогресс в этой обла сти. В то же время развитие численных методов исследования и совер шенствование вычислительной техники позволяют рассматривать все более и более сложные модели, приближаясь к описанию реально существую щих молекулярных систем. Это дает возможность при помощи численных экспериментов не только получить полезную информацию, не прибегая к трудоемким экспериментам на реальных мембранах, но и провести иссле дования процессов, в принципе недоступных экспериментальным методам.
Это касается и детальной микроскопической картины динамики липидов, и деталей процессов взаимодействия отдельных молекул с окружением, и деталей проникновения молекул сквозь толщу бислоя.
В данной работе метод молекулярной динамики (МД), с использова нием полноатомной модели и силового поля OPLS [1, 2], применяется для уточнения микроскопической картины диффузионных процессов в толще липидного бислоя, а также для моделирования проникновения различных молекул и ионов через мембрану. Использование полноатомного прибли жения делает МД-расчеты больших систем, таких как фрагмент липидной мембраны, окруженный водой, достаточно трудоемкими. Это накладыва ет ограничения на максимальные промежутки времени (длину траекторий), доступные в расчетах МД. Поэтому весьма актуальным является разви тие таких методов обработки данных численных экспериментов, которые позволяли бы, основываясь на расчетах разумной длины, корректно вычис лять искомые величины. Так, например, в большинстве работ, в которых методами молекулярной динамики исследовались параметры подвижности молекул липида, для определения коэффициента диффузии используют под ход, дающий корректные результаты только при длине траекторий в десятки наносекунд [3], а при меньших длинах результат оказываются сильно завы шенным [4, 5, 6] Таким образом, актуальной задачей является разработка подходов, поз воляющих вычислять коэффициенты диффузии и другие параметры моле кулярной подвижности, исходя из траекторий длиной в 1 нс. Разработка дополнительных моделей и методов обработки данных, полученных в МД расчетах, позволило с успехом решить эту задачу и получить из относитель но коротких траекторий корректные оценки как микроскопических парамет ров молекулярной подвижности, так макроскопически измеримых величин.
Кроме того, при помощи предложенных новых методов и модификаций ме тодов, использовавшихся ранее, была измерена проницаемость липидного бислоя для частиц с сильно различающимися свойствами. Показано, что на границах гидрофобной области бислоя располагается кинетический барьер, препятствующий проникновению частиц через мембрану.
Целью работы было исследование диффузионной подвижности мо лекул дистеароилфосфатидилхолина на малых временах методами молеку лярной динамики, а также исследование проницаемости липидных мембран для различных молекул и ионов.
Постановка задачи.
1. Построить и отладить устойчивую полноатомную модель фрагмен та бислойной липидной мембраны, окруженного водой. Параметры модельной системы должны полностью соответствовать параметрам реальных мембран. Модель должна быть устойчивой во времени.
2. Разработать протоколы МД расчетов, наилучшим образом подходящие для работы с предложенной моделью. Определить оптимальный шаг интегрирования, а также параметры силового воздействия на молеку лярную систему для проведения расчетов управляемой МД.
3. Разработать методы оценки параметров подвижности молекул липи да. Методы должны учитывать особенности и ограничения, присущие методу МД и давать корректные оценки макроскопически измеримых параметров.
4. Разработать методы исследования подвижности молекул ДСФХ на временах 10 фс – 10 нс. Рассчитать время оседлой жизни и парамет ры тепловых движений. Описать полученные результаты в терминах кинетической теории жидкостей.
5. Разработать метод расчета параметров проницаемости липидных мем бран для молекул и ионов на основании расчетов управляемой МД.
6. Исследовать проникновение молекул воды и кислорода, а также иона натрия через ДСФХ бислой. Выявить характеристики этого процесса, не зависящие от свойств частицы и рассчитать коэффициенты прони цаемости бислоя для этих частиц.
Научная новизна.
Раработаны новые методы оценки коэффициентов диффузии в плоско сти бислоя и времени вращательной корреляции для молекул липида, позволяющие оценивать эти параметры на основании МД расчетов, длиной в несколько наносекунд. Полученные результаты хорошо со гласуются с известными из эксперимента величинами.
Предложен метод оценки времени оседлой жизни, основанный на из мерении средней скорости движения молекулы липида. С помощью этого метода удалось измерить время оседлой жизни молекул ДСФХ и обнаружить коллективные тепловые движения молекул липида.
Для описания движения молекул липида предложена модель неизо тропных двумерных случайных блужданий, что позволило продемон стрировать диффузионную природу тепловых колебаний и оценить такие их характеристики, как амплитуду и частоту. Полученные ре зультаты существенно дополняют физическую картину строения жид костей. В частности показано, что параметры молекулярных движений на временах порядка 1 пс определяют макроскопически наблюдаемый коэффициент диффузии.
Предложен вариант метода оценки внутримембранных коэффициен тов диффузии, основанный на измерении подвижности частицы в тол ще бислоя. Этот метод имеет ряд преимуществ перед предложенным ранее вариантом [6] и позволяет получать согласующиеся с экспери ментальными данными оценки коэффициентов проницаемости.
Получены данные, указывающие на положение внутри бислоя кине тического диффузионного барьера, препятствующего проникновению через бислой небольших частиц. Барьер располагается на границе гид рофобной области мембраны, и его положение не зависит от природы проникающей частицы.
Практическая значимость работы заключается в разработке методов вычисления микроскопических молекулярных параметров, таких как коэф фициенты диффузии, времена оседлой жизни и вращательной корреляции, а также параметры теплового движения молекул. Предложенные в работе подходы могут быть использованы для решения многих биофизических и молекулярно-биологических проблем с использованием молекулярной ди намики и расширяют возможности этого метода.
Апробация работы. Результаты работы были представлены автором на XIV международной конференции студентов, аспирантов и молодых уче ных “Ломоносов” (Москва, 2007), на научно-практической конференции в рамках международной научно-образовательной школы-конференции по биоинженерии и приложениям “Перспективы развития инноваций в биоло гии” (Москва, 2007), на XV научной конференции “Математика компьютер образование” (Дубна, 2008), на XV международной конференции студентов, аспирантов и молодых ученых “Ломоносов” (Москва, 2008). По материалам работы было сделано также несколько докладов на лабораторных семинарах и семинарах кафедры биофизики биологического факультета МГУ.
Публикации. По теме диссертационной работы опубликовано 6 работ, в том числе 1 статья в журнале, рекомендованном ВАК.
Структура и объем диссертации. Работа изложена на 134 страницах, проиллюстрирована 50 рисунками и содержит 5 таблиц. Диссертация состо ит из введения, 4 глав, заключения, выводов и списка литературы, содержа щего 131 источник.
ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ
В первой главе представлен обзор литературы, первая часть которо го посвящена строению липидных мембран. Освещены такие вопросы как разнообразие мембранных липидов и их функции, а также строение и топо логия водно липидных смесей. Во второй части первой главы дан обзор ре зультатов, полученных численными методами для чисто липидных систем.
Рассматриваются различные молекулярные модели липидных систем, ис следования посвященные изучению взаимодействия молекул липидов друг с другом и с водой, изучению подвижности липидов, а также проницаемо сти и слиянию мембран. Дана критическая оценка этих работ и сравнение с результатами, полученными в данной работе.
Во второй главе представлен обзор метода молекулярной динамики.
Дана общая характеристика, рассмотрены особенности и ограничения ме тода. Описано устройство и приемы построения силовых полей, а также файлов молекулярной топологии. Показано, каким образом в модель мо лекулы липида, построенную на основании силового поля OPLS [1, 2], были добавлены недостающие в силовом поле параметры. Для расчета этих параметров использовался пакет квантово-химических программ PC GAMESS [7]. Расчеты велись для малых молекул, в соответствующих кон формациях, представляющих фрагменты молекулы липида, с использова нием DFT и функционала B3LYP5. Подробно обсуждаются возможности корректировки и настройки модели средствами GROMACS [8], а также оп тимизация констант жесткости в случаях сильной связи между внутренними координатами молекулярной системы. В частности, оптимизация констант жесткости, основанная на сравнении спектров собственных частот и Гесси анов молекулярно-механической и квантово-химической систем.
Далее рассмотрена тонкая настройка модели и параметры расчетов. В частности, выбор способа численного решения уравнений движения и под бор оптимального шага интегрирования. Для решения уравнений движения был выбран “стохастически-динамический” (sd) способ, так как этот ал горитм включает в себя алгоритм стохастического термостатирования, что делает его более устойчивым к случайным молекулярным флуктуациям и в значительной степени нивелирует ошибки, вносимые в расчеты термостатом Берендсена [9, 10, 11].
Оптимальный шаг интегрирования был определен, как наибольший из возможных, для которых реальная температура молекулярной системы мало отличается от заданной и не происходит “вымораживания” колеба тельных степеней свободы [10]. Вторым критерием выбора служил неболь шой разброс в координатах атомов в конце 10 пс траекторий, при одинако вых стартовых координатах и различных стартовых скоростях, что отражает устойчивость алгоритмов решения уравнений движения. Чем устойчивее ал горитм, тем медленнее накапливаются ошибки численного интегрирования.
Этим критериям удовлетворяет шаг = 0.25 фс, при значении константы термостатирования 0.5 пс.
Контроль термодинамической достоверности модели производился на основании следующих критериев: распределение энергии системы по сте пеням свободы должно быть нормальным;
средняя энергия, приходящаяся на одну степень свободы, должна составлять 2 kT, что и было достигну то. Расчеты проводились при температуре 300К. Нормальное давление в системе поддерживалось равным 1 атм., тангенциальное – 300 атм/нм, в со ответствии с поверхностным натяжением, характерным для липидных мем бран [12].
Далее описано построение модели фрагмента липидной мембраны.
48 модельных молекул ДСФХ были уложены в бислойную гексагональную структуру (рис. 1). Затем в систему было добавлено 1759 молекул воды. Для воды была использована модель tip4p [13]. После этого был проведен 10 нс релаксационный МД расчет, в ходе которого геометрические и динамиче ские параметры модельной мембраны достигли устойчивых значений, хоро шо совпадающих с таковыми для реальных бислойных мембран, состоящих из ДСФХ и находящихся в жидкокристаллическом (ЖК) состоянии (табл. 1).
Более того, над поверхностью мембраны образовался слой малоподвижной воды, толщиной порядка нескольких, и характеризующийся коэффициен A том самодиффузии 107 см2/с, что на два порядка ниже, чем в объемной фазе (2.1 · 105 см2/с), что соответствует экспериментальным данным [14].
Для построения модели липидной мембраны в гелевом состоянии ис ходная модельная система была подвергнута действию тангенциального дав ления 105 атм. После релаксации молекулярной системы при нормальном давлении (в течении 1 нс), геометрические параметры полученной липид ной мембраны оказались близки к таковым для ДСФХ в гелевом состоя нии (табл. 1). Однако такие параметры оказались неустойчивы и медленно релаксировали к параметрам, характерным для ЖК состояния. Обсуждаются Рис. 1. Общий вид молекулярной модели участка липидной мембраны в жидкокристаллическом состоянии, окруженной водой. Модельная система состоит из 48 молекул дистеароилфосфатидилхолина и 1759 молекул воды.
На рисунке отмечены: a – вода ( 20 ), b – липидная мембрана ( 65 ), A A c – гидрофобная область бислоя ( 35 ) и d – полярные головки ( 15 ).
A A возможные причины неустойчивости модели, в частности ее малые разме ры и возможное некорректное поведение процедуры численного решения уравнений движения.
В конце первой главы рассмотрены методы оценки поступательной и вращательной подвижности молекул липида, предложенные автором дис сертационной работы. Для измерения коэффициента диффузии на основа нии данных численного эксперимента предлагается использовать метод, ос нованный на измерении кинетики роста дисперсии расстояний от центров масс (ЦМ) некоторой группы молекул до их общего ЦМ. Показано, что та кой метод оценки коэффициентов диффузии обладает рядом преимуществ перед широко используемым способом непосредственного его вычисления, как отношения квадрата смещения к времени за которое это смещение про изошло. Преимущества предлагаемого подхода связаны с взаимным вычи танием случайных вкладов, даваемых в оценку коэффициента диффузии Таблица 1. Геометрические и динамические характеристики ДСФХ мем браны в гелевом и жидкокристаллическом состоянии. Представлены вели чины рассчитанные для модели и экспериментальные данные. S – средняя площадь, приходящаяся на липидную молекулу;
L – толщина гидрофобной области бислоя;
– угол наклона ацильных цепей в гелевом состоянии;
D коэффициент латеральной диффузии;
– время вращательной корреляции.
состояние жидкокристаллическое гелевое модель эксперимент модель эксперимент S ( 2 ) A 48.7 52 34.1 L ( ) A 35.3 32 53 (град) – – 10 D (см /сек) 1.46 108 109 (сек) 1010 1010 109 каждой молекулой. Как следствие, для измерения коэффициентов диффузии предлагаемым методом требуются более короткие расчеты, а полученные результаты существенно ближе к экспериментальным (рис. 2, A).
Для определения интенсивности вращательных движений предложен метод оценки времени вращательной корреляции, основанный на анализе эволюции распределения молекул липида по углу поворота относительно начального положения. В начальный момент времени это распределение яв ляется -функцией, с течением времени молекулы поворачиваются, что при водит к переходу распределения в нормальное. На достаточно больших вре менах, когда молекулы полностью “забудут” о своем начальном положении, распределение станет равномерным. В качестве оценки времени вращатель ной корреляции было выбрано время, за которое доверительный интервал ( = 50%) для математического ожидания этого распределения достигает [, ]. То есть это время за которое в среднем половина молекул повернется на угол больший ±, соответственно молекулы в среднем повернуться друг относительно друга на угол (рис. 2, B). Поскольку в основе эксперимен тального измерения времени вращательной корреляции лежат аналогичные процессы дезорганизации, предложенный метод дает оценки хорошо согла Рис. 2. А. Зависимость дисперсии распределения расстояний от ЦМ группы молекул липида, до ЦМ группы. Синим показаны результаты численного эксперимента, красным - линия линейной регрессии, тангенс угла накло на которой соответствует коэффициенту диффузии. В. Положение границ симметричного доверительного интервала (50%) для мат. ожидания рас пределения группы липидов по углу поворота относительно начального положения. Рубеж /4 соответствует времени вращательной корреляции.
Красным показаны данные, описывающие динамику изменения положения границ доверительного интервала для ацильных хвостов, синим - полярных головок.
сующиеся с экспериментальными величинами (табл. 1). Отметим, что время вращательной корреляции, рассчитанное для гидрофобной части молекулы, оказалось большие, чем для полярных головок, что означает, что враща тельная подвижность полярных головок ацетилхолина выше, чем ацильных хвостов. Аналогичный результат был получен ранее [15].
В третьей главе представлено описание исследования подвижности молекул липида, на основании двух новых подходов, позволяющих обраба тывать результаты молекулярно-динамических расчетов. Предложенные ме тодики позволяют описать движение молекул и дать количественные оцен ки параметров этого движения в диапазоне времен от 100 фс до 10 нс.
Вычисленные параметры подвижности молекул хорошо согласуются с экс периментальными данными (в тех случаях, когда эти параметры могут быть измерены экспериментально).
В первой части третьей главы предложен метод измерения времени оседлой жизни. Метод основан на анализе динамики изменения средней скорости движения ЦМ молекулы липида. Допустим молекула находится в состоянии “оседлости”, совершает периодические колебания в полости, об разованной окружением [16]. Тогда средняя скорость ее движения на этом отрезке времени будет невелика. Если же в течении промежутка времени, на котором была измерена средняя скорость, молекула совершит переход меж ду двумя соседними полостями, то измеренная средняя скорость возрастет.
Такие колебания средней скорости движения ЦМ молекулы липида были зафиксированы при анализе МД траекторий (рис. 3, A).
Предложено два альтернативных подхода к анализу зависимости сред ней скорости движения от времени. Первый из них основан на ана лизе автокорреляционных функций динамики изменения средней скоро сти (рис. 3, B).. Очевидно, первый минимум такой функции соответству ет времени оседлой жизни, так как оно соответствует примерно половине полного периода наблюдаемых колебаний. Об этом свидетельствует гар монический вид зависимости средней скорости от времени (рис. 3, A). В среднем для всех липидов, в модели ЖК состояния для времени оседлой жизни получено значение 550 пс, при 300 К.
Второй метод оценки времени оседлой жизни основан на анализе Фу рье спектров динамики изменения средней скорости. Во всех полученных спектрах присутствовали два хорошо выраженных пика - аппаратный, свя занный с конечной длиной траектории, и пик 950 ГГц, соответствующий полному циклу колебаний средней скорости (рис. 4, A). Соответственно, анализ Фурье спектров дает для времени оседлой жизни 530 пс. При анали зе спектров было обнаружено, что их высокочастотная часть имеет быстро затухающую колебательную структуру (рис. 4, B). Это говорит о том, что в молекулярной системе ЦМ молекул липида организованно смещаются на определенные расстояния. Такие коллективные движения были отождеств лены с инфразвуком, стоячими звуковыми волнами, которые генерируются и поддерживаются тепловыми колебаниями молекул во всяком ограниченном куске любого материала.
Во второй части третьей главы обсуждается тепловая подвижность мо лекул липида. В ходе анализа МД траекторий было показано, что отношение квадрата смещения ЦМ молекулы липида к времени за которое это смеще Рис. 3. A. Динамика изменения средней скорости ЦМ липида, вычислен ная с характерным временем 250 пс. Красным показаны истинные значения средней скорости, синим – скользящее среднее (5%). B. Усредненная по всем липидам автокорреляционная функция динамики изменения средней скорости движения центра масс молекул липида. Первый минимум соот ветствует времени оседлой жизни, максимум полному периоду колебаний средней скорости.
ние произошло (Def f ) зависит от этого времени (рис. 5, A). Параметр Def f был назван эффективным коэффициентом диффузии, а временной проме жуток, использованный для его вычисления - временем измерения (t). При вычислении Def f на малых (t 100 пс) временах, случайные тепловые дви жения не успевают усредняться и появляется зависимость Def f от времени измерения (что было также показано [15]), имеющая вид несимметрично го колокола (рис. 5, A). На больших ( 1 нс) временах Def f, очевидно, переходит в коэффициент диффузии и от времени не зависит.
Для описания наблюдаемой зависимости была предложена модель, де тально описывающая движение ЦМ молекул липида. При ее построении весь временной диапазон был разбит на три области. Первая - малые време на, существенно меньшие периода тепловых колебаний. При рассмотрении движения молекулы в таком масштабе времени ЦМ молекулы равномер но движется от одной стенки полости, образованной окружением, к дру гой. Следовательно, вклад такого движения в суммарный Def f имеет вид:
(1) Def f = (Vwarm t) et/t1, где Vwarm – тепловая скорость движения молекулы липида, t0 – характерное время тепловых колебаний, t – время измерения Def f.
При увеличении масштаба измерения этот режим должен смениться другим, связанным с периодическими тепловыми колебаниями в полости, образованной окружением. Для этого в модель введена экспонента, описы вающая затухание первого режима и проявление второго. Соответственно, (2) R2 (1et/t2 ), где R – ампли для второго режима можно записать: Def f = 2t туда тепловых колебаний. При дальнейшем увеличении времени измерения колебательный режим не затухает, а продолжает вносить вклад в суммарное смещение ЦМ молекулы.
Рис. 4. A. Низкочастотная область усредненного по всем липидам Фурье спектра динамики изменения средней скорости движения ЦМ молекул ли пида. Хорошо виден пик, соответствующий полному периоду колебаний средней скорости (второй) и аппаратный пик, соответствующий длине тра ектории (первый). B. Усредненный по всем липидам Фурье спектр динамики изменения средней скорости движений ЦМ молекулы липида. Хорошо вид на периодическая структура высокочастотной области спектра, свидетель ствующая о существовании в модельной системе коллективных движений молекул липида.
При дальнейшем увеличении времени измерения проявляется третий режим, диффузионный, который заключается в периодических случайных перемещениях между потенциальными ямами – диффузионных шагах. Для описания этого режима была использована модель двумерных случайных блужданий. Предполагалось, что все случайные диффузионные шаги име ют одинаковую протяженность и происходят за одинаковое время. Если предположить, что для каждого случайного шага все направления эквива лентны, то для вклада диффузионного режима с суммарный Def f получим:
(3) r2 (1 et/t2 ), где r – характерная длина случайного шага, t0 – харак Def f = t терное время диффузионного шага и время, начиная с которого этот режим становиться значимым. Однако, использование такой модели привело к рас хождению между модельной кривой и данными численного эксперимента в области 1-10 пс (рис. 5, A).
Рис. 5. Зависимость отношения среднего квадрата смещения к времени за которое оно произошло (Def f ) от этого времени. Кружками показаны дан ные, полученные из МД траектории. Красная кривая – результат фитиро вания данных численного эксперимента при помощи изотропной модели движения молекул липида – A и неизотропной модели – B. Хорошо видно, что введение в модель предположение о неравновероятности направлений очередного случайного шага существенно улучшило совпадение между дан ными численного эксперимента и модельной кривой.
Для уменьшения этого расхождения в модель было введено предпо ложение о зависимости направления очередного случайного шага от на правления предыдущего. Введение этого предположения существенно улуч шило совпадение модельной кривой с данными численного эксперимен та (рис. 5, B). В этом случае третий, диффузионный член приобретает вид:
t 1 t 1 2k 3 4k 1 ( 2k ) 2 (3) Def f = r2 0 1 et/t + t (2k 1) t2 2k где k – константа, описывающая неравновероятность направлений очередно го случайного шага, для которой получено значение k = 1.7, что означает большую вероятность возвратных движений. После того, как молекула пе реместиться в соседнюю полость, образованную окружением, вероятность возврата в три раза выше вероятности диффузионного движения в том же направлении (рис. 6).
Рис. 6. Распределение плотности вероятности для азимутального угла очередного случайного диффузионного шага, относительно направления предыдущего шага, построенная в полярных координатах. Плотность по строена в предположении, что предыдущий случайный шаг был совершен в направлении 0o.
Немаловажным результатом является то, что коэффициент диффузии, рассчитанный как предел Def f :
r2 2k D = lim Def f = 0· t2 2k t с полученными в ходе фитирования данных численного эксперимента па раметрами, имеет значение 2.9 · 108см2/сек, что очень хорошо согласуется с экспериментальными данными [12]. Это означает, кроме того, что пред ложенная модель правильно описывает тепловые движения молекул, что макроскопический коэффициент диффузии определяется подвижностью мо лекул на временах 0.1 – 10 пс, то есть на временах тепловых колебаний. А значит величина коэффициента диффузии не связанна с временем оседлой жизни, которое, как показано в первых разделах третьей главы, имеет суще ственно большее значение.
Таблица 2. Оценка параметров молекулярной подвижности молекул ДСФХ в плоскости липидного бислоя. Представлены результаты использования изотропной и неизотропной модели. Температура 310 К.
Параметр Использованная модель изотропная неизотропная Характерное время смены первого и второго режимов (t0, пс) 0.19 ± 5.1 0.16 ± 0. Характерное время диффузионного шага (t0, пс) 1.48 ± 1.22 1.18 ± 0. Характерная длина случайного шага (r, пм) 4.1 ± 2.7 4.5 ± 0. Средняя амплитуда тепловых колебаний (R, пм) 9.6 ± 1.1 8.4 ± 0. Константа анизотропности k – 1.63 ± 1. Коэффициент диффузии (D, см2 /сек) 3.09 · 107 2.86 · В четвертой главе описано исследование проницаемости ДСФХ бис лоя для иона натрия и молекул воды и кислорода. Такой выбор пенетрантов объясняется принципиальными различиями в их строении, что позволит выявить особенности процесса проникновения, не зависящие от природы проникающей частицы.
В первой части четвертой главы рассмотрен метод неравновесной мо лекулярной динамики и его применение для исследования проницаемости липидных мембран. Суть метода заключается в использовании внешних сил, приложенных к различным частям модельной системы. В настоящей работе внешняя сила была приложена к проникающей через мембрану частице, что существенно ускоряло процесс прохождения бислоя. Величина приложенно го силового поля варьировалась в широких пределах. Это было необходимо для того, чтобы убедиться что измеренный коэффициент диффузии не за висит от величины приложенной силы. В среднем время проникновения исследуемой частицы через бислой составляло 500 – 1000 пс.
Далее описан процесс точного вычисления приложенной к проникаю щей частице силы. Под действием слишком большой силы частица разгонит ся, разогреет систему, плотность снизится, следовательно снизится коэффи циент трения, тормозящий частицу. Это приведет к еще большему разгону и разогреву и, в конце концов, к испарению модельной системы. Для того, что бы этого не происходило необходимо использовать небольшие силы. Но в гетерогенных системах, в частности в системе бислой-вода, невозможно выбрать такую силу, которая подходила бы для обеих частей системы.
Рис. 7. A. Нормальная к поверхности мембраны компонента ускорения иона натрия, соответствующая тормозному воздействию термостата Берендсена.
Кружками показаны реальные значения ускорения, красным - сглаженная кривая. Хорошо видно, что на границах мембраны тормозное воздействие ниже, следовательно приложенная сила больше. Серым показан профиль суммарной плотности системы, силовое поле составляло 128 нм/пс2, темпе ратура 300К.
Для решения этой проблемы необходимо, очевидно, использовать си ловое поле, величина которого зависит от силы трения, действующей на частицу. Однако, в явном виде задать такую зависимость нельзя, так как именно тормозное воздействие среды является предметом исследования.
Поэтому для масштабирования силы был использован термостат Берендсе на с маленьким значением (0.001 пс) константы термостатирования. Термо стат, настроенный таким образом, вносит очень существенные возмущения в уравнения движения, активно притормаживая проникающую частицу в тех областях, где коэффициент жидкого трения мал, а скорость движения велика. Соответственно, как только пенетрант достигнет области с высоким коэффициентом трения, а скорость его движения уменьшиться, тормозное воздействие термостата ослабеет, а сила, соответственно, возрастет (рис. 7).
Необходимо отметить, что при использовании такого подхода к масшта бированию внешней силы, исследователь не вносит в модель никаких до полнительных предположений относительно свойств среды, что было бы внесением в модель априорных предположений о итоговом результате ис следования.
Во второй части четвертой главы рассмотрен метод определения внут римембранного коэффициента диффузии и полученные с его помощью ре зультаты. Проницаемость липидных мембран характеризуют коэффициен том проницаемости, связывающим поток вещества (J) с разностью его кон центраций (C) по разные стороны мембраны: J = L · C. Эксперимен тально было установлено, что L связан с внутримембранным коэффициен том диффузии проникающей частицы (D), ее коэффициентом распределе D ния () и толщиной мембраны (l): L = [17]. В наших расчетах толщина l мембраны была принята равной 10 нм. Точное значение коэффициента рас пределения взято из литературы. Поэтому, для определения коэффициента проницаемости L ключевым являлось определение внутримембранного ко эффициента диффузии.
Ранее был предложен метод вычисления среднего внутримембранно го коэффициента диффузии на основании МД расчетов [6]. В основе его лежит соотношение Эйнштейна: D = kT B, k – константа Больцмана, T – температура, а B – подвижность частицы. Подвижность частицы связа на с дрейфовой скоростью v и силой Fout: v = BFout. В вязкой среде на частицу действует сила вязкого трения: Ff ric = v, где - коэффициент вяз кого трения. Если внешняя сила уравновешивается силой трения, то можно Рис. 8. A. Красным показана динамика изменения средней скорости движе ния иона натрия при его движении сквозь бислой под действием внешней силы. Хорошо видно, что на границах гидрофобной области бислоя ско рость движения иона сильно падает. B. Красным показан трансмембранный профиль коэффициента диффузии для иона натрия. Черным на графиках показан профиль суммарной плотности системы. Приложенная сила соот ветствовала ускорению 64 нм/пс2, температура 300К.
приравнять Fout и Ff ric, тогда получим: B =. Для сферических частиц коэффициент вязкого трения можно выразить через вязкость среды – и радиус частицы – r: = 6r. Таким образом, для коэффициента диффузии kT kT получаем: D = [6].
= 6r Коэффициент диффузии может быть вычислен также на основании измерения подвижности частицы. В этом случае нет необходимости при равнивать внешнюю действующую силу и силу вязкого трения, а значит отпадает предположение о равномерности движения. Отметим, что это тре бование не выполняется при движении частицы в такой сильно неоднород ной среде, как липидный бислой (рис. 8). Для проверки работоспособности этого метода были проведены МД расчеты, в которых частицы двигались в воде под действием внешней силы. Полученные коэффициенты диффузии мало зависели от величины приложенной силы, а их значения представлены в табл. 3, откуда видно, что полученные оценки коэффициентов диффузии хорошо согласуются с экспериментальными значениями.
Для всех исследованных частиц процесс проникновения через бислой имел ряд общих черт. Исследуемая частица под действием внешней силы Таблица 3. Коэффициенты диффузии (см2/сек) молекул воды и кислорода, и иона натрия в воде и липидном бислое. Коэффициенты диффузии рас считаны как D = kT B, где B - измеренная в МД подвижность. Истинные коэффициенты диффузии даны по [18, 19, 20].
Коэффициенты Молекула Молекула Ион диффузии кислорода воды натрия в воде расчетный 3.38 ± 1.02 · 105 2.10 ± 0.71 · 105 6.33 ± 2.37 · в воде истинный 1.8 2.6 · 105 2.27 · 105 1 2 · трансмембранный 3.09 ± 1.75 · 106 3.79 ± 3.72 · 107 2.61 ± 4.00 · в водной фазе движется равномерно. Затем сильно замедляется на границе мембраны, соответственно снижается и коэффициент диффузии (рис. 8). В срединной области бислоя коэффициент диффузии и дрейфовая скорость движения возрастают. А на границе с водной фазой, на выходе из мембра ны, вновь уменьшаются (рис. 8). Замедление пенетранта на границах гид рофобной области бислоя наблюдалось во всех проведенных расчетах, что свидетельствует о существовании кинетического барьера, препятствующего проникновению различных агентов, именно в этих областях мембраны, на уровне плотно упакованных начальных участков остатков жирных кислот.
Ранее в МД исследованиях были получены схожие данные [6].
Экспериментально, на реальных мембранах, невозможно измерить профили, аналогичные представленным на рис. 8. Поэтому, необходимо, очевидно, усреднить полученные трансмембранные профили коэффициента диффузии внутри мембраны, чтобы получить средний внутримембранный коэффициент диффузии и рассчитать на его основе коэффициент проницае мости. Подчеркнем, что такое усреднение необходимо только для получения оценок экспериментально измеренных величин, необходимых для верифи кации результатов.
Далее обсуждаются полученные таким образом оценки коэффициен тов проницаемости для исследованных частиц. В частности, коэффици ент проницаемости крупных лецитиновых везикул для натрия Llarge 1010 см/сек [21], в то время как для мелких липосом этот коэффициент существенно ниже Lsmall 1014 см/сек [22, 21]. Это свидетельствует о влиянии кривизны липидных мембран на их проницаемость для натрия.
Следовательно, коэффициент проницаемости полученный в наших расче тах должен быть по порядку величины равен коэффициенту проницае мости крупных везикул или превышать его. Коэффициент проницаемости D = 1.8 · 109 см/сек, рассчитанный исходя из наших расчетов на Lteor = l порядок превышает коэффициент проницаемости крупных везикул Llarge.
Для воды было получено значение коэффициента проницаемости бислоя Lteor = 1.6 · 104 см/сек. Это значение близко к коэффициенту проница емости лецитиновых везикул Llipo = 5 · 104 см/сек [23]. Для кислорода было получено значение коэффициента проницаемости Lteor = 12 см/сек, что несколько меньше проницаемости мембраны клеток млекопитающего Lmembr = 42 см/сек [24], что связано по всей видимости с существованием систем транспорта кислорода через мембраны живых клеток.
Выводы:
1. Построена устойчивая полноатомная модель фрагмента липидной мембраны, окруженного водой. Модель состоит из 48 молекул ДСФХ и 1759 молекул воды. Геометрические и динамические параметры мо дельной системы соответствуют параметрам реальных ДСФХ мем бран.
2. Разработаны и отлажены протоколы проведения МД расчетов, направ ленных на изучение подвижности молекул липида и проницаемости липидных мембран. Определена оптимальная величина шага интегри рования, равная 0.25 фс.
3. Для расчета коэффициента диффузии предложен метод оценки диспер сии центров масс, позволяющий вычислять коэффициенты диффузии более точно, по сравнению с существовавшими ранее методами, и исходя из более коротких траекторий.
4. Для расчета времени вращательной корреляции предложен метод ос нованный на анализе эволюции распределения углов поворота молекул липида относительно начального положения. Предложенный метод да ет результаты, хорошо согласующиеся с экспериментальными данны ми. Подтверждены данные о том, что полярные головки ацетилхолина имеют более высокую вращательную подвижность, чем гидрофобная часть молекулы.
5. Для описания динамики липидов на временах 0.1 – 10 нс предложен метод, основанный на анализе динамики изменения средней скорости движения ЦМ молекул липида. Показано, что периоды относитель ной оседлости молекулы липида сменяются периодами перемещений.
Рассчитано время оседлой жизни ДСФХ.
6. Для описания динамики липидов на временах 0.01 – 10 пс была пред ложена модель, детально описывающая подвижность молекул липида как функцию времени измерения. Рассчитаны параметры тепловых движений ДСФХ.
7. Показано, что характеристики тепловых движений на временах 0.1 – 10 пс полностью определяют величину макроскопического коэффици ента диффузии. Рассчитаны амплитуды и частоты чисто колебатель ных и чисто диффузионных тепловых движений.
8. Предложена модификация метода расчета коэффициента проницаемо сти бислоя, основанного на измерении подвижности частицы в мем бране. Метод дополнен масштабированием внешней силы при помощи высокоэффективного термостатирования проникающей частицы.
9. Показано, что на границе гидрофобной области бислоя существует ки нетический барьер, определяющий диффузионное сопротивление ли пидных мембран. Рассчитаны коэффициенты проницаемости ДСФХ бислоя для иона натрия, а также молекул воды и кислорода.
Список публикаций по теме диссертационной работы:
1. Зленко Д.В., Красильников П.М. “Построение молекулярно динамических моделей молекулярных систем на основе квантово химических расчетов.” Тезисы на международной конференции сту дентов, аспирантов и молодых ученых “Ломоносов”, 2007. Секция “Биоинженерия и биоинформатика”, С. 62-63.
2. Зленко Д.В., Красильников П.М. “Молекулярная модель фрагмента бислойной липидной мембраны.” Тезисы на международной конфе ренции студентов, аспирантов и молодых ученых “Ломоносов”, 2008.
Секция "Биоинженерия и биоинформатика подсекция “Молекулярное моделирование”, С. 6-7.
3. Зленко Д.В. “Молекулярное моделирование липидных систем”. Те зисы на научно-практической конференции в рамках международ ной научно-образовательной школы-конференции по биоинженерии и приложениям “Перспективы развития инноваций в биологии”. 2007.
Москва, МГУ имени М.В. Ломоносова, Биологический факультет.
2007. С. 33-34.
4. Зленко Д.В., Красильников П.М. “Молекулярно моделирование бис лойных липидных мембран”. Тезисы на XV научной конференции “Математика компьютер образование”, 2008. Дубна. Вып. 15, С. 462.
5. Зленко Д.В., Красильников П.М. “Молекулярное моделирование ли пидных мембран.” Сборник научных трудов конференции “Математи ка компьютер образование”. 2008. Вып. 15, Т. 3, С. 61-70.
6. Zlenko D.V., Krasil’nikov P.M. “Construction of a Molecular Model of a Fragment of a Lipid Membrane in the Gel and Liquid Crystalline States”.
Moscow University Physics Bulletin. 2008. V. 63. No. 6. P. 405–409.
Список цитируемой литературы:
1. Jorgensen W.L., Tirado-Rives J. // J. Am. Chem. Soc. 1988. V. 110. P. 1657– 1666.
2. Jorgensen W.L., Maxwell D.S., Tirado-Rives J. // J. Am. Chem. Soc. 1996.
V. 118. P. 11225–11236.
3. Pitman M.C., Suits F., Gawrisch K., Feller S.E. // J. Chem. Phys. 2005. V.
122. P. 244715.
4. Essmann U., Berkowitz M.L. // Biophys. J. 1999. V. 76. P. 2081–2098.
5. Jin B., Hopnger A.J. // Pharm. Res. 1996. V. 13. P. 1786–1794.
6. Шайтан К.В., Антонов М.Ю., Турлей Е.В., Левцова О.В., Терёшкина К.Б., Николаев И.Н., Кирпичников П.М. // Биологические мембраны.
2008. Т. 25(1). С. 66–75.
7. Nemukhin A.V., Grigorenko B.L., Granovsky A.A. // Moscow University Chemistry Bulletin. 2004. V. 45(2). P. 75.
8. Lindahl E., Hess B., van der Spoel D. // J. Mol. Mod. 2001. V. 7. P. 306–317.
9. Berendsen H.J.C., Postma J.P.M., DiNola A., Haak J.R. // Phys. Rev. 1985.
V. 31. P. 1695–1697.
10. Голо В.Л., Шайтан К.В. // Биофизика. 2002. Т. 47(4). С. 611–617.
11. van der Spoel D., Lindahl E., Hess B., el. all. GROMACS user manual. The Netherlands, AG Gronongen, 2000.
12. Р.Б. Генис. Биомембраны. Москва, изд-во Мир, 1997.
13. Jorgensen,W.L., Chandrasekhar, J., Madura, J.D., Impey, R.W., Klein, M.L.
// J. Chem. Phys. 1983. V. 79. P. 926–935.
14. Stern O. // Z. Elektrochem. 1924. V. 30. P. 508–516.
15. Moore P.B., Lopez C.F., Klein M.L. // Biophys. J. 2001. V. 81. P. 2484–2494.
16. Я.И. Френкель. Кинетическая теория жидкостей. Москва, Изд-во Ака демии Наук СССР, 1945.
17. Overton E. // Vierteljahrsschr. Naturforsch. Ges. Zurich. 1899. V. 44. P. 88– 135.
18. Шапошник В.А. // Вестник ВГУ. Серия: Химия. Биология. Фармация.
2003. V. 2. P. 81–85.
19. Subczynski W.K., Hyde J.S. // Biophys. J. 1984. V. 45. P. 743–748.
20. Влаев Л.Т., Гениева С.Д. // Журнал Структурной Химии., 2004. V. 45(5).
– P. 870–876.
21. Л.Б. Марголис, Л.Д. Бергельсон. Липосомы и их взаимодействие с клет ками. Москва, Наука, 1986.
22. Hauser H., Oldani D., Phillips M.C. // Biochemistry. 1973. V. 12. P. 4507– 4517.
23. Ye R., Verkman A.S. // Biochemistry. 1989. V. 28. P. 824–829.
24. Subczynski W.K., Hopwood L.E., Hyde J.S. // J. Gen. Physiol. 1992. V. 100.
P. 69–87.