Biol Cybern (2014) 108:49–60 DOI 10.1007/s00422-013-0577-z

ORIGINAL PAPER

An activity-dependent hierarchical clustering method for sensory organization Jesús Requena-Carrión · Mark Richard Wilby · Ana Belén Rodríguez-González · Juan José Vinagre-Díaz

Received: 29 November 2012 / Accepted: 4 November 2013 / Published online: 19 November 2013 © Springer-Verlag Berlin Heidelberg 2013

Abstract Biological and artificial sensory systems share many features and functionalities in common. One shared challenge is the management setup and maintenance of sensory topological information. In the case of a massive artificial sensory receptor array, this is an extremely complex problem. Biological sensory receptor arrays, such as the visual or tactile system, face the same problem and have found excellent solutions by implementing processes of sensory organization. Not only can biological sensory organization initiate the topological data construction, it can deal with growing systems and repair damaged ones. Importantly, it can use the patterned activity of sensory receptors to extract topological relationships. Using inspiration from these biological processes, we propose an activity-dependent clustering method for organizing large arrays of artificial sensory receptors. We present an algorithm that proceeds hierarchically by building a quadtree description of sensory organization and possesses many qualities of its biological counterpart, namely it can operate autonomously, it uses the

J. Requena-Carrión (B) Department of Signal Theory and Communications, Rey Juan Carlos University, Camino del Molino s/n, 28943 Fuenlabrada, Madrid, Spain e-mail: [email protected] M. R. Wilby · A. B. Rodríguez-González · J. J. Vinagre-Díaz Department of Mathematics Applied to Information Technologies, Universidad Politécnica de Madrid, Ciudad Universitaria s/n, 28040 Madrid, Spain e-mail: [email protected] A. B. Rodríguez-González e-mail: [email protected] J. J. Vinagre-Díaz e-mail: [email protected]

patterned activity of sensory receptors and it is capable of supporting growth and repair. Keywords Sensory systems · Topographic organization · Clustering · Visual system · Sensor networks

1 Introduction Information processing in sensory systems whether they are biological, or artificial, relies on a detailed description of how sensory receptors are spatially distributed. This results in the requirement of the construction of sensory topographic maps, which are database structures containing receptor location information. Considering systems of tens of thousands, or even millions, of sensory receptors, building topographic maps can be an extremely difficult task. If we add that with this number of components, there is a constantly changing distribution of sensory elements caused by, system failures, repairs, location changes or upgrades; any realistic approach dealing with sensory systems of this magnitude must also provide a mechanism to automatically maintain the topographic structure. The approach of many biological systems to the challenge of building and maintaining sensory topographic maps has been to develop two different mechanisms for sensory receptor organization, namely activity-dependent mechanisms, which use sensory patterned activity and activityindependent mechanisms (Goodman and Shatz 1993). Activity-dependent mechanisms are believed to build a fine, detailed description of sensory receptor organization, by using a coarse description of sensory receptor organization built by activity-independent mechanisms. Both mechanisms seem to operate in parallel during development and also in

123

50

adulthood, providing means to adapt to new external environments and to reorganize in case of severe damages in the sensory system itself (Kaas 2000). In extensive artificial sensory systems, given that it is not practical to manage information of this complexity by a purely manual process, automatic systems must be developed. When considering the problem, it is apparent that there are close parallels between artificial and biological systems. Specifically, as we have already mentioned, artificial systems can be considered to grow and undergo repair, and we need to handle these changes in much the same way as biological systems do. In this work, we will use inspiration gained from the study of distributed sensory systems in the biological world, to inform the construction of algorithms and processes that can be applied to organizing artificial massive distributed receptor arrays. We will devise an organizational method that uses the patterned activity of sensory receptors as an input and that proceeds in a hierarchical manner, producing multiple descriptions of sensory organization ranging from coarse to fine grain. Additionally, we will explore the problem of receptor mis-location detection in the context of sensory organization.

Biol Cybern (2014) 108:49–60

2.1 The problem of sensory organization Sensory systems consist of multiple receptors covering a sensory field. When receptors are deployed, they necessarily have a physical location in space and a set of immediate neighbors, in other words, there exists a map M P of their physical positions. However, we initially may not have information about M P . Typically, all we have is every sensed event that the receptors reported and the time when it occurred. Therefore, if we organized the receptor positions with “pseudo” coordinates to estimate a map M E of the sensed field, we would obtain an uncorrelated pattern that bore no resemblance to the original stimulus. In other words, the original topological structure is lost in the mapping M P → M E . In order to reconstruct the original physical map M P in an internal representation map M R , we must apply an appropriate transformation T to M E . These relationships are shown schematically in Fig. 1. The problem of sensory organization can then be simply stated as follows: for a given physical map M P and a specific unordered map M E , find the transformation T , such that T (M E ) = M P .

2.2 Sensory organization in the visual system

In this section, we explore the problem of sensory organization in biological and artificial systems. Firstly, we formulate the problem of sensory organization. Then, we use the visual system as a paradigm for introducing the processes involved in sensory organization in biological systems. Finally, we identify the main elements and strategies in biological sensory organization that can be exported to artificial sensory organization. These elements constitute the basis for the sensory organization method proposed in this study, which will be described in later sections.

Evolution has provided solutions to the problem of sensory organization in biological systems. The most obvious place from which to draw inspiration is that of the visual system (Nassi and Callaway 2009). This system is far from understood, but there are hints within current knowledge that could greatly assist in the design of algorithms to solve this problem for artificial systems. Light receptors develop in the retina at unspecified points and are connected to the central nervous system by retinal ganglion cells (RGC). The RGC can be considered to generate the physical map that corresponds to the map M P defined in Sect. 2.1. Then, the RGC have to be connected into the central nervous system. In mammals, the RGC connect to

