Chaotic Discrete Fractional-Order Food Chain Model and Hybrid Image Encryption Scheme Application

: Using the discrete fractional calculus, a novel discrete fractional-order food chain model for the case of strong pressure on preys map is proposed. Dynamical behaviors of the model involving stability analysis of its equilibrium points, bifurcation diagrams and phase portraits are investigated. It is demonstrated that the model can exhibit a variety of dynamical behaviors including stable steady states, periodic and quasiperiodic dynamics. Then, a hybrid encryption scheme based on chaotic behavior of the model along with elliptic curve key exchange scheme is proposed for colored plain images. The hybrid scheme combines the characteristics of noise-like chaotic dynamics of the map, including high sensitivity to values of parameters, with the advantages of reliable elliptic curves-based encryption systems. Security analysis assures the efﬁciency of the proposed algorithm and validates its robustness and efﬁciency against possible types of attacks.


Introduction
Fractional order differentiation and integration can be considered as a generalization of conventional integer order calculus to noninteger real or even complex valued orders. The history of fractional calculus has started about 300 years ago. Fractional calculus can be used to describe memory dependent behaviors and inherited properties of nonlinear systems. As a mathematical tool, it provides an extra degree of freedom to implement and interpret many real world systems with higher accuracy than the integerorder equivalent [1][2][3][4][5][6]. For example, most physical and engineering systems show complex behaviors such as nonlinear circuits, nanophotonics, viscous systems and laser systems. Compared to the integer-order differential equation, the memory effects are considered in the fractional-order differential equations and allow more accurate description of these natural phenomena [7]. One important aspect of dynamic systems is their chaotic behavior. A great deal of attention to this behavior has been paid in various areas of application where many fruitful results have been achieved over the past decades. For example, the authors of [8] updated Chua's model to include elements of fractional order and demonstrated chaos and other nonlinear behaviors. Several continuous time chaotic fractional order models, such as the fractional-order Lorenz [9][10][11] and the Chen fractional-order systems [12][13][14], have been examined and employed in several interesting fields such as control of chaos [15][16][17], nonlinear circuits [18,19], chaos based encryption [20] and medical applications [21]. The analytical and numerical study of dynamic propagation of light beams in the fractional order Schrödinger equation with a harmonic potential have been presented in [22]. The dynamics in stochastic models of Gaussian-amplitude field, phase-diffusion and chaotic field in laser have been investigated in [23,24].

Preliminaries
In order to overcome the difficulties that arise from dealing with continuous time fractional order system and efficiently capture the memory effects, discrete fractional calculus was introduced [46][47][48][49][50]. The studies of dynamic behaviors and applications of fractional delta difference models attracted increasing interest in the last decade, see [49][50][51][52][53][54] and references therein.

Definition 3.
The order α fractional delta difference equation is defined as [55]: where its associated discrete time integral is expressed as Here, the initial value can be written as Theorem 1 ([55]). The delta fractional difference equation has the following equivalent discrete integral equation

Theorem 2.
The conditions of asymptotic stability of zero equilibrium point to the next discrete fractional-order system for every eigenvalues λ of the discrete fractional system.

Discrete-Time Fractional-Order Food Chain Model
Discrete-time dynamic systems are ideal to describe the population dynamics of nonoverlapping organisms. It is well known that one of the basic population models is the Lotka-Volterra prey-predator model. Holling has introduced the study of more practical predator models since the groundbreaking theoretical works by Lotka [56] and Volterra [57] in the last century, proposing three types of functional responses for different species to model predation dynamics [58].
The discrete time models are known to exhibit more complicated dynamics than their counterpart continuous time models. For example, the chaotic behavior can be induced by the simple one-dimensional logistic map while only exponential growth or decaying of state variable can be observed in logistic differential equation. The main interest of mathematicians and scientists was focusing on continuous time models in mathematical biology and ecology for a long time. The capability of discrete time models to describe some cases of species dynamics efficiently brought them to light quite recently. In 2020, a novel discrete time novel food chain model was introduced [59] as the first discrete model to consider three interacting organisms, having non-overlapping generations, which is subjected to intraspecific competition and strong pressure on preys species. More specifically, we consider a prey population x predated by a first predator species y. The third species is the top predator z that predates on the first predator y and simultaneously interferes with the population growth of prey. The model proposed for studying these ecological interactions is represented by the following nonlinear discrete system [59]: For the function f (n), the delta difference operator is defined by ∆ f (n) = f (n + 1) − f (n).
Based on previous assumptions, we get the following discrete fractional-order food chain model for (1): The fixed points of system (2) satisfy By simple algebraic computations, we obtain four fixed points for the abo E 0 = (0, 0, 0), In the following subsection, the local stability of these fixed points are studied for model (2).

