Численное моделирование формирования и эволюции солитонов в нуклеотидных последовательностях Текст научной статьи по специальности «Физика»

УДК 519.688+519.62+538.9 УДК 519.6

Г. В. Квитко, Н. В. Малышева, С. А. Ишанов, Д. В. Шоть, Е. М. Логвин

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ФОРМИРОВАНИЯ И ЭВОЛЮЦИИ СОЛИТОНОВ В НУКЛЕОТИДНЫХ ПОСЛЕДОВАТЕЛЬНОСТЯХ

В рамках модифицированной модели Холстейна, учитывающей температурные флуктуации среды, численно решается задача Коши о движении избыточного заряда вдоль однородной дискретной цепочки нуклеотидных пар оснований ДНК. Для специально подобранных начальных условий исследуется эволюция особого класса решений этой задачи,

© Квитко Г. В., Малышева Н. В., Ишанов С. А., Шоть Д. В., Логвин Е. М., 2014

Вестник Балтийского федерального университета им. И. Канта. 2014. Вып. 10. С. 20 — 30.

относящихся к типу солитонов. Показано, что существование решений в виде уединенных волн в рамках принятой модели возможно только при сверхнизких температурах термостата.

The Cauchy's problem of the induced charge movement along the discrete chain of nucleotide couples in DNA is solved in the frame of Holstein's model taking account surrounding temperature fluctuations. For the initial conditions chosen specifically the evolution of a special class of solutions of this problem related to the type of soliton is investigated. The existence of solutions in the form of solitary waves in the framework of this model is possible only at very low temperatures the thermostat.

Ключевые слова: солитон, нуклеотидные последовательности, заряд, амплитуда вероятности, задача Коши, температура.

Key words: soliton, nukleotidny sequences, charge, probability amplitude, Cauchy's task, temperature.

Введение

21

В настоящее время в связи с развитием нанобиоэлектроники особое внимание привлекли исследования уникальных проводящих свойств молекул ДНК [1]. Экспериментально было показано [2], что перенос заряда в ДНК возможен на значительные расстояния. Все это сделало актуальной задачу изучения подвижности носителей тока в ДНК. Определенный интерес вызывают работы по исследованию солитонной динамики в одномерных системах [3 — 7]. Процесс формирования солито-на как пример квантовой самоорганизующейся системы уже давно представляет собой самостоятельный научный интерес. В этом смысле следует отметить работу [8] по изучению автолокализованных состояний в деформируемых цепочках, где исследовались стационарные состояния электрона, эквивалентные стоячим уединенным волнам или движущимся солитонам. Солитоны и поляроны в современной науке рассматриваются в качестве главных претендентов на роль переносчиков зарядов в ДНК. Ввиду этого стоит подчеркнуть, что время образования состояний типа уединенной стационарной волны является важной характеристикой, определяющей физические свойства проводящих молекулярных цепочек [6].

1. Модельные уравнения движения

Для решения задачи воспользуемся моделью Холстейна [9], которая была применена для описания движения поляронов в кристаллах. Молекула ДНК во многом похожа на твердое тело. Комплементарные нуклеотидные пары оснований, объединяемые водородными связями, будем их далее называть сайтами [13], уложены в ней как в кристалле. Но ДНК это кристалл линейный и одномерный. Каждая пара сайтов имеет только двух ближайших соседей. Основное отличие кристалла ДНК от обычных — его апериодичность, поскольку последовательность сайтов в нем нерегулярна.

22

В модели введен ряд допущений. В частности, считается, что смещения сайтов из равновесных положений малы. В связи с этим, во-первых, колеблющиеся сайты рассматриваются как гармонические осцилляторы, а во-вторых, предполагается, что энергия электрона en(t) на n-м сайте линейно зависит от амплитуд смещения ük (t) сайтов:

е. (t) = £0 (t).

m

Здесь константа еП — собственная энергия избыточного заряда, внесенного в ДНК и локализованного на n-м сайте (потенциал окисления); a'nm — константа связи между зарядом и сайтами; ün (t) — смещение n-го сайта из положения равновесия в момент времени t . Значками «~» снабжены исходные реальные параметры задачи до проведения последующей процедуры обезразмеривания.

