Dynamics and Synchronization of Weakly Coupled Memristive
  Reaction-Diffusion Neural Networks
A new mathematical model of memristive neural networks described by the
partly diffusive reaction-diffusion equations with weak synaptic coupling is
proposed and investigated. Under rather general conditions it is proved that
there exists an absorbing set showing the dissipative dynamics of the solution
semiflow in the energy space and multiple ultimate bounds. Through uniform
estimates and maneuver of integral inequalities and sharp interpolation
inequalities on the interneuron differencing equations, it is rigorously proved
that exponential synchronization of the neural network solutions at a uniform
convergence rate occurs if the coupling strength satisfies a threshold
condition expressed by the system parameters. Applications with numerical
simulation to the memristive diffusive Hindmarsh-Rose neural networks and
FitzHugh-Nagumo neural networks are also shown.
Yuncheng You, Junyi Tu
# Dynamics and synchronization of weakly coupled memristive reaction-diffusion neural networks

###### Abstract.

A new mathematical model of memristive neural networks described by the partly diffusive reaction-diffusion equations with weak synaptic coupling is proposed and investigated. Under rather general conditions it is proved that there exists an absorbing set showing the dissipative dynamics of the solution semiflow in the energy space and multiple ultimate bounds. Through uniform estimates and maneuver of integral inequalities and sharp interpolation inequalities on the interneuron differencing equations, it is rigorously proved that exponential synchronization of the neural network solutions at a uniform convergence rate occurs if the coupling strength satisfies a threshold condition expressed by the system parameters. Applications with numerical simulation to the memristive diffusive Hindmarsh-Rose neural networks and FitzHugh-Nagumo neural networks are also shown.

Key words and phrases:Memristive neural network, synchronization, reaction-diffusion equations, dissipative dynamics, nonlinear coupling strength 2010 Mathematics Subject Classification: 35B40, 35G50, 35K57, 37N25, 92C20

## 1. Introduction

In this paper we shall consider a new mathematical model of neural networks described by a system of partly diffusive hybrid differential equations with memristors and weak nonlinear interneuron coupling.

Let a network of \(m\) coupled memristive neuron cells be denoted by \(\mathcal{NW}=\{\mathcal{N}_{i}:i=1,2,\cdots,m\}\), where \(m\geq 2\) is a positive integer, which is described by the following model of partly diffusive and memristive equations with nonlinear couplings. Each neuron \(\mathcal{N}_{i},1\leq i\leq m\), in the network is presented by the hybrid differential equations:

\[\begin{split}&\frac{\partial u_{i}}{\partial t}=\eta\Delta u_{i}+f(u_{ i},z_{i})-k\tanh(\rho_{i})u_{i}-Pu_{i}\sum_{j=1}^{m}\Gamma(u_{j}),\\ &\frac{\partial z_{i}}{\partial t}=\Lambda z_{i}+h(u_{i},z_{i}),\\ &\frac{\partial\rho_{i}}{\partial t}=au_{i}-b\rho_{i},\end{split} \tag{1.1}\]
for \(t>0,\ x\in\Omega\subset\mathbb{R}^{n}\ (n\leq 3)\), where \(\Omega\) is a bounded domain with locally Lipschitz continuous boundary \(\partial\Omega\). Here \(\Delta\) is the Laplacian operator with respect to the spatial variable \(x\in\Omega\) and \(\Lambda\) is an \(\ell\times\ell\) constant square matrix. In the nonlinear synaptic coupling term for neurons, the function

\[\Gamma(s)=\frac{1}{1+\exp[-r(s-V)]} \tag{1.2}\]

is a sigmoidal function. The biological meaning of the parameters in this type of nonlinear weak couplings can be seen in [9, 16, 24, 28, 37]. The coefficient \(P>0\) is the coupling strength. The constant \(V\in\mathbb{R}\) is a threshold for neuron bursting, and \(r>0\) shapes the sigmoidal function versus the Heaviside function.

In this system (1.1), for \(1\leq i\leq m\), the transmembrane electric potential \(u_{i}(t,x)\) and the memductance \(\rho_{i}(t,x)\) of the memristor (which is caused by the electromagnetic induction flux across the neuron membrane) are scalar functions, while \(z_{i}(t,x)\) can be an \(\ell\)-dimensional (\(\ell\geq 1\)) vector function whose components represent various ionic currents in the neuron cell. The memristive-potential coupling \(-k\tanh(\rho_{i})u_{i}\) is a nonlinear term.

We impose the homogeneous Neumann boundary condition on the first component function in (1.1):

\[\frac{\partial u_{i}}{\partial\nu}(t,x)=0,\quad\text{for}\ \ t>0,\ x\in \partial\Omega,\quad 1\leq i\leq m. \tag{1.3}\]

The initial states of the system (1.1) will be denoted by

\[u_{i}^{0}(x)=u_{i}(0,x),\ z_{i}^{0}(x)=z_{i}(0,x),\ \rho_{i}^{0}=\rho_{i}(0,x), \ 1\leq i\leq m. \tag{1.4}\]

The scalar function \(f\in C^{1}(\mathbb{R}^{1+\ell},\mathbb{R})\) and the vector function \(h\in C^{1}(\mathbb{R}^{1+\ell},\mathbb{R}^{\ell})\) are continuously differentiable and we assume that

\[\begin{split}& f(s,\sigma)s\leq-\alpha|s|^{4}+\lambda|s|| \sigma|+J,\quad(s,\sigma)\in\mathbb{R}^{1+\ell},\\ &\max\left\{\frac{\partial f}{\partial s}(s,\sigma),\,\left| \frac{\partial f}{\partial\sigma}(s,\sigma)\right|\right\}\leq\beta,\ (s,\sigma)\in\mathbb{R}^{1+\ell},\end{split} \tag{1.5}\]

and the \(\ell\)-dimensional matrix \(\Lambda\) and \(h(s,\sigma)\) satisfy

\[\begin{split}&\langle\Lambda\sigma,\sigma\rangle\leq-\gamma| \sigma|^{2},\quad\sigma\in\mathbb{R}^{\ell},\\ & h(s,\sigma)\sigma\leq q|s|^{2}|\sigma|+L|\sigma|,\quad(s, \sigma)\in\mathbb{R}^{1+\ell},\\ &\left|\frac{\partial h}{\partial s}(s,\sigma)\right|\leq\xi(|s |+1),\ \ \frac{\partial h}{\partial\sigma}(s,\sigma)=0,\quad(s,\sigma)\in\mathbb{R}^{1+ \ell},\end{split} \tag{1.6}\]

where the parameters \(\eta,a,b,k\) in (1.1), \(\alpha,\lambda,J,\beta\) in (1.5), and \(\gamma,q,L,\xi\) in (1.6) are all positive constants. Note that Assumptions (1.5) and (1.6) are satisfied by the
partly diffusive FitzHugh-Nagumo neural networks [33, 34] and the partly diffusive Hindmarsh-Rose neural networks [45, 46], which will be shown in Section 6.

In neurobiology, the ionic currents flowing across a neuron cell's membrane cause changes of the membrane potential over time. The resulting electrical signals propagate through the neuron axon and stimulate the dendrites of neighbor neurons by synapses, which constitute a biological neural network. An excitatory neuron's firing consists of successive spiking followed by relatively long period of quiescence called bursting. Chaotic bursting means the number of spikes per burst is irregular.

Synchronization mechanism of various neural networks revealed by mathematical models and analysis is one of the central topics in the research of neuroscience, medical science, and artificial neural networks [1, 3, 10, 13, 15, 29, 39, 41, 43]. For many neuron models in terms of ODEs with or without time-delay, analysis of Hopf bifurcations, stability by Lyapunov exponents, the energy or Hamiltonian functions, and numerical simulations are the main approaches to show synchronization of neuron ensembles or neural networks [10, 18, 25, 28, 37, 48].

Dynamics and synchronization of cellular neural networks and memristive neural networks modeled by partly diffusive Hindmarsh-Rose equations and FitzHugh-Nagumo equations have been studied in the authors' group recently [26, 27, 33, 34, 45, 46, 47]. These neural network models of hybrid (PDE-ODE) differential equations reflect the structural feature of neuron cells, which contain the short-branch dendrites receiving incoming signals and the long-branch axon transmitting outgoing signals through synapses.

Memristor concept coined by Chua [7] describes the effect of electromagnetic flux on moving charges such as the ionic currents in neuron cells. In advanced biological neuron models and artificial intelligence computing [2, 11, 19, 22, 36, 38, 39, 42], the memristive feature was exhibited as a new type (other than electrical and chemical) synaptic coupling or an ideal component which has the nonvolatile properties and can process dynamically memorized signal information to deal with complex or chaotic behaviors in neural networks. Memristor-based differential equation models now appear in many fields with applications to image encryption, DNA sequences operation, brain criticality, cell physiology, cybersecurity, drift-diffusion models in semiconductor devices, and quantum computers, cf. [20, 22, 23, 31, 35, 36, 40].

Researches on dynamics of memristive neural networks in ODE models with linear interneuron couplings are increasing in the recent years, cf.[2, 14, 11, 21, 25, 30, 37, 38, 39, 44], mainly by computational simulations combined with semi-analytic methods.

Very recently in [45, 46, 47] the authors rigorously proved the dissipative dynamics and the exponential synchronization of the memristive Hindmarsh-Rose neural networks and FitzHugh-Nagumo neural networks with the partly diffusive PDE-ODE models and linear interneuron couplings. Note that the linear interneuron coupling
can be viewed as a strong coupling for neural networks, which is mathematically amenable but may not always well reflect the biological synaptic interactions.

It is an open and challenging problem to theoretically show that partly diffusive and memristive biological neural networks with the nonlinear interneuron coupling shown in (1.1) with (1.2) is dissipative and synchronizable under a threshold condition. Such a nonlinear coupling can be called weak coupling.

In this work we shall prove a sufficient condition on the network coupling strength to ensure an exponential synchronization of the neural networks modeled in (1.1)-(1.2) by the approach of showing dissipative dynamics of the system through sharp and uniform integral estimates. Moreover, the general model (1.1) with the Assumptions (1.5)-(1.6) for neural networks can cover all the typical neuron models including Hindmarsh-Rose equations [16], FitzHugh-Nagumo equations [12], Hodgkin-Huxley equations [17], FitzHugh-Rinzel equations [42], and Morris-Lecar equations [41]. It is worth mentioning that an effective methodology developed here from dissipative dynamics toward synchronization was originated in J.K. Hale's paper in 1997 [15]. This approach can be extended and explored to study many other complex or artificial neural networks in a broad scope.

## 2. **Formulation and Preliminaries**

For a framework to formulate the solution and dynamics problem of the neural network system (1.1)-(1.4), we define two Hilbert spaces of functions:

\[E=[L^{2}(\Omega,\mathbb{R}^{2+\ell})]^{m}\quad\text{and}\quad\Pi=[H^{1}(\Omega )\times L^{2}(\Omega,\mathbb{R}^{\ell+1})]^{m}\]

where \(H^{1}(\Omega)\) is a Sobolev space. One can call \(E\) the energy space and \(\Pi\) the regular space. The norm and inner-product of \(L^{2}(\Omega)\) or \(E\) will be denoted by \(\|\,\cdot\,\|\) and \(\langle\,\cdot,\cdot\,\rangle\), respectively. We use \(|\,\cdot\,|\) to denote a vector norm or a set measure in Euclidean spaces \(\mathbb{R}^{n}\).

The initial-boundary value problem (1.1)-(1.4) can be formulated into an initial value problem of the evolutionary equation:

\[\begin{split}\frac{\partial g}{\partial t}=&\,Ag+F( g),\,\,\,t>0,\\ &\,g(0)=g^{0}\in E.\end{split} \tag{2.1}\]

The unknown function in (2.1) is a column vector \(g(t)=\text{col}\,\,(g_{1}(t),g_{2}(t),\cdots,g_{m}(t))\), where the component subvector

\[g_{i}(t)=\text{col}\,(u_{i}(t,\cdot),\,z_{i}(t,\cdot),\,\rho_{i}(t,\cdot)), \quad\text{for}\,\,1\leq i\leq m,\]

characterizes the dynamics of the neuron \(\mathcal{N}_{i}\), for \(1\leq i\leq m\). The initial data function in (2.1) is

\[g(0)=g^{0}=\text{col}\,\,(g_{1}^{0},\,g_{2}^{0},\cdots,g_{m}^{0})\quad\text{ where}\,\,\,g_{i}^{0}=\text{col}\,(u_{i}^{0},\,z_{i}^{0},\,\rho_{i}^{0}),\,\,\,1\leq i \leq m.\]
The energy norm \(\|g(t)\|\) of a solution \(g(t)\) for the evolutionary equation (2.1) in the space \(E\) is given by

\[\|g(t)\|^{2}=\sum_{i=1}^{m}\|g_{i}(t)\|^{2}=\sum_{i=1}^{m}\left(\|u_{i}(t)\|^{2}+ \|z_{i}(t)\|^{2}+\|\rho_{i}(t)\|^{2}\right).\]

The closed linear operator \(A\) in (2.1) is defined by \(A=\operatorname{diag}\left(A_{1},A_{2},\cdots,A_{m}\right)\), where

\[A_{i}=\begin{pmatrix}\eta\Delta&0&0\\ 0&\Lambda&0\\ 0&0&-bI\end{pmatrix}_{(2+\ell)\,\times\,(2+\ell)}:\mathcal{D}(A)\to E,\quad i=1,2,\cdots,m, \tag{2.2}\]

