# Feedback Control and Parameter Invasion for a Discrete Competitive Lotka–Volterra System

Authors: Li Xu; Shanshan Lou; Panqi Xu; Guang Zhang
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

State feedback is used to stabilize the Turing instability at the unstable equilibrium point of a discrete competitive Lotka–Volterra system. In addition, a regularization method is applied to parameter inversion for the given Turing system and numerical simulation can verify the effectiveness of the algorithm. Furthermore, how less or more sample data and dependence on the initial state affect estimation procedure are tested.


---

## Body

## 1. Introduction

In theoretical ecology, the models governed by difference equations are used to characterize the interactions of species when the size of the population is rarely small [1]. For example, a two-species competitive discrete-time system which reduced from the continuous one with the forward Euler scheme can be obtained as follows [2]:(1)xt+1=r1xt1-a11xt-a12ytyt+1=r2yt1-a21xt-a22ytwhere xt and yt are the quantities of the two species at tth generation, r1>0 and r2>0 are growth rates of the respective species, a11 and a22 represent the strength of the intraspecific competition, and a12 and a21 are the strength of the interspecific competition.It is well known that the role of spatial heterogeneity and dispersal in the dynamics of populations has been the subject of much research, both theoretical and experimental, such as the role of dispersal in the maintenance of patchiness, or spatial population variation. It is also a fact that the motion of individuals is random and isotropic, i.e., without any preferred direction, and the individuals are also absolute ones in microscopic sense, and each isolated individual exchanges materials and information by diffusion with its neighbors [3–5]. Then it is reasonable to consider a 2D spatially discrete reaction-diffusion system as follows [6]:(2)xi,jt+1=r1xi,jt1-a11xi,jt-a12yi,jt+D1∇2xi,jtyi,jt+1=r2yi,jt1-a21xi,jt-a22yi,jt+D2∇2yi,jtfor i,j∈1,2,···,m=[1,m] and t∈1,2,···,n,···, where m,n are positive integers. Here, ∇2 is discrete Laplace operator (3)∇2xijt=xi+1,jt+xi,j+1t+xi-1,jt+xi,j-1t-4xijtand (4)∇2yijt=yi+1,jt+yi,j+1t+yi-1,jt+yi,j-1t-4yijt.This also indicates the coupling or diffusion from the cells to the left (i,j-1) and right (i,j-1) and top (i+1,j) and bottom (i-1,j) respectively.The partial difference systems (space-time discrete systems) for biological patterns resulted from diffusion-driven instability (Turing instability) in plants and animals and the set of equilibrium patterned solutions have been studied in great detail over the last several years; for example, see [6–11]. Models such as discrete competition reaction-diffusion models (2) have also been proposed to explain a wide variety of biological patterning processes and various patterns such as spiral wave, trigger wave, stripes, and chaotic Turing structure can be exhibited in the Turing instability region [6].In some situations, one may wish to recover its stability by means of some ways, in order to move the trajectory towards the desired orbit. Feedback control is an effective one, which is of significance in the control procedure of ecology balance. If one may wish to alter the positions of positive equilibrium and to obtain its stability, to achieve the aim, one of the techniques used is to alter system structurally by introducing “indirect control” variables. Though there are many works on the single species or multispecies competition systems with feedback controls [12–16]. To the best of the authors’ knowledge, there are still no scholars who are investigating the stability property of the 2D spatially discrete reaction-diffusion competitive system with feedback controls; this motivates us to propose such a model as follows:(5)xi,jt+1=r1xi,jt1-a11xi,jt-a12yi,jt-b1ui,jt+D∇2xi,jtyi,jt+1=r2yi,jt1-a21xi,jt-a22yi,jt-b2vi,jt+D∇2yi,jtui,jt+1=1-e1ui,jtt+d1xi,jtvi,jt+1=1-e2vi,jtt+d2yi,jtwhere u and v are feedback control variables and the parameters c1,c2,d1,d2,e1, and e2 are positive constants.It may be a fact that stability analysis, in mathematics, belongs to direct problem which pays attention on the dynamical behavior of the system state. Otherwise, one may be much more concerned to know what reason or what system environment to result into the current state. In mathematics, it can belong to inverse problems of identification and determination of parameters. On the other hand, an important and often difficult step from the viewpoint of testing models against experimental observations is the determination of model parameters from limited data when details of the mechanistic steps involved are not known. Parameter identification is the foundation of state estimation, controller design, diagnosis and fault detection, etc. Therefore, much research on such parameters inversion problems based on equation or reaction-diffusion systems has emerged; for example, see [17–20]. And, there is a few work on parameters identification or estimation for Turing systems [21–23]. Parameter identification for the classic Gierer-Meinhardt reaction-diffusion system is considered in [21], which can result in diffusion-driven instability, and the parameters are extended in time and space and used as distributed control variables. In [22], it is shown how using different combinations of spatial and temporal data can improve parameter estimation in a postulated model and how postprocessing with sensitivity analysis can be used to address the complexity issue. The authors present a Bayesian inference approach to solve both the parameter and the state estimation problem for stochastic reaction-diffusion systems in [23]. So far, to our knowledge, there have been very few corresponding research works focusing on discrete Turing system. In [24], a parameter estimation method called regularization method is applied to estimate the discrete Lotka–Volterra cooperative system and comparison experiments are also done using the regularization method and least square method to confirm the algorithm’s effectiveness. Similarly, a revised parameter estimation method can be used for the 2D spatially discrete reaction-diffusion competitive system (2) and the case to be discussed in this work will hold linear feature.So the paper is organized as follows. After a brief presentation of the model with diffusion for a completely symmetric case of the system (2), local instability conditions can be deduced combining linearization method and inner product technique for a the symmetric system with feedback control in Section 2. A parameter estimation method called regularization method is applied to estimate the discrete Lotka–Volterra competitive system in Section 3, and numerical examples will also support this inference. The final section is the conclusion.


