Статья: К вопросу о численном расчете пространственно-распределенных динамических систем с запаздыванием по времени

Таблица 1. Список реакций транскрипции генов frq и wc. Здесь, и далее k, k1, k-1, k2, k-2, гF , гW, bF, bW - скорости соответствующих реакций.

Известно, что белки FRQ и WCC могут существовать как в форме мономеров, концентрацию которых обозначим F и W, так и в форме димеров 2F и 2W. Обе формы активно взаимодействуют друг с другом - в системе одновременно протекают реакции димеризации и распада димеров соответствующих белков (реакции 1 и 2 в табл.1). Примем во внимание, что во время транскрипции любого гена происходит регуляция этого процесса посредством открытия и закрытия оператор-сайта молекулой димера (отмечен на рис. 4 черным цветом). Наличие такого оператор-сайта можно учесть введением специальной функции , которая принимает значение в случае открытия оператор-сайта и в случае его закрытия. С формальной точки зрения состояние оператор-сайта можно рассматривать как дополнительный реагент в системе. Своеобразие данной системы состоит в том, что гены frq и wc-1,2 регулируют синтез белков друг друга, формируя положительную связь оператор-сайт гена считается открытым, если к нему прикрепляется димер белка-партнёра (реакции 3 в табл.1).

Процессы синтеза белков задаются уравнениями 4 в табл. 1. Они описывают считывание генетической информации в момент времени , если оператор-сайт открыт, РНК-полимераза прикрепляется к началу гена и процесс транскрипции стартует. Конечным продуктом этой реакции являются мономеры белка, которые синтезируются в момент времени . Время запаздывания необходимо системе, чтобы полностью считать информацию от начала до конца гена, а затем синтезировать белок. Хотя времена запаздывания и для синтеза белков FRQ и WCC, следует заметить, разные, мы для простоты полагаем .

Белки могут выводиться из системы с помощью двух механизмов. Один из них предполагает обычную деградацию со временем (реакции 5 в табл. 1). Другой механизм учитывает нелинейное взаимодействие FRQ и WCC, в результате которого формируется общий гетеродимерный комплекс FRQ/WCC (реакция 6 в табл. 1). Мы предполагаем, что комплекс FRQ/WCC сам по себе не участвует в генной регуляции циркадианных ритмов, хотя его образование фактически выводит молекулы белков из системы и этим влияет на общую динамику (см. рис. 4). Эти процессы играют важную роль в системе, формируя петли положительной и отрицательной обратной связи.

Кинетические уравнения для динамики системы, записанные на основании списка реакций (табл.1), после небольших преобразований [12] принимают вид

где , , - оператор Лапласа в двух пространственных измерениях и . В численных расчетах были использованы значения параметров, приведенные в табл.2.

Таблица 2. Параметры модели

Система уравнений (8-9) интегрировалась в квадратной области (, ), на границе которой задавались следующие условия

, ,

, .

Начальным условием являлось случайное распределение концентраций реагентов. В промежутке времени [-6;0] поля концентрации задавались гармоническими функциями со случайной амплитудой, фазой и частотой.

Краевая задача (8-10) решалась с помощью метода конечных разностей. Уравнения и граничные условия аппроксимировались на равномерной сетке 200 на 200 узлов с помощью аппроксимаций второго порядка по пространственным координатам. Использовалась явная схема (3), и для обеспечения ее устойчивости шаг по времени вычислялся по формуле

где - шаг сетки по пространству.

При данных параметрах шаг по времени составляет . Таким образом, в ходе эволюции в пределах диапазона запаздывания () выполняется не менее 200 шагов по времени. Как было указано выше, прямолинейный подход к данной проблеме подразумевает хранение в машинной памяти данных всех временных слоев для каждого реагента двухкомпонентной системы в пределах диапазона запаздывания. Если расчет производится для данных вещественного типа, требуемая память для хранения значения в каждом элементе 10 байт. Соответственно минимальный размер выделяемой оперативной памяти для хранения одного временного слоя одного компонента равен около 390.63 Кбайт, а всего массива данных равен порядка 76.3 Мбайт. Для двухкомпонентной системы необходимый объем памяти удваивается.

3. Анализ результатов расчетов

