On the analysis of a heterogeneous coupled network of memristive Chialvo neurons

We perform a numerical study on the application of electromagnetic flux on a heterogeneous network of Chialvo neurons represented by a ring-star topology. Heterogeneities are realized by introducing additive noise modulations on both the central–peripheral and the peripheral–peripheral coupling links in the topology not only varying in space but also in time. The variation in time is understood by two coupling probabilities, one for the central–peripheral connections and the other for the peripheral–peripheral connections, respectively, that update the network topology with each iteration in time. We have further reported various rich spatiotemporal patterns like two-cluster states, chimera states, coherent, and asynchronized states that arise throughout the network dynamics. We have also investigated the appearance of a special kind of asynchronization behavior called “solitary nodes” that have a wide range of applications pertaining to real-world nervous systems. In order to characterize the behavior of the nodes under the influence of these heterogeneities, we have studied two different metrics called the “cross-correlation coefficient” and the “synchronization error.” Additionally, to capture the statistical property of the network, for example, how complex the system behaves, we have also studied a measure called “sample entropy.” Various two-dimensional color-coded plots are presented in the study to exhibit how these metrics/measures behave with the variation of parameters.


Introduction
Neurons form the fundamental units of the central and peripheral nervous system and supervise the mechanisms of complex information processing and responding to stimuli, by exchange of electrical and chemical signals between them.These complex dynamical behaviors exhibited by neurons can be represented and studied with the help of dynamical systems [25,23] like ordinary differential equations and maps, leading to the science of neurodynamics.Recently, neurodynamics has become an emerging field of research and has attracted a lot of attentions [25] from mathematicians, biologists, computer scientists, to name a few.These dynamical systems-oriented models, which mimic many neuronal behaviors, have been confirmed experimentally [17] too.Examples of models that are represented as continuous dynamical systems include Hodgkin-Huxley model [22], Hindmarsh-Rose system [21], FitzHugh-Nagumo neuron system [14], etc., whereas examples for the ones represented as discrete systems involve Rulkov neuron system [45] and Chialvo neuron system [9].Scarcely any research attention is given to the study of neurodynamical models that are portrayed as discrete maps.Motivated, we focus on an improved model of a network of Chialvo neurons with heterogeneities incorporated, which we believe is a good imitation of a real-world nervous system.It is important to study the corresponding dynamical behaviors to gain an insight about how a nervous system might behave in reality.
One of the striking features exhibited by an ensemble of neurons is the phenomenon of synchronization.Synchronization is a universal concept in dynamical systems that have been studied in fields ranging from biology, to physics, to neuroscience, to economics [30,59,5,42,58].Relevant to our study in neurodynamics, a synchronized state can both pertain to a normal or abnormal cognitive state [63].In the latter case, it is of utmost importance to study and understand the dynamics of neurons that are not completely synchronized, as has been mentioned in the paper [46].Two of the interesting states that represent asynchronicity in the neural nodes are chimera states and solitary nodes.A "chimera" state [29,41,1,57,26], discovered by Kuramoto, is characterized by the simultaneous existence of coherent and incoherent nodes from a specific ensemble after a precise time.Chimera states have been reported in various natural phenomenon like sleeping patterns of animals [43,33,15], flashing of fireflies [19], and many more [68,41,32].On the other hand, nodes falling under the "solitary" regime [47,49,51,31,6,46,53] are the ones which behave completely different from the coherent nodes in a particular ensemble at a specific time.They get separated from the main typical cluster and oscillate with their own particular amplitude.Additionally, traveling waves are reported in various neuron systems [2,12].The authors have reported traveling waves in heterogeneous neuron systems and have claimed that it is responsible for many features of cortical dynamics, spiking time variability, and fluctuations in membrane potential.There have been research works concerning networks of dynamical systems in which the nodal behaviors were found to be rich.It is interesting to note that in a network of neurodynamical systems, the researchers have found chimera states, cluster synchronization, traveling waves, solitary states, etc.A variety of spatiotemporal patterns are reported in the heterogeneous network considered in this study, such as solitary states, traveling waves, synchronized state, asynchronized states, etc.
The fact that the dynamical behaviors portrayed by neurons are remarkably complex, demands the requirement of statistical tools to study and quantify their complexity.Measures like spatial average density, and entropy are effective tools to study the complexity in the dynamics of neurons.Banerjee and Petrovskii [4] applied the spatial average density approach to show that the spatial irregularity in an ecological model is indeed chaotic.Similar approach is considered in the reference [3].Entropy is an ubiquitous concept first introduced by C. E. Shannon in his revolutionary work [52].Since then, it has had a widespread application in various domains of research, including neuroscience [69,60].In another study [67] authors report the study of applying entropy to EEG data from patients and report that Alzheimer's disease could result in complexity loss.More relevant paper [16] related to this topic can also be found.Other applications of entropy in neuroscience include study of topological connectivity in neuronal networks [66,64,24], and estimation of the upper limit on the information content in the action potential of neurons [56].Motivated by these, we utilize the tool of sample entropy [44] to analyze time-series data of the spatially averaged action potential generated from our model and study the complexity of the dynamics.
Empirical evidence says that the functioning of neurons are affected by many external factors such as temperature [8], light [40], electromagnetic radiations [38,36,13] and many more.For example, signature of growth in embryo neuronal cells was observed with the application of electromagnetic flux [27].These evidences motivate us to perform a study on the effect of external electromagnetic flux in a network of Chialvo neurons.The ring-star network topology [37] serves as a promising candidate.Most recently, the dynamical effects of discrete Chialvo neuron map under the influence of electromagnetic flux have been studied [36], where after understanding a single neuron system, analysis on a ring-star network of the neurons has been conducted, making it flexible to work with all the possible scenarios of ring, star, and ring-star topologies.The network that has been mathematically set up in the work, consists of N Chialvo neurons arranged in the formation of a ring-star topology, where the central node is connected to all the N − 1 nodes with the homogeneous coupling strength µ and the peripheral nodes are connected to each of their R neighbors with a coupling strength σ.The nodes follow a periodic boundary condition, meaning that the (N − 1) th node is connected to the 2 nd node to complete the loop.Although the studies in the paper have exhibited some rich dynamics, the topology of the network is too simple to be close to a realworld connection of the nervous system.The study of these neurodynamical systems on this simple network topology gave further motivation to consider even more complex topologies to mimic the real-world connection between the neurons in the brain.Extending on the work, we have tried to tackle this issue by studying the spatiotemporal dynamics on introduction of heterogeneous and time-varying coupling strengths to the network topology of the ring-star Chialvo neurons under the influence of an electromagnetic flux.
In this paper, we mainly catalogue the presence of several interesting temporal phenomena i.e., solitary nodes, chimera, and traveling waves, to name a few, exhibited by the nodes in a heterogeneous network of Chialvo neurons.Heterogeneity is introduced both in space and time.In another study [7], the authors have reported chimera states in time-varying links of a network formed from two coupled populations of Kuramoto oscillators.To quantify what solitary states mean in our system model, we work with a metric called "cross-correlation" coefficient [65].This metric lets us decide the regime a node lies in.
The goals of this paper are as follows: 1. Improve the ring-star topology of the Chialvo neuron network introduced in the original paper [36], by inclusion of heterogeneous links between nodes that vary not only in space but also in time, depending on noise modulations and specific coupling probabilities, 2. Report the appearance of rich spatiotemporal patterns throughout the temporal dynamics of the heterogeneous network, 3. Study the emergence of an important asynchronous behavior called solitary nodes and characterize it using quantitative measures like cross-correlation coefficient and synchronization error that we set up according to our network topology, 4. Use the statistical tool called sample entropy on the time series data of spatially averaged action potential, to gain an intuition on the extent of complexity, and 5. Developing an understanding of the innate mechanism of whether and how the various measures relate by employing numerical simulations.
We organize the paper in the following way: in Sec. 2, we introduce our improved heterogeneous model for a ring-star network of Chialvo neurons, followed by setting up of the quantitative metrics like cross-correlation coefficient and synchronization error in Sec. 3. Next, in Sec.3.3 we establish how we are going to apply the measures of sample entropy and maximum Lyapunov exponents on time series data that get generated from our model dynamics.Then in Sec. 4 we show results from running various numerical experiments (plotting phase portraits, spatiotemporal patterns, recurrence plots, time series plots, several two regime color plots, bifurcation plots for synchronization, etc.), perform time series analysis, and try drawing inferences about the behavior and dynamics of our heterogeneous neuronal network model.Note that all simulations are performed in Python.Finally, in Sec. 5 we provide concluding remarks and future directions.

