Triangulation sensing

Short description: Mathematical theory

Triangulation sensing is a mathematical theory describing the computational steps to recover the location of a source from the fluxes estimated on small windows located on a test surface. The source emits random particles in a medium. The reconstruction steps of the gradient source from the fluxes of diffusing particles arriving to small absorbing receptors are decomposed in several steps:

1. Estimating the arrival of Brownian particles to the small windows located on the test surface.
2. Inverting the Laplace's equation to estimate the position from the fluxes based on a triplet of windows.
3. Reduce the uncertainty by averaging the noise on several triplets.

This theory provides a framework to study neural navigation in the brain.

Computing the diffusion fluxes to small window with respect to the source location

The mathematical formulation consists in considering diffusing cues that have to bind to N narrow windows located on the surface of a three dimensional shaped object, typically a ball (in dimension 3) or a disk in dimension 2. M windows ranging between 10 and 50 are located on a ball B(a) of radius a. Individual cue molecules are described as Brownian particles. The cues are released from a source at position [math]\displaystyle{ x_0 }[/math] outside the ball. The triangulation sensing method consists in reconstructing the source position from the steady-state fluxes estimated at each narrow window for fast binding, i.e. the probability density function has an absorbing boundary condition at the windows.

Steady-state Laplace's equation

The steady-state gradient [math]\displaystyle{ P_0(x) }[/math] is obtained by integrating the diffusion equation from 0 to infinity, which is equivalent to resetting the diffusing particles after they disappear through a window. It is given by the solution of the mixed-boundary value problem (where D is the diffusion coefficient)

[math]\displaystyle{ D\Delta P_0(x) = -Q \delta(x-x_0) \hbox{ for } x\, \in \,\R^3 -\Omega }[/math]

[math]\displaystyle{ \frac{\partial P_0}{\partial n}(x) = 0\hbox{ for } x, \in\, \partial\Omega_r }[/math]

[math]\displaystyle{ P_0(x) = 0 \hbox{ for } x \in \partial \Omega_a. }[/math]

The probability fluxes associated with [math]\displaystyle{ P_0(x) }[/math] at each individual window [math]\displaystyle{ \Omega_a(\epsilon_j) }[/math] of size [math]\displaystyle{ \epsilon_j }[/math] can be computed from the fluxes and depend on the specific window arrangement and the surface of the domain [math]\displaystyle{ \Omega }[/math]. The parameter Q>0 can be calibrated so that there are a fixed number of particles located in a given volume. At infinity, the density [math]\displaystyle{ P_0(x) }[/math] has to tend to zero in three dimensions. More complex domains could be studied if their associated Green's function can be found explicitly.

In dimension 2, although the probability density [math]\displaystyle{ P_0(x) }[/math] diverges when [math]\displaystyle{ x\rightarrow\infty }[/math], the splitting probability between windows is finite because it is the ratio of the steady-state flux at each window divided by the total flux through all windows:

[math]\displaystyle{ J_k= \frac{ \int_{\Omega_{k}} \frac{\partial P_0(x)}{\partial n} dS_{x}}{\sum_{q} \int_{\partial \Omega_{q}} \frac{\partial P_0(x)}{\partial n} dS_{x}}. }[/math]

In two-dimensions, due to the recurrent property of the Brownian motion, the probability to hit a window before going to infinity is one, thus the total flux is one:

[math]\displaystyle{ \sum_{q} \int_{\partial\Omega_{q}} \frac{\partial P_0(x)}{\partial n} dS_{x}=1. }[/math]

The construction of the solution uses an inner and outer solution describe below.

Boundary layer analysis in dimension 2

The inner solution is constructed near each small window by scaling the arclength s and the distance to the boundary [math]\displaystyle{ \eta }[/math] by [math]\displaystyle{ \bar \eta=\frac{\eta}{\epsilon} }[/math] and [math]\displaystyle{ \bar s=\frac{s}{\epsilon} }[/math], so that the inner problem reduces to the classical two-dimensional Laplace equation

[math]\displaystyle{ \Delta w=0 \hbox{ in } \R_+^{2}, \frac{\partial w}{\partial n}=0 \hbox{ for } |\bar s|\gt \frac{1}{2}, \bar \eta=0 , w(\bar s, \bar \eta)=0 \hbox{ for } |\bar s|\lt \frac{1}{2}, \bar \eta=0. }[/math]