Для оценки эффективности разработанного метода был произведен расчет пространственно-временной динамики системы (8-10) с хранением всех опорных слоев в пределах диапазона запаздывания (). В самом начале эволюции, когда хранимых в памяти опорных слоев еще нет, расчет осуществлялся следующим образом если , то для каждого узла сетки возвращалось значение функции, заданной в виде гармонических колебаний, определяемых набором случайных значений частоты, амплитуды и фазы. Этот подход к экстраполяции эволюции системы позволял одновременно задать начальные условия, а также произвести расфазировку в системе. При , в случае совпадения времени со временем одного из хранимых слоев, использовался этот слой. В случае, когда время не совпадало со временем записи ни одного из хранимых слоев, использовалась формула интерполяции Ньютона (7) для определения приближенного значения функции (см. рис. 3).

С точки зрения программной реализации метода хранение опорных слоев организуется с применением динамических массивов по принципу очереди. Таким образом, , где t1, t2 - время первого и второго опорного слоя. Основным критерием эффективности использования данного метода является время расчета модели.

В табл. 3 представлена зависимость времени расчета от параметра . Видно, что в случае, когда в памяти хранится каждый 20-й слой по времени , расчет происходит в 10.5 раз быстрее, причем вне зависимости от технических параметров вычислительной машины.

Таблица 3. Ускорение скорости расчета при использовании разработанного метода

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

Оценка проводилась следующим образом для каждого значения параметра выполнялась серия расчетов с выводом данных с различных узлов сетки. Затем временные ряды для каждого расчета в одинаковых узлах сопоставлялись и находилось среднее значение.

Рис. 5. Погрешность применения метода

Более детально схему оценки результатов методики можно представить следующим образом

Задается фиксированный набор узлов сетки, используемых для оценки (например 10 точек (20;20), (50;50), (80;80), (120;120), (150;150), (180;180), (30;170), (90;130), (140;100), (190;10) ).

Выбирается значение параметра K (например K=1) и производится расчет модели.

По временному ряду каждого узла определяется значение амплитуды, частоты, а также фазы колебаний в конечный момент времени расчета.

Находится среднее значение амплитуды, частоты и фазы для всех временных рядов.

Повтор пунктов 2-5 для остальных значений параметра K (K = [2;5;10;20] и адаптивный метод).

Для результатов пункта 6 вычисляется вносимая методическая погрешность. При этом используется формула относительной погрешности вычислений , где X1 - значение оцениваемого параметра при K=1.

Сводные значения представлены на диаграмме (рис. 5). На основании данных диаграммы можно отметить, что значительное увеличение погрешности вычислений происходит тогда, когда значение параметра K устанавливается менее 5. Если рассматривать колебания концентрации в одном из узлов сетки, то исходя из графика (рис. 6) можно сказать, что наиболее точные результаты могут быть получены при К=2 и с использованием адаптивного метода, однако в первом случае время расчета модели существенно больше.

Рис. 6. Временной ряд для узла [20;20] при использовании различных значений параметра К

Рис. 7. Поле концентрации компонента FRQ в момент времени 1000 часов для различных вариантов расчетов (a) К=1 (хранение 100% слоев в диапазоне запаздывания); (б) К=2 (хранение 50% слоев); (в) К=20 (хранение 5% слоев); (г) адаптивный метод

Представляет интерес, каким образом предлагаемый метод влияет на точность воспроизведения пространственно-временных структур. Как следствие выраженной сильной нелинейности данной системы, даже малые отклонения в периоде, амплитуде и фазе колебаний могут приводить к значительным изменениям в структуре. На рис. 7 представлены поля концентрации белка FRQ после 1000 единиц времени эволюции для четырех различных вариантов расчета. Первый рисунок (рис. 7,а) соответствует стандартному подходу к решению задачи, когда в оперативной памяти компьютера хранятся данные каждого шага по времени в пределах диапазона запаздывания (). При хранении только половины временных слоев () пространственное распределение концентрации реагента практически не изменилось (рис. 7,б). При использовании только каждого 20-го слоя в качестве опорного () результат вычислений меняется качественно (рис. 7,в). Это хорошо видно, если сравнивать выделенные на рисунке участки поля с эталонным полем на рис. 7,а. Вследствие сильной погрешности, вносимой методом, в выделенной области 1 на рис. 7,в значения концентрации белка находятся в противофазе с эталонным значением, а в области 2 вообще не наблюдается развитие спиральной волны, которая наблюдается при эталонном расчете. Схема расчета, использующая адаптивную схему отбора числа опорных слоев показывает превосходный результат - поле концентраций вплоть до мелких деталей совпадает с эталонным (рис. 7,г).

