Bifurcation Patterns and Chaos Control in
Discrete-Time Coral Reef Model
M. Priyankaa, P. Muthukumara111Corresponding Author; Email:
<EMAIL_ADDRESS>Phone: 91-451-2452371, Fax: 91-451-2454466.
<EMAIL_ADDRESS>(M. Priyanka).,
aDepartment of Mathematics,
The Gandhigram Rural Institute (Deemed to be University),
Gandhigram - 624 302, Tamil Nadu, India.

> Abstract: The reduction in coral reef densities, characterized by the
> proliferation of macroalgae, has emerged as a global threat. In this paper,
> we present a discrete-time coral reef dynamical model that incorporates
> macroalgae. We explore all ecologically possible equilibrium points for the
> proposed model. The conditions for the local stability of the interior
> equilibrium point are analyzed, which represents the coexistence of both
> coral and macroalgae. Furthermore, we investigate the model’s behavior using
> the center manifold theorem and bifurcation theory. Our analysis reveals
> that the model undergoes codimension-one bifurcations, specifically period-
> doubling and Neimark-Sacker bifurcations. To address the chaos resulting
> from the emergence of the Neimark-Sacker bifurcation, we apply the OGY
> feedback control method and a hybrid control methodology. Finally, we
> provide numerical simulations not only to validate the obtained results but
> also to demonstrate the complex dynamic behaviors that arise. These
> behaviors include reversal period-doubling bifurcation, period-4, 8, and 24
> bubble bifurcations, as well as chaotic behavior.

Keywords: Bubble bifurcation; Coral reef model; Discrete model; OGY control;
Period-doubling and Neimark-Sacker bifurcation; Stability.

Mathematics subject classification 2020: 37G10; 37G35; 86A08; 92D25.

##  1\. Introduction

The biodiversity of coral reefs plays a crucial role in maintaining the
environment’s health and holds significant economic value. Due to their
ecosystem functions, which include coastal storm defense and the preservation
of fisheries and marine biodiversity, coral reefs are vital to aquatic
environments. Additionally, coral reef-dependent fishing and tourism provide
livelihoods for millions of people [1]. However, reefs are experiencing a
severe decline; many of these populations are already considered seriously
damaged and remain at risk of disappearing soon [2]. Pristine reefs no longer
exist, and despite local efforts to safeguard coral reefs over the past 30
years, regional-scale declines persist.

Scientists believe that millions of unknown organisms thrive in and around
reefs. Future discoveries of new treatments depend on this biodiversity [3,
4]. Coral reef animals and plants are now used to create medicines that
potentially treat diseases such as cancer, arthritis, human bacterial
infections, viruses, and other conditions. Many studies have concentrated on
the effects of macroalgae on coral reefs and the potential role of these
interactions leading to coral demise [5].

The interaction between these species plays a significant role in the oxygen
levels in the atmosphere [6]. Consequently, fluctuations in the population of
these species can contribute to global warming, potentially resulting in
devastating consequences for life on Earth, even risking and leading to the
extinction of various animal species, including humans. Macroalgal forests are
highly productive, providing food and habitat for numerous species, including
fish [7, 8]. Hence, it is crucial to mathematically analyze the interaction
between coral and macroalgae for the continued existence of these species. The
ongoing study of coral reefs and macroalgae presents current and future
threats to humans, along with novel research ideas to support the management
of these essential natural resources.

There are two categories of mathematical models for population dynamics: those
defined by continuous-time differential equations and those defined by
discrete-time difference equations. Generations with overlapping growth
processes, such as the human population, are described by nonlinear
differential equations. In other cases, population growth occurs over discrete
time intervals with non-overlapping generations, as seen in 13-year periodical
cicadas. Nonlinear difference equations mathematically describe this situation
[9].

Numerous recent studies suggest that discrete-time equations describe marine
dynamics more accurately and naturally [10, 11]. Discrete-time models
significantly reduce computing time and retain the essential characteristics
of their related continuous-time counterparts [12, 13]. Moreover, they exhibit
richer dynamics than continuous-time models, leading to chaotic behaviors from
a biological perspective [14, 15]. Therefore, in recent years, many
researchers have focused on discrete-time population models. However, there
remains a research gap in the discrete-time analysis of coral reef models.

The primary purpose of modeling population dynamics is to identify the
controlling parameters and predict the expected outcomes when environmental
parameters change [16]. Model parameters determine the dynamical properties of
equilibriums in nonlinear dynamical systems. Changes in the qualitative
structure of the model corresponding to parameter fluctuations are known as
bifurcation points [17]. Therefore, it is crucial to identify the variables
that can cause instability and unpredictable behavior, such as chaos.

Moreover, detecting codimension-one bifurcations in a model allows one to
anticipate global phenomena, including hysteresis, invariant tori, limit
cycles, homoclinic bifurcations, and chaotic attractors [18]. For this reason,
many authors concentrate on bifurcation analysis and chaos control in
population dynamical models [19, 20].

Several researchers have studied chaotic discrete dynamical systems.
Researchers in [21] employs three chaos control methods to examine the
fractional-ordered discretized two-dimensional Leslie-Gower prey-predator
model. In [22], a non-standard finite difference scheme discretizes a
continuous-time Leslie prey-predator model. Additionally, authors in [13]
comprehensively investigate a discrete food web model enriched with mate-
finding Allee effects and hunting cooperation to explore resilience, chaos,
and bifurcations.

In this paper, we contribute to understanding a discrete-time coral reef
model. To date, there has been no comprehensive qualitative study of a
discrete-time coral reef model that considers the significance of growth rate
on macroalgae. Therefore, this paper focuses on the stability and bifurcation
analysis of a two-dimensional discrete-time coral reef model. This study
employs center manifold theorem [23, 24, 25] and bifurcation theory [26, 23,
24] to analyze the discrete-time coral reef model within the interior of
$\mathbb{R}^{2}_{+}$. We demonstrate the existence of flip and Neimark-Sacker
bifurcation in the proposed model. The primary contribution of this paper is
outlined as follows:

  * •

The proposed coral reef model comprises macroalgae, coral reefs, and algal
turfs. The conditions for the local stability of the model’s steady states are
derived.

  * •

The existence of flip and Neimark-Sacker bifurcation are discussed with the
help of the center manifold theorem and bifurcation theory.

  * •

OGY feedback control method is implemented to control the chaos resulting from
the emergence of Neimark-Sacker bifurcation. Then, the theoretical analysis is
validated using numerical simulations.

  * •

Similar to existing studies [1, 27, 28, 29, 30], we investigated the
interaction between coral reefs and macroalgae with discrete-time. Observing
the model reveals complex dynamic behaviors and demonstrates a delayed
response in bifurcation when subjected to an incremental increase in the
intrinsic growth rate of macroalgae.

Paper organization: Section 2 includes model formation and description.
Section 3 studies the existence and stability of the model’s steady states. In
section 4, we demonstrate that the model undergoes flip and Neimark-Sacker
bifurcation around the interior equilibrium point. In Section 5, the OGY
method and hybrid control procedures are employed to control the chaos due to
the appearance of the Neimark-Sacker bifurcation. Finally, results and
discussions were given in Section 6.

##  2\. Formulation of model with key assumptions

In this section, we start with a basic conceptual framework of the coral reef
model, including corals, macroalgae, and short algal turfs.

$\displaystyle\frac{dM}{dt}$
$\displaystyle=\underbrace{rM\left(1-\frac{M}{k}\right)}_{\textnormal{Logistic
growth}}+\underbrace{aMC}_{\begin{subarray}{c}\textnormal{Coral overgrown}\\\
\textnormal{by
macroalgae}\end{subarray}}-\underbrace{\frac{gM}{M+S}}_{\textnormal{Grazing
rate}}+\underbrace{\gamma MS}_{\begin{subarray}{c}\textnormal{Macroalgae}\\\
\textnormal{spread rate}\end{subarray}},$ $\displaystyle\frac{dC}{dt}$
$\displaystyle=\underbrace{\alpha SC}_{\begin{subarray}{c}\textnormal{Coral
overgrows}\\\ \textnormal{algal
turfs}\end{subarray}}-\underbrace{dC}_{\textnormal{Coral
mortality}}-\underbrace{aMC}_{\begin{subarray}{c}\textnormal{Coral
overgrown}\\\ \textnormal{by macroalgae}\end{subarray}},$
$\displaystyle\frac{dS}{dt}$
$\displaystyle=-\underbrace{rM\left(1-\frac{M}{k}\right)}_{\textnormal{Logistic
growth}}+\underbrace{\frac{gM}{M+S}}_{\textnormal{Grazing
rate}}-\underbrace{\gamma MS}_{\begin{subarray}{c}\textnormal{Macroalgae}\\\
\textnormal{spread rate}\end{subarray}}-\underbrace{\alpha
SC}_{\begin{subarray}{c}\textnormal{Coral overgrows}\\\ \textnormal{algal
turfs}\end{subarray}}+\underbrace{dC}_{\textnormal{Coral mortality}},$ (1)

where $M$, $C$, and $S$ denote the densities of macroalgae, corals, and algal
turfs, respectively. For constructing the mathematical model (2), the
following assumptions are made:

  * (i)

It is assumed that a specified region of the seabed is fully sheltered by
macroalgae $(M)$, coral $(C)$, and algal turfs $(S)$ so that $M+C+S$ is kept
constant over time, by a rescaling, lets assume that $M+C+S=1$. Hence, algal
turfs is defined by $S=1-M-C$ and accordingly $\frac{dS}{dt}$ is given by
$-\frac{dM}{dt}-\frac{dC}{dt}$.

  * (ii)

It is assumed that corals recruit and overgrow algal turfs at a rate $\alpha$,
they have a natural mortality rate of $d$, are overgrown by macroalgae at a
rate $a$, and macroalgae spread vegetatively over algal turfs at a rate
$\gamma$.

  * (iii)

Algal turfs are believed to be the space recolonized after coral mortality.
Additionally, the ratio $\frac{M}{M+S}$ represents the percentage of grazing
$g$ that affects macroalgae.

  * (iv)

Here, we incorporate logistic macroalgae growth into the proposed coral reef
model with an intrinsic growth rate $r$ and a time-varying carrying capacity
$k$.

  * (v)

Contrasting growth dynamics between macroalgae (logistic growth) and algal
turfs (negative logistic growth) reflect competitive interactions within the
ecosystem.

Biological description of the proposed model (2) parameters are given in Table
1. The pictorial representation of the coral reef model is shown in Figure 1.

Figure 1. The schematic representation of coral reef model. Table 1. Model (2) parameter description. Parameter |  Interpretation
---|---
$r$ |  intrinsic growth rate
$k$ |  carrying capacity of macroalgae
$a$ |  rate at which coral overgrows macroalgae
$g$ |  macroalgae grazing rate
$\gamma$ |  spread rate of macroalgae over algal turfs
$\alpha$ |  rate at which coral overgrows algal turfs
$d$ |  natural mortality rate of coral

Therefore, two equations are sufficient to describe the dynamics of this model
(see [1] for more details). From those mentioned above and setting $S=1-M-C$,
the coral dynamics are given by

$\displaystyle\frac{dM}{dt}$
$\displaystyle=rM\left(1-\frac{M}{k}\right)+aMC-\frac{gM}{1-C}+\gamma
M(1-M-C),$ $\displaystyle\frac{dC}{dt}$ $\displaystyle=\alpha(1-M-C)C-dC-aMC.$
(2)

To discretize the continuous-time coral reef model (2), we use the forward
Euler technique, for $n=0,1,2,\cdots,$

$\displaystyle M_{n+1}$
$\displaystyle=M_{n}+\delta\left[rM_{n}\left(1-\frac{M_{n}}{k}\right)+aM_{n}C_{n}-\frac{gM_{n}}{1-C_{n}}+\gamma
M_{n}(1-M_{n}-C_{n})\right],$ $\displaystyle C_{n+1}$
$\displaystyle=C_{n}+\delta\Big{[}\alpha(1-M_{n}-C_{n})C_{n}-dC_{n}-aM_{n}C_{n}\Big{]},$
(3)

such that $M_{0}=M(0)\geq 0$, $C_{0}=C(0)\geq 0,$ and $\delta>0$ is the step
size. We have discretized the model to emphasize the significance of step size
in chaos control analysis for the dynamic behavior of coral reef population
model. The corresponding map representation is then provided as:

$\displaystyle M$ $\displaystyle\longrightarrow
M+\delta\left[rM\left(1-\frac{M}{k}\right)+aMC-\frac{gM}{1-C}+\gamma
M(1-M-C)\right],$ $\displaystyle C$ $\displaystyle\longrightarrow
C+\delta\Big{[}\alpha(1-M-C)C-dC-aMC\Big{]}.$ (4)

##  3\. Stability analysis of two-component coral reef model

This section discusses the existence of fixed points and their stability
analysis for the proposed model (2). The equilibrium points of the proposed
coral reef model are listed below:

  * (i)

Trivial equilibrium $E_{0}(0,0).$

  * (ii)

Axial equilibrium $E_{A_{1}}(M_{A},0),$ where
$M_{A}$=$\frac{(r-g+\gamma)k}{(r+k\gamma)}$.

  * (iii)

Axial equilibrium $E_{A_{2}}(0,C_{A}),$ where $C_{A}=1-\frac{d}{\alpha}$.

  * (iv)

Interior equilibrium, $E^{*}(M^{*},C^{*})$ in the presence of both populations
is the intersection point of two following nontrivial nullclines in the first
quadrant of $\mathbb{R}^{2}$:

$\displaystyle r\left(1-\frac{M}{k}\right)+aC-\frac{g}{1-C}+\gamma(1-M-C)$
$\displaystyle=0,$ $\displaystyle\alpha(1-M-C)-d-aM$ $\displaystyle=0.$

Figure 2. Nullclines indicates the occurrence of two equilibrium points of
model (2). Figure 3. Nullclines indicates the occurrence of unique equilibrium
point of model (2).

Depending on the parameters, the number of feasible interior equilibrium
points for the proposed model is determined by comparing the relative
positions and shapes of the nontrivial nullclines. Numerical analysis suggests
that the potential number of coexistence equilibrium points ranges from zero
to two, as depicted in Figures 2-4, contingent upon the intrinsic growth rate
$r$. Consider the fixed set of parameters $k=0.5,$ $a=0.65,$ $g=0.3,$
$\gamma=0.4,$ $d=0.1,$ and $\alpha=0.6.$ For $r=0.2,$ the algae nullcline
(green curve) and coral nullcline (red curve) cross twice, indicating the
occurrence of two interior equilibrium (0.162603, 0.494576) and (0.369384,
0.063783) in the feasible region shown in Figure 2. For $r=0.112287,$ the
algae nullcline (green curve) and coral nullcline (red curve) touch each other
and give unique equilibrium (0.27, 0.26) is shown in Figure 3. Two nullclines
neither overlap nor contact each other for $r=0.05$, implying no equilibrium
point, as shown in Figure 4.