Stability of Fixed Points
The local stability analysis of the fixed points is established by studying the Jacobian matrix of system (1) at these fixed points. The system (2) can be linearized about any fixed point (x * , y * , z * ) via computing its associated Jacobian matrix which takes the form The characteristic equation of matrix J is computed and expressed in the form 2ab(x * ) 2 + abx * + aby * z * + ab(z * ) 2 − abz * − 2acx * y * − acy * z * − ac(y * ) 2 + acy * + 2ax * + ay * + az * − a + bcx * y * − bx * + bz * − cy * + 1. Thus, the stability analysis of each fixed point is carried out as follows: 3.1.1. Stability Analysis of E 0 For this fixed point, it is found that the eigenvalues of J are given by This implies that Referring to the conditions of asymptotic stability of fixed points which are given in previous section, the fixed point E 0 is locally asymptotically stable if 0 < a < 1. Figure 1 shows the stability region of fixed point E 0 in parameters' plane (a, α).

Stability Analysis of E 1
For this fixed point, it is found that the eigenvalues of J are given by that implies that Referring to the conditions of asymptotic stability of model's steady states which are given in the previous section, the fixed point E 1 satisfies local asymptotic stability in the red colored regions illustrated in Figure 2 which depends on the values of a, b and α.

Stability Analysis of E 2
For this fixed point, it is found that the eigenvalues of J are given by In order to find local asymptotic stability region in space of parameters of model (2), numerical evaluations of values of parameters which satisfy stability conditions can be used. In particular, Figure 3 shows examples of stability regions in (a, b) and (a, c) planes of parameters for distinct values of fractional order α.

Stability Analysis of E 3
First, this point takes values within the feasible space if the following condition is satisfied which means that this condition should be examined along with the aforementioned stability conditions. For the fixed point E 3 , it is found that the eigenvalues of J have very complicated expressions which renders the numerical investigations of asymptotic stability region in space of parameters of model (2). In Figure 4, regions of occurrence of fixed point E 3 in addition to colored stability regions in (a, b) and (a, c) planes of parameters at distinct values of α are depicted.

Numerical Simulations
In this section, results of numerical simulations are shown for integer order model (1) and fractional order model (2). The phase portraits and bifurcation diagrams are employed to determine variation in dynamical behaviors of the two models in terms of variations in parameters in the models. In the following simulations, system (2) is used in the following form: Firstly, the aforementioned conditions of stability of fixed points of model (2) are verified. In particular, the values of parameters which yield a stable fixed point E 0 is from Figure 1 and the associated time series outputs of model (2) are presented in Figure 5 and confirm theoretical results. Similarly, the values of parameters which correspond to stable fixed points E 1 , E 2 and E 3 are presented from Figures 2-4, respectively. Furthermore, the output time series which are illustrated in Figures 6-8 are, respectively, verify stability conditions of fixed points E 1 , E 2 and E 3 . In Figure 6, time series of state variables of (2) show stability of fixed point E 1 at a = 2, b = 1, c = 1 and α = 0.9. Secondly, the bifurcation diagrams are employed to overview the influences of parameter's variations on dynamical behaviors of models (1) and (2). The case of integer order model (1) is presented in Figure 9 where the effects of variations in parameters a and b are shown. The influences of fractional order α along with other parameters are depicted in Figure 10. Finally, some selected examples of phase portraits of the two models (1) and (2) are shown in Figure 11.
It is obvious that increasing predation rate of predators on the prey, the mode can exhibit chaotic dynamics. In particular, the model undergoes a stable Neimarck-Sacker bifurcation followed by period-doubling bifurcations till chaos behavior starts arises. Figure 11. Examples of the two and three dimensional phase portraits of the models (1) and (2)

