domain: math info info.info-ti

text:
Suvadip Mukherjee

Member, IEEE Catalina Gonzalez-Gomez

Lydia Danglot

Thibault Lagache

Jean-Christophe Olivo-Marin

Generalizing the Statistical Analysis of Objects' Spatial Coupling in Bioimaging

Keywords: spatial statistics, colocalization, Ripley's Kfunction, coupling, fluorescence microscopy, level sets

We introduce a novel paradigm for statistical analysis of spatial signals, applied to colocalization studies in bioimaging. Quantitative assessment of the spatial colocalization of different molecules in microscopy provides important functional information on cellular processes. By reducing objects (molecules or cluster of molecules) to their position (points), existing methods are predominantly restricted to point-based analyses, and scale poorly to applications involving large, asymmetric and complexin-shape objects. We address this issue, and propose a statistical model for shape-based colocalization analysis. Our solution is a generalization of the Ripley's K-function for arbitrary shapes, and provides a method to statistically interpret the coupling of molecules to biological objects which are implicitly represented via level set embedding. Compared to the state-of-the-art, our solution is efficient and generic, and establishes the theoretical basis for an one-fits-all approach to tackle the heterogeneous challenges in robust interpretation of molecular association in biology. The efficacy of our method is established via synthetic simulations, and a practical application is also described to measure the statistical accumulation of synaptic molecules' spots near the cell's boundary of cultured neurons.

I. INTRODUCTION

T HE analysis of the spatial distribution of molecules and organelles in bioimaging remains a gold-standard for understanding cellular processes at the molecular level. Colocalization studies are performed across spatial scales which range from standard fluorescence microscopy, to superresolution imaging (see Fig. 1), and the significant variability in the data makes it difficult to design an approach which is robust, yet generic. Traditional intensity-based methods are often sensitive to the parameters of the imaging system [START_REF] Aaron | Image co-localization-cooccurrence versus correlation[END_REF], and have recently led to the object-based techniques which define colocalization in terms of the spatial proximity of the detected objects (biomolecules) [START_REF] Levet | Sr-tesseler: a method to segment and quantify localizationbased super-resolution microscopy data[END_REF]- [START_REF] Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF]. By reducing an object's shape to a spatial coordinate (e.g. center-of-mass), such methods are suitable for small, rounded objects shown in Fig. 1(a) and (b), but their use is rather limited when the shapes are large and complex (Fig. 1(c)). Such large objects can be present when the objects' size is larger than the PSF in standard fluorescence microscopy (e.g. cell edges, tubular organelles such as mithocondria [START_REF] Pagliuso | A role for septin 2 in drp1-mediated mitochondrial fission[END_REF] or synaptic aggregates [START_REF] Dani | Superresolution imaging of chemical synapses in the brain[END_REF]).

We propose a theory to unify both point-and shape-based colocalization problems in a single mathematical framework. †The authors are affiliated to the Bioimage Analysis Unit, Institut Pasteur, Paris, France. L. Danglot is affiliated to Institut de Psychiatrie et Neurosciences de Paris, Inserm U1266, France. In this novel paradigm of shape-based spatial pattern analysis, objects are represented via a level set of a higher dimensional embedding function [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF], which makes this solution attractive for both point-objects (approximated as circular/spherical spots), as well as complex shapes. This allows us to study a generalized notion of spatial proximity, which we call coupling.

A. Background and motivation

Different imaging modalities are used to understand intermolecular association, ranging from optical coherence tomography in histological samples [START_REF] Fujimoto | Optical biopsy and imaging using optical coherence tomography[END_REF] to fluorescence and electron microscopy [START_REF] Lichtman | Fluorescence microscopy[END_REF], [START_REF] Bai | How cryo-em is revolutionizing structural biology[END_REF] in molecular cell biology. The inherent heterogeneity in the imaging data has led to two predominant groups of techniques for quantitative assessment of the spatial coupling of biomolecules. For example, in fluorescence microscopy the colocalization of bio-molecules labeled with different fluorescent markers is traditionally quantified through their signal overlap or correlation [START_REF] Manders | Measurement of co-localization of objects in dual-colour confocal images[END_REF]- [START_REF] Wang | Spatially adaptive colocalization analysis in dual-color fluorescence microscopy[END_REF]. However, the colocalization indices derived via these methods are sensitive to the microscope's point spread function. Moreover, such techniques are not applicable in super-resolution microscopy [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF], [START_REF] Betzig | Imaging intracellular fluorescent proteins at nanometer resolution[END_REF] where the colocalized molecules are either only partially overlapping, or do not overlap at all [START_REF] Lagache | Statistical analysis of molecule colocalization in bioimaging[END_REF]. In this context, object-based methods have been developed which compute colocalization via spatial analysis of the objects' positions (defined by either the single molecule localization or the center of mass of the molecule spot) using a suitable second-order spatial statistic [START_REF] Lagache | Statistical analysis of molecule colocalization in bioimaging[END_REF], [START_REF] Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF], [START_REF] Rossy | Method for co-cluster analysis in multichannel single-molecule localisation data[END_REF]- [START_REF] Yin | Ultrafast data mining of molecular assemblies in multiplexed high-density super-resolution images[END_REF].

