Г.Ю.Ризниченко. Лекции по математическим моделям в биологии

ВВЕДЕНИЕ. МАТЕМАТИЧЕСКИЕ МОДЕЛИ В БИОЛОГИИ Понятие модели. Объекты, цели и методы моделирования. Модели в разных науках. Компьютерные и математические модели. История первых моделей в биологии. Современная классификация моделей биологических процессов. Регрессионные, имитационные, качественные модели. Принципы имитацуионного моделирования и примеры моделей. Специфика моделирования живых систем.

Компьютеры в современном мире стали привычными для человеческой деятельности: в финансовой сфере, в бизнесе, промышленности, образовании, сфере досуга. Благодаря компьютерам западной цивилизации удалось существенно продвинуться в следующих направлениях. 1. Автоматизация трудовой деятельности во всех сферах 2. Информационная революция. Возможность хранить и структурировать огромные и самые разнообразные массивы информации и производить быстрый и эффективный поиск необходимой информации. 3. Прогнозирование. Компьютер позволяет строить имитационные модели сложных систем, проигрывать сценарии и делать прогнозы. 4. Оптимизация. Любая человеческая деятельность, в том числе обыденная жизнь требует постоянной оптимизации действий. В процессе эволюции сформировались биологические системы, которые оказываются оптимальными в том или ином смысле, например, в смысле наиболее экономичного использования энергии. Для того чтобы формализовать целевую функцию, то есть ответить на вопрос, что же является для системы оптимальным, необходимо сформулировать модель оптимизируемого процесса и критерии оптимизации. Компьютер позволяет проектировать и реализовать различные алгоритмы оптимизации.

Компьютер работает не с самой системой, а с моделью. Что же такое МОДЕЛЬ? Наиболее простой и общий ответ на этот вопрос: модель – это копия объекта, в некотором смысле «более удобная», допускающая манипуляции в пространстве и во времени. При моделировании, выборе и формулировке модели, определяющими обстоятельствами являются объект, цель и метод (средства) моделирования. В нашем курсе объектами моделирования будут биологические процессы разного уровня организации. Методами моделирования служат методы динамической теории систем. Средства – дифференциальные и разностные уравнения, методы качественной теории дифференциальных уравнений, компьютерная симуляция. Цели моделирования: 1. Выяснение механизмов взаимодействия элементов системы 2. Идентификация и верификация параметров модели по экспериментальным данным. 3. Оценка устойчивости системы (модели). Само понятие устойчивости требует формализации. 4. Прогноз поведения системы при различных внешних воздействиях, различных способах управления и проч. 5. Оптимальное управление системой в соответствии с выбранным критерием оптимальности.

Примеры моделей. 1. Портрет дамы. Пусть некто заказывает художнику написать портрет любимой женщины. Рассмотрим объект, метод (средства) и цель моделирования. Объектом моделирования является женщина. Метод (средства) – краски, кисти, холст. Эмаль, если портрет будет сделан на медальоне, как это было принято в прошлые века. Фотоаппарат и пленка. Рекламный щит, если некто хочет, чтобы его даму видели все, кто проезжает по оживленной магистрали. Обложка журнала, или экран телевизора. Наконец, сам художник, фотограф или рекламное агентство в лице своих дизайнеров. Цель. При моделировании целью, как правило, является манипуляция с пространством и временем. Сохранить облик дамы во времени. Повесить портрет в гостиной, или медальон с изображением любимой - на шею, как это делали в старину. Чтобы потомки восхищались красотой дамы и своим пращуром, которому удалось запечатлеть такую красоту. Другая цель – воспроизведение изображения (модели) объекта с целью сделать модель доступной некоторому кругу людей. Или многократно тиражировать, если некто хочет, чтобы образ дамы увидели миллионы.

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

Примеры моделей. 3. Аквариум является примером физического моделирования. В аквариуме можно моделировать водную экосистему – речную, озерную, морскую, заселить ее некоторыми видами фито- и зоопланктона, рыбами, поддерживать определенный состав воды, температуру, даже течения. И строго контролировать условия эксперимента. Какие компоненты естественной системы будут воспроизведены, и с какой точностью, зависит от цели моделирования.

Примеры моделей. 4. Выделенные из листьев хлоропласты. На выделенных системах часто изучают процессы, происходящие в живой системе, в этом смысле фрагмент является моделью целой живой системы. Выделение более простой системы позволяет исследовать механизмы процессов на молекулярном уровне. При этом исключается регуляция со стороны более высоких уровней организации, в данном случае, со стороны растительной клетки, листа, наконец, целого растения. В большинстве случаев наблюдать процессы на молекулярном уровне в нативной (ненарушенной) системе не представляется возможным. Говорят, что изученные на выделенном хлоропласте первичные процессы фотосинтеза являются моделью первичных процессов фотосинтеза в живом листе. К сожалению, этот метод фрагментирования приводит к тому, что «…живой ковер жизни распускается по ниточкам, каждая ниточка досконально изучается, но волшебный рисунок жизни оказывается утрачен» (лауреат Нобелевской премии по биохимии Л. Поллинг).

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

Примеры моделей. 6. Популяция дрозофилы, является классическим объектом моделирования микроэволюционного процесса и примером исключительно удачно найденной модели. Еще более удобной моделью являются вирусы, которые можно размножать в пробирке. Хотя не вполне ясно, справедливы ли эволюционные закономерности, установленные на вирусах, для законов эволюции высших животных. В лекции 11 мы увидим, что хорошей моделью микроэволюционных процессов являются также микробные популяции в проточном культиваторе. Из приведенных примеров видно, что любая физическая модель обладает конкретными свойствами физического объекта. В этом ее преимущества, но в этом и ее ограничения.

Компьютерные модели содержат "знания" об объекте в виде математических формул, таблиц, графиков, баз данных и знаний. Они позволяют изучать поведение системы при изменении внутренних характеристик и внешних условий, проигрывать сценарии, решать задачу оптимизации. Однако каждая компьютерная реализация соответствует конкретным, заданным параметрам системы. Наиболее общими и абстрактными являются математические модели. Математические модели описывают целый класс процессов или явлений, которые обладают сходными свойствами, или являются изоморфными. Наука конца 20 века – синергетика, показала, что сходными уравнениями описываются процессы самоорганизации самой разной природы: от образования скоплений галактик до образования пятен планктона в океане.

Если удается сформулировать «хорошую» математическую модель, для ее исследования можно применить весь арсенал науки, накопленный за тысячелетия. Недаром многие классики независимо высказывали одну и ту же мудрую мысль: «Область знания становится наукой, когда она выражает свои законы в виде математических соотношений» С этой точки зрения самая "научная" наука физика. Она использует математику в качестве своего естественного языка. Все физические законы выражаются в виде математических формул или уравнений. В химию математика пришла в тридцатые годы 20 века вместе с химической кинетикой и физической химией. Сейчас книги по химии, в особенности по химической кинетике, физической химии, квантовой химии полны математическими символами и уравнениями. Чем более сложными являются объекты и процессы, которыми занимается наука, тем труднее найти математические абстракции, подходящие для описания этих объектов и процессов. В биологию, геологию и другие «описательные науки» математика пришла по настоящему только во второй половине 20 века.

Первые попытки математически описать биологические процессы относятся к моделям популяционной динамики. Эта область математической биологии и в дальнейшем служила математическим полигоном, на котором "отрабатывались" математические модели в разных областях биологии. В том числе модели эволюции, микробиологии, иммунологии и других областей, связанных с клеточными популяциями. Самая первая известная модель, сформулированная в биологической постановке, знаменитый ряд Фибоначчи, который приводит в своем труде Леонардо из Пизы в 13 веке. Это ряд чисел, описывающий количество пар кроликов, которые рождаются каждый месяц, если кролики начинают размножаться со второго месяца и каждый месяц дают потомство в виде пары кроликов. Ряд представляет последовательность чисел: 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89,....,

Пример Пара, готовая к размножению Пара, неготовая к размножению

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

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

Сейчас биологические журналы полны математическими формулами и результатами компьютерных симуляций. Имеются специальные журналы, посвященные работам в области математических моделей: Theoretical Biology; Biosystems; Mathematical Ecology, Mathematical biology, Biological systems etc. Работы по математическому моделированию печатаются практически во всех российских биологических журналах: Общая биология, Биофизика, Экология, Молекулярная биология, Физиология растений и других. В основном, модели являются инструментом изучения конкретных систем, и работы по моделированию печатают в журналах, посвященных той области биологии, к которой относится объект моделирования. Это означает, что модель должна быть интересна, полезна и понятна специалистам-биологам. В то же время она должна быть, естественно, профессионально сделана с точки зрения математики. Наиболее успешные модели сделаны в содружестве специалистов математиков, или физиков, и биологов, хорошо знающих объект моделирования. При этом наиболее трудная часть совместной работы это формализация знаний об объекте (как правило, в виде схемы) на языке, который может затем быть переформулирован в математическую или компьютерную модель.

Классификация моделей Условно все математические модели биологических систем можно разделить на: Регрессионные Качественные Имитационные

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

Примеры 1. Зависимость между количеством производителей хамсы S и количеством молоди от каждого нерестившегося производителя в Азовском море (используется в большой имитационной модели динамики рыбного стада Азовского моря, Горстко, 1985): S - количество сеголеток (штуки) на каждого нерестившегося производителя. x - количество зашедших весной из Черного моря в Азовское производителей хамсы (млрд штук). - среднеквадратичное отклонение. 2. Скорость поглощения кислорода опадом листьев (Из книги: Д.Джефферс "Введение в системный анализ: применение в экологии", М., 1981) Y поглощение кислорода, измеренное в мкл(0,25 г) -1 ч -1. D - число дней, в течение которых выдерживались образцы, B - процентное содержание влаги в образцах, Т - температура, измеренная в град.С. Эта формула дает несмещенные оценки для скорости поглощения кислорода во всем диапазоне дней, температур и влажностей, которые наблюдались в эксперименте, со средним квадратичным отклонением в поглощении кислорода, равном =

Коэффициенты в регрессионных моделях обычно определяются с помощью процедур идентификации параметров моделей по экспериментальным данным. При этом чаще всего минимизируется сумма квадратов отклонений теоретической кривой от экспериментальной для всех точек измерений. Т.е. коэффициенты модели подбираются таким образом, чтобы минимизировать функционал: Здесь i номер точки измерения, x e 'экспериментальные значения переменных, х t теоретические значения переменных, a 1, a 2... параметры, подлежащие оценке, w i "вес" i-го измерения, N число точек измерения.

Имитационные модели (simulation) По меткому выражению Р. Шеннона (1978) имитационное моделирование это нечто промежуточное между искусством и наукой, направление, появление которого целиком обязано бурному росту возможностей вычислительной техники. Суть имитационного моделирования заключается в исследовании сложной математической модели с помощью вычислительных экспериментов и обработки результатов этих экспериментов. При этом, как правило, создатели имитационной модели пытаются максимально использовать всю имеющуюся информацию об объекте моделирования, как количественную, так и качественную. Грубо говоря, процесс построения имитационной модели можно представить следующим образом. Мы записываем в любом доступном для компьютера формализованном виде (в виде уравнений, графиков, логических соотношений, вероятностных законов) все, что знаем о системе, а потом проигрываем на компьютере варианты того, что может дать совокупность этих знаний при тех или иных значениях внешних и внутренних параметров системы. Если вопросы, задаваемые нами модели, относятся не к выяснению фундаментальных законов и причин, определяющих динамику реальной системы, а к бихевиористскому (поведенческому) анализу системы, как правило, выполняемому в практических целях, имитационная модель оказывается исключительно полезной.

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

Имитационные модели (simulation) Основные этапы построения имитационной модели следующие. Формулируются основные вопросы о поведении сложной системы, ответы на которые мы хотели бы получить. В соответствии с задачами моделирования задается вектор состояния системы. Вводится системное время, моделирующее ход времени в реальной системе. Временной шаг модели также определяется целями моделирования. Производится декомпозиция системы на отдельные блоки, связанные друг с другом, но обладающие относительной независимостью. Для каждого блока определяют, какие компоненты вектора состояния должны преобразовываться в процессе его функционирования. Формулируют законы и гипотезы, определяющие поведение отдельных блоков и связь этих блоков друг с другом. Для каждого блока множество законов функционирования дополняется множеством логических операторов, формализующих опыт наблюдения за динамикой процессов в. системе. При необходимости вводится внутреннее системное время данного блока модели, позволяющее моделировать более быстрые или более медленные процессы. Если в блоке используются случайные параметры, задаются правила отыскания на каждом шаге некоторых их реализаций. Разрабатываются программы, соответствующие отдельным блокам. Каждый блок верифицируется по фактическим данным, и при этом его информационные связи с другими блоками замораживаются. Обычно последовательность действий при верификации блоков такова: часть имеющейся информации используется для оценки параметров модели, а затем по оставшейся части информации сравнением расчетных данных с фактическими проверяется адекватность модели.

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

Имитационные модели (simulation) Таким образом, основные задачи имитационного моделирования: 1. проверка гипотез о взаимодействии отдельных элементов и подсистем; 2. прогноз поведения при изменении внутренних характеристик и внешних условий; 3. оптимизация управления. Ясно, что разработка имитационной модели сложной системы и работа с этой моделью требуют усилий целого коллектива специалистов, как в области машинной математики, так и в предметной области. К настоящему времени в литературе имеются тысячи имитационных моделей биологических систем самого разного уровня, многие модели представлены в ИНТЕРНЕТ.

Примеры Молекулярная динамика. Основные принципы построения моделей и результаты молекулярной динамики представлены на сайте Информационная система Российская биофизика. Биофизическое образование. На протяжении всей истории западной науки стоял вопрос о том, можно ли, зная координаты всех атомов и законы их взаимодействия, описать все процессы, происходящие во Вселенной. Вопрос не нашел своего однозначного ответа. Квантовая механика утвердила понятие неопределенности на микроуровне. В лекциях мы увидим, что существование квазистохастических типов поведения в детерминированных системах делает практически невозможным предсказание поведения некоторых детерминированных систем и на макроуровне. Следствием первого вопроса является второй: вопрос «сводимости». Можно ли, зная законы физики, т.е. законы движения всех атомов, входящих в состав биологических систем, и законы их взаимодействия, описать поведение живых систем. В принципе, на этот вопрос можно ответить с помощью имитационной модели, в которую заложены координаты и скорости движения всех атомов какой-либо живой системы и законы их взаимодействия. Для любой живой системы такая модель должна содержать огромное количество переменных и параметров и практически неосуществима, но попытки моделировать с помощью такого подхода функционирование элементов живых систем биомакромолекул делаются, начиная с 70-х годов.

Примеры 2. Модели систем организма. В настоящее время имеются имитационные модели многих систем организма – сердца, желудочно-кишечного тракта, почек, печени, мозга, и других. 3. Модели продукционного процесса растений. Имитационные модели продукционного процесса растений (агробиоценозов) для разных культур являются одними из первых имитационных моделей. Практическая задача моделирования выбор оптимальной стратегии проведения сельскохозяйственных мероприятий: орошения, полива, внесения удобрений с целью получения максимального урожая. Существует большое число моделей разных культур, как упрощенных, предназначенных для решения конкретных вопросов управления, так и очень подробных, используемых в основном для исследовательских целей. Подробные модели имеют иерархическую блочную структуру. Среди биотических процессов выделяют блок фотосинтеза, блок корневого питания, блок роста и развития, блок почвенной микрофлоры, блок развития болезней сельскохозяйственной культуры и другие. Рассматриваются также геофизические процессы: формирование теплового и водного режима, концентрации и передвижения биогенных и токсических солей, концентрации СО 2 в посеве и других. Методику работы с такими сложными моделями мы рассмотрели выше. Более подробное описание моделей продукционного процесса растений можно найти в книгах: 1.Бондаренко Н.Ф. «Моделирование продуктивности агроэкосистем». Л., 1982; 2.Заславский Б.Г., Полуэктов Р.А. Управление экологическими системами. М Торнли Дж. Математические модели в физиологии растений. Киев, Франс Дж., Торнли Дж. «Математические модели в сельском хозяйстве», М., 1987; 5.Vries de P. Simulation of plant growth and crop production/ Wageningen, Wit C.T. Simulation of assimilation, respiration, and transpiration of crops, Wageningen, 1978 Kниги 3-6 имели несколько более поздних переизданий на Западе.

Примеры 4. Модели водных экосистем. Водная среда гораздо более однородна, чем сухопутные биогеоценозы, и имитационные модели водных систем успешно создаются начиная с 70-х годов 20 века. Описание обменных процессов в водной среде включает описание усвоения азота, фосфора и других биогенных элементов, рост фито- и зоопланктона и детрита. При этом важно учитывать гидробиологические процессы в рассматриваемых водоемах, которые, как правило, являются неоднородными и при моделировании разбиваются на ряд компартментов. С помощью имитационного моделирования решались вопросы выработки стратегии борьбы с эфтрификацией закрытых водоемов, в частности, одного из Великих Американских озер – Озера Эри. Много имитационных моделей посвящено разработке оптимальной стратегии вылова рыбы. Пионерскими в этой области были книги: Меншуткин В.В. Математическое моделирование популяций и сообществ водных животных, Л., 1971 Jorgensen S.E. Lake management. Oxford, 1980 Экологические системы. Адаптивная оценка и управление. (под ред Э.Холлинга), М., 1981 Горстко А.Б., Домбровский Ю.А., Сурков Ф.А. Методы управления эколого- эконоическими. М., 1985 Основные идеи и результаты по моделированию водных систем, так же как и по моделированию продукционного процесса растений изложены в учебном пособии Г.Ю.Ризниченко, А.Б.Рубин «Математические модели биологических продукционных процессов». М., Готовится к печати дополненное и переработанное издание

