The article information
 Huang Shoufang, Zhang Jiqian, Wang Maosheng
 黄守芳, 张季谦, 汪茂胜
 Response Ability to External Signal Enhanced by Biological Spatial Con guration in Coupled HindmarshRose Neural System
 生物构型对耦合HindmarshRose神经元系统响应能力的影响
 Chinese Journal of Chemical Physics, 2016, 29(2): 265270
 化学物理学报, 2016, 29(2): 265270
 http://dx.doi.org/10.1063/16740068/29/cjcp1505088

Article history
 Dated: Received on May 4, 2015
 Accepted on September 29, 2015
Recently,more and more attention has been focused on the dynamic processes in complex networks because of their importance in explicitly mimicking highly complex structure of many realistic systems,especially in biological neural networks [1, 2, 3, 4, 5, 6, 7]. Some researchers reported that the dynamical behavior of the neuron is quite complex,which exhibits subthreshold periodic,chaotic responses,chaotic firing of spikes,as well as modelocked firing. These complex behaviors can be controlled by adjusting different configuration parameters under a periodic stimulation [8]. Furthermore,people found the effects of external period force and noise on the controlled system,and the complexity of system can be suppressed by the slaving principle [9]. These results imply the stimuli from the internal or external environment of neuron system,may play an important role in controlling its collective activities. In recent years,some researchers found the system size and coupling strength also play a vital role in FitzHughNagumo neural network [10]. By controlling bifurcation parameters in the networks of HodgkinHuxley or HindmarshRose (HR) neurons,the transition behavior from spiral wave to other state can be induced [11, 12].
Nowadays,the frontier of the research interest has focused on two kinds of collective behaviors in coupled complex systems. On one hand,the influences of the network topology on dynamical behaviors are investigated,especially the nervous system's dynamic features [13, 14, 15]. For example,the spiking oscillations show the best regularity at optimal network topology randomness [15]. The introduction of smallworld connections can change the secondorder phase transition into a firstorder one in the 2D Ising model [16]. On the other hand,how was collective dynamics of a system controlled by adjusting the parameters or introducing external stimuli. Some fascinating and important research showed the effect of system size on the dynamics behavior of system [17, 18, 19]. Time delays in a network of the HR model can induce phaseflip transitions leading to synchrony or out of synchrony [20, 21].
In the models mentioned above,the biological realistic cellular topology of the neural system was not considered. However,in actual tissue,even in small areas,its topological structure is significantly complicated [22, 23, 24]. It was found that the topology of intercellular connections in tissue has an essential influence on calcium ions signaling [25]. The geometrical factors may play a crucial role in both intra and intercellular calcium ions signaling [26]. An interesting problem is whether the spatial topology is closely related with the periodic bursting behaviors in the response dynamics of the neural system or not. And how does this topology control the efficiency of signal processing in neuronal system? Thus,in this work,by using the HR model as the units,and constructing a 2D regular lattice with the nearest connections,we mainly investigate the effect of spatial topology on the signal propagation behavior of HR coupled system.
Our results really show that in a proper coupling strength region,we can observe the statetostate transition (SST) behaviors between different states,the stimulus intensity and the delay time needed by SST behavior are different in all possible configurations.
Ⅱ. MODELThe paradigmatic model used here is first proposed by Hindmarsh and Rose in 1985 [27]. It is originally introduced to describe the firing behavior of neurons,and qualitatively explain the bursting type with inter spike intervals of real neurons. In this work,we consider the dynamics of a system of $N$ ($N$=4) coupled HR neurons,in which local coupling may correspond to the diffusion process between the nearestneighbours,and the longrange connections to the electrical synaptic connections among nonneighbouring cells. The dynamics of an isolated HR neuron are described by threevariable differential equations as the following (for a single cell $i$=1 and $\gamma$=0.0).
$ \begin{array}{l} {{\dot x}_i} = {y_i}  a{x_i}^3 + b{x_i}^2  {z_i} + {I_{{\rm{ext}}}} + \\ \;\;\;\;\;\frac{\gamma }{N}\sum\limits_{j = 1}^N {{C_{ij}}({x_j}  {x_i})} \end{array} $  (1) 
$ {\dot y_i} = c  d{x_i}^2  {y_i} $  (2) 
$ {\dot z_i} = r\left[{s\left( {{x_i}  {X_0}} \right)  {z_i}} \right] $  (3) 
where $x_i$ is the membrane potential,$y_i$ is a recovery variable,and $z_i$ is a slow adaptation current for the $i$th neuron. The constants $a$,$b$,$c$,$d$,$r$,$s$,and $X_0$ are set as $a$=1.0,$b$=3.0,$c$=1.0,$d$=5.0,$s$=4.0,$r$=0.006,and $X_0$=1.56,respectively. The last term in Eq.(1) denotes the coupling term for the complex neurons. $\gamma$ is the coupling strength,and the indices $i$,$j$ denote the cell number,$C_{ij}$ is the coupling parameter between the two neurons $i$ and $j$. If these two neurons are coupled to each other,$C_{ij}$=1,otherwise,$C_{ij}$=0. In this work,we consider the case of four coupled HR neurons,in which there are ten possible different configurations AJ [23] (see Fig. 1),and all possible plane configurations can be distinguished with different combination of $C_{ij}$.
It is well known,the bifurcation characteristic and the complex vibration modes of HR neurons are closely related to the variable $I_{\mathrm{ext}}$,so it is used as one of the control parameters to generate various firing patterns of neural network. The time series of the variable $x(t)$ is a spike train which can display the most important coded information in neural systems.
The bifurcation property of a single neuron can be obtained by calculating the inter spike intervals of neuronal firing,in which a spike pulse in the $x(t)$ variable is defined when the $i$th cell is excited so that $x(t)$ exceeds a given threshold value $x_0$ (note that a spike occurs when the state variable $x(t)$ is near $x$=1.2,so the excited threshold value $x_0$=1.0). One can notice that from the bifurcation diagram,the complex firing activity of neuron provides a rather rich dynamical behaviors including periodic bursting state and chaotic motions [28].
In the following sections,it is reasonable and feasible to study the effects of topology connectivity on the dynamical behaviors by using the HR model. Numerical integrations of Eq.(1) are carried out using the fourthorder RungeKutta scheme with time step of d$t$=0.001. To obtain each numerical result,those beginning data of simulation are discarded and 10$^6$ time steps are used.
Ⅲ. RESULTS AND DISCUSSIONTo investigate the effects of topological configurations on the dynamics of neural firing patterns in the 4neuron coupling system mentioned above,and show the SST behavior,under the condition of proper coupling strength,we adjust the system parameter as the following scheme: change topological configuration for different critical current $I_{\mathrm{ext}}$,which locals in the righthand side for each bifurcation point [28]. To quantify the SST effects induced by external stimulus,the time series of output signal was analyzed and the values of the inter spike intervals versus the external current $I_{\mathrm{ext}}$ are obtained.
A. Response ability to external weak signal enhanced by topologyFirstly,we consider the case that only the cell 1 is stimulated by external current,and discuss which type of configuration is more advantageous to the signal response in coupled system. To do this,we fix the coupling strength $\gamma$=0.07,and select one type of configurations (such as type H),then set the external current as $I_{\mathrm{ext}}$=1.20,the neural system may keep in the 1period state. One can find,with the increasing of $I_{\mathrm{ext}}$,the external stimulus information could be propagated along this coupled chain,and all cells could turn into 2period state successively,indicating the phenomenon of SST occurs under this condition. Similar synchronous SST behavior could be seen for other configurations. Of cause,the critical values $I_{\mathrm{ext}}$ of the four cells for SST are different under the condition of certain coupling strength,so we select the critical current value of the cell which induces the system to enter the transition state as the critical standard.
To describe the response ability of the system to the external feeble signal enhanced by proper configuration,we introduce the synchronous error $e$ of critical value by the following formula:
$ \bar I = \frac{1}{N}\sum\limits_{i = 1}^N {{I_{i,c}}} ,i = 1,2,\cdots ,N $  (4) 
$ e = \frac{1}{N}\sqrt {\sum\limits_{i = 1}^N {(I_{i,\mathrm{c}}  \bar I)^2 } } $  (5) 
where $I_{i,\mathrm{c}}$ is critical values of $I_{\mathrm{ext}}$,$e$ can provide an indication of firing patterns transition dynamics in neural system. The smaller the current value of $\bar I$ is,the more easily the state transition occurs,also,the smaller $e$ implies the better effect of the synchronous transition between two states.
The critical current values $I_{i,\mathrm{c}}$ and synchronous error $e$ as functions of the different coupling strength for different configuration types are depicted in Fig. 2. In Fig. 2 (a) and (b) we could find with the increment of $\gamma$,the critical values of $\bar I$ and synchronous error $e$ reduce continuously for a certain type of configuration. Of cause,the values of types A and B are larger than others indicating the coupled system could not transit to the 2period state easily in these configurations. This means the 1D chain configuration is not suitable for the transition behavior. On the contrary,if cells are coupled to C,D,and J configurations,the system could jump to the 2periodic state very easily than those in A and B (Fig. 2(c)). Above results imply that the HR system may enhance its response ability to external stimulus through the mode of SST,by using proper spatial topology configuration and adjusting itself parameter.
B. Response time of SST behavior shorten by topologyTo study the effects of spatial configuration on the SST behavior in coupled noisy system,so we introduce an external noise in Eq.(1),i.e.:
$ \begin{array}{l} {{\dot x}_i} = {y_i}  a{x_i}^3 + b{x_i}^2  {z_i} + {I_{{\rm{ext}}}} + \\ \;\;\;\;\;\;\xi (t) + \frac{\gamma }{N}\sum\limits_{j = 1}^N {{C_{ij}}({x_j}  {x_i})} \end{array} $  (6) 
where $\xi (t)$ is Gaussian white noise with $\langle\xi (t)\rangle$=0,$\langle\xi (t)\xi (t')\rangle$=$D^2\delta(tt')$,$D$ is the noise intensity. To discuss the transition between 1period and 2period states,we select the typical configurations B and D as the examples,and fix $\gamma$=0.07,$I_{\mathrm{ext}}$=1.20. Here,the external noise intensity is an adjustable parameter,other parameters are the same as the above,i.e. $a$=1.0,$b$=3.0,$c$=1.0,$d$=5.0,$s$=4.0,$r$=0.006,$X_0$=1.56. One can find that,when noise is introduced,and its intensity is set to be $D$=0.1,the system may transit to the 2period state. We also found the SST happen quickly in cells 2 and 3 when the noise is switched on the cell 1,while cell 4 reach the same state only after a period of time,indicating there is a certain time delay in the propagation of the stimulus along the coupling system,and this time delay depends on noise intensity. Similar phenomenon can also be found in other 9 configurations. We can define this time delay $\tau_{1,j}$ as the time interval among which the SST first appears in first cell and then transmits to the cell $j$ (see Fig. 3(a)). The mean value of time delay $\langle\tau\rangle$ among cells in coupled system reads by the following:
$ \tau _{1,j} = t_j  t_1 , { j = 2,3,4 } $  (7) 
$ \langle\tau\rangle = \frac{1}{{N  1}}\sum\limits_{j = 2}^N {\tau _{1,j } } $  (8) 
The mean delay time $\langle\tau\rangle$ as a function of noise intensity for different configurations are depicted in Fig. 3(b),from this figure,one can find that,for each configuration,the delay time decreases with increment of noise intensity,which means that the signal transmission processing could be accelerated by adjusting the noise intensity.
We also find that the mean value of delay times $\langle\tau\rangle$ is different for ten configurations with the same noise intensity. For example,the mean delay time in configuration A is larger than that in the type D,and especially,with increment of noise intensity the mean delay time $\langle\tau\rangle$ in the type D decreases very faster than that in other types. These results may indicate that in configuration D,the system is more advantageous to improve the response speed and response ability.
After further study,we found that,the delay time not only depends on noise intensity,but also depends on the connection distance (CD,between the stimulated cell and a given cell). For example,in configuration A,the CD between cell 2 and cell 1 is 1,while between cell 3 and cell 1 is 2. The CD between the cell 1 and the another cells in all configurations are listed in Table Ⅰ. These data may provide us some useful clues to analyze the SST behavior. Obviously,in configuration D and G,the average distance is the shortest,while in configuration A the CD is longest,this difference may be the reason for that the D and G configurations are more suitable for the transmission of information than others.
C.The topology influence on selection effect of systemIn this work,we have also discussed the selection effect of the topology to the external signal,the corresponding dynamic equations are as follows:
$ \begin{array}{l} {{\dot x}_i} = {y_i}  a{x_i}^3 + b{x_i}^2  {z_i} + {I_{{\rm{ext}}}} + \\ \;\;\;\;\;A\sin (2\pi ft) + \frac{\gamma }{N}\sum\limits_{j = 1}^N {{C_{ij}}({x_j}  {x_i})} \end{array} $  (9) 
where $A$ is the critical value of amplitude. According to the same scheme used for adjusting noise mentioned above,we first set the system to be in the 1period state,then select the amplitude or frequency as control parameters,respectively. The distribution of inter spike intervals versus the amplitude A for a constant frequency ($f$=0.001) is depicted in Fig. 4. It is found that,in configuration H,cell 1 is in 1period state if $A$<0.003,while if $A$ increases to 0.003,2period state appears,indicating the SST occurs. But in the other three cells,the critical value of $A$ is different when the SST occurs: in cells 2,cell 3,and cell 4,$A$=0.152,0.154,and 0.14,respectively. Thus,to make the whole system appear SST behavior,the amplitude should be set to the maximum critical value,i.e.,$A$=0.154.
Similar results are also observed when we adjust the frequency,and the critical values of amplitude (or frequency) of each configuration for different frequency (or amplitude) of the external signal are depicted in Fig. 4.
In Fig. 4(a),it is observed that in some proper configurations (such as in C,D,and J),with the increment of $f$,the periodic motion in the cell 1 could transmit easily from 1period to 2period,indicating that SST occurs. These results suggest that some proper types of configuration (such as C and J) can enhance effectively the sensitivity of the HR system to the weak signal,implying the selection effect of topology.
We now consider some interesting configurations of four cells in detail. In configuration E and I,the disposition of cell 1,cell 3 and cell 4 is symmetrical relative to cell 2 (Fig. 1 E and I),and the average connection distance of coupled cell system is the same (see Table Ⅰ),but the system shows different selection effect on the external signal and this difference is also manifested in weaker response in cell 3 and cell 4 in configuration E,as compared with configuration I (see Fig. 4(a) and (b)). The topological difference between configurations E and I consists in the presence (E) and absence (I) of a connection between cell 3 and cell 4. This is connected with the fact that cell 3 in configuration E exchanges external signal with cell 4 in addition to cell 2,whereas there is no such exchange in configuration I. That additional exchange decreases the value of $A$ and $f$ in cell 3 and cell 4 when the SST occurs in configuration E.
This interesting phenomenon also takes place in other configurations. Such as in configuration G and H,the disposition of cell 2,cell 3,and cell 4 is symmetrical relative to cell 1 (Fig. 1 G and H),and the average connection distance of coupled cell system is the same (see Table 1),but there is an additional connection between cell 3 and cell 4 in configuration G. Obviously,this additional connection makes the system show different selection effect and decreases the time delay among four cells (see Fig. 3). In other words,additional connection can improve the response speed and ability of coupled cells system.
Ⅳ. CONCLUSIONWe have investigated the effects of topological structures of system on the response ability to external stimulus signal in coupled HR neural system. It is found,when the controlling parameters such as the external current intensity or noise intensity are adjusted to be the critical values,and when the system is coupled to some proper configurations,SST phenomenon could appear in coupled system. This result may suggest that the neural system could improve its response ability and sensitivity to external stimulus signal by choosing suitable configuration.
We hope our results could provide some interesting information for understanding the properties of collective response behavior and mechanisms of information processing in coupled real neurobiological systems. Note that we only consider two dimensional configurations of four cells in this work. It would be interesting to consider the role of topological structure in more cells coupled system and three dimensional tissues,which is a subject of our future research.
Ⅴ. ACKNOWLEDGMENTSThis work was supported by the National Natural Science Foundation and Special Found for the Theoretical Physics of China (No.21103002 and No.11047017),the Special Foundation of Education of Anhui Province for Excellent Young Scientists (No.2011SQRL023),and the Natural Science Funds of Anhui Province of China (No.1508085MA15).
[1]  M. Bucolo, L. Fortuna, and M. LaRosa, Chaos, Solitons Fractals 14, 1059 (2002). 
[2]  F. Han, Q. S. Lu, M. Wiercigroch, and Q. B. Ji, Inter. J. NonLinear Mechanics 44, 298 (2009). 
[3]  H. S. Chen, J. Q. Zhang, and J. Q. Liu, Phys. A 387, 1071 (2008). 
[4]  M. S. Wang, Z. H. Hou, and H. W. Xin, Phys. Lett. A 334, 93 (2008). 
[5]  P. Wang, J. Q. Zhang, and H. L. Ren, J. Chem. Phys. 23, 23 (2010). 
[6]  L. X. Duan, Q. S. Lu, and Q. Y. Wang, Neurocomputing 72, 341 (2008). 
[7]  X. J. Sun, M. Perc, Q. S. Lu, and J. Kurths, Chaos 18, 023102 (2008). 
[8]  E. M. Izhikevich, Int. J. Bifurcation Chaos. 10, 1171 (2000). 
[9]  W. Wang, Y. Wang, and Z. D. Wang, Phys. Rev. E 57, R2527 (1998). 
[10]  D. Q. Wei, X. S. Luo, and Y. L. Zou, Commun. Theor. Phys. 50, 267 (2008). 
[11]  J. Ma, C. N. Wang, W. Y. Jin, and Y. Wu, Appl. Math. Comput. 217, 3844 (2010). 
[12]  C. N. Wang, J. Ma, J. Tang, and Y. L. Li, Commun. Theor. Phys. 53, 382 (2010). 
[13]  O. Kwon and H. T. Moon, Phys. Lett. A 298, 319 (2002). 
[14]  M. Lin and T. L. Chen, Phys. Rev. E 71, 016133 (2005). 
[15]  Q. S. Li and Y. Gao, Phys. Rev. E 77, 036117 (2008). 
[16]  A. D. Sánchez, J. M. López, and M. A. Rodriguez, Phys. Rev. Lett. 88, 048701 (2002). 
[17]  A. Pikovsky, A. Zaikin, and M. A. de la Casa, Phys. Rev. Lett. 88, 050601 (2002). 
[18]  Z. H. Hou, J. Q. Zhang, and H. W. Xin, Phys. Rev. E 74, 031901 (2006). 
[19]  J. Q. Zhang, C. D.Wang, M. S.Wang, and S. F. Huang, Neurocomputing 74, 2961 (2011). 
[20]  M. B. Adhikari, A. Prasad, and M. Dhamala, Chaos 21, 023116 (2011). 
[21]  M. M. Shrii, D. V. Senthilkumar, and J. Kurths, Europhys. Lett. 98, 10003 (2012) 
[22]  J. H. Evans and M. J. Sanderson, Am. J. Physiol. 277, L30 (1999). 
[23]  H. Niessen, H. Harz, P. Bedner, K. Kramer, and K. Willecke, J. Cell Sci. 113, 1365 (2000). 
[24]  M. V. L. Bennett and V. K. Verselis, Semin. Cell. Biol. 3, 29 (1992). 
[25]  I. V. Dokukina, M. E. Gracheva, E. A. Grachev, and J. D. Gunton, Phys. D 237, 745 (2008). 
[26]  K. TsanevaAtanasova, D. I. Yule, and J. Sneyd, Biophys. J. 88, 1535 (2005). 
[27]  J. L. Hindmarsh and R. M. Rose, Proc. R. Soc. London Ser. B 221, 87 (1984). 
[28]  S. F. Huang, J. Q. Zhang, and S. J. Ding, Chin. Phys. Lett. 26, 050502 (2009). 