On the Stability of Critical Point for Positive Systems and Its Applications to Biological Systems

Abstract

The coexistence and extinction of species are important concepts for biological systems and can be distinguished by an investigation of stability. When determining local stability of nonlinear systems, Lyapunov indirect method based on the Jacobian linearization has been widely employed due to its simplicity. Despite such popularity, it is not applicable to singular systems whose Jacobian has at least one eigenvalue that is equal to zero. In such singular cases, an appropriate Lyapunov function should be sought to determine the stability of systems, which is rather difficult and quite involved. In this paper, we seek for a simple criterion to determine stability of the equilibrium that is located at the boundary of the positive orthant, when one of eigenvalues of the Jacobian is zero. The goal of the paper is to present a generalized condition for the equilibrium to attract all trajectories that starting from initial condition in the positive orthant and near the equilibrium. Unlike the Lyapunov direct method, the proposed method requires just a simple algebraic computation for checking the stability of the critical point. Our approach is applied to various biological systems to show the effectiveness of the proposed method.

1. Introduction

Stability is an important concept that can provide powerful insight into qualitative behavior of various biological systems [1, 2] such as predator-prey systems, viral and immune systems, epidemic systems and so on. Stability analysis tells us whether a solution trajectory near an equilibrium will converge toward or move away from it. For predator-prey system [3-15], an investigation of stability distinguishes the coexistence of predators and prey from the extinction of some species. In addition, stability analysis may determine the basin of attraction of each equilibrium, that is, the region of coexistence and the region of extinction. Aside from predator-prey systems, some research works on virus infection of CD4 + cells have been performed in [16-20]. In their works, mathematical models describing infection dynamics have infection-free equilibrium and chronic-infection equilibrium. The infection will die out if the former is stable, while the infection will persist if the latter is stable. For the epidemic model with some vaccination strategy [21-23], an examination of stability can predict whether the vaccination successfully leads to disease eradication. An epidemic outbreak is produced if an endemic equilibrium is stable, whereas the disease will die out if a disease-free equilibrium is stable.

In order to check the local stability of a nonlinear system at its certain equilibrium, Lyapunov indirect method based on the Jacobian linearization has been widely employed due to its simplicity. Despite such popularity, the method cannot be applied to determine stability when one eigenvalue of the Jacobian matrix is equal to zero. In fact, although the research works in [2-12, 15-18] have used the Jacobian matrix to determine stability, most of them have not provided any conclusion on stability when the Jacobian matrix has an eigenvalue at the origin of the complex plane. While the extinction of some species can be predicted by investigating the stability of equilibria with some zero component, stability of such equilibria has not been determined in [2-12, 15-18] since the Jacobian has an eigenvalue at the origin.

In order to analyze stability for such cases, an appropriate Lyapunov function should be constructed as in [24], or some other nonlinear stability theory needs to be employed [25]. Even though Lyapunov stability theory and LaSalle’s invariance principle can be employed to determine stability [26], it is rather difficult and quite involved to construct an appropriate Lyapunov function. An alternative method may be to study the stability of reduced system by invoking the center manifold theory [26, 27]. However, the verification procedure is not convenient in general, because a solution to the partial differential equation for the center manifold needs to be obtained. To avoid the difficulty in solving the partial differential equation, the authors have presented in [28] an alternative technique under a certain assumption of the Jacobian structure.

In this paper, we extend the result of [28], and, from the extension, a variety of ecological systems are investigated about their stability. The motivation of this research lies in the observation that there might exist a boundary equilibrium that attracts all trajectories starting in the positive orthant although it is not stable in the usual sense. A condition is presented for the equilibrium to attract all trajectories that reside in the positive orthant and near the equilibrium (even though it may not attract those trajectories outside the positive orthant). The proposed condition can deal with the case to which the technique of [28] is not applicable (motivated by the models in [6, 9]). We also extend the applicability of the previous technique by relaxing the assumption on the Jacobian structure. To show the effectiveness of the proposed method, we apply the proposed condition to inspect the stability of several different predator-prey models [6, 8, 9] and an HIV infection model [17]. Although the authors of [6, 8, 9, 17] have discussed the stability of various equilibria, they have not determined the stability of boundary equilibria when one of eigenvalues of the Jacobian is zero. Moverover, the result of [28] cannot be used for [6, 9], either. In contrast to those previous results, the proposed method can determine the stability in such singular cases.

The paper is organized as follows: Section II recalls the concept of stability with respect to the positive orthant and presents a simple condition for checking the stability of boundary equilibrium. In Section III, we apply the proposed method to several predator-prey models and an HIV infection model to show the effectiveness of the proposed method. Finally, some concluding remarks are given in Section IV.

