Содержание к диссертации
Введение
Глава 1. Структура и динамика биомембран 6
1.1. Биологические мембраны, строение и функции 6
1.2. Исследование структуры биологических мембран 14
1.3. Динамика липидов 19
1.4. Молекулярная динамика 21
Глава 2. Метод молекулярной динамики 30
2.1. Силовое поле в методе молекулярной динамики 33
2.2. Численное интегрирование 42
2.3. Контроль температуры и давления 45
2.4. Ансамбли в молекулярной динамике 49
2.5. Неравновесная молекулярная динамика 51
Глава 3. Сравнительное изучение молекулярной динамики гидратированных мембран различного липидного состава 52
3.1. Постановка задачи 52
3.2. Протокол молекулярной динамики 53
3.3. Релаксация бислоя 56
3.4. Флуктуации параметров бислоя 58
3.5. Распределения электронной плотности 62
3.6. Радиальные функции распределения 63
3.7. Параметры порядка 65
3.8. Латеральная диффузия липидов 66
Глава 4. Исследование проницаемости мембранных структур для ряда низкомолекулярных лигандов 69
4.1. Постановка задачи 69
4.2. Протокол молекулярной динамики 70
4.3. Результаты и обсуждение 73
Выводы 87
Благодарности 88
Литература 89
Приложение 1 106
Приложение 2 108
- Исследование структуры биологических мембран
- Численное интегрирование
- Протокол молекулярной динамики
- Протокол молекулярной динамики
Введение к работе
Биологические мембраны являются важным структурным компонентом живой клетки. Структурные и кинетические свойства биомембран играют важную роль в организации различных клеточных процессов. Однако такие сильно структурированные анизотропные структуры оказываются весьма сложными для изучения. Исследование биомембранных структур методами молекулярной динамики (МД) приобретает в настоящее время все большую популярность в связи с уникальными возможностями визуализации и детализации на атомарном уровне молекулярных процессов.
Выполнение мембраной барьерной и транспортной функций во многом зависит от ее физико-химических и кинетических параметров, определяемых липидным и белковым составом. Это касается и микроскопической картины массопереноса в сложном мембранном окружении, диффузионных процессов на границе водной и мембранной фаз.
В данной работе, метод молекулярной динамики и метод управляемой молекулярной динамики, в полноатомном приближении используются для изучения динамического поведения ряда мембранных структур различного липидного состава, расчета траекторий движения, определения и сравнительного анализа важных с физической точки зрения макроскопических термодинамических и кинетических параметров.
Интересным и до конца не решенным является вопрос трансмембранного транспорта малых молекул. Изучение этого процесса методами МД сталкивается с определенными трудностями, поскольку характерные времена диффузии намного превосходят достижимые в численном эксперименте. Поэтому ниже развит метод управляемой молекулярной динамики (УМД), в рамках которого эволюция системы стимулируется по определенным степеням свободы, что позволяет дать количественную оценку параметрам, характеризующим механизмы
трансмембранного переноса малых молекул за достижимые времена численного эксперимента.
Исследование структуры биологических мембран
МД эксперимент основан на исследовании молекулярной модели и важным моментом является ее калибровка и проверка соответствия между моделью и экспериментальными данными. Примерами могут служить исследования профилей электронной и атомной плотностей, размеров средней площади приходящейся на молекулу липида, параметров порядка для аклильных цепей липидов, объема связанной воды и электростатических потенциалов. Ниже приведено описание некоторых наиболее характерных методик. Дифракционные методы
Эксперименты по нейтронному и рентгеновскому рассеянию являются одними из наиболее мощных методов для изучения структуры биологических систем с атомным разрешением. К сожалению, липидные системы в жидкокристаллической фазе обладают высокой степенью анизотропности и неупорядоченности и применение этой группы методов встречается с определенными сложностями. Тем не менее, определенный прогресс в понимании структурных особенностей поведения липидных бислоев был достигнут с использованием именно этих методов.
Wiener и White [18] с помощью методов по нейтронному и рентгеновскому рассеянию с достаточно хорошей точностью описали структуру жидкокристаллической фазы ДПФХ. В частности, были получены функции распределения для положения атомов различных типов вдоль оси перпендикулярной плоскости бислоя. Следует отметить, что в этих работах степень гидратации изучаемого бислоя была ниже биологически обоснованной. Степень гидратации липидных голов влияет на формирование структуры бислоя. В ранних работах не была получена также информация о структуре бислоя в его латеральной плоскости, в частности, не была определена величина удельной площади, приходящейся на молекулу липида.
Величина средней площади является важным параметром, по которому производится сравнения результатов различных экспериментов. Данный параметр является важным показателем корректности компьютерной модели при МД моделировании и является ключевым при использовании ансамблей при постоянном объеме или постоянной площади на молекулу липида (NVT, NplAT). В отличие от объема, площадь бислоя может претерпевать сильные изменения в зависимости от условий проведения эксперимента. Значительное влияние на нее оказывают такие факторы как липидный состав мембраны, температура среды, степень гидратации бислоя или наличие примесей в растворе [16].
В более поздних работах Nagle и другие авторы [19] определили структуру полностью гидратированного бислоя ДПФХ при помощи комбинации методов рентгеновского рассеивания и теории, учитывающей волновые колебания бислоя. Основные результаты данной работы включают конструктивные особенности бислоя, которые можно получить с помощью Фурье преобразования профиля электронной плотности, а также набор структурных параметров, включающий среднюю площадь (порядка 62.9 ±1.3 А) и объем приходящиеся на молекулу липида, расстояние между пиками для профиля электронной плотности, количество молекул воды, приходящихся на липид в полностью гидратированном состоянии (порядка 29.1). В этих работах был сделан вывод, что наблюдаемые отклонения в профилях электронной плотности, полученных при изменении степени гидратирования не являются результатом масштабных структурных изменений в бислое, а скорее являются артефактом проявляющемся в результате пренебрежения эффектом волновых колебаний бислоя при интерпретации данных эксперимента. Из этого следовало, что при корректном учете эффекта волновых колебаний результаты экспериментов с низкогидратированными бислоями могут являться вполне обоснованными.
Для изучения движения липидов в липидных бислоях их метят флуоресцентным маркером. Исследуя движение меченных молекул можно получать данные о динамике липидов, в частности, коэффициентах объемной и латеральной диффузии [31]. Подобные исследования происходят на сравнительно продолжительных характерных временах от микросекунд до секунд, что значительно больше характерных времен моделирования доступных в МД эксперименте. Во многих случаях остается неясным как соотносится динамика меченых и немеченых липидов, поскольку сама метка типа родамина может влиять на подвижность липида.
Для исследования различных свойств липидных бислоев применяются и другие экспериментальные методы. Примерами могут служить в частности ЭПР - спектроскопия, измерения проницаемости веществ на черных пленках, ИК/Рамановская спектроскопия [32], в том числе, и для изучения параметров порядка и концевых gawc/ie-дефектов, измерения поверхностных потенциалов мембран, коэффициентов распределения вода-мембрана и др.
Здесь угловыми скобками обозначен средний квадрат отклонения центра массы липида в плоскости бислоя за время т. В методе МД траектория, получаемая в эксперименте, содержит координаты всех атомов, что делает вычисление коэффициента латеральной диффузии тривиальной задачей. Экспериментально измерение коэффициента латеральной диффузии сталкивается с определенными трудностями. Наиболее часто коэффициенты латеральной диффузии измеряют, используя флуоресцентные методы [31], квазиупругое рассеяние нейтронов на бислоях [33, 34] и импульсный ЯМР [35, 36].
Численное интегрирование
Данная система уравнений в общем случае не может быть решена аналитически и применяются численные методы. Существуют различные численные методы решения уравнений, различающиеся точностью и вычислительной сложностью. Компромиссом между точностью процедуры и скоростью её реализации, является широко используемый в МД метод Верле [152]. Важно также обеспечить сохранение постоянного числа частиц (/V), объема (У) и энергии (Е) системы.
Согласно (1.12) определяются координаты и скорости каждого атома на основании значения координат в два предшествующих момента времени. При этом скорость не обязательно вычислять для определения траектории. Однако скорость нужна для вычисления кинетической энергии системы, температуры и тензора давления. Схема имеет порядок точности 0(At ) по координатам и 0(At )по скоростям атомов.
Алгоритм имеет тот же порядок точности, что и оригинальный алгоритм Верле, но обладает существенным плюсом при вычислении скоростей, поскольку не производится деления двух близких к нулю значений, как это делается при вычислении скорости в (2.12). Это при определенных условиях позволяет использовать более длинный шаг интегрирования At. Минусом данной схемы является то, что скорости и координаты вычисляются в разные моменты времени, отличающиеся на At/2, что является определенным неудобством. Обычно скорость в момент времени t приближенно вычисляют как среднее арифметическое значений скоростей в моменты времени t+At/2 и t-At/2.
Существуют схемы интегрирования более высокого порядка точности, например, метод Гира (предиктор-корректор [147]). Существенным минусом таких схем является их большая вычислительная сложность. Более того, при использовании больших значений для шага интегрирования такие схемы в целом оказываются даже хуже, чем простые схемы [100]. По этим причинам эти схемы обычно не используются в молекулярной динамике.
Выбор длины шага интегрирования является принципиальным моментом процедуры интегрирования. Длина шага интегрирования должна выбираться таким образом, чтобы позволять надежно отследить самые быстрые осцилляции, возникающие в системе. Поэтому для полно-атомных систем с нефиксированными длинами связей шаг интегрирования выбирают не более 1-2фс (характерное время периода осцилляции валентной -Н связи порядка 10 фс). Если длины связей фиксированы, длина шага интегрирования может достигать 4 фс (характерное время периода осцилляции валентного угла -Н, порядка 20-35 фс). При использовании тяжело-атомных моделей с фиксированными валентными углами и связями величина шага интегрирования может достигать 1 Офс.
Начальные скорости атомов в нулевой момент времени, обычно выбираются случайным образом, в соответствии с распределением Максвелла при заданной температуре. Этих недостатков лишен используемый в работе столкновительный (коллизионный) термостат [155, 156]. Суть его заключается в том, что атомы системы находятся во взаимодействии со средой виртуальных частиц. Атомы виртуальной среды обладают заданной молекулярной массой и способны упруго взаимодействовать с атомами системы с заданной средней частотой.
Здесь Т - длительность промежутка, р(т) — плотность распределения случайной величины. На каждом шаге интегрирования уравнений движения происходит розыгрыш столкновений виртуальных частиц с атомами системы в соответствии с распределением (2.20). В промежутках между столкновениями частица двигается как в традиционной молекулярной динамике в соответствие с уравнениями (2.9). В случае если атом сталкивается с виртуальной частицей, то новая скорость вычисляется как результат упругого столкновения с виртуальной частицей заданной массы т0 и скорости v в соответствии с (2.19).
Управляя значениями параметров частоты Л и массы виртуальных частиц то можно регулировать энергообмен при столкновении с виртуальными частицами и влиять на динамику флуктуации и скорости релаксационных процессов. При расчетах с использованием виртуальной среды следует учитывать ее вязкостные характеристики. Так, при расчетах в вакууме с использованием столкновительной среды с параметрами массы виртуальных частиц 18 а.е.м; и частоты столкновений 55 пс-1, среда по своим вязкостным свойствам близка к воде,при н.у. [157]
При проведении МД экспериментов, часто бывает удобно использовать ансамбли при постоянном давлении, поскольку именно такие условия чаще всего соответствуют состоянию наблюдаемой системы в эксперименте. В частности это касается экспериментальных работ с бислоями липидов [158].
Протокол молекулярной динамики
Для проведения расчетов МД использовался пакет PUMA [155, 156] с графическим интерфейсом PUMA-GUI. Решение классических уравнений производилось с использованием алгоритма Верле в полноатомном силовом поле Amber-99. В системах, содержащих КЛ, отрицательный заряд компенсировался добавлением в воду ионов натрия. Расчетная ячейка в различных системах содержала 20-32 молекулы липида на монослой. Степень сольватации составляла не менее 32 молекул воды на липид, что является достаточным для полного гидратирования липидных головок [169].
Парциальные заряды атомов в липидах оценивались с использованием неограниченного метода Хартри-Фока, базиса 6-31G и метода Малликена. Использовалась модель воды TIP3P [150]. Валентные связи и валентные углы в молекулах воды и липидов не фиксировались. В стартовой конфигурации молекулы воды помещались на расстоянии не менее 2.3 А от крайних атомов мембраны. (3.1)
Расчет проводился в периодических граничных условиях, при постоянных температуре и объеме (NVT - ансамбль), а также при постоянных температуре и давлении (NPT - ансамбль), в зависимости от стадии эксперимента. Температура в зависимости от стадии эксперимента, выбиралась от 300 К до 2000 К. Для поддержания постоянной температуры среды использовалась виртуальная столкновительная среда (столкновительный термостат [155, 156]), с массой виртуальных частиц 1 а.е.м. и средней частотой столкновений с атомами системы 10 пс"1. Для поддержания условия постоянного давления в NPT-ансамбле, использовался баростат Берендсена [170], с частотой от 0.1 до 1 пс"1. Для учета эффектов поверхностного натяжения бислоя величина давления в латеральной плоскости выбиралась отрицательной, подобранной так, чтобы удельная площадь поверхности, приходящаяся на липид (далее, плотность липидов), соответствовала известным экспериментальным значениям [171].
Система изначально собиралась с помощью оригинального молекулярного конструктора так, чтобы плотность лидидов соответствовала экспериментальным данным и составляла в течение расчета для ПОФХ 62-68 А2, для ДПФХ 59-62 А2 [30, 120], для КЛ 100-142 А2 [172, 173]. В случае бислоев из разных липидов удельная площадь вычислялась как среднее арифметическое площадей участвующих липидов, взвешенных согласно их концентрации. Начальная (стартовая) структура бислоев собиралась с использованием специально разработанного программного обеспечения по оригинальной методике, включающей повороты молекул липидов вокруг длинной оси на случайный угол от 0 до 2л (Рис. 13.-6). Это позволяет минимизировать возможную симметрию в начальном положении различных липидов по отношению друг к другу. На Рис. 13. приведена стартовая конфигурация первой модельной системы.
Чтобы уменьшить упорядоченность в области алкильных хвостов, на первом этапе релаксации проводилось моделирование стартовой конфигурации системы в NVT ансамбле при температуре 1000 К, с фиксированными положениями атомов фосфора в головах липидов, в течение не более 200 пс. Затем, полученная система гидратировалась молекулами воды, и проводилась дальнейшая релаксация бислоя в течение 200-500 пс при температуре 300-1000 К в NPT ансамбле. При использовании диэлектрической проницаемости 1 и изотропного баростатирования с одинаковым во всех направлениях компонентами давления 1 атм наблюдалось сокращение удельной площади поверхности липидов до значений, существенно более низких, чем наблюдаемые в эксперименте. Поэтому к системам прикладывалось отрицательное давление в латеральной плоскости, достаточное для удержания удельной площади приходящейся на липид в соответствии с экспериментальными значениями. Затем производился набор рабочего участка МД траектории в течение не менее 2 не. Определение величины необходимого давления производилось на этапе подготовки бислоя. На рабочем участке траектории параметры баростатирования не изменялись.
Необходимость использования анизотропного баростатирования для поддержания необходимой площади расчетной ячейки обусловлена, по-видимому, эффектами поверхностного натяжения и влияния параметров силового поля. Для компенсации эффектов поверхностного натяжения необходимо использовать отрицательное давление в плоскости мембраны, что неоднократно обсуждалось ранее [96, 145, 171, 174]. Исследование зависимости величины отрицательного давления необходимого для поддержания заданной площади, от используемой величины диэлектрической проницаемости указывает на зависимость эффекта от вклада кулоновских сил (Рис.14.).
Таким образом, влияние на наблюдаемый эффект сжатия мембраны в изотропном баростате оказывают погрешности параметризации силового поля, выбора парциальных зарядов и параметров МД протокола. Отметим, что в настоящее время принципиально не может быть решен вопрос о выборе «правильной» параметризации. Само представление о молекулярной системе с заданными сферическими парными потенциалами изначально содержит погрешности, поэтому задачей является разработка таких протоколов и параметров МД, которые позволяют получать модель в согласии с экспериментальными данными.
Плотность вероятности и ее гауссовская аппроксимация для значений объема расчетной ячейки для систем 1-а, И-б, Ш-в. Использован рабочий участок траектории. В соответствии с известными значениями хт Для мембран, которые могут лежать в диапазоне от ІхІО"10 до 6х10"10Па"1 [175], и известными величинами для сжимаемости воды 4,6х10"10 Па"1, можно говорить о согласии вычисленных и экспериментальных значений хт Для исследуемых систем. Таким образом, полученные разумные значения коэффициента сжимаемости Хт, а также гауссовский характер распределения флуктуации объема свидетельствуют о достижении локального равновесия.
Протокол молекулярной динамики
Для проведения МД расчетов использовался пакет PUMA-B [155, 156, 189] с графическим интерфейсом PUMA-GUI [190] с поддержкой управляемых воздействий. Решение классических уравнений производилось с использованием алгоритма Верле в полноатомном силовом поле Amber-99. В качестве исследуемых мембран, использовались изученные выше бислои систем I и III, состоящие из липидов в различных пропорциях: 1-пальмитоил-2-олеоил-5П-глицеро-3-фосфатидилхолин (ПОФХ), 1,2-дипальмитоил-5п-глицеро-3-фосфатидилхолин (ДПФХ), 1,3-( 1 -стеароил-2-пальмитоил-8п-гл и церо-3-фосфатид ил)- глицерин (кардиолипин, КЛ). В системах, содержащих КЛ, отрицательный заряд компенсировался добавлением в воду ионов натрия. Расчетная ячейка в различных системах содержала 20-32 молекулы липида на монослой. Степень сольватации составляла не менее 32 молекул воды на липид, что является достаточным для полного гидратирования липидных головок [169]. Липидный состав исследуемых систем представлен в Табл. 5.
Парциальные заряды атомов низкомолекулярных соединений оценивались с использованием неограниченного метода Хартри-Фока, базиса 6-31G и метода Малликена. Использовалась модель воды TIP3P [150]. Валентные связи и валентные углы в молекулах воды, пенетрантов и липидов не фиксировались.
Расчет проводился в периодических граничных условиях, при постоянных температуре и давлении (NPT - ансамбль). Для поддержания постоянной температуры среды использовалась виртуальная столкновительная среда (столкновительный термостат [155, 156]), с массой виртуальных частиц 1 а.е.м. и средней частотой столкновений с атомами системы 10 пс"1. Температура термостата составляла 300 К. Для поддержания условия постоянного давления, использовался баростат Берендсена [170]. В соответствии с методикой описанной в предыдущей главе, для учета эффектов поверхностного натяжения бислоя и поддержания величины плотности липидов величина давления в латеральной плоскости выбралась отрицательной. "
Кинетика прохождения молекулы аммиака под действием силы 6 ккал/(А-моль) в направлении нормали к мембране через бислой системы I. (а) - Положение молекулы по оси Y (нормали мембраны). Горизонтальными линиями показаны усредненные координаты атомов фосфора фосфолипидов. (б) - Скорость перемещения молекулы.
Под действием внешней молекулы проникают в мембрану. При этом смещение молекулы определяется в основном движением под воздействием внешней силы, а вкладом диффузии практически можно пренебречь.
Заметим, что имеющиеся данные о микровязкости мембран в настоящее время неполные и противоречивые. Экспериментально известна усредненная величина вязкости поверхностного слоя биомембраны, составляющая от 30 до 190 сПз для различных липидных мембран [191-193]. Для ПОФХ известна также экспериментальная оценка усредненной вязкости порядка 18 сПз [101]. Экспериментальные данные по микровязкости в направлении нормали практически отсутствуют. Поэтому прямое сравнение экспериментальных и расчетных величин затруднительно. Следует также отметить, что приближение сплошной среды для молекулярной системы требует определенных пояснений. В случае липидных бислоев значение эффективной микровязкости в направлении нормали к поверхности мембраны следует понимать как некую количественную- характеристику диссипативных свойств бислоя, выраженную в единицах вязкости.
Наблюдаемая общая закономерность такова: чем больше радиус низкомолекулярного лиганда, тем больше эффективная микровязкость мембраны. Это согласуется и с наблюдаемым здесь эффектом увеличения эффективной микровязкости при уменьшении внешней силы (и, следовательно, замедлении скорости проникновения частицы). Отметим, что при наблюдаемой тенденции увеличения величины эффективной микровязкости при увеличении размера молекул иногда наблюдаются некоторая немонотонность, обусловленная, по-видимому, различием в химической природе пробных молекул. В случае глицерола гидрофильная часть составляет порядка 50% от общего размера молекулы. Эта существенная гидрофильная часть служит, по-видимому, своеобразным якорем, удерживающим молекулу в более гидрофильных областях мембраны (Рис. 26.). Сравнительно низкая эффективная вязкость мембраны в случае этанола, по-видимому, объясняется небольшим размером и внутренней подвижностью молекулы, что приводит к меньшей деформации мембраны при проникновении пенетранта. Полученный коэффициент диффузии для кислорода 2.4x10" см /с для смешанной мембраны и 1.7x10" см /с для мембраны ПОФХ близок к коэффициенту диффузии кислорода в воде (2x10" см /с при 293 К и 3.2x10" см/с при 313 К). Это согласуется с некоторыми экспериментальными данными для коэффициента диффузии кислорода в биологических мембранах, который согласно [195-197] лежит в пределах от 1.5x10" см /с до 4.7х10 5 см2/с.
Кинетика наиболее быстрого (1) и наиболее медленного (2) зафиксированного в серии экспериментов прохождения молекул глицерина (a, Ref=) и этанола (б) через бислой системы III под действием силы 6 ккал/(А-моль) в направлении нормали к мембране. Графиком показано положение молекулы по оси Y (нормали мембраны). Горизонтальными линиями показаны усредненные координаты атомов фосфора фосфолипидов.
Динамика повторных прохождений молекулы аммиака через липидную мембрану системы I. Показано положение центра массы остатка, горизонтальными линиями показаны центры масс атомов фосфора соседних бислоев. Среднее время прохода, после 10 повторных испытаний составило 17 пс. Среднее время первого прохода Т по результатам 4 испытаний составило 13 пс, с величиной среднеквадратичного отклонения 3 пс.
Неоднородность микровязкости а разных слоях по отношению к уровню нормали бислоя вычислялась путем измерения средней скорости, с которой молекула проходила каждый участок бислоя в направлении нормали к мембране с заданным значением силы. Используя формулу Стокса-Эйнштейна, вычислялся локальный коэффициент трения, который далее выражался в терминах микровязкости ((4.1), (4.2)).
