# CryptoLight: An Electro-Optical Accelerator for Fully Homomorphic Encryption

Mengxin Zheng  Indiana University of BloomingtonBloomingtonUSA<EMAIL_ADDRESS>,  Qian Lou  University of Central FloridaOrlandoUSA<EMAIL_ADDRESS>,  Fan
Chen  Indiana University of BloomingtonBloomingtonUSA<EMAIL_ADDRESS>,  Lei Jiang
Indiana University of BloomingtonBloomingtonUSA<EMAIL_ADDRESS>and  Yongxin
Zhu  Shanghai Advanced Research Institute, Chinese Academy of
SciencesShanghaiChina<EMAIL_ADDRESS>
(2022)

###### Abstract.

Fully homomorphic encryption (FHE) protects data privacy in cloud computing by
enabling computations to directly occur on ciphertexts. Although the speed of
computationally expensive FHE operations can be significantly boosted by prior
ASIC-based FHE accelerators, the performance of key-switching, the dominate
primitive in various FHE operations, is seriously limited by their small bit-
width datapaths and frequent matrix transpositions. In this paper, we present
an electro-optical (EO) FHE accelerator, CryptoLight, to accelerate FHE
operations. Its 512-bit datapath supporting 510-bit residues greatly reduces
the key-switching cost. We also create an in-scratchpad-memory transpose unit
to fast transpose matrices. Compared to prior FHE accelerators, on average,
CryptoLight reduces the latency of various FHE applications by ¿$94.4\%$ and
the energy consumption by ¿$95\%$.

††journalyear: 2022††copyright: acmcopyright††conference: 17th IEEE/ACM
International Symposium on Nanoscale Architectures; December 7–9, 2022;
Virtual, OR, USA††booktitle: 17th IEEE/ACM International Symposium on
Nanoscale Architectures (NANOARCH ’22), December 7–9, 2022, Virtual, OR,
USA††price: 15.00††doi: 10.1145/3565478.3572544††isbn:
978-1-4503-9938-8/22/12††ccs: Security and privacy Privacy protections††ccs:
Hardware Emerging optical and photonic technologies

##  1\. Introduction

Data privacy is gaining tremendous importance around the world. This leads a
surge in demand for privacy-preserving computing solutions protecting data
confidentiality while in transit, rest, and in-use. Therefore, data privacy
has become a major barrier to cloud computing, where a client relies on a
server in the cloud. Moving sensitive data, e.g., healthcare reports and
financial records, to untrusted cloud servers is vulnerable to data breaches.

Fully homomorphic encryption (FHE) (Gentry:ACC2012, ) emerges as one of the
most promising solutions to guaranteeing data privacy by allowing computations
to directly happen on ciphertexts. Though trusted execution environments
protect both code and data, they are susceptible to side-channel attacks
(Oleksenko:ATX2018, ). Secure multi-party computation (SMPC)
(Pratyush:USENIX2020, ) protects data privacy in an interactive way. However,
SMPC requires the continuous online presence of all involved parties, each of
which needs a high-bandwidth network connection. For instance, a SMPC-based
private network (Pratyush:USENIX2020, ) has to exchange $2GB$ data for only an
inference on an encrypted CIFAR-10 image. In contrast, FHE enables a client to
encrypt her data and to send only the ciphertexts to an untrusted server,
which can perform computations on the ciphertexts even when the client is
offline. Whenever the client comes back online, she can retrieve and decrypt
the encrypted results. The client and the server do not have to communicate
during FHE computations. Since encrypted intermediate results in the server
are not shared, less is leaked to the client.

FHE operations are extremely time-consuming, i.e., one FHE bootstrapping costs
several seconds on a CPU. Recent ASIC-based hardware accelerators
(Samardzic:ISCA2022, ; Kim:ISCA2022, ) greatly speed up various FHE
operations. However, the performance of prior FHE accelerators is limited by
their small bit-width datapaths and frequent matrix transpositions. An FHE
ciphertext consists of two polynomials of large degrees (e.g., several
thousand) with large integer coefficients (e.g., several hundred bits). To
efficiently compute with polynomials, FHE schemes (e.g., CKKS
(Samardzic:ISCA2022, )) adopt Residue Number System (RNS) (Samardzic:ISCA2022,
), and Number Theoretic Transform (NTT). First, to compute with large integer
coefficients, RNS divides each coefficient into multiple smaller bit-width
(e.g., 60-bit) residues, each of which can be processed by the datapath of
prior FHE accelerators. The latency of an FHE multiplication, rotation, or
bootstrapping is dominated by expensive key-switching (KS) primitives
(Gentry:ACC2012, ) that makes output ciphertexts to be encrypted by the same
secret key as the input ciphertext(s). The computational overhead of KS, i.e.,
the number of NTTs, additions, and multiplications in a KS, greatly increases
with an enlarging number of residues. Second, for two large degree-$N$
polynomials, NTT and inverse NTT (INTT) reduces the time complexity of their
multiplications to $\mathcal{O}(N\log N)$. But it is difficult to perform an
(i)NTT, i.e., NTT or iNTT, on a large degree-$N$ polynomial directly. Prior
FHE accelerators (Samardzic:ISCA2022, ; Kim:ISCA2022, ) place it as an
$n\times n$ matrix, where $N=n^{2}$, perform an (i)NTT on each row, multiply
the matrix with some constants, transpose the matrix, and perform an (i)NTT on
each row again. As a result, frequent matrix transpositions greatly prolong
the latency of KS in various FHE operations by introducing huge volumes of on-
chip memory traffic during (i)NTTs.

We propose an electro-optical (EO) FHE accelerator, CryptoLight, to support a
large bit-width datapath and free its computing units from matrix
transpositions. Our contribution is summarized as:

  * •

A 512-bit EO CU. We propose a 512-bit EO Computing Units (CUs) built upon
ultra-fast EO integer adders and multipliers to process polynomials with large
coefficients. Its 512-bit datapath greatly reduces the number of residues, and
thus the computational overhead of each key-switching primitive.

  * •

An in-SPM TU. We build a low-power eDRAM-based on-chip scratchpad (SPM)
system, where we present an in-SPM Transpose Unit (TU) to transpose matrices
for (i)NTTs without sending matrices having low spatial locality to the CUs
via H-trees of the SPM or NoCs.

  * •

Latency and energy efficiency. Compared to prior ASIC-based FHE accelerators,
CryptoLight reduces the latency of various FHE applications by ¿$94.4\%$ and
the energy consumption by ¿$95\%$.

##  2\. Background

###  2.1. Fully Homomorphic Encryption

Fully Homomorphic Encryption. FHE (Gentry:ACC2012, ) allows computations to
directly occur on ciphertexts. To guarantee security, FHE introduces noise
into ciphertexts. An FHE operation enlarges the noise in the ciphertext which
can accommodate only a limited number of FHE operations (a.k.a.,
multiplicative depth) without a decryption failure. Bootstrappings keep the
noise in check without a security key. But a bootstrapping is slow, e.g.,
several minutes.

CKKS. CKKS (Cheon:ICTACIS2017, ) is one of the widely used FHE schemes in
encrypted machine learning (Kim:JMI2018, ), since it supports complex
arithmetics. CKKS encrypts a vector of complex numbers as a polynomial, whose
coefficients are integers modulo $Q$ and polynomial degree is $N$. For a given
$N$, a ciphertext can encrypt up to $N/2$ complex numbers for SIMD operations,
each of which is an element-wise multiplication or addition between two
ciphertexts. Its critical operations include:

  * •

FMUL: An FHE multiplication (FMUL) is implemented using polynomial
multiplications and additions.

  * •

FROT: An FHE rotation (FROT) rotates the vector encrypted in a ciphertext. ROT
performs an automorphism on the ciphertext.

  * •

FBOT: An FHE bootstrapping (FBOT) reduces the noise in a ciphertext, and
refreshes its consumed multiplicative depth.

Key-Switching. Key-Switching (KS) (Gentry:ACC2012, ) makes the resulting
ciphertext of an FMUL, FROT, and FBOT stay encrypted by the same secret key as
the input ones. KS is the most expensive and dominate kernel in these
operations, i.e., it takes ¿90% (Samardzic:ISCA2022, ; Kim:ISCA2022, ) of all
operations. KS expands the input polynomial to use wider coefficients,
multiplies the polynomial by KS hint matrices (Samardzic:ISCA2022, ), and then
converts the resulting polynomial back to use original coefficients. By using
different KS hint matrices, KS can adjust the trade-off between the noise cost
and the time complexity of an FHE operation.

Number Theoretic Transform. Multiplying 2 polynomials of degree-$N$ has the
time complexity of $\mathcal{O}(N^{2})$. NTT (Riazi:ASPLOS2020, ), a variant
of FTT for modular arithmetic, reduces the polynomial multiplication time
complexity to $\mathcal{O}(N\log N)$. Specifically,
$NTT(\mathbf{a}\cdot\mathbf{b})=NTT(\mathbf{a})\bigodot NTT(\mathbf{b})$,
where $\mathbf{a}$ and $\mathbf{b}$ are two polynomials of degree-$N$, and
$\bigodot$ is element-wise multiplication. Typically, a large degree-$N$
($N=n^{2}$) polynomial is placed as an $n\times n$ matrix. The NTT operation
on the large polynomial can be performed in four steps (Samardzic:ISCA2022, ).
First, $n$ $n$-element NTT operations, each of which occurs on one row of the
matrix, can be done. Second, the matrix is multiplied with some constants.
Third, the matrix is transposed. Finally $n$ $n$-element NTT operations, each
of which happens on one row of the matrix, are conducted. Matrix
transpositions greatly prolong the latency of an NTT by introducing huge
volumes of memory traffic (Samardzic:ISCA2022, ).

Residue Number System. RNS (Riazi:ASPLOS2020, ) denotes an integer by a set of
residues modulo predefined pairwise co-prime moduli. Through RNS, a CKKS
operation working with a large $Q$ is converted to a set of computations on
$Q/W$ residues, where $W$ is the datapath bit-width of the FHE accelerator.
For instance, a 560-bit Q is decomposed to 10 60-bit residues on a CPU, 11
54-bit residues on a FPGA (Riazi:ASPLOS2020, ), or 20 28-bit residues on a GPU
(Palisade:GPU2020, ) and ASIC (Samardzic:ISCA2022, ). However, for a given
$Q$, smaller bit-width residues significantly increase the computational
overhead of KS, i.e., the number of NTTs, multiplications, and additions in
each KS, in various FHE operations.

(a) A 1-bit EO FA (PD: photodetector).

(b) The adder frequency.

Figure 1. A photonic ripple-carry adder.

###  2.2. Electro-optical Ripple-Carry Adders

Optical microdisk (Ying:JSTQE2018, ) emerges as one of the most promising
technologies for photonic computing, due to its CMOS compatibility, low power
consumption, and ultra-fast speed. A 1-bit microdisk-based electro-optical
(EO) full adders (FAs) (Ying:JSTQE2018, ) is shown in Figure 1(a). Via $n$
FAs, an $n$-bit EO ripple-carry adder can be built. For the adder, its sum is
calculated as $S_{n}=C_{n-1}\oplus(A_{n}\oplus B_{n})=C_{n-1}\oplus P_{n}$,
while its carry is computed as $C_{n}=(A_{n}\oplus B_{n})\cdot
C_{n-1}+A_{n}\cdot B_{n}=P_{n}\cdot C_{n-1}+G_{n}$, where $P_{n}$ is the
propagate bit, and $G_{n}$ is the generate bit. Since the carry computation is
on the critical path of a ripple-carry adder, CMOS gates are used to compute
$P_{n}$, $G_{n}$, and $S_{n}$, but microdisks are adopted to compute the
carry. The optical carry signal $C_{n-1}$ is separated into two halves by a
splitter. One half of $C_{n-1}$ is sent to a photodetector and converted to an
electrical signal, which goes through an XOR gate together with $P_{n}$ to
compute the sum $S_{n}$. The other half of $C_{n-1}$ goes through a microdisk
$b$ controlled by $P_{n}$ to produce $P_{n}\cdot C_{n-1}$. A combiner merges
$P_{n}\cdot C_{n-1}$ and $G_{n}$ generated by another microdisk $a$. The
continuous wave $\lambda$ is tuned to ensure the balance of light intensities
at the two arms of combiners. When the bit-width of a ripple-carry adder
increases, as Figure 1(b) shows, the adder frequency significantly decreases
(Dinechin:ICFPLA2010, ). Although a carry look-ahead adder can maintain high
frequency for large bit-widths, the carry look-ahead logic components consume
nontrivial power. Compared to the electrical counterpart, the $n$-bit EO
ripple-carry adder greatly enhances its frequency without introducing large
power overhead. We build a wide datapath operating at high frequency by
$n$-bit EO ripple-carry adders in this paper.

(a) operation # v.s. datapath width.

(b) matrix transpos, BC: base conversion.

Figure 2. The bottlenecks of prior FHE accelerators.

###  2.3. Related Work and Motivation