Полный гамильтониан системы можно записать в виде

H = H quant + H class + H int . (1)

Здесь H quam — гамильтониан квантовой подсистемы, описывающий движение заряда в некотором эффективном одночастичном потенциале; Hclass — гамильтониан классической подсистемы, представляющий собой сумму кинетической (T), потенциальной (U) и диссипативной (D) составляющих энергии сайтов, входящих в рассматриваемую цепочку:

H class = T + U + D . Гамильтониан Hint — это оператор взаимодействия квантовой и классической подсистем, который зависит как от положения заряда, так и смещения сайтов из невозмущенных равновесных положений. Выпишем в соответствии с [13] гамильтонианы в явном виде.

H . =Ve0fl+а +У V (я+а + а+а ).

quant / , nnn / , n ,m\ nm m ni

Здесь операторы а„ и а„ — соответственно операторы рождения и уничтожения заряда на и-м сайте (п = 1, 2,..., N), где N — число сайтов в цепочке); упт — матричные элементы оператора перехода между сайтами п и т.

Составляющие гамильтониана НС1ЯББ, описывающего движение сайтов: Т = 2 I Мп%, и = 1 X К*2, V = 1X ^птйпйт .

2 п 2 п 2 п,т

Здесь Мп — эффективная масса; кп — константа упругости п-го сайта; йпт — коэффициенты дисперсии между и-м и т-м сайтами. Если все коэффициенты йпт = 0, то модель сводится к одночастотной (то есть имеем эйнштейновский осциллятор). В работе [4] оператор энергии смещения сайтов из равновесных положений вдоль направления водородных связей в гармоническом приближении имел вид О = X йп (йп - йп_1)2 /2.

п

Оператор взаимодействия заряда с сайтом: = ХаП,тйт(^ )яИ ®т.

Используя обозначения Дирака, полагаем, что заряд, локализованный на n-м сайте, находится в состоянии | n) — собственного вектора оператора

Hq0 = X^anân , и H0|и) = е°п\и). Поскольку векторы \n) образуют базис в

n

гильбертовом пространстве (ортонормированы и удовлетворяют условию полноты), то волновую функцию всей рассматриваемой системы

представим |Y = Xbn(t)|и), где bn(t) — амплитуда вероятности нахожде-

п

ния заряда на n-м сайте (и = 1, ..., N , N — длина цепочки), для которой в любой момент времени выполняется условие нормировки X |bn (t )|2 = 1.

n

Заметим, что в некоторых работах, например в [6], вместо представления гамильтониана Hquant, описывающего квантовую подсистему путем использования операторов рождения а+ и уничтожения ап, может быть применена эквивалентная и более удобная форма представления гамильтониана через операторы проектирования |и)(и| и операторы

перехода |и)(ш| состояний заряда. В этом случае имеем

^^quant =Х8ПИ(и1+Xvn,Jn)H.

n nïm

Для описания переноса возбуждения будем искать решение линейного уравнения Шредингера

н|*>. (2)

Следующее важное допущение используемой модели — то, что движение сайтов описывается в так называемом приближении ближайших соседей. Математически это означает, что vnm = 0, если т ф и, и +1. Заметим, что vnn =en (t ) — это энергия электрона на n-м сайте. Несложно показать, что осредненный по состоянию оператор (1) системы заряд—сайты с учетом приближения ближайших соседей будет

H\Y) = X VnmbnbH + 2 T(MnU + KnU )+ 1 S dnmUnUm .

n,m=n-1,n,n+1 n n ,m

Так как большинство экспериментов по переносу заряда в ДНК проводится для фрагментов ДНК в (водном) растворе при определенной температуре T, то для учета температурных эффектов в представленном выше варианте модели Холстейна в классическую часть гамильтониана необходимо ввести составляющую типа HStoch(t, T), отвечающую за воздействие температурных флуктуаций термостата на движение сайтов. Следуя [7], воспользуемся для этого методом броуновской динамики, который широко применяется при моделировании движения макромолекул в растворе при заданной температуре и, в частности, для одномерных моделей [10], и представляет собой обобщение метода Ланжевена для описания движения броуновской частицы в растворе.