System Modelling
The two-dimensional discrete map proposed by Chialvo [9] in 1995 that corresponds to the single neuron dynamics is given by where x and y are the dynamical variables representing activation and recovery respectively.Two of the control parameters are a < 1 and b < 1 which indicate the time constant of recovery and the activation dependence of recovery respectively during the neuron dynamics.
Control parameters c and k 0 are the offset and the time independent additive perturbation.
Throughout this study we have kept a = 0.89, b = 0.6, c = 0.28, and k 0 = 0.04.The model was improved [36] to a three-dimensional smooth map by inclusion of electromagnetic flux, realized by a memristor.The improved system of equations is give by where M (φ) = α+3βφ 2 is the commonly used memconductance function [10] in the literature with α, β, k 1 and k 2 as the electromagnetic flux parameters.In (2.3), the k 1 x term denote the membrane potential induced changes on magnetic flux and the k 2 φ term denote the leakage of magnetic flux [36].The parameter k denotes the electromagnetic flux coupling strength with kxM (φ) as the induction current.We note that k can take both positive and negative values.Throughout the paper we have used α = 0.1, β = 0.2, k 1 = 0.1, k 2 = 0.2.We also restrict k in the range [ −1, 4].
Recently [36], a ring-star network of the Chialvo neurons under the effect of electromagnetic flux has been considered where the coupling strength is set to be homogeneous, and they do not vary with time.We further improve the ring-star model with the incorporation of heterogeneous coupling strengths µ m (n) and σ m (n).The improved version is then given by whose central node is further defined as ) having the following boundary conditions: (2.12) where n represents the n th iteration, and m = 1, . . ., N with N as the total number of neuron nodes in the system.For the sake of making the model more complex and closer to a realistic nervous system we introduce heterogeneities to the coupling strengths σ m (n) and µ m (n) both in space and time.In space, the heterogeneities are realised following the application of a noise source with a uniform distribution [48,65] given by where σ 0 and µ 0 are the mean values of the coupling strengths µ m and σ m respectively.Throughout the paper we keep σ 0 ∈ [−0.01, 0.01] and µ 0 ∈ [−0.001, 0.001].The noise sources ξ σ and ξ µ for the corresponding coupling strengths are real numbers randomly sampled from the uniform distribution [−0.001, 0.001].Finally, the D's refer to the "noise intensity" which we restrict in the range [0, 0.1].Note that we have used negative (inhibitory) coupling strengths in this study.They represent a significant proportion of neuronal connectivity in the human nervous system.The authors in the papers [62,61] have mentioned them in their course of simulations of the leaky Integrate-and-Fire (LIF) model.Also, negative coupling strengths are included via a rotational coupling matrix (See Eq. ( 2) in the paper [39]) during simulations of FitzHugh-Nagumo neuronal models.Heterogeneity in time is introduced by considering the network having time-varying links [28,11,35] depending on the two coupling probabilities P µ and P σ , which govern the update of the coupling topology with each iteration n.The probability with which the central node is connected to all the peripheral nodes at a particular n is denoted by P µ .Likewise, the probability with which the peripheral nodes are connected to their R neighboring nodes is given by P σ .We have studied the spatiotemporal dynamics of the improved topological model under the variation of the seven most important parameters which are k, σ 0 , µ 0 , D σ , D µ , P σ and P µ .