Fig. 1 The presence of an object in the sensory field constitutes a stimulus that activates the receptors in the physical map M P . As the actual positions of the receptors are unknown, only an apparently random

pattern can be observed, as shown in the map M E . This has to be transformed into a map M R = T (M E ) where the original sensed pattern can be observed

2 Background

123

Biol Cybern (2014) 108:49–60

the dorsal lateral geniculate nucleus of the thalamus (dLGN). From here, there are radial connections directly into the visual cortex. Obviously, this is a simplified picture, but the important observation is that in both the dLGN and the cortex, the retina is mapped out in topographic structures known as retinotopic maps (Wandell et al. 2007). Retinotopic maps maintain the topology of the retina, though they seldom maintain the metric. The mapping of the retina to the retinotopic maps typify the connection defined by M P → M R . It is important to note that the mappings M P → M R must grow and do not initially exist even in the mature fetus. Therefore, the question arises, what processes of organization make it possible to connect a single cellular receptor to a specific nerve cell in the visual cortex, when the two are separated by distance of several centimeters? In other words, what processes might lead to the organization of retinotopic maps in the visual system? Experimental studies suggest that there exist coarse-grain and fine-grain levels of organization in the nervous system (Goodman and Shatz 1993). In coarse-grain levels of organization, clusters of receptors are mapped together into their correct target regions. They rely on activity-independent mechanisms for pathway and region identification. Based on this coarse mapping of receptors, activity-dependent mechanisms are capable of refining the mapping and create fine-grain levels of organization. Activity-dependent mechanisms use the patterned activity of receptors as a basis for organization and are best described by the Hebb rule (Hebb 1949). This can be stated in a very simple summary rule: neurons that fire together wire together. What it means is that strongly correlated stimuli tend to reinforce their synaptic connections, whereas synaptic connections supporting uncorrelated stimuli tend to weaken and disappear. Despite its simplicity, this rule is an excellent candidate for building retinotopic maps as the underlying process is looking to connect neurons supporting highly correlated stimuli. There is one other feature that should be considered when investigating activity-dependent mechanisms in sensory organization. That is the signals have to travel along axons and this introduces delays into the interaction between neurons. Hence, even when considering the Hebb rule, we could be considering not instantaneous correlations, but correlations that are delayed in time. There has been some investigation into the effect of delay in neural networks (Coombes and Laing 2009). However, here, we only need to note that delay allows us to construct Hebb like rules, but with correlations delayed in time. 2.3 Biological analogy to organization in artificial sensory systems Let us start by defining an artificial sensory system in simple tractable terms. We have a set of receptors covering a sensory

51

field. Each receptor is denoted by i and has a spatial position in the physical map M P . Every time an object enters the coverage region of receptor i, it will produce an event ei . It is worth noting that both the spatial distribution of the receptors across the physical map M P and the dynamics of the objects moving through the sensory field will determine the correlation between the events ei and e j produced by any two receptors i and j. Let us assume that we only have access to the set of events produced by each receptor and that neither the location of the receptors in the physical map M P nor the detailed dynamics of the objects moving through the sensory field are known. Given these constraints, we face the problem of inferring the physical map M P . It is obvious that we can gain inspiration from biological sensory systems to devise sensory organization methods for artificial systems. Firstly, since we can only gain access to the set of events produced by the sensory system, our organizational method needs to use activity-dependent mechanisms accounting for the correlation between the events produced by every receptor. The correlation between the events produced by each receptor allows us to extract patterns of activity in the sensory system and is used as an indirect measure of their spatial proximity in the physical map M P . Secondly, since we have no a priori information about the approximate location of receptors in the physical map M P , we will devise a dynamic sensory organization method that proceeds hierarchically and builds multilevel descriptions of receptor location. Starting from coarse descriptions of organization in which receptors are clustered together within regions, they reach highly specific descriptions by successive organization refinements. 2.4 Sensory organization and retinotopic modeling There has been a considerable body of work developed in modeling retinotopic processes. It would at first sight seem sensible to adapt to our problem of previous algorithms that have been proposed for retinotopic modeling. However, it turns out that this is not so straightforward. In the problem of organizing artificial sensory systems, there is no a priori knowledge about the location of the receptors nor about the dynamics of the objects moving through the sensory field. For example, we do not know which subset of receptors forms the boundary. Given these constraints, the objective is to construct an algorithm that will efficiently find the topographic organization of any distribution of sensory receptors. By contrast, retinotopic modeling has typically focused on the dynamics of nerve cell development in an attempt to explain the underlying biological process. As such, it has had to consider a family of problems and features not present within an artificial sensory system. Generally the approaches attempt to define a dynamical model for the evolution of the initially disordered maps. For

123

52

