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

Authors: Gang Xiong; Xisong Dong; Li Xie; Thomas Yang
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## 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.


---

## Body

## 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 Section2. 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-dimensionalM×N array of cells. Each cell is denoted by C(i,j), where i=1,2,…,M,j=1,2,…,N. The dynamics of each cell is given by the equation:(1)ẋij=-xij+∑k=-rr∑l=-rraklyi+kj+l+∑k=-rr∑l=-rrbklui+kj+l+zij,
where xij,yij,uij are the state, output, and input variables of the cell, respectively. ak,l,bk,l,zij are the elements of the A-template, the B-template, and threshold, respectively. r is the radius of influence sphere. The output yij is the piece-wise linear function given by(2)yi,j=12(|xi,j+1|-|xi,j-1|),i=1,2,…,M;j=1,2,…,N.
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):(3)∂xi∂t=fi(X)+Di(∂2xi∂x2+∂2xi∂y2+∂2xi∂z2),i=1,2,…,n,
where X=(x1,x2,…xn) is state variables, (x,y,z) is spatial coordinates, fi(x1,x2,…,xn) is a coupled nonlinear vector function called the kinetic term, and D1,D2,…,Dn are constants called diffusion coefficients. Replacing the Laplace in above formulation by its discrete version yields(4)∂2xi∂x2+∂2xi∂y2+∂2xi∂z2⟶∇2Xα,β,γ,
where(5)(∇2Xα,β,γ)i=xi(α-1,β,γ)+xi(α+1,β,γ)+xi(α,β-1,γ)+xi(α,β+1,γ)+xi(α,β,γ-1)+xi(α,β,γ+1)-6xi(α,β,γ).
Chua et al. have introduced reaction-diffusion CNN equations:(6)Ẋα,β,γ=f(Xα,β,γ)+D∇2Xα,β,γ,
where D=diag(D1,D2,…Dn), Ẋα,β,γ 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 hasn state variables, but only m (m≤n) state variables couple directly to their nearest neighbors via “reaction-diffusion”. Consequently, each cell has the following state equations:(7)V̇a=fa(Va,Vb)+Ia,V̇b=fa(Va,Vb),
where(8)Va=[V1,V2,…,Vm]T,Vb=[Vm+1,Vm+2,…,Vn]T,fa=[f1(Va,Vb),f2(Va,Vb),…,fm(Va,Vb)]T,fb=[fm+1(Va,Vb),fm+2(Va,Vb),…,fn(Va,Vb)]T,Ia=[D1∇2V1,D2∇2V2,…,Dm∇2Vm]T.The cell equilibrium pointQi=(Vai,Vbi)(∈Rn) of equation (7) can be determined via(9)fa(Va,Vb)=0,fb(Va,Vb)=0.The Jacobian matrix at the equilibrium pointQi has the following form:(10)J(Qi)=[akl(Qi)]n×n=[Aaa(Qi)Aab(Qi)Aba(Qi)Abb(Qi)],
where Akl(Qi) are called cell parameters and(11)Aaa(Qi)=[∂f1∂V1⋯∂f1∂Vm⋮⋱⋮∂fm∂x1⋯∂fm∂xm],Aab(Qi)=[∂f1∂Vm+1⋯∂f1∂Vn⋮⋱⋮∂fm∂xm+1⋯∂fm∂xn],Aba(Qi)=[∂fm+1∂V1⋯∂fm+1∂Vm⋮⋱⋮∂fn∂x1⋯∂fn∂xm],Abb(Qi)=[∂fm+1∂Vm+1⋯∂fm+1∂Vn⋮⋱⋮∂fn∂xm+1⋯∂fn∂xn].
The local state equations at the cell equilibrium point Qi are defined via(12)V̇a=AaaVa+AabVb+Ia,V̇b=AbaVa+AbbVb.Definition 1.
(13)YQ(s)=sI-Aaa-Aab(sI-Abb)-1Aba
is called the admittance matrix at the cell equilibrium point Qi.Lemma 2.
A reaction-diffusion CNN cell is called locally active at the equilibrium pointQi if and only if, its admittance matrix at Qi satisfies at least one of the following four conditions [4].(1)
YQ(s) has a pole in Re[s]>0.(2)
YQH(iω)=Y¯Q(iω)+YQ(iω)<0 for some ω=ω0, where ω0 is any real number.(3)
YQ(s) has a simple pole s=iωp on the imaginary axis, where its associated residue matrix:
(14)k1={lim⁡s→iωp(s-iωp)YQ(s),ifωp<∞lim⁡ωp→∞⁡YQ(iωp)iωp,ifωp=∞
is either a complex number or a negative real number.(4)
YQ(s) has a multiple pole on the imaginary axis.Definition 3.
The cell equilibrium pointQi is called stable if and only if, all the real parts of eigenvalue λi of Jacobianmatrix at the equilibrium point Qi are negative [2].Definition 4.
A “reaction-diffusion” CNN withn state variables and m ports is said to be operating on the “edge of chaos” with respect to an equilibrium point Qi if and only if, Qi is both locally active and stable when Ia=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(15)V̇a=AaaVa+AabVb+Ia,V̇b=AbaVa+AbbVb,
where(16)Va=[V1],Vb=[V2V3V4V5]T,Ia=[I1],Aaa=[a11],Aab=[a12a13a14a15],Aba=[a21a31a41a51]T,Aba=[a21a31a41a51],Abb=[a22a23a24a25a32a33a34a35a42a43a44a45a52a53a54a55].The corresponding CNN cell admittance matrixYQ(s) is given by [1].(17)YQ(s)=sI-Aaa-Aab(sI-Abb)-1Aba=s-a11-T1s3+K1s2+L1s+Δ1s4+Ts3+Ks2+Ls+Δ,
where T,T1,K,K1,L,L1,Δ,Δ1 are the parameters of aij’s.Theorem 5.
A necessary and sufficient condition forYQ(s) to satisfy condition (1) in Lemma 2 is that ∃s, such that g(s)=0(Re[s]>0), and any one of the following conditions holds.(1)
f(s)≠0.(2)
f(s)=0, and m>n, where s is m and n orders zero point of g(s) and f(s), respectively, where f(s)=T1s3+K1s2+L1s+Δ1,g(s)=s4+Ts3+Ks2+Ls+Δ.Proof.
Obviously proved.Denote(18)E=-a11,F=-TT1+K1-a11(-2K+T2),P=LL1-KΔ1-ΔK1,Q=Δ1+KK1-TL1-LT1,G=-a11(2Δ+K2-2LT)-Q,H=-a11(L2-2ΔK)-P,I=-a11Δ2-ΔΔ1,g(Q)=EQ4+FQ3+GQ2+HQ+I,h(λ)=-(TT1-K1)λ3-Qλ2-Pλ-ΔΔ1,λ1,2*=-Q±Q2-3(TT1-K1)P3(TT1-K1),p=-3F216E2+G2E,q=F332E3-FG8E2+H4E,w1,2=-1±i32,D=q24+p327,Aj=(-q2±D1/2)1/3,x1=A1+A2,x2=w1A1+w2A2,x3=w2A1+w1A2,Ωi=xj-F4E,i=1,2,3.Theorem 6.
Let the following parameters be defined as in Theorem5, then YQH(iw)<0 for some w=w0∈R if any one of the following conditions holds.(1)
a11>0.(2)
a11=0,TT1-K1>0.(3)
a11=0,TT1-K1=0,Q>0.(4)
a11=0,TT1-K1=0,Q<0,ΔΔ1>0.(5)
a11=0,TT1-K1=0,Q<0,P≥0,ΔΔ1-P2/Q/4>0,ΔΔ1≤0,(6)
a11=0,TT1-K1=0,Q=0,P>0.(7)
a11=0,TT1-K1=0,Q=0,P≤0,ΔΔ1>0.(8)
a11=0,TT1-K1<0,ΔΔ1>0.(9)
a11=0,TT1-K1<0,ΔΔ1≤0,and λj*≥0,h(λj*)<0, for j=1 or 2.(10)
a11<0,D>0,Ω1>0,g(Ω1)<0.(11)
a11<0,D<0, and Ωj≥0,g(Ωj)<0, for j=1,2 or 3.(12)
a11<0,D=0,p=q=0,g(-F/4E)<0.(13)
a11<0,D=0,q2/4=-p3/27≠0,and Ωj≥0, g(Ωj)<0, for j=1 or 2.Proof.
YQH(iω)=Y¯Q(iω)+YQ(iω)=2Re[YQ(iω)], so YQ(iω) to satisfy condition (2) in Lemma 2 equals to Re[YQ(iω)]<0,
(19)Re[YQ(iω)]=Re[T1(iω)3+K1(iω)2+L1(iω)+Δ1(iω)4+T(iω)3+K(iω)2+L(iω)+Δiω-a11-T1(iω)3+K1(iω)2+L1(iω)+Δ1(iω)4+T(iω)3+K(iω)2+L(iω)+Δ]=Eω8+Fω6+Gω4+Hω2+I(ω4-Kω2+Δ)2+(Lω-Tω2).(1)
Ifa11>0, then Re[YQ(iω)]<0 when ω is large enough (See (1) of Theorem 6).(2)
Ifa11=0, then
(20)Re[YQ(iω)]=-(TT1-K1)ω6+Qω4+Pω2+ΔΔ1(ω4-Kω2+Δ)2+(Lω-Tω2).Let f(λ)=-Qλ2-Pλ-ΔΔ1,(I)
IfTT1-K1>0, then Re[YQ(iω)]<0 when ω is large enough (See (2) of Theorem 6).(II)
IfTT1-K1=0, then
(i)
IfQ>0, then Re[YQ(iω)]<0 when ω is large enough (See (3) of Theorem 6).(ii)
IfQ<0,
(a)
IfΔΔ1>0,∃ω0∈R, such that Re[YQ(iω0)]<0 (See (4) of Theorem 6).(b)
IfΔΔ1≤0, solve f′(λ*)=0, we can get λ*=-0.5P/Q,f(λ*)=0.25P2/Q-ΔΔ1.Then whenP≥0,ΔΔ1-0.25P2/Q>0,∃ω0>λ*, such that Re[YQ(iω0)]<0 (See (5) of Theorem 6).(iii)
IfQ=0, then Re[YQ(iω)]=-Pω2-ΔΔ1.
(a)
IfP>0, then Re[YQ(iω)]<0 when ω is large enough (See (6) of Theorem 6).(b)
IfP≤0,ΔΔ1>0, then ∃ω0, such that Re[YQ(iω0)]=-ΔΔ1<0 (See (7) of Theorem 6).(III)
IfTT1-K1<0, let h(λ)=-(TT1-K1)λ3-Qλ2-Pλ-ΔΔ1,
(i)
IfΔΔ1>0, then ∃ω0, such that Re[YQ(iω0)]<0 (See (8) of Theorem 6).(ii)
IfΔΔ1≤0, solve h(λ*)=0, we can get
(21)λ1,2*=-Q±Q2-3(TT1-K1)P3(TT1-K1).
Then, for i=1,2, if λi*≥0, h(λi*)<0, then ∃ω0, such that Re[YQ(iω0)]<0 (See (9) of Theorem 6).(3)
Ifa11<0, let g(Q)=EQ4+FQ3+GQ2+HQ+I, then g′(Q)=4EQ3+3FQ2+2GQ+H. Let x=Ω+(F/4E), then the above becomes g′(Q)=4E(x3+px+q)=4Ef(x), then xi,i=1,2,3 are the roots of f(x)=0, Ωi are the roots of g′(Ω)=0. If any one of the (10)–(13) of Theorem 6 holds, we can get Re[YQ(iω0)]<0.So, if any one of conditions (1)–(13) holds, Re[YQ(iω0)]<0. YQ(s) Satisfies condition (2) in Lemma 2. This completes the proof.Theorem 7.
Forj=1, or 2, let
(22)wj=(K+2Δ+(-1)jK-2Δ)2,Aj=L-3wj*2,wj*=(K+(-1)jK2-4Δ)2,Bj=2Kwj*-4wj*3,A1j=Δ1-K1wj*2,B1j=L1wj*-T1wj*3.
Then YQ(s) satisfies condition (3) of Lemma 2, if any one of the following conditions holds.(I)
Δ>0,K>2Δ,T=L=0, and any one of the following conditions holds.
(1)
K1w12-Δ1≠0.(2)
K1w12-Δ1=0,(L1-T1w12)(w22-w12)>0.(3)
K1w22-Δ1≠0.(4)
K1w22-Δ1=0,(L1-T1w22)(w12-w22)>0.(II)
K>0,Δ1≠0,Δ=0,L=KT≠0, and any one of the following conditions holds.
(1)
TΔ1<0.(2)
T(K-T1K)-Δ1+KK1≠0.(3)
T(K-T1K)-Δ1+KK1=0,T(Δ1-KK1)+K(L1-T1K)<0.(III)
Δ=0,Δ1L>0.(IV)
Δ<0, orK>0,K2-4Δ>0,and 2L=T(K+K2-4Δ)(K+K2-4Δ>0) or 2L=T(K-K2-4Δ),(K-K2-4Δ>0) and any one of the following conditions holds for j=1, or 2.
(1)
AjB1j-A1jBj≠0.(2)
AjB1j-A1jBj=0,AjA1j-BjB1j>0.Proof.
Letf(s)=T1s3+K1s2+L1s+Δ1,g(s)=s4+Ts3+Ks2+Ls+Δ, obviously, ∞ is not a single pole of YQ(s) on the imaginary axis.
IfYQ(s) has a simple pole s=iω on the imaginary axis, where its associated residue
(23)k1=lim⁡s→iω(s-iω)YQ(s)=lim⁡s→iω(s-iω)f(s)g(s)=f(s)g′(s)|s=iωis either a complex number or a negative real number, then k1≠0, so f(s)≠0, which implies that iω is not a zero point of f(s)=0,iω  is not a removed pole of YQ(s).(I)
IfYQ(s) has four poles s=±iω1,±iω2(ω1≠ω2≠0) on the imaginary axis. In this case,g(s)=(s2+ω12)(s2+ω22)=s4+(ω12+ω22)s2+ω12ω22. Hence we obtain T=L=0, K=ω12+ω22>0, Δ=ω12ω22>0. Then, we can get K+2Δ=(ω1+ω2)2, K-2Δ=(ω1-ω2)2, which implies that K>2Δ, ω1,2=(K+2Δ±K-2Δ)/2. So,
(24)lim⁡s→±iω1⁡(s∓iω1)YQ(s)=w1(L1-T1w12)±i(K1w12-Δ1)2w1(w22-w12),(25)lim⁡s→±iω2⁡(s∓iω2)YQ(s)=w2(L1-T1w22)±i(K1w22-Δ1)2w2(w12-w22).
Then, when condition (I) in Theorem 7 holds, k1 is a complex number or a negative real number. YQ(s) satisfies condition (3) in Lemma 2.(II)
IfYQ(s) has a simple pole s=0 and two conjugate poles ±iω(ω≠0) on the imaginary axis, and another pole is a≠0.In this case, it follows thatΔ=0,Δ1≠0, and g(s) has the form:
(26)g(s)=s(s2+ω2)(s-a)=s4-as3+ω2s2-aω2,
which implies that T=-a,K=ω2>0,L=-aω2=KT,Δ=0,Δ1=0, Therefore,
(27)g(s)=s(s2+K)(s+T)
(1)
The residue ofYQ(s) at s=0 is
(28)lim⁡s→0sYQ(s)=Δ1(s2+K)(s+T)=Δ1KT.
Then, we conclude that if K>0,Δ1≠0,Δ=0,L=KT≠0,TΔ1<0, k1 is a negative real number. YQ(s) satisfies condition (3) in Lemma 2 (See (1) of (II) in Theorem 7).(2)
The residue ofYQ(s) at s=±iK is(29)lim⁡s→±iK(s∓iK)YQ(s)=T(Δ1-KK1)+K(L1-T1K)∓iK(T(L1-T1K)-Δ1+KK1)2K(K+T2).Consequently, we conclude that if (2) or (3) in (II) in Theorem 7 holds, k1 is either an imaginary number or a negative real number. YQ(s) satisfies condition (3) in Lemma 2.(III)
IfYQ(s) has a simple pole s=0 on the imaginary axis, and the other poles are ai,Re[ai]≠0,i=1,2,3, it follows that Δ1≠0,Δ=0, and g(s) has the form
(30)g(s)=s(s-a1)(s-a2)(s-a3)=s(s3-(a1+a2+a3)s2+(a1a2+a1a3+a2a3)s-a1a2a3).
Therefore we obtain that  Δ=0,T=-(a1+a2+a3),K=a1a2+a1a3+a2a3,L=-a1a2a3≠0, hence the reside of YQ(s) at s=0 is
(31)lim⁡s→0sYQ(s)=Δ1a1a2a3=-Δ1L.
Then, when Δ=0,Δ1L>0, k1 is a negative real number. YQ(s) satisfies condition (3) in Lemma 2 (See (III) of Theorem 7).(IV)
IfYQ(s) has two conjugate poles ±iω(ω>0) on the imaginary axis, and the other poles are Re[a]≠0,Re[b]≠0. In this case, g(s) has the form
(32)g(s)=(s-a)(s-b)(s+ω2)=s4-(a+b)s3+(ab+ω2)s2-(a+b)ω2s+abω2.
Therefore, we obtain that T=-(a+b),K=ab+ω2,L=-(a+b)ω2,Δ=abω2≠0. Then, ab=K-ω2,Δ=(K-ω2)ω2⇔ω4-Kω2+Δ=0. Solving it, we have
(33)ωj*=K±K2-4Δ2,j=1,2.
It implies that Δ<0 or K>0,K2-4Δ≥0, and T=-(a+b)=L/ω2. Then, the residue of YQ(s) at s=±ωj* is
(34)lim⁡s→±iωj*⁡(s∓iωj*)YQ(s)=-(A1jAj+B1jBj)±i(A1jBj-AjB1j)Aj2+Bj2.
Hence, if condition (IV) in Theorem 7 holds, k1<0 or it is an imaginary number. YQ(s) satisfies condition (3) in Lemma 2.
Therefore, when any one of conditions (I)–(IV) holds, YQ(s) satisfies condition (3) in Lemma 2. This completes the proof.Theorem 8.
YQ(s) has a multiple pole on the imaginary axis if any one of the following conditions holds.(I)
Δ=L=0,Δ1≠0.(II)
Δ=L=K=0,andΔ1≠0orL1≠0.(III)
Δ=L=K=T=0, and Δ1≠0 or L1≠0 or K1≠0.(IV)
T=L=0,K>0,Δ=(K/2)2, and any one of the following conditions holds.
(1)
2Δ1≠KK1.(2)
2L1≠KT1.Proof.
Letf(s)=T1s3+K1s2+L1s+Δ1,g(s)=s4+Ts3+Ks2+Ls+Δ. Obviously, when the conditions (I)–(III) hold, 0 is the multiply poles of YQ(s).
IfYQ(s) has two multiply nonzero poles ±iω(ω>0), then g(s) has the form:
(35)g(s)=(s2+ω2)2=s4+2ω2s2+ω4,
which implies that T=L=0,K=2ω2>0,Δ=ω4=(K/2)2>0.
If±iω are the multiply poles of YQ(s), thenf(iω)≠0 where
(36)f(±iω)=(Δ1-K1ω2)±iω(L1-T1ω2)=(Δ1-KK12)±iω(L1-T1K2).
Obviously, when any one of (1)-(2) of (IV) in Theorem 8 holds, f(iω)≠0, any one of ±iω is a multiply pole of YQ(s).
So, when any one of condition (I)–(IV) holds, YQ(s) satisfies condition (4) in Lemma 2. This completes the proof.
When any one of Theorems5–8 holds, YQ(s) 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]:
(37)dxdt=λ-dx-bvx-bnxvn,dydt=b(1-e)vx-ay,dvdt=ky-uv,dyndt=bevx+bnxvn-ayn,dvndt=knyn-uvn,
where the five variables—x,y,v,yn,vn 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 dx and become infected at rate bxv by normal virus and infected at rate bnxvn by mutated virus. Infected cells infected by normal and mutated virus are removed at rate ay and ayn, respectively. Normal virus is produced at rate ky and removed at rate uv, mutated virus is produced at rate kyn and removed at rate uvn. e is the rate constant describing the probability of mutation of virus (usual 10-5–10-3), a,b,bn,d,e,k,kn,u,λ 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:(38)dxijdt=λ-dxij-bxijvij-bnxijvnij+D1∇2xij,dyijdt=b(1-e)xijvij-ayij,dvijdt=kyij-uvij,dynijdt=bexijvij+bnxijvnij-aynij,dvnijdt=knynij-uvnij,
where ∇2xij=xi-1j+xi+1j+xij-1+xij+1-4xij.Let equation (38) be zeros (D1=0) and solve it, we can get the two equilibrium points:(39)Q1=(λd,0,0,0,0),(40)Q2=(x0,u(au-bnknux0)(λ-dx0)bkx0(bnknex0+au-knx0),(au-bnknux0)(λ-dx0)bx0(bnknex0+au-knx0),eu(λ-dx0)bnknex0+au-knx0,kne(λ-dx0)bnknex0+au-knx0),
where x0=au/((1-e)bk) and Q1, Q2 stand for the patient’s complete recovery and HBV persistent infection, respectively.Consequently, the Jacobian matrix at the equilibrium pointQi(i=1,2) is(41)J(Qi)=[-d-bv-bnvn0-bx0-bnxb(1-e)v-ab(1-e)x000k-u00bev+bnvn0bex-abnx000kn-u].Takingk,u as variables, and λ=10, a=0.5, b=0.01, bn=0.005, e=0.0001, kn=10, and d=0.01, using Theorems 5–8, we can calculate the bifurcation of the reaction-diffusion CNN model equation (38) at the equilibrium points Q1 and Q2 at k∈[0,40],u∈[0,10], see Figures 1 and 2.Figure 1
Bifurcation diagrams of equation (38) at the equilibrium points Q1 at k∈[0,40],u∈[0,10].Figure 2
Bifurcation diagrams of equation (38) at the equilibrium points Q2 at k∈[0,40],u∈[0,10].In Figures1 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 Q1 does not exist at the edge of chaos domain.Takeλ=10, k=0.01, a=0.5, b=0.01, bn=0.005, kn=10,e=0.0001, andk=1.0,3.0,4.9,5.1,10,24,39, u=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.
No.ukEquilibrium pointEigenvaluesDynamic pattern121.020,0,0,19,9869.4396,29.8638,-0.0090,-33.3646,-71.9897Convergent, divergent223.020,0,0,20,98-0.0097,53.0149,69.4400,-56.5150,-71.9902Convergent, divergent324.920,0,0,20,98-0.0098,68.2343,69.4485,-71.6962,-72.0368Convergent, divergent425.120,20,50,0,0-0.2043±0.3878i,-0.0000,-2.6014,-2.5000Convergent521010,20,99,0,0-2.7130,-0.3935±0.4583i,-0.2192,-2.2808Convergent62244,20,239,0,0-3.2812,-0.8094±0.2709i,-0.3768,-2.1232Convergent72393,20,389,0,0-4.4393,-1.2723,-0.6884,-0.4059,-2.0941Convergent851.050,0,0,19,3867.9709,28.4103,-0.0075,-34.9127,-73.5210Convergent, divergent953.050,0,0,19,38-0.0092,51.5420,67.9713,-58.0426,-73.5215Convergent, divergent1054.950,0,0,19,38-0.0095,66.7552,67.9801,-73.2207,-73.5651Convergent, divergent1155.149,19,19,0,0-0.0920±0.2787i,-0.0092,-5.5160,-5.4909Convergent1251025,19,39,0,0-0.1829±0.3778i,-0.2375,-5.5343,-5.2625Convergent1352410,20,95,0,0-5.5708,-0.4446±0.4784i,-0.3915,-5.1085Convergent145396,20,155,0,0-5.6298,-0.7151±0.4210i,-0.4343,-5.0657Convergent1591.090,0,0,18,2066.0620,26.5823,-0.0055,-37.0867,-75.6121Convergent, divergent1693.090,0,0,18,20-0.0085,49.6418,66.0624,-60.1432,-75.6126Convergent, divergent1794.990,0,0,18,20-0.0091,64.8330,66.0717,-75.3029,-75.6527Convergent, divergent1895.188,18,10,0,0-0.0531±0.2110i,-0.0106,-9.5037,-9.4894Convergent1991045,19,21,0,0-0.1047±0.2973i,-0.2431,-9.5106,-9.2569Convergent2092419,20,52,0,0-0.2481±0.4288i,-9.5339,-0.3897,-9.1103Convergent2193912,20,86,0,0-0.4014±0.4931i,-9.5671,-0.4300,-9.0700ConvergentIn the following discussions, we select some parameters in No. 8-14 andu=5, k=12.5. The simulation results are shown in Figures 3, 4, 5, 6, and 7. During the simulation, we reached a new conclusion.The kinetic trajectories of equation (37) when u=5, k=3.
(a)(b)(c)(d)The kinetic trajectories of equation (37) when u=5,k=4.9.
(a)(b)(c)(d)The kinetic trajectories of equation (37) when u=5,k=5.1.
(a)(b)(c)(d)The kinetic trajectories of equation (37) when u=5,k=12.5.
(a)(b)(c)(d)The kinetic trajectories of equation (37) when u=5,k=24.
(a)(b)(c)(d)From Table1 and Figures 3–7, we can conclude that(I)
whenk is smaller (less than 5),
(1)
these parameters are located in the green domain (the local and unstable domains);(2)
regardless of the value ofu, 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)
Whenk 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 ofu, 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)
Whenk<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 theu 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.



---

*Source: 674243-2012-04-12.xml*