Quantitative metrics and time series analysis
To quantize the solitary nodes and determine the extent of synchronization in the system after a satisfactory number of iterations, we employ two metrics known as the cross-correlation coefficient [51,65] and the synchronization error [34].We also realise the complexity of the time series data that get generated from our simulations using a measure known as ssample entropy.

Cross-correlation coefficient
Keeping in mind that our system has a ring-star topology, we must define the coefficient in such a way that it captures the correct collective behaviour of the network dynamics.The general definition of the cross-coefficient denoted by Γ i,m is given by (3.17) The averaged cross-correlation coefficient over all the units of the network is given by, Throughout the paper, we use Γ 2,m , denoting the degree of correlation between the first peripheral node of the ring-star network and all the other nodes, including the central node.The average is calculated over time with transient dynamics removed and x(n) = x(n)− x(n) refers to the variation from the mean.Note that denotes the average over time.Everywhere in the paper we take 20000 iterations in time, of which we discard the first 10000 points for the dynamical variables, such that the transient dynamics is discarded, and thereafter perform all the calculations and simulations.We use Γ 2,m to characterize all the necessary regimes for our study.When the nodes have Γ Note that these values were selected after running numerous simulations confirming the fact that it is the local dynamics of the network of oscillators that govern the characteristics of solitary nodes from system to system.Γ → 1 implies that almost all the nodes are clustered in the coherent regime after a specific time without any transient dynamics.

Synchronization error
The averaged synchronization-error for the nodes in a system is given by where we again consider node number N = 2 as the baseline for the calculation as has been done in the case of cross-correlation coefficient, and n represents the n th iteration.Note that E → 1 implies that the nodes in the system are moving towards an asynchronous behavior, with E = 1 depicting a complete asynchrony.Similarly E → 0 implies synchronization.Like the first metric, denotes the average over time.

Sample entropy: a measure for complexity
Additionally, we perform a statistical analysis of the dynamics of our system to determine how complex it is.In order to do so, we take the spatial average of all the N nodes at a particular time n and generate a time series data for the action potential, which we utilize to calculate the sample entropy (SE).SE tells us about how complex the time series data is.A high value indicates unpredictability in the behavior, thus offering more complexity [20,18].To calculate SE from the time series data, we utilize an open-source package called nolds [50], that provides us with the function nolds.sampen().

Results
Every time we run a computer experiment, we use a pseudo-random generator to initialize the action potential x between 0 and 1. Additionally y and φ are initialized to 1.This is done for all the N = 100 nodes.In each simulation, the total number of iterations is set to 20000, of which we discard the first 10000 to remove transient dynamics.First, we briefly cover the single neuron dynamics under the fixed parameter values as mentioned in the previous sections, before moving to the analysis of the network.