example, in Cottrell and Fort (1986), a hebbian-based dynamics is defined for a two-dimensional mapping. However, the approach relies on the application of known boundary conditions. In effect, the receptor positions on the boundaries must be known to find the topology of the remaining receptors. The requirement of the a priori knowledge of the boundaries makes the algorithm unsuitable for the problem of sensory organization as defined in Sect. 2.3. A more flexible approach first proposed in one dimension in Häussler and Malsburg (1983), and later extended to an arbitrary number of dimensions (Güßmann et al. 2007), focused specifically on the dynamics of the nerve cell connectivity development. The resultant model presents problems in the form of periodic boundary conditions and the inclusion of a parameter that must be tuned to find a solution. Changing the boundary condition leads to significant changes in the connectivity dynamics. However, the main limitation is related to the parameter dependence and was demonstrated in Zhu (2008). Here, it was shown that there is a very sensitive dependency on the algorithm parameters. The solution “froze” in errors in the form of disordered domains of correctly ordered, but incorrectly located receptors. It did this even for small system sizes. This latter limitation seems to be common to many algorithms we considered. Small regions of correctly ordered receptors are often created; however, they are frequently large distances from their correct location. Relocating them becomes progressively more difficult for the algorithm as system size increases. The more general field of self organizing maps (Kohonen 1990) is also not directly appropriate to the problem of sensory organization as defined in Sect. 2.3 as they implicitly require metric information of the underlying input vector. 3 Sensory organization method In this section, we describe the sensory organization method inspired in the biological processes that were outlined in Sect. 2. Firstly, we define the strength of the connection between two receptors in terms of the probability that they are activated simultaneously. This is an activity-dependent quantity and serves to indirectly infer the proximity between two receptors. Secondly, we propose a hierarchical clustering method for sensory organization that uses the strength of the connection between receptors. This method allows us to organize sensory systems in multiple levels of accuracy, by defining receptor clusters of different sizes. Finally, we discuss how our method could be used to handle system growth and failure. 3.1 Activity-dependent inference of receptor proximity Each receptor i produces a continuous time measurement. That time measurement can be translated into a discrete set of

123

Biol Cybern (2014) 108:49–60

events, E i (t), which contains all the events that have occurred up until time t. Let us denote by eim the m-th event produced at receptor i and by tim the time that it is produced. We can also apply the sensible restriction that we only consider events over a limited time interval ΔT . Hence, we can define a restricted set of events E i (t, ΔT ) as E i (t, ΔT ) := {eim : (eim ∈ E i (t)) and (t > tim ≥ t − ΔT )}

(1)

The set of events E i (t, ΔT ) constitutes the information that will effectively be produced by any receptor i in the sensory system. Note that we represent the events produced by a receptor by their existence, or not, at a particular time instant, and consequently, any details concerning the temporal structure of the events are lost by its projection into this binary state representation. The main advantage of this representation is that once the binary state operator is defined, we need not worry about any temporal variations in the structure of the events, except for failures or malfunctions in the receptor. There are several steps that we need to consider in order to use E i (t, ΔT ) for sensory organization. First, we need to define what we mean by correlation. If an event occurs on receptor i at time t, we wish to see whether an event also occurred on receptor j between times t − t0 and t − t1 , where t0 > t1 . We need a function that given the time of two events can return a value of one if a correlation exists or zero if one does not. We will start by defining the typical step function,  1, if x ≥ 0 θ (x) = (2) 0, otherwise Let tim and t jn be the times of two events, m and n, at receptors i and j, respectively. Using (2), we can build a function r (tim , t jn ) that will return the value one whenever tim − t0 < t jn < tim − t1 and zero otherwise:     r (tim , t jn ) = θ t jn − (tim − t0 ) θ (tim − t1 ) − t jn . (3) All we now need to do is sum all the terms defined by (3) over all the elements in E i (t, ΔT ) and E j (t, ΔT )   Ri j (t0 , t1 ) = r (tim , t jn ). (4) eim ∈E i e jn ∈E j

Equation (4) sums all the correlations. We can turn it into a frequency measure by dividing it by the time interval that defines the event set, ΔT . This gives us, 1 (5) Ri j (t0 , t1 ). ΔT However, to interpret this equation as a frequency that can then be regarded as a probability, a little care must be taken with the units of time. Hence, an interpretation as a probability is problematic, unless we somehow eliminate ΔT from the equation. p(i, j) =

Biol Cybern (2014) 108:49–60

53

With the correct time scale (5) can be considered to be the joint probability that receptors i and j are activated quasisimultaneously, i.e., within the time interval defined by t0 and t1 . Our path now is to find the probability that a given receptor is activated during the time interval ΔT . We will then use this to construct a conditional probability that does not depend on the definition of ΔT . The probability that receptor i is activated is simply found from the cardinality of the event set E i (t, ΔT ) divided by the time interval ΔT  1   (6) E i (t, ΔT ) p(i) = ΔT where |·| denotes set cardinality. This equation like that of (5) can only be considered a full probability, if the measure of time is correctly defined. We can now construct the conditional probability as Ri j (t0 , t1 ) p(i, j) . = p(i| j) =    p( j)  E j (t, ΔT )

that a cluster of physically close receptors can be formed by grouping the receptors that collectively have large conditional probabilities between them. In addition, we want to define clusters of receptors that are compact, i.e. we want clusters whose content to boundary ratio is low. Thus, we need to provide an operative definition of cluster compactness. Let G denote a target cluster in a set of receptors S and G denote its complementary cluster. The latter is the set of all the receptors in S that do not belong to G. We can define two weight functions, WGG and WGG , that will be useful to describe the quality of the target cluster G, as follows:   1 WGG = p(i| j) (8) |G||G| s ∈G si ∈G

WGG

j

  1 = p(i| j) |G||G|

(9)

si ∈G s j ∈G

(7)