FHE operations computing with large polynomials are extremely time-consuming.
Prior work proposes GPU (Palisade:GPU2020, )-, FPGA (Riazi:ASPLOS2020, )-, and
ASIC (Samardzic:ISCA2022, ; Kim:ISCA2022, )-based accelerators to process FHE
operations. Among all, the ASIC-based FHE accelerators, CraterLake
(Samardzic:ISCA2022, ) and BTS (Kim:ISCA2022, ), obtain the state-of-the-art
performance. However, their performance is seriously limited by their narrow
datapaths and intensive matrix transpositions. First, the narrow datapaths of
prior FHE accelerators greatly increase the number of additions (adds),
multiplications (mults), and NTTs in a key-switching (KS), which is the
dominate primitive in FMULs, FROTs, and FBOTs. CraterLake (Samardzic:ISCA2022,
) has a 32-bit modular integer datapath, while BTS (Kim:ISCA2022, ) supports
64-bit modular arithmetics. To compute a KS with a 1024-bit ciphertext modulus
$Q$, 32-bit CraterLake performs 3.2K adds, 3.1K mults, and 192 NTTs, while
64-bit BTS computes 832 adds, 800 mults, and 96 NTTs, as shown in Figure 2(a).
A narrower datapath exponentially increases the number of adds and mults, and
linearly enlarges the number of NTTs in a KS. In contrast, a wider datapath
may greatly decrease the number of adds, mults, and NTTs in each KS. However,
a wider datapath (e.g., 512-bit) also exponentially increases the power
consumption and the chip area of an FHE accelerator. Second, frequent matrix
transpositions greatly prolong the latency of an NTT on prior FHE accelerators
by introducing huge volumes of memory traffic. To transpose a large (e.g.,
$256\times 256$) matrix, all processing elements of prior FHE accelerators
(Samardzic:ISCA2022, ; Kim:ISCA2022, ) have to frequently access scratchpad
memory arrays, and thus cannot focus on the NTT computations. As a result, the
matrix transpositions (NTT-T) during NTTs in various FHE operations become the
largest bottleneck, as highlighted in Figure 2(b), while the computations of
NTTs (NTT) consume only ¡$10\%$ of the FHE operation latency.

Figure 3. The architecture of CryptoLight.

##  3\. CryptoLight

###  3.1. The Architecture of CryptoLight

We propose an EO FHE accelerator, CryptoLight, to process FHE operations. The
architecture of CryptoLight is shown in Figure 3. CryptoLight has 2048 512-bit
computing units (CUs). All CUs are connected to a 512-bank 512MB eDRAM
scratchpad memory (SPM) system by optical crossbar NoCs. In each SPM bank,
there are multiple transpose units (TUs) to transpose matrices for (i)NTTs
without sending the data to CUs via NoCs. CryptoLight has two memory
controllers (MCs) and two HBM2 PHYs to communicate with the off-chip memory.

###  3.2. A 512-bit EO CU

We build a 512-bit EO CU featured by an NTT unit, a modular add/mult unit, an
automorphism unit, and a TRNG unit. Its most important component is the
512-bit EO NTT unit, which has an arithmetic and inversion unit, an address
generation unit, and two butterfly units. A 512-bit EO NTT unit also supports
the kernel of iNTT working in a different data flow. We also use EO adders and
multipliers to construct the other CU components.

An EO NTT Unit. We present an EO NTT unit for the CU. Matrix transpositions
are done by TUs in SPM banks, but the other three steps of the NTT on a large
polynomial are computed by an NTT unit. CryptoLight aims to support
64K-element NTTs, so a NTT unit supports 256-element NTT operations. The
details are summarized as follows.

  * •

Data flow: Like prior FHE accelerators (Samardzic:ISCA2022, ; Kim:ISCA2022, ),
we adopt the Cooley-Tukey data flow (Riazi:ASPLOS2020, ) for (i)NTTs.

  * •

A butterfly unit: We propose an EO butterfly unit (BU) to accelerate radix-2
NTT butterflies, as shown in Figure 4(a). A BU consists of an EO pipelined
integer array multiplier, an EO Montgomery modular reduction unit, and two EO
modular adders. The EO multiplier computes the multiplication between the
input and a twiddle factor $\omega$. The EO Montgomery modular reduction unit
performs modular reduction on the multiplication result. By an EO adder and a
comparator, the EO modular adder performs modular additions and subtractions.
Two EO modular adders can generate the radix-2 butterfly outputs concurrently.

  * •

A pipelined EO array multiplier: Through EO ripple-carry adders
(Ying:JSTQE2018, ), we propose an EO pipelined integer array multiplier. We
show the example of a 4-bit $\times$ 4-bit pipelined array multiplier in
Figure 4(b). An $m$-bit $\times$ $m$-bit pipelined array multiplier consists
of $m$ stages, each of which is an $m$-bit EO ripple-carry adder. Between two
stages, there is an $m$-bit register file to buffer the intermediate result.
The inputs of the first pipeline stage of the multiplier are the results of
AND operations between the corresponding bits of the multiplier inputs. The
input of the other multiplier pipeline stage is the shifted output of the
previous multiplier pipeline stage. Although an $m$-bit $\times$ $m$-bit
multiplication takes $m$ cycles, the pipelined array multiplier generates an
output per cycle. The wide datapaths and in-SPM TUs of CryptoLight make the
(i)NTT kernels memory-friendly, so the (i)NTT kernels become more
computationally intensive. Since each stage, i.e., a $512$-bit EO ripple-carry
adder, can operate at 3GHz (Ying:JSTQE2018, ), the pipeline of a 512-bit
$\times$ 512-bit EO array multiplier also works at 3GHz. Since only the EO
devices on its critical path work at such high frequency, based on our
estimation, the power consumption of the 512-bit $\times$ 512-bit EO array
multiplier is similar to that of its CMOS counterpart operating at
$\sim$800MHz.

  * •

A Montgomery modular reduction unit: As Figure 4(c) shows, we build an EO
Montgomery modular reduction unit (MMRU) in a BU to perform modular reduction
operations. The MMRU implements the modular reduction algorithm
(Riazi:ASPLOS2020, ) shown in Figure 4(d). Besides some logic operations, and
2’s complement conversions, the most intensive operation in a modular
reduction operation is the multiply-add operation (i.e., $T1_{H}+(q_{H}\cdot
T2)+z_{in}$), which can be computed by EO adders and an EO multiplier. The
output of each iteration of the loop can be cached in a register file and used
as the input for the next iteration. A MUX selects one between the outputs of
the register file and an EO adder as the modular reduction output.

Figure 4. An NTT unit in an EO CU, and a TU in a SPM.

A modular add & mult unit. In a CU, we group two modular adders and two
modular multipliers to construct a modular add/mult unit. A modular adder
consists of an EO 512-bit integer adder and a comparator, while a modular
multiplier is composed of an EO 512-bit $\times$ 512-bit integer pipelined
array multiplier and a Montgomery modular reduction unit.

An automorphism unit. In an automorphism unit, we use 128 CMOS 32-bit $\times$
32-bit integer multipliers, 128 binary shifters, and several binary logic
gates to compute the index permutations. But matrix transpositions are done by
TUs.

A TRNG unit. Since key generation and encryption rely on true random number
generation (TRNG), in a CU, we adopt a CMOS TRNG generator (Mathew:JSSC2016,
), which combines the the entropy of multiple independent sources to generate
a TRNG bitstream. It produces $162.5M$ random bits per second and consumes
only $\sim 1.09mW$.

###  3.3. An eDRAM-based SPM System with TUs

Scratchpad Memory. To avoid the large power of SRAM, we present a 512MB eDRAM-
based scratchpad memory (SPM) system consisting of 512 banks, which is large
enough to store ciphertexts and KS hint matrices shared by FHE operations. To
reduce the refresh power of the eDRAM SPM, we skip refreshes on the rows which
are accessed by a TU or a CU in each refresh period. All CUs access SPM banks
via two optical crossbar NoCs.

Transposing a Matrix. Matrix transpositions are heavily used in the (i)NTT and
automorphism kernels of FHE operations. As Figure 4(e) shows, we implement the
recursive algorithm (Samardzic:ISCA2022, ) to transpose a large matrix. An
$E\times E$ matrix is stored in the eDRAM in a row-major order, i.e., the
element$[n,m]$ ($0\leq n,m\leq E-1$) is stored at the address of $E\times
n+m$. After the transposition, the address of the element is $E\times m+n$.
For the transposition of a large $E\times E$ matrix, we divide the matrix into
four $\frac{E}{2}\times\frac{E}{2}$ matrices, i.e., $A$, $B$, $C$, and $D$, at
the top level. Instead of transposing the matrix directly, we compute $A^{T}$,
$C^{T}$, $B^{T}$, and $D^{T}$. By repeating this process recursively, the
$E\times E$ matrix can be transposed.

In-SPM Transpose Unit. We create an in-SPM transpose unit (TU) to transpose a
matrix inside SPM banks without sending it to CUs via NoCs. As Figure 4(f)
shows, along the H-tree of all SPM banks, we hierarchically deploy TUs in all
SPM banks. We assume a SPM bank has four sub-arrays. Each sub-array in a SPM
bank has a TU. One bank has a level-2 TU to enable data movements inside the
bank. All banks share a level-3 TU for inter-bank communication. The level-3
and level-2 TUs perform the recursive matrix transposition algorithm. When the
recursive process reaches the $2\times 2$ matrix level, a TU attached to a
sub-array swaps each element into its new position. The index of an element in
the $2\times 2$ matrix is $0\leq i\leq 4$. The new index ($j$) after the
transposition can be computed by as: if $i=3$, $j=3$; otherwise, $j=(2i\mod
3)$. A TUs share a similar structure shown in Figure 4(g) with different
numbers of FIFOs.

###  3.4. Compiling FHE Operations on CryptoLight

We build a compilation framework to compile C/C++ source code to executable
binaries for CryptoLight. Our compilation framework is shown in Figure 5.

  * •

First, a C/C++-based application can be translated to FHE primitives, e.g.,
FMUL, FROT and FBOT, by an FHE operation compiler Cingulata (Carpov:IWSCC2015,
). The compiler reorders FHE operations to reuse each KS hint as much as
possible to minimize the DRAM traffic.

  * •

Second, we modified the state-of-the-art FHE accelerator simulator, Sapphire-
Sim (Banerjee:TCHES2019, ), by adding an Operation-to-Function-Unit (OtFU)
mapper. The OtFU mapper decomposes each FHE operation to function kernels,
e.g., (i)NTT and SIMD add/mult, converts these function kernels to a data
dependency graph, and then maps the graph to function units of CryptoLight
(specified by the architecture description). Each function unit can be an NTT
unit, a modular add/mult unit, an auto unit, a TRNG unit, or a TU. To simplify
the mapping, the mapper considers the SPM and DRAM as function units, so the
graph includes accesses to on-/off-chip memories. The mapper statically
schedules every node of the graph to a function unit to maximize the operation
concurrency and the utilization of each function unit using the worst-case
latencies of function units and DRAM.

  * •

Third, we modified the cycle-level scheduler of Sapphire-Sim to further refine
the data dependency graph. The scheduler divides a function kernel into
multiple sub-kernels, schedules sub-kernels across all components of a
function unit, and allocates data in the register files of the function unit
to maximize the data locality. This step determines the exact cycles of all
sub-kernels, and produces the binaries all function units. The scheduler also
minimizes the off-chip data movement in multiple passes.

Figure 5. Compilation on CryptoLight. Table 1. The power and area of CryptoLight ($7$nm). Name | Component | Spec | Power | Area
---|---|---|---|---
| (i)NTT |  $\times 2$ BUs | 7.98$mW$ | 28,503$\mu m^{2}$
| Add/Mult |  $\times 2$ mod adds/mults | 5.94$mW$ | 21,852$\mu m^{2}$
CU | auto |  $\times 128$ mults/shifts | 2.92$mW$ | 14,459$\mu m^{2}$
| TRNG | $\times 1$ | 1.09$mW$ | 1,501$\mu m^{2}$
| reg. files | 512KB | 2.15$mW$ | 12,479$\mu m^{2}$
Sub-Total |  |  | 20.09$mW$ | 78,794$\mu m^{2}$
| CU | $\times 2048$ | 161.4$W$ | 41.14$mm^{2}$
| SPM | 512MB, 512-bank | 8.1$W$ | 58.7$mm^{2}$
Cryto- | optical NoC | 2 $32\times 16$ bit-sliced | 28.35$W$ | 16.2$mm^{2}$
Light | HBM2 PHY | $\times 2$ | 29.6$W$ | 31.76$mm^{2}$
| TU |  $\times$ 2,569 | 24.2$W$ | 36.8$mm^{2}$
Total |  |  | 146.2$W$ | 290.1$mm^{2}$
Table 2. Our accelerator baselines (norm. to $7$nm). Name | Bit | Description | SPM | Area | Power
---|---|---|---|---|---
Lake (Samardzic:ISCA2022, ) | 32 | 1GHz, 1.4K PEs, 6 HBMs | 256MB | 276$mm^{2}$ | 151$W$
BTS (Kim:ISCA2022, ) | 64 | 1.2GHz, 2K PEs, 2 HBMs | 512MB | 373$mm^{2}$ | 163$W$

###  3.5. Design Overhead

The power and area overhead of CryptoLight is shown in Table 1. All CMOS logic
units are synthesized by Synopsys design compiler with $7$nm PTM process. The
eDRAM SPM and register files are modeled by CACTI. To simulate photonic
microdisk-based computing components, we used Lumerical FDTD (FDTD, ) and
INTERCONNECT (INTERCONNECT, ). To model the electro-optical full adder, we
adopted optical splitters & combiners, photodetectors, and microdisks from
(Ying:JSTQE2018, ). We set the frequency of electro-optical adders and
multipliers to 3GHz. We use two $32\times 16$ optical crossbars to connect all
CUs and SPM banks. We used two HBM2 PHYs to access off-chip DRAMs. Totally,
CryptoLight occupies $290.1mm^{2}$ and consumes $146.2W$.

##  4\. Experimental Methodology

Simulation. We modeled CryptoLight by a cycle-accurate FHE accelerator
simulator, Sapphire-Sim (Banerjee:TCHES2019, ), which is validated against
several crypto-processor chips. The input of Sapphire-Sim is the FHE operation
schedule generated by our FHE compilation framework. Based on the
architectural description of CryptoLight, Sapphire-Sim simulates the cycle-
level execution and data movement for each FHE operation. At last, Sapphire-
Sim generates the total latency and energy consumption of an FHE application.

Schemes. We compared CryptoLight against the state-of-the-art ASIC-based FHE
hardware accelerators, CrateLake (Lake) (Samardzic:ISCA2022, ) and BTS
(Kim:ISCA2022, ). Their configurations are shown in Table 2. Lake has a 32-bit
datapath operating at 1GHz and consisting of 1.4K NTT processing elements
(PEs), six HBM memory controllers, and a 256MB scratchpad memory (SPM) system.
It occupies $276mm^{2}$ at $7$nm and consumes 151 Watt. In contrast, BTS is a
64-bit 1.2GHz FHE accelerator composed of 2K NTT PEs, two HBM memory
controllers, and a 512MB SPM system. Its chip size is $373mm^{2}$ at $7$nm,
and it consumes 163 Watt. Lake uses a larger off-chip DRAM bandwidth to make
up for its smaller SPM capacity.

