Mixed analytical-stochastic simulation method for the recovery of a Brownian gradient source from probability fluxes to small windows.

Is it possible to recover the position of a source from the steady-state fluxes of Brownian particles to small absorbing windows located on the boundary of a domain? To address this question, we develop a numerical procedure to avoid tracking Brownian trajectories in the entire infinite space. Instead, we generate particles near the absorbing windows, computed from the analytical expression of the exit probability. When the Brownian particles are generated by a steady-state gradient at a single point, we compute asymptotically the fluxes to small absorbing holes distributed on the boundary of half-space and on a disk in two dimensions, which agree with stochastic simulations. We also derive an expression for the splitting probability between small windows using the matched asymptotic method. Finally, when there are more than two small absorbing windows, we show how to reconstruct the position of the source from the diffusion fluxes. The present approach provides a computational first principle for the mechanism of sensing a gradient of diffusing particles, a ubiquitous problem in cell biology.

Introduction

Recovering the source location from incomplete information about the emitting signal is a generic problem in several fields of science, such as finding an emitter in signal processing, the food source by smelling a few molecules and many more. In the context of cell biology, the question of how a cell can sample its environment and decide its final destination remains open, but it starts with the detection of an external gradient concentration that the cell must use to transform cell positional information into its genetic specialization and differentiation [30], [12].

During axonal growth and guidance, the growth cone (which is the tip of a neuronal cell) uses external concentration gradients [10], [25] to decide whether to continue moving or to stop, to turn right or left. Bacteria and spermatozoa can orient themselves in various chemotactical or mechanical gradients [2], [11]. However, most models in the current literature that are concerned with addressing these questions rely on computing the flux to an absorbing or reflecting ball [4], an absorbing or permeable ball [9], [1], or a single receptor sphere [17], all of which is insufficient to differentiate between concentrations to the left or right of the cell. To enable sensing of this difference, the detectors, modeled here as small absorbing windows, should be considered individually.

We compute here in the first part the steady-state fluxes of Brownian particles to small absorbing windows located on the boundary of a an infinite domain. Computing the fluxes of Brownian particles moving inside a bounded domain to small absorbing windows located on a boundary falls into the narrow escape problems [19], [13], [15], [16], [7], [21] and has also been studied numerically [20]. However, the mean passage time to a small hole becomes infinite in an unbounded domain due to long excursions to infinity of Brownian trajectories. This difficulty is resolved here by computing the flux directly using two methods: first, we compute the flux of Brownian particles to small absorbers located on the half-plane, a disk in and in a narrow band. The asymptotic computations are obtained by matched asymptotics of Laplace's equation in infinite domains.

In the second part, we develop a mixed numerical procedure to avoid tracking Brownian trajectories in the entire infinite space. We generate particles near the absorbing windows, computed from the analytical expression of the exit probability on an artificial boundary without introducing any artifacts [23], [24]. This method avoids the costly computation of particle trajectories in the unbounded environment (e.g. extracellular space in the brain or cells moving in two dimensional chamber), containing large excursions away from the cell, thereby allowing direct simulations of Brownian trajectories in the region of interest close to the cell. In the absence of such a procedure, these simulations would be next to impossible to perform due to the aforementioned infinite mean passage time.

We show that the results of both independent methods (Asymptotic and numerical) agree. The local geometry and distribution of windows does matter for the reconstruction of the source position: we show that it is indeed possible to recover the source of a gradient already with three receptors. Finally, the location of the windows might also be critical for the sensitivity of detection: for example, the flux of Brownian particles to small targets depends crucially on their localization [15], [14], [7], [21], [16], [19]. In summary, the manuscript is organized as follows. First, we compute asymptotically the flux of Brownian particles to receptors. Second, we introduce the mixed simulation method. In the third part, we present several applications to various geometry: half-space, a disk and a disk in a narrow band. In the fourth and last section, we apply the methods to reconstruction the source location.

Fluxes of Brownian particles to small targets in an open space

