Description of Trajectories of an Evolution Operator Generated by Mosquito Population Текст научной статьи по специальности «Математика»

Russian Journal of Nonlinear Dynamics, 2024, vol. 20, no. 2, pp. 197-207. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd240401

NONLINEAR PHYSICS AND MECHANICS

MSC 2010: 39A12, 92D25

Description of Trajectories of an Evolution Operator Generated by Mosquito Population

In this paper, we study discrete-time dynamical systems generated by the evolution operator of mosquito population. An invariant set is found and a Lyapunov function with respect to the operator is constructed in this set. Using the Lyapunov function, the global attraction of a fixed point is proved. Moreover, we give some biological interpretations of our results.

Keywords: Lyapunov function, fixed point, limit point, invariant set

1. Introduction

Mosquito control manages the population of mosquitoes to reduce their damage to human health and economies. Mosquito control is a vital public-health practice throughout the world and especially in the tropics because mosquitoes spread many diseases, such as malaria and the Zika virus.

Mathematical models have been formulated in the literature [4, 5, 7] to study the interactive dynamics and control of the wild and sterile mosquito populations. Many concrete models of mathematical biology described by a corresponding nonlinear evolution operator (see [3, 8]). Since there is no general theory of nonlinear operators, for each such operator one has to use a specific method of investigation.

All mosquitoes undergo complete metamorphosis, going through four distinct stages of development during a lifetime: egg, larva, pupa, and adult. Female mosquitoes lay their eggs on either water or moist soil or interior walls of tree holes, cans, buckets, and old tires that are likely to be flooded by water. All mosquito eggs must have water to develop. Most species prefer slow-moving or stagnant water in which to lay their eggs. Mosquitoes can complete their life cycle within 7 to 10 days during the summer months. Eggs hatch within 48 hours and larvae and pupae live in the water for about a week. Adult mosquitoes will emerge from the pupal case

Received November 10, 2023 Accepted February 22, 2024

Zafar S. Boxonov b.hamidouche@univ-djelfa.dz

V. I. Romanovskiy Institute of Mathematics of Uzbek Academy of Sciences Tashkent International University of Financial Management and Technology University st. 9, Tashkent, 100174 Uzbekistan

Z. S. Boxonov

(a process known as eclosion) and rest on the water's surface until the body dries and hardens before flying away. The life cycle stages can easily be seen in the water [12].

To keep our mathematical modeling as simple as possible, due to the fact that the first three stages in a mosquito's life cycle are aquatic, we group the three aquatic stages into the larvae class by x, and divide the mosquito population into the larvae class and the adults, denoted by y. We assume that the density dependence exists only in the larvae stage [4, 5].

We let the birth rate, that is, the oviposition rate of adults be /, and the rate of emergence from larvae to adults be a function of the larvae with the form of a(1 — k(x)), where a > 0 is the maximum emergence rate, 0 ^ k(x) ^ 1, with k(0) = 0, k'(x) > 0, and lim k(x) = 1,

is the functional response due to the intraspecific competition. We further assume a functional response for k(x), as in [4], in the form

k(x) = —-—. v 7 1 + x

We let the death rate of larvae be a linear function, denoted by d0 + d1 x, and the death rate of adults be constant, denoted by ¡. The dynamics for the wild mosquitoes are governed by the following system [4, 5]:

dx ax

-J7= 7T--(do + dxx)x,

dt 1 + x / -i -i \

dy ax

m-

dt 1 + x

In this paper (as in [9-11]) we study the discrete time dynamical systems associated to the system (1.1). Define the operator W: R2 ^ R2 by

{x' = fly--—--(do + dix)x + x,

, ax +x a.2)

v + y

where a > 0, / > 0, ¡i > 0, d0 ^ 0, d1 ^ 0. In the system (1.2), x', y' means the next state relative to the initial state x, y, respectively.

If for each possible state z = (x, y) G R2 describing the current generation, the state z' = = (x', y') G R2 is uniquely defined as z' = W(z), the map W: R2 ^ R2 is called the evolution operator [1, 8].

The main problem for a given operator W and an arbitrarily initial point z(0) = (x(0), y(0)) G G R2 is to describe the limit points of the trajectory {z(m)}where z(m) = Wm (z(0)).

The dynamics of the operator (1.2) were studied in detail under conditions / = ¡ , d0 — d± —