with the domain \(\mathcal{D}(A)=\{g\in[H^{2}(\Omega)\times L^{2}(\Omega,\mathbb{R}^{\ell+1})]^ {m}:\partial u_{i}/\partial\nu=0,1\leq i\leq m\}\). The operator \(A\) is the generator of \(C_{0}\)-semigroup \(\{e^{At}\}_{t\geq 0}\) on the space \(E\) and \(I\) is the identity operator. By the fact that the Sobolev imbedding \(H^{1}(\Omega)\hookrightarrow L^{6}(\Omega)\) is a continuous mapping for space dimension \(n\leq 3\) and according to Assumptions (1.5) and (1.6), the nonlinear mapping

\[F(g)=\begin{pmatrix}f(u_{1},z_{1})-k\tanh(\rho_{1})u_{1}-Pu_{1}\sum_{j=2}^{m} \Gamma(u_{j})\\ h(u_{1},z_{1})\\ au_{1}\\ \vdots\\ f(u_{m},z_{m})-k\tanh(\rho_{m})u_{m}-Pu_{m}\sum_{j=1}^{m-1}\Gamma(u_{j})\\ h(u_{m},z_{m}).\\ au_{m}\end{pmatrix}:\Pi\longrightarrow E \tag{2.3}\]

is a locally Lipschitz continuous mapping.

In this work we shall consider the weak solutions, cf. [6, Section XV.3] and [32, Section 4.2.3], of this initial value problem (2.1).

**Definition 2.1**.: A \((2+\ell)m\)-dimensional vector function \(g(t,x)\), where \((t,x)\in[0,\tau]\times\Omega\), is called a weak solution to the initial value problem of the evolutionary equation (2.1), if the following two conditions are satisfied:

(i) \(\frac{d}{dt}(g(t),\zeta)=(A^{1/2}g(t),A^{1/2}\zeta)+(F(g(t)),\zeta)\) is satisfied for almost every \(t\in[0,\tau]\) and any \(\zeta\in E^{*}=E\).

(ii) \(g(t,\cdot)\in C([0,\tau];E)\cap C^{1}((0,\tau);E)\) and \(g(0)=g^{0}\).

Here \(\mathcal{D}(A^{1/2})=\{g\in\Pi:\partial u_{i}/\partial\nu=0,\,1\leq i\leq m\}\) and \(E^{*}\) is the dual space of \(E\). The bilinear \(E\) vs \(E^{*}\) dual product is in the scalar distribution sense.

The following proposition can be proved by Galerkin spectral approximation method [6] for the first statement on weak solutions and by the compactness property of the
parabolic semigroup \(e^{At}\) in mild solution bootstrap argument [32, Theorem 42.12 and Corollary 42.13] for the second statement on strong solutions when \(t>0\).

**Proposition 2.2**.: _For any given initial state \(g^{0}\in E\), there exists a unique weak solution \(g(t;g^{0}),\,t\in[0,\tau]\), for some \(\tau>0\) may depending on \(g^{0}\), of the initial value problem (2.1) formulated from the memristive neural network equations (1.1). The weak solution \(g(t;g^{0})\) continuously depends on the initial data \(g^{0}\) and satisfies_

\[g\in C([0,\tau];E)\cap C^{1}((0,\tau);E)\cap L^{2}((0,\tau);\Pi). \tag{2.4}\]

_Moreover, for any initial state \(g^{0}\in E\), the weak solution \(g(t;g^{0})\) becomes a strong solution for \(t\in(0,\tau)\), which has the regularity_

\[g\in C((0,\tau];\Pi)\cap C^{1}((0,\tau);\Pi). \tag{2.5}\]

An infinite dimensional dynamical system [6, 32] for time \(t\geq 0\) only is usually called a semiflow. Absorbing set defined below is the key concept to characterize dissipative dynamics of a semiflow on a Banach space.

**Definition 2.3**.: Let \(\{S(t)\}_{t\geq 0}\) be a semiflow on a Banach space \(\mathscr{X}\). A bounded set \(B^{*}\) of \(\mathscr{X}\) is called an _absorbing set_ of this semiflow, if for any given bounded set \(B\subset\mathscr{X}\) there exists a finite time \(T_{B}\geq 0\) depending on \(B\), such that \(S(t)B\subset B^{*}\) for all \(t>T_{B}\). The semiflow is said to be _dissipative_ if there exists an absorbing set.

The Young's inequality in a generic form will be used throughout in this paper. For any two positive numbers \(x\) and \(y\), if \(\frac{1}{p}+\frac{1}{q}=1\) and \(p>1,q>1\), one has

\[x\,y\leq\frac{1}{p}\varepsilon x^{p}+\frac{1}{q}C(\varepsilon,p)\,y^{q}\leq \varepsilon x^{p}+C(\varepsilon,p)\,y^{q},\quad C(\varepsilon,p)=\varepsilon^{ -q/p}, \tag{2.6}\]

where constant \(\varepsilon>0\) can be arbitrarily small. The Gagliardo-Nirenberg interpolation inequalities [32, Theorem B.3] will be exploited in a crucial step to prove the main result on neural network synchronization.

## 3. **Dissipative Dynamics of the Memristive Semiflow**

In this section, first we shall prove the global existence of weak solutions in time for the initial value problem (2.1) and establish a solution semiflow of the memristive neural networks modeled by (1.1). Then we show the existence of an absorbing set of this semiflow in the state spaces \(E\), which exhibits the dissipative dynamics of this memristive neural network semiflow.

**Theorem 3.1**.: _Under the Assumption (1.5), for any initial state \(g^{0}\in E\), there exists a unique global weak solution in time, \(g(t;g^{0})=\operatorname{col}\left(u_{i}(t),v_{i}(t),w_{i}(t),\rho_{i}(t):1\leq i \leq m),t\in[0,\infty)\), to the initial value problem (2.1) of the memristive neural network equations (1.1)._
Proof.: Take the \(L^{2}\) inner-products of the \(u_{i}\)-equation in (1.1) with \(C_{1}u_{i}(t,x)\) for \(1\leq i\leq m\), with a scaling constant \(C_{1}>0\) to be chosen later. Then sum them up. By using the Gauss divergence theorem and the boundary condition (1.3) to treat the Laplacian term and by the Assumprion (1.5), we can get

\[\begin{split}&\frac{C_{1}}{2}\frac{d}{dt}\sum_{i=1}^{m}\|u_{i}(t)\|^ {2}+C_{1}\sum_{i=1}^{m}\eta\|\nabla u_{i}(t)\|^{2}\\ =&\,C_{1}\sum_{i=1}^{m}\int_{\Omega}\left[f(u_{i},z_ {i})u_{i}-k\tanh(\rho_{i})u_{i}^{2}-\sum_{j=1,j\neq i}^{m}\frac{Pu_{i}^{2}}{1+ \exp[-r(u_{j}-V)]}\right]dx\\ \leq&\,C_{1}\sum_{i=1}^{m}\int_{\Omega}\left[-\alpha |u_{i}|^{4}+\lambda|u_{i}||z_{i}|+J-k\tanh(\rho_{i})u_{i}^{2}-\sum_{j=1,j\neq i }^{m}\frac{Pu_{i}^{2}}{1+\exp[-r(u_{j}-V)]}\right]dx\\ \leq&\,C_{1}\sum_{i=1}^{m}\int_{\Omega}\left[-\alpha |u_{i}|^{4}+k|u_{i}|^{2}+\lambda|u_{i}||z_{i}|+J\right]dx\\ \leq&\,-C_{1}\alpha\sum_{i=1}^{m}\int_{\Omega}u_{i}^ {4}(t,x)\,dx+\left(C_{1}k+\frac{C_{1}^{2}\lambda^{2}}{\gamma}\right)\sum_{i=1 }^{m}\|u_{i}\|^{2}+\frac{\gamma}{4}\sum_{i=1}^{m}\|z_{i}\|^{2}+C_{1}mJ\,|\Omega |,\end{split} \tag{3.1}\]

because

\[-\sum_{j=1,j\neq i}^{m}\frac{Pu_{i}^{2}}{1+\exp[-r(u_{j}-V)]}\leq 0,\]

where the Young's inequality (2.6) and the property \(|\tanh(\rho_{i})|\leq 1\) are used. Then sum up the \(L^{2}\) inner-products of the \(z_{i}\)-equation with \(z_{i}(t,x)\) and the \(\rho_{i}\)-equation with \(\rho_{i}(t,x)\) in (1.1), \(1\leq i\leq m\), again by (2.6), we have

\[\begin{split}&\frac{1}{2}\frac{d}{dt}\sum_{i=1}^{m}\left(\|z_{i}(t)\|^ {2}+\|\rho_{i}(t)\|^{2}\right)=\sum_{i=1}^{m}\int_{\Omega}\left(\langle\Lambda z _{i},z_{i}\rangle+h(u_{i},z_{i})z_{i}+au_{i}\rho_{i}-b\rho_{i}^{2}\right)dx\\ \leq&\,\sum_{i=1}^{m}\int_{\Omega}\left[-\gamma|z_{i} |^{2}+q|u_{i}|^{2}|z_{i}|+L|z_{i}|+au_{i}\rho_{i}-b\,\rho_{i}^{2}\right]dx\\ \leq&\,\sum_{i=1}^{m}\frac{q^{2}}{4\gamma}\int_{\Omega }u_{i}^{4}(t,x)\,dx-\gamma\sum_{i=1}^{m}\|z_{i}\|^{2}+\frac{a^{2}}{2b}\sum_{i= 1}^{m}\|u_{i}\|^{2}-\frac{b}{2}\sum_{i=1}^{m}\|\rho_{i}\|^{2}+\frac{mL^{2}}{ \gamma}|\Omega|\\ =&\,\sum_{i=1}^{m}\frac{q^{2}}{4\gamma}\int_{\Omega }u_{i}^{4}(t,x)\,dx-\frac{3\gamma}{4}\sum_{i=1}^{m}\|z_{i}\|^{2}+\frac{a^{2}}{2 b}\sum_{i=1}^{m}\|u_{i}\|^{2}-\frac{b}{2}\sum_{i=1}^{m}\|\rho_{i}\|^{2}+\frac{mL^{2}}{ \gamma}|\Omega|.\end{split} \tag{3.2}\]
Both inequalities (3.1) and (3.2) are valid in the time interval \(I_{max}(g^{0})=[0,T_{max})\) of solution existence for each weak solution \(g(t;g^{0}))\).

Now we add the above two inequalities (3.1) and (3.2) to obtain

\[\begin{split}&\frac{1}{2}\frac{d}{dt}\sum_{i=1}^{m}\left(C_{1}\|u_{i}(t )\|^{2}+\|z_{i}(t)\|^{2}+\|\rho_{i}(t)\|^{2}\right)+C_{1}\eta\sum_{i=1}^{m}\| \nabla u_{i}(t)\|^{2}\\ \leq&-C_{1}\alpha\sum_{i=1}^{m}\int_{\Omega}u_{i}^{4} (t,x)\,dx+\left(C_{1}k+\frac{C_{1}^{2}\lambda^{2}}{\gamma}\right)\sum_{i=1}^{ m}\|u_{i}\|^{2}+\frac{\gamma}{4}\sum_{i=1}^{m}\|z_{i}\|^{2}+C_{1}mJ|\Omega|\\ &+\sum_{i=1}^{m}\frac{q^{2}}{4\gamma}\int_{\Omega}u_{i}^{4}(t,x) \,dx+\frac{a^{2}}{2b}\sum_{i=1}^{m}\|u_{i}\|^{2}-\frac{3\gamma}{4}\sum_{i=1}^{ m}\|z_{i}\|^{2}-\frac{b}{2}\|\rho_{i}\|^{2}+\frac{mL^{2}}{\gamma}|\Omega|\\ =&-\left(C_{1}\alpha-\frac{q^{2}}{4\gamma}\right) \sum_{i=1}^{m}\int_{\Omega}u_{i}^{4}(t,x)\,dx+\left(C_{1}k+\frac{C_{1}^{2} \lambda^{2}}{\gamma}+\frac{a^{2}}{2b}\right)\sum_{i=1}^{m}\int_{\Omega}u_{i}^ {2}(t,x)\,dx\\ &-\frac{\gamma}{2}\|z_{i}\|^{2}-\sum_{i=1}^{m}\frac{b}{2}\|\rho_{ i}\|^{2}+m\left(C_{1}J+\frac{L^{2}}{\gamma}\right)|\Omega|,\quad t\in I_{max}=[0,T_{ max}).\end{split} \tag{3.3}\]

Choose the scaling constant \(C_{1}\) to be

\[C_{1}=\frac{1}{\alpha}\left(1+\frac{q^{2}}{4\gamma}\right)\quad\text{so that}\quad C_{1}\alpha-\frac{q^{2}}{4\gamma}=1. \tag{3.4}\]

With this choice, from (3.3) it follows that

\[\begin{split}&\frac{1}{2}\frac{d}{dt}\sum_{i=1}^{m}\left(C_{1}\|u_{i} \|^{2}+\|z_{i}\|^{2}+\|\rho_{i}\|^{2}\right)+C_{1}\eta\sum_{i=1}^{m}\|\nabla u _{i}\|^{2}\\ &+\sum_{i=1}^{m}\int_{\Omega}u_{i}^{4}(t,x)\,dx-\left(C_{1}k+ \frac{C_{1}^{2}\lambda^{2}}{\gamma}+\frac{a^{2}}{2b}\right)\sum_{i=1}^{m}\int _{\Omega}u_{i}^{2}(t,x)\,dx\\ &+\frac{\gamma}{2}\sum_{i=1}^{m}\|z_{i}\|^{2}+\frac{b}{2}\sum_{i= 1}^{m}\|\rho_{i}\|^{2}\leq m\left(C_{1}J+\frac{L^{2}}{\gamma}\right)|\Omega|, \quad t\in I_{max}=[0,T_{max}).\end{split} \tag{3.5}\]
By completing square, we have