Single neuron dynamics
For a single neuron, we set the parameter values to a = 0.89, b = 0.6, c = 0.28, k 0 = 0.04, α = 0.1, β = 0.2, k 1 = 0.1, and k 2 = 0.2.Additionally, we set k = −1.The corresponding phase portrait and the time series are given in Fig. 2 and Fig. 3 respectively.The sample entropy value is calculated to be 0.041, indicating quite a low complexity.Looking at the phase portrait, we observe that the attractor is a closed invariant curve.Using the QR−factorization method as was done recently [36], it can be noted that the maximum Lyapunov exponent is ∼ 0, exhibiting a quasi-periodic dynamics for the fixed set of parameter values.Keeping this in mind, we analyze the system of N = 100 such neurons arranged in ring-star topology in the next sections.Note that for the above selected local dynamical parameters a, b, c, k 0 , α =, β, k 1 , and k 2 , we perform the whole analysis throughout.In case, these parameters are set different, the dynamical behaviors of the single neuron and the network of neurons are expected to change accordingly.

Phase portraits, spatiotemporal patterns and recurrence plots
In this section, we have numerically plotted and gathered some interesting behaviors exhibited by the network under variations of different parameter values.We showcase a variety of different phase portraits, spatiotemporal patterns exhibited by the heterogeneous Chialov ring-star network in Eq. 2.6 and Eq.2.9.For each simulation, we have shown five separate plots corresponding to different responses.In both Fig. 4  The sample entropy is 0.041, which is quite low.This indicates less disorderedness in the dynamics as can also be seen from not so irregular behaviour in the time series.5. the last plot corresponds to the recurrence plot [38] comparing the distance between the final position of each node with all the other nodes in space after n iterations.
In the first, second, and fourth plots, the nodes lying in the solitary regime are denoted by red dots, the nodes with Γ 2,m ∈ [−0.15, 0.75] are denoted by green dots and the nodes in the coherent domain are denoted by blue dots.Additionally, the second node is denoted by a black dot in the first (phase portrait) plot.Fig. 4(a) displays the solitary nodes and the coherent nodes clustered in their respective domains along with one node that rests in the region where Γ 2,m ∈ [−0.15, 0.75].The blue nodes have cluttered with Γ 2,m ∈ [0.75, 1] and the solitary nodes have accumulated with Γ 2,m ∼ −0.172.We notice that the solitary nodes are distributed over the whole ensemble of the nodes.The fact that there exists almost equal number of nodes in both the clusters is evident from the tiny squares in the recurrence plot.In Fig. 4(b), we see that all the nodes have been clustered and in synchrony except two, which are solitary.Synchronization is also confirmed from the very small value of the normalized synchronization-error, i.e, E = 0.034.Here too the solitary nodes consist of Γ 2,m ∼ −0.172.In the corresponding recurrence plot we notice a deep blue region that covers almost the whole space, visually denoting the fact that the nodes are synchronized.Fig. 4(c) has the similar dynamics as that of Fig. 4(a).The solitary nodes mostly accumulate with Γ 2,m ∼ −0.167.Note that after sufficient time iterations, the blue and the green nodes can rest together in a single cluster (See the fourth plot in Figs.4(b) and (c)).As expected, their recurrence plots also look very similar to each other.In Fig. 5(a) we again see an emergence of clusters with nodes in synchrony at the solitary regime, having Γ 2,m ∼ −0.176 and the remaining clustered around Γ 2,m ∼ 0.5.Here in the recurrence plot, although we see squares denoting multi clusters in the dynamics, they are obviously bigger in area than the ones appearing in Fig. 4(a) and Fig. 4(c), due to the fact that the number of nodes in one of the clusters is much higher.The behavior in Fig. 5(b) refers to the phenomenon of "chimera" where nodes within a particular boundary (approximately in 40 ≤ m ≤ 45 and in 62 ≤ m ≤ 75) are completely asynchronous compared to the other nodes in the space which are convincingly synchronous, and they coexist [49].In Fig. 5(c), we observe the emergence of travelling waves.The recurrence plots in Fig. 5(b) and Fig. 5(c) visually support the fact that we indeed notice chimera and travelling waves, respectively.

Time series analysis
Next, we perform statistical analysis on the time-series data of the action potential x corresponding to the parameter combinations reported in Fig. 4 and Fig. 5.As mentioned in Sec.3.3, we take the spatial average of all the nodes at a particular time n, denoted by x, and this is illustrated in Fig. 6.The plots of the spatial average against time for the parameter combinations in Fig. 4 As seen in Fig. 6, the oscillations are irregular and do not appear to converge to any steady state of the system.To further our analysis on complex dynamics, we use nolds to compute the sample entropy for each of these cases and record them on each time series plot.What we clearly observe is the presence of fairly complex behaviour in all of them, giving us an intuition about the extent of disorderedness.Visually, we also notice irregular oscillatory behaviour in the firing pattern of x, not converging to any stable steady state.Note that Fig. 6(b) which corresponds to the mostly synchronized case (E = 0.034) Fig. 4(b) has the lowest value of sample entropy (SE=0.044)and thus the lowest disorderedness as compared to the other five cases.Furthermore, out of the above six cases, the highest value of synchronization error is observed in Fig. 5(c), having E = 0.136.The corresponding time series Fig. 6(a) too has a very high value of sample entropy, SE=0.139 (second highest).Statistically speaking, it can be expected that an increase in asynchrony leads to an increase in the complexity of the system dynamics.From the color plots and their corresponding E vs. SE plots in the next section, we can infer the above phenomenon.

