Модификация модели Базыкина-Алексеева системы "хищник-жертва" с насыщением

Утюпин Ю.В.
Мирнинский политехнический институт

Аннотация:

В работе рассматривается математическая модель системы хищник-жертва, в которой учитывается действие эффекта насыщения в процессе взаимодействия видов. В отличие от известной модели Алексеева - Базыкина [1,2], в рассматриваемой модели учитывается зависимость скорости прироста численности хищников от условий их питания в течение определенного времени. Для системы трех интегро-дифференциальных уравнений изучаются динамические режимы, которые реализуются при различных значениях параметров.

In this article the mathematical model of a system a predator - prey is considered, in which the operation of effect of a saturation during interaction of species is taken into account. As against known model of Alekseev - Bazykin [1,2] in considered model the association of a velocity of an increase of number of predators on conditions of their feed during the defined time is taken into account. For a system of three integro-differential equations the dynamic conditions are studied which are realized at various values of parameters.

Ключевые слова: модель системы хищник - жертва.

1. Анализу математических моделей системы хищник - жертва посвящено огромное количество работ (см., например, [1-7]). Практически во всех работах рассматриваются различные модификации модели Лотки - Вольтерра [5], учитывающие в ней действие дополнительных регуляторных механизмов. В частности, в [1-4] рассматривались модификации модели Лотки - Вольтерра с "эффектом насыщения" хищников: с ростом численности жертв скорость их потребления хищниками и скорость прироста числа хищников стремились к некоторой конечной величине.

В модели Лотки - Вольтерра [5]

$\begin{displaymath} \frac{dx}{dt}=\alpha x-\gamma xz-\beta x^2, \quad \frac{dz}{dt}=-\alpha_1 z+\gamma_1 xz-\beta_1 z^2, \eqno(1) \end{displaymath}$

где

- численность жертв,

- численность хищников в момент времени

, предполагается, что все параметры постоянны и неотрицательны, $\alpha ,\alpha_1 ,\beta ,\beta_1 ,\gamma ,\gamma_1 \equiv const\geq 0$ . Это означает, что в (1) при фиксированном значении

рост

вызывает неограниченный рост скоростей потребления жертв и прироста численности хищников.

В модификациях [1-4] модели (1) предполагалось, что $\gamma$ и $\gamma_1$ не являются постоянными величинами и зависят от значений :

$\begin{displaymath} \gamma=\frac{a}{1+bx}, \quad \gamma_1=\frac{a_1}{1+bx}, \end{displaymath}$

где

. Следуя [1-4], рассмотрим переменную $\theta$ , соответствующую количеству пищи в пищеварительной системе (ПС) одного хищника. Пусть также

- объем ПС одного хищника,

. Будем предполагать, что $\gamma$ - функция от $\theta$ : $\gamma=\gamma_0 (n-\theta)$ , - и чем больше наполненность ПС, тем меньше скорость уничтожения жертв. Таким образом, первое уравнение системы (1) примет вид:

$\begin{displaymath} \frac{dx}{dt}=\alpha x-\beta x^2-\gamma_0 (n-\theta)xz.\eqno(2) \end{displaymath}$

Учитывая то, что скорость потребления жертв равна $\gamma_0 (n-\theta)xz$ , скорость наполнения ПС хищника равна $\gamma_0 (n-\theta)x$ . Таким образом, динамика переменной $\theta$ определяется выражением

$\begin{displaymath} \frac{d\theta}{dt}=-l\theta+\gamma_0 (n-\theta)x, \eqno(3) \end{displaymath}$

где $l=const>0, l\theta$ - скорость снижения наполненности ПС. В [2] предполагалось, что $\theta$ является быстрой переменной, что, после понижения порядка системы, позволило получить указанную выше модификацию модели (1).

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

$\begin{displaymath} \gamma_1\frac 1 \tau \int\limits_{t-\tau}^t \theta (s)ds. \end{displaymath}$

В наиболее простом случае эта зависимость является прямо пропорциональной. Таким образом получаем уравнение:

$\begin{displaymath} \frac{dz}{dt}=-\alpha_1 z+\gamma_1 z\frac 1 \tau \int\limits_{t-\tau}^t \theta (s)ds-\beta_1 z^2. \eqno(4) \end{displaymath}$