Figure 4. Nullclines indicates the occurrence of no equilibrium point of model
(2).

Considering the permanence of two species, here we concentrate only on the
interior equilibrium $E^{*}(M^{*},C^{*})$. The eigenvalues of the Jacobian
matrices corresponding to $E^{*}(M^{*},C^{*})$ are used to investigate the
local stability of the proposed model (2). Consider the following Jacobian
matrix for model (2) about $E^{*}(M^{*},C^{*})$:

$\displaystyle J=J(M,C)=\begin{bmatrix}a_{11}&a_{12}\\\
a_{21}&a_{22}\end{bmatrix},$ (5)

where $a_{11}=1+\delta r-\frac{2\delta rM}{k}+\delta aC-\frac{\delta
g}{1-C}+\delta\gamma-2\delta\gamma M-\delta\gamma C$, $a_{12}=\delta
aM-\frac{\delta gM}{(1-C)^{2}}-\delta\gamma M$, $a_{21}=-\delta\alpha C-\delta
aC$, and $a_{22}=1+\delta\alpha-\delta\alpha M-2\delta\alpha C-\delta d-\delta
aM.$

Let $F(\lambda)=\lambda^{2}+\mathcal{A}\lambda+\mathcal{B}$ where
$\mathcal{A}=-\textrm{tr}(J)$ and $\mathcal{B}=\det(J)$ be the characteristic
equation of associated Jacobian matrix $J$ evaluated at a fixed point
$E^{*}(M^{*},C^{*})$. Then, the following Lemmas 1 and 2 may be demonstrated
by the relationships between roots of $F(\lambda)=0$ and its coefficients.
This will help us to examine the stability of interior equilibrium
$E^{*}(M^{*},C^{*})$.

######  Lemma 1.

[31] Suppose that $F(1)>0$, $\lambda_{1}$ and $\lambda_{2}$ are two
characteristic roots of $F(\lambda)=0.$ Then

  * 1.

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

  * 2.

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

  * 3.

$|\lambda_{1}|>1$ and $|\lambda_{2}|>1$ if and only if $F(-1)>0$ and
$\mathcal{B}>1,$

  * 4.

$\lambda_{1}=-1$ and $|\lambda_{2}|\neq 1$ if and only if $F(-1)=0$ and
$\mathcal{A}$ is neither 0 nor 2,

  * 5.

$\lambda_{1}$ and $\lambda_{2}$ are complex and $|\lambda_{1}|=1$ and
$|\lambda_{2}|=1$ if and only if $\mathcal{A}^{2}-4\mathcal{B}<0$ and
$\mathcal{B}=1.$

######  Lemma 2.

[31] Suppose $\lambda_{1}$ and $\lambda_{2}$ be characteristics roots of
$F(\lambda)=0$, then the interior equilibrium $E^{*}(M^{*},C^{*})$ of the
model (2) is called

  * 1.

a sink if $|\lambda_{1}|<1$ and $|\lambda_{2}|<1$,

  * 2.

a source if $|\lambda_{1}|>1$ and $|\lambda_{2}|>1$,

  * 3.

a saddle if $|\lambda_{1}|>1$ and $|\lambda_{2}|<1$ (or $|\lambda_{1}|<1$ and
$|\lambda_{2}|>1$),

  * 4.

non-hyperbolic if either $|\lambda_{1}|=1$ or $|\lambda_{2}|=1.$

######  Remark 1.

To establish the necessary and sufficient conditions for the local stability
of the model (2) in the following Theorem 1, it is essential to observe the
requirements given in Lemmas 1 and 2. The characteristic equation of the
Jacobian matrix $J(M,C)$ is given by

$\displaystyle\lambda^{2}+\mathcal{A}(M,C)\lambda+\mathcal{B}(M,C)=0,$

where

$\displaystyle-\mathcal{A}(M,C)=a_{11}+a_{22}=2+\delta\mathcal{U},$

where $\mathcal{U}=u_{1}+u_{2}$ and

$\displaystyle\mathcal{B}(M,C)=a_{11}a_{22}-a_{12}a_{21}=1+\delta\mathcal{U}+\delta^{2}\mathcal{V},$

where $\mathcal{V}=u_{1}u_{2}-\frac{a_{12}a_{21}}{\delta^{2}}.$ Here,
$u_{1}=r-\frac{2rM}{k}+aC-\frac{g}{1-C}+\gamma-2\gamma M-\gamma C$ and
$u_{2}=\alpha-\alpha M-2\alpha C-d-aM$. Then

$\displaystyle F(\lambda)=$
$\displaystyle\lambda^{2}-\textrm{tr}(J)\lambda+\det(J)$ $\displaystyle=$
$\displaystyle\lambda^{2}-(2+\delta\mathcal{U})\lambda+(1+\delta\mathcal{U}+\delta^{2}\mathcal{V}).$

Here, $F(-1)=4+2\delta\mathcal{U}+\delta^{2}\mathcal{V}$ and
$F(1)=\delta^{2}\mathcal{V}>0$. It is clear that from the assumption of Lemma
1, $F(1)>0$ if $\mathcal{V}>0$. Assume $\mathcal{V}>0.$

######  Theorem 1.

Let $E^{*}(M^{*},C^{*})$ be the positive interior equilibrium and Lemma 1 and
2 holds for the model (2). Then $E^{*}(M^{*},C^{*})$ is said to be

  * 1.

a sink if one of the following conditions holds:

    * (i)

$\mathcal{U}=-2\sqrt{\mathcal{V}}$ and
$0<\delta<-\frac{\mathcal{U}}{\mathcal{V}}$,

    * (ii)

$\mathcal{U}<-2\sqrt{\mathcal{V}}$ and
$0<\delta<\frac{-\mathcal{U}-\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}}$.

  * 2.

a source if one of the following conditions holds:

    * (i)

$\mathcal{U}=-2\sqrt{\mathcal{V}}$ and
$\delta>-\frac{\mathcal{U}}{\mathcal{V}}$,

    * (ii)

$\mathcal{U}<-2\sqrt{\mathcal{V}}$ and
$\delta>\frac{-\mathcal{U}+\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}}$,

    * (iii)

$\mathcal{U}\geq 0.$

  * 3.

a saddle if the following conditions holds:

    * $\mathcal{U}<-2\sqrt{\mathcal{V}}$ and $\frac{-\mathcal{U}-\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}}<\delta<\frac{-\mathcal{U}+\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}}$.

  * 4.

a non-hyperbolic if one of the following conditions holds:

    * (i)

$\mathcal{U}<-2\sqrt{\mathcal{V}},$
$\delta=\frac{-\mathcal{U}\pm\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}}$,
and $\delta\neq-\frac{2}{\mathcal{U}},-\frac{4}{\mathcal{V}}$,

    * (ii)

$-2\sqrt{\mathcal{V}}<\mathcal{U}<0$ and
$\delta=-\frac{\mathcal{U}}{\mathcal{V}}.$

##  4\. Bifurcation analysis

This section proposed to derive codimension-one Neimark-Sacker bifurcation
around the interior equilibrium $E^{*}(M^{*},C^{*})$ for the model (2) where
$\delta$ act as a bifurcation parameter.

###  4.1. Flip bifurcation analysis around $E^{*}(M^{*},C^{*})$

Let

$\displaystyle\mathcal{F}_{1}=\left\\{\Big{(}r,k,a,g,\gamma,\alpha,d,S,\delta\Big{)}:\delta=\frac{-\mathcal{U}-\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}},\mathcal{U}<-2\sqrt{\mathcal{V}},r,k,a,g,\gamma,\alpha,d,S,\delta>0\right\\},$

or

$\displaystyle\mathcal{F}_{2}=\left\\{\Big{(}r,k,a,g,\gamma,\alpha,d,S,\delta\Big{)}:\delta=\frac{-\mathcal{U}+\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}},\mathcal{U}<-2\sqrt{\mathcal{V}},r,k,a,g,\gamma,\alpha,d,S,\delta>0\right\\}.$

Now we analyze, the flip bifurcation of $E^{*}(M^{*},C^{*})$ if the parameters
fluctuate in a restricted area around $\mathcal{F}_{1}$ (or
$\mathcal{F}_{2}$).

We begin with a discussion about the flip bifurcation of model (2) at
$E^{*}(M^{*},C^{*})$ when parameters fluctuate in a restricted area around
$\mathcal{F}_{1}$. The other instance $\mathcal{F}_{2}$ can be justified by
the same reasoning. Consider the model (2) with
$(r,k,a,g,\gamma,\alpha,d,S,\delta_{1})\in\mathcal{F}_{1}$, which is given by

$\displaystyle M$ $\displaystyle\longrightarrow
M+\delta_{1}\left[rM\left(1-\frac{M}{k}\right)+aMC-\frac{gM}{M+S}+\gamma
MS\right],$ $\displaystyle C$ $\displaystyle\longrightarrow
C+\delta_{1}\Big{[}\alpha SC-dC-aMC\Big{]}.$ (6)

The interior equilibrium point $E^{*}(M^{*},C^{*})$ of map (4.1), whose
eigenvalues are $\lambda_{1}=-1$, and $\lambda_{2}=3+\mathcal{U}\delta_{1}$
with $|\lambda_{2}|\neq 1$ by Theorem 1.

Since $(r,k,a,g,\gamma,\alpha,d,S,\delta_{1})\in\mathcal{F}_{1}$ and
$\delta_{1}=\frac{-\mathcal{U}-\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}}$,
we examine a perturbation of (4.1) using $\delta_{1}^{*}$ as the bifurcation
parameter as follows:

$\displaystyle M$ $\displaystyle\longrightarrow
M+(\delta_{1}+\delta_{1}^{*})\left[rM\left(1-\frac{M}{k}\right)+aMC-\frac{gM}{M+S}+\gamma
MS\right],$ $\displaystyle C$ $\displaystyle\longrightarrow
C+(\delta_{1}+\delta_{1}^{*})\Big{[}\alpha SC-dC-aMC\Big{]},$ (7)

where $|\delta_{1}^{*}|\ll 1,$ is a small perturbation parameter around
$\delta_{1}$.

######  Theorem 2.

The interior equilibrium point $E^{*}(M^{*},C^{*})$ of model (2) experience
flip bifurcation when $\delta$ changes in a small neighborhood of $\delta_{1}$
if $\Omega_{1}\neq 0$ and $\Omega_{2}\neq 0$ which is defined in Appendix A.
Furthermore, period-2 orbits bifurcates from $E^{*}(M^{*},C^{*})$ is
stable(unstable) if $\Omega_{2}>0(<0)$.

###### Proof.

The Proof of Theorem 2 is given in Appendix A. ∎

######  Remark 2.

In discrete-time models, taking integral step size
$\delta=\delta_{2}=\frac{-\mathcal{U}-\sqrt{\mathcal{U}^{2}-4\mathcal{V}}}{\mathcal{V}}$
(which depends on model’s parameter, see the statement of Theorem 1) as a
bifurcation point can affect the stability and behavior of the model. When the
integral step size is small, the model exhibits stable behavior. However, the
model may exhibit complex behaviors such as chaos as the integral step size
increases. The value of $\delta$ can significantly impact the model’s
behavior, especially near bifurcation points. Changing $\delta$ can alter the
relative strengths of different terms in the model equations, potentially
triggering a qualitative shift in the system’s dynamics (bifurcation).
Therefore, the integral step size is an important parameter to consider when
analyzing the behavior of the discrete-time model. For this reason, the
authors[32, 33, 34, 35, 36] analyzed the behavior of the dynamical model by
considering integral step size, which depends on model parameters as a
bifurcation parameter.

###  4.2. Neimark-Sacker bifurcation analysis around $E^{*}(M^{*},C^{*})$

Consider,

$\displaystyle\mathcal{N}=\Bigg{\\{}\Big{(}r,k,a,g,\gamma,\alpha,d,\delta\Big{)}:\delta=-\frac{\mathcal{U}}{\mathcal{V}},-2\sqrt{\mathcal{V}}<\mathcal{U}<0,r,k,a,g,\gamma,\alpha,d,\delta>0\Bigg{\\}}.$

When parameters fluctuate in a small neighborhood of $\mathcal{N}$, Neimark-
Sacker bifurcation may occur at $E^{*}(M^{*},C^{*})$.

The Neimark-Sacker bifurcation of $E^{*}(M^{*},C^{*})$ is discussed if the
parameters $(r,k,a,g,\gamma,\alpha,d,\delta)$ vary within a small neighborhood
of $\mathcal{N}$. We take parameters $(r,k,a,g,\gamma,\alpha,d,\delta_{2})$
arbitrarily from $\mathcal{N}$ then the model is given by

$\displaystyle x$ $\displaystyle\longrightarrow
x+\delta_{2}\left[rM\left(1-\frac{M}{k}\right)+aMC-\frac{gM}{1-C}+\gamma
M(1-M-C)\right],$ $\displaystyle y$ $\displaystyle\longrightarrow
y+\delta_{2}\Big{[}r(1-M-C)C-dC-aMC\Big{]}.$ (8)

Given that
$(r,k,a,g,\gamma,\alpha,d,\delta_{2})\in\mathcal{N},\delta_{2}=-\frac{\mathcal{U}}{\mathcal{V}}$.
Taking $\delta_{2}^{*}$ as the bifurcation parameter, we take the following
perturbation of (4.2):

$\displaystyle x$ $\displaystyle\longrightarrow
x+(\delta_{2}+\delta_{2}^{*})\left[rM\left(1-\frac{M}{k}\right)+aMC-\frac{gM}{1-C}+\gamma
M(1-M-C)\right],$ $\displaystyle y$ $\displaystyle\longrightarrow
y+(\delta_{2}+\delta_{2}^{*})\Big{[}r(1-M-C)C-dC-aMC\Big{]},$ (9)

where $|\delta_{2}^{*}|\ll 1,$ is a small perturbation parameter.

######  Theorem 3.

The interior equilibrium $E^{*}(M^{*},C^{*})$ of the model (2) experiences
Neimark-Sacker bifurcation when $\delta$ changes in a small neighborhood of
$\delta_{2}$ if $\Psi\neq 0$ which is defined in Appendix B. Furthermore, an
attracting(repelling) invariant closed curve bifurcates from
$E^{*}(M^{*},C^{*})$ for $\delta>\delta_{2}(<\delta_{2})$ if $\Psi<0(>0)$.

###### Proof.

The Proof of Theorem 3 is given in Appendix B. ∎

##  5\. Chaos control

###  5.1. OGY chaos control method

Controlling chaos and bifurcation are believed to be crucial characteristics
of interacting populations. Discrete-time models often exhibit more complex
behavior than continuous ones. Chaos management techniques are required to
protect the population from unpredictable events. In this part, we apply the
OGY approach [37] to the model for guiding an unstable trajectory toward a
stable one. Consider model (2) by taking the integral step size $\delta$ as a
control parameter and applying the OGY approach to find the stable region
around the neighborhood of $\delta$

