Содержание к диссертации
Введение
1 Липидные мембраны, литературный обзор 8
1.1 Строение липидных структур 8
1.1.1 Мембранные липиды 8
1.1.2 Водно-липидные смеси 14
1.2 Молекулярное моделирование липидных мембран 20
1.2.1 Параметры молекулярных моделей 20
1.2.2 Модели липидных структур 23
1.2.3 Исследование взаимодействия липидов 29
1.2.4 Граница бислой/вода 31
1.2.5 Подвижность липидов 32
1.2.6 Проницаемость липидных мембран 34
1.2.7 Слияние мембран 36
2 Молекулярное моделирование 39
2.1 Силовое поле и построение модели молекулы липида 39
2.2 Метод молекулярной динамики 47
2.3 Методы оценки латеральной подвижности 54
2.3.1 Прямой метод оценки 55
2.3.2 Оценка дисперсии центров масс 57
2.4 Оценка времен вращательной корреляции 59
2.5 Построение модели липидной мембраны 63
2.5.1 Модель жидкокристаллического состояния 65
2.5.2 Модель гелевого состояния 68
3 Диффузионная подвижность молекул ДСФХ 71
3.1 Наносекундный временной диапазон 71
3.1.1 Автокорреляционный анализ 73
3.1.2 Анализ Фурье образов 75
3.2 Пикосекундный временной диапазон 79
3.2.1 Изотропная модель случайных блужданий 80
3.2.2 Неизотропная модель случайных блужданий 85
4 Проницаемость липидных мембран 98
4.1 Управляемая молекулярная динамика 100
4.1.1 Внешняя сила переменной величины 101
4.1.2 Расчет коэффициентов проницаемости 104
4.2 Диффузионное сопротивление липидных мембран 108
4.2.1 Общие закономерности 109
4.2.2 Значения коэффициентов проницаемости 111
Заключение 117
Выводы 119
Литература 121
- Водно-липидные смеси
- Метод молекулярной динамики
- Автокорреляционный анализ
- Внешняя сила переменной величины
Введение к работе
Актуальность темы. Молекулярное моделирование липидных бис-лойных мембран представляет значительный интерес в связи с тем, что эти объекты достаточно трудны для экспериментальных исследований. Фундаментальные закономерности микроскопической динамики молекул липида остаются все еще не вполне ясными, несмотря на прогресс в этой области. В то же время развитие численных методов исследования и совершенствование вычислительной техники позволяют рассматривать все более и более сложные модели, приближаясь к описанию реально существующих молекулярных систем. Это дает возможность при помощи численных экспериментов не только получить полезную информацию, не прибегая к трудоемким экспериментам на реальных мембранах, но и провести исследования процессов, в принципе недоступных экспериментальным методам. Это касается и детальной микроскопической картины динамики липидов, и деталей процессов взаимодействия отдельных молекул с окружением, и деталей проникновения молекул сквозь толщу бислоя.
В данной работе метод молекулярной динамики (МД), с использованием полноатомной модели и силового поля OPLS [1, 2], применяется для уточнения микроскопической картины диффузионных процессов в толще липидного бислоя, а также для моделирования проникновения различных молекул и ионов через мембрану. Использование полноатомного приближения делает МД-расчеты больших систем, таких как фрагмент липидной мембраны, окруженный водой, достаточно трудоемкими. Это накладывает ограничения на максимальные промежутки времени (длину траекторий), доступные в расчетах МД. Поэтому весьма актуальным является развитие таких методов обработки данных численных экспериментов, которые позволяли бы, основываясь на расчетах разумной длины, корректно вычислять искомые величины. Так, например, в большинстве работ, в которых методами молекулярной динамики исследовались параметры подвижности
молекул липида, для определения коэффициента диффузии используют подход, дающий корректные результаты только при длине траекторий в десятки наносекунд [3], а при меньших длинах результат оказываются сильно завышенным [4-6]
Таким образом, актуальной задачей является разработка подходов, позволяющих вычислять коэффициенты диффузии и другие параметры молекулярной подвижности, исходя из траекторий длиной в 1 не. Разработка дополнительных моделей и методов обработки данных, полученных в МД расчетах, позволило с успехом решить эту задачу и получить из относительно коротких траекторий корректные оценки как микроскопических параметров молекулярной подвижности, так и слагающихся на их основе на больших временах макроскопически измеримых величин. Кроме того, при помощи предложенных новых и вариантов ранее использовавшихся методик была измерена проницаемость липидного бислоя для частиц с сильно различающимися свойствами. Показано, что на границах гидрофобной области бислоя располагается кинетический барьер, препятствующий проникнове-, нию частиц через мембрану.
Целью работы было исследование диффузионной подвижности мо1 лекул дистеароилфосфатидилхолина на малых временах методами молекулярной динамики, а также исследование проницаемости липидных мембран для различных молекул и ионов.
Постановка задачи.
Построить и отладить устойчивую полноатомную модель фрагмента бислойной липидной мембраны, окруженного водой. Параметры модельной системы должны полностью соответствовать параметрам реальных мембран. Модель должна быть устойчивой во времени.
Разработать протоколы МД расчетов, наилучшим образом подходящие для работы с предложенной моделью. Определить оптимальный шаг интегрирования, а также параметры силового воздействия на молекулярную систему для проведения расчетов управляемой МД.
Разработать методы оценки параметров подвижности молекул липида. Методы должны учитывать особенности и ограничения, присущие
методу МД и давать корректные оценки макроскопически измеримых параметров.
Разработать методы исследования подвижности молекул ДСФХ на временах 10 фс - 10 не. Рассчитать время оседлой жизни и параметры тепловых движений. Описать полученные результаты в терминах кинетической теории жидкостей.
Разработать метод расчета параметров проницаемости липидных мембран для молекул и ионов на основании расчетов управляемой МД.
Исследовать проникновение молекул воды и кислорода, а также иона натрия через ДСФХ бислой. Выявить характеристики этого процесса, не зависящие от свойств частицы и рассчитать коэффициенты проницаемости бислоя для этих частиц.
Научная новизна.
Раработаны новые методы оценки коэффициентов диффузии в плоскости бислоя, времени оседлой жизни и вращательной корреляции для молекул липида, позволяющие оценивать эти параметры на основании ДМ расчетов, длиной в несколько наносекунд. Полученные результаты хорошо согласуются с известными из эксперимента величинами и имеют ряд преимуществ перед предложенными ранее аналогами.
Предложен метод оценки времени оседлой жизни, основанный на измерении средней скорости движения молекулы липида. С помощью этого метода удалось измерить время оседлой жизни молекул ДСФХ и обнаружить коллективные тепловые движения молекул липида.
к для описания движения молекул липида впервые предложена модель неизотропных двумерных случайных блужданий, что позволило продемонстрировать диффузионную природу тепловых колебаний и оценить такие их характеристики, как амплитуду и частоту. Полученные результаты существенно дополняют физическую картину строения жидкостей. В частности показано, что параметры молекулярных движений на временах порядка 1 пс определяют макроскопически наблюдаемый коэффициент диффузии.
* предложен вариант метода оценки внутримембранных коэффициентов диффузии, основанный на измерении подвижности частицы в толще бислоя. Предложенный метод имеет ряд преимуществ перед предложенным ранее вариантом [6] и позволяет получать согласующиеся с экспериментальными данными оценки коэффициентов проницаемости.
-к получены данные, указывающие на положение кинетического диффузионного барьера, препятствующего проникновению через бислой небольших частиц, внутри бислоя. Барьер располагается на границе гидрофобной области мембраны, и его положение не зависит от природы проникающей частицы.
Практическая значимость работы заключается в разработке методов вычисления микроскопических молекулярных параметров, таких как коэффициенты диффузии, времена оседлой жизни и вращательной корреляции, а также параметры теплового движения молекул. Предложенные в работе подходы могут быть использованы для решения многих биофизических и молекулярно-биологических проблем с использованием молекулярной динамики и расширяют возможности этого метода.
Водно-липидные смеси
Как было сказано выше, молекулы липида амфифильны, то есть каждая молекула имеет две части: полярную, обычно называемую головкой, и неполярную - представляющую собой, как правило, два углеводородных "хвоста". За счет высокой гидрофобности углеводородов в водно-липидных смесях молекулы липида стремятся расположиться таким образом, чтобы полярные головки были экспонированы в воду, а гидрофобные участки молекул - наоборот, с водой не контактировали. Этим условиям удовлетворяют несколько принципиально различающихся способов взаимного расположения липидов друг относительно друга. Соответственно, в водно-липидных смесях возможно существование нескольких топологически различных состояний.
В смеси с водой липиды могут образовывать коллоидный раствор, а могут образовывать отдельную фазу, с устойчивым для данных условий составом. Самым простым вариантом сосуществования липидов с водой является коллоидный раствор липида в воде, в котором липиды собраны в мицеллы - сферические частицы, поверхность которых покрыта полярными головками, а середина занята гидрофобным ядром, образованным неполярными участками молекул. Такая форма существования водно-липидной смеси возможна, очевидно, только в том случае, если площадь полярной части молекулы больше площади гидрофобной части. Близким к такому коллоидному раствору по составу состоянием являются неупорядоченные микроэмульсии, отличающиеся неправильной и часто существенно неизо-диаметрической формой скоплений молекул липида, которые могут образо- вывать относительно протяженные структуры. В остальном микроэмульсии похожи на мицеллярный коллоидный раствор.
В том случае, если площадь полярной головки примерно равна площади гидрофобной части молекулы, возможно образование липидных структур с нулевой кривизной.- В таких структурах липиды будут располагаться параллельно друг другу. Очевидно, возможно образование монослойных пленок на границах раздела, например, между водой и воздухом, или водой и неполярной жидкостью. В этом случае полярные головки будут экспонированы в водную фазу, а неполярные ацильные хвосты в воздушную или неполярную. Второй вариант организации не искривленной фазы - это бис-лойные структуры, в которых протяженная неполярная фаза заменена вторым монослоем липидов. Такое строение имеют, как было отмечено выше, мембраны живых клеток.
Липидные бислои могут быть по разному организованы в пространстве друг относительно друга. Расстояние между ними может существенно варьировать, в зависимости от липидного состава мембран и ионной силы растворителя. В одних случаях, как правило между несильно заряженными мембранами, возможно слипание бислоев друг с другом за счет перекрывания диффузных частей двойных электрических слоев. При этом между мембранами остается лишь небольшая (1.5-2 нм) прослойка воды. В англоязычной литературе для обозначения такого состояния используют слово stack (стек). Однако, в случае высокого содержания заряженных липидов или низкой ионной силы расстояние между липидными бислоями возрастает и стек не образуется.
Следует также добавить, что как правило липидные бислои замыкаются сами на себя с образованием "гигантских мультиламмелярных липосом" - замкнутых пузырей, стенки которых образованны бислойной мембраной. Причем в боле крупных пузырях содержаться более мелкие, и уровней вложения может быть произвольное количество. Стенки гигантских мультиламмелярных липосом могут быть собраны в стеки, а могут быть и свободны. Естественно, что эти различия очень существенно влияют на физические свойства водно-липидной смеси. В принципе, мембранные пузыри могут быть и однослойными, тогда говорят об однослойных липосомах.
Сильно варьирует размер липосом. Действием на смесь фосфолипи-дов различного состава с водой ультразвука высокой интенсивности можно получить однослойные сферические липосомы диаметром 40-50 нм, имеющие предельную для липидного бислоя кривизну. Такие везикулы взвешены в водной фазе и, с некоторой натяжкой, могут рассматриваться как вариант коллоидного раствора липида в воде. Такие липосомы из-за небольшого разброса в размерах и способности образовывать коллоидные растворы часто используют на практике, например, в качестве модели для исследования проницаемости мембран и функционирования мембранных белков, а также в фармацевтике и медицине, где липосомы такого строения пытаются использовать для направленной доставки веществ.
Для водно-липидных смесей существуют еще два неламеллярных способа пространственной организации липидов. Первый из них - это гексагональные фазы. Прямая гексагональная фаза Hi представляет собой плотно уложенные в гексагональном порядке тяжи, образованные молекулами липида. Каждый тяж (круглый в поперечном сечении) снаружи покрыт полярными головками, которые экспонированы в непрерывную водную фа зу. В середине каждого тяжа располагаются неполярные участки молекул липида. Обратная гексагональная фаза Нц имеет, соответственно, инвертированное строение (рис. 1.5 - А). Теперь гидрофобные участки молекул липида образуют сплошную фазу, в которой имеются узкие каналы, выстланные полярными головками и заполненные полярным растворителем (водой). Очевидно, на то, какое строение будет иметь гексагональная фаза, прямое или обратное, прежде всего оказывает влияние соотношение площадей полярной головки и гидрофобной части молекулы липида. Для липидов с объемными головами будет характерна прямая гексагональная фаза, а для липидов объемными неполярными участками молекулы - обратная.
Неламеллярные варианты организации водно-липидных смесей. А. Принципиальная схема строения ячейки инвертированной гексагональной фазы Нц; В. Схема, отображающая топологию кубической фазы, имеющей октаэдрическую симметрию; С. Губчатая структура, по сути нерегулярный вариант кубической фазы. А по [8], В и С по [9].
Второй вариант неламеллярной организации видно-липидной смеси это кубические фазы (рис. 1.5 - В). В первом приближении, прямая кубическая фаза Qi представляет собой плотно слипшиеся мицеллы, между которыми произошло слияние в местах контакта. В соответствии с пространственным расположением мицелл друг относительно друга итоговая кубическая фаза может иметь тетраэдрическую, октаэдрическую или иную симметрию. Отметим, что никаких упоминаний о существовании додекаэд-рической или икосаэдрической фаз в литературе обнаружено не было, но и отрицания их существования тоже. Обратная кубическая фаза Qn отличает ся от прямой взаимным расположением водной и углеводородной фаз, как и в случае с гексагональными фазами. Соответственно, в инвертированной кубической фазе углеводородные части молекул липида образуют сплошную фазу, в которой имеются более или менее изодиаметрические полости соединяющиеся друг с другом и заполненные водой.
В соответствии с тем фактом, что геометрия каждого узла кубической фазы может, вообще говоря, отличаться от геометрии соседних узлов, возможно существование так называемых "губчатых структур" - по сути своей нерегулярного аналога кубической фазы (рис. 1.5 - С). Это явление аналогично существованию микроэмульсий, которые можно рассматривать и как нерегулярный мицеллярный раствор, так и как нерегулярный вариант прямой гексагональной фазы, с короткими тяжами и сильно нарушенной гексагональной структурой. Отметим, что для губчатых структур, также как и для кубических фаз характерна обратная, по сравнению с гексагональными фазами кривизна. Иными словами, прямая кубическая фаза образуется из липидов с объемными гидрофобными частями молекул, а обратная - из ли-пидов с крупными полярными головками. В то же время для гексагональных фаз все строго наоборот.
Метод молекулярной динамики
Расчеты МД проводились с использованием пакета GROMACS [85, 86]. Для корректности расчетов МД важным является выбор способа численного решения уравнений движения и подбор оптимального шага интегрирования. Для решения уравнений движения был выбран "стохастически-динамический" (sd) способ, так как этот алгоритм включает в себя алгоритм стохастического термостатирования, что делает его более устойчивым и нивелирует ошибки, связанные с работой термостата Берендсена [87-89].
В ходе численного решения уравнений движения, в связи с ненулевым значением шага интегрирования неизбежно накапливаются ошибки. Отдельные колеблющиеся частицы в конце шага интегрирования оказываются, например, не точно в точке поворота, а чуть дальше. Это означает, что кинетическая энергия соответствующего модельного колебания возросла на некоторую величину. Накапливаясь с течением времени такие ошибки приводят к существенному увеличению температуры всей системы. Надо отметить, что такой самопроизвольный процесс неизбежен при использовании численных методов решения временных уравнений движения. Необходимо также отметить, что уменьшение шага интегрирования влечет за собой ослабление этого негативного эффекта, так как с уменьшением шага интегрирования существенно уменьшаются модули отдельных ошибок.
Для нивелирования описанного выше негативного эффекта используют алгоритмы термостатирования - численные процедуры, отбирающие у системы излишки кинетической энергии. Самым распространенным алгоритмом термостатирования является "термостат Берендсена", в котором в уравнения движения введена специальная функция, описывающая знакопеременное трение. Профиль потенциальной энергии С-С-О-С двугранного угла метил-ацетата. По оси абсцисс отложен угол в градусах, по оси ординат - энергия в кДж/моль. Кружками показан истинный профиль потенциальной энергии, рассчитанный методами квантовой химии. Зеленая кривая - суммарный профиль Кулоновскои и Ван-дер-Ваальсовой потенциальной энергии, полученный методом молекулярной динамики. Квадратами показана разность между истинным профилем потенциальной энергии и результатами молекулярного моделирования. Синим показан результат фитирования разности истинного и молекулярно-динамического профилей потенциальной энергии функцией Рикарта-Беллермана. Красная кривая - сумма профиля молекулярно-динамической Кулоновскои и Ван-дер-Ваальсовой потенциальной энергии и потенциала Рикарта-Беллермана. - = 7 -, где Т - истинная температура, То - референсная температура, г - временная константа, характеризующая скорость убывания отклонений температуры. Технически, в пакете молекулярно-динамических программ GROMACS термостат Берендсена реализован посредством масштабирования скоростей отдельных частиц на каждом шаге интегрирования. При этом коэффициент масштабирования вычисляется как где T{t — At) - отклонение истинной температуры на предыдущем шаге интегрирования от заданной То, At - шаг интегрирования, тт - временная константа, связанная cr = 2CVTT/Ndjk, где Cv - изохорная теплоемкость системы, Ndf - количество степеней свободы системы, к - константа Больц-мана.
Термостат Берендсена обладает рядом преимуществ перед другими схемами термостатирования. В частности этот алгоритм требует относительно небольших вычислительных ресурсов, по сравнению с другими схемами. Есть у этого термостата и существенные недостатки, основным из которых является некорректное перераспределение энергии между степенями свободы молекулярной системы [88]. В частности, использование термостата Берендсена приводит к частичному вымораживанию угловых колебательных степеней свободы. Этого недостатка лишен альтернативный алгоритм термостатирования, известный в рускоязычной литературе как "столкнови-тельный" [90, 91]. Суть процедуры сводится к тому, что в систему вводят газ, состоящий из виртуальных частиц, обладающих Максвелловским распределением по скоростям. При этом в расчетах траектории каждой отдельной частицы не рассматриваются, а учитывается лишь передаваемый ими молекулярной системе импульс. Соответственно, в зависимости от кинетической энергии модельной частицы такая процедура способна как охладить, так и разогреть ее. Основным минусом такого подхода является его чрезмерная ресурсоемкость. Зависимость реальной температуры (К) молекулярной системы от величины шага интегрирования в фс. Синим показаны данные для эффективной температуры системы валентных С-С связей ацильных хвостов молекулы липида. Красным - для валентных С-С-С углов ацильных хвостов молекул липидов. Кружки - экспериментальные данные, сплошные кривые - фитирующие экспоненциальные функции. Референсная температура 310К, г = 0.01 ps.
Для проверки работоспособности комбинированного термостат нами была проведена серия тестовых молекулярно-динамических расчетов. Для этого одиночная модельная молекула липида была помещена в ограниченный объем, после чего проведены расчеты молекулярной динамики с различным шагом интегрирования (рис. 2.5). Из приведенных графиков видно, что при шаге интегрирования в пределах 0.1 - 0.5 фс, реальная температура не сильно отличалась от заданной, что свидетельствует о корректной работе процедуры термостатирования. Однако при увеличении шага интегрирования до 0.5 фс или более, реальная температура угловых степеней свободы существенно снижалась, а температура валентных связей возрастала. Перегрев колебательных степеней свободы валентных связей при больших шагах интегрирования связан с большими по модулю ошибками, возникающими в процессе численного интегрирования уравнений движения. При больших шагах интегрирования велика вероятность того, что реальная температура окажется в несколько раз больше заданной. В этом случае даже такой мощный термостат, как термостат Берендсена оказывается неспособным быстро (за несколько шагов интегрирования) нивелировать перегрев, что неминуемо отразится на средней энергии связей. Это явление хорошо иллюстрирует вид распределения молекулярной системы по длинам С-С связей (рис. 2.6, А и В). Такой же механизм лежит в основе частичного перегрева угловых степеней свободы при больших шагах интегрирования. Таким образом, использование в МД расчетах шагов интегрирования более \ фс некорректно. Полученный результат хорошо согласуется с выводами сделанными другими авторами [53].
Частичная потеря энергии угловыми степенями свободы связанна с некорректной работой термостата Берендсена, которая лучше проявляется при экстремально низких значениях величины шага интегрирования (0.01-0.08 фс, рис. 2.5). С увеличением шага интегрирования за счет возникающих ошибок эффективная температура угловых степеней свободы возрастает и при шаге 0.5 фс угловые степени свободы перегреваются. Вымораживание угловых колебательных степеней свободы хорошо иллюстрируется распределением системы по значениям углов (рис. 2.6, С и D). На приведенных графиках хорошо видно насколько сильно заужается распределение молекулярной системы по угловым колебательным степеням свободы, при использовании малых шагов интегрирования.
Автокорреляционный анализ
Более полное представление о точной величине времени оседлой жизни дают автокорреляционные функции, построенные для динамики изменения средней скорости (рис. 3.3). На графиках представлены автокорреляционные функции, рассчитанные для одиночных молекул липида. На графике (рис. 3.3, В) можно выделить два широких пика, соответствующих временам 1 - 1.1 не и 2 - 2.2 не. Второй пик является обертоном, соответствующим удвоенному времени первого пика. Его появление связанно с самим принципом построения автокорреляционных функций. Следовательно физическим смыслом обладает только первый пик и соответствующее ему время, которое и можно связать с временем оседлой жизни. Стоит лишь оговориться, что время оседлой жизни вдвое меньше, так как периоду оседлости соответствует половина полного периода колебаний средней скорости, представленных на рис. 3.2. Таким образом время оседлой жизни для молекул ДСФХ в бислое в жидкокристаллическом состоянии можно отождествить с первым минимумом на автокорреляционной функции (рис. 3.3, В) и можно оценить как 550 пс. 020 500 1000 Ї500 2000 2500 3000 " 50 500 1000 1500 2000 2500 3000
Автокорреляционные функции, рассчитанные для динамики изменения средней скорости движения центра масс одной молекулы липида. Характерное время вычисления средней скорости 10 пс - А и 250 пс - В, температура 310 К. Первый минимум на рис. В соответствует времени оседлой жизни, которое составило 550 пс.
Конкретный вид экспериментально полученных кривых, описывающих динамику изменения средней скорости движения ЦМ молекул, аналогичных представленным на рис. 3.2, существенно варьирует от молекулы к молекуле. Диффузионные процессы случайны, поэтому зависимости средней скорости от времени для разных молекул дистеароилфосфатидилхолина, как минимум, сдвинуты по фазе друг относительно друга, а на самом деле еще и обладают своими индивидуальными особенностями, связанными с локальными молекулярными флуктуациями. Соответственно, усреднение и иная совместная обработка исходных кривых невозможна. Для получения оценки времени оседлой жизни целесообразно подвергнуть сглаженные зависимости средней скорости от времени Фурье-анализу. Полученные спектры могут быть усреднены по всем однотипным молекулам, что даст более репрезентативный результат. Периоды, соответствующие полученным в ходе Фурье-преобразования частотам, будут соответствовать различным параметрам молекулярной подвижности.
Усредненные по всем однотипным молекулам, нормированные Фурье спектры (ГГц) динамики изменения средней скорости движения центра масс молекул дистеароилфосфатидилхолина. Сначала динамика изменения средней скорости была сглажена методом скользящего среднего с 5% рамкой, затем для каждой сглаженной кривой вычислен Фурье-образ. Полученные спектры были усреднены. Хорошо видны два пика в низкочастотной области. Первый пик - аппаратный, связан с длиной МД траектории, второй пик соответствует периоду колебаний средней скорости движения ЦМ молекулы ДСФХ. Характерное время вычисления средней скорости 10 пс - А и 250 пс - В, температура 310 К. случаях, когда характерное время измерении мало и при визуальной оценке кривая зависимости средней скорости от времени не обнаруживает никакой периодичности (рис. 3.2, А), в усредненных Фурье спектрах можно наблюдать хорошо выраженные пики. Более того, мы наблюдали появление этих пиков в Фурье спектрах полученных в восьми независимых 6 и 10 пико-секундных расчетах, что свидетельствует о неслучайности их появления. Стоит отметить, что при использовании малых ( 10 пс) времен измерения, амплитуда пиков существенно снижается.
В Фурье-спектрах, представленных на рис. 3.4, не показаны первые члены разложения, соответствующие аппаратным характеристикам, таким как нормировочная константа и пик, соответствующий времени 6 не, то есть полной длине траектории. В низкочастотной части спектра хорошо различимы два пика, соответствующие частотам 0.3 - 0.5 и 0.9 - 1.0 ГГц. Низкочастотный пик и связанное с ним время, так же как и в случае с авто корреляционными функциями является обертоном, соответствующим удвоенной частоте основного колебания, наблюдаемого в динамике изменения средней скорости (рис. 3.2).
Положение пика, который мы отождествили с временем оседлой жизни колеблется в зависимости от температуры и конкретного расчета от 0.9 до 1.1 ГГц, что соответствует периоду 0.9 - 1.1 не. Соответственно, время оседлой жизни можно оценить как 0.45 - 0.55 не, так как периоду оседлости соответствует половина полного периода, наблюдаемого в Фурье спектре пика.
Появление дополнительного пика, по всей видимости, связанно с существенным отличием реальной модельной системы от идеализированной жидкости [94]. Таким образом, описание динамики движений молекул модельной системы не может быть сведено к простому чередованию периодов оседлости и периодов трансляционного переноса. Существование дополнительного пика свидетельствует о наличии нескольких различных режимов движения, каждый из которых принципиально соответствует режиму, описанному для идеальной жидкости, но различающихся своими параметрами, в частности характерными временами.
Иными словами, существование дополнительного пика, в частности низкочастотного пика 0.6 - 0.7 ГГц, свидетельствует о том, что молекулы липида не только перемещаются между потенциальными ямами с частотой 1 ГГц. Но также и между потенциальными ямами с другими характеристиками с частотами 0.7 ГГц. На основной режим накладывается существенно менее выраженный и более медленный режим движения.
Обратимся теперь к рассмотрению высокочастотной области спектра. На рис. 3.4, А различимы широкие пики в области 2.5 и 4 ГГц. Амплитуда этих пиков сильно падает с увеличением частоты, поэтому целесообразно рассмотреть те же спектры, построенные в полулогарифмических координатах (рис. 3.5). Нетрудно видеть, что в высокочастотной части спектра наблюдается набор периодически сменяющих друг друга максимумов и минимумов, с сильно затухающей с ростом частоты амплитудой.
Внешняя сила переменной величины
Рассмотрим более подробно проникновение молекулы или иона через бислой под действием внешней силы. В самом начале расчета молекулярной динамики скорость движения исследуемой частицы невелика, но с течением времени приложенное силовое постепенно разгонит частицу. Причем, в том случае если приложенная сила невелика, скорость стабилизируется за счет уравновешивания внешней силы силой вязкого трения. Однако, если приложить достаточно большую силу, то исследуемая частица при столкновениях будет отдавать соседним атомам достаточно кинетической энергии для локального перегрева ее окружения. Причем, как показали наши предварительные исследования, если сила будет достаточно велика, то перегрев не будет полностью компенсироваться термостатом.
Как следствие, в модельной системе начнется неконтролируемый нагрев, повышение давления и уменьшение плотности. Падение плотности, в свою очередь, повлечет за собой уменьшение трения и дальнейшее увеличение скорости движения исследуемой частицы под действием внешней силы, а значит и еще больший разогрев. Такой автокаталитический процесс приведет в конце концов к испарению модельной системы.
Очевидно, что для исследования проницаемости липидного бислоя необходимо использовать такую внешнюю силу, которая не приведет к катастрофическому перегреву, так как в противном случае в процессе МД расчетов произойдет разрушение модельной мембраны. Таким образом, необходимо использовать силу такой величины, чтобы, с одной стороны, вносимые ей возмущения не приводили к нарушению структуры бислоя. А с другой стороны, величина приложенной силы должна быть достаточной для того, чтобы за времена порядка наносекунд частица полностью преодолела бислой. Однако, в том случае если модельная система неоднородна и состоит из областей, в которых коэффициенты трения для исследуемой частицы различаются, как это имеет место в случае липидной мембраны в воде, невозможно указать какое-либо единственное значение приложенной внешней силы, удовлетворяющее поставленным условиям и подходящее для всех областей с различными коэффициентами трения.
Это связанно с тем, что небольшая сила, вполне подходящая для проведения расчетов в одной области модельной системы (в воде), может оказаться неприемлемо маленькой для другой (бислой). Соответственно, время, необходимое для проведения расчетов с такой силой в более плотной области может существенно превысить пороговое значение в 10 нс. Проведение таких длинных расчетов нецелесообразно ввиду их высокой вычислительной стоимости. В то же время сила, удовлетворяющая всем условиям в более вязкой среде, оказывается неприемлемо большой для менее вязкой части системы.
Для преодоления описанных затруднений необходимо использовать переменную силу, зависящую от характеристик среды, в которой находится исследуемая частица в данный момент времени. Однако, в явном виде задать такую зависимость не представляется возможным, так как именно свойства среды, то есть липидной мембраны, являются предметом исследования. Поэтому нельзя ввести в модель никакого дополнительного предположения относительно зависимости сопротивления мембраны проникновению частицы, не внеся тем самым в модель систематическую ошибку или даже основные черты полученного конечного "результата". Следовательно, необходимо проводить расчеты таким образом, чтобы на величину приложенной силы оказывали влияние факторы, возникающие локально, в процессе молекулярной динамики.
В качестве такого фактора мы решили использовать термостат Беренд-сена, процедуру отводящую излишки энергии у перегретых частей модельной системы. Процедура термостатирования была отдельно определена для модельной системы с липидным бислоем и водой и исследуемой частицы к которой прикладывалась внешняя сила. Причем, для пенетрирующей частицы была определена существенно более производительная процедура термостатирования, способная отводить больше лишней энергии, чем для остальной модельной системы. Это необходимо для эффективной компенсации разгона частицы в областях с низким значением коэффициента трения, особенно в случае использования больших по абсолютной величине сил. Технически увеличение производительности термостата было достигнуто уменьшением константы термостатирования до 0.01 пс.
В модельной системе с термостатом, организованным таким образом, в областях модельной системы, в которых коэффициент трения невелик процедура термостатирования будет тормозить излишне разгоняющуюся частицу. В то же время в областях системы с высоким коэффициентом трения и низкими скоростями движения влияние термостата будет меньше, что и обеспечит необходимое регулирование величины внешней приложенной силы. что константа ту -фт связан с тем, что при масштабировании скоростей изменение энергии перераспределяется между кинетической и потенциальной составляющей таким образом, что изменения в температуре системы оказываются меньше, чем изменения в энергии, что не соответствует задаче, поставленной перед термостатом [87, 113].
Проницаемость липидных мембран характеризуют коэффициентом проницаемости, который связывает поток вещества (J) с разностью его концентраций (АС) по разные стороны мембраны: J = L АС. В конце XIX века экспериментально было установлено, что коэффициент проницаемости L связан со средним внутримембранным коэффициентом диффузии проникающей частицы (D), ее коэффициентом распределения (а) и толщиной мембраны (I): L = [114, 115]. В наших расчетах толщина мембраны была принята равной 10 нм. Точное значение коэффициента распределения а может быть получено только в реальных экспериментах. Поэтому, для определения коэффициента проницаемости L в численных экспериментах ключевым является определение среднего внутримембранного коэффициента диффузии.