This eliminates the factor of ΔT ; hence, there is no ambiguity and we can regard (7) as a true conditional probability regardless of the definition of time scale. A further point to note about (7) is that it is implicitly assumed that the stimulus generating the events is statistically homogeneous during the time interval ΔT . This ensures that p(i, j) and p(i) have no time dependency during this interval. One interpretation of p(i| j) is that it measures the strength of the connection between receptors i and j. The larger its value, the larger the probability that the two receptors are adjacent, or at least close. Also, it must be emphasized that the meaning of the conditional probability p(i| j) is determined by the choice of t0 and t1 . If t0 > 0 and t1 < 0, p(i| j) defines the probability that receptor i activates simultaneously with receptor j, given that j activates. Simultaneity here simply means both receptors are activated within the same time window. Whenever t0 > 0 and t1 > 0 (or t0 < 0 and t1 < 0), p(i| j) can be interpreted as the probability that receptor i is activated after (before) receptor j activates. The conditional probability p(i| j) is the activity-dependent quantity that will be used to devise the hierarchical clustering method for sensory organization. 3.2 Hierarchical clustering method 3.2.1 Definition of cluster compactness The first step toward defining a hierarchical clustering method for sensory organization is to devise a method for building clusters of receptors that are physically located close to each other. In reference to the Hebb rule, which would strengthen connections that have a high conditional probability and eliminate those that have low values, we can infer

Since WGG is a function of the conditional probabilities p(i| j) connecting cluster receptors to non-cluster receptors, it can be regarded as the boundary of the cluster G. By contrast, WGG is a function of the conditional probabilities p(i| j) connecting cluster members among themselves, and hence, it can be interpreted as the content of the cluster G. Consequently, the ratio between WGG and WGG can be used to define a measure of cluster compactness Comp{G}: Comp{G} =

WGG WGG

(10)

A maximally compact cluster G will thus have a maximal value of Comp{G}. 3.2.2 Receptor set partitioning The notion of cluster compactness introduced in Sect. 3.2.1 can be used to inform the definition of a cluster of neighboring receptors in a sensory system. Given a set of receptors, a number of options can be devised to group receptors into maximally compact clusters. In this study, we propose a method that builds four equally sized clusters (quadrants) that are maximally compact. This clustering problem is reduced to a traditional graph partitioning problem, whose solution can be obtained by minimizing a cost function, the cut, defined over any partition of the graph, where the conditional probability p(i| j) represents the bond strengths between the vertices of the graph. Several cuts have been defined in the literature in the context of data clustering and image segmentation, such as the cuts proposed in Shi and Malik (2000) and Wu and Leahy (1993). In this work, we define a new cut, the Jcut which is based on the notion of cluster compactness. The Jcut will define a relatively simple partitioning algorithm, and the authors believe much more optimal solutions could be found, especially accounting for the specifics of each problem.

123

54

Biol Cybern (2014) 108:49–60

Given a set of receptors S, our algorithm of quadrant partitioning applies two consecutive bi-partitions. A bi-partition of S consists of two artificial clusters A and B, to which all the receptors are assigned. The clusters are exclusive, i.e., a receptor can only be in one or the other. We also place the added restriction that both clusters must be the same size, |A| = |B|. Given a bi-partition of S into the artificial clusters A and B, we define the cut Jcut (A, B) as follows: 1 1 + Jcut (A, B) = Comp{A} Comp{B} WAB WBA = + W WBB AA  si ∈B,s j ∈A p(i| j) si ∈A,s j ∈B p(i| j) =  + si ∈A,s j ∈A p(i| j) si ∈B,s j ∈B p(i| j) (11) Based on the cut Jcut (A, B), we obtain the following bipartition {H0 , H1 } of the set of receptors S: {H0 , H1 } = arg min Jcut (A, B) {A,B}

(12)

By minimizing Jcut , we seek to obtain two equally sized clusters separated by the smallest possible boundary, i.e. two equally sized clusters that are maximally compact. On average, in a maximally compact bi-partition, two receptors within the same cluster are more likely to generate correlated events than two receptors assigned to different clusters. Hence, since the conditional probability between two receptors can be used as an indirect measure of their spatial proximity, a maximally compact bi-partition constitutes a simplified description of how receptors are spatially distributed. The same bi-partitioning process can be applied to H0 and H1 to obtain four equally sized quadrants Q0 , Q1 , Q2 , and Q3 , such that Q0 ∪ Q1 = H0 , Q2 ∪ Q3 = H1 and Qi ∩ Q j = ∅, for 0 ≤ i, j ≤ 3 and i = j. Alternatively, another bi-partition {H 0 , H 1 } orthogonal to {H0 , H1 } can be obtained, from which the four quadrants Q0 , Q1 , Q2 and Q3 can be defined as the intersections between the two halves defining each bi-partition. Finding two clusters H0 and H1 that minimize the cost function Jcut can be a difficult task, since this cut will in general have multiple local minima. We developed an optimization algorithm combining gradient descent and stochastic sampling which provided satisfactory results in a great variety of problems. Nevertheless, the authors believe other optimization algorithms could be used that improve our results. 3.2.3 Quadtree construction The quadrant partitioning algorithm outlined in Sect. 3.2.2 can be applied to the whole set of receptors in the sensory system S as well as to any subset of receptors of S. Specifi-

123

Fig. 2 The transformation can be defined in layers. The first stage is to build a transformation that coordinates the receptors into quadrants, then order the quadrants into the correct relative position. Finally, apply the same process recursively to each new cluster and order the position accounting for the ordination of neighboring clusters

cally, it can be applied to any of the four quadrants in which S can be partitioned and recursively to any new quadrant that is created. As a result, a hierarchical quadtree-type description of the system can be obtained, in which S constitutes the top-level cluster (Fig. 2). For each level of description l, receptors are organized in maps MlR which consist of 4l quadrants, Qlk . Sensory organization will be coarse in maps MlR generated at the top of the quadtree structure (low l), and it will be detailed for maps MlR generated at the bottom of the structure (high l). It is worth noting that quadrant partitioning by itself will only cluster together spatially close receptors, but it will not necessarily arrange quadrants in the correct position relative to each other within each map. Therefore, an additional mechanism to order quadrants at each level is needed. In