$\displaystyle M_{n+1}$
$\displaystyle=M_{n}+\delta\left[rM_{n}\left(1-\frac{M_{n}}{k}\right)+aM_{n}C_{n}-\frac{gM_{n}}{1-C_{n}}+\gamma
M_{n}(1-M_{n}-C_{n})\right]=f_{11}(M_{n},C_{n},\delta),$ $\displaystyle
C_{n+1}$
$\displaystyle=C_{n}+\delta\Big{[}r(1-M_{n}-C_{n})C_{n}-dC_{n}-aM_{n}C_{n}\Big{]}=f_{12}(M_{n},C_{n},\delta).$
(10)

Here, $\delta$ acts as a control parameter; very small perturbations in
$\delta$ can accomplish the desired chaos control. To do this, we constrain
$\delta$ to lie in a small interval
$\delta\in(\delta_{0}-\epsilon,\delta_{0}+\epsilon)$, $\epsilon>0$ where
$\delta_{0}$ implies that the nominal value belongs to the chaotic region.

We use the stabilizing feedback control approach to direct the trajectory
toward the intended orbit. The following linear map can be used to approximate
model (5.1) in the neighborhood of the unstable equilibrium
$E^{*}(M^{*},C^{*})$ of the model (2) in the chaotic zone created by the
formation of Neimark-Sacker bifurcation.

$\displaystyle\begin{bmatrix}M_{n+1}-M^{*}\\\
C_{n+1}-C^{*}\end{bmatrix}\approx
J(M^{*},C^{*},\delta_{0})\begin{bmatrix}M_{n}-M^{*}\\\
C_{n}-C^{*}\end{bmatrix}+B[\delta-\delta_{0}],$ (11)

where

$\displaystyle J(M^{*},C^{*},\delta_{0})$
$\displaystyle=\begin{bmatrix}\frac{\partial f_{11}}{\partial
M}&\frac{\partial f_{11}}{\partial C}\\\ \frac{\partial f_{12}}{\partial
M}&\frac{\partial f_{12}}{\partial
C}\end{bmatrix}=\begin{bmatrix}\hat{a}&\hat{b}\\\
\hat{c}&\hat{d}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}1\\!+\\!\delta
r\\!-\\!\frac{2\delta rM}{k}\\!+\\!\delta aC\\!-\\!\frac{\delta
g}{1-C}\\!+\\!\delta\gamma\\!-\\!2\delta\gamma M\\!-\\!\delta\gamma C&\delta
aM\\!-\\!\frac{\delta gM}{(1-C)^{2}}\\!-\\!\delta\gamma M\\\ -\delta\alpha
C\\!-\\!\delta aC&1\\!+\\!\delta\alpha\\!-\\!\delta\alpha
M\\!-\\!2\delta\alpha C\\!-\\!\delta d\\!-\\!\delta aM\end{bmatrix}$

and

$\displaystyle B=\begin{bmatrix}\frac{\partial f_{11}}{\partial\delta}\\\
\frac{\partial
f_{12}}{\partial\delta}\end{bmatrix}=\begin{bmatrix}rM\left(1-\frac{M}{k}\right)+aMC-\frac{gM}{1-C}+\gamma
M(1-M-C)\\\ r(1-M-C)C-dC-aMC\end{bmatrix}=\begin{bmatrix}\tilde{a}\\\
\tilde{b}\end{bmatrix}.$

Additionally, model (5.1) is controllable if $C$ is a rank 2 matrix,

$\displaystyle
C=[B:JB]=\begin{bmatrix}\tilde{a}&\hat{a}\tilde{a}+\hat{b}\tilde{b}\\\
\tilde{b}&\hat{c}\tilde{a}+\hat{d}\tilde{b}\end{bmatrix}.$

Take

$\displaystyle[\delta-\delta_{0}]=-H\begin{bmatrix}M_{n}-M^{*}\\\
C_{n}-C^{*}\end{bmatrix},$

where $H=[\rho_{1}~{}\rho_{2}]$ is a gain matrix. Then the model (11) is given
by

$\displaystyle\begin{bmatrix}M_{n+1}-M^{*}\\\
C_{n+1}-C^{*}\end{bmatrix}\approx[J-BH]\begin{bmatrix}M_{n}-M^{*}\\\
C_{n}-C^{*}\end{bmatrix}.$

Furthermore, $E^{*}(M^{*},C^{*})$ is locally asymptotically stable if and only
if all its eigenvalues lie in the unit disk. Consider $J-BH$

$\displaystyle J-BH=\begin{bmatrix}\hat{a}&\hat{b}\\\
\hat{c}&\hat{d}\end{bmatrix}-\begin{bmatrix}\tilde{a}\rho_{1}&\tilde{a}\rho_{2}\\\
\tilde{b}\rho_{1}&\tilde{b}\rho_{2}\end{bmatrix}=\begin{bmatrix}\hat{a}-\tilde{a}\rho_{1}&\hat{b}-\tilde{a}\rho_{2}\\\
\hat{c}-\tilde{b}\rho_{1}&\hat{d}-\tilde{b}\rho_{2}\end{bmatrix}.$

The characteristic equation $J-BH$ is given by

$\displaystyle
P(\lambda)=\lambda^{2}-(\lambda_{1}+\lambda_{2})\lambda+(\lambda_{1}\lambda_{2}),$

where

$\displaystyle\lambda_{1}+\lambda_{2}$
$\displaystyle=\hat{a}+\hat{d}-\tilde{a}\rho_{1}-\tilde{b}\rho_{2},$ (12)
$\displaystyle\lambda_{1}\lambda_{2}$
$\displaystyle=(\hat{a}-\tilde{a}\rho_{1})(\hat{d}-\tilde{b}\rho_{2})-(\hat{c}-\tilde{b}\rho_{1})(\hat{b}-\tilde{a}\rho_{2}).$
(13)

Next, to determine the lines of marginal stability for the associated
controlled model, we consider $\lambda_{1}=\pm 1$ and
$\lambda_{1}\lambda_{2}=1$. Additionally, these limitations guarantee that the
open unit disk contains $\lambda_{1}$ and $\lambda_{2}$. Assume
$\lambda_{1}\lambda_{2}=1$, then (13) became

$\displaystyle\mathcal{L}_{1}:(\hat{a}-\tilde{a}\rho_{1})(\hat{d}-\tilde{b}\rho_{2})-(\hat{c}-\tilde{b}\rho_{1})(\hat{b}-\tilde{a}\rho_{2})-1=0.$

Next, assume that $\lambda_{1}=1$, then from (12) and (13) we have

$\displaystyle\mathcal{L}_{2}:\hat{a}+\hat{d}-\tilde{a}\rho_{1}-\tilde{b}\rho_{2}-1=\hat{a}\hat{d}-\hat{a}\tilde{b}\rho_{2}-\tilde{a}\hat{d}\rho_{1}+\tilde{a}\tilde{b}\rho_{1}\rho_{2}-\hat{c}\hat{b}+\tilde{a}\hat{c}\rho_{2}+\tilde{d}\hat{b}\rho_{1}-\tilde{a}\tilde{b}\rho_{1}\rho_{2}.$

Finally assume $\lambda_{2}=-1$, then from (12) and (13) we have

$\displaystyle\mathcal{L}_{3}:\hat{a}+\hat{d}-\tilde{a}\rho_{1}-\tilde{b}\rho_{2}+1=-[\hat{a}\hat{d}-\hat{a}\tilde{b}\rho_{2}-\tilde{a}\hat{d}\rho_{1}+\tilde{a}\tilde{b}\rho_{1}\rho_{2}-\hat{c}\hat{b}+\tilde{a}\hat{c}\rho_{2}+\tilde{d}\hat{b}\rho_{1}-\tilde{a}\tilde{b}\rho_{1}\rho_{2}].$

Thus, the bounded region obtained by the intersection of lines
$\mathcal{L}_{1}$, $\mathcal{L}_{2}$, and $\mathcal{L}_{3}$ gives the stable
eigenvalues for the model (2).

###  5.2. Feedback control method

We use the hybrid control feedback methodology [38], which was mainly
developed for controlling the period-doubling bifurcation; a similar technique
is used in [39] to manage the chaos caused by the appearance of Neimark-Sacker
bifurcation. Assuming model (2) experiences Neimark-Sacker bifurcation at
equilibrium $E^{*}(M^{*},C^{*})$, the associated controlled model is given by

$\displaystyle M_{n+1}=\zeta f_{11}(M_{n},C_{n},\delta)+(1-\zeta)M_{n},$
$\displaystyle C_{n+1}=\zeta f_{12}(M_{n},C_{n},\delta)+(1-\zeta)C_{n},$ (14)

where $0<\zeta<1.$ The controlled method in (5.2) combines feedback control
$\zeta$ and parameter perturbation $\delta$. Additionally, the Neimark-Sacker
bifurcation of the equilibrium $E^{*}(M^{*},C^{*})$ of the controlled model
(5.2) can be delayed, or perhaps completely removed by making an appropriate
choice of the controlled parameter. From the Jacobian matrix (5), we write the
corresponding Jacobian $J^{*}$ for the controlled model (5.2)

$\displaystyle J^{*}=\begin{bmatrix}\zeta a_{11}+(1-\zeta)&\zeta a_{12}\\\
\zeta a_{21}&\zeta a_{22}+(1-\zeta)\end{bmatrix}.$

Then the characteristic equation is given by

$\displaystyle\lambda^{2}-\textrm{tr}(J^{*})\lambda+\det(J^{*})=0,$

where

$\displaystyle\textrm{tr}(J^{*})$ $\displaystyle=\zeta a_{11}+\zeta
a_{22}+2(1-\zeta),$ $\displaystyle\det(J^{*})$ $\displaystyle=\Big{(}\zeta
a_{11}+(1-\zeta)\Big{)}\Big{(}\zeta
a_{22}+(1-\zeta)\Big{)}-\zeta^{2}a_{12}a_{21}$

where $a_{11},$ $a_{12},$ $a_{21},$ and $a_{22}$ are given previously. Using
[39], the following result identifies the requisites for the positive
equilibrium $E^{*}(M^{*},C^{*})$ of the controlled model (5.2) to be locally
asymptotically stable.

######  Theorem 4.

If the following condition holds,

$\displaystyle|\textrm{tr}(J^{*})|-1<\textrm{det}(J^{*})<1,$

then the positive equilibrium $E^{*}(M^{*},C^{*})$ of the controlled model
(5.2) is locally asymptotically stable.

##  6\. Results and discussions

In this section, to verify the aforementioned theoretical analysis numerically
by exploiting complicated dynamical behaviors such as bifurcation diagrams and
phase portraits for the model (2). Particularly the Macroalgae and coral reef
populations dynamical behaviors are analyzed in marine ecosystems under
intrinsic growth rate.

Figure 5. Reverse period-doubling bifurcation for $\delta$ versus $M$ and $C$.

Consider the following set of parametric value for the proposed model (2),
$k=0.3,$ $a=0.9,$ $g=0.5,$ $\gamma=0.8,$ $d=0.1,$ $\alpha=0.5$ with $\delta$
varies in range $[2.2,2.75]$ and take initial condition as
$(M_{0},C_{0})=(0.04,0.66)$. The model (2) endures full Feigenbaum remerging
tree for $r=0.76$. The term “full Feigenbaum three” refers to the situation
where the sequence of period-doubling leading to chaos is followed by the
reverse process when a parameter varies monotonously. Reverse period-doubling
bifurcation diagrams for $M$ and $C$ for the proposed macroalgae-coral reef
model with macroalgae intrinsic growth rate $r=0.76$ with other parameter
fixed given previously is given in Figure 5.

Figure 6. Period-8 bubble bifurcation for $\delta$ versus $M$ and $C$.

Figure 7. Period-4 bubble bifurcation for $\delta$ versus $M$ and $C$.

If we change the intrinsic growth rate as $r=0.77$ with same fixed set of
parameters and $\delta$ varies in the range $[2.2,2.8]$, the model (2) creates
a primary period-8 bubble which is given in Figure 6. The bifurcation diagram
for the model with the same fixed set of parameters is shown with two
different initial conditions: (0.04, 0.66) plotted in violet and (0.035, 0.59)
plotted in green. As depicted in Figure 5 and 6, it is evident that the model
exhibits multistability. Period-4 bubble is generated for $r=0.78$ with
$\delta$ varies in the range $=[2.2,2.8]$ which is given in Figure 7.
Similarly, a period-24 bubble is developed for $r=0.8$ which is given in
Figure 10. Moreover interesting and exciting bifurcation exist for $r=0.75$
and $r=0.78$ which consists of several reversal period-doubling bifurcation
which is given Figures 8 and 9. When further increasing the values to $r=1$,
the model undergoes period-doubling bifurcation which is given in Figure 11.

######  Remark 3.

The growth rate of macroalgae $r$ has a significant impact on the dynamics of
the model even if we merely change the intrinsic growth rate and leave the
other parameters fixed.

We provided the phase portrait for the above set of parameters, with $\delta$
fixed at $2.8$, and $r$ varied as follows: (a) $r=0.795$, (b) $r=0.78$, (c)
$r=0.77$, and (d) $r=0.765$ (see Figure 12). The figures depict chaotic
attractors with varying values of $r$. Consider the following set of
parameters: $r=0.5$, $k=0.7$, $a=0.65$, $g=0.3$, $\gamma=0.4$, $d=0.1$,
$\alpha=0.6$, and $\delta=3.3$, with different initial conditions: (e) (0.005,
0.6), (f) (0.004, 0.6), and (g) (0.001, 0.6) (see Figure 12). The figures
depict coexisting attracting orbits.

Figure 8. Bifurcation diagram for $r=0.795$.

Figure 9. Bifurcation diagram for $r=0.798$.

######  Remark 4.

In summary, the observation from Figure 5 to Figure 10 underscores the
significant finding that the bifurcation point of the model exhibits a delayed
response when subjected to incremental increases in the intrinsic growth rate,
while all other pertinent parameters remain constant.

######  Case 1.

Consider $r=0.5,$ $k=0.7,$ $a=0.65,$ $g=0.3,$ $\gamma=0.4,$ $d=0.1,$ and
$\alpha=0.6,$ with $\delta$ varies in the range $[2.8,3.32].$

For Case 1 with the given parameters, the model experience Niemark-Sacker
bifurcation with initial condition $(M_{0},C_{0})=(0.06,0.7)$, $\delta$ varies
in the range $[2.8,3.32]$. The bifurcation diagram corresponds to Case 1
parameters in the $(\delta,M)$ and $(\delta,C)$ plane are given in Figure 13.
We fix the value of $\delta=3.1$, then the corresponding phase portrait is
given in Figure 14.

######  Case 2.

