Blinking Networks of Memristor Oscillatory Circuits in the Flux-Charge Domain

Multistability phenomena and complex nonlinear dynamics in memristor oscillators pave the way to obtain efficient solutions to optimization problems by means of novel computational architectures based on the interconnection of single–device oscillators. It is well-known that topological properties of interconnections permit to control synchronization and spatio–temporal patterns in oscillatory networks. When the interconnections can change in time with a given probability to connect two oscillators, the whole network acts as a complex network with blinking couplings. The work of has shown that a particular class of blinking complex networks are able to completely synchronize in a faster fashion with respect to other coupling strategies. This work focuses on the specific class of blinking complex networks made of Memristor–based Oscillatory Circuits (MOCs). By exploiting the recent Flux–Charge Analysis Method, we make clear that synchronization phenomena in blinking networks of memristor oscillators having stochastic couplings, i.e., Blinking Memristor Oscillatory Networks (BMONs), correspond to global periodic oscillations on invariant manifolds and the effect of a blinking link is to shift the nonlinear dynamics through the infinite (invariant) manifolds. Numerical simulations performed on MOCs prove that synchronization phenomena can be controlled just by changing the coupling amongst them.


INTRODUCTION
Synchronization of complex networks with interacting units is an active research topic with applications in many fields such as engineering, computer sciences, neural networks (both biological and not), neuromorphic circuit computing and information security. Arenas et al. (2008) formulated a comprehensive review in 2008, describing the impacts of this field of research. One of the most challenging problem in complex network is the study of spatio-temporal patterns emerging from synchronization phenomena and their relationship with the network topology and interactions' dynamics among the (both oscillatory and/or chaotic) network units (Pecora and Carroll, 1998;Boccaletti et al., 2006).
In literature great attention has been paid to investigate synchronization of complex networks with fixed topology and constant couplings. Only in recent years researchers have begun to focus on time-varying coupling strengths, due to the greater resemblance of this type of networks to real world systems. Different approaches have been adopted to study synchronization properties of these complex dynamical networks. For instance, Wu and Chen (2008) have provided several criteria of global synchronization, in the case of linearly coupled networks with continuous time-varying coupling. Following a complex networked control systems approach, De Lellis et al. (2008) and DeLellis et al. (2009) have considered time-varying adaptive couplings that evolve with respect to the difference among the oscillators. Several other authors (Berner et al., 2019a,b;Menara et al., 2019;Feketa et al., 2020) have focused on cluster synchronization in networks of phase oscillators with adapting coupling based on the dynamical change of the coupling strength between the oscillators. Moreover, Berner et al. (2019a,b) explored how such complex systems behave by "splitting the networks" into different clusters of nonlinear oscillators grouped in such a way that all the oscillators in a poll have the same frequency. For the scope of this work it is also necessary to take into account the dynamics of networks formed by memristors, in particular the work by Ascoli et al. (2015). In their work the authors extensively studied the synchronization dynamics between two memristor-coupled Hindmarsh-Rose oscillatory neural cells through theoretical investigation and numerical verification. On the other hand, by exploiting the mathematical tools from the control theory of switched systems, Stilwell et al. (2006) have studied a network of oscillators with switching topology, in order to model networks where some connections may fail or may create over time. In their work they have demonstrated that complex network can synchronize even when the topology is time-varying and instantaneously disconnected. Their thesis stands, though, on two pillars that serve as a condicio sine qua non: first the different circuits must be able to synchronize in a static time-average fashion; secondly the transition between two different topology must be sufficiently fast. Synchronization criteria in the case of complex network with arbitrary switching topology and strategies for designing switching signals to obtain synchronization have been investigated by Zhao et al. (2009).
A particular class of complex networks with switching topology permits to introduce the concept of blinking complex network Hasler and Belykh, 2005;Hasler et al., 2013a,b) in which the topology evolves in time according to a stochastic process. In particular, blinking networks have the peculiar feature that, starting by a fixed topology (called pristine network), new "long-range links" between two network units are randomly added according to a fixed probability. Once the new topological configuration is set, it will remain constant for a given amount of time, then again new links are randomly created, and so on.  have demonstrated that such time-varying random switching (i.e., blinking) topology induces complete synchronization in the networks. By creating randomly selected couplings, the chance of creating interconnected shortcuts increases, lowering the synchronization threshold sensibly.
The study and the development of blinking complex networks finds its importance especially in boosting functional properties of bio-inspired and neuromorphic circuits acting on the topology of synaptic connections. Neuromorphic circuits made of (1 or 2-d) arrays of silicon neurons exploits its analog nature, among other architectural characteristics, to enhance their computational efficiency and speed, but one of the main drawback is the necessity of a very high number of wired connections between the computational elements (e.g., large numbers of synapses is also used to overcome low resolution in synaptic weights-see also Chicca et al., 2014). Strategies that involve solely the reduction of the required connections, for instance reshaping the network topology adding a sum cell which has inputs and outputs that are not connected to all cells has been exploited by Seiler and Nossek (1993). This approach has permitted to reduce the number of connections for a Winner Take All architecture, however it did not take into account that most of the connections had to cross a wide part of the circuit, thus still rendering it not feasible. A more effective approach to this is to generate long distance connections, activated and deactivated randomly, in such a way that the network computations performance is the same as that of a corresponding locally-connected non-switching (averaged) network. The successful use of blinking/switching networks for solving this kind of issue has been proven, for example in a similar case to the one introduced by Seiler et al., by means of a Winner Take All locally-connected architecture (Devoe and Devoe, 2012). These results highlight the importance of the study of the synchronization dynamics of blinking networks since they can be one representative of the feasibility factors of real world neuromorphic computation.
This work focuses on Blinking Memristor Oscillatory Networks (BMONs) made of memristor-based nonlinear oscillators interacting via blinking couplings. A thorough study of nonlinear dynamics and bifurcation phenomena in memristor oscillators have been recently carried out by means of the Flux-Charge Analysis Method (FCAM). Such circuit methodology permits to show that the state space in the (voltage-current)domain [i.e., (v − i)-domain] can be decomposed into invariant manifolds, where the dynamical behavior of the memristor oscillator may correspond to different complex attractors. In addition, bifurcation phenomena can be induced either by parameters on a fixed manifold (as in the classical theory of dynamical systems) or by changing the invariant manifold but with fixed parameters (a.k.a. bifurcation without parameters).
The ultimate goal of the manuscript is to show that blinking couplings in complex networks of memristor oscillators induce bifurcations without parameters so that the whole dynamics of a BMON takes place on an invariant manifold where periodic oscillations correspond to synchronized states. In this manuscript analysis of nonlinear dynamics and numerical simulations are focused on a specific BMON in which each oscillator is a third-order dynamical system describing the Chua's circuit with a memristor replacing the nonlinear resistor, but the FCAM provides the theoretical framework to grasp the role of network topology in BMONs including Memristor-based Oscillator Circuit (MOCs) of any order. It is also worth noting that bifurcation without parameters, and then the shift of nonlinear dynamics from one manifold to another, can be obtained by applying suitable charges or flux sources via time-varying voltage and/or current pulses with finite time duration (Corinto and Forti, 2017b). This work shows that BMONs exhibit analogous nonlinear phenomena without any external input. Moreover, we show that for memristor oscillators not having all the same initial conditions it is possible to change the invariant manifold on which the nonlinear dynamics takes place only by acting on the coupling topology. Our numerical investigations of BMONs shows how the creation of time-varying long-range random connections results into switching between invariant manifolds and then into synchronization.
At first, a brief introduction to the memristor oscillators is presented showing the internal dynamics and the numerical formulations of such systems in the flux-charge domain. Secondly it will be shown how the change of the topology of interacting memristor oscillators affects invariant manifolds on which nonlinear dynamics occur. Finally, synchronization phenomena in BMONs are analyzed via numerical simulations.