Модели глобальной динамики сыграли особую роль в становлении имитационного моделирования. Именно для этих моделей был разработан формализм представления системы в виде узлов и потоков между ними, который затем в разных видах использовался практически во всех моделях сложных систем. Первая глобальная модель была создана Д. Форрестером и Д. Медоузом с соавторами по заказу Римского клуба в 60 годы 20 века. [J.W.Forrester. World dynamics. Cambridge:Wright-Allen Press, 1972] Полученные с ее помощью результаты были опубликованы в знаменитой переведенной на 35 языков книге «Пределы роста», и впервые послужили предостережением человечеству в том, что Земля – ограниченная система, безудержный прогресс ведет к истощению ее ресурсов, и человечество ждет глобальный экологический кризис. [Donella H.Meadows et.al. The Limits of the Growth. N.-Y. Universe Books. 1972, перевод на русский язык 1991 г.]. Современное состояние проблемы описано в книге Д.Х.Медоуз, Д.Л.Медоуз, Й.Рандерс «За пределами роста» М., Прогресс (Donella H.Meadows et.al Beyond the Limits, (Confronting global collapse. Envisioning a sustainable future.1992) Вторая знаменитая глобальная модель – модель ядерной зимы, была создана под руководством Н.Н. Моисеева в России. Ее результаты наглядно показали, что глобальная ядерная война приведет к уничтожению как побежденных, так и победителей, так как после нее небо над всей Землей закроется тучами и настанет ядерная зима на период в несколько десятков лет. Поэтому победа в такой войне будет быссмысленной. В настоящее время активно разрабатываются глобальные модели, позволяющие рассчитать «парниковый эффект» и другие процессы, протекающие в глобальном масштабе. Ясно, что разработка имитационной модели сложной системы и работа с этой моделью требуют усилий целого коллектива специалистов как в области машинной математики, так и в предметной области. Подробное изучение методологии имитационного моделирования не входит в задачу нашего курса, мы будем заниматься более общими вопросами.

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

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

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

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

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

ЛЕКЦИЯ 2 МОДЕЛИ БИОЛОГИЧЕСКИХ СИСТЕМ, ОПИСЫВАЕМЫЕ ОДНИМ ДИФФЕРЕНЦИАЛЬНЫМ УРАВНЕНИЕМ ПЕРВОГО ПОРЯДКА ЛЕКЦИЯ 2 МОДЕЛИ БИОЛОГИЧЕСКИХ СИСТЕМ, ОПИСЫВАЕМЫЕ ОДНИМ ДИФФЕРЕНЦИАЛЬНЫМ УРАВНЕНИЕМ ПЕРВОГО ПОРЯДКА

Изучение математических моделей биологических систем начнем с систем первого порядка, которым соответствует одно дифференциальное уравнение первого порядка: Если система автономная, то правая часть уравнений не зависит явно от времени и уравнение имеет вид: (2.1) Состояние таких систем в каждый момент времени характеризуется одной единственной величиной – значением переменной x в данный момент времени t. Рассмотрим плоскость t, x. Решениями уравнения (2.1): x(t) являются кривые на плоскости t, x, называемые интегральными кривыми (рис. 2.1) Пусть заданы начальные условия при t=0 или, иначе, пусть на плоскости t, x задана точка с координатами. Если для уравнения (2.1) выполнены условия теоремы Коши, то имеется единственное решение уравнения (2.1), удовлетворяющее этим начальным условиям, и через точку проходит одна единственная интегральная кривая x(t).

Интегральные кривые уравнения (2.1) не могут пересекаться. Решения уравнения (2.1) не могут быть периодическими, они монотонны. Поведение интегральных кривых на плоскости t, x можно установить, не решая в явном виде дифференциального уравнения (2.1), если известен характер движения изображающей точки на фазовой прямой. Рассмотрим плоскость t, x, причем фазовую прямую совместим с осью x. Построим на плоскости t, x точку с абсциссой t и с ординатой, равной смещению изображающей точки по оси x в данный момент времени t. С течением времени в соответствии с уравнением (2.1) изображающая точка будет двигаться по фазовой прямой (рис. 2.2), а на плоскости t, x описывать некую кривую. Это будет интегральная кривая уравнения (2.1). Решения одного автономного дифференциального уравнения либо уходят в бесконечность (чего не бывает в реальных системах), либо асимптотически приближаются к стационарному состоянию.

Стационарное состояние (точка покоя, особая точка, состояние равновесия) В стационарном состоянии значения переменных в системе не меняются со временем. На языке дифференциальных уравнений это означает: (2.2) Если левая часть уравнения равна нулю, значит равна нулю и его правая часть: f(x) = 0. (2.3) Корни алгебраического уравнения (2.3): суть стационарные состояния дифференциального уравнения (2.1). На плоскости (t, x) прямые – асимптоты, к которым приближаются интегральные кривые. На фазовой прямой (рис. 2.2) стационарное состояние – точка, к которой стремится величина x. Реальные биологические системы испытывают многочисленные флуктуации, переменные при малых отклонениях возвращаются к своим стационарным значениям. Поэтому при построении модели важно знать, устойчивы ли стационарные состояния модели.

Каждый имеет интуитивное представление об устойчивости. На рис в обоих положениях (а и б) шарик находится в равновесии, т.к. сумма сил, действующих на него, равна нулю. Попытайтесь ответить на вопрос: «Какое из этих состояний равновесия устойчиво?» Скорее всего, Вы дали правильный ответ. Сказать, как Вы догадались? Вы дали шарику малое отклонение от состояния равновесия. В случае (а) шарик вернулся. В случае (б) покинул состояние равновесия навсегда. Устойчивое состояние равновесия можно определить так: если при достаточно малом отклонении от положения равновесия система никогда не уйдет далеко от особой точки, то особая точка будет устойчивым состоянием равновесия, что соответствует устойчивому режиму функционирования системы. Строгое математическое определение устойчивости состояния равновесия уравнения dx/dt = f(x) выглядит следующим образом :

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

Стационарное состояние называется асимптотически устойчивым, если малые отклонения от него со временем затухают. Пример шарик в ямке в вязкой среде. Стационарное состояние называется неустойчивым, если малые отклонения со временем увеличиваются. Пример: шарик на горке. Устойчивое стационарное состояние представляет собой простейший тип аттрактора. Аттрактором называется множество, к которому стремится изображающая точка системы с течением времени (притягивающее множество). В нашем курсе мы рассмотрим следующие типы аттракторов: ·устойчивая точка покоя; ·предельный цикл режим колебаний с постоянными периодом и амплитудой (начиная с размерности системы 2); ·Области с квазистохастическим поведением траекторий в области аттрактора, например, «странный аттрактор» (начиная с размерности 3).

Аналитический метод исследования устойчивости стационарного состояния (метод Ляпунова). Линеаризация системы в окрестности стационарного состояния. Метод Ляпунова приложим к широкому классу систем различной размерности, точечным системам, которые описываются обыкновенными дифференциальными уравнениями, и распределенным системам, описываемым уравнениями в частных производных, непрерывным и дискретным. Рассмотрим метод линеаризации Ляпунова для одного автономного дифференциального уравнения первого порядка. Пусть стационарное решение уравнения (2.1) Пусть система, первоначально находившаяся в стационарном состоянии, отклонилась от него и перешла в близкую точку с координатой:, причем. Перейдем в уравнении (2.1) от переменной x к переменной, т.е. новой переменной будет отклонение системы от стационарного состояния. Получим:..

Правую часть разложим в ряд Тейлора в точке : Или где Отбросим члены порядка 2 и выше. Останется линейное уравнение: (2.4) которое носит название линеаризованного уравнения или уравнения первого приближения. Интеграл этого уравнения для находится сразу:, (2.5) где, с произвольная постоянная.

В случае одного уравнения вопрос об устойчивости состояния равновесия нетрудно решить, рассматривая график функции f(x). По определению в стационарной точке правая часть уравнения (2.1) функция f(x) обращается в нуль. Здесь возможны три случая (рис. 2.4 а, б, в). Рис Определение устойчивости стационарного состояния по графику функции f(x) a – стационарное состояние устойчиво; б, в стационарное состояние неустойчиво

2. Вблизи состояния равновесия функция f(x) меняет знак с минуса на плюс при возрастании x (рис. 2.4 б). Проведите рассуждения, аналогичные случаю 1. Поместите изображающую точку в область. Теперь в область. В обоих случаях изображающая точка удаляется от состояния равновесия. Стационарное состояние неустойчиво. 3. Вблизи состояния равновесия функции f(x) не меняет знак (рис 2.4 в). Поскольку, это означает, что изображающая точка, помещенная достаточно близко к состоянию равновесия с одной стороны, будет приближаться к нему, помещенная с другой стороны – удаляться. Вопрос. Является ли состояние равновесия в случае 3 устойчивым? Ответ. Нет. По определению устойчивости.

Уравнение Ферхюльста Рассмотрим еще один пример, который относится к классическим моделям математической экологии. Логистическое уравнение было предложено Ферхюльстом в 1838 г. Оно имеет вид: (2.15) Это уравнение обладает двумя важными свойствами. При малых х численность х возрастает, при больших – приближается к определенному пределу К. Уравнение (2.15) можно решить аналитически. Ход решения следующий. Произведем разделение переменных: (2.16) Представим левую часть в виде суммы и проинтегрируем Переходя от логарифмов к переменным, получим: (2.17) Здесь С – произвольная постоянная, которая определяется начальным значением численности x 0 : Подставим это значение С в формулу (2.17): Отсюда получим решение – зависимость численности от времени: (2.18) График функции (2.18) при разных начальных значениях численности популяции представлен на рис. 2.7.

Если начальное значение х 0 К, численность со временем убывает. В приведенных примерах в правой части уравнений стоят полиномы первой и второй степени. Если в правой части более сложная нелинейная функция, алгебраическое уравнение для стационарных значений может иметь несколько корней. Какое из этих решений реализуется в этом случае, будет зависеть от начальных условий. В дальнейшем мы, как правило, не будем искать аналитическое решение для наших моделей. Для более сложных нелинейных уравнений это и невозможно. Однако важные заключения относительно свойств моделей можно сделать и на основании качественного их исследования, в первую очередь путем исследования устойчивости стационарных состояний и типов поведения системы вблизи этих состояний. При этом следует иметь в виду, что с помощью одного автономного дифференциального уравнения могут быть описаны только монотонные изменения переменной, и, следовательно, ни периодические, ни хаотические процессы не могут быть описаны. Для описания более сложного поведения необходимо либо переходить к системам большей размерности (2, 3 порядка и выше), либо вводить время в явном виде в правую часть уравнения. В Лекции 3 мы увидим, что дискретные уравнения и уравнения с запаздыванием могут описать и колебания, и динамический хаос. ПРИМЕР

ЛЕКЦИЯ 3 МОДЕЛИ РОСТА ПОПУЛЯЦИЙ

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

Модель предложена Мальтусом в 1798 г. в его классическом труде «О законе роста народонаселения». Томас Роберт Мальтус ( ) – известный английский демограф и экономист, обратил внимание на тот факт, что численность популяции растет по экспоненте (в геометрической прогрессии), в то время как производствопитания растет со временем линейно (в арифметической прогрессии), из чего сделал справедливый вывод, что рано или поздно экспонента обязательно «обгонит» линейную функцию, и наступит голод. На основании этих выводов Мальтус говорит о необходимости ввести ограничения на рождаемость, в особенности для беднейших слоев общества. «Экономический пессимизм», следующий из прогнозов предложенной им модели, в основу которой положен анализ эмпирических данных, Мальтус противопоставлял модным в начале IXX века оптимистическим идеям гуманистов: Жана - Жака Руссо, Уильяма Годвина и других, предсказывающих человечеству грядущее счастье и процветание. Можно говорить о том, что Мальтус был первым ученым «алармистом», который на основании результатов моделирования «бил тревогу» и предупреждал человечество об опасности следования развитию по используемым ранее сценариям прогресса. Во второй половине XX века такую «алармистскую» роль сыграли работы Римского клуба, и в первую очередь «модель глобального роста» Дж. Форрестера (см. Лекция 1).

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

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

Логистическое уравнение обладает двумя важными свойствами. При малых зачениях х численность возрастает экспоненциально (как в уравнении 3.2) при больших – приближается к определенному пределу К. Эта величина, называемая емкостью экологической ниши популяции, определяется ограниченностью пищевых ресурсов, мест для гнездования, многими другими факторами, которые могут быть различными для разных видов. Таким образом, емкость экологической ниши представляет собой системный фактор, который определяет ограниченность роста популяции в данном ареале обитания. Уравнение (3.3) можно также переписать в виде: (3.4) Здесь - коэффициент внутривидовой конкуренции (за пищевой ресурс, убежища и т. п.) Аналитическое решение уравнения (3.3) мы получили в лекции 2: (3.5)

Формула (3.5) описывает кинетическую кривую, то есть зависимость численности популяции от времени. Ход кинетических кривых для разных начальных условий представлен на рис В случае, если начальная численность меньше половины максимальной, кривая x(t) имеет точку перегиба с координатами Ордината точки перегиба представляет собой половину максимальной численности, а абсцисса зависит как от емкости популяции K, так и от константы собственной скорости роста r чем выше генетические возможности популяции, тем скорее наступает перегиб на кривой численности. Примеры экспериментально наблюдаемой динамики популяций, развивающихся по логистическому закону, приведены на рис. 3.4 а, б. На рис. 3.4 а сплошной линией представлен график функции (3.5), при малых начальных численностях он имеет S- образный характер. При изучении более сложных систем мы не будем искать решение для x (t) в явном виде, а ограничимся исследованием устойчивости их стационарных состояний. Проведем такое исследование и для логистического уравнения. Легко видеть, что уравнение стационарных состояний f ( ) = 0 в данном случае имеет два корня: ПРИМЕР

Посмотрим, будут ли эти корни устойчивыми. Для этого вначале воспользуемся аналитическим методом Ляпунова. Введем новую переменную x, обозначающую отклонение переменной х от ее стационарного значения: x =. Запишем линеаризованное уравнение для (3.4): dx /dt = ax, где а =. Напомним, что знак величины определяет устойчивость соответствующей особой точки : (3.5) Подставив в выражение (3.5) значение первого корня мы получим. Будем считать, что величина r – коэффициент естественной скорости роста популяции положительная; – неустойчивая особая точка. Если же мы подставим в выражение (3.5) то получим – отрицательную величину. Это дает нам право утверждать, что стационарное решение уравнения соответствует устойчивому стационарному режиму существования популяции в ограниченной среде.

Проведем теперь исследование устойчивости стационарных решений этого уравнения, исходя из графика функции правой части. На рисунке 3.5 видно, что при переходе от отрицательных к положительным значениям х в точке = 0 функция f (x) меняет знак с минуса на плюс, т.е. особая точка неустойчива. Наоборот, в точке = К имеет место изменение знака f(x) с ростом x с плюса на минус, следовательно, эта особая точка устойчивая. Несмотря на схематичность положенных в ее основу представлений, логистическая кривая оказалась очень хорошим приближением для описания кривых роста численности многих популяций. В природе внутривидовая конкуренция не удерживает естественные популяции на строго неизменном уровне, но действует в широком диапазоне начальных значений плотности и приводит их к гораздо более узкому диапазону конечных значений, определяя, таким образом, тенденцию к поддержанию плотности в определенных пределах.

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

Рис Кривые пополнения, соответствующие в модели кривым зависимости скорости роста от численности популяции. а, б. - схема расчета кривых пополнения, в, г, д - кривые пополнения для разных видов: в - численность обыкновенного фазана на о. Протекши - Айленд после его интродукции в 1937 г. (Einarsen, 1945); г - экспериментальная популяция плодовой мушки Drosophyla melanogaster (Pearl, 1927); д - оценка численности арктического финвала (Allen, 1972)

Модель популяции с наименьшей критической численностью В рассмотренных моделях прирост численности (биомассы) популяции представлен линейным членом, пропорциональным численности. Строго говоря, это соответствует лишь тем популяциям, размножение которых происходит путем самооплодотворения (микроорганизмы). Если же в основе размножения лежит скрещивание, предполагающее встречи между особями разных полов одного и того же вида, то прирост будет тем выше, чем больше количество встреч между особями, а последнее пропорционально второй степени х. Таким образом, для разнополой популяции в условиях неограниченных ресурсов можно записать (3.6) Уравнение (3.6) хорошо описывает тот факт, что при низких плотностях популяций скорость размножения резко падает, так как вероятность встречи двух особей разных полов уменьшается при понижении плотности популяции пропорционально квадрату плотности. Однако при больших плотностях популяций скорость размножения лимитирует уже не число встреч особей противоположного пола, а число самок в популяции, формула, учитывающая эти оба эффекта, имеет вид (3.7) Графики численности в зависимости от времени (рис. 3.7 а) и скорости размножения как функции численности (рис ) для уравнения (3.7.) представлены на рис В действительности плотность популяции не должна опускаться ниже некоторой критической величины. При падении плотности популяции ниже критической среднее время, в течение которого может состояться оплодотворение, становится больше времени жизни отдельной особи, точнее времени, в течение которого особь способна к размножению. В этом случае популяция вымирает.

Рис а – график зависимости численности от времени и б – скорости размножения как функции численности для уравнения 3.7

Этот эффект может быть учтен, если в формулу (3.7) ввести член, пропорциональный численности и описывающий смертность. Зависимость скорости роста популяции от ее численности при этом примет вид (3.8) Уравнение (3.8) имеет два стационарных решения: =0 и = d ( – d ) =L. Соответствующие графики x(t) и f (x) даны на рис. 3.8, а, б. Из графика 3.8,б видно, что решение = 0 устойчивое, а =L– неустойчивое. При начальных численностях x нач L в соответствии с уравнением (3.8) популяция неограниченно размножается.

Рис.3.8. Зависимость численности популяции от времени (а) и скорости роста от численности (б) для уравнения (3.8). Штриховкой обозначена область вырождения популяции.

Величина нижней критической плотности L различна для разных видов. Наблюдения биологов показали, что это всего лишь одна пара особей на тысячу квадратных километров в случае ондатр и сотни тысяч особей для американского странствующего голубя. Заранее трудно было предугадать, что столь многочисленный вид перешел через критическую границу своей численности и обречен на вырождение. Однако это произошло, несмотря на все усилия по охране этих птиц. Для голубых китов критическая граница общей численности оказалась равной десяткам – сотням. Хищническое истребление этих гигантских животных привело к тому, что их осталось слишком мало в Мировом океане. И хотя охота на них запрещена, надежд на восстановление популяции голубых китов практически нет. Наиболее общая формула, учитывающая как нижнюю границу численности, так и внутривидовую конкуренцию, имеет вид (3.9) Зависимости численности от времени и скорости прироста от численности представлены на рис.3.9 а, б. =0 и =К –устойчивые стационарные состояния, =L – неустойчивое, разделяющее области притяжения устойчивых состояний равновесия. Величины L и К различны для разных популяций и могут быть определены только из наблюдений и экспериментов. Ясно, что их определение представляет значительные трудности. Кривые 1, 2, 3, 4 на рис. 3.9 а соответствуют различным начальным значениям численности популяции.

При любых промыслах особый интерес представляет величина нижней критической границы, при переходе через которую популяция уже не сможет восстановиться. Модель позволяет дать некий методический рецепт определения не самой критической границы, но степени близости к ней численности вида. Обратимся к рис. 3.9 а. Пусть численность вида в начальный момент времени была близка к максимально возможной. При t=0 происходит одноразовое выбивание популяции. Если численность осталась значительно больше критической, восстановление происходит сначала быстро, а затем с монотонным замедлением (кривая 1). Если же оставшаяся популяция близка к критической точке, восстановление происходит сначала очень медленно, численность популяции надолго «застревает» вблизи критической точки, а затем уже, «набрав силы», более быстро приближается к стационарному уровню (кривая 3). Кривая 2 представляет промежуточный случай. Таким образом, наблюдая реакцию системы на возмущение, можно предсказать приближение ее к опасным границам.

ДИСКРЕТНЫЕ МОДЕЛИ ПОПУЛЯЦИЙ. Численность популяции не изменяется непрерывно, а представляет собой дискретную величину, что соответствует экспериментальным данным по переписи реальных популяций. Если предположить, что численность N зависит от численностей в некоторые предшествующие моменты времени, то для описания динамики численности популяций можно применить аппарат разностных уравнений (отображений). Если при этом внешние и внутренние факторы, определяющие развитие популяции, остаются во времени неизменными, то численность популяции в момент времени t может быть описана при помощи разностного уравнения в виде: (3.10) Здесь функция F зависит от численности популяции в k предшествующие моменты времени.

Особенно просто выглядит разностное уравнение в случае, когда численность каждого следующего поколения в популяции N t+1 зависит лишь от предыдущего поколения N t. Это справедливо для многих видов насекомых. Их взрослые особи живут непродолжительное время, достаточное для откладывания яиц, и к моменту появления на свет нового поколения (на стадии взрослой особи) предшествующее поколение прекращает свое существование. То же имеет место для некоторых видов зоопланктона, рыб, птиц. Про эти виды можно сказать, что поколения в них не перекрываются во времени, и уравнение (3.10) может быть записано в виде: (3.11) Опыт показывает, что в таких системах при малых численностях N растет от одной генерации к другой, а при высоких – падает. Это свойство в экономике проявляется как закон «бумов» и «спадов».

В таких случаях функция N – одноэкстремальная функция (рис. 3.10) В зависимости от крутизны графика F (N t ) в системе могут возникать различные режимы: монотонное и колебательное приближение к равновесию, колебательные изменения – циклы разной длины и квазистохастическое поведение хаос. Модели вида (3.11) являются простейшими детерминированными объектами, демонстрирующими хаотическое поведение. Если поколения перекрываются – появляются дискретные возрастные группы (стадии развития). Для их описания обычно используются тоже дискретные, но более сложные матричные модели.

Вспомним логистическое уравнение, которое описывает развитие популяции в непрерывном времени. Заменим dN/dt на DN/ t. Здесь Получим: (3.12) или N t+1 =N t · f(N t ). Однако, уравнение (3.12) биологически некорректно. Если в какой-то момент времени значение численности становится больше определенной величины: N t >K(1+r)/r, то уравнение (3.12) дает отрицательное значение N t+1. Это связано с видом функции f (N t ) (рис а). Этот недостаток отсутствует у непрерывного логистического уравнения. Чтобы исправить положение в дискретном уравнении, в качестве f (N t ) следует взять функцию, асимптотически стремящуюся к нулю при N t. Вид такой функции изображен на рис б.

Экспоненциальная форма такой зависимости была предложена Мораном (1950) для численности насекомых и Рикером (1952) для рыбных популяций: (3.13) Ход решения можно наглядно продемонстрировать графически с помощью диаграммы и лестницы Ламерея. Точка пересечения биссектрисы первого координатного угла N t+1 =N t и функции F(N t ) определяет равновесное состояние системы, аналогичное стационарному состоянию дифференциального уравнения. На рис б показан способ нахождения значений N t в последовательные моменты времени. Пусть в начальный момент времени N=N 0.. F(N 0 )=N 1 задает значение численности в последующий момент времени t=1. Величина N 1 в свою очередь определяет значение F(N 1 )=N 2. И так далее. На рис б изображен случай, когда траектория сходится к равновесному состоянию, совершая затухающие колебания.

Введем для разностных уравнений (отображений) некоторые понятия, соответствующие основным понятиям теории дифференциальных уравнений. Решение – (траектория) – любая последовательность значений {N t }, (t=0,1...), удовлетворяющая данному разностному уравнению при любом t. Разным начальным значениям соответствуют разные решения. Равновесием называется решение вида: N t =const=N *, удовлетворяющее соотношению N * =F(N * ) Как и в случае дифференциальных уравнений для исследования устойчивости применим линейный анализ. Положим: Линеаризуем уравнение (3.11), разлагая F в ряд по степеням x t и отбрасывая члены порядка и выше. Получим: Из условий сходимости геометрической прогрессии следует, что x t 0 при t, если (3.14) и x t при t, если (3.15) Случаи dF/dN =1 и dF/dN =0 требуют дополнительного исследования членов более высокого порядка в разложении.

Обозначим. В случае устойчивого равновесия: при 0

Для уравнения (3.13) равновесие находится из выражения Единственное равновесное значение N * = K>0 существует при любом r. Равновесие устойчиво, если 0

Рис Примеры хаотических решений для уравнения 3.13 Рис Циклическое поведение решения уравнения 3.13 а - двухточечный цикл, б - четырехточечный цикл.

Иногда в качестве дискретного аналога используется более простая формула, в которой численность выражена в безразмерном виде. При N t > 1 численность принимается равной нулю. (3.16) На рис показаны диаграммы Ламерея и динамика численности во времени для разных значений параметра a в уравнении Рис Диаграммы Ламерея и динамика численности во времени для разных значений параметра a в уравнении a - монотонный рост; б – двухточечный цикл; в – квазистохастическое поведение

Матричные модели популяций Детализация возрастной структуры популяций приводит к классу моделей, впервые предложенных Лесли, (1945, 1948). Пусть ресурсы питания не ограничены. Размножение происходит в определенные моменты времени Пусть популяция содержит n возрастных групп. Тогда в каждый фиксированный момент времени популяцию можно охарактеризовать вектор-столбцом(3.18) Вектор X(t 1 ), характеризующий популяцию в следующий момент времени, например, через год, связан с вектором X(t 0 ) через матрицу перехода L: (3.19) Установим вид этой матрицы. Из всех возрастных групп выделим те, которые производят потомство. Пусть их номера будут k, k+1,..., k+p. Предположим, что за единичный промежуток времени особи i-й группы переходят в группу i+1, от групп k, k+1,..., k+p появляется потомство, а часть особей от каждой группы погибает. Потомство, которое появилось за единицу времени от всех групп, поступает в группу 1. (3.20)

Вторая компонента получается с учетом двух процессов. Первый – переход особей, находившихся в момент в первой группе, во вторую. Второй процесс возможная гибель части из этих особей. Поэтому вторая компонента x 2 (t 1 ) равна не всей численности x 1 (t 0 ), а только некоторой ее части, Аналогично получаются третья компонента и все остальные. Предположим, что все особи, находившиеся в момент t 0 в последней возрастной группе к моменту t 1 погибнут. Поэтому последняя компонента вектора X (t 1 ) составляется лишь из тех особей, которые перешли из предыдущей возрастной группы.

Коэффициенты для каждой группы имеют следующий смысл: - коэффициент рождаемости, - коэффициент выживания. Вектор численностей возрастных групп в момент времени t 1 представим в виде: (3.21) Вектор X(t 1 ) получается умножением вектора X(t 0 ) на матрицу (3.22)

По диагонали матрицы стоят нули, под диагональными элементами коэффициенты выживания, на первой строке стоят члены, характеризующие число особей, родившихся от соответствующих групп. Все остальные элементы матрицы равны нулю. (3.23) Таким образом, зная структуру матрицы L и начальное состояние популяции – вектор-столбец X(t 0 ), – можно прогнозировать состояние популяции в любой наперед заданный момент времени. Главное собственное число матрицы L дает скорость, с которой размножается популяция, когда ее возрастная структура стабилизировалась. Пример популяции из трех возрастных групп (Уильямсон, 1967) Пусть возрастная динамика популяции характеризуется матрицей:

Такая запись означает, что исходная популяция состоит из одной самки старшего возраста (вектор столбец в правой части уравнения). Каждое животное старшего возраста, прежде чем умереть, успевает произвести в среднем 12 потомков, каждое животное среднего возраста, прежде чем умереть или перейти в следующий возрастной класс (вероятности этих событий одинаковы) производит в среднем 9 потомков. Молодые животные не производят потомства и с вероятностью 1/3 попадают в среднюю возрастную группу. По прошествии одного временного интервала в популяции будет уже 12 самок младшего возраста: Далее процедуру следует повторять на каждом шаге. Из графика видно, что до некоторого момента времени ( t 10 ), наблюдаются колебания численности, после чего количество самок всех трех возрастов экспоненциально возрастает, причем со–отношение между ними остается постоянным. Главное собственное число 1 при этом равно 2, т.е. размер популяции за каждый временной шаг удваивается.

Наклон графика равен ln 1 – собственной скорости естественного прироста. Соответствующий главному собственному числу собственный вектор отражает устойчивую структуру популяции и в нашем случае равен Этот пример страдает тем же недостатком, что и модель Мальтуса экспоненциального роста: мы допускаем, что популяция может неограниченно расти. Более реалистическая модель должна учитывать, что все элементы матрицы L являются некоторыми функциями размера популяции. Модели с применениями матриц Лесли для крупных возрастных групп могут дать описание колебательных изменений численности популяции. Пример такой модели описание динамики популgfяции овсеца Шелли мелкодерновинного злака северных луговых степей (Розенберг, 1984). Модель позволила описать наблюдаемые в природе явления - старение овсеца и колебания распределений по возрастному спектру в течение ряда лет (рис. 3.19).

Уравнения с запаздыванием В реальных системах всегда имеется некоторое запаздывание в регуляции численности, вызванное несколькими причинами. Развитие любой взрослой особи из оплодотворенного яйца требует определенного времени. Поэтому если какое-нибудь изменение внешних факторов, например, увеличение ресурсов, вызовет повышение продуктивности взрослых особей, то соответствующее изменение численности произойдет лишь по прошествии времени T. Это означает, что уравнение dx/dt=f(x), где x - число взрослых особей в момент времени t следует заменить уравнением dx/dt=f{x(t-T)}.

Логистическое уравнение с запаздыванием может быть записано в виде: (3.24) Более точное уравнение, учитывающее распределение времени запаздывания имеет вид: (3.25) Типичный вид весовой функции (t) изображен на рисунке Такого типа уравнения могут иметь колебательные решения. Это легко проверить для простого линейного уравнения: (3.26) Уравнение (3.26) имеет периодическое решение N(t) = Acos( t/2T) в широком диапазоне значений скоростей роста r и времени запаздывания T.

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

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

Вероятностное описание процессов размножения и гибели. В качестве простейшего примера рассмотрим вероятностное описание процесса роста популяции с учетом только размножения. При детерминистическом подходе мы считали, что существует определенная скорость размножения такая, что численность популяции п за время dt увеличивается на dn = n dt. Это приводит к экспоненциальному закону n = a e t. Здесь а – численность популяции в начальный момент времени. Подойдем к процессу размножения с вероятностной точки зрения. Пусть вероятность появления одного потомка у данной особи в интервале времени dt равна dt. Тогда вероятность появления одной новой особи в целой популяции за время dt равна n dt. Обозначим через p n (t) вероятность того, что в момент t в популяции имеется ровно п особей.

Предположим, что в каждый момент времени может произойти только одно событие, а именно за время dt численность популяции может либо увеличиться на 1, либо остаться неизменной. Размер популяции в момент t можно связать с размером популяции в момент t + dt с помощью следующих рассуждений. Если число особей в момент t + dt равно n, это означает, что либо в момент t их было п – 1 и за время dt появилась еще одна, либо в момент t было п особей и за время dt это число не изменилось. Складывая вероятности, получим соотношение p n (t+dt) = p n-1 (t) (n – 1) dt+ p n (t) (1- ) n dt,(3.27) Или (3.28) (3.29)

Системы дифференциально-разностных уравнений, аналогичных уравнениям (3.27), которые можно рассматривать как динамические уравнения для случайного процесса, обычно бывает трудно разрешить в общем виде. Однако в нашем примере это довольно просто. Проинтегрируем уравнение (3.29) с учетом того обстоятельства, что р а (0) = 1. Это дает В свою очередь этот результат подставляем в последующее уравнение, и весь процесс повторяется. После вычисления нескольких последовательных членов можно записать результат в общем виде: (3.30) Выражение (3.30) определяет распределение вероятностей для любого момента времени, заменяющее при вероятностном описании то единственное значение которое рассматривалось в детерминистической модели. Выражение (3.30) является частным случаем биномиального распределения. Его математическое ожидание и дисперсия записываются следующим образом: m(t) = a e t,(3.31) a 2 (t) = ae 2 t (1 – e - t ).(3.32)

Легко заметить, что математическое ожидание (3.31) совпадает с детерминистическим средним. Таким образом, при большом числе особей детерминистическое описание будет удовлетворительно заменять любую стохастическую модель, в которой основное внимание уделяется нахождению средних значений. Когда же число особей мало, например, когда начальный размер популяции составляет всего лишь несколько единиц, дисперсия, т.е. среднее квадратичное отклонение численности отдельно взятой популяции от математического ожидания, может быть довольно значительной. При этом при t коэффициент вариации величины п, равный / m, стремится к. При рассмотрении какой-либо определенной популяции мы наблюдаем только одно численное значение. График роста обнаруживает значительные колебания. Возникает вопрос: каким образом эти колебания связаны с распределением вероятностей?

Смысл выражения (3.30) состоит в том, что если имеется некоторое большое число популяций и в начальный момент времени t = 0 численность каждой из них равна а, то доля этих популяций, имеющих в момент t численность а, теоретически равна p n (t) с математическим ожиданием m(t) и дисперсией 2 (t). Кривая роста любой отдельно взятой популяции может значительно отклоняться от соответствующей кривой математического ожидания, так что последняя вместе с дисперсией служит показателем случайной флуктуационной изменчивости, характерной для данного процесса. Рассмотрим теперь более сложный процесс – размножение и гибель особей в популяции. Как и ранее, полагаем, что вероятность появления одного потомка у одной особи в интервале времени t равна t, поэтому для всей популяции вероятность увеличения ее численности на единицу равна n t. Допустим также, что вероятность гибели одной особи составляет t. Вероятность того, что размер популяции в момент t+ t составляет n особей, будет в таком случае представлять собой сумму вероятностей трех событий:

1) в момент времени t было n особей, и за время dt это число не изменилось; 2) в момент t было n – 1 особей, за время dt их количество увеличилось на единицу; 3) в момент времени t было n +1 особей, за время dt их количество уменьшилось на единицу. Выражение для p n (t + dt) принимает вид p n (t+dt)=p n-1 (t) (n–l)dt+p n (t)(1– ndt– indt)+p n+1 p(t) (n+l)dt, при n=l, 2,... Эта система уже не решается простым интегрированием, однако применение метода производящей функции (Бейли, 1970) позволяет найти общее решение: (3.33)

