Авторефераты по всем темам  >>  Авторефераты по биологии

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

ТЮТЮНОВ Юрий Викторович ПОСТРОЕНИЕ, ИССЛЕДОВАНИЕ И ПРИЛОЖЕНИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ПРОСТРАНСТВЕННО-ВРЕМЕННОЙ ДИНАМИКИ ПОПУЛЯЦИОННЫХ СИСТЕМ (03.00.02 - Биофизика)

АВТОРЕФЕРАТ

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

Красноярск - 2009

Работа выполнена в отделе математических методов в экономике и экологии Научно-исследовательского института механики и прикладной математики им. Воровича И.И.

Южного федерального университета

Научный консультант: доктор физико-математических наук, профессор МЕДВИНСКИЙ Александр Берельевич

Официальные оппоненты: доктор физико-математических наук, профессор ХЛЕБОПРОС Рем Григорьевич доктор физико-математических наук, профессор СУХИНОВ Александр Иванович доктор физико-математических наук, профессор САДОВСКИЙ Михаил Георгиевич

Ведущая организация: МГУ им М.В.Ломоносова (кафедра биофизики биологического факультета)

Защита диссертации состоится У 20 Ф апреля 2010 г. в 10 часов на заседании диссертационного совета Д 003.007.01 при Институте биофизики СО РАН по адресу: 660036, г. Красноярск, Академгородок, д. 50, стр. 50.

С диссертацией можно ознакомиться в библиотеке Института биофизики СО РАН.

Автореферат разослан У ___ Ф __________ _____ г.

Ученый секретарь диссертационного совета, кандидат биологических наук Л.А. Франк

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

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

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

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

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

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

4. построить и исследовать демо-генетическую модель пространственноЦ временной динамики популяции, позволяющую осуществлять долгосрочный прогноз эволюции генетической структуры диплоидной популяции в условиях неоднородности среды обитания. Обосновать адекватность и преимущества применения демо-генетического подхода к моделированию развития в популяции вредителя устойчивости к токсину генетически модифицированных культур при использовании стратегии Увысокая дозаЦубежищеФ. Применить модель для оценки эффективности подавления кукурузного мотылька (Ostrinia nubilalis Hubner) трансгенной Bt-кукурузой, ткани которой содержат токсичный для вредителя ген земляной бактерии (Bacillus thuringiensis);

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

Материалы и методы исследования. В качестве фактического материала использованы опубликованные в открытой печати данные специалистов Южного отдела института водных проблем РАН, Азовского НИИ рыбного хозяйства, кафедры гидробиологии биологического факультета МГУ, офиса охраны природы кантона Ву Швейцарии, Национального института сельскохозяйственных исследований Франции (INRA), а также другие материалы, опубликованные в работах авторов - специалистов в области сельского хозяйства, популяционной генетики, ихтиологии, гидробиологии и экологии, ссылки на которые приведены в списке литературы диссертации.

При построении моделей пространственноЦвременной динамики популяций использован разнообразный математический аппарат: непрерывные и дискретные модели агрегированной (точечной) динамики популяций, индивидуумЦ ориентированные модели, описывающие пространственное поведение особей, системы дифференциальных уравнения в частных производных типа реакцияЦ адвекцияЦдиффузия, адвективный член которых моделирует таксис. Демогенетические модели, построенные для прогноза генетической структуры популяций насекомых-вредителей, представляют собой системы дифференциальных уравнений в частных производных типа реакцияЦдиффузия, в которых локальная кинетика конкурирующих генотипов вредителя задается модифицированной моделью Костицына (Kostitzin 1936; 1937; 1938а, б, в), а пространственные перемещения насекомых описываются диффузией. Модификации уравнений Костицына заключаются в том, что в качестве приспособленностей генотипов рассматриваются не коэффициенты плодовитости генотипов, а коэффициенты выживаемости личинок, а также предполагается, что все экологические характеристики вредителя одинаковы для различных генотипов за исключением коэффициентов приспособленности к среде (Тютюнов и др. 2006).

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

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

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

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

Для моделей, неявно учитывающих пространственные эффекты, предложена и апробирована на данных по питанию в системе коловраткиЦмикроводоросли новая трофическая функция, обобщающая зависимости АрдитиЦГинзбурга и Холлинга II типа. Для рыбных популяций швейцарских озер и Азовского моря впервые численно решалась стационарная задача долгосрочной оптимизации промысла при минимизации риска чрезмерного падения их численности (квазивымирания).

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

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

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

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

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

Апробация работы. Результаты, полученные в рамках диссертационной работы, докладывались и обсуждались на семинаре французского центра космического мониторинга CLS УИспользование ГИС-технологий и данных спутникового мониторинга в моделировании экологических системФ (CLS, Рамонвиль Сант-Ань, Франция, 2008), на семинаре лаборатории UMR 7625 УЭкология и эволюцияФ (Университет Париж-VI, Париж, Франция, 2008), на XI-XXXVI школах-семинарах УМатематическое моделирование в проблемах рационального природопользования.

Экология. Экономика. ИнформатикаФ (Абрау-Дюрсо, Россия, 1987-2008), на международной конференции УСовременные климатические и экосистемные процессы в уязвимых природных зонах (арктических, аридных, горных)Ф (Азов, 2006), на днях науки лабораторий экологии и эволюционной паразитологии (Париж, Франция, 2006), на 3их междисциплинарных днях науки национального агрономического института ПарижЦГриньон (Гриньон, Франция, 2006), на X Европейском экологическом конгрессе Eureco'05 (Кушадасы, Измир, Турция, 2005), на рабочих совещаниях научной сети CoReV (INAPG, Париж, Франция, 2001-2005), на конференции УВлияние ГМОФ (Институт им. Пастера, Париж, Франция, 2004), на встрече NASA Tropical Team (ESSIC, Университет Мэриленда, Мэриленд, США, 2004); на 32ой международной конференции УДни энтомологовФ (INRA - Университет Ниццы, София-Антиполис, Франция, 2004), на семинаре УМатематическое понимание процессов инвазии в науках о живомФ (CIRM, Люмини, Франция, 2004), на 1ой и 2ой Международных Конференциях в Алкале по Математической Экологии AICME (Алкала, Испания, 1998, 2003), на IX Европейском Экологическом Конгрессе Eureco2002 (Лундский Университет, Швеция, 2002), на международной конференции по новым технологиям и приложениям современных физико-химических методов (ядерный магнитный резонанс, хроматография, масс-спектрометрия, ИК-Фурье, спектроскопия и их комбинации) для изучения окружающей среды (Ростов-на-Дону, Россия, 2001), на конференции УМатематические методы в экологииФ Карельского Научного Центра Российской Академии Наук (Петрозаводск, Россия, 2001), на международной конференции по детерминистическому и стохастическому моделированию биологических взаимодействий DESTOBIO (София, Болгария, 1997), на семинаре УОкуньФ (EAWAG, Кастаньенбаум, Швейцария, 1994) и др.

Благодарности. Автор выражает благодарность и глубокую признательность всем коллегам, оказавшим помощь и поддержку на разных этапах выполнении данной работы, в том числе, сотрудникам НИИМ и ПМ им. И.И. Воровича ЮФУ Ф.А. Суркову, Ю.А. Домбровскому, С.В. Бердникову, Л.И. Титовой, В.Г. Ильичеву, В.В. Селютину, Н.И. Обущенко, С.М. Хартиеву, В.Л. Шустовой, своим ученикам И.Н. Сениной, Н.Ю. Сапухиной, Е.А. Жадановской, А.Д. Загребневой, А.М.

Болдыревой, коллегам по механико-математическому факультету ЮФУ В.Н.

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

Бакаевой, офиса охраны природы швейцарского кантона Ву Б. Буттикеру, швейцарского федерального офиса внешней среды, лесов и ландшафтов Е. Штаубу, национального агрономического института ПарижЦГриньон Р. Ардити, К. Йосту, Т.

Спатаро, М.-О. Ивен, Центра биологии и управления популяциями Национального Института Сельскохозяйственных Исследований Франции Д. Бурге, Института зоологии и экологии животных Университета Лозанны О. Глязо, Д. Кантони, Ю.

Михальскому, Х. Сайа, Тулузского Университета им. П. Сабатье С. Понсард, Института системной биологии и экологии АН Чехии П. Киндлманну, Университета Стони Брук Л.Р. Гинзбургу, Р. Акчакайа, факультета прикладной математики Университета Лестера С.В. Петровскому, а также своему научному консультанту, профессору Института Теоретической и Экспериментальной Биофизики РАН в Пущино А.Б. Медвинскому.

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

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

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

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

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

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

Анализируются основные преимущества и недостатки существующих подходов к моделированию популяционной динамики. Проведено различие между моделями явно (т.е. топографически) описывающими пространство (spatially explicit models) и пространственно-неявными моделями (spatially implicit models). К первой группе относятся непрерывные модели типа клеточных автоматов и систем реакцияЦ диффузияЦадвекция, а также модели, описывающие динамику популяций в системе относительно независимых, связанных между собой миграционными потоками местообитаний (patch abundance models) и индивидуумЦориентированные модели (individual based models), описывающие пространственные перемещения отдельных особей. Ко второй группе - непрерывные и дискретные модели метапопуляций (metapopulations and metacommunities models), а также точечные модели, описывающие агрегированную динамику популяционных систем. Значительное внимание уделено уравнениям типа реакцияЦадвекцияЦдиффузия, позволяющим сочетать развитые методы аналитического исследования с мощной современной техникой вычислительных экспериментов.

Отмечается, что мелкомасштабная мозаичность популяционных систем не может быть объяснена неоднородностью факторов среды обитания. Общепринятые модели реакцияЦадвекцияЦдиффузия способны воспроизвести пространственнонеоднородную динамику микробиологических систем, однако их непосредственное использование для моделирования популяций высокоорганизованных видов, процессы воспроизводства и смертности у которых не столь интенсивны, как у бактерий, некорректно. Формулируется задача построения минимальной' математической модели пищевых миграций в системе хищникЦжертва, которая решается в первой главе на примере моделирования пространственноЦвременной динамики веслоногих раков - гарпактицид (Harpacticoida: Copepoda) - одного из основных компонентов морского бентоса на рыхлых грунтах. Для большинства гарпактицид даже в однородных биотопах характерно выраженное пятнистое (агрегированное) распределение, причем размер скоплений может варьировать от микроагрегаций площадью 0.5-1 см2 до крупных пятен площадью несколько метров (Woods, Tijtjen 1985; Fleeger, Moser 1990; Sun, Fleeger 1991; Sach, Bernem 1996;

Азовский, Чертопруд 2003; Azovsky et al. 2004).

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

Перемещения популяционной плотности гарпактицид могут быть описаны с помощью популяционной модели типа таксисЦдиффузияЦреакция. Обычно в моделях таксиса используется классическое уравнение популяционного потока ПэтлакаЦКеллерЦСегеля, основанное на предположении, что особи находятся в постоянном движении, и характер перемещений определяется частотой тамблинга (tumbling), т.е. частотой с которой особь меняет направление движения при различной концентрации некоторого стимула (Keller, Segel 1971; Okubo, Levin 2001;

Patlak 1953; Kareiva, Odell 1987). Это уравнение потока, однако, не описывает случай, когда особь может какое-то время находиться в относительном покое, для гарпактицид - это периоды, когда особь находится в грунте. Таким образом, строго говоря, использование уравнение потока ПэтлакаЦКеллерЦСегеля при моделировании организмов с подобным типом перемещения нуждается в дополнительном обосновании.

