Dynamic Behaviors of a Stage Structure Commensalism System with

Holling type II Commensalistic Benefits

FENGDE CHEN, ZHONG LI, LIJUAN CHEN

College of Mathematics and Statistics

Fuzhou University

No. 2, wulongjiang Avenue, Minhou County, Fuzhou

CHINA

Abstract: - Noting the fact that commensal species that behave as foragers are subject to the constraints of han-

dling time, a two species commensalism model with Holling type II commensalistic benefits and stage structure is

proposed and studied. We first show that among four possible equilibria, host-only equilibrium and positive equi-

librium are possible asymptotically stable. Next, we establish a powerful lemma on the global stability property

of the single species stage structured model with linear perturbation on mature species. By applying this lemma

and the differential inequalities theory, sufficient conditions which ensure the global attractivity of the host-only

equilibrium and positive equilibrium are obtained, respectively. Our results generalize some known results.

Key-Words: Stage structure; Commensalism; Comparison theorem; Holling II functional response; Global

attractivity

Received: September 30, 2022. Revised: November 3, 2022. Accepted: November 17, 2022. Published: December 12, 2022.

1 Introduction

The aim of this paper is to investigate the dynamic be-

haviors of the following commensalism system with

Holling II functional response and stage structure:

dx1

dt =αx2−βx1−δ1x1,

dx2

dt =βx1−δ2x2−γx2

2+dy

η1+η2yx2,

dt =y(b2−a2y),

(1)

where α, β, δ1,δ2, d, b2, a2, η and γare all positive

constants, x1(t)and x2(t)are the densities of the im-

mature and mature commensal species at time t,y(t)

is the density of the host species at time t. The follow-

ing assumptions are made in formulating the model

(1):

(A1) The birth rate of the immature commensal

species is proportional to the existing mature popula-

tion with a proportionality constant α; for the imma-

ture commensal species, the death rate and transfor-

mation rate of mature are proportional to the existing

immature population with proportionality constants β

and δ1;

(A2) The death rate of the mature commensal species

is proportional to the existing mature population with

a proportionality constant δ2. The mature commensal

species has density restriction, γis density dependent

coefficient;

(A3) The host species satisfies the logistic model;

(A4) The host species only benefits to the mature

commensal species with rate dy

η1+η2y, which is of

Holling II type.

During the last decade, many scholars investigated

the dynamic behaviors of the mutualism model or

commensalism model [1]-[39]. By means of com-

mensal interactions, one species benefits from the

other without either harming or benefiting the lat-

ter. Commensalism is very common in nature, for

example, small plants called epiphytes and the large

tree branches on which they grow. Epiphytes de-

pend on their hosts for structural support but do not

derive nourishment from them or harm them in any

way. Another example is the remora (family Echinei-

dae) that rides attached to sharks and other fishes.

One could refer to https://www.britannica.com/sci-

ence/commensalism for more background and exam-

ples.

Despite commensalism is a very common rela-

tionship in nature, its mathematical modeling lags

far behind. In 2013, for the first time, Sun and

Sun[42] first time proposed a two species commen-

salism model. Since then, many scholars worked on

this direction. Topics such as the influence of har-

vesting [2, 3, 5, 9, 10, 11, 25, 26, 32, 37, 38], the

existence of periodic solution or almost periodic so-

lution [6, 8, 12, 21, 23], the stability of the system

[1, 2, 3, 4, 7, 18, 19, 20, 25, 31, 37, 38, 39, 42], the per-

sistent and extinction property of the system [2, 19],

the influence of stage structure [13], the influence of

Allee effect [14, 15, 16, 17, 22, 29, 30, 34, 35], the in-

fluence of time delays [1, 6, 24], the influence of feed-

back control [20, 31], the bifurcation phenomenon of

the system [26, 28, 29, 30, 34, 35] etc were exten-

sively investigated.

In the real world, almost all animals have the

stage structure of immature and mature. In differ-

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

810

Volume 21, 2022

ent stages, the species have different reactions to

the environment. many scholars investigated the dy-

namic behaviors of the stage structured species, see

[1, 13, 32, 40, 44, 45] and the references cited therein.

Aiello and Freedman [44] first time proposed the fol-

lowing stage-structured single species model.

dx1(t)

dt =αx2(t)−γx1(t)

−αe−γτ x2(t−τ),

dx2(t)

dt =αe−γτ x2(t−τ)−βx2

2(t).

(2)

They showed that the above system admits a unique

positive equilibrium which is globally asymptotically

stable. Such a result is similar to the traditional single

species Logistic model. However, in above system,

the authors did not consider the influence of death rate

of mature species. Chen, Xie, Chen [1] proposed and

studied the following May type stage-structured co-

operation model,

˙x1(t) = b1e−d11 τ1x1(t−τ1)−d12x1(t)

−a11x2

1(t)

c1+f1x2(t)−a12x2

1(t),

˙y1(t) = b1x1(t)−d11y1(t)

−b1e−d11 τ1x1(t−τ1),

˙x2(t) = b2e−d22 τ2x2(t−τ2)−d21x2(t)

−a22x2

2(t)

c2+f2x1(t)−a21x2

2(t),

˙y2(t) = b2x2(t)−d22y2(t)

−b2e−d22 τ2x2(t−τ2).

(3)

In this system, d12 and d21 represent the death rates

of the first and second mature species, respectively.

They showed that despite the cooperation between the

species, the species may still be driven to the extinc-

tion due to the stage structure. Death rate of mature

species plays crucial role on the persistence and ex-

tinction property of the system.

Instead of consider the influence of delay, some

scholars [13], [32], [40], [44] assumed that there

are proportional number of immature species be-

comes mature species. Khajanchi and Banerjee [44]

proposed the following stage structure predator-prey

model with ratio dependent functional response

dx1

dt =αx2(t)−βx1(t)−δ1x1(t),

dx2

dt =βx1(t)−δ2x2(t)−γx2

2(t)

−η(1 −θ)x2(t)y(t)

g(1 −θ)x2(t) + hy(t),

dt =uη(1 −θ)x2(t)y(t)

g(1 −θ)x2(t) + hy(t)−δ3y(t).

(4)

Here the authors assumed that the prey species is stage

structured. The authors only considered the stabil-

ity property of predator-extinction equilibrium and

the positive interior equilibrium. We would like to

point out that Xiao and Lei [12] showed that the death

rate of the mature species plays important role on the

persistence and extinction of the single species stage

structured system.

In our opinion, the commensalism models were

not well studied in the sense that up today, just one

paper [13] considers the influence of the stage struc-

ture. A recent study [43] showed that the relationship

among Brazil nut and three frog species is commen-

salism, and it is well known that frog species are stage

structured species. Hence, it is necessary to propose

some modeling on stage structured commensalism. In

[13], Lei proposed the following two species com-

mensalism model

dx1

dt =αx2−βx1−δ1x1,

dx2

dt =βx1−δ2x2−γx2

2+dx2y,

dt =y(b2−a2y).

(5)

Concerned with the global stability property of the