FHE Benchmarks and Parameters. We first studied the performance and energy of
CKKS FMUL, FROT, and FBOT operations on ciphertexts. We also investigated the
performance and energy of other CKKS operations including
public/private/permutation key generation (KEYGEN), encryption/decryption
(ENC/DEC), addition between two ciphertexts (FADD), addition between a
ciphertext and a plaintext (ADDCP), and multiplication between a ciphertext
and a plaintext (MULTCP). To systematically study FHE performance, we selected
two state-of-the-art CKKS-based applications, i.e., encrypted logistic
regression (Kim:JMI2018, ) (LR) and encrypted neural network inference
(Brutzkus:ICML2019, ) (Lola).

  * •

LR. LR is trained with 2K features. It has 4K samples per batch. Nonlinear
operations in LR are approximated by polynomials. And LR heavily depends on
FMUL, FROT, and FBOT.

  * •

Lola. LoLa is a six-layer convolutional network inferring on the CIFAR-10
dataset. The weights of Lola is not encrypted, and all ReLU activations are
approximated by square functions. So the most intensive operations of Lola are
FROT, MULTCP, and FMUL.

We evaluated the FHE parameters that can support the multiplicative depth of
40 and maintain the 128-bit security, i.e., $N=2^{16}$, and $Q$ is 1536-bit,
since this multiplicative depth is large enough for both FHE applications.
Under this set of FHE parameters, a ciphertext costs $24$MB.

(a) Latency.

(b) Energy.

Figure 6. The latency and energy comparison of FMUL, FROT, and FBOT between
three accelerators (norm. to Lake).

##  5\. Evaluation

In this section, we first report the CKKS operation performance and energy
consumption of CryptoLight, and then discuss the FHE application performance
and power efficiency of CryptoLight. At last, we performed the design space
exploration for CryptoLight.

###  5.1. FHE Operations

FMUL/FROT/FBOT latency. The latency comparison between CryptoLight and various
accelerator baselines is shown in Figure 6(a). Compared to Lake, BTS decreases
the latency of FMUL, FROT and FBOT by 69% on average, due to its larger bit-
width datapath and larger SPM system. The 32-bit datapath of Lake greatly
increases the computational overhead of each key-switching (KS) primitive in
FMUL, FROT and FBOT operations, while the small SPM cannot hold ciphertexts
and KS hint matrices simultaneously, thereby significantly increasing off-chip
memory traffic. Because of the 512-bit EO datapath and the TUs in the SPM,
CryptoLight reduces the latency of FMUL, FROT and FBOT by 96% on average over
BTS. Particularly, the in-SPM TUs of CryptoLight minimize the matrix
transpositions (NTT-T) during (i)NTTs, while the 512-bit datapath of
CryptoLight evenly reduces the latency of the other kernels in each KS
primitive.

FMUL/FROT/FBOT energy. The energy comparison between CryptoLight and various
accelerator baselines is shown in Figure 6(b). Compared to Lake, BTS decreases
the energy consumption of FMUL, FROT and FBOT by 67.5% on average, although it
consumes slightly larger power than Lake, since it uses a much smaller latency
to compute the same operation. This shows that building a 64-bit datapath and
a 512MB SPM actually improves the energy efficiency for FHE accelerators.
Compared to BTS, CryptoLight reduces the energy consumption of FMUL, FROT and
FBOT by 98.8% on average. Instead of power hungry SRAM, CryptoLight uses low-
power eDRAM to construct its SPM. Although operating at higher frequency, the
photonic 512-bit datapath of CryptoLight consumes similar power consumption to
the CMOS 64-bit datapath of BTS. In this way, CryptoLight uses only 89% of the
power consumption of BTS to complete the same FHE operation.

(a) Latency.

(b) Energy.

Figure 7. The latency and energy comparison of other CKKS operations between
three accelerators (norm. to Lake).

Other FHE operations. We also studied the latency and energy improvement of
other FHE operations including KEYGEN, ENC, DEC, FADD, ADDCP, and MULTCP on
CryptoLight in Figure 7. Since these FHE operations do not involve KS, the
latency improvement of these FHE operations achieved by CryptoLight is not as
significant as those of FMUL, FROT and FBOT. On average, CryptoLight reduces
the latency of these FHE operations by $88\%$ over BTS. ENC, ADDCP, and MULCP
still invoke (i)NTT kernels, so their latency improvement is still high, as
shown in Figure 7(a). Moreover, the 512-bit photonic datapath also helps the
FADD operation performing a modular addition between two polynomials. Again,
because other FHE operations has no KS primtive, compared to BTS, as Figure
7(b) highlights, CryptoLight reduces the energy consumption of these FHE
operations by 91% on average.

###  5.2. FHE-based Applications

Latency. The latency comparison of FHE applications between CryptoLight and
various accelerator baselines is shown in Figure 8(a). Since the performance
bottleneck in both LR and Lola is FROT, their latency improvement is heavily
influenced by the speedup of FROT. Compared to Lake, BTS reduces the latency
of LR and Lola by $65.5\%$ and $57.4\%$ respectively, due to its faster FROT
operations. CryptoLight reduces the latency of LR and Lola by $96\%$ and
$93\%$ respectively, compared to BTS. Particularly, the latency improvement of
CryptoLight on Lola is smaller. This is because only square activations on
Lola require FMULs, while the summed products between weights and inputs
involve just MULTCPs.

Energy. The energy comparison of FHE applications between CryptoLight and
various accelerator baselines is shown in Figure 8(b). The energy reduction of
these two applications share a similar trend to the latency reduction.
Compared to Lake, BTS decreases the energy consumption of LR and Lola by
$62.7\%$ and $54\%$ respectively. CryptoLight reduces the energy consumption
of LR and Lola by $96.4\%$ and $93.3\%$ respectively over BTS.

###  5.3. Design Space Exploration

We tried 256-, 512-, and 1K-bit datapaths on our CryptoLight. Compared to
256-bit, the 512-bit CryptoLight reduces the latency of FMUL, FROT and FBOT by
39% on average. However, when the datapath reaches 1K, the latency of these
operations increases by 3%. This is because a longer datapath decreases the
operating frequency of CryptoLight and imposes a larger burden on the small
register file in each CU. We also tried 256MB, 512MB, and 1GB SPMs for
CryptoLight. Compared to 256MB, a 512MB SPM reduces the latency of FMUL, FROT
and FBOT by 27%. However, a 1GB SPM does not reduce their latency obviously.

(a) Latency.

(b) Energy.

Figure 8. The latency and energy of FHE apps (norm. to Lake).

##  6\. Conclusion

In this paper, we propose an EO FHE accelerator, CryptoLight, which has fast
EO 512-bit CUs, a 512MB eDRAM SPM, and in-SPM TUs. CryptoLight greatly reduces
the key-switch cost, while the TUs process matrix transpositions in SPM banks.
On average, for various FHE applications, CryptoLight reduces their latency by
¿$94.4\%$ and their energy consumption by ¿$95\%$ over prior FHE accelerators.

## References

  * (1) U. Banerjee, _et al._ , “Sapphire: A Configurable Crypto-Processor for Post-Quantum Lattice-based Protocols,” _IACR Transactions on Cryptographic Hardware and Embedded Systems_ , (4):17–61, Aug. 2019.
  * (2) A. Brutzkus, _et al._ , “Low Latency Privacy Preserving Inference,” in _International Conference on Machine Learning_ , pages 812–821, 2019.
  * (3) S. Carpov, _et al._ , “Armadillo: A Compilation Chain for Privacy Preserving Applications,” in _3rd International Workshop on Security in Cloud Computing_ , 2015.
  * (4) J. H. Cheon, _et al._ , “Homomorphic encryption for arithmetic of approximate numbers,” in _International Conference on the Theory and Application of Cryptology and Information Security_ , pages 409–437, Springer, 2017.
  * (5) F. de Dinechin, _et al._ , “Pipelined FPGA Adders,” in _International Conference on Field-Programmable Logic and Applications_ , 2010\.
  * (6) C. Gentry, _et al._ , “Homomorphic evaluation of the AES circuit,” in _Annual Cryptology Conference_ , pages 850–867, Springer, 2012.
  * (7) M. Kim, _et al._ , “Secure Logistic Regression Based on Homomorphic Encryption: Design and Evaluation,” _JMIR Medical Informatics_ , 6(2):e19, Apr 2018.
  * (8) S. Kim, _et al._ , “BTS: An Accelerator for Bootstrappable Fully Homomorphic Encryption,” in _ACM International Symposium on Computer Architecture_ , 2022.
  * (9) Lumerical, “FDTD Solutions,” http://www.lumerical.com/
tcad-products/fdtd/.

  * (10) Lumerical, “INTERCONNECT,” http://www.lumerical.com/
tcad-products/interconnect/.

  * (11) S. K. Mathew, _et al._ , “$\mu$RNG: A 300-950 mV, 323 Gbps/W All-Digital Full-Entropy True Random Number Generator in 14nm FinFET CMOS,” _IEEE Journal of Solid-State Circuits_ , 51(7):1695–1704, 2016.
  * (12) P. Mishra, _et al._ , “DELPHI: A cryptographic inference service for neural networks,” in _USENIX Security Symposium_ , 2020.
  * (13) O. Oleksenko, _et al._ , “Varys: Protecting SGX Enclaves from Practical Side-Channel Attacks,” in _USENIX Annual Technical Conference_ , pages 227–240, July 2018.
  * (14) Palisade, “PALISADE Operations using CUDA,” , 2020.
  * (15) M. S. Riazi, _et al._ , “HEAX: An Architecture for Computing on Encrypted Data,” in _Architectural Support for Programming Languages and Operating Systems_ , 2020.
  * (16) N. Samardzic, _et al._ , “CraterLake: A Hardware Accelerator for Efficient Unbounded Computation on Encrypted Data,” in _IEEE/ACM International Symposium on Computer Architecture_ , 2022.
  * (17) Z. Ying, _et al._ , “Electro-Optic Ripple-Carry Adder in Integrated Silicon Photonics for Optical Computing,” _Journal of Selected Topics in Quantum Electronics_ , 2018.


# Global Dynamics of a Discrete Mosquito Model with Allee Effect

Z.S. Boxonov  Z.Boxonov, V.I.Romanovskiy Institute of Mathematics, 9,
University str. 100174, Tashkent, Uzbekistan<EMAIL_ADDRESS>
###### Abstract.

In this paper, we investigate the dynamical behavior of a two-dimensional
discrete mosquito model incorporating an Allee effect. We analyze the system’s
trajectory comprehensively, examining both local and global dynamics. Our
analysis begins by establishing the invariance of the positive trajectory
within the first quadrant. We then investigate the existence of fixed points
within this region. Following this, we perform a stability analysis to
determine both local and global asymptotic stability of the identified fixed
points. Finally, we validate the effectiveness and applicability of our
theoretical results through numerical examples.

###### Key words and phrases:

Discrete mosquito population model; invariant set; cooperative map; local
behavior; global behavior.

###### 2020 Mathematics Subject Classification:

39A12, 39A30, 92D25

##  1\. Introduction

Analyzing the global dynamics and stability of complex systems remains a
significant challenge in discrete dynamical systems. This has spurred recent
interest in the global stability of discrete-time systems, a key area within
dynamic system theory. In this context, we investigate the global dynamics of
a discrete mosquito population model.

Mosquitoes transmit a variety of diseases, including malaria and several
viruses, making their population control a crucial public health practice,
especially in tropical regions.

Mathematical models play a crucial role in understanding mosquito population
dynamics. These models can be continuous-time or discrete-time, as exemplified
by the works in [2]-[5], [13]-[15], [17], and [19].

This work focuses on a discrete-time dynamical system for a simplified
mosquito population model presented in [14]. This model incorporates Allee
effects, which are situations where population growth is limited at low
densities. Mosquitoes undergo complete metamorphosis, with four distinct
development stages: egg, larva, pupa, and adult. To simplify the model, we
combine the three aquatic stages (egg, larva, pupa) into a single larvae class
denoted by $x$. The adult mosquito population is denoted by $y$ [14].
Moreover, we assume that density dependence exists only in the larvae stage.

Denote the birth rate, i.e., the oviposition rate of adults by $\beta(t)$; the
rate of emergence from larvae to adults by a function of the larvae with the
form of $\alpha(1-k(x))$, where $\alpha>0$ is the maximum emergence rate,
$0\leq k(x)\leq 1$, with $k(0)=0,k^{\prime}(x)>0$, and
$\lim\limits_{x\rightarrow\infty}k(x)=1$, is the functional response due to
the intraspecific competition [13]. Moreover, we assume the death rate of
larvae be a linear function, denoted by $d_{0}+d_{1}x$, and the death rate of
adults be constant, denoted by $\mu$. Then, in the absence of sterile
mosquitoes, we get the following system of equations:

$\left\\{\begin{array}[]{ll}\frac{dx}{dt}=\beta(t)y-\alpha(1-k(x))x-(d_{0}+d_{1}x)x,\\\\[8.53581pt]
\frac{dy}{dt}=\alpha(1-k(x))x-\mu y.\\\ \end{array}\right.$ (1.1)

We further assume a functional response for $k(x)$, as in [13], in the form

$k(x)=\frac{x}{1+x}.$

The dynamical system (1.1) was previously studied in [13] for the case where
the birth rate $\beta$ is constant (i.e., $\beta(t)=\beta$ ), in which
mosquito adults have no difficulty finding their mates such that no Allee
effects are concerned. The discrete-time version of this model was then
investigated in [2]-[5].

A component Allee effect describes a decrease in any fitness component (e.g.,
survival, reproduction) as population size or density declines. A prime
example is the decrease in a female’s probability of mating with fewer males
around. This specific phenomenon is commonly referred to as a ”mate-finding
Allee effect” ([8]).

When adult mosquitoes face difficulty finding mates, the model incorporates
Allee effects. This is reflected in the adult birth rate, which is given by:

$\beta(t)=\frac{\beta y(t)}{\gamma+y(t)},$

where $\beta$ is birth rate, $\gamma$ is the Allee effect constant. Then
stage-structured wild mosquito population model is given by:

$\left\\{\begin{array}[]{ll}\frac{dx}{dt}=\frac{\beta
y^{2}}{\gamma+y}-\frac{\alpha x}{1+x}-(d_{0}+d_{1}x)x,\\\\[8.53581pt]
\frac{dy}{dt}=\frac{\alpha x}{1+x}-\mu y.\\\ \end{array}\right.$ (1.2)

This continuous system has been discussed in [14], but its discrete version
has not been investigated as of yet. Discretization of the continuous-time
model offers several advantages for studying mosquito dynamics. It allows
simulations to be conducted on discrete time scales, which are more relevant
to biological processes. This facilitates a closer alignment with actual
mosquito population data and enables the effective integration of time-
sensitive ecological components. Continuous systems can be approached by
several discrete methods, including Euler forward/backward and semi-
discretization. We utilize the Euler forward difference scheme here. In this
paper, we study the discrete time dynamical systems associated with the system
(1.2) by employing Euler’s forward discretization method. Define the operator
$W:{\mathbb{R}^{2}}\rightarrow{\mathbb{R}^{2}}$ by

$\left\\{\begin{array}[]{ll}x^{\prime}=\frac{\beta
y^{2}}{\gamma+y}-\frac{\alpha x}{1+x}-(d_{0}+d_{1}x)x+x,\\\\[8.53581pt]
y^{\prime}=\frac{\alpha x}{1+x}-\mu y+y,\\\ \end{array}\right.$ (1.3)

where $\alpha>0,\beta>0,\gamma>0,\mu>0,\ d_{0}\geq 0,\ d_{1}\geq 0.$ In the
system (1.3), $x^{\prime},y^{\prime}$ means the next state relative to the
initial state $x,y$ respectively.

If for each possible state $\textbf{z}=(x,y)\in\mathbb{R}^{2}$ describing the
current generation, the state
$\textbf{z}^{\prime}=(x^{\prime},y^{\prime})\in\mathbb{R}^{2}$ is uniquely
defined as $\textbf{z}^{\prime}=W(\textbf{z})$ then the map
$W:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ called the evolution operator
[12], [17].

The main problem for a given operator $W$ and arbitrarily initial point
$\textbf{z}^{(0)}=(x^{(0)},y^{(0)})\in\mathbb{R}^{2}$ is to describe the limit
points of the trajectory $\\{\textbf{z}^{(m)}\\}_{m=0}^{\infty}$, where
$\textbf{z}^{(m)}=W^{m}(\textbf{z}^{(0)}).$

In [18] the dynamics of the operator defined by (1.3) is studied under the
condition $d_{0}=d_{1}=0$. In this paper we consider the operator $W$ (defined
by (1.3)) for the case $d_{0}>0$.

##  2\. Preliminaries and Known Results

Consider systems of difference equations of the form

$\left\\{\begin{array}[]{ll}x_{n+1}=f_{1}(x_{n},y_{n}),\\\\[5.69054pt]
y_{n+1}=f_{2}(x_{n},y_{n}),\end{array}n=0,1,2,...\right.$ (2.1)

where $f_{1}$ and $f_{2}$ are given functions and the initial condition
$(x_{0},y_{0})$ comes from some considered set in the intersection of the
domains of $f_{1}$ and $f_{2}$. An operator $W$ that corresponds to the system
(2.1) is defined as $W(x,y)=\left(f_{1}(x,y),f_{2}(x,y)\right)$.

Let $\mathfrak{R}$ be a subset of $\mathbb{R}^{2}$ with a nonempty interior,
and let $W:\mathfrak{R}\rightarrow\mathfrak{R}$ be a continuous operator. When
the function $f_{1}(x,y)$ is increasing in $x$ and decreasing in $y$ and the
function $f_{2}(x,y)$ is decreasing in $x$ and increasing in $y$, the system
(2.1) is called competitive. When the function $f_{1}(x,y)$ is increasing in
$x$ and increasing in $y$ and the function $f_{2}(x,y)$ is increasing in $x$
and increasing in $y$, the system (2.1) is called cooperative. Competitive and
cooperative maps, which are called monotone maps, are defined similarly.
Strongly cooperative systems of difference equations or maps are those for
which the functions $f_{1}$ and $f_{2}$ are coordinate-wise strictly monotone
[1].

Let $\preceq$ be a partial order on $\mathbb{R}^{n}$ with nonnegative cone
$P$. For $z_{1},z_{2}\in\mathbb{R}^{n}$ the order interval
$\mathopen{[\\![}z_{1},z_{2}\mathclose{]\\!]}$ is the set of all $z$ such that
$z_{1}\preceq z\preceq z_{2}$. We say $z_{1}\prec z_{2}$ if $z_{1}\preceq
z_{2}$ and $z_{1}\neq z_{2}$, and $z_{1}\ll z_{2}$ if $z_{2}-z_{1}\in int(P)$.
An operator $W$ on a subset of $\mathbb{R}^{n}$ is order-preserving if
$W(z_{1})\preceq W(z_{2})$ whenever $z_{1}\prec z_{2}$, strictly order-
preserving $W(z_{1})\prec W(z_{2})$ whenever $z_{1}\prec z_{2}$, and strongly
order-preserving if $W(z_{1})\ll W(z_{2})$ whenever $z_{1}\prec z_{2}$.

A point $(x^{*},y^{*})\in\mathbb{R}^{2}$ is called a fixed point of $W$ if
$W\left((x^{*},y^{*})\right)=(x^{*},y^{*})$. The basin of attraction of a
fixed point $(x^{*},y^{*})$ of an operator $W$, is defined as the set of all
initial points $\left(x^{(0)},y^{(0)}\right)$ in the domain of $W$ such that
the sequence of iterates $W^{n}\left(x^{(0)},y^{(0)}\right)$ converges to
$(x^{*},y^{*})$.

The following result is a direct consequence of the Trichotomy Theorem of
Dancer and Hess, (see [6], [10], [11]) proves useful in determining the basins
of attraction for the fixed points.

######  Proposition 1.

If the non-negative cone of a partial ordering $\preceq$ is a generalized
quadrant in $\mathbb{R}^{n}$, and if $W$ has no fixed points in
$\mathopen{[\\![}z_{1},z_{2}\mathclose{]\\!]}$ other than $z_{1}$ and $z_{2}$,
then the interior of $\mathopen{[\\![}z_{1},z_{2}\mathclose{]\\!]}$ is either
a subset of the basin of attraction of $z_{1}$ or a subset of the basin of
attraction of $z_{2}$.

Consider the North-East ordering (NE) on $\mathbb{R}^{2}$ for which the
positive cone is the first quadrant, i.e., this partial ordering is defined by
$(x_{1},y_{1})\preceq_{NE}(x_{2},y_{2})$ if $x_{1}\leq x_{2}$ and $y_{1}\leq
y_{2}$. The South-East (SE) ordering defined as
$(x_{1},y_{1})\preceq_{SE}(x_{2},y_{2})$ if $x_{1}\leq x_{2}$ and $y_{1}\geq
y_{2}$.

A map $W$ on a nonempty set $\mathfrak{R}\subset\mathbb{R}^{2}$ which is
monotone with respect to the $NE$ ordering is called cooperative and a map
monotone with respect to the $SE$ ordering is called competitive.

We define the non-negative quadrant of the plane as
$\mathbb{R}_{+}^{2}=\\{(x,y):x,y\in\mathbb{R},x\geq 0,y\geq 0\\}$. The
operator $W$ is well-defined on $\mathbb{R}^{2}$ except for the lines $x=-1$
and $y=-\gamma$. However, to model a population dynamics system using a
continuous operator, we require both population values (represented by $x$ and
$y$) to be non-negative. Therefore, we restrict our analysis to
$\mathbb{R}^{2}_{+}$ and choose the parameters of the operator $W$ such that
it maps this set onto itself.

######  Lemma 1.

If

$\alpha>0,\ \beta>0,\ \gamma>0,\ 0<\mu\leq 1,\ d_{0}>0,\ \alpha+d_{0}\leq 1,\
d_{1}=0$ (2.2)

then the operator (1.3) maps the set $\mathbb{R}_{+}^{2}$ to itself.

###### Proof.

This can be proved by direct computation. ∎

Lemma 1 provides sufficient conditions for non-negative values of $x$ and $y$.
In this case the system (1.3) becomes

$W_{0}:\left\\{\begin{array}[]{ll}x^{\prime}=\frac{\beta
y^{2}}{\gamma+y}+\left(1-d_{0}-\frac{\alpha}{1+x}\right)x,\\\\[8.53581pt]
y^{\prime}=\frac{\alpha x}{1+x}+(1-\mu)y.\\\ \end{array}\right.$ (2.3)

######  Definition 1.

(see [7]) A fixed point $z$ of an operator $W_{0}$ is called hyperbolic if its
Jacobian $J$ at $z$ has no eigenvalues on the unit circle.

######  Definition 2.

(see [7]) A hyperbolic fixed point $z$ is called:

  * 1)

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

  * 2)

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

  * 3)

a saddle otherwise.

To investigate the local stability analysis of the fixed point of system
(2.3), we present the following lemma [20].

######  Lemma 2.

Let $F(\lambda)=\lambda^{2}+B\lambda+C,$ where $B$ and $C$ are two real
constants. Suppose $\lambda_{1}$ and $\lambda_{2}$ are two roots of
$F(\lambda)=0$. Then the following statements hold.

  1. (i)

If $F(1)>0$ then

  2. (i.1)

$|\lambda_{1}|<1$ and $|\lambda_{2}|<1$ if and only if $F(-1)>0$ and $C<1;$

  3. (i.2)

$\lambda_{1}=-1$ and $\lambda_{2}\neq-1$ if and only if $F(-1)=0$ and $B\neq
2;$

  4. (i.3)

$|\lambda_{1}|<1$ and $|\lambda_{2}|>1$ if and only if $F(-1)<0;$

  5. (i.4)

$|\lambda_{1}|>1$ and $|\lambda_{2}|>1$ if and only if $F(-1)>0$ and $C>1;$

  6. (i.5)

$\lambda_{1}$ and $\lambda_{2}$ are a pair of conjugate complex roots and
$|\lambda_{1}|=|\lambda_{2}|=1$ if only if $-2<B<2$ and $C=1;$

  7. (i.6)

$\lambda_{1}=\lambda_{2}=-1$ if only if $F(-1)=0$ and $B=2.$

  8. (ii)

If $F(1)=0,$ namely, 1 is one root of $F(\lambda)=0,$ then the other root
$\lambda$ satisfies $|\lambda|=(<,>)1$ if and only if $|C|=(<,>)1.$

  9. (iii)

If $F(1)<0$ then $F(\lambda)=0$ has one root lying in $(1;\infty).$ Moreover,

  10. (iii.1)

the other root $\lambda$ satisfies $\lambda<(=)-1$ if and only if
$F(-1)<(=)0;$

  11. (iii.2)

the other root $\lambda$ satisfies $-1<\lambda<1$ if and only if $F(-1)>0.$

##  3\. Analyzing the Existence and Local Stability of Fixed Points

We begin by investigating the existence of fixed points.

The fixed points of the system (2.3) satisfy the following system of
equations:

$\begin{array}[]{ll}x=\frac{\beta
y^{2}}{\gamma+y}+\left(1-d_{0}-\frac{\alpha}{1+x}\right)x,\\\\[8.53581pt]
y=\frac{\alpha x}{1+x}+(1-\mu)y.\\\ \end{array}$ (3.1)

It can be seen that the point $(0,0)$ is a solution to system (3.1). By
solving the second equation of system (3.1) for $y$ and substituting it into
the first equation, we obtain the following equation:

$d_{0}\mu(\alpha+\gamma\mu)x^{2}+\left(\mu\left(d_{0}\gamma\mu+(\alpha+d_{0})(\alpha+\gamma\mu)\right)-\alpha^{2}\beta\right)x+\gamma\mu^{2}(\alpha+d_{0})=0.$
(3.2)

We denote

$\begin{array}[]{ll}D=\left(\mu(d_{0}\gamma\mu+(\alpha+d_{0})(\alpha+\gamma\mu))-\alpha^{2}\beta\right)^{2}-4d_{0}\gamma\mu^{3}(\alpha+d_{0})(\alpha+\gamma\mu),\\\\[8.53581pt]
\nu=\frac{\mu}{\alpha^{2}}\left(\sqrt{d_{0}\gamma\mu}+\sqrt{(\alpha+d_{0})(\alpha+\gamma\mu)}\right)^{2}.\end{array}$
(3.3)

Conditions in (2.2) ensure $d_{0}\mu(\alpha+\gamma\mu)>0$ and
$\gamma\mu^{2}(\alpha+d_{0})>0$. For positive roots in equation (3.2), we
additionally require
$\mu\left(d_{0}\gamma\mu+(\alpha+d_{0})(\alpha+\gamma\mu)\right)-\alpha^{2}\beta<0$
and $D>0$. From here we get $\beta>\nu.$ The number of solutions to equation
(3.2) depends on the relationship between $\beta$ and $\nu$. If $\beta=\nu$,
equation (3.2) has one solution. If $\beta<\nu$, equation (3.2) has no
solution.

The following proposition holds for fixed points of the operator $W_{0}$:

######  Proposition 2.

The fixed points for (2.3) are as follows:

  * •

If $\beta<\nu$ then the operator (2.3) has a unique fixed point $z_{1}=(0,0).$

  * •

If $\beta=\nu$ then mapping (2.3) has two fixed points

$z_{1,1}=(0,0),\ \ z_{1,2}=\left(x^{*}_{1},y^{*}_{1}\right).$

  * •

If $\beta>\nu$ then mapping (2.3) has three fixed points

$z_{2,1}=(0,0),\ \ z_{2,2}=\left(x^{*}_{2},y^{*}_{2}\right),\ \
z_{2,3}=\left(x^{*}_{3},y^{*}_{3}\right),$

where

$\begin{array}[]{llll}x^{*}_{1}=\sqrt{\frac{\gamma\mu(\alpha+d_{0})}{d_{0}(\alpha+\gamma\mu)}},\\\\[8.53581pt]
x^{*}_{2}=\frac{\alpha^{2}\beta-\mu(d_{0}\gamma\mu+(\alpha+d_{0})(\alpha+\gamma\mu))-\sqrt{D}}{2d_{0}\mu(\alpha+\gamma\mu)},\\\\[8.53581pt]
x^{*}_{3}=\frac{\alpha^{2}\beta-\mu(d_{0}\gamma\mu+(\alpha+d_{0})(\alpha+\gamma\mu))+\sqrt{D}}{2d_{0}\mu(\alpha+\gamma\mu)},\\\\[8.53581pt]
y^{*}_{i}=\frac{\alpha x^{*}_{i}}{\mu(1+x^{*}_{i})},\ i=1,2,3.\end{array}$
(3.4)

Now we shall examine the type of the fixed point.

To determine the type of a fixed point $z=(x^{*},y^{*})$ of the operator
defined by (2.3), we can analyze the Jacobian matrix:

$J_{W_{0}}=\left(\begin{array}[]{cc}1-d_{0}-\frac{\alpha}{(1+x)^{2}}&\frac{\beta
y(2\gamma+y)}{(\gamma+y)^{2}}\\\ \frac{\alpha}{(1+x)^{2}}&1-\mu\\\
\end{array}\right).$

The characteristic equation is

$F(\lambda)=\lambda^{2}+B\lambda+C,$ (3.5)

where

$B=\mu+d_{0}+\frac{\alpha}{(1+x^{*})^{2}}-2,$ $C=1-\mu-
d_{0}-\frac{\alpha}{(1+x^{*})^{2}}+\mu
d_{0}+\frac{\alpha}{(1+x^{*})^{2}}\left(\mu-\frac{\beta
y^{*}(2\gamma+y^{*})}{(\gamma+y^{*})^{2}}\right).$

From (3.5) we get

$F(1)=1+B+C=\mu\left(d_{0}+\frac{\alpha}{(1+x^{*})^{2}}\right)-\frac{\alpha}{(1+x^{*})^{2}}\cdot\frac{\beta
y^{*}}{\gamma+y^{*}}\left(1+\frac{\gamma}{\gamma+y^{*}}\right),$
$F(-1)=1-B+C=4-2\left(\mu+d_{0}+\frac{\alpha}{(1+x^{*})^{2}}\right)+\mu\left(d_{0}+\frac{\alpha}{(1+x^{*})^{2}}\right)-$
$\frac{\alpha}{(1+x^{*})^{2}}\cdot\frac{\beta
y^{*}}{\gamma+y^{*}}\left(1+\frac{\gamma}{\gamma+y^{*}}\right).$

(1) Let $\beta<\nu$, i.e., $x^{*}=0,y^{*}=0.$ Then

$F(1)=\mu(\alpha+d_{0})>0,\ C=(1-d_{0}-\alpha)(1-\mu)<1,$ $F(-1)=2(2-\mu-
d_{0}-\alpha)+\mu(d_{0}+\alpha)>0.$

According to item (i.1) of Lemma 2, $|\lambda_{1,2}|<1.$ Therefore, the fixed
point $(0,0)$is attracting.

(2) Let $\beta=\nu$, i.e.,
$x^{*}=x_{1}^{*}=\sqrt{\frac{\gamma\mu(\alpha+d_{0})}{d_{0}(\alpha+\gamma\mu)}}.$
Then

$F(1)=\mu\left(d_{0}+\frac{\alpha}{(1+x^{*})^{2}}\right)-\frac{\alpha}{(1+x^{*})^{2}}\cdot\frac{\beta
y^{*}}{\gamma+y^{*}}\left(1+\frac{\gamma}{\gamma+y^{*}}\right)$
$=\mu\left(d_{0}+\frac{\alpha}{(1+x^{*})^{2}}\right)-\frac{\mu
y^{*}}{(1+x^{*})x^{*}}\cdot\frac{\beta
y^{*}}{\gamma+y^{*}}\left(1+\frac{\gamma}{\gamma+y^{*}}\right)$
$=\frac{\mu}{(1+x^{*})x^{*}}\left(d_{0}x^{*}(1+x^{*})+\frac{\alpha
x^{*}}{1+x^{*}}\right)-\frac{\mu y^{*}}{(1+x^{*})x^{*}}\cdot\frac{\beta
y^{*}}{\gamma+y^{*}}\left(1+\frac{\gamma}{\gamma+y^{*}}\right)$
$=\frac{\mu}{(1+x^{*})x^{*}}\left(d_{0}x{{}^{*}}^{2}+\frac{\beta
y{{}^{*}}^{2}}{\gamma+y^{*}}\right)-\frac{\mu
y^{*}}{(1+x^{*})x^{*}}\cdot\frac{\beta
y^{*}}{\gamma+y^{*}}\left(1+\frac{\gamma}{\gamma+y^{*}}\right)$
$=\frac{\mu}{(1+x^{*})x^{*}}\left(d_{0}x{{}^{*}}^{2}-\frac{\beta\gamma
y{{}^{*}}^{2}}{(\gamma+y^{*})^{2}}\right)=\frac{\mu}{(1+x^{*})x^{*}}\left(d_{0}x{{}^{*}}^{2}-\frac{\beta\gamma\alpha^{2}x{{}^{*}}^{2}}{(\gamma\mu+(\alpha+\gamma\mu)x^{*})^{2}}\right)$
$=\frac{\mu
x{{}^{*}}^{2}}{(1+x^{*})x^{*}}\left(d_{0}-\frac{\beta\gamma\alpha^{2}}{(\gamma\mu+(\alpha+\gamma\mu)x^{*})^{2}}\right)$
$=\frac{\mu
x^{*}}{1+x^{*}}\left(d_{0}-\frac{\gamma\alpha^{2}\mu\left(\frac{\sqrt{d_{0}\gamma\mu}+\sqrt{(\alpha+d_{0})(\alpha+\gamma\mu)}}{\alpha}\right)^{2}}{\left(\gamma\mu+(\alpha+\gamma\mu)\sqrt{\frac{\gamma\mu(\alpha+d_{0})}{d_{0}(\alpha+\gamma\mu)}}\right)^{2}}\right)=\frac{\mu
x^{*}}{1+x^{*}}\left(d_{0}-d_{0}\right)=0,$ $-1<C=1-\mu-
d_{0}-\frac{\alpha}{(1+x^{*})^{2}}+F(1)=1-\mu-
d_{0}-\frac{\alpha}{(1+x^{*})^{2}}<1.$

According to item (ii) of Lemma 2, $\lambda_{1}=1$, $|\lambda_{2}|<1.$ So, the
fixed point $(x_{1}^{*},y_{1}^{*})$ is non-hyperbolic.

(3) Let $\beta>\nu$, i.e., $x^{*}=x_{2,3}^{*}$. Then

$F(1)=\frac{\mu
x{{}^{*}}^{2}}{(1+x^{*})x^{*}}\left(d_{0}-\frac{\beta\gamma\alpha^{2}}{(\gamma\mu+(\alpha+\gamma\mu)x^{*})^{2}}\right)$
$=\frac{\mu
x^{*}}{1+x^{*}}\left(d_{0}-\frac{\beta\gamma\alpha^{2}}{\left(\gamma\mu+(\alpha+\gamma\mu)\cdot\frac{\alpha^{2}\beta-\mu(d_{0}\gamma\mu+(\alpha+d_{0})(\alpha+\gamma\mu))\pm\sqrt{D}}{2d_{0}\mu(\alpha+\gamma\mu)}\right)^{2}}\right)$
$=\frac{\mu
d_{0}x^{*}}{1+x^{*}}\cdot\frac{D\pm\sqrt{D^{2}+4\alpha^{2}\beta\gamma\mu^{2}d_{0}D}}{D+2\alpha^{2}\beta\gamma\mu^{2}d_{0}\pm\sqrt{D^{2}+4\alpha^{2}\beta\gamma\mu^{2}d_{0}}}.$

If $x^{*}=x_{3}^{*}$, then $F(1)>0$, $F(-1)>0$, $C<1.$ According to item (i1)
of Lemma 2, $|\lambda_{1,2}|<1.$ Therefore, the fixed point
$(x_{3}^{*},y_{3}^{*})$ is attracting. If $x^{*}=x_{2}^{*}$, then $F(1)<0$,
$F(-1)<0.$ According to item (iii1) of Lemma 2, $\lambda_{1}>1$,
$\lambda_{2}<-1.$ So, the fixed point $(x_{2}^{*},y_{2}^{*})$ is repelling.

Thus for the type of fixed points the following theorem holds.

######  Theorem 1.

The type of the fixed points for (2.3) are as follows:

  * i)

if $\beta<\nu$, then the operator (2.3) has unique fixed point $(0,0)$ which
is attracting.

  * ii)