PROBLEM SETUP AND BACKGROUND
Let us briefly recall some results obtained in Corinto and Forti (2017a) for a MOC in the flux-charge domain. The nonlinear oscillator consists in a modified version of the well-known Chua's circuit where the nonlinear resistor is replaced by a memristor. The cited paper details and extends the results obtained in Corinto et al. (2011) For further details regarding the circuit model and its analysis please refer to Corinto and Forti (2017a). The State Equations (SEs) in dimensionless form 1 are the following (t ≥ t 0 ): with the initial conditions where ϕ M (t 0 ) is the initial flux across the memristor, ϕ L (t 0 ) is the initial flux across the inductor, whereas q C 2 (t 0 ) is the initial charges of capacitor 2. Moreover, we have and Readers interested in the circuit-theoretic definitions of the state variables, parameters and circuit analysis can refer to the study reported in Corinto and Forti (2017a). Hereinafter and throughout the paper we set the following initial states and parameters: As shown in Corinto and Forti (2017a), the state space of the MOC in the (v, i)-domain can be decomposed in infinitely many three-dimensional manifolds M(Q 0 ) defined, for any Q 0 ∈ R as Furthermore, these manifolds are positive invariant for the dynamics of the MOC in the (v, i)-domain. Therefore, at instant t = t 0 the system is on the manifold M(Q 0 ), with Q 0 = Q(w(t 0 )) and there it stays. The value of Q(w(t 0 )) depends on the memristor initial conditions and the circuit parameters.

