Содержание к диссертации
Введение
1 Литературный обзор 8
1.1 Общие сведения о ДНК 8
1.2 Эксперименты по переносу заряда в ДНК. Основные представления о механизмах переноса 10
1.3 Теоретические подходы к описанию переноса заряда в молекулярных цепочках 15
2 Квантово-механические модели переноса заряда в ДНК 20
2.1 Уотсон-криковские пары как осцилляторы 20
2.2 Обезразмеривание динамических уравнений 24
2.3 Параметры модели и начальные данные 27
2.4 Некоторые свойства модели с детерминированной классической подсистемой 30
2.5 Оценка подвижности заряда при конечной температуре 39
2.6 Различные модели переноса заряда в молекулярных цепочках 40
3 Результаты численного моделирования система без случайной силы (Д = 0) 47
3.1 Моделирование переноса дырки в GTTGGG-фрагменте ДНК 47
3.2 Динамика возбуждения в однородных цепочках 50
3.3 Моделирование переноса заряда в однородных полинуклеотидах с донором и акцептором 54
3.4 Перенос заряда в регулярных и нерегулярных нуклеотидных последовательностях с донором и акцептором 59
3.5 Квазисолитон система со случайной силой (Д ф 0) 64
3.6 Моделирование переноса в GTTGGG-фрагменте при температуре 37 С 64
3.7 Расчеты подвижности дырки в полинуклеотидах при заданной температуре 300К 68
3.8 Расчет подвижности в HSSH-модели ДНК при Т = 300 К 72
3.9 Зависимость подвижности от температуры 73
3.10 Подвижность в случае полярона малого радиуса 76
3.11 Температурный развал стоячей волны 79
3.12 Моделирование динамики переноса заряда в однородном полинуклео-тиде во внешнем электрическом поле 80
4 Вычислительные проблемы 84
4.1 Стандартные схемы 84
4.2 Особенности задачи при = 0 87
4.3 Замечания о численном интегрировании системы 90
4.4 Смешанная схема 93
4.5 Тестирование основного алгоритма 95
Заключение 105
Благодарности 106
Список Литературы
- Эксперименты по переносу заряда в ДНК. Основные представления о механизмах переноса
- Обезразмеривание динамических уравнений
- Моделирование переноса заряда в однородных полинуклеотидах с донором и акцептором
- Замечания о численном интегрировании системы
Введение к работе
В основе множества явлений, происходящих в живых системах, лежат процессы перемещения и транспортировки заряда. Компьютерное моделирование процессов переноса заряда в биологических системах фактически началось в 70-х годах прошлого века и было прежде всего связано с моделированием переноса заряда в белках [1-5]. Начало исследований процессов переноса заряда в ДНК относится к более позднему времени — началу 90-х годов прошлого века, когда в одной из первых экспериментальных работ (Бартон и др. [6]) был продемонстрирован перенос заряда в сравнительно длинном фрагменте ДНК, содержащем 15 пар оснований. С одной стороны это было связано с тем, что стали возможны эксперименты по измерению ультракоротких по времени процессов, обусловленные развитием наносекундной, а затем и фемтосекундной техники, а с другой — с развитием вычислительной техники, позволяющей проводить модельные расчеты таких процессов в сложных биологических системах. Интерес к изучению переноса заряда в ДНК объясняется тем, что перенос является частью механизма таких важных биохимических процессов как репликация, транскрипция, разрушение и репарация ДНК, передвижение радикалов по молекуле ДНК играет существенную роль в процессах мутагенеза и канцерогенеза.
Математическое моделирование процессов переноса заряда в биологических системах связано с использованием дискретных моделей, в которых рассматриваются пути переноса заряда в макромолекулах. Первые дискретные динамические нелинейные модели в биологии были рассмотрены А.С. Давыдовым. В работах [7-9] квантово-классические дискретные модели использовались для описания переноса нелинейных возбуждений в а-спиралях белков. Дальнейшему развитию этого направления посвящены книги и обзоры [10-12]. В представляемой диссертации дискретные динамические модели применяются для описания переноса заряда в поли-нуклеотидных цепочках.
В настоящее время опубликовано множество работ, посвященных моделированию движения заряженной частицы в молекулярных цепочках различного типа (см., например, книги и обзоры [13-15] и ссылки в них). Интерес к этой проблеме связан, с одной стороны, с возросшими экспериментальными возможностями исследования процессов переноса в квазиодномерных молекулярных кристаллах, полимерах и биополимерах, а с другой стороны — с необычными проводящими свойствами таких систем. К новому типу проводящих квазиодномерных молекулярных систем относятся полинуклеотидные цепочки, образующие двойную спираль молекулы ДНК.
Проведенные в последние годы измерения проводимости полинуклеотидных цепочек выявили разброс в проводящих свойствах, простирающийся от изоляторов, полупроводников до проводников и сверхпроводников [16-24].
В связи с этим построение адекватной модели переноса заряда в ДНК является актуальной задачей.
Целью данной работы является исследование модели переноса заряда в ДНК.
Основными задачами, которые были поставлены в ходе исследований, являются разработка теоретической модели переноса заряда в квазиодномерных молекулярных цепочках применительно к ДНК, создание программ и численное исследование динамики переноса заряда в различных пуклеотидных последовательностях, расчет подвижности заряда в полинуклеотидах при различной температуре, а также сопоставление полученных результатов с имеющимися экспериментальными данными.
Научная новизна. В процессе решения поставленных задач разработан новый смешанный алгоритм для численного интегрирования квантово-классическои модели со случайной силой на большие временах счета.
С использованием результатов прямого моделирования и формул Кубо разработан новый способ расчета подвижности заряда в полинуклеотидах.
Впервые теоретически показана возможность переноса заряда в ДНК на большие расстояния.
Впервые рассчитана подвижность дырки для ряда однородных и регулярных последовательностей при комнатной температуре и найдена зависимость подвижности от температуры для однородной синтетической последовательности.
Впервые показано, что при нулевой температуре в однородных полинуклеотидных цепочках, помещенных в постоянное внешнее электрическое поле, возникают блоховские осцилляции.
Практическая значимость. Созданный в ходе выполнения работы комплекс программ позволяет промоделировать перенос заряда вдоль фрагмента ДНК любой заданной последовательности и узнать время переноса, место конечной локализации заряда, подвижность заряда в такой последовательности. В дальнейшем разработанные программы планируется применять для нахождения фрагментов генома, в которых мутация наиболее вероятна, а также при расчетах проводящих свойств "ДНК-нанопроводов" и биочипов на основе ДНК.
Разработанная численная схема может быть использована при исследовании других дискретных моделей квазиодномерных биомакромолекул.
В настоящее время демо-версия (без учета температуры) программы расчета переноса заряда в ДНК доступна на информационно-вычислительном портале "Математическая клетка" (www.mathcell.ru).
Апробация работы. Основные результаты, изложенные в диссертации, докладывались па научных семинарах Института математических проблем биологии РАН (Пущино), на V Международном конгрессе по математическому моделированию (Дубна, 2002г.), на XIV и XV Всероссийских конференциях Теоретические основы и конструирование численных алгоритмов для решения задач математической физики с приложением к многопроцессорным системам" (Дюрсо, 2002 и 2004 гг.), на международной школе-конференции "International School of Crystallography, 38h Course: Structure and Functions of Large Molecular Assemblies" (Erice, Italy, 2006 г.), на I Международной Конференции "Математическая биология и биоинформатика" (Пущино, 2006 г.), на международной конференции EGEE User Forum 01-03 March 2006 (CERN, Switzerland).
Публикации. По результатам диссертации опубликовано 20 работ, в том числе 8 статей в реферируемых научных журналах и две статьи в коллективных монографиях. Перечень статей приведен в конце списка цитируемой литературы.
Структура и содержание диссертации. Диссертация состоит из введения, четырех глав и заключения, содержит 117 страниц, 28 рисунков, 16 таблиц и список цитируемой литературы, включающий 117 наименований.
Во введении дается краткое описание изучаемых в работе задач, обосновывается актуальность темы диссертации, формулируются цели работы, указывается новизна и научная значимость полученных результатов. Описывается структура диссертации и ее краткое содержание по главам.
Первая глава посвящена обзору современного состояния проблемы. Приведены общие сведения о составе и строении дезоксирибонуклеиновой кислоты (ДНК), обсуждаются результаты экспериментов по переносу заряда в ДНК, полученные различными научными группами. Описаны основные дискретные модели квазиодномерных молекулярных цепочек — холстейновская, давыдовская и SSH-модели; приводятся разные способы моделирования температуры системы (термостата).
Во второй главе рассматривается (используемая нами) дискретная квантово-классическая модель переноса заряда в ДНК. Заданная температура моделируется с помощью добавки случайной силы в уравнения движения классических сайтов (уравнения Ланжевена). Приведены значения параметров для исследуемой модели ДНК и описаны варианты задания начальных данных для детерминированной задачи и задачи со случайной силой. Рассмотрены некоторые предельные случаи модели без случайной силы. Выписаны выражения для полной энергии системы, условия локальных минимумов энергии; приведено аналитическое решение "стоячей волны" в однородной цепочке, соответствующее минимуму энергии. Описана схема для оценки подвижности заряда при конечной температуре по формуле Кубо с помощью прямого моделирования. Рассмотрены различные возможные детализации исследуемой модели.
В третьей главе приведены результаты моделирования. Глава разделена на две части; сначала описаны результаты, полученные при численном исследовании детерминированной системы, затем приведены результаты для системы со случайной силой (при конечной температуре). Расчеты проводились для различных фрагментов ДНК; кроме собственно моделирования динамики переноса, проведены расчеты диффузии и подвижности дырки в однородных и регулярных последовательностях и аппроксимирована зависимость подвижности от температуры. Приведены результаты численного моделирования системы при нехарактерных для ДНК параметрах. Найдете численно решение — движущаяся уединенная волна — не является солито-ном в точном смысле, хотя при столкновениях такие волны ведут себя как солитоны. Проведено сравнение зависимости диффузии полученной с помощью прямого моделирования, с аналитическими выражениями Холстейна в случае полярона малого радиуса. Промоделирована динамика стоячей волны при различной температуре, численно найдена температура развала стоячей волны. Рассчитана динамика дырки в однородном полинуклеотиде в постоянном электрическом поле. Показано, что в такой системе движение дырки соответствует блоховским осцилляциям, оценена максимальная температура, при которой такие осцилляции существуют.
В четвертой главе рассмотрен ряд вопросов, связанных с вычислительной задачей. Описаны стандартные методы, которые использовались для численного интегрирования на сравнительно небольших временах счета; показано, что в детерминированной системе время выхода на стационар может быть очень большим. Описана смешанная схема, разработанная для численного интегрирования рассматриваемой системы на больших временах. Приведены тесты смешанного алгоритма для детерминированной системы и для системы со случайной силой.
Эксперименты по переносу заряда в ДНК. Основные представления о механизмах переноса
Основные представления о механизмах переноса Среди макромолекул ДНК занимает особое место. Молекула ДНК похожа на твёрдое тело. Пары оснований уложены в ней как в кристалле. Но это кристалл линейный, как бы одномерный — каждая пара оснований имеет только двух ближайших соседей. Основным отличием кристалла ДНК от обычных является его апериодичность, так как последовательность пар оснований в нём нерегулярна.
Проблема переноса электрона в ДНК давно привлекала внимание исследователей. В работе [27], вышедшей вскоре после открытия Уотсоном и Криком спиральной структуры ДНК, была высказана гипотеза о полупроводниковых свойствах ДНК. Как известно, реакции переноса заряда в ДНК играют важную роль при репликации, транскрипции и репарации молекулы, в процессах мутагенеза и канцерогенеза [28-31]. Однако эксперименты по переносу заряда в ДНК были начаты сравнительно недавно, полтора десятилетия назад. Это связано, с одной стороны, с развитием нано- и фемтосекундной техники, и с другой — с появлением биохимических методов ковалентной привязки молекулярных комплексов, выполняющих роль донора и акцептора, к фрагменту ДНК. В большинстве экспериментов электроны или дырки (катион-радикалы) специальным образом создавались на фрагментах ДНК с известной последовательностью оснований. Скорость электронного транспорта при этом рассчитывалась либо по измерениям тушения флюоресценции, либо сравнением числа разрушений в нуклеотидах, вызванных переносом заряда, в различных участках спирали ДНК. За последние несколько лет была проделана огромная экспериментальная работа по выяснению механизма переноса заряда в ДНК, однако этот вопрос все еще остается открытым.
При описании миграции заряда в молекулярных цепочках обычно рассматриваются три различных механизма переноса заряда.
Первый механизм предполагает некогерентный прыжковый перенос заряда, вызванный температурными флуктуациями [32,33]. В этом случае заряд движется от донора вдоль мостика из сайтов (пуклеотидов либо нуклеотидных пар), совершая прыжки от одного сайта к другому, оставаясь какое-то время на каждом из этих сайтов, продвигаясь, таким образом, от донора к акцептору Таким образом, в случае прыжкового механизма каждый сайт, на который совершается прыжок, является реальным химическим интермедиатом.
Второй механизм — это суперобмен [32,33]. Он соответствует туннелированию заряда от донора к акцептору в один прыжок. В этом случае заряд невозможно обнаружить ни на одном из сайтов мостика, соединяющего донор с акцептором — его присутствие на мостиковых сайтах лишь виртуальное. В предположении, что электрон или дырка передвигаются вдоль фрагмента ДНК со скоростью, экспоненциально падающей с расстоянием (в соответствии с теорией Маркуса скорость реакции к ехр{—/3R} exp{-@NAr}, Аг — расстояние между соседними нуклеотидными парами, N — число пар), результаты экспериментов обычно представлены получен ным значением параметра /?.
И третий механизм — формирование самозахватывающегося полярона, который движется под влиянием температурных флуктуации [34].
Уже первые эксперименты по переносу заряда в ДНК не укладывались в рамки одной модели. Исследования замороженной в воде ДНК показали, что облучение ДНК приводит к миграции зарядов на расстояние менее 8 пар оснований с локализацией на гуанине G [35,36]. Затем было показано существование переноса на расстояние « 20 нм (около 60 пар оснований) [34,37], с эффективностью, зависящей от расстояния переноса.
В экспериментах [38] изучался перенос электрона на расстояние в 10.2,13.6 и 17 А (соответственно 3, 4 и 5 пар оснований между донором и акцептором). В этих экспериментах наблюдалось экспоненциальное убывание скорости переноса электрона с увеличением числа пуклеотидпых пар. Аналогичные результаты были получены в экспериментах [39,40], в которых расстояние переноса составляло 20.5 А(8 нуклео-тидных пар, разделяющих донор и акцептор). Совершенно неожиданными оказались результаты экспериментов [6], в которых изучался перенос электрона на расстояние в 15 нуклеотидных пар, разделяющих донор и акцептор (« 37 А). Несмотря на то, что расстояние переноса в экспериментах [6] почти вдвое больше чем в [38,40], полученная в [6] скорость переноса оказалась на три порядка выше.
В последующих экспериментах было установлено, что в ДНК возможны различные механизмы переноса, а скорость переноса сильно зависит не только от длины нуклеотидной последовательности, вдоль которой происходит перенос, но и от вида последовательности, а также от характеристик доноров и акцепторов, используемых в экспериментах.
Исследования [41] показали, что заряд движется преимущественно вдоль одной нити ДНК. Эксперименты по переносу заряда (катион-радикала), индуцированного светом А = 193 нм, в ДНК в воде [42], продемонстрировали неслучайное распределение однонитевых разрывов и изменений оснований, которое определялось биохимическими методами, в основном на гуанинах в большинстве последовательностей. Если по одной нити ДНК имеется бедная гуанинами область, повреждения неслучайным образом локализуются на аденииах, даже если на комплементарной нити есть один - два гуанина.
Обезразмеривание динамических уравнений
Обозначим характерные величины следующим образом: время т, масса сайта М„, характерный размер смещения сайта Un. Характерный масштаб колебаний Un связан с константой связи а п и может задаваться по-разному. Обозначим безразмерные переменные t = ї/т, ип = un/Un, и перепишем систему (2.14), (2.15)
Получаем следующие соотношения между (почти всеми) размерными и безразмерными параметрами: тап гі/п,„±і 2 т2Кп , Т7„ , 7ч »Ь = Х, 1 = - -, "» = Ж, "п = К, (2-17) и, в соответствии с обозначениями, систему г —г- = VnOn + Vn-i,non-i + Vn+i,nO„+i -\ г—unbn, ( -loj at n Кроме того, случайная сила An(t) должна удовлетворять условиям (2.13). Обозначим т2 Mt)=tnZn{t), (2.20) мпип где
В системе (2.18), (2.19) остались неопределенными выражения a nUnr/h, T2/MnUn и a nT2/MnUn, связь между которыми можно задать по-разному.
На первом этапе мы занимались численным моделированием системы без случайной силы; в дальнейшем оказалось, что для системы со случайной силой (при 0) выбранный нами способ обезразмеривания неудачен, и мы использовали другой способ. Ниже приведены эти способы обезразмеривания. к-способ (при = 0). \ При обезразмеривании потребуем выполнения соотношения a nT2/MnUn = 1, (2.22) т.е. характерная величина колебаний n-го сайта Un = а пт2/Мп зависит от его массы и энергии связи электрона на этом сайте. Система в обезразмеренном виде і - = Vnbn + Vn-i,nbn-i + Vn+i,nbn+i + кпш1ипЬп, (2.23) — = -0, - - -16,,1, (2.24) где через кп обозначена константа связи, (ад) =0, (Zn(t)Zm(s)) = 5nm5(t-s), = 2квТ -. (2.21) .М2 2 т 1тт Т(ап) (2.25)
Обратим внимание на особенность такого выбора безразмерных единиц. В исходных размерных уравнениях движения сайтов (2.15) взаимодействие с зарядом, пропорциональное Ь„2, обращется в ноль при а п = 0. В обезразмерешюй к-способом системе (2.23), (2.24) это взаимодействие не равно нулю ни при каких к для уравнений (2.24). Кроме того, для случайной силы при к-способе обезразмеривания из (2.21) мы получим n = Const\Jijj Jкпш\ (где константа зависит от температуры). Но при этом нельзя провести простейший тестовый расчет среднего квадрата смещения (и2) свободной броуновской частицы из уравнения mil = — fu + A{t) (по формуле Эйнштейна (и2) = 2Dt должна получиться линейная зависимость по времени t), так как для нее ш = О, и температурный коэффициент не определен. Поэтому для задачи со случайной силой более удобным оказался приведенный ниже %-способ обезразме-ривания.
Х-способ. Система со стохастической правой частью ( ф 0). Обозначим характерную температуру Т . Для задачи со случайной силой мы применили следующий способ обезразмери-вания. Вместо (2.22) потребуем выполнения соотношения то есть и Чк х"=а -\1ш„ 2-27)
Обозначим безразмерную температуру через Т = Т/Т ; из формул (2.21) получаем, что действующая на n-ый сайт случайная сила nZn(t), соответствующая этой температуре Т, имеет дисперсию , где „ = 2квТ ф y/bJn. (2.28)
Система в обезразмеренном виде г —Ц- = ЦпК + Vn-i,nbn-i + Vn+i,nbn+i + XnUnb„, (2.29) = -ulun - u/j - Xnlbnl2 + tnZn. (2.30) Мы рассматривали только случай, когда трение ш для всех сайтов одинаковое и, соответственно, п = дая всех гс-Связь между константами связи к и х выражается простой формулой Хп = \Лп #; (2.31) связь безразмерных переменных: $ = $, uMV = uW (2.32)
Отметим, что х-способ более универсален, но при численном интегрировании системы, обезразмеренной к-способом, надо выполнять на одно действие меньше (а именно, умножение: в подсистеме (2.30) член хп&п2, а в (2-24) — просто 6П2)-Связь размерного смещения с безразмерными: для выбранной нами массы сайта у/Нт/Мп и 102.7 А.
Характерное время выбираем соответствующим квантовой подсистеме: г = Ю-14 сек.
Эффективную массу всех сайтов считаем одинаковой, Мп = 10_21гр. Атомная масса нуклеотидной пары равна примерно 600 а.е. « 1.16 10-21гр; эффективная масса Мп = тпітП2/(тпі + тП2) пары оснований с массами тп\ и тп„2 должна быть примерно в четыре раза меньше, но мы дополнительно учитываем молекулы воды, "налипшие" на ДНК [71-73].
Внутримолекулярные колебания в ДНК, отвечающие колебаниям оснований в отдельном сайте, имеют частоты порядка пикосекунд [74,75]. Мы полагаем все частоты и коэффициенты трения на сайтах одинаковыми; соответствующие безразмерные величины лежат в диапазоне шп « 10 2-10 3, ш п « 10-3-10-5 (коэффициенты трения выбираем так, чтобы классическая подсистема не была апериодической).
Для входящих в последовательность оснований прямые измерения потенциалов окисления (энергий сайтов) в настоящее время отсутствуют. Обычно при компьютерном моделировании их принимают равными потенциалам окисления изолированных нуклеотидов в соответствующем полярном растворителе [42,48,76]. Ниже в таблице (Та) приводятся значения потенциалов окисления (ПО) а, полученные с помощью электрохимических измерений [77,78], и соответствующие относительные (за ноль выбирается наименьший ПО) беразмерные значения г]п, которые мы использовали при моделировании переноса заряда в ДНК.
Моделирование переноса заряда в однородных полинуклеотидах с донором и акцептором
В большинстве физических экспериментов по переносу заряда электроны или дырки специальным образом создаются на фрагментах ДНК с известной последовательностью оснований с прикрепленными к ним молекулярными комплексами, например, гуаниновые основания ковалентно связываются с рутений-родиевыми молекулярными комплексами. Комплекс, отдающий заряд, называется донором, а принимающий — акцептором.
Мы рассматривали длинные однородные цепочки с модельными донором и акцептором на концах [А5]. А именно, проводился расчет переноса заряда по цепочкам, состоящим из 150, 200 и 250 сайтов. Первые 10 сайтов моделировали донор заряда, последние 10 — акцептор. Для "мостиковых" сайтов между ними значения параметров т]п = 0 г}п,п±1 = 1-92, и\ = 0.0001, и п = 0.006, кп = 4. Для донора уменьшалась величина потенциала окисления г\в на сайте (остальные параметры — как для сайтов цепочки). Это предположение не противоречит условиям экспериментов но переносу электрона: первоначально заряд находится в ловушке, которая представляет собой потенциальную яму на доноре, и затем переходит на акцептор, параметры которого определяют еще более выгодную с энергетической точки зрения ловушку. Для сайтов акцептора мы выбирали частоты и коэффициенты трения порядка 1, и, как для донора, уменьшали на них потенциал окисления гц. Присвоение сайтам акцептора сравнительно больших значений коэффициента трения ускоряет диссипацию энергии на этих сайтах; в результате возбуждение в процессе перехода углубляет потенциальную яму акцептора, формируя такое распределение смещений, которое понижает энергию возбуждения, локализованного на нем. Расчеты проводились для различных значений параметров донора и акцептора, на времена I порядка ЮОпсек (t « 104). Простая модель При моделировании переноса численно интегрировалась система (3.1), (3.2).
Хотя время необратимого переноса с донора на акцептор различно, качественно картина распространения заряда сходна для различных расчетов. На рисунке 13 показаны графики распределения вероятностей ЬП()2 и суммарных энергий дырок Еп = г}п+Kn nun(t) п0 сайтам цепочки длиной 200 сайтов в одинаковые моменты времени. В нулевой момент времени заряд локализован на доноре (обычно выбиралось h{t = 0) = 1), и сайты находятся в равновесном положении (ип = йп — 0). После нескольких отражений возбуждение практически равномерно размазывается по всей цепочке (рис. 13 в) и остается в таком состоянии сравнительно большое время. Затем вероятность нахождения заряда на акцепторе начинает увеличиваться.
Графики рассчитаны для следующих значениях параметров: rji = ... = Що = %oi = 200 = — 0-4 (потенциал окисления на доноре и акцепторе для наглядности выбран одинаковым), ш\ = 2, и А = 2 (A = 191,...,200), и% = 0.0001, и п = 0.006 (n= 1,...190), кк = 4 (к = 1,...200).
За точку отсчета энергии принято её значение в момент t = 0 на сайтах мостика (п — 11,...190). На доноре и акцепторе в начальный момент энергия отрицательна (на графиках а рис. 13 это соответствует наличию потенциальных ям на концах
цепочки). Поскольку величина шп на донорных и мостиковых сайтах на 2 порядка меньше, чем на акцепторных, вклад смещений в изменение энергии на этих сайтах почти не заметен. Захват дырки на акцепторных сайтах приводит к изменениям в энергии на них. На стадии, когда дырка размазана по всем сайтам цепочки (графики в), профиль энергии почти не отличается от начального. На конечной стадии захвата дырки акцептором происходит углубление потенциальной ямы, соответствующее релаксированному состоянию дырки.
По результатам вычислительных экспериментов можно предположить следующее.
Перенос происходит, если значения щ энергии электрона на доноре и акцепторе близки к значениям на "мостиковых" сайтах. При уменьшении значений гц, на доноре (остальные параметры фиксированы) происходит "частичный перенос", т.е. через достаточно большой промежуток времени вероятность распределяется между донором и акцептором. Например, для параметров г]р = — 3 (D = 1,.. .10), TJA = —0.2, ш\ = 2, и/А = 2 (А = 191,.. .200) значения вероятности на доноре и акцепторе близки к 0.5. При дальнейшем уменьшении энергии электрона на доноре (TJD 10) возбуждение все время расчета локализовано на доноре и нескольких первых мостиковых сайтах.
При прочих фиксированных параметрах существует оптимальная "глубина" ямки энергии электрона на сайтах акцептора г]д, для которой время переноса будет наименьшим. При уменьшении значений г\А время, за которое вероятность на акцепторе будет равна, к примеру, 0.9, увеличивается, затем наблюдается "частичный перенос", и при дальнейшем уменьшении ( —10) возбуждение не захватывается акцептором. При Г\А = 0 возбуждение также не локализуется на акцепторе.
Время переноса уменьшается при увеличении частот сил и коэффициентов трения w A на сайтах акцептора, пока значения частот достаточно близки к значениям матричных элементов 77„іП±і. При значительном изменении частот по отношению к J7n,n±i время увеличивается. В таблицах Т4, Т5 (стр. 57) продемонстрирована зависимость времени переноса от значений этих параметров акцептора (значения остальных параметров приведенных в таблицах расчетов соответствуют параметрам для рис. 13).
Время переноса зависит от начальных условий — оно наименьшее, если возбуждение локализовано на первом сайте, и максимально при равномерном распределении вероятностей по всем сайтам донора.
Замечания о численном интегрировании системы
В рассматриваемой системе ОДУ {4.8), (4.9) существует первый интеграл — сохраняющаяся величина Е = ]n bn()2 = 1. При численном интегрировании с помощью указанных выше алгоритмов этот интеграл сохраняется не точно. Контролем точности численного решения может служить разность є = 1 — Е, которая со временем растет. Следующий прием, позволяющий увеличить времена счета — введение искусственной нормировки, при которой переменные Ьп "подправляются" так, чтобы сумма квадратов их модулей снова равнялась единице. При этом мы, возможно, не возвращаемся на ту же траекторию, а "перескакиваем" на новую (соответствующую траектории точного решения из других начальных данных), но в полной задаче, когда нас интересует среднее по множеству реализаций, индивидуальные решения не важны. В нашем случае искусственная нормировка — это проектирование всех Ьп на единичную сферу в С, т.е. все Ьп делятся на число vE. При этом получается новый метод интегрирования.
Действительно, пусть в момент to мы знаем точное решение В = B(t0). Сделав один расчетный шаг h, мы уйдем от точного решения В = B(to + h) в точку В (либо наружу, либо внутрь сферы). После нормирования мы окажемся в лежащей на сфере точке В. Получаем треугольник, в котором (см. рис. 26). ошибка исходного метода 5orig = 5В , ошибка нового метода 8new = \ВВ\. Если В лежит вне сферы, то треугольник ВВ В тупоугольный, и Snew S ig- Если В лежит внутри сферы, то 5new и 5ОГід(1+5Ігід/2). В любом случае, порядок метода от добавления процедуры нормирования не изменился.
Идея искусственного нормирования далеко не нова. В ряде задач молекулярной динамики, где со временем должна сохраняться полная энергия системы Е = EQ, давно применяются подобные алгоритмы "возвращения" фазовой точки с поверхности Е = EQ + 5Е по нормали на поверхность Е = Еа [87,117]. В этом параграфе подробно описан алгоритм, примененный нами для серийных расчетов. Напомним еще раз про особенности задачи. Нужно численно решить систему
В ней есть две подсистемы осцилляторов — быстрые Ьп (4.20), с характерным временем порядка единицы, и медленные ип (4.21), с характерным временем порядка сотни или тысячи. Подсистемы связаны между собой нелинейно. В быстрой подсистеме есть сохраняющаяся величина (первый интеграл) S = 1. Нас интересует поведение решений на больших временах — миллион и дальше. В медленной подсистеме есть член со случайной величиной, поэтому точно отслеживать индивидуальные траектории не нужно, но при этом Е не должна сильно отклоняться от единицы, иначе решение теряет физический смысл. Поскольку нам затем понадобятся осреднснпые по реализациям траектории, то расчет каждого отдельного решения должен требовать разумного машинного времени.
Два базовых алгоритма, используемых нами для решения полной задачи:
(1) Стандартный метод Рунге-Кутта 4-го порядка для системы (4.20), (4.21) при () = 0. Он позволяет считать с хорошей точностью на большие времена (за счет мельчения шага и значительных затрат машинного времени). Этот метод неприменим для решения систем ОДУ со случайной правой частью.
(2) Метод (4.14) для решения систем ОДУ со случайной силой в правой части (описанный в [68]), который при = 0 совпадает со схемой Рунге-Кутта 2-го порядка; он применим на временах порядка тысячи колебаний переменных Ьп (в дальнейшем нарушение условия нормировки становится существенным).
Идея смешанного алгоритма состоит в следующем. Быструю подсистему (4.20) решаем с помощью более точного алгоритма (1) с мелким шагом т, используя при этом приближенные значения . На некоторые вопросы, связанные с использованием смешанного алгоритма, проще всего оказалось ответить, проведя ряд численных экспериментов. В этом параграфе приведены выборочные тесты, позволяющие оценить величину "производственного" шага т, отношение шагов интегрирования К — h/т, выбор "отклонения нормировки"
Enorm — 1-І 2-і г
Сначала все вопросы рассматривались для системы (4.20), (4.21) в отсутствие случайной силы (при = 0); затем для величин т и К, которые мы сочли оптимальными для серийных расчетов, были проведены тесты при разных значениях . В качестве иллюстраций приведены выборочные результаты тестов.
В тестовой системе мы находим точное решение (проведя расчет методом Рунге-Кутта 4-го порядка с достаточно мелким шагом). Все посчитанные с помощью смешанной схемы результаты сравниваются с точным решением. В таблицах этого параграфа приводятся разности А между точным и расчетным значениями переменных в некоторый момент, и/или разность квадратов модулей Ьп (6„2 для нас более важен, чем фаза Ьп). Величины коэффициентов тестовых систем следующие: ип = 0.01, ш п = 0.006, х = 0.02; значения матричных элементов rjnk соответствуют параметрам ДНК (см. таблицы TQ, Т„ на стр. 27).
Порядок смешанной схемы.
Формально смешанная схема при одинаковых шагах т, h имеет второй порядок точности. По результатам тестов получается, что этот метод имеет четвертый порядок точности по Ъп и второй порядок по ип, vn.
Приведем пример — тест для GAAGG цепочки. При t = 0 xi = 1 (6n = хп + гуп), все остальные начальные данные равны нулю. В таблицах ниже приведены разности Д для переменных первого и третьего сайтов в момент t = 10.