if $\beta=\nu$, then the operator (2.3) has two fixed points $(0,0)$,
$(x_{1}^{*},y_{1}^{*})$, and the point $(0,0)$ is attracting, the point
$(x_{1}^{*},y_{1}^{*})$ is non-hyperbolic (but $(x_{1}^{*},y_{1}^{*})$ is
semi-attracting 111meaning that the second eigenvalue is less than 1 in
absolute value.).

  * iii)

if $\beta>\nu$, then the operator (2.3) has three fixed points $(0,0)$,
$(x_{2}^{*},y_{2}^{*})$, $(x_{3}^{*},y_{3}^{*})$, and the points $(0,0)$ and
$(x_{3}^{*},y_{3}^{*})$ are attracting, the point $(x_{2}^{*},y_{2}^{*})$ is
repelling.

where $(x_{1}^{*},y_{1}^{*})$, $(x_{2}^{*},y_{2}^{*})$ and
$(x_{3}^{*},y_{3}^{*})$ are defined by (3.4).

##  4\. Global Dynamics of the Operator $W_{0}$

In this section for any initial point
$\left(x^{(0)},y^{(0)}\right)\in\mathbb{R}^{2}_{+}$ we investigate behavior of
the trajectories
$\left(x^{(n)},y^{(n)}\right)=W_{0}^{n}\left(x^{(0)},y^{(0)}\right),n\geq 1.$

Let

$\Omega=[0,\omega_{1}]\times[0,\omega_{2}]$ (4.1)

with

$\omega_{1}=\frac{\alpha^{2}\beta}{\mu d_{0}(\alpha+\gamma\mu)},\
\omega_{2}=\frac{\alpha}{\mu}.$

######  Definition 3.

(see [16]) The set $A$ is called the absorbing set for the operator $W_{0}$ if
it satisfies the following properties:

  * 1)

Invariance: the set $A$ is invariant with respect to $W_{0}$, i.e.,
$W_{0}(A)\subset A$;

  * 2)

Eventual absorption: for any initial point
$\left(x^{(0)},y^{(0)}\right)\in\mathbb{R}^{2}_{+}$, there exists
$k_{0}\in\mathbb{N}$ such that $W_{0}^{k}\left(x^{(0)},y^{(0)}\right)\in A$
for all $k\geq k_{0}.$

######  Lemma 3.

(Absorbing set).The set $\Omega$ is an absorbing set for the operator $W_{0}$
given by (2.3).

###### Proof.

We have

$\begin{array}[]{ll}x^{(n)}=\frac{\beta\left(y^{(n-1)}\right)^{2}}{\gamma+y^{(n-1)}}+\left(1-d_{0}-\frac{\alpha}{1+x^{(n-1)}}\right)x^{(n-1)},\\\\[5.69054pt]
y^{(n)}=\frac{\alpha
x^{(n-1)}}{1+x^{(n-1)}}+\left(1-\mu\right)y^{(n-1)}.\end{array}$ (4.2)

Assume that all values of $y^{(n)}$ are greater than $\omega_{2}$. Then

$y^{(n+1)}-y^{(n)}=\frac{\alpha x^{(n)}}{1+x^{(n)}}-\mu y^{(n)}<\alpha-\mu
y^{(n)}=\mu\left(\omega_{2}-y^{(n)}\right)<0.$

So $y^{(n)}$ is a decreasing sequence. Since $y^{(n)}$ is decreasing and
bounded from below we have:

$\lim_{n\to\infty}y^{(n)}\geq\omega_{2}.$ (4.3)

We estimate $y^{(n)}$ by the following:

$y^{(n)}=\frac{\alpha
x^{(n-1)}}{1+x^{(n-1)}}+(1-\mu)y^{(n-1)}<\alpha+(1-\mu)y^{(n-1)}<\alpha+(1-\mu)\left(\alpha+(1-\mu)y^{(n-2)}\right)$
$<\alpha+\alpha(1-\mu)+(1-\mu)^{2}\left(\alpha+(1-\mu)y^{(n-3)}\right)<...$
$<\alpha+\alpha(1-\mu)+\alpha(1-\mu)^{2}+...+\alpha(1-\mu)^{n-1}+(1-\mu)^{n}y^{(0)}=\omega_{2}+(1-\mu)^{n}\left(y^{(0)}-\omega_{2}\right).$

Thus $y^{(n)}<\omega_{2}+(1-\mu)^{n}\left(y^{(0)}-\omega_{2}\right)$.
Consequently

$\lim_{n\to\infty}y^{(n)}\leq\omega_{2}.$ (4.4)

By (4.3) and (4.4) we have

$\lim\limits_{n\to\infty}y^{(n)}=\omega_{2}.$ (4.5)

From (4.5) and $y^{(n)}=\frac{\alpha x^{(n-1)}}{1+x^{(n-1)}}-\mu
y^{(n-1)}+y^{(n-1)}$ it follows

