Brain Connectivity Dynamics and Mittag–Leffler Synchronization in Asymmetric Complex Networks for a Class of Coupled Nonlinear Fractional-Order Memristive Neural Network System with Coupling Boundary Conditions

: This paper investigates the long-time behavior of fractional-order complex memristive neural networks in order to analyze the synchronization of both anatomical and functional brain networks, for predicting therapy response, and ensuring safe diagnostic and treatments of neurological disorder (such as epilepsy, Alzheimer’s disease, or Parkinson’s disease). A new mathematical brain connectivity model, taking into account the memory characteristics of neurons and their past history, the heterogeneity of brain tissue, and the local anisotropy of cell diffusion, is proposed. This developed model, which depends on topology, interactions, and local dynamics, is a set of coupled nonlinear Caputo fractional reaction–diffusion equations, in the shape of a fractional-order ODE coupled with a set of time fractional-order PDEs, interacting via an asymmetric complex network. In order to introduce into the model the connection structure between neurons (or brain regions), the graph theory, in which the discrete Laplacian matrix of the communication graph plays a fundamental role, is considered. The existence of an absorbing set in state spaces for system is discussed, and then the dissipative dynamics result, with absorbing sets, is proved. Finally, some Mittag–Leffler synchronization results are established for this complex memristive neural network under certain threshold values of coupling forces, memristive weight coefficients, and diffusion coefficients.