Таким образом, даже в случае простейшего стохастического процесса размножения и гибели общее выражение для p n (t) оказывается довольно сложным, и выразить его в явном виде, как правило, не удается. Математическое ожидание и дисперсия распределения (3.33) имеют вид (3.34) Как и в случае простого процесса размножения, математическое ожидание совпадает со значением численности в детерминистической модели, а выражение для дисперсии показывает, что имеет место значительная флуктуационная изменчивость.

Рассмотрим случай, когда размножение и гибель уравновешивают друг друга, т.е. когда =. Математическое ожидание и дисперсию находим из формул (3.34), полагая, что, и используя во втором выражении правило Лопиталя для раскрытия неопределенности вида 0/0, получаем m (е) = a, (t) = 2 a t. Первое выражение представляет собой очевидный результат, а именно средний размер популяции сохраняет свое начальное значение. Второе выражение показывает, что дисперсия размера популяции возрастает пропорционально длительности интервала времени, в течение которого протекает процесс. Детерминистическая модель в тех случаях, когда скорость размножения превышает скорость гибели, предсказывает устойчивое экспоненциальное увеличение размера популяции. Однако в вероятностной модели учитывается, что всегда существует определенная вероятность такого большого числа случаев гибели, при котором популяция полностью вымирает. Таким образом, вероятность вымирания является важной характеристикой вероятностной модели. Обозначим через p 0 (t) вероятность того, что в момент времени t не останется ни одной живой особи.

Приравняв правую часть к нулю, из уравнений (3.33) можно найти выражения для этой вероятности в явном виде: (3.35) (3.36) Вероятность того, что рано или поздно произойдет вымирание популяции, можно найти, полагая t. В пределе при t выражения (3.35) и (3.36) для случаев можно записать следующим образом: Следовательно, если скорость размножения не превышает скорости гибели, вымирание рано или поздно обязательно произойдет. Если же скорость размножения выше скорости гибели, то вероятность вымирания составляет (./ ) a.

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

Учет флуктуаций среды. Флуктуации условий среды могут приводить к изменению характера взаимодействий между отдельными особями, т.е. случайному изменению параметров модели. Для модели экспоненциального роста – это случайные изменения коэффициента естественного прироста, учет которых приводит к уравнению dx/dt = ( + y(t))x, m(t)=ae t где 2 – дисперсия y (t). Отсюда

Таким образом, с течением времени колебания численности популяции становятся более резкими; это значит, что детерминистическая система не имеет устойчивого стационарного состояния. Можно показать (Свирежев, Логофет, 1978), что при < 2 вероятность вырождения со временем увеличивается, стремясь в пределе к единице – популяция вероятностно неустойчива, т.е. достаточно длительное воздействие возмущений с большой вероятностью может привести к ее гибели. При > 2 вероятность вырождения уменьшается, и при t стремится к нулю – популяция в этом смысле устойчива. Из полученного результата следуют более жесткие ограничения на коэффициент естественного прироста, чем из детерминистической модели. В самом деле, в последней для невырождения популяции достаточно, чтобы среднее значение коэффициента было положительным, в то время как в стохастической модели этого недостаточно – нужно, чтобы > 2.

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

ЛЕКЦИЯ 4 МОДЕЛИ, ОПИСЫВАЕМЫЕ СИСТЕМАМИ ДВУХ АВТОНОМНЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ.

Наиболее интересные результаты по качественному моделированию свойств биологических систем получены на моделях из двух дифференциальных уравнений, которые допускают качественное исследование с помощью метода фазовой плоскости. Рассмотрим систему двух автономных обыкновенных дифференциальных уравнений общего вида (4.1) P(x,y), Q(x,y) - непрерывные функции, определенные в некоторой области G евклидовой плоскости (x,y декартовы координаты) и имеющие в этой области непрерывные производные порядка не ниже первого. Область G может быть как неограниченной, так и ограниченной. Если переменные x, y имеют конкретный биологический смысл (концентрации веществ, численности видов) чаще всего область G представляет собой положительный квадрант правой полуплоскости: 0 x

Концентрации веществ или численности видов также могут быть ограничены сверху объемом сосуда или площадью ареала обитания. Тогда область значений переменных имеет вид: 0 x < x 0, 0 y< y 0. Переменные x, y во времени изменяются в соответствии с системой уравнений (4.1), так что каждому состоянию системы соответствует пара значений переменных (x, y). Обратно, каждой паре переменных (x, y) соответствует определенное состояние системы. Рассмотрим плоскость с осями координат, на которых отложены значения переменных x,y. Каждая точка М этой плоскости соответствует определенному состоянию системы. Такая плоскость носит название фазовой плоскости и изображает совокупность всех состояний системы. Точка М(x,y) называется изображающей или представляющей точкой. Пусть в начальный момент времени t=t 0 координаты изображающей точки М 0 (x(t 0 ), y(t 0 )). В каждый следующий момент времени t изображающая точка будет смещаться в соответствии с изменениями значений переменных x(t), y(t). Совокупность точек М(x(t), y(t)) на фазовой плоскости, положение которых соответствует состояниям системы в процессе изменения во времени переменных x(t), y(t) согласно уравнениям (4.1), называется фазовой траекторией.

Совокупность фазовых траекторий при различных начальных значениях переменных дает легко обозримый "портрет" системы. Построение фазового портрета позволяет сделать выводы о характере изменений переменных x, y без знания аналитических решений исходной системы уравнений (4.1). Для изображения фазового портрета необходимо построить векторное поле направлений траекторий системы в каждой точке фазовой плоскости. Задавая приращение t>0, получим соответствующие приращения x и y из выражений: x=P(x,y) t, y=Q(x,y) t. Направление вектора dy/dx в точке (x, y) зависит от знака функций P(x, y), Q(x, y) и может быть задано таблицей:

P(x,y)>0, Q(x,y)>0 P(x,y)

Задача построения векторного поля упрощается, если получить выражение для фазовых траекторий в аналитическом виде. Для этого разделим второе из уравнений системы (4.1) на первое: (4.2) Решение этого уравнения y = y(x, c), или в неявном виде F(x,y)=c, где с – постоянная интегрирования, дает семейство интегральных кривых уравнения (4.2) фазовых траекторий системы (4.1) на плоскости x, y.

Метод изоклин Для построения фазового портрета пользуются методом изоклин – на фазовой плоскости наносят линии, которые пересекают интегральные кривые под одним определенным углом. Уравнение изоклин легко получить из (4.2). Положим где А – определенная постоянная величина. Значение А представляет собой тангенс угла наклона касательной к фазовой траектории и может принимать значения от – до +. Подставляя вместо dy/dx в (4.2) величину А получим уравнение изоклин: (4.3)

Уравнение (4.3) определяет в каждой точке плоскости единственную касательную к соответствующей интегральной кривой за исключением точки, где P (x,y) = 0, Q (x,y) = 0, в которой направление касательной становится неопределенным, так как при этом становится неопределенным значение производной: Эта точка является точкой пересечения всех изоклин – особой точкой. В ней одновременно обращаются в нуль производные по времени переменных x и y. Таким образом, в особой точке скорости изменения переменных равны нулю. Следовательно, особая точка дифференциальных уравнений фазовых траекторий (4.2) соответствует стационарному состоянию системы (4.1), а ее координаты – суть стационарные значения переменных x, y.

Особый интерес представляют главные изоклины: dy/dx=0, P(x,y)=0 – изоклина горизонтальных касательных и dy/dx=, Q(x,y)=0 – изоклина вертикальных касательных. Построив главные изоклины и найдя точку их пересечения (x,y), координаты которой удовлетворяют условиям: мы найдем тем самым точку пересечения всех изоклин фазовой плоскости, в которой направление касательных к фазовым траекториям неопределенно. Это – особая точка, которая соответствует стационарному состоянию системы (рис. 4.2). Система (4.1) обладает столькими стационарными состояниями, сколько точек пересечения главных изоклин имеется на фазовой плоскости. Каждая фазовая траектория соответствует совокупности движений динамической системы, проходящих через одни и те же состояния и отличающихся друг от друга только началом отсчета времени.

Таким образом, фазовые траектории системы – это проекции интегральных кривых в пространстве всех трех измерений x, y, t на плоскость x, y (рис.4.3). Если условия теоремы Коши выполнены, то через каждую точку пространства x, y, t проходит единственная интегральная кривая. То же справедливо, благодаря автономности, для фазовых траекторий: через каждую точку фазовой плоскости проходит единственная фазовая траектория.

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

Линейные системы. Рассмотрим систему двух линейных уравнений: (4.4) Здесь a, b, c, d - константы, x, y декартовы координаты на фазовой плоскости. Общее решение будем искать в виде: (4.5) Подставим эти выражения в (4.4) и сократим на e t : (4.6) Алгебраическая система уравнений (4.6) с неизвестными A, B имеет ненулевое решение лишь в том случае, если ее определитель, составленный из коэффициентов при неизвестных, равен нулю:

Раскрывая этот определитель, получим характеристическое уравнение системы: (4.7) Решение этого уравнения дает значения показателя 1,2, при которых возможны ненулевые для A и B решения уравнения (4.6). Эти значения суть (4.8) Если подкоренное выражение отрицательно, то 1,2 комплексно сопряженные числа. Предположим, что оба корня уравнения (4.7) имеют отличные от нуля действительные части и что нет кратных корней. Тогда общее решение системы (4.4) можно представить в виде линейной комбинации экспонент с показателями 1, 2 : (4.9) Для анализа характера возможных траекторий системы на фазовой плоскости используем линейное однородное преобразование координат, которое позволит привести систему к каноническому виду: (4.10)

допускающее более удобное представление на фазовой плоскости по сравнению с исходной системой (4.4). Введем новые координаты ξ, η по формулам: (4.1) Из курса линейной алгебры известно, что в случае неравенства нулю действительных частей 1, 2 исходную систему (4.4) при помощи преобразований (4.11) всегда можно преобразовать к каноническому виду (4.10) и изучать ее поведение на фазовой плоскости ξ, η. Рассмотрим различные случаи, которые могут здесь представиться.

Корни λ 1, λ 2 – действительны и одного знака В этом случае коэффициенты преобразования действительны, мы переходим от действительной плоскости x,y к действительной плоскости ξ, η. Разделив второе из уравнений (4.10) на первое, получим: (4.12) Интегрируя это уравнение, находим:, где.(4.13) Условимся понимать под λ 2 корень характеристического уравнения с большим модулем, что не нарушает общности нашего рассуждения. Тогда, поскольку в рассматриваемом случае корни λ 1, λ 2 – действительны и одного знака, a>1, и мы имеем дело с интегральными кривыми параболического типа. Все интегральные кривые (кроме оси η, которой соответствует ) касаются в начале координат оси ξ, которая также является интегральной кривой уравнения (4.11). Начало координат является особой точкой. Выясним теперь направление движений изображающей точки вдоль фазовых траекторий. Если λ 1, λ 2 – отрицательны, то, как видно из уравнений (4.10), |ξ|, |η| убывают с течением времени. Изображающая точка приближается к началу координат, никогда, однако, не достигая его. В противном случае это противоречило бы теореме Коши, которая утверждает, что через каждую точку фазовой плоскости проходит лишь одна фазовая траектория.