The far field behavior for [math]\displaystyle{ |x|\rightarrow \infty }[/math] and for each hole i=1, 2 is [math]\displaystyle{ w_i(x) \approx A_i \{\log|x-x_i|-\log \epsilon\} +o(1), }[/math] where A_i is the flux [math]\displaystyle{ A_i=\frac{2}{\pi} \int_0^{1/2} \frac{\partial w(0,\bar s)}{\partial \bar \eta}d\bar s. }[/math]

The general solution of equation for n=2 is obtained from the outer solution of the external Neumann-Green's function

[math]\displaystyle{ -\Delta_{x} G(x,y) = \delta(x-y) \hbox{ for } x \in \,\R_{+}^2 }[/math]

[math]\displaystyle{ \frac{\partial G}{\partial n_{x}}(x,x_0) = 0\hbox{ for } x \in {\partial\R_{+}^2}. }[/math]

Given for [math]\displaystyle{ x, x_0 \in \R_{+}^2 }[/math] by [math]\displaystyle{ G(x,x_0)=\frac{-1}{2\pi }\left(\ln |x-x_0| + \ln \left|x-\bar{x_0}\right| \right), }[/math] where [math]\displaystyle{ \bar{x_0} }[/math]is the symmetric image of [math]\displaystyle{ x_0 }[/math]through the boundary axis 0z.

Construction of the general solution

The uniform solution is the sum of inner and outer solution (Neuman-Green's function) [math]\displaystyle{ P(x,x_0)=G(x,x_0) +A_1 \{\log|x-x_1|-\log \epsilon\}+A_2 \{\log|x-x_2|-\log \epsilon\} +C, }[/math] where [math]\displaystyle{ A_1,A_2,C }[/math] are constants to be determined.

To that purpose, the behavior of the solution near each point [math]\displaystyle{ x_i }[/math] is crutial. In the boundary layer, we get

[math]\displaystyle{ P(x,y) \approx A_i \{\log|x-x_i|-\log \epsilon\}. }[/math]

Using this condition on each window, we obtain the two conditions:

[math]\displaystyle{ G(x_1,x_0)+A_2 \{\log|x_1-x_2|-\log \epsilon\} +C=0 }[/math]

[math]\displaystyle{ G(x_2,x_0) +A_1 \{\log|x_2-x_1|-\log \epsilon\} +C=0. }[/math]

Due to the recursion property of the Brownian motion in dimension 2, there are no fluxes at infinity, thus the conservation of flux gives: [math]\displaystyle{ \int_{\partial\Omega_{1}} \frac{\partial P(x,y)}{\partial n}dS_{x} + \int_{\partial\Omega_{2}} \frac{\partial P(x,y)}{\partial n}dS_{x} =-1. }[/math] In the limit of two well separated windows ([math]\displaystyle{ |x_1-x_2| \gg 1) }[/math], using the condition for the flux, we get for each window i=1,2, [math]\displaystyle{ \int_{\partial \Omega_{i}} \frac{\partial P(x,y))}{\partial n}dS_{x}= -\pi A_i }[/math] (the minus sign is due to the outer normal orientation), thus [math]\displaystyle{ \pi A_1+\pi A_2=1. }[/math] We finally obtain the system

[math]\displaystyle{ \begin{cases} \frac{G(x_1,x_0)-G(x_2,x_0)}{\{\log|x_1-x_2|-\log \epsilon\}} +(A_2-A_1) =0 \\ \\ A_1+ A_2=\frac{1}{\pi}. \end{cases} }[/math]

To conclude, the absorbing probabilities are given by

[math]\displaystyle{ \begin{cases} P_1=\pi A_1=\frac12+ \frac{\pi}{2}\frac{G(x_2,x_0)-G(x_1,x_0)}{\{\log|x_1-x_2|-\log \epsilon\}} =\frac12- \frac{1}{4}\frac{\ln \frac{|x_2-x_0|\left|x_2-\bar{x_0}\right|}{|x_1-x_0|\left|x_1-\bar{x_0}\right|}}{\{\log|x_2-x_1|-\log \epsilon\}}. \\\\ P_2=\pi A_2=\frac12+ \frac{\pi}{2}\frac{G(x_1,x_0)-G(x_2,x_0)}{\{\log|x_1-x_2|-\log \epsilon\}} =\frac12- \frac{1}{4}\frac{\ln \frac{|x_1-x_0|\left|x_1-\bar{x_0}\right|}{|x_2-x_0|\left|x_2-\bar{x_0}\right|}}{\{\log|x_1-x_2|-\log \epsilon\}} . \end{cases} }[/math]