Biol Cybern (2014) 108:49–60

order to do so, it is necessary to define a measure of contiguity between clusters. The cut Jcut (A, B) is a good candidate, since it is a function of the connectivity between clusters A and B (through WAB and WBA ) although other measures of contiguity could be used. For the sake of argument, we define S0l = {Ql0 , Ql1 , Ql2 , Ql3 } as the set of the four quadrants defined at level l into which the parent quadrant Ql−1 is partitioned. Consider the values 0 of l = Jcut (Qlk , Ql−1 Pkm m )

where 0 < k < 3, 1 < m < 3. One straightforward interprel is that they measure the connectation of the values of Pkm tivity between the quadrants into which cluster Ql−1 is par0 l−1 l−1 titioned and its neighboring clusters Q1 , Q2 and Ql−1 3 . l−1 in the map M , based Then, if we know the position of Ql−1 0 R l , we can determine the position of each quadrant Ql on Pkm k l = 0, then Ql is not within the map MlR . For example, if P0m 0 l−1 l−1 adjacent to Ql−1 1 , Q2 or Q3 . This uniquely positions the quadrant in the outside corner of Ql−1 0 . Each position in the l ; quadrant structure has a unique pattern of values for Pkm hence, all the quadrants can easily be spatially organized. The described algorithm will function for all levels below the first. Hence, it can be applied recursively from top to bottom, leading to a fully organized quadtree-type structure where top levels provide a coarse description of sensory organization, and bottom levels a fine, detailed description. Figure 2 summarizes the activity-dependent quadtree approach to sensory organization. The nature of the described algorithm is efficient because it measures the strength of the connection between quadrants and regions outside the one being partitioned. The only limitation is that it cannot be applied to the quadrants in the first layer. However, on the assumption that the first layer contains a large number of receptors, a less efficient algorithm comparing the cost function between quadrants in the same layer can be used. Quadrant Q10 can be identified by showing it is adjacent to Q11 and Q12 , but not Q13 , etc.. Add this to the fact that the partitioning process provides partial ordering, the two quadrants found from each half are known to be adjacent, and we can order the first layer. Hence, we can organize the entire structure from the first layer down. 3.3 Application to growth and failure The algorithm, as presented, is used to initiate the organization of the receptor array. However, once the organizational structure is in place it can be used as a real-time tool to rapidly identify the location of newly added receptors. This obviously has applications in terms of growth of the sensory system.

55

Another application of the ongoing use of the algorithm, after the initial ordering is defined, is to check the reliable operation of the existing receptors. If a receptor starts to fail, this will result in the generation of spurious readings, or failure to register actual events. In either case the correlation that the receptor exhibits with its neighbors will be reduced. The constant checking of the correlation relative to its position in the hierarchical structure will allow for the identification of failing receptors. We examine in a later section the effect of mis-locating a receptor and how the algorithm can identify such mislocations.

4 Models of object dynamics The conditional probability p(i| j) defined in Sect. 3.1 captures the patterned activity of the sensory system and constitutes the input data for the organization quadtree algorithm proposed in Sect. 3.2. Its value depends on both the spatial distribution of sensory receptors and the dynamics of the objects moving through the sensory field. Under certain simplifying conditions, it is possible to construct object dynamics, where analytic solutions for p(i| j) can be found. In some cases, these are non-trivial and provide useful models for understanding the system. Therefore, in this section, we will present the asymptotic solution of p(i| j) for three analytic object dynamics, namely infinite plane waves, small diffusing object and gas of diffusing objects. All three analytic models will be used in conjunction with simulated p(i| j) in Sect. 5. 4.1 Infinite plane waves A plane wave, i.e., an excitation wave front that can be described as a straight line, which extends from −∞ to +∞, represents one of the simplest, “realistic” object models. This model leads to a very simple form for p(i| j). Consider two receptors i and j. Any plane wavefront that crosses receptor i will also cross receptor j within a time interval, ±δt, provided the orientation of the wavefront is within a specific range. If this happens, then obviously receptors i and j will activate simultaneously within the time interval ±δt. Let p(i| j) be defined as the probability that receptors i and j are activated simultaneously within the time interval ±δt, given that j is activated. Considering Fig. 3, only wavefronts whose orientation lie between o1 and o2 will provide a contribution to p(i| j). Let θ1 and θ2 be the angles describing the orientation of such wavefronts in a given coordinate system. All we need to know is the number of wavefronts that lie in this range, relative to the total number of wavefronts, and this gives us the value for p(i| j). We can find this number

123

56

Biol Cybern (2014) 108:49–60

Therefore, p(i| j) is finite for cases where i and j are neighbors and zero for all other cases. The value of ai j on a finite receptor field can be found by normalization. In addition, due to the reflecting boundary conditions, the neighboring probabilities will be twice as big, or 4/3 as big, depending on whether the site is a corner site or an edge site, respectively. These boundary effects are dependent on the nature of the dynamics, the choice of boundary conditions and the definition of boundary within the model. 4.3 Gas of diffusing objects Fig. 3 A plane wavefront activates two receptors simultaneously, provided they lie within the angle θi j

by integrating the probability distribution function, P(θ ), for the plane wave orientation between the angles θ1 and θ2 , θ2 P(θ )dθ.

p(i| j) =

(13)

θ1

