METHOD FOR DETECTING SPATIAL COUPLING

Abstract

Method for detecting spatial coupling comprising the steps of: a. providing a set of data, b. identifying and segmenting a first and a second sets of objects of interest, wherein the objects of the second set are assimilated to punctual objects, c. determining, using a level set function, an expected number of objects of the second set present within a specified range of distances to at least one given object of the first set in case there were no interactions between said at least one given object of the first set and the objects of the second set, d. determining, using a level set function, an actual number of objects of the second set within the same range of distances to the at least one given object of the first set, and e. comparing said expected amount and said determined amount.

Claims

1. Method for detecting spatial coupling comprising the steps of: a. providing a set of data, b. identifying and segmenting a first and a second sets of objects of interest, wherein the objects of the second set are assimilated to punctual objects, c. determining, using a level set function, an expected number of objects of the second set present within a specified range of distances to at least one given object of the first set in case there were no interactions between said given object of the first set and the objects of the second set, d. determining, using a level set function, an actual number of objects of the second set within the same range of distances to the at least one given objects of the first set, and e. comparing said expected amount and said determined amount.

2. The method according to claim 1 wherein the set of data is a set of imaging data.

3. The method according to claim 1 wherein the spatial coupling to be identified is a molecular spatial coupling between molecular objects of interest.

4. The method according to claim 3 wherein the provided set of data includes a biological tissues imaging data, preferably medical imaging data.

5. The method according to claim 3 wherein the provided set of data includes microscopy imaging data.

6. The method according to claim 1 wherein steps c and d are performed for several objects of the first set.

7. The method according to claim 6 wherein steps c and d are repeated for every objects of the first set.

8. The method according to claim 1 wherein steps c and d are performed for several given ranges of values.

9. The method according to claim 8 wherein the several given ranges of values do not overlap.

10. The method according to claim 8 wherein the several ranges of values are not consecutive.

11. The method according to claim 1 wherein the data comprise spatial coordinates, preferably bidimensional or tridimensional spatial coordinates, and/or temporal coordinates.

12. The method according to claim 1 wherein the first and second sets of objects of interest comprise punctual objects and/or regions.

13. Computer program product comprising code configured to, when executed by a processor or an electronic control unit, perform the method according to claim 1.

Description

[0050] The invention can be better understood at the reading of the detailed examples below, which constitute non-limitative embodiments of the present invention and at the examining of the annexed drawing, on which:

[0051] FIG. 1 is a schematic representation illustrating how to leverage level-sets for the analysis of coupling between points and larger objects;

[0052] FIG. 2 is a comparison of the results of the method according to the present invention with method of the prior art depending on the shape of the objects considered;

[0053] FIG. 3 shows measurements of the coupling of synaptic molecules to neuron's membrane with wide-field fluorescence microscopy using the method according to the present invention, and

[0054] FIGS. 4 to 6 show the result of three different applications of the method according to the present invention on multicolored images,

[0055] FIG. 7 shows an application of the method to an epidemiological application.

[0056] The annexed drawing includes meaningful colors. Although the present application is to be published in black and white, a colored version of the annexed drawing was filed before the EPO.

[0057] Embedding Complex-Shaped Objects with Level-Sets

[0058] Let us define the domain of analysis to be Ω⊂.sup.2.

[0059] We consider the case of analyzing the co-localization of the elements of a set B to the set A. The set A={s.sub.1(p), . . . , s.sub.n1(p)} of complex-shaped objects is represented with a collection of curves parameterized by p comprised in [0;1]. B={u.sub.1, . . . , u.sub.n2} is a set of points u.sub.j=(x; y) in Ω.

[0060] The set cardinalities are represented as n.sub.1=|A| and n.sub.2=|B|. The method is graphically summarized in FIG. 1.

[0061] The method is designed for analyzing the spatial coupling of points (or small spotty objects, shown in cyan on FIG. 1) 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 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.sup.0.sub.mn above statistical threshold, highlighted with a red dashed line), our analysis provides a statistical map of coupled points.