Notations: A function is said to be of class Ck if it is continuously differentiable k times. For a vector x and a matrix A , the i-th component of x and the i-th row of A are denoted by x (i) and A(i) , respectively, and when there is no confusion, x (i) is abbreviated to xi . For a matrix A , AT denotes a transpose of the matrix A . We denote by ek the column vector [0 0 ··· 0 1 0 ··· 0]T with the entry 1 in the k-th place. The elementary matrix obtained by interchanging the first and k-th rows of identity matrix is denoted by Ek . The n × n identity matrix and the r × 1 zero vector are denoted by In and 0r , respectively, and when there is no confusion, 0r is abbreviated to 0 . When all eigenvalues of a matrix A have negative real parts, A is called a Hurwitz matrix. The maximum and the minimum eigenvalue of matrix A are represented by λmax ( A ) and λmin ( A ) , respectively. For some r > 0 and x 0 ∈ℝn , B (x0 , r) = ｛ x∈ ℝn｜∥x ˗ x0∥ < r ｝ , where ∥x∥ stands for the Euclidean norm of a vector x. For a vector x, we write x >> 0 and x ≥≥ 0 to indicate that every component of x is positive and nonnegative, respectively. Let = ｛ x ∈ ℝn : x ≥≥ 0 ｝ and = ｛ x ∈ ℝn : x ≥≥ 0 ｝ . The order of magnitude notation ο is used as follows: we say f (x) =ο (g (x)) if, for each ε > 0, there exists δ > 0 such that ｜f (x)｜ ≤ε ｜g (x)｜ for ｜x｜ <δ . A set M is said to be a positively invariant set for the system ẋ = f ( x) if the solution x (t ) satisfies x (0)∈M ⇒ x (t )∈M, ∀t ≥ 0.

 

2. Stability Condition for Positive Nonlinear Systems

In this paper, we consider a class of nonlinear system in the absence of an input u, that is, so-called unforced system

where f :ℝn →ℝn is assumed to be C3 and is assumed to be a positively invariant set for system (1). We also assume that there exists an isolated equilibrium x* that is located on the boundary of and satisfies the following assumption.

Assumption JS (Jacobian Structure): Suppose that the Jacobian matrix at x* , has one eigenvalue at the origin and all others with negative real parts. Moreover, there exist a nonsingular matrix M and an integer k (1 ≤ k ≤ n ) such that M(k) ≥≥ 0 and

and

Remark 1: Usually, M is chosen as an identity matrix, for which (2) and (3) become

However, an introduction of the matrix M allows a larger class of systems to satisfy Assumption JS. In fact, as will be illustrated in Section 3.4., the HIV infection model of [17] does not satisfy (4) but assumption JS (if M is chosen as an appropriate matrix).

Although Assumption JS seems to be quite restrictive at first glance, it is not. In fact, the class of systems under consideration includes the following typical form [2-11, 14, 15, 19, 20], which often arises in various biological systems:

Here, we assume that the equilibrium of our interest satisfies that

It is easily seen from (6) that the class of systems (5) satisfies Assumption JS. However, it should be noted that not all systems satisfying Assumption JS can be written as in (5). (As a matter of fact, the HIV infection model of [17] is not of the form (5), but satisfies Assumption JS.)

To illustrate the basic concept of the proposed method, consider the basic 2-species Lotka-Volterra competition model [2]:

The system (7) has several equilibria, but we are interested in the stability of equilibrium (0, 1). The Jacobian matrix at (0, 1) is computed as

which has eigenvalues 0 and 1, and stability of (0, 1) is not determined. In order to see what is going on in this case, we check the phase plane1) in Fig. 1. It is seen that all trajectories starting from initial conditions in converge to the equilibrium (0, 1), while those starting outside do not converge to (0, 1). Thus, this equilibrium is not stable in the usual sense. However, if we restrict our interest to initial conditions in (the shaded region in Fig. 1), we may say that the equilibrium (0, 1) is asymptotically stable with respect to since it is attracting all trajectories in .

Fig. 1.Phase portrait of system (7). In this figure, solid red circles indicate equilibria of system (7). Although phase trajectories starting from initial conditions outside do not converge to the equilibrium (0, 1), those trajectories starting from (shaded region) converge to (0, 1).

Definition 1: The equilibrium point x* is locally stable with respect to (w.r.t.) the set if, for each ε < 0 , there exists δ (ε ) > 0 such that x (0) ∈ B( x* , δ )∩ ⇒ x (t ) ∈ B( x* , ε)∩ .

