# The most probable dynamics of receptor-ligand binding on cell membrane

Xi Chen1, Hui Wang2,∗, Jinqiao Duan3 1 School of Statistics, Xi ’an University
of Finance and Economics, Xi’an 710100,China  2School of Mathematics and
Statistics, Zhengzhou University, Zhengzhou 450001, China  3 Department of
Applied Mathematics, Illinois Institute of Technology, Chicago, IL 60616, USA
<EMAIL_ADDRESS>
###### Abstract

We devise a method for predicting certain receptor-ligand binding behaviors,
based on stochastic dynamical modelling. We consider the dynamics of a
receptor binding to a ligand on the cell membrane, where the receptor and
ligand perform different motions and are thus modeled by stochastic
differential equations with Gaussian noise or non-Gaussian noise. We use
neural networks based on Onsager-Machlup function to compute the probability
$P_{1}$ of the unbounded receptor diffusing to the cell membrane. Meanwhile,
we compute the probability $P_{2}$ of extracellular ligand arriving at the
cell membrane by solving the associated Fokker-Planck equation. Then, we could
predict the most probable binding probability by combining $P_{1}$ and
$P_{2}$. In this way, we conclude with some indication about where the ligand
will most probably encounter the receptor, contributing to better
understanding of cell’s response to external stimuli and communication with
other cells.

  * February 2023

Keywords  binding probability; most probable transition pathway; Onsager-
Machlup function; Fokker-Planck equation

##  1 Introduction

There has been an explosion of interest in the effect of noise in cell biology
in recent years [1, 2, 3, 4, 5, 6]. Stochastic cell biology has attracted the
attention of many researchers in biophysics, statistical physics and
stochastic dynamics [7]. As the cell is an open system that exchanges energy
and matter with the environment, the molecules inside and outside the cell
keep moving by different transport mechanisms. Besides, the cells are in
highly complex and noisy environment, so the motions are stochastic [1].

Many essential cell functions, such as monitoring the environment, responding
to external stimuli and communicating with neighboring cells, are realized by
the cooperation of intracellular and extracellular matters. For the molecules
inside the cell, there are two main transporting forms: self-diffusion and
directional transport. Self-diffusion is passive and depends on the cytoplasm;
while directional transport is active, which can quickly realize material
transport and signal conduction. It is worth mentioning that the receptor and
ligand are an important class of material that transmits information, which
regulates basic tasks essential to life, such as in cell growth, secretion,
contraction, motility, and adhesion and receptor-mediated cell behaviors [4,
8, 9, 10, 11, 12, 13].

The biochemical receptors are proteins which bind the signal molecules and
initiate responses in the target cell. The binding site of the receptor has a
complex structure that is shaped to recognize the signal molecule with high
specificity, helping to ensure that the receptor responds only to the
appropriate signal and not to the many other signaling molecules surrounding
the cell. In this case, the signal molecules are called ligands paired to the
specific receptors. When the ligand binds to the receptor, they become
activated and generate various intracellular signals that alter the behavior
of the cell [14]. As is known, receptor-ligand interactions play a major role
in environmental sensing, signal transduction, and cell-to-cell signaling [15,
16, 17, 18]. Although receptors and ligands have many different forms, they
come in closely matched pairs, with a receptor recognizing just one specific
ligand, and a ligand binding to just one target receptor.

In many cases, receptors are transmembrane proteins on the target cell
surface, while in other cases, the receptor proteins are inside the target
cell which are either located in the cytoplasm and nucleus or in the membrane
of subcellular compartments, such as the endoplasmic reticulum [4]. Receptors
are uniquely able to direct cell behavior by virtue of their ability to sense
the environment, through binding of ligands, and their ability to transmit the
signal to the cell interior, through interaction of their intracellular
domains with enzymes and proteins within the cell. In this paper, we will
focus on the unbounded receptors which bind to ligands on cell menbrane.

An additional influence on receptor-ligand binding on the cell surface is
receptor mobility, or diffusion within the cell. After the recognition of a
receptor on cell surface, the receptor-ligand complex will go into the
endocytic cycle in the cell [19], which governs the recycling to the cell
surface (receptors) and delivery to lysosomes (ligands). Besides, the
receptors are protein molecules which synthesize intracellularly and transport
to the cell surface. It is the combined work that makes the change of
receptors numbers on the cell surface, which undoubtedly affects the cellular
responses to outside signals. As the intracellular transport processes can
alter the number of receptors present on the cell surface, it is necessary to
take the motion of receptors into consideration when examining the binding
dynamics. The intracellular receptors are a class of receptors that locate in
the cytoplasm and nucleus or in the membrane of subcellular compartments. They
can pick up signals from outside the cell, bind to ligands that can be paired,
trigger reactions inside the cell, and produce specific effects [1, 3].
Diffusion effects play a role in ligand binding to receptors [20, 21]. Here we
consider a reduced model by treating a cell as a disk, and the illustration
can be seen in Figure 1.

Figure 1: (Color online) The illustration of receptor and ligand binding on
cell membrane.

Dynamical systems arising in biophysics are often subject to random
fluctuations. The noisy fluctuations may be Gaussian or non-Gaussian, which
are modeled by Brownian motion or $\alpha-$stable Lévy motion, respectively.
There are various studies about Gaussian noise in receptors and ligands [22].
However, as the ligand search for the surface receptor in the form of jump-
diffusion [23], which can not be modeled by purely diffusion any more, using
the Lévy motion with non-Gaussian distribution to depict its motion will be
more appropriate. Indeed? many researchers have begun to consider non-Gaussian
noise in cell biology [1, 25]. We will study the dynamics of ligands under
non-Gaussian fluctuations. It is worth mentioning that noise-induced
transition phenomena were regarded in various dynamical systems, including
stochastic gene expression [26, 27], climatology of the earth [28], chemical
reacting systems [29]. In particular, the most probable transition pathway
attracts a lot of attention [27, 30, 31, 32, 33].

There are several approaches for distinguishing the receptor and ligand
binding events, including performing experiments at lowered temperatures and
utilizing pharmacologic agents. However, the existing studies mostly consider
the binding of ligand and receptor when the receptor is confined to the cell
membrane [23, 34], the dynamical behaviors of the receptor and ligand
themselves are not considered. In this paper, we apply tools from stochastic
dynamical systems to investigate receptor-ligand binding, where receptors and
ligands are randomly moving. We will examine the possibility of the receptor
and ligand binding on the membrane, by calculating the joint probability of
the receptor occurring on the membrane and ligand locating at the same site.
Simultaneously, as a receptor can appear at any location on the membrane, we
will focus on the most probable trajectory of a receptor and find its
probability arrival at the membrane. For the ligand, we examine the
probability of the ligand reaching the domain close to the receptor on the
cell membrane by means of Fokker-Planck equation. Then, we relate the
probability of the receptor diffusing to the cell membrane with the
probability of the ligand arriving at the cell membrane to calculate the most
probable binding probability.

