Reconstructing the gradient source position from steady-state fluxes to small receptors

TL;DR

This paper investigates how a biological cell can determine the position of a diffusing source using steady-state fluxes through small receptors, providing a computational method for gradient source localization.

Contribution

It introduces a novel approach using Narrow Escape Time theory to reconstruct source position from fluxes, requiring at least three receptors within a certain distance.

Findings

Source position can be recovered with at least 3 receptors within 10 cell radii

The method applies in two-dimensional domains and narrow strips

Provides a basis for gradient sensing at the single-cell level

Abstract

Recovering the position of a source from the fluxes of diffusing particles through small receptors allows a biological cell to determine its relative position, spatial localization and guide it to a final target. However, how a source can be recovered from point fluxes remains unclear. Using the Narrow Escape Time approach for an open domain, we compute the diffusion fluxes of Brownian particles generated by a steady-state gradient from a single source through small holes distributed on a surface in two dimensions. {We find that the location of a source can be recovered when there are at least 3 receptors and the source is positioned no further than 10 cell radii away}, but this condition is not necessary in a narrow strip. The present approach provides a computational basis for the first step of direction sensing of a gradient at a single cell level.

Equations

$$\frac{\partial p_t}{\partial t}$$

$$\displaystyle p_{t}(\mbox{\boldmath$x$},\mbox{\boldmath$x$}_{0})$$

$$\displaystyle\displaystyle\frac{\partial p_{t}}{\partial\mbox{\boldmath$n$}}(\mbox{\boldmath$x$},\mbox{\boldmath$x$}_{0})$$

$$\displaystyle p_{t}(\mbox{\boldmath$x$},\mbox{\boldmath$x$}_{0})$$

$$\displaystyle-D\Delta P_{0}(\mbox{\boldmath$x$})$$

$$\displaystyle\displaystyle\frac{\partial P_{0}}{\partial\mbox{\boldmath$n$}}(\mbox{\boldmath$x$})$$

$$\displaystyle P_{0}(\mbox{\boldmath$x$})$$

$$\displaystyle P_{k}=\displaystyle\frac{\displaystyle\int_{\partial\Omega_{k}}\displaystyle\frac{\partial P_{0}(\mbox{\boldmath$x$})}{\partial\mbox{\boldmath$n$}}dS_{\mbox{\boldmath$x$}}}{\displaystyle\sum_{q}\int_{\partial\Omega_{q}}\displaystyle\frac{\partial P_{0}(\mbox{\boldmath$x$})}{\partial\mbox{\boldmath$n$}}dS_{\mbox{\boldmath$x$}}}.$$

$$\displaystyle\sum_{q}\int_{\partial\Omega_{q}}\displaystyle\frac{\partial P_{0}(\mbox{\boldmath$x$})}{\partial\mbox{\boldmath$n$}}dS_{\mbox{\boldmath$x$}}=1.$$

$$\displaystyle P_{2}=\frac{1}{2}+\frac{\pi}{2}\frac{G(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{0})-G(\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})}{\{\log|\mbox{\boldmath$x$}_{1}-\mbox{\boldmath$x$}_{2}|-\log\varepsilon\}},$$

$$\displaystyle G(\mbox{\boldmath$x$},\mbox{\boldmath$y$})=\frac{-1}{2\pi}\left(\ln|\mbox{\boldmath$x$}-\mbox{\boldmath$y$}|+\ln\left|\frac{R^{2}}{|\mbox{\boldmath$x$}|^{2}}\mbox{\boldmath$x$}-\mbox{\boldmath$y$}\right|\right).$$

$$\displaystyle r(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})=\frac{|P_{1}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})-P_{2}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})|}{P_{1}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})+P_{2}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})}$$

$$\displaystyle D_{S}=\{\mbox{\boldmath$x$}_{0}\hbox{ such that }r(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})\geq T_{h}\}.$$

$$\displaystyle f(\mbox{\boldmath$x$}_{0})=\max_{\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2}}|P_{1}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})-P_{2}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{0})|.$$