Для корректной постановки задачи требуется задание начальных данных для функции $\theta$ на отрезке времени $\left [-\tau,0\right]$ . Объединяя уравнения (2) - (4), получаем модель системы хищник - жертва с "эффектом насыщения" с начальными данными $x(0)=x_0>0, z(0)=z_0>0, \theta (t)=\phi (t)$ при $t\in \left [-\tau,0\right]$ , где $\phi (t)\in C_{\left [-\tau,0\right]}$ и $0\leq\phi (t)\leq n$ .

Далее, в пунктах 2-7 мы будем рассматривать модель без саморегуляции хищника (т.е. в уравнении (4) $\beta_1=0$ ), а затем полную модель.

2. Рассмотрим модель без саморегуляции в популяции хищников.

$\begin{displaymath} \frac{dx}{dt}=\alpha x-\beta x^2-\gamma_0 (n-\theta)xz, \end{displaymath}$

$\begin{displaymath} \frac{dz}{dt}=-\alpha_1 z+\gamma_1 z\frac 1 \tau \int\limits_{t-\tau}^t \theta (s)ds, \end{displaymath}$

$\begin{displaymath} \frac{d\theta}{dt}=-l\theta+\gamma_0 (n-\theta)x. \end{displaymath}$

С помощью замены переменных:

$\begin{displaymath} x=\frac l {\gamma_0} u,\quad z=\frac a {\gamma_0 n} v, \quad \theta=np, \end{displaymath}$

число параметров в модели уменьшается и она приобретает вид

$\begin{displaymath} \begin{array}{*{20}c} \frac{du}{dt}=\alpha u(1-bu-(1-p)v),... ...ds\right), \frac{dp}{dt}=l(-p+(1-p)u), \end{array} \eqno(5) \end{displaymath}$

где $b=\frac{\beta l}{\gamma_0\alpha}, \quad m=\frac{\gamma_1 n}{\alpha_1}$ , $u(0)=u_0>0, \quad v(0)=v_0>0, \quad p(t)=\rho (t)$ , где $\rho (t)\in C_{\left [-\tau,0\right]}$ и $0\leq\rho (t)\leq 1$ .

3. Учитывая то, что

$\begin{displaymath} \frac{du}{dt}\vert _{u=0}=0,\quad \frac{dv}{dt}\vert _{v=0}=0, \quad \frac{dp}{dt}\vert _{p=0}=lu\geq 0, \end{displaymath}$

траектории системы (5) не выходят за границы

Заметим, что $\frac{du}{dt}\leq 0$ при $u\geq\frac 1 b$ . Следовательно, все траектории асимптотически входят в полосу $0<u\leq\frac 1 b$ и не могут выйти за ее границы, если $u_0\leq\frac 1 b$ . Также выполнено неравенство $\frac{dp}{dt}\vert _{p=1}<0$ . Следовательно, $\forall t\quad p(t)<1$ . Если $0<u<\frac 1 b$ , то $\frac{dp}{dt}\vert _{p=\frac{1}{b+1}}<0$ . Таким образом, все траектории асимптотически входят в область $0<p<\frac 1 {b+1}$ .

Ограниченность доказывается более сложным образом (см. [10]), но, к сожалению, только ограниченность каждой траектории, но не ограниченность в совокупности.

4. Для доказательства существования и единственности решения задачи приведем (5) к стандартному виду. Для этого примем интеграл во втором уравнении за новую переменную , после чего с помощью замены переменной $t=t'\tau$ (5) приводится к виду:

$\begin{displaymath} \begin{array}{*{20}c} \frac{du}{dt'}=a u(1-bu-(1-p)v),\\ ... ...u),\\ \frac{dR}{dt'}=p(t')-p(t'-1), \end{array} \eqno(6) \end{displaymath}$

где $a=\alpha\tau$ , $c=\alpha_1\tau$ , $d=l\tau$ , с начальными данными $u(0)=u_0>0, v(0)=v_0>0, p(t')=\rho (t')$ , где $\rho (t')\in C_{\left [-1,0\right]}$ , $0\leq\rho (t')\leq 1$ , $R(0)=\int\limits_{-1}^0\rho (s)ds$ . В результате получаем систему (6) ОДУ с запаздывающим аргументом, для которой нетрудно показать существование и единственность решений [10].

5. При любых значениях параметров система (5) имеет две стационарные точки: и $B(\frac 1 b, 0, \frac 1 {b+1} )$ . Кроме того, при условии в существует нетривиальное равновесие $C(\frac{1}{m-1},\frac{m(m-1-b)}{(m-1)^2},\frac{1}{m} )$ .