Brownian molecules are produced by a steady-state source located at position in an open space such as the two-dimensional real space . The steady-state distribution of particles, , is the solution of the Green's function where the parameter measures the injection rate of particles. We study here the flux received by an obstacle Ω containing N-small absorbing windows on its boundary ∂Ω. The fluxes of diffusing particles on the windows can be computed from solving the mixed boundary value problem (we now set ) [16] The reflecting boundary condition accounts for the impenetrable walls and diffusing molecules are reflected on the surface . The absorbing boundary condition on each window represents the extreme case where the binding time of particles is fast and the particle trajectories are terminated.

Although the probability density diverges when , we are interested in the splitting probability between windows, which is the ratio of the steady-state flux at each hole divided by the total flux through all windows: 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: We shall now compute the fluxes asymptotically for three different configurations: 1 – when the windows are distributed on a line at the boundary of the half-plane, 2 – when there are located on a disk in the entire space, and 3 – when the disk is located in a narrow band. We use the Green–Neumann's function and the method of matched asymptotics [19], [21].

Fluxes to small absorbers on a half-plane

We now estimate the fluxes of Brownian particles to two absorbing small holes, and when the source is located at , which is the two-dimensional half-plane (Fig. 1A). Diffusing particles are reflected everywhere on the boundary of half-space, except at the two small targets. The boundary value problem in equation (2) for two windows reduces to We set and derive a solution of equation (5) in the small window limit. We construct an inner and outer solution. The inner solution is constructed near each small window [29] by scaling the arclength s and the distance to the boundary η by and (here we use the same size ), so that the inner problem reduces to the classical two-dimensional Laplace equation The far field behavior for and for each hole is where is the flux The general solution of equation (2) with is obtained from the outer solution of the external Neumann–Green's function given for by where is the symmetric image of through the boundary axis 0z. The uniform solution is the sum of inner and outer solution (Neumann–Green's function) where are constants to be determined. To that purpose, we study the behavior of the solution near each point . In the boundary layer, we get Using this condition on each window, we obtain the two conditions: Due to the recursion property of the Brownian motion in two dimensions, there are no fluxes at infinity, thus the conservation of flux gives: In the limit of two well separated windows (), using the condition for the flux in equation (10) we get for each window (the minus sign is due to the outer normal orientation), thus Using relation (19) and (16), we finally obtain the system of two equations to solve The absorbing probabilities are given by and These probabilities precisely depend on the source position and the relative position of the two windows. When one of the splitting probabilities (either or ) is known and fixed in , recovering the position of the source requires inverting equation (23). For , the position lies on the curve At this stage, we conclude that knowing the splitting probability between two windows is not enough to recover the exact position and distance of the point source , because it leads to a one dimensional curve solution. However the direction can be obtained by simply checking which one of the two probabilities is the highest.

Fig. 1

Brownian fluxes to small windows in different geometries. (A) Two windows of size 2ϵ are placed on the boundary of half-space a distance d apart. Diffusing particles are released from a source at 0 at a distance L = |0| and are absorbed by one of the windows. (B) Two absorbing windows of size 2ϵ are placed on the circumference of a disk with radius R at angles θ1 and θ2 with the x-axes. As before, particles are released at the source position 0. (C) Two windows are placed on a disk as in (B), inside an infinitely long strip with reflecting walls at y = ±a.

Fluxes to small windows on a disk

A similar asymptotic can be derived for the splitting probability when the domain containing the windows is a disk of radius R. The boundary condition are similar: there are no particle fluxes except on and the two windows remain absorbing (Fig. 1B).

We recall that the external Neumann–Green's function of a disk of radius R, solution of the boundary value problem is given explicitly for by It is the sum of two harmonic functions with a singularity at and an image singularity at . A direct computation shows that , where . Following the derivation given for the half-plane above, we can use the Neumann–Green function of the disk (26) directly in expression (22) and obtain the probability to be absorbed on each window: for window 2 and for window 1

Splitting fluxes with many windows

The general solution of equation (2) is given by where are constants to be determined. We derive a matrix equation using the solution behavior near the center of the windows , and obtain the ensemble of conditions for The final equation is given by total flux condition: When the absorbing windows are well separated compared to the distance , a direct computation using (32) gives The ensemble of conditions (34) and (36) is equivalent to a matrix equation where for , , , for , and for , The matrix is symmetric and invertible, but does not have a specific structure, rendering it difficult to compute an explicit solution for a large number of windows in general. However, system (37) can be straightforwardly solved numerically to find the unique solution and the constant C.

