Research article

The security analysis of the key exchange protocol based on the matrix power function defined over a family of non-commuting groups

  • In this paper, we revisited the previously proposed key exchange protocol based on the matrix power function. We prove that the entries of the public key matrices of both parties of the protocol are uniform. Using this result we defined a security game for our protocol and show that the malicious attacker cannot gain any significant advantage in winning this game by applying faithful representation or the linearization approaches. Moreover, we showed that the shared key is computationally indistinguishable from the imitation key if the security parameters are properly chosen.

    Citation: Aleksejus Mihalkovich, Jokubas Zitkevicius, Eligijus Sakalauskas. The security analysis of the key exchange protocol based on the matrix power function defined over a family of non-commuting groups[J]. AIMS Mathematics, 2024, 9(10): 26961-26982. doi: 10.3934/math.20241312

    Related Papers:

    [1] Xiaoming Su, Jiahui Wang, Adiya Bao . Stability analysis and chaos control in a discrete predator-prey system with Allee effect, fear effect, and refuge. AIMS Mathematics, 2024, 9(5): 13462-13491. doi: 10.3934/math.2024656
    [2] Xiongxiong Du, Xiaoling Han, Ceyu Lei . Dynamics of a nonlinear discrete predator-prey system with fear effect. AIMS Mathematics, 2023, 8(10): 23953-23973. doi: 10.3934/math.20231221
    [3] A. Q. Khan, Ibraheem M. Alsulami, S. K. A. Hamdani . Controlling the chaos and bifurcations of a discrete prey-predator model. AIMS Mathematics, 2024, 9(1): 1783-1818. doi: 10.3934/math.2024087
    [4] Ahmad Suleman, Rizwan Ahmed, Fehaid Salem Alshammari, Nehad Ali Shah . Dynamic complexity of a slow-fast predator-prey model with herd behavior. AIMS Mathematics, 2023, 8(10): 24446-24472. doi: 10.3934/math.20231247
    [5] A. Q. Khan, Ibraheem M. Alsulami . Discrete Leslie's model with bifurcations and control. AIMS Mathematics, 2023, 8(10): 22483-22506. doi: 10.3934/math.20231146
    [6] Saud Fahad Aldosary, Rizwan Ahmed . Stability and bifurcation analysis of a discrete Leslie predator-prey system via piecewise constant argument method. AIMS Mathematics, 2024, 9(2): 4684-4706. doi: 10.3934/math.2024226
    [7] Jie Liu, Qinglong Wang, Xuyang Cao, Ting Yu . Bifurcation and optimal harvesting analysis of a discrete-time predator–prey model with fear and prey refuge effects. AIMS Mathematics, 2024, 9(10): 26283-26306. doi: 10.3934/math.20241281
    [8] Wei Li, Qingkai Xu, Xingjian Wang, Chunrui Zhang . Dynamics analysis of spatiotemporal discrete predator-prey model based on coupled map lattices. AIMS Mathematics, 2025, 10(1): 1248-1299. doi: 10.3934/math.2025059
    [9] Ali Al Khabyah, Rizwan Ahmed, Muhammad Saeed Akram, Shehraz Akhtar . Stability, bifurcation, and chaos control in a discrete predator-prey model with strong Allee effect. AIMS Mathematics, 2023, 8(4): 8060-8081. doi: 10.3934/math.2023408
    [10] A. Q. Khan, Ibraheem M. Alsulami . Complicate dynamical analysis of a discrete predator-prey model with a prey refuge. AIMS Mathematics, 2023, 8(7): 15035-15057. doi: 10.3934/math.2023768
  • In this paper, we revisited the previously proposed key exchange protocol based on the matrix power function. We prove that the entries of the public key matrices of both parties of the protocol are uniform. Using this result we defined a security game for our protocol and show that the malicious attacker cannot gain any significant advantage in winning this game by applying faithful representation or the linearization approaches. Moreover, we showed that the shared key is computationally indistinguishable from the imitation key if the security parameters are properly chosen.



    Predator-prey interaction is one of the most important subjects in the bio-mathematical literature. Therefore, many ecologists, mathematicians and biologists have investigated the dynamical behavior of the predator-prey system describing the interaction between prey and predator. The Lotka-Volterra predator-prey system, which is a fundamental population model, was introduced by Lotka [1] and Voltera [2]. The Lotka-Volterra model assumes that "the prey consumption rate by a predator is directly proportional to the prey abundance. This means that predator feeding is limited only by the amount of prey in the environment. While this may be realistic at low prey densities, it is certainly an unrealistic assumption at high prey densities where predators are limited, for example, by time and digestive constraints [3]." Since this model has neglected many real situations and complexities, it has been modified over the years by numerous researchers to ensure an accurate explanation and better understanding [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]. Additionally, many ecological concepts, such as diffusion, functional responses, time delays and the Allee effect, have been added to the predator-prey system to facilitate a more accurate description [7,13,14,15,16,20,23,24,25,26,27,28]. It is known to all that one of the important modifications to the system is immigration. The concept of immigration, which affects the population growth rate, is an external factor expressed for an organism that establishes a particular habitat. Both immigration and migration happen when there is shortage of basic needs due to human activities. Especially, predator-prey interactions may be disrupted due to different reactions to climate change. Sometimes, regular and seasonal migration and immigration drive living organisms to adopt traveling and nomadic lives [29]. Additionally, the process of immigration is a significant phenomenon that helps the ecosystem to achieve a state of stability. Specially, discrete-time systems that take into account seasonal events, such as immigration are more realistic [6,11,12,26,27,30,31,32,33,34,35]. There have also been many studies on the role of immigration and its impact on population dynamics. Sugie et al. studied the existence and uniqueness of limit cycles in a predator-prey system with constant immigration in [26]. Zhu et al. investigated the local and global stability in a delayed predator-prey system with constant-rate immigration [27]. Also, they showed the existence of the global Hopf bifurcation. Tahara et al. analyzed the asymptotic stability of the predator-prey systems by adding an immigration factor to the prey and predator populations in a classical Lotka-Volterra system [36].

    Rana [23] investigated the dynamics of a discrete-time predator-prey system of Holling I type. The author analyzed the existence and local stability of positive fixed points in the system and showed that the system undergoes flip and Neimark-Sacker bifurcations. Also, Rana observed that, when the prey is in a chaotic dynamic state, the predator can tend to extinction or to a stable equilibrium. It is important to consider the effects of the presence of some number of immigrants because most predator-prey systems in nature are not isolated. The immigration effect of the predator-prey system has been rarely studied in literature. Therefore, we will investigate a discrete-time predator-prey system with constant immigration of the prey population, as follows:

    Xt+1=rXt(1Xt)aXtYt+sYt+1=bXtYtdYt, (1.1)

    where Xt and Yt represent the population densities of the prey and predator at time t, respectively. r is the intrinsic growth rate of the prey, a is the per capita searching efficiency of the predator, b is the conversion rate of the predator, d is the death rate of the predator and s is the prey immigration rate. All parameters, r,a,b,d and s, are positive constants. The predator-prey system given by (1.1) assumes that the prey grows logistically with an intrinsic growth rate r and a carrying capacity of one in the absence of predation.

    The contributions of this study are given as follows:

    (1) The intended system consists of two interacting species, where one species is the food source for the other. In this paper, we have analyzed the immigration effect on the prey population in the system.

    (2) The stability of the considered system is analyzed for possible fixed points.

    (3) It has been shown that the proposed system undergoes both flip and Neimark-Sacker bifurcations.

    (4) The OGY method has been applied to the system to control the chaos due to the emergence of the Neimark-Sacker bifurcation.

    (5) Some numerical examples have been presented for our discrete-time predator-prey system with immigration in order to support the accuracy of our theoretical results.

    Our task is to discuss the dynamics of a modified discrete-time predator-prey system with immigration of the prey population. The rest of the paper is organized as follows. In Section 2, the existence conditions and stability of the fixed points of the system are investigated. In Section 3, Neimark-Sacker and flip bifurcation analysis is explored by choosing r as a bifurcation parameter. Next, the directions of both Neimark-Sacker and flip bifurcations are obtained by using normal form theory [21,37]. In Section 4, an OGY feedback control strategy is implemented to control the chaos due to the emergence of Neimark-Sacker bifurcation. In Section 5, some numerical simulations are demonstrated out to support the accuracy of our theoretical finding. In the end, a brief conclusion is given.

    In this section, we will examine the existence and stability of the fixed points of a discrete system with immigration of the prey in the close first quadrant R2+. To find the fixed points of the system given by (1.1), we can rewrite system (1.1) in the following form:

    X=rX(1X)aXY+sY=bXYdY. (2.1)

    Solving the algebraic equation, we get that system (1.1) has the following fixed points:

    E1=((r1)+(r1)2+4rs2r,0):SurvivalofPopulation Xonly,andE2=(d+1b,(d+1)b(r1)r(d+1)2+sb2ab(d+1)):Coexistenceofallpopulations.

    We have obtained the following lemma regarding the existence of the fixed points of system (1.1).

    Lemma 2.1. For system (1.1), the following statements hold true:

    i) System (1.1) always has an axial positive fixed point E1=((r1)+(r1)2+4rs2r,0).

    ii) System (1.1) has a unique positive fixed point E2=(d+1b,(d+1)b(r1)r(d+1)2+sb2ab(d+1)) if (d+1)b(r1)+sb2>r(d+1)2.

    Now, we investigate the stability of the coexistence fixed point E2 of system (1.1) only. The Jacobian matrix of system (1.1) evaluated at the unique positive fixed point E2=(d+1b,(d+1)b(r1)r(d+1)2+sb2ab(d+1)) is given by

    J(E2)=((d+1)(br(d+1))sb2b(1+d)a(d+1)b(d+1)b(r1)r(d+1)2+sb2a(d+1)1).

    Then, the characteristic polynomial of J(E2) is given by

    F(λ)=λ2+((d+1)2r2b(d+1)+sb2(d+1)b)λ+((d+1)d(d2+4d+5)2(d+1)b))r+sbd(d+1)d. (2.2)

    In order to analyze the dynamics of the unique positive fixed point E2=(d+1b,(d+1)b(r1)r(d+1)2+sb2ab(d+1)) of system (1.1), we present Lemma 2.2 [8,18,19]:

    Lemma 2.2. Let F(λ)=λ2+Bλ+C, where B and C are two real constants, and let F(1)>0. Suppose λ1 and λ2 are two roots of F(λ)=0. Then, the following statements hold true:

    i) |λ1|<1 and |λ2|<1 if and only if F(1)>0 and C<1,

    ii) |λ1|>1 and |λ2|>1 if and only if F(1)>0 and C>1,

    iii) ( |λ1|<1 and |λ2|>1) or ( |λ1|>1 and |λ2|<1) if and only if F(1)<0,

    iv) λ1 and λ2 are a pair of conjugate complex roots, and |λ1|=|λ2|=1 if and only if B24C<0 and C=1,

    v) λ1=1 and |λ2|1 if and only if F(1)=0 and C±1.

    Assume that λ1 and λ2 are the roots of the characteristic polynomial at the positive fixed point (x,y). Then, the point (x,y) is called a sink if |λ1|<1 and |λ2|<1 and it is locally asymptotically stable. (x,y) is known as a source or repeller if |λ1|>1 and |λ2|>1; thus, a source is always unstable. The point (x,y) is called a saddle point if (|λ1|<1 and |λ2|>1) or (|λ1|>1 and |λ2|<1). And, (x,y) is called non-hyperbolic if either |λ1|=1 or |λ2|=1.

    Now, we will discuss the topological classification of the unique positive fixed point (d+1b,(d+1)b(r1)r(d+1)2+sb2ab(d+1)) of system (1.1), and we will apply Lemma 2.2 to prove the following lemma. From Eq (2.2), we can obtain  F(1)=(d+1)b(r1)r(d+1)2+sb2b>0 if (d+1)b(r1)+sb2>r(d+1)2.

    Lemma 2.3. Assume that (d+1)b(r1)+sb2>r(d+1)2; then, for a unique positive fixed point E2 of system (1.1), the following holds true.

    i) E2 is a sink-type fixed point if the following conditions hold:

    (1+d)(bd+(2+bd)(1+d)r)+b2ds>0

    and

    (1+d)(bd+(2+bd)(1+d)r)+b2ds<b(1+d).

    ii) E2 is a source (repeller)-type fixed point if the following conditions hold:

    (1+d)(bd+(2+bd)(1+d)r)+b2ds>0

    and

    (1+d)(bd+(2+bd)(1+d)r)+b2ds>b(1+d).

    iii) E2 is a saddle-type fixed point if the following condition holds:

    (1+d)(bd+(2+bd)(1+d)r)+b2ds<0.

    iv) Assume that λ1 and λ2 are the roots of Eq (2.2); then λ1=1 and |λ2|1 if and only if r=r1 and rr2,r2+2bd2bd+3db+2.

    v) The roots of Eq (2.2) are complex with modules of one if and only if

    r=r2 and sb2(d+1)2<r2<(4(d+1)sb)b(d+1)2,

    where r1=(2d+3d2sb(1d))b(d+1)2(bd3) and r2=((d+1)2sbd)b(d+1)2(b(d+2)).

    Example 2.1. For the parameter values a = 3.5, b = 4.5, d = 0.25, s = 0.75, r = 0.5 and the initial condition (x0,y0)=(0.25,0.6), the positive unique fixed point of system (1.1) is obtained as (x,y)=(0.27777778,0.58888889). Figure 1 shows that the fixed point (x,y) of system (1.1) is locally asymptotically stable, with x(t) and y(t) being the prey and predator populations at time t, respectively.

    Figure 1.  Stable fixed points for system (1.1) for a = 3.5, b = 4.5, d = 0.25, s = 0.75, r = 0.5 and the initial condition (x0,y0)=(0.25,0.6).

    In this section, we will investigate flip and Neimark-Sacker bifurcations of the coexistence fixed point E2 by using the bifurcation theory in the sense of [15].

    In this subsection, we investigate the existence conditions for the Neimark-Sacker bifurcation around coexistence fixed point by choosing the parameter r as a bifurcation parameter. Furthermore, the direction of the Neimark-Sacker bifurcation is given by applying bifurcation theory.

    When the term (v) of Lemma 2.3 holds, we say that two eigenvalues of J(E2) are a pair of complex conjugates with a modulus of one. Therefore, we can write the following set:

    NSBE2 =

    {(r,a,b,d,s)R5+:(d+1)b(r1)+sb2>r(d+1)2, rNS=r2, sb2(d+1)2<r2<(4(d+1)sb)b(d+1)2}, (3.1)

    where r2=((d+1)2sbd)b(d+1)2(b(d+2)).

    The Neimark-Sacker bifurcation appears when the parameter r varies in a small neigborhood of NSBE2. The eigenvalues of system (1.1) under these conditions are obtained by

    λ,¯λ=(MT+bK)±iQ2bK, (3.2)

    where

    M=sb2bdb,T=rd2+2rd+r,K=d+1,Q=(M+T)2+2bK(M+T)+b2K2(4Kr1)+4bK2(MT).

    It is easy to see that

    |λ|=|¯λ|=1. (3.3)

    From the transversality condition, we get

    d|λi(r)|dr|r=rNS=Kd2(d+4)(5d2)Kb0, i=1,2. (3.4)

    If the nonresonance condition B=trJE2(rNS)0,1, then

    r2bα,b(α+1K), (3.5)

    which obviously satisfies

    λk(rNS)1 for k=1,2,3,4, (3.6)

    where α=2KsbK2. Assume that q,pC2 are two eigenvectors of J(NSBE2) and the transposed matrix JT(NSBE2) corresponding to λ and ¯λ, respectively. We have

    q(1,12K2a(M+T+bK+iQ))T (3.7)

    and

    p(12K2a(M+T+bK+iQ),1)T. (3.8)

    To achieve the normalization <p,q>=1, where <.,.> means the standard scalar product in C2, we can set the normalized vectors as

    q=(1,12K2a(M+T+bK+iQ))Tp=(12i(M+T+bK)2Q,iK2aQ)T. (3.9)

    From Taylor expansion, substituting xt+x into Xt, and yt+y into Yt, the fixed point E2 of system (1.1) results in the origin (0,0), and system (1.1) converts to

    (xtyt)J(E)(xtyt)+(F1(xt,yt)F2(xt,yt)), (3.10)

    where J(E)=J(x,y), and

    F1(xt,yt)=rx2taytxt+O(X4t),F2(xt,yt)=bxtyt+O(X4t),

    where Xt=(xt,yt)T. System (3.10) can be expressed as

    (xt+1yt+1)=J(E)(xtyt)+12B(xt,xt)+16C(xt,xt,xt)+O(x4t), (3.11)

    where B(x,y)=(B1(x,y)B2(x,y)) and C(x,y,u)=(C1(x,y,u)C2(x,y,u)) are symmetric multi-linear vector functions of x,y,uR2. We define the functions B(x,y) and C(x,y,u) by the following formulas:

    B1(x,y)=2j,k=12F1ξjξk|ξ=0xjyk=2rx1y1ax1y2ax2y1,B2(x,y)=2j,k=12F2ξjξk|ξ=0xjyk=bx1y2+bx2y1,C1(x,y,u)2=j,k,l=13F1ξjξkξl|ξ=0xjykul=0,C2(x,y,u)=2j,k,l=13F2ξjξkξl|ξ=0xjykul=0. (3.12)

    For r near to rNS and zC, we use the transformation X=zq+¯z¯q; then, system (3.11) becomes

    zλ(r)z+g(z,¯z,r), (3.13)

    where λ(r)=(1+φ(r))eiarctan(r), with φ(rNS)=0, and g(z,¯z,r) is a smooth complex-valued function. After Taylor expression of g with respect to (z,¯z), we obtain

    g(z,¯z,r)=k+l21k!l!gkl(r)zk¯zl, with gklk,l=0,1,.... (3.14)

    Given the symmetric multi-linear vector functions, the Taylor coefficients gkl can be expressed by the following formulas:

    g20(rNS)=<p,B(q,q)>,g11(rNS)=<p,B(q,¯q)>,g02(rNS)=<p,B(¯q,¯q)>,g21(rNS)=<p,C(q,q,¯q)>. (3.15)

    The coefficient β2(rNS), which determines the direction of the appearance of the invariant curve in a generic system exhibiting the Neimark-Sacker bifurcation, can be calculated as follows:

    β2(rNS)=Re(eiarctan(rNS)2g21)Re((12eiarctan(rNS))e2iarctan(rNS)2(1eiarctan(rNS))g20g11)12|g11|214|g02|2, (3.16)

    where eiarctan(rNS)=λ(rNS).

    From the above analysis, we have the following theorem:

    Theorem 3.1.1. Suppose that E2 is a positive unique fixed point of system (1.1). If Eq (3.5) holds, β2(rNS)0 and the parameter r changes its value in a small vicinity of NSBE2; then, system (1.1) passes through a Neimark-Sacker bifurcation at the fixed point E2. Moreover, if β2(rNS)<0 (β2(rNS)>0), then the Neimark-Sacker bifurcation of System (1.1) at r=rNS is supercritical (subcritical) and there exists a unique closed invariant curve bifurcation from E2 for r=rNS, which is attracting (repelling).

    In this subsection, we will discuss the existence conditions for and direction of the flip bifurcation at a unique positive fixed point E2 of system (1.1) by choosing r as a bifurcation parameter.

    If Lemma 2.3 (iv) holds, we see that one of the eigenvalues of J(E2) is 1 and the other eigenvalue is neither 1 nor 1. Therefore, there can be a flip bifurcation of the fixed point E2 if the parameter r varies in the small neighborhood of FBE2, where

    FBE2={(r,a,b,d,s)R5+:(d+1)b(r1)+sb2>r(d+1)2, r=rF=r1 , rr2,r2+2bd2bd+3db+2}. (3.17)

    The eigenvalues of system (1.1) under these conditions are given by

    λ1(r1)=1,λ2(r1)=b2s+4d22bs3b2bds+10d+63bdd24d+bd+b3. (3.18)

    In order to be λ2(r1)±1, we have

    bs(b2(1+d))3d2+2b6d+2bd3,5d2+4b14d+4bd9. (3.19)

    Suppose that q,pR2 are two eigenvectors of J(FBE2) and the transposed matrix JT(FBE2), respectively, for λ1(r1)=1. Then, we have

    J(FBE2)q=q and JT(FBE2)=p.

    By direct calculation, we obtain

    q(1,(3d2+6d2bd2bds+3+b22bs2b)ba(7d5d2+bd3+2bd+bd23))T (3.20)

    and

    p(1,a(1+d)2b)T. (3.21)

    To achieve the normalization <p,q>=1, where <,> means the standard scalar product in R2, we can write the normalized vectors p and q as

    q=(1,(3d2+6d2bd2bds+3+b22bs2b)ba(7d5d2+bd3+2bd+bd23))T,p=(2(d24d+bd+b3)5d2+14d4bd2bds+b2s4b2bs+9,a(1+d)(d24d+bd+b3)b(5d2+14d4bd2bds+b2s4b2bs+9))T. (3.22)

    Assume that xt=Xtx, yt=Yty and J(E)=J(x,y). We transform the fixed point E2 of system (1.1) into the origin (0,0). By Taylor expansion, system (1.1) can be taken as

    (xtyt)J(E)(xtyt)+(F1(xt,yt)F2(xt,yt)), (3.23)

    where

    F1(xt,yt)=rFx2taytxt+O(X4t),F2(xt,yt)=bxtyt+O(X4t),

    where Xt=(xt,yt)T. System (3.23) can be expressed as

    (xt+1yt+1)=J(E)(xtyt)+12B(xt,xt)+16C(xt,xt,xt)+O(x4t), (3.24)

    where B(x,y)=(B1(x,y)B2(x,y)) and C(x,y,u)=(C1(x,y,u)C2(x,y,u)) are symmetric multi-linear vector functions of x,y,uR2 that are defined as follows:

    B1(x,y)=2j,k=12F1ξjξk|ξ=0xjyk=2rFx1y1ax1y2ax2y1,B2(x,y)=2j,k=12F2ξjξk|ξ=0xjyk=bx1y2+bx2y1,C1(x,y,u)2=j,k,l=13F1ξjξkξl|ξ=0xjykul=0,C2(x,y,u)=2j,k,l=13F2ξjξkξl|ξ=0xjykul=0. (3.25)

    The sign of the coefficient β2(rFB) determines the direction of the flip bifurcation, and it can be calculated as follows:

    β2(rF)=16p,C(q,q,q)12p,B(q,(JI)1B(q,q)). (3.26)

    From the above analysis, we have the following theorem:

    Theorem 3.2.1. Suppose that E2 is a positive unique fixed point of system (1.1). If the condition given by (3.19) holds, β2(rF)0 and the parameter r changes its value in a small vicinity of FBE2; then, system (1.1) experiences a flip bifurcation at the fixed point E2. Moreover, if β2(rF)>0 (β2(rF)<0)), then there exists stable (unstable) period-2 orbits that bifurcate from E2.

    In dynamical systems, it is desirable to optimize the system according to some performance criteria and avoid chaos. In order to stabilize the chaos in the case of unstable trajectories of the systems, different control techniques are applied. Therefore, chaos control is one of the attractive topics of current studies [8,9,10].

    Chaos control can be obtained by applying various methods to the discrete-time systems. To control the chaos in system (1.1), we chose to explore a feedback control strategy and apply the OGY method to system (1.1) in the sense of [22]. In order to apply OGY techniques to system (1.1), we rewrite it as follows:

    Xt+1=rXt(1Xt)aXtYt+s=f(Xt,Yt,r),Yt+1=bXtYtdYt=g(Xt,Yt,r), (4.1)

    where r is a controlling parameter. Furthermore, r is restricted to lie in some small interval |rr0|<μ, where μ>0 and r0 denotes the nominal value belonging to the chaotic region. We apply a stabilizing feedback control strategy in order to move the trajectory toward the desired orbit. Suppose that (X,Y) is an unstable fixed point of system (1.1) in the chaotic region produced by the emergence of Neimark-Sacker bifurcation; then, system (4.1) can be approximated in the neighborhood of the unstable fixed point (X,Y) by the following linear map:

    [Xt+1XYt+1Y]A[XtXYtY]+B[rr0], (4.2)

    where

    A=[f(X,Y,r0)Xf(X,Y,r0)Yg(X,Y,r0)Xg(X,Y,r0)Y]=[r0d2+2r0d+r+sb2bdbb(d+1)a(d+1)br0bd+r0brd22r0dr0+sb2bdba(d+1)1]

    and

    B=[f(X,Y,r0)rg(X,Y,r0)r]=[d+1b(d+1)2b20].

    On the other hand, system (4.1) can be controlled provided that the following matrix exists:

    C=[B:AB]=[d+1b(d+1)2b2(bd1)(r0d2+2r0d+r+sb2bdb)b30(bd1)(r0bd+r0br0d22r0dr0+sb2bdb)b2a],

    which has a rank of 2. We suppose that [rr0]=K[XtXYtY]; where K=[ρ1 ρ2], then, system (4.2) can be written as follows:

    [Xt+1XYt+1Y][ABK][XtXYtY]. (4.3)

    The corresponding controlled system given by (1.1) can be written as

    Xt+1=(r0ρ1(XtX)ρ2(YtY))Xt(1Xt)aXtYt+s,Yt+1=bXtYtdYt, (4.4)

    where ρ1(XtX)ρ2(YtY) is the control force; ρ1 and ρ2 are defined as the feedback gains. Furthermore, the fixed point (X,Y) of system (4.4) is locally asymptotically stable if and only if both eigenvalues of the matrix ABK lie in an open unit disk. The Jacobian matrix ABK of the controlled system given by (4.4) can be written as follows:

    ABK=[r0d2+2r0d+r+sb2bdbb(d+1)(d+1b(d+1)2b2)ρ1(d+1)(ab+ρ2(bd1))b2r0bd+r0brd22r0dr0+sb2bdba(d+1)1].

    The characteristic equation of the Jacobian matrix ABK is given by

    P(λ)=λ2(tr(ABK))λ+det(ABK)=0. (4.5)

    Let λ1 and λ2 be roots of the characteristic equation given by (4.5); then,

    λ1+λ2=(d+1b(d+1)2b2)ρ1r0d2+2r0d+r+sb2bdbb(d+1)+1,λ1λ2=(d2bd+2db+1b2)ρ1+4r0d25r0d2r0bd+r0bd2+2r0bd+rbr0d3+sb2dbd2b(d+1)+(r02r0bd2+bd2+3r0d2sb2d+db2r0+2bd+3r0db2d4dbr0sb2+sb3+b+b2r0+rb22r0bab2)ρ2 (4.6)

    are valid. In order to obtain the lines of marginal stability, we must solve the equations λ1=±1 and λ1 λ2=1. These constraints ensure that the absolute value of λ1 and λ2 is less than 1. Assume that λ1λ2=1; then, the second part of Eq (4.6) implies that

    L1:=((ad3+(a(3b))d2+(2ab+3a)dab+a)b2a(d+1))ρ1+((d4r0+(2r0b+b+4r)d3+(sb2+b2r0+3b+6r0b26r0b)d2b2a(d+1)+(4r06r0b+3b2sb2+2b2r0+sb32b2)dsb2+sb3+b+(b2+1)r0b(b+2r0))b2a(d+1))ρ2+(ab2+ab2r04abr0)d2+(5abr0+2ab2r02ab2+ab3s)dab22abr0+ab2r0b2a(d+1)=0.

    Moreover, we suppose that λ1=1; then, using Eq (4.6), we get

    L2:=(r0d3+(2r0b+b+3r0)d2+(sb2+b2r+2b+3r0b24r0b)dabr0+ab2r0+ab3sab2ab2)ρ2ρ2d2b+(ab22abr0+ab2r0)dab2+abr0+ab2r0+ab3sab2ab2=0.

    Finally, taking λ1=1 and using Eq (4.6), we get

    L3:=(2(bd+bd22d1)b2)ρ1+(sba+d(r0s1)s1a+2r0d(d+2)+(d+1)22r0ab+rd(d2+3d+3)+r0ab2)ρ2+(d3+(5+b)d2+(2b7)d+b3)r0b(d+1)+bd2+(sb2+2b)dsb2+3bb(d+1)=0.

    Then, stable eigenvalues lie within the triangular region in the ρ1ρ2 plane bounded by the straight lines L1, L2 and L3 for particular parametric values.

    In this section, we will give three examples to illustrate the results obtained in the above sections. Computer algebra systems such as Mathematica, Maple and Matlab were used for the calculations and graphic drawings.

    Example 5.1. With this example, our aim is to present numerical simulations to validate the above theoretical results by choosing r as the bifurcation parameter for system (1.1) around the fixed point E2 . The parameter values were taken from [23]. The Neimark-Sacker bifurcation point rNS was obtained as rNS=0.70339999998. By taking the parameter values (rNS,a,b,d,s)=(0.70339999998,3.5,4.5,0.25,0.9), the positive fixed point of system (1.1) has been evaluated as E2=(0.2777777778,0.7852698412). Using these parameter values, we get the following Jacobian matrix:

    JE2(rNS)=[2.4355555550.97222222233.5337142851].

    The eigenvalues were calculated to be:

    λ1,2=07177777775±0.6962722616i.

    Let q,pC2 be complex eigenvectors corresponding to λ1,2, respectively:

    q(0.19703694300.4861111113i,i)T,

    and

    p(0.71616575381.766857144i,i)T.

    To obtain the normalization p,q=1, we can take the normalized vectors to be

    q=(0.19703694300.4861111113i,i)T

    and

    p=(2.5375951972.313425848109i,1.233553221+0.4999999989i)T.

    To compute the coefficients of normal form, we can transform the fixed point E2 to the point (0,0) by changing the following variables:

    X=x0.2777777778,Y=y0.7852698412.

    Hence, we can compute the coefficients of the normal form of the system by using the formulas given by Eq (3.15), as follows:

    g20(rNS)=1.6458263202.815555559i,g11(rNS)=2.255064549+2.187500003i,g02(rNS)=3.4191588081.559444448i,g21(rNS)=0.

    From Eq (3.16), the critical part is then obtained as β2(rNS)=17.302665<0. Therefore, the Neimark-Sacker bifurcation is supercritical and it shows the correctness of Theorem 3.1.1. The bifurcation diagram, maximum Lyapunov exponents and this phase portraits of system (1.1) are shown in Figures 2 and 3.

    Figure 2.  Bifurcation diagram and Maximum Lyapunov exponents for system (1.1) with values of a = 3.5, b = 4.5, d = 0.25 and s = 0.9, as well as r[0.6, 0.74] and the initial value (x0,y0) = (0.26, 0.8). (a) Bifurcation diagram for Xt; (b) bifurcation diagram for Yt; (c) maximum Lyapunov exponents.
    Figure 3.  Bifurcation diagram and MLE for system (1.1) with values of a = 3.5, b = 4.5, d = 0.25 and r=0.7039999998, as well as s[0.8, 0.94] and the initial value (x0,y0) = (0.26, 0.8). (a) Bifurcation diagram for Xt; (b) bifurcation diagram for Yt; (c) maximum Lyapunov exponents.

    The bifurcation diagrams in Figure 2(a) and (b) show that the stability of E2 happens for r<0.7039999998, but there is a loss of stability at r=0.7039999998, and an attractive invariant curve appears if r>0.7039999998. We computed the maximum Lyapunov exponents to detect the presence of chaos in the system. The existence of chaotic regions in the parameter space is clearly visible in Figure 2(c).

    Now, we present the Neimark-Sacker bifurcation diagrams and maximum Lyapunov exponents for system (1.1) by choosing the immigration parameter s instead of r as a bifurcation parameter.

    The phase portraits for different values of r are displayed in Figure 4, which clearly depicts the process of how a smooth invariant circle bifurcates from the stable fixed point E2=(0.27777778,0.7852698412). When r exceeds 0.7039999998, there appears a circular curve enclosing the fixed point E2, and its radius becomes larger with respect to the growth of r.

    Figure 4.  Phase portraits of system (1.1) for different values of r.

    Example 5.2. Let (rF,a,b,d,s)=(0.1953418483,2.1,2,0.1,1.7). The positive fixed point of system (1.1) has been evaluated as E2=(0.55,1.037529963). Also, the flip bifurcation point rF was obtained as rF=0.1953418483. Using these parameter values, we can obtain the following Jacobian matrix:

    JE2(rF)=[2.1983471071.1552.0750599261.000000000].

    The eigenvalues were calculated to be

    λ1=1,λ2=0.1983471077.

    Let q,pR2 be eigenvectors of J and the transposed matrix JT, respectively, for λ1=1. Then, we have

    J(rFB)q=q and JT(rFB)=p.

    Direct computation gives

    q(0.6939646278,0.7200090940)T,

    and

    p(0.86596925300,0.5000972431)T.

    To obtain the normalization p,q=1, we can take the normalized vectors to be

    q=(0.6939646278,0.7200090940)T

    and

    p=(3.595061266,2.076147879)T.

    Then, the symmetric multi-linear vector functions can be obtained as follows:

    B1(x,y)=2rFx1y12.1(x1y2+x2y1)B2(x,y)=2(x1y2+x2y1),C1(x,y,u)=0,C2(x,y,u)=0.

    From Eq (3.26), the critical part is obtained as β2(rF)=0.6745625580>0. Therefore, unique and stable period-2 orbits bifurcate from E2, and it shows the correctness of Theorem 3.2.1. The bifurcation diagram, maximum Lyapunov exponents and this phase portraits of system (1.1) are shown in Figures 5 and 6.

    Figure 5.  Bifurcation diagram and MLE for system (1.1) with values of a = 2.1, b = 2, d = 0.1 and s = 1.7, as well as the initial value (x0,y0) = (0.91, 0.61). (a) Bifurcation diagram for Xt; (b) bifurcation diagram for Yt; (c) maximum Lyapunov exponents.
    Figure 6.  Bifurcation diagram and MLE for system (1.1) with values of a = 2.1, b = 2, d = 0.1 and r=0.1953418483, and as well as the initial value (x0,y0) = (0.54, 1.02). (a) Bifurcation diagram for Xt; (b) bifurcation diagram for Yt; (c) maximum Lyapunov exponents.

    From Figure 5(a) and (b), we see that the fixed point E2=(0.55,1.037529963) is stable for r<0.1953418483, and loses its stability at  the flip bifurcation parameter value r=0.1953418483. We also observe that, if r>0.1953418483, the flip bifurcation giving a 2-periodic orbit occurs. Moreover, the maximum Lyapunov exponents were computed, and the existence of chaotic regions in the parameter space is clearly visible in Figure 5(c).

    We present the flip bifurcation diagrams and maximum Lyapunov exponents for system (1.1), which were obtained by choosing the immigration parameter s instead of r as a bifurcation parameter.

    In Figure 7, we present the phase portraits of system(1.1) for different values of r.

    Figure 7.  Phase portraits of system (1.1) for different values of r.

    Example 5.3.In order to discuss the OGY feedback control method for system (1.1), we take r0=0.73 and (a,b,d,s)=(3.5,4.5,0.25,0.9). In this case, system (1.1) has unique positive fixed point (X,Y)=(0.2777777778,0.7906349207), which is unstable. Then, the corresponding controlled system is given by

    Xn+1=(0.73ρ1(Xn0.2777777778)ρ2(Yn0.7906349207))Xn(1Xn)aXnYn+s,Yn+1=bXnYndYn, (5.1)

    where K=[ρ1 ρ2] is a matrix and (X,Y)=(0.2777777778,0.7906349207) is an unstable fixed point of system (1.1). We have

    A=[2.4427777780.97222222223.5578571431],
    B=[0.20061728400]

    and

    C=[B:AB]=[0.20061728400.490063443200.7137676369].

    Then, it is easy to check that the rank of Matrix C is 2. Therefore, system (5.1) is controllable. Then, the Jacobian matrix ABK of the controlled system given by (5.1) is

    ABK=[2.4427777780.2006172840ρ10.97222222220.2006172840ρ23.5578571531]. (5.2)

    Moreover, the lines L1, L2 and L3 of marginal stability are given by

    L1=0.016250.2006172840ρ1+0.7137676369ρ2=0,L2=3.459027778+0.7137676369ρ2=0,L3=0.5734722220.4012345680ρ1+0.7137676369ρ2=0.

    Then, the stable triangular region bounded by the marginal lines L1, L2 and L3 for the controlled system given by (5.1) is shown in Figure 8.

    Figure 8.  Triangular stability region bounded by L1, L2 and L3 for the controlled system given by (5.1).

    Since most predator-prey systems in the wild are not isolated, it is important to consider the effects of the presence of some number of immigrants. The immigration effect of the predator-prey system has rarely been studied in the literature. Therefore, we considered a discrete-time predator-prey system with a constant-rate immigration of the prey and investigated the dynamics of modified system (1.1) around coexistence fixed point in this paper.

    In the case of s=0, we obtained the trivial fixed point, non-predatory fixed point for r>1 and coexistence fixed point for r>bbd1. In the case of s0, the origin is not a fixed point. System (1.1) has two fixed points, namely, E1 and E2. Moreover, simple calculations showed that, if (d+1)b(r1)+sb2>r(d+1)2, then system (1.1) has a unique positive fixed point E2. As r is varied, the system exhibits several complicated dynamical behaviors, including the emergence of Neimark-Sacker and flip bifurcations, an invariant circle, a period-2 orbit and chaotic sets. We presented some figures to explain the theoretical analysis in Section 4. In Figure 2 and Figure 5, the bifurcation, maximum Lyapunov exponents and phase portraits are given for some parameter values. Figure 2 and Figure 5 show that, when the bifurcation parameter r passes a critical bifurcation value, the stability of the coexistence fixed point of the system changes from stable to unstable. In Figure 2(c) and Figure 5(c), the computation of the maximum Lyapunov exponents confirm the presence of chaotic behavior in the system. In Figure 4 and Figure 7, the phase portraits of the system for different values of the r parameter are plotted. Moreover, the Neimark-Sacker bifurcation has been successfully controlled with an OGY control strategy. From our numerical investigation, it is clear that an OGY method based on a feedback control strategy can restore the stability. This control method is effective in order to advance or completely eliminate the chaos due to the emergence of Neimark-Sacker bifurcation. In this study, because of consistency with the biological facts, we used the same parameter values that were used in [23]. As a final point, we can say that the parameter r has a strong effect on the stability of system (1.1) for the control of two populations.

    On the other hand, the biological meaning of the term s being included in the system is that a small number of immigrants is added persistently to the prey population in each generation. A very small immigration into the system occurs as a stabilizing factor for the system. A positive immigration factor is enough to change the quality of the dynamics of the system. Specially, adding the term s contributes to asymptotic stability. This means that adding a small immigration s may imply that cyclic population can be stabilized.

    When the immigration effect is applied to a small number of species, the level of the species increases significantly. In the same way, when the harvesting effect is applied to a large number of species, the level of the species is significantly reduced [38]. Thanks to these two external factors, all species in the ecosystem can be preserved to stabilize the unstable system. In a future study, we plan to explore the effects of these factors on the population model.

    The authors declare no conflicts of interest in this paper.



    [1] W. Diffie, M. Hellman, New directions in cryptography, IEEE T. Inform. Theory, 22 (1976), 644–654. http://dx.doi.org/10.1109/TIT.1976.1055638 doi: 10.1109/TIT.1976.1055638
    [2] D. Boneh, V. Shoup, A graduate course in applied cryptography, 0.6 Eds., 2023, unpublished work. Available from: https://toc.cryptobook.us/book.pdf.
    [3] E. Bresson, Y. Lakhnech, L. Mazaré, B. Warinschi, A generalization of DDH with applications to protocol analysis and computational soundness, In: Advances in cryptology-CRYPTO 2007, Berlin: Springer, 2007,482–499. http://dx.doi.org/10.1007/978-3-540-74143-5_27
    [4] R. Rivest, A. Shamir, L. Adleman, A method for obtaining digital signatures and public-key cryptosystems, Commun. ACM, 21 (1978), 120–126. http://dx.doi.org/10.1145/359340.359342 doi: 10.1145/359340.359342
    [5] P. Shor, Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer, SIAM Rev., 41 (1999), 303–332. http://dx.doi.org/10.1137/S0036144598347011 doi: 10.1137/S0036144598347011
    [6] Submission requirements and evaluation criteria for the post-quantum cryptography standardization process, NIST Computer Security Resource Center, 2016. Available from: https://csrc.nist.gov/CSRC/media/Projects/Post-Quantum-Cryptography/documents/call-for-proposals-final-dec-2016.pdf.
    [7] Post-quantum cryptography, selected algorithms 2022, NIST Computer Security Resource Center, 2022. Available from: https://csrc.nist.gov/Projects/post-quantum-cryptography/selected-algorithms-2022.
    [8] Post-quantum cryptography, round 4 submissions, NIST Computer Security Resource Center, 2022. Available from: https://csrc.nist.gov/Projects/post-quantum-cryptography/round-4-submissions.
    [9] K. Ko, S. Lee, J. Cheon, J. Han, J. Kang, C. Park, New public-key cryptosystem using braid groups, In: Advances in cryptology-CRYPTO 2000, Berlin: Springer, 2000,166–183. http://dx.doi.org/10.1007/3-540-44598-6_10
    [10] V. Shpilrain, A. Ushakov, The conjugacy search problem in public key cryptography: unnecessary and insufficient, AAECC, 17 (2006), 285–289. http://dx.doi.org/10.1007/s00200-006-0009-6 doi: 10.1007/s00200-006-0009-6
    [11] E. Sakalauskas, N. Listopadskis, P. Tvarijonas, Key agreement protocol (KAP) based on matrix power function, Proceedings of Sixth International Conference on Information Research and Applications, 2008, 92–96.
    [12] A. Mihalkovich, E. Sakalauskas, Asymmetric cipher based on MPF and its security parameters evaluation, Lietuvos Matematikos Rinkinys, 53 (2012), 72–77. http://dx.doi.org/10.15388/LMR.A.2012.13 doi: 10.15388/LMR.A.2012.13
    [13] J. Liu, H. Zhang, J. Jia, A linear algebra attack on the non-commuting cryptography class based on matrix power function, In: Information security and cryptology, Cham: Springer, 2017,343–354. http://dx.doi.org/10.1007/978-3-319-54705-3_21
    [14] E. Sakalauskas, A. Mihalkovich, Improved asymmetric cipher based on matrix power function resistant to linear algebra attack, Informatica, 28 (2017), 517–524. http://dx.doi.org/10.15388/Informatica.2017.142 doi: 10.15388/Informatica.2017.142
    [15] A. Mihalkovich, M. Levinskas, Investigation of matrix power asymmetric cipher resistant to linear algebra attack, In: Information and software technologies, Cham: Springer, 2019,197–208. http://dx.doi.org/10.1007/978-3-030-30275-7_16
    [16] A. Mihalkovich, E. Sakalauskas, K. Luksys, Key exchange protocol defined over a non-commuting group based on an NP-complete decisional problem, Symmetry, 12 (2020), 1389. http://dx.doi.org/10.3390/sym12091389 doi: 10.3390/sym12091389
    [17] A. Mihalkovich, K. Luksys, E. Sakalauskas, Sigma identification protocol construction based on MPF defined over non-commuting platform group, Mathematics, 10 (2022), 2649. http://dx.doi.org/10.3390/math10152649 doi: 10.3390/math10152649
    [18] M. Durcheva, TrES: tropical encryption scheme based on double key exchange, European Journal of Information Technologies and Computer Science, 2 (2022), 11–17. http://dx.doi.org/10.24018/compute.2022.2.4.70 doi: 10.24018/compute.2022.2.4.70
    [19] X. Jiang, H. Huang, G. Pan, Cryptanalysis of tropical encryption scheme based on double key exchange, Journal of Cyber Security and Mobility, 12 (2023), 205–220. http://dx.doi.org/10.13052/jcsm2245-1439.1224 doi: 10.13052/jcsm2245-1439.1224
    [20] Modular maximal-cyclic group, Groupprops, 2023, Available from: https://groupprops.subwiki.org/wiki/Modular_maximal-cyclic_group.
    [21] H. Grundman, T. Smith, Automatic realizability of Galois groups of order 16, Proc. Amer. Math. Soc., 124 (1996), 2631–2640. http://dx.doi.org/10.1090/S0002-9939-96-03345-X doi: 10.1090/S0002-9939-96-03345-X
    [22] H. Grundman, T. Smith, Realizability and automatic realizability of Galois groups of order 32, Open Math., 8 (2010), 244–260. https://doi.org/10.2478/s11533-009-0072-x doi: 10.2478/s11533-009-0072-x
    [23] H. Grundman, T. Smith, Galois realizability of groups of order 64, Centr. Eur. J. Math., 8 (2010), 846–854. http://dx.doi.org/10.2478/s11533-010-0052-1 doi: 10.2478/s11533-010-0052-1
    [24] A. Mihalkovich, E. Sakalauskas, M. Levinskas, Key exchange protocol based on the matrix power function defined over M16, In: Intelligent computing, Cham: Springer, 2022,511–531. http://dx.doi.org/10.1007/978-3-031-10467-1_32
    [25] Faithful irreducible representation of M16, Groupprops, 2023, Available from: https://groupprops.subwiki.org/wiki/Faithful_irreducible_representation_of_M16.
    [26] J. Faugère, A. Joux, Algebraic cryptanalysis of hidden field equation (HFE) cryptosystems using Gröbner bases, In: Advances in cryptology-CRYPTO 2003, Berlin: Springer, 2003, 44–60. http://dx.doi.org/10.1007/978-3-540-45146-4_3
    [27] A. Kipnis, J. Patarin, L. Goubin, Unbalanced oil and vinegar signature schemes, Advances in Cryptology-EUROCRYPT'99, Berlin: Springer, 1999,206–222. http://dx.doi.org/10.1007/3-540-48910-X_15
    [28] R. Benadjila, T. Feneuil, M. Rivain, MQ on my mind: post-quantum signatures from the non-structured multivariate quadratic problem, Proceedings of IEEE 9th European Symposium on Security and Privacy, 2024,468–485. http://dx.doi.org/10.1109/EuroSP60621.2024.00032
    [29] A. Mihalkovich, J. Zitkevicius, On the decisional problem based on matrix power function defined over non-commutative group, Mathematical Models in Engineering, 10 (2024), 1–9. http://dx.doi.org/10.21595/mme.2024.24071 doi: 10.21595/mme.2024.24071
  • This article has been cited by:

    1. Deniz ELMACI, Figen KANGALGİL, Stability, Neimark-Sacker Bifurcation Analysis of a Prey-Predator Model with Strong Allee Effect and Chaos Control, 2022, 15, 1307-9085, 775, 10.18185/erzifbed.1207680
    2. Karima Mokni, Mohamed Ch-Chaoui, A Darwinian Beverton–Holt model with immigration effect, 2024, 217, 03784754, 244, 10.1016/j.matcom.2023.10.022
    3. Md. Jasim Uddin, Sarker Md. Sohel Rana, Hiroki Sayama, Qualitative Analysis of the Discretization of a Continuous Fractional Order Prey-Predator Model with the Effects of Harvesting and Immigration in the Population, 2024, 2024, 1099-0526, 1, 10.1155/2024/8855142
    4. Md. Jasim Uddin, P. K. Santra, Sarker Md Sohel Rana, G.s. Mahapatra, Chaotic Dynamics of the Fractional Order Predator-Prey Model Incorporating Gompertz Growth on Prey with Ivlev Functional Response, 2024, 6, 2687-4539, 192, 10.51537/chaos.1300754
    5. Karima Mokni, Halima Ben Ali, Bapan Ghosh, Mohamed Ch-Chaoui, Nonlinear dynamics of a Darwinian Ricker system with strong Allee effect and immigration, 2025, 229, 03784754, 789, 10.1016/j.matcom.2024.10.017
    6. Rasha M. Yaseen, May M. Helal, Kaushik Dehingia, Ahmed A. Mohsen, Effect of the Fear Factor and Prey Refuge in an Asymmetric Predator–Prey Model, 2024, 54, 0103-9733, 10.1007/s13538-024-01594-9
    7. Deniz Elmacı, Figen Kangalgil, Complex Dynamics of a Discrete Prey–Predator Model Exposing to Harvesting and Allee Effect on the Prey Species with Chaos Control, 2024, 34, 0218-1274, 10.1142/S0218127424501141
    8. Érika Diz-Pita, Global dynamics of a predator-prey system with immigration in both species, 2024, 32, 2688-1594, 762, 10.3934/era.2024036
    9. Halima Benali, Karima Mokni, Hajar Mouhsine, Mohamed Ch-Chaoui, Unveiling Complexity: A Discrete-Time Prey–Predator Model with Immigration Effects, 2024, 2731-8095, 10.1007/s40995-024-01742-5
    10. Seval Işık, Figen Kangalgil, MD. Jasim Uddin, Sarker MD. Sohel Rana, An Investigation of the Discrete-Time Model Subject to Immigration, Harvesting, and Allee Effect, 2025, 35, 0218-1274, 10.1142/S0218127425500221
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(631) PDF downloads(43) Cited by(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog