Cerebellar Model Controller with new Model of Granule Cell-golgi Cell Building Blocks and Two-phase Learning Acquires Multitude of Generalization Capabilities in Controlling Robot Joint without Exponential Growth in Complexity

ABSTRACT


INTRODUCTION
Artificial neural networks are an important branch of information processing. It started with inspiration from living neuronal circuits present in living organisms, as means for solving problems of technical importance that seem to be handled so easily by them. Attractive was the way of solving problems by mimicking the learning process, without explicitly dealing with the complexity of the problem. Despite initial enthusiasm followed by disappointment and a period of almost abandoned artificial neural network research, it revived and developed to powerful tool for solving wide range of technical problems.
Cerebellum played an important role in this process. It gave hope to decipher internal functional principles of neuronal circuits [1], based on regular and simple layered arrangement of neurons across whole cerebellar cortex, with well defined inputs and outputs, seemingly as uniform repetition of same elementary cerebellar neuronal circuit between groups of neurons. After initial information about anatomy and physiology, target corresponding granule cells. Moreover coding in two sets of mossy fibers is different, being rate coded for those targeting lower dendrite tree. Functionally this input will serve to control multiplicatively activity of granule cell outputs forming circuit with corresponding Granule cell. Existence of these segregated and functionally distinct sets of mossy fibers was later suggested in [22] as explanation for in vivo and in vitro experiments they performed. They also name function of segregated connectivity as "parallel feedforward inhibition", contrasted to function of stereotyped connectivity as "classical feedforward inhibition". This new model puts polynomial processing of granule cell layer as biologically plausible and performed inside the same layer [23]. It may also contribute to sparseness of activity issue at this layer that was recently challenged [24], [25].
In this paper granule cell layer is considered as collection of granule cell-Golgi cell building blocks. Capability of this block in generating multi-dimensional receptive fields modulated by separate input coming to lower dendrite tree of Golgi cell was evaluated with simulations using Simulink single compartment spiking neuronal model. Generalization capabilities of resulting averaging cerebellar model controller with two-phase learning will be explored on a robot joint plant in computed torque configuration.

RESEARCH METHOD
Cerebellum, or little brain, is present in all vertebrates and in most cases weights about 10 percent of total brain weight. On the other side it contains more than half of the total number of neurons in the brain [26]. Functionally it is associated and is considered as part of motor structure. It is important in maintaining posture and balance, in fluent execution of voluntary articulated movements, in motor learning, and based on connections it makes with other parts of the brain it may be involved in many other functions, up to the cognitive level [27], [8].
Spiking model of cerebellar granule cell layer for Simulink was developed to assess the ability of granule cell-Golgi cell building block in generating modulated higher-order receptive fields. Averaging functional equivalent is used to build new averaging cerebellar controller [28], [23]. Robot plant [28] for testing global generalization property to joint speed and joint acceleration is presented.

Neuronal Circuit of the Cerebellum
In contrast to number of neurons contained in the cerebellum [26], number of neuron types is small, common to all vertebrates. It hosts two extreme types of neurons by size, Purkinje cell and granule cell. First is known as largest neuron in the brain, while second one as smallest and most numerous type. All of them are arranged in highly organized three layer structure, named cerebellar cortex. Middle layer one-cell thick contains only Purkinje cells. Next layer in the inner side mainly inhabits granule cells and Golgi cells. Third, outer layer of cerebellar cortex is known as molecular layer. Typical for this layer is presence of huge number of parallel fibers that cross at right angle through flat shaped dendrite trees of Purkinje cells. Parallel fibers are in fact granule cell axons that rise vertically from granule cell layer, through Purkinje cell layer, and when reaching the molecular layer bifurcate in T-shaped form. This description of granule cell axon is from final mature state point of view. Biologically route is opposite [29], from neurogenesis at the temporary external granular layer, through migration on molecular layer, to final state in the granular layer. All circuit formation process is very well orchestrated in time and space, resulting in this sophisticated functional structure, initially thought simple to discern, but it continues to surprise scientific community with finer structural details and molecular diversities [10]. Present in molecular layer are also basket cells and stellate cells.
Simple sketch in Figure 1 (modified from [21]) gives connections between cerebellar neurons. There are only two routes where information can enter cerebellum, mossy fibers and climbing fibers. Mossy fibers form many contacts with granule cells, Golgi cells, and deep cerebellar nuclei. Climbing fiber splits in several branches, and each branch makes very strong connection (composed of hundreds of contacts) with Purkinje cell, but every Purkinje cell in mature state makes exclusive connection with only one climbing fiber. All climbing fibers originate from inferior olivary nucleus residing outside of the cerebellum, but that is reciprocally connected with deep cerebellar nuclei. Contacts between mossy fibers, granule cell dendrites, and Golgi cell axons are done within glomeruli.
Sole output from the cerebellar cortex is through Purkinje cell axons. Most of them, after sending collaterals back to cerebellar cortex, terminate to deep cerebellar nuclei, and from them to other parts of the brain. Remaining part of them make similar neuronal circuits with vestibular nuclei. Groups of cerebellar neurons with microcircuitry as in Figure 1 are organized at higher level in microzones [30], and even higher to zones, that are thought to have specific functions depending on the input and output connections they receive from, and send to other parts of the nervous system. There are very specific attributes of connections between actors of the cerebellum. First, very high divergence of mossy fiber input information, from millions of mossy fibers to billions of the granule cells [27]. Second, convergence from hundreds of thousands parallel fiber Ultrastructure of connections between different actors of the circuit is part of intensive research, and could have important implications for possible attainable functional variations. In this sense mossy fibers, according to [21], are separated into groups that can target differentially specific structures and/or specific locations on them, based on actor and location specific molecular markers [10]. Set mf1 of mossy fibers is common one that target granule cells. They can target Golgi cells as well without consequences for new model. Standard model assumes both targets. Specific to this model is set mf2 and mf3 of mossy fibers that target lower dendrite tree of Golgi cells and deep cerebellar nuclei, respectively, but not the corresponding group of granule cells. Specific functional importance was attributed for set mf2 mossy fibers in [21]. Later, authors of [22] according to their findings suggest that separate sets of mossy fibers may target neighboring granule cells and Golgi cells that are also functionally distinct.