Соответствующими стационарными точками для системы (6) будут , $B(\frac 1 b, 0, \frac 1 {b+1},\frac 1 {b+1} )$ и $C(\frac{1}{m-1},\frac{m(m-1-b)}{(m-1)^2},\frac{1}{m},\frac{1}{m})$ . Исследуем эти состояния равновесия на устойчивость.

Линеаризуя (6) вблизи точки , получаем характеристическое уравнение: $\lambda (\lambda-a) (\lambda+c) (\lambda+d)=0$ . Поскольку это уравнение имеет положительный корень, то, следовательно, точка неустойчива.

Характеристическое уравнение для точки имеет вид:

$\begin{displaymath} \lambda (\lambda+a)\left (\lambda-c\frac{m-b-1}{b+1}\right )\left (\lambda+d\frac{b+1}{b}\right )=0. \end{displaymath}$

Это уравнение имеет положительный корень при

. Отсюда следует, что если в

существует точка

, то

неустойчива.

Характеристическое уравнение для точки имеет вид:

$\begin{displaymath} \lambda^2\left[\left (\lambda+\frac{ab}{m-1}\right )\left ... ...da+\frac{dm}{m-1}\right )- \frac{ad(m-1-b)}{(m-1)^2}\right]+ \end{displaymath}$

$\begin{displaymath} (1-e^{-\lambda})\frac{acd(m-1-b)}{m-1}=0. \eqno(7) \end{displaymath}$

6. Пусть

$\begin{displaymath} A=\frac{ab+dm}{m-1}, \quad B=\frac{ad(bm-m+1+b)}{(m-1)^2}, \quad D=\frac{acd(m-1-b)}{m-1}. \end{displaymath}$

При этих обозначениях уравнение (7) примет вид:

$\begin{displaymath} \lambda^2 (\lambda^2+A\lambda+B)+D(1-e^{-\lambda} )=0. \eqno(8) \end{displaymath}$

Заметим,что

. Из (8) видно, что $\lambda=0$ является корнем. Этот корень возникает при переходе от системы трех уравнений к системе четырех уравнений. Система (6), вообще говоря, имеет бесконечное множество нетривиальных стационарных точек. Но система (6) эквивалентна (5) только при условии $R(0)=\int\limits_{-\tau}^{0}\rho (s)ds$ .

При выполнении этого условия имеется единственное нетривиальное равновесие С. В этом случае появившийся корень $\lambda=0$ не влияет на устойчивость.

Если все корни характеристического уравнения (8) (кроме одного нулевого) лежат слева от мнимой оси, то соответствующее решение является асимптотически устойчивым.

Справедлива следующая теорема [10]:

Теорема. Для того, чтобы все корни квазиполинома (8) (кроме $\lambda=0$ ) лежали слева от мнимой оси, необходимо и достаточно, чтобы выполнялись следующие два условия:

а) функция $G(y)=-\cos y (Ay^3)+\sin y (y^4-By^2+D)$ имеет на интервале $(0,\pi )$ два корня;

б) $1-\frac{D\sin y_1}{Ay^3_1}<0,1-\frac{D\sin y_k}{Ay^3_k}>0\quad \forall k>1$ , где - корни функции . $\Box$

Следствие 1. Если , то не все нули квазиполинома (8) лежат слева от мнимой оси.

Следствие 2. Если $B>\frac D 2\left (D\leq 24\right )$ или $B\geq\frac D 2\left (D>24\right )$ , то $\exists A_0:\forall A>A_0$ все нули квазиполинома (8) лежат слева от мнимой оси.

Следствие 3. Если одно из следующих условий не выполняется:

- функция имеет корень на $(0,\pi )$ ,

- второй положительный корень функции лежит в интервале $(2\pi k+\pi,2\pi k+2\pi )$ ,

- $D<\frac{B\pi^2-\pi^4}2$ , то при $A\to 0$ появляются корни квазиполинома (8), имеющие положительные действительные части.

Замечание. Выполнение всех трех условий последнего следствия не гарантирует отсутствия корней квазиполинома с положительной действительной частью при $A\to 0$ .

**Figure 1:**
$\begin{figure} \begin{center} \epsfxsize=80mm \epsfysize=70mm \epsfbox{ut2.eps} \end{center} \end{figure}$