\[\begin{split}&\sum_{i=1}^{m}\int_{\Omega}u_{i}^{4}(t,x)\,dx-\left(C _{1}k+\frac{C_{1}^{2}\lambda^{2}}{\gamma}+\frac{a^{2}}{2b}\right)\sum_{i=1}^{m} \int_{\Omega}u_{i}^{2}(t,x)\,dx\\ =&\sum_{i=1}^{m}\int_{\Omega}\left(u_{i}^{4}(t,x)- \left(C_{1}k+\frac{C_{1}^{2}\lambda^{2}}{\gamma}+\frac{a^{2}}{2b}\right)u_{i}^ {2}(t,x)\right)dx\\ =&\sum_{i=1}^{m}\int_{\Omega}\left(u_{i}^{2}(t,x)- \frac{1}{2}\left(C_{1}k+\frac{C_{1}^{2}\lambda^{2}}{\gamma}+\frac{a^{2}}{2b}+1 \right)\right)^{2}dx\\ &+\sum_{i=1}^{m}\|u_{i}\|^{2}-\frac{m}{4}\left(C_{1}k+\frac{C_{1 }^{2}\lambda^{2}}{\gamma}+\frac{a^{2}}{2b}+1\right)^{2}|\Omega|\\ \geq&\sum_{i=1}^{m}\|u_{i}\|^{2}-\frac{m}{4}\left(C _{1}k+\frac{C_{1}^{2}\lambda^{2}}{\gamma}+\frac{a^{2}}{2b}+1\right)^{2}|\Omega |.\end{split} \tag{3.6}\]

Substitute (3.6) in (3.5). It yields the inequality

\[\begin{split}&\frac{1}{2}\frac{d}{dt}\sum_{i=1}^{m}\left(C_{1}\|u_{i} \|^{2}+\|z_{i}\|^{2}+\|\rho_{i}\|^{2}\right)+C_{1}\eta\sum_{i=1}^{m}\|\nabla u _{i}\|^{2}\\ &+\sum_{i=1}^{m}\left(\|u_{i}\|^{2}+\frac{\gamma}{2}\|z_{i}\|^{2} +\frac{b}{2}\|\rho_{i}\|^{2}\right)\leq C_{2}m|\Omega|,\quad t\in I_{max}. \end{split} \tag{3.7}\]

Denote by

\[C_{2}=C_{1}J+\frac{L^{2}}{\gamma}+\frac{1}{4}\left(C_{1}k+\frac{C_{1}^{2} \lambda^{2}}{\gamma}+\frac{a^{2}}{2b}+1\right)^{2}. \tag{3.8}\]

We can remove the nonnegative term \(C_{1}\eta\sum_{i=1}^{m}\|\nabla u_{i}\|^{2}\) from (3.7) to obtain the Gronwall-type differential inequality:

\[\begin{split}&\frac{d}{dt}\sum_{i=1}^{m}\left[C_{1}\|u_{i}\|^{2}+\|z_{i}\|^{2}+\| \rho_{i}\|^{2}\right]+\mu\sum_{i=1}^{m}\left[C_{1}\|u_{i}\|^{2}+\|z_{i}\|^{2}+ \|\rho_{i}\|^{2}\right]\\ \leq&\frac{d}{dt}\sum_{i=1}^{m}\left[C_{1}\|u_{i}\|^ {2}+\|z_{i}\|^{2}+\|\rho_{i}\|^{2}\right]+\sum_{i=1}^{m}\left(2\|u_{i}\|^{2}+ \gamma\|z_{i}\|^{2}+b\|\rho_{i}\|^{2}\right)\\ \leq&\,2C_{2}m|\Omega|,\quad\text{for}\,\,t\in I_{ max}=[0,T_{max}),\end{split} \tag{3.9}\]

where

\[\mu=\min\left\{\frac{2}{C_{1}},\;\gamma,\;b\right\}=\min\left\{\frac{8\alpha \gamma}{4\gamma+q^{2}},\;\gamma,\;b\right\}. \tag{3.10}\]
Now solve the differential inequality (3.9) to obtain the following bounding estimate of all the weak solutions on the maximal existence time interval \(I_{max}\),

\[\begin{split}&\|g(t,g^{0})\|^{2}=\sum_{i=1}^{m}\|g_{i}(t,g_{i}^{0}) \|^{2}=\sum_{i=1}^{m}\big{(}\|u_{i}(t)\|^{2}+\|w_{i}(t)\|^{2}+\|\rho_{i}(t)\|^{ 2}\big{)}\\ &\leq\frac{\max\{C_{1},1\}}{\min\{C_{1},1\}}e^{-\mu\,t}\|g^{0}\|^ {2}+\frac{2C_{2}m}{\mu\min\{C_{1},1\}}|\Omega|,\ \ \ \ t\in[0,\infty).\end{split} \tag{3.11}\]

Here it is shown that \(I_{max}=[0,\infty)\) for every weak solution \(g(t,g^{0})\) because the solution will never blow up at any finite time. The uniqueness of any weak solution to the initial value problem (2.1) is shown in Proposition 2.2. Therefore, for any initial data \(g^{0}=(g_{1}^{0},\cdots,g_{m}^{0})\in E\), there exists a unique global weak solution of the initial value problem (2.1) for this memristive reaction-diffusion neural network model (1.1)-(1.4) in the space \(E\) for time \(t\in[0,\infty)\). 

Based on the global existence of weak solutions established in Theorem 3.1, we can define the solution semiflow \(\{S(t):E\to E\}_{t\geq 0}\) of the memristive neural network system (1.1) to be

\[S(t):g^{0}\longmapsto g(t;g^{0})=\operatorname{col}\big{(}u_{i}(t,\cdot),z_{i}( t,\cdot),\rho_{i}(t,\cdot):1\leq i\leq m\big{)},\ \ \ t\geq 0.\]

We shall call this semiflow \(\{S(t)\}_{t\geq 0}\) the _memristive reaction-diffusion neural network semiflow_ generated by the model equations (1.1).

The next theorem shows that the memristive reaction-diffusion neural network semiflow \(\{S(t)\}_{t\geq 0}\) is a dissipative dynamical system in the state space \(E\).

**Theorem 3.2**.: _There exists a bounded absorbing set for the memristive reaction-diffusion neural network semiflow \(\{S(t)\}_{t\geq 0}\) in the state space \(E\), which is the bounded ball_

\[B^{*}=\{g\in E:\|g\|^{2}\leq K\} \tag{3.12}\]

_where the constant_

\[K=1+\frac{2C_{2}m}{\mu\min\{C_{1},1\}}|\Omega|, \tag{3.13}\]

_and the positive constants \(C_{1}\) and \(C_{2}\) are given in (3.4) and (3.8)._

Proof.: This is the consequence of the global uniform estimate (3.11) shown in Theorem 3.1, which implies that

\[\limsup_{t\to\infty}\|g(t,g^{0})\|^{2}=\limsup_{t\to\infty}\,\sum_{i=1}^{m}\|g _{i}(t,g_{i}^{0})\|^{2}<K \tag{3.14}\]
for all weak solutions of (2.1) with any initial data \(g^{0}\) in \(E\). Moreover, for any given bounded set \(B=\{g\in E:\|g\|^{2}\leq\mathcal{R}\}\) in \(E\), there exists a finite time

\[T_{B}=\frac{1}{\mu}\log^{+}\left(\mathcal{R}\,\frac{\max\{C_{1},1\}}{\min\{C_{1 },1\}}\right)\]

such that all the solution trajectories started at the initial time \(t=0\) from the set \(B\) will permanently enter the bounded ball \(B^{*}\) shown in (3.12) for \(t>T_{B}\). Therefore, the bounded ball \(B^{*}\) is an absorbing set in \(E\) for the semiflow \(\{S(t)\}_{t\geq 0}\) so that this memristive neural network semiflow is dissipative. 

## 4. **Higher-Order and Pointwise Ultimate Bounds**

We shall further prove an ultimate uniform bound of the membrane potential functions \(\{u_{i}(t):1\leq i\leq m\}\) for all the weak solutions in the higher-order integrable space \(L^{4}(\Omega)\). Note that Proposition 2.2 and Theorem 3.1 together show that any weak solution \(g(t)\in C((0,\infty),\Pi)\) so that the component function \(u_{i}(t)\in C((0,\infty),H^{1}(\Omega))\subset C((0,\infty),L^{6}(\Omega)) \subset C((0,\infty),L^{4}(\Omega))\).

**Theorem 4.1**.: _There exists a constant \(Q>0\) such that for any initial data \(g^{0}\in E\), the membrane potential components \(u_{i},1\leq i\leq m\), of the weak solution \(g(t,g^{0})=(g_{1}(t),\cdots,g_{m}(t))\) of the initial value problem (2.1) for the memristive reaction-diffusion neural network \(\mathcal{NW}\) are ultimately uniform bounded in the space \(L^{4}(\Omega)\) and_

\[\limsup_{t\to\infty}\,\sum_{i=1}^{m}\|u_{i}(t)\|_{L^{4}}^{4}<Q. \tag{4.1}\]

Proof.: Take the \(L^{2}\) inner-product of the \(u_{i}\)-equation in (1.1) with \(u_{i}^{3}(t,\cdot),1\leq i\leq m\), and sum them up. By the boundary condition (1.3) and Assumption (1.5), we have

\[\begin{split}&\frac{1}{4}\,\frac{d}{dt}\sum_{i=1}^{m}\|u_{i}(t) \|_{L^{4}}^{4}+3\eta\sum_{i=1}^{m}\|u_{i}\nabla u_{i}\|_{L^{2}}^{2}\\ &=\sum_{i=1}^{m}\int_{\Omega}(f(u_{i},x)u_{i}^{3}-k\tanh(\rho_{i })u_{i}^{4})\,dx-\sum_{i=1}^{m}\sum_{j=1}^{m}\int_{\Omega}\frac{Pu_{i}^{4}}{1+ \exp[-r(u_{j}-V)]}\,dx\\ &\leq\sum_{i=1}^{m}\int_{\Omega}\left(-\alpha u_{i}^{6}+\lambda u _{i}^{3}|z_{i}|+Ju_{i}^{3}+ku_{i}^{4}\right)dx,\quad t>0.\end{split} \tag{4.2}\]

By Young's inequality (2.6), it is seen that

\[\lambda u_{i}^{3}|z_{i}|+Ju_{i}^{3}+ku_{i}^{4}\leq\left(\frac{\alpha}{4}u_{i} ^{6}+\frac{\lambda^{2}}{\alpha}z_{i}^{2}\right)+\left(\frac{\alpha}{4}u_{i}^{ 6}+\frac{J^{2}}{\alpha}\right)+\left(\frac{\alpha}{4}u_{i}^{6}+\frac{64}{27 \alpha^{2}}k^{3}\right), \tag{4.3}\]
where the last two terms in a bracket come from

\[ku_{i}^{4}=\left[\frac{3\alpha}{8}u_{i}^{6}\right]^{2/3}\left[\frac{64}{9\alpha^{ 2}}k^{3}\right]^{1/3}\leq\frac{2}{3}\left(\frac{3\alpha}{8}u_{i}^{6}\right)+ \frac{1}{3}\left(\frac{64}{9\alpha^{2}}k^{3}\right)=\frac{\alpha}{4}u_{i}^{6}+ \frac{64}{27\alpha^{2}}k^{3}.\]

Note that

\[u_{i}^{4}\leq\frac{1}{3}+\frac{2}{3}u_{i}^{6}\leq 1+u_{i}^{6}\quad\text{so that} \,-u_{i}^{6}\leq-u_{i}^{4}+1. \tag{4.4}\]

From (4.2) wherein we can use the inequalities (4.3) and (3.14), it follows that

\[\frac{1}{4}\,\frac{d}{dt}\sum_{i=1}^{m}\|u_{i}(t)\|_{L^{4}}^{4}+3 \eta\sum_{i=1}^{m}\|u_{i}\nabla u_{i}\|^{2} \tag{4.5}\] \[\leq \sum_{i=1}^{m}\int_{\Omega}\left[-\alpha u_{i}^{6}+\left(\frac{ \alpha}{4}u_{i}^{6}+\frac{\lambda^{2}}{\alpha}z_{i}^{2}\right)+\left(\frac{ \alpha}{4}u_{i}^{6}+\frac{J^{2}}{\alpha}\right)+\left(\frac{\alpha}{4}u_{i}^{ 6}+\frac{64}{27\alpha^{2}}k^{3}\right)\right]dx\] \[= \sum_{i=1}^{m}\left(-\frac{\alpha}{4}\int_{\Omega}u_{i}^{6}\,dx+ \frac{\lambda^{2}}{\alpha}\|z_{i}(t)\|^{2}\right)+m\left(\frac{J^{2}}{\alpha }+\frac{64k^{3}}{27\alpha^{2}}\right)|\Omega|\] \[\leq \sum_{i=1}^{m}\left(-\frac{\alpha}{4}\int_{\Omega}u_{i}^{4}\,dx+ \frac{\lambda^{2}}{\alpha}\|z_{i}(t)\|^{2}\right)+m\left(\frac{\alpha}{4}+ \frac{J^{2}}{\alpha}+\frac{64k^{3}}{27\alpha^{2}}\right)|\Omega|\] \[< -\frac{\alpha}{4}\sum_{i=1}^{m}\|u_{i}(t)\|_{L^{4}}^{4}+\frac{ \lambda^{2}}{\alpha}K+m\left(\frac{\alpha}{4}+\frac{J^{2}}{\alpha}+\frac{64k^ {3}}{27\alpha^{2}}\right)|\Omega|,\]

