Theorems and Application of Local Activity of CNN with Five State Variables and One Port

Abstract

Coupled nonlinear dynamical systems have been widely studied recently. However, the dynamical properties of these systems are difficult to deal with. The local activity of cellular neural network (CNN) has provided a powerful tool for studying the emergence of complex patterns in a homogeneous lattice, which is composed of coupled cells. In this paper, the analytical criteria for the local activity in reaction-diffusion CNN with five state variables and one port are presented, which consists of four theorems, including a serial of inequalities involving CNN parameters. These theorems can be used for calculating the bifurcation diagram to determine or analyze the emergence of complex dynamic patterns, such as chaos. As a case study, a reaction-diffusion CNN of hepatitis B Virus (HBV) mutation-selection model is analyzed and simulated, the bifurcation diagram is calculated. Using the diagram, numerical simulations of this CNN model provide reasonable explanations of complex mutant phenomena during therapy. Therefore, it is demonstrated that the local activity of CNN provides a practical tool for the complex dynamics study of some coupled nonlinear systems.

1. Introduction

Coupled nonlinear dynamical systems have been widely studied in recent years. However, the dynamical properties of these systems are difficult to deal with. Although the research on emergence and complexity has gained much attention during the past decades, the determination, prediction, and control of the complex patterns generated from high-dimensional coupled nonlinear systems are still far from perfect. Nature abounds with complex patterns and structures emerging from homogeneous media, and the local activity is the origin of these complexities [1, 2]. The cellular neural network (CNN), firstly introduced by Chua and Yang [3] as an implementable alternative to fully connected Hopfield neural network, has been widely studied for image processing, robotic, biological versions, and higher brain functions, and so on [3]. Many of the coupled nonlinear systems can be modeled and studied via the CNN paradigm [4]. The local activity proposed by Chua asserts that a wide spectrum of complex behaviors may exist if the cell parameters of the corresponding CNN are chosen in or nearby the edge of chaos [2, 4]. There have been quite a few new methods developed for complex systems [5–8], and local activity has attracted the attention of many researchers. Now, local activity has been successfully applied to the research of complex patterns generated from several CNNs in physical, biological, and chemical domains, such as Fitzhugh-Nagumo equation [9], Brusselator equation [10], Gierer-Meinhart equation [11], Oregonator equation [12], Hodgkin-Huxley equation [13], Van Der Pol equation [14], the biochemical model [15], coupled excitable cell model [16], tumor growth and immune model [17], Lorenz model [18], advanced image processing [19], Rossler equation [20], images analysis [21, 22], data prediction [23], neutron transport equation [24], vision safety [25], retinomorphic model [26], and theory research [27–30], and so forth.

Although Chua presents the main theorem of local activity at a cell equilibrium point [1, 2], it is actually difficult to “test” directly the complex patterns of the high-dimensional coupled nonlinear systems, since the theorem contains no recipe for finding whether a variable actually exists or not. It is necessary to develop some mathematical criteria according to the numbers of the variables and ports; that is the topic addressed in this paper.

The remaining of this paper is organized as follows. The local activity of CNN is introduced in Section 2. A set of theorems for testing the local activity of reaction-diffusion CNN with five state variables and one port are set up in Section 3. As an application of the theorems, a coupled reaction-diffusion CNN of hepatitis B Virus (HBV) mutation-selection model is introduced, aiming at describing HBV mutation in the therapeutic process. The bifurcation diagrams of this CNN are developed and some numerical simulations are presented in Section 4. Concluding remarks are given in Section 5.

2. Local Activity Theory of CNN

The CNN architecture is composed of a two-dimensional 𝑀×𝑁 array of cells. Each cell is denoted by 𝐶(𝑖,𝑗), where 𝑖=1,2,…,𝑀,𝑗=1,2,…,𝑁. The dynamics of each cell is given by the equation:𝑥̇_𝑖𝑗=−𝑥_𝑖𝑗+^𝑟_𝑘=−𝑟^𝑟_𝑙=−𝑟𝑎_𝑘𝑙𝑦_{𝑖+𝑘𝑗+𝑙}+^𝑟_𝑘=−𝑟^𝑟_𝑙=−𝑟𝑏_𝑘𝑙𝑢_{𝑖+𝑘𝑗+𝑙}+𝑧_𝑖𝑗, (1) where 𝑥_𝑖𝑗,𝑦_𝑖𝑗,𝑢_𝑖𝑗 are the state, output, and input variables of the cell, respectively. 𝑎_𝑘,𝑙,𝑏_𝑘,𝑙,𝑧_𝑖𝑗 are the elements of the A-template, the B-template, and threshold, respectively. 𝑟 is the radius of influence sphere. The output 𝑦_𝑖𝑗 is the piece-wise linear function given by𝑦_𝑖,𝑗=12𝑥_𝑖,𝑗+1−𝑥_𝑖,𝑗−1,𝑖=1,2,…,𝑀;𝑗=1,2,…,𝑁. (2) Clearly, CNN with different template elements may have different functions.

A vast majority of active homogeneous media that are known to exhibit complexity are modeled by a reaction-diffusion partial differential equation (PDE):𝜕𝑥_𝑖𝜕𝑡=𝑓_𝑖(𝑋)+𝐷_𝑖𝜕²𝑥_𝑖𝜕𝑥²+𝜕²𝑥_𝑖𝜕𝑦²+𝜕²𝑥_𝑖𝜕𝑧²,𝑖=1,2,…,𝑛, (3) where 𝑋=(𝑥₁,𝑥₂,…𝑥_𝑛) is state variables, (𝑥,𝑦,𝑧) is spatial coordinates, 𝑓_𝑖(𝑥₁,𝑥₂,…,𝑥_𝑛) is a coupled nonlinear vector function called the kinetic term, and 𝐷₁,𝐷₂,…,𝐷_𝑛 are constants called diffusion coefficients. Replacing the Laplace in above formulation by its discrete version yields𝜕²𝑥_𝑖𝜕𝑥²+𝜕²𝑥_𝑖𝜕𝑦²+𝜕²𝑥_𝑖𝜕𝑧²⟶∇²𝑋_{𝛼,𝛽,𝛾}, (4) where∇²𝑋_{𝛼,𝛽,𝛾𝑖}=𝑥_𝑖𝛼−1,𝛽,𝛾+𝑥_𝑖𝛼+1,𝛽,𝛾+𝑥_𝑖𝛼,𝛽−1,𝛾+𝑥_𝑖𝛼,𝛽+1,𝛾+𝑥_𝑖𝛼,𝛽,𝛾−1+𝑥_𝑖𝛼,𝛽,𝛾+1−6𝑥_𝑖𝛼,𝛽,𝛾. (5) Chua et al. have introduced reaction-diffusion CNN equations:𝑋̇_{𝛼,𝛽,𝛾}=𝑓𝑋_{𝛼,𝛽,𝛾}+𝐷∇²𝑋_{𝛼,𝛽,𝛾}, (6) where 𝐷=diag(𝐷₁,𝐷₂,…𝐷_𝑛), 𝑋̇_{𝛼,𝛽,𝛾} denotes the state variable located at a point in three-dimensional space with spatial coordinates. Chua refers to the process of transforming a PDE into a reaction-diffusion CNN [2].