PROOF.
According to the definition of f in (2), it is easy to see by using (4) and the definition of Q in (5), we deduce: By extending (1), let us consider a one-dimensional array of N diffusively 2 coupled MOCs described in the flux-charge domain by (i = 1, . . . , N): where N i is the sphere of influence of ith MOC and d ik = α R ik with d ik = d ki . Following the same procedure for X 0 in (3), the parameters (4) permit to derive Analogously as in the single MOC, it is possible to prove that the state space of the Network of Memristor-based Oscillatory Circuits (NMOCs) in the (v, i)-domain is foliated by an infinite number of 3N-dimensional manifolds defined, for any From (7) and (8), we easily deduce that In the next sections we exploit NMOCs described by (6) and (7) in the flux-charge domain to make clear the influence of the couplings d ik on synchronization phenomena. In particular, by assuming that d ik describes blinking links added to a pristine complex network (i.e., an underlying networks with fixed topology) then the NMOC acts as a BMON.

COUPLING EFFECTS ON INVARIANT MANIFOLDS EMBEDDING NONLINEAR DYNAMICS
Nonlinear dynamics and global attractors of a single MOC described by (1) with X 0 in (3) can be controlled by applying suitable time-varying voltage and/or current pulses with finite time duration (Corinto and Forti, 2017b). The concept of bifurcation without parameters is the crucial tool to enable the programming of nonlinear dynamics via invariant manifolds. For instance, a sequence of period-doubling bifurcations (without parameters) leading from a periodic oscillation to a chaotic attractor can be induced by means of pulses applied to the MOC; each pulse has the effect of shifting the trajectory (due to the change of X 0 ) from an invariant manifold on which the MOC exhibits a periodic behavior to an invariant manifold with a double-scroll chaotic attractor (Corinto and Forti, 2017b). In other words, bifurcations without parameters make possible to tune nonlinear dynamics of a MOC by selecting via pulses an invariant manifold embedding all the attractors. This section shows how nonlinear dynamics and invariant manifolds of NMOCs described by (6) and (7) are influenced by the coupling terms d ik (i.e., by the network topology).
First of all, it is possible to see that when all the memristors have the same initial conditions, the coupling terms d ik play a role on the dynamics of the network [since they are present on the network Equations (6)], but they do not influence the manifold embedding the dynamics of the whole NMOC.
PROPOSITION 2. Let us consider the NMOC (6). If all the initial conditions of memristors are the same [i.e., ϕ M,i (t 0 ) = ϕ 0 , ∀i = 1, . . . , N], then the attractors of the NMOC belong to an invariant manifold that is independent on the coupling parameters d ik .
PROOF. It is easy to notice that when the initial conditions of each memristor are the same, that is ϕ M,i (t 0 ) = ϕ 0 for all i = 1, . . . , N, then with Q defined as in (5). Therefore, the invariant manifold M c (Q 0 ) is independent on the coupling. Thus, in order to induce bifurcations without parameters by means of the couplings and then change the invariant manifold embedding the attractors of the NMOC, at least one initial condition of the memristors has to be different. Without loss of generality, let us suppose that the memristor initial condition in the first oscillator is different with respect to the others. Therefore, we can prove the following result: PROPOSITION 3. Let us consider the NMOC (6) with the memristor initial conditions: Let us suppose that the global dynamics of the NMOC evolves on the manifold M c (Q 0 ), with Q(w c (t 0 )) = Q 0 = (Q 0,1 , . . . , Q 0,N ) defined as in (8).
Moreover, let us suppose to add a link between the nodes i 1 and j 1 , that are not yet linked. Therefore, we have the following two cases: • if i 1 = 1 and j 1 = 1, as in Figure 1B, the NMOCs Equations (6) are modified but the dynamics of the network still evolves on the manifold M c (Q 0 ) • if i 1 = 1, as in Figure 1C, not only the NMOCs Equations (6) change as in the previous case, but the nonlinear dynamics of the NMOC shift from the invariant manifold M c (Q 0 ) to the PROOF. Let us suppose: In this case, the NMOCs Equation (6) become: Moreover, according to (8), for any t ≥ t 0 the global dynamics of the NMOC (11) evolves on the manifold M c (Q 0 ), with Q(w c (t 0 )) = Q 0 = (Q 0,1 , . . . , Q 0,N ) and For simplicity, we suppose a sphere of influence of radius one and zero boundary conditions. The memristor initial state of the first MOC is different with respect of the other ones. When we add a link between node i 1 and node j 1 = 1, we have two possibilities: either (B) i 1 = 1 or (C) i 1 = 1.
(13) Let us suppose now to add a link between the nodes i 1 and j 1 , that are not yet linked. We can have two different cases (see Figure 1): 1. if i 1 = 1 and j 1 = 1, as in Figure 1B, the NMOCs Equations (6) are modified. Indeed, for the x-components of the MOCs with indexes i 1 and j 1 we have additional coupling terms due to the new link: while the other equations do not change.
It is important to notice that, since ϕ M,i 1 (t 0 ) = ϕ M,j 1 (t 0 ), the addition of the new link does not change the expression of X c 0,1 1 and X c 0,j 1 . Therefore, for all i = 1, . . . , N, the X c 0,i terms are still the same, and from the relation between X c 0,i and Q 0,i in (9), we can deduce that the network still evolves on the manifold M c (Q 0 ), 2. if i 1 = 1, as in Figure 1C, not only the NMOCs Equations (6) change as in the previous case, but also the coefficients X c 0,1 and X c 0,j 1 . In fact, we have the new coefficients: This implies that the nonlinear dynamics of the NMOC shift from the invariant manifold while Q ′ 0,i = Q 0,i does not change for all i different from 1 and j 1 .
The next section shows the application of this analysis to the numerical study of synchronization phenomena in NMOCs and BMONs.

SPATIO-TEMPORAL SYNCHRONIZATION PHENOMENA
This section presents the numerical analysis of spatio-temporal synchronization phenomena in NMOCs in two specific cases of interest: (a) the addition of just one single link in order to make clear the effect of a connection on the manifold shift; (b) the inclusion of blinking interconnections having a given probability to connect two MOCs. As preliminary step, we show the effect of the memristor initial conditions in such NMOC with fixed diffusive topology.

NMOC With Fixed Diffusive Topology
Let us consider N = 4 MOCs forming a NMOC whose SEs are given by (6) and such that each MOC has an influence sphere of radius r = 1, that is N i = {i−1, i, i+1} for i = 1, . . . , 4. Moreover, zero boundary conditions are considered. We suppose that each oscillator has the following circuit parameters: α = 9.5 and β = 15. Moreover, let us consider as the ith MOC's initial conditions: (15) Therefore, here ϕ 0 = −1.5.
Moreover, with this choice of initial conditions, it follows that Hence, the nonlinear dynamics of the NMOC is not on the 3Ndimensional zero-manifold M c (0), but from (3) we deduce that the manifold of the NMOC is M c (Q 0 ), with: In the following, let us consider a fixed topology defined by a diffusive (space-invariant) couplings with d ik = d = 0.05α. In this case, Q 0 is the following: (18) Thus, due to the choice of different memristor initial conditions and the coupling (diffusive topology), the global evolution of the NMOC is not on the zero-manifold and the first three oscillators synchronize nearly in-phase , while the last one is still chaotic (see Figure 2).