Вывод уравнения потока. Разделение двух событий: выхода рака в придонный слой воды и совершения им пространственного перемещения, позволило нам построить уравнение потока популяционной плотности, неявно учитывающее периоды пребывания особи в толще грунта, которое совпадает с уравнением ПэтлакаЦКеллерЦСегеля, но получено на основе гипотез относительно движения особей, отличных от предположений, сделанных в работах (Patlak 1953; Keller, Segel 1971; Okubo, Levin 2001). В случае одномерного местообитания уравнения потока имеет вид:

S J = (S) - (S), (1.1) x x где = (x,t) - плотность гарпактицид в точке x в момент времени t, S = S(x,t) - концентрация некоторого стимула, влияющего на интенсивность вертикальных 2 l l dP(S) миграций индивидуумов, (S) = P(S) - коэффициент диффузии, (S) = - - 2 2 dS коэффициент таксиса, P(S) - вероятность того, что в течение периода времени < особь выйдет из грунта в воду при значении стимула S, l - среднее значение квадрата длины горизонтального перемещения особи (0, ). Если выход рака из грунта в воду является простейшим пуассоновским процессом с переменным параметром (x,t) = f (S(x,t)), где f (S) - зависимость частоты выхода рака в воду от стимула S, то коэффициенты диффузии и таксиса есть:

2 l l df (S) d(S) (S) = f (S) (S) = - = -. (1.2), 2 2 dS dS Полученный результат подтверждает универсальность модели ПэтлакаЦКеллерЦ Сегеля и её применимость к описанию феномена таксиса в различных ситуациях.

ИндивидуумЦориентированная модель, построенная на основе гипотез, использованных при выводе уравнения потока (1.1) иллюстрирует, что рассмотренный механизм положительного таксиса, т.е. убывание частоты выхода рака в воду при возрастании концентрации стимула, приводит к агрегированию организмов в местах с повышенной концентрацией стимула. Как видно на рис. 1, при высоких численностях популяции динамика системы хорошо аппроксимируется эйлеровой непрерывной моделью таксисЦдиффузия (Тютюнов и др. 2009; 2010):

= - J. (1.3) t x Рис. 1. 1D распределение популяции гарпактицид (справа) в средах со стационарным распределением стимула (слева) при t = 325. Ось абсцисс - горизонтальная координата, ось ординат - концентрация стимула (слева) и обилие особей (справа), усл. единицы.

1 - получено с помощью непрерывной модели, 2 - получено с помощью индивидуумЦ ориентированной модели. Параметры: частота выхода в воду f (S) = exp(- S / 2), период времени = 1, количество шагов T = 325, количество особей M = 2000, среднее значение квадрата длины перемещения l = 2 =, распределение стимула в случае непрерывной модели S(x) = 8 exp(- (x - 20) /16). В начальный момент все особи находились в начале координат.

Модель гарпактицидыЦмикроводоросли строится как система хищникЦжертва R R 2 R = rR1- - aR + , R t K x (1.4) = - (S) S + (S) , t x x x где R(x,t) - плотность популяции жертвы (водорослей), r - коэффициент прироста жертвы, K - емкость среды, a - поисковая эффективность хищника, R, (S) - коэффициенты диффузии жертвы и хищников. Поскольку плотность микроводорослей предпочитаемого раками размера невелика (Azovsky et al. 2005) и насыщения рациона гарпактицид не происходит, их трофическая функция аппроксимирована зависимостью ЛоткиЦВольтерра aR. Такая функция типична для многих ракообразных (Jeschke et al. 2004), в том числе - для гарпактицид (De Troch et al. 2007). В модели отсутствует описание рождения и гибели хищников, поскольку демографические процессы в популяции хищника протекают намного медленнее, чем в популяции жертвы. Для замкнутого местообитания это означает, что средняя (по пространству) плотность хищников не зависит от времени, и является параметром системы.

Принимая во внимание (1.2), необходимыми условиями положительного хемотаксисного ответа являются: (S) > 0 для S 0; d(S) dS < 0; и (S) > 0 для , что хорошо согласуется с законом ВебераЦФехнера.

S > 0. Отсюда (S)0S Но какова должна быть природа стимула для того, чтобы в данной модели реализовались неоднородные динамические режимы? Показано, что общепринятая гипотеза о том, что стимулом таксиса хищника является неоднородность распределения жертв, не позволяет получить неоднородные режимы в (1.4), так как однородное распределение устойчиво при любых допустимых значениях параметров модели, для которых оно положительно.

Если стимулом таксиса является некоторое выделяемое жертвой вещество S(x,t), являющееся аттрактантом для хищников, то модель (1.4) принимает вид:

R R 2 R = rR1 - - aR + , R t K x S = - (S) + (S) , (1.5) t x x x x S 2S = R -S + , S t x R S = = = 0, (1.6) x x x x=0,L x=0,L x=0,L где - коэффициент, характеризующий интенсивность выделения химического вещества жертвой, - коэффициент распада химического вещества, S - коэффициент диффузии. Предварительно обезразмерив модель (1.5)-(1.6) таким образом, что r = K = a = = L = 1, и линеаризовав ее в окрестности нетривиального * * однородного режима (R*,, S )= (1-,,(1- ) ), мы определили аналитическое условие колебательной потери устойчивости k -ой моды неоднородного пространственного возмущения k k k (r(x,t), n(x,t), s(x,t)) = e t cos kx, e t cos kx, e t cos kx :

rk nk sk k * k * (S )> (,, k,R, (S ),k )> 0, (1.7) S -* 2 * (,,k,, (S ),S )= ( R*) [(k ) ( + S )( (S*)+ S (S*)+ S + (S )) R R R R 2 * * + (k ) ((R* +)( (S )+ 2(S )( + S )+ S )+ ( + )(S R* + )) R R R S R 2 -1 -* +((R* +) (S )+(R* +)(S R* + )+ ( + S )R*)]+ (k ) (R* +) .

R R * Превышение коэффициентом таксиса (S ) критического значения влечет бифуркацию ПуанкареЦАндроноваЦХопфа: в окрестности колебательно теряющего * * устойчивость стационарного однородного режима (R*,, S ) рождается периодический пространственноЦнеоднородный волновой режим. Возникающий режим может быть довольно сложным, если первой теряет устойчивость мелкомасштабная мода (например, k = 2, 3 или 4 на рис.2 при достаточно высоких значениях коэффициента распада ).

R = 0.01, S = 0.002, (S*) = 0.001, (S*) = 2 4 I 8 II 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 N Рис. 2. Линии, разделяющие области затухания (I) и возбуждения (II) k -ой моды, построенные для первых десяти мод.

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

Однако, возможна и иная, более предпочтительная для изучаемой биологической системы интерпретация модели (1.5)-(1.6). А именно, предположим, что частота выхода раков в воду f (S) и, соответственно, коэффициенты (S) (S) и, зависят от их сытости, характеризующейся средней наполненностью желудка хищника S(x,t).

Тогда динамика системы гарпактицидыЦмикроводоросли описывается моделью:

R R 2R = rR1- - aR + , R t K x S = - (S) + (S) , (1.8) t x x x x S = eaR -S, t где e - коэффициент УусвоенияФ пищи хищником, eaR - количество пищи, которое в единицу времени в среднем попадает в желудок хищника, - коэффициент переваривания, S - количество перевариваемой в единицу времени пищи. В отличие от (1.5) в данной модели имеет смысл положить = 0, поскольку S количество пищи в желудке хищника в точке x не влияет на сытость хищников в соседних точках.

Так как модель (1.8), (1.6) эквивалентна модели (1.5), (1.6) при = ea и = 0, то S для нее справедливы результаты линейного анализа, представленные выше. Обе модели учитывают инерционность (запаздывание) реакции хищника на изменения плотности жертв: хищники реагируют на изменение концентрации стимула (сытости), который изменяется в зависимости от плотности жертв. Именно благодаря этой задержке реакции хищника на градиент распределения жертв, при выполнении условия (1.7) возникают пространственно неоднородные режимы.

Поскольку стимул является потенциалом поля скорости таксиса: v = S, уравнение для S можно представить в виде v t = R -v + v. Таким образом, S фактически, главное отличие предложенного метода моделирования таксиса хищника заключается в гипотезе о зависимости ускорения таксиса, а не скорости, от градиента плотности жертв. При малых надкритичностях модель (1.5)-(1.6) оказывается эквивалентной еще более простой модели (см. Главу 2) с постоянными коэффициентами диффузии и таксиса хищника (Говорухин и др. 2000; Тютюнов и др. 2001; 2002; Arditi et al. 2001).

Волны убеганияЦпреследования. Предложенный подход к моделированию таксиса позволяет описать пространственно-неоднородные режимы в системе хищникЦжертва, индуцируемые исключительно пространственным поведением животных, что невозможно при использовании общепринятых двухкомпонентных моделей кроссЦдиффузии с постоянными коэффициентами (Murray 1993, 2003;

Березовская и др. 1999; Цыганов и др. 2007; Tsyganov et al. 2003, 2004; Tsyganov, Biktashev 2004). Модель (Tyutyunov et al. 2007) имеет вид:

t = -div( v ) + ;

P t = -div(P v ) + P;

P P (1.9) v t = - P + v ;

v vP t = P + vP vP, P n = n = vP = v = 0. (1.10) Здесь и P - плотности жертв и хищников; v и vP - скорости их таксиса; , , P , - диффузионные коэффициенты; , - коэффициенты таксиса. Оператор v vP P Лапласа для скалярных переменных, P имеет вид = div grad, а для векторного поля v : v = grad divv - rot rotv. Отсутствие потоков популяционных плотностей на границе области означает постоянство средних популяционных плотностей и P, которые могут быть масштабированием переменных :=, P := P P приведены к единице. Линейный анализ однородного стационарного режима * ( = 1; P* = 1; v* = 0; v* = 0) дает условие возбуждение k -ой моды разложения в ряд P Фурье малого неоднородного по пространству возмущения:

k ( + )( + )( + )( + )( + )( + ) v vP P vP P v vP v P >. (1.11) P ( + + + )v vP P С ростом произведения моды с номерами k = 1, 2,Е возбуждаются P последовательно. Включение в уравнения скоростей сил трения в виде -v и -vP может изменить порядок возбуждения мод, что делает модель более реалистичной.

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

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

Модель стаеобразования. Предложенный метод позволяет решить проблему, сформулированную в работе (Edelstein-Keshet et al. 1998) и заключающуюся в том, что стандартные модели таксиса не позволяют воспроизвести формирование и перемещение устойчивого скопления животных, индуцированного исключительно неоднородностью распределения популяции и обладающего реалистичным профилем - постоянной плотностью внутри группы и четким границами снаружи.