23

24

После применения процедуры усреднения к уравнению Шредин-гера (2) с последующим обезразмериванием получается следующая система обыкновенных дифференциальных уравнений (ОДУ) для амплитуд вероятностей Ьп (/) обнаружения заряда на и-м сайте и амплитуд смещений ип () сайтов в приближении ближайших соседей [7]:

+ ХпйпЬп, (3)

йп =-(®П + 2С ) йп - Сп (йп-1 + йп+1 ) - Хп \К \2 -<йп +^п (*) . (4)

В этой самосогласованной системе, размерность которой зависит от длины цепочки, учтены процессы взаимного влияния движущегося по сайтам заряда и совершающих колебательное движение сайтов. Здесь / — мнимая единица. Точки над амплитудами — производные по нормированному времени t = I / т, где т — характерный временной масштаб процесса перехода т = 10-14 с. Безразмерная амплитуда смещения сайтов равна йп = йп /ип, где ип — характерный масштаб колебаний. Нормированные элементы матрицы переходов цпт =Уптт / Й, причем

"Лп,т =^т,п, а "Лп =Цп/п =е0т / Й, где Й — постоянная Планка. Константа связи квантовой и классической частей системы хп имеет вид: Хп = а'пт]т3 /(ЙМп), где Мп — эффективная масса п-го сайта. Частота колебаний сайтов юп связана с коэффициентом упругости кп формулой юп = Т\]кп / Мп. Нормированные коэффициенты дисперсии Сп равны Сп = т2йп / М, где йп — реальные коэффициент дисперсии. Параметр а'п =туп /Мп — коэффициент трения или диссипации, задаваемый в диапазоне ю'п и 105 — 103 так, чтобы классическая система не была апериодической. Коэффициенты '%п связаны с температурой среды Т и коэффициентом трения уп формулой = 2квТ*(т / Й)(Т /Т*)уп, где кв — постоянная Больцмана; Т* — характерная температура среды (Т* = 300 К). Отметим, что добавление случайной силы в уравнения движения превращает все динамические переменные в случайные величины.

В приведенных выше уравнениях (3) — (4) мы рассматриваем более упрощенный вариант описанной выше модели, которая, описывает задачу для неоднородной цепочки сайтов. Именно в силу однородности рассматриваемой в данной работе цепочки для каждого п-го сайта большинство параметров задачи — константы. В частности, Мп = М,

Уп = У, £п = ^ Юп =a, Хп = Х, Цп,п+1 = Чп,п-1 . Кроме Tого, исключили

из рассмотрения дисперсионные эффекты, определяемыми составляющей гамильтониана О, иначе говоря, все коэффициенты Сп в (3) равны нулю. Далее, при обсуждении результатов численных расчетов все параметры задачи будут определены количественно.

В уравнении (4) смещение сайтов относительно положения равновесия фактически описывается уравнениями Ланжевена, в которых учтены