Construction of a hybrid analytical-stochastic simulations of Brownian particles to small windows

We present a numerical method to efficiently simulate the splitting probability of Brownian particles arriving at small windows located on the boundary of an obstacle Ω in an infinite two-dimensional domain. The Brownian particles are generated at a single source point .

It is always possible to run naive Brownian trajectories, starting from the source, however, the mean arrival time of a Brownian particle to a target in an infinite two-dimensional domain is infinite, which would render the computational effort prohibitive. Therefore, naive simulations are inefficient, especially when computing average fluxes due the very large excursions of Brownian trajectories before they hit their targets. However the probability for any particle to hit a window is one, hence we dedicate the present section to develop a mixed stochastic simulations for computing the splitting probability.

To resolve the difficulties associated with naive Brownian simulations, we now introduce the simulation procedure. The goal of this procedure is to efficiently produce large ensembles of trajectories for estimating the splitting probability of Brownian particles generated at position and absorbed at small windows located on the surface of a two-dimensional domain. The domains are either a disk of radius R or the boundary of half-space, but any shape is possible. The procedure can be generalized to any obstacle surface in any number of dimensions, where random particles evolve in an unbounded space. We now describe the mixed algorithm consisting of two steps:

Hybrid analytical-stochastic algorithm

We present below the explicit Green's function and the steady-state flux for different geometries of Fig. 2A–C. Because the splitting probability between the windows does not depend on the arrival time, we did not take into account the temporal aspect in the renewal process, when re-injecting Brownian particles.

Fig. 2

(A) hybrid stochastic simulation procedure for two windows on the boundary of half-space. Brownian particles injected at 0 are directly place on semi-circle with radius R according to the exit pdf p (red arrow). Inside the disk, trajectories are generated by the Euler's Scheme (40) until it passes outside the radius R > R, where the trajectory is terminated at point T and restarted at a new position determined by the pdf p. (B) same as in (A) but for a ball. (C) hybrid simulation scheme for windows on a disk in a strip. Brownian particles are injected at the boundary x = d based on the exit probability distribution p. Trajectories with x > d or x < −d are re-injected at x = ±d according to p (same procedure as in (A)). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The first step consists of replacing Brownian paths by repositioning a Brownian particle to the boundary of an imaginary circle with radius (Fig. 2A–B). The position of the particle on is computed from the exit distribution of the steady-state Fokker–Planck equation with a zero absorbing boundary condition on . The exit point probability here is actually the Green's function of the Laplace operator with zero absorbing boundary condition on .

In the second step, we define a larger disk of radius and run Brownian trajectories in the domain after starting on , until they are either terminated because they escape the disk (point T in Fig. 2) or are absorbed on a small window. When a Brownian particle escapes through , we resume the trajectories by assigning a new random initial position on choosing uniformly distributed according to the exit probability .

For part (1), we need to use the explicit Green's function for the given geometry. We consider three examples: a half-disk on the boundary of the two-dimensional half-plane, the exterior of a disk in two dimensions and the two-dimensional half-strip. For part (2), we simulate particle trajectories using the Euler's scheme for the position where is a two-dimensional normal distributed vector with zero mean and variance one and D the diffusion coefficient. The time step Δt is chosen such that the mean square displacement between two time points is smaller than the size of the absorbing window ε.

Construction of the mapping using explicit Green's functions

The first part of the stochastic-analytic hybrid algorithm consists of mapping the source to a point on the artificial circumference . The construction starts with the explicit exit distribution of the Laplace operator with zero absorbing boundary condition on . The initial point is chosen randomly distributed according to the probability , where s is the arclength coordinate.

Hybrid map positioning for the full space

We start with the explicit external Neumann–Green's function in with zero absorbing boundary conditions on a disk of radius R. The solution is constructed by the method of images [22] and given by Thus the probability distribution of exit points on the boundary given that source is located at position is computed by normalizing the flux [26], The flux is computed in polar coordinates , and the angles θ and (with the horizontal axis) of points and respectively: Note that indeed . The probability (44) is used to compute the position of the sequence of points randomly and uniformly chosen, until the trajectory is finally absorbed at one of the windows. Each time a trajectory hits the external circle of radius , the motion is immediately resumed at one of the points (). This procedure disregards the absolute time of the trajectories.

