On the interaction of species capable of explosive growth

One way of solving the logistic population equation (here x = x(t)) (1.1) x′=ax−bx2 x' = ax - b{x^2} is to divide this equation by x², and obtain a linear equation for u=1x u = {1 \over x} . Here a > 0 is the growth rate, and b > 0 is the limitation (or self-limitation) coefficient, both given numbers. We wish to explore the interactions of two species with populations x = x(t) and y = y(t) for which the substitution u=1x u = {1 \over x} and v=1y v = {1 \over y} leads to a linear system. The model we consider is (1.2) x′=ax+x2(by+e)y′=dy+y2(cx+f) , \matrix{ {x' = ax + {x^2}\left( {{b \over y} + e} \right)} \cr {y' = dy + {y^2}\left( {{c \over x} + f} \right){\kern 1pt} ,} \cr } with constants a,b,c,d,e and f. Dividing the first equation by x², the second one by y², and setting u=1x u = {1 \over x} and v=1y v = {1 \over y} , gives a linear system (1.3) −u′=au+bv+e−v′=cu+dv+f . \matrix{ { - u' = au + bv + e} \cr { - v' = cu + dv + f{\kern 1pt} .} \cr } The signs of the coefficients determine the type of interaction, which will include both predator-prey and competing species cases.

Let us compare (1.2) with the classical Lotka-Volterra predator-prey model (1.4) x′=x(a−b y)y′=y(−c+d x) , \matrix{ {x' = x\left( {a - b{\kern 1pt} y} \right)} \cr {y' = y\left( { - c + d{\kern 1pt} x} \right){\kern 1pt} ,} \cr } where the constants a,b,c,d are positive. In (1.4) the species affect each other through the growth rate: the prey, with the number given by x(t), improves the growth rate of the predator, with the number y(t), while the predator decreases the growth rate of the prey. In the model (1.2) the species affect each other through their limitation coefficients. This appears to be more realistic: the presence of foxes is not likely to decrease the fertility of rabbits (new rabbits will be born at the same rate), but will place a limitation on the growth of rabbit population.

Similarly to the Lotka-Volterra model, the proposed model (1.2) predicts oscillatory behavior for predator-prey interaction, and either stable coexistence or competitive exclusion for competing species. Unlike the Lotka-Volterra model, it is possible that the population number of one of the species goes to infinity in finite time, while the number of the other species remains finite and positive. Explosive growth of populations occurs often in nature, in connection with various “outbreaks”, see e.g., D. Ludwig, D.D. Jones and C.S. Holling [4] for a model of spruce budworm. Notice that our analysis leads to some non-standard questions about linear systems. For example, if a solution of (1.3) starts in the first quadrant of the xy-plane, will it stay in the first quadrant for all t?

Using the Floquet theory, we analyze a case of predator-prey interaction with periodic coefficients, and give a condition for the existence of a limit cycle.

In another direction we derive a multiplicity result on the periodic solutions for a class of periodic equations x′(t)=f(t,x(t)) , with f(t+p,x)=f(t,x) . x'(t) = f(t,x(t)){\kern 1pt} ,\;\;\;\;{\kern 1pt} {\rm with}\;f(t + p,x) = f(t,x){\kern 1pt} {\kern 1pt} . The connection to the rest of this paper is both by the method (it uses the same transformation u=1x u = {1 \over x} ), and by ecological applications, particularly to models of harvesting, see e.g., S. Ahmad and A.C. Lazer [1].

Consider the model (2.1) x′=x2(by−1)y′=−y2(dx−1) . \matrix{ {x' = {x^2}\left( {{b \over y} - 1} \right)} \cr {y' = - {y^2}\left( {{d \over x} - 1} \right){\kern 1pt} .} \cr } Here x(t) gives the number of prey, and y(t) the number of predator. If y(t) is small, the prey grows explosively (with x′ behaving like αx², α > 0). If the number of predators y(t) is large, then x′(t) < 0 and x(t) decreases. The number of predators y(t) decreases when x(t) small, and grows explosively for x(t) large. This model corresponds to (1.2), with a = d = 0. The coefficients e and f have been scaled out.

The system (2.1) has a rest point (d,b). Letting X=d/bx X = {{\sqrt {d/b} } \over x} and Y=1y Y = {1 \over y} transforms (2.1) into a perturbed harmonic oscillator (2.2) X′=−bd Y+d/bY′=bd X−1 . \matrix{ {X' = - \sqrt {bd} {\kern 1pt} Y + \sqrt {d/b} } \cr {Y' = \sqrt {bd} {\kern 1pt} X - 1{\kern 1pt} .} \cr } Setting X(t)=ξ(t)+1bd X(t) = \xi (t) + {1 \over {\sqrt {bd} }} , Y(t)=η(t)+1b Y(t) = \eta (t) + {1 \over b} leads to a harmonic oscillator (2.3) ξ′=−bd ηη′=bd ξ , \matrix{ {\xi ' = - \sqrt {bd} {\kern 1pt} \eta } \cr {\eta ' = \sqrt {bd} {\kern 1pt} \xi {\kern 1pt} ,} \cr } so that the solution of (2.2) is (2.4) X(t)=1bd+c1cosbd t−c2sinbd tY(t)=1b+c1sinbd t+c2cosbd t , \matrix{ {X(t) = {1 \over {\sqrt {bd} }} + {c_1}\cos \sqrt {bd} {\kern 1pt} t - {c_2}\sin \sqrt {bd} {\kern 1pt} t} \cr {Y(t) = {1 \over b} + {c_1}\sin \sqrt {bd} {\kern 1pt} t + {c_2}\cos \sqrt {bd} {\kern 1pt} t{\kern 1pt} ,} \cr } which is just a rotation of the point (X(0),Y (0)) around the point (1bd,1b) ({1 \over {\sqrt {bd} }},{1 \over b}) , the rest point of (2.2), on the circle of radius c12+c22 \sqrt {c_1^2 + c_2^2} , which depends on the initial conditions. The solution of (2.1) is then (2.5) x(t)=d/b1bd+c1cosbd t−c2sinbd ty(t)=11b+c1sinbd t+c2cosbd t . \matrix{ {x(t) = {{\sqrt {d/b} } \over {{1 \over {\sqrt {bd} }} + {c_1}\cos \sqrt {bd} {\kern 1pt} t - {c_2}\sin \sqrt {bd} {\kern 1pt} t}}} \cr {y(t) = {1 \over {{1 \over b} + {c_1}\sin \sqrt {bd} {\kern 1pt} t + {c_2}\cos \sqrt {bd} {\kern 1pt} t}}{\kern 1pt} .} \cr } The constants c₁ and c₂ are determined from the initial values (x(0),y(0)): (2.6) c1=d/bx(0)−1bd , and c2=1y(0)−1b . {c_1} = {{\sqrt {d/b} } \over {x(0)}} - {1 \over {\sqrt {bd} }}{\kern 1pt} ,\;\;{\kern 1pt} \;{\rm and}{\kern 1pt} \;\;{c_2} = {1 \over {y(0)}} - {1 \over b}{\kern 1pt} . It is now clear that the rest point (d, b) is a center for (2.1), and we can give a complete description of the behavior of positive solutions.

It remains to prove the lower bounds for the periodic solutions in the first part of the theorem. From (2.4) one sees that the positivity of X(t) and Y (t) implies that X(t)<2bd X(t) < {2 \over {\sqrt {bd} }} and Y<2b Y < {2 \over b} , from which one gets the lower bounds on x(t) and y(t).

Fig. 1

Consider the model (3.1) x′=a x+x2(by−1), x(0)>0y′=d y+y2(cx−1), y(0)>0 , \matrix{ {x' = a{\kern 1pt} x + {x^2}\left( {{b \over y} - 1} \right),\;\;x(0) > 0} \cr {y' = d{\kern 1pt} y + {y^2}\left( {{c \over x} - 1} \right),\;\;y(0) > 0{\kern 1pt} ,} \cr } with positive constants a,b,c and d. Each species grows explosively, if the number of the other one is small, while if the competitor's number is large, the growth is logistic-like. Clearly, the interaction is competitive in nature.

We begin with a simple observation: if x(0) > 0 and y(0) > 0, then x(t) > 0 and y(t) > 0 for all t > 0. Indeed, writing the first equation in the form x′ = A(t)x, with A(t)≡a+x(t)(by(t)−1) A(t) \equiv a + x(t)\left( {{b \over {y(t)}} - 1} \right) , and integrating, obtain x(t)=x(0)e∫0tA(s) ds>0 x(t) = x(0){e^{\int_0^t A(s){\kern 1pt} ds}} > 0 . Similarly, y(t) > 0 for all t > 0. Hence, we can limit our study of (3.1) to the first quadrant of the (x,y) plane.

Setting X=1x X = {1 \over x} and Y=1y Y = {1 \over y} produces a linear system (3.2) X′=−aX−bY+1 , X(0)=1x(0)>0Y′=−cX−dY+1 , Y(0)=1y(0)>0 , \matrix{ {X' = - aX - bY + 1{\kern 1pt} ,\;\;X(0) = {1 \over {x(0)}} > 0} \cr {Y' = - cX - dY + 1{\kern 1pt} ,\;\;Y(0) = {1 \over {y(0)}} > 0{\kern 1pt} ,} \cr } with a unique rest point (X₀,Y₀) given by (3.3) X0=d−bad−bc , Y0=a−cad−bc . {X_0} = {{d - b} \over {ad - bc}}{\kern 1pt} ,\;\;{Y_0} = {{a - c} \over {ad - bc}}{\kern 1pt} . Since x(t) > 0 and y(t) > 0 for all t > 0, we may restrict the system (3.2) to the first quadrant of the (X,Y) plane. The rest point (X₀,Y₀) lies in the first quadrant if either (3.4) d>b and a>c (and then ad>bc) , d > b\;{\kern 1pt} \;{\rm and}{\kern 1pt} \;a > c\;{\kern 1pt} \;({\rm and}\;{\rm then}\;ad > bc){\kern 1pt} {\kern 1pt} , or (3.5) d<b and a<c (and then ad<bc) . d < b\;{\kern 1pt} \;{\rm and}{\kern 1pt} \;a < c\;{\kern 1pt} \;({\rm and}\;{\rm then}\;ad < bc){\kern 1pt} {\kern 1pt} . Letting ξ = X − X₀ and η = Y −Y₀, we translate the rest point to the origin, obtaining the system (3.6) ξ′=−aξ−bηη′=−cξ−dη , \matrix{ {\xi ' = - a\xi - b\eta } \cr {\eta ' = - c\xi - d\eta {\kern 1pt} ,} \cr } with the matrix A=[−a−b−c−d] A = \left[ {\matrix{ { - a} & { - b} \cr { - c} & { - d} \cr } } \right] . The eigenvalues of A are (3.7) λ1,2=12(−a−d±a2−2ad+4bc+d2) . {\lambda _{1,2}} = {1 \over 2}\left( { - a - d \pm \sqrt {{a^2} - 2ad + 4bc + {d^2}} } \right){\kern 1pt} . The corresponding (column) eigenvectors are (3.8) ξ1,2=(−−a+d±a2−2ad+4bc+d22c,1)T . {\xi _{1,2}} = {\left( { - {{ - a + d \pm \sqrt {{a^2} - 2ad + 4bc + {d^2}} } \over {2c}},1} \right)^T}{\kern 1pt} . In case (3.4) holds, both eigenvalues are negative, and the rest point (X₀,Y₀) is a stable node, while in case (3.5) holds, one eigenvalue is negative and the other one is positive, so that (X₀,Y₀) is a saddle.

Fig. 2

We now consider a periodic perturbation of the explosive predator-prey model (4.1) x′=x2(b+β(t)y−1)y′=−y2(d+δ(t)x−1) , \matrix{ {x' = {x^2}\left( {{{b + \beta (t)} \over y} - 1} \right)} \cr {y' = - {y^2}\left( {{{d + \delta (t)} \over x} - 1} \right){\kern 1pt} ,} \cr } with small continuous functions β (t) and δ (t) of period p, so that β (t + p) = β (t) and δ (t + p) = δ (t) for all t. (We make no assumptions on the sign of β (t) and δ (t).) The linear system for X=1x X = {1 \over x} and Y=1y Y = {1 \over y} (4.2) X′=−(b+β(t))Y+1Y′=(d+δ(t))X−1 \matrix{ {X' = - \left( {b + \beta (t)} \right)Y + 1} \cr {Y' = \left( {d + \delta (t)} \right)X - 1} \cr } has p-periodic coefficients. Let F(t) be the normalized fundamental solution matrix (with F(0) = I, the identity matrix) of the corresponding homogeneous system (4.3) X′=−(b+β(t))YY′=(d+δ(t))X . \matrix{ {X' = - \left( {b + \beta (t)} \right)Y} \cr {Y' = \left( {d + \delta (t)} \right)X{\kern 1pt} .} \cr } For small β (t) and δ (t), F(t) is close for t ∈ [0, p] to the normalized fundamental solution matrix F0(t)=[cosbd t−bdsinbd tdbsinbd tcosbd t] {F_0}(t) = \left[ {\matrix{\hfill {\cos \sqrt {bd} {\kern 1pt} t} & \hfill { - \sqrt {{b \over d}} \sin \sqrt {bd} {\kern 1pt} t} \cr \hfill {\sqrt {{d \over b}} \sin \sqrt {bd} {\kern 1pt} t} & \hfill {\cos \sqrt {bd} {\kern 1pt} t} \cr } } \right] of the unperturbed system (4.4) X′=−bYY′=dX . \matrix{ {X' = - bY} \cr {Y' = dX{\kern 1pt} .} \cr } By the continuous dependence of eigenvalues on the coefficients of the matrix, the Floquet multipliers of (4.3), i.e., the eigenvalues of F(p) are close to the eigenvalues ρ₁ and ρ₂ of F₀(p). Clearly, (4.5) ρ1ρ2=1=detF0(p)ρ1+ρ2=2cosbd p=trace F0(p) . \matrix{ {{\rho _1}{\rho _2} = 1 = \det {F_0}(p)} \cr {{\rho _1} + {\rho _2} = 2\cos \sqrt {bd} {\kern 1pt} p = {\rm{trace}}{\kern 1pt} {F_0}(p){\kern 1pt} .} \cr }

We derive next an a priori bound on X_p(t) and Y_p(t), uniform in β (t) and δ (t), provided that |β (t)|+|δ (t)| ≤ c₀, for some constant c₀. Indeed, integrating both equations in (4.2) over (0,t), with t ∈ (0, p), taking absolute values and then adding the corresponding inequalities, obtain (4.6) |Xp(t)|+|Yp(t)|≤a1∫0t(|Xp(s)|+|Yp(s)|) ds+a2 , |{X_p}(t)| + |{Y_p}(t)| \le {a_1}\int_0^t \left( {|{X_p}(s)| + |{Y_p}(s)|} \right){\kern 1pt} ds + {a_2}{\kern 1pt} , for some positive constants a₁ and a₂. The desired bound over (0, p) follows by the Bellman-Gronwall lemma, see e.g., [2].

We claim that X_p(t) > 0 and Y_p(t) > 0 for all t. Setting Xp(t)=ξ(t)+1d {X_p}(t) = \xi (t) + {1 \over d} and Yp(t)=η(t)+1b {Y_p}(t) = \eta (t) + {1 \over b} in (4.2) obtain (4.7) ξ′=−b η−β(t)Yp(t)η′=dξ+δ(t)Xp(t) . \matrix{ {\xi ' = - b{\kern 1pt} \eta - \beta (t){Y_p}(t)} \cr {\eta ' = d\xi + \delta (t){X_p}(t){\kern 1pt} .} \cr } Express [ξ(t)η(t)]=F0(t)[c1c2]+F0(t)∫0tF0−1(s)f(s) ds , \left[ {\matrix{ {\xi (t)} \cr {\eta (t)} \cr } } \right] = {F_0}(t)\left[ {\matrix{ {{c_1}} \cr {{c_2}} \cr } } \right] + {F_0}(t)\int_0^t F_0^{ - 1}(s)f(s){\kern 1pt} ds{\kern 1pt} , with some constants c₁ and c₂, and f(t)=[−β(t)Yp(t)δ(t)Xp(t)] f(t) = \left[ {\matrix{ { - \beta (t){Y_p}(t)} \cr {\delta (t){X_p}(t)} \cr } } \right] . Since the vector [ξ(t)η(t)] \left[ {\matrix{ {\xi (t)} \cr {\eta (t)} \cr } } \right] has period p, and the fundamental solution matrix F₀(t) has period 2πmbd≠p {{2\pi m} \over {\sqrt {bd} }} \ne p , it follows that c₁ = c₂ = 0. The vector f (t) is small by our assumptions, and the a priori estimate (4.6). Both matrices F₀(t) and F0−1(s) F_0^{ - 1}(s) have bounded entries. Then the vector [ξ(t)η(t)] \left[ {\matrix{ {\xi (t)} \cr {\eta (t)} \cr } } \right] is small, so that the trajectory (X_p(t),Y_p(t)) remains near the point (1d,1b) \left( {{1 \over d},{1 \over b}} \right) , and hence it stays in the first quadrant for all t.

The transformation u(t)=1x(t) u(t) = {1 \over {x(t)}} of the preceding sections turns out to be useful for a class of first order equations with periodic coefficients. V.A. Pliss [6] considered the Abel equation: (5.1) x′(t)=a0(t)x3+a1(t)x2+a2(t)x+a3(t) . x'(t) = {a_0}(t){x^3} + {a_1}(t){x^2} + {a_2}(t)x + {a_3}(t){\kern 1pt} . Assuming that the given functions a_i(t), 0 ≤ i ≤ 3, are of period p, and a₀(t) is either positive or negative for all t, he proved that the equation (5.1) has at most three p-periodic solutions. The proof involved a clever combination of the equations that the inverses of solutions satisfy.

What if one changes the a₀(t)x³ term to a₀(t)x²ⁿ⁺¹? In case it is a₀(t)x⁵, the method of V.A. Pliss [6] still gives the same result with a little extra effort. For higher powers things get more involved, and in fact existence of at most three p-periodic solutions was proved by another elegant method in A.A. Panov [5]. It turns out that the following more general result was already known.

This theorem follows from a more general result of A. Sandqvist and K.M. Andersen [7]. They considered the equation (5.2) on the interval (0, p) and called a solution to be closed if x(p) = x(0). Assuming the condition (5.4) holds, they showed that the problem (5.2) has at most three closed solutions, which implies the Theorem 5.1.

A simpler proof of the Theorem 5.1 was found in P. Korman and T. Ouyang [3]. We now simplify the presentation in that paper. The proof will follow from the following three simple lemmas.

The proof is standard, and it can be found in e.g., P. Korman [2], p. 245. The next lemma is crucial.

Turning to the proof of the Theorem 5.1, observe that different solutions of (5.2) do not intersect by the uniqueness theorem. If the equation (5.2) has four p-periodic solutions, let ξ (t) be the smallest one. Then z(t) = x(t) − ξ (t) satisfies (5.7) z′=f(t,z+ξ)−f(t,ξ)≡g(t,z) , z' = f(t,z + \xi ) - f(t,\xi ) \equiv g(t,z){\kern 1pt} , and the equation (5.7) has three positive p-periodic solutions. However, g(t,0) = 0 and g_zzz(t,z) > 0 for z > 0, contradicting the Lemma 5.3.

Equations of the type (5.2) occur often in ecological problems, see e.g., S. Ahmad and A.C. Lazer [1], or P. Korman [2].