Рассмотрим структуру плоскости параметров (рис. 1). Здесь $D=\phi (B)$ - кривая, ограничивающая область существования корней функции $\frac{y^4-By^2+D}D$ на интервале $(0,\frac\pi 2 )$ . Она состоит из части параболы $D=\frac{B^2}4$ и части прямой $D=B{\left (\frac\pi 2\right )}^2-{\left (\frac\pi 2\right )}^4$ .

Можно показать [10], что:

1. Если точка (B,D) лежит в области , ограниченной прямой , линией касания функций $\cos y$ и $\frac{y^4-By^2+D}D$ (пунктирная линия) и прямой $D=\frac{B\pi^2-\pi^4}2$ , то при любом не все нули квазиполинома (8) лежат слева от мнимой оси.

2. Если точка принадлежит областям (рис. 2), то при $A\to\infty$ все нули квазиполинома будут лежать слева от мнимой оси. Если не принадлежит областям , то при $A\to 0$ не все нули квазиполинома будут лежать слева от мнимой оси. Если лежит в области V (заштрихованная область), то при $A\to\infty$ и при $A\to 0$ не все корни квазиполинома лежат слева от мнимой оси.

3. Если точка принадлежит области и

$\begin{displaymath} \begin{array}{*{20}c} A<2\frac {2D-B\left (\frac\pi 2\righ... ...ight )}}{\left (\frac \pi 2\right )^4}, \end{array} \eqno(9) \end{displaymath}$

то не все нули квазиполинома (8) лежат слева от мнимой оси.

4. Если точка принадлежит области и

$\begin{displaymath} A<2\frac {2D-B\left (\frac\pi 2\right )^2-\sqrt{\left (B \... ...ht )^4\left (B^2-4D\right )}}{\left (\frac \pi 2\right )^4}, \end{displaymath}$

то не все нули квазиполинома (8) лежат слева от мнимой оси.

7. Анализ зависимости устойчивости нетривиального стационарного состояния от $\tau$ методом D-разбиений показал [9], что увеличение запаздывания в (5) сужает область устойчивости.

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

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

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

$\begin{displaymath} \begin{array}{*{20}c} \frac{du}{dt}=\alpha u(1-bu-(1-p)v),... ...),\\ \par\frac{dp}{dt}=l(-p+(1-p)u), \end{array} \eqno(10) \end{displaymath}$

где $b=\frac{\beta_l}{\gamma_0\alpha},\quad b_1=\frac{\beta_1\alpha}{\gamma_0\alpha_1 n} \quad m=\frac{\gamma_1 n}{\alpha_1}$ , $u(0)=u_0>0, \quad v(0)=v_0>0, \quad p(t)=\rho (t)$ , где $\rho (t)\in C_{\left [-\tau,0\right]}$ и $0\leq\rho (t)\leq 1$ .

9. Также как и в п. 3 доказывается, что траектории не выходят за границы и все траектории асимптотически входят в область $(0<u\leq\frac 1 b)\times(0<p<\frac 1 {b+1})$ и не могут выйти за ее границы. Кроме того, заметим, что $\frac{dv}{dt}\leq 0$ при $v\geq\frac {1}{b_1}(\frac m{b+1}-1)$ и при $p<\frac 1 {b+1}$ на отрезке $(t-\tau,t)$ . Следовательно, в существует инвариантное компактное притягивающее множество $[0,\frac 1 b]\times\left[0,\frac 1 {b_1}(\frac m{b+1}-1)\right]\times[0,\frac 1 {b+1}]$ .

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

10. Перейдем к исследованию нетривиальных состояний равновесия. Состояния равновесия удовлетворяют следующим уравнениям:

$\begin{displaymath} v^{3}-v^{2}\frac{2}{b_1}\left (m-1\right )+v\frac{1}{b_1}\l... ...ight )^2}{b_1}\right )-\frac m {b_1^2}\left (m-1-b\right )=0, \end{displaymath}$

$\begin{displaymath} p=\frac{1+b_1v} m, \end{displaymath}$

$\begin{displaymath} u=\frac {1-\left (1-p\right )v} b. \end{displaymath}$

Если выполняются условия

$\begin{displaymath} \left (m-1\right )^{2}-3mb_{1}\left (1+b\right )>0 \end{displaymath}$