Two dimensional color coded plots
4.4.1 Parameter space defined by (σ 0 , µ 0 ) Fig. 7 is a collection of two-dimensional color plots of Γ, E, Ns N , and SE in the parameter space defined by (σ 0 , µ 0 ).The parameter space is a 40 × 40 grid with k = −1, P µ = 1, P σ ∼ 0.666, D µ = D σ = 0.005.From the plot that depicts the normalized cross correlation coefficient (Fig. 7(a)), we observe an almost definitive bifurcation boundary within which most of the nodes lie in the coherent regime, i.e, Γ ∼ 1 and outside which the nodes behave incoherently, i.e, Γ < 0.75.A similar kind of bifurcation boundary is observed in the normalized syncronization error color plot (Fig. 7(b)).The region where Γ ∼ 1, the nodes tend to adopt a complete synchronization behavior, making E ∼ 0. Moving on to the sample entropy plot (Fig. 7(d)), interestingly we once again notice an almost similar bifurcation boundary between no complexity (SE∼ 0) and onset of complexity (SE > 0).Whenever, Γ ∼ 1 and E ∼ 0 then SE∼ 0 too.Now in the parameter region where Γ < 1, we notice a transition in the behavior of the nodes from complete synchrony to a degree of asynchrony, such that the normalized synchronization rate lies in the range 0 < E < 0.3 approximately.In an analogous mechanism the sample entropy increases too when Γ < 1 pertaining to a more disordered system dynamics.When Γ < 1, and the mean coupling strength σ 0 ∼ 0, there appears a tinge of violet region in the color plot, implying Γ ∈ (−0.15, −0.38), for which there appear solitary nodes (as evident from the Ns N color plot, Fig. 7(c)) and a further peak in the value of sample entropy (0.3 < SE < 0.75 approximately).Otherwise throughout the parameter space, there exit no solitary nodes.In Fig. 8 we have collected the values of E, Γ, and SE and plotted them against each other.Notice in Fig. 8(a) that it gives a clear inversely proportional trend for E and Γ.In Fig. 8(b), we notice that with increase in E, the sample entropy shows a fairly increasing trend, whereas in Fig. 8(c), with increase in Γ, the sample entropy decreases as expected.
4.4.2Parameter space defined by (σ 0 , D σ ) Fig. 9 is a collection of two-dimensional color plots in the parameter space defined by (σ 0 , D σ ) with µ 0 = −0.001,P µ = 1, P σ ∼ 0.666, D µ = 0.005 and k = −1.Note that for all values of D σ , when −0.01 < σ 0 < −0.005, we see that the nodes prefer to remain distributed throughout the space, with Γ ∼ (0, 0.5) (Fig. 9(a)) and behave in an asynchronous manner as evident from the violet region in the E color plot (Fig. 9(b)) denoting higher asynchronity.As soon as σ 0 > −0.005, the nodes mostly tend towards synchronization in the coherent region as can be seen from the deep yellow and black regimes respectively in the Γ and E color plots respectively, although there might be instances when the nodes do not settle to synchronity in the coherent domain.Interestingly, there exists a very few to no solitary nodes as can be seen from the Ns N color plot (Fig. 9(c)).Now looking into the SE plot (Fig. 9(d)), note that again there is a fair proportional relationship between sample entropy and asynchronicity, with the maximum complex behaviour appearing at σ 0 ∼ 0. It can be seen that with increase in the mean coupling strength σ 0 , from −0.01 to 0.01, sample entropy starts increasing, reaches a peak at around σ 0 ∼ 0 and then the value drops.In Fig. 10 we have collected the values of E, Γ, and SE and plotted them against each other.Notice in Fig. 10   indicates an inversely proportional trend for E and Γ.In Fig. 10(b), we notice that with increase in E, the sample entropy shows a fairly increasing trend, whereas in Fig. 10(c), with increase in Γ, the sample entropy decreases as expected.