Consider $r=0.6,$ $k=0.9,$ $a=0.45,$ $g=0.3,$ $\gamma=0.4,$ and $d=0.1,$ with
$\delta$ varies in the range $[2.5,3.15].$

Figure 10. Period-24 bubble bifurcation for $\delta$ versus $M$ and $C$.

Figure 11. Bifurcation diagram showing the emergence of period-doubling
bifurcation for $r=1.$

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 12. Phase portraits for the above set of parameters are depicted with
fixed $\delta=2.8$ and varying $r$ values: (a) $r=0.795$, (b) $r=0.78$, (c)
$r=0.77$, and (d) $r=0.765$. Phase portraits with initial conditions: (e)
(0.005, 0.6), (f) (0.004, 0.6), and (g) (0.001, 0.6).

For Case 2 with the given parameters, the model experience Niemark-Sacker
bifurcation with initial condition $(M_{0},C_{0})=(0.08,0.68)$, $\delta$
varies in the range $[2.5,3.15]$. The bifurcation diagram corresponds to Case
2 parameters in the $(\delta,M)$ and $(\delta,C)$ plane are given in Figure
15. To illustrate the bifurcation nature of the model (2) with the set of
parameters in Case 2, we fix the value of $\delta$ in the given range such as
$\delta$ =2.95, 3, 3.05, and 3.1. Then the corresponding phase portraits are
given in Figures 16 and 17.

######  Remark 5.

In this remark, we verify the nature of period-2 orbits bifurcates from
$(0.04,0.66).$ For this set of parameters $r=1,$ $k=0.3,$ $a=0.9,$ $g=0.5,$
$\gamma=0.8,$ $d=0.1,$ $\alpha=0.5$, we have eigenvalues $\lambda_{1}=-1$ and
$\lambda_{2}=-12.3527$ for an interior equilibrium point equilibrium point
$(0.04,0.66)$ with $a_{1}=0.0134,$ $a_{2}=0.0666,$ $a_{3}=0,$ $b_{1}=-1.1347,$
$b_{2}=-4.7558,$ and $b_{5}=-0.1134.$ From $b_{1},$ $b_{2},$ and $b_{5}$ we
have two discriminatory values $\Omega_{1}=-9.5115$ and $\Omega_{2}=2.3481$
Clearly, both the values $\Omega_{1}$ and $\Omega_{2}$ are non-zero real
numbers. Hence, from Theorem 2 the model experience flip bifurcation.
Moreover, $\Omega_{2}>0$ which implies that period-2 orbits bifurcates from
the $E^{*}(M^{*},C^{*})$ is stable.

######  Case 3.

Consider $r=0.5,$ $k=0.7,$ $a=0.65,$ $g=0.3,$ $\gamma=0.4,$ $d=0.1,$ and
$\alpha=0.6,$ with $\delta$ varies in the range $[2.8,3.32].$

For Case 1 with the given parameters, the model experiences Niemark-Sacker
bifurcation with initial condition $(M_{0},C_{0})=(0.06,0.7)$, $\delta$ varies
in the range $[2.8,3.32]$. The bifurcation diagram corresponds to Case 1
parameters in the $(\delta,M)$ and $(\delta,C)$ plane are given in Figure 13.
We fix the value of $\delta=3.1$; the corresponding phase portrait is shown in
Figure 14.

######  Case 4.

Consider $r=0.6,$ $k=0.9,$ $a=0.45,$ $g=0.3,$ $\gamma=0.4,$ $d=0.1,$ and
$\alpha=0.6$ with $\delta$ varies in the range $[2.5,3.15].$

For Case 2 with the given parameters, the model experiences Niemark-Sacker
bifurcation with initial condition $(M_{0},C_{0})=(0.08,0.68)$, $\delta$
varies in the range $[2.5,3.15]$. The bifurcation diagram corresponds to Case
2 parameters in the $(\delta,M)$ and $(\delta,C)$ plane are given in Figure
15. To illustrate the bifurcation nature of the model (2) with the set of
parameters in Case 2, we fix the value of $\delta$ in the given range such as
$\delta$ =2.95, 3, 3.05, and 3.1. Then, the corresponding phase portraits are
shown in Figures 16 and 17.

Figure 13. Bifurcation diagram showing the emergence of Neimark-Sacker
bifurcation for Case 1. Figure 14. Phase portrait of Case 1 parameters for
$\delta=3.1.$

Figure 15. Bifurcation diagram showing the emergence of Neimark-Sacker
bifurcation for Case 2.

######  Remark 6.

In this remark, we verify the nature of Neimark-Sacker bifurcation emergence
from the fixed point $(0.06,0.7).$ For the Case 1 parameter, we have
eigenvalue

$\displaystyle\lambda,\bar{\lambda}=1+\frac{\mathcal{U}(\delta_{2}+\delta_{2}^{*})}{2}\pm\frac{i(\delta_{2}+\delta_{2}^{*})}{2}\sqrt{4\mathcal{V}-\mathcal{U}^{2}},$

$\lambda=3.3318$ and $\bar{\lambda}=0.3002$ for an interior equilibrium with
$\varsigma_{1}=4.2199+31.0638i,$ $\varsigma_{2}=-3.7491-16.7558i,$
$\varsigma_{3}=1.4140+3.2607i,$ and $\varsigma_{4}=-8.5369$. From $\lambda,$
$\bar{\lambda},$ $\varsigma_{1},$ $\varsigma_{2},$ $\varsigma_{3},$ and
$\varsigma_{4}$ we have the discriminatory value $\Psi=-125.1820\neq 0.$
Hence, from Theorem 3 the model experience Neimark-Sacker bifurcation.
Moreover, $\Psi<0$ it indicates that repelling invariant closed curve
bifurcates from $E^{*}(M^{*},C^{*}).$

Figure 16. Phase portraits of Case 2 parameters for $\delta=2.95$ and
$\delta=3$, respectively.

Figure 17. Phase portraits of Case 2 parameters for $\delta=3.05$ and
$\delta=3.1$, respectively.

######  Remark 7.

In this remark, we verify the nature of Neimark-Sacker bifurcation emergence
from the fixed point $(0.08,0.68).$ For the Case 2 parameter, we have
eigenvalue $\lambda=3.3710$ and $\bar{\lambda}=0.2966$ for an interior
equilibrium with $\varsigma_{1}=4.2725+26.3698i,$
$\varsigma_{2}=-3.8442-13.3660i,$ $\varsigma_{3}=1.5006+2.4943i,$ and
$\varsigma_{4}=-7.0960$. From $\lambda,$ $\bar{\lambda},$ $\varsigma_{1},$
$\varsigma_{2},$ $\varsigma_{3},$ and $\varsigma_{4}$ we have the
discriminatory value $\Psi=-83.5551\neq 0.$ Hence, from Theorem 3 the model
experience Neimark-Sacker bifurcation. Moreover, $\Psi<0$ it indicates that
repelling invariant closed curve bifurcates from $E^{*}(M^{*},C^{*}).$

It is established that the model experiences Neimark-Sacker bifurcation for
the given set of parameters in Case 2. We obtain the following controlled
model that corresponds to the Case 2 parametric values to apply the OGY
method:

$\displaystyle M_{n+1}$
$\displaystyle=M_{n}+\delta\left[0.6M_{n}\left(1-\frac{M_{n}}{0.9}\right)+0.45M_{n}C_{n}-\frac{0.3M_{n}}{M_{n}+s}+0.4M_{n}s\right],$
$\displaystyle C_{n+1}$
$\displaystyle=C_{n}+\delta\Big{[}0.6s-0.1C_{n}-0.45M_{n}C_{n}\Big{]},$ (15)

where $\delta_{0}=1-\rho_{1}(M_{n}-M^{*})-\rho_{2}(C_{n}-C^{*})$ and
$(M^{*},C^{*})=(0.08,0.28)$ is an unstable equilibrium of the model (2)
corresponds to the given set of parameters.

Now, we have

$\displaystyle J(M^{*},C^{*},\delta_{0})=\begin{bmatrix}1.5050&0.0360\\\
-0.1260&1.2480\end{bmatrix},\quad B=\begin{bmatrix}0.0410\\\
0.0694\end{bmatrix},$

and

$\displaystyle C=[B:JB]=\begin{bmatrix}0.0410&0.0641\\\
0.0694&0.0815\end{bmatrix}.$

Since, $\det(C)=-0.0011\neq 0$, clearly $C$ is a rank 2 matrix. Hence, the
model (2) is controllable. Consider the following Jacobian matrix of the
controlled model:

$\displaystyle
J-BH=\begin{bmatrix}1.5050-0.0409\rho_{1}&0.036-0.0409\rho_{2}\\\
-0.0694\rho_{1}-0.126&1.248-0.0694\rho_{2}\end{bmatrix}.$

The characteristic equation of $J-BH$ is given by

$\displaystyle\lambda^{2}-\textrm{tr}(J-BH)\lambda+\det(J-BH)=0,$

where $\textrm{tr}(J-BH)=2.7529-0.0694\rho_{2}-0.0409\rho_{1}$ and
$\det(J-BH)=1.8827-0.1096\rho_{2}-0.0486\rho_{1}.$ Moreover, the marginal
stability of the controlled model (6) are given by

$\displaystyle\mathcal{L}_{1}=$
$\displaystyle-0.04862\rho_{1}-0.109665\rho_{2}+0.88273,$
$\displaystyle\mathcal{L}_{2}=$ $\displaystyle
0.00765\rho_{1}+0.040226\rho_{2}-0.129767,$ $\displaystyle\mathcal{L}_{3}=$
$\displaystyle-0.089578\rho_{1}-0.179106\rho_{2}+5.635693.$

Figure 18. Stability region of the controlled model (6).

######  Remark 8.

The triangular region covered by the straight lines $\mathcal{L}_{1}$,
$\mathcal{L}_{2}$, and $\mathcal{L}_{3}$ in the $\rho_{1}\rho_{2}-$plane
contains stable eigenvalues of the controlled model by taking the integral
step size $\delta$ as the control parameter with a given parametric value. The
stability region bounded by marginal lines $\mathcal{L}_{1}$,
$\mathcal{L}_{2},$ and $\mathcal{L}_{3}$ of model (6) is given in Figure 18.

######  Remark 9.

We have referenced the work of [1], who examined the stability and Hopf-
bifurcation analysis of the coral reef ODE model. Our study focuses on a
slightly modified version of this model proposed by [1]. Our analysis has
identified codimension-one period-doubling bifurcation and various complex
dynamic behaviors, including reversal period-doubling bifurcation, period-4,
8, and 24 bubble bifurcations, and chaotic behavior. These findings highlight
the richer dynamics exhibited by discrete-time models compared to continuous-
time models, contributing to our understanding of chaotic behaviors in
ecological systems.

## Conclusion

This paper examines the complex dynamics of the discrete-time coral reef
dynamical model incorporating macroalgae. We demonstrate that the model
undergoes codimension-one bifurcation Neimark-Sacker bifurcation. Furthermore,
the proposed model exhibits even more exciting dynamic behavior, including
invariant cycles and chaotic sets. These results illustrate that the discrete
model has far richer dynamics than the continuous model. We employ feedback
control to stabilize the chaotic orbits at an unstable fixed point emerging
through Neimark-Sacker bifurcation. Finally, we verify the theoretical results
numerically with the help of suitable computational software. In the future,
the authors will focus on the bifurcation analysis of a four-dimensional coral
reef mathematical model that incorporates parrotfish interaction.

## Statements and Declarations

Animal Research: This article does not contain any studies with human
participants or animals performed by any of the authors.
Author Contribution: M. Priyanka: Conceptualization, Methodology, Software,
Writing - original draft, Writing - review & editing. P. Muthukumar:
Conceptualization, Validation, Supervision, Formal analysis, Writing - review
& editing.
Code Availability: Not applicable.
Competing Interests: The authors have no competing interests to declare that
are relevant to the content of this article.
Conflict of Interest: The authors declare that they have no conflict of
interest.
Data Availability Statement: Data sharing does not apply to this article as no
datasets were generated or analyzed during the current study.
Financial/Non-financial Interests: The authors have no relevant financial or
non-financial interests to disclose.