$\lim\limits_{n\to\infty}x^{(n)}=+\infty.$ (4.6)

Now we estimate the sequence $x^{(n)}$. The inequality $y^{(n)}\leq y^{(0)}$
holds for all $n$, as the sequence $y^{(n)}$ is decreasing. We denote
$y_{0}^{(0)}=\frac{\beta\left(y^{(0)}\right)^{2}}{\gamma+y^{(0)}}.$ Then

$x^{(n)}=\frac{\beta\left(y^{(n-1)}\right)^{2}}{\gamma+y^{(n-1)}}+\left(1-d_{0}-\frac{\alpha}{1+x^{(n-1)}}\right)x^{(n-1)}<y_{0}^{(0)}+(1-d_{0})x^{(n-1)}$
$<y_{0}^{(0)}+(1-d_{0})\left(y_{0}^{(0)}+(1-d_{0})x^{(n-2)}\right)<...$
$<y_{0}^{(0)}+y_{0}^{(0)}(1-d_{0})+...+y_{0}^{(0)}(1-d_{0})^{n-1}+(1-d_{0})^{n}x^{(0)}$
$=\frac{y_{0}^{(0)}}{d_{0}}+(1-d_{0})^{n}\left(x^{(0)}-\frac{y_{0}^{(0)}}{d_{0}}\right).$

From condition (2.2) it follows $d_{0}\in(0,1).$ Consequently

$\lim_{n\to\infty}x^{(n)}\leq\frac{y_{0}^{(0)}}{d_{0}}.$ (4.7)

Our initial assumption is invalid based on (4.6) and (4.7). Therefore, there
exists $n_{0}$ such that $y^{(n_{0})}$ is not greater than $\omega_{2}$. If
$y^{(n_{0})}\leq\omega_{2}$ then $y^{(n_{0}+1)}<\omega_{2}$. Indeed,

$y^{(n_{0}+1)}=(1-\mu)y^{(n_{0})}+\frac{\alpha
x^{(n_{0})}}{1+x^{(n_{0})}}\leq(1-\mu)\omega_{2}+\frac{\alpha
x^{(n_{0})}}{1+x^{(n_{0})}}$
$=\omega_{2}-\alpha\left(1-\frac{x^{(n_{0})}}{1+x^{(n_{0})}}\right)<\omega_{2}.$

Assume that all values of $x^{(n)}$ are greater than $\omega_{1}$. Then

$x^{(n+1)}-x^{(n)}=\frac{\beta\left(y^{(n)}\right)^{2}}{\gamma+y^{(n)}}-\frac{\alpha
x^{(n)}}{1+x^{(n)}}-d_{0}x^{(n)}<\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}-d_{0}x^{(n)}=d_{0}\left(\omega_{1}-x^{(n)}\right)\leq
0.$

So $x^{(n)}$ is a decreasing sequence. Since $x^{(n)}$ is decreasing and
bounded from below we have:

$\lim_{n\to\infty}x^{(n)}\geq\omega_{1}.$ (4.8)

We estimate $x^{(n)}$ by the following:

$x^{(n)}=\frac{\beta\left(y^{(n-1)}\right)^{2}}{\gamma+y^{(n-1)}}+\left(1-d_{0}-\frac{\alpha}{1+x^{(n-1)}}\right)x^{(n-1)}<\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}+(1-d_{0})x^{(n-1)}$
$<\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}+(1-d_{0})\left(\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}+(1-d_{0})x^{(n-2)}\right)<...$
$<\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}+\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}(1-d_{0})+...+\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}(1-d_{0})^{n-1}+(1-d_{0})^{n}x^{(0)}$
$=\omega_{1}+(1-d_{0})^{n}\left(x^{(0)}-\omega_{1}\right).$

From $d_{0}\in(0,1)$, consequently, the inequality holds.

$\lim_{n\to\infty}x^{(n)}\leq\omega_{1}.$ (4.9)

By (4.8) and (4.9) we have

$\lim\limits_{n\to\infty}x^{(n)}=\omega_{1}.$ (4.10)

Since $x^{(n)}$ converges to $\omega_{1}$, it follows from the second equation
of (4.2) that

$\lim\limits_{n\to\infty}\left(y^{(n)}-(1-\mu)y^{(n-1)}\right)=\frac{\alpha\omega_{1}}{1+\omega_{1}}.$

Since $y^{(n)}$ is bounded, it possesses upper and lower limits, denoted by
$\overline{y}$ and $\underline{y}$, respectively. To arrive at a
contradiction, suppose $y^{(n)}$ has no limit
(i.e.,$\overline{y}\neq\underline{y}$). This would necessitate the existence
of two subsequences: $y^{(n_{i})}$ converging to $\overline{y}$ and
$y^{(n_{j})}$ converging to $\underline{y}$.

$\lim\limits_{i\to\infty}\left(y^{(n_{i})}-(1-\mu)y^{(n_{i}-1)}\right)=\frac{\alpha\omega_{1}}{1+\omega_{1}},\
\Rightarrow\
\lim\limits_{i\to\infty}y^{(n_{i}-1)}=\frac{1}{1-\mu}\left(\overline{y}-\frac{\alpha\omega_{1}}{1+\omega_{1}}\right)\leq\overline{y},\
\Rightarrow\ \overline{y}\leq\frac{\alpha\omega_{1}}{\mu(1+\omega_{1})};$

$\lim\limits_{i\to\infty}\left(y^{(n_{i}+1)}-(1-\mu)y^{(n_{i})}\right)=\frac{\alpha\omega_{1}}{1+\omega_{1}},\
\Rightarrow\
\lim\limits_{i\to\infty}y^{(n_{i}+1)}=(1-\mu)\overline{y}+\frac{\alpha\omega_{1}}{1+\omega_{1}}\leq\overline{y},\
\Rightarrow\ \overline{y}\geq\frac{\alpha\omega_{1}}{\mu(1+\omega_{1})};$

$\lim\limits_{j\to\infty}\left(y^{(n_{j})}-(1-\mu)y^{(n_{j}-1)}\right)=\frac{\alpha\omega_{1}}{1+\omega_{1}},\
\Rightarrow\
\lim\limits_{j\to\infty}y^{(n_{j}-1)}=\frac{1}{1-\mu}\left(\underline{y}-\frac{\alpha\omega_{1}}{1+\omega_{1}}\right)\geq\underline{y},\Rightarrow\underline{y}\geq\frac{\alpha\omega_{1}}{\mu(1+\omega_{1})};$

$\lim\limits_{j\to\infty}\left(y^{(n_{j}+1)}-(1-\mu)y^{(n_{j})}\right)=\frac{\alpha\omega_{1}}{1+\omega_{1}},\
\Rightarrow\
\lim\limits_{j\to\infty}y^{(n_{j}+1)}=(1-\mu)\underline{y}+\frac{\alpha\omega_{1}}{1+\omega_{1}}\geq\underline{y},\
\Rightarrow\ \underline{y}\leq\frac{\alpha\omega_{1}}{\mu(1+\omega 1)}.$

Therefore, we conclude that
$\overline{y}=\underline{y}=\frac{\alpha\omega_{1}}{\mu(1+\omega_{1})}$ i.e.,

$\lim\limits_{n\to\infty}y^{(n)}=\frac{\omega_{1}\omega_{2}}{1+\omega_{1}}.$

Proposition 2 shows that
$\left(\omega_{1},\frac{\omega_{1}\omega_{2}}{1+\omega_{1}}\right)$ is not a
fixed point of system (2.3), contradicting our hypothesis. Therefore, there
exists $m_{0}$ such that $x^{(k_{0})}$ is not greater than $\omega_{1}$. If
$x^{(k_{0})}\leq\omega_{1}$ then $x^{(k_{0}+1)}<\omega_{1}$. Indeed,

$x^{(k_{0}+1)}=\frac{\beta\left(y^{(k_{0})}\right)^{2}}{\gamma+y^{(k_{0})}}+\left(1-d_{0}-\frac{\alpha}{1+x^{(k_{0})}}\right)x^{(k_{0})}<\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}+(1-d_{0})x^{(k_{0})}$
$\leq\frac{\beta\omega_{2}^{2}}{\gamma+\omega_{2}}+(1-d_{0})\omega_{1}=\omega_{1}.$

So for any initial point $\left(x^{(0)},y^{(0)}\right)\in\mathbb{R}^{2}_{+}$,
there exists $k_{0}\in\mathbb{N}$ such that
$W_{0}^{k}\left(x^{(0)},y^{(0)}\right)\in\Omega$ for all $k\geq k_{0}.$
Furthermore, $W_{0}$ maps the set $\Omega$ into itself. ∎

Lemma 3 establishes that $\Omega$ is an absorbing set. Consequently, we can
restrict our analysis to the dynamics of the operator $W_{0}$ within the set
$\Omega$.

######  Proposition 3.

Let $[a,b]$ and $[c,d]$ be intervals of real numbers, and let

$f_{1}:[a,b]\times[c,d]\rightarrow[a,b],\ \ and\ \
f_{2}:[a,b]\times[c,d]\rightarrow[c,d]$

be continuous functions. Consider the system of difference equation

$\begin{array}[]{ll}x_{n+1}=f_{1}(x_{n},y_{n}),\\\\[5.69054pt]
y_{n+1}=f_{2}(x_{n},y_{n}),\end{array}\ \ n=0,1,2,...$ (4.11)

with initial condition $\left(x^{(0)},y^{(0)}\right)\in[a,b]\times[c,d].$
Suppose the following conditions (i) and (ii) are satisfied:

  * (i)

$f_{1}(x,y)$ and $f_{2}(x,y)$ are non-decreasing in both $x$ and $y.$

  * (ii)

If $a_{1},b_{1}\in[a,b]$, $c_{1},d_{1}\in[c,d]$ is a solution of the system of
equations