This integral depends on the form of P(θ ). The simplest assumption to make is that it is a uniform distribution. In this case (13) reduces to 1 p(i| j) = 2π

θ2 dθ = θ1

θi j 1 . (θ2 − θ1 ) = 2π 2π

(14)

where θi j is the angle defined by the coverage field of j, taking the location of i as a reference.

This model constitutes an extension of the small diffusing object dynamics. In this model, we have a collection of independent particles, each obeying the paradigm outlined in Sect. 4.2. After sufficient time, any object will be at any receptor field with the same probability, and therefore, their individual contributions to the activations of each receptor will be the same. By defining p(i| j) in the same sense as in Sect. 4.2, this probability can be found to be:  ai j , for nearest neighbors (16) p(i| j) = ai j , otherwise As opposed to the small diffusing object dynamics, there exist spurious correlations between receptors that are not neighbors, represented by ai j . These spurious correlations arise as a consequence of the existence of multiple particles activating several receptors simultaneously. It is worth noting that although max ai j / max ai j > 1, this quotient approaches one as the number of particles increase. Therefore, our ability to infer spatial proximity between receptors will in general decrease as the number of objects increase.

4.2 Small diffusing object In this model, there exists one single object whose size is smaller than the coverage area of one receptor, which defines the object as a point-like particle. The movement of this particle is described by the probability that it moves from a location in the receptor field to one of its neighboring locations. It is assumed that this probability is memoryless that is the same for every direction and that the boundaries of the receptor field are reflective. This model has a time dependency, but we can assume sufficient time for the transient behavior to decay. Let t¯j be the mean residence time of the object in the coverage field of receptor j. Define p(i| j) as the probability that receptor i activates t¯j ± δt units of time after receptor j is activated. The conditional probability p(i| j) thus defined is given as  ai j , if i and j are nearest neighbors (15) p(i| j) = 0, otherwise

123

5 Results 5.1 Sensory organization We explored the performance of the quadtree organization algorithm proposed in Sect. 3 in a sensory system consisting of 64 receptors arranged in a 8 × 8 square array. This sensory array has an associated quadtree structure consisting of four layers, l = 0, . . . 3, and each layer l in turn consists of 4l groups of 43−l receptors. For example, level 0 consists of one single group of 64 receptors and level 1 consists of four groups of 16 receptors each. The coverage field of each receptor was assumed to be square and non-overlapping, and the sensory system was assumed to completely cover the entire sensory field. We used both the analytical values of p(i| j) obtained in Sect. 4 and simulated ones. The latter were obtained by simulating

Biol Cybern (2014) 108:49–60

57

Table 1 Mis-location errors 12 and 23 for simulated (s) and analytical (a) dynamics Model of object dynamics

log(ΔT )

12

23

Plane wave

→ ∞ (a)

0

0

Small object (N = 1) Gas object (N = 50) Gas object (N = 100)

4 (s)

0

0

→ ∞ (a)

0

0

4 (s)

0.13

0.16

→ ∞ (a)

0

0

4 (s)

0.61

0.48

→ ∞ (a)

0

0

ΔT is the number of time steps in the simulations

the movement through the receptor field of objects obeying the small diffusing object dynamics and the gas of diffusing objects dynamics. We produced sets of activation events E i (t, ΔT ) for time intervals ΔT of different length, and from them, we estimated the conditional probability p(i| j). Therefore, different time intervals ΔT led to different estimations of p(i| j), which converged asymptotically to the analytical expressions for ΔT → ∞. Before feeding the quadtree organization algorithm with p(i| j), receptor labels were randomized to mimic the loss of topological structure in the sensory system as a consequence of the mapping M P → M E . Based on the randomized p(i| j) and assuming that level 1 was perfectly organized, we obtained the quadtree structure by constructing levels 2 and 3. The mis-location error 12 (respectively 23 ) was defined as the total number of receptors incorrectly organized at level 2 (respectively 3), assuming perfect organization at the preceding level 1 (respectively 2). The proposed algorithm achieved an error-free, perfect reconstruction of the spatial organization of receptors when using the analytical values p(i| j) for plane wave dynamics and diffusing object dynamics (Table 1). This was expected trivially for single diffusing object dynamics. Due to the fact that plane wave dynamics and gas-diffusing object dynamics p(i| j) includes spurious, long-range correlations between non-neighboring receptors, the response is less obvious. Despite this, the quadtree strategy constructed correctly the spatial relationship between receptors. When using simulated p(i| j) values, the proposed algorithm showed a performance that depended on the time of simulation ΔT . Figure 4 shows mis-location errors 12 and 23 in a sensory system stimulated by 30 gasdiffusing objects, for increasingly large ΔT . Mis-location errors decreased as the duration of ΔT increased, progressing asymptotically toward zero. Also, the performance of the algorithm degraded as the number of particles increased. Figure 5 shows estimated mis-location errors 12 and 23 versus the number of simulated gas-diffusing objects. The simulation time was intentionally chosen short to allow a bet-

Fig. 4 Estimated mis-location errors 12 and 23 as a function of the duration of the simulation time ΔT . ΔT corresponds to the logarithm of the number of time steps in the simulation

Fig. 5 Estimated mis-location errors 12 and 23 as a function of the number of simulated diffusing objects

ter comparison between the performance of the algorithm at each level. By increasing the number of gas-diffusing objects, the level of spurious correlations increased and accordingly the ability of the proposed algorithm to organize receptors decreased. It is worth noting though that for a small number of particles n ≤ 10 our algorithm reconstructed perfectly the distribution of receptors. For the sake of comparison, Table 1 includes estimated mis-location errors 12 and 23 for ΔT = 104 time steps, and N = 1, 50, 100 diffusing objects. 5.2 Mis-location detection Once a region is ordered, we can use the underlying structure to inform us of the effect of mis-locating of a receptor within this system. This has applications to diagnostic analysis of

