A Density-Driven Method for the Placement of Biological Cells Over Two-Dimensional Manifolds

We introduce a graphical method originating from the computer graphics domain that is used for the arbitrary placement of cells over a two-dimensional manifold. Using a bitmap image whose luminance provides cell density, this method guarantees a discrete distribution of the positions of the cells respecting the local density. This method scales to any number of cells, allows one to specify arbitrary enclosing shapes and provides a scalable and versatile alternative to the more classical assumption of a uniform spatial distribution. The method is illustrated on a discrete homogeneous neural field, on the distribution of cones and rods in the retina and on the neural density of a flattened piece of cortex.


Introduction
e spatial localization of neurons in the brain plays a critical role since their connectivity pa erns largely depends on their type and their position relatively to nearby neurons and regions (short-range or/and long-range connections).Interestingly enough, if the neuroscience literature provides many data about the spatial distribution of neurons in di erent areas and species (e.g.[25] about the spatial distribution of neurons in the mouse barrel cortex, [20] about the neuron spatial distribution and morphology in the human cortex, [4] about the spatial distribution of neurons innervated by chandelier cells), the computational literature exploiting such data is rather scarce and the spatial localization is hardly taken into account in most neural network models (be it computational, cognitive or machine learning models).One reason may be the inherent di culty in describing the precise topography of a population such that most of the time, only the overall topology is described in term of layers, structures or groups with their associated connectivity pa erns (one to one, one to all, receptive elds, etc.).One can also argue that such precise localization is not necessary because for some model, it is not relevant (machine learning) while for some others, it may be subsumed into the notion of cell assemblies [13] that represent the spatiotemporal structure of a group of neurons wired and acting together.Considering cell assemblies as the basic computational unit, one can consider there is actually few or no interaction between assemblies of the same group and consequently, their spatial position is not relevant.However, if cell assemblies allows to greatly simplify models, it also brings implicit limitations whose some have been highlighted in [23].To overcome such limitations, we think the spatial localization of neurons is an important criterion worth to be studied because it could induces original connectivity schemes from which new computational properties can be derived as it is illustrated on gure 2.
Figure 1: Stippling.According to Wikipedia2 , Stippling is the creation of a pa ern simulating varying degrees of solidity or shading by using small dots.Such a pa ern may occur in nature and these e ects are frequently emulated by artists.e pair of boots (le part) have been rst converted into a gray-level image and processed into a stippling gure (right part) using the weighted Voronoi stippling technique by [29] and replicated in [27].Image from [27] (CC-BY license).
However, before studying the in uence of the spatial localisation of neurons, it is necessary to design rst a method for the arbitrary placement of neurons. is article introduces a graphical, scalable and intuitive method for the placement of neurons (or any other type of cells actually) over a twodimensional manifold and provides as well the necessary information to connect neurons together using either an automatic mapping or a user-de ned function.
is graphical method is based on a stippling techniques originating from the computer graphics domain for non-photorealistic rendering as illustrated on gure 1. e lower the density, the shorter the path and the higher the density, the longer the path.On the far le , the shortest path from top to bo om is only 6 connections while this size triples on the far right to reach 19 connections.Said di erently, the le part is the fast pathway while the right part is the slow pathway relatively to some input data that would feed the architecture from the top.C. Due to the asymmetry of cells position, a signal entering on the top side (materialized with small arrows) travels at di erent speeds and will consequently reach the bo om side at di erent times.is represents a spatialization of time.D. Due to the asymmetry of cells position, a signal entering on the le side (materialized with small arrows) slows down while traveling before reaching the right side. is represents a compression of time and may serve as a short-term working memory.

