Delay differential equations as a tool for mathematical modelling of population dynamic

Full Text

ВВЕДЕНИЕ

Математическое моделирование в биологии

В современной науке математическое моделирование играет особую роль. Математические модели – это язык, на котором формулируются наши представления о явлениях в живой и неживой природе [Ризниченко, 2016, с. 5]. Уже, как минимум, несколько десятилетий естествознание характеризуется глубоким проникновением математических методов в различные области биологии [Рубин и др., 1987, с. 3]. Во-первых, биология, химия и физика переплетаются столь тесно, что любое биологическое утверждение (как новое, так и уже известное) нуждается в сопоставлении с законами физики и химии. А для этого необходимо использование математического аппарата [Романовский и др., 1975, с. 10]. Во-вторых, решение вопросов прикладной экологии, агро- и биометеорологии – таких, как повышение продуктивности агроэкосистем, рациональное управление природными ресурсами и охрана окружающей среды, – требует наличия в первую очередь подробно разработанных моделей динамики элементарных компонентов экосистемы [Полуэктов и др., 1980, с. 4]. Растущий интерес к реалистичным и пригодным для применения математическим моделям в популяционной биологии, идет ли речь о росте численности населения земного шара (с возрастным распределением или без него) или о сокращении численности редких видов растений и животных, или о динамике размножения бактерий и вирусов и т.д., объясняется значимостью подобных моделей для понимания происходящих процессов и построения практических прогнозов [Мюррей, 2009, c. 26].

При исследовании сложных биогеоценозов (экологических систем) плодотворными оказались методы общесистемного подхода. К экологическим системам оказался также применим принцип изоморфизма, позволяющий описывать сходными математическими уравнениями системы, разные по своей природе, но одинаковые по структуре и типу взаимодействия между элементами, их составляющими. В данном случае имеется большое сходство систем уравнений химической кинетики и межпопуляционной динамики^[14] [Рубин и др., 1987, с. 241].

Математической основой химической кинетики служит теория дифференциальных уравнений [Романовский и др., 1975, с. 10]. И в экологических исследованиях рабочим математическим аппаратом являются, как правило, дифференциальные уравнения (обыкновенные [Базыкин и др., 1980; Дорофеев и др., 1992; Panikov et al., 1992; Зинченко, 2017] и в частных производных [Орлов и др., 1987; Свирежев, 1987; Мюррей, 2011; Глаголев и др., 2016]). Однако при этом исследователи нередко встречаются с трудностями, связанными с ограниченными возможностями имеющегося математического аппарата, который долгое время развивался под влиянием задач физики, механики и техники [Галкин и др., 1981, с. 5, 7]. Ниже мы рассмотрим в этой связи примеры описания динамики популяций, причем (кроме базового описания единственной популяции) ограничимся частным случаем межпопуляционной динамики – моделями «хищник-жертва».

Математическому моделированию последней задачи посвящена значительная литература (не считая классической монографии Вольтерра [1976], это, например, работы [Leigh, 1968^[15]; Shoemaker, 1977; Базыкин и др., 1980; Bazin, 1981] и многие др.). В большинстве исследований изучались модели, представляющие собой систему двух обыкновенных дифференциальных уравнений [Бурд, 1985] – см., например, [Буриев и Розет, 1985; Казакова, 1993; Leffelaar, 1993, p. 55-56]. При этом рассуждения проводились так, словно будущая эволюция биологической системы полностью определяется знанием ее настоящего состояния (прошедшее влияния не оказывает) [Вольтерра, 1976, с. 162]. Серьезным недостатком упомянутых выше моделей является использование мгновенных значений рождаемости и смертности, полностью определяемых состоянием популяции в данный момент [Мэрди, 1979, c. 113]. Однако в реальных экосистемах это не так. На самом деле всегда имеется некоторое запаздывание в регуляции численности, которое может быть вызвано несколькими причинами [Рубин и др., 1987, с. 251]. Например, согласно R. May^[16], часто на рост популяции оказывает воздействие количество пищи не только в настоящий, но и в предшествующий момент времени [Горбенко и Крышев, 1985, с. 96]. С формальной точки зрения, независимо от причины, запаздывание можно учесть если строить модели не на основе дифференциальных уравнений, а использовать «дифференциальные уравнения с отклоняющимся аргументом» (причем, раз речь идет о запаздывании, то можно ограничиться лишь одним их типом – «дифференциальными уравнениями с запаздыванием»).

Некоторые понятия теории уравнений с отклоняющимся аргументом

Дифференциальными уравнениями с отклоняющимся аргументом  (ДУСОтА)^[17] называются дифференциальные уравнения, в которых неизвестная функция и ее производные входят при различных значениях аргумента. Обычно такие уравнения делятся на три типа – см. табл. [Митропольский и Мартынюк, 1979, с. 6]. Мы в данной статье ограничимся только уравнениями запаздывающего типа (как наиболее важными для математической биологии).

Как видим, в отличие от дифференциальных уравнений, где производная зависит от решения в настоящий момент времени, в системе дифференциальных уравнений с запаздывающим аргументом^[18] производная зависит также от значений решений в предыдущие моменты времени, т.е. уравнения имеют вид

dу(t)/dt = f(t, у(t), у(t-τ₁), у(t-τ₂), ..., у(t-τ_k)), (1)

где запаздывания τ_j > 0 [Шампайн и др., 2009, с. 243]. В общем случае запаздывания могут меняться с течением времени, т.е τ_j = τ_j(t).

Но обширный и очень важный класс образуют системы уравнений с постоянными запаздываниями. Более того С.Т.Н. Baker с соавт.^[19] отмечают, что при моделировании систем на практике в качестве функций запаздывания константы выбираются чаще всего [Шампайн и др., 2009, с. 243].

Поскольку правая часть уравнения содержит члены вида у(t-τ_j), которые соответствуют значениям решения, вычисляемым в моменты времени, предшествующие начальной точке, то заданные начальные данные должны включать не только значение решения у(а) в начальной точке t = а, но и «историю» изменения этого решения. Значения y(t) для всех t Î [а-Т, а]^[20] составляют функцию, которая называется функцией истории [Шампайн и др., 2009, с. 244].

 

Таблица. Классификация дифференциальных уравнений с отклоняющимся аргументом.

Тип	Определение	Пример	Примечание
Уравнения с запаздываю-щим аргумен-том (или запаз-дывающего типа)	Такое уравнение, в котором производная максимального порядка от неизвестной функции входит при одинаковых значениях аргумента и этот аргумент не меньше аргументов неизвестной функции и ее других производных, входящих в уравнение	dу(t)/dt =  = f(t, у(t), у(t-τ(t))),   где τ(t)≥0	Широко используются в математических моделях биологии (особенно – в популяционной динамике).
Уравнения с опережающим аргументом (или опережа-ющего типа)	Такое уравнение, в котором производная максимального порядка от неизвестной функции входит при одинаковых значениях аргумента и этот аргумент не больше аргументов неизвестной функции и ее других производных, входящих в уравнение	dу(t)/dt =  = f(t, у(t), у(t+τ(t))),   где τ(t)≥0	Нам не известно применение уравнений это-го типа в биологии. Вероятно, они могли бы использоваться для моделирования целена-правленного поведения высших животных, если предположить, что животное способно (в определенных границах) действовать в на-стоящий момент времени в соответствии со своим ожидаемым состоянием в будущем. Но поскольку ожидаемое состояние может не совпасть с тем, которое реально удасться достигнуть, то более естественными здесь, по-видимому, будут стохастические уравнения опережающего типа.
Уравнения нейтрального типа	Дифференциальные уравнения с отклоняющимся аргументом, не относящиеся к уравнениям запаздывающего или опережающего типа	dу(t)/dt =  = у(t-2) - dу(t-1)/dt	Задачи, приводящие к таким уравнениям часто появляются при изучении двух и бо-лее колебательных систем с некоторыми связями между ними [Хейл, 1984, с. 13]. Но, хотя в биологии колебательные систе-мы со связями распространены весьма ши-роко, мы не встречали уравнения нейтраль-ного типа в биологических приложениях.

 

Используемые сокращения

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

Уравнение Goudriaan-Roermund (введение запаздывания в уравнение Мальтуса)