## References

  * [1] Li, X., Wang, H., Zhang, Z., Hastings, A.: Mathematical analysis of coral reef models. Journal of Mathematical Analysis and Applications 416(1), 352-373 (2014)
  * [2] Mora, C., Graham, N. A., Nyström, M.: Ecological limitations to the resilience of coral reefs. Coral Reefs 35, 1271-1280 (2016)
  * [3] Doering, T., Wall, M., Putchim, L., Rattanawongwan, T., Schroeder, R., Hentschel, U., Roik, A.: Towards enhancing coral heat tolerance: a “microbiome transplantation” treatment using inoculations of homogenized coral tissues. Microbiome 9(1), 102 (2021)
  * [4] Ranjit, B., Biswas, S., Bhattacharyya, J., Chattopadhyay, J.: Dynamics of Zooplankton-Mediated Disease Outbreak in Coral-reef. Differential Equations and Dynamical Systems 1-29 (2023)
  * [5] Clements, C. S., Hay, M. E.: Disentangling the impacts of macroalgae on corals via effects on their microbiomes. Frontiers in Ecology and Evolution 11, 1083341 (2023)
  * [6] Jokiel, P. L., Morrissey, J. I.: Influence of size on primary production in the reef coral Pocillopora damicornis and the macroalga Acanthophora spicifera. Marine Biology 91, 15-26 (1986)
  * [7] Melis, R., Ceccherelli, G., Piazzi, L., Rustici, M.: Macroalgal forests and sea urchin barrens: Structural complexity loss, fisheries exploitation and catastrophic regime shifts. Ecological Complexity 37, 32-37 (2019)
  * [8] Misra, A. K., Tiwari, P. K., Chandra, P.: Modeling the control of algal bloom in a lake by applying some external efforts with time delay. Differential Equations and Dynamical Systems 29, 539-568 (2021)
  * [9] May, R. M.: Biological populations with nonoverlapping generations: stable points, stable cycles, and chaos. Science 186(4164), 645-647 (1974)
  * [10] Zhang, F., Tian, H., Zhao, H., Zhang, X., Shi, Q.: Spatiotemporal pattern formation in a discrete toxic-phytoplankton-zooplankton model with cross-diffusion and weak Allee effect. International Journal of Bifurcation and Chaos 32(10), 2250156 (2022)
  * [11] Han, X., Lei, C.: Bifurcation and turing instability analysis for a space-and timediscrete predator-prey system with Smith growth function. Chaos, Solitons & Fractals 166, 112910 (2023)
  * [12] Xu, W., Liu, X., Wang, H., Zhou, Y.: Event-based optimal output-feedback control of nonlinear discrete-time systems. Information Sciences 540, 414-434 (2020)
  * [13] Atabaigi, A.: Multiple Bifurcations and Dynamics of a Discrete Time Predator-Prey System with Group Defense and Non-Monotonic Functional Response. Differential Equations and Dynamical Systems 28(1), 107-132 (2020)
  * [14] Ren, J., Yu, L.: Codimension-two bifurcation, chaos and control in a discrete-time information diffusion model. Journal of Nonlinear Science 26, 1895-1931 (2016)
  * [15] Zelinka, I., Diep, Q.B., Snášel, V., Das, S., Innocenti, G., Tesi, A., Schoen, F. and Kuznetsov, N.V.: Impact of chaotic dynamics on the performance of metaheuristic optimization algorithms: An experimental analysis. Information Sciences 587, 692-719 (2022)
  * [16] Murray, J. D.: Mathematical biology: I. An introduction. Interdisciplinary applied mathematics. Mathematical Biology, Springer 17, (2002)
  * [17] Luo, X. S., Chen, G., Wang, B. H., Fang, J. Q.: Hybrid control of period-doubling bifurcation and chaos in discrete nonlinear dynamical systems. Chaos, Solitons & Fractals 18(4), 775-783 (2003)
  * [18] Kuznetsov, Y. A., Kuznetsov, I. A., Kuznetsov, Y.: Elements of applied bifurcation theory. New York: Springer 112, (1998)
  * [19] Nath, B., Kumari, N., Kumar, V., Das, K. P.: Refugia and Allee effect in prey species stabilize chaos in a tri-trophic food chain model. Differential Equations and Dynamical Systems 1-27 (2019)
  * [20] Kumari, N., Kumar, V.: Controlling chaos and pattern formation study in a tritrophic food chain model with cannibalistic intermediate predator. The European Physical Journal Plus 137(3), 1-23 (2022)
  * [21] Singh, A., Deolia, P.: Dynamical analysis and chaos control in discrete-time prey-predator model. Communications in Nonlinear Science and Numerical Simulation 90, 105313 (2020)
  * [22] Akhtar, S., Ahmed, R., Batool, M., Shah, N. A., Chung, J. D.: Stability, bifurcation and chaos control of a discretized Leslie prey-predator model. Chaos, Solitons & Fractals 152, 111345 (2021)
  * [23] Robinson, C.: Dynamical systems: stability, symbolic dynamics, and chaos. CRC press (1998)
  * [24] Layek, G. C.: An introduction to dynamical systems and chaos. New Delhi: Springer (2015)
  * [25] Carr, J.: Applications of centre manifold theory. Springer Science & Business Media (2012)
  * [26] Guckenheimer, J. Holmes, P.: Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Springer Science & Business Media (2013)
  * [27] Fattahpour, H., Zangeneh, H. R., Wang, H.: Dynamics of coral reef models in the presence of parrotfish. Natural Resource Modeling 32(2), e12202 (2019)
  * [28] Blackwood, J. C., Hastings, A.: The effect of time delays on Caribbean coral–algal interactions. Journal of Theoretical Biology 273(1), 37-43 (2011)
  * [29] Blackwood, J. C., Hastings, A, Mumby, P. J.: The effect of fishing on hysteresis in Caribbean coral reefs. Theoretical ecology 5, 105-114 (2012)
  * [30] González-Rivero, M., Yakob, L., Mumby, P. J.: The role of sponge competition on coral reef alternative steady states. Ecological Modelling 222(11), 1847-1853 (2011).
  * [31] Luo, A. C.: Regularity and complexity in dynamical systems. New York: Springer (2012)
  * [32] Salman, S. M., Yousef, A. M., Elsadany, A. A.: Stability, bifurcation analysis and chaos control of a discrete predator-prey system with square root functional response. Chaos, Solitons & Fractals 93, 20-31 (2016)
  * [33] He, Z., Lai, X.: Bifurcation and chaotic behavior of a discrete-time predator–prey system. Nonlinear Analysis: Real World Applications 12(1), 403-417 (2011)
  * [34] Elabbasy, E. M., Elsadany, A. A., Zhang, Y.: Bifurcation analysis and chaos in a discrete reduced Lorenz system. Applied Mathematics and Computation 228, 184-194 (2014)
  * [35] Liu, W., Jiang, Y.: Flip bifurcation and Neimark–Sacker bifurcation in a discrete predator–prey model with harvesting. International Journal of Biomathematics 13(01), 1950093 (2020)
  * [36] Abdelaziz, M. A., Ismail, A, I., Abdullah, F. A., Mohd, M. H.: Bifurcations and chaos in a discrete SI epidemic model with fractional order. Advances in Difference Equations 2018(1), 1-19 (2018)
  * [37] Ott, E., Grebogi, C., Yorke, J. A.: Controlling chaos. Physical review letters 64(11), 1196 (1990)
  * [38] Luo, A. C.: Bifurcation and stability in nonlinear dynamical systems. Springer International Publishing (2019)
  * [39] Yuan, L. G., Yang, Q. B.: Bifurcation, invariant curve and hybrid control in a discrete-time predator–prey system. Applied Mathematical Modelling 39(8), 2345-2362 (2015)

Appendix A

Proof of Theorem 2

Let $x=M-M^{*},$ $y=C-C^{*}$. Then shift the interior equilibrium point
$E^{*}(M^{*},C^{*})$ of map (4.2) into the origin which given by

$\displaystyle\left\\{\begin{array}[]{l}x\longrightarrow
i_{100}x+i_{010}y+i_{200}x^{2}+i_{110}xy+i_{300}x^{3}+i_{101}x\delta_{1}^{*}+i_{011}y\delta_{1}^{*}+i_{201}x^{2}\delta_{1}^{*}\\\
\qquad\quad\\!\\!+i_{111}xy\delta_{1}^{*}+O\Big{(}(|x|+|y|+|\delta_{1}^{*}|)^{4}\Big{)},\\\
y\longrightarrow
j_{100}x+j_{010}y+j_{110}xy+j_{101}x\delta_{1}^{*}+j_{011}y\delta_{1}^{*}+j_{111}xy\delta_{1}^{*}+O\Big{(}(|x|+|y|+|\delta_{1}^{*}|)^{4}\Big{)},\end{array}\right.$
(19)

where $i_{100}=1+\delta r-\frac{2\delta rM}{k}+\delta aC-\frac{\delta
Sg}{(M+S)^{2}}+\delta\gamma S,$ $i_{010}=\delta aM,$ $i_{200}=\frac{-\delta
r}{k}+\frac{\delta Sg}{(M+S)^{3}},$ $i_{110}=\delta a,$
$i_{101}=r-\frac{2rM}{k}+aC-\frac{Sg}{(M+S)^{2}}+\gamma S,$
$i_{300}=\frac{-\delta Sg}{(M+S)^{4}},$ $i_{011}=aM,$
$i_{201}=\frac{-r}{k}+\frac{Sg}{(M+S)^{3}},$ $i_{111}=a,$ $j_{100}=-\delta
aC,$ $j_{010}=1+\delta\alpha S-\delta d-\delta aM,$ $j_{110}=-\delta a,$
$j_{101}=-aC,$ $j_{011}=\alpha S-d-aM,$ $j_{111}=-a,$ and $\delta=\delta_{1}$.
Assume $i_{010}\neq 0.$

The invertible matrix is constructed as,

$\displaystyle P_{1}=\begin{bmatrix}i_{010}&i_{010}\\\
-1-i_{100}&\lambda_{2}-i_{100}\end{bmatrix}$

and by using this translation

$\displaystyle\begin{bmatrix}x\\\
y\end{bmatrix}=P_{1}\begin{bmatrix}\tilde{M}\\\ \tilde{C}\end{bmatrix},$ (20)

which implies

$\displaystyle\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}=\frac{1}{|P_{1}|}\begin{bmatrix}\lambda_{2}-i_{100}&-i_{010}\\\
1+i_{100}&i_{010}\end{bmatrix}\begin{bmatrix}x\\\ y\end{bmatrix},$

where $|P_{1}|=i_{010}(1+\lambda_{2})$. Then

$\displaystyle\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}=\frac{1}{i_{010}(1+\lambda_{2})}\begin{bmatrix}(\lambda_{2}-i_{100})x-i_{010}y\\\
(1+i_{100})x+i_{010}y\end{bmatrix}.$

By map (2), we have

$\displaystyle\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}\longrightarrow\frac{1}{i_{010}(1+\lambda_{2})}\begin{bmatrix}(\lambda_{2}-i_{100})\Big{[}i_{100}x+i_{010}y+i_{200}x^{2}+i_{110}xy+i_{300}x^{3}+i_{101}x\delta_{1}^{*}+i_{011}y\delta_{1}^{*}&\\\
+i_{201}x^{2}\delta_{1}^{*}+i_{111}xy\delta_{1}^{*}+O\Big{(}(|x|+|y|+|\delta_{1}^{*}|)^{4}\Big{)}\Big{]}&\\\
-i_{010}\Big{[}j_{100}x+j_{010}y+j_{110}xy+j_{101}x\delta_{1}^{*}+j_{011}y\delta_{1}^{*}&\\\
\vspace{0.5cm}+j_{111}xy\delta_{1}^{*}+O\Big{(}(|x|+|y|+|\delta_{1}^{*}|)^{4}\Big{)}\Big{]}&\\\
(1+i_{100})\Big{[}i_{100}x+i_{010}y+i_{200}x^{2}+i_{110}xy+i_{300}x^{3}+i_{101}x\delta_{1}^{*}+i_{011}y\delta_{1}^{*}&\\\
+i_{201}x^{2}\delta_{1}^{*}+i_{111}xy\delta_{1}^{*}+O\Big{(}(|x|+|y|+|\delta_{1}^{*}|)^{4}\Big{)}\Big{]}&\\\
+i_{010}\Big{[}j_{100}x+j_{010}y+j_{110}xy+j_{101}x\delta_{1}^{*}+j_{011}y\delta_{1}^{*}&\\\
+j_{111}xy\delta_{1}^{*}+O\Big{(}(|x|+|y|+|\delta_{1}^{*}|)^{4}\Big{)}\Big{]}&\end{bmatrix}.$

Now, the above map (19) can be rewritten as

$\displaystyle\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}\longrightarrow\begin{bmatrix}-1&0\\\
0&\lambda_{2}\end{bmatrix}\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}+\begin{bmatrix}g_{1}(\tilde{M},\tilde{C},\delta_{1}^{*})\\\
g_{2}(\tilde{M},\tilde{C},\delta_{1}^{*})\end{bmatrix},$ (21)

where

$\displaystyle
g_{1}(\tilde{M},\tilde{C},\delta_{1}^{*})=\frac{1}{i_{010}(1+\lambda_{2})}\begin{bmatrix}\Big{[}(\lambda_{2}-i_{100})i_{110}-i_{010}j_{110}\Big{]}xy+\Big{[}(\lambda_{2}-i_{100})i_{101}-i_{010}j_{101}\Big{]}x\delta_{1}^{*}&\\\
+\Big{[}(\lambda_{2}-i_{100})i_{011}-i_{010}j_{011}\Big{]}y\delta_{1}^{*}+\Big{[}(\lambda_{2}-i_{100})i_{111}-i_{010}j_{111}\Big{]}xy\delta_{1}^{*}&\\\
+(\lambda_{2}-i_{100})i_{200}x^{2}+(\lambda_{2}-i_{100})i_{201}x^{2}\delta_{1}^{*}&\\\
+(\lambda_{2}-i_{100})i_{300}x^{3}+O\Big{(}(|x|+|y|+|\delta_{1}^{*}|)^{4}\Big{)}&\end{bmatrix},$

$\displaystyle
g_{2}(\tilde{M},\tilde{C},\delta_{1}^{*})=\frac{1}{i_{010}(1+\lambda_{2})}\begin{bmatrix}\Big{[}(1+i_{100})i_{110}+i_{010}j_{110}\Big{]}xy+\Big{[}(1+i_{100})i_{101}+i_{010}j_{101}\Big{]}x\delta_{1}^{*}&\\\
+\Big{[}(1+i_{100})i_{011}+i_{010}j_{011}\Big{]}y\delta_{1}^{*}+\Big{[}(1+i_{100})i_{111}+i_{010}j_{111}\Big{]}xy\delta_{1}^{*}&\\\
+(1+i_{100})i_{200}x^{2}+(1+i_{100})i_{201}x^{2}\delta_{1}^{*}&\\\
+(1+i_{100})i_{300}x^{3}+O\Big{(}(|x|+|y|+|\delta_{1}^{*}|)^{4}\Big{)}&\end{bmatrix},$

and from (20), we have

$\displaystyle\begin{bmatrix}x\\\
y\end{bmatrix}=\begin{bmatrix}i_{010}\tilde{M}+i_{010}\tilde{C}\\\
-(1+i_{100})\tilde{M}+(\lambda_{2}-i_{100})\tilde{C}\end{bmatrix},$

which implies

$\displaystyle x$ $\displaystyle=i_{010}(\tilde{M}+\tilde{C}),$ $\displaystyle
y$ $\displaystyle=-(1+i_{100})\tilde{M}+(\lambda_{2}-i_{100})\tilde{C},$
$\displaystyle xy$
$\displaystyle=i_{010}\Big{[}-(1+i_{100})\tilde{M}^{2}+(\lambda_{2}-i_{100})\tilde{C}^{2}+(\lambda_{2}-1-2i_{100})\tilde{M}\tilde{C}\Big{]},$
$\displaystyle x^{2}$
$\displaystyle=i_{010}^{2}(\tilde{M}^{2}+\tilde{C}^{2}+2\tilde{M}\tilde{C}),$
$\displaystyle x^{3}$
$\displaystyle=i_{010}^{3}(\tilde{M}^{3}+\tilde{C}^{3}+3\tilde{M}^{2}\tilde{C}+3\tilde{M}\tilde{C}^{2}).$

Next, the center manifold $C_{M}(0,0,0)$ of (21) about the origin in a small
neighborhood of $\delta_{1}^{*}$ is determined. We can conclude that there is
a center manifold $C_{M}(0,0,0)$ from the center manifold theorem given as
follows:

$\displaystyle
C_{M}(0,0,0)=\left\\{(\tilde{M},\tilde{C},\delta_{1}^{*})\in\mathbb{R}^{3},\tilde{C}=\mathfrak{F}(\tilde{M},\delta_{1}^{*}),\mathfrak{F}(0,0)=0,D\mathfrak{F}(0,0)=0\right\\},$

for $\tilde{M}$ and $\delta_{1}^{*}$ sufficiently small. Assume the following
center manifold form:

$\displaystyle\mathfrak{F}(\tilde{M},\delta_{1}^{*})=a_{1}\tilde{M}^{2}+a_{2}\tilde{M}\delta_{1}^{*}+a_{3}{\delta_{1}^{*}}^{2}+O\Big{(}(|\tilde{M}|+|\delta_{1}^{*}|)^{3}\Big{)}.$
(22)

Given $C_{M}(0,0,0)$ must satisfy