4.4.3
Parameter space defined by (µ 0 , D µ ) Fig. 11 is a collection of two-dimensional color plots in the parameter space defined by (µ 0 , D µ ) with σ 0 = 0, P µ = 1, P σ ∼ 0.666, D σ = 0.005 and k = −1.We note that for all values of D µ > 0, probability of all the nodes having Γ < 0.75 is very high, as depicted by reddishyellow region in Fig. 11 (a), with −0.38 < Γ < −0.15 at around µ 0 ∼ 0, as can be seen from the emergence of solitary nodes in Fig. 11 (c) in this region.Otherwise there occurs no solitary nodes for µ 0 > 0. For µ 0 < 0, mostly there exists nodes with Γ ∼ 1 denoting nodes in the coherent domain and synchronized (Fig. 11 (b)), except a certain region of µ 0 < 0, where again solitary nodes emerge.Nodes are synchronized where Γ > 0.75.From the sample entropy plot (Fig. 11 (d)) we observe that for µ 0 > 0, there occurs high disorderdness, i.e, for the paramter combinations where γ < 0.75 and nodes are asynchronous.Again, like Γ and E plot, there exists a region µ < 0, where the sample entropy is > 0, overlapping with the region where solitary nodes emerge.In Fig. 12 we have collected the values of E, Γ, and SE and plotted them against each other.Notice in Fig. 12(a) that it clearly indicates an inversely proportional trend for E and Γ.In Fig. 12(b), we notice that with increase in E, the sample entropy shows a fairly increasing trend, whereas in Fig. 12(c), with increase in Γ, the sample entropy decreases as expected.4.4.4Parameter space defined by (P µ , P σ ) Fig. 13 is a collection of two-dimensional color plots in the parameter space defined by (P µ , P σ ) with σ 0 = 0, µ 0 = −0.001,D σ = D µ = 0.005 and k = −1.From Fig. 13(a), we see that for values close to P µ = 0, for all P σ , there exist regions where Γ < 0.75, where nodes behave asynchronously (Fig. 13(b)), having very high sample entropy value (Fig. 13(d)) and are solitary (Fig. 13(c)).Above this region of P µ where P µ < 0.5, mostly the nodes are synchronized in the coherent domain having very small value of sample entropy.As soon as P µ > 0.5, the region is randomly distributed between both types of extreme behaviors, denoted by the contrasting colored boxes in all the four color coded plots.In Fig. 14 we have collected the values of E, Γ, and SE and plotted them against each other.Notice in Fig. 14(a) that it clearly indicates an inversely proportional trend for E and Γ.In Fig. 14(b), we notice that with increase in E, the sample entropy shows a fairly increasing trend, whereas in Fig. 14(c), with increase in Γ, the sample entropy decreases as expected., we detect a space where almost all the nodes are solitary in the region k ∈ (−0.5, 0.5), σ 0 ∈ (0.005, 0.01).In Fig. 16 we have collected the values of E, Γ, and SE and plotted them against each other.Notice in Fig. 16(a) that it clearly indicates an inversely proportional trend for E and Γ.In Fig. 16(b), we notice that with increase an in E, the sample entropy shows a fairly increasing trend, whereas in Fig. 16(c), with an increase in Γ, the sample entropy decreases as expected.
4.4.6 Parameter space defined by (µ 0 , k) Fig. 17 is a collection of two-dimensional color plots in the parameter space defined by (µ 0 , k) with σ 0 = 0, P µ = 1, P σ ∼ 0.666, D µ = D σ = 0.005.It can be seen that for µ 0 > 0, Γ lies in the range (−0.38, 0.75) (Fig. 17(a)) and the nodes behave in asynchronous fashion as evident from a larger value in the averaged synchronization error (Fig. 17(b)).The value of E is maximum in the region µ 0 > 0, k > 1.5, where the sample entropy (Fig. 17(d)) reaches very high values too.When µ < 0, mostly all the nodes seem to be synchronous (E ∼ 0) and lie in the coherent region (Γ ∼ 1) for −1 < k < 0 and 4 > k > 2 as can be interpreted from the Γ and E color plots.Looking into the SE plot again, for µ 0 < 0 and −1 < k < 0, there exists a region where the sample entropy is very close to 0. For µ > 0, there is a high probability of occurrence of solitary nodes (Fig. 17(c)) spanning over k with 0.2 < Ns N < 0.8.When µ 0 < 0, solitary nodes appear with 0.2 < Ns N < 0.8 within 0.5 < k < 1.5.We also notice an upsurge in the value of the sample entropy at this parameter region.In Fig. 18 we have collected the values of E, Γ, and SE and plotted them against each other.Notice in Fig. 18(a) that it  clearly indicates an inversely proportional trend for E and Γ.In Fig. 18(b), we notice that with increase in E, the sample entropy shows a fairly increasing trend, whereas in Fig. 18(c), with increase in Γ, the sample entropy decreases as expected.

