Study of the dna denaturation based on the peyrard-bishop-dauxois model and recurrence quantification analysis

Abstract


The aim of this work was is to study of the phase transition during thermal denaturation of DNA by the method of recurrence quantification analysis, which is the most modern and efficient approach for the analysis of complex systems. To achieve the goal, the nonlinear dynamics model of the DNA denaturation is implemented, based on the Peyrard-Bishop-Dauxois approach, which allows determining the state of DNA base pairs at an arbitrary point in time. Next, the dynamics of DNA is simulated at different temperatures, and the evolution of hydrogen bond stretching during DNA heating is studied. To analyze the results obtained, the method of recurrence quantification analysis is used, which make it possible to judge the nature of the processes occurring in the system and are of particular interest in relation to the study of nonlinear DNA dynamics, in particular phase transitions, since on the basis of their analysis it is possible to identify states of repetition and fading (laminarity), the occurrence of extreme events, the presence of hidden periodicity and cyclicity. Using this method of analysis, quantitative measures are calculated in a window shifted along the main diagonal of the recurrent plot (RP). These measures quantify fine-scale RP structures. The calculations reveal the most sensitive measures to the analysis of the process under study, with the help of which the precursors of the phase transition are determined.

Full Text

Плавлением ДНК (термической денатурацией) называется инициированный нагревом процесс разрушения комплементарных водородных связей (обусловленных, так называемым стэкинг-взаимодействием) в ДНК при сохранении ковалентных связей. Как показано в ряде экспериментальных и теоретических работ [2; 6; 8; 9; 14; 15; 17; 20; 21], плавление ДНК сопровождается структурными превращениями, имеющими природу фазовых переходов первого или второго рода. Исследования возможности существования фазовых переходов в одномерных системах, к коим относится ДНК, является актуальной проблемой теоретической физики [5; 10]. Впервые наличие фазового перехода в растворе гомогенной ДНК (т.е. молекулах, состоящих из последовательности одинаковых пар оснований, нуклеотидов) экспериментально было установлено при измерении поглощения света с длиной волны 260-268 нм при медленном нагреве [8]. Дальнейшие исследования нативных ДНК молекул показали, что денатурация происходит в несколько этапов и очень чувствительна к составу последовательности [21]. Теоретические исследования, подтвердившие возможность фазовых переходов в ДНК, были выполнены в работах Poland et al. [14], Scheraga et al. [15], Causo et al. [2], Пейраром и соавт. [3; 4]. При этом наиболее эффективной оказалась динамическая модель Пейрара - Бишопа - Доксуа и ее модификации, уточняющие описание процесса денатурации ДНК, многообразие которых подробно описано в обзоре [18]. Метод количественного анализа рекуррентных диаграмм (РД) применяется для изучения нелинейных процессов, характерных в том числе для биологических систем. Рекуррентные диаграммы позволяют судить о характере протекающих в системе процессов и представляют особый интерес применительно к изучению нелинейной динамики ДНК, в частности фазовых переходов, так как на основе их анализа можно выявить состояния «повторения и замирания» (recurrence and laminarity states), совершение экстремальных событий, наличие скрытой периодичности и цикличности. Количественный анализ рекуррентных диаграмм позволяет определить «численные меры», основанные на плотности рекуррентных точек, распределениях длин диагональных, горизонтальных или вертикальных линий, характеризующих эволюцию динамической системы [16; 22]. Целью работы является исследование фазового перехода в ДНК, происходящего при термической денатурации, методом количественного анализа рекуррентных диаграмм. Данная работа состоит из пяти частей. Во введении обсуждаются вопросы актуальности и новизны исследования процесса денатурации ДНК на основе модели Пейрара - Бишопа - Доксуа и количественного анализа рекуррентных диаграмм. Первый раздел посвящен описанию нелинейной динамической модели процесса денатурации ДНК, основанной на подходе Пейрара - Бишопа - Доксуа. Во втором разделе представлены результаты моделирования процесса денатурации ДНК на основе модели Пейрара - Бишопа - Доксуа. Третий раздел посвящен изучению результатов моделирования методом количественного анализа рекуррентных диаграмм. В заключении обсуждаются полученные результаты, делаются выводы и обозначаются перспективы дальнейших исследований. Нелинейная динамическая модель процесса денатурации ДНК Оригинальная модель Пейрара - Бишопа - Доксуа была предложена в работах [3; 4] для описания коллективных движений пар оснований, включая флуктуации большой амплитуды. Модель Пейрара - Бишопа - Доксуа является развитием модели Пейрара - Бишопа, (ПБ-модель), предложенной для описания денатурации ДНК. В ПБ-модели взаимодействие между парами оснований (водородные связи) описывается потенциалом Морзе; стэкинг-взаимодействия вдоль цепи учитываются через гармонический потенциал. Полная энергия системы имеет вид (1) где Ekin - кинетическая энергия, Uhyd - потенциал Морзе, Ustack - потенциал стэкинг-взаимодействия. На рис. 1 приведено схематичное представление динамической модели нелинейной динамики ДНК с указанием водородных связей Uhyd и стэкинг-взаимодействий Ustack. Кинетическая энергия вычисляется согласно следующему соотношению (2) где m - масса пары оснований, pn - импульс пары оснований, n - количество пар оснований. Потенциальная энергия представляет собой сумму потенциала Морзе, который описывает энергию взаимодействия водородных связей при изменении расстояния между парами оснований и стэкинг-взаимодействий между соседними парами оснований вдоль цепи. Потенциал Морзе имеет следующий вид (3) где D - энергия диссоциации пары оснований, α - пространственный масштаб потенциала, yn - переменная, описывающая динамику оснований внутри пары. Стэкинг-взаимодействия между соседними основаниями ДНК описываются гармоническим соотношением (4) где k - параметр гармонического взаимодействия между соседними парами оснований, h - индекс, обозначающий гармоническую связь. Рис. 1. Динамическая модель нелинейной динамики ДНК [3] Для описания движения оснований рассматривается производная от функции Гамильтона: (5) Откуда можно получить соотношение для импульса основания . (6) Для адекватного моделирования необходимо учитывать термические флуктуации ДНК, поскольку ДНК находится в ядрышке клетки при определенной температуре, определяемой окружением. Один из способов учесть температуру - это введение термостата Нозе - Гувера, который будет симулировать реальные условия. Так, вводится новая переменная ξ, представляющая собой термодинамический коэффициент трения ξ = sps/M, где s - дополнительная степень свободы, действующая как внешняя сила на систему из n частиц. Таким образом (7) где М - масса термостата. Производная от сопряженного импульса имеет вид: (8) где - константа Больцмана, Т - температура термостата. Тогда уравнение движения ДНК будет иметь вид (9) Переменная ξ стабилизирует кинетическую энергию реальной системы до канонического математического ожидания . Модель Пейрара - Бишопа позволила сравнить результаты моделирования с экспериментом, но численные расчеты давали значения температуры плавления приблизительно на 150 К выше экспериментальных данных. В связи с этим М. Пейраром, А. Бишопом и Т. Доксуа был введен ангармонический потенциал, описывающий нелинейный характер стэкиг-взаимодействия: (10) где ρ - безразмерный параметр, через который вводится нелинейность процесса стэкинг-взаимодействия, χ - коэффициент затухания, a - индекс, обозначающий ангармоничность потенциала. Следствием введенного «ангармонизма» потенциала явилось резкое увеличение расстояния между молекулами ДНК при возрастании температуры, что является качественной особенностью модели при описании динамики ДНК. Последнее предполагает учет роли термических флуктуаций на развитие нелинейных процессов в ДНК. Один из способов учета температуры - это введение термостата Нозе - Гувера [19]. С учетом ангармонизма потенциала уравнение движения ДНК принимает следующий вид: (11) где - константа Больцмана, Т - температура термостата, М - масса термостата. Таким образом, модель Пейрара - Бишопа - Доксуа представляет собой нелинейную динамическую модель денатурации ДНК, позволившую привести в соответствие результаты расчетов с экспериментальными данными. Моделирование процесса денатурации ДНК Несмотря на то, что существуют модификации модели Пейрара - Бишопа - Доксуа, учитывающие гетерогенность участка ДНК [1; 7; 12], в данной работе рассматривается гомогенный случай. С использованием модели Пейрара - Бишопа - Доксуа проведено численное моделирование процесса нагрева ДНК в диапазоне от 200 до 540 К для 256 пар оснований. На рис. 2 представлена эволюция растяжения водородных связей при нагревании. На данной диаграмме изображена динамика сегмента ДНК. Случайный характер флуктуаций связан с видом начальных условий, а именно равномерным случайным распределением открытых и закрытых состояний пар оснований. Белые пятна на рисунках соответствуют областям, где растяжение между парами оснований на 1,2 Å и более превосходит допустимое значение, равное 10 Å, в течение ограниченного времени, когда пара оснований считается открытой. Рис. 2. Эволюция растяжения водородных связей при нагревании ДНК. Черный цвет соответствует случаю, когда пара оснований закрыта, белый - когда открыта Рис. 3. Диаграмма плавления цепи ДНК Черные области, в свою очередь, обозначают области, где растяжение меньше допустимого, и пара оснований предполагается закрытой. Рассматривая каждую пару оснований отдельно, можно заметить значительные различия амплитуд растяжения. Темные области разделяются белыми областями, и подобное явление объясняется появлением пузырьков денатурации, обнаруженных экспериментально. При увеличении температуры до максимальной становится всё больше белых областей, это означает, что ДНК полностью денатурирует, то есть связи разрываются. На рис. 3 приведена осредненная по 1000 реализаций диаграмма плавления ДНК. Сплошная линия соответствует ангармонической модели, а пунктирная - гармонической. Можно заметить, что температуры плавления, полученные на основе данных моделей, различаются. Данный результат согласуется с результатами эксперимента и оригинальной работой [4], в которой предложена модель. Количественный анализ рекуррентных диаграмм процесса денатурации ДНК Рекуррентный анализ является современным и быстро развивающимся подходом к анализу сложных систем, а также является мощным инструментом для исследования повторяемости в динамических системах. С помощью данного метода можно охарактеризовать поведение системы в фазовом пространстве. Среди использующих данный метод исследований, имеется большое количество работ в области анализа нелинейной динамики биологических систем. С помощью рекуррентных диаграмм (РД) появляется возможность анализа динамических процессов при наличии состояний «повторения и замирания», влияния шума. Количественный анализ РД позволяет сопоставить полученной диаграмме некоторые количественные меры, основанные на распределениях длин линий и плотности, так называемых, рекуррентных точек. Главными преимуществами количественного анализа РД являются наглядность полученных результатов, а также возможность извлечь полезную информацию даже для коротких и нестационарных данных, когда другие методы нелинейного анализа терпят неудачу, например, метод фазовой плоскости, метод марковских процессов или метод кусочно-линейной аппроксимации [12]. Количественный анализ применим к различным видам данных и широко используется при решении инженерных задач, физиологии, геофизики, экономике и медицине [11; 12]. Для количественного анализа рекуррентных диаграмм процесса денатурации ДНК необходимо рассчитать рекуррентные матрицы, являющиеся отображением траектории эволюции ДНК в m-мерном фазовом пространстве (в работе будет использовано одномерное фазовое пространство, соответствующее динамике одной пары оснований) на двумерную двоичную матрицу размером N×N. Единица в ячейке матрицы соответствует проходу траектории в окрестности одной и той же точки фазового пространства при некотором времени i в некоторое другое время j. Обе координатные оси диаграммы являются осями времени. Математически это можно выразить следующим образом (12) где yi - переменная, в случае данного исследования описывающая динамику оснований внутри i-й пары, N - количество рассматриваемых состояний yi, εi - размер окрестности точки yi, Θ - функция Хэвисайда, m - размерность фазового пространства, ||.|| - Евклидова норма. Размер окрестности εi определяет радиус в фазовом пространстве с центром в точке yi. Если точка yj попадает внутрь данной окрестности, то такое состояние считается подобным состоянию yi, и, таким образом, на диаграмме устанавливается Ri,j =1 (i=j). Если значение ε будет малым, то малой будет последовательность точек повторения, и следовательно, «рекуррентная структура» исследуемой системы будет неясна. Если же значение будет выбрано слишком большим, тогда почти каждая точка будет являться соседней для любой другой точки, что приведет к ошибочным представлениям о структуре системы [13]. Поэтому для анализа динамики пары оснований были рассмотрены различные ε, а именно ε ≥ 5σ, ε ≥ 0,5σ, ε = σ, где σ - стандартное отклонение сигнала. На основе качественного анализа рекуррентных диаграмм было выявлено, что наиболее оптимальным является ε = σ. На рис. 4, б, приведена рекуррентная диаграмма динамики одной пары оснований, где по оси абсцисс указано время (относительные временные единицы (о.в.е.)), а по оси ординат - расстояние внутри пары. Динамика растяжения данной пары отражена на рис. 4, а. Для построения данной диаграммы были приняты следующие параметры: N = 100, ε = 0,31. Стоит отметить, на диаграмме присутствует черная диагональная линия, которая называется линией идентичности. Также стоит заметить, рекуррентная диаграмма симметрична относительно линии идентичности, то есть Ri,j ≡ Rj,i. а б Рис. 4. Анализируемая пара оснований: а - динамика растяжения пары оснований; б - рекуррентная диаграмма динамики, по оси абсцисс и ординат - время, в течение которого происходит нагрев (о.в.е.) С точки зрения топологии и текстуры можно выделить несколько видов рекуррентных диаграмм. Если на диаграмме видны белые полосы, это означает нестационарность данных, некоторые состояния редки или далеки от нормы. Если прослеживаются единичные точки, значит в системе происходят сильные колебания. Возникновение отдельных изолированных точек означает, что процесс может быть некоррелированным, случайным или антикоррелированным. Наличие вертикальных/горизонтальных полос показывает, что некоторые состояния не меняются или изменяются медленно в течение некоторого времени; индикация ламинарных состояний [13]. Проанализируем рекуррентную диаграмму (РД), представленную на рис. 4, б. В окрестности 900 относительных временных единиц (о.в.е.) происходит смена режима. Проанализируем часть диаграммы до данной точки с точки зрения топологии и текстуры. На диаграмме четко видны белые и черные вертикальные полосы. Их наличие означает нестационарный характер динамики пары оснований ДНК, система не находится в состоянии равновесия. При этом стоит отметить, что на коротких промежутках времени динамика пары оснований практически не меняет свой характер. Чем ближе к отметке 900 о.в.е., тем заметнее становятся единичные изолированные точки, это означает, что в системе происходят сильные колебания. На основе анализа топологии и текстуры рекуррентной диаграммы можно предположить, что наиболее чувствительными к анализу плавления ДНК являются меры, основанные на вертикальных/горизонтальных структурах, по причине явного выделения на диаграмме вертикальных линий. Для перехода к количественному анализу рекуррентных диаграмм необходимо ввести меры, которые количественно определяют мелкомасштабные структуры рекуррентных диаграмм. Данные меры основаны на плотности точек рекуррентности или структурах линий на диаграммах. Меры могут быть сгруппированы по структурам, на основе которых они рассчитаны. К мерам, основанным на плотности точек, относится мера рекуррентности (RR). Меры, рассчитанные на основе диагональных структур: мера дивергенции (DIV), мера детерминизма (DET), мера отношения (RATIO), мера средней длины диагональной линии (L) и мера энтропии (ENTR). Меры, основанные на вертикальных/горизонтальных линиях: мера замирания (LAM), мера времени запаздывания (TT) и мера максимальной длины вертикальной линии (Vmax). Мера рекуррентности определяется как . (13) На рис. 5 представлен результат расчета меры рекуррентности на основе РД (см. рис. 4, б). Данная мера показывает плотность рекуррентных точек, определяемая их количеством. Рис. 5. Мера рекуррентности, основанная на РД растяжения пары оснований Рассмотрим меры, основанные на диагональных структурах. Мера дивергенции задается следующим соотношением (14) где Lmax - максимальный размер диагональной линии. а б в г Рис. 6. Результаты расчета мер, основанные на РД растяжения пары оснований: а - мера дивергенции; б - мера детерминизма; в - мера отношения; г - мера энтропии Мера детерминизма имеет вид (15) где l - длина диагональной линии, lmin - минимальная длина диагональной линии, P(l) - количество диагональных линий размера l. Меры отношения и энтропии могут быть определены согласно , (16) (17) (18) На рис. 6 представлены результаты расчета мер, основанных на диагональных структурах, по РД динамики пары оснований ДНК (см. рис. 4, б). Мера дивергенции основана на определении размера самой длинной диагональной линии рекуррентной диаграммы. Данная мера связана с расходимостью траектории фазового пространства. Чем больше значение данной меры, тем быстрее расходятся сегменты траектории. Мера детерминизма отражает предсказуемость поведения ДНК. Диагональные структуры показывают время, в течение которого участок фазовой траектории подходит близко к другому участку фазовой траектории. Мера а б в Рис. 7. Результаты расчета мер, основанные на РД растяжения пары оснований: а - мера времени запаздывания; б - мера замирания; в - мера максимальной длины вертикальной линии отношения помогает выявлять переходы в динамике исследуемых молекул. Мера энтропии отражает вероятность p(l) нахождения диагональной линии длины l в рекуррентной диаграмме. Данная мера количественно определяет информативность случайных процессов. Также мера энтропии показывает сложность рекуррентной диаграммы по отношению к диагональным линиям в системе. Для диаграмм, в которых преобладают линии одинаковой длины, значение меры энтропии будет низкое, в ином случае, когда линии разных длин появляются одинаково часто, ENTR будет иметь довольно высокое значение. Мера времени запаздывания, основанная на вертикальных/горизонтальных структурах имеет вид (19) где ν - длина вертикальной линии, νmin - минимальная длина вертикальной линии, P(ν) - количество диагональных линий длины ν. Мера ламинарности определяется как (20) где знаменатель определяет общее количество точек в рекуррентной диаграмме. Максимальная длина вертикальной линии на РД задается (21) где Nv - абсолютное количество вертикальных линий. На рис. 7 представлены результаты расчета мер, основанных на вертикальных структурах, по РД динамики пары оснований ДНК (см. рис. 4, б). С помощью меры времени запаздывания определяется средняя длина вертикальных структур. Также данная мера позволяет оценить среднее время, в течение которого ДНК будет находиться в определенном состоянии, или как долго это состояние будет заблокировано. Мера замирания определяется отношением количества рекуррентных точек, образующих вертикальные линии, к общему количеству рекуррентных точек. Величина LAM характеризует наличие состояний замирания исследуемой молекулы. С помощью меры максимальной а б в Рис. 8. Графики эволюции мер, основанные на РД растяжения первой пары оснований с использованием окна размером S = 20: а - мера времени запаздывания; б - мера максимальной длины вертикальной линии; в - мера замирания длины вертикальной структуры можно определить максимальную длину вертикальной линии. Фазовый переход исследуется при различных значениях температуры, которая варьируется в диапазоне от 200 до 540 К. Расчёт рекуррентных свойств ДНК выпол-няется для центральной пары оснований. По рис. 4 проведен расчет мер в окне размером S = 400, выбор которого установлен экспериментально и вычислен по формуле S ≈ N/5. Важным моментом является выбор размера окна. Слишком малый размер окна может повлечь за собой излишний шум и флуктуации, а слишком большой - даст более гладкие графики мер ввиду незаметного влияния событий во всем объеме данных. Вычисление мер в окне аналогично вычислению мер для всей рекуррентной диаграммы и позволяет проследить динамику изменения мер в зависимости от времени (эпохи). Полученные меры реагируют на смену структуры диаграммы. Действительно, в окрестности координаты 900 о.в.е. на диаграмме видна резкая смена структуры, что говорит о смене режима. Меры RR, DET, RATIO, ENTR и TT отреагировали на это возникновением локальных экстремумов. Вычисление данных мер позволяет установить факт смены режима. Характер перехода лучше всего демонстрируют меры TT, LAM, Vmax, поскольку расчёт данных мер основан на вертикальных структурах. Поскольку при анализе топологии и структуры рекуррентной диаграммы результатов моделирования денатурации ДНК четко видны именно вертикальные структуры, значит именно меры, основанные на них, будут чувствительны к анализу исследуемого процесса. Рассмотрим процесс более детализировано, уменьшив размер окна при расчёте мер, а также рассчитав меры для сигнала других частиц. На рис. 8 приведены графики эволюции мер TT, LAM, Vmax, рассчитанные для динамики одной пары оснований с использованием окна размером S = 20. Затем приведем графики эволюции мер TT, LAM, Vmax, рассчитанные для динамики другой пары оснований с использованием окна с таким же размером, S = 20. На рис. 8, 9 можно заметить явные предвестники фазового перехода. В частности, по графикам эволюции меры ТТ можно увидеть непосредственно перед фазовым переходом снижение значения амплитуды, а затем резкий всплеск. По графикам эволюции меры Vmax - увеличение промежутков времени между локальными экстремумами, а по графикам эволюции меры LAM - резкое снижение амплитуды колебаний. Заключение Целью работы являлось исследование фазового перехода в ДНК, происходящего при термической денатурации, методом количественного анализа рекуррентных диаграмм. Для достижения цели реализована модель Пейрара - Бишопа - Доксуа, которая позволяет определять состояние пар оснований ДНК в произвольный момент времени, проведено моделирование динамики ДНК при различных температурах. Для анализа полученных результатов использован метод количественного анализа рекуррентных диаграмм, благодаря которому появилась возможность судить о характере протекающих в системе процессов. С помощью данного метода анализа произведен расчет мер и представлены их эволюции, исходя из которых был найден момент фазового перехода в системе. а б в Рис. 9. Графики эволюции мер, основанные на РД растяжения второй пары оснований с использованием окна размером S = 20: а - мера времени запаздывания; б - мера максимальной длины вертикальной линии; в - мера замирания По результатам исследования выявлено, что наиболее чувствительными к анализу исследуемого процесса оказались меры, основанные на вертикальных структурах, а именно TT, LAM, Vmax. Данные меры количественно определяют мелкомасштабные структуры рекуррентных диаграмм, изменения в топологии и «восприимчивость» системы к флуктуациям различной природы в окрестности фазового перехода. В отличие от классической теории фазовых переходов, процессы в ДНК демонстрируют качественно новые закономерности, связанные с термодинамикой фазовых переходов, в так называемых одномерных системах. Применительно к термодинамике ДНК это является следствием сочетания в молекуле ДНК свойств, характерных для кристаллической решетки - дальний порядок взаимодействия оснований вдоль цепей ДНК, и «слабых» взаимодействий между цепями (водородные связи), характерные для жидкостей. Это предопределяет низкотемпературный барьер структурных трансформаций в ДНК (локальное плавление), обусловленных термоактивированным разрывом водородных связей, которые также определяются напряженным состоянием (деформированием) ДНК. Энергетический вклад обоих типов взаимодействий отражен в модели Пейрара - Бишопа - Доксуа, которая является базовым представлением при исследовании динамики ДНК и дает уникальную возможность изучения роли коллективных мод дефектов, так называемых открытых комплексов, на термодинамику структурных переходов. Последние определяют фундаментальные закономерности функционирования ДНК (транскрипцию, репликацию) и являются основой при построении термодинамики трансформаций в клетке с учетом роли коллективных степеней свободы, представляющих различные нелинейные моды открытых комплексов. Благодаря возможности выявления фазового перехода, количественный анализ рекуррентных диаграмм является одним из перспективных методов. Данный подход применим не только для исследования денатурации, но и для большинства биологических процессов, например канцерогенез, влияние геомагнитного поля на человека и других.

About the authors

A. S Nikitiuk

Institute of continuous media mechanics UB RAS

O. S Burmistrova

Institute of continuous media mechanics UB RAS

O. B Naimark

Institute of continuous media mechanics UB RAS

References

Statistics

Views

Abstract - 152

PDF (Russian) - 235

Refbacks

  • There are currently no refbacks.

Copyright (c) 2022 Никитюк А.С., Бурмистрова О.С., Наймарк О.Б.