Methods
Blue noise [31] is an even, isotropic yet unstructured distribution of points [22] and has minimal low frequency components and no concentrated spikes in the power spectrum energy [35].Said di erently, blue noise (in the spatial domain) is a type of noise with intuitively good properties: points are evenly spread without visible structure (see gure 3 for the comparison of a uniform distribution and a blue noise distribution).
is kind of noise has been extensively studied in the computer graphic domain and image processing because it can be used for object distribution, sampling, printing, half-toning, etc.One speci c type of spatial blue noise is the Poisson disc distribution that is a 2D uniform point distribution in which all points are separated from each other by a minimum radius (see right part of gure 3).Several methods have been proposed for the generation of such noise, from the best in quality (dart throwing [8]) to faster ones (rejection sampling [6]), see [18] for a review.An interesting variant of the Poisson disk distribution is a non isotropic distribution where local variations follow a given density function as illustrated on gure 1 where the density function has been speci ed using the image gray levels.On the stippling image on the right, darker areas have a high concentration of dots (e.g.boots sole) while lighter areas such as the background display a sparse number of dots.ere exist several techniques for computing such stippling density-driven pa ern (optimal transport [22], variational approach [7], least squares quantization [19], etc.) but the one by [29] is probably the most straightforward and simple and has been recently replicated in [27].

Distribution
e desired distribution is given through a bitmap RGBA image that provides two types of information.e three color channels indicates the identity of a cell (using a simple formula of the type identity = 256 × 256 × R + 256 × G + B for 0 ≤ R, G, B < 256) and the alpha channel indicates the desired local density.is input bitmap has rst to be resized (without interpolation) such that the mean pixel area of a Voronoi cell is 500 pixels.For example, if we want a nal number of 1000 cells, the input image needs to be resized such that it contains at least 500x1000 pixels.For computing the weighted centroid, we apply the de nition proposed in [29] over the discrete representation of the domain and use a LLoyd relaxation scheme.
A ρ(x) More precisely, each Voronoi cell is rasterized (as a set of pixels) and the centroid is computed (using the optimization proposed by the author that allow to avoid to compute the integrals over the whole set of pixels composing the Voronoi cell).As noted by the author, the precision of the method is directly related to the size of the Voronoi cell.Consequently, if the original density image is too small relatively to the number of cells, there might be quality issues.We use a xed number of iterations (n = 50) instead of using the di erence in the standard deviation of the area of the Voronoi regions as proposed in the original paper.Last, we added a threshold parameter that allows to perform a pre-processing of the density image: any pixel with an alpha level above the threshold is set to the threshold value before normalizing the alpha channel.Figure 4 shows the distribution of four populations with respective size 1000, 2500, 5000 and 10000 cells, using the same linear gradient as input.It is remarkable to see that the local density is approximately independent of the total number of cells.