$\displaystyle\mathfrak{F}\Big{(}-\tilde{M}+g_{1}(\tilde{M},\mathfrak{F}(\tilde{M},\delta_{1}^{*}),\delta_{1}^{*}),\delta_{1}^{*}\Big{)}=\lambda_{2}\mathfrak{F}(\tilde{M},\delta_{1}^{*})+g_{2}\Big{(}\tilde{M},\mathfrak{F}(\tilde{M},\delta_{1}^{*}),\delta_{1}^{*}\Big{)}.$
(23)

Substituting (22) into (23), and comparing coefficients of similar powers in
obtained (23), we get that

$\displaystyle a_{1}=$
$\displaystyle\frac{1}{i_{010}(1-\lambda_{2}^{2})}\Big{[}-[(1+i_{100})i_{110}+i_{010}j_{110}]i_{010}(1+i_{100})+(1+i_{100})i_{200}i_{010}^{2}\Big{]},$
$\displaystyle a_{2}=$
$\displaystyle\frac{-1}{i_{010}(1+\lambda_{2})^{2}}\Big{[}[(1+i_{100})i_{101}+i_{010}j_{101}]i_{010}-[(1+i_{100})i_{011}+i_{010}j_{011}](1+i_{100})\Big{]},$
$\displaystyle a_{3}=$ $\displaystyle 0.$

The map (21) that is restricted to the central manifold $C_{M}(0,0,0)$ is
taken into consideration:

$\displaystyle
G:\tilde{M}\longrightarrow-\tilde{M}+b_{1}\tilde{M}^{2}+b_{2}\tilde{M}\delta_{1}^{*}+b_{3}\tilde{M}^{2}\delta_{1}^{*}+b_{4}\tilde{M}{\delta_{1}^{*}}^{2}+b_{5}\tilde{M}^{3}+O\Big{(}(|\tilde{M}|+|\delta_{1}^{*}|)^{4}\Big{)},$
(24)

where

$\displaystyle b_{1}=\frac{1}{i_{010}(1+\lambda_{2})}$
$\displaystyle\Big{[}-[(\lambda_{2}-i_{100})i_{110}-i_{010}j_{110}]i_{010}(1+i_{100})+(\lambda_{2}-i_{100})i_{200}i_{010}^{2}\Big{]},$
$\displaystyle b_{2}=\frac{1}{i_{010}(1+\lambda_{2})}$
$\displaystyle\Big{[}[(1+i_{100})i_{101}+i_{010}j_{101}]i_{010}-[(1+i_{100})i_{011}+i_{010}j_{011}](1+i_{100})\Big{]},$
$\displaystyle b_{3}=\frac{1}{i_{010}(1+\lambda_{2})}$
$\displaystyle\Big{[}[(\lambda_{2}-i_{100})i_{110}-i_{010}j_{110}]i_{010}(\lambda_{2}-1-2i_{100})a_{2}+[(\lambda_{2}-i_{100})i_{101}-i_{010}j_{101}]i_{010}a_{1}$
$\displaystyle+[(\lambda_{2}-i_{100})i_{011}-i_{010}j_{011}](\lambda_{2}-i_{100})a_{1}-[(\lambda_{2}-i_{100})i_{111}-i_{010}j_{111}]i_{010}(1+i_{100})$
$\displaystyle+2(\lambda_{2}-i_{100})i_{200}i_{010}^{2}a_{2}+(\lambda_{2}-i_{100})i_{201}i_{010}^{2}\Big{]},$
$\displaystyle b_{4}=\frac{1}{i_{010}(1+\lambda_{2})}$
$\displaystyle\Big{[}[(\lambda_{2}-i_{100})i_{101}-i_{010}j_{101}]a_{2}i_{010}+[(\lambda_{2}-i_{100})i_{011}-i_{010}j_{011}](\lambda_{2}-i_{100})a_{2}\Big{]},$
$\displaystyle b_{5}=\frac{1}{i_{010}(1+\lambda_{2})}$
$\displaystyle\Big{[}[(\lambda_{2}-i_{100})i_{110}-i_{010}j_{110}]i_{010}(\lambda_{2}-1-2i_{100})a_{1}+2(\lambda_{2}-i_{100})i_{200}i_{010}^{2}a_{1}$
$\displaystyle+(\lambda_{2}-i_{100})i_{300}i_{010}^{3}\Big{]}.$

Two discriminatory values $\Omega_{1}$ and $\Omega_{2}$, must not be zero for
the map (24) to experience a flip bifurcation.

$\displaystyle\Omega_{1}$
$\displaystyle=\left.\left(2\frac{\partial^{2}G}{\partial\tilde{M}\partial\delta_{1}^{*}}+\frac{\partial
G}{\partial\delta_{1}^{*}}\frac{\partial
G}{\partial\tilde{M}}\right)\right|_{(0,0)}=2b_{2},$ $\displaystyle\Omega_{2}$
$\displaystyle=\left.\left(\frac{1}{2}\left(\frac{\partial^{2}G}{\partial\tilde{M}^{2}}\right)^{2}+\frac{1}{3}\left(\frac{\partial^{3}G}{\partial\tilde{M}^{3}}\right)\right)\right|_{(0,0)}=2b_{1}^{2}+2b_{5}.$

Appendix B

Proof of Theorem 3

Let $x=M-M^{*}$, $y=C-C^{*}$. Then shift $E^{*}(M^{*},C^{*})$ of map (4.2)
into the origin which is given by

$\displaystyle\begin{bmatrix}x\\\
y\end{bmatrix}\longrightarrow\begin{bmatrix}i_{100}x+i_{010}y+i_{200}x^{2}+i_{020}y^{2}+i_{110}xy+i_{120}xy^{2}+i_{030}y^{3}+O\Big{(}(|x|+|y|)^{4}\Big{)}\\\
j_{100}x+j_{010}y+j_{020}y^{2}+j_{110}xy+O\Big{(}(|x|+|y|)^{4}\Big{)}\end{bmatrix},$
(25)

where $i_{100}$, $i_{010}$, $i_{200},$ $i_{020},$ $i_{110}$, $i_{120},$
$i_{300},$ $j_{100}$, $j_{010}$, $j_{020},$ $j_{110}$, are given previously by
substituting $\delta$ for $\delta_{2}+\delta_{2}^{*}.$

Consider the following characteristic equation for the linearization of the
map (25) about the origin:

$\displaystyle\lambda^{2}+\mathcal{A}(\delta_{2}^{*})\lambda+\mathcal{B}(\delta_{2}^{*})=0,$

where

$\displaystyle\mathcal{A}(\delta_{2}^{*})$
$\displaystyle=-2-\mathcal{U}(\delta_{2}+\delta_{2}^{*}),$
$\displaystyle\mathcal{B}(\delta_{2}^{*})$
$\displaystyle=1+\mathcal{U}(\delta_{2}+\delta_{2}^{*})+\mathcal{V}(\delta_{2}+\delta_{2}^{*})^{2}.$

Given that the parameters are
$(r,k,a,g,\gamma,\alpha,d,\delta_{2})\in\mathcal{N}$, the eigenvalues of
$(0,0)$ are a pair of complex conjugate numbers called $\lambda$ and
$\bar{\lambda}$ with modulus 1 by Theorem 1, where

$\displaystyle\lambda,\bar{\lambda}=\frac{-\mathcal{A}(\delta_{2}^{*})\pm
i\sqrt{4\mathcal{B}(\delta_{2}^{*})-\mathcal{A}^{2}(\delta_{2}^{*})}}{2}$

and

$\displaystyle|\lambda|_{\delta_{2}^{*}=0}=\sqrt{\mathcal{B}(0)}=1,\quad\left.\frac{d|\lambda|}{d\delta_{2}^{*}}\right|_{\delta_{2}^{*}=0}=-\frac{\mathcal{U}}{2}>0.$

Additionally, it is required that when $\delta_{2}^{*}=0$, $\lambda^{r},$
$\bar{\lambda}^{r}\neq 1$ $(r=1,2,3,4)$ which is equivalent to
$\mathcal{A}(0)\neq-2,0,1,2$. Note that
$(r,k,a,g,\gamma,\alpha,d,\delta_{2})\in\mathcal{N}.$ Thus,
$\mathcal{A}(0)\neq-2,2$, we only need to show that $\mathcal{A}(0)\neq 0,1,$
which implies

$\displaystyle\mathcal{U}^{2}\neq 2\mathcal{V},3\mathcal{V}.$ (26)

When $\delta_{2}^{*}=0$ and (26) is true, the eigenvalues of fixed point
$(0,0)$ of (25) do not lie in the position where the coordinate axes and unit
circle coincide. Then, we investigate the normal form of (25) at
$\delta_{2}^{*}=0$.

Let $\sigma_{1}=1+\frac{\mathcal{U}\delta_{2}}{2}$ and
$\sigma_{2}=\frac{\delta_{2}}{2}\sqrt{4\mathcal{V}-\mathcal{U}^{2}}$. The
invertible matrix is constructed as,

$\displaystyle P_{2}=\begin{bmatrix}i_{010}&0\\\
\sigma_{1}-i_{100}&-\sigma_{2}\end{bmatrix}$

and by using this translation

$\displaystyle\begin{bmatrix}x\\\
y\end{bmatrix}=P_{2}\begin{bmatrix}\tilde{M}\\\ \tilde{C}\end{bmatrix},$ (27)

which implies

$\displaystyle\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}=\frac{1}{|P_{2}|}\begin{bmatrix}-\sigma_{2}&0\\\
i_{100}-\sigma_{1}&i_{010}\end{bmatrix}\begin{bmatrix}x\\\ y\end{bmatrix},$

where $|P_{2}|=-\sigma_{2}i_{010}$. Now,

$\displaystyle\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}=\frac{1}{-\sigma_{2}i_{010}}\begin{bmatrix}-\sigma_{2}x\\\
(i_{100}-\sigma_{1})x+i_{010}y\end{bmatrix}.$

By map (2), we have

$\displaystyle\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}\longrightarrow\frac{1}{-\sigma_{2}i_{010}}\begin{bmatrix}-\sigma_{2}\Big{[}i_{100}x+i_{010}y+i_{200}x^{2}+i_{020}y^{2}+i_{110}xy+i_{120}xy^{2}\\\
\vspace{0.5cm}+i_{030}y^{3}+O\Big{(}(|x|+|y|)^{4}\Big{)}\Big{]}\\\
(i_{100}-\sigma_{1})\Big{[}i_{100}x+i_{010}y+i_{200}x^{2}+i_{020}y^{2}+i_{110}xy+i_{120}xy^{2}\\\
+i_{030}y^{3}+O\Big{(}(|x|+|y|)^{4}\Big{)}\Big{]}\\\
+i_{010}\Big{[}j_{100}x+j_{010}y+j_{020}+j_{110}xy+O\Big{(}(|x|+|y|)^{4}\Big{)}\Big{]}\end{bmatrix}.$

Now, the above matrix can be rewritten as

$\displaystyle\begin{bmatrix}\tilde{M}\\\
\tilde{C}\end{bmatrix}\\!\longrightarrow\\!\frac{1}{-\sigma_{2}i_{010}}\\!\begin{bmatrix}-\sigma_{2}i_{100}&-\sigma_{2}i_{010}\\\
(i_{100}-\sigma_{1})i_{100}+i_{010}j_{100}&(i_{100}-\sigma_{1})i_{010}+i_{010}j_{010}\end{bmatrix}\begin{bmatrix}x\\\
y\end{bmatrix}\\!+\\!\begin{bmatrix}g_{3}(\tilde{M},\tilde{C})\\\
g_{4}(\tilde{M},\tilde{C})\end{bmatrix},$ (28)

where

$\displaystyle g_{3}(\tilde{M},\tilde{C})=\frac{1}{-\sigma_{2}i_{010}}$
$\displaystyle\Bigg{[}-\sigma_{2}\Big{[}i_{200}x^{2}++i_{020}y^{2}+i_{110}xy+i_{120}xy^{2}+i_{030}y^{3}+O\Big{(}(|x|+|y|)^{4}\Big{)}\Big{]}\Bigg{]},$
$\displaystyle g_{4}(\tilde{M},\tilde{C})=\frac{1}{-\sigma_{2}i_{010}}$
$\displaystyle\Bigg{[}(i_{100}-\sigma_{1})\Big{[}i_{200}x^{2}++i_{020}y^{2}+i_{110}xy+i_{120}xy^{2}+i_{030}y^{3}+O\Big{(}(|x|+|y|)^{4}\Big{)}\Big{]}$
$\displaystyle+i_{010}\Big{[}j_{020}y^{2}+j_{110}xy+O\Big{(}(|x|+|y|)^{4}\Big{)}\Big{]}\Bigg{]},$

and from (27), we have

$\displaystyle\begin{bmatrix}x\\\
y\end{bmatrix}=\begin{bmatrix}i_{010}\tilde{M}\\\
(\sigma_{2}-i_{100})\tilde{M}-\sigma_{2}\tilde{C}\end{bmatrix},$

which implies,

$\displaystyle x$ $\displaystyle=i_{010}\tilde{M},$ $\displaystyle y$
$\displaystyle=(\sigma_{1}-i_{100})\tilde{M}-\sigma_{2}\tilde{C},$
$\displaystyle xy$
$\displaystyle=i_{010}(\sigma_{1}-i_{100})\tilde{M}^{2}-\sigma_{2}i_{010}\tilde{M}\tilde{C},$
$\displaystyle x^{2}$ $\displaystyle=i_{010}^{2}\tilde{M}^{2},$ $\displaystyle
y^{2}$
$\displaystyle=(\sigma_{1}-i_{100})^{2}\tilde{M}^{2}+\sigma_{2}^{2}\tilde{C}^{2}-2\sigma_{2}(\sigma_{1}-i_{100})\tilde{M}\tilde{C},$
$\displaystyle xy^{2}$
$\displaystyle=i_{010}(\sigma_{1}-i_{100})^{2}\tilde{M}^{3}+i_{010}\sigma_{2}^{2}\tilde{M}\tilde{C}^{2}-2\sigma_{2}i_{010}(\sigma_{1}-i_{100})\tilde{M}^{2}\tilde{C},$
$\displaystyle y^{3}$
$\displaystyle=(\sigma_{1}-i_{100})^{3}\tilde{M}^{3}-\sigma_{2}^{3}\tilde{C}^{3}+3\sigma_{2}^{2}(\sigma_{1}-i_{100})\tilde{M}\tilde{C}^{2}-3\sigma_{2}(\sigma_{1}-i_{100})^{2}\tilde{M}^{2}\tilde{C},$

and