From Chua and his collaborators’ point, PDEs are merely mathematical abstractions of nature, and the concept of a continuum is in fact an idealization of reality. Even the collection of all electrons in a solid does not form a continuum, because much volume separating the electrons from the nucleus represents a vast empty space [2]. Reaction-diffusion CNNs have been used to model some phenomena with important practical backgrounds, which were described by PDEs.

Generally speaking, in a reaction-diffusion CNN, every cell has 𝑛 state variables, but only 𝑚 (𝑚≤𝑛) state variables couple directly to their nearest neighbors via “reaction-diffusion”. Consequently, each cell has the following state equations:𝑉̇_𝑎=𝑓_𝑎𝑉_𝑎,𝑉_𝑏+𝐼_𝑎,𝑉̇_𝑏=𝑓_𝑎𝑉_𝑎,𝑉_𝑏, (7) where𝑉_𝑎=𝑉₁,𝑉₂,…,𝑉_𝑚^T,𝑉_𝑏=𝑉_𝑚+1,𝑉_𝑚+2,…,𝑉_𝑛^T,𝑓_𝑎=𝑓₁𝑉_𝑎,𝑉_𝑏,𝑓₂𝑉_𝑎,𝑉_𝑏,…,𝑓_𝑚𝑉_𝑎,𝑉_𝑏^T,𝑓_𝑏=𝑓_𝑚+1𝑉_𝑎,𝑉_𝑏,𝑓_𝑚+2𝑉_𝑎,𝑉_𝑏,…,𝑓_𝑛𝑉_𝑎,𝑉_𝑏^T,𝐼_𝑎=𝐷₁∇²𝑉₁,𝐷₂∇²𝑉₂,…,𝐷_𝑚∇²𝑉_𝑚^T. (8)

The cell equilibrium point 𝑄_𝑖=(𝑉^𝑖_𝑎,𝑉^𝑖_𝑏)(∈R^𝑛) of equation (7) can be determined via𝑓_𝑎𝑉_𝑎,𝑉_𝑏𝑓=0,_𝑏𝑉_𝑎,𝑉_𝑏=0. (9)

The Jacobian matrix at the equilibrium point 𝑄_𝑖 has the following form:𝐽𝑄_𝑖=𝑎_𝑘𝑙𝑄_{𝑖𝑛×𝑛}=A_𝑎𝑎𝑄_𝑖A_𝑎𝑏𝑄_𝑖A_𝑏𝑎𝑄_𝑖A_𝑏𝑏𝑄_𝑖, (10) where A_𝑘𝑙(𝑄_𝑖) are called cell parameters andA_𝑎𝑎𝑄_𝑖=𝜕𝑓₁𝜕𝑉₁⋯𝜕𝑓₁𝜕𝑉_𝑚⋮⋱⋮𝜕𝑓_𝑚𝜕𝑥₁⋯𝜕𝑓_𝑚𝜕𝑥_𝑚,A_𝑎𝑏𝑄_𝑖=𝜕𝑓₁𝜕𝑉_𝑚+1⋯𝜕𝑓₁𝜕𝑉_𝑛⋮⋱⋮𝜕𝑓_𝑚𝜕𝑥_𝑚+1⋯𝜕𝑓_𝑚𝜕𝑥_𝑛,A_𝑏𝑎𝑄_𝑖=𝜕𝑓_𝑚+1𝜕𝑉₁⋯𝜕𝑓_𝑚+1𝜕𝑉_𝑚⋮⋱⋮𝜕𝑓_𝑛𝜕𝑥₁⋯𝜕𝑓_𝑛𝜕𝑥_𝑚,A_𝑏𝑏𝑄_𝑖=𝜕𝑓_𝑚+1𝜕𝑉_𝑚+1⋯𝜕𝑓_𝑚+1𝜕𝑉_𝑛⋮⋱⋮𝜕𝑓_𝑛𝜕𝑥_𝑚+1⋯𝜕𝑓_𝑛𝜕𝑥_𝑛. (11) The local state equations at the cell equilibrium point 𝑄_𝑖 are defined via𝑉̇_𝑎=A_𝑎𝑎𝑉_𝑎+A_𝑎𝑏𝑉_𝑏+𝐼_𝑎,𝑉̇_𝑏=A_𝑏𝑎𝑉_𝑎+A_𝑏𝑏𝑉_𝑏. (12)

Definition 1. 𝑌_𝑄(𝑠)=𝑠𝐼−A_𝑎𝑎−A_𝑎𝑏𝑠𝐼−A_𝑏𝑏⁻¹A_𝑏𝑎 (13) is called the admittance matrix at the cell equilibrium point 𝑄_𝑖.

Lemma 2. A reaction-diffusion CNN cell is called locally active at the equilibrium point 𝑄_𝑖 if and only if, its admittance matrix at 𝑄_𝑖 satisfies at least one of the following four conditions [4].(1)𝑌_𝑄(𝑠) has a pole in Re[𝑠]>0.(2)𝑌^𝐻_𝑄(𝑖𝜔)=𝑌_𝑄(𝑖𝜔)+𝑌_𝑄(𝑖𝜔)<0 for some 𝜔=𝜔₀, where 𝜔₀ is any real number.(3)𝑌_𝑄(𝑠) has a simple pole 𝑠=𝑖𝜔_𝑝 on the imaginary axis, where its associated residue matrix: 𝑘₁=lim_{𝑠→𝑖𝜔𝑝}𝑠−𝑖𝜔_𝑝𝑌_𝑄(𝑠),if𝜔_𝑝<∞lim_{𝜔𝑝→∞}𝑌_𝑄𝑖𝜔_𝑝𝑖𝜔_𝑝,if𝜔_𝑝=∞ (14) is either a complex number or a negative real number.(4)𝑌_𝑄(𝑠) has a multiple pole on the imaginary axis.

