Monte-Carlo Simulation of Particle Diffusion in Various Geometries and Application to Chemistry and Biology

1. Introduction

The simulation of systems comprising different types of molecules is of great interest in several fields, notably in chemistry and biological sciences. The conventional approach to simulate biological networks is to write down macroscopic rate equations and solve the corresponding differential equations numerically to obtain the time evolution of molecules concentrations [1-3]. In this method, the system evolution is deterministic, and it is implicitly assumed that the concentrations are large. However, in most biological systems, the concentrations of molecules usually range from nanomolar to micromolar; therefore, the interactions of these molecules are highly stochastic. Recently, various techniques have been developed to take into account the spatial distribution of individual molecules and the stochastic character of reactions between them. Several of these techniques are based on the Green's functions of the Diffusion Equation (DE), and they offer the advantage of being able to follow all the particles in time and space. In this book chapter, this method is applied to simple systems in one, two, and three dimensions. Two applications are discussed. The first is the simulation of ligands molecules near a plane membrane comprising receptors with the possibility of dissociation and initiation of signal transduction. The second application of this theory is the simulation of partially diffusion-controlled chemical reactions.

2. Sampling of the Green's functions of the diffusion equation

The Green's functions of the DE in simple systems under different boundary conditions are known from the theory of heat propagation in solids [4], since the DE and the equation for the propagation of heat are the same. With few assumptions, these Green's functions can be used to model particle systems in specific environments.

2.1. Diffusion equation and boundary conditions

The DE can be written

$$\frac{\partial }{\partial t}p(\text{r}\text{,t}|{\text{r}}_{\text{0}})=D{\nabla }^{\text{2}}p(\text{r}\text{,t}|{\text{r}}_{\text{0}}),$$E1

where r₀ is the initial position vector

Through the book chapter, variables written in bold are vectors.

of a particle and r is its position vector at time t. The diffusion coefficient of a particle (D) is used instead of the thermal diffusivity coefficient in the equation of the propagation of heat. In the model presented here, a particle is located at the position r₀ at t=0. This is expressed mathematically as:

$$p(\text{r}\text{,t}=\text{0}|{\text{r}}_{\text{0}})=\delta (\text{r}-{\text{r}}_{\text{0}}),$$E2

where δ(r-r₀)=δ(x-x₀)δ(y-y₀)δ(z-z₀) is the product of three Dirac's delta functions. To complete the description of the model the boundary conditions at the surfaces are specified. One important case is the radiation boundary condition at the surface X=0:

$$D\frac{\partial p(\text{x,y,z,t}|x_{\text{0}},y_{\text{0}},z_{\text{0}})}{\partial x}{|}_{x=\text{0}}=k_{a}p(\text{0,y,z,t}|x_{\text{0}}{\text{,y}}_{\text{0}},z_{\text{0}}).$$E3

The limit k_a→∞ is known as the absorbing boundary condition. In this case, the boundary condition is written p(0,y,z,t|x₀, y₀, z₀, t₀)=0. The limit k_a→0 is the reflecting boundary condition. The condition k_a>0 implies that particles may bind at the surface; therefore, the Green's function of these systems are generally not normalized to 1. In all systems described in Sections 3 and 4 the boundary conditions on surfaces, if any, are reflecting (k_a=0). In Section 5 two applications of the theory for which the boundary conditions are different will be discussed.

2.2. Sampling of the Green's functions

To simulate the state of a particle system after one time step, the positions of each particle are calculated by sampling the Green's function for the position r. The Green's function is written p(r, t|r₀) to emphasize that it represents the probability of a particle initially located at the position r₀ to be found at the position r at time t. In the systems presented in Sections 3 and 4, there is no binding of particles to boundaries. Therefore, p(r,t|r₀) is normalized as follows:

$${\displaystyle \underset{\Omega }{\int }p({\text{r}}_{\text{1}}{\text{,t}}_{\text{1}}|{\text{r}}_{\text{0}})d{\text{r}}_{\text{1}}}=\text{1},$$E4

where Ω is the domain where the particle can be found. Another important thing to consider is that the simulation of a particle trajectory can be done by using one or several time steps. This can be expressed mathematically as:

$$p\left({\text{r}}_{\text{2}},\Delta t_1+\Delta t_2|{\text{r}}_{\text{0}}\right)={\displaystyle \underset{\Omega }{\int }p\left({\text{r}}_{\text{2}},\Delta t_2|{\text{r}}_{\text{1}}\right)p\left({\text{r}}_{\text{1}},\Delta t_1|{\text{r}}_{\text{0}}\right)d{\text{r}}_{\text{1}}}.$$E5

That is, the probability of a particle to be found at r₂ after the time steps Δt₁ and Δt₂ is the same as going to the intermediate position r₁ after Δt₁ and then the final position r₂ at Δt₁+Δt₂. The integral over the domain Ω covers all intermediate positions. This is the Chapman-Kolmogorov equation, which hold in Markovian systems.

As seen in a previous paper [5], Equations (4) and (5) have to modified in systems where there is a possibility of binding of particles at the boundaries.

3. One-dimensional systems

One dimensional (1D) systems are simpler to study than those which are multi-dimensional. They are also of interest since several multi-dimensional systems can be approached as independent 1D systems [6]. In 1D, the DE is:

$$\frac{\partial p(\text{x,t}|x_{\text{0}})}{\partial t}=D\frac{{\partial }^{\text{2}}}{\partial x^{\text{2}}}p(\text{x,t}|x_{\text{0}}).$$E6

The initial condition is p(x,0|x₀)=δ(x-x₀). Three systems are described in this section: free diffusion, diffusion in the half-space X>0 and diffusion in the space between X=0 and X=L.

3.1. Free diffusion

The Green's function of the DE in 1D for a free particle, noted p_free(x,t|x₀), is well known:

$$p_{\text{free}}\left(\text{x,t}|x_{\text{0}}\right)=\frac{\text{1}}{\sqrt{\text{4}\pi \text{Dt}}}\text{exp}\left[-\frac{{({\text{x-x}}_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right].$$E7

This is a Gaussian distribution, with variance σ²=2Dt. This distribution can be sampled by using, for example, the Box-Muller method [7].

Algorithm 1: Sampling algorithm for p_free(x,Δt|x₀) (free 1D diffusion)

GENERATE N as a standard Normal

$$X\leftarrow x_{\text{0}}+\sqrt{D\Delta t}$$

RETURN X.

The sampling algorithm is used to generate the positions of 1,000,000 particles after one time step. The results are stored in histograms and normalized. The Green's functions and the results of sampling are illustrated

Through this book chapter, the diffusion coefficient is set to D=1 in Figures.

in Figure 1.

Figure 1.

Green's functions p_free(x,t|x₀) at t=1, 2, 4, 8, 16 and 32 time units, for a particle initially at x₀=2. The lines are analytical predictions; the points are the result of sampling from 1,000,000 particle histories.

Because the sampling algorithms use random numbers, the simulation results vary from one simulation to another. Usually, the larger the number of histories, the smaller the statistical fluctuation of the results. The algorithms presented in this book chapter are exact, i.e. the simulation results are expected to converge to the analytical Green's function at the limit of an infinite number of histories. In practice, 1,000,000 particle histories are usually sufficient to converge to the analytical Green's function and to confirm the validity of the algorithm. The time discretization property of the Green's function can also be used to validate the algorithm: since a simulation can be done in several time steps, the positions of the particles obtained after the first time step are used as the initial positions for the next step.

3.2. Diffusion in the half-space x>0

As it will be seen in the applications, this system may be used to represent the motion of particles near a membrane. The reflective boundary condition at x=0 is given by

$$D\frac{\partial p(\text{x,t}|x_{\text{0}})}{\partial x}{|}_{x=\text{0}}=\text{0}.$$E8

The solution of the DE in 1D for this system, noted p_ref(x,t|x₀), is given in ref. [8]:

$$p_{\text{ref}}\left(x,t|x_0\right)=\frac{1}{\sqrt{4\pi Dt}}\left\{\exp \left[-\frac{{(x-x_0)}^2}{4Dt}\right]+\exp \left[-\frac{{(x+x_0)}^2}{4Dt}\right]\right\}.$$E9

This Green's function is the sum of two Gaussian distributions. It can be sampled by the method described in Algorithm 2 (given in ref. [8]):

Algorithm 2: Sampling algorithm for p_ref(x,Δt|x₀) (free 1D diffusion with reflection at x=0)

SET $$N_{\text{0}}\leftarrow (\text{1}/\text{2})\text{Erfc}[-x_{\text{0}}/\sqrt{\text{4DΔt}}]$$
GENERATE uniform [0,1] random variatesU,V
IF (U<N0)	$$X\leftarrow x_{\text{0}}+\sqrt{\text{4DΔt}}{\text{Erfc}}^{\text{-1}}[\text{VErfc}(-x_{\text{0}}/\sqrt{\text{4DΔt}})]$$
ELSE	$$X\leftarrow -x_{\text{0}}+\sqrt{\text{4DΔt}}{\text{Erfc}}^{\text{-1}}[\text{VErfc}(x_{\text{0}}/\sqrt{\text{4DΔt}})]$$
RETURN X

where Erfc(x) is the complementary error function [9]. The Green's functions and the result of sampling are illustrated in Figure 2.

Figure 2.

Green's functions p_ref(x,t|x₀) at t=1, 2, 4, 8, 16 and 32 time units, for a particle in the half-space x>0 initially at x₀=2.5. The lines are analytical predictions; the points are the result of sampling from 1,000,000 particles histories.

3.3. Diffusion of a particle between x=0 and x=L

This system is used to simulate the motion of particles in a confined space. In this case, the reflective boundary conditions at x=0 and x=L are given by

$$\begin{array}{l}D\frac{\partial p(\text{x,t}|x_{\text{0}})}{\partial x}{|}_{x=\text{0}}=\text{0},\\ D\frac{\partial p(\text{x,t}|x_{\text{0}})}{\partial x}{|}_{x=L}=\text{0}.\end{array}$$E10

The Green's function of the DE in 1D for a particle between the boundaries, noted p_conf(x,t|x₀), is [4]:

$$p_{\text{conf}}\left(\text{x,t}|x_{\text{0}}\right)=\frac{\text{1}}{L}{\displaystyle \sum _{n=-\infty }^{\infty }\text{exp}\left(-\frac{{\pi }^{\text{2}}{n}^{\text{2}}\text{Dt}}{L^{\text{2}}}\right)\text{cos}\frac{n\pi x}{L}\text{cos}\frac{n\pi x_{\text{0}}}{L}}.$$E11

This function is complicated by an infinite sum. However, it is possible to use the following algorithm to generate trajectories of particles in this environment [10]. The algorithm is based on a series method, as described in [11]. For this section the functions f(x) and f_n(x) are defined:

$$\begin{array}{l}f(x)=\frac{\text{1}}{L}+\text{2}{\displaystyle \sum _{n=\text{1}}^{\infty }{f}_n(x)},\\ f_n(x)\stackrel{\text{def}}{\text{=}}\frac{\text{1}}{L}\text{exp}\left(-\frac{{\pi }^{\text{2}}{n}^{\text{2}}\text{Dt}}{L^{\text{2}}}\right)\text{cos}\frac{n\pi x}{L}\text{cos}\frac{n\pi x_{\text{0}}}{L}.\end{array}$$E12

The algorithm is:

Algorithm 3A: Generation of random variateX distributed as p_conf(x,Δt|x₀) (particle between x=0 and x=L in 1D)

DEFINE $$H\leftarrow \text{1}/L+\text{1}/\sqrt{\pi D\Delta t}$$
REPEAT
{
	GENERATE U,V uniform on [0,1] and X uniform on [0,L]
	SET Y ← VH
}
UNTIL Y≤f(X)
RETURN X

The verification of Y≤f(X) is dose by using the second part of the algorithm:

Algorithm 3B: Verification of Y<f(X)

SET S ← 1/L+2f1(X), k ← 1 (S holds the approximation sum)
REPEAT FOREVER
{
	T = (L/π^2kDΔt)exp(-π^2k^2DΔt/L^2)
	IF (Y ≥ S+T) THEN RETURN "Y ≥ f(X)" and EXIT
	IF (Y ≤ S-T) THEN RETURN "Y ≤ f(X)" and EXIT
	k ← k+1
	S ← S+2fk(X)
}

The Green's function and the result of sampling for 1,000,000 histories of particles initially located at x₀=5 are shown in Figure 3.

Figure 3.

Green's functions p_conf(x,t|x₀) at t=1, 2, 4, 8, 16 and 32 time units, for a particle initially at x₀=5, allowed to occupy the space between x=0 and x=10. The lines are analytical predictions; the points are the result of sampling of 1,000,000 particle histories.

When time evolves, as illustrated in Figure 3, the particles become uniformly distributed between x=0 and x=10.

4. Two-dimensional systems

Two dimensional (2D) systems are interesting, because the particles trajectories can be visualized in a plane. The DE in 2D can be written:

$$\frac{\partial p(\text{x,y,t}|x_{\text{0}},y_{\text{0}})}{\partial t}=D\left[\frac{{\partial }^{\text{2}}}{\partial x^{\text{2}}}+\frac{{\partial }^{\text{2}}}{\partial y^{\text{2}}}\right]p(\text{x,y,t}|x_{\text{0}},y_{\text{0}}).$$E13

The initial condition is p(x,y,t=0|x₀, y₀)=δ(x-x₀)δ(y-y₀). In several cases, the Green's function can be written as the product of the functions p_x(x,t|x₀) and p_y(y,t|y₀) (see ref. [6]):

$$p(\text{x,y,t}|x_{\text{0}}{\text{,y}}_{\text{0}})=p_x(\text{x,t}|x_{\text{0}})p_y(\text{y,t}|y_{\text{0}}).$$E14

This property of the Green’s function allows the simulation of the motion of particles in the direction x and y independently.

4.1. Free diffusion

The solution of the DE is the product of two Gaussian functions:

$$p_{\text{free}}\left(\text{x,y,t}|x_{\text{0}}{\text{,y}}_{\text{0}}\right)=\frac{\text{1}}{\text{4πDt}}\text{exp}\left[-\frac{{\left({\text{x-x}}_{\text{0}}\right)}^{\text{2}}+{\left({\text{y-y}}_{\text{0}}\right)}^{\text{2}}}{\text{4Dt}}\right]\equiv p_{\text{free}}\left(\text{x,t}|x_{\text{0}}\right)p_{\text{free}}\left(\text{y,t}|y_{\text{0}}\right).$$E15

Therefore, sampling of this Green's function can be done by sampling two Gaussian functions:

Algorithm 4: Sampling algorithm for p_free(x,y,Δt|x₀,y₀) (free 2D diffusion)

GENERATE N1,N2 as a standard Normals

$$X\leftarrow x_{\text{0}}+\sqrt{D\Delta t}{N}_{\text{1}}$$

$$Y\leftarrow y_{\text{0}}+\sqrt{D\Delta t}{N}_{\text{2}}$$

RETURN (X,Y).

Plotting the Green's functions in X and Y would give figures identical to Figure 1. However, in 2D, it is useful to use polar coordinates. That is, the following transformation is used:

$$\begin{array}{l}x=r\text{cos}(\theta ),\\ y=r\text{sin}(\theta ).\end{array}$$E16

The same transformation is done with x₀ and y₀ to obtain r₀ and θ₀. The Green's function may be rewritten:

$$p_{\text{free}}\left(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}}\right)=\frac{\text{1}}{\text{4πDt}}\text{exp}\left[-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}{\text{-2rr}}_{\text{0}}\text{cos}(\theta -{\theta }_{\text{0}})}{\text{4Dt}}\right].$$E17

This form yields an alternate method to sample the Green's function in polar coordinates:

Algorithm 5: Sampling algorithm for r*p_free(r,θ,Δt|r₀,θ₀) (free 2D diffusion)

CALCULATE$$\rho \leftarrow r_{\text{0}}/\sqrt{\text{2D}\Delta t}$$
REPEAT
	GENERATE V,W uniform on [0,1], E1 as standard exponential.
	SET$$\Theta \leftarrow \text{2πV}$$.
	IF $$W(\text{1}+\rho \sqrt{\text{2π}})<\text{1}$$THEN SET$$R\leftarrow \rho +\sqrt{\text{2}{E}_{\text{2}}}$$, where E2 is standard exponential
	ELSE SET R ←ρ+N, where N is standard normal
UNTIL R>0 and log(Max(R,ρ)/R)+Rρ(1-cos(Θ-θ0))<E1 .
RETURN (R,Θ).

Algorithm 5 generates variates r and θ distributed as r*p_free(r,θ,t|r₀,θ₀). The factor r originates from the differential area in polar coordinates: dxdy is replaced by rdrdθ.

It is possible to obtain related distributions such as p(r,t|r₀) and p(θ,t|θ₀). They are found by integrating the Green's function:

$$p\left(\text{r,t}|r_{\text{0}}\right)=\frac{\text{1}}{\text{4πDt}}\text{exp}\left(-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}}{\text{4Dt}}\right){\displaystyle \underset{\text{0}}{\overset{\text{2π}}{\int }}\text{exp}\left[\frac{{\text{2rr}}_{\text{0}}\text{cos}\left(\theta -{\theta }_{\text{0}}\right)}{\text{4Dt}}\right]d\theta =}\frac{\text{1}}{\text{2Dt}}{I}_{\text{0}}\left(\frac{{\text{rr}}_{\text{0}}}{\text{2Dt}}\right),$$E18

where I₀(x) is the Modified Bessel Function of the First Kind of order 0 (see ref. [9]). The angular distribution p(θ,t|θ₀) is obtained by integrating over r:

$$p\left(\theta \text{,t}|{\theta }_{\text{0}}\right)=\frac{\text{1}}{\text{4πDt}}{\displaystyle \underset{\text{0}}{\overset{\infty }{\int }}\text{exp}\left[-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}{\text{-2rr}}_{\text{0}}\text{cos}\left(\theta -{\theta }_{\text{0}}\right)}{\text{4Dt}}\right]\text{rdr}}.$$E19

The factor r is also included in the integral. The result is:

$$p\left(\theta \text{,t}|{\theta }_{\text{0}}\right)=\frac{\text{1}}{\text{2π}}\text{exp}\left(-\frac{r_{\text{0}}^{\text{2}}}{\text{4Dt}}\right)+\text{exp}\left[-\frac{r_{\text{0}}^{\text{2}}{\text{sin}}^{\text{2}}\left(\theta -{\theta }_{\text{0}}\right)}{\text{4Dt}}\right]\frac{r_{\text{0}}\text{cos}\left(\theta -{\theta }_{\text{0}}\right)}{\text{4π}\sqrt{\text{Dt}}}\left(\text{1}+\text{Erf}\left[\frac{r_{\text{0}}\text{cos}\left(\theta -{\theta }_{\text{0}}\right)}{\sqrt{\text{4Dt}}}\right]\right).$$E20

where Erf(x) is the error function [9]. The distributions p(r,t|r₀) and p(θ,t|θ₀) for x₀=2.5 and y₀=2.5 are illustrated in Figure 4.

Figure 4.

Distributions p(r,t|r₀) and p(θ,t|θ₀) for 1,000,000 particles initially located at x₀=2.5 and y₀=2.5 at t=1,2,4,8,16 and 32 units.

Similar results are found using either algorithm 4 or 5. In 2D, it is possible to represent the time evolution of the particle system (Figure 5).

Figure 5.

Positions of 1,000,000 particles initially located at x₀=2.5, y₀=0 at t=1, 2, 4, 8, 16 and 32 time units. The positions are obtained by sampling the Green's function as described in Algorithm 4 or 5.