$\begin{displaymath} \left (\left (m-1\right )^{3}+\frac{27}{2}mb_{1}\left (bm-... ...left (m-1\right )^{2}-3mb_{1}\left (1+b\right )\right )^{3}, \end{displaymath}$

то в

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

$\begin{displaymath} \lambda(\lambda^3+A\lambda^2+B\lambda+C)+(1-e^{-\lambda})D=0,\eqno(11) \end{displaymath}$

где

$\begin{displaymath} A=abu+\alpha_1 b_1 (1-bu)(1+u)+d(1+u), \end{displaymath}$

$\begin{displaymath} B=abu\alpha_1 b_1(1-bu)(1+u)+a b u d(1+u)+\alpha_1 b_1 d (1-bu)(1+u)^2-adu(1-bu), \end{displaymath}$

$\begin{displaymath}C=abud\alpha_1 b_1 (1-bu)(1+u)^2-aud\alpha_1 b_1(1-bu)^2 (1+u), \end{displaymath}$

$\begin{displaymath} D=ad\alpha_1 m \frac{(1-bu)u}{1+u}. \end{displaymath}$

Заметим, что в случае, когда система имеет три стационарные точки (обозначим их за

по порядку возрастания

) точка

всегда является неустойчивой, поскольку в ней

. В проекции фазового пространтва на плоскость OUV эта точка ведет себя подобно седловой точке.

Дальнейший анализ устойчивости сводится к исследованию расположения нулей квазиполинома (11). Для этого мы, как и в пункте 7, воспользовались методом D-разбиений.

Области в пространстве параметров с различным количеством корней квазиполинома (11), имеющих положительные вещественные части, отделяются параметрическими гиперповерхностями

$\begin{displaymath} A=\frac C {y^2} + \frac{D \sin \tau y}{\tau y^3}, \end{displaymath}$

$\begin{displaymath} B=\frac {\tau y^4+D(1-\cos \tau y)}{\tau y^2}. \end{displaymath}$

Из-за большого числа параметров картина D-разбиений получается достаточно сложной для того, чтобы ее приводить здесь. Отметим лишь, что при различных соотношениях параметров C и D характер влияния запаздывания $\tau$ на устойчивость также различен. Увеличение запаздывания может оказывать на стационарную точку как стабилизирующее, так и дестабилизирующее воздействие (в отличие от модели без саморегуляции у хищников, где запаздывание является дестабилизирующим фактором).

С помощью численных экспериментов были выявлены следующие режимы:

- вырождение популяции хищника;

- стабилизация на ненулевом уровне;

- одна неустойчивая стационарная точка, окруженная глобально устойчивым пределным циклом;

- три нетривиальных состояния равновесия, два из которых устойчивы;

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

- три нетривиальных состояния равновесия, одно из которых неустойчиво; все траектории наматываются на один аттрактор - устойчивое состояние равновесия.

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

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

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

Литература

1: Базыкин А.Д. Система Вольтерра и уравнение Михаэлиса - Ментен.// Вопросы математической генетики, ИЦИГ СО АН СССР, Новосибирск, 1974, 103-143.
2: Алексеев В.В. Влияние фактора насыщения на динамику системы "хищник - жертва"// Биофизика, 18, 15, 1973, 922-926.
3: Алексеев В.В. Биофизика сообществ живых организмов.//Успехи физ. наук, 120, 4, 1976, 647-676.
4: Базыкин A.Д. Математическая биофизика взаимодействующих популяций. Наука ,Москва, 1985.
5: Вольтерра В. Математическая теория борьбы за существование. Наука, Москва, 1976.
6: Полуэктов Р.А., Пых Ю.А., Швытов И.А. Динамические модели экологических систем. Гидрометеоиздат, Ленинград, 1980.
7: Свирежев Ю.М., Логофет Д.О. Устойчивость биологических сообществ. Наука, Москва, 1978.
9: Недорезов Л.В., Утюпин Ю.В., Утюпина С.П. Эффект насыщения в модели системы "хищник - жертва".// Сибирский журнал индустриальной математики, IV, 1, 2001, 150-164.

Ваши комментарии

[Головная страница]
[Конференции]
[СО РАН]

© 2001, Сибирское отделение Российской академии наук, Новосибирск
© 2001, Объединенный институт информатики СО РАН, Новосибирск
© 2001, Институт вычислительных технологий СО РАН, Новосибирск
© 2001, Институт систем информатики СО РАН, Новосибирск
© 2001, Институт математики СО РАН, Новосибирск
© 2001, Институт цитологии и генетики СО РАН, Новосибирск
© 2001, Институт вычислительной математики и математической геофизики СО РАН, Новосибирск
© 2001, Новосибирский государственный университет
Дата последней модификации Wednesday, 05-Sep-2001 10:45:36 NOVST