$$\displaystyle f(\mbox{\boldmath$x$}_{0})=\displaystyle\frac{2R}{|\mbox{\boldmath$x$}_{0}|\log\frac{2R}{\varepsilon}}+o\left(\frac{1}{|\mbox{\boldmath$x$}_{0}|}\right),$$

$$\displaystyle r_{3}(\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\mbox{\boldmath$x$}_{3},\mbox{\boldmath$x$}_{0})=\frac{\max\{|P_{1}-P_{2}|,|P_{1}-P_{3}|\}}{P_{1}+P_{2}+P_{3}},$$

Full text

Sensing a molecular gradient made of cue concentration is the first step to transform cell positional information into a genetic specialization and a differentiation signal Wolpert . During axonal growth and guidance, the growth cone (the tip of a neuron) uses the concentration of morphogens Goodhill1 ; reviewREingruber to decide whether or not to continue moving, stop, turn right or left. Bacteria and spermatozoa in particular are able to orient themselves in a chemotaxis gradient Kaupp1 ; Kaupp2 ; Kaupp3 ; Kaupp4 ; Heinrich3 ; Heinrich5 ; Goldstein1 ; Goldstein2 ; Powers . However, how a cell senses an external gradient concentration depends on its ability to estimate the fluxes of cues. These fluxes have been computed assuming that cues are fully or partially absorbed uniformly at the surface of a detecting ball BergPurcell . These computations are used to estimate the sensitivity of the local concentration. This, however, is insufficient to establish the orientation of the source of the gradient. Our aim is to clarify how a cell, which is only a few microns in size, can detect the direction of a source.

The first step of differentiating left from right certainly has to involve the spatial difference in the binding flux of external cues. Our model for cell direction detection uses a reflecting disk covered with small receptors. The receptors are perfect absorbers for Brownian particles (cues), emanating from a point source. Computing the fluxes of Brownian particles to small targets is part of the Narrow Escape Theory PNAS ; PPR2013 ; holcmanschuss2015 , but this theory cannot be applied directly to open and unbounded domains, because the mean passage time of particles to any small target is infinite. To avoid this difficulty, we neglect the receptor binding time of diffusing molecules and consider that cell sensing is possible via direct measurement of the diffusion flux. However, we do not account for any further cellular transduction cascade that translates receptor local activation into an internal signal. A receptor-local memory mechansim is necessary in order to prevent loss of directional information on the gradient due to homogenization of the downstream transduced signal (concentration of a second messenger or surface molecules) inside the cell. Therefore, asymmetric fluxes at the receptor level should lead to an asymmetrical transduction inside the cell. Hence, we do not replace receptors by a homogenized boundary condition that would render measurements of spatial flux differences impossible.

In this letter, we first compute the fluxes of diffusing molecules to small targets ($n=2$ and $3$) located on the surface of a detecting disk. We evaluate the effect of different receptor arrangements and also study the influence of an infinitely long, confining narrow strip. Secondly, we estimate the maximum distance from a source at which a given concentration can be detected with a pre-defined accuracy. Finally, we study how the location of a source can be recovered from the difference of fluxes. We demonstrate that the source position of a gradient can be reconstructed with three receptors, while sensing of the concentration level can be achieved with two only.

Diffusion fluxes through narrow windows. The probability density function p_{t}(\mbox{\boldmathx},\mbox{\boldmathx}_{0}) for a Brownian particle generated at location \mbox{\boldmathx}_{0} that can be absorbed on the boundary of a detecting two-dimensional disk of radius $R$, $\Omega=D(R)$, by windows $\partial\Omega_{1}\cup\ldots\cup\partial\Omega_{N}$ located on the surface $\partial\Omega$ of the disk $\Omega$ satisfies

[TABLE]

The reflecting boundary condition at $\partial\Omega-(\partial\Omega_{1}\cup\ldots\cup\partial\Omega_{N})$ accounts for the impenetrable boundary so that diffusing molecules are reflected on the surface. The absorbing boundary condition on each of the windows $\partial\Omega_{1}\cup\ldots\cup\partial\Omega_{N}$ represents rapid binding with a diffusion limited activation rate. The window sizes are identical and equal to $|\partial\Omega_{1}|=\varepsilon\ll 1$. The steady-state probability density P_{0}(\mbox{\boldmathx}) is computed by solving the mixed boundary value problem for the Laplace equation holcmanschuss2015