Hybrid Image Encryption Scheme
In this section, we propose an encryption scheme that combines elliptic curve key exchange technique with chaotic output of a three-dimensional mapping. Numerical simulations on different color images are used to validate the efficiency of the scheme against differential, statistical in addition to brute-force attacks.
Input: Assume that a colored plain image with size of m × n pixels is given. Then, three values of color components are associated to each pixel in the image. These color components are red, green and blue such that for each pixel with position (i, j), let R(i, j), G(i, j) and B(i, j) refer to the values of red, green and blue color components, respectively. Typi-cally, they take the range [0 : 255]. In addition, assume that the reading of internal clock of encrypting machine is given by T * at the moment when encryption session starts.
Public Keys: Pick one family of elliptic curves standardized by the NIST, its associated group generator and parameters are considered as public keys of the suggested scheme. In the next numerical experiments, we adopt the P-192 curve groups of the following form (2) Update the values of secret parameters according to one of the following rules (3) Apply the updated values of parameters to simulate the model (2) discarding any transient nonchaotic dynamics. The resulting chaotic time series of lengths m × n + 2 × max{m, n} are to be used in the next steps.
(4) Construct the following encrypting sequences The values of S r R , S r G , S r B , S c R , S c G and S c B are arranged in an ascending order to formulate six confusion vectors. Hence, the rows of original matrix of pixels are permuted such that the red components are scrambled according to components of S r R whereas green and blue values in each row are confused by S r G and S r B orders, respectively. By the same way, the columns in the original image are scrambled by utilizing S c R , S c G and S c B . (6) The plain image is reshaped into three vectors each of which has m × n values. These vectors involve the separate values of pixels' color intensity. They are referred as I 1 , I 2 and I 3 .
(7) The bitwise XOR operations are carried out between vectors S 1 , S 2 , S 3 and I 1 , I 2 , I 3 such that the three encrypted components of cipher images are computed as (8) The elliptic curve key exchange technique in the sense of Diffie-Hellman is adopted such that the sender side publishes P t G whereas the receiver side publishes P r G. Thus, the sender and receiver will agree on a common symmetric key P t P r G = P r P t G.
(9) Three perturbation values δ t , t = R, G, B are encoded using the shared secret key. More specifically, El-Gamal scheme can be applied [60] at this step.
(10) The shared secret keys and similar numerical precision settings at both sides, imply that identical versions of chaotic sequence S 1 , S 2 and S 3 are generated at the receiver part. (11) The transmitted ciphered image is deciphered through the aforementioned bitwise XOR operations.
(12) Finally, the deciphered vectors are reshaped to restore the original plain image.

Numerical Simulations
Numerical simulations are performed at a = 2, b = 3.35, c = 9.15, µ k = 0.1250075381 + 10 −3 k and α = 0.985. Figure 12a shows the plain King Tut image, ciphered King Tut image and deciphered King Tut image. The histograms for separate colors components in the pixels of these images are given in Figure 12b. In addition, Figures 13-15 show the results correspond to images of Baboon, pepper and Egyptian pyramids, respectively. It is clear that the distribution of pixels' values in encrypted images is almost flat and uniform for each color which renders the encrypted images invulnerable to different statistical attacks. The uniformity of pixels distribution in ciphered images is quantified using the variance of histogram concept. More specifically, the small values of variances indicate high level uniformity. It can be defined for red, green and blue colors as follows: where h R i , h G i and h B i are the number of pixels with the value of i for red, green and blue values of pixels, respectively. The results obtained for each image are summarized in Table 1. The results show the considerable reduction in variances of histograms for cipher images and thus confirm the uniformity of histogram values.