процессы диссипации (затухания) и температурные флуктуации среды (термостата) в виде случайной силы со специфическими статистическими свойствами: (¥п ^= 0, (¥п ^)¥т (í= 5пт5 (í -1'). Угловыми скобками

обозначено усреднение по времени. Это дельта-коррелированая случайная величина: ее автокорреляционная функция равна дельта-функции. Такой случайный процесс часто называют белым шумом. В последнем уравнении 5пт и 5(£ - £') — соответственно символ Кронекера и дельта-функция Дирака. При условии ¥п (£) = 0 рассматриваемая модель имитирует

перенос заряда вдоль оснований ДНК без учета температуры окружающей среды. Система (3) — (4), дополненная начальными условиями для амплитуд ип (0), йп (0) и Ьп (0), определяет собой задачу Коши.

2. Результаты расчетов

Нас интересуют, только те решения Ьп (Ь) задачи Коши и их эволюция при разных температурах среды, которые имеют вид уединенной волны — солитона. Численное исследование поставленной задачи по формированию и последующему отслеживанию временной динамики солитонов проводилось на основе двух вариантов модели Холстейна: классического, где по определению никак не учтено температурное влияние среды, и модифицированного, учитывающего температуру, путем введения в уравнения Ланжевена (4) случайной силы ¥п (). Формально, можно полагать, что в классической модели ¥п (/) = 0 .

На первом этапе все исследования велись на базе классического варианта модели Холстейна (без учета температуры среды). Численное решение задачи Коши для системы дифференциальных уравнений проводилось по схеме Рунге — Кутты 4-го порядка точности. Точность решений корректировалась требованием нормированности амплитуд

вероятности: |Ьп(0|2 - 1|<е. Поскольку Ьп (/) — это комплексная фун-

п

кция Ьп ^) = хп ^) + 1уп (í), то вместо одного уравнения (3) решались два: для вещественной хп () и мнимой уп (t) компонент.

В [5] было показано, что формирование решений в виде солитонов прежде всего зависит от соответствующего выбора начальных условий Ьп (0), тогда как варьирование ип (0) и ип (0) при фиксированных Ьп (0) практически на результаты вычислений не влияет. Поэтому сначала были проведены расчеты по определению времени tsol выхода солито-на на стационарный установившийся режим в зависимости от выбора начальных условий задачи Коши для плотности вероятности Ьп ^). При

этом все прочие параметры задачи считались заданными. В частности, рассматривалась однородная цепочка типа гуанин—цитозин (С/С) из N = 51 сайтов при фиксированных параметрах системы: к = 2, ^ = 0,8, где ^ = XV / Й, 0, ю2 = ю' = 1. Во всех расчетах начальные амплитуды

25

смещений ип (0) выражались через начальные амплитуды вероятностей Ьп (0) по закону йп (0) = - |Ьп (0)^, а начальные скорости ип (0) = 0 .

В соответствии с [6] рассмотрены два варианта задания начального распределения Ьп (0), удовлетворяющего требованию нормировки: А —

равномерное распределение типа |Ьп (0) = N 1/2 и В — делокализованное:

\Ьп (0)| = (л/2/4)^/кя/\ц\ • сН 1 (кя(п-п0)/4, где а = 0,02 и п0 = N/2. Параметр а «растягивает» функцию обратного гиперболического косинуса на всю длину цепочки сайтов. Так как Ьп (0) = хп (0) + у (0), то определены начальные условия для хп (0) и уп (0). Применены два разных типа условий:

I — негладкие: *я(0) = |Ь„(0)|(-1)" />/2, уп(0) = \Ьп(0)|(-1)"+7ч/2, п = N/2и

II — гладкие: хп (0) = |Ьп (0), уп (0) = 0, п0 = N/2.

На рисунках 1—4 представлены графики решений задачи Коши

без учета температуры термостата для функции |Ьп (/)) в зависимости

от номера п сайта в однородной цепочке (С/С). Они получены для четырех разных конфигураций начальных условий А1, А11, В1, В11. Из рисунков видно, что если начальные условия гладкие, то заряд сначала размазан по цепочке сайтов и только потом уже собирается в солитон.

Рис. 3. График |Ьп (0|2 . Распределение Рис. 4. График |Ьп(1)| . Распределение В1 (делокализованное негладкое) В11 (нормальное гладкое)

Времена формирования солитонов для конфигураций с гладкими распределениями следующие: ЬА11 и 700 (рис. 2) и Ьт и 900 (рис. 4). Для негладких нормальных условий этот процесс происходит быстрее: Ьм и 250 (рис. 1) и Ьш и 50 (рис. 3). Напомним, что реальное время это

I = £ • 10-14 с. Самым лучшим вариантом по скорости формирования со-литона оказался В1 с негладким делокализованным начальным условием, построенным с использованием обратного гиперболического косинуса. Результаты расчетов оказались близки с полученными в [6].

На втором этапе были проведены численные расчеты задачи Коши (без учета температуры) по выявлению зависимости времени образования солитона от различных значений параметра ю ' для наилучшего варианта В1. Связь ю ' с коэффициентом трения у = ю' = ху/М. На рисун-

ках 5 — 6 представлены графики решений для двух вариантов: ю ' = 0,25 и ю ' = 5.

полученные при ю2 = 1

27

Рис. 5. График |Ьп (£)|2. Распределение Рис. 6. График |Ьп(£)|2. Распределение В1 (ю ' = 0,25) В1 (ю ' = 5,0)

2

Из сравнения решений, полученных для разных ю', следует, что ^ (ю' = 0,25) и 10 (рис. 5) < ^ (ю' = 1) и 50 (рис. 3) < ^(ю' = 5) и 150 (рис. 6). Поэтому можно утверждать, что с ростом параметра ю' или коэффициента трения у время стабилизации солитона увеличивается.

Но при учете температуры среды Т численное решение системы (3) —(4) , содержащей случайную силу Рп ^) на основе метода Рунге — Кутты, становится неприемлемым. Поэтому в расчетах применяют другие способы, например, 202ДС -метод (двухшаговый метод второго порядка с одной гауссовой величиной), предложенный Гринзайдом и Гельфандом в работе [11], или его различные модификации типа смешанного алгоритма Фиалко [12]. Следуя последней технологии, были проведены расчеты задачи Коши (3) — (4) с делокализованным начальным условием для плотности вероятности Ьп (0), использующим обратный гиперболический косинус и соответствующим стоячей волне при фиксированных параметрах, совпадающих с заданными в [12]. Начальные условия для смещений сайтов формировались в виде суммы ип (0) = и1п + и2п, где и1п, а также начальная скорость ип (0) определялись

28

из условий равновесного состояния для гармонического осциллятора при фиксированной температуре, а и2п = -х\Ъ„ (0)| /ю2„ — смещение, соответствующее стоячей волне. Расчеты проводились для (С/С) поли-нуклеотида длиной N = 51 сайтов с параметрами хп = 0,02, юп = 0,01, = 0,006, г]пт = 1,276, = 0 на времена порядка г = 5 -104 (реальное

время г = 5 -10-10 с).