An important factor for robust quantitative assessment is to distinguish between genuine inter-object coupling and molecule association due to chance. In [START_REF] Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF], the authors describe a statistical tool, namely SODA, which uses a Ripley-based Fig. 2. Leveraging level-sets for the analysis of coupling between points and larger objects Our method is designed for analyzing the spatial coupling of points (or small spotty objects, shown in cyan) to larger and complex-in-shape objects (shown in red). We leverage level-sets to embed the complex shape of large, red objects and map the distance of all FOV' s points to embedded contours. We then compute the number of points in several, level-set delimited domains ωmn for m < n, and statistically characterize the expected distribution under the hypothesis of (cyan) points' complete spatial randomness. After having determined the level-set domains where there is a significant accumulation of points (normalized number of points K 0 mn above statistical threshold, highlighted with a red dashed line), our analysis provides a statistical map of coupled points.

statistic [START_REF] Ripley | Statistical inference for spatial processes[END_REF] to compute the significance of the measured spatial coupling between individual objects. However, this technique is also primarily geared towards localization-based analysis in super-resolution microscopy.

B. Our contribution

The limitations of the contemporary practices motivate us to develop a generic solution to analyze the statistical coupling of molecular localizations to objects with arbitrary shapes. A novel colocalization statistic is introduced, which generalizes the Ripley's K-function for any closed shapes that are implicitly represented via level set embedding. This also provides a theoretical generalization of our recent development [START_REF] Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF], and extends its relevance beyond point-point data analysis to embrace a wider gamut of object based colocalization studies in bioimaging. In contrast to the cluster-overlap methods [START_REF] Rossy | Method for co-cluster analysis in multichannel single-molecule localisation data[END_REF]- [START_REF] Yin | Ultrafast data mining of molecular assemblies in multiplexed high-density super-resolution images[END_REF], our solution rejects any stochastic association of points and only reflects the statistically significant coupling to the objects. The mathematical details of the algorithm are described next, followed by experimental evaluations and concluding remarks in Sec. III and Sec. IV respectively.

II. METHOD

To simplify notations, the following theory is presented in two-dimensions, but the model extends naturally to three dimensional data. Let us define the domain of analysis to be Ω ⊂ R 2 . We consider the case of analyzing the co-localization of the elements of a set B to the set A. The set A = {s 1 (p), . . . s n1 (p)} of complex-shaped objects is represented with a collection of curves parameterized by p ∈ [0, 1]. These objects are assumed to be delineated from the input image via appropriate segmentation/object detection techniques [START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF], [START_REF] Mukherjee | Region based segmentation in presence of intensity inhomogeneity using legendre polynomials[END_REF]. B = {u 1 , . . . , u n2 } is a set of points (typically centerof-mass of fluorescent spots) u j = (x, y) ∈ Ω. The set cardinalities are represented as n 1 = |A| and n 2 = |B|. We graphically summarize our method in Fig. 2.

A. Embedding complex-shaped objects with level-sets

Since the shapes in A are defined over Ω, one may define a continuous domain image b : Ω → [0, 1] such that b(x, y) = 1 for all pixels (x, y) which are crossed by or enclosed (for closed contours) by at least one contour in A, and zero otherwise. The binary image b(x, y) may be interpreted as the superposition of the (possibly overlapping) segments defined by the curves in A. We propose an implicit representation of the components in b(x, y) by embedding the object boundaries as the zero level shape of a level set function. We introduce the Lipschitz function φ : Ω → R, which is defined as follows:

φ(x, y) = ≤ 0 ∀(x, y) ∈ ω 0 > 0 ∀(x, y) ∈ Ω \ ω 0 (1)
Here ω 0 ⊂ Ω is the area enclosed by the union of the objects in A, and the object boundaries correspond to the zero level set of φ. In general, for any scalar k, the region enclosed by the klevel set of φ is represented by ω k = {(x, y)|φ(x, y) ≤ k}, and the intermediate zone between two levels m and n (m ≤ n) is defined to be ω mn = ω m ∩ ω n .

B. Quantifying the coupling of points to objects

We define the statistic K mn which is proportional to the number of events u ∈ B which occur inside the region ω mn . Mathematically,

K mn = |Ω| n 2 n2 j=1 χ mn (u j ) where u j ∈ B (2)
The region indicator function χ mn (u) = 1 if u ∈ ω mn and zero otherwise.

Proposition. If the spatial distribution of the events u ∈ B follow a homogeneous Poisson distribution in Ω, K mn asymptotically converges to a Gaussian distribution.

If the events u ∈ B are realizations of a homogeneous Poisson spatial process (a model for complete spatial randomness (CSR) [START_REF] Lagache | Analysis of the spatial organization of molecules with robust statistics[END_REF]), it implies that the spatial spread of the events is uniformly random in the domain Ω. Under CSR, the events in B are i.i.d, and therefore this serves as the null hypothesis for (lack of) coupling between the events in B to the regions in A. Under the null hypothesis, χ mn (u) is a Bernoulli random variable with parameter p mn , and consequently K mn is binomially distributed. This parameter may be computed to be p mn = |ω mn |/|Ω| as each event is equally likely to occur in any closed sub-domain due to the CSR assumption.

Furthermore, K mn converges to the Normal law, expressed as K ∼ N (µ mn , σ mn ) for sufficiently large n 2 . Our objective is to estimate the parameters of this distribution in a closed form under the assumption of a spatial randomness of the points in B. The expected value of K mn under the null hypothesis is calculated to be

µ mn = E [K mn ] = |Ω| n 2 n2 j=1 E [χ mn (u j )] = |Ω| n 2 × n 2 |ω mn |/|Ω| = |ω mn | (3)
The variance of the statistic can be expressed as follows:

σ 2 mn = E (K mn -µ mn ) 2 = E K 2 mn -µ 2 mn = |Ω| 2 n 2 2 n2 i=1 E χ 2 mn (u i ) + n2 i=1 n2 j=1 j =i E [χ mn (u i )χ mn (u j )] -n 2 2 p 2 mn ( 4
)
Since χ mn is a Bernoulli random variable with expected value p mn , E χ 2 mn (u) = p mn as well. Furthermore, the complete spatial randomness of the point process ensures that the location of a point

u i ∈ B is independent of another point u j ∈ B ∀i = j, which yields E [χ mn (u i )χ ω (u j )] = E [χ mn (u i )] E [χ mn (u j )]
. Substituting this result in Eq.4, we compute the variance as follows:

σ 2 mn = |Ω| 2 n 2 2 n 2 p mn + n 2 (n 2 -1)p 2 mn -n 2 2 p 2 mn = |ω mn | (|Ω| -|ω mn |) /n 2 (5)
Using vector notations, we have K = K 0 , . . . , K mn , . . . , M = [µ 0 , . . . , µ mn , . . .] and Σ = diag [σ 0 , . . . , σ mn , . . .] such that K ∼ N (M, Σ).

C. Statistical characterization of the spatial coupling

To statistically characterize the coupling of points to objects from the computed K vector, we first reduce it to

K 0 = Σ -1 . [K -M] , (6)
which converges to a standard gaussian vector under the complete spatial randomness hypothesis:

K 0 ∼ N (0 N , 1 N,N )
where N is the length of the vector K, 0 N a 0-vector of length N and 1 N,N a N by N identity matrix. Components of K 0 are independent random variables and we can therefore compute the probability that the maximum component of K 0 , which we denote sup N K 0 , is greater than arbitrary value x is equal to

Pr sup N K 0 > x = 1 - 0<m<n Pr K 0 mn < x = 1 -cdf N (x), (7)
where cdf(

x) = x -∞ exp -u 2 2 √ 2π
du is the cumulative density function of the normal gaussian law. Therefore, we can com- pute the p-value for rejecting the null hypothesis of points' complete spatial randomness (CSR) with