Addition of One Link
We are interested in investigating what happens when we add a link between the i 1 -th and the j 1 MOCs of the NMOC having the fixed diffusive topology analyzed in the previous section. For this purpose, we consider the NMOC with fixed diffusive topology from t = 0 to t = 500, then we add a link at t = 500 with d i 1 ,j 1 = d and analyze the dynamic behavior of the whole NMOC till t = 1, 000. Indeed, according to section 3, the SEs of the NMOC including the additional link change as in (14).
In the choice of i 1 and j 1 , we have here three possibilities: (a) i 1 = 1 and j 1 = 3: since these two MOCs have different memristor initial conditions, at t = 500 the nonlinear dynamics of the whole NMOC shifts from the invariant In this case, the addition of the link does not modify the synchronization properties of the network (see Figure 3). (b) i 1 = 1 and j 1 = 4: also in this case the memristor initial conditions are different, therefore the dynamic evolution of the NMOC for any t ≥ 500 takes place on the invariant As represented in Figure 4, the new link makes the fourth MOC to synchronize with the others, and then the whole NMOC exhibits in-phase synchronization among the memristor-based oscillator circuits. (c) i 1 = 2 and j 1 = 4: in this case the two MOCs have identical memristor initial conditions, therefore the network remains on the manifold M c (Q 0 ). As it is possible to see in Figure 5, the addition of the link gives rise to a (global) periodic state corresponding to a in-phase synchronization among the oscillators in the NMOC.
The three cases (a)-(c) makes clear that adding suitable links in NMOCs with fixed topology can induce synchronization phenomena. In particular, for the case under study the adding of a link connecting the fourth MOC is essential to obtain in-phase synchronization because it was out-of-phase in the MNMOC with fixed diffusive couplings. A more general case is considered in the next section by investigating larger NMOCs with fixed diffusive topology and where new links are added with a given probability and for a fixed time. We refer to such case as BMONs.

