ARTICLE IN PRESS
Physica A 353 (2005) 258–270 www.elsevier.com/locate/physa
Computer simulation of a central pattern generator via Kuramoto model Frederico Alan O. Cruz, Ce´lia Martins Cortez Dept de Cieˆncias Fisiolo´gicas, do Estado do Rio de Janeiro (UERJ), Av. Prof. Manuel de Abreu, 444-51 andar, 20555-170, Rio de Janeiro, RJ, Brazil Received 20 June 2004; received in revised form 11 October 2004 Available online 25 January 2005
Abstract The presence of oscillatory activity in motor command signals was demonstrated, and the possible existence of synchronous oscillatory activity within the sensorimotor system has been suggested. We adopted a set of 1000 coupled oscillators distributed on 10 10 square lattice as a model of a central pattern generator. Each lattice element contained 10 closer neighbor oscillators. We introduced a new method to estimate the equilibrium order parameter to characterize the degree of synchronization of the oscillator system. The numerical scheme was based on a model that is simple, but that contains the essential features of a broad class of cortical neurons. The adopted boundary conditions were based on known characteristics of the synaptic transmission. By means of numerical simulations of the Kuramoto, we simulated the dynamical behavior of a central pattern generator. Our results are in agreement with some experimental data from literature. So, we believe such a method can be used to model dynamic behavior of a central pattern generator with good computational efficiency. r 2005 Elsevier B.V. All rights reserved. Keywords: Computer simulation; Coupling oscillators; Order parameter; Motor control
Corresponding author. Tel.: +55 21 25876407; fax: +55 21 25876129.
E-mail address:
[email protected] (C.M. Cortez). 0378-4371/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2004.12.046
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
259
1. Introduction Biological rhythms have been targets of great interest of many science areas [1–3]. Neurophysiological studies have shown that there is a pattern of local synaptic circuitry in many parts of the brain [4–7]. This pattern consists of excitatory and inhibitory neurons connected to each other so that action potentials (APs) generated earlier excite the latter, which in sequence reciprocally inhibit the earlier, constituting one of the basic mechanisms for the generation of periodic activity in the brain. Such a circuit reproduces an oscillator system, and it has been named neural oscillator [8]. The neural oscillators within one brain structure can be connected into a network because the neurons can have synapse with other distant neurons. The networks have different synaptic organizations, and consequently present diverse dynamical properties. The notion of synaptic organization is intimately related to the anatomicphysiological view of the brain [9,10]. The presence of oscillatory activity has been demonstrated in motor control signals. It has incited interest in probable functions of synchronous oscillatory activity within the sensorimotor system [4,11,12]. Indeed, experimental studies have demonstrated that all rhythmic movements are related to central oscillators functioning as central pattern generators (CPNs) [13,14]. Descending inputs from these generators initiate stereotyped movements, but play relatively slight roles after the movement begins. The wide view of motor control includes a set of central nervous system (CNS) oscillators, local or distributed, controlling individual limbs or body segments. These local oscillators might, if disconnected, display conflicting cycles periods and therefore no stable phase relationships. In normal conditions, these oscillators are coordinated by means of the interplay between neuronal interactions within CNS and peripheral sensory inputs [15]. Sensory inputs can go through the CPG and change oscillatory periods. It could promote mechanical coupling between body segments, or even form a sensory-central oscillatory loop [16]. In past years, physicists have shown great interest in the behavior of systems composed of a large number of coupled oscillators. Some of the most interesting applications include analysis of nonlinear oscillators in neural circuits [10,17–19]. To study the dynamics of a system of oscillators a relatively weak coupling among them together with a random distribution of frequencies is usually assumed. When the coupling is small compared to the extension of natural frequencies, each oscillator runs at its natural frequency and the system behavior is incoherent. The incoherence remains until a certain threshold of coupling is reached, in which a small cluster of oscillators synchronizes. The degree of synchronization increases until all the oscillators become bound in phase and amplitude [20]. A favorable way of trading analytically with several coupled oscillators is through Kuramoto’s model [21], which is one of the simplest models to describe collective behavior of systems in biological sciences. The Kuramoto model considers a collective oscillatory behavior when the natural frequencies are taken from a given frequency distribution, as occurs in a real system. The existence of this behavior with phase and frequency entrain depends on the coupling strength and on the choice of the frequency distribution [22].
ARTICLE IN PRESS 260
F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
The purpose of this paper is to report results obtained from a computer simulation of the dynamical behavior of a CPG, and to introduce a new method to estimate the equilibrium order parameter to characterize the degree of synchronization of the oscillator system. The CPG has been modeled as a set of 1000 coupled oscillators distributed on 10 10 square lattice, where each lattice contains 10 oscillators, trading via numerical simulations of Kuramoto’s model. The adopted boundary conditions have been based on known characteristics of the synaptic transmission.
2. Elements of the sensorimotor circuit The brain, brain stem and spinal cord constitute the CNS. This is responsible for integrating sensory information and responding appropriately. The brain composes the supra-segmental nervous system (SSNS), while the spinal cord and brain stem form the segmental nervous system (SNS). There are two general kinds of tissue in the CNS: gray matter and white matter. In the gray matter masses of the nerve cell bodies, dendrites, each covered with synapses, and axons are found. Neurons in the gray matter organize either in layers, surface layer called the cortex, or in inner neuron clusters called nuclei. The white matter consists mostly of bundles of axons, and its white look is due to the sheath of myelin involving the axons [23]. One of the functions of the SNS is carrying out signals between the brain and the rest of the body. In segmental white matter there are several longitudinal tracts constituted by nerve fibers ascending and descending. The segmental gray matter is the integrative area for the reflexes and other motor functions. It is able to display musculoskeletal reflexes without input from the brain [24–26]. However, segmental reflex responses are modulated by hierarchical motor control that is attributed to serial activation of a motor command chain, because after synapse in segmental gray matter, the sensory fiber collaterals proceed to higher segmental levels and the suprasegmental area [27,28]. Sensory signals enter the CNS through the sensory nerve roots, reaching the supra-segmental area through ascending fibers. The pathway until the brain cortex involves many synapses. The larger the number of synapses in the multisynaptic pathways, the longer the delay between spinal input and output signals, because of the period of time required to complete all synaptic events. This is called synaptic delay, and is normally about 0.5 ms [23]. In the brain, primary somatosensory cortex (PSC) receives somatosensory input from the SNS. PSC contains neurons that register the sense of touch, and is located in the parietal lobe of the brain. Primary motor cortex (PMC) is located in the frontal lobe of the brain. It provides the bulk of the corticospinal tract, but other cortical areas contribute too. The corticospinal tract is an important descending motor pathway. This tract starts in PMC and ends in the gray matter of the spinal cord segment. The other descending tracts involved with motor control originate in several brain stem areas. Both PMC and PSC are highly organized and present specific regions relating to each part of the body [23].
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
261
Most of the computational power of the brain cortex is used for integrating somatosensory inputs and coordinating body movements, both consciously and unconsciously. So, it controls reflexes, muscle tone, posture, and voluntary movements.
3. Adopted model and method The circuit used to model a CPG for motor control is shown in Fig. 1. The signal generated in element R reaches element C via V1 : After data processing, C-outputs follow other elements within CNS. The circuit elements and their basic equations are described as follows. Element R: This element simulates a sensory receptor. R output is represented by Dirac delta function [29] and denotes the APs generated from stimulus transduction. The AP frequency is related to the stimulus intensity RðtÞ ¼
n X
dðt ti Þ ,
(1)
i¼1
where ti denotes the occurrence intents of the ith AP in the sensory element. The interval between the APs is a random Gauss process (truncated for tp0), where the mean is equal to the standard deviation.
SENSORIOMOTOR CORTEX C
O1
O2 ………
On
SNC
V1
R Fig. 1. Schematic representation of the neural circuit adopted as a model of a central pattern generator system to motor control (R ¼ sensory receptor; C ¼ population of neurons in sensorimotor cortex; V1 ¼ pathway between R and C elements; O1 :::On ¼ element C output).
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
262
Element C: It represents one population of excitatory and inhibitory neurons connected to each other within the sensorimotor cortex. The neurons have similar patterns of synaptic connections and dynamic behavior, and constitute one neural oscillator circuit. So element C simulates a CPG. It is able to decode and process signals from element R, as well as integrate all afferent inputs for coordinating movements. Element C consists of a set of 1000 coupled oscillators with long-range coupling, which were distributed on a 10 10 square lattice. Each lattice element contains 10 closer neighbor oscillators. The coupling modifies the frequencies, acting as an external oscillatory force. The degree of synchronization among the individual oscillators of C element defines the response of others system elements. By activity of element C we mean the number of spikes per unit time. For any system of N weakly coupled limit cycle oscillators, we can assume, according to the perturbation theory, that N F ðYj Þ F ðYi Þ dYi KX ¼ F ðYi Þ þ , N iaj raij dt
(2)
where Y is the phase and F ðYi Þ represents the coupling function, k is the coupling strength and rij is the distance between oscillators i and j. For a ¼ 0; the coupling is uniform, and a ! 1 corresponds to very near neighbor coupling. Eq. (2) was applied to the elements of a 10 10 two-dimensional lattice, and we adopted F ðYi Þ ¼ A sinðYÞ; with YðtÞ ¼ Y0 þ ot for p=2potpp=2: Studies using coupled oscillators require the estimation of a certain normalized parameter to characterize the degree of synchronization in the asymptotic regime. Such a parameter was named order parameter by Kuramoto [22]. This author proposed a complex order parameter for phase denoted by sf ¼
N 1X eifj , N j¼1
(3)
where N is the oscillator number within the system, s is the order parameter and fj is the phase of oscillator j: In our study, since we have been interested in the mean behavior of oscillators and cluster formation, we have taken a new method to estimate the equilibrium order parameter (s), which was based on statistical analysis of the phases of oscillators (Y). In this method, s is given by s ¼ expðgÞ ,
(4)
where g is coefficient of variation of the Y-values. So, in Eq. (4) we used the variation coefficient of Y instead of the summation on the phase exponentials used in Eq. (3). This alteration was based on the observation that, for a certain frequency distribution in our model, when phases of oscillators were very near, the order parameter estimated by Eq. (3) already indicated the formation of several clusters or that the system was almost totally isolate. This result could lead us to a false conclusion about phase distribution. The results from this analysis will be discussed in the next section.
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
263
In our method, to find the g-value, firstly it is necessary to evaluate the mean value and the variance of the Y-distribution within the oscillator system. The coefficient of variation of a set of values is defined by ratio variance/mean value of the distribution [30]. From this definition, we can write for oscillators that g¼
sY , ¯ Y
(5)
¯ is the arithmetic mean value of phases. It where sY is the variance of phases and Y follows from the definition of variance and mean value of any distribution [30] that we can also write qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P ¯ 2 N N i¼1 ðYi YÞ g¼ . (6) PN i¼1 Yi The coefficient of variation is usually written as percentage, which becomes null (g ¼ 0) for equal Y values. In general, the following correspondence is used [30]: go10; low variation of values; 20ogp10; intermediate variation of values; 30ogp20; high variation of values; 30pg; very disperse values. In the same way, we could classify our order parameter as follows: s ¼ 1; for coherent system or completely synchronized state; 0:9oso1; for low variation of Y-values or almost synchronized state; sp0:9; for incoherent oscillations or desynchronized state. We used 50 phase measures to calculate the phase mean value. Such large sample size allows us to observe variation of values. At t ¼ 0; we supposed that oscillations are incoherent leading to s ¼ 0: Our model exhibited the two characteristic types of qualitative synaptic behavior: spatial and temporal summation. It could provide an AP in neuronic response to a superthreshold stimulus, but a simple pulse from the sensory fiber was considered subthreshold. A neuron could not fire before its membrane potential had reached firing threshold, but a series of subthreshold excitatory input signals was able to depolarize the neuron membrane gradually. So, only when a sufficient depolarizing effect of inputs was provided, and the plasmatic membrane potential became greater than the threshold potential, the neurons of element C fired. It is known that the transmitter substance released by a single presynaptic terminal is able to generate an excitatory postsynaptic potential no more than 1 mV, but about 15–20 mV is required for the synaptic firing. For instance, when several excitatory presynaptic terminals are simultaneously stimulated, their effects can summate, even though these terminals are distributed over wide areas of the neuron. In this case, the intrasomal potential became more positive by about 1 mV for each added excitatory discharge. When the excitatory postsynaptic potential became great enough, the threshold for firing was reached, and an AP was generated on the neuron axon. Its synaptic characteristic is known as spatial summation. In the same way, successive postsynaptic potentials of a presynaptic terminal occurring rapidly enough can summate. It is called temporal summation, because when a terminal fires, the released transmitter substance opens the membrane channels for 1 ms or so. Since the postsynaptic potential lasts up to 15 ms, a second opening of the same channel
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
264
can increase the postsynaptic potential to a still greater level so that the more rapid the rate of terminal stimulation, the greater the effective postsynaptic potential [23]. A Fortran 77 (Visual Fortran 95 for Windows) based program was developed for the model described above. In this program, numerical values typical of parameters related to the neuron membrane and the interneuronal synapses found in the literature were applied. Then, we examined the ability of our model to reproduce the behavior of CPG for motor control.
4. Results and discussion Our model takes the population of neurons in sensorimotor cortex involved with the motor control as a set of 1000 diffusively coupled nonlinear oscillators distributed on a 10 10 lattice (element C), with each lattice element containing 10 very near neighbor oscillators. The oscillators are connected weakly via a common medium forced by an external input and with a single output. The external inputs model the depolarizing signals generated by muscular receptor activity. In our model, the muscular receptor was represented by element R. Physiological characteristics of digitalling and autopropagability of the AP were represented by Dirac delta function. External inputs were processed by element C, and C-output was dependent of the degree of synchronization of oscillators. The dependence of the coupling strength (k) on system size can be evaluated in Figs. 2(a) and (b). These figures present the variation of effective k, ke ¼ k=Nra (Eq. (2)), in function of distance between number oscillators within the lattice. Consistently with the idea that typical strength of the net local field on each oscillator falls with the distance between the oscillators within the lattice, Fig. 2(a)
α = 0.0 α = 0.3
α = 0.1 α = 0.4
α = 0.2 α = 0.5
α=0 α = 0.3
0.5
0.8
0.4
0.6
0.3
k*
ke
1.0
α = 0.2 α = 0.5
0.4
0.2
0.2
0.1
0.0
0.0 0
(a)
20 40 60 80 Distance between Oscillators (a.u.)
0
100
(b)
25
50 75 100 125 150 175 Number of Oscillators
200
Fig. 2. Effective coupling strength, ke ¼ k=Nra (Eq. (2)), in function of: (a) distance between the oscillators and (b) number oscillators within lattice, for several a-values.
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
265
shows that ke value decreases with the increase in rij ; for aa0: Fig. 2(b) shows the important fall of ke with the increase in the number of oscillators (N), with ke o0:1 for N410: These figures talk on effects from neighbor oscillators on a certain oscillator j: They show that k value required to synchronize the oscillator system is dependent on r and N. The higher the number of oscillators, or the distance between them, then harder the pertubation on the system. We verified the onset of synchronization of the oscillator system as a function of the number of processing time steps (t), by means of numerical simulations using Eqs. (2), (4) and (6). The dependence of s on k for fixed t and different a values can be seen in Fig. 3(a). In this we can see the direct proportionality between the coupling ratio and k, since s is a convenient measure of the extent of synchronization. We may also see that k requiring for s reaches saturation value (s ¼ 1) and decreases with a: It is possible to note in plots of s in function of k for several t values (Fig. 3(b)) the increase in curve inclination with the increase in the number of processing time steps. It reveals an inverse relationship between k and t: In fact, the system was able to selforganize depending on time, for fixed k. Fig. 4 shows the behavior of the oscillator system constituted by a 10 10 square lattice for three s-values, with a ¼ 0: In this figure we can see that the degree of synchronization increases with increase in s; and that the system becomes coherent for s ¼ 1: The uniform regions in graphics correspond to the formation of clusters or local synchronism. The results above show the existence of interrelationship among a; k, t; N, r and s: So, the synchronous oscillatory activity is not dependent exclusively on one or two parameters, but is dependent on all of them. The increase in N requires elevation of k or time for oscillators to synchronize. So, independent of N, the oscillator system can self-organize when its parameters are adjusted. For lattices smaller than 10 10; we observed large variation in values of parameters, indicating the possibility of occurrence of very discrepant values among them, although synchronization could occur.
0.8
0.8
0.6
0.6
0.4 0.2
α=0 α = 0.1 α = 0.3 α = 0.5
Time step 0 5 10 20 30
σ
1.0
σ
1.0
0.4 0.2
0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 (b) (a) k
0.1
0.2
k
0.3
0.4
0.5
Fig. 3. Order parameter in function of k for: (a) different a and fixed t; and (b) several t and a ¼ 0:
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
266
σ = 0.53 α=0
σ=0 α=0
10 9
(a)
10 9 8
7
6
5
4
3
2
(b)
1
8
7
6
5
4
3
2
1
σ=1 α=0
10 9
(c)
8
7
6
5
4
3
2
1
Fig. 4. Behavior of the oscillator system constituted by 10 10 square lattice for three s-values: (a) s ¼ 0; (b) s ¼ 0:25; (c) s ¼ 1; being a ¼ 0 and t ¼ 100: The uniform regions in graphics correspond to the formation of clusters or local synchronism.
Study of s relations using Eqs. (4) and (5) have shown that, despite lattice size, s always have be zero for k ¼ 0; for all used a values (Figs. 3(a) and (b)). In contrast, using Kuramoto’s equation (Eq. (3)) we have observed that initial k only tended to vanish for large N. Fig. 5 presents the temporal evolution of four oscillators of initial phases 4p/9, 5p/18, p/6 and p/2. One can note the low phase coincidence at the stage shown in this figure. In fact, according to Eqs. (4) and (5), the oscillators are far from a synchronized state at this stage, since estimation using these equations gave us gb30
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
267
Fig. 5. Temporal evolution of four oscillators of initial phases: 4p/9 (blue), 5p/18 (red), p/6 (green) and p/2 (orange). One can note the low phase coincidence at the stage shown in the figure, in contrast with the value of the order parameter (s 0:72) given by Eq. (3).
and s ! 0: However, the use of Eq. (3) gave us s 0:72; indicating that oscillators are almost completely synchronized. In Figs. 6(a)–(d) we can see the oscillatory behavior of two oscillators with no neighbor for different s values, according to Eqs. (2), (4) and (6). It is particularly interesting to observe the degree of synchronization for s ¼ 0:75 given by these equations. Based on these results we introduce Eq. (4), which includes a statistical analysis on the phases of oscillators, as a valid method to estimate order parameter. The input spectrum of element C during 100 ms is shown in the Fig. 7(a). Figs. 7(b) and (c) illustrate the firing pattern of element C during 200 ms for two s-values and t ¼ 100 (equivalent to 1 ms). It was the processing time in element C for each input signal, and synaptic delay was taken as 6 ms. One can observe in these last two figures the random firing behavior of the oscillator system with mean frequencies of about 125 Hz for s ¼ 0:53 (Fig. 7(b)) and 45 Hz for s ¼ 0:25 (Fig. 7(c)). This last frequency is in order of magnitude of that observed in experimental studies [4,5,9,12] which that demonstrated that neurons within the sensorimotor cortex of monkeys and humans have synchronous oscillatory activity around 15–44 Hz. This activity range is known to influence descending motor commands to contralateral muscles [5,12]. It is interesting to observe the sigmoidal feature of Figs. 3(a) and (b). It is known that activation functions used in neural networks generate sigmoidal feature curves [31,32]. So, we can say that our model gives a response similar to that used in neural networks of the kind multi-layer perceptron (MLP), where the adopt transference function is typically sigmoidal.
ARTICLE IN PRESS 268
F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
Fig. 6. Oscillatory behavior of two oscillators with no neighbor for different s values estimated by use of Eqs. (2), (4) and (6): (a) s ¼ 0; (b) s ¼ 0:25; (c) s ¼ 0:75; (d) s ¼ 1:0:
In summary, we have presented an efficient model for simulating central pattern generator and have introduced a new method to estimate the equilibrium order parameter to characterize the degree of synchronization of the oscillator system. We have focused on the ability of our model to provide numerical results in order of magnitude of that observed in experimental studies with sensorimotor neurons. The numerical scheme was based on a model that is simple and yet contains the essential features of a broad class of cortical neurons. The model equations were able to describe the behavior of each element of the adopted circuit. Thus, the model presented here offers much in terms of simplicity, speed, and wide applicability. The simulation environment can be easily expanded to include additional neuronal geometry; cell populations can be added to the model as well. So, we believe that such method can be used to model dynamic behavior of a central pattern generator with good computational efficiency.
ARTICLE IN PRESS F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
269
σ = 0.53 1
Spike
Spike
1
0
0 0
10
(a)
20
30
40 50 60 Time (ms)
70
80
0
90 100
25
50
(b)
75
100 125
150 175
200
Time (ms)
σ = 0.25
Spike
1
0 0
(c)
25
50
75
100 125
150 175
200
Time (ms)
Fig. 7. (a) Spectrum of element C input during 100 ms. (b) Firing pattern of element C output during 200 ms. Adopted time constant was 10 ms, k ¼ 0:035; t ¼ 100 and s ¼ 0:25:
Acknowledgements This work was supported by FAPERJ - Fundac- a˜o Carlos Chagas Filho de Amparo a` Pesquisa do Estado do Rio de Janeiro.
References [1] S.H. Jezzini, A.A. Hill, P. Kuzyk, R.L. Calabrese, Detailed model of intersegmental coordination in the timing network of the leech heartbeat central pattern generator, J. Neurophysiol. 91 (2004) 958. [2] M. Carro-Juarez, S.L. Cruz, G. Rodriguez-Manzo, Evidence for the involvement of a spinal pattern generator in the control of the genital motor pattern of ejaculation, Brain Res. 975 (2003) 222. [3] T. Ichikawa, Mutual coupling among insect neurosecretory cells with an ultradian firing rhythm, Neuroscience 299 (2001) 73. [4] S.N. Baker, E.M. Pinches, R.N. Lemon, Synchronization in Monkey motor cortex during a precision grip task. II. Effect of oscillatory activity on corticospinal output, J. Neurophysiol. 89 (2003) 1941. [5] S.N. Baker, E. Olivier, R.N. Lemon, Coherent oscillations in monkey motor cortex and hand muscle EMG show task-dependent modulation, J. Physiol. 501 (1997) 225. [6] S. Salenius, K. Portin, M. Kajola, R. Salmelin, R. Hari, Cortical control of human motoneuron firing during isometric contraction, J. Neurophysiol. 77 (1997) 3401.
ARTICLE IN PRESS 270
F.A.O. Cruz, C.M. Cortez / Physica A 353 (2005) 258–270
[7] B.A. Conway, D.M. Halliday, S.F. Farmer, U. Shahani, P. Maas, A.I. Weir, J.R. Rosenberg, Synchronization between motor cortex and spinal motoneuronal pool during the performance of a maintained motor task in man, J. Physiol. 489 (1995) 917. [8] W.O. Friesen, G.D. Block, What is a biological oscillator?, Am. J. Physiol. 246R (1984) 847. [9] J.P. Donoghue, J.N. Sanes, N.G. Hatsopoulos, G. Gaal, Neural discharge and local field potential oscillations in primate motor cortex during voluntary movements, J. Neurophysiol. 79 (1998) 159. [10] F.C. Hoppensteadt, E.M. Izhikevich, Synaptic organizations and dynamical properties of weakly connected neural oscillators, Biol. Cybern. 75 (1996) 117. [11] J.M. Kilner, S. Salenius, S.N. Backer, A. Jackson, R. Hari, R.N. Lemon, Task-dependent modulations of cortical oscillatory activity in human subjects during a bimanual precision grip task, NeuroImage 18 (2003) 67. [12] V.N. Murthy, E.E. Fetz, Coherent 25-Hz to 35-Hz oscillations in the sensorimotor cortex of awake behaving monkeys, Proc. Natl. Acad. Sci. USA 89 (1996) 5670. [13] K. Pearson, Motor systems, Curr. Biol. 10 (2000) 649. [14] E. Marder, Neural modulation: following your own rhythm, Curr. Biol. 6 (1996) 119. [15] J.E. Skinner, C.D. Yingling, Regulation of slow potential shifts in nucleus reticularis thalami, Electroenceph. Clin. Neurophys. 40 (1976) 288. [16] W.O. Friesen, J. Cang, Sensory and central mechanism control intersegmental coordinator, Curr. Opin Neurobiol. 11 (2001) 678. [17] P.G. Vaidya, Monitoring and speeding up chaotic synchronization, Chaos Soliton. Fract. 17 (2003) 433. [18] D.V.R. Reddy, A. Sen, G.L. Johnston, Time delay induced death in coupled limit cycle oscillators, Phys. Rev. Lett. 80 (1998) 5109. [19] G. Ermentrout, Oscillator death in population of ‘‘All to All’’ coupled nonlinear oscillators, Physica D 41 (1990) 219. [20] S.H. Strogatz, Exploring complex networks, Nature 410 (2001) 268. [21] T. Po¨schel, J.A.C. Gallas, Synchronization in the dynamical behaviour of elevators, Phys. Rev. E 50 (1994) 2654. [22] Y. Kuramoto, Chemical Oscillation, Waves and Turbulence, Springer, Berlin, 1984. [23] A.C. Guyton, J.E. Hall, Human Physiology and Mechanisms of Disease, sixth ed., W.B. Saunders Co., Philadelphia, USA, 1997. [24] T. Chang, S.J. Schiff, T. Sauer, J.P. Gossard, R.E. Burke, Stochastic versus deterministic variability in simple neuronal circuits: I. Monosynaptic spinal cord reflexes, Biophys. J. 67 (1994). [25] R.E. Burke, B. Walmsley, J.A. Hodgson, HRP anatomy of group Ia afferent contacts on alpha motoneurons, Brain Res. 160 (1979) 347–3452. [26] L.M. Mendell, E. Henneman, Terminals of single Ia fibres: location, density and distribution within a pool of 300 homonymous motoneurons, J. Neurophysiol. 34 (1971) 171. [27] A. Bu¨schges, A. El Manira, Sensory pathways and their modulation in the control of locomotion, Curr. Opin. Neurobiol. 8 (1998) 733. [28] R.J. Nelson, Interactions between motor commands and somatic perception in sensorimotor cortex, Curr. Opin. Neurobiol. 6 (1996) 801. [29] M. Khamsi, Impulse functions: Dirac function. S.O.S. Mathematics, 1997. http://www.math.utep.edu/sosmath/diffeq/laplace/dirac/dirac.html. [30] B.R. Martin, Statistics for Physicists, Academic Press Inc., London, 1971. [31] C. Godin, M.B. Gordon, J.D. Muller, SpikeCell: a deterministic spiking neuron, Neural Networks 15 (2002) 873. [32] P. Poirazi, T. Brannon, B.W. Mel, Pyramidal neuron as two-layer neural network, Neuron 37 (2003) 989.