Definition 3. The cell equilibrium point 𝑄_𝑖 is called stable if and only if, all the real parts of eigenvalue 𝜆_𝑖 of Jacobian matrix at the equilibrium point 𝑄_𝑖 are negative [2].

Definition 4. A “reaction-diffusion” CNN with 𝑛 state variables and 𝑚 ports is said to be operating on the “edge of chaos” with respect to an equilibrium point 𝑄_𝑖 if and only if, 𝑄_𝑖 is both locally active and stable when 𝐼_𝑎=0 [4].

Using the above lemma and definitions, the bifurcation of CNN with respect to an equilibrium point can be divided into three parts: the edge of chaos domains (the locally active and stable domains), the locally active and unstable domains, and the locally passive domains. Numerical simulations indicated that many complex dynamical behaviors, such as oscillatory patterns, chaotic patterns, or divergent patterns, may emerge if the selected cell parameters are located in or nearby the edge of chaos domains.

3. Analytical Criteria for Local Activity of CNN with Five State Variables and One Port

For the reaction-diffusion CNN with five state variables and one port, its local state equations have the form𝑉̇_𝑎=A_𝑎𝑎𝑉_𝑎+A_𝑎𝑏𝑉_𝑏+𝐼_𝑎,𝑉̇_𝑏=A_𝑏𝑎𝑉_𝑎+A_𝑏𝑏𝑉_𝑏, (15) where𝑉_𝑎=𝑉₁,𝑉_𝑏=𝑉₂𝑉₃𝑉₄𝑉₅^𝑇,𝐼_𝑎=𝐼₁,A_𝑎𝑎=𝑎₁₁,A_𝑎𝑏=𝑎₁₂𝑎₁₃𝑎₁₄𝑎₁₅,A_𝑏𝑎=𝑎₂₁𝑎₃₁𝑎₄₁𝑎₅₁^𝑇,A_𝑏𝑎=𝑎₂₁𝑎₃₁𝑎₄₁𝑎₅₁,A_𝑏𝑏=𝑎₂₂𝑎₂₃𝑎₂₄𝑎₂₅𝑎₃₂𝑎₃₃𝑎₃₄𝑎₃₅𝑎₄₂𝑎₄₃𝑎₄₄𝑎₄₅𝑎₅₂𝑎₅₃𝑎₅₄𝑎₅₅. (16)

The corresponding CNN cell admittance matrix 𝑌_𝑄(𝑠) is given by [1].𝑌_𝑄(𝑠)=𝑠𝐼−A_𝑎𝑎−A_𝑎𝑏𝑠I−A_𝑏𝑏⁻¹A_𝑏𝑎=𝑠−𝑎₁₁−𝑇₁𝑠³+𝐾₁𝑠²+𝐿₁𝑠+Δ₁𝑠⁴+𝑇𝑠³+𝐾𝑠²,+𝐿𝑠+Δ (17) where 𝑇,𝑇₁,𝐾,𝐾₁,𝐿,𝐿₁,Δ,Δ₁ are the parameters of 𝑎_𝑖𝑗’s.

Theorem 5. A necessary and sufficient condition for 𝑌_𝑄(𝑠) to satisfy condition (1) in Lemma 2 is that ∃𝑠, such that 𝑔(𝑠)=0(Re[𝑠]>0), and any one of the following conditions holds.(1)𝑓(𝑠)≠0.(2)𝑓(𝑠)=0, and 𝑚>𝑛, where 𝑠 is 𝑚 and 𝑛 orders zero point of 𝑔(𝑠) and 𝑓(𝑠), respectively, where 𝑓(𝑠)=𝑇₁𝑠³+𝐾₁𝑠²+𝐿₁𝑠+Δ₁,𝑔(𝑠)=𝑠⁴+𝑇𝑠³+𝐾𝑠²+𝐿𝑠+Δ.

Proof. Obviously proved.

Denote𝐸=−𝑎₁₁,𝐹=−𝑇𝑇₁+𝐾₁−𝑎₁₁−2𝐾+𝑇²,𝑃=𝐿𝐿₁−𝐾Δ₁−Δ𝐾₁,𝑄=Δ₁+𝐾𝐾₁−𝑇𝐿₁−𝐿𝑇₁,𝐺=−𝑎₁₁2Δ+𝐾²−2𝐿𝑇−𝑄,𝐻=−𝑎₁₁𝐿²−2Δ𝐾−𝑃,𝐼=−𝑎₁₁Δ²−ΔΔ₁,𝑔(𝑄)=𝐸𝑄⁴+𝐹𝑄³+𝐺𝑄²ℎ(𝜆)+𝐻𝑄+𝐼,=−𝑇𝑇₁−𝐾₁𝜆³−𝑄𝜆²−𝑃𝜆−ΔΔ₁,𝜆^∗_1,2=−𝑄±𝑄²−3𝑇𝑇₁−𝐾₁𝑃3𝑇𝑇₁−𝐾₁,𝑝=−3𝐹²16𝐸²+𝐺𝐹2𝐸,𝑞=³32𝐸³−𝐹𝐺8𝐸²+𝐻,𝑤4𝐸_1,2=−1±𝑖√32,𝑞𝐷=²4+𝑝³27,𝐴_𝑗=−𝑞2±𝐷^1/21/3,𝑥₁=𝐴₁+𝐴₂,𝑥₂=𝑤₁𝐴₁+𝑤₂𝐴₂,𝑥₃=𝑤₂𝐴₁+𝑤₁𝐴₂,Ω_𝑖=𝑥_𝑗−𝐹4𝐸,𝑖=1,2,3. (18)