Moreover, it is locally asymptotically stable w.r.t. if it is stable w.r.t. and δ can be chosen such that x (0) ∈ B( x* , δ )∩ ⇒

Now, we provide a simple condition for the stability w.r.t. that does not require the existence of an appropriate Lyapunov function. From Assumption JS, we obtain (MAM˗1)(k) [0 0 … 0],

which implies

for some matrix A21 ∈ℝ (n˗1)×1 and A2 ∈ℝ (n˗1) × (n˗1) . (Recall that premultiplication and postmultiplication by an elementary matrix results in elementary row and column operation, respectively.) Note that A2 is a Hurwitz matrix because of Assumption JS.

Lemma 1: Under Assumption JS, the matrix

Satisfies

Proof: Let

Then, we obtain T = Ek M and

Moreover, it follows from (8) and (10) that

which results in

With matrix T of Lemma 1, let g1 : ℝn →ℝ1 and g2 : ℝn →ℝn˗1 be functions

where =[z2, z3,···, zn ]T . In addition, for nonnegative integers m, k and 2 ≤ i1,···, ik ≤ n, let

Theorem 1: For an equilibrium x* of system (1), suppose that x* and (1) satisfy Assumption JS. Then, we have the followings:

(Case 1) x* is locally asymptotically stable w. r. t. (respectively, unstable) if P[2,0] < 0 (respectively, P[2,0] > 0 ).

(Case 2) Suppose that P[2,0] = 0. Then, x* is locally asymptotically stable (respectively, unstable) if (respectively, P[3, 0] + > 0 ), where

Remark 2: At first glance, it may seem rather complicated to apply Theorem 1. But, as will be illustrated in Section 3, the application of the theorem is not difficult since it only requires numerical computations. Moreover, when M = I, it is convenient to use , where

Proof: (Case 1): Let x ∈ and

Then, we have

and

Since is C3 and = 0, (14) can be represented as where is C3 and

With matrix T of Lemma 1, the change of variables transforms (14) into where

From (9), the system is rewritten as

where g1 and g2 are C3 functions such that g = [g1, ]T and, for i = 1, 2,

Thus, g1 (z1, ) = g1(z1, ···, zn ) can be written as g1(z1, ···, zn)

Moreover, according to center manifold theorem (e.g., see Theorem 8.1 of [26]), there exists δ > 0 and a C1 functionπ : ℝ1 → ℝn˗1, defined for all │ z1│ < δ1, such that =π (z1 ) =: ［π2 (z1 ), ···, πn (z1 ) ］T is a center manifold for (16), i.e.,

and

From (19), π can be written as

for some . Therefore, we get

where and

Now, another change of variables transforms (16) into

Subtracting (21) into (22) yields and subtracting (18) from (23) gives

Hence, we can rewrite (22) and (23) as

where

It is easily seen that N1 and N2 are C2 and

and, by virtue of (17),

Thus, there exists δ2 > 0 such that ∥[ z1, wT ]T∥ <δ2 implies ∥Ni (z1, w)∥ ≤ ki∥w∥, i = 1,2, where k1 and k2 can be made arbitrary small by choosing δ2 small enough.

Now, suppose that P[2,0] < 0 and consider a Lyapunov function candidate for full system (24) where P0 is the positive definite solution of

Then, there exists δ3 < min ( δ1, δ2 ) such that ∥[ z1, wT ]T∥ < δ3 implies

Thus, the derivative of V along the trajectories of the system (24) is given by

On the other hand, from Lemma 1, we obtain which implies that z1 ≥ 0 since x ∈ and M(k) ≥≥ 0. Therefore, using the standard Lyapunov stability techniques, it can be shown that ( z1 ( t ), w( t )) →0 as t →∞. Since π (0) = 0, it follows that ( z1 ( t ), z2 ( t )) →0 as t →∞ and as a consequence, x (t )→ x* as t →∞ for x (0)∈ B( x* ,δ ) ∩ for some δ > 0.

On the contrary, suppose that P[2,0] > 0. For the convenience, system (24) is rewritten as

in which the origin corresponds to the equilibrium x* of (1). For given , let be the solution with Then, cannot be kept within a small neighborhood of the origin because and ˗C2 = P[2,0] > 0. (In fact, > 0 when is positive and sufficiently small.) But, observe that (, 0) is also a solution of system (27). This implies that the origin is unstable since some initial condition ( , 0) with > 0 moves away from the origin. Therefore, x* is unstable.

(Case 2): From (17), g2 (z1, ) can be written as where C12 and C22 are some matrices of appropriate dimensions. Since