Такая особая точка, через которую проходят интегральные кривые, подобно тому, как семейство парабол проходит через начало координат, носит название узла (рис. 4.5) v Такая особая точка, через которую проходят интегральные кривые, подобно тому, как семейство парабол проходит через начало координат, носит название узла (рис. 4.5)

Состояние равновесия типа узел при λ 1, λ 2 0, то |ξ|, |η| возрастают с течением времени и изображающая точка удаляется от начала координат. В этом случае особая точка – неустойчивый узел. На фазовой плоскости x, y общий качественный характер поведения интегральных кривых сохранится, но касательные к интегральным кривым не будут совпадать с осями координат. Угол наклона этих касательных будет определяться соотношением коэффициентов α, β, γ, δ в уравнениях (4.11).

Корни λ 1, λ 2 – действительны и разных знаков. Преобразование от координат x,y к координатам ξ, η опять действительное. Уравнения для канонических переменных снова имеют вид (4.10), но теперь знаки λ 1, λ 2 различны. Уравнение фазовых траекторий имеет вид: где,(4.14) Интегрируя (4.14), находим (4.15) Это уравнение определяет семейство кривых гиперболического типа, где обе оси координат – асимптоты (при a=1 мы имели бы семейство равнобочных гипербол). Оси координат и в этом случае являются интегральными кривыми – это будут единственные интегральные кривые, проходящие через начало координат. Каждая из них состоит из трех фазовых траекторий: из двух движений к состоянию равновесия (или от состояния равновесия) и из состояния равновесия.

Все остальные интегральные кривые – суть гиперболы, не проходящие через начало координат (рис. 4.6) Такая особая точка носит название «седло». Линии уровня вблизи горной седловины ведут себя подобно фазовым траекториям в окрестности седла.

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

и система имеет своими состояниями равновесия все точки прямой: Остальные интегральные кривые представляют собой семейство параллельных прямых с угловым коэффициентом, по которым изображающие точки либо приближаются к состоянию равновесия, либо удаляются от него в зависимости от знака второго корня характеристического уравнения λ 2 = a+d. (Рис.4.7) В этом случае координаты состояния равновесия зависят от начального значения переменных. Рис Фазовый портрет системы, один из характеристических корней которой равен нулю, а второй отрицателен.

Корни λ 1, λ 2 – комплексные сопряженные В этом случае при действительных x и y мы будем иметь комплексные сопряженные ξ, η (4.10). Однако, вводя еще одно промежуточное преобразование, можно и в этом случае свести рассмотрение к действительному линейному однородному преобразованию. Положим: (4.16) где a,b, и u,v – действительные величины. Можно показать, что преобразование от x,y к u,v является при наших предположениях действительным, линейным, однородным с детерминантом, отличным от нуля. В силу уравнений (4.10, 4.16) имеем: откуда (4.17) Разделив второе из уравнений на первое, получим:

которое легче интегрируется, если перейти к полярной системе координат (r, φ). После подстановки получим, откуда: (4.18) Таким образом, на фазовой плоскости u, v мы имеем дело с семейством логарифмических спиралей, каждая из которых имеет асимптотическую точку в начале координат. Особая точка, которая является асимптотической точкой всех интегральных кривых, имеющих вид спиралей, вложенных друг в друга, называется фокусом (рис.4.8). Рис Фазовый портрет системы в окрестности особой точки типа фокус на плоскости координат u, v.

Рассмотрим характер движения изображающей точки по фазовым траекториям. Умножая первое из уравнений (4.17) на u, а второе на v и складывая, получаем: где Пусть a 1 < 0 (a 1 =Reλ). Изображающая точка тогда непрерывно приближается к началу координат, не достигая его в конечное время. Это означает, что фазовые траектории представляют собой скручивающиеся спирали и соответствуют затухающим колебаниям переменных. Это – устойчивый фокус. В случае устойчивого фокуса, как и в случае устойчивого узла, выполнено не только условие Ляпунова, но и более жесткое требование. Именно, при любых начальных отклонениях система по прошествии времени вернется как угодно близко к положению равновесия. Такая устойчивость, при которой начальные отклонения не только не нарастают, но затухают, стремясь к нулю, называют абсолютной устойчивостью. Если в формуле (4.18) a 1 >0, то изображающая точка удаляется от начала координат, и мы имеем дело с неустойчивым фокусом. При переходе от плоскости u,v к фазовой плоскости x,y спирали также останутся спиралями, однако будут деформированы. Рассмотрим теперь случай, когда a 1 =0. Фазовыми траекториями на плоскости u, v будут окружности которым на плоскости x,y соответствуют эллипсы:

Таким образом, при a 1 =0 через особую точку x=0, y=0 не проходит ни одна интегральная кривая. Такая изолированная особая точка, вблизи которой интегральные кривые представляют собой замкнутые кривые, в частности, эллипсы, вложенные друг в друга и охватывающие особую точку, называется центром. Таким образом, возможны шесть типов состояния равновесия в зависимости от характера корней характеристического уравнения (4.7). Вид фазовых траекторий на плоскости x, y для этих шести случаев изображен на рис. 4.9.

Бифуркационная диаграмма Введем обозначения: (4.11) Тогда характеристическое уравнение запишется в виде: (4.12) Рассмотрим плоскость с прямоугольными декартовыми координатами, и отметим на ней области, соответствующие тому или иному типу состояния равновесия, который определяется характером корней характеристического уравнения (4.13) Условием устойчивости состояния равновесия будет наличие отрицательной действительной части у 1 и 2. Необходимое и достаточное условие этого – выполнение неравенств > 0, > 0.

На диаграмме (4.15) этому условию соответствуют точки, расположенные в первой четверти плоскости параметров. Особая точка будет фокусом, если 1 и 2 комплексны. Этому условию соответствуют те точки плоскости, для которых, т.е. точки между двумя ветвями параболы 2 = 4. Точки полуоси = 0, >0, соответствуют состояниям равновесия типа центр. Аналогично, 1 и 2 - действительны, но разных знаков, т.е. особая точка будет седлом, если

Если коэффициенты линейной системы a, b, c, d зависят от некоторого параметра, то при изменении этого параметра будут меняться и величины,. При переходе через границы характер фазового портрета качественно меняется. Поэтому такие границы называются бифуркационными – по разные стороны от границы система имеет два топологически различных фазовых портрета и, соответственно два разных типа поведения. На диаграмме видно, как могут проходить такие изменения. Если исключить особые случаи – начало координат, – то легко видеть, что седло может переходить в узел, устойчивый или неустойчивый при пересечении оси ординат. Устойчивый узел может перейти либо в седло, либо в устойчивый фокус, и т.д. Отметим, что переходы устойчивый узел – устойчивый фокус и неустойчивый узел – неустойчивый фокус не являются бифуркационными, так как топология фазового пространства при этом не меняется. Более подробно мы поговорим о топологии фазового пространства и бифуркационных переходах в лекции 6. При бифуркационных переходах меняется характер устойчивости особой точки. Например, устойчивый фокус через центр может переходить в неустойчивый фокус. Эта бифуркация называется бифуркацией Андронова-Хопфа по именам исследовавших ее ученых. При этой бифуркации в нелинейных системах происходит рождение предельного цикла, и система становится автоколебательной (см. лекцию 8).

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

Пусть биологическая система описывается системой двух автономных дифференциальных уравнений второго порядка общего вида: (5.1) Стационарные значения переменных системы определяются из алгебраических уравнений: (5.2) Стационарные состояния соответствуют особым точкам дифференциального уравнения первого порядка, определяющего интегральные кривые: (5.3)

Математический анализ поведения траекторий этой системы на фазовой плоскости связан с именами французского математика Анри Пуанкаре и русского математика и механика Александра Михайловича Ляпунова ( ). Ляпунов показал, что в большом числе случаев анализ устойчивости стационарного состояния нелинейной системы можно заменить анализом устойчивости системы, линеаризованной в окрестности стационарного состояния. Рассмотрим характер поведения переменных при некотором небольшом отклонении системы от состояния равновесия. Введем вместо переменных x, y новые независимые переменные,, определив их как смещения относительно равновесных значений переменных (5.4) Подставив эти выражения в (5.1), получим: (5.5), так как - координаты особой точки.

Предположим, что функции P и Q непрерывны и имеют непрерывные производные не ниже первого порядка. Тогда мы можем разложить правые части уравнений (5.5) в ряд Тейлора по переменным,. (5.6) где (5.7) Учтем, что по определению особой точки и отбросим в уравнениях (5.6) нелинейные члены. Получим систему линейных уравнений с постоянными коэффициентами систему первого приближения: (5.8)

Решение этой системы было рассмотрено в Лекции 4. Оно определяется корнями характеристического уравнения системы: (5.9) Ляпунов показал, что в случае, если оба корня уравнения (5.9): (5.10) имеют отличные от нуля действительные части, исследование уравнений первого приближения (5.8) всегда дает правильный ответ на вопрос о типе устойчивости состояния равновесия в системе (5.1). А именно: · если оба корня имеют отрицательную действительную часть и, следовательно, все решения уравнений первого приближения (5.8) затухают, то состояние равновесия устойчиво; · если хотя бы один корень имеет положительную действительную часть, то есть система (5.8) имеет нарастающие решения, то состояние равновесия неустойчиво. Если действительные части обоих корней характеристического уравнения равны нулю или если один корень равен нулю, а другой отрицателен, то уравнения (5.8) не дают ответа на вопрос об устойчивости состояния равновесия, и необходимо рассматривать члены более высокого порядка малости в разложении в ряд Тейлора правых частей уравнений (5.6).

В случае, когда оба корня характеристического уравнения имеют отличные от нуля действительные части (грубые системы), уравнение первого приближения определяют не только устойчивость стационарного состояния, но и характер фазовых траекторий в достаточно малой его окрестности. Как и в случае линейных уравнений (Лекция 4) здесь возможны пять типов грубых состояний равновесия: устойчивый узел, неустойчивый узел, устойчивый фокус, неустойчивый фокус и седло. Для исследования типов состояний равновесий удобно пользоваться диаграммой, изображенной на рис Для системы (5.1): (5.11) (5.12) Грубым состояниям равновесия соответствуют все точки плоскости параметров,, лежащие вне оси =0 и полуоси =0, >0. Точкам оси = 0 и полуоси = 0, >0 соответствуют негрубые состояния равновесия (негрубые особые точки). Их свойства могут быть изменены сколь угодно малыми изменениями правых частей уравнений (5.1) за счет сколь угодно малых изменений функций P(x,y), Q(x,y) и их производных. Поэтому характер негрубых состояний равновесия (в частности, устойчивость) уже не определяется значениями коэффициентов в правых частях уравнений первого приближения (5.8). В отличие от линейных систем, уже при небольших изменений в правых частях содержащихся там нелинейных членов может произойти качественное изменение фазового портрета бифуркация.

Примеры 1. Кинетические уравнения Лотки Лоткой была исследована гипотетическая химическая реакция: Модель очень простая и служит хорошей иллюстрацией применения исследования устойчивости стационарного состояния системы методом линеаризации. Пусть в некотором объеме находится в избытке вещество А. Молекулы А с некоторой постоянной скоростью превращаются в молекулы вещества X (реакция нулевого порядка). Вещество X может превращаться в вещество Y, причем скорость этой реакции тем больше, чем больше концентрация вещества Y – реакция второго порядка. В схеме это отражено обратной стрелкой над символом Y. Молекулы Y в свою очередь необратимо распадаются, в результате образуется вещество B (реакция первого порядка). Запишем систему уравнений, описывающих реакцию: (5.13) Здесь x, y, B - концентрации химических компонентов. Первые два уравнения этой системы не зависят от B, поэтому их можно рассматривать отдельно. Рассмотрим стационарное решение системы:

Из этих условий получим систему алгебраических уравнений, связывающих равновесные концентрации : (5.14) Координаты особой точки: Исследуем устойчивость этого стационарного состояния методом Ляпунова. Введем новые переменные,, характеризующие отклонения переменных от равновесных концентраций : Линеаризованная система в новых переменных имеет вид: (5.15)

Отметим, что величины отклонений от стационарных значений переменных, могут менять знак, в то время как исходные переменные x, y, являющиеся концентрациями, могут быть только положительными. Запишем характеристическое уравнение системы (4.3): или Корни характеристического уравнения: Фазовый портрет системы (5.13) изображен на рис. 5.1.

Рис Фазовый портрет системы а – устойчивый фокус, б – устойчивый узел.

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

Соотношение параметров соответствует изменению типа особой точки системы уравнений (5.13). Рассмотрим плоскость параметров, где по оси абсцисс отложены значения константы k 2, а по оси ординат – произведение k 0 k 1. Парабола k 0 k 1 = 4 k 2 2 делит изображенную на рис. 5.2 плоскость параметров на две области – устойчивых узлов и устойчивых фокусов. Задавая те или иные значения параметров, можно получить колебательный и бесколебательный режимы изменения концентраций веществ x и y, и фазовый портрет системы, соответственно, будет собой представлять фокус (а) или узел (б), изображенные соответственно на рис 5.1а, и 5.1б. Если в системе установятся стационарные концентрации веществ x и y, это приведет к установлению постоянной скорости прироста концентрации вещества В в третьем уравнении системы (5.13): Ясно, что в действительности такая система реализоваться не может, так как в ней при t концентрация вещества В стремится к бесконечности. Однако система, подобная системе реакций Лотки, может представлять собой фрагмент более сложной химической системы. Исследованные нами уравнения правильно описывают поведение компонентов x и y, если приток вещества x (скорость его постоянна и равна k 0 ) осуществляется из большого «резервуара», а отток вещества y – в большой «резервуар» (значение В очень велико). При этих предположениях на малых промежутках времени (по сравнению с временем существенного изменения заполненности емкости B) наше рассмотрение является вполне правомерным.

2. Модель Вольтерра В качестве второго примера рассмотрим классическую модель взаимодействия видов, которая впервые была предложена В. Вольтерра в тридцатые годы XX века для объяснения периодических изменений числа особей, так называемую вольтерровскую модель «хищник-жертва». Более подробно модели взаимодействия видов мы рассмотрим в Лекции 9. Пусть в некотором замкнутом районе живут хищники и жертвы, например, зайцы и волки. Зайцы питаются растительной пищей, имеющейся всегда в достаточном количестве. Волки могут питаться лишь зайцами. Обозначим число зайцев (жертв) x, а число волков (хищников) – y. Так как количество пищи у зайцев неограниченно, мы можем предположить, что они размножаются со скоростью, пропорциональной их числу: (5.16) Если рождаемость зайцев превышает их смертность, 0. Выражение (5.16) соответствует автокаталитической реакции первого порядка. Пусть убыль зайцев пропорциональна вероятности встречи зайца с волком, т.е. пропорциональна произведению численностей xy. Можно предположить по аналогии с бимолекулярными реакциями, где вероятность появления новой молекулы пропорциональна вероятности встречи двух молекул, что и количество волков нарастает тем быстрее, чем чаще происходят их встречи с зайцами, а именно, пропорционально xy. Кроме того, имеет место процесс естественной смертности волков, причем скорость смертности пропорциональна их количеству.

Эти рассуждения приводят к системе уравнений для изменений численности зайцев-жертв x и волков-хищников y. (5.17) Покажем, что система уравнений (5.17) имеет на фазовой плоскости переменных xy ненулевую особую точку типа центр. Координаты этой особой точки легко найти, приравняв правые части уравнений системы (5.17) нулю. Это дает стационарные ненулевые значения: Так как все параметры положительны, точка расположена в положительном квадранте фазовой плоскости. Линеаризация системы вблизи этой точки дает: (5.18)

Рис Фазовый портрет системы Особая точка типа «центр». а – параметры системы: x = 4, xy = 0,3, y = yx = 0,4 б – параметры системы: x =2, xy = 0,3, y = yx = 0,4

Здесь, - отклонения численностей от их стационарных значений: Характеристическое уравнение системы (5.18): Корни этого уравнения чисто мнимые: Таким образом, исследование системы показывает, что траектории вблизи особой точки являются концентрическими эллипсами, а сама особая точка – центром. Расcматриваемая модель Вольтерра и вдали от особой точки имеет замкнутые траектории, хотя форма этих траекторий уже отличается от эллипсоидальной, и определяется параметрами системы (рис. 5.3).

Изменения численности жертвы и хищника во времени представляют собой колебания, причем колебания численности хищника отстают по фазе от колебаний жертв. Как мы уже отмечали в Лекции 4, особая точка типа центр устойчива по Ляпунову, но не асимптотически. Покажем на данном примере, в чем это проявляется. Пусть колебания x(t) и y(t) происходят таким образом, что изображающая точка движется по фазовой траектории 1 (рис 5.3). В момент, когда точка находится в положении М 1, в систему добавляется извне некоторое число особей y такое, что изображающая точка переходит скачком из точки M 1 в точку M 2. Если после этого систему предоставить самой себе, колебания x(t), y(t) уже будут происходить с большими амплитудами, чем прежде, и изображающая точка будет двигаться по траектории 2. Это и означает, что колебания в системе неустойчивы: они навсегда изменяют свои характеристики при внешнем воздействии.

Рис Кривые численности зайца и рыси в Канаде (по К. Вилли, В. Детье, 1974)

В дальнейшем мы рассмотрим модели, описывающие устойчивые колебательные режимы, и покажем, что на фазовой плоскости такие асимптотически устойчивые периодические движения описываются предельными циклами. На рис. 5.4 кривые колебаний численности пушных зверей по данным компании Гудзонова залива о числе заготовленных шкурок. Во всех классических учебниках в течение многих лет колебательный характер этих изменений приводили как подтверждение гипотез, положенных в основу модели Вольтерра, которую мы только что рассмотрели. Действительно, периоды колебаний численности зайцев (жертв) и рысей (хищников) примерно одинаковы и составляют порядка 9 – 10 лет. При этом максимум численности зайцев опережает, как правило, максимум численности рысей на один год. Можно полагать, что мы видим регулярные колебания, осложненные случайными факторами, связанными с погодой и проч. Однако возможна и другая интерпретация этих данных наблюдений на основе моделей детерминированного хаоса. О дискретных моделях такого типа мы уже говорили в Лекции 3. Непрерывные модели популяционной динамики, приводящие к детерминированному хаосу, мы рассмотрим в Лекции 9.

Серьезным недостатком рассмотренной модели Вольтерра является неустойчивость решений по отношению к малым случайным воздействиям, приводящим к изменению переменных. Кроме того, в силу «негрубости» этой системы произвольно малое изменение вида правых частей уравнений (величин параметров системы) приведет к изменению типа особой точки, и, следовательно, к изменению характера фазовых траекторий. Поскольку природные системы подвергаются огромному количеству случайных воздействий, реалистическая модель должна быть по отношению к ним устойчивой. Поэтому негрубые системы не могут давать адекватное описание природных явлений. Различные модификации рассмотренной нами системы, изученные самим Вольтерра и другими авторами, лишены этих недостатков. Наиболее широко известные из них будут рассмотрены в Лекции 9. Здесь мы остановимся на модели, которая учитывает самоограничение в росте обеих популяций. На ее примере видно, как может меняться характер решений при изменении параметров системы. Итак, рассмотрим систему:

Система (5.19) отличается от ранее рассмотренной системы наличием в правых частях членов: Эти члены отражают тот факт, что численность популяции жертв не может расти до бесконечности даже в отсутствие хищников в силу ограниченности пищевых ресурсов, ареала существования и проч. Такие же «самоограничения» накладываются на популяцию хищников. Система имеет два стационарных решения: нулевое и ненулевое. Анализ показывает, что нулевое решение представляет собой неустойчивый узел. Рассмотрим систему алгебраических уравнений, решение которых дает координаты ненулевого стационарного состояния. (5.20)

Стационарное решение: Корни характеристического уравнения системы, линеаризованной в окрестности особой точки:

Из выражения для характеристических чисел видно, что если выполнено условие то численности хищников и жертв совершают во времени затухающие колебания. Система имеет особую точку – устойчивый фокус. Рис Фазовый портрет системы 5.19 а – устойчивый фокус, параметры системы: x = 2, xy = 18, x =1, y = 3, yx = 5, y =1 б – устойчивый узел, параметры системы: x = 2, xy = 1, x =1, y = 3, yx = 1, y =1

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

МЕТОД ФУНКЦИЙ ЛЯПУНОВА ИССЛЕДОВАНИЯ УСТОЙЧИВОСТИ СТАЦИОНАРНОГО СОСТОЯНИЯ. При аналитическом исследовании устойчивости стационарного состояния часто используется метод подбора функции, линии уровня которой представляют собой замкнутые траектории – «ловушки» для фазовых траекторий системы типа (5.1) Этот метод применим к автономной системе уравнений n-го порядка (5.21) где f i (0,0,…, 0) = 0, (i = 1,…, n). Он состоит в непосредственном исследовании устойчивости ее стационарного состояния при помощи подходящим образом подобранной функции Ляпунова. Метод основан на двух теоремах.