Заключение

Представленный метод позволяет существенно сократить время численного моделирования пространственно-распределенных систем с запаздыванием по времени. В работе рассмотрено применение данного метода к численной реализации модели цепочки химических реакций двух взаимодействующих реагентов, описывающих процесс транскрипции генов. Время расчета данной модели с использованием разработанного метода сокращается на 87%. Наряду с этим полученные результаты имеют вносимую относительную погрешность не более 2% по периоду колебаний и не более 0.6% амплитуды колебаний. Это свидетельствует об эффективности применения данного метода для расчета пространственно-распределенных систем с запаздыванием по времени.

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

Мюррей Дж. Математическая биология Т.1 Введение. Москва-Ижевск Изд-во ИКИ-РХД, 2009.

Ризниченко Г.Ю., Рубин А.Б. Математические модели биологических продукционных процессов. М. Изд-во МГУ, 1993.

Трубецков Д.И., Мчедлова Е.С., Красичков Л.В. Введение в теорию самоорганизации открытых систем. М. Физматлит, 2002.

Гурли С. А., Соу Дж. В.-Х., Ву Дж. Х. Нелокальные уравнения реакции-диффузии с запаздыванием биологические модели и нелинейная динамика // Тр. Междунар. конф. по дифференциальным и функционально-дифференциальным уравнениям ICM-2002 (Москва, МАИ, 11-17 августа, 2002 г.). М. МАИ, 2003. Ч.1. С. 84-120.

Bratsun D., Volfson D., Hasty J., Tsimring L. Delay-induced stochastic oscillations in gene regulation // PNAS. 2005a. Vol.102, №41. P.14593-14598.

Янушевский Р.Т. Управление объектами с запаздыванием. М. Наука, 1978.

Брацун Д.А., Зюзгин А.В., Половинкин К.В., Путин Г.Ф. Об активном управлении равновесием жидкости в термосифоне // ПЖТФ. 2008. Т.34, вып.15. С.36-42.

Нигматулин Р.И. Основы механики гетерогенных сред. М. Наука, 1978.

Мохов И.И., Елисеев А.В., Хворостьянов Д.В. Эволюция характеристик межгодовой климатической изменчивости, связанной с явлениями Эль-Ниньо/Ла-Нинья // Изв. РАН. Физика атмосферы и океана. 2000. Т. 36, № 6. С. 741-751.

Bratsun D.A. Effect of unsteady forces on the stability of non-isothermal particulate flow under finite-frequency vibrations // Microgravity Sci. Technol. 2009. Vol.21. P.153-158.

Брацун Д.А., Захаров А.П. Моделирование пространственно-временной динамики циркадианных ритмов Neurospora crassa // Компьютерные исследования и моделирование. 2011. Т.3, № 2. С.191-213.

Dunlap J. C. Molecular bases for circadian clocks // Cell. 1999. Vol. 96, №. 2. P.271-90.

Loros J.J., Dunlap J.C. Genetic and molecular analysis of circadian rhythms in Neurospora // Annu. Rev. Physiol. 2001. Vol. 63. P.757-794.

Аннотация

Предложен новый алгоритм оптимизации хранения промежуточных полей при численном расчете эволюции пространственно-распределенных запаздывающих систем методом конечных разностей. Алгоритм предполагает хранение в памяти не всех, а только некоторых опорных временных слоев и последующую интерполяцию данных при восстановлении промежуточных слоев. Применение данной методики позволяет производить численные расчеты без использования вычислительных систем с большим объемом оперативной памяти. Эффективность предложенного алгоритма продемонстрирована на примере численного моделирования процессов транскрипции белков, определяющих циркадианные ритмы в клетках.

Ключевые слова пространственно-распределенные динамические системы; запаздывание; метод конечных разностей.

A novel algorithm for optimizing storage of intermediate space fields in the numerical simulation of evolution of spatially extended time-delayed systems by finite difference method is proposed. The algorithm involves storing in the memory not all, but only some selected time layers, and the subsequent interpolation of the data in the construction of intermediate layers. This approach allows the numerical calculations to be implemented on computer systems without large memory. The effectiveness of the developed algorithm is demonstrated by the example of numerical simulation of the transcription of proteins that determine circadian rhythms in cells.