A Geometrical Method for the Smoluchowski Equation on the Sphere

The diffusion of a tracer particle on the surface of a two-dimensional sphere has served as a simple model to describe many processes observed in nature. On the one hand, such a tracer might correspond to the tip of a vector performing random rotations, where the vector may represent the spin of a molecule, the axis of a rigid body or the direction of motion of a swimming bacteria. On the other, it describes the random motion of particles on the surface of a sphere. In this context, relevant to biophysics, it approximately describes the motion of certain proteins or phospholipids (PIP2 Fujiwara et al. (2002); Metzler et al. (2016)), on the cell membrane, which are crucial messengers in cell signaling Haastert and Devrotes (2004); Ma et al. (2004); van Meer et al. (2008); Meyers et al. (2006); Hancock (2010).

While the mathematical framework that approaches the isotropic rotations of a rigid-body orientation has been developed Furry (1957); Favro (1960); Ivanov (1964); Hubbard (1972); Valiev and Ivanov (1973); McClung (1980), the Brownian motion of a particle on the surface of a sphere has been analyzed as an extension to general manifolds Graham (1977); van Kampen (1986); Risken (1989) of the theory of translational Brownian motion in Euclidean spaces and, at the turn of the 21st century, the improvement on single-particle tracking methods triggered a resurgence of the analysis the random motion of biological tracers on curved surfaces Faraudo (2002).

Advances on this field have been made for the case when the underlying surface is (nearly) spherical in shape (for instance, the pseudopod formation in cell membranes can be inhibited with a specific hormone Janetopoulos et al. (2004)). In this case, two mathematical frameworks have been considered to take into consideration the holonomic constraint of moving on the sphere’s surface, namely, the Fokker-Planck equation and the Langevin equation, which have also become the theoretical frameworks for studying systems far from thermal equilibrium such as active matter Castro-Villarreal and Sevilla (2018), or reaction-diffusion systems Garduño et al. (2019).

The free and isotropic motion (without obstacles or external forces) of a particle diffusing on the surface of the two dimensional sphere is well-known, however, the necessity to devise numerical algorithms to generate trajectories on the sphere to investigate the transport processes associated has recently been pointed out Castro-Villarreal et al. (2014); Yang and Li (2019); Novikov et al. (2020). A much less investigated case corresponds to that one when the motion is under the influence of external forces. This would be the situation, for instance, when the orientational dynamics of the swimming direction of a bacteria are subject to chemotaxis Schnitzer (1993) or when the rotational Brownian motion is anisotropic Ford et al. (1979).

The general framework to describe Brownian motion constrained to the two dimensional sphere of radius $r$, denoted by $\mathbb{S}^{2}_{r}$, that is, under a holonomic constraint, is based on the Fokker-Planck equation and has been presented in Refs. N.G.Van.Kampen.1986; Risken (1989). In few words, one starts from a formulation of the Fokker-Planck equation in Cartesian coordinates $x^{k}$ ($k=1,2,3$). This is transformed by choosing an appropriate coordinates system (see Appendix A), that allows for the implementation of the constraint. In Cartesian coordinates, the probability density $P(x^{k},t)$ of finding a particle at the coordinates $x^{k}$ at time $t$ satisfies the Fokker-Planck equation

where $A^{k}$ plays the role of the $k$-th component of a velocity field generated by a force field in the overdamped limit and $B^{kl}=D\,\delta^{kl}$ denotes the diffusion tensor with $D$ the diffusion constant. $[\cdot],_{k}$ and $[\cdot],_{kl}$ denote $\frac{\partial}{\partial x^{k}}[\cdot]$ and $\frac{\partial^{2}}{\partial x^{k}x^{l}}[\cdot]$, respectively, and Einstein’s summation convention must be understood on repeated indexes. The particle motion is constrained to $\mathbb{S}^{2}_{r}$ of radius $r$ by requiring that at each instant $x^{k}x_{k}=r^{2}$. Such a constraint is naturally implemented by the use of spherical coordinates: $(x,y,z)\to(r\sin{\theta}\cos{\phi},r\sin{\theta}\sin{\phi},r\cos{\theta})$, which induces the transformation $P(x^{k},t)dx^{1}dx^{2}dx^{3}=r^{2}\sin\theta P(\theta,\phi,t)drd\theta d\phi$ and under the constriction we get $P(x^{k},t)\delta\bigl{(}\sqrt{x^{k}x_{k}}-r\bigr{)}\,dxdydz=\,r^{2}\sin\theta\,P(\theta,\phi,t)d\theta d\phi$. By defining $\overline{P}(\theta,\phi,t)=r^{2}\sin\theta\,P(\theta,\phi,t)$, this satisfies