Теорема 2 Если существует дифференцируемая функция V (x 1,…,x n ), удовлетворяющая в окрестности начала координат следующим условиям: a) V(x 1,…,x n ) = 0 и сколь угодно близко от начала координат имеются точки, в которых V(x 1,…,x n ) > 0; б) причем лишь при x 1 =…= x n = 0, то точка покоя системы (5.21) неустойчива. С доказательством этих теорем можно познакомиться в книге Л.Э. Эльсгольц «Теория дифференциальных уравнений» или в других учебниках по теории дифференциальных уравнений. Общего методы построения функции Ляпунова не существует. Однако для линейных автономных систем ее следует искать в виде: и т.п., подбирая надлежащим образом коэффициенты a > 0, b > 0. Для нелинейных систем a и b могут быть произвольных знаков.

Примеры 1. Рассмотрим линейную систему: Выберем функцию Ляпунова: V = x 2 +y 2. Тогда Это выражение всегда отрицательно при х 0, т.к. в скобках стоят четные степени x. Следовательно, точка (0, 0) устойчива.

2. Рассмотрим систему уравнений, описывающую конкуренцию видов, численности которых x и y. Каждый из видов размножается в соответствии с логистическим законом, а при встрече (произведения в правых частях уравнений), численность как одного, так и другого вида уменьшается. (5.22) Исследуем стационарное состояние, соответствующее сосуществованию видов ( x, y) – ненулевое для x и y. Его координаты: (5.23)

В. Вольтерра показал, что стационарное состояние (5.23) устойчиво для параметров системы a > 0, b 1, построив функцию Ляпунова: Ее производная равна и отрицательна при малых значениях коэффициентов a, b и x, y > 0. Доказательство приведено в книге В. Вольтерра. «Математическая теория борьбы за существование» (М., 1976)

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

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

Благодаря фотосинтезу образуется органическое вещество из углекислого газа и воды с использованием энергии солнечного света и неорганических веществ из почвы и воды. Фотосинтез также служит источником земного кислорода, необходимого для дыхания всех аэробных организмов. Степень подробности моделирования изучаемых явлений зависит от цели моделирования. Однако в любом случае задача моделирования заключается в том, чтобы построить модель явления, содержащую возможно меньшее число переменных и произвольных параметров, и в то же время правильно отражающую свойства явления. Учет временной иерархии процессов позволяет сократить число дифференциальных уравнений. «Совсем медленные» переменные не меняются на временах рассматриваемых процессов, и их можно считать постоянными параметрами. Для «быстрых» переменных можно вместо дифференциальных уравнений записать алгебраические уравнения для их стационарных значений, поскольку «быстрые» переменные достигают своих стационарных значений практически мгновенно по сравнению с «медленными».

Средние, быстрые и медленные времена. Пусть имеется три группы переменных с различными характерными временами: Переменные изменяются с разными характерными временами, причем Пусть мы наблюдаем за переменной y, характерное время изменения которой T y. Тогда за время T y «совсем медленная» переменная z практически не будет изменяться, и ее можно считать постоянным параметром, обозначим его z *. Система дифференциальных уравнений с учетом этого обстоятельства будет содержать два уравнения и может быть записана в виде:

Отметим, что z * не является истинно стационарным значением, «медленная» переменная z будет продолжать меняться и «вести» за собой более быстрые переменные x и y. В этом смысле медленная переменная является ведущей, или «параметром порядка». Рассмотрим теперь уравнение для x. Эта «быстрая» переменная изменяется значительно быстрее, чем y, и за время T y успеет достичь своего стационарного значения. Значит, для переменной x дифференциальное уравнение можно заменить алгебраическим: P(x, y, z * ) =0 Или Таким образом, благодаря учету иерархии времен, исходную систему из трех дифференциальных уравнений удается свести к одному дифференциальному уравнению для переменной y: В химической кинетике метод такой редукции системы был впервые предложен Боденштейном и носит название метода квазистационарных концентраций (КСК).

Обычно он применяется для систем химических реакций, промежуточные продукты которых являются частицами с высокой реакционной способностью. К ним относятся каталитические процессы, свободно радикальные и цепные реакции. В процессах с участием активных промежуточных частиц разность скоростей образования v о и расхода v р этих частиц мала по сравнению с этими скоростями. Режим называется квазистационарным, а отвечающие ему концентрации активных промежуточных веществ квазистационарными концентрациями. Дифференциальные уравнения для промежуточных соединений: можно заменить алгебраическими: Из l алгебраических уравнений можно выразить l квазистационарных концентраций промежуточных химических соединений. По мере расходования исходных веществ, квазистационарные концентрации промежуточных соединений будут меняться, но если время установления квазистационарного режима мало, он не будет нарушаться в течение всего процесса. Конечно, такое рассмотрение не правомерно для начальных стадий процесса, когда R i меняются от нуля до своих квазистационарных значений. Этот период носит название периода индукции.

Теорема Тихонова Математически строгое обоснование применения метода квазистационарных концентраций (редукции системы в соответствии с иерархией времен) и формулировка условий его применимости дана в работе А.Н. Тихонова (1952). Рассмотрим простейший случай двух дифференциальных уравнений (6.1) Пусть y медленная, а x быстрая переменная. Это означает, что отношение приращений y и x за короткий промежуток времени t много меньше единицы: y/ x1. Первое уравнение системы можно представить в виде:

Разделив левую и правую часть уравнения на А и обозначив =1/A, получим полную систему уравнений, тождественную исходной: (6.2) где

Скорость движения по горизонтальным участкам траектории dx/dt 1/ =A, т.е. очень велика по сравнению со скоростью движения в окрестности кривой F(x,y)=0. Поэтому общее время достижения некоего состояния на кривой F(x,y) определяется лишь характером движения вдоль этой кривой, т.е. зависит лишь от начальных значений медленной переменной y и не зависит от начальных значений быстрой переменной x. Отметим, что квазистационарные значения быстрых переменных являются функциями не окончательных стационарных значений медленных переменных, а лишь их мгновенных значений. В этом смысле говорят о том, что быстрая переменная «подчинена» медленной.

Теорема Тихонова устанавливает условия редукции системы дифференциальных уравнений с малым параметром (условия замены дифференциальных уравнений для быстрых переменных алгебраическими). Запишем систему N уравнений, часть из которых содержит малый параметр перед производной. (6.4) (6.5) Назовем систему (6.4) присоединенной, а систему (6.5) вырожденной.

Решение полной системы ( ) стремится к решению вырожденной системы (6.5) при 0, если выполняются следующие условия: а) решение полной и присоединенной системы единственно, а правые части непрерывны; б) решение представляет собой изолированный корень алгебраической системы (в окрестности этого корня нет других корней); в) решение устойчивая изолированная особая точка присоединенной системы (6.4) при всех значениях ; г) начальные условия попадают в область влияния устойчивой особой точки присоединенной системы.

Число начальных условий вырожденной системы меньше, чем полной: начальные значения быстрых переменных не используются в вырожденной системе. Согласно теореме Тихонова, если выполняется условие в), результат не зависит от начальных условий для переменных присоединенной системы. Таким образом, необходимым условием редукции является наличие малого параметра в уравнениях (6.4). Представляет интерес система двух дифференциальных уравнений вида (6.2), в которой особая точка расположена на неустойчивой ветви кривой F(x,y)=0. Такая система совершает релаксационные колебательные движения. Вопрос о релаксационных колебаниях мы обсудим в лекции 8. Теорема Тихонова явно или неявно применяется при исследовании практически любых моделей биологических систем, в этом мы убедимся в дальнейшем (лекции 7-12).

Фермент-субстратная реакция Михаэлиса-Ментен Классическим примером является модель базовой ферментативной реакции, предложенная Михаэлисом и Ментен в 1913 г. Схема реакции может быть представлена в виде: E + S ES, ES P + E (6.6) Схема означает, что субстрат S соединяется с ферментом E в комплекс ES, в котором происходит химическое превращение, и который затем распадается на фермент E и продукт P. По закону действующих масс, скорость реакции пропорциональна произведению концентраций. Обозначив концентрации реагентов малыми буквами: s=[S], e=[E], c=[ES], p=[P], получим систему уравнений: (6.7)

В системе (6.7.) учтены следующие процессы: · Субстрат S расходуется, образуя комплекс ES (бимолекулярная реакция), и его концентрация увеличивается при распаде комплекса; · Фермент E расходуется на образование комплекса ES, его концентрация увеличивается при распаде комплекса. · Комплекс ES образуется из фермента E и субстрата S (бимолекулярная реакция) и распадается на субстрат S и фермент E. · Продукт P образуется при распаде комплекса. Для полной математической формулировки задачи Коши необходимо задать начальные условия: s 0 (0)=s 0, e(0)=e 0, c(0)=0, p(0)=0.(6.8) Уравнения (6.7) не являются независимыми. Кроме того, последнее уравнение отделяется от первых трех. Если система первых трех уравнений решена, концентрация продукта может быть рассчитана по формуле: В соответствии со схемой реакций ( ) общее количество фермента, свободного и связанного в комплекс, сохраняется:

е(t) + с(t) = e 0. Это условие позволяет одно из дифференциальных уравнений системы (6.7) заменить алгебраическим, и модель сводится к двум дифференциальным уравнениям: (6.9) c начальными условиями: s 0 (0)=s 0, c(0)=0. Введем безразмерные переменные и параметры: (6.10) Запишем уравнения (6.9) в безразмерном виде: (6.11) Из (6.10) следует, что (К )>0.

Поскольку реакция превращения фермент- субстратного комплекса необратима, уже из схемы реакций (6.6) ясно, что с течением времени весь субстрат будет превращен в продукт, и в стационарном состоянии концентрации и субстрата и комплекса станут равны нулю: x=0, y=0. Систему (6.7) нельзя решить аналитически. Проанализируем качественно, как ведут себя x(t) и y(t). Вблизи =0 dx/d 0, y растет от y=0 до величины y=x/(x+K), при которой правая часть уравнения для dy/dt обращается в нуль. После этого величина y будет уменьшаться до нуля. Таким образом, концентрация фермент-субстратного комплекса y проходит через максимум. В это время величина x (концентрация субстрата) монотонно уменьшается. Относительная концентрация свободного фермента e/e 0 сначала убывает а затем снова возрастает до величины e/e 0 =1, поскольку с течением времени субстрат исчерпывается, и все меньшая доля фермента оказывается связанной. Кинетические кривые изображены на рис. 6.2.

Предположим, что концентрация субстрата значительно превышает концентрацию фермента: s 0 >>e 0. Тогда из соотношений (6.10) следует, что

В размерном виде это – классическая формула Михаэлиса - Ментен для кинетики изменения субстрата в ферментативной реакции: (6.14) Таким образом, формула (6.14) верно отражает изменение концентрации субстрата, но ничего не может сказать об изменении концентраций свободного фермента и фермент-субстратного комплекса, которые на малых временах ведут себя немонотонно (см. рис. 6.2). Величина K m называется константой Михаэлиса и имеет размерность концентрации. При s

Бифуркации динамических систем. Мы рассматриваем динамические модели биологических процессов, то есть считаем, что система может быть описана системой дифференциальных уравнений: (6.15) Здесь x – вектор переменных, - вектор параметров. Пусть – стационарное решение – особая точка системы, координаты которой представляют собой решение системы алгебраических уравнений: (6.16) Зафиксируем некоторое = *, и рассмотрим фазовые портреты системы при данном значении параметра, а также при > * и < *. Фазовые портреты топологически эквивалентны, если существует невырожденное непрерывное преобразование координат, которое переводит все элементы одного фазового портрета в элементы другого. Для того чтобы представить себе такое преобразование на поверхности, представим себе, что поверхность резиновая, ее можно сжимать и изгибать, но нельзя перекручивать. При таких преобразованиях все начальные точки будут однозначно переходить в точки деформированной «резиновой» поверхности, незамкнутые кривые будут переходить в незамкнутые, замкнутые – в замкнутые, связность множеств не будет нарушаться. Такое преобразование происходит с фазовыми кривыми при невырожденном непрерывном преобразовании координат. Недаром говорят, что топология – это «резиновая геометрия»

Если фазовые портреты при значениях > * и < * топологически не эквивалентны, это означает, что при происходит качественная перестройка системы. Тогда говорят, что * бифуркационное значение параметра. Простейший пример бифуркационного значения параметра – нулевое значение собственной константы скорости роста в уравнении экспоненциального роста (2.7): При r>0 стационарное значение x=0 – неустойчиво, при r

Для оценки «сложности» бифуркации вводится понятие «коразмерности». Коразмерность k совпадает с числом параметров, при независимой вариации которых эта бифуркация происходит. В системе происходит бифуркация коразмерности k (codim k, dimension размерность), если в ней выполняются k условий типа равенств. Значение k=0 соответствует отсутствию бифуркации в данной точке. На рис линии представляют собой бифуркации коразмерности 1, а начало координат – бифуркацию коразмерности 2. Бифуркации разделяют на локальные и нелокальные. Все рассмотренные нами ранее бифуркации, а также другие бифуркации смены устойчивости или исчезновения предельного множества в результате слияния с другим предельным множеством (как мы это увидим при параметрическом переключении триггера в лекции 7) – локальные. Они диагностируются с помощью линейного анализа ляпуновских показателей (собственных чисел). Нелокальные бифуркации нельзя определить на основе линейного анализа окрестности стационарного состояния, здесь требуется нелинейный анализ системы. К нелокальным бифуркациям относятся образование сепаратрисных петель, касание аттрактором сепаратрисных кривых или поверхностей. Бифуркации аттракторов принято подразделять на мягкие (внутренние) бифуркации и кризисы (жесткие бифуркации). Внутренние бифуркации приводят к топологическим изменениям самих притягивающих множеств, не затрагивая их бассейнов притяжения – областей, из которых фазовые траектории сходятся к данному аттрактору.

Кризисы – бифуркации аттракторов, сопровождающиеся качественной перестройкой границ областей притяжения (бассейнов) аттракторов. Пример бифуркация слияния устойчивого узла с седлом, в результате чего аттрактор исчезает (рис. 6.5). Часто кроме бифуркационных диаграмм для наглядности строят фазопараметрические диаграммы. В этом случае по одним координатным осям откладывают значения параметров, а по другим – динамические переменные или связанные с ними величины. Получают некоторую гиперповерхность, точки которой соответствуют определенным динамическим режимам, меняющимся с изменением параметров. Бифуркации на таких диаграммах могут проявляться в образовании складок поверхности или в расщеплении ее на несколько частей. Резкие значительные изменения переменных состояния динамической системы, вызванные малыми возмущениями в правых частях уравнений, в частности, малыми изменениями параметров, часто называют катастрофами. Теория катастроф была развита топологом Рене Тома (Thom R. Structural Stability and Morphogenesis. N.Y., 1972). В основу ее была положена разработанная ранее теория особенностей Уитни. Показано, что существует небольшое количество элементарных катастроф, с помощью которых можно локально описать поведение системы.

Модельные системы Для описания событий, происходящих вблизи бифуркационной границы удобно использовать системы самых простых уравнений, обычно – полиномиальных, которые описывают качественные особенности процесса. Такие системы называются модельными и активно используются в теории бифуркаций и в теории катастроф. Например, для системы, которая может быть описана одним автономным дифференциальным уравнением, модельная система имеет вид: Условием вырождения (бифуркации) является равенство нулю коэффициента a, то есть отсутствие в правой части линейного члена. В качестве модельной системы, описывающей бифуркацию коразмерности k, обычно выступает полиномиальная система l k уравнений, зависящая от k малых параметров. При нулевых значениях параметров в системе возникает вырождение, а при вариации параметров происходит бифуркация. В простейшем случае в качестве параметров выступают вещественные части собственных чисел. Размерность модельной системы l совпадает с количеством собственных чисел, вещественные части которых обращаются в нуль при бифуркационном значении параметра. Рассмотрим основные бифуркации – катастрофы.

Седло-узловая бифуркация (складка). Пусть в системе при < * существуют два состояния равновесия: устойчивый узел Q и седло S (Рис. 6.5, а). При = * происходит слияние узла и седла с образованием негрубого состояния равновесия, называемого седло-узлом. (рис. 6.5, б). При > * положение равновесия исчезает (рис. 6.5, в). Переменная x с течением времени стремится к бесконечности. Поскольку в результате бифуркации аттрактор (узел) исчезает, границы бассейнов должны качественно перестроиться. Следовательно, данная бифуркация является кризисом (катастрофой). Простейшая модельная система, описывающая данную бифуркацию, имеет вид: (6.17) Уравнение (6.17) имеет два стационарных состояния Линеаризуем уравнение (6.17) в окрестности стационарного состояния. Собственные значения

Таким образом x 1 устойчивое состояние, x 2 неустойчивое. При = 0 имеем x 1 = x 2 = 0, и собственное значение в этой точке равно нулю. Бифуркация имеет коразмерность 1, так как выделяется одним условием ( )=0. На рис. 6.6 а изображена фазопараметрическая диаграмма системы (6.17). Если бифуркация седло-узел происходит в двупараметрической системе, то в фазопараметрическом пространстве ей соответствует особенность (катастрофа) типа складки вдоль линии l* на плоскости параметров.

Поясним, как можно пользоваться образами теории катастроф при изучении математических моделей на примере модели второго порядка, содержащей переменные x и u. u – это фактически управляющий параметр, бифуркационному значению которого = * соответствует u=0. Пусть x – «быстрая» переменная, но исключить ее нельзя, поскольку система не удовлетворяет условиям Теоремы Тихонова (см. выше), так как быстрый процесс не везде устойчив. «Складка» соответствует модели (6.18) Здесь >> 1, характерное время изменения переменной x будем считать порядка единицы. Изоклина P=0 имеет устойчивую ветвь – аттрактор в форме «складки». При медленном уменьшении u в соответствии с первым уравнением (6.18) при достижении u=0 произойдет срыв изображающей точки, которая либо уйдет на, либо перескочит на другой устойчивый аттрактор. Заметим, что в реальных моделях такая устойчивая ветвь всегда присутствует. Катастрофа типа «складки» появляется в моделях, описывающих релаксационные колебания, «ждущие» режимы и триггерные системы (параметрическое переключение). В распределенных моделях (2 том лекций) модели, имеющие «складки», используются при описании автоволновых процессов и диссипативных структур.

Трехкратное равновесие (сборка) Бифуркация состоит в слиянии трех состояний равновесия – узлов Q 1, Q 2 и седла Q 0 между ними (в рождении двух устойчивых узлов из седла) (рис. 6.7, 6.8) Рис Трансформации фазового портрета при бифуркации «рождение двух узлов из седла». а – фазовый портрет в незаштрихованной области (рис. 6.8 а); б – фазовый портрет на границе l 1 ; в – фазовый портрет на границе l 2 ; в – фазовый портрет в заштрихованной области представлен двумя устойчивыми узлами и седлом между ними. Рис Бифуркация трехкратного равновесия (катастрофа – сборка). а бифуркационная диаграмма, б – фазопараметрическая диаграмма

Бифуркация имеет коразмерность 2 и требует для своего описания как минимум двух параметров. Модельной системой для нее служит уравнение: (6.19) Система имеет три особых точки. Линейный анализ показывает, что при 2 >0 и любом 1 система имеет единственное состояние равновесия Q 0 с отрицательным собственным значением, то есть асимптотически устойчивое. При 2

Границы области бистабильности образованы линиями l 1 и l 2, соответствующими бифуркациям седло-узел, на которых два из состояний равновесия сливаются и исчезают. Линии l 1 и l 2 сходятся в точке А ( 1 = 2 = 0), где одновременно выполняются два условия: ( 1, 2 ) в точках Q 1 и Q 2 одновременно равны нулю, поэтому бифуркация в этой точке, называемая трехкратным равновесием, имеет коразмерность 2. Для уравнения в точке А фазовый портрет представляет собой седло В фазопараметрическом пространстве (рис. 6.8, б) имеет место структура, называемая сборкой. Верхний и нижний лист сборки соответствуют устойчивым состояниям равновесия, а средний – неустойчивому. На ребрах сборки имеют место катастрофы типа складки. Модели, содержащие катастрофу типа сборки, используются при описании релаксационных автоколебаний малой амплитуды, колебательных режимов со смещением средней точки и диссипативных структур ступенчатого типа. Слияние четырех или пяти особых точек приводит к катастрофам типа «ласточкин хвост» (рис. 6.9) и «бабочка». (Арнольд В.И. Теория катастроф. М., Знание, 1983). Фазовые пространства при этом – четырех- и пятимерные. Отметим важное различие катастроф типа складки и сборки. «Складка» не описывает поведение системы на больших временах. Изображающая точка уходит из рассматриваемой области фазового пространства, где справедлива формула (6.18). Катастрофа складка не локализуема, то же относится к катастрофе «ласточкин хвост» с четной коразмерностью (рис. 6.9).

Рис Бифуркация «Ласточкин хвост» В случае сборки изображающая точка остается вблизи прежнего стационарного состояния. Cборка локализуема, как и катастрофа «бабочка» с четной коразмерностью. Бифуркации, приводящие к возникновению незатухающих колебаний и квазистохастических режимов, мы рассмотрим в лекциях 8 и 10, соответственно.