[TABLE]

Although the density P_{0}(\mbox{\boldmathx}) is non-normalizable in two dimensions, we are only interested in the splitting probability between windows, i.e. the normalized steady-state flux at window $k$,

[TABLE]

Due to the recurrent property of the Brownian motion in two dimensions, the probability to hit a window before going to infinity is one, thus

[TABLE]

The fluxes for $N=2$ windows can be computed using the Green’s function G(\mbox{\boldmathx},\mbox{\boldmathy}) of the domain using matched asymptotic expansion Coombs ; Ward and involves a Green’s function Matrix in general. However, using identity Eq. (4), it is sufficient to compute only one probability and we get

[TABLE]

where the external Neumann-Green’s function of a disk $\Omega=D(R)$ of radius $R$, for \mbox{\boldmathx},\mbox{\boldmathy}\in\hbox{\bb R}^{2}-B(R) is

[TABLE]

To evaluate how the probability $P_{2}$ changes with the distance of the source \mbox{\boldmathx}_{0} and the relative position of the windows, we compare Brownian simulations with the analytical expression (5) for a disk (Fig. 1A). For the simulations, we generated Brownian trajectories near the disk (on the surface of a disk of radius $R_{e}$) according to the exit point distribution of a process from an internal disk schuss2013 . Interestingly, already at a distance of $L=10R$, the absolute difference between the fluxes $\Delta P=|P_{1}-P_{2}|$ is within $5\%$, making it almost impossible to determine source direction or concentration differences in a noisy environment. The results are independent of the window positions and is qualitatively the same in dim 2 and 3 (compared Fig. 1B-C-D to E-F-G). Moreover, $\Delta P\to 0$ as $L$ increases, see Fig. 1B-D.

In contrast, when the disk is located in a narrow strip (Fig. 2A), the difference of fluxes between the two windows converges asymptotically to a finite difference $\Delta P$ depending on the strip width $a$, even for large source distances $L\geq 100R$ (see Fig. 2B-D.). Indeed, the fluxes hardly show any dependence on source distance $L$. The narrow funnel PPR2013 between the strip and the disk prevents Brownian particles to reach a window located on the opposite side of the disk, leading to the observed effects.

To further investigate how the window positions could influence the recovery of the source location, we estimated the maximum distance between the source and the disk containing two absorbing windows located at position \mbox{\boldmathx}_{1},\mbox{\boldmathx}_{2} that gives a significant difference of probability flux. For that purpose, we define the sensitivity ratio as

[TABLE]

(note that here P_{1}(\mbox{\boldmathx}_{1},\mbox{\boldmathx}_{2},\mbox{\boldmathx}_{0})+P_{2}(\mbox{\boldmathx}_{1},\mbox{\boldmathx}_{2},\mbox{\boldmathx}_{0})=1). The domain of sensitivity for a threshold $T_{h}$ is the interior of the two-dimensional region

[TABLE]

We plotted the boundary of the region $D_{S}$ for two absorbing windows symmetrically positioned (Fig. 3A) and when the angle is $\theta_{12}=\pi/2$ (Fig. 3B) for $T_{h}=1\%,5\%$ and $10\%$. The region $D_{S}$ consists of two connected components and no detection (for a threshold below $T_{h}$) is possible outside $D_{S}$. Interestingly, with a $1\%$ precision, the domain is around 20$\times$ the size of the detecting disk. Beyond this distance, no detection is possible.

The detection sensitivity for a given source location \mbox{\boldmathx}_{0} and optimal window placement is defined by

[TABLE]

The maximum is achieved for a window configuration aligned with the position of the source and symmetric with respect to the center of the disk centered at the origin. An explicit computation with \mbox{\boldmathx}_{2}=-\mbox{\boldmathx}_{1}, |\mbox{\boldmathx}_{1}|=|\mbox{\boldmathx}_{2}|=R gives

[TABLE]

where $R$ is the radius of the disk. Hence, the detection sensitivity decreases as the reciprocal of the distance to the source L=|\mbox{\boldmathx}_{0}|. With three windows, detection is possible if at least one of the difference between the splitting probability is higher than the threshold $T_{h}$. We thus define the new sensitivity using