This paper is organized as follows. In Section 2, we present the mathematical
model and compute the probability of the receptor diffusing to the cell
membrane and the probability of the ligand arriving at the cell membrane
respectively. In Section 3, we calculate the most probable binding
probability. Finally, we end with some discussions in Section 4.

##  2 Binding of a ligand to a receptor on cell membrane

The binding event can be regarded as the encounter of a ligand and a receptor
on the cell membrane. To investigate the binding dynamics, we consider the
intracellular motion of the receptor and the extracellular motion of the
ligand respectively. As the receptors are large protein molecules diffusing in
cell in the form of vesicles and the intracellular environment is crowded, the
motion of receptors can be treated as Brownian motion on a certain time scale
[35, 36, 37].

###  2.1 The probability of a cell-surface receptor diffusing to the cell
membrane

The intracellular motion of the receptor can be described by a Langevin
equation with Brownian motion [1]. By assuming the cell to be a 2-dimensional
disk, we can rewrite the motion equation of receptor in the form:

$\begin{cases}dx=b_{1}(x,y)dt+\sqrt{2D_{1}}dB_{1}(t),\\\
dy=b_{2}(x,y)dt+\sqrt{2D_{2}}dB_{2}(t),\end{cases}$ (1)

where $b_{1}(x,y),b_{2}(x,y)$ represent the external force, $D_{1}$ and$D_{2}$
are diffusion coefficients with respect to the Brownian motion, and
$B_{t}=(B_{1}(t),B_{2}(t))$. Transforming Eq.(1) from Cartesian coordinates to
polar coordinates by setting $x=r\cos\theta$ and $y=r\sin\theta$, we have

$\begin{cases}dx=-r\sin\theta d\theta+\cos\theta dr,\\\ dy=r\cos\theta
d\theta+\sin\theta dr.\end{cases}$

Thus Eq.(1) can be rewritten as

${\left(\begin{array}[]{cc}\cos\theta&-r\sin\theta\\\
\sin\theta&r\cos\theta\end{array}\right)}{\left(\begin{array}[]{c}dr\\\
d\theta\end{array}\right)}={\left(\begin{array}[]{c}b_{1}(x,y)dt\\\
b_{2}(x,y)dt\end{array}\right)}+{\left(\begin{array}[]{cc}\sqrt{2D_{1}}&0\\\
0&\sqrt{2D_{2}}\end{array}\right)}{\left(\begin{array}[]{c}dB_{1}(t)\\\
dB_{2}(t)\end{array}\right),}$ (2)

which leads to an evolution equation about $r$ and $\theta$ read as

$\left(\begin{array}[]{c}dr\\\
d\theta\end{array}\right)=B\left(\begin{array}[]{c}b_{1}(x,y)dt\\\
b_{2}(x,y)dt\end{array}\right)+B\left(\begin{array}[]{cc}\sqrt{2D_{1}}&0\\\
0&\sqrt{2D_{2}}\end{array}\right)\left(\begin{array}[]{c}dB_{1}(t)\\\
dB_{2}(t)\end{array}\right),$ (3)

with $B$ denoting the inverse of the coefficient matrix in the left hand sides
of Eq.(2), which has forms

$B={\left(\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta/r&\cos\theta/r\end{array}\right).}$

Eq.(3) gives us the Langevin equation in the polar coordinates. For
simplicity, we further assume the intracellular force $f$ in polar coordinates
to be $\mathbf{r}=(r,\theta)\triangleq f$ [38]. The explicit choice of $f$
above is only to simplify many technicalities; the method used will apply to
more general but intracellular force. So the Langevin equation Eq.(3) can be
reduced as