Keyspace Analysis
The proposed hybrid encryption technique has three initial conditions, four parameters of the model, i.e., a, b, c and α, three image-depending parameters µ i s. Suppose that the IEEE 754 floating-point format is used. Therefore, the secret keys of the proposed scheme will have space size equals 2 530 . This does not include the parameters of elliptic curve. Note that 2 100 space size is the minimum necessary size in order to break brute-force attacks [58,59]. Therefore, our encryption system has a big enough keyspace to make brute force attacks useless.

Analysis of Key Sensitivity
Another essential requirement for any reliable encryption scheme is to have high sensitivity to tiny perturbations in the values of secret keys. To investigate this point, the value of each secret key is perturbed by 10 −14 at time and then the generated chaotic time series are used to decrypt the encrypted images. Then, the sensitivity to mismatch in parameters is investigated. Figure 16 shows an example when the value of b is perturbed by 10 −14 . The deciphered images are shown in the figure. It is clear that when very small differences in b occur, we successfully decrypt any encrypted image. Numerical simulations illustrate that similar conclusions are inferred about the remaining secret keys in the system.
Finally, in order to measure the sensitivity to mismatches in parameters, Table 2 shows the original values of secret keys employed in encryption process, the mismatch of secret keys during decryption process, and the average difference percentage between the correct and incorrect decrypted images.

Analysis of Information Entropy
In order to investigate the degree of randomness and uncertainty of encrypted images, the information entropy is utilized. The large values of information entropy indicates high randomness of encrypted information. The information entropy has units in bits and it defined for an input information as follows [60] where s i denotes specific form i of input data, N is the number of possible different forms of input data and p(s i ) refers to the probability of data symbol s i . Generally, the optimal value of H(s) of encrypted images is to be near to the value of 8. Table 3 illustrates the values of H(s) in the produced encrypted images. The computed values are almost equal to eight which shows the efficiency of the suggested technique.

Differential Attacks Analysis
For image encryption scheme to be efficient, it should be also very sensitive to teeny variations in plain images along with secret keys of encryption process. This shows that the very small perturbations that added to input data should make significant alternations in the resulting cipher data and hence the chaos based encryption scheme become more immune to different differential attacks [61][62][63].
For quantifying sensitivity to teeny changes in plain images, two measures can be successfully employed. The first of them is the Number of Pixels Change Rate (NPCR) that measures the percentage of unlike pixels between two encrypted images for one pixel difference in their associated plain images. The second measure is the Unified Average Changing Intensity (UACI), which represents the mean of intensity differences between two encrypted images for a single pixel value difference in the two input plain images.
Assume that C 1 and C 2 are two cipher images whose plain images are identical except for one pixel only. Let C 1 (i, j) and C 2 (i, j) refer to the specific color components of the pixel at position (i, j) in each of the two images. Thus, the NPCR is defined by [64,65] The second measure is given by Table 4 presents the values of NPCR and UACI for difference in one color component in random single pixel of the two original images.