Модели с запаздыванием впервые появились в экологии, где и сейчас применяются наиболее широко [Бабский и Мышкис, 1983, с. 383]. По-видимому, одними из первых ввели в рассмотрение дифференциальные уравнения с отклоняющимся аргументом в качестве моделей биологических процессов с запаздыванием А.J. Lotka и F.R. Sharpe^[21]. Поэтому наш обзор моделей с запаздыванием мы начнем именно с экологических приложений.

В свое время Одум подчеркивал, что «простая ситуация, когда сопротивление среды возрастает линейно при увеличении плотности, по-видимому, имеет место только в популяциях организмов с простым жизненным циклом». В популяциях высокоорганизованных растений и животных существенное влияние на динамику плотности оказывает внутренняя структура (демографическая и генетическая) и эффекты запаздывания в реакции популяции на изменения абиотических компонентов экосистемы. Последнее означает, что изменение количества доступных ресурсов, условий среды обитания или внутрипопуляционных характеристик сказывается на рождаемости, смертности и миграции не мгновенно, а только через некоторый промежуток времени. Следовательно, и в моделях необходимо учитывать эффекты запаздывания [Полуэктов и др., 1980, с. 22, 112]. В частности, имеется промежуток времени между зачатием и рождением. В природе наличие такого запаздывания по времени легко проследить. Так, например, у млекопитающих имеются длительные периоды созревания, вследствие чего число созревающих особей зависит от численности популяции в предыдущий период размножения, и поэтому на рождаемость в момент времени t влияют условия, существовавшие в момент времени t-τ [Мэрди, 1979, c. 113]. Если какое-то изменение в окружающей среде, например увеличение ресурсов вызовет внезапное повышение продуктивности взрослых особей, то соответствующее изменение численности взрослых особей произойдет лишь по прошествии времени τ₂ [Рубин и др., 1987, с. 251].

Согласно [Смит, 1976, с. 53; Рубин и др., 1987, с. 251; Ризниченко, 2016, с. 91] это означает, что уравнение

dу(t)/dt = f(у(t)), (2)

где у – численность взрослых особей, следует заменить уравнением dх(t)/dt = f(х(t-τ₂)), здесь х(t-τ₂) – численность половозрелых особей в момент времени t-τ₂. Если в качестве (2) рассматривать закон Мальтуса (dy/dt = μ·y), то, казалось бы, мы придем к уравнению^[22] dх(t)/dt = μ·(х(t-τ₂)), где μ – удельная скорость роста.

Но следует помнить, что в уравнении Мальтуса, μ представляет собой разность рождаемости (b₀) и смертности (d₀). Если раскрыть это уравнение, то в правой части мы будем иметь разность скоростей рождения (b₀·y) и гибели (d₀·y) [Полуэктов и др., 1980, с. 110-111; Мюррей, 2009, с. 27]. То, что скорость рождения популяции зависит от ее состояния в момент времени t-τ₂, обсуждалось выше. Но разве скорость гибели будет зависеть тоже от него же?! Мэрди [1979, c. 113] полагал, что с меньшей ошибкой будет связано предположение о зависимости смертности от состояния популяции в данный момент, поэтому уравнение Мальтуса превращается вот в такое [Goudriaan and Roermund, 1993, p. 89]:

dу(t)/dt = b₀·у(t-τ₂) - d₀·у(t). (3)

 

Итак, при учете запаздывания, обусловленного промежутком времени между зачатием и рождением, уравнение Мальтуса заменяется уравнением Goudriaan-Roermund. Впрочем, если отвлечься от биологического смысла коэффициентов, то область применения этого уравнения оказывается шире, поскольку к нему сводятся некоторые другие уравнения популяционной динамики.

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

 

dХ(t)/dt = 0.5·k·s·w - d₀·Х(t) - 0.5·m·k·s·Х(t-τ₂),

 

где d₀ – удельная скорость гибели взрослых особей (она составляет некоторую постоянную величину, т.е. не учитывается старение, а также возможности повышения смертности взрослых особей при увеличении численности у и сокращении количества корма на одну взрослую особь); m – скорость поступления корма, необходимая для поддержания жизни каждой взрослой особи (m = const); s –вероятность того, что каждое данное яйцо выживет и из него получится взрослая особь (допущение о постоянстве этой вероятности при наличии неограниченного количества корма для личинок представляется вполне оправданным); Х(t) – численность взрослых особей Lucilia в момент времени t; w – постоянная скорость выдачи корма этим особям (следовательно, за время dt количество корма, приходящееся на одну особь, равно w·dt/Х); τ₂ – время, необходимое для того, чтобы из яйца развилась взрослая особь. При выводе уравнения предполагалось, что половозрелые самки превращают избыток корма в яйца с некоторой постоянной эффективностью k, т.е. число отложенных яиц за время dt, приходящееся на одну самку, составит k·(w/Х - m)·dt, а общее число яиц, отложенных в популяции будет равно 0.5·Х·k·(w/Х - m)·dt = 0.5·k·(w - m·Х)·dt (при условии, что число самцов и самок одинаково). После ряда несложных преобразований данное уравнение приводится к (3), где b₀ = ‑0.5·m·k·s; у = Х/Х_С - 1, здесь Х_С = 0.5·w/[d₀/(k·s) + 0.5·m] – равновесная плотность популяции, т.е. такая, что если Х = const = Х_С в течение времени τ₂, то и впоследствии всегда Х = Х_С [Смит, 1976, с. 55-57; Рубин и др., 1987, с. 253-254].

Однако «выше головы не прыгнешь» – рассмотренное нами уравнение (3) является линейным и оно оказывается неприменимо для описания нелинейных взаимодействий в популяции, даже таких простых, которые описывает одно из самых, пожалуй, известных в динамической экологии уравнений – логистическое. Так что теперь вполне естественным будет посмотреть, как преобразуется логистическое уравнение в случае учета запаздывания.

Уравнения Хатчинсона и Мэрди (введение запаздывания в логистическое уравнение)

В реальных популяциях интенсивность таких процессов как отравление среды продуктами метаболизма, каннибализм и т.п., зависит от численности взрослых особей, т.е. отрицательное влияние на коэффициент естественного прироста оказывают особи предыдущего поколения. С учетом этих явлений логистическое уравнение (dN/dt = N·[b₀ - δ·N])^[23] перепишется в виде [Рубин и др., 1987, с. 252] так называемого уравнения Хатчинсона^[24]:

 

dN(t)/dt = N(t)·[b₀ - δ·N(t-τ₁)]. (4)

 

Мы ввели здесь новое обозначение N для численности всей популяции, оставив прежнее у – для взрослых особей^[25]. Вообще говоря, как это часто бывало в истории, научное сообщество связало название уравнение вовсе не с первооткрывателями, а назначило его довольно случайным образом.

Задолго до 1948 г., когда появилась работа Хатчинсона, два экономиста – Фриш и Холм – изучали такое же уравнение в связи с анализом цикличности деловой активности. Только в 1955 г. математики наконец открыли уравнение Хатчинсона-Райта. Райт^[26] изучал его в связи с распределением простых чисел [Хэссард и др., 1985, с. 134], причем в 1945 г. писал об этом уравнении: «Лорд Черуэлл привлек мое внимание к одному уравнению…, которое он встретил, занимаясь приложением вероятностных методов к задаче о распределении простых чисел» [Хайрер и др., 1990, с. 311]. И вот почему-то Фриша и Холма (1935 г.) забыли, лорда Черуэлла и Каннингема^[27] тоже забыли, а запомнили – Хатчинсона (1948 г.) и иногда вспоминают Райта (публикация 1955 г., т.е. ровно через 20 лет после первооткрывателей).