Theorem 6. Let the following parameters be defined as in Theorem 5, then 𝑌^𝐻_𝑄(𝑖𝑤)<0 for some 𝑤=𝑤₀∈𝑅 if any one of the following conditions holds.(1)𝑎₁₁>0. (2)𝑎₁₁=0,𝑇𝑇₁−𝐾₁>0. (3)𝑎₁₁=0,𝑇𝑇₁−𝐾₁=0,𝑄>0. (4)𝑎₁₁=0,𝑇𝑇₁−𝐾₁=0,𝑄<0,ΔΔ₁>0. (5)𝑎₁₁=0,𝑇𝑇₁−𝐾₁=0,𝑄<0,𝑃≥0,ΔΔ₁−𝑃²/𝑄/4>0,ΔΔ₁≤0, (6)𝑎₁₁=0,𝑇𝑇₁−𝐾₁=0,𝑄=0,𝑃>0. (7)𝑎₁₁=0,𝑇𝑇₁−𝐾₁=0,𝑄=0,𝑃≤0,ΔΔ₁>0. (8)𝑎₁₁=0,𝑇𝑇₁−𝐾₁<0,ΔΔ₁>0. (9)𝑎₁₁=0,𝑇𝑇₁−𝐾₁<0,ΔΔ₁≤0,and 𝜆^∗_𝑗≥0,ℎ(𝜆^∗_𝑗)<0, for 𝑗=1 or 2.(10)𝑎₁₁<0,𝐷>0,Ω₁>0,𝑔(Ω₁)<0. (11)𝑎₁₁<0,𝐷<0, and Ω_𝑗≥0,𝑔(Ω_𝑗)<0, for 𝑗=1,2 or 3.(12)𝑎₁₁<0,𝐷=0,𝑝=𝑞=0,𝑔(−𝐹/4𝐸)<0. (13)𝑎₁₁<0,𝐷=0,𝑞²/4=−𝑝³/27≠0,and Ω_𝑗≥0, 𝑔(Ω_𝑗)<0, for 𝑗=1 or 2.

Proof. 𝑌^𝐻_𝑄(𝑖𝜔)=𝑌_𝑄(𝑖𝜔)+𝑌_𝑄(𝑖𝜔)=2Re[𝑌_𝑄(𝑖𝜔)], so 𝑌_𝑄(𝑖𝜔) to satisfy condition (2) in Lemma 2 equals to Re[𝑌_𝑄(𝑖𝜔)]<0, Re𝑌_𝑄(𝑖𝜔)=Re𝑖𝜔−𝑎₁₁−𝑇₁(𝑖𝜔)³+𝐾₁(𝑖𝜔)²+𝐿₁(𝑖𝜔)+Δ₁(𝑖𝜔)⁴+𝑇(𝑖𝜔)³+𝐾(𝑖𝜔)²()=+𝐿𝑖𝜔+Δ𝐸𝜔⁸+𝐹𝜔⁶+𝐺𝜔⁴+𝐻𝜔²+𝐼𝜔⁴−𝐾𝜔²+Δ²+𝐿𝜔−𝑇𝜔². (19)(1)If 𝑎₁₁>0, then Re[𝑌_𝑄(𝑖𝜔)]<0 when 𝜔 is large enough (See (1) of Theorem 6).(2)If 𝑎₁₁=0, then Re𝑌_𝑄(𝑖𝜔)=−𝑇𝑇₁−𝐾₁𝜔⁶+𝑄𝜔⁴+𝑃𝜔²+ΔΔ₁𝜔⁴−𝐾𝜔²+Δ²+𝐿𝜔−𝑇𝜔². (20)Let 𝑓(𝜆)=−𝑄𝜆²−𝑃𝜆−ΔΔ₁,(I) If 𝑇𝑇₁−𝐾₁>0, then Re[𝑌_𝑄(𝑖𝜔)]<0 when 𝜔 is large enough (See (2) of Theorem 6).(II)If 𝑇𝑇₁−𝐾₁=0, then(i)If 𝑄>0, then Re[𝑌_𝑄(𝑖𝜔)]<0 when 𝜔 is large enough (See (3) of Theorem 6).(ii)If 𝑄<0,(a)If ΔΔ₁>0,∃𝜔₀∈R, such that Re[𝑌_𝑄(𝑖𝜔₀)]<0 (See (4) of Theorem 6).(b)If ΔΔ₁≤0, solve 𝑓′(𝜆^∗)=0, we can get 𝜆^∗=−0.5𝑃/𝑄,𝑓(𝜆^∗)=0.25𝑃²/𝑄−ΔΔ₁.Then when 𝑃≥0,ΔΔ₁−0.25𝑃²/𝑄>0,∃𝜔₀>𝜆^∗, such that Re[𝑌_𝑄(𝑖𝜔₀)]<0 (See (5) of Theorem 6).(iii)If 𝑄=0, then Re[𝑌_𝑄(𝑖𝜔)]=−𝑃𝜔²−ΔΔ₁.(a)If 𝑃>0, then Re[𝑌_𝑄(𝑖𝜔)]<0 when 𝜔 is large enough (See (6) of Theorem 6).(b)If 𝑃≤0,ΔΔ₁>0, then ∃𝜔₀, such that Re[𝑌_𝑄(𝑖𝜔₀)]=−ΔΔ₁<0 (See (7) of Theorem 6).(III)If 𝑇𝑇₁−𝐾₁<0, let ℎ(𝜆)=−(𝑇𝑇₁−𝐾₁)𝜆³−𝑄𝜆²−𝑃𝜆−ΔΔ₁,(i)If ΔΔ₁>0, then ∃𝜔₀, such that Re[𝑌_𝑄(𝑖𝜔₀)]<0 (See (8) of Theorem 6).(ii)If ΔΔ₁≤0, solve ℎ(𝜆^∗)=0, we can get 𝜆^∗_1,2=−𝑄±𝑄²−3𝑇𝑇₁−𝐾₁𝑃3𝑇𝑇₁−𝐾₁. (21) Then, for 𝑖=1,2, if 𝜆^∗_𝑖≥0, ℎ(𝜆^∗_𝑖)<0, then ∃𝜔₀, such that Re[𝑌_𝑄(𝑖𝜔₀)]<0 (See (9) of Theorem 6).(3)If 𝑎₁₁<0, let 𝑔(𝑄)=𝐸𝑄⁴+𝐹𝑄³+𝐺𝑄²+𝐻𝑄+𝐼, then 𝑔′(𝑄)=4𝐸𝑄³+3𝐹𝑄²+2𝐺𝑄+𝐻. Let 𝑥=Ω+(𝐹/4𝐸), then the above becomes 𝑔′(𝑄)=4𝐸(𝑥³+𝑝𝑥+𝑞)=4𝐸𝑓(𝑥), then 𝑥_𝑖,𝑖=1,2,3 are the roots of 𝑓(𝑥)=0, Ω_𝑖 are the roots of 𝑔′(Ω)=0. If any one of the (10)–(13) of Theorem 6 holds, we can get Re[𝑌_𝑄(𝑖𝜔₀)]<0.So, if any one of conditions (1)–(13) holds, Re[𝑌_𝑄(𝑖𝜔₀)]<0. 𝑌_𝑄(𝑠) Satisfies condition (2) in Lemma 2. This completes the proof.