Bifurcation diagrams for synchronization
Next, we separately plot the bifurcation diagrams for the last instances of x m against the parameters σ 0 , µ 0 , P σ , P µ , and k.See Fig. 19.The first plot Fig. 19(a) depicts the bifurcation plot against σ 0 .The parameters have been set to P µ = 1, P σ ∼ 0.666, k = −1, µ 0 = −0.001,D µ = 0.005 and D σ = 0.005.Keeping aside the values of k and σ 0 , note that we have utilized the same parameter combinations in Fig. 15.Going back to the bifurcation diagram Fig. 19(a) we observe several cavities of synchronous states and asynchronous states almost evenly distributed in the range −0.01 < σ 0 < 0.01.For example see the fourth plot in Fig. 4(a) where we have used the same parameter combinations.Although there exists two clusters (one coherent, one solitary) and a single green node with Γ 2,m ∈ [−0.15, 0.75], they are very close to each other, as evident from the x−scale in the Y −axis and obviously from the synchronization error value E = 0.073, which is quite low.Similar kind of patterns have been spotted by the authors in the references [11].identical arguments can be made for the bifurcation diagram where we vary µ 0 in Fig. 19(b).Parameter values are P µ = 1, P σ ∼ 0.666, k = −1, σ 0 = 0, D µ = 0.005 and D σ = 0.005.In the bifurcation plot shown in Fig. 19(b), we observe that for an approximate range of µ 0 from −0.00045 to −0.0002, there is a gap where the nodes are synchronized.In the plot Fig. 19(e), we can thus state that the dynamics mostly exhibits asynchronity in the approximate range of k from 0.01 to around 3. Outside this range, the dynamics portray synchronity.Similar inferences

Conclusion
In this article, we have made a substantial improvement in the ring-star network of Chialvo neurons under the influence of an electromagnetic flux by considering a heterogeneous topology of the network.The motivation behind this was to study the spatiotemporal dynamics of an ensemble of neurons which mimics a realistic nervous system.Heterogeneity is realized not only in space but also in time.We have introduced a noise modulation to incorporate heterogeneity in space and a time-varying structure of the links between neurons that update probabilistically.Noise sources are sampled from a uniform distribution.How the network dynamics change on induction of other complicated distributions opens a future scope of study.Exploring various computer-generated diagrams like phase portraits, spatiotemporal plots and recurrent plots, we observe rich dynamical behaviors like two-cluster states, solitary nodes, chimera states, traveling waves, and coherent states.
One of the purposes of the paper was to study the appearance of this special kind of asynchronous behavior called solitary nodes and characterize them using a metric called cross-correlation coefficient.It was observed that the cross-correlation coefficient indeed characterizes the solitary nodes in an efficient manner.Two more measures that we implemented here were the synchronization error to study how synchronized the nodes behave, and sample entropy to study the complexity of the network dynamics.Under different pairwise combinations of the model parameters we have studied the response of each of these three measures and have tried to conclude how efficiently they relate.Asynchronization and  sample entropy, for example, portray a fairly convincing proportional relationship, whereas asynchronization and cross-correlation coefficient exhibit an inverse relation.These open a research direction to study how exactly they relate and if it is possible to establish an analytical relationship between them.Note that the trend obtained was a bit noisy.We hypothesize that this may be due to the heterogeneity introduced by noise modulations and time-varying links in the system.It would be interesting to see a noise-free relationship in a noise-free system which we plan to address in the future.Also, it would be very interesting to analytically and numerically study the existence of chaos in the ring-star network dynamics (by studying the Lyapunov exponent of the network system and methods to compute it).
Finally, we have also studied the one-parameter bifurcation diagrams of the synchronization of the nodes by plotting the final state of all the nodes against one of the parameters of the system.Interestingly, We observe some windows in the parameter regimes where the probability of all the nodes to be in a synchronized state is much higher than being asynchronous.
The fact that a ring-star network of chialvo neurons exhibits quite a rich dynamics, this study raises the question of how identical or contrasting the reported dynamical behaviors would be if we had considered different topologies and/or perturbations.Keeping this in mind, we plan to study the behaviors of an ensemble of Chialvo neurons under different topological arrangements: heterogeneous but quenched (no change in topology with time); has coupling strengths that are dependent on a normalized distance between a pair of neurons; is a multiplex heterogeneous network (See the recent work [55]); is perturbed by either temperature or photosensitivity.The fact that solitary nodes oscillate in a completely different phase from the main coherent ensemble, raises a future direction of exploring anti-phase synchronization with different coupling forms such as attractive and/or nonlinear repulsive  couplings.A recent similar study [54] with Van der Pol oscillators has also been published.

Figure 1 :
Figure 1: Heterogeneous ring-star network of Chialvo neurons.The nodes are numbered from 1, . . ., N .The star and ring coupling strengths are denoted by µ m and σ m for each node m = 1, . . ., N respectively.Different colors in the ring-star topology signify a range of heterogeneous values of µ m and σ m .
2,m ∼ 1, they lie in the coherent domain.Due to the noise modulations, there might be node behaviors that are realized by their Γ 2,m values lying in the range [−0.15, 0.75].Finally, the solitary regime is characterized by the domain where the nodes have Γ 2,m ∈ [−0.38, −0.15].