Hybrid map positioning for a half-space

The Neumann–Green's function for the half-space with zero absorbing boundary condition on a half a disk of radius R is the solution of the boundary value problem The solution is obtained by the method of image charges using the Green's function for the absorbing disk in free space computed in section in eq. (3.2.1). The Green's function for the half-space is then constructed by symmetrizing with respect to the reflecting z-axis: where is the mirror reflection of on the vertical axis. The exit probability distribution is the flux through the absorbing half disk boundary where the length in polar coordinates are , and the angles θ and of and are given with respect to the horizontal axis respectively.

Green's function for the semi-strip

Finally, we summarize the Neumann–Green's function for a semi-strip of width . The normalized flux is the distribution of exit points [26]. A zero absorbing boundary condition is imposed on the boundary and a reflecting boundary condition on the rest of the strip (Fig. 2C). The function is solution of the boundary value problem The exit probability distribution is given explicitly (see appendix) by the flux through the artificial boundary

Results of the hybrid algorithm

Computing the splitting probability in half a space

To illustrate the mixed stochastic-analytical algorithm described in the previous section, we computed the splitting probability for two windows located on the y-axis, where Brownian particles are released at position in the half space (Fig. 3A). We compared the splitting probability computed analytically (formula (22)) with the results of the hybrid stochastic-analytical simulations and found perfect agreement (Fig. 3B) when we vary the distance of the source from the origin O. When the source is located more than 10 times the distance between the two windows, the difference in fluxes between the two windows decays to less than , suggesting that small fluctuations of the same order of magnitude render it impossible to measure the direction of the source from the steady state fluxes. We obtained similar results for the case of the disk, but not for a disk located in a band with reflecting walls, where the detection sensitivity extends much further [8].

Fig. 3

Diffusion fluxes to small windows on the boundary of a half-plane. (A) Particles are released at the source 0 at a distance L = |0| from the origin forming an angle θ with the x-axis. (B) Splitting probability (normalized flux) at window 2 as a function of θ for different L. The exact analytical solution given in equation (23) (lines) is compared to hybrid simulations (markers).

In order to study the sensitivity with respect to the relative positions of the windows and the source, we previously introduced the sensitivity ratio as the difference between the fluxes to the two windows [8] We then used this ratio to define the domain of sensitivity as consisting of all possible source locations that yield a ratio larger than a pre-defined threshold T. Similarly, we defined the maximum detection threshold function as [8] For a disk of radius R in free space, this maximum is indeed obtained 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. Hence, we obtain and and In particular, a Taylor expansion of for large source distance , leads to the decay of the maximum detection threshold function [8]

For two windows located on the boundary of half-space, the sensitivity ratio (51) can only be influenced by the spacing between the windows . Therefore, we do not need to find the optimal arrangement and can directly compute the sensitivity ratio where θ is the angle between the x-axis and the vector from the origin O to the source location . A Taylor expansion for of the logarithmic term yields to where the maximum of the detection threshold is similar to the one of the disk in equation (54) with and .

We conclude that for a disk and the half-plane, the detection threshold decays algebraically with the distance. In a biological context, which involve two different types of absorbing windows, each accepting only one of two types of Brownian particles, the splitting probabilities are independent. In this case, we would define the sensitivity as the product of each particle's sensitivity function by Interestingly, this formula would predict for that case a decay of the splitting probability of with respect to the source position.

Recovering the position of the source from the fluxes to several windows

To reconstruct the location of a source from the measured fluxes, at least three windows are needed. Indeed, with two windows only, a source located on the line perpendicular to the one of the connecting windows would, for example generate the same splitting probability , leading to a one dimensional curve degeneracy for the reconstructed source positions .

