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

74 Вестник СамГУ — Естественнонаучная серия. 2009. № 6(72)

МЕХАНИКА

УДК 519.6

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

© 2009 Л.И. Громова, С.Е. Курушина^ А.А. Иванов2

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

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

1. Предварительные сведения

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

Общепринятой моделью для описания возбудимых сред является система нелинейных параболических уравнений типа ”реакция-диффузия” (см., например [1, 2]).

= ¡к(хк, г) + БХк Ахк, к = 1, 2, 3 ... (1.1)

В отсутствие шумов система типа (1.1) аналитически исследовалась в [3-5]. В ней были обнаружены ДС, которые при существенно различных

1Громова Лидия Ивановна, Курушина Светлана Евгеньевна ([email protected]), кафедра физики Самарского государственного аэрокосмического университета имени акад. С.П. Королева, 443086, Россия, г. Самара, Московское шоссе, 34.

2Иванов Андрей Александрович ([email protected]), студент 5 курса Самарского государственного аэрокосмического университета.

коэффициентах диффузии реагентов приобретают пичковый характер. Численное моделирование ДС при различных граничных условиях проведено в [3, 6, 7].

Однако известно, что под действием случайных полей поведение системы типа ”реакция-диффузия” может качественно меняться. Теоретическое изучение динамики систем вида (1.1) при этом чрезвычайно усложняется, и в этой области существует ограниченное число работ, обычно касающееся линейной стадии образования ДС, например [8, 9]. В работе [10] исследовались уравнения типа (1.1) с нелинейной функцией ¡к(хк,г) (к=3) при условии малости отклонения переменных от стационарных значений, определенных ограничениях на интенсивность шума и ряде других ограничений, серьезно сужающих область применимости полученных результатов.

В данной работе исследуется влияние аддитивного однородного изотропного поля гауссовых флуктуаций динамических переменных на образование ДС при мягком режиме их возбуждения в отсутствие указанных выше ограничений в модели морфогенеза Гирера — Майнхардта [4]. Скорости изменения концентраций автокаталитической и демпфирующей переменных (активатора и ингибитора, соответственно) описываются системой уравнений:

да = р + к а2 - ца + БаАа + Рх(т,г) ( )

Ц = еа2 — pH + БнАН + Е2(т, г). '

Здесь р — скорость образования автокаталитической переменной, ц — ее постоянная распада, Оа — коэффициент диффузии активатора; V — постоянная скорости распада демпфирующей переменной, член са 2 задает скорость образования ингибитора, — коэффициент диффузии ингибитора. Слагаемое ка2/Н учитывает, что увеличение концентрации ингибитора снижает скорость автокаталитического процесса. Случайные функции Рг(т,г) определяют пространственные и временные гауссовы флуктуации скорости изменения концентраций активатора и ингибитора, соответственно, с корреляционными функциями {Рг(г,г)Рг(г',г)) = 9ifi(\f — г^) 5(Ь — т) и нулевыми средними значениями. Взаимная корреляция случайных функций Р^т,г) отсутствует. Величины 9i характеризуют интенсивность флуктуаций, Д(\г — т*\) — некоторая функция, определяющая пространственную зависимость корреляций однородного и изотропного случайного поля [11]; ¿-коррелированность во времени фактически означает, что время корреляции случайного поля гораздо меньше всех характерных времен задачи.

2. Система уравнений для амплитуд неустойчивых мод

Анализ системы (1.2), линеаризованной вблизи однородного стационарного состояния, показывает, что в данной системе существует набор неустойчивых мод [3], которые, нарастая, переходят на нелинейную стадию

развития. Для описания нелинейной стадии эволюции этих мод представим систему (1.2) в виде

д

- — К (А)) р = д + р,

где

К(А) =

К — 1

-у

(р+1)2

± (р +1) —1 + БА

\Р =

91

д2

; д1 — а — ао, д2 — Н — Но,

д1, д2 описывают малые отклонения концентраций относительно

стационарного состояния Но=(р + 1)2 /ц2; ао=(р + 1)/ц. Вектор д

содержит квадратичные и кубические нелинейности, полученные

разложением в ряд правой части (1.2), и определяется так: д=

= ((д1/Но — (9292 + 2ао9192) /Н0+ (2аод1д| + а0д2) /Но — а0д!/Н0), д2). Будем, следуя [3], искать решения системы (1.2) в виде суперпозиции Р(г, г) =

= ^2,к ' О(^)С^?) (г)ё1^кг, где О(^)— собственные векторы линейного оператора

К(А) \ з =1,2 ). р р