equilibria of system (5), by constructing suitable Lya-

punov function, the author obtained the following re-

sult (Theorem 3.1 and 3.2 in [13]):

Theorem A. (1) if

(β+δ1)δ2−db2

a2−αβ > 0(6)

hold, then A2(0,0,b2

a2)is globally asymptotically sta-

ble. If

(β+δ1)δ2−db2

a2−αβ < 0(7)

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

811

Volume 21, 2022

hold, then A4(x1, x2, y)is globally asymptotically

stable, where

x1=αx2

β+δ1

x2=αβ −δ2−dy(β+δ1)

(β+δ1)γ,

y=b2

(8)

Traditional two species Lotka-Volterra coopera-

tion system, which takes the form

dt =x(r1−a11x+a12y),

dt =y(r2−a21y+a22x),

is criticized by many biologist. The reason lies in the

unrealistic assumption that the benefits of the interac-

tion were unlimited and increased in direct proportion

to the density of the mutualistic partner. To overcome

this, and inspired by the fact that mutualists that be-

have as foragers are also subject to the constraints of

handling time, for example, ( for a species of solitary

bee visiting a flower species. The rate of collection

of pollen is limited by the handling time per plant),

Wright[43] proposed the following two species mu-

tualism model

dt =Nr1(1 −c1N) + baM

1 + aT11M,

dt =Mr2(1 −c2M) + cdN

1 + dT22N,

where the author used Holling type II functional re-

sponse baM