To study the reconstruction of a source location from the splitting probabilities, we need to invert system (34). The general solution is given by where are constants to be determined. Following the steps of section 2.3, the three absorbing boundary conditions for give The normalization condition for the fluxes is and the solution is Using the determinant: and the general notation for any we get and These equation uniquely solve the problem of determining the source location from the fluxes. Indeed, choosing and such that , the position of the source is located at the intersection of the two curves: Due to the normalization condition (17), the flux condition (relation (69)) on window 3 is redundant. Analytically solving the system (70) remains difficult, hence we investigate the position of the source numerically by inverting system (70) using expression (67)–(68). The result is shown in Fig. 4B. We conclude that for three and more receptors, it is always possible to reconstruct the source location. Here, we positioned the source 8 times the distance between the effective receptors, which could also represent clusters of receptors. Indeed, this organization would correspond to receptors located on the diameter of a round cell. Hence, we interpret the present result as recovering a source located 8 times the diameter of a cell (see [8] for an application).

Fig. 4

Reconstruction of the source position in the half-plane. (A) Two windows placed a distance d = 1 apart allow the recovery of the source position up to a curve. Three different flux configurations are shown. (B) Three windows positioned at y = −1,0,1 yield two independent curves, the intersection of which is the position of the source at 0 = (−2,8). The redundant third condition is shown for completeness. (C) The shaded areas indicate the uncertainty resulting from fluctuations in the fluxes with an amplitude of η = 0.005. The resulting sensitivity of the reconstructed source position (overlapping shaded areas 1 and 2) is highly inhomogeneous. (D) Five windows with nearest-neighbor distance d = 1 yield four independent curves. The point of intersection of all curves is the recovered source position.

We next tested the effect of possible uncertainty in the steady state fluxes on the recovery of the source , by adding a small perturbation to the fluxes, so that with in (25). Using numerical solutions, we found that the resulting uncertainty in has a highly non-linear spatial dependency, as shown by the relative sizes of the areas labelled 1 and 2 in Fig. 4C.

Finally, we studied the consequence of adding more windows. These additional windows allows to refine the reconstruction of the source. We increased the number to 5 (Fig. 4D) and indeed found that the source is precisely located at the intersection of all curves for a given set of fluxes. There are other points at which two curves intersect, however, there is only one location where more than two curves (all of them) intersect, which corresponds to the source position. We conclude that having several windows could reduce the area of the uncertainty region when the fluxes contains some steady-state fluctuations.

In order to further investigate the sensitivity of the recovered source position to small fluctuations in the flux, we numerically solved system (37) for three windows. Fig. 5A shows how the recovered distance depends the fluxes and for windows on the boundary of the half-plane, while Fig. 5B displays the same result for a disk in .

Fig. 5

Distance to the source as a function of the flux. (A) Three receptors arranged on the boundary of the half-plane, and (B) three receptors equally spaced on the circumference of a disk. The color shading indicates the logarithmically scaled distance to the source as a function of the two independent fluxes P1 and P3. The grey area indicates invalid flux combinations due to the condition ∑P = 1. The top inset displays the arrangement of windows and the lower inset shows the relationship between the fluxes and the source distance along the blue line. Note that in both cases, the distance increases nonlinearly as a function of P1/3, hence the recovery of the source position is robust against flux noise only at short distances. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Discussion and conclusion

In this paper, we studied the steady-state distribution of fluxes across absorbing windows located on a surface of a disk embedded in the plane and narrow strip and on the boundary of the half-plane. Interestingly, we found that with three absorbing windows, it is possible to recover the location of a Brownian source of particles from the splitting probability in two dimensions. The analytical computations are based on matched asymptotics to construct the Green's function used in the analysis. We developed a new mixed efficient algorithm to compute numerically the fluxes by generating truncated Brownian trajectories. Both the analytical results and the simulation procedure presented here rely on the Green's function of the domain of interest. The motivation of the present work is the problem that a cell has to face for navigation: finding a gradient source inside a tissue, which is often a two dimensional rather than three dimensional problem, as cells are moving along other cells, thereby reducing the dimensionality. In addition, many chemotaxis experiments occurring in microfluidic chambers are almost two dimensional. The methods we developed here are however also applicable without restrictions in three dimensions. In particular, the same procedure can be used to simulate sensing in three dimensional spaces.

Other methods use Green's function for the simulations of stochastic particle trajectories (reaction–diffusion method [31]), applied to transient receptor binding [17]. However, these kinetic simulation methods are quite different from the direct time step propagation method that we have developed here. The model is a molecular gradient, generated by a fixed source emitting Brownian particles. We simplified the cell geometry to a round disk containing small, fast absorbing targets on its surface (receptors). A diffusing molecule may find one of these receptors, leading to its activation. We neglected the binding time. Furthermore, we only focus on the steady-state regime in which the external gradient is already established and thus there is no intrinsic time scale. This is in contrast to a transient regime, that could represent the regime shortly after the source first starts to emit, in which the time scale can be defined as the first passage time of particles to reach the cell.