Simulink Spiking Model for Granule Cells and Golgi Cells
To analyze information processing at granule cell layer, spiking models for two most wide spread inhabitants of this layer was developed in Simulink, i.e. for granule cells and Golgi cells. Adopted model is single compartment modified integrate-and-fire that computes membrane potential from [31], with parameters from the same publication This is a first order differential equation that accounts for chemical synapses (for AMPA, NMDA, and GABA receptors) and leaking (resting conductivity Grest). C m is membrane capacitance. Nature of each member in (1) is current and represents postsynaptic currents caused from ions passing through corresponding ion channels. g AMPA , g NMDA and g GABA are conductances, and E AMPA , E NMDA and E GABA are reversing potentials for AMPA, NMDA and GABA receptors, correspondingly. E resr is resting membrane potential. Conductance of NMDA channel is dependent on membrane voltage [32] and is given as  (2) with constants equal to α=62V -1 , [Mg 2+ ]=1.2mM and β=3.57mM. Ion channel conductances are only time dependant, like in spike response model [33], and are given as follows where, t s denotes spike arrival time at presynaptic site. ICC stands for type of ion channel conductance, for AMPA, NMDA, or GABA. Conductances are modeled as falling exponential functions with maximum values g AMPA (t s ), g NMDA (t s ) and g GABA (t s ), and with time constants τ AMPA , τ NMDA and τ GABA , correspondingly. First two conductances will result with excitatory and the last one with inhibitory postsynaptic current. Data used for parameters in equations for granule cell and Golgi cell model are the one given in Table 1 and Table 3 of [31]. Authors of the paper have collected these data from many sources, as mentioned at their Table  1. Figure 2 shows Simulink model for granule cell. Golgi cell Simulink model is similar, but it has two sets of AMPA receptors, for upper and lower dendrite tree synapses, lack of NMDA receptors, and presence of leak current and spontaneous oscillation [31,34]. At left, first two blocks model excitatory postsynaptic currents, i_AMPA and i_NMDA, genarated from presynaptic pulse arriving at port 1 (Exc) from mossy fiber. Third block models inhibitory currents i_GABA from pulse at port 2 (Inh). Forth block models leaking current i_leak. Model of somatic action is at the bottom right, with submodels for integration of currents, generation of spike, and generation of refractory period Model was built with objects that will allow processing of multidimensional inputs signals. Each channel conductance is modeled according to upper simulation diagram. Four blocks at left model four sources of synaptic currents according to (1) submodels: for integration, spike generation and for refractory time. Spike is idealized as rectangular pulse of unitary amplitude and with width (Sw) equal to 0.1 ms. Gain of ion channel conductance is scaled according to these values. Propagation of spikes in the axon is not modeled. These models are encased inside two types of neuron subsystem blocks with different processing of excitatory synaptic currents, Figure 3. First one, GrC_ns, is model of a single neuron with n synaptic inputs. Second one, GrC_nn, is model for n neurons. Summation of spikes targeting dendrite tree for later case is modeled out of neurons block. Number of synapses in the first model and number of neurons in the second model will be determined by dimensionality of corresponding input signal. Golgi cell will use model for single neuron with many synapses resulting to GoC_ns subsystem model. It will have two excitatory inputs, from upper dendrite tree (udt) and lower dendrite tree (ldt). Inhibitory inputs to Golgi cell [35] are not modeled. Figure 3. Two forms of granule cell subsystem blocks with different processing of excitatory currents i_AMPA and i_NMDA to i_Exc: (a) processing that will generate model of single neuron with many napses, subsystem GrC_ns; (b) processing that will generate model for many neurons. Number of synapses n the first model and number of neurons in the second model will be determined by dimensionality of i_Exc