with $A^{\theta}=r^{-1}[\cos\theta\cos\phi\,A^{x}+\cos\theta\sin\phi\,A^{y}-\sin\theta\,A^{z}]$, $A^{\phi}=(r\sin\theta)^{-1}[-\sin\phi\,A^{x}+\cos\phi A^{y}]$. The ‘spurious’ or ‘noise-iduced’ drift, $(D/r^{2})\,\cot{\theta}\,\,\overline{P}$ Ryter and Deker (1980), appears not only as a consequence of the coordinates transformation, but as part of the total drift, whose effects are physically observable (notice that is proportional to the diffusion constant). For instance, it modifies the ‘hopping’ rate over the potential barriers depending on the particular latitude on the sphere. Indeed, above the equator, the particles are affected by an extra ‘push’ toward the south pole, proportional to the diffusion coefficient. Similarly, below the equator, the extra ‘push’ is toward the north pole. So, curvature will manifest (indirectly) in these kind of physical measurements; for example, if a chemical reaction depends on a potential energy function, the reaction rate will depend on the curvature of the ambient space in which the reaction is taking place.

Exact analytical solutions for (2) are seldom available, except for the case of free diffusion ($A^{\theta}=A^{\phi}=0$), this is why it is important to count on validated numerical methods either to solve it explicitly or implicitly by simulations of such solutions. This becomes relevant especially in the context of applications.

In this work, we devise a numerical algorithm to generate ensemble trajectories of Brownian particles under the effects of an external field, which properly samples the conditional probability density $\bar{P}(\theta,\phi,t|\theta_{0},\phi_{0},t_{0})$, solution of equation (2). We validate the algorithm by comparing its predictions against analytical results in two different scenarios: one considering the time evolution, the other in the stationary regime. In contrast to the algorithm presented in Carlsson et al. (2010), our algorithm is based on the projection over the spherical surface, after a proper rescale of the Brownian step computed on the tangent plane. These two steps, rescaling and projecting, implement in an exact manner the constriction that maintains the particle motion on the sphere surface, and describe equivalently the diffusion process that underlies equation (2). Furthermore, there are two main advantages of this method. Firstly, we avoid the implementation of rotations so we avoid as many trigonometric evaluations as possible; secondly, there is no need of prior knowledge of Riemannian geometry or stochastic calculus to understand its essence, but only elementary geometry.

The paper is organized as follows: in section II we derive the numerical algorithm that we designed to tackle the problems just defined mathematically in this section. In section III we present our numerical results as well as the analytical solutions to free diffusion and the solutions to the stationary distributions when there is an external interaction, putting emphasis in the differences with respect to the solutions of flat spaces. We explicitly derive from the associated Fokker-Planck equations the analytic solution to the stationary distributions in the appendix B. In section IV we present our concluding remarks and we leave to appendix A within the appendices, a brief deduction of the transformation rules for the elements of the Fokker-Planck equation that justifies equation (2).

We present a numerical algorithm that generates an ensemble of trajectories of non-interacting Brownian particles that diffuse on the surface of the sphere $\mathbb{S}^{2}_{r}$, in the overdamped limit, and under the effects of an arbitrary external potential. The algorithm is implemented in three dimensions without making use of the standard, but sometimes singular, Langevin-like evolution equations of the spherical angles $\theta$ and $\phi$.

The instantaneous particle’s position at time $t$, with respect to the sphere center, is denoted with $\boldsymbol{r}(t)=r\hat{\boldsymbol{n}}(t)$ with $\hat{\boldsymbol{n}}=\bigl{(}\sin\theta(t)\cos\phi(t),\sin\theta(t)\sin\phi(t),\cos\theta(t)\bigr{)}$ $\in\mathbb{R}^{3}$. At this point on the surface, the tangent plane $\text{T}_{\boldsymbol{r}(t)}\mathbb{S}^{2}$ is used to approximate the particle position after a time interval $dt$ namely, $\boldsymbol{r}(t+dt)=\boldsymbol{r}(t)+d\boldsymbol{r}(t)$, i.e. the new approximated particle position is computed locally by calculating the net change $d\boldsymbol{r}(t)=d\boldsymbol{r}_{n}(t)+d\boldsymbol{r}_{f}(t)$, where $d\boldsymbol{r}_{n}(t)$ is the change due to the effects of noise, and $d\boldsymbol{r}_{f}(t)$ the change caused by all of the forces on the particle as is shown schematically in figure 1.