4.2. The rectangular corner x>0, y>0

As in section 4.1, the solution of the DE is the product of two independent 1D reflective Green's functions:

$$\begin{array}{l}p(\text{x,y,t}|x_{\text{0}}{\text{,y}}_{\text{0}})=\\ =\frac{\text{1}}{\text{4}\pi \text{Dt}}\left\{\text{exp}\left[-\frac{{(x-x_0)}^{\text{2}}}{\text{4Dt}}\right]+\text{exp}\left[-\frac{{(x+x_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right]\right\}\left\{\text{exp}\left[-\frac{{(y-y_0)}^{\text{2}}}{\text{4Dt}}\right]+\text{exp}\left[-\frac{{(y+y_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right]\right\}\end{array}.$$E21

The positions of the particles can be sampled by using two applications of Algorithm 2 (for X and Y). An example of the time evolution of 1,000,000 particles is illustrated in Figure 6.

Figure 6.

Positions of 1,000,000 particles initially located at x₀=10, y₀=10 at t=1, 2, 4, 8, 16 and 32 time units, assuming reflective boundaries at X=0 and Y=0. The positions are obtained by sampling the Green's function as described in Algorithm 2.

4.3. The rectangle 0<x<a, 0<y<b

The solution of the DE is also the product of two independent 1D Green's functions:

$$\begin{array}{l}p(\text{x,y,t}|x_{\text{0}}{\text{,y}}_{\text{0}})=\\ =\frac{\text{1}}{\text{ab}}\left[{\displaystyle \sum _{m=-\infty }^{\infty }\text{exp}\left(-\frac{{\pi }^{\text{2}}{m}^{\text{2}}\text{Dt}}{a^{\text{2}}}\right)\text{cos}\frac{m\pi x}{a}\text{cos}\frac{m\pi x_{\text{0}}}{a}}\right]\left[{\displaystyle \sum _{n=-\infty }^{\infty }\text{exp}\left(-\frac{{\pi }^{\text{2}}{n}^{\text{2}}\text{Dt}}{b^{\text{2}}}\right)\text{cos}\frac{n\pi y}{b}\text{cos}\frac{n\pi y_{\text{0}}}{b}}\right]\end{array}.$$E22

Therefore the positions of the particles can be sampled by using Algorithm 3 for the coordinates x and y.

Figure 7.

Positions of 1,000,000 particles initially located at x₀=2.5, y₀=0 in the rectangle of 10 units height and 15 units width at t=1, 2, 4, 8, 16 and 32 time units. The positions are obtained by sampling the Green's function as described in Algorithm 3.

The particle system depicted in Figure 7 illustrates how particles initially located in a given position will diffuse and eventually become uniformly distributed in the rectangular area.

4.4. The particles outside a disk of radius R

This problem is of interest to simulate how particles move and/or bind to curved cell membranes comprising receptors. In this section, only the case of the reflective boundary condition is discussed. In polar coordinates, the DE can be written

$$D\left[\frac{\text{1}}{r}\frac{\partial }{\partial r}\left(r\frac{\partial }{\partial r}\right)+\frac{\text{1}}{r^{\text{2}}}\frac{{\partial }^{\text{2}}}{\partial {\theta }^{\text{2}}}\right]p\left(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}}\right)=\frac{\partial }{\partial t}p\left(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}}\right).$$E23

The reflective boundary condition at r=R is:

$$D\frac{\partial p\left(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}}\right)}{\partial r}{|}_{r=R}=\text{0}.$$E24

The Green's function for a particle initially located at (r₀,θ₀) is:

$$p\left(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}}\right)=\frac{\text{1}}{\text{2π}}{\displaystyle \sum _{n=-\infty }^{\infty }\text{cos}[n(\theta -{\theta }_{\text{0}})]{\displaystyle \underset{\text{0}}{\overset{\infty }{\int }}{e}^{{\text{-u}}^{\text{2}}\text{Dt}}{\text{uC}}_n(\text{u,r})C_n({\text{u,r}}_{\text{0}})\text{du}}},$$E25

where

$$C_n(u,r)=\frac{J_n(ur)Y_n^{\text{'}}(uR)-J_n^{\text{'}}(uR)Y_n(ur)}{\sqrt{J_n^{\text{'}}{(uR)}^{\text{2}}+Y_n^{\text{'}}{(uR)}^{\text{2}}}}.$$E26

J_n(x) and Y_n(x) are the Bessel Functions of the First and Second Kind, respectively (see ref. [9]). The parameters of this distribution are R>0, 0≤θ₀≤2π, r₀>R and Dt>0. The variates of this distributions are r>R and 0≤θ≤2π.

The normalization of this function is:

$${\displaystyle \underset{\text{0}}{\overset{\text{2π}}{\int }}{\displaystyle \underset{R}{\overset{\infty }{\int }}p(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}})\text{rdrd}\theta =\text{1}}}.$$E27

We have not been able so far to develop a sampling algorithm for the distribution function p(r,θ,t|r₀, θ₀). An approximate algorithm to sample the positions can be done by using free diffusion and reflecting particles which are inside the disk (r<R). The results given by this simulation are illustrated in Figure 8.

Figure 8.

Positions of 1,000,000 particles initially located at x₀=2.5, y₀=2.5 near a reflective boundary at R=1 at t=1, 2, 4, 8, 16 and 32 time units.

In this example, since an approximate method is used to simulate the motion of particles, the distributions of r and θ do not verify the time discretization equations, as an exact sampling algorithm of p(r,θ,t|r₀, θ₀) would (results not shown). Therefore, this example illustrate that the time discretization property of the Green’s function is a very powerful way to verify the validity of the sampling algorithms.

It is interesting to note that for R→0, C_n(u,r)→J_n(ur); in this case, the integral can be evaluated by using Weber's formula (Watson P. 395 [12], Gradsteyn 6.633, Equation 2 [13]):

$$p_{\text{ref}}\left(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}}\right)=\frac{\text{1}}{\text{2π}}{\displaystyle \sum _{n=-\infty }^{\infty }\text{cos}[n(\theta -{\theta }_{\text{0}})]\frac{\text{1}}{\text{2Dt}}\text{exp}\left(-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}}{\text{4Dt}}\right)I_n\left(\frac{{\text{rr}}_{\text{0}}}{\text{2Dt}}\right)}.$$E28

The series which is obtained can be evaluated (http://functions.wolfram.com/ 03.02.23.0007.01):

$$p_{\text{ref}}\left(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}}\right)=\frac{\text{1}}{\text{2π}}\frac{\text{1}}{\text{2Dt}}\text{exp}\left(-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}}{\text{4Dt}}\right)\text{exp}\left[\frac{{\text{rr}}_{\text{0}}}{\text{2Dt}}\text{cos}(\theta -{\theta }_{\text{0}})\right]\equiv p_{\text{free}}\left(\text{r,}\theta \text{,t}|r_{\text{0}},{\theta }_{\text{0}}\right).$$E29

Thus, for R→0, p_ref(r,θ,t|r₀,θ₀) → p_free(r,θ,t|r₀,θ₀), as expected.

5. Three-dimensional systems

Three dimensional systems are more complex and harder to visualize. In 3D, the DE can be written in cartesian coordinates as:

$$\frac{\partial }{\partial t}p(\text{x,y,z,t}|x_{\text{0}}{\text{,y}}_{\text{0}}{\text{,z}}_{\text{0}})=D\left[\frac{{\partial }^{\text{2}}}{\partial x^{\text{2}}}+\frac{{\partial }^{\text{2}}}{\partial y^{\text{2}}}+\frac{{\partial }^{\text{2}}}{\partial z^{\text{2}}}\right]p(\text{x,y,z,t}|x_{\text{0}}{\text{,y}}_{\text{0}}{\text{,z}}_{\text{0}}).$$E30

The initial condition is p(x,y,z,t=0|x₀, y₀, z₀)=δ(x-x₀)δ(y-y₀)δ(z-z₀). As in 2D, in several cases, the Green's function can be written as the product of the functions p_x(x,t|x₀), p_y(y,t|y₀) and p_z(z,t|z₀):

$$p(\text{x,y,z,t}|x_{\text{0}}{\text{,y}}_{\text{0}}{\text{,z}}_{\text{0}})=p_x(\text{x,t}|x_{\text{0}})p_y(\text{y,t}|y_{\text{0}})p_z(\text{z,t}|z_{\text{0}}).$$E31

5.1. Free diffusion

The solution of the DE for free particles is the product of three Gaussian functions:

$$\begin{array}{l}{p}_{\text{free}}\left(\text{x,y,z,t}|x_{\text{0}}{\text{, y}}_{\text{0}}{\text{,z}}_{\text{0}}\right)=\\ =\frac{\text{1}}{{(\text{4πDt})}^{\text{3}/\text{2}}}\text{exp}\left[-\frac{{({\text{x-x}}_{\text{0}})}^{\text{2}}+{({\text{y-y}}_{\text{0}})}^{\text{2}}+{({\text{z-z}}_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right]\equiv p_{\text{free}}\left(\text{x,t}|x_{\text{0}}\right)p_{\text{free}}\left(\text{y,t}|y_{\text{0}}\right)p_{\text{free}}\left(\text{z,t}|z_{\text{0}}\right)\end{array}.$$E32

Since the Gaussian functions are independent, the positions x, y and z can be obtained by sampling the Gaussian distribution for each variateX, Y and Z:

Algorithm 6: Sampling algorithm for p_free(x,y,z,Δt|x₀,y_0,z₀) (free 3D diffusion)

GENERATE N1, N2, N3 as a standard Normals

$$X\leftarrow x_{\text{0}}+\sqrt{D\Delta t}{N}_{\text{1}}$$

$$Y\leftarrow y_{\text{0}}+\sqrt{D\Delta t}{N}_{\text{2}}$$

$$Z\leftarrow z_{\text{0}}+\sqrt{D\Delta t}{N}_{\text{3}}$$

RETURN (X,Y,Z).

In 3D, it is useful to use spherical coordinates (r,θ,ϕ)

The following notation is used: (radial, azimutal, polar). r: [0,∞),θ: [0,π],ϕ: [0,2π).

. That is, the following transformation is used:

$$\begin{array}{l}x=r\text{sin}(\theta )\text{cos}(\varphi ),\\ y=r\text{sin}(\theta )\text{sin}(\varphi ),\\ z=r\text{cos}(\theta ).\end{array}$$E33

The differential volume element is given by

$$dV=r^{\text{2}}\text{sin}(\theta )drd\theta d\varphi .$$E34

The DE can be written in spherical coordinates as:

$$\frac{\partial p\left(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}}\right)}{\partial t}=D\left[\frac{{\partial }^{\text{2}}}{\partial r^{\text{2}}}+\frac{\text{2}}{r}\frac{\partial }{\partial r}+\frac{\text{1}}{r^{\text{2}}{\text{sin}}^{\text{2}}\theta }\frac{{\partial }^{\text{2}}}{\partial {\varphi }^{\text{2}}}+\frac{\text{1}}{r^{\text{2}}\text{sin}\theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial }{\partial \theta }\right)\right]p\left(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}}\right).$$E35

Using the transformation to spherical coordinates, the Green's function of the DE for free particles can be rewritten as:

$$\begin{array}{l}{p}_{\text{free}}(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}})=\\ =\frac{\text{1}}{{(\text{4πDt})}^{\text{3}/\text{2}}}\text{exp}\left\{-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}-{\text{2rr}}_{\text{0}}\left[\text{cos}\theta \text{cos}{\theta }_{\text{0}}+\text{sin}\theta \text{sin}{\theta }_{\text{0}}\text{cos}(\varphi -{\varphi }_{\text{0}})\right]}{\text{4Dt}}\right\}\end{array}.$$E36

The angle between the vectors (r,θ,ϕ) and (r₀,θ₀,ϕ₀) is given by cos(γ):

$$\text{cos}(\gamma )=\cos \theta \cos {\theta }_0+\sin \theta \sin {\theta }_0\cos (\varphi -{\varphi }_0).$$E37

Therefore, Eq. (36) can be written in the compact form:

$$p_{\text{free}}\left(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}}\right)=\frac{\text{1}}{{(\text{4πDt})}^{\text{3}/\text{2}}}\text{exp}\left[-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}-{\text{2rr}}_{\text{0}}\text{cos}\gamma }{\text{4Dt}}\right].$$E38

In spherical coordinates, as in 2D, it is possible to obtain related distributions such as p(r,t|r₀). It is obtained by integrating the Green's function over the angles:

$$p\left(\text{r,t}|r_{\text{0}}\right)={\displaystyle \underset{\text{0}}{\overset{\text{2π}}{\int }}{\displaystyle \underset{\text{0}}{\overset{\pi }{\int }}{r}^{\text{2}}}}{p}_{\text{free}}\left(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}}\right)\text{sin}(\varphi )d\theta d\varphi =\frac{r}{\sqrt{{\text{4πr}}_{\text{0}}^{\text{2}}\text{Dt}}}\text{exp}\left(-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}}{\text{4Dt}}\right)\text{sinh}\left(\frac{{\text{rr}}_{\text{0}}}{\text{2Dt}}\right).$$E39

5.2. The corner x>0, y>0, z>0

The Green's function is the product of three independent 1D Green's functions:

$$\begin{array}{l}p\left(\text{x,y,z,t}|x_{\text{0}}{\text{,y}}_{\text{0}}{\text{,z}}_{\text{0}}\right)=\frac{\text{1}}{{(\text{4πDt})}^{\text{3/2}}}\left[e^{-{({\text{x-x}}_{\text{0}})}^{\text{2}}/\text{4Dt}}+e^{-{(x+x_{\text{0}})}^{\text{2}}/\text{4Dt}}\right]\\ \left[e^{-{({\text{y-y}}_{\text{0}})}^{\text{2}}/\text{4Dt}}+e^{-{(y+y_{\text{0}})}^{\text{2}}/\text{4Dt}}\right]\left[e^{-{({\text{z-z}}_{\text{0}})}^{\text{2}}/\text{4Dt}}+e^{-{(z+z_{\text{0}})}^{\text{2}}/\text{4Dt}}\right]\end{array}.$$E40

As it is the case in 2D, the coordinates x,y and z can be treated independently. Therefore, the positions of the particles can be generated by using Algorithm 2 for each of the coordinates.

5.3. The parallelipiped 0<x<a, 0<y<b, 0<z<c

The Green's function is also the product of three independent 1D Green's functions:

$$\begin{array}{l}p(\text{x,y,z,t}|x_{\text{0}}{\text{,y}}_{\text{0}}{\text{,z}}_{\text{0}})=\\ \frac{\text{1}}{\text{abc}}{\displaystyle \sum _{\text{m,n,l}=-\infty }^{\infty }\text{exp}\left[-{\pi }^{\text{2}}\text{Dt}\left(\frac{m^{\text{2}}}{a^{\text{2}}}+\frac{n^{\text{2}}}{b^{\text{2}}}+\frac{l^{\text{2}}}{c^{\text{2}}}\right)\right]\text{cos}\frac{m\pi x}{a}\text{cos}\frac{m\pi x_{\text{0}}}{a}}\text{cos}\frac{n\pi y}{b}\text{cos}\frac{n\pi y_{\text{0}}}{b}\text{cos}\frac{l\pi z}{c}\text{cos}\frac{l\pi z_{\text{0}}}{c}\end{array}.$$E41

Therefore, the positions of particles can be generated by using Algorithm 3 for each of the coordinates. This sampling algorithm can be used to simulate the trajectories of particles confined into a box of sides a, b and c. At equilibrium, as it is the case in 2D (Figure 7), the particles are uniformly distributed within the volume.

6. Applications

Two applications of the approach described in the previous sections are presented. The first is the simulation of particles near a membrane comprising receptors. The second one is the simulation of chemical reactions. These applications will eventually be used to study the evolution of the chemical species created by ionizing radiations in biological media [14].

6.1. Particle near a 1D membrane comprising receptors

To apply the theory described in the chapter to cell culture simulations, it is necessary to make some assumptions. The details of the model presented here and more results are given in ref. [5]. Similar models have been used to study autocrine and paracrine signals in cell culture assays [15], ligand accumulation in autocrine cell cultures [16] and the spatial range of autocrine signaling [17].

6.1.1. Model and assumptions

A particle is assumed to be located in the proximity of a plane membrane comprising receptors, corresponding to the plane X=0. The particle can move freely in the Y and Z directions, and can bind reversibly to receptors located on the membrane and initiate signal transduction. This is illustrated in Figure 9:

Figure 9.

Ligands near a surface comprising receptors. The rate constants for association, dissociation and initiation of signal transduction are k_a, k_d and k_e, respectively.

The binding rate constant (k_a) is obtained from the ligand-receptor association rate constant and the density of receptors on the surface, which are assumed to be uniformly distributed on the membrane [16]. The ligand-receptor complex has the rate constants k_d for dissociation and k_e for initiation of signal transduction. As a first approximation, this is equivalent to the problem

$$A+B\underset{k_d}{\overset{k_a}{\rightleftarrows }}{(AB)}^{}\stackrel{k_e}{\to }AB.$$E42

Since the particles are moving freely in the Y and Z direction, the Green's function can be written as the product of the independent functions p_x(x,t|x₀), p_y(y,t|y₀) and p_z(z,t|z₀). But p_y(y,t|y₀) and p_z(z,t|z₀) are Gaussians functions, which can be sampled as described in section 4.1; therefore, only the sampling of p_x(x,t|x₀) will be discussed here. For convenience, the index "x" will be omitted in this section.

A particle can be free at position x>0, in a reversible bound state (*) at x=0 or activating signal transduction (**). The boundary condition at x=0 can be written:

$$D\frac{\partial p(\text{x,t}|x_{\text{0}})}{\partial x}{|}_{x=\text{0}}=k_{a}p(\text{0,}t|x_{\text{0}})-k_{d}p(\text{*,}t|x_{\text{0}}),$$E43

where p(*,t|x₀) is the probability of a particle initially located at x₀ to be in the reversible bound state at time t. The time evolution of p(*,t|x₀) is:

$$\frac{dp(\text{*,t}|x_{\text{0}})}{dt}=k_{a}p(\text{0,}t|x_{\text{0}})-(k_d+k_e)p(\text{*,}t|x_{\text{0}}).$$E44

Similarly, the time evolution of p(**,t|x₀) is:

$$\frac{dp(\text{**,t}|x_{\text{0}})}{dt}=k_{e}p(\text{*,}t|x_{\text{0}}).$$E45

6.1.2. Green's functions

This problem can be solved analytically by using Laplace Transforms [18]. The solutions are expressed using three coefficients α, β and γ. They are the roots of a cubic polynomial

At least one of the roots of a third-order polynomial is real, the two other roots being either real or complex conjugates.

, for which the coefficients are given by the rate constants k_a, k_d and k_e as follows:

$$\begin{array}{l}\alpha +\beta +\gamma =k_a/D,\\ \alpha \beta +\beta \gamma +\gamma \alpha =\left(k_e+k_d\right)/D,\\ \alpha \beta \gamma =k_e{k}_a/D^{\text{2}}.\end{array}$$E46

The mathematical expressions of the Green's function for a particle in this system are quite long. To simplify the following expressions, the variables $${\rho }_x=\text{x/}\sqrt{\text{4Dt}}$$, $${\rho }_{\text{0}}=x_{\text{0}}/\sqrt{\text{4Dt}}$$, $$\alpha \text{'}=\alpha \sqrt{\text{Dt}}$$, $$\beta \text{'}=\beta \sqrt{\text{Dt}}$$ and $$\gamma \text{'}=\gamma \sqrt{\text{Dt}}$$ are defined.