Simulink Spiking Model for Granule Cell-Golgi Cell Circuit
Spiking models of granule cell and Golgi cell are used to build Simulink spiking model of granule cell-Golgi cell circuit shown in Figure 4, where GrC_nn is model for n granule cell neurons, GoC_ns is model for single Golgi cell with n synapses. Summation block with four inputs, x1_L1 to x1_L4, represents granule cell dendrites receiving inputs from four groups of population coded mossy fibers. Lower mossy fibers correspond to one rate coded mossy fiber v1_R that connects to lower dendrite tree ldt of Golgi cell. Synaptic strengths are assumed to be constant, with nominal value [31] set at corresponding neuron model. Gain blocks S_MF-GrC, S_MF_GoC, and S_GoC-GrC, currently set to 1, represent relative synaptic strength. They can be used to explore the range and effect of variable synaptic strengths. If synaptic strengths of different inputs to granule cell are different, each input (of four of them in Figure 4) may have its own gain block S_MF-GrCi, and by removing S_MF-GrC block after summation block. Synaptic strength S_GrC-GoC is used to correct effective synaptic strength [31] between granule cell axon and upper dendrite tree of Golgi cell in accordance with number of granule cells targeting Golgi cell. Parallel fibers contain processed multidimensional receptive fields (in case of two-dimensional input space they will be two-dimensional receptive fields, shown later in Figure 5). Generation of information with population coded receptive fields that drive mossy fibers is thought not to be part of the cerebellum and corresponds to pre-cerebellar structures.

Mossy Fibers and Information Entering the Cerebellum
There are two routes where information can enter the cerebellum, mossy fibers and climbing fibers. In most cerebellar models later one are assumed to carry error information. Initial model and most of model variations following assume that information carried by mossy fibers is population coded [6,18]. In this case unique input states are represented with unique activity in a subset of mossy fibers. Simplest one used is in form of several one-dimensional layers composed from nonoverlapping receptive fields, with their activity described by rectangular basis functions. In this case each layer will have only one of them active for certain range of input values, where active means for output being equal to one and zero for inactive case. This form of coding was proposed by Albus [6] and layers are known as Albus overlays, Figure 5. Later layers were joined to a single layer of overlapping receptive fields, with extension to many shapes of basis functions describing receptive field activity. Another view is thinking of each receptive field as one layer. Regardless of how we view them, each signal (usually one-dimension, but not necessarily) will be coded with a group of mossy fibers where subset of them (one or more) will be more or less active in a specific range (typically form 0 to 1). Mossy fibers carrying rate coded information is common to some later versions of cerebellar [19,20] or hybridized fuzzy-cerebellar models [18]. Generation of population coded spiking mossy fibers from input signal is done with two processing stages, upper part in Figure 6. First stage from input signal will generate a number of new signals with selected basis function, where square, triangular and Gaussian being most common.
Each signal is then converted to stream of pulses (spikes) at second stage, where pulse frequency is proportional to amplitude of input signal. Rate coded spiking mossy fibers are generated by skipping the first stage, i.e. by directly feeding input signal to second processing stage, lower part in Figure 6. As functional model of the second stage a modified spiking neuron model from Section 2.3 was used. Modifications consist on keeping only somatic action model, removing refractory period subsystem, and treating the input signal as excitatory synaptic current with proper scaling, resulting in PulseGen (pulse generator) subsystem block. Parameters were set to have mean output pulse frequency equal to control value at input. Internally a noise generator was added, with parameters adjustable through mask of the corresponding block. Each pulse generator block will have different seed for noise generator for uncorrelated noise between signals. Processing of the first stage is done by passing input signal through center shifted group of selected basis function, block 1L_SBF for one layer of square basis functions in Figure 6.
Specific to this cerebellar model is that population coded mossy fibers (mf1 in Figure 1) will target granule cells, while separate set of rate coded mossy fibers (mf2 in Figure 1) will target lower dendrite tree of Golgi cell [21,28,23]. Input signals x1 and v1 represent any signal to be processed by cerebellum. Typically for servomotor driven joint x1 will be joint position, while v1 will be joint speed or joint acceleration. Receptive field width is 3 quanta. For dimension x a layers are ADG, BEH and CFI, and layers for dimension x b are adg, beh and cfi. Input at x a dimension will make active receptive fields D, E, and F, whereas input at dimension x b will make active receptive fields g, h, and f. For active fields corresponding to input (x a , x b ) there will be 9 two-dimensional receptive fields, Df, Dg, Dh, Ef, Eg, Eh, Ff, Fg dhe Fh. Positions of the two-dimensional receptive fields have been shifted a little from their real position determined from corresponding one-dimensional receptive fields to make them more distinctive Figure 6. Model for generating population and/or rate coded spiking mossy fibers. Upper part generates population coded spiking signal x1_L1 from input signal x1. Block 1L_SBF generates 6 signals with square basis function. Amplitude of these signals is converted to corresponding stream of pulses (spikes) with pulse generator block PulseGen L1. Lower part generates rate coded spiking signal v1_R from input signal v1, by directly feeding input signal to the input of pulse generator block PulseGen R