Главный вывод, к которому привели численные расчеты в рамках взятой модели, что тепловые флуктуации в термостате с ростом температуры быстро приводят к разрушению начального состояния, заданного в виде уединенной стационарной волны — солитона.

Ниже приведены графики распределения вероятностей нахождения заряда на сайтах со временем Ьп (г) для трех вариантов сверхнизких температур термостата: Т = 0,1 К, Т = 15 К и Т = 25 К. Они иллюстрируют динамку эволюции решений задачи Коши для системы ОДУ (3) — (4), описывающей движение избыточного заряда по сайтам цепочки. Расчеты показывают, что при температуре Т = 0,1 К (рис. 7) сформированный солитон практически не меняет свою форму на протяжении всего времени счета (\ = 0,5 нс). При Т = 15 К видно (рис. 8), что волна со временем меняет свой вид, ее классическая форма искажается.

Рис. 7. График

Рис. 8. График )| . Т = 15 К

При Т = 25 К (рис. 9) со-литон быстро распадается и к концу времени счета полностью расплывается по всей цепочке сайтов. Такая картина с разрушением солитона естественно наблюдается и при всех температурах выше Т = 25 К. Просто с дальнейшим ростом температуры время существования солитона быстро стремится к нулю. Это говорит о том, что при реальных нормальных температурах солитон в однородной цепочке практически мгновенно разрушается или же попросту не образуется.

Разумеется, такой вывод мы получаем, решая поставленную задачу в рамках выбранной нами математической модели Холстейна. Как было замечено в работе [12], «в используемой простой модели не учиты-

Рис. 9. График \Ьп(*)|2. Т = 25 К

ваются сольватационные эффекты. Для ДНК в водном окружении энергия связи заряда с сайтом сильно возрастает и можно ожидать, что при биологически значимых температурах Т = 300 К солитонные состояния могут оказаться устойчивыми». В работе [3] авторы предлагают отказаться от линейных квантово-механических уравнений системы в целом и перейти к нелинейной квантовой механике подсистем.