Простое обоснование возможности колебательного режима для частного случая уравнения Хатчинсона дано в [Марри, 1983, с. 185], а срогое качественное исследование поведения его решений приводится в монографии Ю.М. Свирежева и О.Д. Логофета^[28] (основные свойства решений кратко обсуждаются в [Полуэктов и др., 1980, с. 158^[29]; Косолапова и Ковров, 1988, с. 10]). Численные решения для уравнения Хатчинсона в безразмерных переменных при различных значениях параметра^[30] приведены в [Хайрер и др., 1990, с. 311-313]. Большая группа работ, посвященных уравнению Хатчинсона, выполнена Ю.С. Колесовым и его сотрудниками. Кроме того Ю.С. Колесов показал, что это уравнение достаточно хорошо описывает динамику изменения численности многих массовых видов животных, так как у таких видов средняя продолжительность жизни обычно лишь немногим больше τ₁. Если же это условие нарушается, то, согласно указанному исследователю, в простейшем случае можно перейти к так называемому «обобщенному уравнению Хатчинсона» (в безразмерном времени)

dN(t)/dt = b₀·N(t)·[1 - {α₁·N(t-τ₁) + α₂·N(t-τ₁-1)}/К],

где α₁·α₂ ≥ 0, α₁ + α₂ = 1. Заметим, что при моделировании популяционной динамики исследователи предлагали различные обобщения уравнения Хатчинсона. Так Braddock и van den Driessche  изучали существование и устойчивость периодических решений (а также их области притяжения) для уравнения

dN(t)/dt = b₀·N(t)·[1 - β·N(t) - γ·N(t-τ₁)],

Mufti изучал уравнение

dN(t)/dt = b₀·N(t)·[1 - {N(t-τ₁)/К}^θ] + р,

а Брауэр – уравнение

dN(t)/dt = N(t)·f(N(t-τ₁)) - E

 

с точки зрения устойчивости и временных характеристик его решения в зависимости от Е, т.е. от постоянного снятия урожая [Бабский и Мышкис, 1983, с. 384-385]. Наконец, объединяя идеи моделей Goudriaan-Roermund и Хатчинсона, Мэрди [1979, c. 113] дает общее популяционное уравнение с запаздыванием:

dN(t)/dt = b₀·[1 - N(t-τ₁)/K]·N(t-τ₂) - d₀·N(t), (5)

где τ₁ и τ₂ – времена запаздывания, а b₀ и d₀ – соответственно «основные значения»^[31] рождаемости и смертности. То, что все вышеприведенные уравнения являются обобщениями уравнения Хатчинсона, совершенно очевидно, ибо они сводятся к этому уравнению в частных случаях: уравнение Мэрди – при d₀ = τ₂ = 0; уравнение Брауэра – при f(N(t-τ₁)) = b₀ - δ·N(t-τ₁) и E = 0; уравнение Mufti – при θ = 1 и р = 0; уравнение Braddock-van den Driessche – при β = 0.

Исследовалась также задача о влиянии неоднородности среды обитания на динамику численности вида при учете запаздывания, в частности, уравнение Хатчинсона с диффузией и ряд других [Бабский и Мышкис, 1983, с. 387-388; Колесов, 2001; Алешин и др., 2017]. Такие модели представляют собой уравнения в частных производных с запаздыванием, а не обыкновенные ДУСОтА, которыми мы решили ограничиться в данной статье.

В целом предмет теории уравнений с запаздывающим аргументом весьма обширен, как и приложения этой теории к экологическим моделям. Если говорить именно о приложениях к экологии, то кроме уже рассмотренной нами литературы нельзя не упомянуть книгу Мэя^[32], исследование Стирзакера (модель с двумя запаздываниями, описывающую популяцию полевок Microtus agrestis; на русском языке биологические механизмы, лежащие в основании этой модели, кратко описаны в [Хэссард и др., 1985, с. 126, 139]) и модель хищник-жертва Макдональда^[33], а также [Бабский и Мышкис, 1983; Мюррей, 2009, с. 42-53]^[34]. Отметим, что большой вклад в приложение ДУСОтА к экологии внесли отечественные ученые, среди которых упомянем Юрия Серафимовича Колесова и Сергея Александровича Кащенко (см., например, [Колесов и Швитра, 1979; Kolesov et al., 2010; Кащенко, 2017; 2017а]). К моделям экологии весьма близки и уравнения демографии, которые также могут иметь форму ДУСОтА – см., например, [Пэнтл, 1979, с. 75]. Более того, модели с запаздыванием не ограничиваются только сферой экологии. И, хотя это не является основной целью нашей статьи, мы кратко упомянем некоторые другие сферы применения ДУСОтА в биологии в одном из разделов ниже.

 

Запаздывание в моделях межпопуляционных взаимодействий

При рассмотрении уравнений динамики биологической системы с учетом теории встреч традиционно предполагалось, что небезразличные встречи между индивидуумами разных видов оказывают влияние посредством немедленного изменения численности индивидуумов. Однако это законно для истребляемого вида и совсем не так для вида, который извлекает выгоду из встречи. Очевидно, что благоприятное действие встречи может проявляться только с некоторым запаздыванием^[35]. Состояние биологической системы в данный момент должно зависеть от встреч, имевших место в течение более или менее длительного периода, предшествующего этому моменту [Вольтерра, 1976, с. 163]. Можно предположить, что увеличение числа хищников после поедания жертв увеличивается спустя некоторое время τ. Эту инерционность позволяет учесть запаздывание аргумента^[36] [Дьяконов, 2008, c. 429]. По-видимому, одними из первых, кто предожил модель «хищник-жертва» в виде системы ДУСОтА, были Вангерски и Каннингем^[37]. Их модель имела вид^[38]:

 

 dN₁(t)/dt = r₁·N₁(t)·[m-N₁(t)]/m - α·N₂(t)·N₁(t),    dN₂(t)/dt = -d₀·N₂(t) + r·N₁(t-τ₁)·N₂(t-τ₁),

 

где N₁ и N₂ – количество особей в популяциях, соответственно, жертвы и хищника, а τ₁ имеет смысл осредненного интервала времени между моментом гибели одной особи жертвы и моментом соответствующего увеличения числа взрослых хищников [Бабский и Мышкис, 1983, с. 386].

В 1979 г. Колесовым (в предположении, что сопротивление внешней среды для каждого вида есть его внутренняя характеристика) была предложена модель

 

dN₁(t)/dt = r₁·{1 + α·[1 - N₂(t)/K₂] - N₁(t-τ₁)/K₁}·N₁(t),    dN₂(t)/dt = r₂·{N₁(t)/K₁ - N₂(t-τ₂)/K₂}·N₂(t),

 

где К₁ и К₂ – средние размеры популяций, соответственно, жертвы и хищника; r₁ и r₂ – мальтузианские коэффициенты роста^[39]; α – параметр, характеризующий давление популяции хищника на популяцию жертвы, причем предполагается, что при изменении α коэффициенты r₁ и К₁ меняются по законам r₁ = r_о/(1+α) и К₁ = К_о/(1+α), здесь r_о и К_о – соответственно, коэффициент роста жертвы и средний размер ее популяции при отсутствии хищника. Поведение решений вышеприведенных уравнений оказалось весьма сложным [Колесов и Кубышкин, 1980]!

Примерно к этому же времени оносится проникновение в математическое моделирование популяционной динамики так называемого «метода самоорганизации»^[40], в котором функция выбора вида модели передана ЭВМ^[41]. В основе построения самоорганизующихся моделей лежат два главных принципа – принцип «внешнего дополнения» и селекции уравнений по некоторому критерию. Согласно первому принципу, лучшим считается не то уравнение, которое минимизирует некоторый критерий (например минимум ошибки в методе наименьших квадратов) на последовательности точек, используемых для построения этого уравнения, а то, которое дает наименьшее значение критерия на проверочной последовательности (т.е. не участвовавшей в построении уравнения). А принцип селекции уравнений способствует сокращению полного перебора различных комбинаций исходных параметров. Принципы построения самоорганизующихся моделей воплотились, в частности, в методе группового учета аргументов, в котором в качестве функций чаще всего используется квадратичный полином [Галкин и др., 1981, с. 88].

Этот метод был применен для описания нелинейными ДУСОтА (оптимальной сложности) хорошо известных данных о численности рысей и зайцев в Канаде, собранных Компанией Гудзонова залива за 200 лет. В результате, выдающимся российским экологом – чл.-корр. РАН Розенбергом^[42] (в соавторстве с П.М. Брусиловским) было построено две модели, представлявших собой системы из двух уравнений вида