Наиболее известной моделью агрегирования организмов является модель КеллерЦ Сегеля (Keller, Segel 1971). Постулируя пропорциональность скорости таксиса градиенту выделяемого животными аттрактанта, v = S, она позволяет описать направленное перемещение организмов, но в случае, если размерность пространства превышает единицу, прогнозирует взрыв популяционной плотности [возможность взрыва в одномерном случае была исключена сравнительно недавно (Osaki, Yagi 2001; Hillen, Potapov 2004)]. Могильнером и Эдельштейн-Кешет (Mogilner, EdelsteinKeshet 1999) предложена более реалистичная модель, представляющая собой систему интегро-дифференциальных уравнений, в которой взаимное притяжение и отталкивание организмов в каждой точке определяется нелокальной информацией о плотности популяции.

1500 150 1500 1РЯД 1, x РЯД 1, y 1000 10 120 1R = 0.79 R = 0.500 590 0 -5-5-10-10-1500 -1500 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 800 8РЯД 2, y РЯД 2, x 1144R = 0.11R = 0.-4-4-800 -800 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 Рис. 3. Средние значения ускорения dv dt () и плотность роящихся мошек () в проекции на горизонтальную плоскость (x, y). УРяд 1Ф и УРяд 2Ф соответствует данным (Okubo et al. 1977).

Начало координат соответствует центу масс роя, единичное деление равно стандартному отклонению координат мошек. Пунктир показывает приближение данных уравнением (1.12).

R - коэффициент множественной регрессии.

Нами (Tyutyunov et al. 2004) предложена альтернативная локальная модель, отличающаяся гипотезой о зависимости ускорения (а не скорости) направленного перемещения организмов (аутотаксиса) от градиента их пространственного распределения. Данная гипотеза базируется на результатах кинематического анализа движения роящихся мошек (Okubo et al. 1977). Как видно на рис. 3, среднее ускорение движения мошек максимально на краю роя (где градиент плотности популяции максимален) и равно нулю в центре роя (где градиент нулевой).

Простейшим соотношением, связывающим ускорение и градиент плотности, является линейная зависимость dv dt = P. А для того, чтобы учесть силы отталкивания, сменяющие силы притяжения, когда плотность организмов превышает некоторое предпочитаемое ими значение Pc, следует использовать не Плотность популяции, P Среднее ускорение d v dt, см /с постоянный, а зависящий (в простейшем случае линейно) от плотности коэффициент (1- P Pc ). Получаем уравнение:

dv dt = (1- P Pc )P, (1.12) хорошо аппроксимирующее данные наблюдений Окубо и имеющее лишь два параметра.

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

P t + div(Pv) = P, P (1.13) v t = (1- P Pc -v + v.

)P v n P = n v = 0, (1.14) x x Показано, что если коэффициент аутотаксиса превышает критическое значение 2 c c,mn = P(vkmn +)P(PP, (1.15) - P) c то происходит монотонная потеря устойчивости однородного режима (P* = P, v* = 0) относительно моды малого неоднородного возмущения с волновым числом kmn. Данный результат получен для областей с различной геометрией:

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

А) t = Б) t = 2.2 Pmin = 0.0 Pmin = 510-Pmax = 25.Pmax = 19.3В) Г) t = 5.7t Pmin = 1.910- Pmin = 110-Pmax = 19.4Pmax = 19.Рис. 4. Численный эксперимент моделирующий процесс формирования и дальнейшего перемещения скопления организмов в случае сингулярного возмущения однородного распределения. Параметры: = 3.7, = 0.5, v = 0.05, = 0, Pc = 10, P = 1, L1 = L2 = 1.5.

P Как показывает рис. 4, несмотря на свою простоту, модель (1.13)-(1.14) не только позволяет реалистично описать образование плотных скоплений организмов, но и демонстрирует свойства, согласующиеся с такой особенностью пространственного поведения многих животных, как следование вдоль стены [wall-following behaviour, см.: (Creed, Miller 1990; Jeanson et al. 2003)].

В работе обсуждается связь предложенного метода моделирования таксиса с классической моделью КеллерЦСегеля. В частности, отмечается, что система (1.13) может быть использована для моделирования различных механизмов образования скоплений: и бактериального хемотаксиса, и аутотаксиса высокоразвитых видов.

Глава 2. Таксис в моделях трофических систем Модель трофотаксиса хищника. При небольших надкритичностях бифуркации, происходящей при выполнении условия (1.7), когда и отклонения переменных модели (1.5)-(1.6) от равновесных значений, и их градиенты невелики, она может быть приближена еще более простой моделью трофотаксиса в системе хищникЦжертва (Говорухин и др. 1999, 2000; Тютюнов и др. 2001; 2002; Arditi et al.

* 2001), с постоянными коэффициенты диффузии и таксиса хищника: = (S ), P * = (S ):

t = r (1- K)- a P + , (2.1) P t + div(Pv) = P, (2.2) P v t = + vv, (2.3) n v = n = n P = 0, x . (2.4) Для одномерного (отрезок [0, 1]) и двумерного (прямоугольник [0, L][0, 1]) местообитаний исследованы как минимальная модель (2.1)-(2.4), так и ее обобщения, в том числе, случай учета трения в уравнении скорости таксиса:

v t = -v + vv, (2.5) а также случай неконсервативной модели, с описанием процессов воспроизводства и смертности в уравнении баланса плотности популяции хищника:

P t + div(Pv) = ea P - P + P, (2.6) P где e - коэффициент эффективности хищничества, - естественная смертность.

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

Утверждение 2.1. Если = 0, т.е. процесс поиска жертв является случайным, то * ненулевой однородный стационарный режим модели трофотаксиса (, P*, v* 0) устойчив к малым пространственным возмущениям при любых допустимых значениях параметров (в том числе при любых значениях диффузионных коэффициентов).

Утверждение 2.2. Если > 0, т.е. хищник способен к направленному поиску и преследованию жертв, то для любого набора параметров существует критическое * значение коэффициента трофотаксиса > 0, превышение которого приводит к колебательной потере устойчивости однородного стационарного режима * (,P*,v*) и возникновению пространственно-неоднородных волновых решений.

k(a) (б) kkkkkkkk10 kКоэффициент эффективности поиска, a Рис. 5. Критические кривые устойчивости первых пяти мод однородного стационарного режима * модели трофотаксиса (,P*,v*), полученные аналитически и рассчитанные для L = 5, = 0.01, r = K = e = 1, = 0, v = 0.001. а) = 0.01, = 0.6 ; б) = 0.005, = 0.001.

P P На рис. 5 приведены критические кривые устойчивости однородного равновесия * (,P*,v*) на плоскости параметров (a, ) для двумерной неконсервативной модели с трением (2.1), (2.6), (2.5), (2.4). Точки, лежащие ниже критической кривой моды возмущения с волновым вектором k = (n L, m ), соответствуют устойчивости nm однородного режима, выше - неустойчивости. В случае рис. 5а моды возбуждаются в порядке возрастания длины волнового вектора k. Заметим, что при уменьшении nm коэффициентов диффузии порядок возбуждения мод может измениться (рис. 5б).

Потеря устойчивости однородного режима может сопровождаться возбуждением как квазиодномерной ( n = 0 или m = 0 ), так и истинно двумерной моды (n, m 0).

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

Одним их важных результатов, полученных для консервативной модели (2.1)(2.4), является то, что при установлении автоколебательного пространственнонеоднородного режима, диапазон сосуществования хищника и жертвы расширяется за счет пространственной активности хищника. Другими словами, даже при значениях параметров модели, соответствующих области устойчивости нулевого равновесия переменной плотности жертвы в точечной модели (например, при P > r a ), пространственная модель (2.1)-(2.4) демонстрирует стабильное сосуществование двух видов в неоднородном неравновесном режиме (Говорухин и др. 1999; Arditi et al. 2001). Таким образом, модель объясняет миграционный Коэффициент трофотаксиса механизм преодоления дефицита пищевого ресурса, основанный на периодическом посещении хищником пятен с высокой плотностью популяции жертвы, играющий важную роль в поддержании гомеостаза природных экосистем, в частности - в рассмотренной в Главе 1 литоральной трофической системе гарпактицидыЦмикроводоросли (Azovsky et al. 2005 a,b).

Возникновение интерференции хищников. Бифуркация, происходящая при * превышении коэффициентом таксиса критического значения , имеет место даже при использовании простейшей трофической функции ЛоткиЦВольтерра g( ) = a.

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

T def 1 < g(, P)P > < g(, P)P > Q = lim dt = < >>, (2.7) T T < > < > показал, что в системе возникает эффект интерференции хищников, проявляющейся на популяционном уровне как зависимость среднего рациона хищника от их численности.

Рис. 6. Зависимость агрегированной смертности жертв обусловленной хищничеством от размера популяции хищников, рассчитанная при L = , v = 0.00001, = 0.03, = 0.02 в сравнении с P теоретическими кривыми точечных моделей хищникЦжертва с трофическими функциями ЛоткиЦ Вольтерра (пунктир) и АрдитиЦГинзбурга (сплошная линия).

Обнаружено (рис. 6), что зависимость Q(P), где P - агрегированная по пространству плотность хищника, является линейной лишь в однородном случае, * т.е. при < (P). При установлении и продолжении по параметру P пространственноЦнеоднородного режима модели (2.1)-(2.4), зависимость Q(P) деформируется, принимая вид нелинейной функции, насыщающейся с увеличением P. Подобным образом насыщается, в частности, зависимость Q(P), рассчитанная для RD (ratioЦdependent) трофической функции АрдитиЦГинзбурга a P g( P) = (Arditi, Ginzburg 1989), представляющей один из возможных 1+ ah P способов учета зависимости рациона хищников от их численности в точечных моделях (см. Главу 3).

Парадокс биологического контроля. Благодаря обусловленному пятнистостью популяций эффекту интерференции, предложенный метод явного моделирования активного пространственного поведения хищника позволяет разрешить известный Упарадокс биологического контроляФ, заключающийся в невозможности при помощи классических методов моделирования трофических систем воспроизвести наблюдаемую в природе стабильную динамику системы хищникЦвредитель при низкой численности популяции вредителя (Luck 1990; Arditi, Berryman 1991).

Количественный критерий адекватности модели хищникЦжертва предложен в работе (Beddington et al. 1978): модель должна воспроизводить стабильную динамику численности вредителя на уровне, не превышающем 2.5% от численности, достигаемой в отсутствие хищника. В рамках точечных моделей, описывающих агрегированную динамику систем хищникЦжертва, разрешение парадокса биоконтроля возможно путем модификации классических, зависящих лишь от плотности жертв, трофических функций g( ) таким образом, чтобы учесть в них эффект оказываемый интерференцией хищников на их рацион g(, P). Это позволяет неявно учесть в точечных моделях пространственные эффекты, стабилизирующие природные системы. В частности, проблема решается при использовании упомянутой выше трофической функции АрдитиЦГинзбурга.

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

Для системы хищникЦвредитель построена модель трофотаксиса, локальная кинетика в которой описывается классической системой РозенцвейгаЦМакАртура с логистическим воспроизводством жертв и трофической функцией Холлинга II типа:

a P = (1- )- + ;

t 1+ ah P a P = - P - div(Pv)+ P; (2.8) P t 1+ ah v = + vv;

t n v = n = n P = 0, x . (2.9) Изучена обусловленная трофотаксисом возможность образования пространственных структур из устойчивого в отсутствие пространственных эффектов однородного по пространству стационарного режима:

a(1- h)- 1+ h * (, P*,v*)=, 0, где < a <, (2.10) a(1- h), 1- h h(1- h) a2(1- h) а также механизм возникновения пространственной неоднородности из устойчивого в отсутствие пространственных эффектов однородного периодического режима C0, являющемся аттрактором модели РозенцвейгаЦМакАртура при a > (1+ h) (h - h2).

При больших a режим C0 характеризуется высокой амплитудой - численность жертв периодически достигает значений близких емкости среды. Именно наличие такого режима в классических моделях послужило причиной возникновения парадокса биоконтроля.

Сформулируем условие устойчивости равновесия для случая 1D модели:

Теорема 2.1 (условие устойчивости). Однородный стационарный режим * (, P*, v* ) модели (2.8)-(2.9) устойчив к малым пространственно неоднородным возмущениям вида (n(x,t), p(x,t),v(x,t)) = et cos kx, pk et cos kx, et sin kx тогда nk vk k k k и только тогда, когда для каждого волнового числа k = n L, n 0, выполняется неравенство 1 a + - a 0 < (v + )(( + v )k -a11)+, (2.11) P P 2 P* k k где a11 < 0, a21 > 0 - элементы Якобиана системы (2.8).

Следствием данной теоремы являются утверждения, аналогичные приведенным выше для несколько более простой модели трофотаксиса (2.1)-(2.6), о * возникновении адвективной неустойчивости режима (, P*, v* ) при достаточно высоком значении .

Рис. 7. Неустойчивый однородный периодический режим C0 и устойчивый неоднородный периодический режим R, индуцируемый = 0.5, представлены на плоскости осредненных по пространству численностей популяций. Параметры модели: a = 0.7, h = 2.857, L = 5, = 0.01, = 0.01, = 0.6, v = 0.001.

P Влияние трофотаксиса на устойчивость однородного периодического режима C0, существующего при a > (1+ h) (h - h2), исследовалось численно, путем определения положения мультипликаторов оператора монодромии на комплексной плоскости при последовательном увеличении коэффициента таксиса. Таким образом, были * определены бифуркационные значения , приводящие к адвективной потере устойчивости однородного режима C0 и возникновению неоднородных по пространству колебаний R.

Пространственное осреднение Локальные колебания a б в Рис. 8. Фазовые траектории и локальная динамика (x = 20) режимов модели (2.8)-(2.9):

(а) = 0.5 - неоднородный периодический режим; (б) = 1 - неоднородный квазипериодический режим; (в) = 1.5 - неоднородный хаотический режим со вспышками плотности жертв. При = 0 устойчив пространственно-однородный периодический режим C0 : < >[0.00001; 0.62], < P >[0.51; 3.06]. Параметры приведены в подписи рис. 7.

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

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

(а) (б) Рис. 9. Распределение плотности популяции жертвы на двумерном местообитании для (а) периодического режима, реализующегося при возбуждении моды с волновым вектором s( = 0.00069273), (б) пространственноЦвременного хаоса, возникающего при продолжении режима по параметру в надкритическую область ( = 0.0009). Другие параметры модели:

Lx = 4, Ly = 3, a = 0.06, h = 10, = 0.01, = 0.001, = 0.001, v = 0.001.

P Интересно, что и более простая модель (2.1), (2.6), (2.5), (2.4), отличающаяся от системы (2.8)-(2.9) трофической функцией, описывающей локальные взаимодействия жертв с хищниками (линейная зависимость ЛоткиЦВольтерра вместо насыщающейся функции Холлинга II типа, а), позволяет получить сложные динамические режимы, качественно совпадающие с представленными на рис. (Тютюнов и др. 2002). Это объясняется тем, что реальный контакт хищников и жертв имеет место на границах пятен их популяций, где плотность жертв не высока и насыщения функции Холлинга в большинстве случаев не происходит.

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

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

Таблица 1. Примеры трофических функций, с учетом и без учета интерференции хищников.

Название Выражение Источник (авторы) ЛоткиЦВольтерра (LV) g( ) = a Lotka (1925), Volterra (1926) Холлинга II-го типа g( ) = a (1+ ah ) Holling (1959) Ивлева g( ) = gmax(1- e-s ) Ивлев (1955), Ivlev (1961) ХасселаЦВарлей (HV) g(, P) = Pm Hassell, Varley (1969) Pm ХасселаЦВарлейЦ Sutherland (1983), g(, P) = Холлинга (HVH) Arditi, Akakaya (1990) 1+ h Pm БеддингтонаЦ Beddington (1975), g(, P) = a (1+ awP + ah ) ДеАнжелиса (BDA) DeAngelis et al. (1975) АрдитиЦГинзбурга (RD, Гинзбург, Гольдман, Раилкин P g(, P) = g( P) = 'ratio-dependent' (1971), Arditi, Ginzburg (1989) 1+ h P зависимость) Contois (1959) P Гибридная модель 1- g(, P) = (1- ) Trn (2008) Удележа жертвФ Трана P Гибридная модель g(, P) = (1- e-P ) Trn (2008) Увыедания жертвФ Трана P a Обобщенная RD g(, P) = Tyutyunov et al. (2008) P P0 + exp(- P P0)+ ah зависимость (GRD) a g(, P) = Вариант GRD функции Тютюнов и др. (2010) P P0 +1 (1+ P P0 )+ ah Несмотря на значительный объем наблюдений обратного влияния численности хищников на их рацион, пока нет согласия относительно того, какова должна быть общая форма трофической функции, учитывающей это влияние. Практическое значение данной теоретической проблемы очевидно - выбор трофической функции в корне меняет динамические свойства моделей. Зависимость рациона от численности хищников позволяет, неявно учитывая в точечных моделях эффекты пространственной неоднородности и поведения хищников, стабилизировать их динамику и повысить адекватность прогноза. Существующее разнообразие моделей интерференции, однако, не способствует ясности в понимании ее роли в природных системах.

Качественный анализ точечных моделей КолмогороваЦГаузе вида d dt = F( ) - g(, P)P;

(3.1) dP dt = eg(, P)P - P, для случаев мальтузианского и логистического роста популяции жертв и двух вариантов трофической функции, BDA и HVH (см. табл. 1), позволил сделать вывод о том, что как слишком низкие, так и слишком высокие значения параметров интерференции в этих моделях ( w и m ) оказывают скорее неблагоприятное действие на динамику трофической системы, либо дестабилизируя нетривиальное равновесие, либо сильно уменьшая равновесное значение численности хищников (Arditi et al. 2004). Высказывается гипотеза, что эволюционным преимуществом могут обладать виды с умеренной интерференцией, согласующаяся с тем, что эмпирические оценки параметра m функции HVH в работе (Arditi, Akakaya 1990) для многих трофических систем оказались близкими к 1.

Оценка трофической функции. Явные пространственные модели поведения хищников могут объяснить и проиллюстрировать механизм возникновения эффекта интерференции в распределенных трофических системах. В Главе 2, при помощи модели трофотаксиса (2.1)-(2.4), было продемонстрировано, как на популяционном (макро-)уровне функционирования системы хищникЦжертва проявляется влияние численности хищников на их средний рацион (агрегированную трофическую функцию). Однако, поскольку в построенной нами динамической модели численность жертв постоянно менялась, непосредственная оценка значений агрегированной трофической функции g(, P) при различных комбинациях численностей жертв и хищников была невозможна - наличие интерференции определено косвенно, путем вычисления агрегированной смертности жертв (2.7), обусловленной хищничеством. Непостоянство популяционных численностей затрудняет определение рациона хищников и в теоретических (Arditi et al. 2001), и в лабораторных (Arditi, Saah 1992; Jost, Ellner 2000; Skalski, Gilliam 2001; Tully et al.

2005), и в полевых исследованиях (Jost, Arditi 2000; Jost et al. 2005). ИндивидуумЦ ориентированная модель, явно описывающей перемещения и взаимодействие особей в системе хищникЦжертва (Tyutyunov et al. 2008), позволяет преодолеть эту трудность и изучить влияние пространственного поведения на вид трофической функции хищника.

В начальный момент случайным образом задаются пространственные координаты жертв и P хищников внутри прямоугольного местообитания = [0, Lx ][0, Ly]. Изменение положения j -го хищника xP за временной шаг j,t описывается уравнением:

P xP = xP + + vP, ( j = 1,Е, P).

j,t +1 j,t j,t j,t P Здесь - случайная величина, равномерно распределенная внутри круга радиуса j,t P , характеризующего интенсивность ненаправленных перемещений, а скорость P P P P таксиса vP вычисляется как vP = St (xP )+ vP, где 0 и [0;1] - j,t j,t j,t j,t-коэффициенты таксиса и инерции, а St (xP ) - градиент УзапахаФ жертв, j,t являющегося аттрактантом для хищников. Предполагается, что запах каждой жертвы распределен нормально:

x - xi,t Si,t (x) = exp-, 2 ( 2 ) а запах всей популяции вычисляется как сумма St (x) = Si,t (x). Пространственные i=перемещения жертв описываются аналогичным образом, с возможностью включения в модель их отрицательного таксиса против градиента запаха хищников.

Жертва, оказавшаяся в R -окрестности хищника съедается с вероятностью 1- .

При различных предположениях относительно пространственной активности видов (различных соотношениях их коэффициентов , и ), перемещения индивидуумов отслеживались в течение T итераций (t = 1, 2,Е,T; T = 200 ) с подсчетом среднего количества жертв, пожираемых хищником в единицу времени.

Эта процедура выполнялась 500 раз для всевозможных сочетаний количества жертв и хищников (, P = 1, 2,Е,25 ). Постоянство численностей популяций обеспечивалось тем, что съеденная жертва сразу воскресала' в случайно выбранной точке.

Аналогичный прием с заменой жертвы используется и в лабораторных экспериментах по определению рациона насекомых (Tully et al. 2005). Результатом счета является таблица средних рационов G = (g ), представляющих собой,P,P=значения трофической функции g(, P), которая не задается a-priory, а является результатом динамики системы.

Эксперименты показали, что для случайно перемещающихся животных трофическая функция не зависит от численности хищников, однако, при способности хищников к активному преследованию случайно перемещающихся P жертв ( > 0, = 0 ) наблюдается уменьшение рациона хищников с увеличением размера их популяции - типичное проявление интерференции (см. рис. 10).

Полученные значения g аппроксимировались теоретической зависимостью HVH,P (см. табл. 1). Изучено влияние стратегии пространственного поведения хищника P P (параметров , ) на эффективность поиска , время обработки жертвы h и степень интерференции m в функции HVH.

Полученные результаты соответствуют результатам непрерывной модели (Глава 2), и подтверждают гипотезу о том, что эффект интерференции проявляется лишь при определенных соотношениях популяционных плотностей хищников и жертв, а именно, - как при низких плотностях хищников, так и при недостатке жертв, рацион хищников адекватной трофической функции не должен зависеть от P (Abrams, Ginzburg 2000; Berdnikov et al. 1999; Ginzburg, Jensen 2008). Изотрофы такой зависимости будут почти вертикальны при малых P и. В отличие от BDA-, HVH- и RD-функции, данным свойством обладает предложенная нами GRD-функция (см.

табл. 1), обобщающая RD-зависимость АрдитиЦГинзбурга на случай низких популяционных плотностей, и хорошо приближающая результаты имитационных экспериментов. Для данных лабораторных экспериментов Бакаевой (1999), по определению рационов коловраток (Brachionus calyciflorus и Philodina acuticornis) в монокультуре микроводорослей (Chlorella, Scenedesmus и Synechocistis), GRDзависимость также обеспечивает наилучшую в сравнении с вышеназванными функциями аппроксимацию (Тютюнов и др. 2010). Вместе с тем показано, что оригинальная RD-функция, имеющая всего два параметра, тоже достаточно хорошо описывает трофические взаимодействия в данной системе.