$\left\\{\begin{aligned} &dr=rdt+\sqrt{2D_{1}}\cos\theta
dB_{1}(t)+\sqrt{2D_{2}}\sin\theta dB_{2}(t),\\\ &d\theta=\theta
dt-\sqrt{2D_{1}}\sin\theta/rdB_{1}(t)+\sqrt{2D_{2}}\cos\theta/rdB_{2}(t).\end{aligned}\right.$
(4)

In the practical biological processes, the receptor can locate many possible
sites on the cell membrane. This indicates there could exist a mass of
pathways for the receptor transiting to these sites. We are interested in the
most probable transition pathways among all the possible trajectories. Based
on Eq.(4), we will get the most probable pathway for a receptor starting from
the center$(r,\theta)=(0,0)$ to the membrane $(r,\theta)=(1,\hat{\theta})$.
The probability of receptor occurring at site $\theta$ on membrane can be
obtained by the Onsager-Machlup (OM) action functional [40, 41]

$S(z,\dot{z})=\frac{1}{2}\int_{0}^{T}[\dot{z}-f(z)][\sigma\cdot\sigma^{T}(z)]^{-1}[\dot{z}-f(z)]+\bigtriangledown_{x}\cdot
f(z)dt,$ (5)

with the Lagrangian

$\begin{split}L(z,\dot{z})&=\frac{1}{2}[\dot{z}-f(z)][\sigma\cdot\sigma^{T}(z)]^{-1}[\dot{z}-f(z)]+\frac{1}{2}\bigtriangledown_{x}\cdot
f(z)\\\
&=\frac{1}{2}|D^{-1}[\dot{z}-f(z)]|^{2}+\frac{1}{2}\bigtriangledown_{x}\cdot
f(z),\end{split}$ (6)

where $D$ denoting the matrix in Eq.(3)

$D=B\left(\begin{array}[]{cc}\sqrt{2D_{1}}&0\\\
0&\sqrt{2D_{2}}\end{array}\right)=\left(\begin{array}[]{cc}\sqrt{2D_{1}}\cos\theta&\sqrt{2D_{2}}\sin\theta\\\
-\sqrt{2D_{1}}\sin\theta/r&\sqrt{2D_{2}}\cos\theta/r\end{array}\right).$

Thus the Lagrangian becomes [42]

$L(z,\dot{z})=\frac{1}{2}(\frac{\cos^{2}\theta}{2D_{1}}+\frac{\sin^{2}\theta}{2D_{2}})(r-\dot{r})^{2}+\frac{1}{2}(\frac{r^{2}\sin^{2}\theta}{2D_{1}}+\frac{r^{2}cos^{2}\theta}{2D_{2}})(\theta-\dot{\theta})^{2}+1.$
(7)

According to the Euler-Lagrange equation

$\frac{d}{dt}\frac{\partial}{\partial\dot{z}}L(z,\dot{z})=\frac{\partial}{\partial
z}L(z,\dot{z}),$ (8)

we get the following second order differential equations for $r$ and $\theta$

$\left\\{\begin{aligned}
&\ddot{r}=r+r\theta^{2}(\frac{\sin^{2}\theta}{D_{1}}+\frac{\cos^{2}\theta}{D_{2}})/(\frac{\cos^{2}\theta}{D_{1}}+\frac{\sin^{2}\theta}{D_{2}}),\\\
&\ddot{\theta}=\theta+\sin\theta\cos\theta(\frac{1}{D_{2}}-\frac{1}{D_{1}})(1-\theta^{2})/(\frac{\sin^{2}\theta}{D_{1}}+\frac{\cos^{2}\theta}{D_{2}}).\end{aligned}\right.$
(9)

This is rewritten as the following four dimensional system

$\left\\{\begin{aligned} &\dot{r}=w,\\\ &\dot{\theta}=v,\\\
&\dot{w}=r+r\theta^{2}(\frac{\sin^{2}\theta}{D_{1}}+\frac{\cos^{2}\theta}{D_{2}})/(\frac{\cos^{2}\theta}{D_{1}}+\frac{\sin^{2}\theta}{D_{2}}),\\\
&\dot{v}=\theta+\sin\theta\cos\theta(\frac{1}{D_{2}}-\frac{1}{D_{1}})(1-\theta^{2})/(\frac{\sin^{2}\theta}{D_{1}}+\frac{\cos^{2}\theta}{D_{2}}).\end{aligned}\right.$
(10)

As we aim at getting the most probable transition pathway in a fixed period of
time and its corresponding Onsager-Machlup functional action value
$S(z,\dot{z})$, as well as the probability, the next step will be reduced to
solve a two-point boundary problem Eq.(10) with initial location
$(r(0),\theta(0))=(0,0)$ and final location
$(r(T),\theta(T))=(1,\hat{\theta})$ where $\hat{\theta}\in(0,2\pi)$. Here $T$
is the transition time for a receptor from the cytoplasm to cell membrane
which we consider as a fixed period of time and can be set $T=5$.

Assuming that the receptor starting inside the cell can reach any site on the
membrane, we divide the membrane (the circle) into 8 equal segments:
$(0,\frac{\pi}{4}),(\frac{\pi}{4},\frac{\pi}{2}),...$, and find the most
probable transition pathway on each segment as well as the minimal action
value. The key lies in solving Eq.(10), but the only known conditions are the
boundary value so for with the initial velocity still unknown. So here we
adopt a neural network approach to handle the tricky problem. The main idea is
to construct a neural network and approximate the mapping from the final point
to the initial velocity, such that the solution of Eq.(10) with this initial
velocity at time $T$ meets the target given final point. This method has been
used in solving differential equations with boundary problems which is called
neural shooting [43]. Generally, the algorithm can be summarized into the
following steps:

Step 1: Generate the data sets to supply the training basis for this neural
network, by simulating the Eq.(10) given a series of initial velocities with
setting the noisy intensity parameters $D_{1}=D_{2}=0.01$. Then we get a set
of final points at time $T$ which corresponds to a set of initial velocities.
This dataset consists of final point-velocity pairs.

Step 2: Adopt the neural network “2-8-16-8-2” constructed by Tensorflow to
approximate the mapping from the initial velocity to the final point at time
$T$ by using the dataset generated in Step 1.

Step 3: Use the learned model in Step 2 to predict the initial velocities
given the possible final points on the surface, and insert the initial
velocity into the system to compute the transition pathway as well as its
corresponding Onsager-Machlup action functional value.

The Onsager-Machlup action functional value of each pathway on the divided
segments is calculated in Table 1-8 and the minimum action functional value is
selected with bold font, that is the most probable transition pathway we are
looking for.

Table 1:  The action functional value of pathways terminated on $(0,\pi/4)$
with predicted initial velocity $(w_{0},v_{0})$

.  $\hat{\theta}$ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 $\pi/4$ $w_{0}$ 0.1030 0.1024
0.1016 0.1005 0.0991 0.0974 0.0956 0.0935 $v_{0}$ 0.0096 0.0196 0.0297 0.0398
0.0498 0.0599 0.0700 0.0801 $S$ 3.99 3.92 3.88 3.89 3.97 4.16 4.65 7.12

Table 2:  The action functional value of pathways terminated on
$(\pi/4,\pi/2)$ with predicted initial velocity $(w_{0},v_{0})$

.  $\hat{\theta}$ 0.9 1.0 1.1 1.2 1.3 1.4 1.5 $\pi/2$ $w_{0}$ 0.0882 0.0855
0.0827 0.0797 0.0766 0.0733 0.0700 0.0666 $v_{0}$ 0.0934 0.1034 0.1134 0.1234
0.1334 0.1434 0.1534 0.1634 $S$ 28.83 348.86 1568.5 498.76 99.31 63.05 216.42
123.02

Table 3:  The action functional value of pathways terminated on $(\pi/2,3\pi/4)$ with predicted initial velocity $(w_{0},v_{0})$. $\hat{\theta}$ | 1.7 | 1.8 | 1.9 | 2.0 | 2.1 | 2.2 | 2.3 | $3\pi/4$
---|---|---|---|---|---|---|---|---
$w_{0}$ | 0.0617 | 0.0582 | 0.0548 | 0.0513 | 0.0478 | 0.0444 | 0.0411 | 0.0378
$v_{0}$ | 0.1711 | 0.1814 | 0.1917 | 0.2020 | 0.2123 | 0.2226 | 0.2329 | 0.2432
$S$ | 41.86 | 28.54 | 373.05 | 17.34 | 435.13 | 16.48 | 177.57 | 305.81
Table 4:  The action functional value of pathways terminated on $(3\pi/4,\pi)$ $\hat{\theta}$ | 2.4 | 2.5 | 2.6 | 2.7 | 2.8 | 2.9 | 3.0 | $\pi$
---|---|---|---|---|---|---|---|---
$w_{0}$ | 0.0454 | 0.0423 | 0.0393 | 0.0363 | 0.0333 | 0.0306 | 0.0279 | 0.0253
$v_{0}$ | 0.2408 | 0.2509 | 0.2611 | 0.2711 | 0.2811 | 0.291 | 0.301 | 0.3108
$S$ | 197.84 | 320.53 | 208.36 | 804.55 | 1055.6 | 161.68 | 224.56 | 1906.5
Table 5:  The action functional value of pathways terminated on $(-\pi/4,0)$ $\hat{\theta}$ | $-\pi/4$ | -0.6 | -0.5 | -0.4 | -0.3 | -0.2 | -0.1
---|---|---|---|---|---|---|---
$w_{0}$ | 0.0957 | 0.0974 | 0.0989 | 0.1002 | 0.1012 | 0.102 | 0.1024
$v_{0}$ | -0.0692 | -0.0593 | -0.0494 | -0.0396 | -0.0297 | -0.0197 | -0.0099
$S$ | 12.75 | 7.96 | 6.22 | 5.36 | 4.84 | 4.52 | 4.29
Table 6:  The action functional value of pathways terminated on $(-\pi/2,-\pi/4)$ $\hat{\theta}$ | $-\pi/2$ | -1.5 | -1.4 | -1.3 | -1.2 | -1.1 | -1.0 | -0.9
---|---|---|---|---|---|---|---|---
$w_{0}$ | 0.0747 | 0.0773 | 0.0799 | 0.0824 | 0.0849 | 0.0872 | 0.0896 | 0.0917
$v_{0}$ | -0.1595 | -0.1494 | -0.1393 | -0.1293 | -0.1192 | -0.1092 | -0.0992 | -0.0892
$S$ | 5193.3 | 235.1 | 204.79 | 996.73 | 1114.8 | 774.8 | 435.49 | 7503.1
Table 7:  The action functional value of pathways terminated on $(-3\pi/4,-\pi/2)$ $\hat{\theta}$ | $-3\pi/4$ | -2.3 | -2.2 | -2.1 | -2.0 | -1.9 | -1.8 | -1.7
---|---|---|---|---|---|---|---|---
$w_{0}$ | 0.0522 | 0.0547 | 0.0573 | 0.0599 | 0.0625 | 0.0652 | 0.0679 | 0.0706
$v_{0}$ | -0.2314 | -0.2312 | -0.221 | -0.2107 | -0.2005 | -0.1904 | -0.1802 | -0.17
$S$ | 104.15 | 7134.9 | 50.41 | 4051.5 | 68.02 | 387.07 | 148.92 | 131.47
Table 8:  The action functional value of pathways terminated on $(-\pi,-3\pi/4)$ $\hat{\theta}$ | $-\pi$ | -3.1 | -3.0 | -2.9 | -2.8 | -2.7 | -2.6 | -2.5
---|---|---|---|---|---|---|---|---
$w_{0}$ | 0.0342 | 0.0547 | 0.0384 | 0.0406 | 0.043 | 0.0455 | 0.048 | 0.0507
$v_{0}$ | -0.3272 | -0.317 | -0.3068 | -0.2966 | -0.2864 | -0.2761 | -0.2658 | -0.2556
$S$ | 76.65 | 103.4 | 522.58 | 58.17 | 294.69 | 452.05 | 66.9 | 633.63

According to the minimal action functional value on each segment, we can get
the most probable transition pathways from the origin to the corresponding
points on the circle. However, not all the pathways could reach exactly to the
circle, because the pathways are obtained from the given initial velocities,
which are estimated from the neural network algorithm. Even we have expanded
the sample size and limited the predict error, there still exist pathways that
did not satisfy the judgment condition whose final points were not exactly
located on the circle. For these pathways, we call them “unreachable”. Tables
1-8 give us indication to find the most probable transition paths from the
origin to the segments on surface and we show them in Figures 2-3.

Figure 2:  The most probable transition pathways (with minimal action
functional value) terminated on the upper half segments of the circle are
shown in the figure. The blue pathways denote the reachable transitions, while
the red ones fail to transit namely unreachable pathways.  Figure 3:  The most
probable transition pathways (with minimal action functional value) terminated
on the half bottom of the circle.

According to the Onsager-Machlup theory of stochastic dynamical systems, the
probability of the transition pathway on a small tube can be estimated as

$P_{1}({||X-z||_{T}<\delta}|z(0)=(0,0),z(T)=(r,\theta))\approx
c(T,\delta)e^{-S(z,\dot{z})}.$ (11)

When the parameters $\delta$ and $T$ are fixed, it can be seen that $P_{1}$ is
only dependent on the action functional value, and the coefficient
$c(T,\delta)$ can be regarded as a constant which does not vary with the final
location on the membrane. As we are aiming at qualitatively exploring the
difference of the probabilities transiting to each segment, here we get the
approximated probability of the most probable transition pathway on each
segment by ignoring the constant $c(T,\delta)$ in the following table:

Table 9:  The probability $P_{1}$ ($\propto e^{-S(z,\dot{z})}$) of most probable transition pathway on each segment $Segment$ | $(0,\pi/4)$ | $(\pi/4,\pi/2)$ | $(\pi/2,3\pi/4)$ | $(3\pi/4,\pi)$
---|---|---|---|---
$P_{1}$ | $0.02$ | $3.02\times 10^{-13}$ | $7.10\times 10^{-8}$ | $6.07\times 10^{-71}$
$Segment$ | $(-\pi/4,0)$ | $(-\pi/2,-\pi/4)$ | $(-3\pi/4,-\pi/2)$ | $(-\pi,-3\pi/4)$
$P_{1}$ | $0.01$ | $1.15\times 10^{-89}$ | $1.28\times 10^{-22}$ | $5.46\times 10^{-26}$

It should be noted that the most probable transition pathways shown above are
different from the real trajectories, as the calculated result is smooth
curves obtained from solving the Euler-Lagrange equation, but it
quantitatively points out the most possible reaching sites on the membrane,
and there exist at least one trajectory around the most probable transition
pathway.

###  2.2 The probability of the ligand arriving at the cell membrane

The intercellular motion of the ligand can be regarded as a jump-diffusion
process along the cell membrane [23, 24]. Ignoring the external forces of the
external environment, which can be characterized by the Lévy motion obeying
$L^{\alpha}_{t}\sim S_{\alpha}(|t|^{1/\alpha},0,0)$ . The model can be
represented as [30]:

$dX_{t}=\sigma dB_{t}+\epsilon dL_{t}.$ (12)

The corresponding Fokker-Planck equation is

$\frac{\partial p(x,t)}{\partial t}=\sigma^{2}\Delta
p(x,t)+\epsilon^{\alpha}\frac{\partial^{\alpha}}{\partial|x|^{\alpha}}p(x,t).$
(13)

The ligand locates on the site $\theta$ with probability

$P_{2}(\theta)\triangleq\int_{\theta-\delta}^{\theta+\delta}p(x,t)ds,$ (14)

here $P_{2}(\theta)$ can be calculated by integrating the density function
$p(x,t)$ over$(\theta-\delta,\theta+\delta)$, $s$ is the arc length between
$(\theta-\delta)$ and $(\theta+\delta)$, where $p(x,t)$ satisfies the Fokker-
Planck equation. For computational convenience, we treat small segments on the
cell membrane as one dimensional. In this case, there is no symmetry for the
probability of the ligand arriving at the cell membrane. The fractional order
derivation can be defined as follows [30]:

$\displaystyle\frac{\partial^{\alpha}}{\partial|x|^{\alpha}}p(x,t)$
$\displaystyle=-(-\Delta)^{\frac{\alpha}{2}}p(x,t)$ (15)
$\displaystyle=\int_{\mathbb{R}\backslash\\{0\\}}[p(x+y,t)-p(x,t)-I_{\\{|y|<1\\}}(y)yp^{\prime}(x)]\nu_{\alpha}(dy).$

This nonlocal equation can be numerically solved by a similar finite
difference method as in [44]. We divide the cell membrane into the same parts
as that in subsection 2.1 and list some computing results in Table 10, Table
11, Table 12. We also observe the evolution of $P_{2}$ with related parameters
as shown in the following figures.

Figure 4:  (a) Compare Gaussian and non-Gaussian noise. (b) Evolution with
different non-Gaussian noise intensity $\epsilon$: T=5, $\sigma=0$,
$\alpha=0.5$.

In Figure 4, we present the probability $P_{2}$ in the cases of Brownian
motion and Lévy motion with different noise intensities. we draw it by setting
$T=5$, $\sigma=1$ in $[0,\pi]$ domain. It shows that the probability $P_{2}$
with Gaussian noise is much bigger than that with non-Gaussian noise, which
means that the ligand with purely diffusion is more likely to arrive at a
given site on the surface than with jump-diffusion motion, and the larger the
noise intensity of the Lévy motion is, the less possibly the ligand can find
the site. That can be explained by the difference between the two motions: the
particles with Brownian motion find the binding site one by one and can
exactly locate a site on the memnrane, while in the jump-diffusion case, the
particles could jump over the site and wander around that, which will reduce
the possibility of correct bindings.

Figure 5: Evolution with different non-Gaussian index of Lévy process
$\alpha$: T=5, $\sigma=0.25$, $\epsilon=0.25$. (a) $[0,\pi]$. (b)$[\pi,2\pi]$.

To explore the influence of non-Gaussian index on probability $P_{2}$, we
paint Figure 5 by setting $T=5$, $\sigma=0.25$ and $\epsilon=0.25$ and
consider that in $[0,\pi]$ and $[\pi,2\pi]$ domains separately. In $[0,\pi]$,
$P_{2}$ first increases then decreases with respect to time $t$ and the non-
Gaussian index of Lévy process $\alpha$ has a significant effect on increasing
the probability $P_{2}$, which also tells us that as time goes, it becomes
difficult for the ligand to bind the given site, and the strong non-
Gaussianity leads to a same result. In $[0,\pi]$ and $[\pi,2\pi]$, $P_{2}$ has
the similar results in time.

Figure 6:  Evolution with different Gaussian noise intensity $\sigma$: T=5,
$\epsilon=0.25$. (a) $[0,\pi]$. (b)$[\pi,2\pi]$.

Figure 7:  Evolution with different non-Gaussian noise intensity $\epsilon$:
T=5, $\sigma=0.15$. (a) $[0,\pi]$. (b)$[\pi,2\pi]$.

In Figure 6, the ligand locating probability $P_{2}$ is computed when T=5 and
$\epsilon=0.25$. We can see that in both $[0,\pi]$ and $[\pi,2\pi]$ domains,
$P_{2}$ increases monotonically with respect to non-Gaussian index of Lévy
process $\alpha$. Besides, the noise intensity is helpful to increase the
binding probability $P_{2}$. We present the ligand locating probability
$P_{2}$ with T=5 and $\sigma=0.15$ in Figure 7. For small Lévy noise
intensity, $P_{2}$ increases monotonically with respect to non-Gaussian index
of Lévy process $\alpha$; while for large Lévy noise intensity, $P_{2}$
increases first and then decreases with $\alpha$. That points to the existence
of an optimal $\alpha$ which will lead to the max probability for the ligand
to reach the binding site.

Table 10:  The probability $P_{2}$ on each segment with $\epsilon=0.25$, $\sigma=0.1$. $Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
---|---|---|---|---
$\alpha=0.5$ | $6.76\times 10^{-187}$ | $9.46\times 10^{-108}$ | $4.90\times 10^{-50}$ | $9.36\times 10^{-14}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=0.5$ | $6.52\times 10^{-298}$ | $2.42\times 10^{-195}$ | $3.32\times 10^{-114}$ | $1.69\times 10^{-54}$
$Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
$\alpha=1.5$ | $1.26\times 10^{-187}$ | $1.76\times 10^{-108}$ | $9.13\times 10^{-51}$ | $1.74\times 10^{-14}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=1.5$ | $1.21\times 10^{-298}$ | $4.51\times 10^{-196}$ | $6.20\times 10^{-115}$ | $3.15\times 10^{-55}$
Table 11:  The probability $P_{2}$ on each segment with $\epsilon=0.25$, $\sigma=0.3$. $Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
---|---|---|---|---
$\alpha=0.5$ | $9.03\times 10^{-187}$ | $1.26\times 10^{-107}$ | $6.54\times 10^{-50}$ | $1.25\times 10^{-13}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=0.5$ | $8.71\times 10^{-298}$ | $3.23\times 10^{-195}$ | $4.44\times 10^{-114}$ | $2.26\times 10^{-54}$
$Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
$\alpha=1.5$ | $1.39\times 10^{-187}$ | $1.95\times 10^{-108}$ | $1.01\times 10^{-50}$ | $1.93\times 10^{-14}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=1.5$ | $1.34\times 10^{-298}$ | $4.99\times 10^{-196}$ | $6.86\times 10^{-115}$ | $3.49\times 10^{-55}$
Table 12:  The probability $P_{2}$ on each segment with $\epsilon=0.4$, $\sigma=0.3$. $Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
---|---|---|---|---
$\alpha=0.5$ | $6.83\times 10^{-187}$ | $9.56\times 10^{-108}$ | $4.95\times 10^{-50}$ | $9.45\times 10^{-14}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=0.5$ | $6.58\times 10^{-298}$ | $2.45\times 10^{-195}$ | $3.36\times 10^{-114}$ | $1.71\times 10^{-54}$
$Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
$\alpha=1.5$ | $6.43\times 10^{-188}$ | $9.00\times 10^{-109}$ | $4.66\times 10^{-51}$ | $8.91\times 10^{-15}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=1.5$ | $6.21\times 10^{-299}$ | $2.30\times 10^{-196}$ | $3.17\times 10^{-115}$ | $1.61\times 10^{-55}$

##  3 The Most Probable Binding Probability

The last section gives the quantitative illustration about the most probable
probability of receptors and ligands location on the cell surface
respectively. Based on that we could further get the most probable binding
site on the surface by combining $P_{1}$ and $P_{2}$ and calculating the most
probable binding probability $P=P_{1}P_{2}$. Indeed, it is hard to figure out
the probabilities of all the possible parameters, so we just base on the
probability of the partial location what we have figured out and provide the
corresponding most probable binding probability.

The binding probability for receptors and ligands at different site with
different parameters was shown in Tables 13-15. We find that the binding
probability is related with diffusion coefficient $\sigma$, non-Gaussnian
noise coefficient $\epsilon$ and non-Gaussnian index $\alpha$. Note that
$P_{1}$ has a common estimation factor, but we ignore it in the following
calculations.

Table 13:  The binding probability $P$ on each segment with $\epsilon=0.25$, $\sigma=0.1$. $Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
---|---|---|---|---
$\alpha=0.5$ | $1.35\times 10^{-188}$ | $2.85\times 10^{-120}$ | $3.47\times 10^{-57}$ | $5.68\times 10^{-84}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=0.5$ | $3.55\times 10^{-323}$ | $3.09\times 10^{-217}$ | $3.81\times 10^{-203}$ | $1.69\times 10^{-56}$
$Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
$\alpha=1.5$ | $2.52\times 10^{-189}$ | $5.31\times 10^{-121}$ | $6.48\times 10^{-58}$ | $1.05\times 10^{-84}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=1.5$ | $6.60\times 10^{-324}$ | $5.77\times 10^{-218}$ | $7.13\times 10^{-204}$ | $3.15\times 10^{-57}$
Table 14:  The binding probability $P$ on each segment with $\epsilon=0.25$, $\sigma=0.3$. $Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
---|---|---|---|---
$\alpha=0.5$ | $1.80\times 10^{-188}$ | $3.80\times 10^{-120}$ | $4.64\times 10^{-57}$ | $7.58\times 10^{-84}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=0.5$ | $4.75\times 10^{-323}$ | $4.13\times 10^{-217}$ | $5.10\times 10^{-203}$ | $2.26\times 10^{-56}$
$Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
$\alpha=1.5$ | $2.78\times 10^{-189}$ | $5.88\times 10^{-121}$ | $7.17\times 10^{-58}$ | $1.17\times 10^{-84}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=1.5$ | $7.31\times 10^{-324}$ | $6.38\times 10^{-218}$ | $7.88\times 10^{-204}$ | $3.49\times 10^{-57}$
Table 15:  The binding probability $P$ on each segment with $\epsilon=0.4$, $\sigma=0.3$. $Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
---|---|---|---|---
$\alpha=0.5$ | $1.36\times 10^{-188}$ | $2.88\times 10^{-120}$ | $3.51\times 10^{-57}$ | $5.73\times 10^{-84}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=0.5$ | $3.59\times 10^{-323}$ | $3.06\times 10^{-217}$ | $3.86\times 10^{-203}$ | $1.71\times 10^{-56}$
$Segment$ | $(0,\frac{\pi}{4})$ | $(\frac{\pi}{4},\frac{\pi}{2})$ | $(\frac{\pi}{2},\frac{3\pi}{4})$ | $(\frac{3\pi}{4},\pi)$
$\alpha=1.5$ | $1.28\times 10^{-189}$ | $2.17\times 10^{-121}$ | $3.30\times 10^{-58}$ | $5.40\times 10^{-85}$
$Segment$ | $(\pi,\frac{5\pi}{4})$ | $(\frac{5\pi}{4},\frac{3\pi}{2})$ | $(\frac{3\pi}{2},\frac{7\pi}{4})$ | $(\frac{7\pi}{4},2\pi)$
$\alpha=1.5$ | $3.39\times 10^{-324}$ | $2.94\times 10^{-218}$ | $3.64\times 10^{-204}$ | $1.61\times 10^{-57}$

We have picked out the most probable binding location corresponding to the
largest binding probability in each case in bold font. According to the
results, we find that the most probable binding location with largest
probability is in $(\frac{7\pi}{4},2\pi)$, both in $\alpha=0.5$ and
$\alpha=1.5$ cases, with our choosen noise intensities $\sigma$ and
$\epsilon$. In this way, we can make predictions in the cases of other non-
Gaussian index $\alpha$ and noise intensities $\sigma$ and $\epsilon$. This
receptor-ligand binding prediction method is based on a simple diffusion model
for the receptor and a jump-diffusion model for the ligand.

##  4 Conclusion

In this paper, we have considered stochastic dynamics for an cell surface
receptor binding to an extracellular ligand on the cell membrane. The receptor
is modeled as a usual diffusion, while the ligand is regarded as a pure jump-
diffusion. To explore the binding probability, we divided the cell surface
into eight segments, and calculate the largest probabilities of receptor and
ligand reaching each segment. As for the receptor, we used the Onsager-Machlup
functional to compute the most probable transition pathway from the cell
center to the cell membrane. To solve the Euler-Lagrange equation, we proposed
a neural shooting method for the most probable transition pathway. We further
on obtained the action functional value as well as the transition probability.
The motion of extracellular ligand is quantized by a fractional Fokker-Planck
equation. We used a fast discretization method to solve this equation and got
the probability of the ligand arriving at each segment on the membrane. We
also investigated the influence of the system parameters on the ligand
arriving probability on cell membrane, and found the Brownian noise had
advantage on locating the closer site on the membrane over Lévy noise, while
the phenomenon flips when the ligand arrived at a distant site. This result is
consistent with that discussed in [45, 46].

Based on the probabilities of the receptor diffusing to the cell membrane and
the ligand arriving at the cell membrane, we further compute the largest
binding probability and indicate the most probable binding site. In our study,
the binding event happens once the receptor and ligand arrive at the same
location on the cell membrane, which is a rare event. Our investigation
provides a quantitative analysis for the extreme statistics. Besides, our work
helps understand the receptor-ligand binding on the membrane, which also sheds
light on understanding of cell’s response to external stimuli and
communication with other cells.

To sum up, we have devised a method for predicting the likelihood for
receptor-ligand binding on cell surface. Specifically,

(a)We use neural networks based on Onsager-Machlup function to compute the
most probable transition pathway of the receptor diffusion to the cell
membrane.

(b)We compute the probability of the receptor diffusion to a site on the cell
membrane and the probability of ligand locating the same site respectively. We
also compare the impact of non-Gaussnian noise with Gaussnian noise in the
ligand dynamics.

(c)We combine the probability of the receptor diffusion to a site on the cell
membrane and the probability of ligand locating the same site to obtain the
most probable binding site.

## Acknowledgments

We would like to thank Jianyu Chen for helpful discussions. This work is
supported by the National Natural Science Foundation of China (NSFC) (Grant
No.11901536 and No.12101473 ).

## Data Availability

The data that support the findings of this study are openly available in
Github https://github.com/StringsLi/binding-code.

## Appendix

###  Lévy motion

A scalar stable Lévy motion $L_{t}^{\alpha}$, for $0<\alpha<2$, is a non-
Gaussian stochastic process with the following properties [30, 47, 48, 49]:
(i) $L_{0}^{\alpha}=0$, almost surely;
(ii) $L_{t}^{\alpha}$ has independent increments;
(iii)$L_{t}^{\alpha}$ has stationary increments:
$L_{t}^{\alpha}-L_{s}^{\alpha}$ has probability distribution
$S_{\alpha}((t-s)^{\frac{1}{\alpha}},0,0)$ for $s\leq t$; in particular,
$L_{t}^{\alpha}$ has distribution $S_{\alpha}(t^{\frac{1}{\alpha}},0,0)$;
(iv) $L_{t}^{\alpha}$ has stochastically continuous sample paths, i.e.,
$L_{t}^{\alpha}\rightarrow L_{s}^{\alpha}$ in probability, as $t\rightarrow
s$.

Here $S_{\alpha}(\sigma,\beta,\mu)$ is the so-called stable distribution [30,
49] and is determined by four parameters, non-Gaussianity index
$\alpha(0<\alpha<2)$, skewness index $\beta(-1\leq\beta\leq 1)$, shift index
$\mu(-\infty<\mu<+\infty)$ and scale index $\sigma(\sigma\geq 0)$.

The stable Lévy motion $L_{t}^{\alpha}$ has the jump measure

$\nu_{\alpha}(dy)=C_{\alpha}|y|^{-(1+\alpha)}\,dy,$

where the coefficient

$C_{\alpha}=\frac{\alpha}{2^{1-\alpha}\sqrt{\pi}}\frac{\Gamma(\frac{1+\alpha}{2})}{\Gamma(1-\frac{\alpha}{2})}.$

Note that the well-known Brownian motion $B_{t}$ is a special case
corresponding to $\alpha=2$.

### Nonlocal Fokker-Planck equation and numerical methods

Consider a scalar stochastic differential equation with Lévy noise

$dX_{t}=f(X_{t})dt+\varepsilon dL_{t},\;\;X_{0}=x_{0},$ (16)

where $f$ is a given vector field (or drift) and $\varepsilon$ is the noise
intensity.

The generator for this stochastic differential equation is

$A\varphi(x)=f(x)\varphi^{\prime}(x)+\int_{\mathbb{R}^{1}\backslash\\{0\\}}[\varphi(x+\varepsilon
y)-\varphi(x)]\nu_{\alpha}(dy).$ (17)

Let $z=\varepsilon y$. The generator becomes

$A\varphi(x)=f(x)\varphi^{\prime}(x)+\varepsilon^{\alpha}\int_{\mathbb{R}^{1}\backslash\\{0\\}}[\varphi(x+z)-\varphi(x)]\nu_{\alpha}(dz).$

The Fokker-Planck equation for this stochastic differential equation, i.e.,
the probability density $p(x,t)$ for the solution process $X_{t}$ with initial
condition $X_{0}=x_{0}$ is [30]

$p_{t}=A^{*}p,\;\;p(x,0)=\delta(x-x_{0}),$ (18)

where $A^{*}$ is the adjoint operator of the generator $A$ in Hilbert space
$L^{2}(R^{1})$, as defined by

$\int_{\mathbb{R}^{1}\backslash\\{0\\}}A\varphi(x)u(x)dx=\int_{\mathbb{R}^{1}\backslash\\{0\\}}\varphi(x)A^{*}u(x)dx.$

Then via integration by parts, we get the adjoint operator for $A$

$A^{*}u(x)=\int_{\mathbb{R}^{1}\backslash\\{0\\}}\varepsilon^{\alpha}[u(x-y)-u(x)]\;\nu_{\alpha}(dy).$
(19)

Thus we have the nonlocal Fokker-Planck equation

$p_{t}=-(f(x)p(x,t))_{x}+\int_{\mathbb{R}^{1}\backslash\\{0\\}}\varepsilon^{\alpha}[p(x-y,t)-p(x,t)]\;\nu_{\alpha}(dy).$
(20)

We use a numerical finite difference method developed in Gao et al. [44] to
simulate the nonlocal Fokker-Planck equation (20).

## References

## References

  * [1] Bressloff P C 2021 Stochastic Processes in Cell Biology volumes I and II (Springer, New York)
  * [2] Schuss Z, Holcman D 2017 100 years after smoluchowski: stochastic processes in cell biology _J. Phys. A: Math. Theor._ 50 093002
  * [3] Bialek W 2012 Biophysics:Searching for Principles. (Princeton University Press, Princeton)
  * [4] Lauffenburger D A, Linderman J J 1996 Receptors: Models for Binding, Trafficking, and Signaling ( Oxford University Press, Oxford )
  * [5] Qian H, Ge H 2021 Stochastic Chemical Reaction Systems in Biology( Springer, New York)
  * [6] Anderson D F, Kurtz T G 2015 Stochastic Analysis of Biochemical Systems (Springer, New York)
  * [7] Logan J D 2013 Applied Mathematics ( John Wiley and Sons, New York)
  * [8] Wiley H S 1992 Receptors: topology, dynamics, and regulation. In Fundamentals of Medical Cell Biology 5A
  * [9] Lawley S, Miles C 2019 How Receptor Surface Diffusion and Cell Rotation Increase Association Rates. SIAM J. Appl. Math. 79 1124-1146
  * [10] Handy G, Lawley 2021 Revising Berg-Purcell for finite receptor kinetics. Biophys. J.120(11) 2237-2248
  * [11] Lawley S, Lindsay A, Miles C 2020 Receptor organization determines the limits of single-cell source location detection. Phys. Rev. Lett. 125(1), 018102
  * [12] Weikl T R, Asfaw M, Krobath H, Rozycki B, Lipowsky R 2009 Adhesion of membranes via receptor-ligand complexes: domain formation, binding cooperativity, and active processes. Soft Matter. 5 3213-3224
  * [13] Abney J R, Scalettar B A, Owicki J C 1989 Self diffusion of interacting membrane proteins. Biophys. J. 55 817-833
  * [14] Alberts B, Johnson A, Lewis J, Raff M, Walter K R 2008 Molecular Biology of the Cell (Garland, New York )
  * [15] Zheng X, Cheung S L, Schroeder J A , Jiang L, Zohar Y 2011 Cell receptor and surface ligand density effects on dynamic states of adhering circulating tumor cells. Lab on A Chip11(20) 3431-3439
  * [16] Hake J, Lines G T 2008 Stochastic binding of Ca2+ ions in the dyadic cleft; continuous versus random walk description of diffusion. Biophys. J. 94 4184-4201
  * [17] Lin J, Low-Nam S T, Alfieri K N, McAffee D B, Fay N C, Groves J T 2019 Mapping the stochastic sequence of individual ligand-receptor binding events to cellular activation: T cells act on the rare events. Science Signaling 12(564), eaat8715
  * [18] Thewes N, Thewes A, Loskill P, Peisker H, Bischoff M, Herrmann M, Santen L , Jacobs K 2015 Stochastic binding of staphylococcus aureus to hydrophobic surfaces. Soft Matter 11(46), 8913-8919
  * [19] Pastan I and Willingham M C 1985 _Endocytosis_.( New York: Plenum Press, 326 )
  * [20] Dobramysl U, Holcman D 2020 Triangulation sensing to determine the gradient source from diffusing particles to small cell receptors. Phys. Rev. Lett. 125 14, 148102
  * [21] Linn S, Lawley S 2022 Extreme hitting probabilities for diffusion. J. Phys. A: Math. Theor. 55 345002
  * [22] Basnayake K, Schuss Z, Holcman D 2019 Asymptotic formulas for extreme statistics of escape times in 1, 2 and 3-dimensions. J. Nonlinear Sci. 29 461-499
  * [23] Vauquelin G and Packeu A 2009 Ligands, their receptors and plasma membranes.Mol. Cell Endocrinol. 311(1-2) 1-10
  * [24] Reynolds A M, Geritz S A H 2015 Tsallis distributions, Lévy walks and correlated-type anomalous diffusion result from state-dependent diffusion. Physica A: Statistical Mechanics and its Applications  424 317-321
  * [25] Lomholt M A, Koren T, Metzler R, Klafter J 2008 Lévy strategies in intermittent search processes are advantageous. Proc. Natl. Acad. Sci. USA 105 11055
  * [26] Cheng X, Wang H, Wang X, Duan J, Li X 2019 Most probable transition pathways and maximal likely trajectories in a genetic regulatory system. Physica A 531 121779
  * [27] Hu J, Chen J 2021 Transition pathways for a class of high dimensional stochastic dynamical systems with Lévy noise. Chaos 31 063138
  * [28] Zheng Y, Yang F, Duan J, Sun X, Fu L, Kurths J 2020 The maximum likeli hood climate change for global warming under the influence of greenhouse effect under Lévy noise. Chaos 30(1) 013132
  * [29] Dykman M, Mori E, Ross J, Hunt P 1994 Large fluctuations and optimal paths in chemical kinetics. J. Chem. Phys. 100 5735-5750
  * [30] Duan J 2015 An Introduction to Stochastic Dynamics (Cambridge: Cambridge University Press)
  * [31] Li Y, Duan J, Liu X, Zhang Y 2020 Most probable dynamics of stochastic dynamical systems with exponentially light jump fluctuations. Chaos 30 063142
  * [32] Lu Y and Duan J 2020 Discovering transition phenomena from data of stochastic dynamical systems with Lévy noise. Chaos 30 093110
  * [33] Cai R, Chen X, Duan J, Kurths J, Li X 2017 Lévy noise-induced escape in an excitable system. J. Stat. Mech. 063503
  * [34] Hu J, Xu G K, Lipowsky R, Weikl T R 2015 Binding kinetics of membrane-anchored receptors and ligands: Molecular dynamics simulations and theory. J. Chem. Phys.  143(24) 243137
  * [35] Bouchaud J P, Georges A 1990 Anomalous diffusion in disordered media-statistical mechanisms, models and physical applications. Phys. Rep. 195 127-293
  * [36] Saffman P G, Delbruck M 1975 Brownian motion in biological membranes. Proc. Natl. Acad. Sci. USA 72 3111-3113
  * [37] Havlin S, ben-Avraham D 2002 Diffusion in disordered media. Adv. Phys. 51 187-292
  * [38] Holcman D , Schuss Z 2015 _Stochastic Narrow Escape in Molecular and Cellular Biology_( New York: Springer)
  * [39] Wang H, Duan J, Geng X , Chao Y 2020 Mean escape time for randomly switching narrow gates in a steady flow. Comput. Math. Appl. 79(10) 2795-2804
  * [40] Fujita T, Kotani S-i 1982 The Onsager-Machlup function for diffusion processes. Journal of Mathematics of Kyoto University 22(1) 115-130
  * [41] Capitaine M 1995 Onsager-Machlup functional for some smooth norms on Wiener space. Probability theory and related fields 102(2) 189-201
  * [42] Chao Y, Duan J 2019 The Onsager-Machlup function as Lagrangian for the most probable path of a jump-diffusion process. Nonlinearity  32 3715
  * [43] Ibraheem K I, Khalaf B M 2011 Shooting neural networks algorithm for solving boundary value problems in odes. Appl. Applied Math. 6(11) 1927-1941
  * [44] Gao T, Duan J, Li X 2016 Fokker-Planck equations for stochastic dynamical systems with symmetric Lévy motions. _Appl. Math. Comput._ 278, 1-20
  * [45] Palyulin V V, Chechkin A V, Metzler R 2014 Lévy flights do not always optimize random blind search for sparse targets. _Proc. Natl. Acad. Sci. USA_ 111 2931-2936
  * [46] Chen X, Cheng X, Kang Y and Duan J 2019 Target search of a protein on DNA in the presence of position-dependent bias. J. Stat. Mech. 033501
  * [47] Applebaum D 2009 Lévy Processes and Stochastic Calculus 2nd Edition. (Cambridge University Press, New York,USA)
  * [48] Sato K 1999 Lévy Processes and Infinitely Divisible Distributions (Cambridge University Press, New York,USA)
  * [49] Samorodnitsky G, Taqqu M S 1994 Stable Non-Gaussian Random Processes (Chapman and Hall, New York, USA)