## 2. Feedback Control and Its Stability

In this section, a completely symmetric discrete Lotka–Volterra competitive system can be given as follows:(6)xi,jt+1=rxi,jt1-xi,jt-ayi,jt+D∇2xi,jtyi,jt+1=ryi,jt1-axi,jt-yi,jt+D∇2yi,jtwith periodic boundary conditions(7)xi,0t=xi,mt,xi,1t=xi,m+1t,x0,jt=xm,jt,x1,jt=xm+1,jtand(8)yi,0t=yi,mt,yi,1t=yi,m+1t,y0,jt=ym,jt,y1,jt=ym+1,jtfor i,j∈1,2,···,m=[1,m] and t∈1,2,···,n,···, where m,n are positive integers.From [6], the above system is a diffusion-driven unstable one at the nontrivial coexistence point ((r-1)/ra+1,(r-1)/ra+1) when the condition Dkls2>3-r holds for some positive number D and kls2, where (9)kls2=4sin2⁡l-1πm+sin2⁡s-1πmforl,s∈1,m.In order to stabilize the orbit at an unstable equilibrium point of system (2), we use the state feedback control method and indirect control variables are added; then we can get the system(10)xi,jt+1=rxi,jt1-xi,jt-ayi,jt-b1ui,jt+D∇2xi,jtyi,jt+1=ryi,jt1-axi,jt-yi,jt-b2vi,jt+D∇2yi,jtui,jt+1=1-e1ui,jtt+d1xi,jtvi,jt+1=1-e2vi,jtt+d2yi,jtwhere u and v are feedback control variables and the parameters c1,c2,d1,d2,e1, and e2 are positive constants.According to the definition of fixed points, the fixed points of map (10) are solved by direct calculation yielding four fixed points,(11)E00,0,0,0,E10,a2r-1a2+b2c2r,0,a2c2r-1a2+b2c2a2r,E2a1r-1a1+b1c1r,0,a1c1r-1a1+b1c1a1r,0,E3x∗,y∗,u∗,v∗,where(12)x∗=r-1a1a2+b2c2-a2ara1+b1c1a2+b2c2-a1a2a2,y∗=r-1a2a1+b1c1-a1ara1+b1c1a2+b2c2-a1a2a2,u∗=c1a1x∗,v∗=c2a2y∗.There exists positive fixed point if and only if r>1 and(13)a1+b1c1a2+b2c2a1a2a2>1,a<min⁡1+b1c1a1,1+b2c2a2,or(14)a1+b1c1a2+b2c2a1a2a2<1,a>max⁡1+b1c1a1,1+b2c2a2hold. And we only care about the positive fixed point in this paper.For the reaction-diffusion system (10), we linearize about the steady state, to get(15)wijt+1=Awijt+D∇2wijt,with the periodic boundary conditions (16)wi,0t=wi,mt,wi,1t=wi,m+1t,w0,jt=wm,jt,w1,jt=wm+1,jt,where (17)A=fxfyfufvgxgygugvd101-e100d201-e2x∗,y∗,u∗,v∗,wijt=xijt-x∗yijt-y∗uijt-u∗vijt-v∗=xijtyijtuijtvijt,and (18)D=DD0·0In order to study instability of (15), we firstly consider eigenvalues(19)∇2Xij+μXij=0,with the periodic boundary conditions(20)Xi,0=Xi,m,Xi,1=Xi,m+1,X0,j=Xm,j,X1,j=Xm+1,j.In view of [11], the eigenvalue problems (19)-(20) have the eigenvalues(21)μl,s=4sin2⁡l-1πm+sin2⁡s-1πm=kls2forl,s∈1,m.Then taking the inner product of (15), respectively, with the corresponding eigenfunction Xtsij of the eigenvalue μl,s, we see that(22)∑i,j=1mXlsijxijt+1=fx∑i,j=1mXlsijxijt+fy∑i,j=1mXlsijyijt+D∑i,j=1mXlsij∇2xijt∑i,j=1mXtsijyijt+1=gx∑i,j=1mXlsijxijt+gy∑i,j=1mXlsijyijt+D∑i,j=1mXlsij∇2yij∑i,j=1mXtsijuijt+1=1-e1∑i,j=1mXtsijuijt+1+d1∑i,j=1mXlsijxijt∑i,j=1mXtsijvijt+1=1-e1∑i,j=1mXtsijvijt+1+d1∑i,j=1mXlsijyijt.LetXt=∑i,j=1mXlsijxij, Yt=∑i,j=1mXtsijuij, Ut=∑i,j=1mXlsijuij, and Vt=∑i,j=1mXtsijvij and use the periodic boundary conditions (20); then we have(23)Xt+1=fxXt+fyYt+fuUt+fvVt-Dkls2UtYt+1=gxXt+gyYt+guUt+gvVt-Dkls2VtUt+1=1-e1Ut+d1XtVt+1=1-e2Vt+d1Yt,or(24)Xt+1=fx-DCXt+fyYt+fuUt+fvVtYt+1=gxXt+gy-Dkls2Yt+guUt+gvVtUt+1=1-e1Ut+d1XtVt+1=1-e2Vt+d1Yt,which has the eigenvalue equation (25)λ4+R1λ3+R2λ2+R3λ+R4=0,where(26)R1=2Dkls2-fx-gy+e1+e2-2,R2=Dkls2-fxDkls2-gy+e1+e2-22Dkls2-fx-gy+e1-1e2-1-gxfy-d1fu,R3=Dkls2-fxDkls2-gye1+e2-2+2Dkls2-fx-gye1-1e2-1-gxfye1+e2-2+d2fv-d1fygu-fuDkls2+e2-gy-1,R4=Dkls2-fxDkls2-gye1-1e2-1+d2fve1-1-fygxe1-1e2-1-d1e1-1fygu-d1fuDkls2-gye2-1+d1d2fugv-fvgu.and f=rx(1-x-ay-b1u), g=ry(1-ax-y-b2v), fx=fx(x∗,y∗,u∗,v∗), fy=fy(x∗,y∗,u∗,v∗), fu=fu(x∗,y∗,u∗,v∗), fv=fv(x∗,y∗,u∗,v∗), gx=gx(x∗,y∗,u∗,v∗), gy=gy(x∗,y∗,u∗,v∗), gu=gu(x∗,y∗,u∗,v∗), and gv=gv(x∗,y∗,u∗,v∗).According to the Routh–Hurwitz criterion, we can draw the following conclusion.Theorem 1.
The positive homogeneous steady stateE is stable if the following conditions are satisfied:(27)R1>0,R1R2-R3>0,R3R1R2-R3-R12R4>0,R4>0When conditions (27) are not satisfied, the positive steady state is unstable and bifurcations may occur.