$\displaystyle
g_{3\tilde{M}\tilde{M}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}-\sigma_{2}[2i_{200}i_{010}^{2}+2i_{020}(\sigma_{1}-i_{100})^{2}+2i_{110}i_{010}(\sigma_{1}-i_{100})+i_{120}[6i_{010}(\sigma_{1}-i_{100})^{2}\tilde{M}$
$\displaystyle\qquad\qquad\qquad\quad-4\sigma_{2}i_{010}(\sigma_{1}-i_{100})\tilde{C}]+i_{030}[6(\sigma_{1}-i_{100})^{3}\tilde{M}-6\sigma_{2}(\sigma_{1}-i_{100})^{2}\tilde{C}]]\Big{]},$
$\displaystyle
g_{3\tilde{C}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}-\sigma_{2}[2i_{020}\sigma_{2}^{2}+2i_{120}i_{010}\sigma_{2}^{2}\tilde{M}+i_{030}[-6\sigma_{2}^{3}\tilde{C}+6\sigma_{2}^{2}(\sigma_{1}-i_{100})\tilde{M}]]\Big{]},$
$\displaystyle
g_{3\tilde{M}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}-\sigma_{2}[-2i_{020}\sigma_{2}(\sigma_{1}-i_{100})-i_{110}i_{010}\sigma_{2}+i_{120}[2i_{010}\sigma_{2}^{2}\tilde{C}-4\sigma_{2}i_{010}(\sigma_{1}-i_{100})\tilde{M}]$
$\displaystyle\qquad\qquad\qquad\quad+i_{030}[6\sigma_{2}^{2}(\sigma_{1}-i_{100})\tilde{C}-6\sigma_{2}(\sigma_{1}-i_{100})^{2}\tilde{M}]]\Big{]},$
$\displaystyle
g_{3\tilde{M}\tilde{M}\tilde{M}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}-\sigma_{2}[6i_{120}i_{010}(\sigma_{1}-i_{100})^{2}+6i_{030}(\sigma_{1}-i_{100})^{3}]\Big{]},$
$\displaystyle
g_{3\tilde{M}\tilde{C}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}-\sigma_{2}[2i_{120}i_{010}\sigma_{2}^{2}+6i_{030}\sigma_{2}^{2}(\sigma_{1}-i_{100})]\Big{]},$
$\displaystyle
g_{3\tilde{M}\tilde{M}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}-\sigma_{2}[-4i_{120}\sigma_{2}i_{010}(\sigma_{1}-i_{100})-6i_{030}\sigma_{2}(\sigma_{1}-i_{100})^{2}]\Big{]},$
$\displaystyle
g_{3\tilde{C}\tilde{C}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}-\sigma_{2}[-6i_{030}\sigma_{2}^{3}]\Big{]},$
$\displaystyle
g_{4\tilde{M}\tilde{M}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}(i_{100}-\sigma_{1})[2i_{010}^{2}i_{200}+2i_{020}(\sigma_{1}-i_{100})^{2}+2i_{110}i_{010}(\sigma_{1}-i_{100})$
$\displaystyle\qquad\qquad\qquad\quad+i_{120}[6i_{010}(\sigma_{1}-i_{100})^{2}\tilde{M}-4\sigma_{2}i_{010}(\sigma_{1}-i_{100})\tilde{C}]$
$\displaystyle\qquad\qquad\qquad\quad+i_{030}[6(\sigma_{1}-i_{100})^{3}\tilde{M}-6\sigma_{2}(\sigma_{1}-i_{100})^{2}\tilde{C}]]$
$\displaystyle\qquad\qquad\qquad\quad+i_{010}[2j_{020}(\sigma_{1}-i_{100})^{2}+2j_{110}i_{010}(\sigma_{1}-i_{100})]\Big{]},$
$\displaystyle
g_{4\tilde{C}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}(i_{100}-\sigma_{1})[2i_{020}\sigma_{2}^{2}+2i_{120}i_{010}\sigma_{2}^{2}\tilde{M}-6i_{030}\sigma_{2}^{3}\tilde{C}+6i_{030}\sigma_{2}^{2}(\sigma_{1}-i_{100})\tilde{M}]$
$\displaystyle\qquad\qquad\qquad\quad+i_{010}[2j_{020}\sigma_{2}^{2}]\Big{]},$
$\displaystyle
g_{4\tilde{M}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}(i_{100}-\sigma_{1})[-2i_{020}\sigma_{2}(\sigma_{1}-i_{100})-i_{110}i_{010}\sigma_{2}+2i_{120}i_{010}\sigma_{2}^{2}\tilde{C}$
$\displaystyle\qquad\qquad\qquad\quad-4i_{120}\sigma_{2}i_{010}(\sigma_{1}-i_{100})\tilde{M}+6i_{030}\sigma_{2}^{2}(\sigma_{1}-i_{100})\tilde{C}-6i_{030}\sigma_{2}(\sigma_{1}-i_{100})^{2}\tilde{M}]$
$\displaystyle\qquad\qquad\qquad\quad+i_{010}[-2j_{020}\sigma_{2}(\sigma_{1}-i_{100})-j_{110}i_{010}\sigma_{2}]\Big{]},$
$\displaystyle
g_{4\tilde{M}\tilde{M}\tilde{M}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}(i_{100}-\sigma_{1})[6i_{120}i_{010}(\sigma_{1}-i_{100})^{2}+6i_{030}(\sigma_{1}-i_{100})^{3}]\Big{]},$
$\displaystyle
g_{4\tilde{C}\tilde{C}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}(i_{100}-\sigma_{1})[-6i_{030}\sigma_{2}^{3}]\Big{]},$
$\displaystyle
g_{4\tilde{M}\tilde{M}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}(i_{100}-\sigma_{1})[2i_{120}i_{010}\sigma_{2}^{2}+6i_{030}\sigma_{2}^{2}(\sigma_{1}-i_{100})]\Big{]},$
$\displaystyle
g_{4\tilde{M}\tilde{C}\tilde{C}}=\frac{-1}{\sigma_{2}i_{010}}\Big{[}(i_{100}-\sigma_{1})[-4i_{120}\sigma_{2}i_{010}(\sigma_{1}-i_{100})-6i_{030}\sigma_{2}(\sigma_{1}-i_{100})^{2}]\Big{]}.$

The following discriminatory value $\Psi$, must not be zero for the map (28)
to experience Neimark-Sacker bifurcation:

$\displaystyle\Psi=\left.\left[-\textrm{Re}\left(\frac{(1-2\lambda)\bar{\lambda}^{2}}{1-\lambda}\varsigma_{1}\varsigma_{2}\right)-\frac{1}{2}|\varsigma_{2}|^{2}-|\varsigma_{3}|^{2}+\textrm{Re}(\bar{\lambda}\varsigma_{4})\right]\right|_{\delta_{2}^{*}=0},$

where

$\displaystyle\varsigma_{1}$
$\displaystyle=\frac{1}{8}\Big{[}(g_{3\tilde{M}\tilde{M}}-g_{3\tilde{C}\tilde{C}}+2g_{4\tilde{M}\tilde{C}})+i(g_{4\tilde{M}\tilde{M}}-g_{4\tilde{C}\tilde{C}}-2g_{3\tilde{M}\tilde{C}})\Big{]},$
$\displaystyle\varsigma_{2}$
$\displaystyle=\frac{1}{4}\Big{[}(g_{3\tilde{M}\tilde{M}}+g_{3\tilde{C}\tilde{C}})+i(g_{4\tilde{M}\tilde{M}}+g_{4\tilde{C}\tilde{C}})\Big{]},$
$\displaystyle\varsigma_{3}$
$\displaystyle=\frac{1}{8}\Big{[}(g_{3\tilde{M}\tilde{M}}-g_{3\tilde{C}\tilde{C}}-2g_{4\tilde{M}\tilde{C}})+i(g_{4\tilde{M}\tilde{M}}-g_{4\tilde{C}\tilde{C}}+2g_{3\tilde{M}\tilde{C}})\Big{]},$
$\displaystyle\varsigma_{4}$
$\displaystyle=\frac{1}{16}\Big{[}(g_{3\tilde{M}\tilde{M}\tilde{M}}+g_{3\tilde{M}\tilde{C}\tilde{C}}+g_{4\tilde{M}\tilde{M}\tilde{C}}+g_{4\tilde{C}\tilde{C}\tilde{C}})+i(g_{4\tilde{M}\tilde{M}\tilde{M}}+g_{4\tilde{M}\tilde{C}\tilde{C}}-g_{3\tilde{M}\tilde{M}\tilde{C}}-g_{3\tilde{C}\tilde{C}\tilde{C}})\Big{]}.$


# Ghost Image Processing

Harry Penketh  William L. Barnes  and Jacopo Bertolotti  Department of Physics
and Astronomy, University of Exeter, Exeter, Devon, EX4 4QL, UK
<EMAIL_ADDRESS>
###### Abstract

In computational ghost imaging the object is illuminated with a sequence of
known patterns, and the scattered light is collected using a detector that has
no spatial resolution. Using those patterns and the total intensity
measurement from the detector, one can reconstruct the desired image. Here we
study how the reconstructed image is modified if the patterns used for the
reconstruction are not the same as the illumination patterns, and show that
one can choose how to illuminate the object, such that the reconstruction
process behaves like a spatial filtering operation on the image. The ability
to measure directly a processed image, allows one to bypass the post-
processing steps, and thus avoid any noise amplification they imply. As a
simple example we show the case of an edge-detection filter.

††journal: osajournal††articletype: Research Article

##  1 Introduction

Ghost imaging relies on the combination of two signals which individually are
insufficient for image formation [1, 2, 3]: the sequence of patterns
illuminating the object, and the transmitted (or scattered) light, measured
with a single element (bucket) detector [4, 5, 6, 7, 8, 9, 10, 11].
In computational ghost imaging one has a great deal of control over the choice
of the projected patterns, allowing one to tailor them based on a knowledge of
the nature of the object. In Principle Component Analysis the illuminating
wavefront is designed to match the principle components of the object [12,
13], and in other adaptive imaging works the spatial resolution of the
wavefront is enhanced locally in response to the detection of high-frequency
regions of the object [14, 15]. It is however less common to see the
illumination basis modified in response to the way in which the image is to be
processed [16].
In this paper we show that any post-processing step which can be described by
a matrix multiplication with the image, such as convolution with an image
filter, can be incorporated into the illumination basis. Doing so enables one
to avoid image noise amplification by the filtering process, at the cost of
increasing the complexity of the projected patterns. We demonstrate this
technique experimentally for a basic edge-detection filter in a modified
raster basis, and compare the resulting signal-to-noise ratios (SNRs) with
those obtained via post-processing with the same filter. We also discuss a
theoretical method to predict the performance of an arbitrary filter.

##  2 Changing the Illumination Basis

The ghost imaging measurement process can be described as follows in Bra-Ket
notation, with $N\times N$ pixel projection patterns or images represented as
$N^{2}\times 1$ element column vectors for convenience. The reconstructed
image $\ket{I}$ of the object $\ket{O}$ can be written

$\ket{I}=\sum_{j}\braket{\psi_{j}}{O}\ket{\psi_{j}},$ (1)

where $\ket{\psi_{j}}$ is the $j$th pattern in the basis $\Psi$ illuminating
the object. The inner product $\braket{\psi_{j}}{O}$, which becomes the
weighting coefficient of $\ket{\psi_{j}}$ in the reconstruction, measures the
spatial overlap between the projected pattern and object and is recorded with
a bucket detector.
Typically the illumination basis $\Psi$ is the same basis in which the image
is reconstructed, but this need not be the case. A change from an illumination
basis $\Psi$ to a new basis $\Phi$ can be written as $\Phi=B\,\Psi$, where B
is the matrix that performs the basis change. If one makes this substitution
for the illuminating basis in eq.1

$\ket{I}=\sum_{j}\braket{\phi_{j}}{O}\ket{\psi_{j}}=\sum_{j}\braket{\psi_{j}}{B^{T}}{O}\ket{\psi_{j}}=\sum_{j}\left\langle\psi_{j}\right\rvert\left(B^{T}\left\lvert
O\right\rangle\right)\ket{\psi_{j}},$ (2)

one can see by comparison with eq.1 that the reconstructed object effectively
becomes $(B^{T}\ket{O})$. The matrix $B$ can then be chosen such that it
performs any desired operation, as long as it can be expressed as a matrix
multiplication, directly during the measurement process.
To demonstrate this equivalence, we chose as our operation convolution with
the edge-detection filter kernel $K$:

$K=\begin{bmatrix}0&-1&0\\\ -1&0&1\\\ 0&1&0\end{bmatrix}.$ (3)

A 1D convolution between two discrete signals can be written in the form of a
matrix multiplication by converting one signal to the appropriate circulant (a
subclass of Toeplitz) matrix [17]. The matrix $B$ that performs the
convolution with a given kernel $K$ can be constructed by extending this
process to 2D signals with additional zero-padding and flattening steps.
Examples of the resulting illumination patterns ($\bra{\psi_{j}}B^{T}$), are
shown in Fig.1 for the cases where the initial basis $\Psi$ is either the
canonical or Hadamard basis.

The modified illumination patterns are the result of convolving the original
patterns with the filter kernel $K$, providing a more intuitive method of
generating the matrix $B$. As with any image convolution a choice of boundary
conditions is required, for which we have chosen cyclic conditions which wrap
the values at opposite edges [18].
In the following sections we demonstrate this method experimentally for the
edge-detection kernel (eq.3) in the canonical basis and compare with the
results obtained when the filtering operation is instead applied after
reconstruction.

Figure 1: Examples of illumination patterns before and after multiplication
with $B$, in the canonical basis (a) then (b) and Hadamard basis (c) then (d)
respectively. The matrix $B$ is such that it performs the edge-detection
operation of eq.3, with cyclic boundary conditions. The patterns shown are the
85${}^{\text{th}}$ in 16$\times$16 resolution bases.

##  3 Method

###  3.1 Experimental Setup

The experimental setup used to compare post-processed ghost images with those
generated with a modified illumination basis is shown in Fig.2. A 455 nm
fibre-coupled LED is collimated by plano-convex lens L1 (f = 35 mm) and
illuminates the 1080p resolution digital micromirror device (DMD). A
beamsplitter allows the LED output to be monitored by a photodiode (PD1). The
patterned light from the DMD is then imaged onto the object at reduced
magnification by plano-convex lens L2 (f = 200 mm) and biconvex lens L3 (f =
35 mm). Finally, the light transmitted by the object is collected by biconvex
lens L4 (f = 25.4 mm) and focused onto the photodiode PD2. The signal from PD2
can then be divided by the signal from PD1, compensating for fluctuations in
the LED output.

Figure 2: Schematic of the experimental setup. A fibre coupled blue LED is
collimated and illuminates a DMD. The DMD is imaged at reduced magnification
by lenses L2 and L3 onto a planar transmissive object (Obj). The transmitted
light is then collected by lens L4 and focused onto photodiode 2 (PD2). The
photodiode 1 (PD1) and beam-splitter combination allow for compensation for
fluctuations in the light source intensity.

###  3.2 Performing a Fair Comparison