Receptor activation can mediate cellular transduction that transform a signal from the external environment into a cellular biochemical activation cascade. When a cell has to differentially compare the fluxes arriving at either of its sides sides, the local transduction of the signal at the scale of a receptor must not be homogenized throughout the rest of the cell domain, such that the local information about the gradient directionality is preserved. The internal transduced signal can be carried by the concentration of second messenger or diffusing surface molecules. Hence, receptor activations need to be localized inside the cell, leading to an asymmetrical response. Therefore, we studied the flux to stationary and localized receptors, and do not replace receptors with a homogenized boundary condition that is unable to preserve flux differences across cells and directional information.

Studies on the estimation of the fluctuations in the number of receptor-ligand molecules reaching a cell can be found in [4], [10], [3]. These models are generally based on homogenization of the boundary condition, rendering it impossible to recover any directional information since they assign the same flux to the entire boundary. In order to find the fluxes to each window separately, we based our analysis on the narrow escape theory [16]. Interestingly, we find that in two dimensions, the difference in the probability flux decays algebraically with , where L is the distance to the source. Although we focus the present investigation only on two or three windows, the results would be very similar for clusters of windows [14]. The low number of relevant windows located on a neurite must be involved in detecting a gradient concentration, its direction to turn, its forward or retracting motion [5], [27]. The model we have studied here is equivalent to fast binding [18] (without rebinding). It remains a challenge to explain how bacteria [2], sperm [28] or neurite growth [6] localize the source of a chemotactic gradient. Sensing the fluxes to receptor across the cell body is certainly the first step and the present study shows that this information is sufficient to reconstruct the location of a source. The optimal distribution of receptors could also vary from a uniform to a cluster distribution of receptors, a problem that should also be studied.

Appendix

Stability of the hybrid-simulations

We evaluated the stability of the computational method described in section 3.1 by varying the inner radius where Brownian particles are injected. The scheme of the algorithm is presented in Fig. 6A. Varying the inner radius from to 3 has no impact on the measured steady-state fluxes, as shown in Fig. 6B. This method neglects the return of far away trajectories, which in principle occurs with probability 1 in dimension 2, due to the recurrent properties of the Brownian motion.

Fig. 6

Simulation results are independent of the choice of entering radius R and outer radius R. (A) Two windows are placed on a disk. Brownian particles originating from the source at x0 are injected at the circumference of the circle with radius R. Trajectories leaving the region with radius R are restarted at a radius R. (B) Varying the exit radius R and R does not change the splitting probability (flux) from simulations. The box plot shows the distribution of errors in the flux to window 2. The error is defined as the deviation of the flux from the mean over all observed radii combinations (R,R) and all angles θ2, but separately for all source distances L. Note that there are no systematic deviations when changing either R or R.

Explicit Green's function in a band

The hybrid algorithm is based on the exact expression of the Neumann–Green's function for the semi-strip where . The normalized flux is the distribution of exit points [26]. We impose zero absorbing boundary condition on the boundary and reflecting boundary condition on the rest of the strip . The boundary value problem is We compute by expanding in eigenfunctions following the classical method of [22] (p. 80). To start, we write the following Ansatz with , the components of and , the components of , which solves the inhomogeneous diffusion equation Because is bounded to between 0 and a, we can write u and f in terms of a Fourier series along with By inserting this expression for u into equation (73), we arrive at the following ODE for the For , the fundamental solutions to the homogeneous equation are and . Therefore, the inhomogeneous problem is solved by where due to the boundedness condition on the solution as and due to the absorbing boundary at .

For , the fundamental solutions to the homogeneous system are given by . Hence, Thus, the complete solution to equation (73) reads where we inserted the Fourier coefficients for f from equation (75). By inspection, we arrive at the expression for the Green's function Using the identity (see [22] p. 84) we can further simplify to get with . The exit probability distribution is again given by the flux through the boundary