123

58

Biol Cybern (2014) 108:49–60

Fig. 6 Weights associated with a receptor at different points within its own region and its neighboring region for plane wave dynamics

the receptor operation as well as in the growth and evolution of the receptor array. For the sake of the argument, let us focus on the level 1 of the quadtree structure. Our aim is to measure the degree of fitness of a given receptor within any of the four quadrants constructed during partitioning. As in the case of the organization algorithm, p(i| j) can be used for defining a measure of receptor fitness within a quadrant. Given a receptor, i, and a quadrant, Qk , let us define the mean weight wi (Qk ), the weight squared deviation σi2 (Qk ) and the maximum weight maxi (Qk ) as: wi (Qk ) =

1  p(i| j), |Qk |

(17)

1  ( p(i| j) − wi (Qk ))2 , |Qk |

(18)

s j ∈Qk

σi2 (Qk ) =

s j ∈Qk

maxi (Qk ) = max( p(i| j))Qk .

(19)

These weighting functions can be used to measure how strongly receptor i fits within quadrant Q. Figure 6 shows the values of wi (Q0 ), σi2 (Q0 ) and maxi (Q0 ) as a function of the true position of the receptor, relative to the top-left quadrant Q0 for a 64 × 64 array of receptors and assuming random incident plane waves. The x-axis is defined such that x = 0 corresponds to a receptor having a position that is not displaced from the central point in the top-left quadrant. All subsequent values are normalized to this first one, to enable us to display all three weights on the same graph. Two paths of receptor positions are considered: a vertical path connecting the central points of the top-left quadrant Q0 and bottom-left quadrant Q1 (this is the path shown in red in the insert of Fig. 6), and a diagonal path connecting the centers of the top-left quadrant Q0 and

123

bottom-right quadrant Q3 (the blue path). Each increment of x displaces the receptor location by one vertical lattice spacing in the case of the red path and one vertical and one horizontal lattice spacing in the case of the blue path. The gray dotted vertical line indicates the position of the boundary between quadrants. The red curves indicate the values for the vertical path and the blue curves correspond to the diagonal path. In the case of maxi (Qk ), only the vertical path is used and it is marked in a slightly darker red to indicate this difference. It is clear that there is a distinct change in the value of the weights as we move across the boundary between quadrants. Both the standard deviation and the maximum weights exhibit a step-like function behavior. The mean value, wi (Q0 ), shows smoother sigmoid-like transition, but again they show a significant difference between receptors that belong to Q0 and receptors that belong to other quadrants. Therefore, based on the values of wi (Q0 ), σi2 (Q0 ) and maxi (Q0 ), we can assess how well receptor i fits within the predefined cluster Q0 . There is an observable difference between the values of the mean and the standard deviation for each of the two path types. However, although there are some slight changes caused by the differences between corners and edges, the majority of this difference is a consequence of the different lengths associated with each step of the √ path. The diagonal path is obviously longer, by a factor of 2, but for convenience, we have plotted it along the same scale. Note, for the case of maximum value, there is very little difference between the results for the two paths, even considering the length scale difference, which is why only one path is shown. As the quadrants become smaller, then the effective sharpness of the transition becomes less pronounced, but even for the 8×8 case, the transition is clear. As we have used random incident plane waves to construct this graph, we have a weakening of the resolution for smaller clusters. This is because plain waves have built-in correlations along the length of the plane wave front (which is from one end of the cluster to the other). Plane waves are excellent for constructing large clusters, the correlations build in longer-range information into the conditional probability matrix, which helps distinguish positions on a larger-length scale. By contrast, the use of a correlation matrix derived from fine detailed stimulus, such as a diffusive gas, provides much better discrimination for very small clusters. However, such correlations are useless for discriminating on a large scale.

6 Discussion and conclusions Biological sensory systems build topographically organized maps for encoding and representing external information. Topographic maps represent information in an orderly and

Biol Cybern (2014) 108:49–60

detailed manner, and this representation is surprisingly constant both across time and across individuals (Kaas 2000). Hence, pre-programmed mechanisms must be in place during development, for organizing sensory projections into topographic maps which mimic the topology of sensory organs. Recent studies have shown that the ability to organize sensory projections is not a feature restricted to development stages, but is also present in mature sensory systems. Specifically, biological sensory systems can be reorganized both structurally and functionally in adulthood in response to severe damages, such as spinal cord injury (Ghosh et al. 2009) or stroke (Nelles et al. 1999; Johnson-Frey 2004). In general, there exists a high consensus among neuroscientists that all levels and types of sensory systems, including the somatosensory, the visual and the auditory systems, are mutable both during development and in adulthood (Kaas 2000). By making an analogy with biological sensory systems, we have developed a robust activity-dependent process for organizing large arrays of receptors. The process proceeds hierarchically, from coarse levels of organization to fine ones, and uses the correlation between detected events to infer the relative locations of receptors. The process can be used to configure a system from unknown initial arrangement, or it can be used as an ongoing process to compensate damage and physical changes in the array structure, or incorporate new receptors to the array. Although the aim of the present work was not developing optimal partitioning algorithms, the obtained results demonstrate that the proposed strategy functions well all the way down to the placement of single receptors. More sophisticated strategies could be devised. In addition, once the orientation structure is found (the receptor topology), it would be feasible to devise a method that use time difference analysis to assign relative separation distances (the receptor metric). It is also worth noting that although the organization algorithm was presented as a hierarchical structure, it could operate on multiple scales simultaneously. The proposed organization algorithm represents a generalized process. It does not depend upon specifically tailored input signals to produce the desired organization. On the contrary, it exploits any sensed data collected in the usual operation mode of the sensory system, analyzing their inherent correlations. This makes it a truly nonparametric approach. It relies upon the detection of the actual stimulus that the receptors are designed for, and therefore, it does not require the addition of systems and extra receptors dedicated to the management of location. Due to its generality, the proposed method can be applied to a variety of artificial systems. For example, it is ideal for use in large-scale sensor networks and opens up the possibility of different approaches to management and maintenance of such systems. It can equally be considered as a control and management system to support the construction of large-