Connection
Most computational models need to de ne the connectivity between the di erent populations that compose the model.
is can be done by specifying projections between a source population and a target population.Such projections correspond to the axon of the source neuron making a synaptic contact with the dendritic tree of the target neuron.In order to de ne the overall model connectivity, one can specify each individual projection if the model is small enough (a few neurons).However, for larger models (hundreds, thousands or millions of neurons), this individual speci cation would be too cumbersome and would hide any structure in the connectivity scheme.Instead, one can use generic connectivity description [11] such as one-to-one, one-to-all, convergent, divergent, receptive elds, convolutional, etc.For such connectivity scheme to be enforced, it requires either a well structured populations (e.g.grid) or a simple enclosing topology [26] such as a rectangle or a disc.In the case of arbitrary shapes as shown on gure 5, these methods cannot be used directly.However, we can use an indirect mapping from a reference shape such as the unit disc and take advantage of the Riemann mapping theorem that states (de nition from [5]): Riemann mapping theorem (from [5]).Let Ω be a (non empty) simply connected region in the complex plane that is not the entire plane.en, for any z 0 ∈ Ω, there exists a bianalytic (i.e.biholomorphic) map f from Ω to the unit disc such that f (z0) = 0 and f (z0) > 0.
Such mapping is conformal, that it, it preserves angles while isometric mapping preserves lengths (developable surfaces) and equiareal mapping preserves areas.Kerzman and Trummer [17] introduced a method to compute the Riemann mapping function using the Szegö kernel that is numerically stable while Trefethen [30] introduced numerical methods for solving the more speci c conformal Schwarz-Christo el transformation (conformal transformation of the upper half-plane onto the interior of a simple polygon).Furthermore, a Matlab toolkit is available in [12] as well as a Python translation (https://github.com/AndrewWalker/cmtoolkit)that has been used to produce the gure 5 that shows some examples of arbitrary shapes and the automatic mapping of the polar and Cartesian domains.However, even if automatic, this mapping can be perceived as not intuitive.Provided the shape are not too distorted, we'll see in the results section that ad-hoc mapping can also be used.

Visualization
Having now a precise localization for each cell of each population, we have several ways of visualizing the activity within the model.e most straightforward way is to simply draw the activity of a cell at its position using a disc of varying color (a.k.a.colormap) or varying size, correlated with cell activity.is requires the total number of cells to be not too large or the display would be clu ered.For a moderate 1000 cells

Results
We'll now illustrate the use of the proposed method on three di erent cases.

Case 1: Retina cells
e human retina counts two main types of photoreceptors, namely rods, and cones (L-cones, M-cones and S-cones).
ey are distributed over the retinal surface in an non uniform way, with a high concentration of cones (L-cones and M-cones) in the foveal region while the rods are to be found mostly in the peripheral region with a peak density at around 18-20 • of foveal eccentricity.Furthermore, the respective size of those cells is di erent, rods being much smaller than cones.e distribution of rods and cells in the human retina has been extensively studied in the literature and is described precisely in a number of work [10,1].Our goal here is not to t the precise distribution of cones and rods but rather to give a generic procedure that can be eventually used to t those gures, for a speci c region of the retina or the whole retina.e main di culty is the presence of two types of cells having di erent sizes.Even though there exist blue-noise sampling procedures taking di erent size into account [35], we'll use instead the aforementioned method using a two stage procedure as illustrated on gure 6.
A rst radial density map is created for the placement of 25 cones and the stippling procedure is applied for 15 steps to get the nal positions of the 25 cones.A linear rod density map is created where discs of varying (random) sizes of null density are created at the position of the cones.ese discs will prevent the rods to spread over these areas.Finally, the stippling procedure is applied a second time over the built density map for 25 iterations.e nal result can be seen on gure 6C where rods are tightly packed on the le , loosely packed on the le and nicely circumvent the cones.A. e density map for cones placement (n=25) is a circular and quadratic gradient with highest density in the center.B. e density map for rods placement (n=2500) is built using the rods distribution.Starting from a linear density, "holes" with di erent sized are created at the location of each cone, preventing rods to spread over these areas during the stippling procedure.C. Final distribution of cones and rods.Cones are represented as white blobs (splines) while rods are represented as Voronoi regions using random colors to be er highlight the covered area.

Case 2: Neural eld
Neural elds describe the dynamics of a large population of neurons by taking the continuum limit in space, using coarse-grained properties of single neurons to describe the activity [34,33,2,9].In this example, we consider a neural eld with activity u that is governed by an equation of the type: w(x, y)f (u(y, t))dy + I(x) + h e lateral connection kernel w is a di erence of Gaussian (DoG) with short range excitation and long range inhibition and the input I(x) is constant and noisy.In order to solve the neural eld equation, the spatial domain has been discretized into 40×40 cells and the temporal resolution has been set to 10ms.
On gure 7A, one can see the characteristic Turing pa erns that have formed within the eld.e number and size of clusters depends on the lateral connection kernel.Figure 7B shows the discretized and homogeneous version of the DNF where each cell has been assigned a position on the eld, the connection kernel function and the parameters being the same as in the continuous version.e result of the simulation shown on gure 7B is the histogram of cell activities using 40 × 40 regular bins.One can see the formation of the Turing pa erns that are similar to the continuous version.On gure 7C however, the position of the cells have been changed (using the proposed stippling method) such that there is a torus of higher density. is is the only di erence with the previous model.While the output can still be considered to be Turing pa erns, one can see clearly that the activity clusters are precisely localized onto the higher density regions.Said di erently, the functional properties of the eld have been modi ed by a mere change in the structure.is tends to suggest that the homogeneous condition of neural elds (that is the standard hypothesis in most works because it facilitates the mathematical study) is actually quite a strong limitation that constrains the functional properties of the eld.

Case 3: Basal ganglia
e basal ganglia is a group of sub-cortical nuclei (striatum, globus pallidus, subthamalic nucleus, subtantia nigra) associated with several functions such as motor control, action selection and decision making.ere exists a functional dissociation of the ventral and the dorsal part of the striatum (caudate, putamen and nucleus accumbens) that is believed to play an important role in decision making [24,3,21] since these two regions do not receive input from the same structures.For a number of models, this functional dissociation results in the dissociation of the striatum into two distinct neural groups even though such anatomical dissociation does not exist per se (see [14]).Without any proper topography of the striatal nucleus, it is probably the most straightforward way to proceed.However, if each group would possess its own topography, it would become possible to distinguish the ventral from the dorsal part of the BG, as illustrated on gure 8 on a coronal view of the BG.We do not pretend this simpli ed view is su cient to give account on all the intricate connections between the di erent nuclei composing the basal ganglia, but it might nonetheless help to have be er understanding of the structure because it becomes possible to link external input to speci c part of this or that structure (eg.ventral or dorsal part of the striatum).is could lead to di erential processing in di erent part of the striatum and may reconcile di erent theories regarding the role of the ventral and the dorsal part.

Discussion
We've introduced a graphical, scalable and intuitive method for the placement and the connection of biological cells and we illustrated its use on three use-cases.We believe this method, even if simple and obvious, might be worth to be considered in the design of a new class of model, in between symbolic model and realistic model.Our intuition is that such topography may be an important aspect that needs to be taken into account and studied in order for the model to bene t from structural functionality.Furthermore, the proposed speci cation of the architecture as an SVG le associated with the scalability of the method could guarantee to some extent the scalability of the properties of the model.
Notes: All gures were produced using the Python scienti c stack, namely, SciPy [16], Matplotlib [15] and NumPy [32].All sources are available on GitHub  A. Scalable Vector Graphic (SVG) source le de ning each structure in terms of border (solid black lines), major and minor axis (dashed lines), input (red line) and output (blue line).Local density is given by the alpha channel and structure identity is given by the color.In this coronal view of the basal ganglia, the Caudate is red (RGB=(0.83,0.15,0.15)), the GPe is blue (RGB=(0.12,0.46,0.70))and the GPi is green (RGB=(0.17,0.62,0.17

Figure 2 :
Figure 2: In uence of spatial distribution on signal propagation.A. A k-nearest neighbours (k=5) connectivity pa ern shows mid-range connection lengths in low local density areas (le part) and short-range connection lengths in high density areas (right part).B. Shortest path from top to bo om using a k-nearest neighbours connectivity pa ern (k=5).elower the density, the shorter the path and the higher the density, the longer the path.On the far le , the shortest path from top to bo om is only 6 connections while this size triples on the far right to reach 19 connections.Said di erently, the le part is the fast pathway while the right part is the slow pathway relatively to some input data that would feed the architecture from the top.C. Due to the asymmetry of cells position, a signal entering on the top side (materialized with small arrows) travels at di erent speeds and will consequently reach the bo om side at di erent times.is represents a spatialization of time.D. Due to the asymmetry of cells position, a signal entering on the le side (materialized with small arrows) slows down while traveling before reaching the right side. is represents a compression of time and may serve as a short-term working memory.

Figure 3 :
Figure 3: Centroidal Voronoi Tesselation.A. Voronoi diagram of a uniform distribution (n=256) where black dots represent the uniform distribution and white circles represent the centroids of each Voronoi cells.B. Centroidal Voronoi diagram where the point distribution matches the centroid distribution.

Figure 4 :
Figure4: Non-uniform distribution (linear gradient).Di erent population distribution (size of 1000, 2500, 5000 and 10000 cells) using the same linear gradient as input have been computed.Each distribution has been split into four equal areas and the respective proportion and number of cells present in the area is indicated at the bo om of the area.e proportion of cells present in each areas is approximately independent (±2.5%) of the overall number of cells.

Figure 5 :
Figure 5: Conformal mappings.Examples of conformal mappings on arbitrary spline shapes using the conformal Riemann mapping via the Szegö kernel [17].Top line shows conformal mapping of the polar domain, bo om line show conformal mapping of the Cartesian domain.

Figure 6 :
Figure 6: Cones and rods distribution.A. e density map for cones placement (n=25) is a circular and quadratic gradient with highest density in the center.B. e density map for rods placement (n=2500) is built using the rods distribution.Starting from a linear density, "holes" with di erent sized are created at the location of each cone, preventing rods to spread over these areas during the stippling procedure.C. Final distribution of cones and rods.Cones are represented as white blobs (splines) while rods are represented as Voronoi regions using random colors to be er highlight the covered area.

Figure 7 :
Figure 7: Non-homogeneous discrete neural eld. A. Turing pa erns resulting from a continuous and homogeneous neural eld with constant and noisy input.B. Turing pa erns resulting from a discrete and homogeneous neural eld with constant and noisy input.White dots indicate the position of the cells.Mean activity is compute from the histogram of cells activity using 40 × 40 bins.C. Localized Turing pa erns resulting from a discrete and non-homogeneous neural eld with constant and noisy input.White dots indicate the position of the cells.Mean activity is computed from the histogram of cells activity using 40 × 40 bins.

Figure 8 :
Figure 8: Coronal view of the basal ganglia.A. Scalable Vector Graphic (SVG) source le de ning each structure in terms of border (solid black lines), major and minor axis (dashed lines), input (red line) and output (blue line).Local density is given by the alpha channel and structure identity is given by the color.In this coronal view of the basal ganglia, the Caudate is red (RGB=(0.83,0.15,0.15)), the GPe is blue (RGB=(0.12,0.46,0.70))and the GPi is green (RGB=(0.17,0.62,0.17)).B. Distribution of 2500 neurons respecting the local density and structural organization (Caudate: 1345 cells, GPe: 884 cells, GPi: 271 cells).Neurons receiving input are drawn in red, neurons sending output are drawn in blue.Each neuron possesses two set of coordinates: one global Cartesian coordinate set and a local curvilinear coordinates set de ned as the distances to the major and the minor axis of the structure the neuron belongs to.C. Mean activity histogram of the di erent structures using 32x32 bins and a bi-cubic interpolation lter.Each bin includes from zero to several neurons.D. Cell activities represented using the dual Voronoi diagram of the cell position.Each Voronoi region is painted according to the activity of the corresponding centroid (i.e.neuron).
Figure 8: Coronal view of the basal ganglia.A. Scalable Vector Graphic (SVG) source le de ning each structure in terms of border (solid black lines), major and minor axis (dashed lines), input (red line) and output (blue line).Local density is given by the alpha channel and structure identity is given by the color.In this coronal view of the basal ganglia, the Caudate is red (RGB=(0.83,0.15,0.15)), the GPe is blue (RGB=(0.12,0.46,0.70))and the GPi is green (RGB=(0.17,0.62,0.17)).B. Distribution of 2500 neurons respecting the local density and structural organization (Caudate: 1345 cells, GPe: 884 cells, GPi: 271 cells).Neurons receiving input are drawn in red, neurons sending output are drawn in blue.Each neuron possesses two set of coordinates: one global Cartesian coordinate set and a local curvilinear coordinates set de ned as the distances to the major and the minor axis of the structure the neuron belongs to.C. Mean activity histogram of the di erent structures using 32x32 bins and a bi-cubic interpolation lter.Each bin includes from zero to several neurons.D. Cell activities represented using the dual Voronoi diagram of the cell position.Each Voronoi region is painted according to the activity of the corresponding centroid (i.e.neuron).