One important consequence of modifying the illumination basis is that it might
increase the number of unique intensity values required to generate the new
projection patterns. This can be seen in Fig.1. To experimentally generate
such patterns with a digital micromirror device, which is only capable of
binary amplitude modulation, one needs to project and measure multiple times
per desired pattern. The number of projectable sub-patterns needed is such
that the desired pattern can be formed from a linear combination of sub-
patterns. This results in a factor of 2 increase in the total projected
patterns when comparing the canonical basis with its edge-detection
counterpart (e.g. Fig.1(a) and (b) respectively). We compensate for this in
our comparison by repeating each raster pattern twice and using their mean,
keeping the total number of projections for each method consistent. Note that
the normalisation measurement from the secondary photodiode is performed once
per pair of measurements.

###  3.3 Quantifying the Signal-to-Noise Ratio

We characterize the quality of the reconstructed images from each method and
over a range of detector integration times by means of the signal-to-noise
ratio (SNR). We calculate the image signal-to-noise ratios by comparing the
intensity values in defined peak signal ($\langle I_{P}\rangle$) and
background ($\langle I_{B}\rangle$) regions as

$\text{SNR}=\frac{\langle I_{P}\rangle-\langle I_{B}\rangle}{\sigma_{B}},$ (4)

where $\langle.\rangle$ denotes the spatial average over a region and
$\sigma_{B}$ the standard deviation of the intensity in the background region.
The background region is selected manually and is shown in green in Fig.3(b).
We define the peak signal locations as those with the top 10$\%$ of intensity
values, taken from a high SNR experimental image whilst excluding values at
the borders.

Figure 3: (a) USAF 1951 negative resolution target with the region imaged
marked by a red dashed rectangle. (b) The region highlighted in (a), with the
green dashed shape indicating the region defined as the background for later
SNR calculations.

##  4 Experimental Results

The object, a USAF target as shown in Fig.3, was imaged in the experimental
configuration shown in Fig.2 for a range of integration times. In Fig.4 we
show the main result of this paper, a comparison between the images obtained
when projecting the modified illumination basis (‘basis-processed’, left
column) and when using a raster basis which is then post-processed by
convolution with $K$ (eq.3) (‘post-processed’, right column).

Figure 4: Comparison between the ghost images obtained using a modified
projection basis (‘basis-processed’, left column) and those measured with a
raster basis and then convolved with the edge detection kernel $K$ (‘post-
processed’, right column). The three rows show varying detector integration
times, increasing from bottom to top as 20, 100 and 220 ms. The images are 64
$\times$ 64 resolution.

The experimental comparison between post-processing raster images and using a
modified illumination basis, shown in Fig.4, clearly demonstrates that the
modified patterns perform the desired spatial filtering operation. When the
integration times and thus SNRs of the images are high (Fig.4(a) and (b)), the
difference in image quality is subtle, whilst in the low SNR regime (Fig.4(e)
and (f)) the improvement offered by using a modified illumination basis
becomes quite clear.
The visual differences apparent in Fig.4 have two contributing factors. First,
the SNR is higher in the basis-processed case. This is quantified in Fig.5 as
a factor of 2 enhancement over a range of integration times. The second
difference is that the spatial character of the noise has been detrimentally
modified in the post-processing case. A post-processing filter creates
correlated noise due to the convolution theorem, whilst in the basis-
processing approach the noise remains uncorrelated (white) and is less likely
to be misinterpreted as part of the object signal.

Figure 5: Comparison of calculated SNRs from experimental images acquired
using the basis vs post-processing methods. The error bars are calculated from
the variation in three repeat measurements. Images were 64 $\times$ 64
resolution. Lines are a guide to the eye.

##  5 Noise Amplification Model

Finally we consider a simple theoretical framework for predicting the
difference in signal-to-noise ratios for the two methods and for any
convolution kernel $K$. We use additive white Gaussian noise as a model for
the detector noise in our measurement system, which we assume dominates. If an
image is corrupted by additive white noise, a given filter kernel $K$ will
increase the standard deviation of the noise by a factor of $\sqrt{E_{k}}$,
where $E_{K}$ is the energy of the filter defined as [19, 20]:

$E_{K}=\int_{\Re^{2}}|K(x,y)|^{2}\mathrm{d}x\mathrm{d}y.$ (5)

For $K$ equal to the edge-detection filter of eq.3 used in the presented
experiments, $E_{K}$ = 4 and the filter increases the standard deviation of
the noise by a factor of 2. The additive Gaussian detector noise introduced
with each of the two-step measurements of the object $O$ is denoted
$\sigma_{1}$ and $\sigma_{2}$ of equal statistics.
We represent the power normalisation step with a factor $A$, which tracks
gradual changes in the lamp intensity with time. The normalisation will
generally be imperfect due to two problems: first, the measurement of $A$ is
noisy in itself, and introduces a new detector noise $\sigma_{3}$. Second,
that either of the photodiode measurements may contain a background, e.g. due
to stray light. These are denoted $B_{m}$ and $B_{n}$ for the measurement and
normalisation signals respectively.
The signal for the post-processing method $S_{P}$, comprised of two repeat
measurements can be written

$\displaystyle S_{P}=$
$\displaystyle\,K*\left[\frac{1}{2}\left(\frac{AO+B_{m}+\sigma_{1}}{A+B_{n}+\sigma_{3}}+\frac{AO+B_{m}+\sigma_{2}}{A+B_{n}+\sigma_{3}}\right)\right]$
(6) $\displaystyle=$
$\displaystyle\,K*\left[\frac{2AO+2B_{m}+\sqrt{2}\sigma}{2A+2B_{n}+2\sigma_{3}}\right]=\,K*\left[\frac{O+B_{m}/A+\frac{\sigma}{\sqrt{2}{A}}}{1+B_{n}/A+\frac{\sigma_{3}}{A}}\right]=K*\left[\frac{O+\frac{1}{A}\left(B_{m}+\frac{\sigma}{\sqrt{2}}\right)}{1+\frac{1}{A}\left(B_{n}+\sigma_{3}\right)}\right],$

where $*$ denotes the 2D convolution and independent noise ($\sigma$) terms
can be combined in quadrature. In the case where $A$ is large compared to both
the noise and the background (i.e. there is a decent amount of signal compared
with the artefacts):

$\displaystyle S_{P}$ $\displaystyle\simeq
K*\left[\left(O+\frac{1}{A}\left(B_{m}+\frac{\sigma}{\sqrt{2}}\right)\right)-\left(O+\frac{1}{A}\left(B_{m}+\frac{\sigma}{\sqrt{2}}\right)\right)\left(\frac{1}{A}\left(B_{n}+\sigma_{3}\right)\right)\right]$
(7)
$\displaystyle=K*\left[O\left(1-\frac{1}{A}\left(B_{n}+\sigma_{3}\right)\right)+\frac{1}{A}\left(B_{m}+\frac{\sigma}{\sqrt{2}}\right)-\frac{1}{A^{2}}\left(B_{m}+\frac{\sigma}{\sqrt{2}}\right)\left(B_{n}+\sigma_{3}\right)\right]$
$\displaystyle\simeq
K*\left[O\left(1-\frac{1}{A}\left(B_{n}+\sigma_{3}\right)\right)+\frac{1}{A}\left(B_{m}+\frac{\sigma}{\sqrt{2}}\right)-\frac{1}{A^{2}}\left(B_{m}B_{n}+B_{m}\sigma_{3}+B_{n}\frac{\sigma}{\sqrt{2}}\right)\right].$

In the final step we apply the filter kernel $K$, denoting the edge-processed
version of the object $O_{E}$ and neglecting the spatial dependence of the
normalisation factor $A$, as it varies slowly compared to the size of the
filter kernel.

$S_{P}~{}\approx
O_{E}\left(1-\frac{B_{n}}{A}\right)-\frac{2O\sigma_{3}}{A}+\frac{\sqrt{2}\sigma}{A}-\frac{2B_{m}\sigma_{3}+B_{n}\sqrt{2}\sigma}{A^{2}},$
(8)

where the last term is negligible if both $B_{m}$ and $B_{n}$ are independent
from $A$. For the ‘basis-processed’ signal $S_{B}$ we can proceed similarly

$\displaystyle S_{B}$
$\displaystyle=\frac{(K_{1}*[AO+B_{m}])+\sigma_{1}-(K_{2}*[AO+B_{m}])+\sigma_{2}}{A+B_{n}+\sigma_{3}}$
(9)
$\displaystyle=\frac{K*[O+\frac{B_{m}}{A}]+\frac{\sqrt{2}\sigma}{A}}{1+\frac{B_{n}+\sigma_{3}}{A}}\simeq\left(O_{E}+\frac{\sqrt{2}\sigma}{A}\right)\left(1-\frac{B_{n}+\sigma_{3}}{A}\right)$
$\displaystyle=O_{E}\left(1-\frac{B_{n}+\sigma_{3}}{A}\right)+\frac{\sqrt{2}\sigma}{A}-\frac{B_{n}\sqrt{2}\sigma}{A^{2}},$

where $K*B_{m}$ is zero if $B_{m}$ is flat and the last term is negligible if
$B_{n}$ is independent from $A$.
Any comparison between eq.8 and eq.9 requires one to make some assumptions
about which term will be negligible and which will dominate. If the signal is
much stronger than any stray light or noise, we predictably get $S_{B}=S_{P}$.
If we keep the terms divided by $A$ (but neglect those divided by $A^{2}$), we
see that $S_{P}$ is twice as susceptible to the effect of $\sigma_{3}$ than
$S_{B}$. But in general the exact ratio between the two SNRs is very dependent
on the exact parameters of the experiment, and we expect different outcomes
depending on the specific filter kernel chosen.

##  6 Conclusions

At its very core ghost imaging is the choice of basis patterns (the
illuminations), the measurement of the coefficients (the intensity measured by
the bucket detector), and the reconstruction of the image using the two. A big
advantage of the ghost imaging method is the total freedom in the choice of
the basis used, and we have shown that the freedom of measuring the
coefficient of the expansion in a different basis from the one used for the
reconstruction allows for even more freedom. In particular one can use this
freedom to recover any linear map of the image, thus skipping the post-
processing step. We have shown this experimentally for the simple example of
an edge-detection filter.
We also studied how the signal-to-noise ratios compare between performing a
post-processing filtering, and directly measuring the processed image using
ghost imaging. Whether the ghost imaging approach results in a smaller amount
of noise depends on the linear map/filter used, and on the experimental
details (e.g. the noise in the intensity normalization), but in our experiment
we found that the ghost imaging SNR was a factor 2 better than the equivalent
post-processed image.

## Funding

This work was supported by the Leverhulme Trust’s Philip Leverhulme Prize. We
acknowledge support from the Engineering and Physical Sciences Research
Council (EPSRC) of the United Kingdom, via the EPSRC Centre for Doctoral
Training in Metamaterials (Grant No. EP/L015331/1).

## Acknowledgements

The authors acknowledge useful discussions with David B Phillips.

## Disclosures

The authors declare no conflicts of interest.

## Data Availability

Data created during this research are openly available from [21]

## References

  * [1] T. B. Pittman, Y. H. Shih, D. V. Strekalov, and A. V. Sergienko, “Optical imaging by means of two-photon quantum entanglement,” Physical Review A 52 (1995).
  * [2] R. S. Bennink, S. J. Bentley, and R. W. Boyd, “Two-Photon” Coincidence Imaging with a Classical Source,” Physical Review Letters 89, 9–12 (2002).
  * [3] A. Gatti, E. Brambilla, M. Bache, and L. A. Lugiato, “Ghost imaging with thermal light: Comparing entanglement and classical correlation,” Physical Review Letters 93, 1–4 (2004).
  * [4] J. H. Shapiro, “Computational ghost imaging,” Physical Review A - Atomic, Molecular, and Optical Physics 78, 1–4 (2008).
  * [5] Y. Bromberg, O. Katz, and Y. Silberberg, “Ghost imaging with a single detector,” Physical Review A - Atomic, Molecular, and Optical Physics 79, 1–4 (2009).
  * [6] B. I. Erkmen, “Computational ghost imaging for remote sensing,” Journal of the Optical Society of America A 29, 782 (2012).
  * [7] A. Valencia, G. Scarcelli, M. D’Angelo, and Y. Shih, “Two-photon imaging with thermal light,” Physical Review Letters 94, 1–4 (2005).
  * [8] F. Ferri, D. Magatti, A. Gatti, M. Bache, E. Brambilla, and L. A. Lugiato, “High-resolution ghost image and ghost diffraction experiments with thermal light,” Physical Review Letters 94, 2–5 (2005).
  * [9] P. Zerom, Z. Shi, M. N. O’Sullivan, K. W. C. Chan, M. Krogstad, J. H. Shapiro, and R. W. Boyd, “Thermal ghost imaging with averaged speckle patterns,” Physical Review A - Atomic, Molecular, and Optical Physics 86, 1–6 (2012).
  * [10] M. J. Padgett and R. W. Boyd, “An introduction to ghost imaging: Quantum and classical,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375 (2017).
  * [11] J. H. Shapiro and R. W. Boyd, “The physics of ghost imaging,” Quantum Information Processing 11, 949–993 (2012).
  * [12] M. A. Neifeld and P. Shankar, “Feature-specific imaging,” Applied Optics 42, 3379–3389 (2003).
  * [13] M. Liang, Y. Li, H. Meng, M. A. Neifeld, and H. Xin, “Reconfigurable array design to realize principal component analysis (PCA)-based microwave compressive sensing imaging system,” IEEE Antennas and Wireless Propagation Letters 14, 1039–1042 (2015).
  * [14] M. Aßmann and M. Bayer, “Compressive adaptive computational ghost imaging,” Scientific Reports 3, 1–5 (2013).
  * [15] D. B. Phillips, M. J. Sun, J. M. Taylor, M. P. Edgar, S. M. Barnett, G. M. Gibson, and M. J. Padgett, “Adaptive foveated single-pixel imaging with dynamic supersampling,” Science Advances 3, 1–11 (2017).
  * [16] P. del Hougne, M. F. Imani, A. V. Diebold, R. Horstmeyer, and D. R. Smith, “Learned Integrated Sensing Pipeline: Reconfigurable Metasurface Transceivers as Trainable Physical Layer in an Artificial Neural Network,” Advanced Science 7, 1–8 (2020).
  * [17] R. M. Gray, “Toeplitz and circulant matrices: A review,” Foundations and Trends in Communications and Information Theory 2, 155–239 (2006).
  * [18] W. Burger and M. J. Burge, _Digital Image Processing - An Algorithmic Introduction Using Java_ (Springer, London, 2008), 2nd ed.
  * [19] M. Jacob and M. Unser, “Design of steerable filters for feature detection using Canny-like criteria,” IEEE Transactions on Pattern Analysis and Machine Intelligence 26, 1007–1019 (2004).
  * [20] R. N. Strickland and M. Y. Aly, “Image Sharpness Enhancement Using Adaptive 3X3 Convolution Masks,” Optical Engineering 24 (1985).
  * [21] https://doi.org/10.5281/zenodo.5779444.