$d{N}_1\left(t\right)/dt = \sum _{i=1}^m(a_i·{у}_i), d{N}_2\left(t\right)/dt = \sum _{i=1}^n\left(b_i·z_i\right).$

В первой модели m = 30, n = 8, а члены у_i и z_i имеют вид, соответственно, N₁(t-1)^α·N₁(t-2)^β·N₂(t-2)^γ и N₁(t-2)^β·N₂(t-1)^γ, где α, β и γ – могут принимать целые значения от 0 до 4 включительно. Во второй модели m = 29 и n = 32, а члены у_i и z_i имеют вид W(t-1)^δ·W(t-2)^γ·N₁(t-1)^α·N₁(t-2)^β, где W(t-i) – значение солнечной активности в числах Вольфа c запаздыванием на i лет; δ может принимать целые значения от 0 до 4 включительно. Вторая модель (учитывающая солнечную активность) дала точность аппроксимации 12.9% для плотности популяции зайцев и 10.7% – рысей, а точность прогноза, соответственно, 29.5% и 20%, что лучше и аппроксимации, и прогноза 5 других моделей, проверенных авторами [Галкин и др., 1981, с. 84, 88-91].

К сожалению, крайне ограниченный объем журнальной статьи не позволяет нам подробно остановиться на характере и свойствах решений ДУСОтА. Для многих дифференциальных уравнений с запаздыванием характерен достаточно богатый спектр решений, которые весьма сложно зависят как от параметров модели, так и от начальных условий (начальной истории). Часто встречаются периодические и квазипериодические решения и многое другое, что привлекает исследователей – специалистов по математической биологии. Видя аналогию между поведением реальных процессов и получаемыми решениями, они пытаются нащупать механизмы, обуславливающие такую аналогию. Подробнее с характером и свойствами решений уравнений с запаздыванием, являющихся математическими моделями биологических процессов, заинтересованный читатель может познакомиться по богатой литературе, из которой мы упомянем [Колесов и Кубышкин, 1980; Бабский и Мышкис, 1983; Мюррей, 2009, с. 44-66; Хидиров, 2014]. Более полное и строгое изложение всего спектра решений ДУСОтА (правда, уже без связи с моделями биологии) можно найти в математической литературе, в частности в той, которую мы будем цитировать ниже (в разд. «Как практически работать с дифференциально-разностными моделями?»).

Важно, однако, не стать жертвой заблуждения о корректности и обоснованности модели только на основании того, что некоторые ее решения хорошо согласуются с данными эксперимента [Мюррей, 2009, с. 47]. Необходимо проверить модель на соответствие здравому смыслу в каких-то очевидных ситуациях. К сожалению, обе модели Розенберга-Брусиловского не выдерживают такой проверки, поскольку оказываются в значительной степени лишенными биологического смысла. Действительно, как нетрудно заметить, во второй из них скорость изменения плотности популяции ни для жертв, ни для хищников не зависит от количества хищников – соответствующее уравнение вообще не содержит членов, включающих в себя N₁(t), N₁(t-1) или N₁(t-2). Таким образом, даже если хищников не было вовсе, то они, согласно второму уравнению второй модели, невесть откуда возникнут, причем интенсивность этого «самозарождения» будет определяться только количеством жертв и солнечной активностью (но сколько бы хищников ни возникло, это не окажет никакого воздействия на жертвы, динамика численности которых определяется, согласно первому уравнению, также лишь солнечной активностью и численностью популяции жертв в прошлом). Итак, вторая модель представляет собой апофеоз абсурда (не удивительно, что применив ее, авторы приходят к закономерному выводу: «удовлетворительного описания получить так и не удалось» [Галкин и др., 1981, с. 91]). Но и первая модель немногим лучше! Хотя в ней скорости изменения численности популяций и жертв, и хищников, зависят, как это имеет место в действительности, и от количества жертв, и от количества хищников, но пресловутое самозарождение окопалось и в этой модели^[43].

Краткий обзор моделей с запаздыванием в других областях биологии

На широкую применимость в биологии уравнений, содержащих информацию о прошлом, указали Bellman и Danskin в своей монографии^[44], выполненной в Rand Corporation [Хейл, 1984, с. 10]. Разумеется, здесь мы не сможем полностью очертить сферу применения этих уравнений (да это и не является нашей целью), но приведем несколько известных примеров. Во-первых, это, конечно, иммунология и тесно связанная с ней эпидемиология.

При исследовании заболеваемости гонореей Cooke и Yorke^[45] в 1972 г. изучали уравнение

$dI\left(t\right)/dt=g\left(I\right(t-{\tau }_1\left)\right)-g\left(I\right(t-{\tau }_2\left)\right),$

где I представляет число носителей инфекции, а g – неотрицательная функция, обращающаяся в нуль вне компактного интервала. Годом позже, описывая распространение кори в городской среде, London и Yorke^[46] встретились с уравнением

$dS\left(t\right)/dt=-\beta \left(t\right)·S\left(t\right)·[2·\gamma +S\left(t-14\right)-S\left(t-12\right)]+\gamma ,$

где S(t) – количество восприимчивых к инфекции индивидов в момент t; γ – скорость, с которой восприимчивые индивиды включаются в популяцию; β(t) – функция, характеризующая популяцию; индивид, зараженный в момент времени t, будет заразным на  протяжении интервала времени от t+12 до t+14 [Хейл, 1984, с. 12]. Понятно, что запаздывания, вводимые эмпирически в модели эпидемиологии, обусловлены имунными процессами в организмах восприимчивых к инфекции индивидуумов, носителей и заболевших. Поэтому неудивительно, что приблизительно в эти же годы активно стали развиваться математические модели имуннологии.

При формировании даже простейшей модели имунной защиты используется следующее положение иммунологии. Считается, что антитело связывает антиген, образуя комплексы антитело-антиген. Пропорционально количеству таких комплексов в организме через время τ формируются плазматические клетки, которые осуществляют массовую выработку антител (каскадный процесс образования клона плазматических клеток длится от нескольких часов до нескольких дней)^[47] [Марчук, 1991, с. 20, 29-30]. Отсюда естественным образом и возникает член с запаздыванием в уравнениях математических моделей иммунологии.

Базовая математическая модель инфекционного заболевания, представляющая собой систему нелинейных обыкновенных дифференциальных уравнений с запаздывающим аргументом^[48], для описания и исследования закономерностей, общих для всех инфекционных заболеваний, была построена Г.И. Марчуком в 1975 г. Эта базовая модель оказалась полезной для механизмов развития смешанных инфекций, а также изучения общей картины течения инфекций и интерпретации некоторых результатов наблюдений (в частности, по биостимуляции иммунитета при хронических инфекциях). С формальной точки зрения эта модель является системой ДУСОтА из 4 уравнений [Романюха, 2012, с. 47-49, 133], но фактически она может быть представлена в виде системы из 3 взаимосвязанных уравнений и еще одного, которое можно проинтегрировать в дальнейшем, имея результаты интегрирования этой сокращенной системы. В [Хайрер и др., 1990, с. 317-318] приведен пример численного интегрирования указанной сокращенной системы (3 уравнений).

Но задача количественного моделирования вирусных и бактериальных инфекций определила необходимость построения более детальных математических моделей. Для этого в работе [Марчук, 1991] была сформулирована математическая модель противовирусного имунного ответа, а на ее основе и модель противобактериального имунного ответа, также представлявшие собой системы нелинейных ДУСОтА [Романюха, 2012, с. 47-49, 133]. Модель (в форме системы ДУСОтА) лекарственной терапии ВИЧ-инфекции см. в [Мюррей, 2009, с. 493-495], и еще одну модель аутоиммунного заболевания – в [Fall et al., 2002, p. 259-260]. Вообще же, достаточно полное представление о состоянии исследований по математическому моделированию в иммунологии до 80-х гг. ХХ в. включительно, можно получить из специализированных статей, монографий и сборников^[49] [Романовский и др., 1975, с. 180-181, 338; Марчук, 1991, с. 58-59], а после 80-х гг. – см. [Романюха, 2012] и ссылки на литературу там. На этом мы закончим рассмотрение моделей имуннологии и эпидемиологии. Какие же еще модели следует упомянуть, говоря о применении ДУСОтА в биологии?