Blinking Memristor Oscillatory Networks (BMONs)
Several works in literature (Berner et al., 2019a,b;Menara et al., 2019;Feketa et al., 2020) have investigated networks of phase oscillators with adapting coupling based on the dynamical change of the coupling strength between the oscillators in order to obtain cluster synchronization. In the following we adopt the approach of , Hasler and Belykh (2005), and Hasler et al. (2013a,b) and we consider Blinking Memristor Oscillatory Networks, that is networks with a fixed topology (the so-called pristine network) and where new long-range links are randomly added and removed at each time period τ . The Blinking Memristor Oscillatory Networks is defined as follows (i = 1, . . . , N): Following , Hasler and Belykh (2005), Hasler et al. (2013a,b), once divided the time axis into intervals of length τ , we consider the coupling coefficients as follows: N), d = 0.075α, for all i, and s m ik equal to 1 with probability p and equal to 0 with probability 1 − p. Also in this case, zero-boundary conditions are considered. Therefore, it means that the corresponding pristine network (whose links are always present for each t) is a diffusive MOC network where each oscillator has a sphere of influence of radius 1. In addition, new long-range links are added with a probability p and during time intervals of length τ . In the following, we consider a switching probability of p = 0.03 and τ = 0.1. It is possible to see, as shown also , that adding long-range links improve the synchronizations properties of the network. Here we have considered N = 30. Moreover, we suppose that each oscillator has the following circuit parameters: α = 9.5 and β = 15. As in the previous case, we suppose that all the oscillators have identical initial conditions except the first one: Thus, ϕ 0 = −1.5. The numerical simulations of the pristine (i.e., fixed topology) diffusive network, reported in Figure 6, show that the network dynamics is such that the majority of the oscillators are on the large limit cycle surrounding the double-scroll except few of them. Moreover, the oscillators are not synchronized, as it is possible to see in the phase portraits of Figure 7. It is interesting to notice that in this case the last MON is on the large limit cycle while in the previous diffusive network of four MOCs, the last one was chaotic. Since the boundary conditions for both the cases are the same, we presume that the different behavior is due to the number of oscillators in the two networks. Figure 8 shows a realization of the BMON (19). Thanks to the addition of time-dependent on-off long-range links the blinking network, unlike the pristine one, exhibits in-phase synchronization.
Let us investigate how the blinking topology affects the invariant manifold on which the BMON evolves. According to (17), the pristine network at time t = 0 is on the manifold M c (Q 0 ) with (21) As seen in section 3, if a new link is added between two oscillators with different memristor initial conditions, the whole BMON changes of invariant manifold on which nonlinear attractors occur. Therefore, in our BMON, at each time we add a new long-range link, depending if it involves the first MOC or not, the whole BMON can switch from one manifold to another. In Figure 9, the components Q i (t) = Q i (w c (t)) (i = 1, . . . , 30) of Q(w c (t)) as function of time are represented. It is worth noting that these components are piece-wise constant since each manifold is positive invariant, so the system stays on one manifold till a new link (between MOCs with different memristor initial conditions) is added or removed. Clearly, the component Q 1 is the one who switches the most, since it changes each time we add a link with the first MOC as edge destination. On the contrary, for i ≥ 3, the components Q i change only if the new links involve the first MOC and the ith MOC. As we have pointed out in section 3, if the two edges of the new link are different from 1, the BMON evolves on the same invariant manifold. Finally, the component Q 2 does not change at all, since the link from the first and the second MOC (the only one who could make change this component) is already present in the pristine network. All the possible new links that involve the second MOC are among it and a j-th oscillator (with j ≥ 4), but, since these MOCs have the same memristor initial conditions, the manifold embedding the evolution of the BMON does not change.
For further information we have also considered the average system whose SEs are the same as (19), with (for i = k)    (22) It is interesting to notice that in this case the average system does not exhibit in-phase synchronization because of the last MOC (see Figure 10). Indeed, the fact that the average BMON synchronizes is a sufficient but not necessary condition in order to have the synchronization of the blinking network . Moreover, in our case, due to the different memristor initial conditions, the MOCs are not identical as in . Thus, as stated in the beginning of this work, we show that for memristor oscillators not having all the same initial conditions it is possible to change the invariant manifold on which the nonlinear dynamics takes place only by acting on the coupling topology. The numerical procedures that exploit the FCAM analysis method, of BMONs shows how the creation of time-varying long-range random connections results into switching between invariant manifolds which leads them to synchronization. Our results show that nonlinear dynamical effects generated by the blinking action of couplings and the consequent switching among invariant manifolds is needed in order to achieve a in-phase synchronization for the whole BMON.

CONCLUSION
In this paper we have considered a network of memristor nonlinear oscillators in the Flux-Charge Analysis Method (FCAM) framework, in order to show that it is possible to tune nonlinear dynamics of complex networks of memristor-based oscillators by simply adding new links among the oscillators, thus by acting on the blinking coupling topology. This has been done under the hypothesis of different memristor initial conditions. The FCAM permits to interpret the dynamical evolution of the whole blinking network of memristor oscillators as the periodic/chaotic attractors occurring on different invariant manifolds. Numerical simulation show that the synchronization properties can be improved by adding blinking long-range links and consequently by randomly switching of the invariant manifold.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
VL and FC supervised and developed the mathematical model described in the paper. JS performed the simulations. FC wrote abstract and section 1. VL and JS wrote sections 3-7. All authors contributed in the analysis of the results and in writing the manuscript.

FUNDING
This work was supported by the Ministero degli Affari Esteri e della Cooperazione Internazionale (MAECI) under the project n. PGR00823.