ЛЕКЦИЯ 7 МУЛЬТИСТАЦИОНАР НЫЕ СИСТЕМЫ

Важная особенность биологических систем – переключение из одного режима функционирования в другой. Приведем простые примеры переключения процессов в живых системах: · Сон и бодрствование – это разные типы метаболизма. Переключение происходит периодически и синхронизируется геофизическим ритмом. · У большинства насекомых один и тот же организм может существовать в виде гусеницы, куколки, бабочки. Переключение происходит последовательно в соответствии с генетической программой. · Дифференцировка тканей – клетки получаются путем деления из одного типа клеток, но впоследствии каждая выполняет свои функции. На фазовой плоскости триггерной системе в простейшем случае соответствует два или несколько устойчивых стационарных решений, разделенных сепаратрисами. Напомним, что все особые точки (устойчивые и седло) лежат на пересечении главных изоклин – изоклин вертикальных и горизонтальных касательных (см. Лекция 4). На рис. 7.1 представлен относительно простой фазовый портрет триггерной системы, описывающей явление конкуренции двух одинаковых видов:

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

Соответствующая система уравнений имеет вид: (7.1) Такая система имеет четыре стационарных решения: 1. x 1 =x 2 =0 – неустойчивый узел; 2. – седло; 3. – устойчивый узел; 4. – устойчивый узел.

Система описывается уравнениями: (7.2) В такой системе три стационарных состояния. Состояние x=y=0 – седло. Два других стационарных состояния – устойчивые фокусы. Вблизи этих стационарных состояний траектории представляют собой закручивающиеся спирали. Вдали от стационарных состояний области притяжения имеют слоистую структуру. Толщина слоев уменьшается при уменьшении параметра a. Как видно из приведенных выше примеров, в триггерных системах и поведение во времени и стационарное решение зависит не только от параметров, но и от начальных условий. Способы переключения триггера Слово триггер означает переключатель. Встает вопрос, как можно переключить триггер из одного в другое стационарное состояние?

Рассмотрим фазовый портрет системы, обладающей двумя устойчивыми стационарными состояниями (рис. 7.3). Здесь a,c – устойчивые стационарные состояния, b – седло. Рис Триггерная система. Жирными линиями показаны главные изоклины. Пунктирной линией обозначена сепаратриса, отделяющая области влияния двух устойчивых стационарных состояний a и с Двойная стрелка показывает процесс силового переключения триггера. Если начальное положение изображающей точки расположено левее сепаратрисы седла (пунктирная линия), система находится в области притяжения особой точки a и со временем стремится к этому устойчивому стационарному состоянию. Из точек, лежащих правее сепаратрисы, система будет двигаться к особой точке c. Рассмотрим возможные способы переключения системы из режима a в режим c.

1. Силовое переключение. Можно изменить значения концентраций (например, добавить определенное количество вещества x 1, так что система «перепрыгнет» через сепаратрису, например в некоторую точку c 1, которая находится по правую сторону сепаратрисы в области влияния устойчивого стационарного состояния с, к которому система перейдет сама с течением времени. На фазовом портрете рис. 7.3 силовое (специфическое) переключение показано двойной стрелкой. Кинетика переменных при таком переключении показана на рис. 7.4.

Рис Поведение переменных во времени при силовом переключении после добавления в систему вещества x в количестве, достаточном для переключения системы из режима a в режим c (смотри рис. 7.3).

2. Параметрическое переключение. Другой – неспецифический способ переключения показан на рис. 7.5, 7.6. При таком способе переключения непосредственному воздействию подвергаются не переменные, а параметры системы. Это может быть достигнуто разными способами, например, изменением скорости поступления субстрата, температуры, рН. Сущность такого способа переключения состоит в использовании зависимости фазового портрета системы от некоторого управляющего параметра. Пусть с изменением этого параметра фазовый портрет претерпевает последовательность превращений, показанных на рис. 7.5 (а – г). На стадии (в) устойчивый узел (а) и седло (b) сливаются в одну полуустойчивую точку седло узел. На стадии (с) в системе остается лишь одно устойчивое стационарное состояние, к которому и сходятся все фазовые траектории.

Рис Параметрическое переключение триггера. Последовательные стадии трансформации фазового портрета. Стрелками обозначено направление фазовых траекторий. Рис Параметрическое переключение триггера. Последовательные стадии трансформации фазового портрета. Стрелками обозначено направление фазовых траекторий.

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

Модели отбора 1000 тыс. летКАЙНОЗОЙЭволюция человека МЕЗОЗОЙПоявление млекопитающих ПАЛЕОЗОЙПервые многоклеточные 2000 – 3000 тыс. летПРОТЕОЗОЙБиологическая эволюция 4000 тыс. летАРХЕЙМикроископаемые 5000 тыс. лет Образование земли 6000 тыс. лет Возникновение солнечной системы Как мы видели выше, в триггерной системе изображающая точка выбирает (в зависимости от параметров и начальных условий) стационарный режим функционирования. Триггерные модели могут быть использованы при описании процесса отбора, и потому применимы к процессам эволюции. Геохронологическая таблица

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

Структурирование в пространстве (тип 1) обычно предшествует конкурентному отбору (тип 2). До образования изотопов химических элементов (водород и гелий) должны были возникнуть "сгустки материи" зародыши галактик и звезд. До возникновения макромолекул должна была образоваться планетная структура и атмосфера Земли, доступная для солнечных лучей. Человеческие языки возникают в замкнутых коллективах и проч. Наоборот, процессы отбора ведут к возможности образования новых, более сложных структур. На базе разнообразия макромолекул идет становление живого организма. Биогеоценозы формируются из разных видов живых существ и т.д.

Главный вопрос эволюции: «как появилось свойство авторепродукции?» включает в себя несколько вопросов. Вот основные из них: · Как возникли комплексы белка и полинуклеотидов? · Как образовался единый генетический код? (т.е. соответствие между последовательностями нуклеотидов в ДНК и аминокислот в белках?) · Почему именно эти три нуклеотида кодируют данную аминокислоту (кодон). Действительно, существующий генетический код не связан с физико-химическими свойствами аминокислот и кодонов. Число равноправных кодов 20! а реализован только один. Вероятность случайного возникновения именно существующего кода крайне мала. На вопрос «Как же произошел отбор»? возможно несколько ответов: Кастлер: начальный код возник случайно, другие комбинации не успели возникнуть. Эйген: возникло несколько разных кодов, но отобрались наилучшие. Чернавский: произошел отбор одного из равноправных.

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

Отбор одного из равноправных Общая модель такого отбора имеет вид: (7.3) Здесь a - эффективный коэффициент репродукции, - вероятность гибели в результате встречи. Пусть N= 2, X 1 = x, X 2 = y. Система уравнений имеет вид: (7.4) Стационарные решения находятся из алгебраических уравнений, полученных приравниванием правых частей нулю. (7.5) Система имеет два стационарных решения: (7.6)

Для второго – нетривиального симметричного стационарного состояния характеристический определитель системы имеет вид: (7.7) Характеристическое уравнение: или Выражения для характеристических чисел находятся из уравнения: (7.8) Корни положительны и разных знаков. Это означает, что симметричное стационарное состояние представляет собой седло. Аналогичный анализ показывает, что нулевая особая точка представляет собой неустойчивый узел. Изоклины горизонтальных касательных: y=0 – ось абсцисс и вертикальная прямая x=a/ ; изоклины вертикальных касательных: x=0 – ось ординат – и горизонтальная прямая y=a/. Все траектории уходят на бесконечность, так как самоограничение роста популяции в данной модели не учитывается.

Рис Фазовый портрет системы (7.4), описывающей отбор одного из двух равноправных видов в отсутствие ограничений роста. а (начало координат) – неустойчивый узел, b – седло

Биологический смысл модели. Модель (7.4) демонстрирует принципиальную возможность отбора в системе равноправных видов, где симметричное состояние сосуществования является неустойчивым. Вот один из примеров такой системы. Известно, что сахара и аминокислоты являются оптически активными соединениями, причем сахара – левовращающие плоскость поляризации света, аминокислоты – правовращающие. Противоположные изомеры не только не встречаются в живых организмах и не усваиваются ими, но являются ядами. В этом заключается одна из сложностей искусственного синтеза. Ясно, что «зеркальные» организмы не лучше и не хуже. В неживой природе распространены рацемические смеси, содержащие равное количество зеркальных изомеров. То же – при небиологическим синтезе. По-видимому, и первичный бульон был рацемической смесью. Рассмотренная модель описывает выживание одних и уничтожение других. Условие, которое обеспечивает при этом отбор одного вида, заключается в том, что при встрече они взаимно отравляются и гибнут. Причина отбора здесь – не преимущество одного из видов, а их взаимный антагонизм.

Однако модель (7.4) не может описывать реальную систему, так как описывает неограниченный рост биомассы с течением времени. Этот недостаток может быть исправлен несколькими способами. Один их них – введение самоограничения численности вида в виде ферхюльстовских членов. Тогда мы придем к модели (7.1). Другой способ – ввести в модель переменную, описывающую поступающий в систему с определенной скоростью питательный ресурс, общий для обоих видов. Учтем ограниченность питательных ресурсов. Пусть S -лимитирующий субстрат (световая энергия, минеральное питание и т.п.). Сам субстрат не является оптически активным, но преобразуется в оптически активные продукты. Выразим скорость роста каждой популяции a через S в соответствии с формулой Моно (7.9.). График этой функции приведет на рис (7.9) Пусть интенсивность притока субстрата. Расход субстрата пропорционален поглощению его организмами, т.е. сумме их концентраций. Уравнение для скорости изменения концентрации субстрата во времени имеет вид: (7.10)

Здесь а > 1 - экономический коэффициент - указывает, сколько субстрата идет на образование единицы биомассы. Уравнения для концентраций объектов типа x и y: (7.11) Введем безразмерные переменные: (7.12)

Система в безразмерном виде: (7.13) Пусть процессы поглощения субстрата существенно более быстрые, чем процессы репродукции. В этом случае может быть использован метод квазистационарных концентраций (лекция 6), и дифференциальное уравнение для быстрой переменной z (S) – концентрации субстрата – заменено алгебраическим. Тогда субстрат на интересующих нас временах достигнет квазистационарной концентрации: dz/dt=0. Отсюда (7.14) В итоге получается система двух безразмерных уравнений (7.15)

Построим фазовый портрет системы. (Рис.7.8) Изоклины вертикальных касательных: x=0 (ось ординат) и кривая ; Изоклины горизонтальных касательных: y=0 (ось абсцисс) и 0 – (1+x)(x+y) = 0 или: – гипербола. Переменные x и y симметричны, поэтому изоклина вертикальных касательных симметрична изоклине горизонтальных касательных. Система имеет четыре особые точки: 1) x=0, y=0 - неустойчивый узел; 2) x=0, y= 0 - устойчивый узел; 3) x= 0, y=0 - устойчивый узел; 4) и, наконец, симметричную точку - седло (7.16) В такой системе выживет один из видов: x или y. Его стационарная концентрация определяется скоростью притока субстрата и экономическим коэффициентом. Как и в предыдущей системе (7.4) здесь причина отбора – неустойчивость симметричного состояния.

Рис Фазовый портрет системы 7.15, описывающей отбор одного из двух равноправных видов когда субстрат поступает в систему с постоянной скоростью. а (начало координат) – неустойчивый узел, b – седло, c, d –устойчивые узлы.

Генетический триггер Жакоба и Моно Рассмотрим модель биохимической регуляции белкового синтеза, предложенную Жакобом и Моно в 1964 г. и математически разработанную Д.С. Чернавским в 1967 г. Эта модель показывает принципиальные возможности триггерных систем. Она легла в основу целой серии более подробных и конкретных моделей. Подробный вывод модели описан в книге Романовский Ю.М., Степанова Н.В., Чернавский Д.С. «Математическое моделирование в биофизике. М., Схема взаимной регуляции двух систем синтеза ферментов изображена на рис Ген-регулятор каждой системы синтезирует неактивный репрессор. Этот репрессор, соединяясь с продуктом противоположной системы, образует активный комплекс, который обратимо реагируя с участком структурного гена, называемым опероном, блокирует синтез mRNK. Таким образом, продукт первой системы P 1 является корепрессором второй системы, а продукт второй системы P 2 - корепрессором первой. При этом в процессе корепрессии могут принимать участие одна, две и более молекул продукта.

Рис Схема синтеза двух ферментов по Жакобу и Моно

Очевидно, что при таком характере взаимодействий при интенсивной работе первой системы вторая заблокирована, и наоборот. Простейшая система уравнений, описывающая такой тип взаимодействий, имеет вид (7.17) Здесь P 1, P 2 - концентрации продуктов, величины A 1, A 2, B 1, B 2, выражаются через параметры своих систем. Показатель степени m показывает, сколько молекул активного репрессора (соединений молекул продукта с молекулами неактивного репрессора, который предполагается в избытке) соединяются с опероном для блокировки синтеза mRNK. Введем безразмерные переменные: (7.18)

Опустив штрих у времени, перепишем систему в безразмерном виде: (7.19) Исследование системы (7.19) показало, что при m = 0 фазовый портрет имеет одну устойчивую особую точку в первом квадранте фазовой плоскости (рис. 7.10а) и не может описывать процессов переключения в системе. При m 2 и определенных значениях отношения L 1 / L 2 > система приобретает триггерные свойства. На фазовой плоскости такая система имеет две устойчивые особые точки, между которыми расположено седло (рис б). Значение параметра является бифуркационным, причем бифуркация имеет триггерный характер (образуется седло). Отношение L 1 /L 2 служит управляющим параметром, изменение значения L 1 /L 2 > которого может привести к смене стационарного режима в системе, как это было описано выше при рассмотрении параметрического способа переключения системы. Величина параметров L 1, L 2 зависит от многих биохимических характеристик: скорости снабжения субстратами, активности ферментов, времени жизни ферментов, mRNK и продуктов и проч.

Рис Главные изоклины на фазовой плоскости системы При m = 1. система имеет единственное устойчивое стационарное состояние (а). При m = 2 в системе три стационарных состояния, два из которых (а и с) – устойчивые узлы, а третье (b) – седло.

Фазовый портрет триггерной системы 7.19 приведен на рис Таким образом, триггерные модели могут описывать процессы отбора и дифференцировки. Подобные механизмы взаимодействия в распределенной системе (при учете пространственной неоднородности и процессов переноса) могут описывать процессы морфогенеза (формообразования). Эту модель – «распределенный генетический триггер» мы рассмотрим во второй части лекций.

ЛЕКЦИЯ 8 КОЛЕБАНИЯ В БИОЛОГИЧЕСКИХ СИСТЕМАХ

Для биологических систем характерно периодическое изменение различных характеристик. Период этих колебаний может быть связан с периодическими изменениями условий жизни на Земле смена времен года, смена дня и ночи. Существуют и другие геофизические ритмы –солнечные, лунные, связанные с периодами атмосферных явлений. Геофизические и биологические ритмы сопоставлены на рис Но многие периодические процессы имеют частоту изменения, не связанную очевидным образом с внешними геокосмическими циклами. Это так называемые «биологические часы» различной природы, начиная от колебаний биомакромолекул, биохимических колебаний, вплоть до популяционных волн. Внутриклеточные колебания задают эндогенные биологические ритмы, которые свойственны всем живым системам. Именно они определяют периодичность деления клеток, отмеряют время рождения и смерти живых организмов. Модели колебательных систем используются в ферментативном катализе, теории иммунитета, в теории трансмембранного ионного переноса, микробиологии и биотехнологии. С некоторыми из типов периодических движений мы уже имели дело при рассмотрении особых точек типа центр и затухающих или нарастающих колебаний в случае устойчивого и неустойчивого фокуса. Однако «биологические часы» имеют свойство, отличающее их от рассмотренных типов колебаний неизменность во времени периода и амплитуды таких колебаний, означающую стационарность и устойчивость колебательного режима.

Рис Космофизические, геофизические и биологические ритмы. Справа – шкала периодов, слева – шкала частот

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

Предельный цикл. В фазовом пространстве такому типу поведения соответствует притягивающее множество (аттрактор), называемое предельным циклом. Предельный цикл есть изолированная замкнутая кривая на фазовой плоскости, к которой в пределе при t стремятся все интегральные кривые. Предельный цикл представляет стационарный режим с определенной амплитудой, не зависящий от начальных условий, а определяющийся только организацией системы. Существование предельного цикла на фазовой плоскости есть основной признак автоколебательной системы. Очевидно, что при автоколебательном процессе фаза колебаний может быть любой. Остановимся на общих характеристиках автоколебательных систем. Рассмотрим систему уравнений общего вида: (8.1) Если T (T > 0) наименьшее число, для которого при всяком t то изменение переменных x = x(t), y = y(t) называется периодическим изменением с периодом T.

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

Например, для системы (8.2) Траектория является предельным циклом. Его параметрические уравнения будут: а уравнения всех других фазовых траекторий запишутся в виде: Значениям постоянной интегрирования С > 0 соответствуют фазовые траектории, накручивающиеся на предельный цикл изнутри (при t ), а значениям – 1 < C < 0 траектории, накручивающиеся снаружи.

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

Неустойчивые предельные циклы, подобные изображенному на рис. 8.2 б, такие, что все траектории с одной стороны (например, изнутри) приближаются к ним, а с другой стороны (например, извне) удаляются от них при t, называют «полуустойчивыми» или двойными. Последнее название связано с тем, что обычно такие циклы при подходящем изменении параметра системы расщепляются на два, один из которых устойчив, а другой неустойчив. А.М. Ляпунов показал, что для исследования устойчивости периодического движения x = (t), y = (t) можно идти по пути линеаризации уравнений, подобно тому, как мы это делали при исследовании устойчивости состояний равновесия.

Если положить подставить эти выражения в уравнения (8.1), разложить правые части этих уравнений функции в ряды по степеням и и отбросить нелинейные члены, то мы получим линейные уравнения (уравнения первого приближения) для координат возмущения и : Коэффициенты в правой части:

Это система линейных дифференциальных уравнений с периодическими коэффициентами периода T, поскольку a, b, c, d суть функции от, периодических функций времени с периодом T. Общий вид ее решения Здесь некоторые периодические функции с периодом T. От показателей и которые носят название «характеристических показателей», зависят свойства решений для отклонений от стационарного периодического решения и. А именно, знаки их действительных частей определяют, являются ли эти решения нарастающими или затухающими. Можно показать, что в силу автономности исходной системы (8.1) один из характеристических показателей равен нулю, а другой равен h.

где x = (t), y = (t) любое периодическое решение, соответствующее рассматриваемому предельному циклу, T период решения. Таким образом, устойчивость предельного цикла (и устойчивость в смысле Ляпунова соответствующих периодических движений) определяется знаком характеристического показателя. Предельный цикл устойчив, если h 0. Если же h = 0, уравнения первого приближения не решают вопроса об устойчивости периодического движения. Для нахождения предельных циклов не существует таких простых аналитических методов, как для нахождения стационарных точек и исследования их устойчивости. Однако, исследование фазовой плоскости системы позволяет ответить на вопрос, есть в данной системе предельный цикл, или нет. Сформулируем несколько теорем, определяющих наличие предельного цикла по топологическому строению фазовой плоскости. Они могут быть полезны как при аналитическом, так и при компьютерном анализе системы.

Теорема 1. Пусть на фазовой плоскости существует область, из которой фазовые траектории не выходят, и в которой нет положений равновесия (особых точек). Тогда в этой области обязательно существует предельный цикл, причем все остальные траектории обязательно наматываются на него. На рис изображена такая область G, из которой фазовые траектории не выходят. Это означает, что фазовые траектории либо входят, пересекая границу, внутрь области, либо сама граница является траекторией. Легко видеть, что такая область не может быть односвязной. Поскольку траектория наматывается на предельный цикл изнутри, это означает, что внутри этого предельного цикла на фазовой плоскости существует либо неустойчивая особая точка, либо неустойчивый предельный цикл, очевидно, не принадлежащие рассматриваемой области G. Таким образом, если найти на фазовой плоскости такую двусвязную область, что направления фазовых траекторий на всей границе обращены внутрь этой области, то можно утверждать, что внутри этой области имеется предельный цикл.

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

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

Рождение предельного цикла. Бифуркация Андронова-Хопфа. Существование предельных циклов возможно лишь в системе типа (8.1), правые части которой представлены нелинейными функциями. На бифуркационной диаграмме 4.11 мы видели, что при пересечении оси абсцисс происходит смена устойчивости фокуса. Нулевым значениям действительной части характеристических чисел (ляпуновских показателей) соответствует особая точка типа центр. В нелинейной системе, где возникает неустойчивый фокус, при этом возможно рождение предельного цикла. Такой переход легко проследить в «модельной» системе: (8.3) Схематически возникновение предельного цикла в системе (8.3) изображено на фазопараметрической диаграмме на рис. 8.6.

Рис Закритическая (суперкритическая) бифуркация Андронова-Хопфа. Мягкое возбуждение. При с>0 возникают автоколебания, амплитуда которых растет с увеличением с.