Introduction and Mathematical Setting of the Problem
Complex networks of dynamical systems appear naturally in real-world systems such as biology, physics (e.g., plasma, laser cooling), intelligent grid technologies (e.g., power grid networks, communications networks), neuromorphic engineering, social networks, as well as neuronal networks.Of particular interest are complex memristive neural networks in the brain network.The macroscopic anatomical brain network, which is composed of a large number of neurons (≈10 11 ) and their connections (≈10 15 ), is a complex largescale network system that exhibits various subsystems at different spatial scales (micro and macro) and timescales, yet is capable of integrated real-time performance.These subsystems are huge networks of neurons, which are connected with each other and are modified based on the activation of neurons.The communication, via a combination of electrical and chemical signals between neurons, occurs at small gaps called synapses (see, e.g., [1]).This brain network structure implements mechanisms of regulation at different scales (from microscopic to macroscopic scales).On the microscopic level, neural plasticity regulates the formation and behavior of synaptic connections between individual neurons in response to new information.The association matrix of pair-wise synaptic weights between neurons will be of the order of 10 100 .In view of the scale of the network and the specific difficulties of accurately measuring all synaptic connections, an accurate description of a complete network diagram of brain connections is an ongoing challenge and a task of great difficulty (see, e.g., [2]).
However, the rise of functional neuroimaging and related neuroimaging techniques, such as Electroencephalography (EEG), Magnetoencephalography (MEG), functional Magnetic Resonance Imaging (fMRI), diffusion-weighted Magnetic Resonance Imaging (dMRI), and Transcranial Magnetic Stimulation (TMS), has led to good mapping and deeper understanding of the large-scale network organization of the human brain.The fMRI technique can be used to determine which regions of the brain are in charge of crucial functions, and assess the consequences of stroke, malignant brain tumors, abscess, other structural alterations, or direct brain therapy.For this last one, see, e.g., [3] for the used of online sensor information to automatically adjust, in real time, the brain tumor drugs through MRIcompatible catheters, via nonlinear model-based control techniques and a mathematical model describing tumor-normal cell interaction dynamics.
The EEG and MEG signals measure, respectively, electrical activity and magnetic fields induced by the electrical activity, from various brain regions with a high temporal resolution (but with limited spatial coverage).However, fMRI measures whole-brain activity indirectly (by detecting changes associated with blood flow for each network, over a specified interval of time) with a high spatial resolution (but with limited temporal resolution).Consequently, in order to improve both spatial and temporal resolution, EEG and MEG signals are often associated with fMRI (see, e.g., [4]).
This whole-brain connectomics approach, which relies on macroscopic measurements of structural and functional connectivity, has notably favored a major development in the identification and analysis of effects of brain injury or neurodegenerative and psychiatric diseases on brain systems related to cognition and behavior, for better diagnosis and treatments (see, e.g., [5][6][7]).
The dynamic interaction between neuronal networks and systems, which takes into account the dynamic flows of information that pass through different interconnected, widely distributed, and functionally specialized brain regions, is crucial for normal brain function (see, e.g., [8,9]).Measuring electrical activity (from, e.g., dMRI connectivity, magnetoencephalography or electroencephalography recordings) has allowed researchers to point out the existence of oscillations characterized by their frequency, amplitude, and phase (see, e.g., [8,[10][11][12]).This phenomenon is considered to result from oscillatory neuronal (local) synchronization and long-range cortical synchronization (which is linked to many cognitive and memory functions).Moreover, several studies have established that the activity pattern of cortical neurons depends on the history of electrical activities (e.g., caused by the long-range interaction of ionic conductances), as a result of changes in synaptic strength or shape of synaptic plasticity (see, e.g., [13][14][15]).In addition, the diffusion terms play an important role in dynamics and stability of neural networks when, e.g., electrons are moving in asymmetric electromagnetic fields (see, e.g., [16][17][18]).So, diffusion phenomena cannot be neglected.Understanding mechanisms behind these synchronized oscillations and their alterations is very important for better diagnosis and treatments of neurological disorders.In particular, the relationship between stability of synchronization and graph theory is established.This relationship characterizes the impact of network topology on the disturbances.Disturbances or alteration of such synchronized networks, taking into account memory characteristics of neurons and their past history, play an important role in several brain disorders, such as neuropsychiatric diseases, epileptic seizures, Alzheimer's disease, and Parkinson's disease (see, e.g., [19][20][21][22][23] and references therein).
Moreover, noninvasive brain stimulation is attracting considerable attention due to its potential for safe and effective modulation of brain network dynamics.In particular, in the context of human cognition and behavior, the targeting of cortical oscillations by brain stimulation with electromagnetic waveforms has been widely used, whether it is to ensure a safe treatment, to improve quality of life, or to understand and explain the contribution The concept of a memristor, which is a passive and nonlinear circuit element, was first introduced by Chua [25].The author estimates that the memristor was to be considered as basic as the three classical passive electronic elements: the resistor, the capacitor, and the inductor.In the resistor, there is a relation between the instantaneous value of voltage and the instantaneous value of current.Unlike the resistor, the memristance depends on the amount of charge passing through it and consequently, the memristor can remember its past dynamic history.It is a natural nonvolatile memory.
Memristor-based neural network models can be divided into (see, e.g., [26,27]) (a) memristive synaptic (autapse or synapse) models, in which a memristor is used as a variable connection weight for neurons; (b) neural network models affected by electromagnetic radiation, in which a memristor is used to simulate the electromagnetic induction effects.The memristive synaptic model uses the bionic properties of the memristor to realistically mimic biological synaptic functions, while the neural network model under electromagnetic radiation employs a magnetic flux-controlled memristor to imitate the electromagnetic induction effect on membrane potential.
Recently, considerable efforts have been made to mimic significant neural behaviors of the human brain by means of artificial neural networks.Due to the attractive property of this new type of information storage and processing device, it can be implemented using synaptic weights in artificial neural networks.It can also be an ideal component for simulating neural synapses in brain networks due to its nano-scale size, fast switching, passive nature, and remarkable memory characteristics (see, e.g., [28][29][30] and references therein).In recent years, the dynamical characteristics of memiristor-based neural networks have been extensively analyzed and, in this context, several studies have been reported as regards chaos, passivity, stability, and synchronization (see, e.g., [31][32][33][34][35][36][37] and references therein).
Currently, fractional calculus is particularly efficient, when compared to the classical integer-order models, for describing the long memory and hereditary behaviors of many complex systems.Fractional-order differential systems have been studied by many researchers in recent years; the genesis of fractional-order derivatives dates back to Leibniz.Since the beginning of the 19th century, many authors have addressed this problem and have devoted their attention to developing several fractional operators to represent nonlin-ear and nonlocal phenomena (such as Riemann, Liouville, Hadmard, Littlewood, Hilfer, and Caputo among others).Fractional integrals and fractional derivatives have proved to be useful in many real-world applications; in particular, they appear naturally in a wide variety of biological and physical models (see, for instance, Refs.[38][39][40][41][42][43] and references therein).Both the Riemann-Liouville and Caputo derivative operators are the most common and widely used way of defining fractional calculus.Unlike this Riemann-Liouville operator, when solving a fractional differential equation, we use the Caputo fractional operator [44], for which there is no need to define the fractional-order initial conditions.One of the most important characteristics of fractional operators is their nonlocal nature, which accounts for the infinite memory and hereditary properties of underlying phenomena.Recently, in [39], we proposed and analyzed a mathematical model of the electrical activity in the heart through the torso, which takes into account cardiac memory phenomena (this phenomenon, also known as the Chatterjee phenomenon, can cause dynamical instabilities (as alternans) and give rise to highly complex behavior including oscillations and chaos).We have shown numerically the interest of modeling memory through fractional-order derivatives, and that, with this model, we are able to analyze the influence of memory on some electrical properties, such as the duration of action potentials (APD), action potential morphology, and spontaneous activity.
On the other hand, synchronization of neural networks plays a significant role in the activities of different brain regions.In addition, compared with the concept of stability, the synchronization mechanism (with possible control), for two or more apparently independent systems, is being paid increasing attention in neuroscience research and medical science, because of its practicality.The study of dynamical systems and the synchronization of biological neural networks, with different types of coupling, have attracted a large amount of theoretical research, and consequently, the literature in this field is very extensive, especially in the context of integer-order differential systems.Since the literature on this topic has been receiving a significant amount of attention, it is not our intention to comment in detail here on all the works cited.For a general presentation of the synchronization phenomenon and its mathematical analysis, we can cite, e.g., [45][46][47][48].Concerning problems associated with integer-order partial differential equations various methods and technique, we can refer to, e.g., [49][50][51][52][53][54].Finally, for problems associated with fractional-order partial differential equations, various methods have recently been studied in the literature, such as, for example, [55], which considers the synchronization control of a neural network's equation via Pinning Control, ref. [56], which investigates the dissipativity and synchronization control of memristive neural networks via feedback controller, and [57], which explores the stability and pinning synchronization of a memristive neural network's equation.
Motivated by the above discussions, to take into account noninvasive brain stimulation and the effect of memory in the propagation of brain waves, together with other critical brain material parameters, we propose and analyze a new mathematical model of fractionalorder memristor-based neural networks with multidimensional reaction-diffusion terms, by combining memristor with fractional-order neural networks.The models with time fractional-order derivatives integrate all the past activities of neurons and are capable of capturing the long-term history-dependent functional activity in a network.The diffusion can be seen as a local connection (at a lower scale), like, e.g., in [58], whereas the coupling topology relates to dynamical properties of network dynamical systems, corresponding to physical or functional connections at the upper scale.
Thus, the derived brain neural network model is precisely the system (1)-( 5) (see further), which is a nonlinear coupled reaction-diffusion system in the shape of a set of fractional-order differential equations coupled with a set of fractional-order partial differential equations (interacting via a complex network).
In the present work, we are interested in the synchronization phenomenon in a whole network of diffusively nonlinear coupled systems which combines past and present interactions.First, we will impose initial data on a closed and bounded spatial domain, and analyze some complex dynamical property of the long-time behavior of a derived fractional-order large-scale neural network model.After a rigorous investigation of dissipative dynamics, different synchronization problems of the developed complex dynamical networks are studied.

Formulation of Memristive-Based Neural Network Problem
We shall consider a network of m coupled neurons denoted by N G = {N i : i = 1, 2, . . ., m}, where the network size m ≥ 2 is a positive integer.Our model of a memristorautapse-based neural network with external electromagnetic radiation can be depicted as in Figure 2. In this paper, motivated by the above discussions (see Section 1), we introduce the following new fractional-order memristive neural network of coupled neurons, modeled by the following Caputo-fractional system, including the magnetic flux coupling (on each neuron N i of the network, for i = 1, m): where Q ∞ = Ω×]0, +∞[ and the spatial domain Ω is a bounded open subset and its boundary denoted by Γ = ∂Ω is locally Lipschitz continuous.Here, ∂ α 0 + denotes the forward Caputo fractional derivative with α a real value in ]0, 1].The state variables ϕ i , u i , and w i describe, respectively, the membrane potential, the ionic variable, and the magnetic flux across the membrane of the i-th neuron N i , for i = 1, m.The term F is the nonlinear activation operator.The functions f ex and g ex represent the external forcing current and the external field electromagnetic radiation, respectively.The parameters ξ h ≥ 0, ξ f ≥ 0, ξ e ≥ 0 are the coupling strength constants with ξ e ξ f ̸ = 0.The values a and σ can be any nonzero number constants and all the parameters b, c 1 , c 2 , κ, ν 1 , and ν 2 can be any positive constants.
The fractional parameter c α is c α = κ s C α > 0, where C α is the membrane pseudocapacitance per unit area and κ s is the surface area-to-volume ratio (homogenization parameter).The membrane is assumed to be passive, so the pseudo-capacitance C α can be assumed to be constant.Moreover, since the electrical restitution curve (ERC) is affected by the action potential history through ionic memory, we have represented the memory via u (respectively, via w) by a time fractional-order dynamic term c β ∂ α 0 + u (respectively, by c γ ∂ α 0 + w), where the positive parameters c β and c γ are assumed to be constants.The fractional parameters c α , c β , and c γ depend on the fractional-order α.