letting yields

Thus, since , we get from (21)

Now, suppose that P[3, 0] + < 0, i.e., c3 > 0. Then, with the same Lyapunov function candidate as in the proof of (Case 1), there exists a δ4 > 0 such that ∥[z1, wT]T∥ <δ4 leads to

This implies that there exists a δ > 0 such that x (t )→ x* as t →∞ for x (0) ∈ B( x* ,δ ).

 

3. Applications to Biological Systems

To show the effectiveness of the proposed method we apply our result to several predator-prey models and an HIV model in the ecological literatures.

3.1 A predator-prey model with hawk and dove tactics

In this subsection, we consider the predator-prey model that incorporates individual behavior of the predators [6]. It is assumed that individual predators can use two tactics when fighting with other predators to keep a captured prey, the hawk and dove tactics. While the hawk fights in any case, the dove is never aggressive. When a hawk meets a dove, there is no fighting and the hawk keeps the prey and dove gets nothing. When two hawks encounter, they fight and both of them suffer injuries. After the fighting the winner keeps the prey. When two doves encounter, there is no fighting and they share the prey. The model consists of two parts: a fast part that represents the change of tactics of predators and a slow part that expresses the predator-prey interactions. Let n, PH, and PD be the population of prey, hawk predator, and dove predator, respectively. Then, the system model is given by [6]

where and

with G = an. Taking advantage of the two-time scales, the complete system has been studied from the system model of reduced dimension. When , , the aggregated model is given by [6]

where p represents the total population of predators, i.e. p = pH + pD. The system (30) has five equilibria, E1(0, 0), E2 (k, 0), E3(, ), E4(, ), E5(, ), where , and =

We first consider the case where αC > 8μ. It has been shown in [6] that, when αC > 8μ , E2 is stable (respectively, unstable) if K < or < K < (respectively, < K < ). But, the stability of E2 has not been determined when K = or K = . Now, using the proposed theorem, we are going to determine the stability of 2 E for such cases. Since the Jacobian matrix at E2 when K = is computed as

Assumption JS is satisfied with M = I2 and k = 2. More-over, with D1 := Cα a ˗ , we obtain

which leads to

since K = < . Therefore, E2 is locally asymptotically stable w.r.t. and, as a result, the predator goes to extinction when K = . On the other hand, for the case where K = , we can show that where D2 = Cα a + . Therefore, E2 is un-stable, and hence predator and prey coexist when K = .

Next, we consider the stability of E2 when αC = 8μ for which = = . Suppose that we need to check the stability of E2 when K = . Again, the stability has not been determined in [6] since the Jacobian has one eigenvalue at the origin. It should be also noted that the result of [28] cannot be used to determine the stability, while Theorem 1 can be used. In contrast to the case where αC > 8μ , we need to rely on (Case 2) of Theorem 1 since αC = 8μ implies p[2,0] = 0. Since And

we obtain

Thus, when αC = 8μ and K = , E2 is locally asymptotically stable w.r.t. .

3.2 A partial-dependent predator-prey system

In this subsection, we consider a partial-dependent predator-prey system that has been discussed in [9]. The model is given by

where r, s, K, L, α, β, and α1 are some positive constants.

The state variables x, and y represent the population density of prey and predator, respectively. The system (32) has several equilibria such as E O (0,0), E C (K,0), and so on. It has been shown in [9] that the system (32) experiences a saddle-node bifurcation at the equilibrium E P when r = α L and α1 Ks − s −β K ≠ 0, and a pitchfork bifurcation at the equilibrium E P when r =α L and α1 Ks − s −β K = 0. Although it has been shown in [9] that E P is a stable node (respectively, a saddle) r <α L if (respectively, r >α L ), the stability of E P has not been determined when r =α L. Now, using the proposed theorem, we are going to determine the stability ofE P when r =α L. The Jacobian matrix at E P is computed as

which shows that Assumption JS is satisfied with M = I2 and K = 1. Moreover, we obtain so that

Therefore, EP is locally asymptotically stable (respectively, unstable) w.r.t. if α1 Ks − s −β K < 0 (respectively, α1 Ks − s −β K > 0 ).(See Fig. 2 and Fig. 3 for each cases.)

Fig. 2.Phase portrait of system (32) with α = 1, K = 2, s = 3, L = 2, β = 1, r = 2, and α1 = 1/2, for which r =α L and α1 Ks − s −β K < 0 . In this figure, solid red circles indicate equilibria of system (32). Although not all trajectories move toward to EP (0, 2), all trajectories residing in (shaded region) converge to EP.