for time \(t\) sufficiently large, where the constant \(K\) is given in (3.13), which is valid for all the weak solutions. Consequently, with the nonnegative gradient terms removed, the differential inequality (4.5) shows that

\[\frac{d}{dt}\sum_{i=1}^{m}\|u_{i}(t)\|_{L^{4}}^{4}+\alpha\sum_{i= 1}^{m}\|u_{i}(t)\|_{L^{4}}^{4} \leq\frac{4\lambda^{2}}{\alpha}K+m\left(\alpha+\frac{4J^{2}}{ \alpha}+\frac{256k^{3}}{27\alpha^{2}}\right)|\Omega|. \tag{4.6}\] \[<\frac{4\lambda^{2}}{\alpha}K+m\left(\alpha+\frac{4J^{2}}{\alpha }+\frac{10k^{3}}{\alpha^{2}}\right)|\Omega|,\]

for \(t>T(g^{0})\), where \(T(g^{0})>0\) is a finite time when the solution trajectory started from the initial state \(g^{0}\) be absorbed into the absorbing set \(B^{*}\) in the state space \(E\), as shown in Theorem 3.2.

By the parabolic regularity stated in Proposition 2.2, for any weak solution \(g(t,g^{0})\) one has \(u_{i}(T(g^{0}))\in H^{1}(\Omega)\subset L^{4}(\Omega)\) for \(1\leq i\leq m\). Then the second statement in Proposition 2.2 shows that any weak solution has the regularity

\[\sum_{i=1}^{m}u_{i}(t)\in C([T(g^{0}),\infty),H^{1}(\Omega))\subset C([T(g^{0} ),\infty),L^{4}(\Omega)).\]
Apply the Gronwall inequality to (4.6). It results in the bounding estimate of all the \(u_{i}\) components, \(1\leq i\leq m\), in the space \(L^{4}(\Omega)\) as follows:

\[\begin{split}\sum_{i=1}^{m}\|u_{i}(t)\|_{L^{4}}^{4}&<e^ {-\alpha(t-T(g^{0}))}\sum_{i=1}^{m}\|u_{i}(T(g^{0}))\|_{L^{4}}^{4}\\ &+\frac{4\lambda^{2}}{\alpha^{2}}K+m\left(1+\frac{4J^{2}}{\alpha ^{2}}+\frac{10k^{3}}{\alpha^{3}}\right)|\Omega|,\quad\text{for}\,\,\,t\geq 0, \end{split} \tag{4.7}\]

Therefore (4.1) is proved with

\[Q=1+\frac{4\lambda^{2}}{\alpha^{2}}K+m\left(1+\frac{4J^{2}}{\alpha^{2}}+\frac{ 10k^{3}}{\alpha^{3}}\right)|\Omega|. \tag{4.8}\]

which is a constant independent of any initial data. 

The pointwise estimation in the following theorem will be used to deal with the nonlinear weak coupling toward exponential synchronization featured in this work.

**Theorem 4.2**.: _There exists a constant \(G>0\) such that for any initial data \(g^{0}\in E\), the membrane potential component \(u_{i}(t,x),1\leq i\leq m\), of the weak solution \(g(t,g^{0})=(g_{1}(t),\cdots,g_{m}(t))\) of the initial value problem (2.1) for the memristive reaction-diffusion neural network \(\mathcal{NW}\) is ultimately uniform bounded in the space \(\mathbb{R}\) and_

\[\limsup_{t\to\infty}\,\sum_{i=1}^{m}|u_{i}(t,x)|_{\mathbb{R}}<G,\quad\text{for} \,\,\,x\in\Omega. \tag{4.9}\]

Proof.: Similar to (3.1) in the proof of Theorem 3.1, we can multiply the \(u_{i}\)-equation in (1.1) by \(C_{1}u_{i}(t,x)\) for \(1\leq i\leq m\), where \(C_{1}\) is the same constant given in (3.4), and sum them up to get

\[\begin{split}&\frac{C_{1}}{2}\frac{d}{dt}\sum_{i=1}^{m}|u_{i}(t,x) |^{2}+C_{1}\sum_{i=1}^{m}\eta|\nabla u_{i}(t,x)|^{2}\\ \leq&\sum_{i=1}^{m}\left[-C_{1}\alpha\,u_{i}^{4}(t,x )+\left(C_{1}k+\frac{C_{1}^{2}\lambda^{2}}{\gamma}\right)u_{i}^{2}(t,x)+\frac {\gamma}{4}z_{i}^{2}(t,x)\right]+C_{1}mJ.\end{split} \tag{4.10}\]

Similar to (3.2), by multiplication and summation but without integration, we have

\[\begin{split}&\frac{1}{2}\frac{d}{dt}\sum_{i=1}^{m}\left(|z_{i}(t,x )|^{2}+|\rho_{i}(t,x)|^{2}\right)\\ \leq&\sum_{i=1}^{m}\left[\frac{q^{2}}{4\gamma}u_{i}^ {4}(t,x)-\frac{3\gamma}{4}z_{i}^{2}(t,x)+\frac{a^{2}}{2b}u_{i}^{2}(t,x)-\frac {b}{2}\rho_{i}^{2}(t,x)\right]+\frac{mL^{2}}{\gamma}.\end{split} \tag{4.11}\]
Then parallel to the steps from (3.3) through (3.9) in the proof of Theorem 3.1, one can reach the pointwise differential inequality

\[\begin{split}&\frac{d}{dt}\sum_{i=1}^{m}\big{[}C_{1}u_{i}^{2}(t,x)+z _{i}^{2}(t,x)+\rho_{i}^{2}(t,x)\big{]}\\ &+\mu\sum_{i=1}^{m}\big{[}C_{1}u_{i}^{2}(t,x)+z_{i}^{2}(t,x)+\rho _{i}^{2}(t,x)\big{]}\leq 2C_{2}m,\quad t>0,\ x\in\Omega,\end{split} \tag{4.12}\]

where \(C_{2}\) and \(\mu\) are two universal positive constants given in (3.8) and (3.10) respectively. It follows that

\[\sum_{i=1}^{m}\big{(}u_{i}^{2}(t,x)+z_{i}^{2}(t,x)+\rho_{i}^{2}(t,x)\big{)} \leq\frac{\max\{C_{1},1\}}{\min\{C_{1},1\}}e^{-\mu\,t}|g^{0}(t,x)|^{2}+\frac{2 C_{2}m}{\mu\min\{C_{1},1\}}, \tag{4.13}\]

for \(t>0,x\in\Omega\), which implies that (4.9) is valid with a uniform constant

\[G=\left[1+\frac{2C_{2}m}{\mu\min\{C_{1},1\}}\right]^{1/2}, \tag{4.14}\]

which is independent of any initial data. 

## 5. **Synchronization of Memristive Reaction-Diffusion Neural Networks**

In this section, we shall prove the main result on the synchronization of the memristive reaction-diffusion neural networks described by (1.1) in the state space \(E\). This result provides a quantitative threshold condition for the interneuron coupling strength to reach the neural network synchronization.

**Definition 5.1**.: For a model evolutionary equation of a general neural network called _GNW_, such as (2.1) formulated from the memristive reaction-diffusion equations (1.1), we define the asynchronous degree of this neural network in a state space (as a Banach space) \(W\) to be

\[deg_{s}\left(\text{\emph{GNW}}\right)=\sum_{1\,\leq i\,<j\,\leq\,m}\left\{ \sup_{g_{i}^{0},\,g_{j}^{0}\,\in\,W}\ \left\{\limsup_{t\to\infty}\,\|g_{i}(t;g_{i}^{0})-g_{j}(t;g_{j}^{0})\|_{W} \right\}\right\}\]

where \(g_{i}(t)\) and \(g_{j}(t)\) are any two solutions of this model evolutionary equation with the initial states \(g_{i}^{0}\) and \(g_{j}^{0}\) respectively for two neurons \(\mathcal{N}_{i}\) and \(\mathcal{N}_{j}\), \(1\leq i,j\leq m\), in the network. The neural network is said to be asymptotically synchronized if

\[deg_{s}\left(\text{\emph{GNW}}\right)=0.\]
If the asymptotic convergence to zero of the difference norm for any two neurons in the network admits a uniform exponential rate, then the neural network is called exponentially synchronized.

Introduce the neuron difference functions: For \(i,j=1,\cdots,m\), define

\[U_{ij}(t,x) =u_{i}(t,x)-u_{j}(t,x),\] \[Z_{ij}(t,x) =z_{i}(t,x)-z_{j}(t,x),\] \[R_{ij}(t,x) =\rho_{i}(t,x)-\rho_{j}(t,x).\]

Given any initial state \(g^{0}=\operatorname{col}\left(g^{0}_{1},\cdots,g^{0}_{m}\right)\) in the space \(E\), the difference between any two solutions of (2.1) associated with two neurons \(\mathcal{N}_{i}\) and \(\mathcal{N}_{j}\) in the network is what we consider:

\[g_{i}(t,g^{0}_{i})-g_{j}(t,g^{0}_{j})=\operatorname{col}\left(U_{ij}(t,\cdot),Z _{ij}(t,\cdot),R_{ij}(t,\cdot)\right),\quad t\geq 0.\]

By subtraction of the three governing equations for the \(j\)-th neuron from the corresponding governing equations for the \(i\)-th neuron in (1.1), we obtain the following differencing reaction-diffusion equations. For \(i,j=1,\cdots,m\),

\[\begin{split}\frac{\partial U}{\partial t}=&\, \eta\Delta U+f(u_{i},z_{i})-f(u_{j},z_{j})-k(\tanh(\rho_{i})u_{i}-\tanh(\rho_{ j})u_{j})\\ &-P\left[u_{i}\sum_{\nu=1}^{m}\Gamma(u_{\nu})-u_{j}\sum_{\nu=1}^{ m}\Gamma(u_{\nu})\right].\\ \frac{\partial Z}{\partial t}=&\,\Lambda Z+h(u_{i}, z_{i})-h(u_{j},z_{j}),\\ \frac{\partial R}{\partial t}=&\,aU-bR.\end{split} \tag{5.1}\]

Here and after, for any given \(i\) and \(j\), we shall simply write \(U(t,x)=U_{ij}(t,x),Z(t,x)=Z_{ij}(t,x),R(t,x)=R_{ij}(t,x)\) as a notational convenience.

The following exponential synchronization theorem is the main result of this paper.

**Theorem 5.2**.: _For memristive reaction-diffusion neural networks \(\mathcal{NW}\) with the model (1.1)-(1.3) and the Assumptions (1.5)-(1.6), if the following threshold condition is satisfied by the coupling strength coefficient \(P\),_

\[\begin{split} P>&\frac{1+\exp[r(G+|V|)]}{m}\times\\ &\times\left[\beta+\frac{\beta^{2}+\xi^{2}}{\gamma}+k+\frac{a^{2} }{b}+2\sqrt{2Q}\,C^{*}\left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right)+ \frac{64Q^{2}C^{*4}}{\eta^{3}}\left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma} \right]^{4}\right]\end{split} \tag{5.2}\]
_where the positive constant \(Q\) and \(G\) are given in (4.8) and (4.14) respectively and \(C^{*}\) is a coefficient in the Gagliardo-Nirenberg interpolation inequality (5.9), then the memristive neural network \(\mathcal{NW}\) is exponentially synchronized in the state space \(E\) at a uniform exponential rate_

\[\delta(P)=\min\left\{b,\,\gamma,\,2\left(\frac{mP}{1+\exp[r(G+|V|)]}-\kappa \right)\right\}, \tag{5.3}\]

_where the positive constant \(\kappa\) is_

\[\kappa=\beta+\frac{\beta^{2}+\xi^{2}}{\gamma}+k+\frac{a^{2}}{b}+2\sqrt{2Q}\,C ^{*}\left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right)+\frac{64Q^{2}C^{*4}}{ \eta^{3}}\left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right]^{4}. \tag{5.4}\]

Proof.: The proof will go through three steps.

Step 1. Take the \(L^{2}\) inner-products of the first equation in (5.1) with \(U(t)\), the second equation in (5.1) with \(Z(t)\), and the third equation in (5.1) with \(R(t)\). Then sum them up and use the Assumptions (1.5)-(1.6) to get

\[\frac{1}{2}\frac{d}{dt}(\|U(t)\|^{2}+\|Z(t)\|^{2}+\|R(t)\|^{2})+ \eta\|\nabla U(t)\|^{2}+\gamma\,\|Z(t)\|^{2}+b\|R(t)\|^{2} \tag{5.5}\] \[+P\int_{\Omega}\left[u_{i}\sum_{\nu=1}^{m}\Gamma(u_{\nu})-u_{j} \sum_{\nu=1}^{m}\Gamma(u_{\nu})\right]U(t,x)\,dx\] \[= \int_{\Omega}(f(u_{i},z_{i})-f(u_{j},z_{j}))U\,dx-\int_{\Omega}k (\tanh(\rho_{i})u_{i}-\tanh(\rho_{j})u_{j})U]\,dx\] \[+\int_{\Omega}(h(u_{i},z_{i})-h(u_{j},z_{j}))Z\,dx+\int_{\Omega }aUR\,dx\] \[\leq \int_{\Omega}\frac{\partial f}{\partial s}\left(\zeta u_{i}+(1- \zeta)u_{j}\right)U^{2}\,dx+\int_{\Omega}\frac{\partial f}{\partial\sigma} \left(\varsigma z_{i}+(1-\varsigma)z_{j}\right)UZ\,dx\] \[+\int_{\Omega}\frac{\partial h}{\partial s}\left(\epsilon u_{i} +(1-\epsilon)u_{j}\right)UZ\,dx\] \[-k\int_{\Omega}\left[\operatorname{sech}^{2}(\varepsilon\rho_{i} +(1-\varepsilon)\rho_{j})R\,u_{i}U+\tanh(\rho_{j})U^{2}\right]dx+\int_{ \Omega}aUR\,dx\] \[\leq \int_{\Omega}\left(\beta(U^{2}+|UZ|)+\xi(|u_{i}|+|u_{j}|+1)|UZ|+k (|u_{i}|RU+U^{2})+aUR\right)dx.\] \[\leq \int_{\Omega}\left(\beta+\frac{\beta^{2}+\xi^{2}}{\gamma}+k+ \frac{a^{2}}{b}\right)U^{2}(t,x)dx+\frac{\gamma}{4}\|Z(t)\|^{2}+\frac{b}{4}\| R(t)\|^{2}\] \[+\int_{\Omega}\xi(|u_{i}|+|u_{j}|)|UZ|\,dx+\int_{\Omega}k|u_{i}| RU\,dx,\quad t>0.\]
where the mean value theorem in differentiation and the hyperbolic function properties \(|\tanh(\rho_{j})|\leq 1,\operatorname{sech}^{2}(\alpha\rho_{i}+(1-\alpha)\rho_{j})\leq 1\) are used and the numbers \(\zeta,\varsigma,\epsilon,\varepsilon\in[0,1]\).