а б P P Рис. 10. Изотрофы g(, P) = const (вверху), полученные в экспериментах с (а) = 0.2, = 0.P P (оцененные коэффициенты HVH-функции - = 0.1; m = 1.0 ; h = 25.7 ) и (б) = 0.9, = ( = 0.006 ; m = 0.31; h = 18.4 ). Внизу представлены типичные для выбранных коэффициентов P P , распределения запаха жертв, а также 25 жертв () и 25 хищников ().

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

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

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

В обзорной части главы рассматриваются три подхода к описанию изучаемой биологической системы: (i) сложные имитационные пространственные модели, использующие чрезвычайно подробные предположения о популяционной генетике и жизненном цикле видов насекомых (Peck et al. 1999; Guse et al. 2002; Ives, Andow 2002; Storer et al. 2003; Heimpel et al. 2005), (ii) модели, получаемые путем формального добавления диффузионного члена к классическим уравнениям математической генетики (модель Фишера-Холдена-Райта), концентрирующиеся на процессах, происходящих на генетическом уровне, но пренебрегающие экологией насекомых (Alstad, Andow 1995; Vacher et al. 2003; Cerda, Wright 2004;

Tabashnik et al. 2004), (iii) демо-генетические диффузионные модели, учитывающие и демографию, и генетику изучаемых видов (Hillier, Birch 2002а, б; Richter, Seppelt 2004; Жадановская и др. 2006; Тютюнов и др. 2006; Medvinsky et al. 2008).

Отмечается, что оценки развития Bt-устойчивости насекомых-вредителей, получаемые с помощью детальных имитационных моделей, выглядят достаточно реалистично, но огромное число параметров [до 41000 в работе (Onstad 1988)] затрудняет их идентификацию и снижает ценность. Подход, основанный на расширении уравнений Фишера-Холдена-Райта описанием диффузионного распространения алели устойчивости, приводит, по крайней мере, к двум проблемам (Абросов, Боголюбов 1988): (i) область адекватного применения таких уравнений ограничена биологическими видами с чередованием гаплоидных и диплоидных поколений, к которым не относятся насекомые (в том числе, кукурузный мотылек), (ii) центральные гипотезы отсутствия мутаций, отбора и миграций в популяции противоречат характеру a-priori неравновесного переходного процесса развития Btустойчивости.

Предложенный подход к прогнозированию эволюции Bt-устойчивости в популяциях насекомых-вредителей основан на явном описании пространственного распространения плотностей конкурирующих генотипов вредителя, локальная динамика которых моделируется демо-генетическими уравнениями Костицына (Kostitzin 1936; 1937; 1938а, б, в; Свирежев, Пасеков 1982; Абросов, Боголюбов 1988). Предполагается, что ген, отвечающий за наличие признака Bt-устойчивости у отдельной особи вредителя, может находиться в одном из двух состояний (называемых аллелями): в состоянии Bt-восприимчивости (s аллель) или Btустойчивости (r аллель). Эти две аллели формируют три генотипа вредителя: Btвосприимчивые генотипы ss и rs (при рецессивности гена Bt-устойчивости) и Btустойчивый генотип rr. Модель - система УреакцияЦдиффузияФ:

ss t = Wssb( ss + rs 2)2 - ss - ss + ss, ) t = Wrs 2b( + 2)( 2 + - - + , rs ss rs rs rr rs rs rs (4.1) t = Wrrb( 2 + - - + , ) rr rs rr rr rr rr n = n = n = 0, x , ss rs rr где = (x,t) - плотность генотипа ij в точке x в момент времени t ( или i, j = r ij ij s), = + + - общая плотность популяции; параметры b, и - ss rs rr коэффициенты плодовитости, смертности и внутривидовой конкуренции вредителя.

При моделировании стратегии Увысокая дозаЦубежищеФ полагается, что ареал вредителя состоит из участков, засеянных либо Bt-кукурузой (Bt ), либо обыкновенной кукурузой (Ref ). Приспособленности генотипов (выживаемость личинок), согласно (Bourguet et al. 2000; Vacher et al. 2003) представляется в виде:

1- hc c, x ref ; 1, x ref ;

Wrr = 1- c, на всем ; Wrs (x) = Wss(x) = (4.2) 1- + h ( - c), x Bt ; 1- , x Bt, где - коэффициент отбора по признаку Bt-устойчивости; с - цена, которую платит генотип, имеющий ген Bt-устойчивости, за преимущество, проявляющееся на Btполях; h - уровень доминантности отбора по признаку устойчивости к Bt-токсину;

hc - уровень доминантности цены c. Параметры , c, h, hc [0;1].

В случае пространственной однородности области , суммирование уравнений (4.1) дает уравнение логистического роста популяции. Представление частот аллелей через плотности генотипов ps = ( + 2), pr = ( + 2), ss rs rr rs позволяет доказать теорему о сохранении постоянства частот. Также показано, что для равновесий демо-генетической модели выполнено условие ХардиЦВайнберга ussurr = urs 4 где uij = - частота генотипа ij. Данные результаты подтверждают ij адекватность демо-генетической модели (1.4) в пространственно однородном случае.

В случае неоднородности области (сочетания Bt-участков и убежищ) справедлива следующая теорема:

Теорема 4.1. В биологической системе, описываемой пространственной демогенетической моделью (4.1), помимо диффузионного распространения аллельных частот в пространственно-неоднородной среде возникает адвективный поток аллельных частот, индуцируемый неоднородностью пространственного распределения общей плотности популяции и аллельных частот pr и ps:

pr t = bpr (Wr -W )+ pr + 2 pr ln ;

(bW ( ))+ t = - + ;

2 ps + pr = 1; Wr = Wrs ps + Wrr pr ; W = Wss ps + 2Wrs ps pr + Wrr pr ; (4.3) pr n = n = 0, x ..

Модель (4.3) эквивалентна (4.1) и отличается от классических уравнений математической генетики, формально объединенных с диффузией (пространственной модели ФишераЦХолденаЦРайта), только членом 2 pr ln, который влияет на распространение Bt-устойчивой аллели r и может быть интерпретирован как адвективный' член, описывающий направленный поток аллельной частоты pr со скоростью - 2 ln. Эта адвекция индуцируется неоднородностью пространственного распределения и pr и исчезает, если распределение или pr однородно.

Следствие. Диффузионная модель ФишераЦХолденаЦРайта адекватно описывает пространственное распространение гена Bt-устойчивости лишь в случае однородного распределения общей плотности в пределах всего моделируемого поля .

Теорема 4.2. При эволюции генетической структуры популяции, описываемой пространственной демо-генетической моделью (4.1), отклонение системы от равновесия ХардиЦВайнберга = ussurr - urs 4 (Свирежев, Пасеков 1982) описывается дифференциальным уравнением:

2 t = b(ps pr (Wss + Wrr - 2Wrs )- W)+ + 2 ln + 2 pr. (4.4) Согласно (4.4), система (4.1) стремится к равновесию ХардиЦВайнберга ( 0, t ) только, если одна из аллельных частот ( pr или ps ) стремится к нулю. В противном случае система (4.1) эволюционирует за пределами равновесия Харди-Вайнберга к полиморфному состоянию. Кроме того, отклонение увеличивается за счет неоднородности распределения аллельных частот.

ref, Bt, убежище Bt поле 0 Lx xref Рис. 11. Шаблон убежища для одномерного ареала = [0, Lx ].

Для нахождения стационарных режимов пространственной демогенетической модели (4.1) на одномерном (1D) ареале = [0, Lx], состоящем из двух участков: Bt-поля и убежища (рис. 11), решалась стационарная краевая задача следующего вида:

d 2 ss 1 b rs Wss ss + - ss - ss ;

= dx2 d rs 1 2b rs rs = - Wrs ss + 2 2 + rr - rs - rs ;

dx2 (4.5) rr d = - 1 b + - - Wrr rs rr , rr rr dx2 d dx = d dx = d dx = 0, x [0, Lx].

ss x=0,Lx rs x=0,Lx rr x=0,Lx Для обоснования корректности постановки задачи (4.5) были сформулированы и доказаны теоремы о существовании ее решений:

Теорема 4.3. Всегда существуют стационарные однородные по пространству * * * решения системы (4.1) вида (x)= 0, (x)= 0, (x) = K = (bWrr - ) , x [0, Lx], ss rs rr где Wrr = 1- c согласно формулам (4.2).

Теорема 4.4 (необходимое условие). Если существуют стационарные неоднородные по пространству решения системы (4.1) для конфигурации поля, * * * представленной на рис. 11, такие, что x [0, Lx] (x)(0, K), а (x) = (x) 0, и ss rs rr при этом Wss(x) = 0 для x Bt и Wss(x) = 1 для x ref, то выполняется условие:

2 (K - 2 3) < b , (4.6) 0 0 ref где >, = (0), = (xref ) ( xref (0, Lx ) - точка, разграничивающая Bt-поле 0 ref 0 ref и убежище).

В ходе доказательства был определен вид и способ построения решения (рис. 12): т.к. коэффициент приспособленности Wss в (4.5) определяет качественное отличие динамики системы в убежище (ref ) и на Bt-поле (Bt ), задача решалась на каждом из этих участков, затем два решения УсшивалисьФ в точке xref (точка В на рис. 13).

* Рис. 12. Стационарное решение системы (4.1) такое, что x [0, Lx] (x)(0, K), а ss * * (x) = (x) = 0. xref - точка перегиба, разграничивающая убежище (20%) и Bt-поле (80%).

rs rr Рис. 13. Фазовые портреты редуцированной системы (4.5): (x) = (x) x [0, Lx]. Здесь ss * ref * Bt Bt ref ref Bt ref, ; (б) = (Wss b - ) , = (Wss b - ) , Wss < Wss. (а) Wss > b Wss < b Wss > b, 1 Bt Wss > b. Черные линии - фазовые траектории системы в убежище, серые линии - на Bt-поле.

Пунктирные линии - сепаратрисы. Дуга AB соответствует участку замкнутой фазовой траектории, вдоль которой образующая точка движется в убежище. Дуга BC соответствует участку незамкнутой фазовой траектории гиперболического типа, вдоль которой образующая точка движется на Bt-поле. Точка А соответствует левому краевому условию задачи (4.5), точка С - правому. Точка В - точка УсшивкиФ решения в убежище с решением на Bt-поле.

Верна и более общая теорема:

Теорема 4.5 (необходимое условие). Если существуют стационарные неоднородные по пространству решения системы (4.1) для конфигурации поля, * ref представленной на рис. 11, такие, что x [0, Lx] (x)(0, (Wss b - ) ), а ss * * ref Bt ref ref (x)= (x)= 0, и при этом Wss > Wss и Wss > b (Wss = Wss(x), rs rr xref Bt Wss = Wss (x) ), то выполняется условие:

xBt ref ref Bt - 2 b(Wss -Wss ), 2 Wss b - < (4.7) 0 0 ref 3 где >, = (0), = (xref ) ( xref (0, Lx ) - точка, разграничивающая Bt-поле 0 ref 0 ref и убежище) (рис. 12).

Для численного исследования пространственно-временной динамики популяции кукурузного мотылька, использовались следующие параметры модели (4.1):

b = 9.94 год-1, = 6.84 год-1, K = 147.4 106 особей/км2. Начальное распределение популяции мотылька предполагалось однородным ( = 2.948106 особей/км2), а частота Bt-устойчивого гена - малой pr = 0.0015 (Andow et al. 2000; Bourguet et al.

2003).

Бифуркационной анализ демоЦгенетической модели показал, что и при малых, и при высоких значениях коэффициента диффузии , бифуркация в (4.1) происходит жестко (рис. 14а). При высокой подвижности вредителя ( 1 км2год-1) в системе (4.1) имеет место замкнутый гистерезисный цикл, который не позволяет мягко управлять системой за счет небольших вариаций размеров однополосного убежища: например, если при =10 км2год-1 размер убежища станет немного меньше P1, то Bt-устойчивая аллель r полностью вытеснит Bt-восприимчивую аллель s (вся популяция будет состоять только из Bt-устойчивых особей rr генотипа), и снова подавить ген Bt-устойчивости можно только при увеличении доли убежища свыше P2, т.е. при достаточно больших размерах убежища.

Фактически речь может идти о полном прекращении использования Bt-кукурузы на полях. При низкой подвижности вредителя ( 0.1 км2год-1) можно плавно управлять системой, варьируя размер убежища: например, при = 0.1 км2год-переходя критическую точку P2 при уменьшении размера убежища, всегда можно вернуться к устойчивому режиму, когда в популяции нет Bt-устойчивого гена ( pr = 0 ), слегка увеличив размер убежища. Таким образом, подвижность вредителя играет определяющую роль в эффективности стратегии Увысокая доза-убежищеФ.

Рис. 14. Частота Bt-устойчивой аллели pr в зависимости от размера убежища (в %) и коэффициента диффузии для демо-генетической модели (4.1) (а) и модели Фишера-Холдена-Райта (б). Lx = 16 км. Сплошными линиями отмечены ветви устойчивых равновесий систем, пунктирными линиями - ветви неустойчивых равновесий. Заданные генетические параметры соответствуют стратегии Увысокая дозаЦубежищеФ: = 1,, h = 0., hc = 0.2 c = 0.1 (особи, обладающие геном Bt-устойчивости, платят за свою приспособленность к токсичной среде 10% от коэффициента выживаемости Wrr ).

Прогноз диффузионной модели Фишера-Холдена-Райта (рис. 14б) только в случае нулевой и бесконечно большой подвижности вредителя ( = 0 и ) совпадает с демо-генетической моделью (4.1). Согласно фишеровской модели, эффективное применение стратегии Увысокая доза - убежищаФ с целью предотвращения адаптации мотылька к Bt-кукурузе ( pr = 0 ) требует слишком больших размеров убежищ, близких к 100% (98% от всего поля для = 0.1 км2год-1), в то время как в демоЦгенетической модели (4.1) полностью предотвратить развитие Btустойчивости у мотылька удается при гораздо меньших размерах убежища (21% от всего поля для = 0.1 км2год-1).

а) б) адвекция адвекция диффузия pr pr диффузия хref 0 хref Lx 0 Lx Рис. 15. Схема пространственного распределения частоты и общей плотности в демоpr генетической модели (4.1) на 1D ареале [см. (4.3)], когда (а) плотность Bt-устойчивых особей вредителя на Bt-поле очень низка ( 0 ) и (б) плотность Bt-устойчивых особей вредителя на Btполе возросла в результате распространения в популяции Bt-устойчивого гена.

