автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Моделирование динамики распределенных систем с запаздывающей обратной связью

<' ' На правах рукописи

Захаров Андрей Павлович

МОДЕЛИРОВАНИЕ ДИНАМИКИ РАСПРЕДЕЛЕННЫХ СИСТЕМ С ЗАПАЗДЫВАЮЩЕЙ ОБРАТНОЙ СВЯЗЬЮ

05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

1 3 КАР 2014

Пермь-2014

005545869

Работа выполнена в ФГБОУ ВПО «Пермский государственный гуманитарно-педагогический университет»

Научный руководитель:

доктор физико-математических наук, доцент, зав. кафедрой теоретической физики и компьютерного моделирования, Брацун Дмитрий Анатольевич

Официальные оппоненты:

Русаков Сергей Владимирович, доктор физико-математических наук, профессор, зав. кафедрой прикладной математики и информатики ФГБОУ ВПО «Пермский государственный национальный исследовательский университет»

Ризниченко Галина Юрьевна,

доктор физико-математических наук, профессор,

зав. сектором информатики и биофизики сложных систем

биологического факультета ФГБОУ ВПО «Московский

государственный университет имени М.В. Ломоносова»

Ведущая организация:

ФГБОУ ВПО «Вятский государственный университет», г. Киров

Защита состоится «22 » апреля 2014 года в 16-00 на заседании диссертационного совета Д 212.188.08 при ФГБОУ ВПО «Пермский национальный исследовательский политехнический университет» по адресу: 614990, г. Пермь, Комсомольский пр., 29, ауд. 345.

С диссертацией можно ознакомиться в библиотеке Пермского национального исследовательского политехнического университета.

Автореферат разослан «2Г» февраля 2014 года

Ученый секретарь

диссертационного совета Д 212.188.08, кандидат физико-математических наук

А.И. Швейкин

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность работы. Если состояние динамической системы в данный момент времени зависит от всей её предыдущей эволюции (или некоторого интервала предыстории, в т.ч. нулевой меры), то такие системы называют наследственными (или запаздывающими). В последнее время системы такого типа вызывают все больший интерес у исследователей в связи с многочисленными приложениями, которые включают популяционную динамику и социальные процессы, нелинейные химические реакции, процессы генной регуляции, поведение систем с автоматическим управлением, механику деформируемых твердых тел с памятью, механику реологических сред и т.д. Запаздывание может быть обусловлено самыми различными причинами: ограниченностью скорости распространения сигнала (например, электромагнитной волны), растянутостью изучаемого процесса во времени (воспроизводство популяции), наличием инерционности некоторых элементов (в теории управлении с обратной связью), существованием длинных цепочек многоэтапных реакций с известным результатом в конце (в процессах генной регуляции) и т.д.

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

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

Цель данной работы заключается в разработке и применении эффективного вычислительного метода для математического моделирования пространственно-распределенных динамических систем с запаздывающей обратной связью.

Для достижения поставленной цели потребовалось решение следующих основных задач:

1. Разработка и программная реализация нового адаптивного численного мето-

да для расчета динамики систем с наследственной или запаздывающей обратной связью.

2. Тестирование разработанного численного метода на задачах, имеющих точные решения или решения, полученные на основе известных численных методов. Анализ и сравнение результатов тестов.

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

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

Создан программный комплекс для проведения вычислительных экспериментов.

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

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

Эффективность предложенного метода продемонстрирована на примере модели реакции-диффузии белков, ответственных за биоритмы в живой ткани. Впервые исследованы вопросы пространственной синхронизации циркадианных ритмов в рамках как непрерывной, так и дискретно-непрерывной модели живой ткани.

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

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

Исследования, представленные в диссертационной работе, выполнены при поддержке Российского Фонда Фундаментальных Исследований (14-01-96021

р_урал_а, 14-01-96022 р_урал_а), Министерства образования и науки РФ (ГК № 1.3103.2011), Министерства образования и науки Пермского края (С-26/244) и Пермского государственного гуманитарно-педагогического университета (проект № 031-Ф Программы стратегического развития университета).

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

Основные положения, выносимые на защиту:

1. Адаптивный метод расчета пространственно-распределенных систем с запаздывающей обратной связью.

2. Анализ точности и эффективности предложенного метода.

3. Модель циркадианных ритмов в живой ткани, основанная на запаздывании в процессах генной регуляции.

4. Применение программного комплекса на основе адаптивного метода для исследования синхронизации циркадианных ритмов в живой ткани.