Averaging Cerebellar Model Controller with Granule Cell-Golgi Cell Building Blocks
Final averaging model of the cerebellum will have in the first processing layer several granule cell-Golgi cell building blocks. Each of them generates a group of higher-order receptive fields by using multiplication operator, for all possible combinations from a set of lower-order input signals. Blocks have ability of modulating group output activity by separate modulating input signal [21]. Input signals are functional equivalent of mossy fibers. Set of mossy fibbers mf1 that target directly granule cells is expected to be with population code. Other set mf2 that targets Golgi cell lower dendrite tree (modulating input signal) is expected to be rate coded. Standard CMAC is special case of this cerebellar model with modulating input set to constant, typically 1.
Outputs of building blocks will be functional equivalent of parallel fibers. Next stage is composed from one or more Purkinje cells blocks, modeled as standard perceptron with linear activation function, and a number of synaptic weight that is determined from dimensionality of signal (parallel fibers) entering the block. Purkinje cell blocks receive teaching signals (equivalent of climbing fibers) that can be same or different for all blocks, or in any combination. Each output will have one neuron of deep cerebellar nuclei as last stage, modeled as simple summation, receiving outputs from corresponding Purkinje cells. All signals may have positive and negative values, contrary to biological equivalent, but common to most artificial neural networks. Other sites of synaptic connection that may be adjustable [36], beside parallel fiber-Purkinje cell, are not modeled. For biologically consistent structural model for some of them adjustability would be needed to bring model at the operating point, or some interdependency fulfilled. With new functionality assigned to Golgi cell, including new functionalities to stellate cell and basket cell from [21] this cerebellar model was referred to as CMACgbs [28], [23].
Simulation model of CMACgbs is given in the Figure 7 for single output. Model contains preprocessing blocks for generating population coded mossy fibers of selected shape (square, triangular, or Gaussian) from input signals, transformation common for most of these models. Functionally these blocks are thought to reside outside of the cerebellum. These blocks and granule cell-Golgi cell blocks are sharable components [23]. If teaching signal is same for all Purkinje cell blocks they can be joined to single Purkinje cell block, making deep cerebellar nuclei unnecessary. Standard CMAC follows this route. Gaining knowledge with better generalization may require different teaching signals for different Purkinje cells, similar to biological case where single climbing fiber (teaching signal) targets only several Purkinje cells and each mature Purkinje cell is targeted from single climbing fiber. CMACgbs preserves this structural property. Granule cell-Golgi cell blocks have four inputs matching typical number of dendrites that granule cells have. If receptive fields of lower dimensionality are required unused inputs are tied to constant 1, upper block in Figure 7. If pure rate coded signal is to be passed to Purkinje cell all inputs can be tied to constant, choice offered by switches for two next blocks. In simulations there is no need to have these two blocks at all, but biologically there is no direct route from mossy fibers to Purkinje cells.