Theorem 7. For 𝑗=1, or 2, let 𝑤_𝑗=𝐾+2√Δ+(−1)^𝑗𝐾−2√Δ2,𝐴_𝑗=𝐿−3𝑤^∗_𝑗²,𝑤^∗_𝑗=𝐾+(−1)^𝑗√𝐾²−4Δ2,𝐵_𝑗=2𝐾𝑤^∗_𝑗−4𝑤^∗_𝑗³,𝐴_1𝑗=Δ₁−𝐾₁𝑤^∗_𝑗²,𝐵_1𝑗=𝐿₁𝑤^∗_𝑗−𝑇₁𝑤^∗_𝑗³. (22) Then 𝑌_𝑄(𝑠) satisfies condition (3) of Lemma 2, if any one of the following conditions holds.(I)Δ>0,𝐾>2√Δ,𝑇=𝐿=0, and any one of the following conditions holds.(1)𝐾₁𝑤²₁−Δ₁≠0. (2)𝐾₁𝑤²₁−Δ₁=0,(𝐿₁−𝑇₁𝑤²₁)(𝑤²₂−𝑤²₁)>0. (3)𝐾₁𝑤²₂−Δ₁≠0. (4)𝐾₁𝑤²₂−Δ₁=0,(𝐿₁−𝑇₁𝑤²₂)(𝑤²₁−𝑤²₂)>0. (II)𝐾>0,Δ₁≠0,Δ=0,𝐿=𝐾𝑇≠0, and any one of the following conditions holds.(1)𝑇Δ₁<0. (2)𝑇(𝐾−𝑇₁𝐾)−Δ₁+𝐾𝐾₁≠0. (3)𝑇(𝐾−𝑇₁𝐾)−Δ₁+𝐾𝐾₁=0,𝑇(Δ₁−𝐾𝐾₁)+𝐾(𝐿₁−𝑇₁𝐾)<0. (III)Δ=0,Δ₁𝐿>0. (IV)Δ<0, or𝐾>0,𝐾²−4Δ>0,and 2𝐿=𝑇(𝐾+√𝐾²−4Δ)(𝐾+√𝐾²−4Δ>0) or 2𝐿=𝑇(𝐾−√𝐾²−4Δ),(𝐾−√𝐾²−4Δ>0) and any one of the following conditions holds for 𝑗=1, or 2.(1)𝐴_𝑗𝐵_1𝑗−𝐴_1𝑗𝐵_𝑗≠0. (2)𝐴_𝑗𝐵_1𝑗−𝐴_1𝑗𝐵_𝑗=0,𝐴_𝑗𝐴_1𝑗−𝐵_𝑗𝐵_1𝑗>0.