p-value[CSR] = 1 -cdf N sup N K 0 . ( 8
)
To detect significant components of K 0 where points accumulate around objects, we use the threshold T (N ) = 2 log(N ) which computes the significant components of a N -length vector corrupted with unit-variance, white noise [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF]. Thus, for each component K 0 mn > T (N ) we estimate the number of coupled points C mn between level sets 0 < m < n, i.e. statistically above the expected number of points under complete spatial randomness, with

C mn = χ {K 0 mn >T (N )} n 2 |Ω| (K mn -|ω mn |) (9)
The indicator function χ {K 0 mn >T (N )} = 1 if K 0 mn > T (N ) and zero otherwise. Finally, for each point's position u j ∈ B, for 1 ≤ j ≤ n 2 , we can compute the probability p C (u j ) that this point is coupled to objects of set

A p C (u j ) = |Ω| n 2 0<m<n χ mn (u j ) C mn K mn (10)
The mean coupling distance d C [B → A] of points of set B to objects of set A is therefore given by

d C [B → A] = 1 j p C (u j ) 1≤j≤n2 p C (u j )d [u j → A] (11)
where the distance d [u j → A] of point u j to set A is equal to the minimum distance of u j to any pixel x, y such that the level set function Φ(x, y) = 0:

d [u j → A] = min x,y s.t. Φ(x,y)=0 d [u j → (x, y)].

III. EXPERIMENTS

In this section we evaluate the proposed methodology using both simulation studies and examples in bioimaging. For different levels of coupling α = 0, 0.1 . . . 1.0 and distance µ c = 1 pixel and σ c = 0.3 pixels, we compare the simulated and the measured percentage of coupling for increasing objects' length. Ideally, the measured percentage of coupling should be close to the simulated percentage, which corresponds to the black dashed line in Fig. 3. We also compare the accuracy of our method with Ripley-based analysis SODA [START_REF] Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF]. We observe that, while level-sets and Ripley-based methods produce almost similar results for point objects (L = 1), the efficacy of level-set method is significantly better for more elongated structures.

B. Coupling of synaptic spots to neuronal dendrites

To demonstrate an application in bioimaging, we use widefield fluorescence imaging (Fig. 4) and measure the spatial coupling of pre-synaptic spots (Synapsin, labeled in green) with the dendrite of a post-synaptic neuron labeled with MAP2 (blue). Projecting axons of pre-synaptic neurons are known to be apposed to post-synaptic dendrites, and we expect to find a significant accumulation of Synapsin in close proximity to the cell boundary. We validate this phenomenon quantitatively using our proposed technique. Conventional colocalization methods are not applicable here, as the neuron dendrites are large, asymmetrical objects and the Synapsin spots are physically apposed, but not entirely superposed to the cell shapes.

Quantitative evaluations are performed on a set of 18 different images by designing a protocol in ICY [START_REF] De Chaumont | Icy: an open bioimage informatics platform for extended reproducible research[END_REF]. For spatial analysis, the set of objects (set A) is obtained by segmenting the MAP2 response using k-mean thresholding algorithm in ICY [START_REF] De Chaumont | Icy: an open bioimage informatics platform for extended reproducible research[END_REF]. Centroids of the Synapsin spots, automatically extracted via a wavelet-based spot detector [START_REF] Olivo-Marin | Extraction of spots in biological images using multiscale products[END_REF] constitute the point set B. Over 92% of the detected spots were found to be accumulated within the domains ω mn (0 ≤ m < n ≤ 9) defined by the level sets of the segmented cells, of which more that 89% spots were found to be statistically apposed to the cell shape with non-zero coupling probability (Eq. 10). Statistically significant accumulation of spots (p-value = 10 -178±18 ) were observed at an average coupling distance of 2.1 pixels, with a peak at a distance less than 100nm from cell boundary. Additionally, our method also predicts the uncertainty in the spatial process via the average coupling probability pC = 1 n2 n2 j=1 p C (u j ) = 56.4%. This relatively low value is expected, as the dendrites span a considerable area which increases the risk of observing Synapsin accumulation due to chance. The distance of accumulation of synaptic vesicles are found to be in agreement with the typical size of synaptic buttons [START_REF] Wilhelm | Composition of isolated synaptic boutons reveals the amounts of vesicle trafficking proteins[END_REF], and the quantitative results indeed suggest a strong association between the two markers used in this positive control.