Robot Plant
Robot plant with one rotary link, mounted on at arbitrary positioning and orienting base, was used as control plant for testing generalization capabilities of CMACgbs. General dynamics model [37] describing this robot in disaggregated form [21] is given by following equation [28] , where: g x0 , g y0 , and g z0 are projections of gravity vector to 0-th coordinate system; G 1r1 (q 1 ), G 1r1 (q 1 ) and G 1r1 (q 1 ) represent gravity loadings cause from robot link and gripper; d 11r represents inertial loadings, constant in this case; F 1d and F 1s represent dynamic and static friction; m L is the load mass. Other similar terms in (4) have similar meaning, but are related to load mass. G and d terms are in general nonlinear functions dependent on joint position. If load is constant and its effects are accounted within link and/or gripper, or there is no load, four last terms can be omitted. Dynamics of base movement and rotations are not included, assumed to be slow enough to be treated as static. Dimensionality of the input space depending on the dynamical effects that will be accounted for ranges from 1 to 7, being 1 for link and gripper gravity loading only, 2 for standard robot in fixed base without friction and 3 including friction effects, 5-6 for previous case with different base rotation configurations, and 7 for full model including variable mass. For pendulum-like robot with rod-like link, gripper modeled as point mass, and including viscous friction (4) becomes , ) sin( 2 1 3 where m S and l s are link mass and length. Simulation model of the robot [23] was developed in SimMechanics. As main controller Proportional-Derivative (PD) controller was used.

RESULTS AND ANALYSIS
First processing of information entering the cerebellum by mossy fibers is done in the granular layer, with circuits formed by granule cells and Golgi cells, and result is passed to molecular layer by granule cells axons for further processing. Granule cells axons are composed of two distinctive parts: ascending axon and parallel fiber, aa and PF in Figure 1. Vertical ascending axon part starts from granule cell and ascends up to some level of the molecular layer, where it bifurcates in form of letter T and creates parallel fibers that travel some distance in both directions. Outputs from group of about 5000 granule cells will make synaptic connections with one Golgi cell [27,23], which in turn will inhibit same group of granule cells. This feedback inhibition is done through synaptic connections inside glomeruli, Gl in Figure 1, where set mf1 of mossy fibers make excitatory synaptic connections with granule cell dendrites. Second route with functionally distinct information comes by set mf2 of mossy fibers to Golgi cell through its lower dendrite tree ldt, with feedforward inhibitory effect on granule cells. This would be one idealized elementary granule cell-Golgi cell building block.

Modulated Higher Order Receptive Fields with Granule Cell-Golgi Cell Building Blocks
Lower dimensional information at the input of the cerebellum carried by mossy fibers is first processed by granule cell layer. It will expand the information space by several orders of magnitude by generating higher-order receptive fields that serve as associations for next processing steps, termed by Albus as expansion recoding [3]. Typically in simulations inputs are one-dimensional information and processed information carried by axons of granule cells (ascending axons and parallel fibers) at the output of granule cell layer are multi-dimensional.
How effective is spiking single granule cell-Golgi cell elementary building block in generating higher-order receptive field from one-dimensional inputs is explored. Simulations were done with model from Figure 4. All granule cells in the model have same values for parameters [31], and each of them has four inputs (basal dendrites). According to idealized model, only one Golgi cell is used. It receives excitatory synaptic inputs from all granule cells in the block and sends inhibitory synaptic outputs to the same granule cells. All but one relative synaptic strength gain blocks in the model are set to 1.
The one different is gain S_GrC-GoC set to 1/(number of granule cells targeting Golgi cell) for proper scaling of synaptic strength [31]. Input space is two dimensional composed from two one-dimensional input signals. Two layers of square receptive fields are used for each input signal, similar to the one shown in Figure  5. In this way there will be four groups of nonoverlapping receptive fields that carry information for inputs of granule cells. Layers are composed from 2 receptive fields 2 quanta wide. Four combinations, from 9 possible, were implemented targeting 4 granule cells. Model from Figure 6 is used to generate spiking mossy fibers from receptive field activity. To explore whole input space two sawtooth signals are used. Frequency of x1 signal (horizontal axis) is several times higher than that of x2 signal (vertical axis), depending on required vertical resolution. This simulation used 1 Hz and 0.01 Hz frequencies. Input space is scanned in [-8, 8] range, but effective used space is in [-4, 2] range, for both dimensions. resulting two-dimensional receptive fields are marked with squares on upper-left plot. Green square is center of active receptive field. It can be seen that activity of this receptive field is almost suppressed at centers of other receptive fields marked with red squares. At centers of two-dimensional receptive fields total output activity is dominated mainly from activity of single corresponding receptive field, as is shown in four central plots of Figure 8. Number 1 in squares denotes that total output activity will be dominated from one receptive field. Other numbers, 2 and 4, indicate that total output activity will be shared between 2 or 4 neighboring active two-dimensional receptive fields, but each contributing receptive field being proportionally less active.
Low-pass filtered output spike activity of a single granule cell at similar working conditions is given in Figure 9 in two views, top and oblique view (axis scaling for this figure is different from the one in Figure 8). Intensity is color coded, with red being most active and dark blue being inactive state. Looking at excitation and inhibition that target granule cell, response in Figure 9 is qualitatively similar to center surround response [38], where at the center excitation dominates the response while the surrounding is dominated by inhibition, mostly caused by activity of other granule cells in the block. Effective used space is marked with white square in Figure 9(a).
Relative suppression of activity between active receptive fields is handled by feedback (upper dendrite tree) inhibition from Golgi cell, controlled by total activity of all granule cells in the block. This can be viewed as gain control functionality [3,21] that tries to keep total output activity of the group to constant level. Artificial neural networks often use normalization of excitation, either by selecting proper receptive field basis functions having this property inbuilt [39], or by having it as additional processing [23]. Normalization aids learning and would be equivalent to gain control. Other view would be that total activity at the output of the granule cell-Golgi cell building block will move the firing threshold of granule cell accordingly, leading to suppression of less active and enhancing more active granule cells. These suppressed less active granule cells are treated as noise that should be removed from signal, i.e. granule cells with higher activity. If Golgi cell inhibition is not present activity at the output of single granule cell would result with higher-order receptive field typical for neurons with summation processing of incoming inputs, shown with inset at left in Figure 8. It has high activity at side lobes, horizontal and vertical, where only one or two inputs of granule cells are active. These side lobes are suppressed inside effective used space where there are other more active neurons at corresponding positions of input space. Outside of the effective used space side lobe activity reappears, seen in four central plots of Figure 8 and Figure 9.
Typical operator used for creating higher-dimensional receptive fields from lower-dimensional one in case of artificial neural networks is multiplication. Square receptive fields at inputs would result with hypercube (or hyper-parallelepiped) of corresponding dimensionality. In case of square receptive fields with 0/1 activity same result would be obtained with AND operator. For given excitation in this simulation result would be uniform activity inside square denoting center of the resulting receptive field (seen in Figure 8, or only red zone in Figure 9), and zero otherwise. Activity would code only small portion of the input space, resulting in high number of neuron to cover whole input space, possibly leading to combinatorial explosion. Wider square receptive field would relax this problem, but resolution would become worse. Wider receptive fields with processing proposed by Albus [6], shown in Figure 5 for full overlaid two-dimensional CMAC case, will give better generalization and preserve the resolution. Different shapes of receptive fields for input information, like triangular or Gaussian, is other way for tackling the problem [39]. Combination of both is common to later modified models of CMAC [18][19][20]. Processing with granule cell-Golgi cell building block with feedback from upper dendrite tree for same input excitation gives more graded response. It resembles collective response of all corresponding two-dimensional Albus overlays (seen in Figure 5), or to stepped approximation of two-dimensional receptive field that would result from triangular receptive fields at input processed with multiplication operator, or min operator in fuzzy networks case. Higher order receptive fields at input will result with smother high-dimensional receptive fields. Intensity is color coded. Effective used space is marked with white square How activity of receptive fields changes in time for a possibly real input trajectory encoding is shown in Figure 10. Input signals used are x1=4·sin(2π·1·t) and x2=3·sin(2π·2·t), that generate trajectory for drawing number 8. Spike activities at the output of four granule cells are marked with colored dots. Signal amplitudes were selected such that part of trajectory will travel through space not covered by two-dimensional receptive fields. Numbers inside squares denote how many receptive fields will dominate encoding. At parts of the input space that are not covered by corresponding receptive fields, side lobes will appear caused by lucking inhibition from missing receptive fields. These activities will send incorrect associations to other processing layers. Two such cases are encircled in Figure 9, with arrows pointing to correct association. Another example for onedimensional trajectory can be found in [23]. Figure 10. Spiking output activity of 4 two-dimensional receptive fields at 4 granule cells with "number 8" trajectory as input: x1=4·sin(2π·1·t) and x2=3·sin(2π·2·t). Activity of two receptive fields (blue and violet) is shifted to the right for 0.3 units for better visibility, by avoiding overlapping. Encircled are wrong output associations that lie outside of the effective used space that have equivalent associations inside that space. Numbers inside squares denote how many receptive fields will dominate encoding. Dot of the specific color denotes presence of spike at corresponding output, space, and time. Four central plots in Figure 8 are generated without activity at mossy fibers targeting lower Golgi cell dendrite tree. Increasing excitation on this input will increase inhibition to all granule cells of this elementary circuit, and multi-dimensional receptive field activity will get proportionally lower, and vanish for some higher input level. Inset at the bottom-right in Figure 8 shows activity at the output of one granule cell when lower dendrite tree of the Granule cell is active. We notice that center activity is still present, but mean activity at any place of the input space became lower. Golgi cell upper dendrite tree (feedback) and lower dendrite tree (feedforward) inhibition of one-dimensional trajectory can be seen in Figure 5 of [23]. In idealized case signal targeting lower dendrite tree of Golgi cell (mf2 in Figure 1) is assumed to modulate multiplicatively amplitude of resulting set of higher-order receptive fields [21,28,23]. If inhibition through lower dendrite tree is not used (not active), or it follows stereotypical configuration being connected to same group of mossy fibers that target granule cells of the same building block, output activity would be determined solely from activity of the mf1 group of mossy fibers as shown in Figure 1. Under any of these conditions model reverts to classical one, generating unmodulated higher-order receptive fields.
Functional result of processing in granule cell-Golgi cell building block will be a group of multidimensional receptive fields, whose activity can be modulated with input at lower dendrite tree of Golgi cell. Without activity in the modulatory input, activity of resulting multidimensional receptive fields would be maximal, and determined from mean total activity at group mf1 of mossy fibers. This would produce high activity on Purkinje cell, that finally would inhibit deep cerebellar nuclei, and with that output from cerebellum. This could be balanced and brought to operating point with same group mf1 of mossy fibers targeting also deep cerebellar nuclei. This condition would produce output results in function of synaptic strength between parallel fibers and Purkinje cell, assumed as only site with synaptic weights where learning can take place. In addition to that, current activity for certain point in the input space can be directly modulated with input at lower dendrite tree of Golgi cell. Increased input would inhibit granule cell activity that in turn would decrease Purkinje cell activity, causing less inhibition to deep cerebellar nuclei, resulting with increased final output activity. From this inhibition of inhibitory neuron, effective constant of proportionality is positive.

CMACgbs with Two-Phase Learning as Robot Controller and Generalization
Pendulum-like robot plant, with fixed base including viscous friction effects, described in Section 2.6 was used as control plant for exploring generalization capabilities of averaging cerebellar controller with new model of granule cell layer, composed from several granule cell-Golgi cell building blocks, given in Figure 7. Parameters of the robot for following simulations are: link mass m S = 1 kg, link length l S = 1 m, and gripper mass m G = 0.2 kg. Proportional-Derivative (PD) controller, with K p and K v set to 30 and 4, was used as main controller. Cerebellar controller will serve as secondary learning controller for robot dynamics, with one granule cell-Golgi cell building block for each term of (4) that will be accounted by cerebellar controller. It can be in feedforward or in computed torque configuration [40], depending on the source of position and speed information.
In following computed torque configuration is used. Actual values of joint position and joint speed, and desired joint acceleration are used as inputs to the cerebellar controller, signals q_a, dq_a and d2q_d in Figure 7. Since base is fixed gravity projection inputs will be constant, and grav_x input in Figure 7 will be tied to constant equal to 9.81. This configuration of the controller will set dimensionality of the input space to 3. Only position information will be transformed to population coded information with 73 triangular receptive fields of 4 quanta wide, and will serve as mossy fibers that target mf inputs of granule cell-Golgi cell building blocks. Block 1L_1D_TBF (1-Layer of 1-Dimensional Triangular Basis Functions) in Figure 7 performs this function. Two other input signals, joint speed and joint acceleration, will be rate coded and serve as mossy fibers that will be tied to mf_ldt (Golgi cell lower dendrite tree input) of granule cell-Golgi cell building blocks. Since controller will be implemented in digital form, calculated output value will be held constant until new calculation is done with zero-order-hold (block ZOH_CMACgbs in Figure 7), where cerebellar controller update period can be set, being 0.01 second for simulations in this paper. Outputs of the PD controller and cerebellar controller are summed together and result is used to drive simulated robot.
CMACgbs controller was trained on sinusoidal trajectory with cycling from 0 → 2π → 0 → -2π → 0 in selected time. Training was done in two phases. First phase was done with very slow motion, 100 s / cycle, during which dynamics is dominated with gravitational loading. During this phase only gravity term is trained by connecting error signal (output of the PD controller) to teaching input CF (climbing fiber) of Purkinje cell block PC_Gx as shown in Figure 7. Teaching for two other Purkinje cell blocks is disabled by feeding 0 to their CF inputs. Progress of learning and final result is given in Figure 11. Position error in Figure 11(a) shows significant decrease of error during first cycle of learning (0-100 s), and some minor adjustments during second cycle (100-200 s). Remaining error of cosinusoidal form is joint speed related and cannot be learned with only joint position related term, and is caused from viscous friction torques. In fine scale error will change in rhythm

4305
of weight training (same as cerebellar controller update period, i.e. 100 trainings per second) and with amount controlled by learning rate (equal to 0.0001), but it is hard to be noticed in a present scale of Figure 11. Average update during one cycle is close to zero, causing error to be cycle-by-cycle similar during subsequent cycles. High value of learning rate may lead to smaller error when learning is on, but with a risk instability (especially at the beginning of learning when errors may be quite large) or of ending with less accurate cerebellar model when learning is set to off, similar to behavior of adaptive controllers. Inset in the Figure 11(a) shows testing error for one cycle with learning turned off, as measure of success in acquiring plant model [41]. Note vertical axis scaling (position error axis in the range from -0.01 to 0.01) to make details visible. Learned synaptic weights for gravitational term are pictured as heat map in Figure 11(b), with apparent sinusoidal relation of weights from joint position. During learning this term there could be interference from inertial effects, but considering maximal values for inertial and gravitational torques to be 0.0132 Nm and 6.867 Nm, it is expected to be negligible. These weights are not adjusted any more, since other circumstances with higher speed are less informative regarding this term. Resulting receptive fields from granule cell-Golgi cell block GrC-GoC_Gx are common unmodulated ones, since modulating signal is tied to constant, as in standard CMAC model [42,28]. During second phase cycling was speeded up to 10 s/cycle, and only inertial and friction terms were adapted. Since gravity loading is accounted properly from results of first phase of learning, current control error will be caused by inertial or viscous friction torques. These effects will be compensated with two additional subsystems shown in Figure 7, each composed of one granule cell-Golgi cell and one Purkinje cell block, with proper learning. Two cases will be explored, first one when terms are not joint position dependent and second one when they are. First case is configured with two switches at current position in Figure 7, resulting in only one weight for each term [21].
Second case is for the other position of switches, resulting in 73 weights for each term. There are three blocks for generating population coded receptive fields in Figure 7 for each term, but since input for all three of them is joint position only one can be used and then sharing mossy fibers between three GrC-GoC building blocks, similar to biological mossy fibers that distribute information to several granule cells in reach of their length. Functionally both cases will have modulated receptive field/fields in parallel fibers targeting corresponding Purkinje cell blocks. Output of the PD controller will serve as learning signal for both terms, while for gravitational term learning signal will be tied to zero. Process of learning during 20 cycles, and convergence of weights to ideal values for the first case is given in Figure 12(a). Middle plot shows trend of decreasing joint position error during learning. Both weights converge to ideal values. Process of weight convergence for both terms in second case is shown in Figure 12(b). There are 73 weights for each term. Effective activity of the receptive fields is 2 (triangular receptive fields 4 quanta wide) causing weights to converge to half the ideal value. Since information content of error signal is joint position dependent and term dependent, different weights follow different convergence path, but group trend converges to proper value. Heat map of all weights after second phase of learning is pictured in Figure 12(c). Stripes for inertial and friction synaptic weights in ideal case should be uniform and positive. Stripe for gravity synaptic weights is same one shown in Figure 11(b).

Global Generalization to Joint Speed and Joint Acceleration Inputs
Convergence of weights to ideal values is important for proper global generalization property of cerebellar controller with granule cell-Golgi cell building blocks by rate coded inputs targeting Golgi cell lower dendrite tree. These inputs will be targeted by joint speed and joint acceleration [21,28,23], as described in previous section. Two-phase learning was used to achieve learning of three controller terms without interference between them, as one form of tackling credit assignment problem. Simultaneous learning of thee terms will efficiently decrease error, but may fail to converge to ideal values [28]. Generalization to population coded inputs remains local as in standard CMAC, determined by receptive field width and basis function, and layered processing. Trained cerebellar dynamics from previous section was assessed at same sinusoidal trajectory cycling from 0 → 2π → 0 → -2π → 0 at different cycle times. During this phase learning is tuned off and all terms will participate in control process, to illustrate quality of learned model [21]. Figure 13 gives results of joint position error during four test runs with cycle times equal to 100 s, 10 s, 4 s and 2 s. Fifth dashed plot in the same figure corresponds to test run at 2 s cycle time that was subject to additional second phase of learning at same 2 s cycle time. Reference error is the one with 10 s notation (black line), obtained after learning from Figure 12(a). Other lines give errors for other cycling times. Error increases with shorter cycle times, but it remains at very low values. Additional second phase of training at high speed halved joint position error at same speed, with additional improvements to lower cycle times (not shown in the figure). Since different sinusoidal trajectories will be subject to different joint speeds and joint accelerations from reference one. Seen in the three dimensional input state space with axis as joint position, joint speed, and joint acceleration they are almost completely different trajectories [28]. Standard CMAC generalization would provide almost negligible contribution to all new test trajectories, as shown in [28] for two times faster trajectory than reference. Results of tracking trajectories that span the range from 10 times slower to 5 times faster than reference show very good generalization property for this version of cerebellar controller. Errors could be made even smaller by increasing update rate of the controller, but update period of 0.01 s is selected for practical reasons when realtime control is to be implemented.

Int J Elec & Comp Eng
ISSN: 2088-8708  Figure 13. Position error in tracking sinusoidal trajectory, cycling from 0 → 2π → 0 → -2π → 0, with PD+CMACgbs controller. Reference error is the one with 10 s notation (black line). It was obtained with 5 cycles of learning with 100 s / cycle for gravitational loading terms only, followed with 20 cycles of learning with 10 s / cycle for other two terms. Testing was done with the same trajectory (normalized to 2 s) at different speeds of cycling. Generalization is very good for wide range of changes in cycle time, showing good generalization over different joint speeds and joint accelerations. Error during 2 s cycle was halved with one additional second phase of learning at same 2 s cycle time (dashed blue line) Extending dimensionality of this cerebellar model controller with other terms in (4) can be done adding one granule cell-Golgi cell building block and one Purkinje cell for each new term [23]. Purkinje cells with same learning signal can be joined into one. For robots of up to four dimensions, extension is straightforward: precerebellar blocks that create one-dimensional receptive fields are added. They will be connected to other inputs of granule cell-Golgi cell building block (tied to constant 1 if not used), and resulting receptive fields will be of corresponding dimensionality.
One-dimensional receptive fields for all dimensions can be shared between any number of granule cell-Golgi cell building blocks. If dimensionality of input space is higher than four, results of [43] for approximation of frequency band-limited high-dimensional receptive fields can be approximated with full overlaid limited connectivity (up to four in granule cell-Golgi cell building blocks) multi-dimensional receptive fields.

CONCLUSION
Single compartment spiking neuronal model for the granule cell-Golgi cell building block illustrated creation of multi-dimensional receptive fields from one-dimensional receptive fields at granule cell axons (ascending axon and parallel fibers). Separate inputs to the lower dendrite tree of the Golgi cell can be used to modulate group of higher dimensional receptive fields at outputs of the granule cells. Cerebellar model with granule cell layer composed of collection of idealized granule cell-Golgi cell building blocks was used as dynamics controller for robot joint plant, in computed torque configuration. Cerebellar model of this structure, being biologically plausible, spans generalization capabilities to global for rate coded inputs that target lower dendrite tree of the Golgi cell, while they remain local as in conventional cerebellar models for inputs targeting granule cells. Simulations for robot sinusoidal trajectory tracking task over wide range of cycling time was used for exploring generalization capabilities of controller. New cerebellar model controller with two-phase learning, with separate training of gravitational term on very slow cycle time (100 s / cycle) and other terms to some mid cycle time (10 s / cycle), was able to acquire robot plant dynamics model. Testing phase with learning turned off, over wide range of unlearned slower and faster trajectories with cycling times (from 100 s / cycle to 2 s / cycle) performed with very low tracking error. Model can be extended to other dynamical effects, like adaption to any orientation, variable mass load, and contact with environment, without exponential growth in complexity.