Список литературы

1. Lakhno V. D. DNA nanobioelectronics // International Journal of Quantum Chemistry. 2008. № 108. P. 1970-1981.

2. Porath D., Bezryadin A., de Vries S. et al. Direct measurement of electrical transport through DNA molecules // Nature. 2000. Vol. 403. P. 635-638.

3. Коршунова А. Н., Лахно В. Д. Моделирование образования самозахваченного состояния в полинуклеотидной цепочке // Нелинейная динамика. 2008. Т. 4, № 2. С. 193-214.

4. Лахно В. Д., Фиалко Н. С. Механизм остановки движущегося солитона в однородной молекулярной цепочке // Компьютерные исследования и моделирование. 2009. Т. 1, № 1. С. 93 - 99.

5. Квитко Г. В., Малышева Н. В. Исследование влияния начальных условий на время образования уединенной волны в нуклеотидных последовательностях ДНК // Высокопроизводительные вычисления — математические модели и алгоритмы : материалы II Международной конференции, посвященной Карлу Якоби. Калининград, 3 — 5 октября 2013. C. 135-137.

6. Коршунова А. Н., Лахно В. Д. Моделирование образования солитона в однородной цепочке // Математическое моделирование. 2007. Т. 19, № 1. С. 3 — 13.

7. Тихонов Д. А., Соболев Е. В., Лахно В. Д. и др. Адиабатическое приближение при расчетах подвижности заряда в холстейновской модели ДНК // Математическая биология и биоинформатика. 2011. Т. 6, № 2. С. 264—272.

8. Давыдов А. С. Солитоны в квазиодномерных молекулярных структурах // УФН. 1983. Т. 138, № 12. С. 605—643.

9. Holstein T. Studies of polaron motion: Part I. The molecular-crystal model // Annals of Physics. 1959. Vol. 8. P. 325-342.

10. Helfand E. Brownian dynamics study of transitions in a polymer chain of bistable oscillators // J. Chem. Phys. 1978. Vol. 69. P. 1010-1018.

11. Greenside H. S., Helfand E. Numerical integration of stochastic differential euations-II // The Bell System Tech. J. 1981. Vol. 60. P. 1927-1940.

12. Фиалко Н. С. Смешанный алгоритм расчета динамики переноса заряда в ДНК на большие временные интервалы // Компьютерные исследования и моделирование. 2010. Т. 2, № 1. С. 63 — 72.

13. Фиалко Н. С., Лахно В. Д, Динамика переноса заряда в упорядоченных и неупорядоченных нуклеотидных последовательностях // Математическая биология и биоинформатика. 2006. Т. 1, № 1. С. 58—65.

Об авторах

Геннадий Васильевич Квитко — канд. физ.-мат. наук, доц., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: [email protected]

29

Надежда Васильевна Малышева — асп., Балтийский федеральный университет им. И. Канта, Калининград.

E-mail: [email protected]

Г. В. Квитко, Н. В. Малышева, С. А. Ишанов, Д. В. Шоть, Е.М. Логвин

30

Сергей Александрович Ишанов — д-р физ.-мат. наук, проф., Балтийский федеральный университет им. И. Канта, Калининград. E-mail: [email protected]

Дмитрий Владимирович Шоть — ассист., Балтийский федеральный университет им. И. Канта, Калининград. E-mail: [email protected]

Евгений Михайлович Логвин — асп., Балтийский федеральный университет им. И. Канта, Калининград. E-mail: [email protected]

About the authors

Dr Gennady Kvitko — Ass. Prof., I. Kant Federal University, Kaliningrad. E-mail: [email protected]

Nadegda Malysheva — PhD student, I. Kant Federal University, Kaliningrad. E-mail: [email protected]

Prof. Sergey Ishanov — I. Kant Federal University, Kaliningrad. E-mail: [email protected]

Dmitry Shott — Ass., I. Kant Federal University, Kaliningrad. E-mail: [email protected]

Evgeny Logvin — PhD student, I. Kant Federal University, Kaliningrad. E-mail: [email protected]