## 3. Parameter Inversion

From above section, it is clear that if all the system parameters are given, we can solve the concentration distribution with time and space, which is possible to generate patterns of species distribution with the system evolution. It can be called direct problem. However, some parameters cannot be determined in advance or measured directly. Thus, we need to estimate the parameters via mathematical algorithms by means of data which can be measured. It can be so-called inverse problem, namely, parameter identification. The purpose of the section is to determine the parameters which best fit the simulations to the measurements.We only consider the above system (6), which can be denoted as follows:(28)xi,jt+1yi,jt+1=xi,jt-xi,jt2-xi,jtyi,jt∇2xi,jtyi,jt-yi,jt2xi,jtyi,jt∇2yi,jtrarDor(29)xi,jt+1yi,jt+1=xi,jt-xi,jt2-xi,jtyi,jt∇2xi,jtyi,jt-yi,jt2xi,jtyi,jt∇2yi,jtrADwhere A=ar.Let(30)U1,n=xnijm×m,U2,n=ynijm×m,H1U1,n,U2,n=h1xnij,ynijm×m=xi,jt-xi,jt2,-xi,jtyi,jt,∇2xi,jtm×m,H2U1,n,U2,n=h2xnij,ynijm×m=yi,jt-yi,jt2,-xi,jtyi,jt,∇2yi,jtm×m,C=r,A,DT,where C=(cj)3=(r,A,D)T.System (29) can be represented as the following form:(31)Un+1=U1,n+1U2,n+1=H1U1,n,U2,nH2U1,n,U2,nC=HUnC.Although we obtain the above linear form,xi,jt-(xi,jt)2,-xi,jtyi,jt,and ∇2xi,jt (or yi,jt-(yi,jt)2,-xi,jtyi,jt, and ∇2yi,jt) will possess serious colinearity, which will result in the fact that many traditional parameter identification methods, proposed in the past, just like least square method, maximum likelihood method, etc., are not effective. To deal with the problem, some revised parameter identification methods, such as regularization method, have been put forward [24–26]. Here, the regularization method will be used to deal with the parameter identification.Suppose thatUn+1-Un+1∗≤h,H(Un)-H∗(Un)≤δ,and h,δ≥0, where (H∗(Un),Un+1∗) are real value. The least-squares problem associated with (31) can be written as(32)minC∈R3⁡HUnC-Un+1,which is generally ill-posed based on Hadamard [25].To overcome the ill-posedness, we can employ the regularization approach to get a best fitted solution. By using the Tikhonov regularization, we can convert problem (32) into the following form:(33)minC∈R3⁡MθHUn,Un+1,C=minC∈R3⁡HUnC-Un+12+θC2,where θ>0 is regularization parameter. The solution of (33) is given by(34)Cθ=HTUnHUn+θI-1HTUnUn+1,where the superscript T denotes transposition. The value of the regularization parameter θ determines how well Cθ approximates C and how sensitive Cθ is to the error η=(h,δ) in the available data. Because of the importance of being able to determine a suitable value of the regularization parameter, a large number of parameter choice rules have been developed. When the error level η=(h,δ) is known, the most well known of these rule will be Morozov discrepancy principle [27, 28]. When no accurate estimate of η=(h,δ) is known, the rules include the L-curve criterion, generalized cross validation, error estimation methods based on extrapolation, and Regińska’s method; for example, see [29–31] for discussions and comparisons of the performance of these and several other parameter choice rules. This paper will apply L-curve criterion for choosing the parameter. The detail of this method can be seen in [25].Then, the proposed algorithm to obtain parametersC consists of the following steps:(1) By using the L-curve approach, select an appropriate regularization parameterθ∗ and compute the value of Cθ∗ by using relation (34).(2) Check parameters positive or negative.The convergence and feasibility of the Tikhonov regularization have been studied well. Therefore, we no longer discuss the details. Next we will make some numerical experiments to test the effectiveness for the parameter estimation method, namely, generating data computationally for a the model with known parameters and then testing our recovery of parameters from the data which is available as panel data.The test data can be generated by means of iteration computation whenr=2.98,a=0.5,D=0.1, and the small amplitude random perturbation is 1% around the steady state. Firstly, we select the 100000 interaction values of x,y. To confirm the above algorithm’s effectiveness, comparison experiments are done using the regularization method (RM) and least square method (LS) and the numerical result corresponding to the real value (RV) can be found in Table 1, with the regularization parameter θ∗=1.0875×10-6. Furthermore, to intuitively reflect the fit goodness, patterns generated from the system whose parameters are obtained by using the regularization method and least square method can be shown in Figure 1.Table 1
The estimation result when the sample values are generated from rand initial value fromt=1 to t=100000.
r a D RV 2.98 0.5 0.1  LS 2.9839 0.5211 0.0947 RM 2.9758 0.5051 0.0963Figure 1
Selective patterns resulting from different parameters values and fixed initial value whent=100000. (a) r=2.98,a=0.5,D=0.1. (b) r=2.9839,a=0.5211,D=0.0947. (c) r=2.9758,a=0.5051,D=0.0963.
(a) (b) (c)It may be true that patterns formation can depend on not only fluctuations of system parameters but also variation of the initial conditions. Initial distribution of an immobile reactive species can affect pattern formation. Then we make some numerical experiments to check out the estimation effective when the system parameters are fixed but the initial values are different. We still chooser=2.98,a=0.5,D=0.1, and special initial value is given by (35)xi,j0=rand×0.01+r-1ra+1i=j=m20other,and(36)yi,j0=rand×0.01+r-1ra+1i=j=m20other,Then comparison experiments are also done using the regularization method (RM) and least square method (LS); the numerical result can be found in Table2. Furthermore, patterns generated from the system whose parameters are obtained by using the regularization method and least square method can be shown in Figure 2. Although there exist interesting spiral structures in Figure 2(a), the symmetry breaking around the fixed point does not emerge in Figure 2(b) and a stable pattern of rhombi shapes can be observed from Figure 2(c). Comparing the results of two experiments, it may be a fact that factors that might affect estimation procedure depend on the initial state which results in different time evolution process and less sample data.Table 2
The estimation result when the sample values are generated from special initial value fromt=1 to t=10000.
r a D RV 2.98 0.5 0.1  LS 2.9892 0.6735 0.0906 RM 2.9710 0.5201 0.0924Figure 2
Selective patterns resulting from different parameters values and fixed initial value whent=50000. (a) r=2.98,a=0.5,D=0.1. (b) r=2.9892,a=0.6735,D=0.0906. (c) r=2.9710,a=0.5201,D=0.0924.
(a) (b) (c)


## 4. Conclusion

We give several concluding remarks in this section.(1) In this work, the local stability of a positive interior equilibrium for a discrete competitive Lotka–Volterra system model with feedback controls is investigated by means of inner product technique and eigenvalue analysis. Otherwise, the global stability still is not been obtained and the method of global Lyapunov functions may be applied in our further work.(2) The numerical results obtained demonstrate the feasibility and potential advantages of applying a regularization method to parameter estimation in discrete competitive Lotka–Volterra systems. Some facts can also show that more sample data can increase estimation accuracy relieve colinearity, minute estimation error can result in significant effect, and different patterns can emerge, which may prove the fact that the nonlinear system remarkably depends on the system parameters and initial state.



---

*Source: 7473208-2018-08-16.xml*