[0062] Since the shapes in A are defined over Ω, one may define a continuous domain image b: Ω.fwdarw.[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 by embedding the object boundaries as the zero level shape of a level set function. We introduce the Lipschitz function Φ: Ω.fwdarw.R, which is defined as follows:

[00003]$\begin{array}{cc}\varphi \left(x,y\right)=\{\begin{array}{c}\le 0\forall (x,y\}\in \omega \\ >0\forall (x,y\}\in \Omega /\omega \end{array}& \left(1\right)\end{array}$

[0063] Here ω represents the contours or the interior zone enclosed by the union of the closed curves in A. The advantage of embedding the shapes via the level set function is that it allows an implicit representation of the curves as the zero level set Φ. Therefore, the region enclosed by the m.sup.th and n.sup.th level set of Φ may be represented as ω.sub.mn={(x,y)|x∈H(n−ϕ)−H(m−ϕ)}, where n>m≥0 and H(Φ) is the standard Heaviside function.

[0064] Quantifying the Coupling of Points to Objects

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

[00004]$\begin{array}{cc}{K}_{\mathrm{mn}}=\frac{.Math.\Omega .Math.}{n_2}\underset{j=1}{\overset{n_2}{.Math.}}{\chi }_{\mathrm{mn}}\left(u_j\right)\mathrm{where}{u}_j\in B.& \left(2\right)\end{array}$

[0066] The region indicator function χ.sub.mn(u)=1 if u belongs to ω.sub.mn and zero

[0067] otherwise. As χ.sub.mn(u) is a Bernoulli random variable with parameter p.sub.mn, K.sub.mn obeys the Binomial distribution. When the spatial points in B follow a homogeneous Poisson distribution and exhibit complete spatial randomness, this parameter may be computed to be p.sub.mn=|ω.sub.mn|/|Ω|.

[0068] Furthermore, K converges to the Normal law, expressed as K˜N (μ.sub.mn, σ.sub.mn) for sufficiently large n2. 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. In a stochastic setup, we assume K to be a random variable with an associated distribution function. Therefore, the expected value of this random variable may be computed as

[00005]$\begin{array}{cc}\begin{array}{c}{\mu }_{\mathrm{mn}}=𝔼\left[K_{\mathrm{mn}}\right]=\frac{.Math.\Omega .Math.}{n_2}\underset{j=1}{\stackrel{n_2}{.Math.}}𝔼\left[{\chi }_{\mathrm{mn}}\left(u_j\right)\right]\\ =\frac{.Math.\Omega .Math.}{n_2}\times n_2\left(.Math.{\omega }_{\mathrm{mn}}.Math./.Math.\Omega .Math.\right)=.Math.{\omega }_{\mathrm{mn}}.Math.\end{array}& \left(3\right)\end{array}$

[0069] The variance of the statistic can be expressed as follows:

[00006]$\begin{array}{cc}\begin{array}{c}{\sigma }_{\mathrm{mn}}^2=𝔼\left[{\left(K_{\mathrm{mn}}-{\mu }_{\mathrm{mn}}\right)}^2\right]=𝔼\left[K_{\mathrm{mn}}^2\right]-{\mu }_{\mathrm{mn}}^2\\ =\frac{{.Math.\Omega .Math.}^2}{n_2^2}(\underset{i=1}{\overset{n_2}{.Math.}}𝔼\left[{\chi }_{\mathrm{mn}}^2\left(u_i\right)\right]+\\ \underset{i=1}{\overset{n_2}{.Math.}}\underset{\begin{array}{c}j=1\\ j\ne i\end{array}}{\overset{n_2}{.Math.}}𝔼\left[{\chi }_{\mathrm{mn}}\left(u_i\right){\chi }_{\mathrm{mn}}\left(u_j\right)\right]-n_2^2{p}_{\mathrm{mn}}^2)\end{array}& \left(4\right)\end{array}$

Since χ.sub.mn is a Bernoulli random variable with expected value p.sub.mn, E[χ.sup.2.sub.mn(u)]=p.sub.mn as well. Furthermore, the complete spatial randomness of the point process ensures that the location of a point u.sub.i of B is independent of another point u.sub.j of B for every i≠j, which yields E[χ.sub.mn(ui)χ.sub.ω(uj)]=E[χ.sub.mn(ui)]E[χ.sub.mn(uj)].

[0070] Substituting this result in Eq. 4, we compute the variance as follows:

[00007]$\begin{array}{cc}\begin{array}{c}{\sigma }_{\mathrm{mn}}^2=\frac{{.Math.\Omega .Math.}^2}{n_2^2}\left(n_2{p}_{mn}+n_2\left(n_2-1\right)p_{\mathrm{mn}}^2-n_2^3{p}_{\mathrm{mn}}^2\right)\\ =.Math.{\omega }_{\mathrm{mn}}.Math.\left(.Math.\Omega .Math.-.Math.{\omega }_{\mathrm{mn}}.Math.\right)/n_2\end{array}& \left(5\right)\end{array}$

Using vector notations, we have K=[K.sup.0, . . . , K.sub.mn, . . . ]′, M=[μ.sub.0, . . . , μ.sub.mn, . . . ]′ and Σ=diag([σ.sub.0, . . . , σ.sub.mn, . . . ]′) such that K˜N(M,Σ).

[0071] Statistical Characterization of the Spatial Coupling

[0072] To statistically characterize the coupling of points to objects from

[0073] the computed K vector, we first reduce it to

K.sup.0=Σ.sup.−1.Math.[K−M],  (6)

[0074] which converges to a standard gaussian vector under the complete spatial randomness hypothesis: K.sup.0˜N(0.sub.N,1.sub.N,N) where N is the length of the vector K, 0.sub.N a 0-vector of length N and 1.sub.N,N a N by N identity matrix. Components of K.sup.0 are independent random variables and we can therefore compute the probability that the maximum component of K.sup.0, which we denote sup.sub.N[K.sup.0], is greater than arbitrary value x is equal to

[00008]$\begin{array}{cc}Pr\left\{\underset{N}{\sup }\left[K^0\right]>x\right\}=1-\underset{0<m<n}{.Math.}Pr\left\{K_{\mathrm{mn}}^0<x\right\}=1-{\mathrm{cdf}}^N\left(x\right).& \left(7\right)\end{array}$

[0075] where cdf(x) is the cumulative density function of the normal gaussian law. Therefore, we can compute the p-value for rejecting the null hypothesis of points' complete spatial randomness (CSR) with

[00009]$\begin{array}{cc}p-\mathrm{value}\left[\mathrm{CSR}\right]=1-{\mathrm{cdf}}^N\left(\underset{N}{\sup }\left[K^0\right]\right).& \left(8\right)\end{array}$

[0076] The null-hypothesis framework ensures interpretability of the result. The fact that there are no tunable parameters and that no training data is required allows the process to operate fully unsupervised. To detect significant components of K.sup.0 where points accumulate around objects, we use the universal threshold T(N)=√(2 log(N)) introduced by Donoho and Johnston, and used for computing the significant components of a N-length vector corrupted with unit-variance, white noise. Thus, for each component K.sup.0.sub.mn>T(N) we estimate the number of coupled points C.sub.mn between level-sets 0<m<n, i.e. statistically above the expected number of points under complete spatial randomness, with

[00010]$\begin{array}{cc}{C}_{\mathrm{mn}}={\chi }_{\left\{K_{\mathrm{mn}}^0>T\left(N\right)\right\}}\left[\frac{n_2}{.Math.\Omega .Math.}\left(K_{\mathrm{mn}}.Math.{\omega }_{\mathrm{mn}}.Math.\right)\right]& \left(9\right)\end{array}$

[0077] The indicator function χ=1 if K.sup.0.sub.mn>T(N) and zero otherwise. Finally, for each point's position uj of B, for 1≤j≤n2, we can compute the probability pc(uj) that this point is coupled to objects of set A

[00011]$\begin{array}{cc}{p}_C\left(u_j\right)=\underset{0<m<n}{.Math.}{\chi }_{\mathrm{mn}}\left(u_j\right)\frac{C_{\mathrm{mn}}}{K_{\mathrm{mn}}}.& \left(10\right)\end{array}$

[0078] The mean coupling distance d.sub.C[B.fwdarw.A] of points of set B of objects of set A is therefore given by

[00012]$\begin{array}{cc}{d}_C\left[B.fwdarw.4\right]=\frac{1}{{.Math.}_{1\le j\le n_2}{p}_C\left(u_j\right)}\underset{1\le j\le n_2}{.Math.}{p}_C\left(u_j\right)d\left[d_j.fwdarw.A\right],& \left(11\right)\end{array}$

where the distance d [uj.fwdarw.A] of point uj to set A is equal to the minimum distance of uj to any pixel x; y such that the level set function Φ(x; y)=0: d [uj.fwdarw.A]=.sub.minx;y s.t. Φ(x;y)=0d[uj.fwdarw.(x,y)].

[0079] Validation with Simulations

[0080] To test the accuracy and robustness of our method, we generated synthetic sequences in Icy (http://icy.bioimageanalysis.org) with different levels of coupling between points and elongated objects.

[0081] Simulating Coupling with Elongated Objects

[0082] To simulate elongated objects with different shapes and length, we designed a stochastic algorithm where, for each object we simulate an open contour sj(p) with p comprised in [0; 1]: For each object 1≤j≤n1, we start by drawing a random initial position s(0) in the FOV Ω, a 256_256 pixels square here. Then, for 1=1::L, with L the length of objects, we compute iteratively the curvilinear position as sj(p=1/L)=sj((1−1)/L)+re.sup.iθ1 (in complex form). In our simulations r=1 and 01 follows a uniform random variable over [0;π].

[0083] We highlight that restricting the range of each independent random variable θ to [0; π] ensures that each object in the set A is rather elongated. Finally, from the simulated set of objects contours A for discretized curvilinear coordinate p=0, . . . 1/L, . . . , 1 we compute the pixelized level-set function Φ(x; y) such that Φ(x; y)=0 for each pixel (x; y) containing at least one curvilinear position [sj(1/L)] for 1≤j≤n1 and 0≤l≤L. To simulate the coupling of a given proportion 0≤α≤1 of n2 points of the set B to previously simulated objects of the set A, we simulate a Thomas process: we first randomly draw (1−α)n2 points in Ω. Then, for each remaining coupled point uk for 1≤k≤αn2, we choose randomly one simulated curvilinear position sj (l/L), and simulate the position u.sub.k.

[0084] Finally, we generate the green B sequence by binarizing the pixels that contain at least one position u.sub.k, 1≤k≤n2.

[0085] FIG. 2a shows regions for {m; n}={2; 4} (blue),{6; 8} (purple) and {10; 12} pixels (red). For each level of coupling (α=0; 0.1; . . . ; 1, n2=100 points, coupling distance˜Nμc=1 pixel; σc=0:3 pixel)), FIG. 2b compares the level of coupling computed with Ripley-based analysis (in cyan) and the present level-set method (in red). Dashed black line corresponds to the ideal method that would measure exactly the same level of coupling than the simulated one. Standard errors for n=10 independent simulations per coupling percentage are shown

[0086] Results

[0087] For different levels of points 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 (L=1 corresponding to points, L=10 and L=30).

[0088] Ideally, the measured percentage of coupling should be close to the simulated percentage, which corresponds to the black dashed line in FIG. 2. We also compare the accuracy of the method of the prior art with Ripley based analysis. We observed that, while level-sets and Ripley-based methods give similar results for pointillistic objects (L=1) and are both very close to ground truth (black dashed line), the accuracy of level-set method is maintained for longer objects whereas accuracy of Ripley-based analysis is degraded as objects' length increases. Indeed, Ripley-based analysis is based on the distance between the center-of-mass of elongated objects of set A and coupled points B, and is therefore less accurate when objects' length L increases. We highlight that our method slightly underestimates the simulated percentage of points' coupling, and that underestimation is more pronounced for longer A objects. This is due to the fact that, as the length of objects A increases, the area |ω.sub.mn| contained between level sets 0<m<n increases together with the expected number of points inside |ω.sub.mn| under complete spatial randomness. Therefore, while the total number of putative coupled points, i.e. the points that are inside level-set domains where a significant accumulation of points is detected, is close to the actual simulated number of coupled points, the coupling probability of each coupled point (Eq. 10) is decreased because the expected number of points inside |ω.sub.mn| under complete spatial randomness increases with the area of |ω.sub.mn|.

[0089] Using Wide-Field Fluorescence Microscopy to Measure the Coupling of Synaptic Spots to Neuron's Membrane

[0090] To demonstrate an application in bioimaging, we use wide-field fluorescence imaging (FIG. 3) 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.

[0091] Quantitative evaluations are performed on a set of 18 different images by designing a protocol in Icy. For spatial analysis, the set of objects (set A) is obtained by segmenting the MAP2 response using k-mean thresholding algorithm in Icy. Centroids of the Synapsin spots, automatically extracted via a wavelet-based spot detector 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 100 nm from cell boundary. Additionally, our method also predicts the uncertainty in the spatial process via the average coupling probability

[00013]$\stackrel{\_}{p_C}=\frac{1}{n_2}.Math.p_C\left(u_i\right)=56.4\%.$

This 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, and the quantitative results indeed suggest a strong association between the two markers used in this positive control.

[0092] To analyze the colocalization between points (or small spots) and complex-shape objects, we leveraged the versatility of level-sets to map objects contours and the distance of any point inside the FOV to the closest contour. After having characterized the expected points' distribution inside level-set domains under complete spatial randomness, we measured any significant accumulation of points and quantified the spatial coupling of points to objects. We validated our method with synthetic simulations, and showed that it outperforms standard Ripley-based analysis, especially for elongated objects. Finally, we highlight that the restriction of the level set function to the region of interest by using suitable boundary conditions, eliminates the need to explicitly correct for edge-artifacts. Therefore, our method provides a generic, and robust tool to study colocalization for a variety of problems commonly encountered in fluorescent microscopy, super resolution imaging, and histopathology.

Using the Method on Multicolor Images

[0093] Application 1

[0094] As shown on FIG. 4, the method according to the invention was used to determine the coupling between the CY3 channel (red) and the DAPI channel (green) response. The CY3 intensity image was reduced to binary spots (using the Icy plugin spot detector) and this serves as the first set for the method according to the invention.

[0095] The DAPI channel response was converted to point localizations by detecting the binary spots (using the Icy plugin spot detector) and estimating the positions of the center of masses of the spots.

[0096] Colocalization analysis on the above-mentioned channels using the method according to the present invention found statistically significant coupling (p-value=0) between the two channels, with a coupling index of 0.82, average coupling distance of 0.79 pixels, and an average coupling probability of 0.84. These findings support the qualitative observation of the merged channel (last column) which shows a significant response in yellow, thereby signifying co localization.

[0097] Application 2

[0098] As shown on FIG. 5, the method according to the present invention was used to determine the coupling between the FITC channel (red) and the RHOD channel (green) response. The FITC intensity image was reduced to binary spots (using the Icy plugin spot detector) and this serves as the first set for the method.

[0099] The RHOD channel response was converted to point localizations by detecting the binary spots (using the Icy plugin spot detector) and estimating the positions of the center of masses of the spots.

[0100] The method according to the invention allowed finding statistically significant coupling (p-val=1e-19) between the channels, although the coupling index is lower than the previous experiment (0.46). The couples were present at an average distance of 2.7 pixels, which explains the absence of yellow response in the merged channels, although the distribution of the individual channels indicate spatial correlation. Average Coupling probability was determined to be 0.61.

[0101] Application 3: Local Investigation Inside a Predefined Region of Interest

[0102] A local region of interest (ROI) was selected on the previous wide-field image of FIG. 5 so as to analyze the colocalization of the green channel response on the red channel. This ROI is depicted on FIG. 6. Qualitatively, it appears that the two channels are un-correlated. This is established by using the method according to the present invention which finds an average coupling of 0.03 between the red and green channels with a p-value of 0.01 and average coupling distance of 3.5 pixels. It is interesting to note that the method according to the present invention can be used on a portion of interest of an image and need not to be performed on a full image.

[0103] It is understood that the described embodiments are not restrictive and that it is possible to make improvements to the invention without leaving the framework thereof.

[0104] Unless otherwise specified, the word “or” is equivalent to “and/or”. Similarly, the word ‘one’ is equivalent to ‘at least one’ unless the contrary is specified. Unless otherwise specified, all percentages are weight percentages.

[0105] Application for Spatio-Temporal Data Analysis

[0106] In order to develop and validate predictive models for early detection of epidemics or monitoring thereof, the present document propose a method for statistical analysis of spatio-temporal point cloud data. Such spatio-temporal point-clouds may correspond to GPS coordinates of individuals, or may represent physical geolocations such as local clusters or epidemiological hot-spots. In contrast to existing techniques that measure interaction by simply counting local associations between two spatial sets, the present document provide a robust mechanism to measure statistically significant interaction between spatio-temporal data sets. Given a set of reference spatial localizations (such as a group of individuals, geographical hot-spots or disease containment zones), the method provides a temporal index of association between the reference and an observed spatial dataset realized over time.

[0107] To illustrate said method, an epidemiological case study is simulated in FIG. 7. In FIG. 7 (a), we analyze the association between reference geolocations (in green, representing, for example, disease containment zones) to synthetic spatial geo-locations of individuals (in black). Each temporal snapshot corresponds to the spatial spread of infection with respect to the containment zones. Ideally, the epidemic is considered well contained when the infected individuals (black dots) are statistically coupled to the reference. In FIGS. 1(b) and (c), we plot infection containment indices over time, corresponding to a linear, and a quadratic model respectively. The computed temporal index accurately predicts the rate of infection containment, or spread, over time. It could, therefore, be used as a robust quantitative tool to monitor disease outbreaks, and to take informed decisions based on the epidemiological findings.