Таблица 2. Время T10 (годы), за которое частота Bt-устойчивой аллели в популяции мотылька достигает 10% в демо-генетической модели (4.1) при разных размерах убежища и значениях коэффициента диффузии мотылька. В скобках даны значения времени для диффузионной модели Фишера-Холдена-Райта. Размер поля: 16 км 16 км. Генетические параметры соответствуют стратегии Увысокая дозаЦубежищеФ: = 1, c = 0,, hc = 0.2.

h = 0. Размер убежища (%) Диффузия, (км2год-1) 5 10 15 20 25 30 0 24 24 24 24 24 24 0.02 24 (8) 24 (9) 24 (9) 25 (9) 25 (9) 25 (9) 26 (9) 0.04 24 (8) 25 (8) 485 (9) 720 (9) 979 (9) 1262 (9) 1567 (9) 0.05 92 (8) 266 (8) 473 (9) 700 (9) 947 (9) 1213 (9) 1496 (9) 0.1 67 (8) 205 (8) 382 (8) 574 (9) 776 (9) 988 (9) 1208 (9) 0.5 20 (7) 69 (8) 140 (8) 229 (8) 329 (8) 435 (8) 545 (9) 1 13 (7) 37 (7) 78 (8) 130 (8) 192 (8) 261 (8) 335 (9) 5 8 (5) 12 (7) 20 (8) 32 (8) 48 (8) 66 (9) 86 (9) 7 7 8 8 9 9 Если нет цены за приспособленность к Bt-токсину ( c = 0 ), то предотвратить адаптацию вредителя к Bt-растениям невозможно, поэтому оценивался период времени, на который можно задержать ее развитие (табл. 2). Установлено, что при средней подвижности вредителя демо-генетическая модель (4.1) прогнозирует задержку адаптации мотылька к Bt-кукурузе в течение сотен и даже тысяч лет - результат, который не может быть получен с фишеровской диффузионной моделью, в которой время развития Bt-устойчивости при любых сочетаниях значений коэффициента диффузии и размера убежища не превышает 24 лет. Заметим, что полученные для демо-генетической модели значения времени T10 в табл. 2 вполне способны объяснить, почему, несмотря на продолжительное (более 10 лет) и интенсивное выращивание трансгенной кукурузы в США и других странах, Bt устойчивые особи мотылька до сих пор не были выявлены в природных условиях.

Таким образом, несмотря на достаточную простоту предложенной демоЦ генетической модели (4.1), она, в отличие от дифузионной модели ФишераЦ ХолденаЦРайта, позволяет воспроизвести и объяснить эффективное воздействие убежища на задержку развития Bt-устойчивости в популяции вредителя.

Демо-генетическая модель Модель Фишера-Холдена-Райта (a) t = 11 лет (д) t = 1 год K K pr = 0.0(б) t = 574 года (е) t = 9 лет K K pr = 0.0 (ж) t = 10 лет (в) t = 576 лет K K pr = 0.(з) (г) t = 45 лет t = 615 лет K K pr = 0 8 16 км 0 8 16 км ss генотип rs генотип rr генотип Рис. 16. Плотности распределения генотипов мотылька в моменты времени (годы), когда частота Bt-устойчивой аллели последовательно достигает 0.002, 0.1, 0.7 и 1, для 20% убежища (рис. 11) pr в демо-генетической и фишеровской моделях. Коэффициент диффузии = 0.1.

Единственным отличием сравниваемых моделей является наличие в системе (4.3) адвективного слагаемого 2 ln pr, возникающего в демо-генетической модели (4.1) за счет пространственной неоднородности среды. Результаты вычислительных экспериментов показывают, насколько велико может быть влияние этого члена. Именно адвективный поток, пока Bt-устойчивых особей на Bt-поле мало (рис. 15а), сдерживает диффузионное распространение Bt-устойчивого гена в популяции вредителя (рис. 16а). После того как число Bt-устойчивых особей становится достаточно большим (рис. 15б), тот же самый адвективный поток, меняя свое направление, способствует ускорению пространственного распространения гена Bt-устойчивости (рис. 16б). В диффузионной фишеровской модели эффективность убежищ не столь высока (рис. 16д-з).

Шаблоны убежищ для поля: 16 км 16 км 646 3(a) (б) 3.3.0.18 0.(г) 1(в) квази-вымирание 12.0.46 квази-вымирание (е) 202 1(д) 2.84 1.0.44 0.- кол-во лет, в течение которого pr < 0.0.- особи/растение на (ж) 0.- особи/растение на Bt Рис. 17. Время (годы), за которое частота Bt-устойчивой аллели достигает 10% на всем поле T, и уровни заражения мотыльком (особи/растение) на всем поле и на Bt-поле, пока Bt (в течение ) в демо-генетической диффузионной модели (4.1). - участки убежищ, - pr < 0.1 TBt-участки. Размер убежища - 20% от общей площади поля для всех шаблонов. Коэффициент диффузии.

= 0.07 км2год-Влияние пространственной конфигурации убежищ на эффективность стратегии Увысокая доза-убежищеФ исследовано для фиксированных значений процента убежища и подвижности вредителя в случае прямоугольного ареала вредителя = [0, Lx][0, Ly] (рис. 17). Численные эксперименты с демо-генетической моделью показали, что для ареала мотылька 16 км на 16 км однополосный шаблон убежища (рис. 17а) оказывается наиболее эффективным из всех рассматриваемых шаблонов (рис. 17а-ж). Расположение единственной полосы убежища в центре поля (рис. 17б) приблизительно вполовину сокращает время T10, вычисленное для шаблона на рис.17а, способствуя значительному росту уровня заражения Bt-поля мотыльком. Разбивка одной полосы убежища на несколько полос также снижает эффективность убежища (рис. 17в-г). Аналогичный эффект имеет место и для убежищ квадратной формы (рис. 17д-ж). Очевидно, данный результат зависит от размера области и коэффициента , и при меньшей подвижности вредителя (или при рассмотрении проблемы в большем, региональном масштабе) пятнистые структуры убежищ могут оказаться эффективнее агрегированных убежищ.

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

pr t = bpr (Wr -W )+ pr + 2 ln pr ;

t = - + a P + ;

(bW ( )) (4.8) P t = ea P - P P + P;

P ps + pr = 1; pr n = n = P n = 0, x , в которой учитывается как генетическуя структура популяции вредителя, так и трофическая динамика. Аналогичным образом может быть детализировано рассмотрение состояния растительного ресурса (см. Тютюнов и др. 2007).

Глава 5. Агрегированные модели динамики промысловых популяций В обзорной части главы отмечается, что при долгосрочном прогнозе динамики промысловых популяций, как правило, нет необходимости явного моделирования мелкомасштабных (быстрых) процессов, связанных с пространственными перемещениями животных. Однако, нелинейности агрегированных моделей, неявно учитывающих пространственную неоднородность и взаимосвязь миграционного и жизненного циклов популяций, могут быть причиной сложных (хаотических) динамических режимов даже в точечных системах небольшой размерности (Ricker 1954; Ашихмина и др. 2004; Ильичев 2009). Ситуация дополнительно осложняется необходимостью учета в моделях промысла размерной и возрастной структур, а также открытостью популяционных систем внешним, в том числе - стохастическим воздействиям (Домбровский и др. 1991; Абакумов 2001). Характерными чертами возникающих при оптимизации долгосрочного промысла задач являются стохастичность и многокритериальность - модели должны обеспечить возможность поиска управляющих воздействий (как вылова, так и искусственного воспроизводства), максимизирующих урожай и, в то же время, стабилизирующих динамику и минимизирующих риск вымирания популяции (Тютюнов и др. 1996;

2000; Сенина и др. 1996). Для решения таких задач применим подход, предложенный Вильфредо Парето (Моисеев 1981; Vincent, Grantham 1981), позволяющий находить множество неулучшаемых точек, в рассматриваемом случае представляющих собой компромиссные стратегии эксплуатации популяции, которые являются материалом для дальнейшей экспертной оценки и принятия решения. Стохастичность моделей обуславливает необходимость использования метода МонтеЦКарло.