The noisy term $d\boldsymbol{r}_{n}(t)$ is computed by generating two pseudo-random numbers: one normally distributed using the Mersenne Twister method, and the other uniformly distributed in $[0,2\pi]$. In the finite time interval $\Delta t$ we have that

where $\mathcal{A}^{\theta}$ and $\mathcal{A}^{\phi}$ are the components of the deterministic forces along the directions given by the unit vectors $\hat{\boldsymbol{\theta}}_{\boldsymbol{r}}(t)=\bigl{(}\cos\theta(t)\cos\phi(t),\cos\theta(t)\sin\phi(t),-\sin\theta(t)\bigr{)}$ and $\hat{\boldsymbol{\phi}}_{\boldsymbol{r}}(t)$ $=\bigl{(}-\sin\phi(t),\cos\phi(t),0\bigr{)}$ at position $\boldsymbol{r}(t)$, respectively. $\zeta$ in equation (3b) denotes the damping constant that emerges from the coupling between the particle motion and the thermal bath, which gives rise to the fluctuating motion of the particle.

Since $\boldsymbol{r}(t+\Delta t)=\boldsymbol{r}(t)+d\boldsymbol{r}$ does not lie on $\mathbb{S}^{2}_{r}$, but on $\text{T}_{\boldsymbol{r}(t)}\mathbb{S}^{2}$, it must to be projected back to the sphere surface (see figure 1). This can be accomplished by the implementation of the following algebraic scheme that allows an accurate calculation of the particle displacement on the sphere. First, once $d\boldsymbol{r}$ is computed, it is scaled by a factor $\alpha$ in such a way that the arc length of the geodesic $\gamma$ between $\boldsymbol{r}(t)$ and the projection of $\boldsymbol{r}(t)+\alpha d\boldsymbol{r}$ on the sphere (see figure 1b), denoted with $\boldsymbol{r}_{{}_{\Pi}}(t+\Delta t)$, coincides with $||d\boldsymbol{r}||$. From simple geometry it can be shown that $\alpha=r\tan(||d\boldsymbol{r}/r||)/||d\boldsymbol{r}||$. Thus, the particle position is updated according to

which forms the basis of our numerical algorithm that takes into account the effects of the external forces. In (4), the symbol $\prod_{\mathbb{S}^{2}}(\mathbf{x})$ denotes the projection operator $\boldsymbol{x}(\boldsymbol{x},\boldsymbol{x})^{-1/2}$, with $(\,,\,)$ the natural scalar product in $\mathbb{R}^{3}$, used to bring back the particles from $\mbox{T}_{\boldsymbol{x}}\mathbb{S}^{2}\subset\mathbb{R}^{3}$ onto $\mathbb{S}^{2}$; and $\hat{d\boldsymbol{r}}=d\boldsymbol{r}/||d\boldsymbol{r}||$. In the absence of noise, equation (4) reduces to the integration of the deterministic part of the dynamics, i.e.

The method is of first order in $\Delta t$ in the context of stochastic differential equations, however it could be possible to develop an analogous form of the the Miltein’s method Higham (2001), which is based on an Itô-Taylor expansion. Some of these topics are discussed in Gardiner (1997).

In the absence of external forces, the vector

is the one that takes into consideration the particle’s Brownian motion on the sphere, where $d\boldsymbol{r}_{n}$ is given in equation (3a). As such, and as our numerical analysis shows afterwards, the underlying geometry of this updating rule sheds light over the geometric meaning of the three terms proportional to $D$ in the Smoluchowski equation (2), since these same terms fully describe the diffusive aspects of Brownian motion on the sphere. We want to emphasize that our algebraic scheme, the combination of scaling and projecting after updating in the tangent plane, is equivalent to the exact updating rotation of $\boldsymbol{r}(t)$ and therefore no additional systematic error is introduced in such procedure (see SupMat (2021)).

In this section, we compare the results obtained from an statistical analysis based on an ensemble of trajectories that start at $\theta(t=0)=0$ (north pole) computed from the algorithm described in the last section. Firstly, we perform a comparison in the context of free diffusion on $\mathbb{S}^{2}_{r}$, for which an analytic solution for the marginal probability distribution $\overline{P}(\theta,t|0,0)=2\pi\overline{P}(\theta,\phi,t|0,0,0)$ is available at all times, namely Roberts and Ursell (1960); Hobson (1965); Asmar (2005)