1 + aT11Mto describe the feeding rate (items

per unit of time, T11 is the handling time, bis a coeffi-

cient converting Mto new units of N. The model

is based on the traditional logistic equation with a

term added to include the per capita benefits of in-

teracting with the population of the mutualist partner.

The model may have none, one or two positive equi-

librium, i.e., with the introduction of functional re-

sponse, the dynamic behaviors of the system becomes

more complex.

Probably inspired by the works of Wright [43] or

some similar works, several scholars [7],[15],[16],

[17],[23],[28],[31],[39] argued that the the relation-

ship of two commensalism model should be described

by the suitable functional response. Li, Lin and

Chen [23] for the first time adopt the idea of func-

tional response of predator prey system, they pro-

posed the following discrete commensalism model

with Holling II functional response.

x1(k+ 1) = x1(k)exp na1(k)−b1(k)x1(k)

+c1(k)x2(k)

e1(k) + f1(k)x2(k)o,

x2(k+ 1) = x2(k)exp {a2(k)−b2(k)x2(k)}.

They studied the positive periodic solution of the sys-

tem.

Wu [7] argued that a relationship of nonlinear type

between two species is more feasible, and she estab-

lished the following two species commensal symbio-

sis model

dt =xa1−b1x+c1yp

1 + yp,

dt =y(a2−b2y),

(9)

where ai, bi, i = 1,2, p and c1are all positive con-

stants, p≥1. The results of [7] is then generalized

by Lei [23] and Wu, Li and Lin [16] to the commen-

salism model with Allee effect.

Recently, Jawad [26] proposed the following com-

mensalism model with Michaelis-Menten type har-

vesting and Holling type II functional response:

dt =ru1−u

k+βuv

α+u−qEu

cE +lu ,

dt =sv1−v

m−dv,

where u(t)and v(t)denote the densities of the first

and second species at time t, respectively. Topics

such as permanence, saddle node bifurcation were

discussed in [26].

Chen, Chong and Lin [31] proposed the follow-

ing commensal symbiosis model with Holling type II

functional response and feedback controls:

dt =xb1−a11x+a12y

a13 +a14y−α1u1,

dt =y(b2−a22y−α2u2),

du1

dt =−η1u1+a1x,

du2

dt =−η2u2+a2y,

(10)

where x(t)and y(t)denote the density of the first and

second species at time t, and u1and u2are feedback

control variables. By developing some new analytical

technique, the authors showed that the unique positive

equilibrium is globally attractive.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

812

Volume 21, 2022

Xu, Xue, Xie and Lin [39] proposed and studied

the commensalism model with Crowley-Martin func-

tional response. i.e.,

dt =x−a1−b1x+c1y

d1+e1x+f1y+g1xy ,

dt =y(a2−b2y).

For autonomous case, i.e., when all the coefficients of

the system are positive constants, authors investigated

the local and global dynamic behaviors of the sys-

tem. For nonautonomous case, authors investigated

the persistent and extinction properties of the system.

Li and Wang [28] argued that a suitable model

should include some past state of the system, and they

investigated the dynamic behaviors of the following

commensalism system

dt =xa1−b1x+c1y

1 + e1x+f1y,

dt =y(a2−b2y(t−τ)).

As was shown above, there are several kinds of

functional response used in modelling the commen-

salism model, however, it seems that the discussion

of Wright [43] is one of the most reasonable, and

the biological explanation is plausible. Stimulated by

the above works, specially stimulated by the work of

Wright [43], we propose the system (1). As far as sys-

tem (1) is concerned, since it seems that the system is

similar to system (5), only with the cooperation term

dx2yin system (5) changed to the term with func-

tional response dy

η1+η2yx2in system (1). One may

expect the analysis method used in Lei [13] could be

directly applied to system (1), i.e., one could investi-

gate the stability property of the positive equilibrium

by constructing the Lyapunov function

V(x1, x2, y)

=k1x1−x∗∗

1−x∗∗

1ln x1

x∗∗

1

+k2x2−x∗∗

2−x∗∗

2ln x2

x∗∗

2

+k3y−y∗∗ −y∗∗ ln y

y∗∗ ,

where x∗∗

1, x∗∗

2, y∗∗ are defined in (15). However, it

is very difficult to deal with the nonlinear term to en-

sure D+V(t)<0. By constructing Lyapunov func-

tion as above, one could obtain some sufficient con-

ditions to ensure the global stability of the solution to

system (1), however, the condition, generally speak-

ing, is not a good one, since in dealing with the term

η1+η2yx2, some additional conditions are needed,

which could not reflect the essential characteristic of

system (1). We mention here that in [7], Wu investi-

gated the global stability of the equilibrium of system

(9) by using the Dulac criterion, which could only be

applied to the two dimensional systems, and could not

be applied to three dimensional systems. Recently, in

the study of dynamic behaviors of the system (10),

Chen, Chong and Lin [31] developed some analytical

technique, more precisely, they essentially grasped

the characteristics of the commensalism systems, and

applied the differential inequalities theory, to obtain

some interesting results about system (10). We will

try to develop the analytic idea of [31] to investigate

the dynamic behaviors of system (1).

The paper is arranged as follows. We investigate

the existence and locally stability property of the equi-

libria of system (1) in Section 2. In Section 3, we

establish a stability result about the single species

stage structured system. In Section 4, by using the

differential inequalities theory and the Lemma estab-

lished in Section 3, we provide conditions which en-

sure the global attractivity property of the equilibria.

In Section 5, we present some numerical simulations

to show the feasibility of the main results. We end this

paper by a brief discussion.

2 Local stability

In this section, we will investigate the existence and

local stability property of system (1).

The equilibria of system (1) is determined by the

following system

αx2−βx1−δ1x1= 0,

βx1−δ2x2−γx2

2+dy

η1+η2yx2= 0,

y(b2−a2y) = 0.

(11)

The system always admits two boundary equilibria:

A1(0,0,0),A2(0,0,b2

), if, in addition

αβ > δ2(β+δ1),(12)

then the system admits another boundary equilibrium

A3x∗

1, x∗

2,0, where

x∗

1=αx∗

β+δ1

, x∗

2=αβ −δ2(β+δ1)

γ(β+δ1).(13)

αβ − δ2−dy∗∗

η1+η2y∗∗ !(β+δ1)>0,(14)

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

813

Volume 21, 2022

then system (1) admits a unique positive equilibrium

A4x∗∗

1, x∗∗

2, y∗∗ , where

x∗∗

1=αx∗∗

β+δ1

x∗∗

αβ − δ2−dy∗∗

η1+η2y∗∗ !(β+δ1)

(β+δ1)γ,

y∗∗ =b2

(15)

Obviously, x∗∗

1, x∗∗

2and y∗∗ satisfies the equations

αx∗∗

2−βx∗∗

1−δ1x∗∗

1= 0,

βx∗∗

1−δ2x∗∗

2−γ(x∗∗

2)2+dy∗∗

η1+η2y∗∗ x∗∗

2= 0,

b2−a2y∗∗ = 0.

(16)

We shall now investigate the local stability prop-

erty of the above equilibria.

The variational matrix of system (1) is

J(x1, x2, y)

=





−β−δ1α0

β W1W2

0 0 −2a2y+b2







where

W1=−δ2−2γx2+dy

η1+η2y,

W2=dx2

η1+η2y−dη2x2y

(η1+η2y)2.

Theorem 2.1 A1(0,0,0) is unstable.

Proof. The Jacobian matrix of the equilibrium point

A1(0,0,0) is given by







−β−δ1α0

β−δ20

0 0 b2







The characteristic equation of the above matrix is

(λ−b2)λ2+(δ1+δ2+β)λ+βδ2+δ1δ2−αβ= 0.

Hence, it has one positive characteristic root λ1=b2,

consequently, A1(0,0,0) is unstable. This ends the

proof of Theorem 2.1.

Theorem 2.2 Assume that

(β+δ1) δ2−dy∗∗

η1+η2y∗∗ !−αβ > 0,(17)

then A2(0,0,b2

a2)is locally asymptotically stable. As-

sume that

(β+δ1) δ2−dy∗∗

η1+η2y∗∗ !−αβ < 0,(18)

then A2(0,0,b2

a2)is unstable.

Proof. The Jacobian matrix of the system about the

equilibrium point A2(0,0,b2

a2)is given by







−β−δ1α0

βdy∗∗

η1+η2y∗∗ −δ20

0 0 −b2







The characteristic equation of the above matrix is

(λ+b2)λ2+A1λ+A2= 0.(19)

where

A1=δ1+δ2+β−dy∗∗

η1+η2y∗∗ ,

A2= (β+δ1)δ2−dy∗∗

η1+η2y∗∗ −αβ.

Hence, it has one negative characteristic root λ1=

−b2<0, the other two characteristic roots are deter-

mined by the equation

λ2+A1λ+A2= 0.(20)

Noting that the two characteristic roots of equation

(20) satisfy

λ2+λ3=−A1,

λ2λ3=A2,

(21)

under the assumption (18), λ2λ3<0,hence, λ2

and λ3should be one positive and the other negative,

which means that one characteristic root is positive,

consequently, A2(0,0,b2

a2)is unstable. Instead, as-

sumption (17) implies that

δ2>dy∗∗

η1+η2y∗∗ ,

and so, from (21) and (17), one has λ2+λ3<

0, λ2λ3>0.Hence, λ2<0, λ3<0.That is, un-

der the assumption (17), three characteristic roots of

the equation (19) are all negative, hence, A1(0,0,b2

a2)

is locally asymptotically stable. This ends the proof

of Theorem 2.2.

Theorem 2.3 A3(x∗

1, x∗

2,0) is unstable.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

814

Volume 21, 2022

Proof. Obviously, x∗

1, x∗

2(see (13)) satisfy the equa-

tions

αx∗

2−βx∗

1−δ1x∗

1= 0,

βx∗

1−δ2x∗

2−γ(x∗

2)2= 0.

Consequently, the Jacobian matrix at the equilibrium

point A3(x∗

1, x∗

2,0) is given by







−β−δ1α0

β−δ2−2γx∗

dx∗

η1

0 0 b2







The characteristic equation of the above matrix is

(λ−b2) λ2+B1λ+B2!= 0,(22)

where

B1=δ1+ 2γx∗

2+β+δ2,

B2= (β+δ1)(δ2+ 2γx∗

2)−αβ.

Equation (22) has at least one positive root λ1=b2,

consequently, A3(x∗

1, x∗

2,0) is unstable. This ends

the proof of Theorem 2.3.

Theorem 2.4 Assume that (14) holds, then

A4(x∗∗

1, x∗∗

2, y∗∗ ), where x∗∗

1, x∗∗

2, y∗∗ are expressed

in (15), is locally asymptotically stable.

Proof. The Jacobian matrix of the system about the

equilibrium point A4(x∗∗

1, x∗∗

2, y∗∗ )is given by







−β−δ1α0

β V1V2

0 0 −2a2y∗∗ +b2







where

V1=−δ2−2γx∗∗

2+dy∗∗

η1+η2y∗∗ ,

V2=dx∗∗

η1+η2y∗∗ −dη2x∗∗

(η1+η2y∗∗)2.

Noting that

−2a2y∗∗ +b2=−2a2

+b2=−b2,

and, from the second equation of (16), the definitions

of x∗∗

1and x∗∗

2(see (15)), we have

−δ2−2γx∗∗

2+dy∗∗

η1+η2y∗∗

=−βx∗∗

x∗∗

−γx∗∗

2=−αβ

β+δ1

−γx∗∗

So, the characteristic equation of above matrix is

λ+b2hλ2+C1λ+C2i= 0,

where

C1=β+δ1+βα

β+δ1

+γx∗∗

C2= (β+δ1)βα

β+δ1

+γx∗∗

2−αβ.

Hence, it has one negative characteristic root λ1=

−b2<0, the other two characteristic roots are deter-

mined by the equation

λ2+C1λ+C2= 0.(23)

Noting that from the expression of x∗∗

2(see (15)) and

condition (14), the two characteristic roots of equation

(23) satisfy

λ2+λ3

=−C1<0,

λ2λ3

= (β+δ1)βα

β+δ1

+γx∗∗

2−αβ

= (β+δ1)βα

β+δ1

αβ −δ2−dy∗∗

η1+η2y∗∗ (β+δ1)

β+δ1







−αβ

=αβ − δ2−dy∗∗

η1+η2y∗∗ !(β+δ1)>0,

and so, λ2<0, λ3<0.Therefore, all of the

three characteristic roots are negative, consequently,

A4(x∗∗

1, x∗∗

2, y∗∗ )is locally asymptotically stable.

This ends the proof of Theorem 2.4.

Remark 2.1. System (1) admits four equilibria,

moreover, the local stability property of this equilib-

ria is similar to the equilibria of system (5). Assume

that the inequality (17) holds, then A2(0,0,b2

a2)

is locally asymptotically stable, and assume that

(14) holds, then the positive equilibrium is locally

asymptotically stable. However, A1(0,0,0) and

A3(x∗

1, x∗

2,0) are all unstable, which means that the

second species could not be driven to extinction.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

815

Volume 21, 2022

Remark 2.2. Theorem 2.4 shows that if the positive

equilibrium exists, it’s locally asymptotically stable.

Remark 2.3. Noting that from (13) and (15), one

could easily see that x∗∗

2> x∗

2,x∗∗

1> x∗

1, that is,

commensalism increases the final density of the first

species.

Remark 2.4. When η1= 1, η2= 0,system (1) is

degenerate to system (5), and Theorems 2.1-2.4 is de-

generate to Theorems 2.1-2.4 of Lei [13], thus, we

generalize the main results of Lei [13] to the nonlinear

case.

3 Dynamic Behaviors of Single

Species Stage System with

Perturbation

This section will focus on the dynamic behaviors of a

single species stage structured system.

Now let’s consider the system

dx1

dt =αx2−βx1−δ1x1,

dx2

dt =βx1−δ2x2−γx2

(24)

where α, β, δ1, δ2and γare all positive constants,

x1(t)and x2(t)are the densities of the immature and

mature members of the species at time t. One could

refer to system (1) for more detailed explanation of

the biological meaning of those coefficients. From

Theorem 4.1 and 4.2 in [32] by Xiao and Lei, we

have

Lemma 3.1. Assume that

αβ < δ2β+δ1(25)

holds, then the boundary equilibrium O(0,0) of sys-

tem (24) is globally stable. Assume that

αβ > δ2β+δ1(26)

holds, then the positive equilibrium B(x∗

1, x∗

2)of sys-

tem (24) is globally stable, where x∗

1and x∗

2are de-

fined by (13).

Now let’s consider the single species system with

linear perturbation

dx1

dt =αx2−βx1−δ1x1,

dx2

dt =βx1−δ2x2−γx2

2+εx2,

(27)

where εis positive constant. Obviously, if ε < δ2, the

dynamic behaviors of the system (27) is similar to that

of the system (24). However, if ε≥δ2, then we could

not directly apply Lemma 3.1 to system (27). In this

case, set δ0

2=ε−δ2, system (27) could be rewritten

as follows:

dx1

dt =αx2−βx1−δ1x1,

dx2

dt =βx1+δ0

2x2−γx2

(28)

whether system (28) has the similar dynamic behav-

iors as that of system (24) is still unknown. In (28),

δ0

2can not be explained as the death rate of the mature

species, indeed it is a sum of the death rate of mature

species and the perturbation coefficient, the dynamic

behaviors of this system, to the best of our knowledge,

is not investigated by the scholars. So, for the sake

of completeness, we give here the complete proof of

Lemma 3.2 below.

By simple computation, the system admits a posi-

tive equilibrium B2(x∗

1, x∗

2), where

x∗

1=α(αβ +βδ0

2+δ1δ0

γ(β+δ1)2, x∗

2=αβ +βδ0

2+δ1δ0

γ(β+δ1).

Concerned with the stability property of this posi-

tive equilibrium, we have the following result.

Lemma 3.2. The positive equilibrium B2(x∗

1, x∗

2)of

system (28) is globally asymptotically stable.

Proof. We will prove Lemma 3.2 by adapting the idea

of Xiao and Lei [32]. More precisely, we will prove

Lemma 3.2 by constructing a suitable Lyapunov func-

tion as follows. Set

V(x1, x2) = k1x1−x∗

1−x∗

1ln x1

x∗

1

+k2x2−x∗

2−x∗

2ln x2

x∗

2,

where k1, k2are some positive constants determined

later.

One could easily see that the function Vis zero at

the equilibrium B2(x∗

1, x∗

2)and is positive for all other

positive values of x1and x2. The time derivative of

Valong the trajectories of (28) is

D+V(t)

=k1

x1−x∗

x1αx2−(β+δ1)x1

+k2

x2−x∗

x2βx1+δ0

2x2−γx2

2.

(29)

Since αx∗

2−βx∗

1−δ1x∗

1= 0,

βx∗

1+δ0

2x∗

2−γ(x∗

2)2= 0,

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

816

Volume 21, 2022

then

αx2−(β+δ1)x1

=α

x∗

1x2x∗

1−x1x∗

2+αx1

x∗

−β+δ1x1

=α

x∗

1−x2(x1−x∗

1) + x1(x2−x∗

2),

(30)

and

βx1+δ0

2x2−γx2

=β

x∗

2x1x∗

2−x2x∗

1+βx2

x∗

+δ0

2x2−γx2

=β

x∗

2x1x∗

2−x1x2+x1x2−x2x∗

1

+αβ

β+δ1

+δ0

2x2−γx2

=β

x∗

2x1(x∗

2−x2) + x2(x1−x∗

1)

+γx∗

2x2−γx2

(31)

Applying (30) and (31) to (29), and choosing k2=

1, k1=βx∗

x∗

2α,we finally obtain

D+V(t)

=k1

x1−x∗

x∗

1−x2(x1−x∗

1) + x1(x2−x∗

2)

+k2

x2−x∗

x∗

2x1(x∗

2−x2) + x2(x1−x∗

1)

+k2

x2−x∗

x2γx∗

2x2−γx2

2

=−β

x∗

2hrx2

(x1−x∗

1)−rx1

(x2−x∗

2)i2

−γx2−x∗

22.

that is, D+V(t)<0strictly for all x1, x2>0with

the exception of the positive equilibrium B2(x∗

1, x∗

2),

where D+V(t) = 0. Thus, V(x1, x2)satisfies Lya-

punov’s asymptotic stability theorem, and the posi-

tive equilibrium B2(x∗

1, x∗

2)of system (28) is globally

asymptotically stable.

This completes the proof of Lemma 3.2.

4 Global attractivity

We had shown in Section 2 that A2(0,0,b2

a2)and

A4(x∗∗

1, x∗∗

2, y∗∗ )could be locally asymptotically

stable under some suitable assumptions. One in-

teresting issue is to investigate the global stability

property of the equilibria. In this section we will

try to obtain some sufficient conditions to ensure

the global attractivity of the equilibria A2and A4of

system (1).

Theorem 4.1 Assume that

(β+δ1) δ2−dy∗∗

η1+η2y∗∗ !−αβ > 0(32)

holds, then A2(0,0,b2

a2)is globally attractive.

Proof. For ε > 0enough small, condition (32) im-

plies that

(β+δ1)δ2−dy∗∗ +ε)

η1+η2(y∗∗ +ε)−αβ > 0.(33)

Since the third equation of system (1) is a Logistic

equation, then

lim

t→+∞

y(t) = b2

.(34)

In corresponding to the above chosen ε > 0, there

exists a T > 0such that

y(t)<b2

+εdef

=y∗∗ +εfor all t > T.

From above inequality and the first and second equa-

tions of system (1), also, noting that the function

η1+η2yis monotonically increasing, for t > T , we

have

dx1

dt =αx2−βx1−δ1x1,

dx2

dt =βx1−δ2x2−γx2

2+dy

η1+η2yx2

≤βx1−δ2x2−γx2

+dy∗∗ +ε

η1+η2(y∗∗ +ε)x2.

(35)

Now let’s consider the system

dv1

dt =αv2−βv1−δ1v1,

dv2

dt =βv1−δ2v2−γv2

+dy∗∗ +ε

η1+η2(y∗∗ +ε)v2.

(36)

It follows from (33) and Lemma 3.1 that the boundary

equilibrium O(0,0) of system (36) is globally stable.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

817

Volume 21, 2022

That is, for any positive solution (v1(t), v2(t)) of the

system (36), one has

lim

t→+∞

v1(t) = 0,lim

t→+∞

v2(t) = 0.

Let (x1(t), x2(t), x3(t)) be any positive so-

lution of system (1) with initial condition

(x1(T), x2(T), y(T)) = (x10, x20, y0), and let

(v1(t), v2(t)) be the positive solution of system

(36) with the initial condition (v1(T), v2(T)) =

(x10, x20), then it follows from (35), (36) and the

differential inequality theory that

xi(t)≤vi(t)for all t≥T, i = 1,2,

and therefore that

lim sup

t→+∞

xi(t)≤lim

t→+∞

vi(t) = 0, i = 1,2.

Thus, from the positivity of the solution of system (1),

it immediately follows that

0≤lim inf

t→+∞

xi(t)≤lim sup

t→+∞

xi(t)≤0, i = 1,2.

Therefore

lim

t→+∞

xi(t) = 0, i = 1,2.(37)

(34) together with (37) shows that A2(0,0,b2

a2)is

globally attractive.

This completes the proof of Theorem 4.1.

Theorem 4.2 Assume that

αβ − δ2−dy∗∗

η1+η2y∗∗ !(β+δ1)>0(38)

holds, then A4(x∗∗

1, x∗∗

2, y∗∗ ), where x∗∗

1, x∗∗

2, y∗∗ are

defined in (15), is globally attractive.

Proof. For ε > 0enough small, without loss of gener-

ality, we may assume that ε < b2

2a2, so condition (38)

implies that

αβ − δ2−d(y∗∗ +ε)

η1+η2(y∗∗ +ε)!(β+δ1)>0(39)

and

αβ − δ2−d(y∗∗ −ε)

η1+η2(y∗∗ −ε!(β+δ1)>0(40)

hold. The third equation of system (1) is a Logistic

equation, thus

lim

t→+∞

y(t) = b2

.(41)

In correspondence to the above chosen ε > 0, there

exists a T > 0such that

−ε < y(t)<b2

+εfor all t > T, i = 1,2.

(42)

From (42) and the first and second equation of system

(1), also, noting that the function y

η1+η2yis strictly

increasing, for t > T , we have

dx1

dt =αx2−βx1−δ1x1,

dx2

dt ≤βx1−δ2x2−γx2

+dy∗∗ +ε

η1+y∗∗ +εx2.

(43)

Now let’s consider the system

dv1

dt =αv2−βv1−δ1v1,

dv2

dt =βv1− δ2−dy∗∗ +ε

η1+η2(y∗∗ +ε)!v2

−γv2

(44)

There are two subcases.

(i)

δ2> d y∗∗ +ε

η1+η2(y∗∗ +ε).

In this case, it follows from (39) and Lemma 3.1 that

the system (44) admits a unique positive equilibrium

which is globally asymptotically stable.

(ii)

δ2≤d

+ε

η+b2

+ε

In this case, it follows from Lemma 3.2 that the sys-

tem (44) admits a unique positive equilibrium which

is globally stable.

Hence, in any case, system (44) admits a unique

positive equilibrium (v1(ε), v2(ε)), which is globally

asymptotically stable, where

v1(ε) = αv2(ε)

β+δ1

v2(ε) =

αβ − δ2−dy∗∗ +ε

η1+η2(y∗∗ +ε)!(β+δ1)

(β+δ1)γ.

(45)

Therefore, let (v1(t), v2(t)) be any positive solution

of the system (44), one has

lim

t→+∞

v1(t) = v1(ε),lim

t→+∞

v2(t) = v2(ε).(46)

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

818

Volume 21, 2022

Let (x1(t), x2(t), x3(t)) be any positive so-

lution of system (1) with initial condition

(x1(T), x2(T), y(T)) = (x10, x20, y0), and let

(v1(t), v2(t)) be the positive solution of system

(44) with the initial condition (v1(T), v2(T)) =

(x10, x20), it then follows from (43), (44) and the

differential inequality theory that

xi(t)≤vi(t)for all t≥T, i = 1,2.(47)

The positivity of the solution of system (1) and (47)

lead to

lim sup

t→+∞

xi(t)≤lim

t→+∞

vi(t) = vi(ε), i = 1,2.(48)

From (42) and the first and second equation of system

(1), since y

η1+η2yis strictly increasing, for t > T ,

we have

dx1

dt =αx2−βx1−δ1x1,

dx2

dt ≥βx1−δ2x2−γx2

+dy∗∗ −ε

η1+η2(y∗∗ −ε)x2.

(49)

Now let’s consider the system

dw1

dt =αw2−βw1−δ1w1,

dw2

dt =βw1− δ2−dy∗∗ −ε

η1+η2(y∗∗ −ε)!w2

−γw2

(50)

There are two subcases.

(i)

δ2> d y∗∗ −ε

η1+η2(y∗∗ −ε).(51)

In this case, it follows from (40) and Lemma 3.1 that

the system (50) admits a unique positive equilibrium

which is globally asymptotically stable.

(ii)

δ2≤dy∗∗ −ε

η1+η2(y∗∗ −ε).

In this case, it follows from Lemma 3.2 that the sys-

tem (50) admits a unique positive equilibrium which

is globally stable.

Hence, in any case, system (50) admits a unique

positive equilibrium (w1(ε), w2(ε)), which is glob-

ally asymptotically stable, where

w1(ε) = αw2(ε)

β+δ1

w2(ε) =

αβ − δ2−dy∗∗ −ε

η1+η2(y∗∗ −ε)!(β+δ1)

(β+δ1)γ.

(52)

Therefore, let (w1(t), w2(t)) be any positive solution

of the system (50), one has

lim

t→+∞

w1(t) = w1(ε),lim

t→+∞

w2(t) = w2(ε).(53)

Let (x1(t), x2(t), y(t)) be any positive so-

lution of system (1) with initial condition

(x1(T), x2(T), y(T)) = (x10, x20, y0), and

let (w1(t), w2(t)) be the positive solution

of system (50) with the initial condition

(w1(T), w2(T)) = (x10, x20), it then follows

from (49), (50) and the differential inequality theory

that

xi(t)≥wi(t)for all t≥T, i = 1,2.(54)

The positivity of the solution of system (1) and (53),

(54) lead to

lim inf

t→+∞

xi(t)≥lim

t→+∞

wi(t) = wi(ε), i = 1,2.

(55)

Inequality (48) together with (55) leads to

wi(ε) = lim

t→+∞

wi(t)

≤lim inf

t→+∞

xi(t)

≤lim sup

t→+∞

xi(t)

≤lim

t→+∞

vi(t) = vi(ε), i = 1,2.

(56)

From (45) and (52) one can easily see that

wi(ε)→x∗∗

i, vi(ε)→x∗∗

i,as ε →0, i = 1,2.

(57)

Noting that εis any enough small positive constant,

letting ε→0in (56), (57) leads to

lim

t→+∞

xi(t) = x∗∗

i, i = 1,2.(58)

So, (41) together with (58) shows that

A4(x∗∗

1, x∗∗

2, y∗∗ )is globally attractive.

This completes the proof of Theorem 4.2.

Remark 4.1. Theorem 4.1 and 4.2 depicts a very intu-

itive biological phenomenon. From Zhang, Chen and

Neumann [40], for single species stage structured sys-

tem (24), we can regard α

δ2

as a relative birth rate of

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

819

Volume 21, 2022

the fist mature species, β

β+δ1

as a relative transfor-

mation rate of the first immature species. Then con-

dition (32) and (38) are equivalent to

δ2−dy∗∗

η1+η2y∗∗

β+δ1

<1,(59)

and α

δ2−dy∗∗

η1+η2y∗∗

β+δ1

>1,(60)

respectively. Hence, with the help of the host species,

the relative birth rate of the mature commensal

species is increasing from α

δ2

to α

δ2−dy∗∗

η1+η2y∗∗

thus finally increasing the chance of the survival of

the first species.

Remark 4.2. When η1= 1, η2= 0,system (1) is de-

generate to system (5), and Theorem 4.1-4.2 are de-

generate to Theorem 3.1-3.2 of Lei[13], respectively.

Thus, we generalize the main results of Lei [13] to the

nonlinear case.

5 Examples

Example 5.1. Consider the following stage structured

commensalism system

dx1

dt = 3x2−x1−x1,

dx2

dt =x1−x2−x2

2+dy

η1+η2yx2,

dt =y(1 −y).

(61)

Here, corresponding to system (1), we take α=

3, β =δ1=δ2=γ=b2=a2= 1.Note that

αβ = 3 >2 = δ2(β+δ1).(62)

In this case, without the commensalism of the second

species, the first species is globally stable. For any

positive constant d,η1and η2, inequality (38) holds,

so it follows from Theorem 4.2 that A4(x∗∗

1, x∗∗

2, y∗∗ )

is globally attractive. Take d=η1=η2= 1,

then the positive equilibrium A4(3

2,1,1) is globally

attractive, noting that lim

t→+∞

y(t) = 1, hence, we only

give the numeric simulations of immature species

x1and mature species x2. Also, here we fixed

x1(0) + x2(0) = 2.5, y(0) = 2. Fig. 1 shows that

lim

t→+∞

x1(t) = 3

2, though the solutions are finally

approaching to 3

2, moreover for the initial conditions

x1(0) <3

2, the solution is strictly increasing at first,

then, after it reaches the extremum, which is larger

than 3

2, after that, the solutions becomes strictly

decreasing. On the contrary, for the initial conditions

x1(0) >3

2, the solution is strictly decreasing at first,

then, after it reach the extremum, which is smaller

than 3

2, then it becomes strictly increasing. Fig.2

shows that lim

t→+∞

x2(t) = 1. However, the dynamic

behaviors is relative simple, for solutions with initial

conditions larger than 1.5, the solutions is strictly

decreasing, and for solutions with initial conditions

smaller than 1, the solutions is strictly increasing.

Figure 1: Dynamic behaviors of the first com-

ponents x1in system (61) with the initial condi-

tion (x1(0), x2(0), y(0)) = (1,1.5,2),(0.5,2,2),

(1.5,1,2) and (2,0.5,2), respectively.

Figure 2: Dynamic behaviors of the second com-

ponents x2in system (61) with the initial condi-

tion (x1(0), x2(0), y(0)) = (1,1.5,2),(0.5,2,2),

(1.5,1,2) and (2,0.5,2), respectively.

Example 5.2. Consider the following stage structured

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

820

Volume 21, 2022

commensalism system

dx1

dt =x2−x1−x1,

dx2

dt =x1−x2−x2

2+3y

η1+η2yx2,

dt =y(1 −y).

(63)

Here, corresponding to system (1), we take α=β=

δ1=δ2=γ=b2=a2= 1, d = 3.For

η1= 1, η2= 0, the system degenerates to the case 2

in system 4.1 of Lei [13]. From [13], without the com-

mensal of the second species, the first species, which

satisfies the equations

dx1

dt =x2−x1−x1,

dx2

dt =x1−x2−x2

(64)

will be driven to extinction. Also, with the commen-

sal of the second species, the system admits a unique

positive equilibrium A4(3

4,3

2,1), which is globally

asymptotically stable.

Noting that from the third equation of (63), we

have y∗∗ = 1. By simple computation, if

(β+δ1)δ2−3

η1+η2−αβ

= 2 1−3

η1+η2−1>0,

(65)

which is equivalent to η1+η2>6holds true,

then it follows from Theorem 4.1 that A2(0,0,1)

is globally attractive. Let’s take η1= 9, η2= 1.

Numeric simulations reported in Fig. 3 and 4 sup-

ports this assertion. By a similarly discussion, if

η1+η2<6holds, then it follows from Theorem

4.2 that A4(x∗∗

1, x∗∗

2,1) is globally attractive. Let’s

take η1= 1, η2= 1, then A4(0.5,1,1) is globally

attractive. Numeric simulations reported in Fig. 5

and 6 supports this assertion. In Fig. 3-6, we choose

the initial conditions x1(0) + x2(0) = 2, y(0) = 1.

6 Conclusion

Stimulated by the work of Wright [41], we argued that

in the commensal relationship between two species,

the commensal species is also subjected to the con-

straints of handling time. For example, consider a

species of solitary bee visiting a flower species. The

rate of collection of pollen is limited by the handling

time per plant. This finally motivated us to propose a

two species commensalism model with Holling type

II functional response and stage structure. If η1=

Figure 3: Dynamic behaviors of the first com-

ponent x1in system (64) with η1= 9, η2= 1,

and the initial condition (x1(0), x2(0), y(0)) =

(1,1,1),(1.8,0.2,1) and (0.2,1.8,1), respec-

tively.

Figure 4: Dynamic behaviors of the second com-

ponent x2in system (64) with η1= 9, η2= 1,

and the initial condition (x1(0), x2(0), y(0)) =

(1,1,1),(1.8,0.2,1) and (0.2,1.8,1), respec-

tively.

Figure 5: Dynamic behaviors of the first com-

ponent x1in system (64) with η1= 1, η2= 1,

and the initial condition (x1(0), x2(0), y(0)) =

(1,1,1),(1.8,0.2,1) and (0.2,1.8,1), respec-

tively.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

821

Volume 21, 2022

Figure 6: Dynamic behaviors of the second com-

ponent x2in system (64) with η1= 1, η2= 1,

and the initial condition (x1(0), x2(0), y(0)) =

(1,1,1),(1.8,0.2,1) and (0.2,1.8,1), respec-

tively.

1, η2= 0, then our model is degenerate to the model

considered by Lei [13].

We first showed that the system can admit four

equilibria, however, only two of them can be locally

asymptotically stable.

From Lemma 3.1 we know that for the stage struc-

tured single species system, depending on the rela-

tionship of the coefficients, the species may be driven

to extinction or survival in the long run. Theorem 4.2

shows that if the species without commensalism of the

second species could survive in the long run, then for

the commensalism system, two species could coex-

ist in a stable state. If the first species will be driven

to extinction without the commensalism, then, Theo-

rem 4.1 shows that limited commensalism still could

not avoid the extinction of the first species. However,

Theorem 4.2 shows that if the commensalism effect is

large enough, then two species can coexist in a stable

state.

One can easily see that if η1= 1, η2= 0, then

Theorems 2.1-2.4 degenerate to Theorems 2.1-2.4,

Theorems 4.1-4.2 degenerate to Theorems 3.1-3.2 in

Lei [13], hence, we generalize the main result of

Lei [13].

What we really concern is the influence of Holling

type II response, which, from the point of Wright [43],

can be explained as the handling time of commensal

time. Example 5.2 shows that handling time has nega-

tive effect on the persistent property of the commensal

species. If η1+η2y∗∗ is large enough, the influence

of the host will be reduced, and with the influence

of stage structure, despite the commensalism of host

species, the commensal species will still be driven to

extinction.

We mention here that the following two aspects

need to be studied. The first one is about the influ-

ence of delay. A more plausible model should include

some past state of the system [28], this leads to the

delayed modelling. A recent study [28] shows that

delay could change the dynamic behaviors of the sys-

tem greatly. Up today still no scholars propose and

study the delayed stage structured commensal model,

we will try to do some work in this direction. The sec-

ond one is to investigate the influence of Allee effect.

Already, several scholars [14], [15], [16], [17], [34],

[35] studied the commensalism system with Allee ef-

fect, however, all of those works did not consider the

stage structure of the species. We will try to do some

works on this direction in the future.

References:

[1] Chen F. D., Xie X. D. and Chen X. F. , Dy-

namic behaviors of a stage-structured coopera-

tion model, Commun. Math. Biol. Neurosci., Vol.

2015, 2015, 19 pages.

[2] Zhu Z., Chen F., Lai L., et al, Dynamic behav-

iors of a discrete May type cooperative system

incorporating Michaelis-Menten type harvesting,

IAENG Int. J. Appl. Math. Vol. 50, No.1, 2020,

pp. 1-10.

[3] Xie X. D., Chen F. D., Xue Y. L., Note on

the stability property of a cooperative system in-

corporating harvesting, Discrete Dyn. Nat. Soc.,

Vol.2014, No.1, 2014, 5 pages.

[4] Xue Y. L., Chen F. D. , Xie X. D. , et al. Dynamic

behaviors of a discrete commensalism system,

Annals of Applied Mathematics, Vol.31, No.4,

2015, pp. 452-461.

[5] Zhu Z. , Wu R., et al, Dynamic behaviors of a

Lotka-Volterra commensal symbiosis model with

non-selective Michaelis-Menten type harvesting,

IAENG Int. J. Appl. Math. Vol.50, No.1, 2020,

pp.396-404.

[6] Miao Z. S. , Xie X. D., Pu L. Q., Dynamic be-

haviors of a periodic Lotka-Volterra commensal

symbiosis model with impulsive, Commun. Math.

Biol. Neurosci. Vol.2015, No.1, 2015, 15 pages.

[7] Wu R. X. , Lin L. , Zhou X. Y., A commensal

symbiosis model with Holling type functional re-

sponse, J. Math. Computer Sci., Vol.16, No.1,

2016, pp.364-371.

[8] Xie X. D., Miao Z. S., Xue Y. L. , Positive pe-

riodic solution of a discrete Lotka-Volterra com-

mensal symbiosis model, Commun. Math. Biol.

Neurosci. Vol.2015, No.1, 2015, 10 pages.

[9] Chen B. , The influence of commensalism on a

Lotka-Volterra commensal symbiosis model with

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

822

Volume 21, 2022

Michaelis-Menten type harvesting, Advances in

Difference Equations, Vol. 2019, No.1, 2019, Ar-

ticle ID 43.

[10] Liu Y., Xie X., Lin Q., Permanence, par-

tial survival, extinction, and global attractivity

of a nonautonomous harvesting Lotka-Volterra

commensalism model incorporating partial clo-

sure for the populations, Advances in Difference

Equations, 2018, Article ID 211.

[11] Deng H., Huang X., The influence of partial clo-

sure for the populations to a harvesting Lotka-

Volterra commensalism model, Commun. Math.

Biol. Neurosci., Vol.2018, No.1, 2018: Article ID

10.

[12] Xue Y., Xie X., Lin Q., Almost periodic solu-

tions of a commensalism system with Michaelis-

Menten type harvesting on time scales, Open

Mathematics, Vol.17, No. 1, 2019, 1503-1514.

[13] Lei C., Dynamic behaviors of a stage-structured

commensalism system, Adv. Differ. Equ., Vol.

2018, 2018, Article ID 301.

[14] Lin Q., Allee effect increasing the final den-

sity of the species subject to the Allee effect in

a Lotka-Volterra commensal symbiosis model,

Adv. Differ. Equ., Vol. 2018, 2018, Article ID 196.

[15] Chen B., Dynamic behaviors of a commensal

symbiosis model involving Allee effect and one

party can not survive independently, Adv. Differ.

Equ. Vol.2018, 2018, Article ID 212.

[16] Wu R., Li L., Lin Q., A Holling type commen-

sal symbiosis model involving Allee effect, Com-

mun. Math. Biol. Neurosci., Vol.2018, 2018, Ar-

ticle ID 6.

[17] Lei C., Dynamic behaviors of a Holling type

commensal symbiosis model with the first species

subject to Allee effect, Commun. Math. Biol.

Neurosci., Vol.23, No. 1, 2019, Article ID 3.

[18] Vargas-De-León C. , Gómez-Alcaraz G., Global

stability in some ecological models of commen-

salism between two species, Biomatemática,

Vol.23, No.1, 2013, pp. 139-146.

[19] Chen F., Xue Y., Lin Q., et al, Dynamic be-

haviors of a Lotka-Volterra commensal symbio-

sis model with density dependent birth rate, Adv.

Differ. Equ. Vol.2018, 2018, Article ID 296.

[20] Han R., Chen F., Global stability of a commensal

symbiosis model with feedback controls, Com-

mun. Math. Biol. Neurosci., Vol.2015, 2015, Ar-

ticle ID 15.

[21] Chen F., Pu L., Yang L., Positive periodic solu-

tion of a discrete obligate Lotka-Volterra model,

Commun. Math. Biol. Neurosci., Vol.2015, 2015:

Article ID 14.

[22] Guan X., Chen F., Dynamical analysis of a

two species amensalism model with Beddington-

DeAngelis functional response and Allee effect

on the second species, Nonlinear Analysis: Real

World Applications, Vol. 48, No.1, 2019, pp.71-

93.

[23] Li T., Lin Q., Chen J., Positive periodic solution

of a discrete commensal symbiosis model with

Holling II functional response, Commun. Math.

Biol. Neurosci., Vol.2016, No.1, 2016, Article ID

22.

[24] Ji W., Liu M., Optimal harvesting of a stochas-

tic commensalism model with time delay, Phys-

ica A: Statistical Mechanics and its Applications,

Vol.527, No.1, 2019, Article ID: 121284.

[25] Puspitasari N., Kusumawinahyu W. M.,

Trisilowati T., Dynamic analysis of the symbiotic

model of commensalism and parasitism with

harvesting in commensal populations, JTAM

(Jurnal Teori dan Aplikasi Matematika), Vol. 5,

No. 1, 2021, pp. 193-204.

[26] Jawad S., Study the dynamics of commensalism

interaction with Michaels-Menten type prey har-

vesting, Al-Nahrain Journal of Science, Vol.25,

No. 1, 2022, pp. 45-50.

[27] Kumar G. B., Srinivas M.N., Influence of spa-

tiotemporal and noise on dynamics of a two

species commensalism model with optimal har-

vesting, Research Journal of Pharmacy and

Technology, Vol. 9, No. 10, 2016, pp. 1717-1726.

[28] Li T., Wang Q., Stability and Hopf bifurcation

analysis for a two-species commensalism system

with delay, Qualitative Theory of Dynamical Sys-

tems, Vol. 20, No.3, 2021, pp. 1-20.

[29] Chen L., Liu T., et al, Stability and bifurcation

in a two-patch model with additive Allee effect,

AIMS Mathematics, Vol. 7, No. 1, 2022, pp. 536-

551.

[30] Zhu Z. , Chen Y., et al. Stability and bifurca-

tion in a Leslie-Gower predator-prey model with

Allee effect, International Journal of Bifurca-

tion and Chaos, Vol. 32, No. 03, 2022, Artical ID:

2250040.

[31] Chen F., Chong Y., Lin S., Global stability

of a commensal symbiosis model with Holling

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

823

Volume 21, 2022

II functional response and feedback controls,

Wseas Trans. Syst. Contr. Vol. 17, No. 1., 2022,

pp.279–286.

[32] Xiao A., Lei C. Q., Dynamic behaviors of a non-

selective harvesting single species stage structure

system incorporating partial closure for the pop-

ulations, Advances in Difference Equations, Vol.

2018, 2018, Article ID 245.

[33] Yue Q., Stability property of the prey free equi-

librium point, Open Mathematics, 17 (1) (2019)

646-652.

[34] Wei Z., Xia Y., Zhang T., Stability and bifurca-

tion analysis of a commensal model with additive

Allee effect and nonlinear growth rate, Internat.

J. Bifur. Chaos Vo. 31, No.1, Artical ID: 2150204.

[35] He X., Zhu Z., Chen J., et al., Dynamical analy-

sis of a Lotka Volterra commensalism model with

additive Allee effect, Open Mathematics. Vol.20,

No.1, 2022, pp. 646-665

[36] Zhou Q., Lin S., Chen F., Wu R., Positive pe-

riodic solution of a discrete Lotka-Volterra com-

mensal symbiosis model with Michaelis-Menten

type harvesting, WSEAS Trans. Math. Vol. 21,

No.1, 2022, pp. 515-523.

[37] Xu L., Xue Y., Lin Q., Lei C., Global attractiv-

ity of symbiotic model of commensalism in four

populations with Michaelis-Menten type harvest-

ing in the first commensal populations, Axioms,

Vol. 11, No. 1, 2022, Article ID 337.

[38] Chen F., Zhou Q., Lin S., Global stability of

symbiotic model of commensalism and para-

sitism with harvesting in commensal populations,

WSEAS Trans. Math. Vol. 21, No.1, 2022, pp.

424-432.

[39] Xu L., Xue Y., Xie X., Lin Q., Dynamic behav-

iors of an obligate commensal symbiosis model

with Crowley-Martin functional responses, Ax-

ioms, Vol. 11, No. 6, 2022. Article ID 298.

[40] Zhang X., Chen L., Neumann A. U., The

stage-structured predator-prey model and optimal

harvesting policy, Mathematical Biosciences,

Vol.168, No. 2, 2000. pp. 201-210.

[41] Wright D. H. , A simple, stable model of mutual-

ism incorporating handling time, The American

Naturalist, Vol. 134, No.4, 1989, pp. 664-667.

[42] Sun G. C., Sun H., Analysis on symbiosis model

of two populations, Journal of Weinan Normal

University, Vol.28, No. 1, 2013. pp. 6-8.

[43] Lemes P., Barbosa F. G., Naimi B., Araujo

M. B. Dispersal abilities favor commensalism in

animal-plant interactions under climate change,

Science of The Total Environment, Vol. 835, No.

1, 2022, Article ID: 155157.

[44] Khajanchi S., Banerjee S., Role of constant prey

refuge on stage structure predator-prey model

with ratio dependent functional response, Applied

Mathematics and Computation, Vol.314, No. 1,

2017, pp.193-198.

[45] Aiello W. G., Freedman H. I., A time-delay

model of single-species growth with stage struc-

ture, Mathematical Biosciences, Vol. 101, No. 2,

1990. pp.139-144.

Contribution of individual authors to

the creation of a scientific article

(ghostwriting policy)

Author Contributions: Please, indicate the role

and the contribution of each author:

Example

John Smith, Donald Smith carried out the simulation

and the optimization.

George Smith has implemented the Algorithm 1.1

and 1.2 in C++.

Maria Ivanova has organized and executed the exper-

iments of Section 4.

George Nikolov was responsible for the Statistics.

Follow: www.wseas.org/multimedia/contributor-

role-instruction.pdf

Sources of funding for research

presented in a scientific article or

scientific article itself

Report potential sources of funding if there is any

Creative Commons Attribution

License 4.0 (Attribution 4.0

International , CC BY 4.0)

This article is published under the terms of the

Creative Commons Attribution License 4.0

https://creativecommons.org/li-

censes/by/4.0/deed.en_US

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2022.21.93

Fengde Chen, Zhong Li, Lijuan Chen

E-ISSN: 2224-2880

824

Volume 21, 2022