Proof. Let 𝑓(𝑠)=𝑇₁𝑠³+𝐾₁𝑠²+𝐿₁𝑠+Δ₁,𝑔(𝑠)=𝑠⁴+𝑇𝑠³+𝐾𝑠²+𝐿𝑠+Δ, obviously, ∞ is not a single pole of 𝑌_𝑄(𝑠) on the imaginary axis.
If 𝑌_𝑄(𝑠) has a simple pole 𝑠=𝑖𝜔 on the imaginary axis, where its associated residue 𝑘₁=lim_{𝑠→𝑖𝜔}(𝑠−𝑖𝜔)𝑌_𝑄(𝑠)=lim_{𝑠→𝑖𝜔}()𝑓(𝑠)𝑠−𝑖𝜔𝑔(𝑠)=𝑓(𝑠)𝑔^′(𝑠)_{𝑠=𝑖𝜔} (23)is either a complex number or a negative real number, then 𝑘₁≠0, so 𝑓(𝑠)≠0, which implies that 𝑖𝜔 is not a zero point of 𝑓(𝑠)=0,𝑖𝜔  is not a removed pole of 𝑌_𝑄(𝑠).(I) If 𝑌_𝑄(𝑠) has four poles 𝑠=±𝑖𝜔₁,±𝑖𝜔₂(𝜔₁≠𝜔₂≠0) on the imaginary axis. In this case,𝑔(𝑠)=(𝑠²+𝜔²₁)(𝑠²+𝜔²₂)=𝑠⁴+(𝜔²₁+𝜔²₂)𝑠²+𝜔²₁𝜔²₂. Hence we obtain 𝑇=𝐿=0, 𝐾=𝜔²₁+𝜔²₂>0, Δ=𝜔²₁𝜔²₂>0. Then, we can get 𝐾+2√Δ=(𝜔₁+𝜔₂)², 𝐾−2√Δ=(𝜔₁−𝜔₂)², which implies that 𝐾>2√Δ, 𝜔_1,2=(𝐾+2√Δ±𝐾−2√Δ)/2. So, lim_{𝑠→±𝑖𝜔1}𝑠∓𝑖𝜔₁𝑌_𝑄(𝑠)=𝑤₁𝐿₁−𝑇₁𝑤²₁±𝑖𝐾₁𝑤²₁−Δ₁2𝑤₁𝑤²₂−𝑤²₁,lim_{𝑠→±𝑖𝜔2}𝑠∓𝑖𝜔₂𝑌_𝑄(𝑠)=𝑤₂𝐿₁−𝑇₁𝑤²₂±𝑖𝐾₁𝑤²₂−Δ₁2𝑤₂𝑤²₁−𝑤²₂. (24)(25) Then, when condition (I) in Theorem 7 holds, 𝑘₁ is a complex number or a negative real number. 𝑌_𝑄(𝑠) satisfies condition (3) in Lemma 2.(II)If 𝑌_𝑄(𝑠) has a simple pole 𝑠=0 and two conjugate poles ±𝑖𝜔(𝜔≠0) on the imaginary axis, and another pole is 𝑎≠0.In this case, it follows that Δ=0,Δ₁≠0, and 𝑔(𝑠) has the form: 𝑔(𝑠)=𝑠s²+𝜔²(s−𝑎)=s⁴−𝑎𝑠³+𝜔²s²−𝑎𝜔², (26) which implies that 𝑇=−𝑎,𝐾=𝜔²>0,𝐿=−𝑎𝜔²=𝐾𝑇,Δ=0,Δ₁=0, Therefore, 𝑔(𝑠)=𝑠s²+𝐾(s+𝑇) (27)(1)The residue of 𝑌_𝑄(𝑠) at 𝑠=0 is lim_𝑠→0𝑠𝑌_𝑄(𝑠)=Δ₁s²)=Δ+𝐾(s+𝑇𝐾𝑇¹. (28) Then, we conclude that if 𝐾>0,Δ₁≠0,Δ=0,𝐿=𝐾𝑇≠0,𝑇Δ₁<0, 𝑘₁ is a negative real number. 𝑌_𝑄(𝑠) satisfies condition (3) in Lemma 2 (See (1) of (II) in Theorem 7).(2)The residue of 𝑌_𝑄(𝑠) at 𝑠=±𝑖√𝐾 islim_{√𝑠→±𝑖𝐾}𝑠∓𝑖√𝐾𝑌_𝑄(𝑠)=𝑇Δ₁−𝐾𝐾₁+𝐾𝐿₁−𝑇₁𝐾∓𝑖√𝐾𝑇𝐿₁−𝑇₁𝐾−Δ₁+𝐾𝐾₁2𝐾𝐾+𝑇². (29)Consequently, we conclude that if (2) or (3) in (II) in Theorem 7 holds, 𝑘₁ is either an imaginary number or a negative real number. 𝑌_𝑄(𝑠) satisfies condition (3) in Lemma 2.(III)If 𝑌_𝑄(𝑠) has a simple pole 𝑠=0 on the imaginary axis, and the other poles are 𝑎_𝑖,Re[𝑎_𝑖]≠0,𝑖=1,2,3, it follows that Δ₁≠0,Δ=0, and 𝑔(𝑠) has the form 𝑔(𝑠)=𝑠s−𝑎₁s−𝑎₂s−𝑎₃=𝑠𝑠³−𝑎₁+𝑎₂+𝑎₃𝑠²+𝑎₁𝑎₂+𝑎₁𝑎₃+𝑎₂𝑎₃𝑠−𝑎₁𝑎₂𝑎₃. (30) Therefore we obtain that  Δ=0,𝑇=−(𝑎₁+𝑎₂+𝑎₃),𝐾=𝑎₁𝑎₂+𝑎₁𝑎₃+𝑎₂𝑎₃,𝐿=−𝑎₁𝑎₂𝑎₃≠0, hence the reside of 𝑌_𝑄(𝑠) at 𝑠=0 is lim_𝑠→0𝑠𝑌_𝑄(𝑠)=Δ₁𝑎₁𝑎₂𝑎₃Δ=−¹𝐿. (31) Then, when Δ=0,Δ₁𝐿>0, 𝑘₁ is a negative real number. 𝑌_𝑄(𝑠) satisfies condition (3) in Lemma 2 (See (III) of Theorem 7).(IV)If 𝑌_𝑄(𝑠) has two conjugate poles ±𝑖𝜔(𝜔>0) on the imaginary axis, and the other poles are Re[𝑎]≠0,Re[𝑏]≠0. In this case, 𝑔(𝑠) has the form 𝑔(𝑠)=(s−𝑎)(s−𝑏)s+𝜔²=𝑠⁴−(𝑎+𝑏)𝑠³+𝑎𝑏+𝜔²𝑠²−(𝑎+𝑏)𝜔²𝑠+𝑎𝑏𝜔². (32) Therefore, we obtain that 𝑇=−(𝑎+𝑏),𝐾=𝑎𝑏+𝜔²,𝐿=−(𝑎+𝑏)𝜔²,Δ=𝑎𝑏𝜔²≠0. Then, 𝑎𝑏=𝐾−𝜔²,Δ=(𝐾−𝜔²)𝜔²⇔𝜔⁴−𝐾𝜔²+Δ=0. Solving it, we have 𝜔^∗_𝑗=𝐾±√𝐾²−4Δ2,𝑗=1,2. (33) It implies that Δ<0 or 𝐾>0,𝐾²−4Δ≥0, and 𝑇=−(𝑎+𝑏)=𝐿/𝜔². Then, the residue of 𝑌_𝑄(𝑠) at 𝑠=±𝜔^∗_𝑗 is lim_{𝑠→±𝑖𝜔∗𝑗}𝑠∓𝑖𝜔^∗_𝑗𝑌_𝑄(𝑠)=−𝐴_1𝑗𝐴_𝑗+𝐵_1𝑗𝐵_𝑗±𝑖𝐴_1𝑗𝐵_𝑗−𝐴_𝑗𝐵_1𝑗𝐴²_𝑗+𝐵²_𝑗. (34) Hence, if condition (IV) in Theorem 7 holds, 𝑘₁<0 or it is an imaginary number. 𝑌_𝑄(𝑠) satisfies condition (3) in Lemma 2.
Therefore, when any one of conditions (I)–(IV) holds, 𝑌_𝑄(𝑠) satisfies condition (3) in Lemma 2. This completes the proof.

Theorem 8. 𝑌_𝑄(𝑠) has a multiple pole on the imaginary axis if any one of the following conditions holds.(I)Δ=𝐿=0,Δ₁≠0. (II)Δ=𝐿=𝐾=0,andΔ₁≠0or𝐿₁≠0. (III)Δ=𝐿=𝐾=𝑇=0, and Δ₁≠0 or 𝐿₁≠0 or 𝐾₁≠0.(IV)𝑇=𝐿=0,𝐾>0,Δ=(𝐾/2)², and any one of the following conditions holds.(1)2Δ₁≠𝐾𝐾₁. (2)2𝐿₁≠𝐾𝑇₁.

Proof. Let 𝑓(𝑠)=𝑇₁𝑠³+𝐾₁𝑠²+𝐿₁𝑠+Δ₁,𝑔(𝑠)=𝑠⁴+𝑇𝑠³+𝐾𝑠²+𝐿𝑠+Δ. Obviously, when the conditions (I)–(III) hold, 0 is the multiply poles of 𝑌_𝑄(𝑠).
If 𝑌_𝑄(𝑠) has two multiply nonzero poles ±𝑖𝜔(𝜔>0), then 𝑔(𝑠) has the form: 𝑔(𝑠)=𝑠²+𝜔²²=𝑠⁴+2𝜔²𝑠²+𝜔⁴, (35) which implies that 𝑇=𝐿=0,𝐾=2𝜔²>0,Δ=𝜔⁴=(𝐾/2)²>0.
If ±𝑖𝜔 are the multiply poles of 𝑌_𝑄(𝑠), then𝑓(𝑖𝜔)≠0 where 𝑓(±𝑖𝜔)=Δ₁−𝐾₁𝜔²±𝑖𝜔𝐿₁−𝑇₁𝜔²=Δ₁−𝐾𝐾₁2±𝑖𝜔𝐿₁−𝑇₁𝐾2. (36)
Obviously, when any one of (1)-(2) of (IV) in Theorem 8 holds, 𝑓(𝑖𝜔)≠0, any one of ±𝑖𝜔 is a multiply pole of 𝑌_𝑄(𝑠).
So, when any one of condition (I)–(IV) holds, 𝑌_𝑄(𝑠) satisfies condition (4) in Lemma 2. This completes the proof.
When any one of Theorems 5–8 holds, 𝑌_𝑄(𝑠) satisfies Lemma 2, which implies that the reaction-diffusion CNN with five state variables and one port at the equilibrium point is active.
These theorems can be implemented by a computer program for calculating the bifurcation diagram of the general corresponding CNN to determine emergence of complex dynamic patterns of the corresponding CNN.

4. Analysis and Simulations of Reaction-Diffusion CNN of HBV Mutation-Selection Model

Life systems consist of locally coupled homogeneous media. Mostly, dynamics of life systems are suitable to be described via locally connected reaction-diffusion CNNs. It may be expected that reaction-diffusion CNN will become a promising candidate for modeling life phenomena.

In Chapter 11 “Timing the emergence of resistance” (Page 110) of the book “Virus dynamic: mathematical principles of immunology and virology” (Oxford university press), Nowak et al. proposed a mathematical model which describes the mutation selection of HBV infection during the therapy [31]: 𝑑𝑥𝑑𝑡=𝜆−𝑑𝑥−𝑏𝑣𝑥−𝑏_𝑛𝑥𝑣_𝑛,𝑑𝑦()𝑑𝑡=𝑏1−𝑒𝑣𝑥−𝑎𝑦,𝑑𝑣𝑑𝑡=𝑘𝑦−𝑢𝑣,𝑑𝑦_𝑛𝑑𝑡=𝑏𝑒𝑣𝑥+𝑏_𝑛𝑥𝑣_𝑛−𝑎𝑦_𝑛,𝑑𝑣_𝑛𝑑𝑡=𝑘_𝑛𝑦_𝑛−𝑢𝑣_𝑛, (37) where the five variables—𝑥,𝑦,𝑣,𝑦_𝑛,𝑣_𝑛 represent the numbers of uninfected cells, infected cells infected by normal virus, normal virus, infected cells infected by mutated virus, and mutant viruses, respectively. 𝜆 is the rate of reproduction of uninfected cells. Uninfected cells die at rate 𝑑𝑥 and become infected at rate 𝑏𝑥𝑣 by normal virus and infected at rate 𝑏_𝑛𝑥𝑣_𝑛 by mutated virus. Infected cells infected by normal and mutated virus are removed at rate 𝑎𝑦 and 𝑎𝑦_𝑛, respectively. Normal virus is produced at rate 𝑘𝑦 and removed at rate 𝑢𝑣, mutated virus is produced at rate 𝑘𝑦_𝑛 and removed at rate 𝑢𝑣_𝑛. 𝑒 is the rate constant describing the probability of mutation of virus (usual 10⁻⁵−10⁻³), 𝑎,𝑏,𝑏_𝑛,𝑑,𝑒,𝑘,𝑘_𝑛,𝑢,𝜆 are positive constants. The model was briefly analyzed in Nowak’s book.

The reaction-diffusion CNN of HBV mutation selection of model has the form:𝑑𝑥_𝑖𝑗𝑑𝑡=𝜆−𝑑𝑥_𝑖𝑗−𝑏𝑥_𝑖𝑗𝑣_𝑖𝑗−𝑏_𝑛𝑥_𝑖𝑗𝑣_𝑛𝑖𝑗+𝐷₁∇²𝑥_𝑖𝑗,𝑑𝑦_𝑖𝑗()𝑥𝑑𝑡=𝑏1−𝑒_𝑖𝑗𝑣_𝑖𝑗−𝑎𝑦_𝑖𝑗,𝑑𝑣_𝑖𝑗𝑑𝑡=𝑘𝑦_𝑖𝑗−𝑢𝑣_𝑖𝑗,𝑑𝑦_𝑛𝑖𝑗𝑑𝑡=𝑏𝑒𝑥_𝑖𝑗𝑣_𝑖𝑗+𝑏_𝑛𝑥_𝑖𝑗𝑣_𝑛𝑖𝑗−𝑎𝑦_𝑛𝑖𝑗,𝑑𝑣_𝑛𝑖𝑗𝑑𝑡=𝑘_𝑛𝑦_𝑛𝑖𝑗−𝑢𝑣_𝑛𝑖𝑗, (38) where ∇²𝑥_𝑖𝑗=𝑥_𝑖−1𝑗+𝑥_𝑖+1𝑗+𝑥_𝑖𝑗−1+𝑥_𝑖𝑗+1−4𝑥_𝑖𝑗.

Let equation (38) be zeros (𝐷₁=0) and solve it, we can get the two equilibrium points:𝑄₁=𝜆𝑑𝑄,0,0,0,0,₂=𝑥₀,𝑢𝑎𝑢−𝑏_𝑛𝑘_𝑛𝑢𝑥₀𝜆−𝑑𝑥₀𝑏𝑘𝑥₀𝑏_𝑛𝑘_𝑛𝑒𝑥₀+𝑎𝑢−𝑘_𝑛𝑥₀,𝑎𝑢−𝑏_𝑛𝑘_𝑛𝑢𝑥₀𝜆−𝑑𝑥₀𝑏𝑥₀𝑏_𝑛𝑘_𝑛𝑒𝑥₀+𝑎𝑢−𝑘_𝑛𝑥₀,𝑒𝑢𝜆−𝑑𝑥₀𝑏_𝑛𝑘_𝑛𝑒𝑥₀+𝑎𝑢−𝑘_𝑛𝑥₀,𝑘_𝑛𝑒𝜆−𝑑𝑥₀𝑏_𝑛𝑘_𝑛𝑒𝑥₀+𝑎𝑢−𝑘_𝑛𝑥₀, (39)(40) where 𝑥₀=𝑎𝑢/((1−𝑒)𝑏𝑘) and 𝑄₁, 𝑄₂ stand for the patient’s complete recovery and HBV persistent infection, respectively.

Consequently, the Jacobian matrix at the equilibrium point 𝑄_𝑖(𝑖=1,2) is𝐽𝑄_𝑖=−𝑑−𝑏𝑣−𝑏_𝑛𝑣_𝑛0−𝑏𝑥0−𝑏_𝑛𝑥𝑏(1−𝑒)𝑣−𝑎𝑏(1−𝑒)𝑥000𝑘−𝑢00𝑏𝑒𝑣+𝑏_𝑛𝑣_𝑛0𝑏𝑒𝑥−𝑎𝑏_𝑛𝑥000𝑘_𝑛−𝑢. (41)

Taking 𝑘,𝑢 as variables, and 𝜆=10, 𝑎=0.5, 𝑏=0.01, 𝑏_𝑛=0.005, 𝑒=0.0001, 𝑘_𝑛=10, and 𝑑=0.01, using Theorems 5–8, we can calculate the bifurcation of the reaction-diffusion CNN model equation (38) at the equilibrium points 𝑄₁ and 𝑄₂ at 𝑘∈[0,40],𝑢∈[0,10], see Figures 1 and 2.

Figure 1 

Bifurcation diagrams of equation (38) at the equilibrium points 𝑄1 at 𝑘∈[0,40],𝑢∈[0,10].

Figure 2 

Bifurcation diagrams of equation (38) at the equilibrium points 𝑄2 at 𝑘∈[0,40],𝑢∈[0,10].

In Figures 1 and 2, the domains are coded as follows: edge of chaos (locally active and stable) domain (shown red), locally active and unstable domain (shown green) and locally passive domain (shown blue). From Figure 1(a), we can see that the bifurcation at equilibrium point 𝑄₁ does not exist at the edge of chaos domain.

Take 𝜆=10, 𝑘=0.01, 𝑎=0.5, 𝑏=0.01, 𝑏_𝑛=0.005, 𝑘_𝑛=10,𝑒=0.0001, and𝑘=1.0,3.0,4.9,5.1,10,24,39, 𝑢=2,5,9, we model the dynamic trajectories of equation (37) using MATLAB, see Table 1.

Table 1 

Cell parameters and correspongding dynamic properties of the reaction-diffusion CNN of HBV mutation-selection of HBV infection.

In the following discussions, we select some parameters in No. 8-14 and 𝑢=5, 𝑘=12.5. The simulation results are shown in Figures 3, 4, 5, 6, and 7. During the simulation, we reached a new conclusion.

(a)

(b)

(c)

(d)

Figure 3 

The kinetic trajectories of equation (37) when 𝑢=5, 𝑘=3.

(a)

(b)

(c)

(d)

Figure 4 

The kinetic trajectories of equation (37) when 𝑢=5,𝑘=4.9.

(a)

(b)

(c)

(d)

Figure 5 

The kinetic trajectories of equation (37) when 𝑢=5,𝑘=5.1.

(a)

(b)

(c)

(d)

Figure 6 

The kinetic trajectories of equation (37) when 𝑢=5,𝑘=12.5.

(a)

(b)

(c)

(d)

Figure 7 

The kinetic trajectories of equation (37) when 𝑢=5,𝑘=24.

From Table 1 and Figures 3–7, we can conclude that(I) when 𝑘 is smaller (less than 5),(1)these parameters are located in the green domain (the local and unstable domains);(2)regardless of the value of 𝑢, the dynamic pattern of equation (37) is convergent or divergent depending on initial values;(3)the No. 1, No. 4, and No. 5 variables in equation (37) increase and the No. 2 and No. 3 variables in equation (37) decrease to 0. This means the numbers of the mutant virus and of infected cells infected by mutant virus both increase, and the numbers of normal virus and of the infected cells infected by normal virus both decrease, even to near zero. Also, the No. 1 variable in equation (37), that is, the number of uninfected cells increases as compared to the initial number. All these indicate that the potency is perfect except for some virus mutation. The potency is ideal.(II)When 𝑘 is larger (greater than 5) but less than a threshold value (according to initial values and parameters, for example 12.5 in Figure 6), we can conclude the following.(1)These parameters are located in the red domain (edge of chaos).(2)Regardless of the value of 𝑢, the dynamic pattern of equation (37) is convergent.(3)The No. 1, No. 2, and No. 3 variables in equation (37) increase and No. 4 and No. 5 variables decrease to 0. This means that the number of uninfected cells, the numbers of the normal virus, and of the cells infected by normal virus all increase. Meanwhile, the numbers of the mutant virus and of the cells infected by mutant virus both decrease, even to near zero. All these imply that the drug cannot clean the normal virus, but can destroy the mutant virus and increase the infection cells. The potency is also ideal.(III)When 𝑘<40 and greater than a threshold value (according to initial values and parameters),(1)these parameters are located in the red domain (edge of chaos);(2)regardless of the 𝑢 value, the dynamic pattern of equation (37) is convergent;(3)The No. 2 and No. 3 variables in equation (37) increase and No. 1, No. 4, and No. 5 variables decrease, which means the numbers of the mutant virus and of the uninfected cells decrease, and the number of normal virus increases. These imply that although the drug can prevent the mutation of the HBV effectively, it may destroy uninfected cells and the liver. The potency is not ideal.

5. Conclusions

The local activity of CNN has provided a powerful tool for studying the emergence of complex patterns in a homogeneous lattice formed by coupled cells. Based on the local activity principle, the analytic criteria for the local activity in reaction-diffusion CNN with five state variables and one port are set up. The analytical criteria include four theorems, which provide the inequalities involving the parameters of the CNN. The inequalities can be used for calculating the bifurcation diagram to determine emergence of complex dynamic patterns of the reaction-diffusion CNN. As an application example, a reaction-diffusion CNN of HBV mutation-selection model is analyzed and simulated, and the bifurcation diagrams are calculated. Numerical simulations show this CNN model may explain certain complex mutant conditions during the therapy. We conclude that the local activity theory provides a practical tool for the study of the complex dynamics of certain coupled nonlinear systems.

Acknowledgment

This work is supported in part by NSFC 70890084, 60921061, 90920305.