Диссертация (1149183), страница 6
Текст из файла (страница 6)
По этой же29причине в системе используется модель воды TIP3P [112]. Для описаниямоновалентных и двухвалентных ионов были использованы стандартные параметрыполя Amber99. Так как используемое для вычислений программное обеспечениеGromacs 4.5.6 не включает последние версии Amber99, параметры Lipid14 иparmbsc0 были имплементированы в программный пакет.Ко всем системам были применены широко используемые и проверенныеметоды и подходы. Каждая система была помещена в ящик с периодическимиграничными условиями. Размеры ящика вдоль осей X и Y были заданы согласноразмерам рассматриваемого бислоя, помещенного в центре ящика перпендикулярнооси Z (~ 6 и 10 нм ).
Размер ящика по Z координате для системы “ДНК-бислой”варьировался от 10,5 до 13 нм. Cистемы были изучены в изотермо-изобарическомансамбле при температуре 303 K и давлении 1 Бар. Температура поддерживалась припомощимодифицированноготермостатаБерендсенасдополнительнойстохастической поправкой к кинетической энергии, обеспечивающей ее правильноераспределение (velocity-resсaling) [113]. Согласно методу Берендсена поддержаниетемпературы T0 системы описывается формулойdT T0 TdtT(2)где τT – параметр термостата, t – время.Таким образом, отклонение температуры T от T0 убывает экспоненциально.
Так кактемпература определяется кинетической энергией, то корректируются не саматемпература, а скорости частиц. На каждом шаге скорость частицы умножается нафактор λ 1t T0 1 T T (3)В случае модифицированного термостата Берендсена, используемого в работе,30кинетическая энергия корректировалась согласно следующей формулеdK K 0 K dtT2KK 0 dWN f T(4)где K – кинетическая энергия молекул в момент времени t, K0 – заданнаякинетическая энергия, Nf – число степеней свободы системы, dW – винеровскийпроцесс.Давление в системе контролировалось при помощи анизотропного баростатаБерендсена [114].
Данный метод поддерживает давление P0 в системе аналогичноуравнению (2)dP P0 PdtP(5)где τP – параметр баростата. Давление P является тензорной величиной, котороеопределяется какP2 K V(6)где V – объем ящика, Ξ – вириал сил. Как и в случае с термостатоммасштабированию подвергалось не само давление P, а координаты и линейныеразмеры ящика.ВзаимодействияЛеннарда-Джонсабылиобрезанына1нм,дляэлектростатических взаимодействий был применен современный метод particle-meshEwald (PME) [115]. Шаг по времени составлял 2 фс.С целью проверки корректности имплементации силового поля Amber99parmbsc0 было проведено моделирование молекулы ДНК в водном растворе вприсутствии противоионов натрия. Для проведения моделирования молекулы ДНК вводном растворе была использована схема, схожая с методикой подготовки уравновешенного состояния ДНК, отработанной исследователями в течение более 10 лет31[79, 116]. Методика, применённая в данной работе, включала в себя последовательную минимизацию энергии, нагревание и уравновешивание отдельных компонентовсистемы (ДНК и молекул воды) независимо друг от друга.
Процесс уравновешиваниясостоял из следующих шагов: 1) 20 пс моделирование при T = 100 K; 2) минимизацияэнергии воды; 3) минимизация энергии всей системы; 4) моделирование всей системы в течение 20 пс с наложением ограничений на положения атомов ДНК (алгоритмposition restraints) с силовой константой k = 41840 кДж*мол-1 нм-2 при T = 100 К; 5) тоже, что и шаг 4, но время моделирования 50 пс, T = 303 К; 6) то же, что и шаг 5, при k= 20920 кДж*мол-1 *нм-2; 7) то же, что и шаг 5, но k = 10460 кДж*мол-1 *нм-2; 8) тоже, что и шаг 6, но k = 2092 кДж*мол-1 *нм-2; 9). то же, что и шаг 6, но k = 84кДж*мол-1 *нм-2. Затем система была промоделирована в течение 50 пс без использования ограничений в изотермо-изобарическом (NPT) ансамбле. Конечная конфигурация была использована в качестве стартовой для моделирования в течение микросекунды (production run).Одной из главных характеристик проверки надежности силового поля являетсясреднеквадратичное отклонение (СКО) от начальной структуры (Рис.
2). Для додекамера в водном растворе было получено среднее значение СКО = 0,26 ± 0,05 нм длятраектории в 1000 нс, что хорошо согласуется с ранее полученными результатами[16].Что касается липидной мембраны, в работе использовался липидный бислой,состоящий из 128/288 молекул фосфатидилхолина (ПОФХ), являющегося одним изосновных компонентов биомембраны. Как было уже указано выше, молекула ПОФХсостоит из 2 гидрофобных хвостов и гидрофильной головной группы (Рис. 1,б). Дляпроверки имплементации силового поля Amber Lipid14 в программный пакетGromacs 4.5.6 липидный бислой был уравновешен в течение 200 нс, и в результатебыл измерен главный показатель правильной работы поля — площадь, приходящаяся32на один липид. Значение площади было равно 0.65 ± 0.01 нм2, что хорошо согласуется с ранее полученными результатами [17].Рис.
2. Среднеквадратичное отклонение додекамера ДНК от начальной В-формы.Таким образом, имплементация и валидация силовых полей в Amber99оказалась успешной и их можно применять для описания межмолекулярныхвзаимодействий системы ДНК-липидный бислой.Для расчета профиля свободной энергии (ПСЭ) в работе использовался методзонтичной выборки, основанный на приложении гармонического потенциала междуцентрами масс ДНК и бислоя вдоль координаты Z. Для этого додекамер ДНК былпомещен параллельно поверхности бислоя (перпендикулярно оси Z).
Начальноерасстояние между центрами масс липидного бислоя и ДНК было равно 4,8 нм. Нацентр масс ДНК действовала сила, направленная к центру масс липидного бислоявдоль координаты Z (силовая константа k = 1000 кДж·моль-1·нм-2). В результатетакого моделирования было получено 35 конфигураций системы на расстоянии 0.133нм друг от друга (таким образом, расстояние между центрами масс варьировалось от4,8 нм до 1,4 нм с шагом 0,1 нм вдоль оси Z).
Каждая конфигурация служиланачальной для независимых моделирований методом зонтичной выборки в течение100 нс (силовая константа k = 3000 кДж·моль-1·нм-2), первые 20 нс из которыхрассматривались как уравновешивание системы и не использовались для анализа.Профиль свободной энергии вдоль оси Z был вычислен при помощи метода анализавзвешенныхгистограмм(WeightedHistogramAnalysisMethod,WHAM),имплементированного в пакет Gromacs [117].Что касается анализа полученных траекторий, одной из важных характеристикявляется расчет количества контактов (ненормированное координационное число).Количество контактов было сосчитано согласно методике, описанной в статье [48].Сначала рассчитывалась радиальная функция распределения для пары атомов.Следующие пары были рассмотрены: Nфх (атомы азота холиновых групп ФХ) −Pднк (атомы фосфатов ДНК), Ca − Pфх (атомы фосфатов ФХ головных групп), Ca −Pднк, Na − Pфх, и Na − Pднк.
Далее были определены положения первых минимумоврадиальной функции распределения; эти минимумы обозначили радиус первойкоординационной сферы. Были получены следующие радиусы: 0.45 нм для пар Ca −Pфх и Ca − Pднк, 0.42 нм для пар Na − Pфх и Na − Pднк, и 0.6 нм для пары Nфх −Pднк. (Рис. 3). Количество контактов между атомами A и B рассчитывалось поколичеству атомов B в первой координационной сфере атомов A. Считалось, чтоионы адсорбировали на мембрану, когда они попадали в первую координационнуюсферу атома Pфх липида. Кроме того, считалось, что “Ca-мостик” (“Na-мостик”)между фосфатными группами ДНК и липидных молекул образуется, когда оба атомаPфх и Pднк находятся внутри координационной сферы одного и того же иона Ca2+(Na+).34а)б)Рис 3. Радиальная функция распределения в зависимости от расстояния междуатомами: а) Ca-Pфх и Ca-Pднк; б) Nфх-Pднк.35Анализ динамических и структурных характеристик липидного бислоя включал в себя расчет таких характеристик как площадь, приходящаяся на один липид,ориентация вектора, соединяющего атомы Pфх и Nфх (PN вектор) по отношению квектору нормали к поверхности мембраны, параметр порядка, коэффициент латеральной диффузии.
Параметр порядка SCD для углеводородных групп гидрофобныххвостов липидов дает информацию об изменении структуры внутренней областибислоя и рассчитывается следующим образом:SCD 13 cos 2 1 ,2(7)где θ – угол между CH связью и нормалью к мембране.Для изучения мобильности липидов использовался коэффициент латеральнойдиффузии DL, который был рассчитан согласно описанной в статье [118] методике:1[r (t )]2t 4tDL lim(8)2где [r (t )] – среднеквадратичное смещение, t – время.В данной работе последние 100 нс моделирования были разбиты на 10 частей.Для каждого участка рассчитывалось среднеквадратичное смещение, коэффициентдиффузии DL определялся как наклон прямой среднеквадратичного отклонения. Дляувеличения точности коэффициент диффузии был измерен отдельно для каждого 10нс участка траектории.
Конечное значение коэффициента диффузии липидов былополучено усреднением коэффициентов каждого участка на всей 100-нс траектории.Длявизуализациисистемыиспользовалсясовременный программный пакет VMD [119].многофункциональный36ГЛАВА 3. Взаимодействие ДНК с липидным бислоем ПОФХДля изучения взаимодействия молекулы ДНК с липидным бислоем первымшагом было проведение компьютерного моделирования системы без двухвалентныхионов, в которой фрагмент ДНК помещен в раствор параллельно поверхностибислоя, состоящего из липидных молекул ПОФХ. Начальное расстояние между ДНКи липидной мембраной определялось как расстояние между ближайшими атомамиДНК и липидов вдоль вектора нормали к поверхности мембраны.