Во-вторых, это некоторые модели (начиная, как минимум, с 1974 г.) регуляторных механизмов клеточного деления (начиная, по-видимому, с модели Tsanev-Sendov^[50]). Здесь члены с запаздываниями возникают в силу того, что временные взаимоотношения этих регуляторных механизмов можно учесть при помощи латентного периода реализации генетической информации (ДНК®белок), интервала времени, необходимого для перехода регуляторных белков из внешней среды в цитоплазму и других параметров запаздывания [Хидиров, 2014, с. 9].

Некоторые другие модели молекулярной биологии, в частности, простейшая модель регуляции синтеза белка, основанная на репрессорной регуляции активности генов с учетом временных взаимоотношений, также представляют собой системы ДУСОтА [Хидиров, 2014, с. 19; Смит, 2005, с. 139]. Гудвин в 1965 г. предложил модели синтеза белка с механизмами генетического управления, простейший из которых состоял в следующем. Регулирующий ген G производит информационную рибонуклеиновую кислоту (мРНК). Затем мРНК взаимодействует с рибосомами, и образуются молекулы фермента. Разумно, вообще говоря, ожидать некоторого запаздывания (τ₁), вызванного, например, диффузией в цитоплазме, т.е. временного интервала между моментом образования мРНК и моментом достижения ею рибосом, активно участвующих в синтезе фермента. Фермент, в свою очередь, катализирует реакцию с некоторым субстратом. Часть продукта этой реакции представляет собой репрессор. Репрессор поступает назад к регулирующему гену, подавляя его активность, так что мРНК не может производиться [Марри, 1983, с. 173, 183-184]. Похожую модель в 1975 г. предложил Хидиров [2014, с. 19-24]. Главное ее отличие от предыдущей состоит в том, что учитывается еще одно запаздывание: предполагается, что не только концентрация информационных РНК в цитоплазме в момент времени t пропорциональна концентрации информационной РНК в ядре в момент времени (t‑τ₁), но и концентрация регуляторных белков в ядре (место их действия) в момент времени t пропорциональна концентрации белков в цитоплазме (место их образования) в момент времени (t-τ₂). Аналогичная модель рассматривается и в [Fall et al., 2002, p. 247-250].

В дальнейшем эти работы привели к формулировке весьма общей модели: пусть в некотором ограниченном объеме существует N взаимосвязанных элементов – регуляторов, способных к восприятию и синтезу сигналов определенной природы; и пусть взаимосвязь между регуляторами осуществляется посредством указанных сигналов со средним временем прохождения обратной связи τ. Примерами регуляторов могут являться клетки, молекулы, гены, рибосомы и т.д. Величины Х_i(t), характеризующие количества синтезируемого сигнала, соответствующего i-му элементу (1 ≤ i ≤ N) в момент времени t в рамках этой модели удовлетворяют системе ДУСОтА (с запаздыванием τ; в случае необходимости можно взять вместо среднего τ функции от Х) [Хидиров, 2014, с. 33-35].

Сюда же примыкают и работы в области ферментативной и микробиологической кинетики. В [Хайрер и др., 1990, с. 315-316] рассматривается одна из моделей ферментативной кинетики, в которой конечный продукт линейной последовательности ферментативных реакций ингибирует первую реакцию в этой последовательности. Молекула ингибитора должна перемещаться в место нахождения регулирующего энзима посредством таких процессов, как диффузия или активный транспорт. Рассматривая эти медленные процессы как причину запаздывания, приходят не к ОДУ, как это обычно имеет место в классической ферментативной кинетике, а к системе ДУСОтА.

Понятно, что подобные механизмы должны отразится и в запаздывании отклика роста биомассы на изменение питания. Однако в клетке микроорганизма идет огромное количество ферментативных реакций, причем не только последовательных, но и формирующих биохимические циклы (например, цикл трикарбоновых кислот и др.). Подробно описать реальную систему ферментативных реакций микробной клетки вряд ли возможно (по крайней мере, это не было возможным в 60-х гг. ХХ в.). Потому в исследовании Ноака^[51] для учета конечного времени изменения скорости роста микроорганизмов при изменении уровня питания (S) удельная^[52] скорость роста (μ) выражалась просто запаздывающей функцией времени:

$\mu \left(t\right)={\mu }_m·S(t-\tau )/[K_S+S(t-\tau \left)\right].$  (6)

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

$dx\left(t\right)/dt=\sum _{i=0}^N[A_i·x(t-{\tau }_i\left)\right]$

является подходящей моделью для описания разведения красителя в центральной емкости, когда окрашенная вода циркулирует через большое количество трубок. Его приложение к распределению в человеческом теле меченного альбумина, проникающего из кровяного русла в тканевые жидкости и обратно, обсуждали Bailey и Reeve^[53] [Хейл, 1984, с. 11-12]. Еще две физиологические модели (дыхание Чейна-Стокса и регуляция гемопоэза) рассмотрены в [Мюррей, 2009, с. 54-66].

И, наконец, в-четвертых, это отдельные модели продукционных процессов. Например, Башалханов и др. [1988] рассматривают довольно сложную математическую модель продукционного процесса огурца, в которой скорости изменения максимума удельной листовой поверхности (S₁) и высоты растения (Н) в момент времени t обратно пропорциональны массе надземной части растения (y) и прямо пропорциональны скорости ее изменения не в тот же самый момент времени, а в отстоящий от него на одни сутки: dS₁/dt ~ [Н·y(t-1)]^-1·dy(t-1)/dt, dН/dt ~ [S₁·y(t-1)]^-1·dy(t-1)/dt.

Всегда ли модели с запаздыванием имеют биологический смысл?

ДУСОтА, в особенности с запаздывающим аргументом, описывающие процессы с последействием, находят много приложений: в теории автоматического управления, в теории автоколебательных систем (пример: колебания молоточка электромагнитного прерывателя), при изучении проблем, связанных с горением в ракетном двигателе, ряда экономических проблем и многих других [Норкин, 1965, с. 7, 11]. Но, как видим, приложения эти преимущественно технические или вообще связанные с деятельностью человека. В таких приложениях запаздывание возникает совершенно естественно и имеет ясный физический смысл.

Например, в системах автоматического регулирования запаздыванием является принципиально всегда имеющийся промежуток времени, который нужен системе для реагирования на входной импульс [Норкин, 1965, с. 7]. Но в некоторые области биологии ДУСОтА были, по нашему мнению, привнесены несколько искусственно. Удивительным образом, вместо осмысления этой искусственности, произошло едва ли не повальное увлечение биологов моделями, содержащими запаздывание. И это затушевывает многочисленные факты того, что рассматриваемые модели, подчас, приводят к биологически бессмысленным результатам. Причем сейчас мы говорим не о той фундаментальной биологической бессмысленности, которая встретилась выше в моделях Брусиловского-Розенберга (она проистекала вовсе не из-за введения запаздывания в модель), а именно о бессмысленности, порождаемой запаздыванием. К сожалению, это досадное обстоятельство практически не рассматривается в литературе (по крайней мере, нам не известны публикации, где бы этот вопрос обсуждался), а потому мы считаем нужным привлечь к нему внимание вдумчивого читателя. Проиллюстрируем указанную биологическую бессмыслицу на примере модели эксплуатируемой популяции.

Простейшая модель эксплуатируемой популяции, построенная на основе уравнения Мальтуса, как известно, имеет вид dy/dt = μ·y - с, где с ≥ 0 – интенсивность уничтожения особей. Совершенно очевидно, что за конечное время можно добиться того, что численность популяции станет равной нулю [Заславский и Полуэктов, 1988, с. 112]. После этого все время в дальнейшем будет оставаться y = 0 – действительно, если популяция уничтожена, то она не сможет из ничего возникнуть вновь (а члены, описывающие какие-либо процессы, способные возродить популяцию естественным путем, например – миграционные процессы, в данном простейшем уравнении отсутствуют). Будем считать, что с = const > 0 до тех пор, пока количество особей в популяции не упадет до нуля, после чего отлов навсегда прекращается (т.е. становится с = 0). Это совершенно не сложно реализовать в некоторых искусственных экосистемах в случае достаточно мелких животных (примером такой экосистемы может быть аквариум с рыбками).