Выполнению условия Re 1,2 = 0, причем Im 1,2 0, соответствует бифуркация Андронова–Хопфа или бифуркация рождения (исчезновения) предельного цикла. Бифуркация впервые была исследована А.А. Андроновым для случая N = 2 и обобщена Е. Хопфом на системы с произвольной размерностью. (Андронов А.А., Витт А.А., Хайкин С.Э. Теория колебаний. М., Наука, 1981; Hopf E., 1942) Существуют два типа бифуркации Андронова–Хопфа. Только что мы рассмотрели суперкритическую бифуркацию (мягкое возбуждение автоколебаний). Возможна также субкритическая бифуркация (жесткое возбуждение автоколебаний). В этом случае при бифуркационном значении параметра устойчивый фокус теряет устойчивость из-за «влипания» в него неустойчивого предельного цикла (рис. 8.7). Фокус становится неустойчивым, а аттрактором при этом может стать предельный цикл большой амплитуды.

«Модельной» системой (см. лекция 6), описывающей рождение предельного цикла при жестком возбуждении, является система: (8.4) Приравняв правую часть первого уравнения нулю, получим стационарные значения r: Ветвь r = 0 устойчива при c 0.

При с > –1 стационарное решение устойчивый предельный цикл. При –1 < с < 0 стационарное решение неустойчивый предельный цикл. Рассмотрим, что произойдет, если двигаться по параметру с, начиная с отрицательных значений (Рис.8.8).

Первоначально имеется единственное устойчивое стационарное состояние r = 0, колебаний нет. При c > –1 существует также устойчивый предельный цикл, но система не покидает своего устойчивого стационарного состояния. Однако после того как с становится положительным, стационарное состояние становится неустойчивым, и происходит резкий скачок к устойчивому предельному циклу. В системе начинаются колебания сразу большой амплитуды. Если двигаться от положительных значений с к отрицательным, колебания большой амплитуды сохраняются до тех пор, пока с не станет меньше –1, а затем внезапно исчезнут. Таким образом при –1 < с < 0 могут существовать два различных типа поведения. Какой из них реализуется, зависит от предыстории системы. Такой феномен называется эффектом гистерезиса.

При увеличении параметра с и его переходе через ноль скачком возникают устойчивые автоколебания конечной амплитуды и частоты. Для промежуточных значений параметра с существуют два типа устойчивого поведения (два аттрактора) устойчивое стационарное состояние и устойчивый предельный цикл. Винфри (Winfree A.T.) назвал области, в которых возможны два режима: устойчивая точка покоя и предельный цикл, ­ черной дырой (рис. 8.8 б). В этой области параметров можно так приложить возмущение к колебательной системе, что она попадет в область притяжения точки покоя, что приведет к прекращению колебаний. В частности, это показано для уравнений Ходжкина– Хаксли, моделирующих проведение нервного импульса (см. ниже).

Брюсселятор. Простейшим классическим примером существования автоколебаний в системе химических реакций является тримолекулярная модель «Брюсселятор», предложенная в Брюсселе Пригожиным и Лефевром (1965). Основной целью при изучении этой модели было установление качественных типов поведения, совместимых с фундаментальными законами химической и биологической кинетики. В этом смысле блюсселятор играет роль базовой модели, такую же как гармонический осциллятор в физике, или модели Вольтерра в динамике популяций. Во 2-й части лекций мы остановимся на пространственно временных свойствах распределенной системы, локальным элементом которой является брюсселятор. Здесь мы рассмотрим свойства брюсселятора как автоколебательной системы. Брюсселятор содержит простейшую реализацию кубической нелинейности посредством химической реакции 2X + Y 3X(8.5)

Хотя тримолекулярная стадия в химической кинетике не столь распространена, как бимолекулярные процессы, выражения для скорости ряда биохимических реакций в определенных случаях можно свести к кубическому виду. В качестве примера приведем следующую последовательность ферментативных реакций: X + E EX EX + Y XY EXY + X EX 2 Y Здесь предполагается что фермент E имеет по крайней мере три каталитических центра, способных одновременно фиксировать две молекулы X и одну молекулу Y. Если образующиеся комплексы распадаются с достаточно большой скоростью, а ферменты присутствуют в небольших количествах, легко показать, что всю последовательность реакций можно свести к одной стадии, дающей нелинейный член типа X 2 Y в выражении для скорости реакции. Брюсселятор представляет собой следующую схему гипотетических химических реакций:

Здесь А, В исходные вещества, C, R продукты, X, Y промежуточные вещества. Пусть конечные продукты С и R немедленно удаляются из реакционного пространства. Это означает, что обратные константы k 3 = k 4 = 0. Если субстрат A находится в избытке, k 1 = 0. Предположим также, что k 2 = 0. Значения остальных констант положим равными единице. Тогда схема реакций 9.2 (в случае точечной системы) описывается системой уравнений: (8.6) Модель (8.5) имеет одну особую точку с координатами: (8.7)

Исследуем стационарное решение (8.6) на устойчивость по методу Ляпунова. Введем переменные, характеризующие отклонения от особой точки: Линеаризованная система имеет вид: Характеристическое уравнение или 2 + (A B) + A 2 = 0 имеет корни: (8.7)

Напомним, что особая точка является устойчивой, если действительные части корней характеристического уравнения отрицательны. Из выражения (8.7) видно, что при B < 1 + A 2 особая точка (8.6) устойчива. Если же B > 1 + A 2 особая точка становится неустойчивой, и у системы (8.5) появляется устойчивый предельный цикл. Значение B = 1 + A 2 является бифуркационным. Если величина B лишь немного превосходит бифуркационный порог, автоколебания в системе носят квазигармонический характер. Таким образом, брюсселятор при выполнении условияявляется автоколебательной системой. Фазовый портрет брюсселятора при разных значениях параметров изображен на рис. 8.9 Здесь мы приведем краткий обзор нескольких «успешных» моделей колебательных биологических процессов. Более подробно некоторые колебательные процессы будут рассмотрены в лекциях 9, 11, 12.

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

Известно, что в условиях смены дня и ночи интенсивность фотосинтеза, то есть скорость выделения кислорода и поглощения СО 2 изменяется периодически (рис.8.10 а). Если растение поместить в условия непрерывной освещенности, то периодичность в интенсивности фотосинтеза с периодом несколько часов сохраняется достаточно длительное время. По-видимому, растение имеет свой внутренний ритм, синхронизованный с периодическим внешним воздействием. Напомним, что в процесс фотосинтеза входят световой и темновой циклы химических реакций. Первый включает поглощение энергии квантов света и, через ряд промежуточных стадий, приводит к образованию высокоэнергетических восстановленных химических соединений и богатых энергией молекул АТФ. Эти вещества употребляются в темновом цикле (цикле Кальвина), в котором свет непосредственно не участвует. Здесь происходит восстановление углекислоты СО 2 с помощью веществ, богатых энергией, и доноров водорода, полученных в световом цикле, и превращение ее в углеводы фруктозу и глюкозу (рис. 8.11).

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

Предполагалось, что прибыль тяжелых сахаров с 6 может осуществляться за счет соединения двух легких с 3. Их убыль, так же как и убыль тяжелых сахаров, происходит в результате бимолекулярного взаимодействия тяжелых и легких сахаров. Имеет место также приток продукта с 3 в сферу реакции за счет биохимически сходных процессов (гликолиза, дыхания). Эти предположения приводят к системе уравнений:(8.9)

Переменные представляют собой нормированные концентрации легких (x) и тяжелых (y) сахаров. В положительном квадранте имеется одно состояние равновесия с координатами (1,1). Изоклины горизонтальных касательных определяются из уравнения а изоклины вертикальных касательных – из уравнения На рис изображены фазовые портреты системы. При это устойчивый фокус (рис а). При неустойчивый фокус, окруженный предельным циклом (рис б).

Колебания в гликолизе. Классическим примером колебательной биохимической реакции является гликолиз. В процессе гликолиза осуществляется распад глюкозы и других сахаров, при этом соединения, содержащие шесть молекул углерода, превращаются в трикарбоновые кислоты, включающие три молекулы углерода. За счет избытка свободной энергии в процессе гликолиза на одну молекулу шестиуглеродного сахара образуются две молекулы АТФ. Основную роль в генерации наблюдаемых колебаний концентраций компонентов реакции: фруктозо-6-фосфата, фруктозо-1,6- фосфата и восстановленного НАД (никотинаминадениндинуклеотид) играет ключевой фермент гликолитического пути фосфофруктокиназа (ФФК). Упрощенная схема реакций представлена на рис. 8.13:

На схеме [Гл] - глюкоза, Ф6Ф - фруктозо-6-фосфат субстрат ключевой реакции, ФДФ продукт этой реакции, который является субстратом в следующей стадии. Обе реакции катализируются ферментами. В безразмерных переменных система описывающих реакции уравнений может быть записана в виде: Здесь зависимости скоростей реакций от переменных записаны в форме Михаэлиса-Ментен (Моно), как это было представлено в уравнении (8.10).Если выполняются условия: K mx >>x, K my >>y, можно произвести замену переменных

Опустив штрихи у новых переменных, получим систему в безразмерном виде : (8.10) где Кинетика изменений переменных и фазовые портреты системы (8.10) при разных значениях параметров представлены на рис. 8.14

Рис Модель гликолиза (8.10). Кинетика изменений концентраций фруктозо 6 фосфата (х) и фруктозодифосфата (у) (справа) и фазовый портрет системы (слева) при разных значениях параметров системы, в – колебания с постоянной амплитудой и фазой, близкие к гармоническим (предельный цикл на фазовой плоскости), α = 6; r = 0.2. г – релаксационные колебания с постоянной амплитудой и фазой, (предельный цикл почти треугольной формы на фазовой плоскости ) α = 8; r = 0.5 Интересно, что колебательные реакции в системе гликолиза были сначала предсказаны на математической модели (Higgins, 1964), и лишь после этого зарегистрированы экспериментально с помощью метода дифференциальной спектрофотометрии в лаборатории Б. Чанса (1966).

Внутриклеточные колебания кальция. Во многих типах живых клеток наблюдаются колебания внутриклеточной концентрации кальция, период которых может варьировать от 0,5 до 10 мин. Простейшая схема процессов, приводящих к гармонально обусловленным колебаниям кальция, основой которых служит кальций индуцированный выходакальция из клетки, приведена на рис Такие колебания впервые наблюдались Эндо с соавторами (1970) на клетках скелетных мышц, Фабиато (1975) на клетках саркоплазматического ретикулума сердца быка, и позднее многими другими исследователями. Схема и модель процессов, предложена и описана в [Dupont and Goldbetter (1989, 1994)]. Рассматриваются приток и отток кальция в клетку через плазматическую мембрану (константы скоростей v 1 и v 2, соответственно); гормонально активируемое освобождение кальция из пула (скорость v 3 ); активный транспорт цитозольного кальция в пул, (v 4 ), освобождение кальция из пула, активируемое цитозольным кальцием (v 5 ); свободный отток кальция из пула в цитозоль (v 6 ).

Модель состоит из двух дифференциальных уравнений (8.11) Здесь S 1 концентрация кальция в цитозоле, S 1 концентрация кальция в гормонально чувствительном пуле. Выражения для величин скоростей были предложены в Simogyi, Stuckin (1991): (8.12) Модель предсказывает колебания концентрации кальция во времени, по форме близкие к экспериментальным (рис. 8.16).

Клеточные циклы. В процессе жизненного цикла клетка удваивает свое содержимое и делится на две. В организме млекопитающего для поддержания жизни производятся ежесекундно миллионы новых клеток. Нарушение регуляции пролиферации клеток проявляется как онкологическое заболевание. Этим вызван большой интерес к изучению и моделированию механизмов регуляции клеточного деления. Схема клеточного цикла изображена на рис Клеточный цикл состоит из двух периодов: митоз (М фаза) включает разделение предварительно удвоенного ядерного материала, деление ядра и деление самой клетки цитокинез и занимает около часа. Значительно более длительный период между двумя митозами занимает интерфаза, включающая стадию роста G 1, фазу репликации ДНК (S), фазу подготовки к делению G 2. Клеточный цикл (рис.8.17) регулируется генами и белками-ферментами двух основных классов.

Циклин-зависимые протеин киназы (Cdk) индуцируют последовательность процессов путем фосфорилирования отдельных белков. Циклины, которые синтезируются и деградируют при каждом новом цикле деления, связываются с молекулами Cdk и контролируют их способность к фосфорилированию, без циклина Cdk не активны. Количество этих молекул регуляторов различно в разного вида клетках. В делении дрожжевой клетки основные роли играют один Cdk и девять циклинов, которые образуют девять разных циклин Cdk комплексов. У гораздо более сложно организованных млекопитающих изучено шесть Cdk и полтора десятка циклинов. Контроль выхода клетки из G 1, и G 2 фаз осуществляют промотор фактор S фазы (SPF) и промотор фактор M фазы (MPF), представляющие собой гетеродимеры. Cуществует особая контрольная точка клеточного цикла (Start), с которой заканчивается рост (G 1 фаза) и начинается процесс синтеза ДНК. Простая модель процесса предложена Тайсоном (Tyson, 1995). Постулируется существование фактора транскрипции SBF, который может быть в активной S a и пассивной S i форме.

Он переходит в активную форму под действием циклина Cln (N) и Start киназы (Cdc28-Cln3) (A) и инактивируется другим веществом (Е). Циклин продуцируется путем активации SBF и деградирует. SBF активируется Chu и Start-киназой и инактивируется фосфатазой. Безразмерная модель процессов имеет вид: (8.13) Модель имеет одно или три стационарных решения (два устойчивых) в зависимости от значений параметров, и при увеличении параметра (в процессе роста клетки) описывает переключение системы из G 1 в S фазу. Добавление двух уравнений сходного вида позволяет описать также переключение из G 2 в фазу митоза М. Полная модель, учитывающая и другие регуляторные ферменты в фосфорилированной и дефосфорилированной форме содержит 9 нелинейных уравнений (Novak, Tyson 1993) и хорошо описывает кинетику деления ооцитов Xenopus. При соответствующем подборе параметров она применима к описанию деления других типов клеток.

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

ЛЕКЦИЯ 9 МОДЕЛИ ВЗАИМОДЕЙСТВИЯ ДВУХ ВИДОВ

Основателем современной математической теории популяций справедливо считается итальянский математик Вито Вольтерра, разработавший математическую теорию биологических сообществ, аппаратом которой служат дифференциальные и интегро- дифференциальные уравнения. (Vito Volterra. Lecons sur la Theorie Mathematique de la Lutte pour la Vie. Paris, 1931). В последующие десятилетия популяционная динамика развивалась, в основном, в русле высказанных в этой книге идей. Русский перевод книги Вольтерра вышел в 1976 г. под названием: «Математическая теория борьбы за существование» с послесловием Ю.М. Свирежева, в котором рассматривается история развития математической экологии в период гг. Книга Вольтерра написана так, как пишут книги по математике. В ней сначала сформулированы некоторые предположения о математических объектах, которые предполагается изучать, а затем проводится математическое исследование свойств этих объектов. Системы, изученные Вольтерра, состоят их двух или нескольких видов. В отдельных случаях рассматривается запас используемой пищи. В основу уравнений, описывающих взаимодействие этих видов, положены следующие представления.

Гипотезы Вольтерра 1. Пища либо имеется в неограниченном количестве, либо ее поступление с течением времени жестко регламентировано. 2. Особи каждого вида отмирают так, что в единицу времени погибает постоянная доля существующих особей. 3. Хищные виды поедают жертв, причем в единицу времени количество съеденных жертв всегда пропорционально вероятности встречи особей этих двух видов, т.е. произведению количества хищников на количество жертв. 4. Если имеется пища в ограниченном количестве и несколько видов, которые способны ее потреблять, то доля пищи, потребляемой видом в единицу времени, пропорциональна количеству особей этого вида, взятому с некоторым коэффициентом, зависящим от вида (модели межвидовой конкуренции). 5. Если вид питается пищей, имеющейся в неограниченном количестве, прирост численности вида в единицу времени пропорционален численности вида. 6. Если вид питается пищей, имеющейся в ограниченном количестве, то его размножение регулируется скоростью потребления пищи, т.е. за единицу времени прирост пропорционален количеству съеденной пищи.

Аналогии с химической кинетикой Эти гипотезы имеют близкие параллели с химической кинетикой. В уравнениях популяционной динамики, как и в уравнениях химической кинетики, используется принцип соударений, когда скорость реакции пропорциональна произведению концентраций реагирующих компонентов. Действительно, согласно гипотезам Вольтерра, скорость процесса отмирания каждого вида пропорциональна численности вида. В химической кинетике это соответствует мономолекулярной реакции распада некоторого вещества, а в математической модели – отрицательным линейным членам в правых частях уравнений. Согласно представлениям химической кинетики, скорость бимолекулярной реакции взаимодействия двух веществ пропорциональна вероятности столкновения этих веществ, т.е. произведению их концентрации. Точно так же, в соответствии с гипотезами Вольтерра, скорость размножения хищников (гибели жертв) пропорциональна вероятности встреч особей хищника и жертвы, т.е. произведению их численностей. И в том и в другом случае в модельной системе появляются билинейные члены в правых частях соответствующих уравнений.

Наконец, линейные положительные члены в правых частях уравнений Вольтерра, отвечающие росту популяций в неограниченных условиях, соответствуют автокаталитическим членам химических реакций. Такое сходство уравнений в химических и экологических моделях позволяет применить для математического моделирования кинетики популяций те же методы исследований, что и для систем химических реакций. В соответствии с гипотезами Вольтерра взаимодействие двух видов, численности которых x 1 и x 2, могут быть описаны уравнениями: (9.1) Здесь параметры a i константы собственной скорости роста видов, c i константы самоограничения численности (внутривидовой конкуренции), b ij константы взаимодействия видов, (i, j=1,2). Знаки этих коэффициентов определяют тип взаимодействия.

В биологической литературе обычно классифицируют взаимодействия по участвующим в них механизмам. Разнообразие здесь огромно: различные трофические взаимодействия, химические взаимодействия, существующие между бактериями и планктонными водорослями, взаимодействия грибов с другими организмами, сукцессии растительных организмов, связанные в частности, с конкуренцией за солнечный свет и с эволюцией почв и т.д. Такая классификация кажется необозримой. E. Одум, учитывая предложенные В.Вольтерра модели, предложил классификацию не по механизмам, а по результатам. Согласно этой классификации, оценивать взаимоотношения следует как положительные, отрицательные или нейтральные в зависимости от того, возрастает, убывает или остается неизменной численность одного вида в присутствии другого вида. Тогда основные типы взаимодействий могут быть представлены в виде таблицы.

ТИПЫ ВЗАИМОДЕЙСТВИЯ ВИДОВ СИМБИОЗ++b 12,b 21 >0 КОММЕНСАЛИЗМ +0b 12,>0, b 21 =0 ХИЩНИК-ЖЕРТВА + b 12,>0, b 21

Уравнения КОНКУРЕНЦИИ: Как мы уже видели в лекции 6, уравнения конкуренции имеют вид: (9.2) Стационарные решения системы: (1). Начало координат, при любых параметрах системы представляет собой неустойчивый узел. (2). (9.3)

Cтационарное состояние (9.3) представляет собой седло при a 1 >b 12 / с 2 и устойчивый узел при a 1 b 21 /c 1 и устойчивый узел при a 2< b 21 /c 1 (4). (9.5) Стационарное состояние (9.5) характеризует сосуществование двух конкурирующих видов и представляет собой устойчивый узел в случае выполнения соотношения:

Отсюда следует неравенство: b 12 b 21 b 12,, c 1 > b 21. Эти неравенства показывают, что увеличение численности одного из конкурентов сильнее подавляет его собственный рост, чем рост другого конкурента. Если численность обоих видов ограничивается, частично или полностью, различными ресурсами, приведенные выше неравенства справедливы. Если же оба вида имеют совершенно одинаковые потребности, то один из них окажется более жизнеспособным и вытеснит своего конкурента.

Поведение фазовых траекторий системы дает наглядное представление о возможных исходах конкуренции. Приравняем нулю правые части уравнений системы (9.2): x 1 (a 1 –c 1 x 1 – b 12 x 2 ) = 0 (dx 1 /dt = 0), x 2 (a 2 –b 21 x 1 – c 2 x 2 ) = 0 (dx 2 /dt = 0), При этом получим уравнения для главных изоклин системы x 2 = – b 21 x 1 / c 2 +a 2 /c 2, x 2 = 0 – уравнения изоклин вертикальных касательных. x 2 = – c 1 x 1 / b 12 +a 1 /b 12, x 1 = 0 – уравнения изоклин вертикальных касательных. Точки попарного пересечения изоклин вертикальных и горизонтальных касательных систем представляют собой стационарные решения системы уравнений (9.2.), а их координаты суть стационарные численности конкурирующих видов.