В главе представлен унифицированный подход к построению имитационных моделей, апробированный для трех объектов рыболовного промысла: (i) популяции окуня, (ii) пелагического сообщества азовских тюльки и хамсы, а также (iii) азовского судака. Подход основан на описании ключевых процессов (воспроизводства, смертности, роста, вылова) при помощи концептуальных моделей с известными свойствами (модель воспроизводства Риккера, модель возрастной структуры Лесли, модель роста Берталанфи и др.). Несмотря на относительную простоту, модели позволили сделать ряд содержательных рекомендаций относительно промысловой политики.

Модель популяции окуня (Perca fluviatilis L.) создана по заказу офиса охраны фауны кантона Ву Швейцарии и предназначена для исследования последствий изменения стратегии промысла, в частности, изменения размера ячеи жаберных сетей (Tyutyunov et al. 1993). Данный вид живет до 8 лет, характеризуется коротким циклом воспроизводства, высокой плодовитостью, непродолжительным периодом нереста и тенденцией к каннибализму, который является естественным компенсаторным механизмом (Thorpe 1977). В Швейцарии наиболее представительный объём данных получен для популяции оз. Констанс (Staub et al.

1987; Staub, Kramer 1991). Динамика уловов и запаса окуня сложна и с трудом поддается прогнозированию. Ряды наблюдений характеризуются высокими амплитудами: максимальные и минимальные значения могут отличаться на одиндва порядка.

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

С учетом пола и 7-и возрастных групп, вектор состояния системы включает в себя 14 переменных. Модель работает с полумесячным шагом t (24 шага в год).

Внутригодовая динамика представляет собой экспоненциальную функцию выживаемости с учетом естественной и промысловой смертности. Во время нереста (первая половина мая в оз. Констанс), каждая самка возраста производит bp количество икринок, рассчитываемое по формуле E( ) = ap(Lf ( )), где Lf ( ) - средняя длина тела самки, вычисляемая по модели Берталанфи. Воспроизводство описывается формулой Риккера. Так как прикладная задача, для которой создавалась модель, связана с выбором размера ячеи жаберных сетей, моделирование роста особей принципиально необходимо, как и учет половой структуры популяции окуней, самки которых крупнее самцов. Соответственно, промысловая смертность вычисляется с учетом селективности орудий лова относительно длины тела рыбы.

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

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

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

Так как рыбаки стремятся уменьшить размер ячеи, получая сиюминутный выигрыш, контроль размера ячейки жаберных сетей со стороны природоохранных организаций, прежде всего, подразумевает его ограничение снизу. Отметим, что служба рыбоохраны в Швейцарии обладает реальными возможностями контроля орудий лова; нелегальный промысел минимален. Известны и примеры официального временного разрешения использования мелких сетей в некоторых озерах. Такая практика преследует цель регулирования нежелательных колебаний запасов окуня путем интенсивного облова так называемых УсильныхФ многочисленных когорт (Pattay 1978). Построенная модель - один из инструментов прогноза последствий такого рода управления популяциями.

Модель динамики пелагических промысловых рыб Азовского моря: тюльки (она же азовская килька) (Clupeonella delicatula delicatula Nordm.) и хамсы (анчоуса) (Engraulis encrasicholus maeoticus Pus) несколько сложнее (Сенина и др. 1996; Senina et al. 1999). Ихтиологи рассматривают эти два вида как пищевых конкурентов, динамика популяций которых в значительной степени определяется факторами среды (Бронфман и др. 1979; Луц 1986). Это обусловливает необходимость создания неавтономной модели системы конкурирующих популяций в виде:

2 -S Qt -Q St -,0 -,0 S Q t t t t+1 t 2 = a e-b e- (K1 +K2 )e e ; (5.1) 1 t +1 t t = (p + q ) (1 - h t ); (5.2) 2 1 2 -SK Qt -QK St -,0 -,0 S Q t t t t+1 t K K K K 1 2 K1 = aK K2e-b K2e- ( + )e e ; (5.3) t t t t K2+1 = (pK K1 + qK K2) (1- hK ), (5.4) где, - численности младшей и старшей возрастных групп хамсы; K1, K2 - 1 тюльки; a, aK - плодовитости при условии оптимальности внешних факторов; b, bK и , K - коэффициенты внутривидовой и межвидовой конкуренции соответственно; p, pK и q, qK - коэффициенты выживаемости годовиков и взрослых; h, hK - доли промыслового изъятия; S и Q - среднегодовые значения S солености и речного стока; S0, Q0 - оптимальные значения внешних факторов; , Q - толерантности популяций к солености и стоку соответственно. Шаг модели - год.

Используя натурные данные по изучаемым популяциям за 1953-1981 гг., была проведена идентификация параметров детерминированной модели (5.1)-(5.4).

Аналитическое исследование устойчивости равновесий модели было проведено при постоянных значениях внешних факторов, зафиксированных на оптимальном уровне. Поскольку модель с полученными параметрами не демонстрирует хаотического поведения, то ее динамика при установлении существенно не зависит от начальных условий. Это свидетельствует в пользу высказываемой ихтиологами гипотезы о том, что наблюдаемые флуктуации численности возникают в основном благодаря влиянию внешних факторов. Речной сток и соленость моря выбраны как наиболее важные факторы, оказывающие критическое влияние на агрегированную динамику популяций, ввиду их прямого воздействия на состав кормовой базы пелагических видов (Бронфман др. 1979, Луц 1986). В целях неявного учета влияния других внешних факторов, таких как температура, ветровая активность, продуктивность планктона и др., был применен подход, ранее предложенный академиком РАН И. И. Воровичем с соавт. (1989) - разбиение множества рассматриваемых лет по значениям функции воспроизводства на урожайные и неурожайные годы. Применение такой усложненной модели со стохастическим воспроизводством обосновано при помощи критерия максимального правдоподобия. Для всех оцененных параметров модели рассчитаны их стандартные ошибки и проведен анализ чувствительности.

Имитационная модель, идентифицированная для популяций хамсы и тюльки, использовалась для расчета траекторий среднегодовых численностей и уловов. При этом, стохастические ряды внешних факторов генерировались по моделям Ратковича (1977), Бронфмана и Суркова (1976) и Суркова с соавт. (1985). Используя метод Монте-Карло, рассчитывались 50 репликаций, а полученные траектории использовались для вычисления осредненных характеристик, необходимых для численного решения стационарной задачи оценки риска и поиска оптимальных стратегий промысла:

max {C K = C + CK };

t t 0

w h T t t C = (p + q );

1 T t=wK hK T t t CK = (pK K1 + qK K2), T t =где w и wK - вес одной особи хамсы и тюльки (рыночные цены на тюльку и хамсу t t одинаковы). P = P{ } и PK = P{K2 Kкр} - вероятности квазивымирания 2 кр (падения численности хотя бы одной из популяций ниже зафиксированного критического уровня в течение рассматриваемого периода времени T ), рассчитанные для каждого вида (Ginzburg et al. 1982). Критические значения и кр Kкр, при которых наступает квазивымирание, были положены равными 5% от средней наблюденной численности взрослой возрастной группы.

Результаты численного решения оптимизационной задачи представлены в виде изоплетных диаграмм (рис. 19) на плоскости промысловых усилий (h, hK ) и поверхностей Парето для различных оптимизационных критериев (средний вылов, вариация вылова, средняя и минимальная численность, риск вымирания). Анализ результатов показал, что текущая (для рассматриваемого периода времени) стратегия (вылов 20% от взрослой части популяции хамсы и 30% от взрослой тюльки) могла быть улучшена без ущерба для популяций. Согласно имитационной модели, даже при удвоении промысловых усилий, популяции остаются на относительно низком уровне риска. С выбранными порогами квазивымирания вероятность падения численности ниже или Kкр в пределах рассматриваемого кр периода времени T = 35 лет не превышает 0.07. Однако, дальнейшее увеличение промысловых усилий приведет к сильному перелову и, как следствие, резкому увеличению риска квазивымирания, который может быть даже больше, если критические численности (, Kкр ) зафиксировать на более высоком уровне.

кр Рис. 19. Изоплетная диаграмма суммарного годового вылова (в тыс. тонн) в зависимости от коэффициентов промыслового изъятия.

Рассматривая сообщество двух популяций, немаловажно рассмотреть вопрос:

как конкуренция за пищевой ресурс влияет на выводы относительно оптимальной промысловой стратегии (Clark 1976). Как видно на рис. 19, модель показывает, что перелов эндемичной тюльки приводит к сильному падению общего вылова, тогда как практически полное изъятие мигрирующей хамсы лишь увеличивает критерий C K за счет роста численности тюльки в отсутствие конкурента.

Модель азовского судака (Lucioperca lucioperca L.) построена по аналогичной схеме, с учетом структурирования популяции по трем возрастным группам:

2 t -S0 T -T0 St - - t t T t+1 t t S 1 = (a1 + a2 )e-b( + )e e ; (5.5) 0 1 t+1 t t = p0 (1 - h1 ); (5.6) 1 t +1 t t t = ( p1 + p2 )(1- h2 ), (5.7) 2 1 где a1, a2 - плодовитости годовиков (1+) и взрослых особей (2+) судака; b - коэффициент внутривидовой конкуренции, p0, p1, p2 - коэффициенты выживаемости сеголеток (0+), годовиков и взрослых, h1, h2 - доли промыслового изъятия, S и T - среднегодовые значения солености и температуры, S0, T0 - оптимальные значения внешних факторов, , - толерантности популяции к солености и температуре S T воды соответственно.

Идентификация модельных параметров, показавшая хорошее приближение натурных данных моделью (для переменных, и вычисленные значения F- 0 1 статистики высоко значимы (p < 0.005) и равны, соответственно, F = 4.22, F = 4.01 и F = 4.40), проводилась в два этапа. На первом шаге были получены оценки значений параметров a1, a2, b, , при помощи метода наименьших квадратов, путем S T приближения траектории нелинейного уравнения воспроизводства (5.13) к натурным значениям численности сеголеток, при использовании наблюденных рядов,, S и T в качестве независимых внешних данных. Полученные таким 1 образом значения параметров затем были использованы в качестве начальных значений на втором шаге идентификации параметров. На этом этапе значения параметров a1, a2, b, ,, p0, p1 калибровались, путем приближения к данным S T t t t наблюдений траекторий,, полной модели (5.5) - (5.7). Начальные 0 1 численности возрастных групп соответствовали данным 1949 года. Для определения ошибок оценки параметров использовалась техника бутстрапинга (Efron, Tibshirani 1993).

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

кружочек соответствует оптимальной (в случае единственного экономического критерия) стратегии промысла. Звездочкой отмечен пример Парето-оптимальной стратегии ( h1 = 0.0; h2 = 0.07 ).