Step 2. We treat the last two integral terms on the right-hand side of the inequality (5.5). By the Holder inequality,

\[\begin{split}&\int_{\Omega}k|u_{i}|RU\,dx\leq k\int_{\Omega} \left(\frac{b}{4k}R^{2}(t,x)+\frac{k}{b}u_{i}^{2}(t,x)U^{2}(t,x)\right)dx\\ \leq&\,\frac{b}{4}\|R(t)\|^{2}+\frac{k^{2}}{b}\left[ \int_{\Omega}u_{i}^{4}(t,x)\,dx\right]^{1/2}\left[\int_{\Omega}U^{4}(t,x)\,dx \right]^{1/2}\\ =&\,\frac{b}{4}\|R(t)\|^{2}+\frac{k^{2}}{b}\|u_{i}( t)\|_{L^{4}}^{2}\|U(t)\|_{L^{4}}^{2},\quad t>0.\end{split} \tag{5.6}\]

Similarly we have

\[\begin{split}&\int_{\Omega}\xi(|u_{i}|+|u_{j}|)|UZ|\,dx\leq\xi \int_{\Omega}\left(\frac{\gamma}{4\xi}Z^{2}(t,x)+\frac{2\xi}{\gamma}(u_{i}^{2} +u_{j}^{2})U^{2}(t,x)\right)dx\\ \leq&\,\frac{\gamma}{4}\|Z(t)\|^{2}+\frac{2\xi^{2}} {\gamma}\left(\left[\int_{\Omega}u_{i}^{4}\,dx\right]^{1/2}+\left[\int_{ \Omega}u_{j}^{4}\,dx\right]^{1/2}\right)\left[\int_{\Omega}U^{4}(t,x)\,dx \right]^{1/2}\\ =&\,\frac{\gamma}{4}\|Z(t)\|^{2}+\frac{2\xi^{2}}{ \gamma}\left(\|u_{i}(t)\|_{L^{4}}^{2}+\|u_{j}(t)\|_{L^{4}}^{2}\right)\|U(t)\|_ {L^{4}}^{2},\quad t>0.\end{split} \tag{5.7}\]

Substitute the term estimates (5.6) and (5.7) into the differential inequality (5.5). We obtain

\[\begin{split}&\frac{1}{2}\frac{d}{dt}(\|U(t)\|^{2}+\|W(t)\|^{2}+\|R(t)\|^{2})+ \eta\|\nabla U(t)\|^{2}+\frac{\gamma}{2}\|Z(t)\|^{2}+\frac{b}{2}\|R(t)\|^{2} \\ &+P\int_{\Omega}\left[u_{i}\sum_{\nu=1}^{m}\Gamma(u_{\nu})-u_{j} \sum_{\nu=1}^{m}\Gamma(u_{\nu})\right]U(t,x)\,dx\\ \leq&\,\left[\beta+\frac{\beta^{2}+\xi^{2}}{\gamma} +k+\frac{a^{2}}{b}\right]\|U(t)\|^{2}\\ &+\left[\left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right)\|u_ {i}(t)\|_{L^{4}}^{2}+\frac{2\xi^{2}}{\gamma}\|u_{j}(t)\|_{L^{4}}^{2}\right]\| U(t)\|_{L^{4}}^{2},\quad t>0.\end{split} \tag{5.8}\]

The challenge is to handle the last two terms of \(L^{4}\)-norm products on the right-hand side of the above inequality (5.8). We exploit the Gagliardo-Nirenberg interpolation inequalities [32, Theorem B.3] and [4]. It states that Sobolev embedding

\[H^{1}(\Omega)\subset L^{4}(\Omega)\subset L^{2}(\Omega)\]
implies

\[\begin{split}&\|U(t)\|_{L^{4}}^{2}\leq C^{*}\|U(t)\|_{H^{1}}^{2 \theta}\|U(t)\|^{2(1-\theta)}\\ \leq&\,C^{*}(\|U(t)\|+\|\nabla U(t)\|)^{2\theta}\|U(t )\|^{2(1-\theta)}\\ \leq&\,C^{*}2^{2\theta}(\|U(t)\|^{2\theta}+\|\nabla U (t)\|^{2\theta})\|U(t)\|^{2(1-\theta)}\\ =&\,2\sqrt{2}C^{*}\|U(t)\|^{2}+2\sqrt{2}C^{*}\|\nabla U (t)\|^{3/2}\|U(t)\|^{2(1-3/4)}\end{split} \tag{5.9}\]

where an inequality in [5, Theorem 4.7] is used and the coefficient \(C^{*}(\Omega)>0\) only depends on the spatial domain \(\Omega\). Here the interpolation index \(\theta=3/4\) is determined by

\[-\frac{n}{4}\leq\theta\left(1-\frac{n}{2}\right)-(1-\theta)\,\frac{n}{2}, \quad\text{for }1\leq n=\dim\Omega\leq 3,\]

and the equality holds for \(n=3\). The interpolation inequality (5.9) shows that

\[\|U(t)\|_{L^{4}}^{2}\leq 2\sqrt{2}C^{*}\|U(t)\|^{2}+2\sqrt{2}C^{*}\|\nabla U(t )\|^{3/2}\|U(t)\|^{1/2}. \tag{5.10}\]

According to Theorem 3.2 and Theorem 4.1, we know that \(\limsup_{t\to\infty}\|U(t)\|^{2}<K\) and \(\limsup_{t\to\infty}\sum_{i=1}^{m}\|u_{i}(t)\|_{L^{4}}^{4}<Q\). Thus for any given initial state \(g^{0}\in E\) there exists a finite time \(T(g^{0})\geq 0\) such that

\[\sum_{i=1}^{m}\|u_{i}(t)\|_{L^{4}}^{2}<Q^{1/2},\quad\text{for all}\;\;t>T(g^{0}).\]

Therefore, by (5.10) and Young's inequality (2.6), we achieve the estimate

\[\begin{split}&\left[\left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma} \right)\|u_{i}(t)\|_{L^{4}}^{2}+\frac{2\xi^{2}}{\gamma}\|u_{j}(t)\|_{L^{4}}^{2 }\right]\|U(t)\|_{L^{4}}^{2}\\ =&\left[\frac{k^{2}}{b}\|u_{i}(t)\|_{L^{4}}^{2}+\frac {2\xi^{2}}{\gamma}\left(\|u_{i}(t)\|_{L^{4}}^{2}+\|u_{j}(t)\|_{L^{4}}^{2} \right)\right]\|U(t)\|_{L^{4}}^{2}\\ \leq&\,\left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma} \right)Q^{1/2}2\sqrt{2}C^{*}(\|U(t)\|^{2}+\|\nabla U(t)\|^{3/2}\|U(t)\|^{1/2} )\\ \leq&\,\left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma} \right)Q^{1/2}2\sqrt{2}C^{*}\|U(t)\|^{2}+\eta\|\nabla U(t)\|^{(3/2)\times(4/3) }\\ &\,+\frac{1}{\eta^{3}}\left[\left(\frac{k^{2}}{b}+\frac{2\xi^{2} }{\gamma}\right)Q^{1/2}2\sqrt{2}C^{*}\|U(t)\|^{1/2}\right]^{4}\\ =&\,\eta\|\nabla U(t)\|^{2}+\left[2\sqrt{2Q}\,C^{*} \left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right)+\frac{64Q^{2}C^{*4}}{ \eta^{3}}\left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right)^{4}\right]\| U(t)\|^{2},\end{split} \tag{5.11}\]

for \(t>T(g^{0})\).
Substitute (5.11) in (5.8) and then cancel the gradient terms \(\eta\|\nabla U(t)\|^{2}\) on two sides of that inequality. It follows that

\[\begin{split}&\frac{1}{2}\frac{d}{dt}(\|U(t)\|^{2}+\|Z(t)\|^{2}+\|R(t )\|^{2})+\frac{\gamma}{2}\|Z(t)\|^{2}+\frac{b}{2}\|R(t)\|^{2}\\ &+P\int_{\Omega}\left[u_{i}\sum_{\nu=1}^{m}\Gamma(u_{\nu})-u_{j} \sum_{\nu=1}^{m}\Gamma(u_{\nu})\right]U(t,x)\,dx\\ &\leq\left[\beta+\frac{\beta^{2}+\xi^{2}}{\gamma}+k+\frac{a^{2}}{ b}+2\sqrt{2Q}\,C^{*}\left(\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right)+\frac{64Q^{2} C^{*4}}{\eta^{3}}\left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right]^{4} \right]\|U(t)\|^{2}.\end{split} \tag{5.12}\]

Step 3. Another challenge is to handle the nonlinear difference term of the weak coupling on the left-hand side of the inequality (5.12). For any given \(1\leq i\neq j\leq m\), we have

\[\begin{split}& P\int_{\Omega}\left[u_{i}\sum_{\nu=1}^{m}\Gamma(u_{ \nu})-u_{j}\sum_{\nu=1}^{m}\Gamma(u_{\nu})\right]U(t,x)\,dx\\ =&\,P\int_{\Omega}\sum_{\nu=1}^{m}\frac{u_{i}-u_{j}} {1+\exp[-r(u_{\nu}-V)]}U(t,x)\,dx\\ =&\,P\int_{\Omega}\sum_{\nu=1}^{m}\frac{U^{2}(t,x)} {1+\exp[-r(u_{\nu}-V)]}\,dx.\end{split} \tag{5.13}\]

By Theorem 4.2 and (4.9), for each solution trajectory \(g(t,g^{0})\) there exists a finite time \(\tau(g^{0})>0\) such that \(\sum_{i=1}^{m}|u_{i}(t,x)|_{\mathbb{R}}<G\) for \(t>\tau(g^{0})\). Hence it holds that

\[\frac{1}{1+\exp\left[-r(u_{\nu}(t,x)-V)\right]}\geq\frac{1}{1+\exp\left[r(G+| V|)\right]},\quad t>\tau(g^{0}), \tag{5.14}\]

for all \(1\leq\nu\leq m\). Now substitute (5.13) and (5.14) in the differential inequality (5.12) on the left-hand side. We obtain
\[\frac{1}{2}\frac{d}{dt}(\|U(t)\|^{2}+\|Z(t)\|^{2}+\|R(t)\|^{2})+ \frac{\gamma}{2}\|Z(t)\|^{2}+\frac{b}{2}\|R(t)\|^{2}\] \[+\frac{mP}{1+\exp\left[r(G+|V|)\right]}\|U(t)\|^{2}\] \[= \frac{1}{2}\frac{d}{dt}(\|U(t)\|^{2}+\|Z(t)\|^{2}+\|R(t)\|^{2})+ \frac{\gamma}{2}\|Z(t)\|^{2}+\frac{b}{2}\|R(t)\|^{2}\] \[+\frac{P}{1+\exp\left[r(G+|V|)\right]}\int_{\Omega}\sum_{\nu=1}^ {m}U^{2}(t,x)\,dx\] \[\leq \frac{1}{2}\frac{d}{dt}(\|U(t)\|^{2}+\|Z(t)\|^{2}+\|R(t)\|^{2})+ \frac{\gamma}{2}\|Z(t)\|^{2}+\frac{b}{2}\|R(t)\|^{2}\] \[+P\,\sum_{\nu=1}^{m}\int_{\Omega}\frac{1}{1+\exp\left[-r(u_{\nu} (t,x)-V)\right]}U^{2}(t,x)\,dx\] \[= \frac{1}{2}\frac{d}{dt}(\|U(t)\|^{2}+\|Z(t)\|^{2}+\|R(t)\|^{2})+ \frac{\gamma}{2}\|Z(t)\|^{2}+\frac{b}{2}\|R(t)\|^{2}\] \[+P\int_{\Omega}\left[u_{i}\sum_{\nu=1}^{m}\Gamma(u_{\nu})-u_{j} \sum_{\nu=1}^{m}\Gamma(u_{\nu})\right]U(t,x)\,dx\] \[\leq \left[\beta+\frac{\beta^{2}+\xi^{2}}{\gamma}+k+\frac{a^{2}}{b}+2 \sqrt{2Q}\,C^{*}\left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right]+\frac{64 Q^{2}C^{*4}}{\eta^{3}}\left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right]^{4} \right]\|U(t)\|^{2} \tag{5.15}\]

for \(t>\tau(g^{0})\). From (5) and by the threshold condition (5.2) stated in this theorem, it results in the following Gronwall-type inequality:

\[\frac{d}{dt}(\|U(t)\|^{2}+\|Z(t)\|^{2}+\|R(t)\|^{2})+\delta(P)(\| U(t)\|^{2}+\|Z(t)\|^{2}+\|R(t)\|^{2})\] \[\leq\frac{d}{dt}(\|U(t)\|^{2}+\|W(t)\|^{2}+\|R(t)\|^{2})+\gamma\| Z(t)\|^{2}+b\|R(t)\|^{2}+2\left[\frac{mP}{1+\exp\left[r(G+|V|)\right]}\right.\] \[-\left[\beta+\frac{\beta^{2}+\xi^{2}}{\gamma}+k+\frac{a^{2}}{b}+ 2\sqrt{2Q}\,C^{*}\left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right]+\frac{64 Q^{2}C^{*4}}{\eta^{3}}\left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right]^{4} \right]\right]\|U(t)\|^{2}\] \[\leq 0,\quad\text{for}\;\;t>\tau(g^{0}). \tag{5.16}\]
Denote by

\[\kappa=\beta+\frac{\beta^{2}+\xi^{2}}{\gamma}+k+\frac{a^{2}}{b}+2\sqrt{2Q}\,C^{*} \left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right]+\frac{64Q^{2}C^{*4}}{\eta^ {3}}\left[\frac{k^{2}}{b}+\frac{2\xi^{2}}{\gamma}\right]^{4}.\]

Finally we can solve this linear Gronwall inequality (5.16) to reach the exponential synchronization result: For any initial state \(g^{0}\in E\) and any two neurons \(\mathcal{N}_{i}\) and \(\mathcal{N}_{j}\) in this memristive reaction-diffusion neural network model (1.1), their difference function \(g_{i}(t;g_{i}^{0})-g_{j}(t;g_{j}^{0})\) converges to zero in the state space \(E\) exponentially at a uniform convergence rate \(\delta(P)\) shown below. Namely, for any \(1\leq i\neq j\leq m\),

\[\begin{split}\|g_{i}(t)-g_{j}(t)\|_{E}^{2}&=\|u_{i} (t)-u_{j}(t)\|^{2}+\|z_{i}(t)-z_{j}(t)\|^{2}+\|\rho_{i}(t)-\rho_{j}(t)\|^{2}\\ &=\|U_{ij}(t)\|^{2}+\|Z_{ij}(t)\|^{2}+\|R_{ij}(t)\|^{2}\\ &\leq e^{-\delta(P)\,t}\left\|g_{i}^{0}-g_{j}^{0}\right\|^{2} \to 0,\ \ \text{as}\,\,t\to\infty.\end{split} \tag{5.17}\]

Here the constant convergence rate in (5.17) is

\[\delta(P)=\min\left\{b,\,\gamma,\,2\left(\frac{mP}{1+\exp\left[r(G+|V|)\right] }-\kappa\right)\right\},\]

which is exactly (5.3)-(5.4) stated in the threshold condition (5.2) of this theorem. Hence it is proved that

\[deg_{s}(\mathcal{N}\mathcal{W})=\sum_{1\,\leq\,i\,\neq\,j\,\leq\,m}\left\{\sup _{g^{0}\,\in\,E}\,\left\{\limsup_{t\to\infty}\|g_{i}(t)-g_{j}(t)\|_{E}^{2} \right\}\right\}=0. \tag{5.18}\]

The proof of this theorem is completed. 

As a meaningful extension of Theorem 5.2, we can also prove the exponential synchronization of memristive reaction-diffusion neural networks denoted by \(\mathbb{N}\mathbb{W}=\{N_{i}:i=1,2,\cdots,m\}\) with the following model equations, cf. [3, 9, 24, 28],

\[\begin{split}&\frac{\partial u_{i}}{\partial t}=\eta\Delta u_{i}+f(u_{ i},z_{i})-k\tanh(\rho_{i})u_{i}-\sum_{j=1}^{m}\frac{P(u_{i}-u_{e})}{1+\exp[-r(u_{ j}-V)]},\\ &\frac{\partial z_{i}}{\partial t}=\Lambda z_{i}+h(u_{i},z_{i}), \\ &\frac{\partial\rho_{i}}{\partial t}=au_{i}-b\rho_{i},\end{split} \tag{5.19}\]

where the weak coupling terms involve a constant \(u_{e}\in\mathbb{R}\) called the reversal potential, on a bounded spacial domain \(\Omega\) and satisfy the same boundary conditions as specified in Section 1.
**Theorem 5.3**.: _Assume that the nonlinear terms \(f(s,\sigma)\) and \(h(s,\sigma)\) in the memristive neural network model (5.19) are respectively scalar and vector polynomials and satisfy the same Assumptions (1.5) and (1.6), then there exists a positive constant \(\Psi>0\) which depends only on the parameters including \(u_{e}\) but independent of any initial data, such that if the threshold condition_

\[P>\Psi \tag{5.20}\]

_is satisfied, then the solution semiflow of the memristive reaction-diffusion neural network \(\mathbb{NW}\) will be exponentially synchronized in the same state space \(E\) at a uniform convergence rate._

Proof.: We just briefly sketch the proof. Make the variable changes to denote \(\tilde{u}_{i}=u_{i}-u_{e}\), \(1\leq i\leq m\). Then the system (5.19) becomes

\[\frac{\partial\tilde{u}_{i}}{\partial t}=\eta\Delta\tilde{u}_{i}+f(\tilde{u}_ {i}+u_{e},z_{i})\,-k\tanh(\rho_{i})(\tilde{u}_{i}+u_{e})-\sum_{j=1}^{m}\frac{P \,\tilde{u}_{i}}{1+\exp\left[-r(\tilde{u}_{j}+u_{e}-V)\right]}\,,\]

\[\frac{\partial z_{i}}{\partial t}=\Lambda z_{i}+h(\tilde{u}_{i}+u_{e},z_{i}),\]

\[\frac{\partial\rho_{i}}{\partial t}=a(\tilde{u}_{i}+u_{e})-b\rho_{i}.\]

This system of equations can be written as

\[\frac{\partial\tilde{u}_{i}}{\partial t}=\eta\Delta\tilde{u}_{i} +\tilde{f}(\tilde{u}_{i},z_{i},\rho_{i})-\sum_{j=1}^{m}\frac{P\,\tilde{u}_{i}}{ 1+\exp[-r(\tilde{u}_{j}-\tilde{V})]}, \tag{5.21}\] \[\frac{\partial z_{i}}{\partial t}=\Lambda z_{i}+\tilde{h}(\tilde {u}_{i},z_{i}),\] \[\frac{\partial\rho_{i}}{\partial t}=a\tilde{u}_{i}-b\rho_{i}+au_{ e}.\]

where \(\tilde{V}=V-u_{e}\) and the two new functions \(\tilde{f}\) and \(\tilde{h}\) are

\[\tilde{f}(\tilde{u}_{i},z_{i},\rho_{i}) =f(\tilde{u}_{i}+u_{e},z_{i})-k\tanh(\rho_{i})(\tilde{u}_{i}+u_{e}), \tag{5.22}\] \[\tilde{h}(\tilde{u}_{i},z_{i}) =h(\tilde{u}_{i}+u_{e},z_{i}).\]

By expansion of the scalar and vector polynomials \(f(s+u_{e},\sigma)\) and \(h(s+u_{e},\sigma)\) and \(|\tanh(\rho_{i})|\leq 1\), using Young's inequality (2.6) appropriately, it follows from the
Assumptions (1.5) and (1.6) that the new functions \(\tilde{f}\) and \(\tilde{h}\) possess the properties

\[\begin{split}&\tilde{f}(s,\sigma,\rho)s\leq-\tilde{\alpha}|s|^{4}+ \tilde{\lambda}|s||\sigma|+\tilde{J},\quad(s,\sigma,\rho)\in\mathbb{R}^{2+\ell}, \\ &\max\left\{\frac{\partial\tilde{f}}{\partial s}(s,\sigma,\rho), \left|\frac{\partial\tilde{f}}{\partial\sigma}(s,\sigma,\rho)\right|\right\} \leq\tilde{\beta},\;(s,\sigma,\rho)\in\mathbb{R}^{2+\ell},\\ &\left|\frac{\partial\tilde{f}}{\partial\rho}(s,\sigma,\rho) \right|\leq k|s+u_{e}|,\;(s,\sigma,\rho)\in\mathbb{R}^{2+\ell},\end{split} \tag{5.23}\]

and

\[\begin{split}&\tilde{h}(s,\sigma)\sigma\leq\tilde{q}|s|^{2}|\sigma| +\tilde{L}|\sigma|,\quad(s,\sigma)\in\mathbb{R}^{1+\ell},\\ &\left|\frac{\partial\tilde{h}}{\partial s}(s,\sigma)\right|\leq \tilde{\xi}(|s|+1),\;\;\frac{\partial\tilde{h}}{\partial\sigma}(s,\sigma)=0, \quad(s,\sigma)\in\mathbb{R}^{1+\ell},\end{split} \tag{5.24}\]

where the positive constants \(\tilde{\alpha},\tilde{\lambda},\tilde{J},\tilde{\beta}\) in (5.23) and \(\tilde{q},\tilde{L},\tilde{\xi}\) in (5.24) are the new parameters for the new model equations (5.21), which may also depend on the constant \(u_{e}\).

We notice the structural similarity and the new parameters between the Assumptions (1.5)-(1.6) and the properties (5.23)-(5.24) and we see the maneuverable differences in the term \(-k\tanh(\rho_{i})u_{i}\) of (5.22) in (5.23) and a spillover constant \(au_{e}\) of the third equation of (5.21).

Then we can conduct _a priori_ estimates parallel to the steps shown in Section 3 and Section 4 in the same formulated framework. It can be shown that the weak solutions of this neural network model (5.21) exist globally in time and the solution semiflow has an absorbing set \(\tilde{B}^{*}\) in the same state space \(E\). Specifically, we have

\[\tilde{B}^{*}=\{\tilde{g}\in E:\|\tilde{g}\|^{2}\leq K^{*}\}\]

where \(g=\operatorname{col}\left(u_{1},z_{1},\rho_{1},\cdots,u_{m},z_{m},\rho_{m}\right)\) and

\[K^{*}=1+\frac{2C_{4}m}{\tilde{\mu}\min\{C_{3},1\}}|\Omega|, \tag{5.25}\]

in which

\[\begin{split}&\tilde{\mu}=\min\left\{\frac{2}{C_{3}},\;\gamma, \;b\right\},\quad C_{3}=\frac{1}{\tilde{\alpha}}\left(1+\frac{\tilde{q}^{2}} {4\gamma}\right),\\ & C_{4}=C_{3}\tilde{J}+\frac{\tilde{L}^{2}}{\gamma}+\frac{a^{2}u _{e}}{b}+\frac{1}{4}\left(\frac{C_{3}^{2}\,\tilde{\lambda}^{2}}{\gamma}+\frac {a^{2}}{2b}+1\right)^{2}.\end{split} \tag{5.26}\]
Moreover, the ultimate bound property holds:

\[\limsup_{t\to\infty}\,\sum_{i=1}^{m}\|\tilde{u}_{i}(t)\|_{L^{4}}^{4}<Q^{*}\]

where

\[Q^{*}=1+\frac{4\tilde{\lambda}^{2}}{\tilde{\alpha}}K^{*}+m\left(\tilde{\alpha}+ \frac{4\tilde{J}^{2}}{\tilde{\alpha}}\right)|\Omega|. \tag{5.27}\]

We can also get the pointwise ultimate bound:

\[\limsup_{t\to\infty}\,\sum_{i=1}^{m}|\tilde{u}_{i}(t,x)|_{\mathbb{R}}<G^{*}\]

where

\[G^{*}=\left[1+\frac{2C_{4}m}{\tilde{\mu}\min\{C_{3},1\}}\right]^{1/2}. \tag{5.28}\]

Finally we define the neuron difference functions: For \(1\leq i,j\leq m\,(i\neq j)\),

\[\tilde{U}_{ij}=\tilde{u}_{i}-\tilde{u}_{j},\qquad Z_{ij}=z_{i}-z_{j},\qquad R_ {ij}=\rho_{i}-\rho_{j}.\]

They satisfy the following differencing reaction-diffusion equations:

\[\begin{split}\frac{\partial\tilde{U}_{ij}}{\partial t}=& \,\eta\Delta\tilde{U}_{ij}+\tilde{f}(\tilde{u}_{i},z_{i},\rho_{i})- \tilde{f}(\tilde{u}_{j},z_{j},\rho_{j})-\sum_{\nu=1}^{m}\frac{P(\tilde{u}_{i}- \tilde{u}_{j})}{1+\exp\left[-r(\tilde{u}_{\nu}-\tilde{V})\right]}\\ \frac{\partial Z_{ij}}{\partial t}=&\,\Lambda Z_{ij} +\tilde{h}(\tilde{u}_{i},z_{i})-\tilde{h}(\tilde{u}_{j},z_{j}),\\ \frac{\partial R_{ij}}{\partial t}=&\,a\tilde{U}_{ ij}-bR_{ij}.\end{split} \tag{5.29}\]

Parallel to the steps in the proof of Theorem 5.2, one can show that if the threshold condition (5.20) is satisfied, where the threshold constant

\[\begin{split}&\Psi=\frac{1+\exp\left[r(G^{*}+|\tilde{V}|)\right]}{m }\times\\ &\left[\tilde{\beta}+\frac{\tilde{\beta}^{2}+\tilde{\xi}^{2}}{ \gamma}+k+\frac{a^{2}}{b}+2\sqrt{2Q^{*}}C^{*}\left(\frac{k^{2}}{b}+\frac{2 \tilde{\xi}^{2}}{\gamma}\right)+\frac{64(Q^{*})^{2}(C^{*})^{4}}{\eta^{3}} \left[\frac{k^{2}}{b}+\frac{2\tilde{\xi}^{2}}{\gamma}\right]^{4}\right]\end{split} \tag{5.30}\]