$\left\\{\begin{array}[]{ll}a_{1}=f_{1}(a_{1},c_{1}),\\\\[5.69054pt]
c_{1}=f_{2}(a_{1},c_{1}),\end{array}\right.\left\\{\begin{array}[]{ll}b_{1}=f_{1}(b_{1},d_{1})\\\\[5.69054pt]
d_{1}=f_{2}(b_{1},d_{1})\end{array}\right.$ (4.12)

then $a_{1}=b_{1}$ and $c_{1}=d_{1}$.

Then there exists exactly one fixed point $(x^{*},y^{*})$ of system (4.11),
and trajectory of any initial point $\left(x^{(0)},y^{(0)}\right)$ in
$[a,b]\times[c,d]$ converges to $(x^{*},y^{*})$.

###### Proof.

By the Brouwer fixed-point theorem, the operator

$W:[a,b]\times[c,d]\rightarrow[a,b]\times[c,d]$

defined by

$W(x,y)=\left(f_{1}(x,y),f_{2}(x,y)\right)$

has a fixed point $(x^{*},y^{*})$ (cf. with proof of Theorem 1.16 of [9]).
Since $(x^{*},y^{*})$ also satisfies system (4.11), we can conclude that
system (4.11) has at least one fixed point.

Consider the sequence $\left(x^{(n)},y^{(n)}\right)$ generated by the operator
$W$ for any initial point $\left(x^{(0)},y^{(0)}\right)$ within the rectangle
$[a,b]\times[c,d]$. It suffices to show that

$\lim\limits_{n\to\infty}\left(x^{(n)},y^{(n)}\right)=(x^{*},y^{*}).$

We denote $a_{1}^{(0)}=a,b_{1}^{(0)}=b,c_{1}^{(0)}=c,d_{1}^{(0)}=d$, and for
each $i\geq 0$, let

$\begin{array}[]{ll}a_{1}^{(i+1)}=f_{1}\left(a_{1}^{(i)},c_{1}^{(i)}\right),\
\ b_{1}^{(i+1)}=f_{1}\left(b_{1}^{(i)},d_{1}^{(i)}\right),\\\\[8.53581pt]
c_{1}^{(i+1)}=f_{2}\left(a_{1}^{(i)},c_{1}^{(i)}\right),\ \
d_{1}^{(i+1)}=f_{2}\left(b_{1}^{(i)},d_{1}^{(i)}\right).\end{array}$ (4.13)

Then

$\begin{array}[]{ll}a_{1}^{(0)}=a\leq
f_{1}\left(a_{1}^{(0)},c_{1}^{(0)}\right)\leq
f_{1}\left(b_{1}^{(0)},d_{1}^{(0)}\right)\leq b=b_{1}^{(0)},\\\\[8.53581pt]
c_{1}^{(0)}=c\leq f_{2}\left(a_{1}^{(0)},c_{1}^{(0)}\right)\leq
f_{2}\left(b_{1}^{(0)},d_{1}^{(0)}\right)\leq d=d_{1}^{(0)},\end{array}$

and from this, we can see that

$a_{1}^{(0)}\leq a_{1}^{(1)}\leq b_{1}^{(1)}\leq b_{1}^{(0)}$ and
$c_{1}^{(0)}\leq c_{1}^{(1)}\leq d_{1}^{(1)}\leq d_{1}^{(0)}.$

Following the same approach, we have

$\begin{array}[]{ll}a_{1}^{(1)}=f_{1}\left(a_{1}^{(0)},c_{1}^{(0)}\right)\leq
f_{1}\left(a_{1}^{(1)},c_{1}^{(1)}\right)\leq
f_{1}\left(b_{1}^{(1)},d_{1}^{(1)}\right)\leq
f_{1}\left(b_{1}^{(0)},d_{1}^{(0)}\right)=b_{1}^{(1)},\\\\[8.53581pt]
c_{1}^{(1)}=f_{2}\left(a_{1}^{(0)},c_{1}^{(0)}\right)\leq
f_{2}\left(a_{1}^{(1)},c_{1}^{(1)}\right)\leq
f_{2}\left(b_{1}^{(1)},d_{1}^{(1)}\right)\leq
f_{2}\left(b_{1}^{(0)},d_{1}^{(0)}\right)=d_{1}^{(1)},\end{array}$

and from this, we can see that

$a_{1}^{(0)}\leq a_{1}^{(1)}\leq a_{1}^{(2)}\leq b_{1}^{(2)}\leq
b_{1}^{(1)}\leq b_{1}^{(0)}$ and $c_{1}^{(0)}\leq c_{1}^{(1)}\leq
c_{1}^{(2)}\leq d_{1}^{(2)}\leq d_{1}^{(1)}\leq d_{1}^{(0)}.$

Note that

$a_{1}^{(0)}=a\leq x^{(n)}\leq b=b_{1}^{(0)}$ and $c_{1}^{(0)}=c\leq
y^{(n)}\leq d=d_{1}^{(0)}$ for all $n\geq 0.$

For all $n\geq 0$, we have

$\begin{array}[]{ll}a_{1}^{(1)}=f_{1}\left(a_{1}^{(0)},c_{1}^{(0)}\right)\leq
f_{1}\left(x^{(n)},y^{(n)}\right)\leq
f_{1}\left(b_{1}^{(0)},d_{1}^{(0)}\right)=b_{1}^{(1)},\\\\[8.53581pt]
c_{1}^{(1)}=f_{2}\left(a_{1}^{(0)},c_{1}^{(0)}\right)\leq
f_{2}\left(x^{(n)},y^{(n)}\right)\leq
f_{2}\left(b_{1}^{(0)},d_{1}^{(0)}\right)=d_{1}^{(1)}.\end{array}$

Hence, it follows that

$a_{1}^{(1)}\leq x^{(n)}\leq b_{1}^{(1)}$ and $c_{1}^{(1)}\leq y^{(n)}\leq
d_{1}^{(1)}$ for all $n\geq 1.$

For all $n\geq 1$, we have

$\begin{array}[]{ll}a_{1}^{(2)}=f_{1}\left(a_{1}^{(1)},c_{1}^{(1)}\right)\leq
f_{1}\left(x^{(n)},y^{(n)}\right)\leq
f_{1}\left(b_{1}^{(1)},d_{1}^{(1)}\right)=b_{1}^{(2)},\\\\[8.53581pt]
c_{1}^{(2)}=f_{2}\left(a_{1}^{(1)},c_{1}^{(1)}\right)\leq
f_{2}\left(x^{(n)},y^{(n)}\right)\leq
f_{2}\left(b_{1}^{(1)},d_{1}^{(1)}\right)=d_{1}^{(2)}.\end{array}$

This implies that

$a_{1}^{(2)}\leq x^{(n)}\leq b_{1}^{(2)}$ and $c_{1}^{(2)}\leq y^{(n)}\leq
d_{1}^{(2)}$ for all $n\geq 2.$

It follows by induction that for $i\geq 0$, the following statements are true:

(i) $a=a_{1}^{(0)}\leq a_{1}^{(1)}\leq...\leq a_{1}^{(i-1)}\leq
a_{1}^{(i)}\leq b_{1}^{(i)}\leq b_{1}^{(i-1)}\leq...\leq b_{1}^{(1)}\leq
b_{1}^{(0)}=b.$

(ii) $c=c_{1}^{(0)}\leq c_{1}^{(1)}\leq...\leq c_{1}^{(i-1)}\leq
c_{1}^{(i)}\leq d_{1}^{(i)}\leq d_{1}^{(i-1)}\leq...\leq d_{1}^{(1)}\leq
d_{1}^{(0)}=d.$

(iii) $a_{1}^{(i)}\leq x^{(n)}\leq b_{1}^{(i)}$ for all $n\geq i.$

(iv) $c_{1}^{(i)}\leq y^{(n)}\leq d_{1}^{(i)}$ for all $n\geq i.$

Assertions (i) and (ii) guarantee that the sequences $a_{1}^{(i)}$,
$b_{1}^{(i)}$, $c_{1}^{(i)}$ and $d_{1}^{(i)}$ converge. Let their limits be
denoted by $a_{1}$, $b_{1}$, $c_{1}$, and $d_{1}$, respectively. Then

$a\leq a_{1}\leq b_{1}\leq b$ and $c\leq c_{1}\leq d_{1}\leq d.$

By the continuity of $f_{1}$ and $f_{2}$, we have

$\left\\{\begin{array}[]{ll}a_{1}=f_{1}(a_{1},c_{1}),\\\\[5.69054pt]
c_{1}=f_{2}(a_{1},c_{1}),\end{array}\right.\left\\{\begin{array}[]{ll}b_{1}=f_{1}(b_{1},d_{1})\\\\[5.69054pt]
d_{1}=f_{2}(b_{1},d_{1})\end{array}\right.$

and hence $a_{1}=b_{1}$ and $c_{1}=d_{1}$ from which the proof follows. ∎

The following theorem describes the trajectory of any initial point
$(x^{(0)},y^{(0)})$ in $\Omega$.

######  Theorem 2.

Assume that (2.2) holds and let $\beta<\nu$. Then the trajectory converges to
the unique fixed point zero for any initial point
$\left(x^{(0)},y^{(0)}\right)\in\Omega$, i.e.,

$\lim\limits_{n\to\infty}W_{0}^{n}\left(x^{(0)},y^{(0)}\right)=\left(0,0\right).$

###### Proof.

Let

$f_{1}(x,y)=\frac{\beta
y^{2}}{\gamma+y}+x\left(1-d_{0}-\frac{\alpha}{1+x}\right),$
$f_{2}(x,y)=\frac{\alpha x}{1+x}+(1-\mu)y.$

Since functions $f_{1}(x,y)$ and $f_{2}(x,y)$ satisfy all the conditions of
Proposition 3, the trajectory of any initial point within the set $\Omega$
converges to zero. ∎

###  4.1. Dynamics on invariant sets.

A set $A$ is called invariant with respect to $W_{0}$ if $W_{0}(A)\subset A$.
We denote

$\Omega_{1}(x^{*},y^{*})=\left\\{(x,y)\in\Omega,\ \ 0\leq x\leq x^{*},\ \
0\leq y\leq y^{*}\\}\setminus\\{(x^{*},y^{*})\right\\},$
$\Omega_{2}(x^{*},y^{*})=\left\\{(x,y)\in\Omega,\ \ x^{*}\leq
x\leq\omega_{1},y^{*}\leq
y\leq\omega_{2}\right\\}\setminus\left\\{(x^{*},y^{*})\right\\},$

where $(x^{*},y^{*})$ is the positive fixed point defined by $(\ref{fix123})$.

######  Lemma 4.

The sets $\Omega_{1}$ and $\Omega_{2}$ are invariant with respect to $W_{0}.$

###### Proof.

1) Let $0\leq x\leq x^{*},0\leq y\leq y^{*}.$ Then

$x^{\prime}-x^{*}=\frac{\beta
y^{2}}{\gamma+y}+x(1-d_{0}-\frac{\alpha}{1+x})-x^{*}\leq\frac{\beta
y^{*2}}{\gamma+y^{*}}+x^{*}(1-d_{0}-\frac{\alpha}{1+x})-x^{*}$ $=\frac{\beta
y^{*2}}{\gamma+y^{*}}-d_{0}x^{*}-\frac{\alpha x^{*}}{1+x}=0,$
$y^{*}-y^{\prime}=y^{*}-(1-\mu)y-\frac{\alpha x}{1+x}\geq
y^{*}-(1-\mu)y^{*}-\frac{\alpha x^{*}}{1+x^{*}}=\mu y^{*}-\frac{\alpha
x^{*}}{1+x^{*}}=0.$

So, $(x^{\prime},y^{\prime})\in W_{0}(\Omega_{1})\subset\Omega_{1}$

2) Let $x^{*}\leq x\leq\omega_{1},y^{*}\leq y\leq\omega_{2}.$ Then

$x^{*}-x^{\prime}=x^{*}-\frac{\beta
y^{2}}{\gamma+y}-x(1-d_{0}-\frac{\alpha}{1+x})\leq
x^{*}-\frac{\beta{y^{*}}^{2}}{\gamma+y^{*}}-x^{*}(1-d_{0}-\frac{\alpha}{1+x})$
$=-\frac{\beta y^{*2}}{\gamma+y^{*}}+d_{0}x^{*}+\frac{\alpha x^{*}}{1+x}=0,$
$y^{\prime}-y^{*}=\frac{\alpha x}{1+x}+(1-\mu)y-y^{*}\geq\frac{\alpha
x^{*}}{1+x^{*}}+(1-\mu)y^{*}-y^{*}=0.$

Therefore, $(x^{\prime},y^{\prime})\in W_{0}(\Omega_{2})\subset\Omega_{2}.$ ∎

######  Lemma 5.

The operator $W_{0}$ defined in (2.3) is monotone with respect to the $NE$
ordering, i.e., is a cooperative map.

###### Proof.

By (2.3) we have

$\begin{array}[]{ll}x^{\prime}_{1}-x^{\prime}_{2}=\frac{\beta(y_{1}-y_{2})\left(\gamma(y_{1}+y_{2})+y_{1}y_{2}\right)}{(\gamma+y_{1})(\gamma+y_{2})}+(x_{1}-x_{2})\left(1-d_{0}-\frac{\alpha}{(1+x_{1})(1+x_{2})}\right),\\\
y^{\prime}_{1}-y^{\prime}_{2}=\frac{\alpha(x_{1}-x_{2})}{(1+x_{1})(1+x_{2})}+(1-\mu)(y_{1}-y_{2}).\end{array}$
(4.14)

Take $(x_{1},y_{1})\preceq_{NE}(x_{2},y_{2})\in\mathbb{R}^{2}_{+}$ then
$x_{1}\leq x_{2}$ and $y_{1}\leq y_{2}$. By the last equalities from (4.14)
one can see that $x^{\prime}_{1}\leq x^{\prime}_{2}$ and $y^{\prime}_{1}\leq
y^{\prime}_{2}$, i.e.,

$(x^{\prime}_{1},y^{\prime}_{1})=W_{0}(x_{1},y_{1})\preceq_{NE}(x^{\prime}_{2},y^{\prime}_{2})=W_{0}(x_{2},y_{2}).$

This completes the proof. ∎

The following theorem describes the trajectory of any point
$(x^{(0)},y^{(0)})$ in invariant sets $\Omega_{1}(x_{1}^{*},y_{1}^{*})$ and
$\Omega_{2}(x_{1}^{*},y_{1}^{*})$.

######  Theorem 3.

Assume that (2.2) holds and let $\beta=\nu$. For the operator $W_{0}$ and any
initial point $(x^{(0)},y^{(0)})$, the following hold

$\lim_{n\to\infty}x^{(n)}=\left\\{\begin{array}[]{ll}0,\ \ \ \ \ \ \mbox{if}\
\ (x^{(0)},y^{(0)})\in\Omega_{1}(x_{1}^{*},y_{1}^{*}),\\\\[5.69054pt]
x_{1}^{*},\ \ \mbox{if}\ \
(x^{(0)},y^{(0)})\in\Omega_{2}(x_{1}^{*},y_{1}^{*}),\end{array}\right.$
$\lim_{n\to\infty}y^{(n)}=\left\\{\begin{array}[]{ll}0,\ \ \ \ \ \ \mbox{if}\
\ (x^{(0)},y^{(0)})\in\Omega_{1}(x_{1}^{*},y_{1}^{*}),\\\\[5.69054pt]
y_{1}^{*},\ \ \ \ \ \mbox{if}\ \
(x^{(0)},y^{(0)})\in\Omega_{2}(x_{1}^{*},y_{1}^{*}),\end{array}\right.$

where $(x_{1}^{*},y_{1}^{*})$ is the fixed point defined by $(\ref{fix123})$.

###### Proof.

For $\beta=\nu$ (see Theorem 1), $z_{1,1}=(0,0)$ is attracting, while
$z_{1,2}=(x_{1}^{*},y_{1}^{*})$ is a non-hyperbolic. By applying Proposition 1
to the ordered interval
$\mathopen{[\\![}z_{1,1},z_{1,2}\mathclose{]\\!]}_{NE}$, we can conclude that
the trajectory of an arbitrary initial point in
$\Omega_{1}(x_{1}^{*},y_{1}^{*})$ converges to zero. Proposition 3 implies
that the trajectory of any initial point selected from set
$\Omega_{2}(x_{1}^{*},y_{1}^{*})$ converges to a limit point
$(x^{*}_{1},y^{*}_{1})$. See Fig. 1(a) for numerical values and the visual
illustration. ∎

The following theorem describes the trajectory of any point
$(x^{(0)},y^{(0)})$ in invariant sets $\Omega_{1}(x_{2}^{*},y_{2}^{*})$ and
$\Omega_{2}(x_{2}^{*},y_{2}^{*})$.

######  Theorem 4.

Assume that (2.2) holds and let $\beta>\nu$. For the operator $W_{0}$ and any
initial point $(x^{(0)},y^{(0)})$, the following hold

$\lim_{n\to\infty}x^{(n)}=\left\\{\begin{array}[]{ll}0,\ \ \ \ \ \ \mbox{if}\
\ (x^{(0)},y^{(0)})\in\Omega_{1}(x_{2}^{*},y_{2}^{*}),\\\\[5.69054pt]
x_{3}^{*},\ \ \mbox{if}\ \
(x^{(0)},y^{(0)})\in\Omega_{2}(x_{2}^{*},y_{2}^{*}),\end{array}\right.$
$\lim_{n\to\infty}y^{(n)}=\left\\{\begin{array}[]{ll}0,\ \ \ \ \ \ \mbox{if}\
\ (x^{(0)},y^{(0)})\in\Omega_{1}(x_{2}^{*},y_{2}^{*}),\\\\[5.69054pt]
y_{3}^{*},\ \ \ \ \ \mbox{if}\ \
(x^{(0)},y^{(0)})\in\Omega_{2}(x_{2}^{*},y_{2}^{*}),\end{array}\right.$

where $(x_{2}^{*},y_{2}^{*})$, $(x_{3}^{*},y_{3}^{*})$ are the fixed point
defined by $(\ref{fix123})$.

###### Proof.

Under the condition $\beta>\nu$ (see Theorem 1), the fixed points
$z_{2,1}=(0,0)$ and $z_{2,3}=(x_{3}^{*},y_{3}^{*})$ are attracting, while
$z_{2,2}=(x_{2}^{*},y_{2}^{*})$ is a repelling fixed point. Then, by applying
Proposition 1 to the ordered interval
$\mathopen{[\\![}z_{2,1},z_{2,2}\mathclose{]\\!]}_{NE}$, we can conclude that
the trajectory of an arbitrary initial point in
$\Omega_{1}(x_{2}^{*},y_{2}^{*})$ converges to zero.

Recall that $\omega_{1}$ and $\omega_{2}$ are defined by (4.1). We divide the
set $\Omega_{2}(x_{2}^{*},y_{2}^{*})$ into four parts:

$\Omega^{1}(x_{3}^{*},y_{3}^{*})=\left\\{(x,y)\in\Omega_{2}(x_{2}^{*},y_{2}^{*}),\
\ x_{3}^{*}\leq x\leq\omega_{1},\ \ y_{3}^{*}\leq y\leq\omega_{2}\right\\},$
$\Omega^{2}(x_{3}^{*},y_{3}^{*})=\left\\{(x,y)\in\Omega_{2}(x_{2}^{*},y_{2}^{*}),\
\ x_{2}^{*}\leq x\leq x_{3}^{*},\ \ y_{3}^{*}\leq y\leq\omega_{2}\right\\},$
$\Omega^{3}(x_{3}^{*},y_{3}^{*})=\left\\{(x,y)\in\Omega_{2}(x_{2}^{*},y_{2}^{*}),\
\ x_{2}^{*}\leq x\leq x_{3}^{*},\ \ y_{2}^{*}\leq x\leq y_{3}^{*}\right\\},$
$\Omega^{4}(x_{3}^{*},y_{3}^{*})=\left\\{(x,y)\in\Omega_{2}(x_{2}^{*},y_{2}^{*}),\
\ x_{3}^{*}\leq x\leq\omega_{1},\ \ y_{2}^{*}\leq x\leq y_{3}^{*}\right\\}.$

By applying Proposition 1 to the ordered interval
$\mathopen{[\\![}z_{2,2},z_{2,3}\mathclose{]\\!]}_{NE}$, we can conclude that
the trajectory of an arbitrary initial point in
$\Omega^{3}(x_{3}^{*},y_{3}^{*})$ converges to $(x^{*}_{3},y^{*}_{3})$.
Proposition 3 implies that the trajectory of any initial point in
$\Omega^{1}(x_{3}^{*},y_{3}^{*})$ converges to $(x^{*}_{3},y^{*}_{3})$.

Assume that $\left(x^{(0)},y^{(0)}\right)\in
int\left(\Omega^{2}(x_{3}^{*},y_{3}^{*})\right)$ or
$\left(x^{(0)},y^{(0)}\right)\in
int\left(\Omega^{4}(x_{3}^{*},y_{3}^{*})\right)$. Then, there exists
$\left(\overline{x^{(0)}},\overline{y^{(0)}}\right)\in
int\left(\Omega^{3}(x_{3}^{*},y_{3}^{*})\right)$ such that

$\left(\overline{x^{(0)}},\overline{y^{(0)}}\right)\preceq_{NE}\left(x^{(0)},y^{(0)}\right)$

and $\left(\widetilde{x^{(0)}},\widetilde{y^{(0)}}\right)\in
int\left(\Omega^{1}(x_{3}^{*},y_{3}^{*})\right)$ such that

$\left(x^{(0)},y^{(0)}\right)\preceq_{NE}\left(\widetilde{x^{(0)}},\widetilde{y^{(0)}}\right).$

By monotonicity of $W_{0}$, we have

$W_{0}^{n}\left(\overline{x^{(0)}},\overline{y^{(0)}}\right)\preceq_{NE}W_{0}^{n}\left(x^{(0)},y^{(0)}\right)\preceq_{NE}W_{0}^{n}\left(\widetilde{x^{(0)}},\widetilde{y^{(0)}}\right)$

which implies
$W_{0}^{n}\left(x^{(0)},y^{(0)}\right)\rightarrow(x^{*}_{3},y^{*}_{3})$ as
$n\rightarrow\infty$. This completes the proof. See Fig. 1(a) for numerical
values and the visual illustration. ∎

$\omega_{1}$$x$$x^{*}_{1}$$\omega_{2}$$y$$y^{*}_{1}$$0$$(a)$$\omega_{1}$$x$$x^{*}_{2}$$x^{*}_{3}$$\omega_{2}$$y$$0$$y^{*}_{2}$$y^{*}_{3}$$(b)$

Figure 1.(a) Visual illustration of Theorem 3 when $\beta=\nu$,
$\alpha=0.4,\beta=\frac{9}{4}(2+\sqrt{3}),\gamma=1,d_{0}=0.5,\mu=0.6$, in the
case when $(x_{1}^{*},y_{1}^{*})$ is a non-hyperbolic point. (b) Visual
illustration of Theorem 4 when $\beta>\nu$,
$\alpha=0.4,\beta=\frac{9}{4}(2+\sqrt{3}),\gamma=1,d_{0}=0.5,\mu=0.6$. The
case when $(x_{2}^{*},y_{2}^{*})$ is a repeller and $(0,0)$ and
$(x_{3}^{*},y_{3}^{*})$ are attracting points.

###  4.2. On the set
$\Omega\setminus\left(\Omega_{1}(x^{*},y^{*})\bigcup\Omega_{2}(x^{*},y^{*})\right)$

In the following examples, we show trajectories of initial points from the set
$\Omega\setminus\left(\Omega_{1}\bigcup\Omega_{2}\right)$.

Let us consider the operator $W_{0}$ with parameter values
$\alpha=0.4,\beta=\frac{9}{4}(2+\sqrt{3}),\gamma=1,d_{0}=0.5,\mu=0.6,$
satisfying the condition $\beta=\nu$. Then by (3.4), we get
$\left(x_{1}^{*},y_{1}^{*}\right)=\left(\frac{3\sqrt{3}}{5},9-5\sqrt{3}\right)$.

######  Example 1.

For initial points $z_{1}^{(0)}=(0.1,0.5)$ and $z_{2}^{(0)}=(2.5,0.1)$, the
trajectory of the system (2.3) is shown in Fig.2(a). In this case, the system
converges to the origin, i.e.,

$\lim_{n\to\infty}x^{(n)}=0,\ \lim_{n\to\infty}y^{(n)}=0.$

For initial points $z_{3}^{(0)}=(0.1,0.6)$ and $z_{4}^{(0)}=(3,0.1)$, the
trajectory of the system (2.3) is also shown in Fig.2(a) (refer to the same
figure for both cases). Here, the first coordinate of the trajectory converges
to $\frac{3\sqrt{3}}{5}$, while the second coordinate has a limit point of
$9-5\sqrt{3}$:

$\lim_{n\to\infty}x^{(n)}=\frac{3\sqrt{3}}{5},\
\lim_{n\to\infty}y^{(n)}=\frac{\alpha}{\mu}=9-5\sqrt{3}.$

We consider the operator $W_{0}$ with the following parameter values:
$\alpha=0.4,\beta=9,\gamma=1,d_{0}=0.5,$ and $\mu=0.6$. Notably, these
parameters satisfy the condition $\beta>\nu$. By (3.4), we obtain the fixed
points $\left(x_{2}^{*},y_{2}^{*}\right)=\left(\frac{3}{5},\frac{1}{4}\right)$
and $\left(x_{3}^{*},y_{3}^{*}\right)=\left(\frac{9}{5},\frac{3}{7}\right)$.

######  Example 2.

Figure 2(b) illustrates the trajectories of system (2.3) for initial points
$z_{1}^{(0)}=(0.1,0.3)$ and $z_{2}^{(0)}=(1,0.1)$. In both cases, the system
converges to the origin.

Considering initial points $z_{3}^{(0)}=(0.5,0.5)$, $z_{4}^{(0)}=(0.5,0.66)$,
$z_{5}^{(0)}=(2,0.1)$, and $z_{6}^{(0)}=(4,0.3)$, the trajectory of system
(2.3) is shown in Fig. 2(b). In this case, the first coordinate of the
trajectory appears to converge towards $\frac{9}{5}$, while the second
coordinate appears to have a limit point of $\frac{3}{7}$.

$\omega_{1}$$x$$x^{*}_{1}$$\omega_{2}$$y$$y^{*}_{1}$$0$$(a)$$z_{1}^{(0)}$$z_{2}^{(0)}$$z_{3}^{(0)}$$z_{4}^{(0)}$$\omega_{1}$$x$$x^{*}_{2}$$x^{*}_{3}$$\omega_{2}$$y$$0$$y^{*}_{2}$$y^{*}_{3}$$z_{1}^{(0)}$$z_{2}^{(0)}$$z_{3}^{(0)}$$z_{4}^{(0)}$$z_{5}^{(0)}$$z_{6}^{(0)}$$(b)$

Figure 2.(a) Visual illustration of Example 1 when $\beta=\nu$, and initial
points $z_{1}^{(0)}=(0.1,0.5)$, $z_{2}^{(0)}=(2.5,0.1)$,
$z_{3}^{(0)}=(0.1,0.6)$, $z_{4}^{(0)}=(3,0.1)$. (b) Visual illustration of
Example 2 when $\beta>\nu$, and initial points $z_{1}^{(0)}=(0.1,0.3)$,
$z_{2}^{(0)}=(1,0.1)$, $z_{3}^{(0)}=(0.5,0.5)$, $z_{4}^{(0)}=(0.5,0.66)$,
$z_{5}^{(0)}=(2,0.1)$, $z_{6}^{(0)}=(4,0.3)$.

##  5\. Discussion and conclusion

We note that for continuous time system (1.2) the following results are known
(see [14]): All fixed (equilibrium) points are found and their types are
determined. Moreover, local behavior of the dynamical system in the
neighborhood of the fixed point $(0,0)$ is studied. In our discrete-time case,
for $d_{0}>0$ we also determined types of all fixed points. Besides this we
have been able to study global behavior of the system in neighborhood of
$(0,0)$. Also, we have found invariant sets and studied the dynamical system
on the sets. The last results are not known for the continuous time.

## References

  * [1] Bilgin A., Brett A., Kulenovic M.R.S., Pilav E. (2018), Global Dynamics of a Cooperative Discrete System in the Plane, Inter. Jour. of Bif. and Chaos, Vol. 28, No. 7, DOI:10.1142/S0218127418300227
  * [2] Boxonov Z.S., Rozikov U.A. (2021), A discrete-time dynamical system of stage-structured wild and sterile mosquito population, Nonlinear studies, Vol. 28, No. 2, p.413-425.
  * [3] Boxonov Z.S., Rozikov U.A. (2021), Dynamical system of a mosquito population with distinct birth-death rates, Journal of Applied Nonlinear Dynamics, Vol. 10, No. 4, p.807-816. DOI:10.5890/JAND. 2021.12.015
  * [4] Boxonov Z.S., (2023), A discrete-time dynamical system of mosquito population, Journal of Difference Equations and Applications, Vol. 29, No. 1, p.67-83.doi.org/10.1080/10236198.2022.2160245
  * [5] Boxonov Z.S., (2024), Description of trajectories of an evolution operator generated by mosquito population. Russian Journal of Nonlinear Dynamics, Vol. 20, No. 2, pp. 197–207. DOI: 10.20537/nd240401
  * [6] Dancer E., Hess P., (1991) Stability of fixed points for order preserving discrete-time dynamical systems, J.Reine Angew. Math. 419, 125–139.
  * [7] Devaney R.L., (2003), An Introduction to Chaotic Dynamical System, Westview Press.
  * [8] Fauvergue X., (2012), A review of mate-finding Allee effects in insects: from individual behavior to population management, https://doi.org/10.1111/eea.12021
  * [9] Grove E.A., Ladas G., (2005), Periodicities in Nonlinear Difference Equations, CHapman Hall/CRC, 377 pp.
  * [10] Hess P., (1991) Periodic-Parabolic Boundary Value Problems and Positivity, Pitman Research Notes in Mathematics Series, Vol. 247 (Longman Scientific and Technical, Harlow; John Wiley and Sons, Inc., NY).
  * [11] Kulenovic, M.R.S., Merino O., (2009) Global bifurcation for competitive systems in the plane, Discr. Contin. Dyn. Syst. Ser. B 12, 133–149.
  * [12] Lyubich Y.I. (1992), Mathematical structures in population genetics, Springer-Verlag, Berlin.
  * [13] Li J., (2011), Malaria model with stage-structured mosquitoes, Math. Biol.Eng., 8, 753-768.
  * [14] Li J., Cai L., Li Y., (2017), Stage-structured wild and sterile mosquito population models and their dynamics, Journal of Biological Dynamics, 11(2) , 79-101.
  * [15] Lutambi A. M., Penny M.A., Smith T., Chitnis N. (2013), Mathematical modelling of mosquito dispersal in a heterogeneous environment, Mathematical Biosciences. 241(2) 198-216.
  * [16] Niu L., Wang S., (2024) Global stability and period-doubling bifurcations of a discrete Kolmogorov predator-prey model with Ricker-type prey growth, arXiv:2405.11806v1
  * [17] Rozikov U.A., (2020), Population dynamics: algebraic and probabilistic approach, World Sci. Publ. Singapore., 460 pp.
  * [18] Rozikov U.A., Boxonov Z.S., (2024) A discrete-time dynamical system of wild mosquito population with Allee effects, Discontinuity, Nonlinearity, and Complexity, Vol.13, No. 3, p. 495-506. DOI:10.5890/DNC.2024.09.008
  * [19] Rozikov U.A., Velasco M.V. (2019), A discrete-time dynamical system and an evolution algebra of mosquito population, Jour. Math. Biology. 78(4) 1225-1244.
  * [20] Wang Ch., Li X., (2014), Stability and Neimark-Sacker bifurcation of a semi-discrete population model, Journal of Applied Analysis and Computation, Vol.4, No.4, 419-435