Как и модель пелагических рыб (5.1)-(5.4), имитационная модель (5.5)-(5.7) со стохастически моделируемыми внешними факторами (среднегодовой соленостью и температурой моря), использовалась для генерации осредненных по 50-ти репликациям 35-летних траекторий среднегодовых численностей возрастных групп и последующего решения задачи двухкритериальной оптимизации. В качестве стратегии промысла рассматривались комбинации промысловых усилий на годовиков и взрослых особей:

t t t t t max {Ct = w1h1 p0 + w2h2(p1 + p2 )} 0 1 t t 0

Результаты решения оптимизационной задачи подтверждают критическое влияние солености Азовского моря на продуктивность популяции судака. Кроме того, продемонстрирована настоятельная необходимость ограничения вылова в целях снижения риска квазивымирания, который даже при гипотетическом прекращении промысла, остается высоким, близким к 0.7 (при критической численности, зафиксированной на уровне 10% от наблюденной численности взрослой части популяции). Показано также, что непреднамеренный вылов годовиков нежелателен при любой интенсивности промысла, однако не так опасен для популяции, как перелов взрослых особей. Результаты моделирования обосновывают необходимость контроля промысла судака, развития искусственного воспроизводства, рационализации сезонного распределения стока Дона и Кубани, и совершенствования селективности орудий лова для предотвращения попадания в сети молоди судака.

ОСНОВНЫЕ ЗАЩИЩАЕМЫЕ ПОЛОЖЕНИЯ

И ВЫВОДЫ ДИССЕРТАЦИИ:

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

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

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

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

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

4. Продемонстрирована эффективность учёта адвективного потока генов, связанного с неоднородностью ареала обитания популяции. А именно, построена и исследована демоЦгенетическая модель пространственноЦвременной динамики диплоидной популяции с учётом адвективного потока генов. На основании этой модели предложено объяснение значительной (десятки и сотни лет) задержки развития устойчивости кукурузного мотылька Ostrinia nubilalis Hubner к токсину, продуцируемому трансгенной Bt кукурузой. В результате укрепляется теоретический базис стратегии Увысокая доза - убежищеФ, широко используемой в сельскохозяйственном производстве с применением генетически модифицированных инсектицидных культур;

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

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

Список основных работ, в которых изложены результаты диссертации:

1. Тютюнов Ю.В., Титова Л.И., Сурков Ф.А., Бакаева Е.Н. Трофическая функция коловратокЦ фитофагов (Rotatoria). Эксперимент и моделирование // Журнал общей биологии, 2010, 71, № 1, С. 52-62.

2. Тютюнов Ю.В., Загребнева А.Д., Сурков Ф.А., Азовский А.И. Микромасштабная пятнистость распределения веслоногих рачков как результат трофически-обусловленных миграций // Биофизика, 2009, 54, № 3, С. 508-514.

3. Тютюнов Ю.В., Загребнева А.Д., Сурков Ф.А., Азовский А.И. Моделирование потока популяционной плотности организмов с периодическими миграциями // Океанология, 2010, 50, № 1, С. 1-10.

4. Tyutyunov Yu., Zhadanovskaya E., Bourguet D., Arditi R. Landscape refuges delay resistance of the European corn borer to Bt-maize: a demo-genetic dynamic model // Theoretical Population Biology, 2008a, 74: 138-146.

5. Tyutyunov Yu., Titova L., Arditi R. Predator interference emerging from trophotaxis // Ecological Complexity, 2008, 5: 48-58.

6. Medvinsky A.B., Gonik M.M., Tyutyunov Y.V., Li, B.-L., Rusakov A.V., Malchow H. Insecticidal Bt crops under massive Bt-resistant pest invasion: Mathematical simulation. In: Aspects of Mathematical Modelling (ed. by R. Hosking and E. Venturino). Birkhauser, Basel. 2008. pp. 81-94.

7. Tyutyunov Yu., Titova L., Arditi R. A minimal model of pursuitЦevasion in a predatorЦprey system // Mathematical Modelling of Natural Phenomena, 2007, 2(4): 122-134.

8. Тютюнов Ю.В., Жадановская E.A., Ардити Р., Медвинский A.Б. Пространственная модель развития устойчивости насекомых-вредителей к трансгенной инсектицидной сельскохозяйственной культуре на примере кукурузного стеблевого мотылька // Биофизика, 2007, Т. 52, Выпуск 1, С. 95-113.

9. Жадановская E.A., Тютюнов Ю.В., Ардити Р. Моделирование стратегии Увысокая доза - убежищеФ при использовании генетически модифицированной кукурузы для подавления кукурузного стеблевого мотылька // Известия ВУЗов, Северо-Кавказский регион.

Естественные науки, 2006. - Приложение, № 11. - С. 5-9.

10. Тютюнов Ю.В., Сапухина Н.Ю., Сенина И.Н., Ардити Р. Таксис как фактор, стабилизирующий трофическую систему // Обозрение прикладной и промышленной математики, М.: Научное изд-во ТВП, 2005, Т. 12, Выпуск 4, С. 810-814.

11. Tyutyunov Yu., Senina I., Arditi R. Clustering due to acceleration in the response to population gradient: a simple self-organization model // The American Naturalist, 2004, 164(6): 722-735.

12. Arditi R., Callois J.-M., Tyutyunov Yu., Jost. C. Does mutual interference always stabilize predatorЦ prey dynamics? A comparison of models // Comptes Rendus Biologies, 2004, 327: 1037-1057.

13. Сапухина Н.Ю., Тютюнов Ю.В. Явная модель пространственной динамики трофического сообщества растительная культураЦвредительЦхищник // Компьютерное Моделирование.

Экология: Выпуск 2 (под. ред. Г.А. Угольницкого). М.: Вузовская книга, 2004. C. 66-91.

14. Sapoukhina, N., Tyutyunov, Yu., Arditi, R. The role of preyЦtaxis in biological control: a spatial theoretical model // The American Naturalist, 2003, 162(1): 61-76.

15. Tyutyunov Yu., Senina I., Jost C., Arditi R. Risk assessment of the harvested pike-perch population of the Azov Sea // Ecological Modelling, 2002, 149: 297-311.

16. Сенина И.Н., Тютюнов Ю.В. Моделирование стаеобразования как следствия автотаксиса // Журнал общей биологии, 2002, Т. 63, № 6. С. 483-488.

17. Тютюнов Ю.В., Сапухина Н.Ю., Сенина И.Н., Ардити Р. Явная модель поискового поведения хищника // Журнал общей биологии, 2002, Т. 63, № 2. С. 137-148.

18. Сенина И.Н., Тютюнов Ю.В., Кузнецов А.В. Моделирование миграционного цикла пелагических рыб Азовского моря. // Экосистемные исследования Азовского моря и побережья (под. ред. акад. РАН Г.Г. Матишова). Апатиты: изд. Кольского научного центра РАН, 2002, Т. 4. C. 327-333.

19. Сапухина Н.Ю., Тютюнов Ю.В., Сенина И.Н., Ардити Р. Влияние кормовых миграций рыб на возникновение пятнистости распределения морских сообществ // Экосистемные исследования Азовского моря и побережья (под. ред. акад. РАН Г.Г. Матишова). Апатиты: изд. Кольского научного центра РАН, 2002, Т. 4. C. 310-326.

20. Тютюнов Ю.В., Сенина И.Н., Сурков Ф.А., Ардити Р. Модели оценки риска сокращения численности промысловых популяций рыб // Среда, Биота и Моделирование Экологических Процессов в Азовском море (под. ред. акад. РАН Г.Г. Матишова). Апатиты: изд. Кольского научного центра РАН, 2001a, 413. C. 380-396.

21. Arditi R., Tyutyunov Yu., Morgulis A., Govorukhin V., Senina I. Directed movement of predators and the emergence of density-dependence in predatorЦprey models // Theoretical Population Biology, 2001, 59(3): 207-221.

22. Тютюнов Ю.В. Сапухина Н.Ю., Моргулис А.Б., Говорухин В.Н. Математическая модель активных миграций как стратегии питания в трофических сообществах // Журнал общей биологии, 2001, Т. 62, № 3. С. 253-262.

23. Сапухина Н.Ю., Тютюнов Ю.В. Влияние миграций насекомых на регуляторную роль энтомофага при биологическом контроле популяции вредителя // Компьютерное Моделирование. Экология (под. ред. Г.А. Угольницкого). М.: Вузовская книга, 2000. С. 36-57.

24. Тютюнов Ю.В., Сенина И.Н., Титова Л.И. Экономический и экологический критерии оптимизации промысла Азовского судака. Опыт имитационного моделирования // Компьютерное Моделирование. Экология (под. ред. Г.А. Угольницкого). М.: Вузовская книга, 2000. С. 58-78.

25. Говорухин В.Н., Моргулис А.Б., Тютюнов Ю.В. Медленный таксис в модели хищникЦжертва // Доклады академии наук, 2000, T. 372, № 6. C. 730-732.

26. Senina I., Tyutyunov Yu., Arditi R. Extinction risk assessment and optimal harvesting of anchovy and sprat in the Azov Se // Journal of Applied Ecology, 1999, 36(2): 297-306.

27. Говорухин В.Н., Моргулис А.Б., Сенина И.Н., Тютюнов Ю.В. Моделирование активных миграций пространственно распределенной популяции // Обозрение прикладной и промышленной математики, М.: Научное изд-во ТВП, 1999, Т. 6, Выпуск 2. C. 271-295.

28. Тютюнов Ю.В., Сенина И.Н. Трофическая функция как результат пространственного поведения // Проблемы проектирования и управления экономическими системами:

инвестиционный аспект. Ростов-на-Дону: РГЭА, 1998, С. 132-135.

29. Сенина И.Н., Домбровский Ю.А., Тютюнов Ю.В., Обущенко Н.И., Титова Л.И. Имитационное моделирование и оптимизация промысла сообщества двух конкурирующих популяций // Деп.

в ВИНИТИ, 1996, №2563-В96, 57 с.

30. Тютюнов Ю.В., Домбровский Ю.А., Обущенко Н.И. Оптимальное управление эксплуатируемой популяцией при минимизации риска ее вымирания в условиях стохастичности среды обитания // Обозрение прикладной и промышленной математики, М.:

Научное изд-во ТВП, 1996, Т. 3, Выпуск 3, С. 412-433.

31. Tyutyunov Yu., Dombrovsky Yu., Arditi R., Surkov F. The influence of dispersal behaviour on metapopulation viability // Journal of Biological Systems, 1996, 4(2): 277-290.

32. Tyutyunov Yu., Arditi R. Application du modele REEL a des donnees recentes. In: Rapport de la Commission internationale pour la protection des eaux du Leman contre la pollution. Campagne 1993, (1994), pp. 243-255.

33. Tyutyunov Yu., Arditi R., Buttiker B., Dombrovsky Yu, Staub E. Modelling fluctuations and optimal harvesting in perch populations // Ecological Modelling, 1993, 69: 19-42.

34. Домбровский Ю.А., Обущенко Н.И., Тютюнов Ю.В. Рыбные Популяции в Стохастической Среде: Модели Управления и Выживаемости. Ростов-на-Дону: Издательство Ростовского госуниверситета, 1991, 160 c.

Авторефераты по всем темам  >>  Авторефераты по биологии