Resistance against Other Attacks
The Kerckhoff's principle [66] states that through evaluation of security strength of a given encryption system, it should be assumed that the attacker has knowledge about the design and detailed structure of encryption system. The exception is the unknown values of secret keys. Taking this principle into account, there are four possible types of attacks, known as, known plaintext, ciphertext only, chosen plaintext (CP) and chosen ciphertext (CC) attacks, which can be utilized against the proposed system. Most importantly, the opponent in chosen ciphertext attack (CCA) can get access to the decryption system while in chosen plaintext attack (CPA), he secretly can establish access to the encryption system itself. Now, if the suggested cryptography system is immune to the aforementioned CCA and CCP, then the resistance against the other two types are confirmed [67,68].
The presented encryption scheme involves setting up the values of time-dependent and plain image-dependent parameters, scrambling of pixels, bit XOR of pixel values perturbing chaotic signal and elliptic curve key exchange as follows. Firstly, some key features are read out from the plain image to modify the original values of secret keys of the system. This indicates that different plain images will induce different cipher images even for the cases where input images having very small differences. Moreover, feeding the encryption machine with the same plain image, but at distinct times, will produce separated encrypted images because secret keys are time dependent on the moments that images are supplied to the system. In addition, assume that the attackers try to defeat the elliptic curve key exchange step to unveil the values of secret keys via applying the well known Pollard's Rho attack or Baby Step, Giant Step attack [69][70][71][72]. The attacker will find that it is practically infeasible to achieve his goal. The interpretation is that although this attacks can obviously reduce the computations complexity of discrete logarithmic problem, it still requires approximately √ P EC = 7.9 × 10 28 of operations so as to break the encryption system. More specifically, it approximately requires more than 10 13 years to complete the attacking process utilizing 16 GB RAM and Intel Core i7-8550U CPU. Subsequently, the present hybrid technique is reliable against the known-plaintext and chosen-plaintext attacks [67,68] in addition to EC attacks.

Discussion and Conclusions
Analytical and numerical frameworks are presented to analyze a proposed discrete fractional-order food chain model. The present model exhibits rich dynamics and a variety of nonlinear phenomena is exhibited by its state variables.
It is the first model to consider discrete fractional order three-dimensional food chain model while the other models in literature are exclusively continuous time or integer-order discrete time models. The fractional order which represents the memory influences is found to have a significant impact on the stability of nontrivial fixed points of food chain model and it should not be ignored. More specifically, when the value of fractional order decreases from one, the stability of fixed point may be changed. For example, the fixed point E 1 is stable at a = 2.9 and b = 1 when fractional order is very close to one and it loses its stability and becomes unstable when α = 0.8. A second example is the fixed point E 2 which is stable at a = 3.05 and b = 2.4, i.e., when the predator species has increased value of its benefit from predation on preys while simultaneously the fractional order value is very close to one. The fixed point E 2 losses its stability when α is decreased to 0.8. A third example is the coexistence fixed point E 3 which is stable at α > 0.88, a = 3.05, b = 3 and c = 4.35 while it losses its stability when α is decreased. Figures 1-4 depict stability regions of fixed points in the space of parameters in more details. From biological point of view, the above discussion indicates that it possible to find two separate food chain systems which have identical values of parameters, i.e., same species for prey, predators and top-predators in the two systems, whilst the two systems undergo different stable equilibrium points. This can be therefore interpreted by referring to memory impacts or the different values of fractional order differences possessed by the two systems.
Employing multi-dimensional chaotic discrete fractional difference equations in encryption application is a very recent approach in cryptography systems. This approach overcomes numerical errors of continuous fractional order equations, has the advantages of extending the key space size by fractional order secret key and it is more appropriate for being implemented on digital platforms than the continuous time fractional systems. In addition, combining elliptic curve scheme for key exchange with the chaos source of three-dimensional fractional mapping in a novel hybrid encryption system, for first time, inherits the reliability and efficiency of both systems. From security point of view, the secret keys of the hybrid encryption system are made simultaneously time-varying and plain data-dependent that render the system capable to defeat the powerful CCA and CPA attacks along with other attacks, as discussed above. The realization on digital appliances can be conducted via using Arduino boards, such as ARDUINO NANO 33 BLE or ARDUINO MEGA 2560 REV3, Digital Signal Processors (DSPs), such as TMS320C6452 or ADSP-2136x, field-programmable gate array (FPGAs) or microcontrollers, such as nRF52840 or ATmega640. Future work may include adopting mixed functional response for the model to better reflect biological factors. In addition, the more advanced supersingular isogeny elliptic curves can be employed to further improve security strength of the proposed hybrid encryption scheme against risks of post quantum computers era.