Remark 1.
According to the expression of the Caputo derivative, we can obtain that the unit for the dimension of c α C is s α−1 , with s the unit for the dimension of time and C the capacitance c α for α = 1 (this result remains valid for c β and c γ ).
The functions a ij , for i, j = 1, m, represent the memristor's synaptic connection weight coefficient (an example with three neurons is depicted in Figure 3).The term ∑ m j=1 a ij (x, t)H j (ϕ j ) could be regarded as an emulation of a neurological disease, e.g., epileptic seizures, in which the nonlinear operators H j are some activation functions.We assume that (for i, j = 1, m) where (a ij , a ij ) are in IR 2 , and we denote The operator Ψ = q 1 Ψ 1 + q 2 Ψ 2 is the memory conductance of the flux-controlled memristor, where Ψ 1 (w) = δ 1 + η 11 w + η 12 w 2 , Ψ 2 (w) = δ 2 + ζ 2 tanh(w), all the parameters δ 1 , δ 2 , η 11 , η 12 , ζ 2 can be any positive constant, and q 1 q 2 ̸ = 0, with q k ≥ 0, for k = 1, 2. The magnetic flux coupling κϕΨ(w) could be regarded as an additive induction current on the membrane and represents the dynamic effect of electromagnetic induction on neurological diseases (examples with three neurons are depicted in Figure 4).For simplicity, we can write The operator G i , which contains the information of network topology, is defined by with the coupling (or connectivity) matrix G = (G ij ) 1≤i,j≤m (−G is called the Laplacian matrix of the graph), in which G ij are the coefficients of connection from the i-th to the j-th neuron, satisfying the assumption G ik , i.e., the matrix G has vanishing row and column sums and non-negative off-diagonal elements.
Then, G i (v 1 , . . ., v m ) can be written as In graph theory, the Laplacian matrix −G, also called the graph Laplacian or Kirchhoff matrix, defines the graph topology with m the number of vertices/nodes N i in graph G(N G , E G ) (with N G set of vertices and E G set of edges/links).The diagonal matrix D = − diag(G 11 , . . ., G mm ) is called the degree matrix of graph and A := D + G is called the adjacency matrix of the graph.
The state variable (ϕ i , u i , w i ), for i = 1, m in this network, is coupled with the other neurons {N j : j ̸ = i} in the network through the state equation by G i and/or through the boundary conditions as follows (fully connected network on boundary): where n is the outward normal to Γ and p f > 0, p e > 0 are the coupling strength constants on boundary.The tensors K f and K e are the effective diffusion tensors that describe the heterogeneity of brain tissue and local anisotropy of cell diffusion.The initial conditions of (1) to be specified will be denoted by (for i = 1, m) The rest of the paper is organized as follows.In the next section, we give some preliminary results useful in the sequel.In Section 4, we shall prove the existence, stability, and uniqueness of weak solutions of the derived model, under some hypotheses for data and some regularity of nonlinear operators.An important feature of the uniqueness of the solution is the semiflow physical property of the system; starting the system at time t 0 , letting it run until time s = t 0 + τ, and then restarting it and letting it run from time s to the final time of T f amounts to running the system directly from time t 0 to final time T f .In Section 5, the existence of an absorbing set in state spaces for the system is discussed, an estimate of the solutions is derived when time is large enough, and then the dissipative dynamics result, with absorbing sets, is proved.In Section 6, synchronization phenomena are discussed and some Mittag-Leffler synchronization criteria for such complex dynamical networks are established in different situations.Precisely, some sufficient conditions for the synchronization are obtained first for the complete (or identical) synchronization (which refers to the process by which two or more identical dynamical systems adjust their motion in order to converge to the same dynamical state as time approaches infinity) problem and then for the master-slave synchronization problem via appropriate pinning feedback controllers and adaptive controllers.In Section 7, conclusions are discussed.In Appendix A, the full proof of well-posedness of the derived system is shown, and in Appendix B, a brief introduction to some definitions and basic results in fractional calculus in the Rieman-Liouville sense and Caputo sense is given.