$$\begin{array}{l}p\left(\text{x,t}|x_{\text{0}}\right)=p_{\text{ref}}\left(\text{x,t}|x_{\text{0}}\right)+\frac{\alpha (\gamma +\alpha )(\alpha +\beta )}{(\gamma -\alpha )(\alpha -\beta )}W\left({\rho }_x+{\rho }_{\text{0}},\alpha \text{'}\right)+\\ +\frac{\beta (\alpha +\beta )(\beta +\gamma )}{(\alpha -\beta )(\beta -\gamma )}W\left({\rho }_x+{\rho }_{\text{0}},\beta \text{'}\right)+\frac{\gamma (\beta +\gamma )(\gamma +\alpha )}{(\beta -\gamma )(\gamma -\alpha )}W\left({\rho }_x+{\rho }_{\text{0}},\gamma \text{'}\right)\end{array}.$$E47

where p_ref(x,t|x₀) is the Green's function for a reflective membrane, given in Section 3.2. The functions W(x,y) and Ω(x) are defined as:

$$\begin{array}{l}\Omega (x)=\text{exp}(x^{\text{2}})\text{Erfc}(x),\\ W(\text{x,y})\equiv \text{exp}(\text{2xy}+y^{\text{2}})\text{Erfc}(x+y)=\text{exp}(-x^{\text{2}})\Omega (x+y)\end{array}$$E48

It is possible for two of the roots to be complex conjugates. Nevertheless, it is still possible to use the Green's function given in Equation (47) by implementing the numerical evaluation of W(x,y) for complex numbers [19].

The probability of survival for a free particle is obtained by integrating p(x,t|x₀) on [0,∞):

$$\begin{array}{l}Q\left(t|x_{\text{0}}\right)=\text{1}-\frac{(\gamma +\alpha )(\alpha +\beta )}{(\gamma -\alpha )(\alpha -\beta )}W\left({\rho }_{\text{0}},\alpha \text{'}\right)-\\ -\frac{(\alpha +\beta )(\beta +\gamma )}{(\alpha -\beta )(\beta -\gamma )}W\left({\rho }_{\text{0}},\beta \text{'}\right)-\frac{(\beta +\gamma )(\gamma +\alpha )}{(\beta -\gamma )(\gamma -\alpha )}W\left({\rho }_{\text{0}},\gamma \text{'}\right)-\text{Erfc}\left({\rho }_{\text{0}}\right)\end{array}.$$E49

The probability of a particle initially at x₀>0 to be reversibly bound at time t is:

$$p\left(\text{*,t}|x_{\text{0}}\right)=\frac{\alpha (\alpha +\beta +\gamma )}{(\gamma -\alpha )(\alpha -\beta )}W\left({\rho }_{\text{0}},\alpha \text{'}\right)+\frac{\beta (\alpha +\beta +\gamma )}{(\alpha -\beta )(\beta -\gamma )}W\left({\rho }_{\text{0}},\beta \text{'}\right)+\frac{\gamma (\alpha +\beta +\gamma )}{(\beta -\gamma )(\gamma -\alpha )}W\left({\rho }_{\text{0}},\gamma \text{'}\right).$$E50

A particle initially at x₀>0 can also be found in an irreversibly bound state at time t. The probability to be found in this state is given by p(**,t|x₀)=1-p(*,t|x₀)-Q(t|x₀), which can be written as:

$$\begin{array}{l}p\left(\text{**,t}|x_{\text{0}}\right)=\frac{\beta \gamma }{(\gamma -\alpha )(\alpha -\beta )}W\left({\rho }_{\text{0}},\alpha \text{'}\right)+\\ +\frac{\alpha \gamma }{(\alpha -\beta )(\beta -\gamma )}W\left({\rho }_{\text{0}},\beta \text{'}\right)+\frac{\alpha \beta }{(\beta -\gamma )(\gamma -\alpha )}W\left({\rho }_{\text{0}},\gamma \text{'}\right)+\text{Erfc}\left({\rho }_{\text{0}}\right)\end{array}.$$E51

To obtain the Green's function of a particle in a reversibly bound state, the material balance condition k_ap(x,t|*)=k_dp(*,t|x) (given in ref. [18,20]) is used. This yields the Green's function for the dissociation of a bound particle:

$$\frac{D}{k_d}p\left(\text{x,t}|*\right)=\frac{\alpha }{(\gamma -\alpha )(\alpha -\beta )}W\left({\rho }_x,\alpha \text{'}\right)+\frac{\beta }{(\alpha -\beta )(\beta -\gamma )}W\left({\rho }_x,\beta \text{'}\right)+\frac{\gamma }{(\beta -\gamma )(\gamma -\alpha )}W\left({\rho }_x,\gamma \text{'}\right).$$E52

The probability of dissociation of a particle that was initially in a reversibly bound state Q(t|*) is found by integrating p(*,t|x) over [0,∞):

$$\frac{D}{k_d}Q\left(t|*\right)=-\frac{\text{1}}{(\gamma -\alpha )(\alpha -\beta )}\Omega \left(\alpha \text{'}\right)-\frac{\text{1}}{(\alpha -\beta )(\beta -\gamma )}\Omega \left(\beta \text{'}\right)-\frac{\text{1}}{(\beta -\gamma )(\gamma -\alpha )}\Omega \left(\gamma \text{'}\right).$$E53

The probability of a particle initially in a reversibly bound state to stay in this state is given by (ref. [18]):

$$p\left(\text{*,t}|*\right)=\frac{\alpha (\beta +\gamma )}{(\gamma -\alpha )(\alpha -\beta )}\Omega \left(\alpha \text{'}\right)+\frac{\beta (\gamma +\alpha )}{(\alpha -\beta )(\beta -\gamma )}\Omega \left(\beta \text{'}\right)+\frac{\gamma (\alpha +\beta )}{(\beta -\gamma )(\gamma -\alpha )}\Omega \left(\gamma \text{'}\right).$$E54

Finally, the probability of a particle initially in a reversibly bound state to activate signal transduction is:

$$p\left(\text{**,t}|*\right)=\text{1}+\frac{\text{1}}{\alpha +\beta +\gamma }\left[\frac{\beta \gamma (\beta +\gamma )}{(\gamma -\alpha )(\alpha -\beta )}\Omega \left(\alpha \text{'}\right)+\frac{\alpha \gamma (\alpha +\gamma )}{(\alpha -\beta )(\beta -\gamma )}\Omega \left(\beta \text{'}\right)+\frac{\alpha \beta (\alpha +\beta )}{(\beta -\gamma )(\gamma -\alpha )}\Omega \left(\gamma \text{'}\right)\right].$$E55

For some particular values of α, β and γ, the Green's functions of this section have singularities (0/0). In these cases, the Green's functions take slightly different forms, but it is still possible to use the approach described here (see ref. [5]).

6.1.3. Time discretization equations

A simulation is usually done by dividing the final simulation time in a finite number of time steps. If the simulation can be done in two time steps such as t=Δt₁+Δt₂, a particle initially located at position x₀ can:

1. go to position x₁ during Δt₁ and then go to its final x position during Δt₂,

2. go to position x₁ during Δt₁ and reversibly bind to the membrane during Δt₂,

3. go to position x₁ during Δt₁ and activate signal transduction during Δt₂,

4. reversibly bind to the membrane during Δt₁ and stay in this state during Δt₂,

5. reversibly bind to the membrane during Δt₁ and dissociate during Δt₂,

6. reversibly bind to the membrane during Δt₁ and activate signal transduction during Δt₂,

or

1. activate signal transduction during Δt₁.

From this the following time discretization equations are obtained:

$$\begin{array}{l}p\left(\text{x,t}|x_{\text{0}}\right)={\displaystyle \underset{\Omega }{\int }p\left(\text{x,}\Delta t_{\text{2}}|x_{\text{1}}\right)p\left(x_{\text{1}},\Delta t_{\text{1}}|x_{\text{0}}\right)}d{x}_1+p\left(x,\Delta t_{\text{2}}|*\right)p\left(\text{*,}\Delta t_{\text{1}}|x_{\text{0}}\right),\\ p\left(\text{*,t}|x_{\text{0}}\right)=p\left(\text{*,}\Delta t_{\text{2}}|*\right)p\left(\text{*,}\Delta t_{\text{1}}|x_{\text{0}}\right)+{\displaystyle \underset{\Omega }{\int }p\left(\text{*,}\Delta t_{\text{2}}|x_{\text{1}}\right)p\left(x_{\text{1}},\Delta t_{\text{1}}|x_{\text{0}}\right)}d{x}_1,\\ p\left(\text{**,t}|x_{\text{0}}\right)=p\left(\text{**,}\Delta t_{\text{1}}|x_{\text{0}}\right)+p\left(\text{**,}\Delta t_{\text{2}}|*\right)p\left(\text{*,}\Delta t_{\text{1}}|x_{\text{0}}\right)+{\displaystyle \underset{\Omega }{\int }p\left(\text{**,}\Delta t_{\text{2}}|x_{\text{1}}\right)p\left(x_{\text{1}},\Delta t_{\text{1}}|x_{\text{0}}\right)}d{x}_1.\end{array}$$E56

The first of these equations is similar to the Chapman-Kolmogorov equation, to which a term has been added to account for the possibility of dissociation. Each term of these equations can be associated with the possibilities i-vii. Similarly, a particle initially in a reversibly bound state can:

1. dissociate to position x₁ during Δt₁ and go to its final position x during Δt₂,

2. dissociate to position x₁during Δt₁ and re-bind reversibly during Δt₂,

3. dissociate to position x₁during Δt₁ and initiate signal transduction during Δt₂,

4. stay reversibly bound during Δt₁ and dissociate to position x during Δt₂,

5. stay reversibly bound during Δt₁ and Δt₂,

6. stay reversibly bound during Δt₁ and initiate signal transduction during Δt₂, or

7. initiate signal transduction during Δt₁.

This yields the time discretization equations for the bound particle:

$$\begin{array}{l}p\left(\text{x,t}|*\right)={\displaystyle \underset{\Omega }{\int }p\left(\text{x,}\Delta t_{\text{2}}|x_{\text{1}}\right)p\left(x_{\text{1}},\Delta t_{\text{1}}|*\right)}d{x}_1+p\left(x,\Delta t_{\text{2}}|*\right)p\left(*,\Delta t_{\text{1}}|*\right),\\ p\left(*\text{,t}|*\right)=p\left(*,\Delta t_{\text{2}}|*\right)p\left(*,\Delta t_{\text{1}}|*\right)+{\displaystyle \underset{\Omega }{\int }p\left(*,\Delta t_{\text{2}}|x_{\text{1}}\right)p\left(x_{\text{1}},\Delta t_{\text{1}}|*\right)}d{x}_1,\\ p\left(*\text{,t}|*\right)=p\left(*,\Delta t_{\text{1}}|*\right)+p\left(*,\Delta t_{\text{2}}|*\right)p\left(*,\Delta t_{\text{1}}|*\right)+{\displaystyle \underset{\Omega }{\int }p\left(*,\Delta t_{\text{2}}|x_{\text{1}}\right)p\left(x_{\text{1}},\Delta t_{\text{1}}|*\right)}d{x}_1.\end{array}$$E57

These equations are proven in the supplemental material of ref. [5].

6.1.4. Brownian dynamics algorithm

Using the Green's functions, it is possible to construct a Brownian Dynamics algorithm to simulate the motion and binding of the ligand and the possible initiation of signal transduction. First, it should be noted that p(x,t|x₀) is a sub-density

A sub-density means that $$P={\displaystyle \underset{\text{0}}{\overset{L}{\int }}f(x)dx\le \text{1}}$$, and f(x)≥0. In this case, a random variate with the density f/p is generated with probability P.

. Since p(x,t|x₀)≤p_ref(x,t|x₀), the rejection method can be used [21]. This yields the following algorithm:

Algorithm 7a: sampling of p(x,Δt|x₀) (free particles)

REPEAT
{
	GENERATE U uniform on [0,1], X distributed as pref(x,Δt|x0)
}
UNTIL p(X,Δt|x0)≤U*pref(X,Δt|x0)
RETURN X

There are also particles which are bound to the membrane. These particles may dissociate with probability Q(Δt|*). If this is the case, the position x of the particle is distributed as p(x,Δt|*), which can be written with a Gaussian factor in evidence:

$$\begin{array}{l}p\left(\text{x,}\Delta t|*\right)=\frac{k_d}{D}\text{exp}\left(-\frac{x^{\text{2}}}{\text{4DΔt}}\right)[\frac{\alpha }{(\gamma -\alpha )(\alpha -\beta )}\Omega \left(\frac{x}{\sqrt{\text{4DΔt}}}+\alpha \sqrt{D\Delta t}\right)\\ +\frac{\beta }{(\alpha -\beta )(\beta -\gamma )}\Omega \left(\frac{x}{\sqrt{\text{4DΔt}}}+\beta \sqrt{D\Delta t}\right)+\frac{\gamma }{(\beta -\gamma )(\gamma -\alpha )}\Omega \left(\frac{x}{\sqrt{\text{4DΔt}}}+\gamma \sqrt{D\Delta t}\right)]\end{array}.$$E58

This form is convenient because it is the product of a Gaussian function by a function g(x) with three terms in Ω(x). Since 0≤Ω(x)≤1, p(x,Δt|*) can be sampled by a rejection method:

Algorithm 7b: sampling of p(x,Δt|*) (reversibly bound particles)

CALCULATE H = |α/(γ-α)(α-β)|+|β/(α-β)(β-γ)|+|γ/(β-γ)(γ-α)|
REPEAT
{
	GENERATE U uniform on [0,1], $$X=\left|N\right|\sqrt{\text{2D}\Delta t}$$, where N is a standard Normal
}
UNTIL UH≤g(X)
RETURN X

6.1.5. Simulation results

A simulation has been done by using k_a=6, k_d=10 and k_e=1. The results are shown in Figure 10. It should also be noted that the time discretization equations (56) and (57) were used to validate the simulation results. In this case, it is possible to verify that the simulation results for free particles corresponds to p(x,t|x₀), but also that the number of bound particles can be accurately predicted after several time steps.

Figure 10.

Left) Probability distribution of a particle near a partially absorbing and reflecting boundary with back reaction for x₀=2.5, k_a=6, k_d=10, k_e=1 and D=1 (α=1, β=2 and γ=3) at t=1, 2, 4, 8 and 16 time units. The lines are the analytical predictions p(x,t|x₀). The dots are given by the simulation of 1,000,000 particle histories either with one or multiple time steps. (Right) Probability of a particle in these conditions to be free, reversibly bound or irreversibly bound as a function of time.

6.2. Bimolecular reactions

Bimolecular reactions are an important application of the theory described in the book chapter. Let the particles A and B initially at the positions r_A0 and r_B0. They can react with the reaction rate k_a:

$$A+B\stackrel{k_a}{\to }C.$$E59

6.2.1. Green's functions for the particles

The probability distribution is solution of the following equation [2]:

$$\frac{\partial }{\partial t}p({\text{r}}_{\text{A}},{\text{r}}_{\text{B}}\text{,t}|{\text{r}}_{{\text{A}}_{\text{0}}},{\text{r}}_{{\text{B}}_{\text{0}}})=\left[D_A{\nabla }_A^2+D_B{\nabla }_B^2\right]p({\text{r}}_{\text{A}},{\text{r}}_{\text{B}}\text{,t}|{\text{r}}_{{\text{A}}_{\text{0}}},{\text{r}}_{{\text{B}}_{\text{0}}}).$$E60

It is convenient at this point to use a transformation of coordinates, similar to a center of mass transformation:

$$\begin{array}{l}\text{R}=\sqrt{D_B/D_A}{\text{r}}_{\text{A}}+\sqrt{D_A/D_B}{\text{r}}_{\text{B}},\\ \text{r}={\text{r}}_{\text{B}}-{\text{r}}_{\text{A}}.\end{array}$$E61

The same transformation is done to r_A0 and r_B0 to obtain r₀ and R₀. The operators ∇_R=∂/∂R and ∇r=∂/∂r are also defined. Using the transformation, Equation (60) can be rewritten as:

$$\frac{\partial }{\partial t}p(\text{R},\text{r}\text{,t}|{\text{R}}_{\text{0}},{\text{r}}_{\text{0}})=D\left[{\nabla }_{\text{R}}^2+{\nabla }_{\text{r}}^2\right]p(\text{R},\text{r}\text{,t}|{\text{R}}_{\text{0}},{\text{r}}_{\text{0}}).$$E62

where D=D_A+D_B, the sum of the diffusion coefficients of the particles. Therefore, p(R,r,t|R₀,r₀) can be factorized as follows:

$$p(\text{R},\text{r}\text{,t}|{\text{R}}_{\text{0}},{\text{r}}_{\text{0}})=p^{\text{R}}(\text{R}\text{,t}|{\text{R}}_{\text{0}})p^{\text{r}}(\text{r}\text{,t}|,{\text{r}}_{\text{0}}).$$E63

Using this factorization, two DE are obtained in r and R:

$$\begin{array}{l}\frac{\partial }{\partial t}{p}^{\text{R}}(\text{R}\text{,t}|{\text{R}}_{\text{0}})=D{\nabla }_{\text{R}}^2{p}^{\text{R}}(\text{R}\text{,t}|{\text{R}}_{\text{0}}),\\ \frac{\partial }{\partial t}{p}^{\text{r}}(\text{r}\text{,t}|{\text{r}}_{\text{0}})=D{\nabla }_{\text{r}}^{\text{2}}{p}^{\text{r}}(\text{r}\text{,t}|{\text{r}}_{\text{0}}).\end{array}$$E64

This equation describes two independent random processes - diffusion in the coordinates r and R. The initial condition for Equation (64a) is p^R(R,t|R₀)=δ(R-R₀), and the boundary condition is p^R(|R|→∞,t|R₀)=0. As seen in Section 5.1, Equation (64a) represents free diffusion in 3D:

$$p^{\text{R}}\left(\text{R}\text{,t}|{\text{R}}_{\text{0}}\right)=\frac{\text{1}}{{(\text{4}\pi \text{Dt})}^{\text{3/2}}}\text{exp}\left[\frac{{(\text{R}-{\text{R}}_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right].$$E65

This function can be sampled for R as in Algorithm 3, by using three Gaussian random numbers. Equation (64b) is the diffusive motion of the separation vector between particles. The initial condition is p^r(r,t|r₀,t₀)=δ(r-r₀), and the outer boundary condition is p^r(|r|→∞,t|r₀,t₀)=0. The possibility of a chemical reaction is introduced as the inner boundary condition

$$\text{4}\pi {\sigma }^{2}D\frac{\partial p^{\text{r}}(\text{r}\text{,t}|{\text{r}}_{\text{0}})}{\partial \text{r}}{|}_{\text{r}=\sigma }=k_a{p}^{\text{r}}(|\text{r}|=\sigma \text{,t}|{\text{r}}_{\text{0}}),$$E66

where σ is the reaction radius. As it is the case for free particles, Equation (66) is the radiation boundary condition. The limit k_a→0 is the reflecting boundary condition (no reactions between particles), whereas the limit k_a→∞ is the absorbing boundary condition. In the latter case, the boundary condition is written p^r(|r|=σ,t|r₀)=0.

The exact solution for this problem is given in ref. [4]:

$$p\left(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}}\right)=\frac{\text{1}}{\text{4π}\sqrt{{\text{rr}}_{\text{0}}}}{\displaystyle \sum _{n=\text{0}}^{\infty }(\text{2n}+\text{1})P_n(\text{cos}\gamma ){\displaystyle \underset{\text{0}}{\overset{\infty }{\int }}{e}^{{\text{-u}}^{\text{2}}\text{Dt}}\text{u}{F}_{n+\text{1/2}}(\text{u,r})F_{n+\text{1/2}}({\text{u,r}}_{\text{0}})\text{du}}}.$$E67

where P_n(x) are the Legendre polynomials,

$$F_{\upsilon }(\text{u,r})=\frac{({\text{2σk}}_a+\text{1})\left[J_{\upsilon }(\text{ur})Y_{\upsilon }(u\sigma )-Y_{\upsilon }(\text{ur})J_{\upsilon }(u\sigma )\right]-\text{2uσ}\left[J_{\upsilon }(\text{ur})Y_{\nu }^{\text{'}}(u\sigma )-Y_{\upsilon }(\text{ur})J_{\nu }^{\text{'}}(u\sigma )\right]}{\sqrt{{\left[({\text{2σk}}_a+\text{1})J_{\upsilon }(u\sigma )-{\text{2uσJ}}_{\nu }^{\text{'}}(u\sigma )\right]}^{\text{2}}+{\left[({\text{2σk}}_a+\text{1})Y_{\upsilon }(u\sigma )-{\text{2uσY}}_{\nu }^{\text{'}}(u\sigma )\right]}^{\text{2}}}}.$$E68

and γ is the angle formed by (r,θ,ϕ), the origin and (r₀,θ₀,ϕ₀). In this expression, it is possible to evaluate P_n(cosγ), knowing that

$$P_n(\text{cos}\gamma )=P_n(\text{cos}\theta )P_n(\text{cos}{\theta }_{\text{0}})+\text{2}{\displaystyle \sum _{m=\text{1}}^n\frac{(\text{n-m})!}{(n+m)!}{P}_n^m(\text{cos}\theta )P_n^m(\text{cos}{\theta }_{\text{0}})\text{cos}m(\varphi -{\varphi }_{\text{0}})}.$$E69

where P_n^m(x) are the associated Legendre polynomials (ref. [9]). For θ₀=0 and ϕ₀=0, P_n(cosγ)=P_n(cosθ). It is also interesting to note that for k_a=0 and σ→0, F_n+1/2(u,r) →J_n+1/2(ur). In this case, p(r,θ,ϕ,t|r₀,,θ₀,ϕ₀) can be written:

$$p\left(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}}\right)=\frac{\text{1}}{\text{4π}\sqrt{{\text{rr}}_{\text{0}}}}{\displaystyle \sum _{n=\text{0}}^{\infty }(\text{2n}+\text{1})P_n(\text{cos}\gamma ){\displaystyle \underset{\text{0}}{\overset{\infty }{\int }}{e}^{{\text{-u}}^{\text{2}}\text{Dt}}\text{u}{J}_{n+\text{1/2}}(\text{ur})J_{n+\text{1/2}}({\text{ur}}_{\text{0}})\text{du}}}.$$E70

The integral can be evaluated using Weber's formula (ref. [13], section 6.633, eq 2):

$$p\left(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}}\right)=\frac{\text{1}}{\text{4π}\sqrt{{\text{rr}}_{\text{0}}}}{\displaystyle \sum _{n=\text{0}}^{\infty }(\text{2n}+\text{1})P_n(\text{cos}\gamma )\frac{\text{1}}{\text{2Dt}}}\text{exp}\left(-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}}{\text{4Dt}}\right)I_{n+\text{1/2}}\left(\frac{{\text{rr}}_{\text{0}}}{\text{2Dt}}\right).$$E71

This equation can be further simplified by using $$I_{\upsilon }(z)=e^{-i\upsilon \pi /\text{2}}{J}_{\upsilon }(e^{i\pi /\text{2}}z)$$ (ref. [12], p. 77; ref. [13] section 8.406), $$J_{\upsilon }(e^{im\pi }z)=e^{im\pi \upsilon }{J}_{\upsilon }(z)$$ (ref. [12], p. 75; ref. [13], section 8.476) and the identity

$$\text{exp}\left(\text{izcos}\gamma \right)=\sqrt{\frac{\pi }{\text{2z}}}{\displaystyle \sum _{n=\text{0}}^{\infty }(\text{2n}+\text{1})i^n{J}_{n+\text{1/2}}(z)P_n(\text{cos}\gamma )}.$$E72

(ref. [12], p. 368; ref. [13], section 8.511). This yields the following equation:

$$\begin{array}{l}p(\text{r,}\theta ,\varphi \text{,t}|r_{\text{0}},{\theta }_{\text{0}},{\varphi }_{\text{0}})=\\ =\frac{\text{1}}{\text{8πDt}\sqrt{{\text{rr}}_{\text{0}}}}\text{exp}\left(-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}}{\text{4Dt}}\right)\sqrt{\frac{{\text{rr}}_{\text{0}}}{\pi \text{Dt}}}\text{exp}\left(\frac{{\text{rr}}_{\text{0}}\text{cos}\gamma }{\text{2Dt}}\right)\equiv \frac{\text{1}}{{(\text{4πDt})}^{\text{3}/\text{2}}}\text{exp}\left[-\frac{r^{\text{2}}+r_{\text{0}}^{\text{2}}-{\text{2rr}}_{\text{0}}\text{cos}\gamma }{\text{4Dt}}\right]\end{array}.$$E73

Therefore, p(r,θ,ϕ,t|r₀,θ₀,ϕ₀) is equivalent to free 3D diffusion (Equation 38) for k_a=0 and σ=0.

We were not able to develop a sampling algorithm for p(r,θ,ϕ,t|r₀,θ₀,ϕ₀). However, it is possible to use the radial component of the DE in 3D:

$$\frac{\partial p(\text{r,t}|r_{\text{0}})}{\partial t}=\frac{D}{r^{\text{2}}}\frac{\partial }{\partial r}\left[r^{\text{2}}\frac{\partial }{\partial r}p(\text{r,t}|r_{\text{0}})\right],$$E74

with the boundary condition given by Equation 66. The solution is much simpler in this case:

$${\text{4πrr}}_{\text{0}}p\left(\text{r,t}|r_{\text{0}}\right)=\frac{\text{1}}{\sqrt{\text{4πDt}}}\left\{\text{exp}\left[-\frac{{(r-r_0)}^{\text{2}}}{\text{4Dt}}\right]+\text{exp}\left[-\frac{{(r+r_{\text{0}}-\text{2σ})}^{\text{2}}}{\text{4Dt}}\right]\right\}+\alpha W\left(\frac{r+r_{\text{0}}-\text{2σ}}{\sqrt{\text{4Dt}}},-\alpha \sqrt{\text{Dt}}\right),$$E75

where α=-(k_a+4πσD)/(4πσ²D). Therefore, this approximation is widely used in chemistry and biophysics [22]. The survival probability of a pair of particles, Q(t|r₀), is well known and can be calculated by integrating p(r,t|r₀) over space:

$$Q\left(t|r_{\text{0}}\right)={\displaystyle \underset{\sigma }{\overset{\infty }{\int }}{\text{4πr}}^{\text{2}}p\left(\text{r,t}|r_{\text{0}}\right)\text{dr}}=\text{1}+\frac{\sigma \alpha +\text{1}}{r_{\text{0}}\alpha }\left[W\left(\frac{r_{\text{0}}-\sigma }{\sqrt{\text{4}Dt}},\alpha \sqrt{Dt}\right)-Erfc\left(\frac{r_{\text{0}}-\sigma }{\sqrt{\text{4}Dt}}\right)\right].$$E76

For this system, the time discretization equations takes the form:

$$\begin{array}{l}p\left(r_{\text{2}},\Delta t_{\text{1}}+\Delta t_{\text{2}}|r_{\text{0}}\right)={\displaystyle \underset{\sigma }{\overset{\infty }{\int }}{\text{4πr}}_{\text{1}}^{\text{2}}}p\left(r_{\text{2}},\Delta t_{\text{2}}|r_{\text{1}}\right)p\left(r_{\text{1}},\Delta t_{\text{1}}|r_{\text{0}}\right)\text{dr},\\ p\left(\text{*,}\Delta t_{\text{1}}+\Delta t_{\text{2}}|r_{\text{0}}\right)=p\left(\text{*,}\Delta t_{\text{1}}|r_{\text{0}}\right)+{\displaystyle \underset{\sigma }{\overset{\infty }{\int }}{\text{4πr}}_{\text{1}}^{\text{2}}}p\left(\text{*,}\Delta t_{\text{2}}|r_{\text{1}}\right)p\left(r_{\text{1}},\Delta t_{\text{1}}|r_{\text{0}}\right)\text{dr,}\end{array}$$E77

where p(*,t|r₀)=1-Q(t|r₀).

6.2.2. Sampling algorithm

The objective is to generate random variates from f(r)≡4πr²p(r,t|r₀) (ref. [23]). It can be noted that f(r) is composed by the sum of two Gaussians and a negative term, since k_a≥0. Thus, f(r)≤h(r), where

$$h(r)=\frac{r}{r_{\text{0}}\sqrt{\text{4πDt}}}\left\{\text{exp}\left[-\frac{{(r-r_0)}^{\text{2}}}{\text{4Dt}}\right]+\text{exp}\left[-\frac{{(r+r_{\text{0}}-\text{2σ})}^{\text{2}}}{\text{4Dt}}\right]\right\}.$$E78

Since f(r) is bound by the sum of two Gaussian functions, it can be sampled by a rejection method. To do so, h(r) is written as the sum of four terms, which, after rearrangement and truncation to the positive ranges, are:

$$\begin{array}{l}{h}_{\text{1}}(r)=\frac{{\text{r-r}}_{\text{0}}}{r_{\text{0}}\sqrt{\text{4πDt}}}\text{exp}\left[-\frac{{({\text{r-r}}_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right]\times {\text{1}}_{\left[r\ge r_{\text{0}}\right]},\\ h_{\text{2}}(r)=\frac{r+r_{\text{0}}-\text{2σ}}{r_{\text{0}}\sqrt{\text{4πDt}}}\text{exp}\left[-\frac{{(r+r_{\text{0}}-\text{2σ})}^{\text{2}}}{\text{4Dt}}\right]\times {\text{1}}_{\left[r\ge \sigma \right]},\\ h_{\text{3}}(r)=\frac{r_{\text{0}}}{r_{\text{0}}\sqrt{\text{4πDt}}}\text{exp}\left[-\frac{{({\text{r-r}}_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right]\times {\text{1}}_{\left[r\ge \sigma \right]},\\ h_{\text{4}}(r)=\frac{{\text{2σ-r}}_{\text{0}}}{r_{\text{0}}\sqrt{\text{4πDt}}}\text{exp}\left[-\frac{{(r+r_{\text{0}}-\text{2σ})}^{\text{2}}}{\text{4Dt}}\right]\times {\text{1}}_{\left[r\ge \sigma \right]}\times {\text{1}}_{\left[r_{\text{0}}\le \text{2σ}\right]}\end{array}$$E79

where 1_[condition] takes the value 1 when the condition is true, and 0 when it is false. A related function is also defined:

$$h_{\text{4}}^{\text{'}}(r)=\frac{r_{\text{0}}}{r_{\text{0}}\sqrt{\text{4πDt}}}\text{exp}\left[-\frac{{({\text{r-r}}_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right]\times {\text{1}}_{\left[r\le \sigma \right]}.$$E80

Thus, $$f(r)\le h(r)\stackrel{def}{=}{\displaystyle \sum h_i(r)}$$. It is possible to generate a random variate with density proportional to h(r), since h(r) is a mixture of known probability distribution functions. To do so, the weight of the contributions of h₁, h₂ and h₃+h₄ are needed:

$$\begin{array}{l}{p}_{\text{1}}={\displaystyle \int h_{\text{1}}(y)\text{dy}=\sqrt{\frac{\text{Dt}}{\pi r_{\text{0}}^{\text{2}}}}},\\ p_{\text{2}}={\displaystyle \int h_{\text{2}}(y)\text{dy}=\sqrt{\frac{\text{Dt}}{\pi r_{\text{0}}^{\text{2}}}}}\text{exp}\left[-\frac{{(\sigma {\text{-r}}_{\text{0}})}^{\text{2}}}{\text{4Dt}}\right],\\ p^{\text{'}}={\displaystyle \int (h_{\text{3}}(y)+h_{\text{4}}(y))\text{dy}=}\Phi \left(\frac{r_{\text{0}}-\sigma }{\sqrt{\text{2Dt}}}\right)+{\text{1}}_{\left[r_{\text{0}}\le \text{2σ}\right]}\frac{\text{2σ}-r_{\text{0}}}{r_{\text{0}}}\Phi \left(\frac{\sigma {\text{-r}}_{\text{0}}}{\sqrt{\text{2Dt}}}\right).\end{array}$$E81

where Φ(x) is the normal distribution function:

$$\Phi (x\sqrt{\text{2}})=\text{1}-\text{Erfc}(x)/\text{2}.$$E82

From this algorithm 8a is obtained:

Algorithm 8a

CALCULATE p1 , p2 and p'. SET s=p1+p2+p‘

GENERATE a uniform [0,1] random variateU

IF sU∈ [0,p1], SET$$X\leftarrow r_{\text{0}}+\sqrt{\text{4DtE}}$$, where E is standard exponential.

IF sU∈ (p1,p2], SET$$X\leftarrow \text{2σ}-r_{\text{0}}+\sqrt{{(r_{\text{0}}-\sigma )}^{\text{2}}+\text{4DtE}}$$, where E is standard exponential.

IF sU∈ [p1+p2,s], use algorithm 8b to generate X

RETURN X

The second part is given by Algorithm 8b:

Algorithm 8b

REPEAT

GENERATE N as standard normal, U uniform on [0,1]

Set $$X\leftarrow r_{\text{0}}+\sqrt{\text{2Dt}}N$$

UNTIL X>σ, or jointly X<σ, Dt<2σ^2 and U≤(2σ-r0)/r0

In the former case, RETURN X

In the latter case, RETURN 2σ-X

Algorithm 8 generates numbers with density probability proportional to h(r). However, random numbers distributed as f(r) are needed. Since f(r) is bound by h(r), a rejection method is used as the final step, i.e. the overall algorithm generate pairs of random numbers (X,U) with X having density distribution proportional to h and U uniform on [0,1] until Uh(X)≤f(X) and returns X. Results obtained by using this algorithm are shown in Figure 11.

Figure 11.

a) Green's function 4πr²p(r,t|r₀) for D=1, R=1, r₀=1.5 and α=-2 at t=1, 2, 4, 8, 16 and 32. Analytical functions: (―); Result of sampling: (■). b) Survival probability Q(t|r₀) as function of time for D=1, R=1, r₀=1.5 and k_a=0.1 x 4πRD (―), 4πRD (---) and 10 x 4πRD ( ).

The time discretization property of the Green’s functions (Equations 77) are also used to validate the simulation results and the sampling algorithm.

This algorithm is used to sample the length of the separation vector r. To fully caracterize the vector its direction should be specified by the angles θ and ϕ in a spherical coordinate system. An exact answer including the angles could be obtained by sampling the multi-variate Green's function given by Equation 67. An approximate method to obtain the deflection angles, which can be used here, is given by [24].

After sampling of the vectors r and R, the new positions of the particles r_A and r_B are obtained by inverting the transformation given by Equation (61):

$$\begin{array}{l}{\text{r}}_{\text{A}}=\frac{\sqrt{D_A/D_B}\text{R}-D_A\text{r}}{D_A+D_B},\\ {\text{r}}_{\text{B}}=\frac{\sqrt{D_A/D_B}\text{R}+D_B\text{r}}{D_A+D_B}.\end{array}$$E83

6.2.3. Application to radiation chemistry

Ionizing radiation interacts with water molecules and creates radiolytic species, which react with cellular components and eventually lead to biological effects [14]. The precise location of the created radiolytic species are highly dependent of the type and energy of ionizing radiation. Therefore, considerable effort has been devoted to the Monte-Carlo simulation of ionizing radiation tracks, which provides such information [25]. After their formation, the radiolytic species react with each other and with other cellular molecules. Hence, radiation chemistry is a very important link between the radiation track structure and cellular damage and biological effects [26].

The simulation of radiation chemistry is difficult for several reasons. One important problem is the calculation time. Since the number of interactions between N particles is N(N-1)/2, the calculation time increases as ~N². Therefore, simulation of systems comprising over 10000 particles is difficult, even in modern computers. For example, Uehara and Nikjoo [27] reported that the simulation of time-dependent yields of chemical species produced by five 1-MeV electron tracks comprising about 1,500 particles from 10^-12 to 10^-6 s requires about 250 h of calculation time. One promising approach is the new hardware called general purpose graphic processing units (GPGPU), which are special cards comprising a large number of processors designed to perform a large number of calculations simulaneously. Regarding this, the sampling algorithms proposed in this book chapter do not require much memory and could be implemented on a GPGPU, which could allow the simulation of large particle systems.

We conclude this section by showing a simulation result of water radiolysis by a ¹²C⁶⁺ ion of 25 MeV/amu [28] in Figure 12. Briefly, the track structure and the position of all the radiolytic species (^.OH, H^., e^-_aq, H₂, H₂O₂,..) at ~10^-13 s are calculated by the code RITRACKS, described in [25] and references therein. The time evolution of the radiolytic species are calculated by using a step-by-step Monte-Carlo simulation code of radiation chemistry [28-30]. Each dot in Figure 12 represents a radiolytic species. The changes in the colors of the dots indicates that chemical reaction occurs. The number of chemical reactions (the radiolytic yields) which are calculated by the program have been validated with experimental data [30].

Figure 12.

Time evolution of the radiolytic species created by a ¹²C⁶⁺ ion track, 25 MeV/amu (LET~75 keV/μm) at 10^-13, 10^-9, 10^-7 and 10^-6 s.

7. Conclusion

The simulations based on the Green's function of the DE are able to accurately describe the evolution of particles in time and space in various geometries, even if the number of particles is small. In fact, this approach is now used in several fields. Of particular interest in radiobiology is the study of radiation chemistry, biochemical networks and bimolecular reactions.

The Green's functions and their respective sampling algorithms have been presented for simple systems in one, two and three dimensions. Obviously, it is difficult to develop sampling algorithm even for simple systems because of the complex mathematical form of the Green’s functions. However, it is very useful to develop sampling algorithms for systems for which analytical solutions are known, because they can be used to validate future models which will be based on numerical solutions of the DE. The DE can be solved easily by numerical methods such as finite-difference schemes, but these techniques may have restrictions for computational stability and the boundary conditions are often difficult to implement [31]. The systems with complex boundary conditions are particularly relevant to study biological networks. However, as seen in the applications, since the Green's functions of systems with complex boundary conditions are often subdensities of the Green's functions for the corresponding system with reflective boundary conditions, the rejection method can be used.

The algorithms developed in this chapter are already used to simulate the radiation chemistyoccuring in liquid water [28-30]. In the future, they will allow the validation and verification of future calculations, which will use the numerical solution of the diffusion equation. Eventually, the approach depicted in this book chapter will be used to simulate the time evolution of radiolytic species created by ionizing radiation and their reactions with molecules in biological media [23].

Acknowledgments

This work was support by the NASA Space Radiation Risk Assessment Program. The help of Dr. Luc Devroye (McGill University) with sampling algorithms is greatly acknowledged.