Очевидно, модификация модели Goudriaan-Roermund для случая эксплуатируемой популяции (популяции с отловом) будет иметь вид:

$dу\left(t\right)/dt=b_0·у(t-{\tau }_2)-d_0·у\left(t\right)-с\left(t\right),$

где скорость отлова с(t) = с₀ = const > 0 при t < T и с(t) = 0 при t ≥ T (здесь T – наименьшее положительное значение времени, при котором у = 0). Чтобы не возникало лишних проблем с непрерывностью, выберем такую функцию истории, которая будет гарантировать непрерывность первой производной. Зададим историю в виде у(t) = а₀·t + у₀, где у₀ – количество особей в популяции в нулевой момент времени. Реализовать такую динамику на практике вполне возможно (по крайней мере, в искусственных экосистемах): начиная с момента времени t = -τ₂ и до момента t = 0 экспериментатор должен вводить все новые и новые особи в экосистему со скоростью^[54]

${а}_0= (b_0·{у}_0-d_0·{у}_0-{с}_0)/(1+ b_0·{\tau }_2).$

Решение уравнения Goudriaan-Roermund с отловом приведено на рис. 1 (оно получено при помощи MATLAB-функции dde23 – соответствующую программу см. в Приложении – но уравнение это столь простое, что его можно было решить и аналитически). Как видим, согласно полученному решению, после полного уничтожения популяции она немедленно возрождается… из ничего! И в дальнейшем растет довольно бодро.

Кстати, справедливости ради отметим, что, на первый взгляд, как раз для многих рыб описанная ситуация, может показаться имеющей биологический смысл: они отложили икру, их самих всех выловили (т.е. плотность популяции рыб упала до нуля), но из икринок вылупились новые рыбки-мальки и популяция возродилась. Хорошо,  пусть это будут не рыбки, а мышки, например – думаем, основная идея нашей критики вполне ясна. И заодно напомним, что, обосновывая уравнение (5), частным случаем которого (при K ® ¥) является уравнение Goudriaan-Roermund, Мэрди [1979, c. 113] взывал вовсе не к рыбкам, а говорил о длительных периодах созревания у млекопитающих. Да, собственно говоря, и рыбки-то мечут икру не каждый день, а ведь решение уравнения (как это ясно видно из рисунка и как это можно строго показать аналитически) начинает возрастать от нуля сразу же после прекращения отлова^[55]. В общем, рассматриваемое уравнение если и пригодно для описания какой-либо популяции, то только лишь для популяции мифических фениксов.

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

 

Рис. 1. Решение уравнения Goudriaan-Roermund с отловом. Значения коэффициентов в данном примере какого-то особого смысла не имеют, они приняты чисто условно – лишь для иллюстрации: b0=1, τ2=1, d0=0.1, с0=1.1; у0=1.

 

Уже завершая обсуждение использования ДУСОтА, сразу после рассмотрения уравнения (5), Мэрди [1979, с. 113] отмечает: «приведенные выше соображения относятся к видам с дискретными поколениями (например, в случаях, когда размножение происходит только в данное время года)». Но если размножение происходит только в какое-то определенное время года, то член уравнения, ответственный за размножение, т.е. b₀·[1 - N(t-τ₁)/K]·N(t-τ₂), должен работать только в это определенное время. Очевидно, что это не так: указанный член вычисляется в каждый момент времени и в каждый момент времени вносит вклад в скорость роста популяции.

Откуда могла возникнуть ошибочная интерпретация?

Все биологические системы являются сложными многокомпонентными, пространственно-структурированными, их элементы обладают индивидуальностью. При моделировании таких систем возможны два подхода. Первый – агрегированный, феноменологический. В соответствии с этим подходом выделяются определяющие характеристики системы (например, общая численность видов). Такой подход является исторически наиболее древним и свойственен динамической теории популяций. Другой подход – подробное рассмотрение элементов системы и их взаимодействий, построение имитационной модели, параметры которой имеют ясный физический и биологический смысл^[56]. Такая модель не допускает аналитического исследования, но при хорошей экспериментальной изученности фрагментов системы может дать количественный прогноз ее поведения при различных внешних воздействиях [Ризниченко, 2016, с. 10].

В теоретической биофизике возникли свои приемы и принципы. Согласно «принципу простоты» биологическая система (и, следовательно, описывающая ее математическая модель) должна быть сконструирована максимально просто (при условии выполнения заданной функции) [Романовский и др., 1975, с. 11-12]. Вот тут-то и кроется опасность! Мы здесь не будем оспаривать то, что биологическая система должна быть сконструирована максимально просто (скажем только, что не все авторы данной статьи разделяют это беспочвенное утверждение). Но даже если она устроена именно так, то эта «максимальная простота», может оказаться весьма сложной, чтобы действительно выполнять заданные (подчас чрезвычайно сложные) функции. Поэтому адекватная модель, если ее аккуратно построить, также окажется весьма сложной (оставаясь «максимально простой» – любое дальнейшее упрощение приведет к утрате каких-то функций). Однако велико искушение еще больше упростить модель, в частности, перейти от подробного рассмотрения элементов системы и их взаимодействий (т.е. от модели, имеющей ясный физический и биологиеский смысл) к феноменологическому описанию агрегированных переменных, которое, как правило, оказывается математически проще, но часто приводит к отсутствию физического и биологиеского смысла поведения модели в ряде ситуаций. В качестве первого примера рассмотрим обсуждение Романовским и др. [1975, с. 160-163] модели Ноака (6). Сначала указанные авторы рассматривают физически осмысленную картину биологической инерционности. Она такова.

Удельная скорость роста биомассы (μ) определяется, кроме концентрации субстрата, многими различными факторами. Во-первых, концентрациями ферментов и скоростями ферментативных реакций в «узком месте» (наиболее медленном звене) катаболической цепи. Во-вторых, скоростью синтеза мРНК на матрице ДНК и, наконец, скоростью синтеза белков и нуклеотидов из мономеров – аминокислот и мононуклеотидов. В этом последнем звене существенную роль играют рибосомы, на которых и происходит сборка готовых молекул белка. Из всех перечисленных звеньев цепи переработки субстрата в клеточную биомассу самым инерционным должно быть то, которое связано с наиболее крупными органеллами – рибосомами. Опыты показываю, что рибосомальный аппарат клетки в стационарных условиях достаточно загружен и при повышении уровня питания начинается интенсивный синтез новых рибосом. Поэтому при изменении концентрации субстрата скорость роста клеток изменяется не сразу, а лишь спустя некоторое время [Романовский и др., 1975, с. 160].

На этой основе строится математическая модель, которая, как и должно быть в соответствии с законами биохимической кинетики, оказывается системой обыкновенных дифференциальных уравнений, а не ДУСОтА!. Но данная модель, содержащая 3 уравнения (для концентраций биомассы, питательного субстрата и рибосом), способна описать конечное время адаптации клеток к новому уровню питания [Романовский и др., 1975, с. 160-162], хотя это время нигде явно не вводится. В модели же Ноака только 2 уравнения – для концентраций биомассы и субстрата, а вместо уравнения для концентрации рибосом вводится в явном виде время запаздывания. Таким образом, формально, система упростилась: мы получили два уравнения вместо трех, но биологический смысл в некоторой степени редуцировался. Вместо понятного процесса синтеза рибосом, зависящего от текущей концентрации субстрата, но определяющего концентрацию рибосом в будущем и соответственно, скорость синтеза биомассы тоже в будущем, а не сейчас, в модели Ноака мы имеем просто «данный свыше закон»: удельная скорость роста почему-то определяется концентрацией субстрата в прошлом. А если мы сделаем следующий эксперимент: в некоторый момент времени t₁ увеличим в среде, на которой растут микроорганизмы, концентрацию субстрата (до очень большой величины), а потом в момент времени (t₁+0.5·τ) уменьшим ее до нуля^[57]?  Вот тут-то и выяснится принципиальная разница между моделью основанной на процессах (она предсказывает, что клетки не будут расти, как это и есть в действительности – по крайней мере, они не могут расти с той интенсивностью, как росли в богатой питательной среде) и горе-моделью Ноака. Согласно последней, до момента (t₁+1.5·τ) микробы будут расти весьма активно, причем (в дистиллированной воде – вообще без питательного субстрата!!!) в интервале времени от (t₁+0.5·τ) до (t₁+1.5·τ) они будут расти интенсивнее, чем росли на питательной среде до момента времени t₁!