Возможное расположение главных изоклин в системе (9.2) изображено на рис.9.1. Рис. 9.1 а соответствует выживанию вида x 1, рис. 9.1 б – выживанию вида x 2, рис. 9.1 в – сосуществованию видов при выполнении условия (9.6). Рисунок 9.1 г демонстрирует триггерную систему. Здесь исход конкуренции зависит от начальных условий. Ненулевое для обоих видов стационарное состояние (9.5) – неустойчивое. Это – седло, через которое проходит сепаратриса, отделяющая области выживания каждого из видов. Рис Расположение главных изоклин на фазовом портрете вольтерровской системы конкуренции двух видов (9.2) при разном соотношении параметров. Пояснения в тексте.

Для изучения конкуренции видов ставились эксперименты на самых различных организмах. Обычно выбирают два близкородственных вида и выращивают их вместе и по отдельности в строго контролируемых условиях. Через определенные промежутки времени проводят полный или выборочный учет численности популяции. Регистрируют данные по нескольким повторным экспериментам и анализируют. Исследования проводили на простейших (в частности, инфузориях), многих видах жуков рода Tribolium, дрозофиллах, пресноводных ракообразных (дафниях). Много экспериментов проводилось на микробных популяциях (см. лекция 11). В природе также проводили эксперименты, в том числе на планариях (Рейнольдс), двух видах муравьев (Понтин) и др. На рис изображены кривые роста диатомовых водорослей, использующих один и тот же ресурс (занимающих одну и ту же экологическую нишу). При выращивании в монокультуре Asterionella Formosa выходит на постоянный уровень плотности и поддержвает концентрацию ресурса (силиката) на постоянно низком уровне. Б. При выращивании в монокультуре Synedrauina ведет себя сходным образом и поддерживает концентрацию силиката на еще более низком уровне. В. При совместном культивировании (в двух повторностях) Synedrauina вытесняет Asterionella Formosa. По-видимому, Synedra выигрывает конкуренцию благодаря своей способности к более полному использованию субстрата (см. также Лекцию 11).

Рис Конкуренция у диатомовых водорослей. а при выращивании в монокультуре Asterionella Formosa выходит на постоянный уровень плотности и поддерживает концентрацию ресурса (силиката) на постоянно низком уровне. б при выращивании в монокультуре Synedrauina ведет себя сходным образом и поддерживает концентрацию силиката на еще более низком уровне. в при совместном культивировании (в двух повторностях) Synedruina вытесняет Asterionella Formosa. По-видимому, Synedra выигрывает конкуренцию благодаря своей способности к более полному использованию субстрата (см. также Лекцию 11).

Широко известны эксперименты по изучению конкуренции Г. Гаузе, продемонстрировавшие выживание одного из конкурирующих видов и позволившие ему сформулировать «закон конкурентного исключения». Закон гласит, что в одной экологической нише может существовать только один вид. На рис приведены результаты экспериментов Гаузе для двух видов Parametium, занимающих одну экологическую нишу (рис. 9.3 а, б) и видами, занимающими разные экологические ниши (рис в). Рис а - Кривые роста популяций двух видов Parametium в одновидовых культурах. Черные кружки – P Aurelia, белые кружки – P. Caudatum б - Кривые роста P Aurelia и P. Caudatum в смешанной культуре. По Gause, 1934

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

Система ХИЩНИК+ЖЕРТВА Для взаимоотношений типа хищник жертва или паразит хозяин система уравнений (9.1) принимает вид: (9.7) Здесь, в отличие от (9.2) знаки b 12 и b 21 - разные. Как и в случае конкуренции, начало координат (9.8) является особой точкой типа неустойчивый узел.

Три других возможных стационарных состояния: (9.9) (9.10) (9.11) Таким образом, возможно выживание только жертвы (9.10), только хищника (9.9) (если у него имеются и другие источники питания) и сосуществование обоих видов (9.11). Последний вариант уже был рассмотрен нами в лекции 5. Возможные типы фазовых портретов для системы хищник-жертва представлены на рис. 9.4.

Рис Расположение главных изоклин на фазовом портрете вольтерровской системы хищник-жертва (9.7) при различных соотношениях параметров. Стрелками указано направление фазовых траекторий. Пояснения в тексте.

Изоклины горизонтальных касательных представляют собой прямые x 2 = – b 21 х 1 /c 2 + a 1 /c 2, х 2 = 0, а изоклины вертикальных касательных – прямые x 2 = ­– c 1 х 1 /b 12 + a 2 /b 12, х 1 = 0. Стационарные точки лежат на пересечении изоклин вертикальных и горизонтальных касательных.

Из рис. 9.4 видно следующее. Система хищник – жертва (9.7) может иметь устойчивое положение равновесия, в котoром популяция жертв полностью вымерла () и остались только хищники (точка 2 на рис. 9.4 а). Очевидно, такая ситуация может реализоваться лишь в случае, если кроме рассматриваемого вида жертв х 1 хищник х 2 – имеет дополнительные источники питания. Этот факт в модели отражается положительным членом в правой части уравнения для х 2. Особые точки (1) и (3) (рис. 9.4 а) являются неустойчивыми. Вторая возможность – устойчивое стационарное состояние, в котором популяция хищников полностью вымерла и остались одни жертвы – устойчивая точка (3) (рис ). Здесь особая точка (1) – также неустойчивый узел. Наконец, третья возможность – устойчивое сосуществование популяций хищника и жертвы (рис. 9.4 в), стационарные численности которых выражаются формулами (9.11).

Как и в случае одной популяции (см. Лекция 3), для модели (9.7) можно разработать стохастическую модель, но для нее нельзя получить решение в явном виде. Поэтому мы ограничимся общими рассуждениями. Допустим, например, что точка равновесия находится на некотором расстоянии от каждой из осей. Тогда для фазовых траекторий, на которых значения x 1, x 2 остаются достаточно большими, вполне удовлетворительной будет детерминистическая модель. Но если в некоторой точке фазовой траектории какая– либо переменная не очень велика, то существенное значение могут приобрести случайные флюктуации. Они приводят к тому, что изображающая точка переместится на одну из осей, что означает вымирание соответствующего вида.

Таким образом, стохастическая модель оказывается неустойчивой, так как стохастический дрейф рано или поздно приводит к вымиранию одного из видов. В такого рода модели хищник в конечном счете вымирает, это может произойти либо случайно, либо вследствие того, что сначала элиминируется популяция его жертвы. Стохастическая модель системы хищник – жертва хорошо объясняет эксперименты Гаузе (Гаузе, 1934), в которых инфузория Paramettum candatum служила жертвой для другой инфузории Didinium nasatum – хищника. Ожидавшиеся согласно детерминистическим уравнениям (9.7) равновесные численности в этих экспериментах составляли примерно всего по пять особей каждого вида, так что нет ничего удивительного в том, что в каждом повторном эксперименте довольно быстро вымирали либо хищники, либо жертвы (а за ними и хищники) Результаты экспериментов представлены на рис. 9.5.

Рис Рост Parametium caudatum и хищной инфузории Dadinium nasutum. Из: Gause G.F. The struggle for existence. Baltimore, 1934

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

Обобщенные модели взаимодействия двух видов Было предложено большое число моделей, описывающих взаимодействие видов, правые части уравнений которых представляли собой функции численностей взаимодействующих популяций. Решался вопрос о выработке общих критериев, позволяющих установить, какого вида функции могут описать особенности поведения временного численности популяции, в том числе устойчивые колебания. Наиболее известные из этих моделей принадлежат Колмогорову (1935, переработанная статья 1972) и Розенцвейгу (1963). А.Н.Колмогоров рассмотрел обобщенную модель взаимодействия биологических видов типа хищни жертва или паразит хозяин. Модель представляет собой систему двух уравнений общего вида (9.12)

В модель заложены следующие предположения: 1) Хищники не взаимодействуют друг с другом, т.е. коэффициент размножения хищников k 2 и число жертв L, истребляемых в единицу времени одним хищником, не зависит от y. 2) Прирост числа жертв при наличии хищников равен приросту в отсутствие хищников минус число жертв, истребляемых хищниками. Функции k 1 (x), k 2 (x), L(x), непрерывны и определены на положительной полуоси x, y 0. 3) dk 1 /dx0, k 2 (0) < 0 < k 2 ( ). С ростом численности жертв коэффициент размножения хищников монотонно убывает с возрастанием численности жертв, переходя от отрицательных значений, (когда нечего есть) к положительным. 5) Число жертв, истребляемых одним хищником в единицу времени L(x)>0 при N>0; L(0)=0. Возможные типы фазовых портретов системы (9.12) представлены на рис. 9.6:

Рис Фазовые портреты системы Колмогорова (9.12), описывающей взаимодействие двух видов при разных соотношениях параметров. Пояснения в тексте.

Стационарные решения (их два или три) имеют следующие координаты: (1). x=0; y=0. Начало координат при любых значениях параметров представляет собой седло (рис. 9.6 а-г). (2). x=A, y=0. (9.13) A определяется из уравнения: k 1 (A)=0. Стационарное решение (9.13) седло, если BA. Это устойчивый узел. Последний случай, который соответствует гибели хищника и выживанию жертвы, представлен на рис. 9.6 в. (3). x=B, y=C.(9.14)

Величина С определяется из уравнений: Точка (9.14) – фокус (рис.9.6 а) или узел (рис.9.6 г), устойчивость которых зависит от знака величины 2 = – k 1 (B) – k 1 (B)B+L(B)C. Если >0, точка устойчива, если

Модель взаимодействия двух видов насекомых (MacArthur, 1971) Модель, которую мы рассмотрим ниже, использовалась для решения практической задачи борьбы с вредными насекомыми с помощью стерилизации самцов одного из видов. Исходя из биологических особенностей взаимодействия видов, была написана следующая модель (9.16) Здесь x,y - биомассы двух видов насекомых. Трофические взаимодействия описываемых в этой модели видов – весьма сложные. Этим определяется вид многочленов в правых частях уравнений. Рассмотрим правую часть первого уравнения. Насекомые вида х поедают личинок вида у (член + k 3 y), но взрослые особи вида у поедают личинок вида х при условии высокой численности видов х или у или обоих видов (члены – k 4 xy, – у 2 ). При малых х смертность вида х выше, чем его естественный прирост (1–k 1 +k 2 x–x 2 < 0 при малых х). Во втором уравнении член k 5 отражает естественный прирост вида y; –k 6 y – самоограничение этого вида, –k 7 x – поедание личинок вида у насекомыми вида х, k 8 xy – прирост биомассы вида у за счет поедания взрослыми насекомыми вида у личинок вида х. На рис. 9.7 представлен предельный цикл, являющийся траекторией устойчивого периодического решения системы (9.16).

Рис 9.7. Фазовый портрет модели (9.16) (MacArthur, 1971) Жирная замкнутая кривая – предельный цикл, штрих-пунктир – изоклины вертикальных и горизонтальных касательных

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

Управление может сводиться к кратковременному, скачкообразному изменению величин численности х и у. Такой способ отвечает методам борьбы типа однократного уничтожения одной или обеих популяций химическими средствами. Из сформулированного выше утверждения видно, что для совместимых популяций этот метод борьбы будет малоэффективным, поскольку с течением времени система опять выйдет на стационарный режим. Другой способ изменение вида функций взаимодействия между видами, например, при изменении значений параметров системы. Именно такому, параметрическому, способу отвечают биологические методы борьбы. Так при внедрении стерилизованных самцов уменьшается коэффициент естественного прироста популяции. Если при этом мы получим другой тип фазового портрета, такой, где имеется лишь устойчивое стационарное состояние с нулевой численностью вредителя, управление приведет к желаемому результату – уничтожению популяции вредного вида. Интересно отметить, что иногда воздействие целесообразно применить не к самому вредителю, а к его партнеру. Какой из способов более эффективен, в общем случае сказать нельзя. Это зависит от имеющихся в распоряжении средств управления и от явного вида функций, описывающих взаимодействие популяций.

Модель А.Д.Базыкина Рассмотрим одну из изученных в этой книге моделей типа хищник- жертва. (9.17) Система (9.17) является обобщением простейшей модели хищник- жертва Вольтерра (5.17) с учетом эффекта насыщения хищников. В модели (5.17) предполагается, что интенсивность выедания жертв линейно растет с ростом плотности жертв, что при больших плотностях жертв не соответствует реальности. Для описания зависимости рациона хищника от плотности жертв могут быть выбраны разные функции. Наиболее существенно, чтобы выбранная функция с ростом x асимптотически стремилась к постоянному значению. В модели (9.6) использовалась логистическая зависимость. В модели Базыкина в роли такой функции выбрана гипербола x/(1+px). Вспомним, что такой вид имеет формула Моно, описывающая зависимость скорость роста микроорганизмов от концентрации субстрата. Здесь в роли субстрата выступает жертва, а в роли микроорганизмов – хищник.

Система (9.17) зависит от семи параметров. Число параметров можно уменьшить с помощью замены переменных: x (A/D)x; y (A/D)/y; t (1/A)t; =c/A; =PD/A; =E/D; =M/B В новых безразмерных переменных система имеет вид (9.18) и зависит от четырех параметров.

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

Рис Параметрический портрет системы (9.18) при фиксированных и малых

Поведение системы при различных соотношениях параметров может быть существенно различным (рис.9.9). В системе возможны: 1) одно устойчивое равновесие (области 1 и 5); 2) один устойчивый предельный цикл (области 3 и 8); 3) два устойчивых равновесия (область 2) 4) устойчивый предельный цикл и неустойчивое равновесие внутри него (области 6, 7, 9, 10) 5) устойчивый предельный цикл и устойчивое равновесие вне его (область 4). В параметрических областях 7, 9, 10 область притяжения равновесия ограничивается неустойчивым предельным циклом, лежащим внутри устойчивого. Наиболее интересно устроен фазовый портрет, соответствующий области 6 на параметрическом портрете. Детально он изображен на рис Область притяжения равновесия В 2 (заштрихована) представляет собой улитку, скручивающуюся с неустойчивого фокуса В 1. Если известно, что в начальный момент времени система находилась в окрестности В 1, то судить о том, придет ли соответствующая траектория в равновесие В 2 или на устойчивый предельный цикл, окружающий три точки равновесия С (седло), В 1 и В 2 можно лишь на основе вероятностных соображений.

Рис Набор фазовых портретов системы (9.18), возможных в конечной части первого квадранта и соответствующих областям 1 10 параметрического портрета рис. 9.8 (Базыкин, 1985) Рис Фазовый портрет системы 9.18 для параметрической области 6. Область притяжения В 2 заштрихована

На параметрическом портрете (9.7) имеются 22 различные бифуркационные границы, которые образуют 7 различных типов бифуркаций. Их изучение позволяет выявить возможные типы поведения системы при изменении ее параметров. Например, при переходе из области 1 в область 3 происходит рождение малого предельного цикла, или мягкое рождение автоколебаний вокруг единственного равновесия В. Аналогичное мягкое рождение автоколебаний, но вокруг одного из равновесий, а именно B 1, происходит при пересечении границы областей 2 и 4. При переходе из области 4 в область 5 устойчивый предельный цикл вокруг точки B 1 лопается на петле сепаратрис и единственной притягивающей точкой остается равновесие B 2 и т.д.

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

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

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

Приложение Математическая модель изменения концентрации токсина в крови

Пусть [мл] - весь объем крови пациента во время t. Тогда общая сумма ядовитого вещества в крови во время t будет равна произведению * [мг/мл], где обозначает концентрацию ядовитого вещества во время t. Предположим, что некоторая очистка почек пациента [мл/минута] происходит так, что существует постоянный поток [мг/минута] ядовитого вещества из крови через почки в мочу, которая имеет форму: (2.1) Далее предполагаем, что, в пределах данного временного интервала [0, T] пациент привязан к диализеру (название искусственной почки) в течение временного интервала [0, ], где

где [мг/мл] обозначает концентрацию ядовитого вещества в жидкости во время t. Временным изменением концентрации ядовитого вещества управляет дифференциальное уравнение:, (2.3) где параметр F(t) дан в формуле (2.2), и G [мг/минута] является средней скоростью ядовитого вещества. Если предполагаем, что параметр известен для t Є [0,T], то двумя неизвестными функциями в (2.3) являются и. Предположим впредь, что ультрафильтрацией можно пренебречь, то есть. Это подразумевает, что объем крови больше не зависит от времени t и является постоянным [мл]. Поэтому (2.3) принимает форму: (2.4) где обозначает производную относительно параметра t.

Пусть: Решение (2.4) для > 0 тогда будет равна: (2.5) для и (2.6) для Если, то решением (2.4) будет: (2.7) для и (2.8) для

Заданные величины и T временное развитие параметра будет приближаться к устойчивому положению, которое выражено тем, что - это T-периодическая функция времени t. Для этого необходимо и достаточно, чтобы равенство (2.9) сохранялось. Далее принимаем =0, чтобы облегчить математические рассмотрения. Тогда (2.7) и (2.8): и равенство (2.9) оказывается эквивалентным уравнению: (2.10)

Пусть = 150 [мг/минута] и = [мл]. Тогда получаем: График изменения концентрации токсина в крови строится двумя способами. В первом случае график строится по формуле: Во втором случае по формуле:

Диализ выполняется через день, поэтому значение параметра T = 48часов(h) = 2880 минут. В работе для диаметра отверстия дается значение 150 [мг/мин] и для объема очищаемой крови [мл]. Также фиксированные значения представлены в таблице 1 и в таблице 2: Таблица 1 - Изменение концентрации токсина в крови Таблица 2 - Изменение концентрации токсина в крови ПРИМЕР 4h5h6h7h8h h5h6h7h8h

Математическая модель определения эффекта диализа

Эффект диализа в работе [4] также определяется двумя формулами. В первом случае для графика дана формула: (2.11) и определена таблица 3: Таблица 3 – Эффект диализа Во втором случае для построения графика дана формула: (2.12) ПРИМЕР 4h4h5h5h6h6h7h7h8h8h

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

В дифференциальных уравнениях (4.36), (4.37) работы [4] количества, и G доступны для измерения или экспериментального определения. Времена и T могут быть также измерены. Это, однако, не возможно для начальной величины. Можно доказать, что для каждого выбора параметров и в системах (4.36), (4.37) работы [4] существует одна пара решений с условиями: для всех и (2.13) Будем использовать этот факт, чтобы определить неизвестный параметр для мочевины внутри пределов. С этой целью предполагаем, что. Сначала рассмотрим системы (4.36), (4.37) работы [4] во временном интервале. Легко показать, что функция различия равна: и удовлетворяет дифференциальному уравнению: (2.14) гдеи

Для получения величины из решения (2.14) на промежутке считаем, что Пусть [мл] и [мл]. Выбираем T=2880 и [минут]. Считается, что параметр для мочевины, располагается между 700 и 900 [мл/минута]. Если возьмем, то получаем, и. Это подразумевает и, в связи с условием (2.13) из T-периодичности, отношение (2.15)

То же самое считается верным для маленьких величин. Например, для =300 [мл/минута], получаем и который также приводит к (2.15). Чтобы определить параметр, выбираем величину для : =1.5 [мг/мл]. Тогда каждый экспериментально выбирает величину для параметра и определяет скорость токсина G таким образом, что данная величина и соответствующая ей величина, полученные из (2.15), оказываются начальными величинами решения T-периода формул (4.36), (4.37). Для определения диаметра отверстий мембраны для мочевины используется функция: для. где [мл/мин]. ПРИМЕР :

Математическая модель лечения лекарствами

Есть различные пути, которыми лекарства могут быть введены в организм. Лекарства могут быть введены в виде инъекций, вводимых внутримышечно или внутривенно. Также лекарства бывают в виде таблеток, которые, проходя через желудок, поглощаются через стены кишечников. В первом случае лекарство прибывает непосредственно в ту часть тела, где это необходимо. Этот процесс может быть смоделирован моделью c одним отделением. В модели с одним отделением предполагается, что лекарство дано в дозировках в разное время, для с. После администрации дозировки предполагается, что количество лекарства разлагается по экспоненте в пределах временного интервала для (где ) согласно закону: для (2.29) Это желательно, чтобы для интервала всего времени [0,T] определенный уровень лекарства был поддержан. Это конечно не возможно. Графически получается следующая ситуация, показанная на рисунке 1: ПРИМЕР ПРИМЕР …ПРИМЕР Рисунок 1 – Временное изменение количества лекарства …

Математическая модель лечения диабета

Диабет вызван фактом, что гормональный инсулин произведен в недостаточном количестве в поджелудочной железе, и поэтому нужно позаботиться о процессе, в котором глюкоза транспортируется от крови в соматические клетки или сохраняется в печени в форме гликогена. Это имеет такой эффект, что уровень глюкозы в крови становится слишком высоким и это может вызвать смертельную кому. Обозначим концентрацию глюкозы в крови во время t через G(t) и гормональную концентрацию через H(t). Тогда временное развитие G(t) и H(t) может быть описано следующими дифференциальными уравнениями: (2.30) где p=p (t) обозначает увеличение концентрации глюкозы потреблением сахара.

Предполагаем, что без потребления сахара, то есть,, две постоянные состояния величин, и станут равными: и (2.31) Далее предполагаем, что G и H не сильно отличаются от и и определяют функции различия: (2.32) Линеаризуя (2.30) получаем дифференциальные уравнения: (2.33) где ПРИМЕР