Fig. 3.Phase portrait of system (32) with α = 1, K = 8, s = 3, L = 2,β = 1, r = 2, and α1 = 1/2, for which r =α L and α1 Ks − s −β K > 0. All trajectories starting from initial conditions in move away from EP .

When α1 Ks − s −β K = 0 as well as r =α L, the determination of stability of E P becomes more complicated. Although the stability of E P cannot be determined by the result of [28], Theorem 1 (Case 2) can be applied. To this end, we compute And which imply that

Thus, when r =α L and α1 Ks − s −β K = 0, E P is locally asymptotically stable (respectively, unstable) if s > K β (respectively, s

Fig. 4.Phase portrait of system (32) with α = 1, K = 3, s = 6, L = 2,β = 1, r = 2, and α1 = 1/2, for which r =α L, α1 Ks − s −β K = 0, and s < Kβ . All trajectories move toward to E P(0,2), which implies that EP is asymptotically stable in the usual sense.

Fig. 5.Phase portrait of system (32) with α = 1, K = 6, s = 3, L = 2,β = 1, r = 2, and α1 = 1/2, for which r =α L, α1 Ks − s −β K = 0, and s < Kβ . Most trajectories move away from EP (0, L) in the direction of unstable manifold, which implies that E P is unstable.

3.3 A predator-prey model with habitat destruction

In this subsection, we consider the predator-prey model that has been discussed in [8], which was developed to study the consequences of habitat destruction. The model is given by

where cx , ex , μ , cy , ey , φ , and D are some positive constants. The state variables x and y represent the proportion of patches occupied by prey and predators, respectively.

Four possible equilibria exist for the model (36):

E 0 (0,0) , , and E 3 ( x*, y* ), where and It has been shown in [8] that E 1 is stable (respectively, unstable) if D > Dx (respectively, D < Dx ), where Dx := However, when D = Dx , the stability of E 1 has not been determined since the Jacobian matrix has one eigenvalue at the origin.

Now, using the proposed theorem, we are going to show that E 1 is indeed stable when D = Dx . The Jacobian matrix at E 1 is computed as

from which Assumption JS is satisfied with M = I 2 and k = 2. Moreover, we obtain which implies that

Thus, we get which leads to

Therefore, E 1 is locally asymptotically stable w.r.t. when D = Dx .

3.4 HIV dynamics model

In this subsection, we consider the HIV model that has been discussed in [17], which depicts the interaction between T cells and virus particles. The model is given by

where k, β , γ , and N are some positive constants, i can be either 1 or 0, and f is a smooth function that satisfies f (T) > 0, 0 ≤ T < , f ( ) = 0, f ' ( ) < 0, and f (T ) < 0, T > .

The state variables T, T* , and V represent the population of uninfected T cells, productively infected T cells, and free virus particles, respectively. The system (38) has two equilibria, E 0 ( , 0 , 0 )and E e (T e , , V e ), where , and. It has been shown in [17] that E 0 is stable (respectively, unstable) if R 0 := (respectively, R 0 > 1 ). But, the stability of E 0 has not been determined when R 0 = 1.

Now, using the proposed theorem, we are going to determine the stability of E 0 when R 0 = 1, that is, γ = kN ˗ki. The Jacobian matrix, evaluated at E 0 , is

which shows that Assumption JS is not satisfied with M = I 3 . If we choose

then

from which Assumption JS is satisfied with k = 3. With some computation, we obtain

which is negative since f ' ( ) < 0 and N is typically large. Therefore, according to Theorem 1, E 0 is locally asymptotically stable w.r.t. , which suggests that the virus population will decline and die out.

 

4. Conclusions

We have proposed a simple test for checking local stability of an equilibrium that is located on the boundary of the positive orthant. While Laypunov’s indirect method based on the Jacobian linearization cannot draw any conclusion on the stability when the Jacobian matrix has an eigenvalue at the origin, the proposed method is able to determine the stability. Since our approach is based on the approximate solution to center manifold equation, it only guarantees the local result. To the contrary, if the global behavior needs to be studied, rather elaborate tools such as Lyapunov direct method or LaSalle’s invariance principle should be resorted to. Nonetheless, we believe that our results are quite attractive in that it requires just a simple algebraic computation but may provide a clue to the prediction of global behavior. One might think that phase portrait is sufficient to determine the stability. But, the phase portrait is hardly applicable to higher dimensional systems (e. g., ℝn with n ≥ 3 ), while the proposed method can be easily employed. Moreover, we believe that applying our result requires less time and effort than sketching a complete phase portrait.