Assumptions, Notations, and Some Fundamental Inequalities
Some basic definitions, notations, fundamental inequalities, and preliminary lemmas are introduced and other results are developed.
Our study involves the following fundamental inequalities, which are repeated here for review: (i) Hölder's inequality: (ii) Young's inequality (∀a, b > 0 and ϵ > 0): ab ≤ ϵ p a p + ϵ −q/p q b q , f or p, q ∈]1, +∞[ and (iii) Minkowski's integral inequality (for p ∈]1, +∞[ and t > 0): Finally, we denote by | 0 | the Lebesgue measure of 0, by sgn(x) the sign of the scalar x, and by L(A; B) the set of linear and continuous operators from a vector space A into a vector space B. The operator R * stands for the adjoint to linear operator R between Banach spaces.
From now on, we assume that the following assumptions hold for nonlinear operators, matrix coupling, and tensor functions appearing in our model on Ω.
First we introduce the following spaces: H = L 2 (Ω) and V = H 1 (Ω) (endowed with their usual norms).We will denote by V ′ the dual of V. We have the following continuous embeddings (p ′ is such that and the injection V ⊂ H is compact, where p ≥ 2 if d = 2 and 2 ≤ p ≤ 6 if d = 3.
For tensor functions, we assume that the following assumptions hold.
(H1) We assume that the conductivity tensor functions K θ ∈ W 1,∞ (Ω), θ ∈ {e, f }, are symmetric, positive definite matrix functions and that they are uniformly elliptic, i.e., there exist constants 0 < d f ≤ r f and 0 < d e ≤ r e such that (∀v The operators F and (H i ) i=1,m , which describe behavior of the system, are supposed to satisfy the following assumptions.
(H2) The operators F and (H i ) i=1,m are Carathéodory functions from Ω × IR into IR.Furthermore, the following requirements hold.
(H2) 2 The nonlinear scalar activation functions H i are bounded with H i (0) = 0 and satisfy l i -Lipschitz condition, i.e., For the operators (G i ) i=1,m , which are defined from the matrix coupling G, we introduce first the following notations: the matrix S = (ϵ ij ) 1≤i,j≤m , where Then, the matrix S is symmetric, the matrix A is antisymmetric, and G can be represented uniquely as G = S + A. It is evident that both S and A have zero row sums (i.e., Lemma 3.For (ϕ, w, u) in L 4 (Ω) × L 2 (Ω) × L 2 (Ω) and ( f ex , g ex , ρ 0 ) ∈ (L 2 (Ω)) 3 , we have the following inequalities: where and υ > 0 (υ is any positive constant, to be chosen appropriately).
The estimates of the previous lemma are needed to construct a priori estimates (to establish the existence result).Lemma 4. For all (v i ) i=1,m ∈ IR m , we have the following relations (with and then the relation (i).For the relation (ii), we have (we proceed in a similar way as in [62]) Since δ jk = 0, we can deduce that We can deduce that m ∑ j,k=1;k̸ =j For the term I, we have that and then ϵ jk ψ 2 jk (since S is symmetric and has zero row sums). ( This completes the proof. (b) The matrix (µ ij ) i,j is symmetric and satisfies ∑ 1≤i,j≤m µ ij = 0.

Remark 4.
Let A be an arbitrary antisymmetric matrix.Then, for any vector w, we have w t Aw = 0.
To end this section, we give the following lemmas and definitions.From, e.g., [63,64], we can deduce, respectively, the following Lemmas.Lemma 5. (A generalized Gronwall's inequality) Assume γ > 0, h is a nonnegative function locally integrable on (0, T) (some T ≤ +∞) and b is a nonnegative, bounded, nondecreasing continuous function defined on [0, T).Let f be a nonnegative and locally integrable function on (0, T) with, for t ∈ (0, The used function E θ 1 ,θ 2 is the classical two-parametric Mittag-Leffler function (usually The invertibility of f (t) = E α,β (−t), t > 0 follows from the complete monotonicity property of E α,β .As shown in [65], this function is completely monotone if and only if 0 < α ≤ 1 and β ≥ α.Since E α,β (0) = Γ(β) −1 and lim Definition 2 (see, e.g., [66]).The initial-boundary value problem (1)-( 6) is called dissipative in a Banach space L(Ω) if there is a bounded set, in L(Ω), K S = B(0, r) for some positive constants r > 0 such that for any bounded subsets K I (of L(Ω)) and all initial data (ϕ 0i , u 0i , w 0i ) i=1,m belonging to K I , there exists a time T K , depending on K I , such that the corresponding solution (ϕ i (t), u i (t), w i (t)) i=1,m is contained in K S for all t ≥ T K .The set K S is called an absorbing set and r is called a radius of dissipativity.

Remark 5. (i)
We can also express the previous definition as follows.For any bounded set K I , there exists a finite time T K such that any solution started with initial condition in K I remains within the bounded ball K S for all time t ≥ T K .
(ii) The study of dissipative dynamics opens the way to the analysis of the synchronization problem.Definition 3 (see, e.g., [67]).The initial-boundary value problem (1)-( 6) is (locally) complete synchronized if  6) is locally complete Mittag-Leffler synchronized if there exist X * > 0 and λ * > 0 such that for any t > t * (∀i, j = 1, m), whenever the initial state belongs to an appropriate bounded open set.The value λ * is called Mittag-Leffler synchronization rate (or degree).

Well-Posedness of the System
This section concerns the existence and uniqueness of a weak solution to problem (1), under Lipschitz and boundedness assumptions on the non-linear operators.We now define the following bilinear forms: Under hypothesis (8), it can easily be shown that the forms A k (for k = e, f ) are symmetric coercive and continuous on V. We can now write the weak formulation of the initialboundary value problem ( 1)-( 6) (for all v, ϑ, ρ in V and a.e.t ∈ (0, T), with T > 0): The first main theorem of this paper is then the following result.
Theorem 1.Let assumptions (H1) and (H2) be fulfilled and T > 0. Assume that there exists ν2 > 0 such that Then, for any (ϕ 0i , Proof.To establish the existence of a result of a weak solution to system (1), we proceed as in [39] by applying the Faedo-Galerkin method, deriving a priori estimates, and then passing to the limit in the approximate solutions using compactness arguments.The uniqueness result can be evaluated in the classical way.The full proof is given in Appendix A.

Dissipative Dynamics of the Solution
In this section, first we prove that all the weak solutions of initial value problem (1) exist for time t ∈ [0, ∞).Then, we show that there exists an absorbing set in space Theorem 2. Under the assumptions of Theorem 1, there exists a unique global weak solution (ϕ i (t), w i (t), u i (t)) i=1,m in the space H m , for time t ∈ [0, ∞), of problem (1)- (6).Moreover, for any given ε 0 > 0, there exists an absorbing set for the semiflow for problem (1)-( 6) in space H m , which is Proof.Taking the scalar product of the equations of system (1) by ϕ i , u i , and w i , respectively, and adding these equations for i = 1, m, we obtain, according to Lemma 4 and assumption (H1): By using similar arguments to derive (A9), we can deduce where The values ϵ > 0 and υ > 0 are chosen appropriately so that λ 0 > 0 and υ 0 > 0, i.e., ϵ > 4c 2 2 bλ and υ < 5b 4ϵ .In particular, we have We can solve this Caputo fractional differential inequality (17) to obtain the following estimate in maximal existence time interval (from Lemma 6): where ) and ).Since the solution will never blow up at any finite time, then the maximal existence time interval is [0, ∞).Consequently, the solution of initial-boundary value problem ( 1)-( 6) exists in space H m , for any time t ∈ [0, ∞).Then we have the existence of solution semiflow for (1)-( 6), that is, a mapping m ] for any s, t ≥ 0.Moreover, from (18), we can deduce that for any ε 0 > 0 (in view of the asymptotic property of the Mittag-Leffer function): where and then there exists a finite T ρ , given by T ρ = (− 1 , such that the solution trajectories that started at initial time t = 0 from the set B(0, ρ) will permanently enter the set B(0, D), for all t ≥ T ρ .So, B(0, D) is an absorbing set in H m for the semiflow Σ S and this semiflow is dissipative.This completes the proof.

Synchronization Phenomena
Before analyzing the synchronization problems, we will examine the uniform boundedness of the solution in

Uniform Boundedness in K m
In this section, we prove a result on the ultimately uniform boundedness of solution (ϕ i , u i , w i ) i=1,m in K m Lemma 7.For all (v i ) i=1,m ∈ IR m , we have the following results.
For (ii), we have Since (from Young's inequality and the fact that This completes the proof. ) and the assumptions of Theorem 2 hold.If there exists a constant ν1,2 > 0 such that where Proof.As g 1 : y −→ y 3 is an increasing function, then the primitive g 1 is convex and we can have, from [68] (since g 1 is independent on time), 1 4 Consequently, by taking the scalar product of the first and third equations of system (1) by ϕ 3  i and w 3 i , respectively, and adding these equations for i = 1, m, we can deduce (according to Lemma 7) where Since ∀v ∈ IR, we have | tanh(v) |≤ 1 and v 4 ≤ 1 + v 6 ; then (by using Young's inequality), From ( 20), ( 22) and ( 23), we can deduce with C 0 > 0 to be chosen appropriately.In particular, we can deduce where ), with θ (i) 2 ̸ = 0 for an appropriate Using a similar argument to derive (A8), we can deduce that If (ϕ 0i , u 0i , w 0i ) i=1,m is in a ball B(0, ρ) of K m , then (ϕ 0i , u 0i , w 0i ) i=1,m is in a ball B(0, ρ) of H m and then, from (19), (ϕ i (t), u i (t), w i (t)) i=1,m is bounded by some constant D > 0 in H m .Consequently, where 2 |.Then, as in the proof of Theorem 2, we can deduce that if (ϕ 0i , u 0i , w 0i ) i=1,m is in a ball B(0, ρ) of K m , then there exists t * ≥ t 0 , such that for t ≥ t * , (ϕ i (t), u i (t), w i (t)) i=1,m is in a ball B(0, D ρ ) in K m .Therefore, B(0, D ρ ) is an absorbing set in K m for (1)-( 6) (and then (1)-( 6) is a dissipative dynamical system in K m ).The proof is complete.
In the sequel, we assume that (ϕ 0i , u 0i , w 0i ) i=1,m is in a bounded set B(0, ρ) of K m , and then, from Theorem 3, there exists t * > t 0 , such that for t ≥ t * , (ϕ i (t), u i (t), w i (t)) i=1,m is in a ball B(0, D ρ ) of K m , where D ρ > 0 depend on some given parameters of the problem.

Local Complete Synchronization
In this section, we assume that ξ h = 0 and we consider the local synchronization solutions of (1), whether the synchronous state is robust to perturbations, whenever the initial conditions belong to some appropriate open and bounded set.
with the boundary conditions Theorem 4.Under the assumptions of Theorem 3, if there exist appropriate constants ℵ ij > 0 and ℶ ij > 0, i, j = 1, m with i ̸ = j, such that where and 0 < γ 0,ν < 1 (which can be appropriately and arbitrarily selected), then the response system (1) is local complete (Mittag-Leffler) synchronized in H m at a uniform Mittag-Leffler rate.
Proof.Multiply (28) by the function (Φ ij , U ij , W ij ) and integrate the resulting system over all of Ω to obtain (according to (29) and Lemma 4) Since we have (because, for any (y, z) | tanh(y) |≤ 1 and sech 2 (y) ≤ 1), according to assumption (H1): Then (from Young's inequality), where 0 < γ 0,ν < 1 is any positive parameter (to be chosen appropriately).Consequently (by summing the first and the second inequalities of (33)), We take ϵ such that −ϵ α .From Gagliardo-Nirenberg inequalities, we have that there exists where Thus, the previous relations with (34) yield the following inequalities: By summing (for all 1 ≤ i, j ≤ m), we can deduce that (according to Lemma 4) By adding the above two inequalities, we can deduce that From the Poincaré-Steklov inequality, we can deduce that and then, Since, from (30), there exist appropriate constants ℵ ij > 0 and ℶ ij > 0 such that it holds that (for t ≥ t * ) where Consequently (from Lemma 6), where λ 3 is the Mittag-Leffler synchronization rate and X * is given by Finally, Corollary 1.We assume that there exists appropriate positive constants ℵ ij and ℶ ij , for i, j = 1, m with i ̸ = j, such that where and 0 < γ 0,ν < 1 is any positive parameter (to be chosen appropriately).Then, the response system (1) is local complete (Mittag-Leffler) synchronized in H m at a uniform Mittag-Leffler rate.

Master-Slave Synchronization via Pinning Control
The goal of pinning control is to synchronize the whole of the memristor-based neural network by controlling a select part of neurons of the network.Without losing any generality, we may assume that the first q (1 < q < m) neurons would be pinning.Controlled synchronization refers to a case when the synchronization phenomenon is artificially induced by using a suitably designed control law.In order to explore the synchronization behavior via pinning control, we introduce the corresponding slave system to the master system (1) by (for i = 1, m) where, for i, j = 1, m, Ξ i are reasonable controllers and ãij satisfies the assumption (2).The initial and boundary conditions of ( 38) are ( φi , ũi , wi )(0) = ( φ0i , ũ0i , wi0 ), on Ω. ( Introduce now the symmetric matrix P = (p ij ) 1≤i,j≤m such that p ij = 1 if i ̸ = j and Remark 6.We can prove easily that (for Set Thus, from ( 38), ( 39), (1), ( 5), and (6), we can deduce that (for i = 1, m) with the initial and boundary conditions of (38) (ϕ i (0), u i (0), w i (0)) = (ϕ 0i , u 0i , w i0 ), on Ω. (43)

Feedback Control
In this section, we assume that the control is a feedback law.
Combining the first and second inequalities of (49), we obtain By taking ϵ such that −ϵ α 3 8 + c 2 1 b = 0 and by using similar arguments to derive (34), we can deduce that ϵ 6 with Adding the above inequality, we can write From ( 41), we can deduce that According to Poincaré-Steklov, we obtain (using the expression of p ij ) Let us put Then (with Φ = (ϕ i ) i=1,m , W = (w i ) i=1,m ) From ( 44), we can deduce (with P ϖ = diag(ϖ 11 , . . .., ϖ 1m )) According to assumption (45), we can deduce that where Θ = min( min ).Then, by Lemma 6, we can deduce that where . Hence, the error system ( 42) is globally asymptotically stable and the system (1) and response system (38) are globally (Mittag-Leffler) synchronized in H m under the feedback controllers (Ξ i ) i=1,m (at a uniform Mittag-Leffler rate).This completes the proof of Theorem.
Corollary 2. Assume that the assumptions of Theorem 3 hold.Suppose that the pinning feedback control functions Ξ i (for i = 1, m) satisfy assumption (44).
m, and if there exist diagonal matrices N f and N w with strictly positive diagonal terms such that Ω P − D w − N w are positive-defined matrices, where and 0 < γ 0,ν < 1 any positive parameter (to be chosen appropriately), the master-slave systems (1) and (42) can achieve Mittag-Leffler synchronization via the pinning feedback control functions (Ξ i ) i=1,m .
Example 1.The pinning feedback control functions (Ξ i ) i=1,m can have, for example, the following forms.
For this example, we have if Consequently, we obtain that
Proof.According to (58) and ( 59), we can have that ∂ α t + * S(t) ≤ I (from the expression of S) with Consequently, we derive (according to (55)) and then with R φ = diag( κ1i φ 1i , . . .., κ1m φ 1m ).According to assumption (61), we can deduce that where Θ = min( min Then, by Lemma 4.3 of [33], there exists a T S > t * such that (for where Hence, the error system (42) is globally asymptotically stable and the system (1) and the response system (38) are Mittag-Leffler synchronized under the feedback controllers (Ξ i ) i=1,m .This completes the proof of the theorem.
We end this analysis by the following corollary (of Theorem 6).Corollary 3. When the assumptions of Theorem 3 are satisfied, the master-slave system can achieve Mittag-Leffler synchronization via pinning feedback controllers (58) (for all i = 1, m) and if there exist diagonal matrices N f and N w with strictly positive diagonal

Conclusions
Recent studies in neuroscience have shown the need to develop and implement reliable mathematical models in order to understand and effectively analyze the various neurological activities and disorders in the human brain.In this article, we mainly investigate the long-time behavior of the proposed fractional-order complex memristive neural network model in asymmetrically coupled networks, and the Mittag-Leffler synchronization problems for such coupled dynamical systems when different types of interactions are simultaneously present.A new mathematical brain connectivity model, taking into account the memory characteristics of neurons and their past history, the heterogeneity of brain tissue, and the local anisotropy of cell diffusion, is developed.This developed model is a set of coupled nonlinear Caputo fractional reaction-diffusion equations, in the shape of a fractional-order differential equation coupled with a set of time fractional-order partial differential equations, interacting via an asymmetric complex network.The existence and uniqueness of the weak solution as well as regularity result are established under assumptions on nonlinear terms.The existence of some absorbing sets for this model is established, and then the dissipative dynamics of the model (with absorbing sets) are shown.Finally, some synchronization problems are investigated.A few synchronization criteria are derived to obtain some Mittag-Leffler synchronizations for such complex dynamical networks.Precisely, some sufficient conditions are obtained first for the complete synchronization problem and then for the master-slave synchronization problem via appropriate pinning feedback controllers and adaptive controllers.
The developed analysis in this work can be further applied to extended impulsive models by using the approach developed in [16].Moreover, it would be interesting to extend this study to synchronization problems for fractional-order coupled dynamical networks in the presence of disturbances and multiple time-varying delays.For predicting and acting on phenomena and undesirable behavior (as opposed to a synchronous state) occurring in brain network dynamics, the established synchronization results can be applied in the study of robustness behavior of uncertain fractional-order neural network models by considering the approach developed in [61].
The future objective is to simulate and validate numerically the developed theoretical results.These studies will be the subject of a forthcoming paper.According to the stochastic nature of synapses (and then the presence of noise), it would also be interesting to investigate the stochastic processes in the brain's neural network and their impact on the synchronization of network dynamics.

Appendix A
Proof of Theorem 1.To establish the existence result of a weak solution to system (1), we proceed as in [39] by applying the Faedo-Galerkin method.
Let (ϖ k ) k≥1 be a Hilbert basis and orthogonal in L 2 of V.For all N ∈ IN * , we denote by W N = span(ϖ 1 , • • • , ϖ N ) the space generated by (ϖ k ) N≥k≥1 , and we introduce the orthogonal projector L N on the spaces W N .For each N, we would like to define the approximate solution (ϕ N i , w N i , u N i ) i=1,m of problem (1).Setting where (ϕ N,k i , w N,k i , w N,k i ) i=1,m;k=1,N are unknown functions, and replacing (ϕ i , w i , u i ) i=1,m by (ϕ N i , w N i , u N i ) i=1,m in (1), we obtain ∀l = 1, N and a.e.t ∈ (0, T), the system of Galerkin equations (∀i = 1, m): where (L N ϕ 0i , L N w 0i , L N u 0i ) satisfies (by construction) Step 1.We show first that, for every N, the system (A1) admits a local solution.The system (A1) is equivalent to an initial value for a system of nonlinear fractional differential equations for functions (ϕ N,l i , w N,l i , w N,l i ) i=1,m;l=1,N , in which the nonlinear term is a Carathéodory function.The existence of a local absolutely continuous solution on interval [0, T], with T ∈]0, T] is insured by the standard FODE theory (see, e.g., [69,70]).Thus, we have a local solution of (A1) on [0, T].
Step 2. We next derive a priori estimates for functions (ϕ N i , w N i , w N i ) i=1,m , which entail that T = T, by applying iteratively step 1.For simplicity, in the next step, we omit the "˜" on T. Now, we set where ϑ N,k , π N,k and ζ N,k are absolutely continuous coefficients.
Then, from (A1), the approximation solution satisfies the following weak formulation: Take (h N , v N , e N ) = (ϕ N i , w N i , u N i ) and add these equations for i = 1, m.We obtain, according to Lemma 4 and assumption (H1), From Lemma 3 and assumption (H2), we can deduce that where θ 0 = η 1 2η 2 , θ w = κ(η 2 θ 2 0 + ζ − δ) and υ > 0 to be chosen appropriately.Consequently (from the assumption (13)), b > 0 and υ 0 = 5b 8 − ϵυ 2 > 0, we can deduce, by summing the previous three inequalities : where and then we have (from (A7)) where In particular, we have From Lemma 6 and the uniform boundedness of (ϕ N i , u N i , w N i ) i=1,m (t = 0) (from (A2)), we can deduce that (for all t ∈ (0, T)) and then Ψ N (t) is uniformly bounded with respect to N. This ensures that, for i = 1, m, Moreover, from the inequality of (A9), relations (A32) and (A2), we have that for t ∈ (0, T) We can deduce that Now we estimate the fractional derivative ) and using uniform coercivity of forms A f and A e , we can deduce According to assumptions (H1)-(H2) and the boundedness of function tanh, we obtain from (A16) that Then, from (A12)-(A15), the regularity of (ρ 3 , f ex , g ex ), and the continuous embedding of H 1 in L 6 , we can derive (according to Young and Hölder inequalities) Moreover, since −F 1 is an increasing function, then the primitive −E F of −F 1 is convex.Hence, from [68], we can deduce (since −F 1 is independent on time) Now, we estimate the following term ) and then (according to the regularity of ρ 4 and (A12)) Consequently, According to (A12)-(A15), the continuous embedding of H 1 in L 6 , and the uniform boundedness of quantities ϕ N i (0), w N i (0), and w N i (0) in H 1 , we obtain from the third equation of (A19) that the sequence (w N i ) is uniformly bounded in L ∞ (0, T; H 1 (Ω)) and we can derive the following estimate (for all t ∈ (0, T)): This implies (since and then (using Minkowski inequality and the continuous embedding of Using (A12), (A14) and uniform boundedness of quantities u N i (0 This estimate enables us to say that In order to prove that the local solution can be extended to the whole interval (0; T), we use the following process.We suppose that a solution of (A1) on [0, T k ] has already been defined and we shall derive the local solution on [T k , T k+1 ] (where 0 < T k+1 − T k is small enough) by making use of the a priori estimates and fractional derivative ∂ α T + k with beginning point T k .So, by this iterative process, we can deduce that Faedo-Galerkin solutions are well defined on the interval (0, T).So, we omit details.
Then, by the weak convergence and Lebesgue's dominated convergence arguments, we have (A24) Consequently, φi = ∂ α 0 + ϕ i in the weak sense.In the same way, we prove, in the weak sense, that (∂ α 0 + w i , ∂ α 0 + u i ) = ( wi , ũi ).Consider ℵ ∈ D(]0, T[) and (h N , v N , e N ) ∈ W 3 N .According to (A1) , we can deduce that According to (A14), (A23) and to density properties of the space spanned by (ϖ l ), and using similar arguments as used to obtain relation (A24), passing to the limit when N goes to infinity is easy for the linear terms.For passing to the limit (N → ∞) in nonlinear terms, which requires hypothesis (H2), we can use the standard technique, which consists of taking the difference between the sequence and its limit in the form of the sum of two quantities such that the first uses the strong convergence property and the second uses the weak convergence property.So we omit details.

Appendix B
The purpose of what follows is to recall some basic definitions and results of fractional integrals and derivatives in the Riemann-Liouville sense and Caputo sense.We start from a formal level and, for a given Banach space X, we consider a sufficiently smooth function with values in X, f : t ∈ [a, +∞) → f (t) (with −∞ < a < b ≤ +∞).Each fractional-order parameter γ is assumed to be in ]0, 1].Definition A1.For each fractional-order parameter γ, the forward and backward γth-order Riemann-Liouville fractional integrals are defined, respectively, by (t ∈ [a, +∞) and b > a) where Γ(z) = ∞ 0 e τ τ z−1 dτ is the Euler Γ-function.
For each fractional-order parameters γ 1 and γ 2 , the following equality for the fractional integral holds for an L q -function f (1 ≤ q ≤ ∞).

Figure 1 .
Figure 1.Schematic design of a biological neuron.

Figure 2 .
Figure 2. Concept map of the coupled neural network.

Figure 3 .
Figure 3. Abridged general view of synaptic connection.

Figure 4 .
Figure 4. Two examples of connection topology of the neural network with three neurons (based on memristor-autapse and under electromagnetic radiation effects).

δ
jk = 0) and 2δ ii = m ∑ k=1 G ki .Then, we can now derive the following two lemmas.

Definition 4 .
whenever the initial state belongs to an appropriate bounded open set.The initial-boundary value problem (1)-(