and the mathematical coefficient \(C^{*}\) remains the same as in Theorem 5.2, then the solutions of this memristice neural network \(\mathbb{NW}\) is exponentially synchronized in the state space \(E\) at a uniform convergence rate.
## 6. **Examples and Numerical Simulation**

In this section we shall provide two typical and most useful mathematical models of biological neural networks with memristors to illustrate the applications of the achieved exponential synchronization result in Theorem 5.2.

To avoid notational overlap or confusion, the parameters in the following two subsections will be attached with subscript \(1\) and subscript \(2\), respectively.

Numerical simulation for these two types of memristive neural networks will be performed to show the synchronization convergence behavior with a relatively higher threshold and possibly lower convergence rate due to the nonlinear weak coupling of the solution trajectories in the depicted curves of their \(L^{2}\)-norms.

### Diffusive Hindmarsh-Rose Equations with Memristor

Consider a model of memristive diffusive Hindmarsh-Rose neural networks with memristor [16, 26, 45]:

\[\begin{split}&\frac{\partial u_{i}}{\partial t}=\eta_{1}\Delta u_{i} +a_{1}u_{i}^{2}-b_{1}u_{i}^{3}+v_{i}-w_{i}-k_{1}\tanh(\rho_{i})u_{i}-Pu_{i} \sum_{j=1}^{m}\Gamma(u_{j}),\\ &\frac{\partial v_{i}}{\partial t}=\alpha_{1}-\beta_{1}u_{i}^{2}- v_{i},\\ &\frac{\partial w_{i}}{\partial t}=q_{1}u_{i}-r_{1}w_{i},\\ &\frac{\partial\rho_{i}}{\partial t}=c_{1}u_{i}-\delta_{1}\rho_{ i},\end{split} \tag{6.1}\]

for \(t>0,\ x\in\Omega\subset\mathbb{R}^{n}\ (n\leq 3)\), where \(1\leq i\leq m\) and \(\Omega\) is a bounded domain up to three dimension with locally Lipschitz continuous boundary. The nonlinear function \(\Gamma(s)\) is the same as in (1.2).

In this system (6.1), the variable \(u_{i}(t,x)\) refers to the membrane electric potential of a neuron cell, the variable \(v_{i}(t,x)\) represents the transport rate of the ions of sodium and potassium through the fast channels and can be called the spiking variable, while the variables \(w_{i}(t,x)\) called the bursting variable represents the transport rate across the neuron membrane through slow channels of calcium and some other ions. All the involved parameters \(a_{1},b_{1},c_{1},\eta_{1},\alpha_{1},\beta_{1},q_{1},r_{1},\delta_{1},k_{1}\) and the coupling strength coefficient \(P\) can be any positive constants.

We impose the homogeneous Neumann boundary conditions for the \(u\)-component, \(\frac{\partial u}{\partial\nu}\left(t,x\right)=0\), \(x\in\partial\Omega\), and the initial conditions of the components are denoted by

\[u_{i}^{0}(x)=u_{i}(0,x),\ v_{i}^{0}(x)=v_{i}(0,x),\ w_{i}^{0}(x)=w_{i}(0,x),\ \rho_{i}^{0}=\rho_{i}(0,x),\ \ 1\leq i\leq m.\]

For illustrating the synchronization result Theorem 5.2, we can simply check all the Assumptions in (1.5) and (1.6) are satisfied by this model of memristive Hindmarsh-Rose equations (6.1). The vector functions \(z_{i}(t,x)\) in the general model (1.1) in this
case is

\[z_{i}(t,x)=\begin{pmatrix}v_{i}(t,x)\\ w_{i}(t,x)\end{pmatrix}\]

and correspondingly the vector \(\sigma=\operatorname{col}\left(\sigma_{v},\sigma_{w}\right)\) has two components.

Verify the Assumptions (1.5) and (1.6): In this model, we have the scalar function

\[f(s,\sigma)=a_{1}s^{2}-b_{1}s^{3}+\sigma_{v}-\sigma_{w},\]

the 2-dimensional square matrix and the vector function

\[\Lambda=\begin{pmatrix}-I&0\\ 0&-r_{1}I.\end{pmatrix},\quad h(s,\sigma)=\begin{pmatrix}\alpha_{1}-\beta_{1}s ^{2}\\ q_{1}s\end{pmatrix}.\]

We can verify that

\[\begin{split} f(s,\sigma)s&=s(a_{1}s^{2}-b_{1}s^{3}+\sigma_{v}- \sigma_{w})=a_{1}s^{3}-b_{1}s^{4}+s(\sigma_{v}-\sigma_{w})\\ &\leq\left(\frac{3b_{1}}{4}|s|^{4}+\frac{a_{1}^{4}}{4b_{1}^{3}} \right)-b_{1}|s|^{4}+|s|(|\sigma_{v}|+|\sigma_{w}|)\\ &\leq-\frac{b_{1}}{4}|s|^{4}+\sqrt{2}|s||\sigma|+\frac{a_{1}^{4}}{4b _{1}^{3}},\quad\text{for}\,\,\,(s,\sigma)\in\mathbb{R}^{3},\end{split} \tag{6.2}\]

and

\[\begin{split}&\max\left\{\frac{\partial f}{\partial s}(s,\sigma ),\left|\frac{\partial f}{\partial\sigma}(s,\sigma)\right|\right\}=\max\left\{ 2a_{1}s-3b_{1}s^{2},1\right\}\\ &\leq\max\,\left\{\frac{a_{1}^{2}}{2b_{1}}+(2b_{!}-3b_{1})s^{2},1 \right\}\leq\max\,\left\{\frac{a_{1}^{2}}{2b_{1}},\,1\right\},\quad\text{for }\,\,(s,\sigma)\in\mathbb{R}^{3}.\end{split} \tag{6.3}\]

Therefore Assumption (1.5) is satisfied. Moreover,

\[\begin{split}\langle\Lambda\sigma,\sigma\rangle&=- \sigma_{v}^{2}-r_{1}\sigma_{w}^{2}\leq-\min\left\{1,\,r_{1}\right\}|\sigma|^{2 },\quad\text{for}\,\,\,\sigma\in\mathbb{R}^{2},\\ h(s,\sigma)\sigma&=(\alpha_{1}-\beta_{1}s^{2}) \sigma_{v}+q_{1}s\sigma_{w}\leq\beta_{1}s^{2}|\sigma_{v}|+q_{1}s^{2}|\sigma_ {w}|+(\alpha_{1}|\sigma_{v}|+q_{1}|\sigma_{w}|)\\ &\leq(\beta_{1}+q_{1})s^{2}|\sigma|+(\alpha_{1}+q_{1})|\sigma|, \quad\text{for}\,\,\,(s,\sigma)\in\mathbb{R}^{3},\end{split} \tag{6.4}\]

and

\[\begin{split}&\frac{\partial h}{\partial s}(s,\sigma)\leq| \operatorname{col}\left(-2\beta_{1}s,\,q_{1}\right)|\leq\max\left\{2\beta_{1},q _{1}\right\}(|s|+1),\\ &\frac{\partial h}{\partial\sigma}(s,\sigma)=0,\quad\text{for} \,\,\,(s,\sigma)\in\mathbb{R}^{3}.\end{split} \tag{6.5}\]
Therefore Assumption (1.6) is also satisfied. We can record the specific parameters in (1.5) and (1.6) for this memristive Hindmarsh-Rose neural network model as follows:

\[\begin{gathered}\alpha=\frac{b_{1}^{4}}{4},\quad\lambda=\sqrt{2}, \quad J=\frac{a_{1}^{4}}{4b^{3}},\quad\beta=\max\,\left\{\frac{a_{1}^{2}}{2b_{1 }},\,\sqrt{2}\right\},\\ \gamma=\max\,\left\{1,\,r_{1}\right\},\quad q=\beta_{1}+q_{1}, \quad L=\alpha_{1}+q_{1},\quad\xi=\max\,\{2\beta_{1},q_{1}\}.\end{gathered} \tag{6.6}\]

Apply the proved synchronization Theorem 5.2 to this memristive diffusive Hindmarsh-Rose neural network model (6.1). Then we reach the following result.

**Theorem 6.1**.: _For memristive diffusive Hindmarsh-Rose neural networks with the model (6.1), if the threshold condition (5.2) with the parameters in (6.6) is satisfied by the coupling strength coefficient \(P\), then the neural network is exponentially synchronized in the state space \(E=[L^{2}(\Omega,\mathbb{R}^{4})]^{m}\) at a uniform exponential convergence rate \(\delta(P)\) shown in (5.3) with the parameters given in (6.6)._

We numerically solve the differential equations of the memristive Hindmarsh-Rose neural network with the model (6.1) in a two-dimensional square domain. And we use the finite difference method for the numerical scheme programmed in Python.

Choose the following parameters in the model (6.1):

\[\begin{gathered} m=4,\ \ \eta_{2}=5,\ \ a_{1}=1,\ \ b_{1}=2,\,k_{1}=0.3,\ \ V=0.5,\ \ r=0.1,\\ \alpha_{1}=0.4,\ \ \beta_{1}=0.06,\ \ q_{1}=0.2,\ \ r_{1}=4,\ \ c_{1}=1,\ \ \delta_{1}=7.\end{gathered}\]

Take time-step to be \(0.00025\) and spatial-step to be \(1\) on a \(32*32\) membrane. We compute and show the \(L^{2}\)-norm curves of the neuron membrane potential variable \(u_{i}\), the spiking variable \(v_{i}\), the bursting variable \(w_{i}\) and the memductance variable \(\rho_{i}\) in Figure 1 to Figure 4. Lastly the pairwise difference of vector solutions \(g_{i},\ i=1,2,3,4\) in the energy space \(E\) is shown in Figure 5.

In Figure 1 to Figure 4, with a comparison between results after \(666\) iterations and results after \(2000\) iterations, one can observe the synchronization tendency of the four characterizing variables \((u_{i},v_{i},w_{i},\rho_{i})\) among the neurons in the simulated mimristive Hindmarsh-Rose neural network. From Figure 5, we observe that the \(L^{2}\)-norms of pairwise differences \(\|g_{i}-g_{j}\|\) tend to \(0\).

We also calculate the following key constants involved in Theorem 5.2 based on our selection of parameters, rounding up to \(2\) digits.

\[\begin{gathered} C_{1}=0.25,\quad C_{2}=0.44,\quad\mu=4.0,\quad K =3630.45,\quad Q=23719.02,\\ G=2.12,\quad C^{*}=0.4,\quad\kappa=16.69,\quad P=19.60,\quad\delta=4. 0.\end{gathered}\]

The constant \(C^{*}\) from Gargliardo-Nirenberg inequality is chosen to be \(0.4\) based on [4].

Table 1 to Table 4 list the sampled values of the four components \(u_{i},v_{i},w_{i}\), and \(\rho_{i}\) of the simulated solution \(g_{i}\) at one same point in the domain at the initial time
\(t=0\), at the 200 and 2000 time-step. It is seen that with a big difference on the initial values, after a certain time, the values of \(u_{i}\), \(v_{i}\), \(w_{i}\), and \(\rho_{i}\) tend to be close to each other between various neurons.

The synchronization result rigorously proved in this work is illustrated by the example with sample selections of the system parameters and a randomized set of initial data. Our numerical simulation also exhibits that the neuron potentials \(u_{i}\) seem to be synchronized fastest within a limited time partly due to the memristor-potential coupling, while it takes much longer time to observe the synchronization on the other three variables \(v_{i}\), \(w_{i}\) and \(\rho_{i}\).

\begin{table}
\begin{tabular}{c|c|c|c}  & Initial Value & At the 200 time step & At the 2000 time step \\ \hline \(u_{1}\) & 0.0075524 & 0.00160041 & 0.00441385 \\ \(u_{2}\) & 0.038242690 & 0.01277219 & 0.00990575 \\ \(u_{3}\) & 0.064522400 & 0.02294659 & 0.01506631 \\ \(u_{4}\) & 0.098338950 & 0.03892756 & 0.02388300 \\ \end{tabular}
\end{table}
Table 1. Comparison of the \(u_{i}\) at the point \(x=10,\ y=10\)

Figure 1. The \(L^{2}\) norm of the neurons component \(u_{i}\) after 666 iterations (upper figure) and after 2000 iterations (lower figure)
### Diffusive FitzHugh-Nagumo Equations with Memristor

Consider a model of memristive neural networks described by the diffusive FitzHugh-Nagumo equations

\begin{table}
\begin{tabular}{c|c|c|c}  & Initial Value & At the 200 time step & At the 2000 time step \\ \hline \(v_{1}\) & 0.01623901 & 0.03495750 & 0.16725151 \\ \(v_{2}\) & 0.35418338 & 0.35641666 & 0.37220956 \\ \(v_{3}\) & 0.61091695 & 0.60062472 & 0.52791116 \\ \(v_{4}\) & 1.30404545 & 1.25993808 & 0.94827405 \\ \end{tabular}
\end{table}
Table 2. Comparison of the \(v_{i}\) at the point \(x=10,\ y=10\)

\begin{table}
\begin{tabular}{c|c|c|c}  & Initial Value & At the 200 time step & At the 2000 time step \\ \hline \(w_{1}\) & 0.00289204 & 0.00239707 & 0.00053399 \\ \(w_{2}\) & 0.05901589 & 0.04850523 & 0.00841518 \\ \(w_{3}\) & 0.08612489 & 0.07083853 & 0.01235578 \\ \(w_{4}\) & 0.14579938 & 0.11989174 & 0.02088769 \\ \end{tabular}
\end{table}
Table 3. Comparison of the \(w_{i}\) at the point \(x=10,\ y=10\)