59

scale sensor arrays for artificial skin (Shih et al. 2010) or compound visual detectors (Horisaki et al. 2010; Hornsey et al. 2004). In this context, it is ideal for processing information from sensor arrays that are likely to be deployed in an irregular pattern. Its ability to reconfigure its data model as a consequence of damage, or component failure (or even the case of growth, or component addition), adds to its flexibility and possible applications. Finally, its ability to incorporate new receptors (deal with system growth) without modification makes it an extremely flexible approach.

References Coombes S, Laing C (2009) Delays in activity-based neural networks. Philos Trans R Soc A Math Phys Eng Sci 367:1117–1129. doi:10. 1098/rsta2008.0256 Cottrell M, Fort JC (1986) A stochastic model of retinotopy: a selforganizing process. Biol Cybern 53:405–411 Ghosh A, Sydekum E, Haiss F, Peduzzi S, Zörner B, Schneider R, Baltes C, Rudin M, Weber B, Schwab ME (2009) Functional and anatomical reorganization of the sensory-motor cortex after incomplete spinal cord injury in adult rats. J Neurosci 29(39):12,210–12,219 http:// www.ncbi.nlm.nih.gov/pubmed/19793979 Goodman CS, Shatz CJ (1993) Developmental mechanisms that generate precise patterns of neuronal connectivity. Cell 72(Suppl):77–98 Güßmann M, Pelster A, Wunner G (2007) Synergetic analysis of the Häussler-von der Malsburg equations for manifolds of arbitrary geometry. Annalen der Physik 519:379–394. doi:10.1002/andp. 200610243 Häussler AF, von der Malsburg C (1983) Development of retinotopic projections—an analytical treatment. J Theor Neurobiol 2:47–73 Hebb DO (1949) The organization of behavior: a neuropsychological theory. Wiley, New York Horisaki R, Kagawa K, Nakao Y, Tanida J (2010) Irregular lens arrangement design to improve imaging performance of compound-eye imaging systems. Appl Phys Express 3(2):022–501. doi:10.1143/ APEX.3.022501 Hornsey R, Thomas P, Wong W, Pepic S, Yip K, Krishnasamy R (2004) Electronic compound-eye image sensor: construction and calibration. In: Blouke MM, Sampat N, Motta RJ (eds) Society of photooptical instrumentation engineers (SPIE) conference series, society of photo-optical instrumentation engineers (SPIE) conference series, vol. 5301, pp 13–24 doi:10.1117/12.526811 Johnson-Frey SH (2004) Stimulation through simulation? motor imagery and functional reorganization in hemiplegic stroke patients. Brain and Cognition 55(2):328–331 http://www.ncbi.nlm.nih.gov/ pubmed/15177807 Kaas JH (2000) The reorganization of sensory and motor maps after injury in adult mammals. In: Gazzaniga MS (ed) The new cognitive neurosciences, 2nd edn. The MIT Press, Cambridge, pp 223–236 Kohonen T (1990) The self-organizing map. Proc IEEE 78:1464–1480 Nassi JJ, Callaway EM (2009) Parallel processing strategies of the primate visual system. Nature reviews. Neuroscience 10(5):360–372. doi:10.1038/nrn2619 Nelles G, Spiekermann G, Jueptner M, Leonhardt G, Mueller S, Gerhard H, Diener HC (1999) Reorganization of sensory and motor systems in hemiplegic stroke patients. A positron emission tomography study. Stroke J Cereb Circul 30(8):1510–1516 http://www.ncbi.nlm.nih. gov/pubmed/10436092 Shi J, Malik J (2000) Normalized cuts and image segmentation. IEEE Trans Pattern Anal Mach Intell 22(8):888–905. doi:10.1109/34. 868688

123

60 Shih WP, Tsao LC, Lee CW, Cheng MY, Chang C, Yang YJ, Fan KC (2010) Flexible temperature sensor array based on a graphitepolydimethylsiloxane composite. Sensors 10(4):3597–3610. doi:10. 3390/s100403597 http://www.mdpi.com/1424-8220/10/4/3597/ Wandell B, Dumoulin S, Brewer A (2007) Visual field maps in human cortex. Neuron 56:366–383 Wu Z, Leahy R (1993) An optimal graph theoretic approach to data clustering: theory and its application to image segmentation. IEEE Trans Pattern Anal Mach Intell 15(11):1101–1113. doi:10.1109/34. 244673

123

Biol Cybern (2014) 108:49–60 Zhu J (2008) Synaptic formation rate as a control parameter in a model for the ontogenesis of retinotopy. In: Proceedings of the 18th international conference on Artificial Neural Networks, Part II, ICANN ’08, pp. 462–470. Springer, Berlin, Heidelberg. doi:10. 1007/978-3-540-87559-848