В качестве второго примера возьмем математические модели иммунологии. Выше мы отмечали, что большой вклад в построение математических моделей иммунологии (именно в форме ДУСОтА) внес Гурий Иванович Марчук. Но сам акад. Марчук [1991, с. 164-165] ясно указывал: «Запаздывание… обусловлено тем, что… чтобы соответствующие иммунокомпетентные клетки стали производить антитела… им необходимо пройти несколько стадий деления, на что затрачивается время τ… Но мы можем рассмотреть процесс прохождения клетками различных стадий… Тогда процесс образования плазматических клеток можно представить системой… без запаздывания,… модель представляет собой систему обыкновенных дифференциальных уравнений».

Вернемся теперь к моделям популяционной динамики и рассмотрим третий пример. Вспомним сказанное относительно уравнения Хатчинсона Рубиным и др. [1987, с. 252]: что в реальных популяциях интенсивность таких процессов, как отравление среды продуктами метаболизма зависит от численности взрослых особей, т.е. отрицательное влияние на коэффициент естественного прироста оказывают особи предыдущего поколения.

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

$dN\left(t\right)/dt=N\left(t\right)·[b_0-\alpha ·N\left(t\right)-\underset{0}{\overset{t}{\int }}N\left(s\right)·f\left(t-s\right)ds].$

Физический смысл функции f очевиден – она определяет степень влияния предыстории на динамику сообщества в звисимости от отдаленности этой предыстории во времени от рассматриваемого момента [Полуэктов и др., 1980, с. 156]. Похожее уравнение

$dN\left(t\right)/dt=b_0·N\left(t\right)·[1-K^{-1}·\underset{-\yen }{\overset{t}{\int }}N\left(s\right)·w(t-s)ds],$ (7)

(где w – это весовой коэффициент, который показывает, насколько численность популяции в предшествующие годы важна для определения доступности ресурсов сегодня) рассматривает Мюррей [2009, с. 44], а Братусь и др. [2010, c. 204] записывают его в виде

$dN\left(t\right)/dt=b_0·N\left(t\right)·[1-K^{-1}·\underset{-\yen }{\overset{0}{\int }}N(t+s)·w(-s)ds].$

Как видим, и в случае самоотравления популяции продуктами метаболизма, и в случае исчерпания доступных ресурсов в среде мы имеем не ДУСОтА, а интегро-дифференциальные уравнения. И это более правильно.

Действительно, например, в последнем случае эффект запаздывания должен был бы являться усреднением по численности за предшествующие годы и выражаться, таким образом, интегро-дифференциальным уравнением [Мюррей, 2009, с. 44], поскольку среднее значение функции представляет собой интеграл от нее, нормированный некоторым образом. Совершенно то же самое можно сказать и по поводу самоотравления популяции: необходимо провести усреднение по численности за предшествующие годы, поскольку все особи вносили вклад в выработку продуков метаболизма и, следовательно, в самоотравление. В каком же случае эти уравнения превратятся в уравнение Хатчинсона?

 

Рис. 2. Типичная весовая функция w(t) для обобщенного эффекта запаздывания в модели (7) [Мюррей, 2009, с. 45].

 

На практике w(t) обычно стремится к нулю при больших положительных и отрицательных t, а максимума достигает при некотором значении времени Т. В типичном случае w(t) имеет вид, представленный на рис. 2. Если w(t) дает более острый пик в том смысле, что область вокруг Т более узкая, то тогда в предельном случае мы можем считать, что w(t) аппроксимируется функцией Дирака δ(t - Т), а интегро-дифференциальное уравнение Мюррея в этом случае упрощается до уравнения Хатчинсона, поскольку [Мюррей, 2009, с. 44]

$\underset{-\yen }{\overset{\yen }{\int }}N\left(s\right)·\delta (t - Т - s)ds=N(t - Т).$

Совершенно аналогично в уравнение Хатчинсона вырождается и уравнение Вольтерра – если в качестве весовой функции f выступит δ-функция Дирака. По биологическому смыслу использование δ-функции в качестве весовой функции для w и f означает, что только особи строго одного возраста могли, соответственно, потреблять доступный ресурс или отравлять среду. Очевидно, что это – чепуха!

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

Наконец, очень кратко обсудим ДУСОтА, при помощи которых моделируется то, что размножается не вся популяция, а только особи определенного возраста – ур. типа (3) и (5) при условии τ₂ > 0. Вообще говоря, аккуратное описание динамики популяции, имеющей непрерывное распределение по возрасту, приводит к уравнению в частных производных (см., например, [Романовский и др., 1975, с. 171; Корзухин и Семевский, 1992; Кузнецов и др., 2005, с. 55-56]), в котором одна производная берется по физическому времени, а вторая – по возрасту. Кузнецов и др. [2005, с. 57-59] показали, что это уравнение можно свести к ДУСОтА только в двух частных случаях, которые довольно далеки от законов смертности, прироста и внутривидовой борьбы в реальной природе. Однако авторы тут же пишут: «Полезность проведенных выкладок состоит в том, что они дают подсказку как можно решать проблему устойчивости… Разобранная здесь возможность такого сведéния дает нам право опереться на мощные результаты, полученные для уравнений с запаздывающим аргументом». Вот, оказывается, в чем дело: поскольку аналитическая теория для ДУСОтА разработана математиками довольно хорошо, есть искушение свести более сложные (но и более реалистичные!) уравнения к ним, чтобы была возможность применить эту теорию. Причем свести любой ценой – даже ценой потери биологического смысла. В общем, как это часто бывает, «ищут там где светлее, а не там, где потеряли». Более того, перейдя к уравнениям, содержащим не слишком много биологического смысла, пытливые исследователи довольно быстро обнаруживают, что и эти уравнения оказываются слишком сложны, чтобы к ним возможно было применить «мощные результаты, полученные для уравнений с запаздывающим аргументом». Поэтому, на той же самой странице, где Кузнецов и др. [2005, с. 59] уповали на «мощные результаты», через 3 абзаца, видимо осознав сложность проблемы, они пишут: «Правда, простого экспоненциального полинома можно и не получить… На этом мы не будем останавливаться, …сложные случаи исследуем численными методами». И слава Богу! Ибо численными методами (в качестве которых избираются неявные схемы бегущего счета [Кузнецов и др., 2005, с. 60]) в данном случае исследуется классическое уравнение динамики популяции, имеющей непрерывное распределение по возрасту, т.е. уравнение в частных производных, а не выводимые из него для частных случаев ДУСОтА.

Итак, во многих случаях запаздывание вводится как временная характеристика второстепенных или малоизученных процессов, которые на данном этапе построения модели в нее не включаются. Это может быть, например, время транспорта молекул от места их синтеза к месту их включения в систему реакций; время формирования клеток определенного типа, участвующих в имунной реакции; длительность реакции части популяции на лимитирующие факторы окружающей среды и т.д. [Бабский и Мышкис, 1983, с. 383]. Таким образом, использование ДУСОтА относится, в основном, к первому (феноменологическому) из двух указанных в начале этого раздела модельных подходов. Поэтому модели в виде дифференциально-разностных уравнений не стоит абсолютизировать^[58].

Тем не менее, мы не предлагаем совершенно отказаться от математических моделей в виде ДУСОтА, а призываем лишь к осторожности в их использовании – использовать их можно, но только до той поры, пока не возникает какая-то биологически бессмысленная ситуация. А раз использовать их, все-таки, можно, то на первый план выходит вопрос о практической реализации данных моделей – о математическом аппарате и программном обеспечении для их исследования.

Как практически работать с дифференциально-разностными моделями?

Хотя дифференциальные уравнения с отклоняющимся аргументом появились в литературе еще в XVIII в., но до недавнего времени не были сформулированы основные теоремы общей теории ДУСОтА, более того, в литературе не было даже четкой постановки начальной задачи. Это впервые было сделано в 1949-1950 гг. в диссертации А.Д. Мышкиса^[59] [Норкин, 1965, c. 8]. В романтический период оформления новой области математического анализа, именуемой «Дифференциальные уравнения с отклоняющимся аргументом», почти каждое слово, сказанное тогда в этой области оказывалось новым. Как изменилось положение с тех пор [Мышкис, 1972, c. 6]!