Figure 2. The \(L^{2}\) norm of the neurons component \(v_{i}\) after 666 iterations (upper figure) and after 2000 iterations (lower figure)
cf. [12, 33, 47] with nonlinear weak coupling:

\[\frac{\partial u_{i}}{\partial t} =\eta_{2}\Delta u_{i}+\alpha_{2}u_{i}(u_{i}-\beta_{2})(1-u_{i})- \gamma_{2}\,w_{i}-k_{2}\tanh(\rho_{i})u_{i}-Pu_{i}\sum_{j=1}^{m}\Gamma(u_{j}),\] \[\frac{\partial w_{i}}{\partial t} =a_{2}u_{i}+c_{2}-b_{2}w_{i},\] \[\frac{\partial\rho_{i}}{\partial t} =q_{2}u_{i}-r_{2}\rho_{i}, \tag{6.7}\]

\begin{table}
\begin{tabular}{c|c|c|c}  & Initial Value & At the 200 time step & At the 2000 time step \\ \hline \(\rho_{1}\) & 0.00516599 & 0.00377392 & 0.00065141 \\ \(\rho_{2}\) & 0.05720875 & 0.04118323 & 0.00309032 \\ \(\rho_{3}\) & 0.09721523 & 0.07001327 & 0.00511662 \\ \(\rho_{4}\) & 0.13192033 & 0.09538789 & 0.00755887 \\ \end{tabular}
\end{table}
Table 4. Comparison of the \(\rho_{i}\) at the point \(x=10,\ y=10\)

Figure 3. The \(L^{2}\) norm of the neurons component \(w_{i}\) after 666 iterations (upper figure) and after 2000 iterations (lower figure)
for \(t>0,\,x\in\Omega\subset\mathbb{R}^{n}\) (\(n\leq 3\)), where \(1\leq i\leq m\) and \(\Omega\) is a bounded domain with locally Lipschitz continuous boundary \(\partial\Omega\). All the involved parameters are positive constants. The nonlinear function \(\Gamma(s)\) is the same as in (1.2).

In this system, the fast excitatory variable \(u_{i}(t,x)\) refers to the transmembrane electrical potential of a neuron cell and the slow recovering variable \(w_{i}(t,x)\) represents the integrated ionic current across the neuron membrane. The memductance \(\rho_{i}(t,x)\) of the memristor caused by the electromagnetic induction flux across the neuron membrane is a scalar function. We impose the homogeneous Neumann boundary condition is \(\frac{\partial u_{i}}{\partial\nu}(t,x)=0,\,t>0,\,x\in\partial\Omega,1\leq i\leq m\), and the initial states of the system are denoted by

\[u_{i}^{0}(x)=u_{i}(0,x),\quad w_{i}^{0}(x)=w_{i}(0,x),\quad\rho_{i}^{0}=\rho_{ i}(0,x),\quad 1\leq i\leq m.\]

As an application of the synchronization result Theorem 5.2, here we just check all the Assumptions in (1.5) and (1.6) are satisfied by this model of memristive FitzHugh-Nagumo neural networks (6.7).

Figure 4. The \(L^{2}\) norm of the neurons component \(\rho_{i}\) after 666 iterations (upper figure) and after 2000 iterations (lower figure)
In this case the generic functions in the Assumptions (1.5) and (1.6) are

\[f(s,\sigma) =\alpha_{2}s(s-\beta_{2})(1-s)-\gamma_{2}\,\sigma, \tag{6.8}\] \[\Lambda =-b_{2},\] \[h(s,\sigma) =a_{2}s+c_{2},\quad(s,\sigma)\in\mathbb{R}^{2}.\]

Check the Assumptions (1.5) and (1.6): We can verify that

\[f(s,\sigma)s =-\alpha_{2}s^{4}+\alpha_{2}(1+\beta_{2})s^{3}-\alpha_{2}\beta_{ 2}s^{2}-\gamma_{2}\,s\,\sigma \tag{6.9}\] \[\leq-\alpha_{2}\left(s^{4}-\frac{3}{4}s^{4}-\frac{1}{4}(1+\beta_{ 2})^{4}\right)+\gamma_{2}|s||\sigma|\] \[=-\frac{1}{4}\alpha_{2}s^{4}+\gamma_{2}|s||\sigma|+\frac{1}{4} \alpha_{2}(1+\beta_{2})^{4},\quad(s,\sigma)\in\mathbb{R}^{2},\]

and

Figure 5. The \(L^{2}\) norm of pairwise differences between neural network solutions after 666 iterations (upper figure) and after 2000 iterations (lower figure)
\[\begin{split}&\max\left\{\frac{\partial f}{\partial s}(s,\sigma), \left|\frac{\partial f}{\partial\sigma}(s,\sigma)\right|\right\}=\max\,\left\{-3 \alpha_{2}s^{2}+2\alpha_{2}(1+\beta_{2})s-\alpha_{2}\beta_{2},\,\gamma_{2} \right\}\\ \leq&\max\,\left\{-3\alpha_{2}s^{2}+\alpha_{2}s^{2}+ (1+\beta_{2})^{2}-\alpha_{2}\beta_{2},\,\gamma_{2}\right\}<\max\,\left\{(1+ \beta_{2})^{2},\,\gamma_{2}\right\}.\end{split} \tag{6.10}\]

Therefore the Assumption (1.5) is satisfied. Moreover,

\[\begin{split}&\langle\Lambda\sigma,\,\sigma\rangle=-b_{2}|\sigma |^{2},\\ & h(s,\sigma)\sigma=(a_{2}s+c_{2})\sigma\leq\frac{1}{4}a_{2}|s|^{ 2}|\sigma|+(a_{2}+c_{2})|\sigma|,\\ &\frac{\partial h}{\partial s}(s,\sigma)=a_{2},\quad\frac{ \partial h}{\partial\sigma}(s,\sigma)=0,\end{split} \tag{6.11}\]

for \((s,\sigma)\in\mathbb{R}^{2}\). Therefore Assumption (1.6) is also satisfied. We record the specific parameters in (1.5) and (1.6) for this memristive FitzHugh-Nagumo neural network model as follows:

\[\begin{split}\alpha=\frac{1}{4}\alpha_{2},&\lambda =\gamma_{2},\quad J=\frac{1}{4}\alpha_{2}(1+\beta_{2})^{4},\quad\beta=\max\, \left\{(1+\beta_{2})^{2},\,\gamma_{2}\right\},\\ &\gamma=b_{2},\quad q=\frac{1}{4}a_{2},\quad L=a_{2}+c_{2},\quad \xi=a_{2}.\end{split} \tag{6.12}\]

Apply the proved synchronization Theorem 5.2 to this memristive diffusive FitzHugh-Nagumo neural network model. We also reach the following result.

**Theorem 6.2**.: _For memristive diffusive FitzHugh-Nagumo neural networks with the model (6.7), if the threshold condition (5.2) with the parameters in (6.12) is satisfied by the coupling strength coefficient \(P\), then the neural network is exponentially synchronized in the state space \(E=[L^{2}(\Omega,\mathbb{R}^{3})]^{m}\) at a uniform exponential convergence rate \(\delta(P)\) shown in (5.3) with the parameters given in (6.12)._

Now we numerically solve the differential equations of the memristive FitzHugh-Nagumo neural network model (6.7) in a two-dimensional square domain. We use the finite difference method for the numerical scheme and programmed in Python.

Choose the following parameters in the model (6.7):

\[\begin{split} m=4,\,\,\eta_{2}=10,\,\,\alpha_{2}=0.5,\beta_{2}=0. 1\gamma_{2}=0.05,\,\,k_{2}=0.1,\\ a_{2}=0.3,\,\,b_{2}=3,\,\,c_{2}=1,\,\,q_{2}=0.2,\,\,r_{2}=10,\\ V=0.5,\,\,\,\,r=0.1.\end{split}\]

Take the time-step to be \(0.00025\) and spatial-step to be \(1\) on a \(32*32\) membrane. We compute the \(L^{2}\) norm of the neuron membrane potential \(u_{i}\), the recovering variable \(w_{i}\), the memductance \(\rho_{i}\), and also the vector solution \(g_{i}\) of the model equations (6.7) in the energy space \(E\). The plotted curves are shown in Figure 6 to Figure 9.
In comparison between the results after 333 iterations and after 1000 iterations in Figure 6 to Figure 8, one can observe the synchronization tendency of the three characterizing variables \((u_{i},w_{i},\rho_{i})\) among the neurons in the simulated mimristive FitzHugh-Nagumo neural network. From Figure 9, we observe that the pairwise differences of the \(L^{2}\)-norm \(\|g_{i}-g_{j}\|\) tend to 0.

We can calculate the following constants involved in Theorem 5.2 based on our selection of parameters, rounding up to 2 digits.

\[C_{1}=8.01,\quad C_{2}=2.89,\quad\mu=0.25,\quad K=94714.73,\quad Q =15101.69,\] \[G=9.67,\quad C^{*}=0.4,\quad\kappa=15.49,\quad P=19.58,\quad\delta =3.\]

The constant \(C^{*}\) from Gargliardo-Nirenberg inequality is chosen to be 0.4 based on [4].

Table 5 to Table 7 list the sampled values of the three components \(u_{i},w_{i},\) and \(\rho_{i}\) of the simulated solution \(g_{i}\) at one same point in the domain at \(t=0\), at the 100 and 10000 time-step. It is seen that with a big difference on the initial values, after a certain time, the curves of \(u_{i}\), \(w_{i}\), and \(\rho_{i}\) tend to be close to each other among various neurons.

Figure 6. The \(L^{2}\) norm of the neurons component \(u_{i}\) after 333 iterations (upper figure) and after 1000 iterations (lower figure)
The synchronization result rigorously proved in this work is illustrated by the presented example with sample selections of the system parameters and a randomized

\begin{table}
\begin{tabular}{c|c|c|c}  & Initial Value & At the 100 time step & At the 1000 time step \\ \hline \(w_{1}\) & 0.00280643 & 0.02672859 & 0.17726598 \\ \(w_{2}\) & 0.02656693 & 0.04889201 & 0.18859465 \\ \(w_{3}\) & 0.01219399 & 0.03569609 & 0.18192625 \\ \(w_{4}\) & 0.09637967 & 0.11393263 & 0.22179753 \\ \end{tabular}
\end{table}
Table 6. Comparison of the \(w_{i}\) at the point \(x=10,\ y=10\)

Figure 7. The \(L^{2}\) norm of the neurons component \(w_{i}\) after 333 iterations (upper figure) and after 1000 iterations (lower figure)
set of initial data. Our numerical simulation also exhibits that the neuron potentials \(u_{i}\) seem to be synchronized fastest within a limited time, while it takes much longer time to observe the synchronization on the other two variables \(w_{i}\) and \(\rho_{i}\).

This observation actually enhances the neurodynamical conjecture that adding a nonlinear memristor coupling in the neuron potential equation would accelerate the synchronization for the main variable of neuron membrane potential. On the other hand, it also hints that although the main Theorem 5.2 confirmed the exponential synchronization has a uniform but maybe small convergence rate, each of the three components may have a different synchronization rate, which turns out to be a new interesting problem for further research.

\begin{table}
\begin{tabular}{c|c|c|c}  & Initial Value & At the 100 time step & At the 1000 time step \\ \hline \(\rho_{1}\) & 0.00822771 & 0.02672859 & 0.00067401 \\ \(\rho_{2}\) & 0.01216352 & 0.04889201 & 0.00101167 \\ \(\rho_{3}\) & 0.08205714 & 0.03569609 & 0.00674807 \\ \(\rho_{4}\) & 0.08914909 & 0.11393263 & 0.00734463 \\ \end{tabular}
\end{table}
Table 7. Comparison of the \(\rho_{i}\) at the point \(x=10,\ y=10\)

Figure 8. The \(L^{2}\) norm of the neurons component \(\rho_{i}\) after 333 iterations (upper figure) and after 1000 iterations (lower figure)
**Conclusions**. We summarize the new contributions of the results in this paper.

1. In this paper we propose and study a general mathematical framework that can cover many typical and useful partial-ordinary differential equation models to characterize spatiotemporal dynamics of biological neural networks with memristors and weak synaptic coupling, which is a challenging and open problem in mathematical neuroscience and potentially in complex artificial learning dynamics.

2. The advancing contributions of this work are in three aspects.

First it is proved in Section 3 and Section 4 that the solution semiflow of these memristive neural networks exhibits dissipative dynamics in common and admits ultimately uniform bounds in multiple norms.

Second and more important is the exponential synchronization Theorem 5.2 and Theorem 5.3, in which we rigorously proved an explicit threshold condition in terms of the involved biological parameters and one mathematical parameter to ensure a synchronization at a uniform exponential rate in the \(L^{2}\) energy norm.

Third we provide an effective analytic approach and a significant methodology to pursue the synchronization investigation through scaled _a priori_ estimates, leverage of dynamic integral inequalities, and sharp interpolation inequalities (such as the

Figure 9. The \(L^{2}\) norm of pairwise differences between neural network solutions after 333 iterations (upper figure) and after 1000 iterations (lower figure)
crucial Gagliardo-Nirenberg inequalities on Sobolev spaces) to tackle and control the memristive effect and nonlinearity by the weak synaptic coupling strength only.

3. Two illustrative applications of the main result on synchronization are presented by the memristive diffusive Hindmarsh-Rose neural networks and FitzHugh-Nagumo neural networks.

It is expected that the mathematical framework and approach presented in this work and related computational simulations can be further generalized to a broader field and integrated with more applications in neurodynamics and network dynamics.
