Численно-аналитическое исследование математических моделей популяционной динамики (
На правах рукописи
Апонина Елена Александровна ЧИСЛЕННО-АНАЛИТИЧЕСКОЕ ИССЛЕДОВАНИЕ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ПОПУЛЯЦИОННОЙ ДИНАМИКИ (03.00.02 – Биофизика)
Автореферат диссертации на соискание учёной степени кандидата физико-математических наук
Пущино - 2008
Работа выполнена в Институте математических проблем биологии РАН
Научный консультант: доктор физико-математических наук Базыкин Александр Дмитриевич
Официальные оппоненты: доктор физико-математических наук, профессор Ризниченко Галина Юрьевна доктор биологических наук, профессор Фрисман Ефим Яковлевич Ведущая организация – Вычислительный центр им. А.А. Дородницына РАН. г. Москва
Защита состоится _ октября 2008 г. в на заседании диссертационного совета Д 501.001.96 при Московском государственном Университете имени М.В. Ломоносова по адресу: 119991, г. Москва, Ленинские горы, МГУ, биологический факультет, кафедра биофизики, аудитория “новая”.
С диссертацией можно ознакомиться в библиотеке биологического факультета МГУ имени М.В. Ломоносова.
Автореферат разослан сентября 2008 г.
Учёный секретарь диссертационного совета доктор биологических наук, профессор Т.Е. Кренделева Построение и исследование математических моделей является одним из наиболее распространённых методов научного познания. Сегодня математическое моделирование становится действенным инструментом исследования почти в каждой области биологической науки. С помощью этого метода получены интересные и важные результаты в биофизике, биохимии, микробиологии, в популяционной генетике, в экологии, нашедшие отражение в многочисленных статьях и монографиях (Романовский Ю.М., Степанова Н.М., Чернавский Д.С., 1971, 1975, 1984;
Свирежев Ю.М., Елизаров Е.Я., 1972;
Гимельфарб А.А., Гинзбург Л.Р., Полуэктов Р.А., Пых Ю.А., Ратнер В.А., 1974;
Алексеев В.В., 1976;
Рубин А.Б., 1976, 1984, 2004;
Фрисман Е.Я., Шапиро А.П., 1977;
Иваницкий Г.Р., Кринский В.И., Сельков Е.Е, 1978;
Свирежев Ю.М., Логофет Д.О., 1978;
Скалецкая Е.И., Фрисман Е.Я., Шапиро А.П., 1979;
Свирежев Ю.М., Пасеков В.П., 1982;
Пых Ю.А., 1983;
Рубин А.Б., Шинкарёв В.П., 1984;
Базыкин А.Д., 1985, 2003;
Свирежев Ю.М., 1987;
Заславский Б.Г., Полуэктов Р.А., 1988;
Молчанов А.М., 1992;
Алексеев В.В., Крышев И.И., Сазыкина Е.Г., 1992;
Ризниченко Г.Ю., Рубин А.Б., 1993, 2004;
Ризниченко Г.Ю., 2002, 2003).
В последние годы в связи с возрастающей мощностью современных ЭВМ и суперкомпьютеров открываются новые возможности для построения имитационных моделей биологических систем. Однако довольно часто возникает ситуация, когда результаты численных экспериментов, проведённых на таких моделях, с трудом поддаются объяснению и зачастую совсем непонятны.
Опыт моделирования сложных биосистем показывает, что новое знание о моделируемом объекте, углублённое понимание его свойств возникает при построении и исследовании иерархии “вложенных” друг в друга моделей увеличивающейся сложности (Полетаев И.А., 1971;
Гильманов Т.Г., 1978;
Галицкий В.В., Тюрюканов А.Н., 1981;
Базыкин А.Д., 1990). В связи с этим, наряду с разработкой имитационных моделей, становится актуальным выделение и исследование базовых моделей, т.е.
аналитических моделей, относящихся к нижним этажам создаваемой сегодня иерархии математических моделей биологических систем.
Следует также отметить, что исследование базовых моделей представляет интерес не только с теоретической точки зрения.
Аналитические модели с небольшим числом переменных и параметров успешно применяются при решении прикладных задач, например, при математическом описании и оптимизации процессов лабораторного и промышленного культивирования микроорганизмов (Monod J., 1950;
Novick A., Szilard L., 1950;
Степанова Н.В., Романовский Ю.М., Иерусалимский Н.Д., 1965;
Романовский Ю.М., Степанова Н.В., Чернавский Д.С., 1971, 1975, 1984;
Свирежев Ю.М., Елизаров Е.Е., 1972;
Печуркин Н.С., Терсков И.А., 1973, 1975;
Рубин А.Б., Пытьева Н.Ф., Ризниченко Г.Ю., 1977;
Печуркин Н.С., 1978;
Перт С.Д., 1978;
Полуэктов Р.А., Пых Ю.А., Швытов И.А., 1980;
Станишкис Ю., 1984;
Гуревич Ю.Л., 1984;
Бирюков В.В., Кантере В.М., 1985;
Варфоломеев С.Д., Калюжный С.В., 1990;
Печуркин Н.С., Брильков А.В., Марченкова Т.В., 1990;
Вавилин В.А., Васильев В.Б., Рытов С.В., 1993;
Ризниченко Г.Ю., Рубин А.Б., 1993, 2004;
Ризниченко Г.Ю., 2002;
Минкевич И.Г., 2005).
Цель работы. Построить и исследовать математические модели динамики численности популяций, взаимодействующих по принципу «продуцент – консумент – хищник» и «две конкурирующие жертвы – один хищник», а также математическую модель динамики численности плазмид в бактериальных популяциях. Исследование моделей довести до построения параметрических и фазовых портретов. Дать содержательную биологическую интерпретацию параметрическим и фазовым портретам.
Методы исследования. При исследовании математических моделей использовались аналитические методы качественной теории дифференциальных уравнений, теории устойчивости и бифуркаций.
Использовались также вычислительные алгоритмы и компьютерные программы для исследования бифуркаций сепаратрис (Апонин Ю.М., Апонина Е.А., 1976;
Кузнецов Ю.А., 1983;
Govaerts W., Kuznetsov Yu.A.
and Sijnave B., 2000), отыскания предельных циклов и слежения за предельным циклом при изменении параметра (Хибник А.И., 1979;
Khibnik A.I., Kuznetsov Yu.A., Levitin V.V., Nikolaev E.V., 1993).
Основными задачами
работы являются:
а) анализ роли нижнего порога численности продуцента в формировании типов динамического поведения трёхзвенных трофических цепей «продуцент – консумент – хищник»;
б) исследование вольтерровской модели системы трёх популяций, взаимодействующих по принципу «хищник – две конкурирующие жертвы»;
в) математическое моделирование популяционной динамики гибридных плазмид, исследование области устойчивости генно инженерных плазмидных штаммов при длительном непрерывном культивировании.
Научная новизна.
1. Построена и впервые исследована базовая математическая модель экосистемы «продуцент – консумент – хищник» при наличии у продуцента нижней критической плотности популяции. Показано, что существование нижней критической плотности популяции продуцента порождает разнообразие пороговых эффектов при кратковременных воздействиях на экосистему, вызывающих резкие изменения численностей составляющих её видов. На математической модели продемонстрирована роль хищника в поддержании жизнеспособности экосистемы. Установлена возможность устойчивого сосуществования продуцента, консумента и хищника при невозможности сосуществования продуцента и консумента в отсутствие хищника.
2. Впервые установлено существование бесконечного множества предельных циклов и сложного динамического поведения (незатухающих хаотических колебаний) в вольтерровской модели системы трёх взаимодействующих популяций.
3. Впервые описан механизм возникновения сложного динамического поведения в системе трёх популяций, взаимодействующих по принципу «хищник – две конкурирующие жертвы». Показано, что возникновение хаотических колебаний численностей популяций при увеличении промысловой нагрузки на экосистему является критерием приближения к опасной параметрической границе, за которой сосуществование не только трёх, но и двух популяций становится невозможным.
4. Построена и впервые исследована математическая модель непрерывного культивирования микроорганизмов, содержащих гибридные плазмиды, с учётом структурной и репликационной (сегрегационной) нестабильности плазмид. Модель позволяет рассчитывать область значений управляемых параметров процесса культивирования, при которых популяция плазмидосодержащих клеток устойчиво поддерживается на протоке. Впервые проанализирована зависимость этой области от значений кинетических параметров плазмидного штамма и его дериватов, возникающих вследствие структурной и репликационной нестабильности плазмид. Впервые выделена область стационарных значений управляемых параметров, при которых исходный плазмидный штамм поддерживается на протоке в режиме незатухающих колебаний его численности.
Эти основные результаты и выносятся на защиту.
Практическая значимость.
1. Рассмотренные в диссертации математические модели представляют собой нелинейные системы трёх – четырёх обыкновенных дифференциальных уравнений. Исследование каждой модели доведено до построения параметрического и соответствующих фазовых портретов.
Поэтому все рассмотренные модели применимы как базовые модели для построения и исследования последующих уровней иерархии математических моделей экологии и биотехнологии.
2. Математическая модель популяционной динамики плазмид использовалась для объяснения колебательной динамики численности популяции клеток в непрерывной культуре генно-инженерного штамма Bacillus stearothermophilus, несущих рекомбинантную плазмиду pZAM (Koizumi J.-I., Aiba S., 1988). В последние годы эта модель используется при построении и исследовании более сложных моделей с учётом копийности и конъюгационного переноса плазмид (Апонин Ю.М., Апонина Е.А., 2008).
3. Математическая модель популяционной динамики плазмид применима при решении прикладных задач направленного генно инженерного конструирования гибридных плазмид (методами генетической инженерии) и прогнозирования эффективности промышленного использования штаммов, их содержащих.
4. Результаты исследования базовых моделей указывают критерии приближения к опасным границам (в параметрическом и фазовом пространстве) функционирования экологических и биотехнологических систем.
Апробация работы и публикации. Материалы диссертации докладывались на Всесоюзной конференции «Теория и практика программирования на ЭВМ серии Мир» (Душанбе, 1974), Всесоюзной конференции по асимптотическим методам в теории сингулярно возмущённых уравнений (Алма-Ата, 1979), Годовой конференции НИВЦ (Пущино, 1981), I Всесоюзном биофизическом съезде (Москва, 1982), Шестой международной конференции «Математика. Компьютер.
Образование» (Пущино, 1999), Четырнадцатой международной конференции «Математика. Компьютер. Образование» (Пущино, 2007).
По теме диссертации опубликовано 29 печатных работ, из них статьи в российских журналах, 2 статьи в международных журналах.
Структура и объём диссертации. Диссертация состоит из введения, четырёх глав, основных результатов и выводов, списка цитируемой литературы, списка публикаций по теме диссертации. Работа изложена на 132 страницах машинописного текста. Список цитируемой литературы содержит 199 наименований.
Глава 1. Иерархия моделей математической биологии и численно аналитические методы их исследования Эта глава носит обзорный характер. В ней рассматриваются некоторые общие понятия и концепции математической биологии:
многомоделие, интермодельные отношения, иерархия моделей, минимальная модель, базовая модель, максимальная модель и др.
Характерной особенностью математической биологии является разнообразие математических моделей биосистем. Очень часто для описания функционирования заданного биологического объекта используется не одна модель, а целый ряд математических моделей различной степени сложности. В первой главе обсуждается проблема иерархии интермодельных отношений. Выделены основные механизмы генерации простых моделей из сложных: замена переменных и асимптотическая декомпозиция, введение агрегированных переменных, принцип сведения и др. Рассматривается и обратный процесс построения сложной модели путём поэтапного усложнения некоторого исходного набора простых моделей, а также развёртывание иерархии моделей различной степени сложности.
Строительство иерархической сети математических моделей может начинаться, так сказать, высоко в воздухе с некоторой исходной модели и разветвляться как по направлению вверх, так и вниз. Вверх разветвления могут расходиться в разные стороны к различным более сложным моделям. Необходимость продвижения вниз возникает в связи с потребностью более глубокого исследования исходной модели, например, путём асимптотической декомпозиции её на более простые модели меньшей размерности.
В зависимости от цели моделирования можно построить и другую сеть моделей того же самого объекта реальности, выбрав в качестве исходной другую модель, другой исходный узел разветвлений. Однако каждая такая иерархическая сеть представляет собой лишь фрагмент создаваемой сегодня иерархии моделей математической биологии.
Иерархия моделей – это не застывшее, а постоянно развивающееся образование. На каждом этапе своего развития её можно рассматривать как частично упорядоченное множество, в котором существуют непустые подмножества как минимальных, так и максимальных элементов. Эти элементы называются соответственно минимальными и максимальными моделями.
В первой главе даётся также краткий обзор численно-аналитических методов качественной теории дифференциальных уравнений, теории устойчивости и бифуркаций. Подчёркивается важная роль аналитических методов исследования математических моделей. Применение аналитических методов начинается с обезразмеривания уравнений математической модели, выявления малых безразмерных параметров. Во многих случаях удаётся получить явные аналитические формулы для координат состояний равновесия и вывести простые аналитические условия локальных бифуркаций. Критерии Бендиксона и Дюлака позволяют установить отсутствие предельных циклов на фазовой плоскости. Аналитическое исследование состояний равновесия в конечной части фазовой плоскости и на бесконечности помогает иногда положительно решить вопрос о существовании предельного цикла и установить расположение сепаратрис. При аналитическом исследовании двумерных математических моделей оказываются полезными понятия циклов и кривых без контакта, свойства поворота векторного поля (Андронов А.А., Леонтович Е.А., Гордон И.И., Майер А.Г., 1966, 1967;
Баутин Н.Н., Леонтович Е.А., 1990).
Аналитические методы позволяют получить интересные и важные результаты и при исследовании многомерных математических моделей.
Например, в работе [3] рассматривалась математическая модель функционирования антибиотика нигерицина в биологических мембранах.
Соответствующая кинетическая модель включала в себя 32 элементарные реакции, а список реагирующих веществ содержал 18 компонентов ( переменных математической модели). На основании аналитического исследования установлены существование и единственность состояния равновесия модели и предложен эффективный алгоритм для вычисления координат этого равновесия. С учётом результатов аналитического исследования разработана компьютерная программа для проведения массовых расчётов на ЭВМ интенсивности осуществляемого нигерицином K+/ H+ – обмена через бислойную липидную мембрану в зависимости от концентрации нигерицина и некоторых управляемых параметров.
Проведённые расчёты показали, что модель описывает все наблюдаемые в эксперименте типы нелинейных и колеблющихся зависимостей трансмембранных потоков ионов калия и водорода от концентрации нигерицина в мембране и концентраций некоторых буферных соединений в окружающих её растворах [3].
Особую роль в математической биологии играет понятие базовой математической модели (Романовский Ю.М., Степанова Н.В., Чернавский Д.С., 1984;
Чернавский Д.С., Старков Н.И., Щербаков А.В., 2002;
Ризниченко Г.Ю., 2003;
Рубин А.Б., 2004;
Чернавский Д.С., Чернавская Н.М., Малков С.Ю., Малков А.С., 2005;
Малинецкий Г.П., 2005). Так называются математические модели, допускающие детальное исследование аналитическими методами и, возможно, с использованием доказательных вычислений на ЭВМ. Например, математические модели, представленные нелинейными системами обыкновенных дифференциальных уравнений, исследование которых доведено до построения параметрического и фазовых портретов, являются базовыми. В этом смысле математические модели, рассмотренные в главах 2, 3 и настоящей диссертации, являются базовыми. Исследование каждой из этих моделей доведено до построения детального параметрического портрета с указанием соответствующего набора трёхмерных фазовых портретов.
Глава 2. Модель экосистемы трёх трофических уровней с учётом существования нижней критической плотности популяции продуцента В этой главе проводится исследование математической модели экосистемы, состоящей из трёх популяций, относящихся к различным трофическим уровням и взаимодействующих между собой по принципу «продуцент – консумент – хищник». Для популяции продуцента было принято допущение о существовании нижней критической её плотности.
Известно, что многие биологические популяции не могут возрастать и даже сохранять свою численность, если их размер опускается ниже некоторого критического значения. При падении плотности популяции ниже этой пороговой величины из-за неблагоприятных условий или в результате хищнического промысла восстановление популяции становится невозможным и она обречена на вымирание (Тимофеев-Ресовский Н.В., Яблоков А.В., Глотов Н.В., 1973;
Дрё Ф., 1976;
Леме Ж., 1976;
Фрисман Е.Я., Худолей Ю.И., 1976, 1977;
Базыкин А.Д., Березовская Ф.С., 1978, 1979;
Базыкин А.Д., 1985, 2003;
Ризниченко Г.Ю., 2002, 2003).
Изложение в главе 2 начинается с обсуждения вопроса о причинах существования нижней критической плотности. Рассматривается общее уравнение динамики отдельной популяции с перекрывающимися поколениями:
dx = r ( x ) x, где r ( x ) = f ( x ) g ( x ). (2.1) dt Здесь x – плотность популяции, r – удельная скорость изменения плотности, f и g – функции рождаемости и смертности.
Отмечается, что для ряда популяций, допускающих описание вида (2.1), выполняются следующие условия:
• графики функций f и g пересекаются в двух точках с абсциссами x = L и x = K, причём 0 L K;
• для любого x 0 неравенство f(x) g(x) истинно только при x (L, K).
Предположим, что эти условия выполнены. Тогда график функции r = f – g устроен качественно как перевёрнутая парабола, пересекающая ось абсцисс в точках x = L и x = K. В таком случае при аппроксимации полиномом второй степени результатов лабораторных или полевых измерений функции r(x) получается следующая зависимость удельной скорости изменения плотности популяции r = x / x от самой этой плотности x:
x a r ( x ) = a ( x L) 1 = rm ( x xm ) 2, (2.2) K K где K+L a K L xm = rm =,. (2.3) 2 K Здесь xm – плотность популяции, при которой удельная скорость r(x) достигает своего максимального значения rm, rm = r(xm). В экологической литературе величина rm называется биотическим потенциалом популяции.
С учётом (2.2) уравнение (2.1) принимает следующий вид:
x x = a x ( x L) 1. (2.4) K В работах Базыкина А.Д. & Березовской Ф.С. (1978, 1979) уравнение (2.4) рассматривалось как эталонное (наиболее простое) уравнение динамики популяций с критическим порогом плотности. В настоящее время это уравнение и его модификации применяются в общей теории роста гомеостатических систем (Кряжимский Ф.В., 2007) и при математическом моделировании динамики роста лесных насаждений (Исаев А.С., Суховольский В.Г., Овчинникова Т.М., 2008).
При заданных значениях L и K параметр а связан с биотическим потенциалом rm линейной зависимостью, см. вторую формулу в (2.3).
Поэтому для популяций, подчиняющихся уравнению (2.4), величину а естественно называть константой скорости роста популяции.
Интерпретация параметров L и K очевидна: L – нижняя критическая плотность популяции, K – значение плотности в устойчивом состоянии равновесия, которое достигается при заданных условиях среды обитания.
Как и в логистическом уравнении, величину K естественно называть «ёмкостью местообитания» по отношению к данной популяции (Одум Ю., 1975).
Во второй главе диссертации уравнение (2.4) используется как базовая модель динамики популяции продуцента при построении более сложной математической модели сообщества «продуцент – консумент – хищник»:
x x = a x ( x L ) 1 g1 x y, K y = h1 x y g 2 y z c1 y, (2.5) z = h2 y z c2 z.
Здесь x, y и z – соответственно плотность популяций продуцента, консумента и хищника;
g1 и g2 – соответственно константы скорости выедания продуцента консументом и консумента хищником;
h1 и h2 – константы скорости роста популяции консумента и хищника;
с1 и с2 – эффективная смертность соответственно консумента и хищника, которая может иметь три аддитивные составляющие: естественную смертность, промысловую смертность (отлов, отстрел) и интенсивность эмиграции.
Параметры a, K, L имеют тот же смысл, что и в уравнении (2.4).
aK hK, z = 1, t = x = Ku, y = Заменой переменных g1 g2 aK система уравнений (2.5) приводится к следующей безразмерной системе:
u = u ( u l )(1 u ) u, = 1 ( u m1 ), (2.6) = 2 ( m2 ), где u,, – безразмерные переменные (плотности популяций соответственно продуцента, консумента и хищника в выбранном масштабе их измерения). Здесь точка над переменной обозначает её дифференцирование по безразмерному времени ;
l, 1, 2, m1, m2 – безразмерные параметры:
L h h c cg, 1 = 1, 2 = 2, m1 = 1, m2 = 2 1.
l= K a g1 h1 K a h2 K Параметр m1 (m2) представляет собой выраженную в безразмерном виде эффективную смертность в популяции консумента (хищника). Значения m и m2 можно варьировать за счёт изменения интенсивности промысла.
Поэтому величины m1 и m2 естественно рассматривать как управляемые параметры, характеризующие «сопротивление среды» (Одум Ю., 1975) по отношению к популяции соответственно консумента и хищника.
Основными результатами проведённого в главе 2 исследования системы (2.6) являются её параметрический и фазовые портреты. На параметрическом портрете, построенном в плоскости управляемых параметров (m1, m2), выделено 12 областей. Для каждой из этих областей предъявлен соответствующий трёхмерный фазовый портрет системы (2.6).
На основании анализа параметрического и фазовых портретов математической модели (2.6) сделаны следующие выводы.
1. Существует область значений управляемых параметров m1 и m2, при которых в системе сосуществуют вместе сразу три аттрактора. Два из них представляют собой устойчивые состояния равновесия: равновесие R, расположенное внутри первого октанта, и расположенное в начале координат равновесие О, соответствующее гибели всех составляющих экосистемы. Третий аттрактор лежит в плоскости (u, ) и в зависимости от параметров является либо состоянием равновесия P, либо предельным циклом CP, родившимся из равновесия P в результате бифуркации Андронова – Хопфа. Каждый из этих аттракторов имеет свою область притяжения в пространстве состояний.
2. Множественность аттракторов и их бифуркации порождают гистерезисные эффекты, которые могут наблюдаться при плавных изменениях управляемых параметров. Например, зависимость установившейся плотности популяции хищника от интенсивности его промысла (от параметра m2) оказывается неоднозначной, содержит точки разрыва и петлю гистерезиса. Такого типа гистерезисные зависимости от параметра m2 реализуются при фиксированных значениях эффективной смертности консумента m1, принадлежащих некоторому интервалу.
Верхняя граница этого интервала вычислена аналитически (m1 = (1 + l) / 2), а его нижняя граница (m1 = mсп) соответствует бифуркации слияния сепаратрис и легко вычисляется на ЭВМ с помощью разработанных компьютерных программ [19].
3. Существование нижней критической плотности популяции продуцента порождает разнообразие пороговых эффектов в динамических реакциях экосистемы в ответ на кратковременные внешние воздействия, вызывающие резкие изменения численности составляющих её видов. Эти эффекты объясняются чашеобразной формой области притяжения равновесия R и своеобразием строения сепаратрисных поверхностей, разделяющих области притяжения альтернативных аттракторов в пространстве состояний.
4. Математическая модель указывает на важную роль хищника в поддержании жизнеспособности экосистемы. Например, при низких значениях эффективной смертности консумента m1 (точнее, при m1 mcn) его существование в отсутствие хищника оказывается невозможным: в отсутствие хищника в системе «продуцент – консумент» обречены на вымирание оба составляющих её вида. Однако интродукция хищника в этих условиях спасает систему от гибели: хищник закрепляется в экосистеме и оказывается возможным сосуществование продуцента, консумента и хищника в устойчивом состоянии равновесия R.
Глава 3. Исследование математической модели трёхвидового сообщества “хищник – две жертвы” В этой главе рассматривается вольтерровская система трёх обыкновенных дифференциальных уравнений, описывающая динамику трёхвидового биологического сообщества «хищник – две конкурирующие жертвы». Первые исследования такого рода модельных систем были предприняты в связи с работой Paine (1966), в которой сделан вывод о том, что хищники могут играть главную роль в поддержании видового разнообразия нижестоящего трофического уровня. Изучая сообщество беспозвоночных обитателей литоральной (приливной) зоны тихоокеанского побережья США в штате Вашингтон, Paine (1966) в результате ряда проведённых экспериментов обнаружил, что удаление из сообщества хищника (морской звезды Pisaster ochraceus) приводит к существенному снижению числа видов жертв (моллюсков и ракообразных), которыми питался этот хищник.
Parrish & Saila (1970), Cramer & May (1972) попытались воспроизвести данные наблюдений Paine (1966) на вольтерровской модели сообщества «две жерты – один хищник». Исследуя эту модель, Cramer & May (1972) установили, что присутствие хищника может стабилизировать сообщество двух конкурирующих видов жертв даже в том случае, когда в отсутствие хищника сосуществование этих видов невозможно. Несколько позже, Fujii (1977), исследуя ту же модель, показал, что стабилизирующая роль хищника может проявляться не только в стационарном, но и в колебательном режиме: при определённых значениях параметров модели внутри первого октанта существует устойчивый предельный цикл. Vance (1978) численно обнаружил сложный колебательный режим поведения системы, который он назвал «квазициклическим».
Заметим, что во всех упомянутых выше работах исследование вольтерровской модели сообщества «две жертвы – один хищник» не было доведено до построения бифуркационной диаграммы соответствующей системы обыкновенных дифференциальных уравнений (т.е. построения её параметрического и фазовых портретов).
В работе [8] показано, что в некоторой области значений параметров в вольтерровской системе «две жертвы – один хищник» существуют устойчивые предельные циклы сложной формы, бесконечное множество седловых предельных циклов и более сложные траектории нерегулярных (хаотических) автоколебаний. В главе 3 диссертации излагаются, с некоторыми уточнениями, результаты работы [8]. Рассматривается вольтерровская модель системы из трёх популяций, взаимодействующих по принципу «две жертвы – один хищник»:
N1 = N1 ( N1 a1 N 2 b1 P ), N 2 = N 2 ( N 2 a2 N1 b2 P ), (3.1) P = P ( 1 + d1 N1 + d 2 N 2 P ).
Здесь N1 и N2 – соответственно численности популяций первого и второго вида жертв, P – численность популяции хищника. Параметры, характеризующие интенсивность внутривидовой конкуренции во всех трёх популяциях, и скорость вымирания хищника в отсутствие жертв посредством нормировки (т. е. выбора единиц измерения численностей популяций и времени) сделаны равными единице.
В отличие от модели, изучавшейся в работах Parrish & Saila (1970), Cramer & May (1972), Fujii (1977) и Vance (1978), в модели (3.1) учитывается и внутривидовая конкуренция в популяции хищника.
Значения параметров межпопуляционных взаимодействий предполагаются фиксированными:
a1 = 6, a2 = 1, b1 = 4, b2 = 10, d1 = 0.25, d 2 = 4. (3.2) Эти значения подобраны таким образом, чтобы при некоторых значениях параметров и в системе (3.1) существовал сепаратрисный контур, проходящий через три состояния равновесия, одно из которых является седлофокусом, см. рис. 3.1.
Рис. 3.1. Схематическое изображение сепаратрисного контура системы (3.1) – (3.2), проходящего через седлоузлы А5, А2 и седлофокус А4. Такой контур реализуется при значениях и на некоторой линии L в плоскости параметров (, ). При этом на оси N1 существует устойчивое состояние равновесия А1. Внутри первого октанта существует также неустойчивое состояние равновесия А6 седлового типа, которое не показано, чтобы не усложнять рисунок.
В плоскости параметров (, ) построен параметрический портрет системы (3.1) – (3.2). На параметрическом портрете выделено 19 областей, для каждой из которых приведён соответствующий трёхмерный фазовый портрет. Бифуркационное множество системы (3.1) – (3.2) включает в себя линию сепаратрисного контура L: при (, ) L в системе реализуется сепаратрисный контур, см. рис. 3.1.
Бифуркационные явления, происходящие в окрестности сепаратрисного контура могут быть весьма сложными. При исследовании этих бифуркаций использовались некоторые общие математические результаты, полученные в работах Шильникова Л.П. (1965, 1970) и Белякова Л. А. (1981, 1984).
Показано, что с уменьшением параметра при подходе к линии сепаратрисного контура L могут происходить бесконечные серии бифуркаций предельных циклов. Например, существует бесконечная серия бифуркаций рождения однообходных предельных циклов из уплотнения траекторий. Предельные циклы рождаются парами. В каждой паре один цикл седловой, а другой – устойчивый, причём с уменьшением параметра устойчивый цикл почти сразу после рождения становится седловым, отщепляя от себя устойчивый двухобходный цикл. В результате при подходе к линии сепаратрисного контура L происходит накопление (до бесконечного числа) седловых однообходных предельных циклов.
С уменьшением параметра при переходе через линию L бифуркационный процесс развивается в обратном порядке: некоторые из однообходных седловых циклов превращаются в устойчивые, захватывая соответствующие устойчивые двухобходные циклы, после чего они исчезают «аннигилируя» с другими однообходными седловыми циклами.
Процесс завершается слиянием и уничтожением последней пары однообходных циклов (устойчивого и седлового).
В численных экспериментах на модели (3.1) – (3.2) при фиксированном = 2,4 с уменьшением параметра наблюдалась также серия бифуркаций потери устойчивости предельным циклом с отщеплением от него устойчивого предельного цикла удвоенного периода.
Несколько предельных циклов различной обходности показано на рис. 3.2.
В некотором диапазоне значений параметра из результов численных экспериментов создаётся впечатление, что траектория системы всюду плотно заполняет сложно устроенную ограниченную «поверхность» в пространстве состояний (см. рис. 3.2 г). При этих значениях параметра в системе реализуются незатухающие хаотические колебания, что соответствует известным математическим результатам о существовании в окрестности сепаратрисного контура нетривиального притягивающего множества (Шильников Л.П., 1965).
Рис. 3.2.. Результаты численного эксперимента по анализу механизма возникновения сложного динамического поведения в системе (3.4);
x1 = N1, x2 = P, x3 = N2 ;
a) = 2,4, = 1,7537, x1(0) = 0, 40906, x2(0) = 0,10157, x3(0) = 0,32860;
б) = 2,4, = 1,7533, x1(0) = 0, 35304, x2(0) = 0,10367, x3(0) = 0,36352;
в) = 2,4, = 1,7532, x1(0) = 0,32070, x2(0) = 0,10451, x3(0) = 0,38733;
г) = 2,4, = 1.7531, x1(0) = 0,30155, x2(0) = 0,10480, x3(0) = 0,40307.
Перечислим кратко основные полученные в главе 3 результаты и выводы.
1. В плоскости параметров (, ) построен параметрический портрет вольтерровской модели сообщества трёх популяций, взаимодействующих по принципу «хищник – два конкурирующих вида жертвы». Параметр () характеризует разность между рождаемостью и эффективной смертностью в популяции первого (второго) вида жертвы. При этом эффективная смертность наряду с естественной может аддитивно включать в себя и промысловую смертность, а также параметры, характеризующие интенсивность эмиграции.
2. Выделено 19 областей параметрического портрета, для каждой из которых предъявлен соответствующий трёхмерный фазовый портрет математической модели.
3. На математической модели продемонстрирована роль хищника в поддержании видового разнообразия биологических сообществ. Согласно модели (3.1) – (3.2) сосуществование двух конкурирующих видов жертв в отсутствие хищника невозможно ни при каких значениях параметров и. Однако в широкой области значений этих параметров оказывается возможным устойчивое сосуществование трёх видов (двух видов жертвы и одного вида хищника).
4. Все три вида могут сосуществовать в стационарном или автоколебательном режиме в отсутствие каких-либо внешних дестабилизирующих возмущений. Вблизи границы области (в пространстве параметров и ) сосуществования трёх видов они могут сосуществовать также в режиме хаотических колебаний.
5. Одноразовые внешние воздействия, вызывающие резкие изменения численности видов, могут приводить к деградации сообщества (к исчезновению хищника и одного из видов жертвы). Например, существует область значений параметров и, при которых трёхвидовое сообщество оказывается неустойчивым к инвазии первого вида жертвы, вызывающей резкое увеличение её численности выше некоторого порога. При этом хищник и второй вид жертвы обречены на вымирание и в системе остаётся только первый вид жертвы. Аналогичное поведение возможно и в ответ на резкое увеличение численности хищника выше некоторого порога при однократной интродукции его в систему.
6. В плоскости (, ) вблизи границы области существования трёхвидового сообщества могут происходить бесконечные серии бифуркаций предельных циклов и реализовываться режимы хаотических колебаний. Эти бифуркации и хаотические режимы объясняются существованием, при определённых значениях параметров и, сепаратрисного контура, включающего в себя три особые точки, одна из которых является седлофокусом.
7. Установлено, что возникновение хаотических колебаний является критерием приближения в плоскости (, ) к опасной границе, за которой сосуществование не только трёх, но и двух видов становится невозможным ни в каком режиме. При переходе через эту границу популяция хищника и популяция одного из видов жертв обречены на вымирание.
Глава 4. Математическое моделирование процессов непрерывного культивирования микроорганизмов, содержащих нестабильные гибридные плазмиды В этой главе построена и исследована математическая модель непрерывной культуры микроорганизмов, полезная биосинтетическая активность которой (то есть синтез полезного продукта, сокращённо ПБА) поддерживается за счёт генно-инженерного штамма с рекомбинантными плазмидами. В модели учитывается два механизма снижения и потери ПБА культуры: репликационная и структурная нестабильность рекомбинантных плазмид.
Рассматривается хемостатная культура микроорганизмов, проявляющая некоторую ПБА, причём уровень ПБА определяется концентрацией клеток с активными, то есть направляющими сверхсинтез полезного продукта, плазмидами. Вследствие структурной и репликационной нестабильности в культуральной жидкости присутствуют также неактивные, то есть утратившие способность к суперпродукции, клетки двух типов: бесплазмидные клетки и клетки с мутантными (структурно изменёнными) плазмидами. Концентрации бесплазмидных, неактивных плазмидных и активных плазмидных клеток обозначим соответственно через N0, N1 и N2.
Предполагается, что превращение активных плазмидных клеток в неактивные клетки с изменёнными плазмидами аналогично реакции первого порядка N2 N k с константой скорости k. Предполагается также, что процессы деления клеток и потери плазмид можно охарактеризовать следующими реакциями (1 i ) µi Ni Ni + Ni µ Ni Ni + N i i где i = 1, 2, i – вероятность появления бесплазмидной клетки в одном акте деления клетки, несущей плазмиду. Здесь и далее через µi обозначается удельная скорость роста популяции клеток i-го типа (i = 0, 1, 2).
Пусть S – концентрация лимитирующего субстрата в ферментёре.
Тогда при сделанных выше предположениях изменение концентраций S, N0, N1 и N2 во времени описывается следующей системой дифференциальных уравнений:
S = D ( S0 S ) l0 µ0 N 0 l1 µ1 N1 l2 µ2 N N 0 = ( µ0 D ) N 0 + 1 µ1 N1 + 2 µ 2 N (4.1) N1 = ( µ1 D ) N1 1 µ1 N1 + k N N = µ D N µ N k N, (2 ) 2 22 2 где S0 – концентрация лимитирующего субстрата в свежей питательной среде, поступающей в ферментёр извне;
D – скорость протока;
li – стехиометрический коэффициент популяции клеток i-го типа (i = 0, 1, 2) ;
точкой сверху обозначено дифференцирование по времени t.
Относительно удельных скоростей роста принимается обычное для теории непрерывного культивирования предположение (Monod J., 1942, 1950):
µim S µi = µi ( S ) = i = 0, 1, 2,, (4.2) Ki + S где Ki – константа полунасыщения, µim – максимальная удельная скорость роста для популяции клеток i-го типа.
С учётом (4.2) математическая модель (4.1) представляет собой замкнутую систему четырёх дифференциальных уравнений, зависящую от четырнадцати параметров. В четвёртой главе диссертации проводится численно-аналитическое исследование математической модели (4.1) – (4.2) в наиболее важном частном случае, когда популяция бесплазмидных клеток по основным параметрам роста µm, K и l слабо отличается от популяции клеток с мутантными плазмидами. Точнее, предполагается, что µ0 m = µ1m, K 0 = K1, l0 = l1. (4.3) В силу (4.3) утратившие способность к ПБА бесплазмидные клетки и клетки с мутантными плазмидами можно рассматривать как единую (неактивную) популяцию, которая характеризуется параметрами роста µ1m, K1 и l1. Популяция клеток с интанктными плазмидами (активная популяция) характеризуется, вообще говоря, другими значениями этих параметров: µ2 m, K 2, l2.
Введём суммарную плотность микробной популяции N = N 0 + N1 + N 2 и положим = N 2 N. Тогда безразмерная переменная характеризует уровень обеспечиваемой плазмидами ПБА культуры. Введём также безразмерные переменные, x и безразмерное время t :
S l = t = µ1m t.
x = 1 N,, K1 K При условиях (4.3) из (4.1) – (4.2) следует, что переменные, x, удовлетворяют замкнутой системе дифференциальных уравнений:
= ( 0 ) (1 ) + x + 1+ x = (1 ) + x (4.4) 1+ + = (1 ), + 1+ + где = 2 ;
,,,, 0, – безразмерные параметры:
µ S K l k D = 2m, = 2, = 2, =, =, 0 = 0.
µ1m µ1m µ1m K1 l1 K Приведём основные результаты аналитического исследования системы (4.4). Установлено, что максимальное число состояний равновесия системы (4.4) равно четырём. При этом внутри первого октанта фазового пространства (x,, ) существует не более одного состояния равновесия, обозначим его через B. Остальные состояния равновесия лежат на границе первого октанта. При изменении параметров могут происходить бифуркации слияния состояний равновесия, причём возможно только четыре типа бифуркаций слияния коразмерности один.
Плоскость управляемых параметров (0, ) разбивается на области линиями бифуркаций слияния состояний равновесия. Эти области различаются по соответствующим им типам фазового портрета системы (4.4). Установлено, что при изменении параметров,,, могут реализовываться только три невырожденных случая этого разбиения (см.
рис. 4.1;
подчеркнём, что это разбиение не зависит от параметра ).
Рис. 4.1. Три возможных типа разбиения плоскости ( 0, ) на области линиями бифуркаций слияния положений равновесия системы (4.4). Отвечающие этим областям фазовые портреты системы (4.4) соответствующим образом пронумерованы на рис. 4.2.
Заштрихована область значений управляемых параметров 0 и, при которых в системе (4.4) существует положение равновесия B, то есть область устойчивости исходного плазмидного штамма.
Соответсвующие фазовые портреты системы (4.4) представлены на рис. 4.2.
На рис. 4.1 штриховкой выделена область значений управляемых параметров 0 и, при которых штамм с интактной плазмидой (активная популяция) не вымывается из ферментёра. Эта область называется областью устойчивости данного штамма. Она совпадает с областью значений параметров 0,, при которых внутри первого октанта пространства состояний (x,, ) сущесвует равновесие B.
Установлено, что при изменении параметров 0 и в области устойчивости исходного плазмидного штамма (показана штриховкой на рис. 4.1) может происходить бифуркация смены устойчивости равновесия В (бифуркация Андронова– Хопфа). Соответствующая бифуркационная линия в плоскости управляемых параметров (0, ) называется линией нейтральности. На линии нейтральности характеристическое уравнение в точке В имеет один отрицательный и пару чисто мнимых корней. При пересечении линии нейтральности, с вхождением внутрь ограниченной ею области, равновесие В теряет устойчивость и из него рождается устойчивый предельный цикл.
Рис. 4.2. Фазовые портреты системы (4.4) для семи областей параметрического портрета, соответствующим образом пронумерованных на рис. 4.1. На некоторых траекториях, примыкающих к точке B, направление движения не указано, так как характер устойчивости положения равновесия B не определяется однозначно по разбиению плоскости ( 0, ), представленному на рис. 4.1.
На основании аналитического исследования системы (4.4) разработана программа, посредством которой для любых заданных значений параметров,,,, на экран ЭВМ выводятся все линии, отвечающие бифуркациям состояний равновесия системы (4.4). Используя эту программу, можно проследить эволюцию параметрического портрета системы (4.4) в плоскости управляемых параметров (0, ) при любых изменениях параметров роста,, и параметров нестабильности,.
На рис. 4.3 представлены построенные численно параметрические портреты системы (4.4) при различных значениях параметра. Некоторые из областей пронумерованы в соответствии с нумерацией, принятой на рис. 4.1. При выбранных значениях структурных параметров область является настолько узкой, что оказывается практически неразличимой. На рис.4.3 область устойчивости исходного, активного (т.е. обеспечивающего ПБА культуры) плазмидного штамма состоит из двух частей, отделённых Рис. 4.3. Эволюция параметрического портрета системы (4.4) при изменении параметра ;
= 7, = 16, = 0,25, = 4,5. Штриховкой отмечена область устойчивости положения равновесия B, двойной штриховкой – ограниченная линией нейтральности область автоколебаний (т.е. область значений параметров 0 и, при которых в системе (4.4) существует устойчивый предельный цикл).
друг от друга линией нейтральности. Штриховкой отмечена та часть области устойчивости, в которой положение равновесия B устойчиво, то есть штамм с интактной плазмидой поддерживается на протоке в устойчивом стационарном состоянии. Во второй части области устойчивости, отмеченной двойной штриховкой, равновесие B неустойчиво, то есть исходный плазмидный штамм хотя и поддерживается на протоке, но стационарный режим невозможен и в системе устанавливаются незатухающие колебания (при 0 const, const). При этом ПБА культуры периодически то падает, то вновь нарастает.
Иллюстрирующий пример таких незатухающих колебаний представлен на рис. 4.4.
В заключительной части главы 4 обсуждаются некоторые работы, в которых немонотонная колебательная динамика микробных популяций наблюдалась в экспериментах при хемостатном культивировании плазмидосодержащих штаммов.
Рис. 4.4. Результаты численного интегрирования системы (4.4) при = 60,4, = 27, = 0,01, = 5, = 0, 0 = 21,2, = 2,5. Начальные условия: (0) = 4, x(0) = 80, (0) = 0,3;
Т – безразмерное время: Т = µ 1m t.
Математическая модель (4.1) – (4.2) использовалась в работе Koizumi J.- I. & Aiba S. (1989) для объяснения колебаний плотности популяции клеток B. stearothermophilus CU21 с рекомбинантной плазмидой pZAM26, которые наблюдались в течение 80-ти часов хемостатного культивирования. В ходе эксперимента регистрировались суммарная плотность популяции бактерий N = N0 + N1 + N2 и доля плазмидных клеток в этой популяции p = (N1 + N2) / N. Наблюдались колебания переменной N с периодом 19 часов, при этом переменная p практически не изменялась. В эксперименте обнаружены следующие характерные особенности этих колебаний: а) фаза убывания переменной N, была приблизительно в два раза продолжительнее фазы её возрастания, б) переменная p практически не изменялась и принимала значения, близкие к единице. Эти две особенности популяционных колебаний хемостатной культуры рекомбинантного штамма и его дериватов наблюдались и в серии вычислительных экспериментов на математической модели (4.1) – (4.2), см. [4].
Основные результаты и выводы 1. Установлено, что существование нижней критической плотности популяции продуцента порождает разнообразие пороговых эффектов при кратковременных воздействиях на экосистему, вызывающих резкие изменения численностей составляющих её видов. Эти эффекты прослежены с помощью трёхмерных фазовых портретов модельной экосистемы «продуцент – консумент – хищник» и объясняются множественностью аттракторов в пространстве состояний и особым расположением сепаратрисных поверхностей, разделяющих области притяжения этих аттракторов.
2. На математической модели продемонстрирована роль хищника в поддержании жизнеспособности экосистемы «продуцент – консумент – хищник». Установлено существование нижнего порога численности хищника. Промысловое изъятие хищника, понижающее его численность ниже этого порога, приводит к катастрофическим последствиям для экосистемы в целом. После такого воздействия обречены на вымирание все составляющие экосистемы – и продуцент, и консумент, и хищник.
3. Проанализированы возможные типы сосуществования трёх видов в вольтерровской модели, описывающей динамику популяций двух конкурирующих видов жертвы и популяции хищника. В зависимости от значений параметров модели эти три вида могут сосуществовать в стационарном режиме, в режиме периодических колебаний или в режиме хаотических колебаний. Механизм возникновения хаотических колебаний при изменении значений параметров связан с бифуркациями, происходящими в окрестности сепаратрисного контура, включающего в себя три особые точки, одна из которых является седлофокусом.
Установлено, что возникновение хаотического поведения при увеличении промысловой нагрузки на экосистему может быть одним из критериев приближения к опасной параметрической границе, за которой сосуществование не только трёх, но и двух видов становится невозможным ни в каком режиме.
4. Исследована математическая модель непрерывного культивирования плазмидного штамма микроорганизмов с учётом двух основных факторов нестабильности штамма – потерь плазмид при делении клеток и структурных перестроек плазмид. Модель позволяет рассчитывать область устойчивости популяции плазмидосодержащих клеток, т.е. область значений управляемых параметров процесса культивирования, при которых эта популяция устойчиво поддерживается на протоке, обеспечивая полезную биосинтетическую активность непрерывной культуры микроорганизмов.
5. Показано, что по скорости протока область устойчивости ограничена не только сверху, но и снизу некоторым не равным нулю значением скорости протока. При значениях скорости протока ниже этой пороговой величины популяция клеток исходного плазмидного штамма с течением времени вымывается из ферментёра. Существование нижней критической скорости протока подтверждается экспериментальными исследованиями хемостатного культивирования некоторых штаммов E. coli K12, содержащих рекомбинантные плазмиды.
6. На основании исследования математической модели установлено, что при постоянных условиях непрерывного культивирования популяция плазмидных клеток может поддерживаться на протоке не только в стационарном режиме, но и в режиме периодических колебаний её численности. Этот результат использовался Koizumi & Aiba (1988) для объяснения колебательной динамики численности популяции клеток в непрерывной культуре генно-инженерного штамма Bacillus stearothermophilus, несущих рекомбинантную плазмиду pZAM26.
7. Математические модели популяционной динамики, рассмотренные в диссертации, представляют собой нелинейные системы трёх обыкновенных дифференциальных уравнений. Исследование каждой из этих моделей доведено до построения параметрического портрета с указанием соответствующего набора трёхмерных фазовых портретов, качественно объясняющих наблюдаемые в эксперименте особенности динамического поведения взаимодействующих популяций. Детальное исследование этих моделей позволяет рассматривать их как базовые в иерархии моделей популяционной динамики. Эти базовые модели можно использовать при построении более сложных моделей, количественно описывающих данные экспериментов.
В заключение считаю своим долгом отметить, что автору посчастливилось быть среди учеников Александра Дмитриевича Базыкина (1940 – 1994) и мне хотелось бы выразить здесь постоянное чувство благодарности за всё, чему он меня научил. Автор глубоко признателен А.М. Молчанову и Э.Э. Шнолю за интерес к работе, В.Д. Лахно и М.Н. Устинину за поддержку работы на её завершающей стадии, Ю.М. Апонину за многолетнее и плодотворное сотрудничество, а также коллегам из лаборатории математического моделирования Института математических проблем биологии Ф.С. Березовской, Ю.А. Кузнецову и А.И. Хибнику.
Список публикаций по теме диссертации 1. Bazykin A.D., Khibnik A.I., Aponina E.A. A model of evolutionary appearance of dissipative structure in ecosystems // J. Math. Biology, 1983, v.18, № 1, p. 13 – 23.
2. Kuznetsov Yu.A., Antonovsky M.Ya., Biktashev V.N., Aponina E.A.
A cross-diffusion model of forest boundary dynamics // J. Math. Biology, 1994, v. 32, p. 219 – 232.
3. Апонин Ю.М., Антоненко Ю.Н., Апонина Е.А., Ковбаснюк О.Н., Ягужинский Л.С. Неэлектрогенный K+/H+ – обмен, включающий стадию обмена протонами между переносчиком и кислотно-основными группировками фосфолипидов на поверхности мембраны. Сравнение теории и эксперимента // Биологические мембраны, 1996, т. 13, № 3, с. 258 – 270.
4. Апонин Ю.М., Апонина Е.А. Математическое моделирование эволюции бактериальной популяции в непрерывной культуре с учётом немутационной изменчивости генома // Биофизика, 2008 (в печати).
5. Апонин Ю.М., Апонина Е.А. Иерархия моделей математической биологии и численно-аналитические методы их исследования // Математическая биология и биоинформатика, 2007, т. 2, №2, с. 347 – 360, http://www.matbio.org/downloads/Aponin2007(2_347).pdf 6. Апонин Ю.М., Апонина Е.А., Кузнецов Ю.А. Математическое моделирование пространственно-временной динамики возрастной структуры популяции растений // Математическая биология и биоинформатика, т. 1, № 1, 2006, с. 1 – 16. http://www.
matbio.org/downloads/Aponin 2006(1_1).pdf 7. Базыкин А.Д., Апонина Е.А. Модель экосистемы трёх трофических уровней с учётом существования нижней критической плотности популяции продуцента // Проблемы экологического мониторинга и моделирования экосистем. Л.: Гидрометеоиздат, 1981, т. 4, с. 186 – 203.
8. Апонина Е.А., Апонин Ю.М., Базыкин А.Д. Анализ сложного динамического поведения в модели хищник-две жертвы // Проблемы экологического мониторинга и моделирования экосистем. Л.:
Гидрометеоиздат, 1982, т.5, с. 163 – 180.
9. Апонин Ю.М., Апонина Е.А. Бифуркации в обобщённой модели Вольтерра экосистемы двух трофических уровней // Математика.
Компьютер. Образование. Сб. научн. трудов. Том 2. М. – Ижевск: НИЦ «Регулярная и хаотическая динамика», 2007, с. 131 – 138.
10. Апонин Ю.М., Апонина Е.А. О некоторых условиях устойчивого поддержания нестабильных плазмид в микробных популяциях при длительном непрерывном культивировании // Исследования по математической биологии. Пущино, 1996, с. 32 – 48.
11. Базыкин А.Д., Хибник А.И., Апонина Е.А. Нейфельд А.А. Модель эволюционного возникновения диссипативной структуры в экологической системе // Факторы разнообразия в математической экологии и популяционной генетике. Пущино, 1980, с. 33 – 47.
12. Антоновский М.Я., Апонина Е.А., Кузнецов Ю.А. Пространственно временная структура границы разновозрастного леса: простейшая математическая модель // Проблемы экологического мониторинга и моделирования экосистем. Л.: Гидрометеоиздат, 1991, т.13, с. 189 – 199.
13. Кузнецов Ю.А., Бикташев В.Н., Антоновский М.Я., Апонина Е.А.
Кросс-диффузионная модель динамики границы леса // Проблемы экологического мониторинга и моделирования экосистем. Санкт Петербург: Гидрометеоиздат, 1996, т. 16, с. 213 – 227.
14. Апонин Ю.М., Апонина Е.А., Вельков В.В. Математическое моделирование процессов непрерывного культивирования микроорганизмов, содержащих нестабильные гибридные плазмиды.
Препринт. ОНТИ НЦБИ АН СССР, Пущино, 1984, 21 с.
15. Апонина Е.А., Апонин Ю.М., Вельков В.В. Кинетические коэффициенты плазмид и методология конструирования рекомбинантных ДНК. Препринт. ОНТИ НЦБИ АН СССР, Пущино, 1984, 11 с.
16. Апонин Ю.М., Апонина Е.А., Ванякин Е.Н. Математическое моделирование процессов непрерывного культивирования с учётом гетерогенности микробных популяций. Препринт. ОНТИ НЦБИ АН СССР, Пущино, 1989, 31 с.
17. Апонин Ю.М., Апонина Е.А., Кузнецов Ю.А. Математическое моделирование пространственно-временной динамики возрастной структуры популяции растений. Препринт. ОНТИ ПНЦ РАН. Пущино, 2003, 23 с.
18. Апонин Ю.М., Апонина Е.А., Крейцер Г.П., Шноль Э.Э. Избранные алгоритмы и программы для ЭВМ МИР-2. Предельные циклы системы двух дифференциальных уравнений. Материалы по математическому обеспечению ЭВМ. Пущино, 1974, 45 с.
19. Апонин Ю.М., Апонина Е.А. Избранные алгоритмы и программы для ЭВМ МИР-2. Сепаратрисы системы двух дифференциальных уравнений. Материалы по математическому обеспечению ЭВМ.
Пущино, 1976, 36 с.
20. Antonovsky M.Ya., Aponina E.A., Kuznetsov Yu.A. Spatial-temporal structure of mixed-age forest boundary: the simplest mathematical model.
WP-89-54. Laxenburg, Austria: IIASA, 1989.
21. Antonovsky M. Ya., Aponina E.A., Kuznetsov Yu.A. On the stabiliti analysys of the standing forest boundary. WP-91-010. Laxenburg, Austria:
IIASA, 1991.
22. Апонин Ю.М., Апонина Е.А., Крейцер Г.П., Шноль Э.Э. Об использовании ЭВМ МИР-2 для качественного исследования обыкновенных дифференциальных уравнений. Теория и практика программирования на ЭВМ серии МИР. Душанбе, 1974, с.17 – 18.
23. Апонин Ю.М., Апонина Е.А. Вычисление асимптотики для периода предельного цикла при вырождении в петлю сепаратрисы. Всесоюзная конференция по асимптотическим методам в теории сингулярно возмущённых уравнений, ч. 2. Алма-Ата, 1979, с. 139 – 140.
24. Апонина Е.А., Апонин Ю.М., Базыкин А.Д. Сложное поведение в вольтерровской модели сообщества “две жертвы – один хищник”.
Тезисы докладов стендовых сообщений I Всесоюзного биофизического съезда, т. 3. Москва, 1982, с. 19.
25. Апонина Е.А., Апонин Ю.М. Предельные циклы в вольтерровской модели чисто конкурентного сообщества трёх видов. Тезисы докладов стендовых сообщений I Всесоюзного биофизического съезда, т. 3.
Москва, 1982, с. 19.
26. Апонина Е.А. Анализ математической модели непрерывного культивирования микроорганизмов, содержащих нестабильные гибридные плазмиды. Математика. Компьютер. Образование. Вып. 6.
Тезисы докладов 6-ой международной конференции. Москва, 1999, с. 264.
27. Апонин Ю.М., Апонина Е.А. Математическое моделирование пространственно-временной динамики популяции растений.
Математика. Компьютер. Образование. Тезисы докладов 11-ой международной конференции. Москва – Ижевск: НИЦ “Регулярная и хаотическая динамика”, 2004, с. 248.
28. Апонин Ю.М., Апонина Е.А. Ограниченность полутраекторий и точка максимального вырождения обобщённой вольтерровской модели экосистемы двух трофических уровней. В кн.: Математика. Компьютер.
Образование. Тезисы. Выпуск 14. Москва–Ижевск, НИЦ «Регулярная и хаотическая динамика», 2007, с. 28.
29. Апонин Ю.М., Апонина Е.А. Математическое моделирование эволюции бактериальной популяции с учётом немутационной изменчивости генома. В кн.: Математика. Компьютер. Образование. Тезисы. Выпуск 15. Москва–Ижевск, НИЦ «Регулярная и хаотическая динамика», 2008, с. 153.