Triangulation sensing

Triangulation sensing is a theory describing the computational steps of a cell containing on its surface small windows, to recover the location of a source emitting random particles in a medium. This is particularly relevant for neuron navigation in the Brain .Reconstructing the source location allows the cell to triangulate its position. The reconstruction steps of the gradient source from the fluxes of diffusing particles arriving to small absorbing receptors are 1-arrival of the Brownian particles to the small windows 2- counting of the particles 3-inversion of the Laplace's equation to estimate the position from the fluxes 4- possible noise reduction by applying the same procedure to several triplets Dobramysl, U., & Holcman, D. (2018). Reconstructing the gradient source position from steady-state fluxes to small receptors. Scientific reports, 8(1), 1-8. .
== Computing the diffusion fluxes to small window with respcet to the source location ==
The mathematical formulation consists inconsidering 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. For a ball B(a) of radius a. M can range between 10 and 50, but can be much higher for other receptor types. Individual cue molecules are described as Brownian particles. The cues are released from a source at position x_0 outside the ball. The triangulation sensing method consists in reconstructing the source position from estimated steady-state fluxes at each narrow window for fast binding (i.e. the probability density has an absorbing boundary condition at the windows).Shukron, O., Dobramysl, U., & Holcman, D. (2019). Chemical Reactions for Molecular and Cellular Biology. Chemical Kinetics: Beyond The Textbook, 353.
The steady-state gradient P_0(x) is obtained by integrating the diffusion equation from 0 to infinity, which is equivalent to resetting a particle after it disappears through a window. It is given by the solution of the mixed-boundary value problem (where D is the diffusion coefficient)
D\Delta P_0(x) = -Q \delta(x-x_0) \hbox{ for } x\, \in \,\R^3 -\Omega 
\frac{\partial P_0}{\partial n}(x) = 0\hbox{ for } x, \in\, \partial\Omega_r
P_0(x) = 0 \hbox{ for } x \in \partial \Omega_a.
The probability fluxes associated with P_0(x) at each individual window \Omega_a(\epsilon_j) can be computed from the fluxes and depend on the specific window arrangement and the domain \Omega. 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 P_0(x) has to tend to zero in three dimensions. More complex domains could be studied if their associated Green's function can be found explicitly.
Note that in dimension, although the probability density P_0(x) diverges when x\rightarrow\infty, 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:
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}}.
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:
\sum_{q} \int_{\partial\Omega_{q}} \frac{\partial P_0(x)}{\partial n} dS_{x}=1.
The construction of the solution uses an inner and outer solution. The inner solution is constructed near each small window by scaling the Arc length s and the distance to the boundary \eta by \bar \eta=\frac{\eta}{\epsilon} and \bar s=\frac{s}{\epsilon}, so that the inner problem reduces to the classical two-dimensional Laplace equation
\Delta w=0 \hbox{ in } \R_+^{2}, 
\frac{\partial w}{\partial n}=0 \hbox{ for } |\bar s|>\frac{1}{2}, \bar \eta=0 , 
w(\bar s, \bar \eta)=0 \hbox{ for } |\bar s|<\frac{1}{2}, \bar \eta=0.
The far field behavior for |x|\rightarrow \infty and for each hole i=1, 2 is w_i(x) \approx A_i \{\log|x-x_i|-\log \epsilon\} +o(1), where A_i is the flux A_i=\frac{2}{\pi} \int_0^{1/2} \frac{\partial w(0,\bar s)}{\partial \bar \eta}d\bar s.
The general solution of equation for n=2 is obtained from the outer solution of the external Neumann-Green's function
-\Delta_{x} G(x,y) = \delta(x-y) \hbox{ for } x \in \,\R_{+}^2
\frac{\partial G}{\partial n_{x}}(x,x_0) = 0\hbox{ for } x \in {\partial\R_{+}^2}.
Given for x, x_0 \in \R_{+}^2 by G(x,x_0)=\frac{-1}{2\pi }\left(\ln |x-x_0| + \ln \left|x-\bar{x_0}\right| \right), where \bar{x_0}is the symmetric image of x_0through the boundary axis 0z. The uniform solution is the sum of inner and outer solution (Neuman-Green's function) 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, where A_1,A_2,C are constants to be determined. To that purpose, we study the behavior of the solution near each point x_i. In the boundary layer, we get
P(x,y) \approx A_i \{\log|x-x_i|-\log \epsilon\}.
Using this condition on each window, we obtain the two conditions:
G(x_1,x_0)+A_2 \{\log|x_1-x_2|-\log \epsilon\} +C=0
G(x_2,x_0) +A_1 \{\log|x_2-x_1|-\log \epsilon\} +C=0.
Due to the recursion property of the Brownian motion in dimension 2, there are no fluxes at infinity, thus the conservation of flux gives: \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. In the limit of two well separated windows (|x_1-x_2| \gg 1), using the condition for the flux, we get for each window i=1,2, \int_{\partial \Omega_{i}} \frac{\partial P(x,y))}{\partial n}dS_{x}= -\pi A_i (the minus sign is due to the outer normal orientation), thus \pi A_1+\pi A_2=1. We finally obtain the system
\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}.
Finally, the absorbing probabilities are given by 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\}}andP_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\}}.
These probabilities precisely depend on the source position x_0 and the relative position of the two windows. When one of the splitting probabilities (either P_1 or P_2) is known and fixed in [0,1], recovering the position of the source requires inverting the previous equations. For P_2=\alpha \in [0,1], the position x_0 lies on the curve
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) \}. To conclude knowing the splitting probability between two windows is not enough to recover the exact the source x_0 , 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 probleme: computing 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 P_1=P_2, leading to a one dimensional curve degeneracy for the reconstructed source positions x_0. Reconstructing the source location x_0 from the splitting probabilities, requires to invert a system equation. Starting from the general solution as given above, we get for three windows:
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,
where A_1,A_2,A_3,C are constants to be determined. The three absorbing boundary conditions forP(x,x_0)leads to the system of equations
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
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
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,
with the normalization condition for the fluxes \pi A_1+\pi A_2+\pi A_3=1. With \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}, and using the notations d_{ij}=|x_i-x_j|, we obtain the solution
A_1=\frac{1}{\pi}-A_2 -A_3.
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}}.
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}}
To conclude, we can specify the flux values \alpha>0 and \beta>0 such that\alpha+\beta<1, with \alpha=\pi A_1=\int_{\partial \Omega_{1}} \frac{\partial P(x,y))}{\partial n}dS_{x}and \beta=\pi A_2=\int_{\partial \Omega_{2}} \frac{\partial P(x,y))}{\partial n}dS_{x}. We now provide the final step which is to express the source x_0with respect to the flux coefficients A_1,A_2,A_3.With three windows in a x-y plane in the confirguration x_1=(0, 0, 0), x_2=(d, 0, 0)\, \hbox{and} \, x_3=(e, f, 0) (window 1 is at the origin and window 2 is on the x-axis). Then, using the leading order from the expansion of the fluxes, the location of the source is the solution of the three non-linear equations
\gamma_1^2 = (x_0^{(1)})^2 + (x_0^{(2)})^2 + (x_0^{(3)})^2
\gamma_2^2 = (d-x_0^{(1)})^2 + (x_0^{(2)})^2 + (x_0^{(3)})^2
\gamma_3^2 = (e-x_0^{(1)})^2 + (f-x_0^{(2)})^2 + (x_0^{(3)})^2,
where \gamma_i=\frac{2\epsilon}{\pi\Phi_i}. Solving for the coordinates of \x_0 and requiring that x_0^{3}>0, leads to the analytical solutionx_0^{(1)}=\frac{d^2+\gamma_1-\gamma_2}{2d}
x_0^{(2)}= \frac{1}{2df}\left[d(e^2+f^2+\gamma_1-\gamma_3)-e(d^2+\gamma_1-\gamma_2)\right]
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}.
Thus the position can be recovered from the steady-state fluxesDobramysl, U., & Holcman, D. (2019). Triangulation sensing: how cells recover a source from diffusing particles in three dimensions. arXiv preprint arXiv:1911.02907. Dobramysl, U., & Holcman, D. (2018). Mixed analytical-stochastic simulation method for the recovery of a Brownian gradient source from probability fluxes to small windows. Journal of computational physics, 355, 22-36. DOBRAMYSL, Ulrich et HOLCMAN, David. Triangulation Sensing to Determine the Gradient Source from Diffusing Particles to Small Cell Receptors. Physical Review Letters, 2020, vol. 125, no 14, p. 148102.. 
In general, with N windows, no explicit inverse is easily available, hence numerical procedures are used to find the position of the source x_0 at any order in \epsilon. To find the position of the source x_0 from the measured fluxes \Phi_i, i=1...N,we need to invert the previous Eqs. Each of these equations describes a non-planar surface S_{i}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 S_{i}and S_{j} intersect, forming three-dimensional curves C_{ij}and all of these curves intersect at the location of the source x_0. In the 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. Any combination leads to the same source position.