= 0 in [9] and under conditions f3 / ¡j,, d0 = d1 = 0 in [10]. In [11] under condition f3 < ¡j, + we showed that any trajectory converges to a unique fixed point. Thus, the remaining case is / ^ ^ ^ (l + a1)' which has not been studied yet. In this paper, we study the case ¡3 ^ + and give an analysis of limit points of a trajectory generated by the operator W given by (1.2).

2. Preliminaries and known results

Consider systems of difference equations of the form

ixn+1 = f(xn,yn) n = 0, 1, 2,... (2.1)

[ yn+1 = g(xn, yn),

where f and g are given functions and the initial condition (x0, y0) comes from some considered set in the intersection of the domains of f and g.

Let R be a subset of R2 with nonempty interior, and let W: R ^ R be a continuous operator. When the function f (x, y) is increasing in x and decreasing in y and the function g(x, y) is decreasing in x and increasing in y, the system (2.1) is called competitive. When the function f (x, y) is increasing in x and increasing in y and the function g(x, y) is increasing in x and increasing in y, the system (2.1) is called cooperative. An operator W that corresponds to the system (2.1) is defined as W(x, y) = (f (x, y), g(x, y)). Competitive and cooperative operators, which are called monotone operators, are defined similarly. Strongly cooperative systems of difference equations or operators are those for which the functions f and g are coordinate-wise strictly monotone [6].

Let R+ = {(x, y): x, y e R, x ^ 0, y ^ 0}.

Note that the operator W given by (1.2) is well defined on R2 \ {(x, y): x = —1}. But to define a discrete-time dynamical system of continuous operator as population, we assume x ^ 0 and y ^ 0. Therefore, we choose parameters of the operator W to guarantee that it maps R+ to itself.

Lemma 1 (see [9]). If

a > 0, ((> 0, 0 < p < 1, d0 > 0, a + d0 < 1, dl = 0, (2.2)

then the operator (1.2) maps the set R+ to itself.

Lemma 1 provides sufficient conditions for nonnegative values of x and y.

In this case the system (1.2) becomes

x' = fiy- + do ~ i) x,

Wo : < V1+X J (2.3)

0 ax

y> = — + (l-,)y.

A point z G R+ is called a fixed point of W0 if W0(z) = z.

Definition 1 (see [3]). A fixed point z of the operator W0 is called hyperbolic if its Jaco-bian J at z has no eigenvalues on the unit circle.

Definition 2 (see [3]). A hyperbolic fixed point z is called :

(i) attracting if all the eigenvalues of the Jacobi matrix J(z) are less than 1 in absolute value ;

(ii) repelling if all the eigenvalues of the Jacobi matrix J(z) are greater than 1 in absolute value ;

(iii) a saddle otherwise.

The following proposition is known about the fixed points of the operator W0 and their type.

Proposition 1 ([11]). The types of the fixed points for (2.3) are as follows :

i) if ¡3 ^ ¿t + then the operator (2.3) has a unique fixed point, (0, 0), the point,

attracting if ¡3 < fi ( 1 + —

<°-0) ^ 1 nonkyperbolic if ¡3 = ¡j, i 1 H—-

ii) if ¡3 > At ^ 1 + then the operator has two fixed points (0, 0), (x*, y*), and the point, (x*, y*) is attracting, the point

repelling if ¡3 > a (1 + — ) + a*,

a

(0, 0) = saddle if ^ + ^ < (3 < n (l + ^j+a*,

nonkyperbolic if ¡3 = a (1 H—- ) + a*,

a

where

a* = -(4-2 (a + Li + d0)), a0

= a{(3 - At) _ 1 * = a{(3 - At) - Atrip >2 ^

The following theorem describes the trajectory of any initial point

y(0)) in R+.

Theorem 1 ([11]). For the operator W0 given by (2.3) (i. e., under condition (2.2)) and for any initial point (x(0), y(0)) G R+ the following hold:

(i) If y^ > j- for any natural number n, then lim x^ = +oo, lim y^ =

(ii) If there exists n0 number such that y^n^ ^ j- and (3 < AMI + tt). then lim x^ = 0,

" V a J n—tt

lim y(n) = 0,

y(n)

where (x(n), y(n)) = W^ (x(0), y(0)), with Wq is being the nth iteration of W0.

Denote

tt = \ (x, y) G R+, 0 < x < O^y^ -1.

I + ¡¡d0 ¡¡)

From Theorem 1 one can get that for any initial point z(0) = (x(0), y(0)) taken from R+, in a trajectory, if all y^ are greater than then the limit point is ^+oo, j^j. If is less than for some n0, then y^ is less than ^ for all n > n0. One can see that, if y^ ^ then x(n) < ^ [11]. Therefore, it suffices to study the dynamics of the operator W0 on the set Q

3. Main part

Let us consider the dynamics of the operator W0 given by (2.3) in the sets Q under condition ¡3 ^ At (l + • We divide the set. Q into four parts as follows:

Q1 = {(x, y) G Q, 0 < x < x*, 0 < y < y*},

afî^ a

= { (x, y) G Q, x* < x < y*

= y*

( a

= \ {x, y) e 0 ^ X < x*, y* <y ^ -I p

where x*, y* is the fixed point defined by (2.4).

3.1. Invariant sets

A set A is called invariant with respect to W0 if W0(A) C A. Lemma 2. The sets Q17 Q2 and Q are invariant with respect to W0. Proof. (1) Let (x, y) e Q1, i.e., 0 ^ x ^ x*, 0 ^ y ^ y*. Then

x' — x* = fiy + x (1 — d0--—— ) — x* ^ py* + x* (1 — dQ--—— ) — x*

V 1 + xj \ 1 + x J

ax* ax*

= f3y* - d0x* - < f3y* - d0x* - 1 +'= f3y* - d0x* - fj,y* =

= (ß- ß) • a(/j 7/t} "^ -d0-(^^ - 1 ) = 0, ß(ß - ß) V ßdo

aax^, . j, aax

Thus, (x', y') e Q1, i.e., W0(Q1) C Q1.

(2) Let (x, y) e Q2, i.e., a;* < .r < y* < y < fr Then

x* — x' = x* — fiy — x (1 — d0--—— ) ^ x* — py* — x* (1 — d0--—— ) =

V 1+ xj \ 1 + x J

ax* ax*

= d0x* + —--fly* < d0x* + " - f3y* = d0x* + fiy* - f3y* = 0,

1 + x 1 + x*

So (x', y') G i.e., W0(^2) C

(3) Let (x, |/)efi, fr Then

' aß f a \ aß aß aß ( a \ aß

x--— = ßy + x [ 1 — dQ----— ^--1--— 1 — dQ----—

ßd0 \ 1 + x ) ßd0 ß ßd0 \ 1 + x ) ßd0

aß aß i a \ a2 ß

ß ßd,0 ^ 1+x*) ßd0(l + x) <

a i aa , , aax aa , , a aax a

--y =--(1 - a)y--^--(1 - a)---=- > 0.

ß ß 1 + x ß v > 1+ x 1+ x

Thus, (x',y') G Q, i.e., W0(Q) C Q. □

3.2. Periodic points

A point z in W0 is called a periodic point of W0 if there exists p so that Wp(z) = z. The smallest positive integer p satisfying Wq (z) = z is called the prime period or least period of the point z.

Theorem 2. For p ^ 2 the operator (2.3) does not have any p-periodic point in the set R+.

Proof. Let us first describe periodic points with p = 2 on R+, in this case the equation is W0(W0(z)) = z. That is,

ax . \ 1-he i ax

a (/ty - T^ + U - do>)

a (f3y -+ (I - d,0)x) / ax \

= 1+- ff, + (1 - + (1"(TTi + (1 - ■

(3.1)

By adding the first and second equations of the system (3.1) one can find y:

= x(d,0(2 - d,0)( 1 + x) - a(f3 -fi + d,0)) y (l+x)№-v,-do)-№-»)) ■

If /3 ^ fj, (l + , then /3(2 - fi - d0) > ¿u(2 - fj,). Clearly, if d0(2 - d0) ^ a(f3 - ¡j, + d0),

then x ^ 0 y ^ 0. If d0{ 2 - d0) < a{(3 + d0), then x ^

We obtain the following equation by substituting the y defined in (3.2) into the second equation of the system (3.1):

where

x BX + B2x2 + B3x + B4) = 0, (3.3)

B1 = A0A3,

B2 = 2Aq A3 + Ai A3 + AA - a^0,

B3 = A0A3 + A1A3 + A3 + Aq A4 + Ai A4 - 2aAo - A2, B4 = A3 + A4 - aAo - A2,

/3do(2 - do)

Aq = 1 - do +

3(2 - i - do) - i(2 - i)'

A =1 aP(P-fj, + d0) (3.4)

1 - fj, - dQ) - fj,(2 - fj,)'

An = —a2 (l +

3(2 - i - do) - i(2 - i))'

= n(2 - fj,)do(2 - dQ)

3 /5(2 - n - do) - M2 -

A4 = - a(1 - i) +

0

ai(2 - i)(3 - i + d0)

3(2 — i — d0) — ¡(2 — ¡), Obviously, the fixed points x = 0 and x = x* (see (2.4)) are solutions of Eq. (3.3). Thus, dividing the left-hand side of Eq. (3.3) by x (x — x*), we get the following quadratic equation:

B1x2 + (B2 + B1 x*) x + B3 + B2x* + B1 x*2 = 0. (3.5)

r

Denote

= a(3 - fj, + d0) _ ° d0(2-d0)

If B0 ^ 0, then

(2 - f)(2 - d0) - a(2 + 3 - f) ^ 0. (3.6)

Indeed, (2 - f)(2 - d0) - a(2 + / - f) = (2 - f)(2 - d0) - 2a - a(3 - f) ^ (2 - f)(2 - d0) -- 2a - d0(2 - d0) + ad0 = (2 - d0)(2 - f - d0 - a) ^ 0. We write the coefficients of Eq. (3.5) as follows:

B1 — Ao A3,

B2 + B1x* —

ßdo(ß - ß)(2 - ß)(2 - do)((2 - ß)(2 - do) - a(2 + ß - ß))

2 ^ 1 f/^O _ M _ rl \ _ _ „^2

B3 + B2 x* + B1x*2 —

(ß(2 - ß - do) - ß(2 - ß))2 ' (3.7)

2 _ ßdo((2 - ß)(2 - do) - a(2 + ß - ß))

ß(2 - ß - do) - ß(2 - ß)

It can be seen from (2.2) and (3.6) that, if B0 ^ 0, then the coefficients of Eq. (3.5) are positive. Therefore, there is no positive solution of Eq. (3.5).

Let us consider the case B0 > 0. For y to be nonnegative, it must be x ^ B0. So, we solve Eq. (3.5) in [B0, +rc>). In Eq. (3.5) we put x + B0 instead of x and get:

B1x2 + (B2 + Blx* + 2B0B1) x + B3 + B2 (x* + B0) + Bl (x*2 + B02 + B0x*) = 0. (3.8)

We need to find positive solutions (3.8) instead of looking for all solutions (3.5) in [B0, +rc>). If B0 > 0, then

(3 - x)(1 - d0) >d0. (3.9)

Indeed, (3 - f)(1 - d0) - d0 ^ a(3 - f) - d0 > d0(2 - d0 - a) - d0 = d0(1 - d0 - a) ^ 0. If (3 - f)(1 - d0) > d0, then

(3 - f)2(1 - d0) > (3 - f + 1)d2. (3.10)

Indeed, (3 - f)2(1 - d0) - (3 - f + 1)d§ > (3 - f)d0 - (3 - f)d2 - d2 = (3 - f)d0(1 - d0) -- d0 > d0 - d0 = 0.

We write the coefficients of Eq. (3.8) in the following form:

B2 + B1 x* + 2Bo B1 — BoB1 +

Bi — Ao A3,

ß(2 - ß)((do(2 - do) + a((ß - ß)(1 - do) - do))

01 011 ß(2-ß-d0)-ß(2-ß) B3 + S2 (,* + S0) + (x*2 + 2* + B0X*) = (311)

+

do(2 - do)(ß(2 - ß - do) - ß(2 - ß))

aßa{2 - ß)((ß - ß)2(l - d0) - (3 - ß + d0(2 - d0)(ß(2 - ß - d0) - ß(2 - ß))

From (2.2), (3.9) and (3.10) one can see that all coefficients of Eq. (3.8) are positive, i.e., there is no positive solution. Thus, there is no positive solution for Eq. (3.8).

Finally, the operator W0 does not have two periodic points in the set R+.

The Jacobian matrix of the operator W0 has a sign configuration at all points z

f+ +\

sign v+ +

Therefore, the system (2.3) is strictly cooperative. Since the system (2.3) is strictly cooperative, Sharkovskii's ordering holds for periodic points [2] and so the nonexistence of a two-period point would imply the absence of periodic points of all other periods. □

3.3. Lyapunov functions

By using Lyapunov functions, one can handle an w-limit point. Recall the definition of a Lyapunov function.

Definition 3. A continuous functional y: Q ^ R is called a Lyapunov function for a operator W if y(W(z)) ^ y(z) for all z (or y(W(z)) ^ y(z) for all z).

Proposition 2. For any ¡3, / (under condition (2.2)), the function

y(z) = /x + ¡y is a Lyapunov function for the operator W0 defined by (2.3).

Remark 1. If /3 = ¿t (l + ^j, then z = (x, y) e Q, and if /3 > fj, (l + ^j, then z = (x, y) e e u Q2.

Proof. Let 3 = ^ (l + • Then

^(Wo(z)) = pa! + M = px + = <p(z) ~

1 | X 1 | x

Thus, y(W0(z)) ^ y(z) for all z = (x, y) e Q, that is, the function y(z) is a Lyapunov function for the operator W0.

Let f3 > ¡j, + . Let us look at the expression

(XX

ip(W0(z)) - ip(z) = fix' + /V - (fix + py) = (/3- - d^LX =

d>ol-ix (a(/3 - fj) _ 1 _ \ = d,0Lix _

1 + x \ d0 / J 1 + x

If z = (x, y) e Q^ then y(W0(z)) ^ y(z), if z = (x, y) e Q2, then y (W0(z)) ^ y(z).

Hence, the function y(z) is increasing on the set Q1 and decreasing on the set Q2. The function y(z) is a Lyapunov function for the operator W0. □

3.4. The w-limit set

The problem of describing the w-limit set of a trajectory is of great importance in the theory of dynamical systems.

The following theorem describes the trajectory of any initial point (x(0), y(0)) in Q.

Theorem 3. For the operator W0 given by (2.3) (i. e, under condition (2.2)) the following

hold:

(i) If ¡3 = fi + , then for any initial point, ye Q \ Fix W0

lim x(n) = 0, lim y(n) = 0;

n—tt n—tt

(ii) If /3 > ¡j, + , then for any initial point, ye Q \ Fix W0

lim x(n) = x*, lim y(n) = y*;

n—tt n—tt

where (x(n), y(n)) = W01 (x(0), y(0)), with Wq being the nth iteration of W0 and (x*, y*) is the fixed point defined by (2.4).

Proof. We have

ax(n-1) ax(n-1)

x{n) = py(n~l) - {n „ y(n) = T^TIT -w(n-l)+y(n~l)- (3.12)

1+ x(n-1) 1+ x(n-1)

First, we prove assertion (i). Let f3 = /x^l + ^j. The function ip(z) is a decreasing Lyapunov function by Proposition 2 and bounded from below. This implies the existence of lim ^>(z(n)):

n—tt

2

a(

hm fz{n)) = hm L - a[ß~ß) = hm -

n—œ \ ) n—œ I V ) 1+ X(n-l> I n—œ \ )

(x(n-1))2 (x(n-1))2 . ,

-a(P-fi) lim --^rr lim —--——r = 0 lim = 0.

n—tt 1 + x(n-1) n—tt 1 + x(n-1) n—tt

By (3.12) and lim x(n) =0 we have lim y(n) = 0.

n—tt n—tt

Let us prove assertion (ii). We consider the trajectory of any initial point (x(0), y(0)) in the sets Q1 and Q2. The function ^>(z) is increasing on the set Q1 and is bounded from above. The function ^>(z) is decreasing in the set Q2 and is bounded from below. This implies the existence of

lim p (z(n)) = 0, (3.13)

n—œ V /

lim <p(zM) = lim f <p fz(n-V) + d°Mx('2 fx*-x^)] = lim <p(z(n~») +

n—œ \ J n—œ\ V J 1+ x(n—^ V Jj n—œ

x(n—1) , , x(n— 1) , ,

+doß lim , fx* - x{n~l) ) lim , fx* - X(ra"1} ) = 0

n—œ 1+ x(n—^ V J n—œ 1+ X(n—1) V J

=> lim x(n—1) =0 or lim x(n—1) = x*.

Assume lim x(n) = 0. Then, from (3.12) we have lim y(n) = 0. Hence, lim p (z(n)) = 0. This

n—œ n—œ n—œ

is a contradiction to (3.13). So, lim x(n) = x*.

n—œ

By (3.12) and lim x(n) = x* we have lim y(n) = y*

n^tt

Now, we consider the trajectory of an initial point (x(0), y(0)) in the sets Q3 and Q4. Suppose (x(n), y(n)) e for all n e N and for any initial point (x(0), y(0)) taken from the set Q3. Then the sequence x(n) is decreasing and bounded from below, the sequence y(n) is increasing and bounded from above. Indeed, from x* < x^ ^ 0 ^ y^ ^ y* we get:

^(n+1) - = _ _ doX(n) ^ f3y* _ _ doX* = 0>

v(n+1) _ (n) = ax{n)__In) > ax*__* = o

y y 1 + x(n) 1 + x*

Hence, the sequences x(n), y(n) have their limits. From (3.12) it follows that lim x(n) = x*,

n^tt

lim y(n) = y*.

n^tt

Assume (x(n), y(n)) e 04 for all n e N and for any initial point (x(0), y(0)) taken from the set Q4. Then the sequence x(n) is increasing and bounded from above, the sequence y(n) is decreasing and bounded from below. Indeed, from 0 ^ x^ < x*, y* < y^ ^ jL we get:

sCH-U _ = f5y[n) _ _ dox{n) ^ f3y* _ _ doX* = 0)

1 + x(n) 0 1 + x

y(n+l) - y(n) = T^T - W(n) < T^-7 - W* = 01 + x(n) 1 + x*

So, the sequences x(n), y(n) have their limits. Since (3.12), we have lim x(n) = x*, lim y(n) =

n^tt n^tt

= y*. □

4. Biological interpretations

In biology a population biologist is interested in the long-term behavior of the population of a certain species or collection of species. Namely, what happens to an initial population of members. Does the population become arbitrarily large as time goes on? Does the population tend to zero, leading to extinction of the species [8]? In this section we briefly give some answers to these questions related to our model of the mosquito population.

Each point (vector) z = (x, y) e R+ can be considered as a state (a measure) of the mosquito population.

Let us give some interpretations of our main results (interpretation of Theorem 3):

(a) Let ¡3 = ^{l + . Under this condition on ¡3 (i. e., the birth rate of adults), the mosquito population dies;

(b) If the inequality ¡3 > ^{l + holds for /3, then the mosquito population tends to the equilibrium state (x*, y*) with the passage of time;

(c) Let [3 < + . Under this condition on /3, the mosquito population dies (interpretation of Theorem 1).

Acknowledgment

The author is grateful to Professor U.A.Rozikov for his comments and suggestions.

Conflict of interest

The author declares that he has no conflict of interest.

References

[1] Lyubich, Yu.I., Mathematical Structures in Population Genetics, Biomathematics, vol. 22, Berlin: Springer, 1992.

[2] Wang, Y. and Jiang, J., The General Properties of Discrete-Time Competitive Dynamical Systems, J. Differential Equations, 2001, vol. 176, no. 2, pp. 470-493.

[3] Devaney, R.L., An Introduction to Chaotic Dynamical Systems, 3rd ed., Boca Raton, Fla.: CRC, 2022.

[4] Li, J., Malaria Model with Stage-Structured Mosquitoes, Math. Biosci. Eng., 2011, vol. 8, no. 3, pp. 753-768.

[5] Li, J., Cai, L., and Li, Y., Stage-Structured Wild and Sterile Mosquito Population Models and Their Dynamics, J. Biol. Dyn, 2017, vol. 11, suppl. 1, pp. 79-101.

[6] Bilgin, A., Brett, A., Kulenovic, M.R. S., and Pilav, E., Global Dynamics of a Cooperative Discrete System in the Plane, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 2018, vol. 28, no. 7, 1830022, 17 pp.

[7] Rozikov, U.A. and Velasco, M.V., A Discrete-Time Dynamical System and an Evolution Algebra of Mosquito Population, J. Math. Biol., 2019, vol. 78, no. 4, pp. 1225-1244.

[8] Rozikov, U. A., Population Dynamics-Algebraic and Probabilistic Approach, Hackensack, N.J.: World Sci., 2020.

[9] Boxonov, Z. S. and Rozikov, U.A., A Discrete-Time Dynamical System of Stage-Structured Wild and Sterile Mosquito Population, Nonlinear Stud., 2021, vol. 28, no. 2, pp. 413-425.

[10] Boxonov, Z. S. and Rozikov, U. A., Dynamical System of a Mosquito Population with Distinct Birth-Death Rates, J. Appl. Nonlinear Dyn., 2021, vol. 10, no. 4, pp. 791-800.

[11] Boxonov, Z. S., A Discrete-Time Dynamical System of Mosquito Population, J. Differ. Equ. Appl., 2023, vol. 29, no. 1, pp. 67-83.

[12] Mosquito Life Cycle, https://www.baycounty-mi.gov/MosquitoControl/MosquitoLifeCycle.aspx