IV. CONCLUSION This paper describes a mathematically rigorous approach to analyze inter-molecular coupling in quantitative bioimaging. The method is robust, generic and addresses the limitations of the state-of-the-art. Experimental validations establish the efficacy of our algorithm over standard Ripley-based techniques, especially for elongated objects. Although this work primarily addresses colocalization studies in bioimaging, the described principles are suitable for generic spatial data analysis with applications in histopathology, geoscience and remote sensing.

Fig. 1 .

 1 Fig. 1. Examples of colocalization studies in bioimaging. A two-channel fluorescent microscopy image is shown in (a), where the fluorescent spots in the green and red channel correspond to the individual protein responses. In (b), 3D single molecule localization of the proteins Synapsin and VGLUT are shown. Finally, (c) illustrates a widefield fluorescent microscopy image of a rodent neuron, where the challenge is to quantitatively analyze the association of the synaptic molecules (green spots) to the neuronal dendrites (in blue).

Fig. 3 .

 3 Fig. 3. Validation with synthetic images a-Using the stochastic algorithm described in sub-section III-A, we generate elongated objects (n 1 = 100) with different lengths, and the corresponding regions ωmn are shown in different colors. b-For each colocalization level α, we compare the predicted coupling with Ripley-based analysis SODA [4] (in cyan) and our method (in red). Dashed black line corresponds to the ideal estimation.

Fig. 4 .

 4 Fig. 4. Measuring the coupling of synaptic molecules to neuronal dendrites with wide-field fluorescence microscopy. a-Imaging of synaptic molecules (Synapsin in green) and neuronal dendrites (MAP2 marker in blue) in rodents' hippocampal neurons (culture). Upper-right corner is zoomed. b-The cell cytoplasm is automatically segmented, and the cell boundary corresponds to zero level-set ω 0 (shown in white). We also show other level-set enclosed domains ωmn, for 0 ≤ m < n ≤ 9 in different colors. Here 1 pixel 100 nanometers. Normalized accumulation of spots are shown for 18 cells.A. Validation with simulationsTo simulate elongated objects with different shapes and length, for each object we generate an open contour s j (p), p ∈ [0; 1]. For each object 1 ≤ j ≤ n 1 , we start by drawing random initial position s j (0) in the FOV Ω. Then, for l = 1..L, with L the length of objects, we compute iteratively the position ass j p L = s j l-1 L re iθ .In our simulations r = 1 and l ∼ U [0; π]. Restricting the range of each θ l to [0; π] ensures that each object in the set A is rather elongated. Finally, the set of objects A = {s 1 (p), . . . s n1 (p)} are embedded as the zero level of a level set function computed via signed distance transform[START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF]. To simulate spatial coupling of a given proportion 0 ≤ α ≤ 1 of n 2 for the points in B to the objects of the set A, we simulate a Thomas process[START_REF] Lagache | Statistical analysis of molecule colocalization in bioimaging[END_REF] to randomly sample (1 -α)n 2 points in Ω. Then, for each remaining αn 2 coupled point u k , one simulated curvilinear position s j l L is selected at random, and the spatial coordinate is computed asu k = s j l L +r c k e iθ c k , where r c k ∼ N (µ c , σ c ) and θ c k ∼ U [0; 2π].For different levels of coupling α = 0, 0.1 . . . 1.0 and distance µ c = 1 pixel and σ c = 0.3 pixels, we compare the simulated and the measured percentage of coupling for increasing objects' length. Ideally, the measured percentage of coupling should be close to the simulated percentage, which corresponds to the black dashed line in Fig.3. We also compare the accuracy of our method with Ripley-based analysis SODA[START_REF] Lagache | Mapping molecular assemblies with fluorescence microscopy and object-based spatial statistics[END_REF]. We observe that, while level-sets and Ripley-based methods produce almost similar results for point objects (L = 1), the efficacy of level-set method is significantly better for more elongated structures.

ACKNOWLEDGMENT

This work was partially supported by grants from the Region IdF DIM ELICIT program, Labex IBEID (ANR-10-LABX-62-IBEID), France-BioImaging infrastructure (ANR-10-INBS-04), and the program PIA INCEPTION (ANR-16-CONV-0005).