where the spherical addition theorem Jackson (1999) is used to write the sum of the product of spherical harmonics as a Legendre polynomial $P_{n}(\cos\theta)$. We focus on the standard parameters that characterize the particle diffusion, namely, the mean square displacement and the position autocorrelation function $\langle\hat{\boldsymbol{n}}(t)\cdot\hat{\boldsymbol{n}}(0)\rangle=\langle\cos\theta(t)\rangle=\langle P_{1}(\theta)\rangle$. Furthermore, we calculate the complete histogram of the particle positions and the mean value of the polar angle, $\langle\theta\rangle$.

Although the quantities of interest have an analytical expression obtained from the solution (7) of the Fokker-Planck equation (2) in the absence of external fields, only the position autocorrelation function can be written in a closed form, namely

For the calculation of the mean polar angle $\langle\theta(t)\rangle=\int_{0}^{\theta}d\theta\,\theta\overline{P}(\theta,t|0,0)$ and the mean square displacement $\langle\boldsymbol{r}^{2}(t)\rangle$ (which under the initial conditions coincides with $\langle\theta^{2}(t)\rangle=\int_{0}^{\theta}d\theta\,\theta^{2}\overline{P}(\theta,t|0,0)$) we used a numerical evaluation of the involved integrals using (7).

Afterwards, we use our numerical algorithm to analyze the of diffusion of particles on $\mathbb{S}^{2}$ in the presence of an external field which in general can be expanded in the spherical harmonics $Y_{l}^{m}(\theta,\phi)$. For the sake of a clear analysis, we consider an external potential $U_{\lambda}$ that depends only on $\cos\theta$ (or equivalently that depends only on the quantity $z/(x^{2}+y^{2}+z^{2})^{1/2}$), thus focusing on the cases with azimuthal symmetry. $\lambda$ denotes the external potential strength whose physical dimensions are of energy. In such a case we have that the components of the force field in spherical coordinates are given by

That being so, $U_{\lambda}(\theta)$ can be expanded in the Legendre polynomials $P_{l}(\cos\theta)$. We chose three specific cases for $U_{\lambda}(\theta)$, namely

which are depicted in figure 2.

The (minus) sign in (10b) and (10c) is chosen in order to have a potential with local minima at $\theta=0,\,\pi$ and one local maximum at $\theta=\pi/2$ in the first case. For the second case, we have a local minimum at $\theta=0$ and a maximum at $\theta=\pi$ separated by a local minimum at $\theta=\arccos\left(1/\sqrt{5}\right)$.

As has been pointed out before, there are no available analytical solutions to the time dependent problem, but there are for the stationary solutions. The solution corresponding to (10c) cannot be expressed in terms of elementary functions, so we will evaluate it numerically. We use as the characteristic length the radius of the sphere, such that when combined with the diffusion coefficient it defines the time scale $\tau=r^{2}/D$. Except when is stated otherwise, we have used $D=0.1$, $\lambda\approx 2.24$, and $\Delta t=\ln{2}\times 10^{-3}\tau$ for the time step.

In this paper, we presented a numerical algorithm to generate the trajectories of Brownian particles diffusing on the surface of the two-dimensional sphere. The algorithm implemented is based on three-dimensional geometry, thus avoiding the complications introduced by the interpretation of the multiplicative Gaussian-white noise involved in the corresponding stochastic differential equation, and takes into account the effects of an external potential. Our numerical method is weakley convergent, that is, it statistically converges to the solution of the Smoluchowski equation (2).

The algorithm was validated in two scenarios: First, in the absence of external potential, by computing the time dependence histogram of the polar angle, as well as the first two moments, and the time correlation function of the particle positions, from an ensemble of trajectories generated by our algorithm. We observed an excellent agreement at all times between our results and the corresponding quantities obtained from the solution to the diffusion equation on the sphere (7). Second, in the presence of a external potential, for which we calculated the stationary probability distribution when the external potential depends only on $\cos{\theta}$. Again, excellent agreement was found between our numerical results and those given by the stationary solutions of the Smoluchoswki equation (2) for three different potentials. In addition, we have also calculated the time correlation of the particle positions for each of these three potentials. Although no analytical solution are available to our knowledge, the traits of correctness are observed. Although we have considered simple potentials, our algorithm can be applied to an arbitrary but physically motivated potential.

An alternative to our numerical algorithm might consider either the numerical integration of the differential stochastic equations that determine the time evolution of $\theta(t),\varphi(t)$ or the numerical solution of the Smoluchowski equation (2). However, in any case, this procedure requires to deal with singular terms that diverge at the poles. Although this can be solved by changing from one set of coordinates (chart) to another one, this is however, computationally time consuming.