These probabilities precisely depend on the source position [math]\displaystyle{ x_0 }[/math] and the relative position of the two windows. When one of the splitting probabilities (either [math]\displaystyle{ P_1 }[/math] or [math]\displaystyle{ P_2 }[/math]) is known and fixed in [0,1], recovering the position of the source requires inverting the previous equations. For [math]\displaystyle{ P_2=\alpha \in [0,1], }[/math] the position [math]\displaystyle{ x_0 }[/math] lies on the curve

[math]\displaystyle{ S_{source}=\{ x_0 \hbox{ such that} \, \frac{|x_1-x_0|\left|x_1-\bar{x_0}\right|}{|x_2-x_0|\left|x_2-\bar{x_0}\right|}=\exp \left((4\alpha-2)\{\log|x_1-x_2|-\log \epsilon\} \right) \}. }[/math]

To conclude knowing the splitting probability between two windows is not enough to recover the exact the source [math]\displaystyle{ x_0 }[/math] , because it lies on one dimensional curve solution. However the direction can be obtained by simply checking which one of the two probability is the highest.

Inverse problem: computing the source location by triangulation

To reconstruct the location of a source from the measured fluxes, at least three windows are needed. With two windows only, a source located on the line perpendicular to one of the connecting windows would, for example generate the same splitting probability [math]\displaystyle{ P_1=P_2 }[/math], leading to a one dimensional curve degeneracy for the reconstructed source positions [math]\displaystyle{ x_0 }[/math].

Reconstructing the source location [math]\displaystyle{ x_0 }[/math] from the splitting probabilities, requires to invert a system equation. Starting from the general solution as given above, we get for three windows:

[math]\displaystyle{ P(x,x_0)=G(x,x_0) +A_1 \{\log|x-x_1|-\log \epsilon\}+A_2 \{\log|x-x_2|-\log \epsilon\}+A_3 \{\log|x-x_3|-\log \epsilon\} +C, }[/math]

where [math]\displaystyle{ A_1,A_2,A_3,C }[/math] are constants to be determined. The three absorbing boundary conditions for[math]\displaystyle{ P(x,x_0) }[/math]leads to the system of equations

[math]\displaystyle{ G(x_1,x_0) +A_2 \{\log|x_1-x_2|-\log \epsilon\}+A_3 \{\log|x_1-x_3|-\log \epsilon\} +C=0 }[/math]

[math]\displaystyle{ G(x_2,x_0) +A_1 \{\log|x_2-x_1|-\log \epsilon\} +A_3 \{\log|x_2-x_3|-\log \epsilon\}+C=0 }[/math]

[math]\displaystyle{ G(x_3,x_0) +A_1 \{\log|x_1-x_3|-\log \epsilon\}+A_2 \{\log|x_2-x_3|-\log \epsilon\} +C=0, }[/math]

with the normalization condition for the fluxes [math]\displaystyle{ \pi A_1+\pi A_2+\pi A_3=1. }[/math] With [math]\displaystyle{ \Delta_{123}=\left(\log\frac{d_{13}d_{12}}{d_{32}\epsilon}\right)^2-4\log \frac{d_{12}}{\epsilon} \log \frac{d_{13}}{\epsilon} }[/math], and using the notations [math]\displaystyle{ d_{ij}=|x_i-x_j| }[/math], we obtain the solution

[math]\displaystyle{ A_1=\frac{1}{\pi}-A_2 -A_3. }[/math]

[math]\displaystyle{ A_2=\frac{ \log\frac{d_{13}d_{12}}{d_{32}\epsilon}(G_{30}-G_{10} +\frac{1}{\pi} \log\frac{d_{13}}{\epsilon})-( G_{1 0}-G_{2 0}+\frac{1}{\pi}\log\frac{d_{12}}{\epsilon})\log\frac{d_{13}^2}{\epsilon^2 })}{\Delta_{123}}. }[/math]