Хотя теория дифференциальных уравнений с отклоняющимся аргументом принадлежит к числу сравнительно молодых и бурно развивающихся разделов теории ОДУ ⁶⁰, но теории ДУСОтА и ее многообразным приложениям уже посвящено большое количество журнальных статей [Норкин, 1965, c. 5] (добрая сотня работ, представляющих математический интерес, появляется ежегодно в этой области). Среди них имеется ряд обзоров^[61] [Мышкис, 1972, c. 6] и монографий, целиком или частично посвященных различным аспектам данной теории [Норкин, 1965, c. 5]. Укажем прежде всего монографии [Эльсгольц и Норкин, 1971; Мышкис, 1972; Солодов и Солодова, 1980; Хейл, 1984]. Очень простой и понятный анализ линейного уравнения дан в [Хайрер и др., 1990, с. 309-311]. Особо отметим [Хэссард и др., 1985, с. 126-153], где рассматривается приложение теории бифуркации рождения цикла к уравнениям с запаздыванием, что имеет важное значение для математических моделей биологии, поскольку колебательные режимы в биологических системах встречаются весьма часто. Надеемся, заинтересованный читатель вполне сможет освоить аналитические основы теории уравнений с запаздыванием по вышеприведенной литературе.

Усилиями математиков разных стран теория ДУСОтА (в основном, запаздывающего и в несколько меньшей степени нейтрального типов) глубоко разработана в различных направлениях, найдены естественные постановки задач, установилась адэкватная терминология [Мышкис, 1972, c. 6]. Однако аналитическое исследование систем с запаздыванием в отличие от одиночных уравнений с запаздывающим аргументом развито хуже, и общих теоретических результатов здесь получено немного^[62] [Марри, 1983, с. 186]. В силу сложности математических моделей реальных биологических процессов, чаще всего они представляют собой именно системы уравнений, так что аналитическое исследование их если и возможно, то сопряжено с огромными трудностями и потому требует очень большого времени. В связи с этим далее мы будем говорить лишь о численном подходе к решению ДУСОтА.

В настоящее время теория численного интегрирования ДУСОтА развита довольно хорошо (для знакомства с ней см., например, [Холл и Уатт, 1979, с. 256-268; Хайрер и др., 1990, с. 304-308; Bellen and Zennaro, 2013]), однако первые компьютерные программы в этой области по ряду причин были неудовлетворительными^[63]. Поэтому в настоящее время не следует брать наугад случайно подвернувшуюся под руку программу. Впрочем, из уважения к памяти и интеллекту предшественников упомянуть их стоит (хотя вряд ли сейчас надежно и однозначно удастся решить вопрос о приоритете, ибо с момента создания первых программа для численного решения ДУСОтА прошло несколько десятков лет).

В СССР о машинных экспериментах на БЭСМ-6 (с использованием языка Алгол) над системой нелинейных дифференциально-разностных уравнений в 1974 г. сообщил Б.Н. Хидиров (в 71-ом выпуске издававшихся в Ташкенте «Вопросов кибернетики»: с. 129-136). Поскольку эта система была весьма сложной (содержала 14 уравнений) [Хидиров, 2104, с. 9-10], очевидно, что она не могла быть проинтегрирована аналитически. Следовательно, речь идет о какой-то программе приближенного (скорее всего – численного) интегрирования. К сожалению, больше о ней авторам ничего не известно.

Но даже если указанная программа была действительно первой, она, оставшись по эту сторону «железного занавеса», к сожалению, не оказала какого-либо влияния на развитие мировой науки. А в свободном мире, по мнению Л.Ф. Шампайна и др. [2009, с. 248], первой процедурой численного решения дифференциальных уравнений с запаздыванием была DMRODE^[64], в которой для получения аппроксимаций решения на интервале интегрирования вне точек сетки применялась кубическая эрмитова интерполяция. Во-первых, этот подход не может быть признан удовлетворительным, поскольку размер шага интегрирования, выбираемый с учетом необходимости обеспечения заданной точности интегрирования, может оказаться недостаточно малым для точной интерполяции.

Во-вторых, многие ранние версии численных процедур решения ДУСОтА были написаны на языке Fortran^[65], в реализации которого отсутствовала возможность динамического выделения компьютерной памяти. Это существенно усложняло интерфейс с пользователем и нередко приводило к тому, что выделенного объема оказывалось недостаточно для полного решения задачи. Обоих указанных недостатков лишена современная MATLAB-функция dde23. В частности, динамическое выделение компьютерной памяти, предусмотренное в MATLAB при использовании структур, позволило упростить и усовершенствовать интерфейс пользователя численной процедуры dde23 [Шампайн и др., 2009, с. 249].

Численное решение ДУСОтА в системе MATLAB при помощи стандартной функции dde23 описано, например, в [Ануфриев и др., 2005, с. 323-330], кроме того, пример использования dde23 (для решения системы Лотки-Вольтерра с запаздывающим аргументом) приведен еще и в [Дьяконов, 2008, с. 429-430]^[66]. Но особенно подробно тонкости решения задач для ДУСОтА различной сложности рассмотрены в [Шампайн и др., 2009, с. 252-281]^[67]. Здесь мы не будем вдаваться в многочисленные математические подробности, но не можем не отметить одно важное обстоятельство.

Численные методы решения задач с начальными условиями для дифференциальных уравнений с запаздывающим аргументом предполагают, что несколько первых производных соответствующих решений непрерывны. К сожалению, нарушения непрерывности производных малого порядка нередко имеют место, т.к. первая производная функции истории почти всегда отличается от первой производной решения в начальной точке [Шампайн и др., 2009, с. 244]. Почему же авторы моделей не озадачиваются построением такой функции истории, которая обеспечивала бы непрерывность решения? По-видимому, пока не предложено общепринятого естественного способа построить ее таким образом, чтобы, с одной стороны, она имела бы биологический смысл, а, с другой – производные в начальной точке были бы непрерывны (чаще всего на практике в качестве функции истории берется некоторая разумная константа, тогда как решение уравнения константой не является; в результате, при t®-а lim dу/dt = lim dconst/dt = 0, но при t®+а lim dу/dt ¹ 0).

БЛАГОДАРНОСТИ

Авторы выражают признательность одному из анонимных рецензентов, замечания которого позволили существенно улучшить статью, а часть текста его рецензии была использована нами в конце разд. «Запаздывание в моделях межпопуляционных взаимодействий» и в подстрочной сноске № 27. Другой рецензент указал нам на то, что в первоначальном варианте статьи получили недостаточное отражение работы к.ф.-м.н. Ю.С. Колесова, а публикации д.ф.-м.н. С.А. Кащенко вообще не были упомянуты. Мы благодарны за указание на эти недопустимые упущения.

About the authors

Mikhail V. Glagolev

Lomonosov Moscow State University; Yugra State University; Tomsk State University; Institute of Forest Science of the Russian Academy of Sciences; Water Problems Institute of the Russian Academy of Sciences

Author for correspondence.
Email: m_glagolev@mail.com
Russian Federation, 1, Leninskie gory, Moscow, 119991 ; 16, Chehova street, Khanty-Mansiysk, 628012 ; 34a, Lenina prospect, Tomsk, 634050 ; 21, Sovetskaya street, Uspenskoe village, Moscow region, 143030 ; 3, Gubkina street, Moscow 119333

Aleksandr F. Sabrekov

Yugra State University; Tomsk State University; Institute of Forest Science of the Russian Academy of Sciences; Water Problems Institute of the Russian Academy of Sciences

Email: sabrekovaf@gmail.com
16, Chehova street, Khanty-Mansiysk, 628012 ; 34a, Lenina prospect, Tomsk, 634050 ; 21, Sovetskaya street, Uspenskoe village, Moscow region, 143030 ; 3, Gubkina street, Moscow 119333

Vladimir M. Goncharov

Moscow State University

Email: m_glagolev@mail.com
1, Leninskie gory, Moscow, 119991

References

Supplementary files