[TABLE]

with $P_{1}$, $P_{2}$ and $P_{3}$ defined above depend on \mbox{\boldmathx}_{1}, \mbox{\boldmathx}_{2}, \mbox{\boldmathx}_{3} and \mbox{\boldmathx}_{0}. Note that this definition illustrates the range of sensitivity and it might be possible that cells perform using biochemical reactions this computation. The detectability region is now completely surrounding the detecting disk $\Omega$ and the boundaries are much larger than in the two-windows case. Numerical simulations show the region $D_{S}$ defined in Eq. (8) for the function $r_{3}$ with windows positioned at the corners of an equilateral triangle is now connected and seems to extend to 40 times the size of the detecting disk (Fig. 3C). Finally, reconstructing the location of the source in the detectable region $D_{S}$ from the steady-state probability fluxes, requires inverting Eq. (5) and to find \mbox{\boldmathx}_{0} when the fluxes are known. With two windows located on a detecting disk and using the expression of the Green’s function (6), we obtain a one dimensional curve (Fig. 4A). At least three windows are required to recover the point source, located at the intersection of two curves (Fig. 4B). Indeed, the recurrent Brownian motion in dimension 2 schuss2013 implies that $P_{1}+P_{2}+P_{3}=1$, thus we only need to estimate P1 and P2 from Eq. (2), by the matched asymptotic method, except that now we have to invert a three-dimensional matrix equation (see SI). We use the expression Eq. (6) to invert the equations and the intersection point is found numerically via MINPACK’s multidimensional nonlinear root finding method hybrj hybrj . Interestingly, fluctuations in the probability fluxes $P_{i}$ (implemented by changing $P_{i}\to[1\pm\eta]P_{i}$) yields a nonlinear and spatially inhomogeneous uncertainty in the reconstruction of the source position (overlap of the shaded areas in Fig. 4B).

Discussion and Conclusion. The ability of a cell to sense a gradient concentration is mediated by the binding of Brownian molecules to receptors BergPurcell ; HS2005 . Computing the diffusion fluxes via homogenizing over local receptor positions Endres ; Goodhill1 ; Berez , renders a recovery of any directional information impossible, as it assigns the same flux to the entire boundary of the detecting disk. Based on Narrow Escape Theory Ward ; Coombs ; holcmanschuss2015 , we estimated the probability fluxes on each individual receptor window separately and found that in two dimensions (i.e. a flat environment), the direction of the source can be recovered. This mechanism requires a comparison between the fluxes of at least three boundary receptors. In addition, the Green’s function approach allows us to estimate the boundary of the region of sensitivity characterized by a difference of fluxes between receptors being larger than a predefined sensitivity threshold. The sensitivity decays with the reciprocal of the distance to the source according to Eq. (10). Furthermore, we evaluated the effect of the external geometry on the threshold of sensitivity: although a disk of radius R with absorbing windows cannot sense the source position located beyond ten cell radii, in a narrow band, the detection is possible due to narrow passages for the diffusing molecules PPR2013 , as long as one detecting window is facing away from the source.

The present method can be used for several applications, such as growth-cone navigation inside the developing brain. Neurons have to travel millimeters to centimeters to find the correct cortical regions and to form synaptic connections reviewREingruber ; Goodhill1 ; Yogev . We propose that narrow extracellular tubes formed of neurons and glial cells probably allows the genesis of shallow gradients detected by receptors located on the growth cone, although it is not known precisely how these receptors are organized. The dynamics of the growth cone including moving protrusions is certainly an additional mechanism worthy of further investigation in context of direction sensing, but should dominate for short-range distances only.

Although we focussed on two and three windows only, the results would be similar for two or three receptor clusters holcmanschuss2015 . The model we used here was developed in the fast binding limit without rebinding Kaupp1 . It remains a challenge to apply the present theory to understand how bacteria Heinrich3 , sperm Kaupp2 or neurite growth Yogev localize the source of their cues. We limited the present approach to the initial level of source detection, however the asymmetry of receptor detection needs to persist in the downstream transduction, which certainly is another key question to investigate.