[math]\displaystyle{ A_3=\frac{ \log\frac{d_{13}d_{12}}{d_{32}\epsilon}(G_{20}-G_{10} +\frac{1}{\pi} \log\frac{d_{12}}{\epsilon})-( G_{1 0}-G_{3 0}+\frac{1}{\pi}\log\frac{d_{13}}{\epsilon}) \log\frac{d_{12}^2}{\epsilon^2})}{\Delta_{123}} }[/math]

To conclude, we can specify the flux values [math]\displaystyle{ \alpha\gt 0 }[/math] and [math]\displaystyle{ \beta\gt 0 }[/math] such that [math]\displaystyle{ \alpha+\beta\lt 1 }[/math], with [math]\displaystyle{ \alpha=\pi A_1=\int_{\partial \Omega_{1}} \frac{\partial P(x,y))}{\partial n}dS_{x} }[/math]and [math]\displaystyle{ \beta=\pi A_2=\int_{\partial \Omega_{2}} \frac{\partial P(x,y))}{\partial n}dS_{x}. }[/math]

To conclude, when two fluxes are know, the flux coefficients [math]\displaystyle{ A_1,A_2,A_3 }[/math]can be computed and the position [math]\displaystyle{ x_0 }[/math] can be reconstructed from the three equations.

Generalization to three dimensions

With three windows in a x-y plane in the confirguration [math]\displaystyle{ x_1=(0, 0, 0), x_2=(d, 0, 0)\, \hbox{and} \, x_3=(e, f, 0) }[/math] (window 1 is at the origin and window 2 is on the x-axis). Using the leading order from the expansion of the fluxes, the location of the source is the solution of the three non-linear equations

[math]\displaystyle{ \gamma_1^2 = (x_0^{(1)})^2 + (x_0^{(2)})^2 + (x_0^{(3)})^2 }[/math]

[math]\displaystyle{ \gamma_2^2 = (d-x_0^{(1)})^2 + (x_0^{(2)})^2 + (x_0^{(3)})^2 }[/math]

[math]\displaystyle{ \gamma_3^2 = (e-x_0^{(1)})^2 + (f-x_0^{(2)})^2 + (x_0^{(3)})^2, }[/math]

where [math]\displaystyle{ \gamma_i=\frac{2\epsilon}{\pi\Phi_i} }[/math]. Solving for the coordinates of [math]\displaystyle{ x_0 }[/math] and requiring that [math]\displaystyle{ x_0 }[/math].x >0, leads to the analytical solution

[math]\displaystyle{ x_0^{(1)}=\frac{d^2+\gamma_1-\gamma_2}{2d} }[/math]

[math]\displaystyle{ x_0^{(2)}= \frac{1}{2df}\left[d(e^2+f^2+\gamma_1-\gamma_3)-e(d^2+\gamma_1-\gamma_2)\right] }[/math]

[math]\displaystyle{ x_0^{(3)}= \frac{1}{2df}\left[(e^2+f^2)(\{\gamma_1-\gamma_2\}^2-d^4)+2de(e^2+f^2+\gamma_1-\gamma_3)(d^2+\gamma_1-\gamma_2)\right. \left.-d^2(e^4+f^4+[\gamma_1-\gamma_3]^2+2e^2[f^2+2\gamma_1-\gamma_2-\gamma_3]-2f^2[\gamma_2+\gamma_3])\right]^{1/2}. }[/math]Thus the position can be recovered from the steady-state fluxes.

Conclusion

In general, with N windows, no explicit inverse is easily available, hence numerical procedures are used to find the position of the source [math]\displaystyle{ x_0 }[/math] at any order in [math]\displaystyle{ \epsilon }[/math]. To find the position of the source [math]\displaystyle{ x_0 }[/math] from the measured fluxes [math]\displaystyle{ \Phi_i, i=1...N, }[/math]we need to invert a system of equations. Each of these equations describe a non-planar surface [math]\displaystyle{ S_{i} }[/math] in three dimensions, corresponding to window i and intersecting the half-plane (in the case of the windows located on the half-plane) or the unit ball (in the case of the windows located on the ball). Each pair of surfaces [math]\displaystyle{ S_{i} }[/math]and [math]\displaystyle{ S_{j} }[/math] intersect, forming three-dimensional curves [math]\displaystyle{ C_{ij} }[/math]and all of these curves intersect at the location of the source [math]\displaystyle{ x_0 }[/math] . In case of N>3 windows, the system is overdetermined and we can simply choose any combination k, l and m of three fluxes from the N available windows. Any combination lead to the same source position.