Figure 3 :
Figure3: Time series plot for the corresponding neuron.The sample entropy is 0.041, which is quite low.This indicates less disorderedness in the dynamics as can also be seen from not so irregular behaviour in the time series.
(a) Coherent and solitary nodes in two clusters (b) Mostly synced with two solitary nodes (c) Two clusters
Fig.9is a collection of two-dimensional color plots in the parameter space defined by (σ 0 , D σ ) with µ 0 = −0.001,P µ = 1, P σ ∼ 0.666, D µ = 0.005 and k = −1.Note that for all values of D σ , when −0.01 < σ 0 < −0.005, we see that the nodes prefer to remain distributed throughout the space, with Γ ∼ (0, 0.5) (Fig.9(a)) and behave in an asynchronous manner as evident from the violet region in the E color plot (Fig.9(b)) denoting higher asynchronity.As soon as σ 0 > −0.005, the nodes mostly tend towards synchronization in the coherent region as can be seen from the deep yellow and black regimes respectively in the Γ and E color plots respectively, although there might be instances when the nodes do not settle to synchronity in the coherent domain.Interestingly, there exists a very few to no solitary nodes as can be seen from the Ns N color plot (Fig.9(c)).Now looking into the SE plot (Fig.9(d)), note that again there is a fair proportional relationship between sample entropy and asynchronicity, with the maximum complex behaviour appearing at σ 0 ∼ 0. It can be seen that with increase in the mean coupling strength σ 0 , from −0.01 to 0.01, sample entropy starts increasing, reaches a peak at around σ 0 ∼ 0 and then the value drops.In Fig.10we have collected the values of E, Γ, and SE and plotted them against each other.Notice in Fig.10(a) that it clearly

Figure 7 :
Figure 7: Collection of two-dimensional color plots of (a) Γ, (b) E, (c) Ns N , and (d) SE in the parameter space defined by the mean coupling strengths (σ 0 , µ 0 ).The parameter space is a 40 × 40 grid with k = −1, P µ = 1, P σ ∼ 0.666, D µ = D σ = 0.005.An almost definitive bifurcation boundary within which most of the nodes lie in the coherent regime, is synchronous and has low sample entropy, is observed.Solitary nodes appear when the mean coupling strengths are σ 0 ∼ 0 and µ 0 > 0, corresponding the region where the nodes are asynchronous and have high sample entropy.

Figure 8 :
Figure 8: Comparison plots for the various measures: (a) E vs. Γ, (b) SE vs. E, and (c) SE vs. Γ, collected from Fig. 7. Figures (a) and (c) show an inverse trend whereas figure (b) shows a proportional trend.

4. 4 . 5
Fig 15 is a collection of two-dimensional color plots in the parameter space defined by (σ 0 , k) with µ 0 = −0.001,P µ = 1, P σ ∼ 0.666, D µ = D σ = 0.005.We observe that the Γ colored plot (Fig 15(a)) is mostly dominated by values that lie in (−0.38, 0.75), for which the nodes behave in a fairly asynchronous manner as depicted from the corresponding synchronization error plot E (Fig 15(b)).There exist two moderately distinct regions in the Γ plot where the nodes all lie in the coherent region and they are almost completely synchronous (see the deep yellow regions and the deep black regions in the Γ and E plots respectively).Comparing the synchronization error plot and the sample entropy plot, we observe again a fairly increasing

Figure 13 :Figure 14 :
Figure 13: Collection of two-dimensional color plots of (a) Γ, (b) E, (c) Ns N , and (d) SE in the parameter space defined by the coupling probabilities (P σ , P µ ).The parameter space is a 40 × 40 grid with σ 0 = 0, µ 0 = −0.001,D σ = D µ = 0.005 and k = −1.For values close to P µ = 0, for all P σ , there exist regions where Γ < 0.75, and the nodes behave asynchronously, having very high sample entropy value and are solitary.As soon as P µ > 0.5, the region is randomly distributed between both types of extreme behaviors.

Figure 15 :Figure 16 :
Figure 15: Collection of two-dimensional color plots of (a) Γ, (b) E, (c) Ns N , and (d) SE in the parameter space defined by the mean coupling strength and the electromagnetic flux coupling (σ 0 , k).The parameter space is a 40 × 40 grid with µ 0 = −0.001,P µ = 1, P σ ∼ 0.666, D µ = D σ = 0.005.The nodes mostly behave in a fairly asynchronous manner having Γvalues that lie in (−0.38, 0.75).There exist two moderately distinct regions in the where the nodes all lie in the coherent region and they are almost completely synchronous.The sample entropy is the highest at σ 0 ∼ 0 and k ∈ (0, 1.5).there exist some regions where all the nodes are solitary within k ∈ (−0.5, 0.5), σ 0 ∈ (0.005, 0.01)

Figure 18 :
Figure 18: Comparison plots for the various measures: (a) E vs. Γ, (b) SE vs. E, and (c) SE vs. Γ, collected from Fig. 17.Figures (a) and (c) show an inverse trend whereas figure (b) shows a proportional trend.