Полагая, что неустойчивость возникает при к = кс, получим систему уравнений для амплитуд неустойчивых мод £(1)(г) с точностью до слагаемых третьего порядка малости включительно:

сЦ'

(1)

— А1(кс)«<1) = Ек.Гм е(кЖ’С'т + кС — кс)+

(1) (1)

М ~ ,к’с' Л с!^к'с ^кЦ

+

+

кС ,к

*(кс, Щ¿(кС + кс + кш — кс) +

кс кс кс

Ь(1)Л—^ | Г /О(1)о_Жт:

Х(кс,к)[ /ОЦе-*кгFl(r,t)dr+ ]' О2 ке +

£(1 )б(к'с + к — кс) +

+е(кс,к)^1012le-iГrfFl(r,t)dr + / 022ке- Г?F2(r,t)dr

+ / ОЦ e.-iГifFl(r, + / О^ е-к?F2(f,

(2.1)

Значение критического волнового числа определено в [3]. В уравнениях (2.1) введены следующие обозначения:

е( с) =

ц2051

(1+ Р)2

2, кс 2

(1 + р)

+

х(кс, к) = —

А1( к)

+0Й. (0 1

2^2°(1)с (0( 1 ) 0( 1 ) +

( 1+р)2

—ц+р)0 1к02) + 2021О 1к0 1 %;

е(кс,к =— Аш «о(2к + 021р% —

—(1+Р) (0112с 022к+0212с 012к)) + 20211с 0И 012

01,к с 01 к +

о!1) о!1) -

1, к (1+р)2 2> кс 2,к

) к

С

1

(К к к) — — д °1,кс (б(1М б(1) +___Д-___(б(1М______________— (б(1) ^

а(Кс,К) — (1+р)4 °2,кс + (1+р)2 (1+р) \ °2,кс)

+ (і+рг (х(кс, + е(кс, к)б1к) X

0(1П2 , д2 ^бап2_ б(1) б(1) ) +

б1 к- + С1_1_„12 б2 С1_1_^ б1 к_ б2 +

+

'1, к^ "Г" (1+р)2 уб2, к^ (1+р) б1 ,ксб2 ,кс

(о(1к.с)2 + еОсАоЮ (о<1кс)2.

Здесь бПк и б^П'к (І =1,2) — компоненты собственных векторов линейного оператора К(А) задачи (1.2) и сопряженного к нему, соответственно:

б «> — ( 2<Д„ 1; аи — ( (-1),1+^

д(1+)к2+Л^') / У (-іу

д(1+)к2+Лі)(1+)к2+Л2) 2(р+1)(Лі-Л2)

Здесь

а12 = ОР ±\[Щ^-т;

а(к) = -(Б + 1)к2 + 2ц/(р +1) - ц - 1; в (к) = (к2 + ^)(Бк2 + 1) - 2^Бк /(р + 1).

Функции 5(к' + к" -к) и 5(к' + к" + к"' -к) дают ’’правила отбора”, которые важны для определения возникающих пространственных структур. Моды (^(Ъ) служат параметрами порядка. Их кооперация или конкуренция определяют вид возникающих структур.

Будем полагать, что реакция протекает в двумерном слое, размеры которого ¡1 и ¡2, и примем периодические граничные условия. Волновой вектор определяется как

к = 2к(^т1/1^) , п’т = 0, ±1, ±2...

Заменим вектор кс его модулем кс и углом ф, который этот вектор образует с фиксированной осью. Этот угол изменяется от 0 до п. В первой сумме уравнений (2.1) может остаться только слагаемое, для которого к!с + к'С = кс. Это возможно, если рассматриваемая тройка векторов образует правильный треугольник. Аналогично во второй сумме уравнений (2.1) к'с + Щ + к'” = кс. Следовательно, кС = -к!с;к'<! = кс. Модуль

и угол в волнового вектора устойчивых мод к легко выражаются через кс, ф и ф (ф — угол, который вектор к’с образует с выбранной осью): к = \/2кс\/1 - со§(ф' - ф); 1дв = -сЬд(ф' + ф)/2.

В результате систему уравнений (2.1) для двумерного слоя можно переписать так:

^ - А1(кс)^кС,¥ = с(кс)(кс!^+п/3(кс^-п/3 +

+ £„ «ф - ф'\)('£г «к!!/ + «к!>1(,1(|ф - ин (2‘2)

ъкс,р’

+е ^ ^ ¿2 ~ (^2(\<р - ^\)+б1кс ^(і)+б2і г2,кс(ь)’

где ^1(\ф - ф\)=011х(\ф — ф\)+э^£(ф — ф^ ^(\ф — ф 1)=02,1х(\ф — ф\^

+ 0^£(\ф — ф'\); (^) = I е-г^гЕ, (т,1)д,т — компоненты случайного век-

торного поля £(1), ] и к — индексные аргументы этого поля. Корреляционные функции для компонент поля г(1) имеют вид: (^,к(^)г]к'(т)} =

= 26,Ф,(к) 5(к + к')5(1 — т). Очевидно, ^к^р являются функционалами компонент поля г(1).

Усредним уравнения системы (2.2) [12] и, принимая во внимание ее формальные решения и раскрывая коррелятор к(^)^ с помощью много-

мерного обобщения формулы Фурутцу — Новикова [11] , окончательно получим:

)(2)

¿а

(1) к о,р

¿Ь

Л^к ^ + С(кс){ ^Ц^+п/З ^кС^-п/З) +

(1) к с,ф'

+

(2.3)

+ Ш(^2(1ф — ф'\)61 Ф1(к)+ и2(\ф — ф\)62Ф2(к)) (С'

Здесь

Л = \1(кс) + 6^2(0)Ф1(0) + 62 ^ (0)Ф2(0), (2.4)

Ф12 (0) > 0; <х>2 2 (0) в общем случае комплексные величины. Знак о над суммой означает, что в нее не входит слагаемое с ф = ф. Из выражения (2.4) для собственных чисел Л неустойчивых мод видно, что их действительная часть пропорциональна интенсивности флуктуаций.

а б

Рис. 2.1. Зависимость действительной части Л от волнового числа: а — при разных интенсивностях флуктуаций (к/ = 1); б — при различных

радиусах корреляции (6 = 0,1)

Зависимость действительной части Л от волнового числа к численно исследовалась при различных интенсивностях шума 6 и радиусах корреляции rf. Пусть для простоты 61 = 62 = 6, fг = exp(—kf \г — ■Р'!). Здесь kf — величина, обратная радиусу корреляции. Результаты численных расчетов для двумерной среды в случае, когда становится неустойчивой только одна мода кс (маргинальная) при 6 = 0, и для 61 <62 <63 < 64 приведены на

2

рис. 2.1, а. Из рисунка видно, что при увеличении интенсивности флуктуаций концентраций активатора и ингибитора область неустойчивых мод, для которых Яв(Л) > 0, увеличивается.

На рис. 2.1, б приведена зависимость действительной части Л от волнового числа к при различных радиусах корреляции и интенсивности шума 6 = 0,1. Видно, что при увеличении радиуса корреляции область неустойчивых мод нелинейно увеличивается.

3. Результаты численного моделирования

Для двумерной среды было проведено численное моделирование эволюции системы (1.2) с граничными условиями (ГУ) вида [2]:

(I) тип

да да да да

дх х=0,у дх х=Іі,у ду х,у=0 ду х,у=І2

ОН ОН он ОН

дх х=0,у дх х=Іі,у ду х,у=0 ду х,у=І2

(II) тип

а(0, у, і) = а(ї\,у, і); а(х, 0, і) = а(х, 12,і); Н(0, у, і) = Н(ї\,у, і); Н(х, 0, і) = Н(х, 12, і);

да да да да

дх х=0,у дх 5 х=Іі,у ду х,у=0 ду х,у=І2

ОН ОН ОН ОН

дх х=0,у дх 5 х=Іі,у ду х,у=0 ду х,у=І2

Формирование пичков при наличии шума различной интенсивности и разных граничных условиях представлено на рис. 3.1 и 3.2. Из рис. 3.1 и 3.2, а, б видно, что при увеличении интенсивности флуктуаций процесс формирования ДС начинается раньше вдали от точки бифуркации. Кроме того, из сравнения рис. 3.1, а и 3.1, б можно заметить, что при наличии флуктуаций формируются ДС других периодов, т. е. флуктуации играют роль фактора отбора, выбирая на стадии формирования ДС наиболее устойчивые. Последний результат согласуется с выводами работы [2].

Результаты моделирования вблизи неустойчивости Тьюринга показаны на рис. 3.2, б, в. При выбранных параметрах модели В и ц (см. рис. 3.2, б, в) критическая скорость образования автокаталитической переменной рс= 0,814. Численно исследовались процессы изменения концентраций реагентов при р = 0,8 в отсутствие и при наличии флуктуаций.

Оценим времена разрушения однородного состояния и образования ДС. Пусть и— распределение концентраций реагентов в ДС, ик € Кпхп —

распределение концентраций активатора или ингибитора в момент времени к, а Є Я, г,і = !...,и — концентрация активатора или ингибитора в точке с координатой (і/]) в момент времени к. Введем нормы = = ||и^ — ик|| , зависимость которой от времени определяет близость состояния системы к ДС, и С = ||ик+1 — ик|| , которая показывает бли-

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

п

а 6 = 0,01 б 6 = 0,1

Рис. 3.1. Формирование ДС концентрации активатора и ингибитора при наличии флуктуаций вдали от точки бифуркации (ГУ II). Параметры модели: р = 0,1; 1 = 1; В = 400; kf = 1; а — 6 = 0,01, б — 6 = 0,1. Слой 200x200 единиц длины. Левый вертикальный ряд — концентрация автокаталитической переменной, правый — демпфирующей. В середине указаны соответствующие моменты времени

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

Рис. 3.2. Формирование ДС концентрации активатора при наличии флуктуаций. Слой 50x50 единиц длины (ГУ I): а - вдали от точки бифуркации. Левый вертикальный ряд: параметры модели: р = 0,1; р = 1; Б = 400; kf = 5, 0 = 0,01. Правый вертикальный ряд: р = 0,1; р = 1; Б = 400; kf = 5, 0 = 0,1; б — вблизи точки бифуркации: р = 0,8; р = 1; Б = 400; к/ = 5; 0 = 0,2; в — формирование ДС концентраций реагентов при наличии флуктуаций вблизи точки бифуркации при р = 0,8; р = 1; Б = 400; kf = 5; 0 = 0,2. Слой 200x200 (ГУ II). Левый вертикальный ряд — концентрация автокаталитической переменной, правый — демпфирующей. В середине указаны соответствующие моменты времени

14000 ад 1 і і і активатор _

12000

10000 еооо \ 6=0,05 -

4000 \в=02 V.

2000 0 | \_ . ^ и

50 100 150 200 250 300 380 400

Гг-1

ингибитор

\ ' \в=0,05 .

. \в=0,2 •

50 1 СО 150 200 250 300

Рис. 3.3. Зависимость нормы £к матрицы разности распределений концентрации активатора (ингибитора) от времени при разных интенсивностях флуктуаций:

- вдали точки бифуркации Тьюринга; б — вблизи точки бифуркации

а

а

' 1 1— I 1 активатор

1.5 -

1 ^ 0=0.2,0.1 Ш 0=0.2,Л,= 1 "я \ 0=0 \ ^ ——, Ь

Т> й « 0=0.2, Л;=10 60 80 100 120

..«к* . .

ингибитор в- АН -

:1 1 20 40 \ *к ——— 1 1- * і а 80 100 120

Рис. 3.4. Зависимость нормы Бк матрицы разности распределений концентрации активатора (ингибитора) от времени при разных радиусах корреляции вдали точки бифуркации Тьюринга

Рис. 3.5. Зависимость нормы Ск матрицы разности распределений концентрации активатора (ингибитора) от времени при разных интенсивностях флуктуаций: а — вдали точки бифуркации Тьюринга; б — вблизи точки бифуркации

Зависимость нормы Ск матрицы разности распределений концентрации реагентов от времени позволяет также оценить время формирования ДС по величине временного интервала между началом формирования первого пика и окончанием формирования последнего пика этой зависимости. Например, из рис. 3.5, б видно, что при интенсивности шума 0,2 подобный интервал составляет ~ 90 единиц модельного времени, а при 9=0,05 ~ 70 единиц. Следовательно, вблизи точки бифуркации Тьюринга при росте интенсивности шума время формирования ДС увеличивается.

Заключение

Таким образом, в модели Гирера — Майнхардта вблизи маргинального состояния флуктуации динамических переменных приводят к увеличению области неустойчивых мод (рис. 2.1). Показано, что действительная часть собственных чисел неустойчивых мод пропорциональна интенсивности флуктуаций и некоторой степени радиуса корреляции, зависящей от размерности пространства реакционно диффузионной системы (см. (2.4)). Проведено численное моделирование эволюции рассматриваемой системы при различных граничных условиях. Показано, что вдали от точки бифуркации Тьюринга при увеличении параметров случайной компоненты динамических переменных процессы разрушения однородного состояния и формирования ДС происходят быстрее; вблизи точки бифуркации Тьюринга флуктуации замедляют разрушение однородного состояния и образование ДС. Проведена оценка времен изменения состояния системы путем расчета норм матриц разности распределений концентрации реагентов.

Литература

[1] Cohen D.S., White A.B. Sharp fronts due to diffusion and viscoelastic

relaxation in polymers // SIAM J. Appl. Math. 1991. V. 51. P. 472-483.

[2] Романовский Ю.М., Степанова Н.В., Чернавский Д.С. Математическое моделирование в биофизике (введение в теоретическую биофизику). М.; Ижевск: ИКИ, 2004. 472 с.

[3] Хакен Г. Синергетика. М.: Мир, 1980. 406 с.

[4] Meinhardt H., Gierer A. Generation and regeneration of sequences of structures during morphogenesis // J. Theor. Biol. 1980. V. 85. P. 429.

[5] Keener I.P. Activaters and ingibitors in pattern formation // Stuadies and Applied Mathematics. 1978. V. 59. P. 1.

[6] Gierer A., Meinhardt H. Biological pattern formation involving lateral

inhibition // Lectures on Mathematics in the Life Sciences. 1974. V. 7.

P. 163.

[7] Meinhardt H., Gierer A. Applications of a theory of biological pattern formation based on lateral inhibition // Journ. Cell. Sci. 1974. V. 15. P. 321.

[8] Белинцев Б.Н. Динамические коллективные свойства развивающихся систем: дис. ... канд. физ.-мат. наук. М.: МФТИ, 1979.

[9] Соляник Г.И., Чернавский Д.С. Математические модели морфогенеза // Препринт ФИАН. 1980. № 8.

[10] Михайлов А.С., Упоров И.В. Критические явления в средах с размножением, распадом и диффузией // УФН. 1984. Т. 144. Вып. 1. С. 79-114.

[11] Кляцкин В.И. Стохастические уравнения глазами физика. М.: Физмат-лит, 2001. 528 с.

[12] Свирижев Ю.М., Логофет Д.О. Устойчивость биологических сообществ. М.: Наука, 1978. 354 с.

Поступила в редакцию 29/У/2009; в окончательном варианте — 29/У/2009.

MATHEMATICAL MODELLING OF THE CONTRAST DISSIPATIVE STRUCTURES IN THE FIELD OF FLUCTUATIONS OF DYNAMICAL VARIABLES

© 2009 L.I. Gromova, S.E. Kurushina3 A.A. Ivanov4

The influence of additive homogeneous isotropic field of Gauss fluctuations of dynamical variables of Gierer-Meinhardt model on the formation of dissipative structures in the soft mode of their excitation was investigated. The system of equations for the description of undamped modes interaction was received. It was shown that fluctuations of dynamic variables widen the instability region. The numerical modelling of evolution of the considered model with different boundary conditions was performed.

Key words: dissipative structures, undamped modes, field of fluctuations of dynamical variables, numerical modelling.

Paper received 29/У/2009. Paper accepted 29/V/2009.

3Gromova Lydiya Ivanovna, Kurushina Svetlana Evgenjevna ([email protected]), Dept. of Physics, Samara State Aerospace University, Samara, 443086, Russia.

4Ivanov Andrey Alexandrovich ([email protected]), Samara State Aerospace University, student.