Построение, исследование и приложения математических моделей пространственно-временной динамики популяционных систем (
На правах рукописи
ТЮТЮНОВ Юрий Викторович ПОСТРОЕНИЕ, ИССЛЕДОВАНИЕ И ПРИЛОЖЕНИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ПРОСТРАНСТВЕННО-ВРЕМЕННОЙ ДИНАМИКИ ПОПУЛЯЦИОННЫХ СИСТЕМ (03.00.02 – Биофизика)
АВТОРЕФЕРАТ
диссертации на соискание ученой степени доктора физико-математических наук
Красноярск – 2009 2
Работа выполнена в отделе математических методов в экономике и экологии Научно-исследовательского института механики и прикладной математики им. Воровича И.И.
Южного федерального университета
Научный консультант: доктор физико-математических наук, профессор МЕДВИНСКИЙ Александр Берельевич
Официальные оппоненты: доктор физико-математических наук, профессор ХЛЕБОПРОС Рем Григорьевич доктор физико-математических наук, профессор СУХИНОВ Александр Иванович доктор физико-математических наук, профессор САДОВСКИЙ Михаил Георгиевич
Ведущая организация: МГУ им М.В.Ломоносова (кафедра биофизики биологического факультета)
Защита диссертации состоится “ 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 = (S ) µ (S ), (1.1) J x x где = (x, t ) – плотность гарпактицид в точке x в момент времени t, S = S (x, t ) – концентрация некоторого стимула, влияющего на интенсивность вертикальных 2 миграций индивидуумов, µ (S ) = l P (S ) – коэффициент диффузии, (S ) = l dP ( S ) – 2 2 dS коэффициент таксиса, P(S ) – вероятность того, что в течение периода времени особь выйдет из грунта в воду при значении стимула S, l 2 – среднее значение квадрата длины горизонтального перемещения особи (0, ). Если выход рака из грунта в воду является простейшим пуассоновским процессом с переменным параметром (x, t ) = f (S (x, t )), где f (S ) – зависимость частоты выхода рака в воду от стимула S, то коэффициенты диффузии и таксиса есть:
l 2 df (S ) d µ (S ) l µ (S ) = f ( S ), (S ) =. (1.2) = 2 2 dS dS Полученный результат подтверждает универсальность модели Пэтлака–Келлер– Сегеля и её применимость к описанию феномена таксиса в различных ситуациях.
Индивидуум–ориентированная модель, построенная на основе гипотез, использованных при выводе уравнения потока (1.1) иллюстрирует, что рассмотренный механизм положительного таксиса, т.е. убывание частоты выхода рака в воду при возрастании концентрации стимула, приводит к агрегированию организмов в местах с повышенной концентрацией стимула. Как видно на рис. 1, при высоких численностях популяции динамика системы хорошо аппроксимируется эйлеровой непрерывной моделью таксис–диффузия (Тютюнов и др. 2009;
2010):
(1.3) = J.
t x Рис. 1. 1D распределение популяции гарпактицид (справа) в средах со стационарным распределением стимула (слева) при t = 325. Ось абсцисс – горизонтальная координата, ось ординат – концентрация стимула (слева) и обилие особей (справа), усл. единицы.
1 – получено с помощью непрерывной модели, 2 – получено с помощью индивидуум– ориентированной модели. Параметры: частота выхода в воду f (S ) = exp ( S / 2 ), период времени = 1, количество шагов T = 325, количество особей M = 2000, среднее значение квадрата длины перемещения l 2 = 2 = 1, распределение стимула в случае непрерывной модели ( ) (x ) = 8 exp (x 20 )2 / 16. В начальный момент все особи находились в начале координат.
S Модель гарпактициды–микроводоросли строится как система хищник–жертва R 2R R = rR1 aR + R 2, t x K (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 ) 0 0, что хорошо согласуется с законом Вебера–Фехнера.
S Но какова должна быть природа стимула для того, чтобы в данной модели реализовались неоднородные динамические режимы? Показано, что общепринятая гипотеза о том, что стимулом таксиса хищника является неоднородность распределения жертв, не позволяет получить неоднородные режимы в (1.4), так как однородное распределение устойчиво при любых допустимых значениях параметров модели, для которых оно положительно.
Если стимулом таксиса является некоторое выделяемое жертвой вещество S ( x, t ), являющееся аттрактантом для хищников, то модель (1.4) принимает вид:
R 2R R = rR1 aR + R 2, t x K S = (S ) + µ (S ) (1.5), t x x x x S 2S = R S + S 2, 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 -ой моды неоднородного пространственного возмущения (r ( x, t ), n( x, t ), s( x, t )) = rk e k t cos kx, nk e k t cos kx, s k e k t cos kx :
() ( () ) k S * k,, k,, µ S *,k 0, (1.7) R S (,, k,, µ (S ), ) = ( R ) [(k ) ( + )( µ (S ) + µ (S ) + + µ (S )) * 1 * * * 2 * R S R S R S R S + (k ) ((R + )(µ (S ) + 2 µ (S )( + ) + ) + ( + )( R + )) 2 * 2 * * * R S R S R S S R + ((R + ) µ (S ) + (R + )( R + ) + ( + )R )]+ (k ) (R + ).
2 1 * * * * * * S R R S Превышение коэффициентом таксиса (S ) критического значения влечет * бифуркацию Пуанкаре–Андронова–Хопфа: в окрестности колебательно теряющего устойчивость стационарного однородного режима (R *, *, S * ) рождается периодический пространственно–неоднородный волновой режим. Возникающий режим может быть довольно сложным, если первой теряет устойчивость мелкомасштабная мода (например, k = 2, 3 или 4 на рис.2 при достаточно высоких значениях коэффициента распада ).
R = 0.01, S = 0.002, µ(S* ) = 0.001, (S* ) = 2 4 I 5 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 2R R = rR1 aR + R 2, t x K S = (S ) + µ (S ) (1.8), t x x x x S = eaR S, t где e – коэффициент “усвоения” пищи хищником, eaR – количество пищи, которое в единицу времени в среднем попадает в желудок хищника, – коэффициент переваривания, S – количество перевариваемой в единицу времени пищи. В отличие от (1.5) в данной модели имеет смысл положить S = 0, поскольку количество пищи в желудке хищника в точке x не влияет на сытость хищников в соседних точках.
Так как модель (1.8), (1.6) эквивалентна модели (1.5), (1.6) при = ea и S = 0, то для нее справедливы результаты линейного анализа, представленные выше. Обе модели учитывают инерционность (запаздывание) реакции хищника на изменения плотности жертв: хищники реагируют на изменение концентрации стимула (сытости), который изменяется в зависимости от плотности жертв. Именно благодаря этой задержке реакции хищника на градиент распределения жертв, при выполнении условия (1.7) возникают пространственно неоднородные режимы.
Поскольку стимул является потенциалом поля скорости таксиса: v = S, уравнение для S можно представить в виде v t = R v + S v. Таким образом, фактически, главное отличие предложенного метода моделирования таксиса хищника заключается в гипотезе о зависимости ускорения таксиса, а не скорости, от градиента плотности жертв. При малых надкритичностях модель (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 v P t = P + v v P, P (1.10) P n = n = v P =v = 0.
Здесь и P – плотности жертв и хищников;
v и v P – скорости их таксиса;
, P, v, v – диффузионные коэффициенты;
, P – коэффициенты таксиса. Оператор P Лапласа для скалярных переменных, P имеет вид = div grad, а для векторного поля v : v = grad divv rot rotv. Отсутствие потоков популяционных плотностей на границе области означает постоянство средних популяционных плотностей и P, которые могут быть масштабированием переменных, P := P P := приведены к единице. Линейный анализ однородного стационарного режима ( * = 1;
P * = 1;
v * = 0;
v *P = 0) дает условие возбуждение k -ой моды разложения в ряд Фурье малого неоднородного по пространству возмущения:
k 4 ( v + v P )( P + v P )( P + v )( + v P )( + v )( + P ) P. (1.11) ( v + v P + + P ) С ростом произведения P моды с номерами k = 1, 2,… возбуждаются последовательно. Включение в уравнения скоростей сил трения в виде v и v P может изменить порядок возбуждения мод, что делает модель более реалистичной.
Заметим, что уравнения скоростей могут быть заменены на эквивалентные уравнения баланса стимулов (выделяемых жертвами и хищниками веществ феромонов, аттрактанта для хищников и репеллента для жертв). При этом коэффициенты трения “превращаются” в коэффициенты распада феромонов.
Бифуркация возникновения и дальнейшего развития волновых режимов убегания– преследования наблюдалась в численных экспериментах с сеточной аппроксимацией модели (1.9)-(1.10).
Модель стаеобразования. Предложенный метод позволяет решить проблему, сформулированную в работе (Edelstein-Keshet et al. 1998) и заключающуюся в том, что стандартные модели таксиса не позволяют воспроизвести формирование и перемещение устойчивого скопления животных, индуцированного исключительно неоднородностью распределения популяции и обладающего реалистичным профилем – постоянной плотностью внутри группы и четким границами снаружи.
Наиболее известной моделью агрегирования организмов является модель Келлер– Сегеля (Keller, Segel 1971). Постулируя пропорциональность скорости таксиса градиенту выделяемого животными аттрактанта, v = S, она позволяет описать направленное перемещение организмов, но в случае, если размерность пространства превышает единицу, прогнозирует взрыв популяционной плотности [возможность взрыва в одномерном случае была исключена сравнительно недавно (Osaki, Yagi 2001;
Hillen, Potapov 2004)]. Могильнером и Эдельштейн-Кешет (Mogilner, Edelstein Keshet 1999) предложена более реалистичная модель, представляющая собой систему интегро-дифференциальных уравнений, в которой взаимное притяжение и отталкивание организмов в каждой точке определяется нелокальной информацией о плотности популяции.
1500 150 1500 РЯД 1, x РЯД 1, y 1000 120 R = 0.79 R = 0. dvdt, см /с - - 30 - - Плотность популяции, P -1500 -1500 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 Среднее ускорение 800 РЯД 2, y РЯД 2, x R = 0. R = 0.71 90 - - -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, следует использовать не постоянный, а зависящий (в простейшем случае линейно) от плотности коэффициент (1 P Pc ). Получаем уравнение:
dv dt = (1 P Pc )P, (1.12) хорошо аппроксимирующее данные наблюдений Окубо и имеющее лишь два параметра.
С учетом замкнутости области и сопротивления среды полная модель имеет вид:
P t + div(Pv ) = P P, (1.13) v t = (1 P Pc )P v + v v.
(1.14) n P = n v = 0, x x Показано, что если коэффициент аутотаксиса превышает критическое значение c,mn = P ( v k mn + ) Pc, (1.15) P (Pc P ) то происходит монотонная потеря устойчивости однородного режима (P = P, v * = 0) относительно моды малого неоднородного возмущения с волновым * числом k mn. Данный результат получен для областей с различной геометрией:
отрезка длины L, прямоугольной области L1 L2, и диска радиуса L. Во всех случаях первой возбуждается мода, имеющая максимальный пространственный масштаб.
А) t = 0 Б) t = 2. Pmin = 0.092 Pmin = 5 10 Pmax = 25.91 Pmax = 19. В) t = 5.798 Г) t Pmin = 1.9 10 7 Pmin = 1 10 Pmax = 19.447 Pmax = 19. Рис. 4. Численный эксперимент моделирующий процесс формирования и дальнейшего перемещения скопления организмов в случае сингулярного возмущения однородного распределения. Параметры: = 3.7, P = 0.5, v = 0.05, = 0, Pc = 10, P = 1, L1 = L2 = 1.5.
Как показывает рис. 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), с постоянными коэффициенты диффузии и таксиса хищника: P = µ (S * ), = (S * ) :
t = r (1 K ) a P +, (2.1) (2.2) P t + div( Pv ) = P P, v t = + v, (2.3) v n v = n = n P = 0, x. (2.4) Для одномерного (отрезок [0, 1] ) и двумерного (прямоугольник [0, L] [0, 1] ) местообитаний исследованы как минимальная модель (2.1)-(2.4), так и ее обобщения, в том числе, случай учета трения в уравнении скорости таксиса:
v t = v + v v, (2.5) а также случай неконсервативной модели, с описанием процессов воспроизводства и смертности в уравнении баланса плотности популяции хищника:
P t + div( Pv ) = ea P µP + P P, (2.6) где e – коэффициент эффективности хищничества, µ – естественная смертность.
Линейный анализ устойчивости однородных равновесий рассматриваемых моделей к малым пространственным возмущениям позволил сделать следующие общие выводы:
Утверждение 2.1. Если = 0, т.е. процесс поиска жертв является случайным, то ненулевой однородный стационарный режим модели трофотаксиса ( *, P *, v * 0) устойчив к малым пространственным возмущениям при любых допустимых значениях параметров (в том числе при любых значениях диффузионных коэффициентов).
Утверждение 2.2. Если 0, т.е. хищник способен к направленному поиску и преследованию жертв, то для любого набора параметров существует критическое значение коэффициента трофотаксиса * 0, превышение которого приводит к колебательной потере устойчивости однородного стационарного режима ( ) ***, P, v и возникновению пространственно-неоднородных волновых решений.
Коэффициент трофотаксиса (a) (б) k k01 k k k k k20 k k10 k Коэффициент эффективности поиска, a Рис. 5. Критические кривые устойчивости первых пяти мод однородного стационарного режима ( ) ***, P, v, полученные аналитически и рассчитанные для L = 5, модели трофотаксиса µ = 0.01, r = K = e = 1, = 0, v = 0.001. а) = 0.01, P = 0.6 ;
б) = 0.005, P = 0.001.
На рис. 5 приведены критические кривые устойчивости однородного равновесия ( ), P *, v * на плоскости параметров (a, ) для двумерной неконсервативной модели с * трением (2.1), (2.6), (2.5), (2.4). Точки, лежащие ниже критической кривой моды возмущения с волновым вектором k nm = (n L, m ), соответствуют устойчивости однородного режима, выше – неустойчивости. В случае рис. 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.
Очевидно, возникающая при этом пространственная неоднородность оказывает влияние на агрегированную динамику исследуемой системы. Численный анализ агрегированной по времени и пространству смертности жертв обусловленной хищничеством:
1 g (, P) P def g (, P) P T dt =, (2.7) Q = lim T T показал, что в системе возникает эффект интерференции хищников, проявляющейся на популяционном уровне как зависимость среднего рациона хищника от их численности.
Рис. 6. Зависимость агрегированной смертности жертв обусловленной хищничеством от размера популяции хищников, рассчитанная при L =, v = 0.00001, = 0.03, P = 0.02 в сравнении с теоретическими кривыми точечных моделей хищник–жертва с трофическими функциями Лотки– Вольтерра (пунктир) и Ардити–Гинзбурга (сплошная линия).
Обнаружено (рис. 6), что зависимость Q(P ), где P – агрегированная по пространству плотность хищника, является линейной лишь в однородном случае, т.е. при * (P ). При установлении и продолжении по параметру P пространственно–неоднородного режима модели (2.1)-(2.4), зависимость Q(P ) деформируется, принимая вид нелинейной функции, насыщающейся с увеличением P. Подобным образом насыщается, в частности, зависимость Q (P ), рассчитанная для RD (ratio–dependent) трофической функции Ардити–Гинзбурга aP g( P) = (Arditi, Ginzburg 1989), представляющей один из возможных 1 + ah P способов учета зависимости рациона хищников от их численности в точечных моделях (см. Главу 3).
Парадокс биологического контроля. Благодаря обусловленному пятнистостью популяций эффекту интерференции, предложенный метод явного моделирования активного пространственного поведения хищника позволяет разрешить известный “парадокс биологического контроля”, заключающийся в невозможности при помощи классических методов моделирования трофических систем воспроизвести наблюдаемую в природе стабильную динамику системы хищник–вредитель при низкой численности популяции вредителя (Luck 1990;
Arditi, Berryman 1991).
Количественный критерий адекватности модели хищник–жертва предложен в работе (Beddington et al. 1978): модель должна воспроизводить стабильную динамику численности вредителя на уровне, не превышающем 2.5% от численности, достигаемой в отсутствие хищника. В рамках точечных моделей, описывающих агрегированную динамику систем хищник–жертва, разрешение парадокса биоконтроля возможно путем модификации классических, зависящих лишь от плотности жертв, трофических функций g ( ) таким образом, чтобы учесть в них эффект оказываемый интерференцией хищников на их рацион g (, P ). Это позволяет неявно учесть в точечных моделях пространственные эффекты, стабилизирующие природные системы. В частности, проблема решается при использовании упомянутой выше трофической функции Ардити–Гинзбурга.
Однако, являясь феноменологическими, точечные модели с интерференцией не объясняют механизма стабилизации системы хищник–жертва. Данное обстоятельство определяет актуальность построения пространственных моделей, способных явно описать пространственное поведение животных.
Для системы хищник–вредитель построена модель трофотаксиса, локальная кинетика в которой описывается классической системой Розенцвейга–МакАртура с логистическим воспроизводством жертв и трофической функцией Холлинга II типа:
aP t = (1 ) 1 + ah + ;
P aP µP div(Pv ) + P P;
(2.8) = t 1 + ah v t = + v v;
n v = n = n P = 0, x. (2.9) Изучена обусловленная трофотаксисом возможность образования пространственных структур из устойчивого в отсутствие пространственных эффектов однородного по пространству стационарного режима:
a (1 µh ) µ µ µ 1 + µh ( ), P*, v* =, 0, где *, (2.10) a, a (1 µh ) a (1 µh ) h(1 µh ) 1 µh а также механизм возникновения пространственной неоднородности из устойчивого в отсутствие пространственных эффектов однородного периодического режима C0, являющемся аттрактором модели Розенцвейга–МакАртура при a (1 + µh ) (h µh 2 ).
При больших a режим C 0 характеризуется высокой амплитудой – численность жертв периодически достигает значений близких емкости среды. Именно наличие такого режима в классических моделях послужило причиной возникновения парадокса биоконтроля.
Сформулируем условие устойчивости равновесия для случая 1D модели:
Однородный стационарный режим Теорема 2.1 (условие устойчивости).
(, P, v ) модели (2.8)-(2.9) устойчив к малым пространственно неоднородным * * * возмущениям вида (n( x, t ), p( x, t ), v( x, t ) ) = nk e t cos kx, p k e t cos kx, v k e t sin kx тогда k k k и только тогда, когда для каждого волнового числа k = n L, n 0, выполняется неравенство µa 1 a ( ) + P 11 ( v + P ) ( + v )k a 11 + 221, 0 (2.11) P*µ k 2 k где a11 0, a21 0 – элементы Якобиана системы (2.8).
Следствием данной теоремы являются утверждения, аналогичные приведенным выше для несколько более простой модели трофотаксиса (2.1)-(2.6), о возникновении адвективной неустойчивости режима ( *, P *, v * ) при достаточно высоком значении.
Рис. 7. Неустойчивый однородный периодический режим C 0 и устойчивый неоднородный периодический режим R, индуцируемый = 0.5, представлены на плоскости осредненных по пространству численностей популяций. Параметры модели: a = 0.7, h = 2.857, L = 5, µ = 0.01, = 0.01, P = 0.6, v = 0.001.
Влияние трофотаксиса на устойчивость однородного периодического режима C 0, существующего при a (1 + µh ) (h µh 2 ), исследовалось численно, путем определения положения мультипликаторов оператора монодромии на комплексной плоскости при последовательном увеличении коэффициента таксиса. Таким образом, были определены бифуркационные значения *, приводящие к адвективной потере устойчивости однородного режима C 0 и возникновению неоднородных по пространству колебаний R.
Пространственное осреднение Локальные колебания a б в Рис. 8. Фазовые траектории и локальная динамика ( x = 20 ) режимов модели (2.8)-(2.9):
(а) = 0.5 – неоднородный периодический режим;
(б) = 1 – неоднородный квазипериодический режим;
(в) = 1.5 – неоднородный хаотический режим со вспышками плотности жертв. При = 0 устойчив пространственно-однородный периодический режим C 0 : [0.00001;
0.62], P [0.51;
3.06]. Параметры приведены в подписи рис. 7.
Как видно на рис. 7 и рис. 8, амплитуда колебаний численностей популяций режима R значительно меньше амплитуды исходного пространственно однородного периодического режима C 0. Устойчивая по Лагранжу динамика наблюдается для широкого диапазона значений, пока слишком высокая активность хищника не дестабилизирует систему.
Для случая прямоугольного местообитания результаты аналогичны. Рис. показывает примеры периодического и хаотического режимов двумерной модели.
(а) (б) Рис. 9. Распределение плотности популяции жертвы на двумерном местообитании для (а) периодического режима, реализующегося при возбуждении моды с волновым вектором s ( = 0.00069273 ), (б) пространственно–временного хаоса, возникающего при продолжении режима по параметру в надкритическую область ( = 0.0009 ). Другие параметры модели:
L x = 4, L y = 3, a = 0.06, h = 10, µ = 0.01, = 0.001, P = 0.001, v = 0.001.
Интересно, что и более простая модель (2.1), (2.6), (2.5), (2.4), отличающаяся от системы (2.8)-(2.9) трофической функцией, описывающей локальные взаимодействия жертв с хищниками (линейная зависимость Лотки–Вольтерра вместо насыщающейся функции Холлинга II типа, а), позволяет получить сложные динамические режимы, качественно совпадающие с представленными на рис. (Тютюнов и др. 2002). Это объясняется тем, что реальный контакт хищников и жертв имеет место на границах пятен их популяций, где плотность жертв не высока и насыщения функции Холлинга в большинстве случаев не происходит.
Таким образом, явное моделирование пространственной активности хищников, согласно предложенному методу, позволяет значительно упростить описание локальной динамики трофических сообществ. При этом, на макро-уровне, в системе возникает интерференция хищников, а агрегированная динамика модели удовлетворяет критериям адекватности Беддингтона, разрешая парадокс биоконтроля. В качестве примера работоспособности метода, в работе проводится исследование модели трехуровневой трофической системы растительный ресурс – вредитель – хищник в котором и вредитель и его естественный враг (хищник, паразитод) способны к трофотаксису. Результаты главы интерпретированы в терминах прикладных задач биологических методов защиты растений и обоснования свойств видов насекомых, используемых в качестве агентов биоконтроля (в частности, божьих коровок как естественных врагов тли).
Глава 3. Интерференция хищников как проявление пространственных эффектов Обзорная часть главы посвящена фундаментальному понятию теории моделирования трофических систем – трофической функции, количественно характеризующей зависимость среднего рациона хищника от плотности популяций хищников и их жертв (см. табл. 1).
Таблица 1. Примеры трофических функций, с учетом и без учета интерференции хищников.
Название Выражение Источник (авторы) g( ) = a Лотки–Вольтерра (LV) Lotka (1925), Volterra (1926) g ( ) = a (1 + ah ) Холлинга II-го типа Holling (1959) g ( ) = g max (1 e ) s Ивлева Ивлев (1955), Ivlev (1961) g(, P) = Хассела–Варлей (HV) Hassell, Varley (1969) Pm P m Хассела–Варлей– Sutherland (1983), g(, P) = Холлинга (HVH) Arditi, Akakaya (1990) 1 + h P m Беддингтона– 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 g(, P) = 1 (1 ) Trn (2008) “дележа жертв” Трана P (1 e ) Гибридная модель g(, P) = P Trn (2008) “выедания жертв” Трана P Обобщенная RD a 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, L x ] [0, L y ]. Изменение положения j -го хищника x P t за временной шаг j, описывается уравнением:
x P,t +1 = x P,t + P,t + v P,t, ( j = 1,…, P ).
j j j j Здесь P t – случайная величина, равномерно распределенная внутри круга радиуса j,, характеризующего интенсивность ненаправленных перемещений, а скорость P таксиса v P t вычисляется как v P,t = P S t (x P,t ) + P v P,t 1, где P 0 и P [0;
1] – j, j j j коэффициенты таксиса и инерции, а S t (x j,t ) – градиент “запаха” жертв, P являющегося аттрактантом для хищников. Предполагается, что запах каждой жертвы распределен нормально:
xx 1 exp S i,t (x ) = i,t, ( ) 2 2 а запах всей популяции вычисляется как сумма S t (x ) = i =1 S i,t (x ). Пространственные перемещения жертв описываются аналогичным образом, с возможностью включения в модель их отрицательного таксиса против градиента запаха хищников.
Жертва, оказавшаяся в R -окрестности хищника съедается с вероятностью 1.
При различных предположениях относительно пространственной активности видов (различных соотношениях их коэффициентов, и ), перемещения индивидуумов отслеживались в течение T итераций ( t = 1, 2,…, T ;
T = 200 ) с подсчетом среднего количества жертв, пожираемых хищником в единицу времени.
Эта процедура выполнялась 500 раз для всевозможных сочетаний количества жертв и хищников (, P = 1, 2,…,25 ). Постоянство численностей популяций обеспечивалось тем, что съеденная жертва сразу воскресала в случайно выбранной точке.
Аналогичный прием с заменой жертвы используется и в лабораторных экспериментах по определению рациона насекомых (Tully et al. 2005). Результатом счета является таблица средних рационов G = (g, P )25,P =1, представляющих собой значения трофической функции g (, P ), которая не задается a-priory, а является результатом динамики системы.
Эксперименты показали, что для случайно перемещающихся животных трофическая функция не зависит от численности хищников, однако, при способности хищников к активному преследованию случайно перемещающихся жертв ( P 0, = 0 ) наблюдается уменьшение рациона хищников с увеличением размера их популяции – типичное проявление интерференции (см. рис. 10).
Полученные значения g, P аппроксимировались теоретической зависимостью HVH (см. табл. 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-функция, имеющая всего два параметра, тоже достаточно хорошо описывает трофические взаимодействия в данной системе.
а б Рис. 10. Изотрофы g (, P ) = const (вверху), полученные в экспериментах с (а) P = 0.2, P = 0. (оцененные коэффициенты HVH-функции – = 0.1 ;
m = 1.0 ;
h = 25.7 ) и (б) P = 0.9, P = ( = 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. Модель – система “реакция–диффузия”:
t = Wss b( 2) µ + +, ss ss rs ss ss ss t = Wrs 2b( 2 )( ) µ + + 2+, (4.1) rs ss rs rs rr rs rs rs t = Wrr b( rr ) µ + 2+ rr, rr rs rr rr n = n = n = 0, x, ss rs rr где ij = ij (x, t ) плотность генотипа ij в точке x в момент времени t ( i, j = r или = ss + rs + rr общая плотность популяции;
параметры b, µ и – s), коэффициенты плодовитости, смертности и внутривидовой конкуренции вредителя.
При моделировании стратегии “высокая доза–убежище” полагается, что ареал вредителя состоит из участков, засеянных либо Bt-кукурузой ( Bt ), либо обыкновенной кукурузой ( Ref ). Приспособленности генотипов (выживаемость личинок), согласно (Bourguet et al. 2000;
Vacher et al. 2003) представляется в виде:
1 hc c, x ref ;
1, x ref ;
Wrr = 1 c, на всем ;
Wrs (x ) = W ss (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) дает уравнение логистического роста популяции. Представление частот аллелей через плотности генотипов p s = ( ss + rs 2), pr = ( rr + rs 2), позволяет доказать теорему о сохранении постоянства частот. Также показано, что для равновесий демо-генетической модели выполнено условие Харди–Вайнберга частота генотипа ij. Данные результаты подтверждают u ss urr = urs 4 где uij = ij адекватность демо-генетической модели (1.4) в пространственно однородном случае.
В случае неоднородности области (сочетания Bt-участков и убежищ) справедлива следующая теорема:
Теорема 4.1. В биологической системе, описываемой пространственной демо генетической моделью (4.1), помимо диффузионного распространения аллельных частот в пространственно-неоднородной среде возникает адвективный поток аллельных частот, индуцируемый неоднородностью пространственного распределения общей плотности популяции и аллельных частот pr и ps:
p r t = bp r (Wr W ) + p r + 2 p r ln ;
t = (bW (µ + )) + ;
p s + p r = 1;
Wr = Wrs p s + Wrr p r ;
W = Wss p s2 + 2Wrs p s p r + Wrr p r2 ;
(4.3) p r n = n = 0, x..
Модель (4.3) эквивалентна (4.1) и отличается от классических уравнений математической генетики, формально объединенных с диффузией (пространственной модели Фишера–Холдена–Райта), только членом 2 pr ln, который влияет на распространение Bt-устойчивой аллели r и может быть интерпретирован как адвективный член, описывающий направленный поток аллельной частоты pr со скоростью 2 ln. Эта адвекция индуцируется неоднородностью пространственного распределения и pr и исчезает, если распределение или pr однородно.
Следствие. Диффузионная модель Фишера–Холдена–Райта адекватно описывает пространственное распространение гена Bt-устойчивости лишь в случае однородного распределения общей плотности в пределах всего моделируемого поля.
Теорема 4.2. При эволюции генетической структуры популяции, описываемой пространственной демо-генетической моделью (4.1), отклонение системы от равновесия Харди–Вайнберга = u ss u rr u rs 4 (Свирежев, Пасеков 1982) описывается дифференциальным уравнением:
( ) t = b p s2 p r2 (Wss + Wrr 2Wrs ) W + + 2 ln + 2 p r. (4.4) Согласно (4.4), система (4.1) стремится к равновесию Харди–Вайнберга ( 0, t ) только, если одна из аллельных частот ( p r или p s ) стремится к нулю. В противном случае система (4.1) эволюционирует за пределами равновесия ХардиВайнберга к полиморфному состоянию. Кроме того, отклонение увеличивается за счет неоднородности распределения аллельных частот.
ref, Bt, убежище Bt поле Lx 0 xref Рис. 11. Шаблон убежища для одномерного ареала = [0, L x ].
Для нахождения стационарных режимов пространственной демо генетической модели (4.1) на одномерном (1D) ареале = [0, Lx ], состоящем из двух участков: Bt-поля и убежища (рис. 11), решалась стационарная краевая задача следующего вида:
d 2 1 b = Wss ss + rs ss µ ss ;
ss dx 1 2b d rs ss + rs rs + rr rs µ rs ;
= Wrs 2 2 (4.5) dx 2 d rr = 1 W b rs + rr µ rr, dx 2 rr rr 2 d ss dx x =0, L = d rs dx x = 0, L = d rr dx x = 0, L = 0, x [0, L x ].
x x x Для обоснования корректности постановки задачи (4.5) были сформулированы и доказаны теоремы о существовании ее решений:
Теорема 4.3. Всегда существуют стационарные однородные по пространству решения системы (4.1) вида ss (x ) = 0, rs (x ) = 0, rr (x ) = K = (bWrr µ ), x [0, Lx ], * * * где Wrr = 1 c согласно формулам (4.2).
Теорема 4.4 (необходимое условие). Если существуют стационарные неоднородные по пространству решения системы (4.1) для конфигурации поля, представленной на рис. 11, такие, что x [0, Lx ] ss (x ) (0, K ), а rs (x ) = rr (x ) 0, и * * * при этом Wss (x ) = 0 для x Bt и Wss (x ) = 1 для x ref, то выполняется условие:
(K 2 3) b, 2 (4.6) 0 0 ref где 0 ref, 0 = (0), ref = (xref ) ( xref (0, Lx ) – точка, разграничивающая Bt-поле и убежище).
В ходе доказательства был определен вид и способ построения решения (рис. 12): т.к. коэффициент приспособленности Wss в (4.5) определяет качественное отличие динамики системы в убежище ( ref ) и на Bt-поле ( Bt ), задача решалась на каждом из этих участков, затем два решения “сшивались” в точке xref (точка В на рис. 13).
Рис. 12. Стационарное решение системы (4.1) такое, что x [0, L x ] ss ( x ) (0, K ), * а (x ) = rr (x ) = 0. xref точка перегиба, разграничивающая убежище (20%) и Bt-поле (80%).
* * rs ( x ) = ss (x ) x [0, Lx ]. Здесь Рис. 13. Фазовые портреты редуцированной системы (4.5):
( ) ( ) * 2 = Wss b µ, Wss Wss. (а) Wss * 1 = Wss b µ, µ b, Wss µ b ;
(б) Wss µ b, Bt Bt ref ref ref Bt ref Wss µ b. Черные линии – фазовые траектории системы в убежище, серые линии – на Bt-поле.
Bt Пунктирные линии – сепаратрисы. Дуга AB соответствует участку замкнутой фазовой траектории, вдоль которой образующая точка движется в убежище. Дуга BC соответствует участку незамкнутой фазовой траектории гиперболического типа, вдоль которой образующая точка движется на Bt-поле. Точка А соответствует левому краевому условию задачи (4.5), точка С – правому. Точка В – точка “сшивки” решения в убежище с решением на Bt-поле.
Верна и более общая теорема:
Теорема 4.5 (необходимое условие). Если существуют стационарные неоднородные по пространству решения системы (4.1) для конфигурации поля, ss ( x ) (0, ( ss b µ ) ), представленной на рис. 11, такие, что x [0, Lx ] * W ref а и Wss µ b ( Wss = Wss (x ) x, rs ( x ) = rr ( x ) = 0, * * и при этом Wss Wss ref Bt ref ref ref = Wss (x ) x ), то выполняется условие:
Bt W ss Bt ( ), Wss b µ 2 b Wss Wss ref ref Bt 2 (4.7) ref где 0 ref, 0 = (0), ref = (x ref ) ( x ref (0, L x ) – точка, разграничивающая Bt-поле и убежище) (рис. 12).
Для численного исследования пространственно-временной динамики популяции кукурузного мотылька, использовались следующие параметры модели (4.1):
-1 -1 b = 9.94 год, µ = 6.84 год, K = 147.4 10 6 особей/км. Начальное распределение популяции мотылька предполагалось однородным ( 0 = 2.948 10 6 особей/км2), а частота Bt-устойчивого гена – малой pr0 = 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-устойчивого гена ( p r = 0 ), слегка увеличив размер убежища. Таким образом, подвижность вредителя играет определяющую роль в эффективности стратегии “высокая дозаубежище”.
Рис. 14. Частота Bt-устойчивой аллели pr в зависимости от размера убежища (в %) и коэффициента диффузии для демо-генетической модели (4.1) (а) и модели ФишераХолденаРайта (б). L x = 16 км. Сплошными линиями отмечены ветви устойчивых равновесий систем, пунктирными линиями – ветви неустойчивых равновесий. Заданные генетические параметры соответствуют стратегии “высокая доза–убежище”: = 1, h = 0.14, hc = 0.2, c = 0.1 (особи, обладающие геном Bt-устойчивости, платят за свою приспособленность к токсичной среде 10% от коэффициента выживаемости Wrr ).
Прогноз диффузионной модели ФишераХолденаРайта (рис. 14б) только в случае нулевой и бесконечно большой подвижности вредителя ( = 0 и ) совпадает с демо-генетической моделью (4.1). Согласно фишеровской модели, эффективное применение стратегии “высокая доза – убежища” с целью предотвращения адаптации мотылька к Bt-кукурузе ( p r = 0 ) требует слишком больших размеров убежищ, близких к 100% (98% от всего поля для = 0.1 км2год-1), в то время как в демо–генетической модели (4.1) полностью предотвратить развитие Bt устойчивости у мотылька удается при гораздо меньших размерах убежища (21% от всего поля для = 0.1 км2год-1).
а) б) адвекция адвекция диффузия pr pr диффузия хref 0 хref Lx Lx Рис. 15. Схема пространственного распределения частоты p r и общей плотности в демо генетической модели (4.1) на 1D ареале [см. (4.3)], когда (а) плотность Bt-устойчивых особей вредителя на Bt-поле очень низка ( 0 ) и (б) плотность Bt-устойчивых особей вредителя на Bt поле возросла в результате распространения в популяции Bt-устойчивого гена.
Таблица 2. Время T10 (годы), за которое частота Bt-устойчивой аллели в популяции мотылька достигает 10% в демо-генетической модели (4.1) при разных размерах убежища и значениях коэффициента диффузии мотылька. В скобках даны значения времени для диффузионной модели ФишераХолденаРайта. Размер поля: 16 км 16 км. Генетические параметры соответствуют стратегии “высокая доза–убежище”: = 1, c = 0, h = 0.05, hc = 0.2.
Размер убежища (%) Диффузия, (км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-устойчивой аллели p r последовательно достигает 0.002, 0.1, 0.7 и 1, для 20% убежища (рис. 11) в демо-генетической и фишеровской моделях. Коэффициент диффузии = 0.1.