Влияние катионов на структурные и электрические свойства липидного бислоя. Молекулярно-динамическое исследование (1102545), страница 10
Текст из файла (страница 10)
Во-вторых, механизм его срабатывания вслучае с мембранами позволяет достигать равновесия в липидном бислое и в двойномэлектрическом слое на более коротких временах. Полуассиметричный баростат, поддерживающий поверхностное натяжение (σ), был выбран нами из-за несовершенствасилового поля для липидов: для получения адекватных экспериментальным значенийплощади на липид, на мембране в некоторых случаях поддерживалось небольшое по-– 44 –стоянное поверхностное натяжение, величина которого подбиралась в зависимостиот высоты вычислительной сцены. Тиелеман и Берендсен анализировали поведениебислоя из DPPC при различных граничных условиях, NσT (постоянное поверхностноенатяжение), NpT (постоянное давление) и NVT (постоянный объем); авторы сделаливывод о допустимости использования NσT и NpT условий [5].
Отметим, что данныеусловия не являются естественными для мембран в ненапряженном состоянии [85],однако расчет в таких условиях допустим, когда механические свойства мембраныне являются предметом наблюдения. В расчетах, направленных на изучение механохимических свойств липидного бислоя, использовались граничные условия NpST(стандартные условия NPT и постоянная площадь на липид).Название параметра расчетаВид уравненийИспользование голономных ограниченийШагКоррекция дальнедействующего дисперсионного взаимодействияЭлектростатические взаимодействияВан-дер-ваальсовы взаимодействияТермостатированиеБаростатированиеМодель водыЗначение параметраСтохастическая динамика, sdбез ограничений либо фиксировалисьвалентные связи с атомами водорода0.5 или 2 фсвключена для U и p [86]PME [87]switch-функция с длинами 9 и 11 ÅТермостат Нозе-Хувера [88]полуассиметричныйбаростатБерендсена [56]TIP4P («гибкая») [89]Таблица 2.2.1: Стандартные параметры расчета молекулярной динамики, используемые в работе приполучении равновесных траекторий системы липидный бислой–вода.Вопрос о предпочтительности использования дальнодействующей электростатики(PME) является спорным: с одной стороны в ряде работ убедительно показано, чтообрезание электростатических взаимодействий приводит к различного рода артефактам [58], с другой стороны PME усиливает искусственную периодичность системы,что само по себе может приводить к актефактам другого рода [90, 91].
В большинстве наших расчетов использовалось суммирование PME, и мы полагаем, что вкомбинации с достаточно крупным размером системы нам удалось избежать артефактовпериодичности.Протокол работы.Как было описано выше, сборка начальной геометрии липидногобислоя осуществлялась при помощи скрипта memgen. После этого запускалась минимизация геометрии до достижения максимальных сил порядка 102 кДж·моль−1 ·нм−1 .Далее в свободный объем добавлялись молекулы растворителя. На этой стадии частьмолекул растворителя заменялась на ионы для достижения нужной ионной силы и– 45 –AB500,/1000050010000204060,801000204060,80100Рис 2.2.1: Изменение потенциальной энергии при уравновешивании системы DPPC:DPPS 4:1 из исходнособранного состояния (А) и в ходе продуктивного расчета той же системы.снова производилась минимизация геометрии.Достижение стационарного состояния фиксировалось по выходу на плато следующих параметров: потенциальной энергии системы, площади на липид, компоненттензора давления (см.
рис. 2.2.1А). Это состояние интерпретировалось здесь идалее как состояние термодинамического равновесия. Особенности баростата приводят к появлению некоторых небольших осцилляций плотности упаковки липидов,происходящих вокруг равновесного положения (рис. 2.2.1B). Эти небольшие осцилляции иногда влияют на высоту липидного бислоя, что потребовало использованияпослойного суммирования на этапе обработки траекторий (см.
раздел “Привязка котносительным координатам”).Объекты вычислительного эксперимента.В качестве объектов использовались бис-лойные мембраны из DPPC/DPPS в различных соотношениях — 1:0, 4:1, 3:2, 2:3 и0:1 — размером 90 липидов на монослой, с водной прослойкой порядка 6–7 нм скаждой стороны. В некоторых вычислительных экспериментах использовались бислоиразмером 25 липидов на монослой и порядка 2 нм водной прослойки с каждой стороны.Моделирование жидко-кристаллической (Lα ) фазы в таких системах производилось притемпературе 340 K с постоянным «поверхностным натяжением» порядка -20 дин/см.L′β -фазы в таких системах моделировались при температуре 300 К в “полуизотропном”баростате, поддерживающем давление 1 Бар в латеральном и нормальном направле-– 46 –ниях. В ходе вычислительных экспериментов по адсорбции лизина, БЛМ из DOPS (90липидов на монослой) моделировались в полуизотропном баростате (1 Бар) при температуре 300 K.
При исследовании механических свойств бислоя DMPS (90 липидовна монослой), моделирование производилось при температуре 300 К при постояннойплощади, давление поддерживалось только в направлении нормали к поверхности.2.3. Методы обработки траекторийМетоды обработки траекторий молекулярной динамики являются одним из ключевых элементов настоящего исследования. Разработанные методы хорошо подходят для самых различных вычислительных экспериментов, объектом которых является бислойная липидная мембрана или любой другой относительно плоский объект.Реализация утилит обработки траектория выполнена с использованием библиотекMDAnalysis [92], а также в некоторых случаях — библиотек GROMpy9 .2.3.1.
Привязка к относительным координатамКак было сказано выше, мембрана может претерпевать некоторые осцилляциивокруг равновесного положения, которые могут влиять на распределение измеряемыхвеличин (например плотностей) в абсолютных координатах. Очевидно, это явлениеперестает играть роль, если перейти к координатам относительно поверхностимембраны. Поскольку, в основном, нас интересовали распределения величин вдольнормали к мембране, то точку отсчета необходимо выбрать вдоль оси Z. В работеза относительный нуль выбрано среднее положение Cα атома глицерина фосфолипида,как в работе [4]. В простом случае некоторая величина Ω усреднялась в каждомi-слое по формуле:Ωi = ⟨Ω(zCα + i∆z; zCα + (i + 1)∆z)⟩f(2.3.1)где i – номер слоя, ∆z – толщина слоя по z, f означает суммирование по фреймамтраектории.
Специальное усреднение производилось в случае работы с «рифленой»P-фазой, когда мембрана имела неровную форму в ZX или ZY-плоскости. В этомслучае XY-плоскость разбивалась на сетку, в каждой ячейке которой вычислялосьзначение zCα , далее координаты всех атомов над каждой ячейкой корректировалисьв соответствии с этим значением — таким образом достигалось математическое«выпрямление» поверхности мембраны.9Библиотеки GROMpy представляют python-интерфейс к C-библиотекам GROMACS и на данный момент неподдерживаются.– 47 –Поскольку, в данной работе рассматривались симметричные бислои, параллельнос (2.3.1) производилось усреднение по двум монослоям мембраны.2.3.2.
Вычисление коэффициентов диффузииВ работе [93] обсуждаются различные подходы к вычислению значения коэффициента самодиффузии воды, и делается вывод о наиболее адекватном методе расчетачерез формулу Грина-Кубо, использующую автокорреляционную функцию (АКФ) скорости:1D=N∫+∞⟨v(0)v(τ )⟩ dτ0где v — скорость молекулы, N — размерность пространства диффузии. Формально интегрирование производится до бесконечности, что недоступно в численномэксперименте и не является необходимым.
Также формула симметрична относительно нижнего предела интегрирования, поэтому для численного эксперимента можновыписать уточненную статистическую формулу:⟨D=1N∫t+∆t⟩⟨vi (t)vi (t + τ )⟩ dτt,t,iгде индекс i — означает усреднение по множеству частиц, t — начальное время, откоторого строится АКФ, ∆t — соответствует интервалу, на котором вычисляется этафукнция. Параметры подбираются эмпирически в зависимости от скорости затуханияАКФ.
Однако для такого расчета необходим достаточно маленький шаг записи скоростей (в работе использовалась 1 фс), тогда как продуктивные траектории мембрандостаточно большие и зачастую выписываются с низкой частотой. Кроме того, надиффузию воды различные факторы могут влиять с большим характерным временем. Например, молекулы могут захватываться поверхностью, потом снова выходить в объем.Их коэффициент диффузии, таким образом будет снижаться на больших временах.Вторым методом, обсуждаемом в работе [93], широко применяющимся при обработки различных вычислительных экспериментов, является расчет по формуле Эйнштейна:D=|r(τ ) − r(0)|21lim,2N τ →∞τгде r(τ ) — радиус-вектор наблюдаемой частицы.
По теореме об эргодичности мы можемперейти от среднего по времени к среднему по ансамблю и, вместо того, чтобынаблюдать за удалением одной молекулы бесконечное время, наблюдать за смещением– 48 –множества молекул за конечное время:⟨D = limτ →∞|ri (t + τ ) − ri (t)|22N τ⟩(2.3.2)t,iТакой предел должен сходиться существенно быстрее.
Однако определение локальногокоэффициента самодиффузии невозможно при наблюдении на очень больших временах,поскольку молекулы уходят из наблюдаемого локуса. Поэтому для анализа подвижности воды в слоях, находящихся на разном удалении от поверхности, мы используемформулу:⟨Dτ =|ri (t + τ ) − ri (t)|22N τ⟩;τ = const(2.3.3)t,i:z0 <zi (t)<z1Значение коэффициента локальной самодиффузии Dτ оказывается несколько заниженным, чем коэффициента самодиффузии, вычисленного по формуле (2.3.2).
Однако,этот параметр хорошо описывает изменение подвижности воды в ДЭС на различномудалении от внешней плоскости Гельмгольца. Поэтому в настоящем исследованииприменялся именно этот метод расчета.2.3.3. Вычисление электростатического потенциалаЛокальный электростатический потенциал вычислялся методом двойного интегрирования зарядовой плотности, построенной в относительных координатах (2.3.1).Формула для вычисления потенциала в системе СГС:∫ψ(z) = ψ(−∞) −zdz−∞′∫z′−∞ρ(z ′′ ) ′′dzε(2.3.4)Термодинамический потенциал ионов, имеющий смысл свободной энергии ∆G,может быть рассчитан как потенциал средней силы (ΨP M F ) исходя из распределенияконцентрации (2.3.5). Этот потенциал включает в себя по мимо электростатического вклада энергию дисперсионного взаимодействия, силы изображения, эффекты изменения ионной атмосферы и изменения состояния воды вблизи поверхности,привнесенные ионами [24].
В вычислительном эксперименте оказывается возможнымвычислить полный потенциал взаимодействия, который имеет смысл энтальпии (∆H).Этот потенциал можно найти, рассчитывая непосредственно энергию взаимодействия(в основном, электростатического) каждого иона с соседними ионами, молекуламиводы и липидами (2.3.6).Важно отметить, что согласно определению этих потенциалов, потенциал взаимодействия отличается от потенциала средней силы тем, что в– 49 –первом случае не учитывается энтропийная составляющая. Приведенные ниже формулывыражают эти потенциалы в единицах электрического потенциала.zRTC(x)lnFC(∞)(2.3.5)z(Uion−atmosphere + Uion−water + Uion−membrane )F(2.3.6)ΨP M F =Ψint =Здесь z — зарядовое число иона, x — расстояние от выбранного нами начала координат, C — локальная концентрация иона в точке, C(∞) — концентрация иона вобъемной фазе, а R,T,F имеют общепринятое значение.Техническе детали вычисления потенциала взаимодействия.Вычисления потенциала взаимодей-ствия иона с окружением вычислялось по следующей схеме: для N ионов подготавливалось N различныхindex-файлов, в которых были прописаны группы [ion] (в ней находился один i-й ион) и [atmosphere](в ней находились все ионы, кроме i-ого).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