Апробация работы. Результаты исследования докладывались на VI Всероссийской научной конференции «Математическое моделирование развивающейся экономики, экологии и биотехнологий» (Киров, 2011); ХХ-ой Всероссийской школе-конференции «Математическое моделирование в естественных науках» (Пермь, 2011); European Conference on Complex System (Brussels, Belgium, 2012); Всероссийской научно-практическая конференция с международным участием «Актуальные проблемы механики, математики, информатики — 2012» (Пермь, 2012); XX Международной конференции «Математика. Компьютер. Образование» (Пущино, 2013); Международной школе-конференции «Теоретическая биофизика. Анализ и моделирование», (Пущино, 2013); Interdisciplinary Symposium on Complex Systems (Prague, Czech Republic, 2013); European Conference on Complex System (Barcelona, Spain, 2013). Полностью диссертационная работа докладывалась и обсуждалась на семинарах: кафедры теоретической физики и компьютерного моделирования и кафедры общей и экспериментальной физики Пермского государственного гуманитарно-педагогического университета (рук. проф. В.Г. Козлов), Института механики сплошных сред УрО РАН (рук. акад. РАН В.П. Матвеенко), кафедры математического моделирования систем и процессов (рук. проф. П.В. Трусов), кафедры механики композиционных материалов и конструкций (рук. проф. Ю.В. Соколкин) Пермского национального исследовательского политехнического университета.

Публикации. Основные положения диссертации опубликованы в 20 публикациях, в том числе в 7 публикациях в изданиях, рекомендованных ВАК для публикации материалов диссертации.

Личный вклад автора. Разработка метода и анализ полученных результатов произведены совместно с научным руководителем; адаптивный метод хранения данных предложен автором; создание программного комплекса, его тестирование и проведение вычислений осуществлены автором.

Структура и объём работы. Диссертационная работа состоит из введения, трех глав, заключения, списка использованной литературы и приложения.

Работа изложена на 155 страницах, содержит 35 иллюстраций и 2 таблицы. Список литературы включает 126 наименований.

ОСНОВНОЕ СОДЕРЖАНИЕ РАБОТЫ

Во введении обоснована актуальность исследования, определены цели и задачи исследования, показана научная и практическая значимость работы.

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

Вторая глава посвящена описанию и тестированию нового эффективного численного метода для расчета динамики пространственно-распределенных систем с наследственной или запаздывающей обратной связью.

В п. 2.1. выполнена постановка задачи и описание разработанного численного метода. Для этого рассмотрена система безразмерных нелинейных кинетических уравнений, описывающих поведение пространственно-распределенной динамической системы диффузионного типа:

эи_,

д!

о

+£>Ли, (1)

где и - вектор, содержащий набор физических (химических, биологических и т.д.) переменных в «-мерном пространстве, который характеризует нелинейную кинетику системы; В - матрица коэффициентов пространственной диффузии; Д - оператор Лапласа. Нелинейная функция / может содержать слагаемые, зависящие только от текущего времени /, запаздывающие слагаемые, а также наследственные интегральные слагаемые, обусловленные всей предыдущей эволюцией системы. Для замыкания система уравнений (1) должна быть дополнена граничными и начальными условиями.

Для простоты в работе рассматривается явная конечно-разностная схема, хотя это не является обязательным для предлагаемого метода. Для проведения численного исследования в рассмотрение вводится трехмерная пространственная сетка (используется равномерное разбиение):

х,=т, (/ = МГ); и = 0Жу); гк=ки, (* = МГ), (2)

где к - шаг сетки; , Ыу, - количество узлов вдоль х, у, г. С помощью (2) уравнения (1) могут быть записаны в конечно-разностной форме, где используются центральные разности для пространственных производных и разность вперед для производной по времени с шагом по времени в:

и::, +^/(1^,и- -кт;: »+

+ии* -би"м)\ (3)

1+и* ил |,./-и ' 1..М+1 1

Согласно (3) для вычисления поля на (и + 1)-м шаге по времени требуется информация о поле не только в момент времени и, но и в момент и - X в про-

шлом (при запаздывающей обратной связи) либо о всех предыдущих моментах времени от 1-го до и-го (при наследственной обратной связи). Если общее время эволюции / (или время запаздывания т) достаточно велико, то велико и I, поскольку необходимо сохранять в памяти компьютера данные для всех временных слоев, участвующих в построении решения. Ситуация усугубляется еще и тем, что при использовании явной схемы, которая более предпочтительна для расчета эволюционирующих систем, очень важно контролировать шаг по времени в\ для сохранения устойчивости схемы он должен быть меньше заданного значения в * - числа Куранта:

Ь2

= (4)

2 тах и[

К примеру, если размер интегрируемой области / = 50 у.е. по пространственным переменным х,у,г и рассматривается однородная сетка 100x100x100, а время запаздывания в системе т = 10, то должно выполняться неравенство в < в *, где в* = 0.125 (при условии тах|и| = 1, что довольно часто имеет место, например, в задачах реакции-диффузии). Тогда в диапазон времени запаздывания укладывается не менее 80 фиксированных шагов по времени. Если каждый элемент сеточной функции определяется вещественным числом, занимающим 8 байт памяти, то для хранения данных потребуется выделить порядка 610,4 Мб памяти. Организация хранения данных на внешнем носителе неэффективна, в силу наличия значительных затратах времени на доступ к этим устройствам. Выгрузка данных в кэш-память процессора также не возможна, так как эта память имеет ограниченный объем. Таким образом, ключевым фактором, влияющим на скорость проведения расчетов, становится оптимизация записи и считывания из оперативной памяти больших объемов данных на каждом шаге по времени.

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

г

Рисунок 1 - Схема выделения опорных временных слоев (отмечены аббревиатурой ОС): / - текущее время расчета, т - время запаздывания, в - шаг по времени

Вводится в рассмотрение важный для метода параметр К, задающий количество опорных слоев в интервале запаздывания. Значение этого параметра устанавливается в начале расчета.

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

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

Для оценки гладкости функции U введем функцию G, которая описывает относительную величину максимального изменения U в момент времени t + в:

+ (5)

Пусть, для определенности, в пространственно-распределенной системе (1) имеется только запаздывающая обратная связь. Адаптивность метода обеспечивается введением в рассмотрение уже не фиксированного числа опорных слоев, а функции времени K(t):

К=

0(1 + exp(20G))_

(6)

Функция К принимает только целочисленные значения, которые означают, что в потоке данных опорным становится каждый К -й слой. Поэтому при использовании выражения (6), полученного эмпирически для определенного класса дифференциальных операторов, рассматриваемых далее, предполагается округление значения К в меньшую сторону.

Значение К = 1 соответствует отказу от применения адаптивного метода и переходу к прямолинейному запоминанию всего потока данных. Соотношение (6) позволяет задать степень детализации информационных полей, хранимых в памяти компьютера, в зависимости от скорости изменения вектора-функции и, которая определяется переменной С(/) (рисунок 2).

г

Рисунок 2 - Зависимость количества опорных слоев К от гладкости функции G(t) в случае использования адаптивного метода.

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

В п.2.3. выполнена апробация и оценка погрешности метода на основе его применения для моделирования систем трех различных типов.

П.2.3.1. содержит описание результатов тестового расчета системы функционально-дифференциальных уравнений с постоянным запаздыванием, имеющей точное решение (Zubik-Kowal В., Vandewalle S. Waveform relaxation for functional-differential equations // SIAM J. Sei. Comput. 1999. V.21, P.207-226): Э U(t,x) 1 d2U(t,x)

+ e" U(t-l,x) + g(t,x)

(7)

Э/ 10 + 40/2 Э

доя (е [0,10], л-е [-5,5]. Функция g(t,x) и начальные условия подобраны таким образом, чтобы точное решение уравнений (7) определялось как

[/(?,*) = (8) Выполнена оценка методической погрешности (рисунок 3):

£(/) = шах, и""" (I, ¡К) - и(1, Ш) , (9)

где С/"1""0,г/г) - численное значение функции в момент времени I в узле г; II(г, ¿к) - соответствующее точное значение функции в момент времени (.

Ш

0.01

0.001

0.0001

1S-005

.........Г 1 ' Адаптивный метод •• ».....

K=50 - • -

K=10 - -

Zubik-Kowal (Jac/Pic, К-.20) ,. ж . . ,.

• Zubik-Kowai {G-S/Pic, k»15) ■ ш ■ ■ ---------

i ; * .......- - - - ».....- - •• г... • ......» • -

г ' - , -Ol v" £3- —Q ~ - 0- - - * - ~ ... :* • • » • • « ■ * ■ : m- '- - - 4

— - ■ - - > • •••«•- -...................-i.........................................*..............i,..............................................■...... . ........

10

{

Рисунок 3 - Методическая ошибка Е(Г) как функция времени для уравнения (7)

Показано, что погрешность разработанного адаптивного метода не превышает погрешность вычислений произведенных в работе Зубик-Коваль с приме-

нением техники волновой релаксации Якоби-Пикара и Гаусса-Пикара в функциональном пространстве; относительная погрешность метода составила не более 0.01%.

В п. 2.3.2. выполнено тестирование метода на примере численного решения системы уравнений типа «реакция - диффузия», содержащего нелокальное слагаемое с дискретным запаздыванием. Модель предложена в работе Ли, Зу (Li J., Zou X. Modeling spatial spread of infectious diseases with a fixed latent period in a spatially continuous domain // Bull. Math. Biol. 2009. V.71. P.2048-2079) для исследования распространения инфекционных заболеваний с фиксированным латентным периодом.

Система уравнений записана следующим образом:

dS(t,x) n d2S(t,x) JO/ ч

+ —- , - - dS(t, х) - rf(t, x)S(t, x),

dt dl(t,x) dt

= D,

dx2 d2I(t,x)

(10)

dx2

pl{t, X) + £ J rl(t - T, y)S(t - r, y)fa (x-y)dy, (11)

где t>0, xe9t,/„(*) =

1 -

ju, r ,/j ,£ ,Ds,Dt - параметры модели.

л/4 па

Уравнения (10,11) имеют волновое решение (рисунок 4), при определенном наборе параметров. В работе Ли, Зу проблема начальных и граничных условий была преодолена при помощи метода прямых с использованием специальной математической техники, позволяющей избавиться от нелокальности и применимой только для данной системы уравнений.

Рисунок 4 - Волновое решение задачи (10,11), полученное при помощи адаптивного метода.

Показано, что адаптивный метод нахождения численного решения для системы (10,11) позволяет получить результат, идентичный результатам работы Ли, Зу. Установлено, что относительная погрешность метода, в частности, для показателя скорости движения волны, составила не более 0.91%. Кроме того, если численный метод авторов модели применим только к системам с одной пространственной переменной, в настоящей работе задача (10,11) была переформулирована для двух пространственных переменных и численно получено волновое решение задачи на плоскости.

В п. 2.3.3. выполнена проверка эффективности метода на примере численного расчета системы уравнений типа «реакция - диффузия» с дискретным запаздыванием на плоскости (подробнее - в главе 3):

Э ^ _ 1

3/ 1 +

ЭРГ _ 1

э7 _1+4л:Гг

А№+к1

к^к'Ж-и-т) к^к^'а-т)

+ ОА

-БА Ж.

(12)

(13)

1 + «^2(/--г)

где ^, й7 обозначены концентрации белков и ^СС в момент времени / в точке с координатами (х,у); К^, кЕ Я,, уг)г, к - параметры задачи; Дху - оператор Лапласа в двух пространственных измерениях х и у, О - коэффициент диффузии. Модельные уравнения (12,13) записаны в естественных единицах.

Проведена оценка эффективности разработанного метода (Таблица 1). Для этого произведен расчет пространственно-временной динамики концентраций белков в системе (12, 13) с хранением всех данных в пределах диапазона запаздывания (К = 1). Затем выполнено сравнение с результатами, полученными при других значениях параметра К и при использовании адаптивного метода. Таблица 1 - Ускорение расчета для (12,13) при использовании различных способов хранения

Способ хранения данных

Послойный с значением К Адаптивный метод

1 2 5 10 20

Полное расчетное время, с 7012 2320 913 720 665 909

Ускорение процедуры расчета, раз 1 3,02 7,68 9,74 10,54 7,71

Выполнена оценка вносимой методической погрешности (рисунок 5). Проведен сравнительный анализ пространственно-временных структур, получаемых при разных количествах опорных слоев.

В Ошоснтшьная погрешность периода колебаний, % ■ Ошоситепьная погрешность амплитуды колебаний, % 24.2

□ Сдвиг фазы колебаний, часы __

12.29

I 2 5 1(1 20 Адаптивный

метод

Рисунок 5 - Погрешность применения метода при различных степенях детализации хранимых слоев; графы 1, 2, 5, 10, 20 соответствуют фиксированным значениям К; адаптивный метод подразумевает изменение К согласно выражениям (5-6)

Третья глава работы посвящена моделированию пространственно распределенных систем с запаздыванием по времени.

В разделе 3.1. выполнено численное моделирование процессов реакции-диффузии белков, отвечающих за биоритмы в клетках живых организмов. В п.3.1.1. предложена модель циркадианных ритмов (12,13), которая является дальнейшим развитием существующих в литературе моделей, однако построена не феноменологически, а напрямую выведена из цепочки биохимических реакций, происходящих при транскрипции генов. Модель основывается на нелинейном взаимодействии белков FRQ и WCC, кодируемых генами frequency и white collar, и включает в себя как положительную, так и отрицательную петлю обратной связи. Главным элементом механизма колебаний является эффект запаздывания в биохимических реакциях транскрипции генов.

В п.3.1.2. представлены результаты моделирования. Проведено исследование пространственно-временной динамики системы, которое позволило выявить большое разнообразие и сложность возникающих колебательных режимов. На рисунке 6 показаны результаты моделирования системы при воздействии базального фактора транскрипции А (наличие потока извне одного из белков). Обнаружен эффект возникновения пространственно-временной синхронизации колебаний концентраций реагентов при локальном воздействии данного фактора на систему (рисунок 7).

Рисунок 6 - Поля концентраций белка Р в системе (12,13) после 5000 единиц времени эволюции, возникающие в зависимости от значений параметра базального фактора транскрипции. Расчет производился при случайных начальных условиях.

Рисунок 7

синхронизация

системе (12,

воздействии

транскрипции.

соответствуют

t = 5000, 7200,

- Пространственная колебаний белка /•' в 13) при локальном базального фактора Кадры эволюции моментам времени 11000, 20000 ч.

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

группе клеток может происходить, в том числе, через диффузионные потоки между клетками (рисунок 8). Изучен вопрос о влиянии толщины межклеточного пространства на процессы синхронизации ритмов.

а)

б)

и к т ' ! в. . 1 т 1 ! 1 т ЙГ-VO

ivy® Ж ж . . , т? ^ г

Ш ш ш. "♦.л* г Ш ш I щ ш

1 nvjiv О - ^nnciivi Г11ЧС4 lipuv ipunv 1 uvuuui't v f ..«..-V----f--.---- -------- - \ -

для случая взаимодействия четырех «клеток», а) Процессы транскрипции/трансляции происходят только в одной «клетке», б) Все четыре «клетки» функционируют, вырабатывая белок. Последовательно слева направо представлены моменты времени t=200, 5000, 10000. 20000 ч.

В п.3.2. рассмотрено применение разработанного метода для расчета дискретной распределенной системы, моделирующей одновременное взаимодействие более 1600 клеток. В п.3.2.1 выполнена постановка задачи. Для описания механизма обмена сигналами между клетками использована хемомеханическая модель растущей ткани эпителия, предложенная в работе Салма-Письмена (Salm М., Pismen L.M. Chemical and mechanical signaling in epithelial spreading // Phys. Biol. 2012. V.9. P.026009-026023). Модель эпителия дополнена расчетом биоритмов в каждой клетке сообщества согласно следующим уравнениям:

dF __ _1_ dt ~(\+4K{'F)

dW 1

\ + K"K';W-(t-T)

-yFF-kFW

kur

-y^W-kF.W,

j

i = \,W) (15)

л 0+4к['ю^г \+к';к:;р'{1-т) где /Г, Щ обозначают концентрации белков и IVСС в момент времени г в /-ой клетке пространства; , кГ №, уГ№, к - параметры задачи;, а -

коэффициент переноса белка Г через клеточную мембрану между г -ой и у -ой клеткой, длиной ¿,;. Клетки составляют везде непрерывную ткань эпителия.

Приведены результаты моделирования. Обнаружен эффект кластеризации - появления в пространстве клеток доменов сложной формы, внутри каждого из которых происходит выравнивание фаз колебаний (рисунок 9).

мдши шинысычи

Рисунок 9 - Эффект кластеризации циркадианных ритмов в пространстве более 1600 клеток. Слева направо последовательно представлены моменты времени 1=0, 25, 50, 75, 100, 125 ч. для концентрации белка РЯ() в системе (14,15).

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

СО)

Рисунок 10 - Эффект бегущей волны концентрации белка в системе (14,15).в

пространстве, состоящего из более 1600 клеток. В качестве начальных условий использовано возмущение концентрации белка в одной клетке в центре сообщества. Последовательно представлены моменты времени 1=0, 20, 40, 60. 80, 100 ч.

В заключении сформулированы основные результаты работы.

В приложении приведен скрипт текста программы, реализующей адаптивный метод.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ И ВЫВОДЫ

1. Разработан новый метод оптимизации хранения информации при численном расчете пространственно-распределенных систем с наследственной или запаздывающей обратной связью. Метод предполагает хранение в памяти данных не со всех, а только с некоторых, опорных временных слоев, и последующую интерполяцию для восстановления данных с промежуточных слоев.

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

3. Представлено три численных примера применения метода адаптивного хранения наследственных данных:

- решение уравнения с постоянным дискретным запаздыванием, имеющего точное решение;

- расчет системы интегрально-дифференциальных уравнений с пространственным нелокальным запаздыванием;

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

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

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

5. Выполнено численное исследование предложенной модели пространственно-распределенных циркадианных ритмов. Обнаружены сложные колебательные режимы, отражающие процессы синхронизации ритмов в живой ткани.

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

СПИСОК ОСНОВНЫХ РАБОТ, ОПУБЛИКОВАННЫХ АВТОРОМ ПО ТЕМЕ ДИССЕРТАЦИИ Публикации в изданиях, рекомендованных ВАК

1. Брацун Д.А., Захаров, А.П. Моделирование пространственно-временной динамики циркадианных ритмов N. Crassa // Компьютерные исследования и моделирование. - 2011. - Т.З, №2. - С.191-213. (из перечня ВАК)

2. Захаров А.П., Брацун Д.А. Адаптивный алгоритм хранения полей при расчете динамики сплошной среды с наследственной или запаздывающей обратной связью // Вычислительная механика сплошных сред - Computational continuum mechanics. - 2013. - Т.6, №2. - С. 198-206. (из перечня ВАК)

3. Захаров А.П., Брацун Д.А. Синхронизации циркадианных ритмов в масштабах гена, клетки и всего организма // Компьютерные исследования и моделирование. - 2013. - Т.5, №2. - С.255-270. (из перечня ВАК)

4. Bratsun D.A., Zakharov А.Р. Modeling of spatially extended delay-induced cir-cadian oscillations synchronized by cell-to-cell communications // Proceedings of the European Conference on Complex Systems 2012. Springer International Publishing. - 2014. - P.445-452. (международная база цитирования Springer)

5. Bratsun D., Zakharov A. Adaptive numerical simulations of reaction-diffusion systems with history and time-delayed feedback // ISCS 2013: Interdisciplinary Symposium on Complex Systems: Emergence, Complexity and Computation. Springer International Publishing. - 2014. - Vol.8. - P.70-81. (международная база цитирования Springer)

6. Bratsun D., Zakharov A. Deterministic modeling spatio-temporal dynamics of delay-induced circadian oscillations in N. crassa. // ISCS 2013: Interdisciplinary Symposium on Complex Systems: Emergence, Complexity and Computation. Springer International Publishing. - 2014. - Vol.8. - P.82-90. (международная база цитирования Springer)

7. Zakharov A., Bratsun D. Synchronization of circadian rhythms at scale of gene, cell and whole organism // ISCS 2013: Interdisciplinary Symposium on Complex Systems: Emergence, Complexity and Computation. Springer International Publishing. - 2014. - Vol.8. - P. 162-173. (международная база цитирования Springer)

Публикации в прочих изданиях

8. Bratsun D., Zakharov A. Modeling of spatially extended delay-induced circadian oscillations synchronized by cell-to-cell communications // European conference on complex system: book of abstracts. - Belgium, Brussels: Université Libre de Bruxelles, 2012.-P.67.

9. Брацун Д.А., Захаров А.П. К вопросу о численном расчете пространственно-распределенных динамических систем с запаздыванием по времени // Вестник Пермского университета: Математика. Механика. Информатика. - 2012.

— №4. — С.32-42.

10.Bratsun D.A., Zakharov А.Р. Adaptive numerical simulations of reaction-diffusion systems with time-delayed feedback // Interdisciplinary Symposium on Complex Systems: book of abstracts. - Czech Republic, Prague: Czech Technical University Press, 2013. - P.21.

11.Zakharov A., Bratsun D. Modeling dynamics of biorhythms at different scales of biological system // European Conference on Complex System: book of abstracts.

- Spain, Barcelona, 2013. - P.60-61.

Подписано в печать 18.02.2014. Формат 60x90/16. Усл. печ. л. 1,0. Тираж 100 экз. Заказ № 625/2014.

Отпечатано с готового оригинал-макета в типографии издательства Пермского национального исследовательского политехнического университета. Адрес: 614990, г. Пермь, Комсомольский пр-т, 29, к. 113. Тел. (342) 219-80-33.