From molecular dynamics to Brownian dynamics

Three coarse-grained molecular dynamics (MD) models are investigated with the aim of developing and analysing multi-scale methods which use MD simulations in parts of the computational domain and (less detailed) Brownian dynamics (BD) simulations in the remainder of the domain. The first MD model is formulated in one spatial dimension. It is based on elastic collisions of heavy molecules (e.g. proteins) with light point particles (e.g. water molecules). Two three-dimensional MD models are then investigated. The obtained results are applied to a simplified model of protein binding to receptors on the cellular membrane. It is shown that modern BD simulators of intracellular processes can be used in the bulk and accurately coupled with a (more detailed) MD model of protein binding which is used close to the membrane.


Introduction
Brownian dynamics (BD) simulations have been used for the modelling of a number of spatio-temporal processes in cellular and molecular biology in recent years, including models of intracellular calcium dynamics [1], the MAPK pathway [2] and signal trasduction in Escherichia coli chemotaxis [3]. In these applications, trajectories and interactions between key biomolecules (e.g. proteins) are calculated using BD methods, while other components of the system (e.g. solvent molecules), which are of no special interest to a modeller, are not explicitly included in the simulation, but contribute to the dynamics of Brownian particles collectively as a random force. This reduces the dimensionality of the problem, making BD less computationally intensive than the corresponding molecular dynamics (MD) simulations.  Denoting the position of a Brownian particle by X = [X 1 , X 2 , X 3 ] and its diffusion constant by D, a simple model of Brownian motion is given by the (overdamped) Langevin equation where W i , i = 1, 2, 3, are three independent Wiener processes [4]. BD approaches which are based on (1.1) have been implemented in a number of software packages designed for spatiotemporal modelling in systems biology, including Smoldyn [5], MCell [6], Green's function reaction dynamics [7] and the first-passage kinetic Monte Carlo method [8]. The software package Smoldyn discretizes (1.1) using a fixed time step t, i.e. it computes the time evolution of the position X ≡ X(t) of each molecule by where [ξ 1 , ξ 2 , ξ 3 ] is a vector of normally distributed random numbers with zero mean and unit variance. A different BD approach is implemented in the Green's function reaction dynamics [2] which evolves time using a variable time step. It approximately computes the time when the next reactive event happens. This means that trajectories of molecules which are not surrounded by other reactants can be simulated over longer time steps.
Although the BD models are becoming a popular choice for stochastic modelling of intracellular spatio-temporal processes, several difficulties prevent the use of BD for some systems. First of all, detailed BD models are often more computationally intensive than coarser spatio-temporal models which are written for concentrations of biochemical species. In some applications (e.g. intracellular calcium dynamics [1] or actin dynamics in filopodia [9]), individual trajectories (computed by BD) are important only in certain parts of the computational domain, while in the remainder of the domain a coarser, less detailed, method can be used. In these applications, the computational intensity of BD simulations can be decreased by using multiscale methods which efficiently and accurately combine models with a different level of detail in different parts of the computational domain [10,11].
Another difficulty of BD simulations in cell and molecular biology is that detailed BD models require more parameters than coarser (macroscopic) models. In some studies, macroscopic parameters are used to infer BD parameters [5,12]. For example, knowing the macroscopic reaction rate k of a bimolecular reaction A + B → C and diffusion constants of reactants, one can calculate a (microscopic) reaction radius of BD simulations which gives the corresponding macroscopic parameters in the limit of many particles. In the classical Smoluchowski limit [13], a bimolecular reaction occurs whenever the distance of reactants is less than the reaction radius where D A (resp. D B ) is the diffusion constant of reactant A (resp. B). Although this approach is commonly applied in stochastic reaction-diffusion models, it is not the most satisfactory, because different microscopic models can lead to the same macroscopic process and parameters [12,14]. For example, the simplest Smoluchowski model (1.3) assumes that all collisions are reactive but, in reality, many non-reactive collisions of molecules happen before a reactive collision occurs. Therefore, some algorithms postulate that molecules only react with a certain rate (probability) when the distance between reactants is less than a modified reaction radius (which is larger than ). Other methods discretize the Langevin equation with time step t and substitute the Smoluchowski formula (1.3) (which is valid for an infinitely small time step) by a tabulated function computed numerically [5]. However, all of these approaches are verified by considering the macroscopic limit (of many reactants) and showing that the reaction occurs with the given rate k in this limit. A different approach to parametrize BD models is to use a more detailed description written in terms of MD. In this paper, we investigate connections between BD and MD models with the aim of developing and analysing multi-scale methods which couple BD and MD simulations. We consider a (computationally intensive) MD simulation in a domain Ω which is either one or three dimensional, i.e. Ω ⊂ R or Ω ⊂ R 3 . Our main goal is to design and analyse multiscale methods which can compute spatio-temporal statistics with an MD level of detail in the subdomain Ω D ⊂ Ω. We define (1.4) where Ω D and Ω are open sets, the (open) set Ω C is the complement of Ω D and I is the shared interface (boundary) between Ω D and Ω C . In the multi-scale set-up (1.4), we use a detailed MD model in Ω D and a coarser BD model in Ω C .
In this paper, we focus on a simple MD approach which is introduced in §2 and §3. A few (heavy) particles with mass M and radius R are coupled with a large number of light point particles with masses m M. The collisions of particles are without friction, which means that post-collision velocities can be computed using the conservation of momentum and energy [15,16]. We will introduce and study three MD models which make use of elastic collisions. They will be denoted as MD models [ 1). One can also show that these models converge to the Langevin description [15][16][17]

One-dimensional molecular dynamics model [A]
The MD model [A] is described in terms of positions x i and velocities v i , i = 1, 2, 3, . . ., of heat bath particles, and positions X i and velocities V i , i = 1, 2, . . . , N, of heavy particles of mass M m, where m is the mass of a heat bath particle [15]. In our computer implementations, we will consider a finite number of heat bath particles. However, we formulate the MD model in terms of (countably) infinitely many heat bath particles which are initially distributed along the real line according to the Poisson distribution with density where μ is given by (1.5), and D and γ are positive constants. This means that the probability that there are j particles in a subinterval [a, b] ⊂ R, a < b, is equal to . Initial velocities of heat bath particles are given by the normal distribution Let us consider a model with a single heavy particle, i.e. N = 1. Then its location X 1 and velocity V 1 will be denoted as X and V to simplify our notation. Whenever the heavy particle collides with the light particle with velocity v i , their velocities are updated using the conservation of mass and momentum byṼ where tildes denote post-collision velocities. Using (1.5), the equations (2.3) and (2.4) can be rewritten asṼ The following result can be shown for the above MD model [A].
Lemma 2.1. Let γ > 0 and D > 0. Let us consider the heavy particle of mass M with initial position X μ (0) = X 0 and initial velocity V μ (0) = V 0 which is subject to elastic collisions (2.5) with heat bath particles of mass m whose initial positions and velocities are distributed according to (2.1) and (2.2). Then X μ and V μ converges (as μ → ∞) in distribution to the solution X and V of equations where X(0) = X 0 and V(0) = V 0 . That is, X and V solve the one-dimensional version of equations (1.6) and (1.7).
Proof. This lemma can be proved using the main theorem in Holley [15], where it is shown that a similar process converges to the Ornstein-Uhlenbeck process (1.7) for velocities. Although our function f μ does not satisfy all assumptions of the main theorem of Holley [15], a simple rescaling of our parameters leads to a process which is covered by Holley's theorem. In §4 of this paper, we also rederive this result as one of the consequences of multi-scale analysis (see (4.11)).
As the goal of this paper is to study the behaviour of computational algorithms, we formulate the MD model . Let b ∈ R and t > 0. Let us assume that heat bath particles are distributed according to the Poisson distribution with density (2.1) in the interval (−∞, b). Their initial velocities are given according to (2.2) and there are no particles in the interval (b, ∞) at time t = 0. Then the average number of particles in the interval (b, ∞) at time t = t is The positions x and velocities v of these particles are distributed according to where H(·) is the Heaviside step function. In particular, the positions of the particles at point x ∈ (b, ∞) are distributed at time t = t according to where erfc(z) = 2/ √ π ∞ z exp(−s 2 ) ds is the complementary error function.
Proof. Particles which are at point x ∈ (b, ∞) at time t = t were previously at point x − v t at time t = 0. In particular, there will be non-zero heat bath particles with velocity v at point x at time t = t provided that x − v t < b, which implies (2.8). Consequently, the density of particles which are at point Thus, we proved (2.9). Integrating this formula over x in interval (b, ∞), we obtain the average number of particles which are in the interval (b, ∞) at time Substituting (2.1) for λ μ and (2.2) for σ μ , we obtain (2.7).
We use lemmas 2.1 and 2.2 to design a computational test for multi-scale methods. As the number n ≡ n(t) of heat bath particles in [−L, L] is much larger than the number N of large particles, we will focus on models of a single large particle, i.e. N = 1, which is described by its position X and velocity V. We choose a small time step t. One iteration of the MD algorithm is presented in table 1. We first compute the positions of all particles at time t + t in step [A1] by assuming that particles do not interact. Then we use (2.5) to incorporate collisions in step [A2]. As all heat bath particles have the same mass, the collisions between them result in exchange of colliding particles' positions and velocities. In particular, step [A2] can be implemented by sorting the heat bath particles during every iteration. If the large particle collides with a heat bath particle, then their positions at time t + t will be obtained in . We assume that t is chosen so small that (2.7) is much smaller than 1. Then (2.7) can be interpreted as a probability of introducing one particle from the left (resp. right) during one time step. Using lemma 2.2, the new particle will be introduced at the (left boundary) position which is sampled according to the probability distribution proportional to (x; t, −L) in step [A4]. To sample from this probability distribution, we scale and shift a random number sampled from the complementary error function distribution

[A1]
Compute 'free-flight positions' of heat bath particles and the large particle at time t + t bŷ [A2] Compute post-collision velocities by (2.5) for every pair of particles which collided. Compute their postcollision positions x i (t + t) and X(t + t) by updating their 'free-flight positions'x i (t + t) and X(t + t).
[A3] Terminate trajectories of heat bath particles which left the domain [−L, L]. Update n accordingly.
[A4] Generate a random number r 1 uniformly distributed in (0, 1). If r 1 < γ (μ + 1) t/8, then increase n by 1, and introduce a new heat bath particle at a position sampled according to the probability distribution proportional to (x; t, −L). Its velocity is sampled according to the probability distribution proportional to [A5] Generate a random number r 2 uniformly distributed in (0, 1). If r 2 < γ (μ + 1) t/8, then increase n by 1, and introduce a new heat bath particle at position x n (t + t) with velocity v n (t + t) which are sampled according to probability distributions (2.11) and (2.12).
[A6] Continue with step [A1] using time t = t + t.  Table 2. The acceptance-rejection algorithm which is used to sample random numbers which are distributed according to the probability distribution π erfc(z), where z ∈ (0, ∞). In our simulations, we use constants a 1 and a 2 given by (2.10). -Generate a random number ζ 1 uniformly distributed in (0,1).
-If ζ 1 ζ 3 < a 2 erfc(ζ 2 ), then choose ζ 2 as a sample from the probability distribution π erfc(z). Otherwise, repeat the algorithm. π erfc(z), where z ∈ (0, ∞). An acceptance-rejection algorithm for sampling random numbers from π erfc(z) is given in table 2. We use it with the constants a 1 and a 2 given by a 1 = 0.532 and a 2 = 0.814. (2.10) The values of constants a 1 and a 2 were numerically estimated to maximize the total acceptance probability a 2 /(a 1 √ π ) of the acceptance-rejection algorithm in table 2. Using (2.10), we accept 86% of proposed numbers ζ 2 . A particle introduced close to the right boundary in step [A5] will have its position sampled according to the probability distribution where C 1 is a normalization constant. The probability distribution (2.11) is proportional to (−x; t, −L) and can be justified using the same argument as lemma 2.2. To sample from the probability distribution (2.11), we again use the acceptance-rejection algorithm in table 2 with parameters a 1 and a 2 given by (2.10). In step [A5], we also sample the velocity v ∈ R of the new particle using the truncated Gaussian distribution   figure 1b as the solid line. To illustrate the limiting result in lemma 2.1, we also plot the mean squared displacement corresponding to the limiting solution X of (2.6). It can be analytically computed as where E[·] denotes the expected value. It is plotted as the dashed line in figure 1b. We also plot the mean squared displacement corresponding to the overdamped limit (1.1), i.e. √ 2Dt, as the dot-dashed line in figure 1b. If we neglect the exponential terms in (2.13), we obtain This approximation is plotted in figure 1b as the dotted line. We will use (2.14) later in §6 to couple the overdamped BD model (1.1) with MD simulations.

Three-dimensional molecular dynamics models [B] and [C]
. . , N, of heavy particles of mass M m, where m is the mass of a heat bath particle [16]. We again define μ by (1.5). We will denote by R the radius of a heavy particle.
where v j is the velocity of the heat bath molecule which collided with the ith heavy molecule, tildes denote post-collision velocities, superscripts ⊥ denote projections of velocities on the line through the centre of the molecule and the collision point on its surface, and superscripts denote tangential components.
Let the velocities of heat bath particles be distributed according to  [16]. Their theorem expresses the limiting equation of a process with given λ μ and f μ (v) in terms of moments of f μ . These moments can be analytically evaluated to verify the statement of lemma 3.1. We will also rederive this result in §5 as a consequence of the analysis of multi-scale methods.
Lemma 3.1 can be viewed as a different formulation of theorem 2.1 in [16]. They were interested in the limit m → 0, which is equivalent to μ → ∞. Considering the scaling m 3/2 f (vm 1/2 ) of the velocity distribution of heat bath particles (with density scaled as λ/m 1/2 ), they derived formulae for γ and D in terms of moments of f and λ. To formulate lemma 3.1, we inverted their results by deriving the appropriate distributions (3.3) and (3.4) which lead to the limiting BD model with a given D and γ .

(b) Molecular dynamics model [C]
In lemma 3.1, we used the normal distribution for velocities (3.4). Another option is to use heat bath particles with fixed speed as is done in the following lemma 3.2. We denote the resulting MD model as the MD model [C]. Let the velocities of heat bath particles be distributed according to and δ is the Dirac delta function. Let us consider one heavy molecule in this heat bath at position X μ with velocity V μ . Then X μ and V μ converge (in the sense of distribution) to the solution of (1.6) and (1.7) in the limit μ → ∞.
Lemma 3.2 can again be proved using theorem 2.1 in [16] that is applicable to any spherically symmetric velocity distribution which has at least four finite moments.

(c) Boundary conditions for molecular dynamics models [B] and [C]
Next, we generalize lemma 2.2 to the three-dimensional case. This will help us to specify boundary conditions for simulations which use the MD models [B] and [C] in finite domains.

7)
and the average number of particles in the semi-infinite cuboid (b, ∞) × (0, 1) 2 at time t = t is Proof. Formula (3.7) is a generalization of formula (2.8) in lemma 2.2 and can be justified using the same arguments. To prove (3.8), we will distinguish two cases.
First, let us consider that λ μ and f μ (v) are given by ( which is formula (3.8).

From one-dimensional molecular dynamics model [A] to Brownian dynamics
In the case of one-dimensional MD model [A], the situation is schematically shown in figure 2a-c, where we consider only one large (heavy) particle, i.e. N = 1. The large particle can either be in The large particle covers the interval (X(t) − R, X(t) + R). Let us consider that the large particle intersects the interface I, as is shown in figure 2b. Then I ⊂ (X(t) − R, X(t) + R), which is equivalent to X(t) ∈ (−R, R). The heat bath particles are simulated in Ω D using the MD model [A]. Let us choose t so small that the probability of two collisions happening in the time interval (t, t + t) is negligible. As we do not explicitly simulate heat bath particles in Ω C , we will consider an additional correction of the velocity of the heavy particle in the form whereṼ(t + t) is the post-collision velocity of the heavy particle at time t + t, which only takes into account collisions with the heat bath particles from the left. It is either equal to V(t) or computed by (2.5) if a collision with a heat bath particle occurred in Ω D . Equation (4.1) is adding both the drift term α(V(t)) t and the noise term β(V(t)) √ t ξ , where ξ is a normally distributed random number with zero mean and unit variance. The drift and noise terms implicitly take into account collisions at the right boundary (X(t) + R) of the heavy particle. Passing t → 0, we observe that the contributions of the collisions at the right boundary are given by the Itō stochastic differential equation If we explicitly modelled heat bath particles in Ω C , then they would be distributed according to the Poisson distribution with density λ μ in the interval (X(t) + R, ∞). Their initial velocities would be given according to (2.2). Thus, using lemma 2.2 and (2.5), we can estimate the drift coefficent of the stochastic differential equation (4.2) to get where H(·) is the Heaviside step function. Using (2.2), we obtain where λ μ and σ μ are given by (2.1) and (2.2). In the limit μ → ∞, we have V/ √ μ + 1 → 0. Thus, we use the Taylor expansion in (4.3) to get The noise term in (4.2) can be computed by Using (2.2), we obtain Using (2.1), (2.2) and the Taylor expansion, we obtain Step [M4] is the same as [A4], which introduces heat bath particles that have entered Ω D through its left boundary x = −L during the time interval (t, t + t). The boundary at x = 0 is treated in step [M5] if the heavy particle does not intersect with this boundary. We assume that t is chosen so small that (2.7) is much smaller than 1. Then (2.7) can be interpreted as a probability of introducing one particle from the left (resp. right) during one time step. A particle introduced close to the right boundary of Ω D in step [M5] will have its position sampled according to the probability distribution  Terminate trajectories of heat bath particles which left the subdomain Ω D = (−L, 0). Update n accordingly.
[M4] Implement the influx of heat bath particles through the boundary x = −L using step [A4].
If X(t) ∈ (−R, R), then generate a random number r 2 uniformly distributed in (0, 1). If r 2 < γ (μ + 1) t/8, then increase n by 1, and introduce a new heat bath particle at position x n (t + t) with velocity v n (t + t) which are sampled according to probability distributions (4.6) and (4.7).
If X(t) ∈ [R, L), then update the velocity of the heavy particle using (4.8).
[M8] Continue with step [M1] using time t = t + t. distribution (4.6), we again use the acceptance-rejection algorithm in table 2 with parameters a 1 and a 2 given by (2.10). In step [M5], we also sample the velocity v ∈ R of the new particle using the truncated Gaussian distribution where C 2 is a normalization constant. To sample random numbers according to the truncated normal distributions in steps [M4] and [M5], we again use the acceptance-rejection algorithm which is derived as proposition 2.3 in [18]. If the heavy particle does intersect with the boundary I, then step [M6] is executed. It uses (4.1) to incorporate collisions of heat bath particles from the right. If the particle does not intersect with Ω D , then we simulate it in step [M7] using the discretized version of (1.7) given by , is plotted in figure 3a. It is compared with the distribution obtained by the limiting BD model (2.6), which is, for t γ −1 , given by [19] 1 4π , where t * = 3 2γ . The stochastic differential equation (4.2) was derived for collisions from the right. Using the same argument, we can also derive a stochastic differential equation which is approximating the effect of collisions from the left. We obtain (4.10) Adding (4.2) and (4.10) and using the independence of noise terms in (4.2) and (4.10), we can approximate collisions from both sides by the following stochastic differential equation for the velocity of the heavy particle:

From three-dimensional molecular dynamics models [B] and [C] to Brownian dynamics
We use a simple multi-scale geometry where domain Ω = R 3 is divided into two half spaces. Boundary conditions for heat bath particles at the interface I = {0} × R 2 can be specified using lemma 3.3. As in §4, we need to analyse the behaviour of a heavy molecule when it intersects with the interface I. Such a molecule is subject to the collisions with heat bath particles on the part of its surface which lies in Ω D . This has to be compensated by using a suitable random force from Ω C , so that the overall model is equivalent to (1.6) and (1.7) in the BD limit. To simplify the presentation of the algorithm, we use the same time step in Ω D and Ω C . In §6, we present coupling of threedimensional MD models with the BD model (1.1), which will make use of different time steps in different parts of the computational domain.
The heavy particle is the ball with centre X = [X 1 , X 2 , X 3 ] with velocity V = [V 1 , V 2 , V 3 ] and radius R. It intersects the interface I if X 1 (t) ∈ (−R, R). Let us consider that heat bath particles are simulated in Ω D using the MD model [B] or the MD model [C]. Let us choose t so small that the probability of two collisions happening in the time interval (t, t + t) is negligible. As we do not explicitly simulate the heat bath particles in Ω C , we will consider an additional correction of the velocity of the heavy particle in the form

V(t + t) =Ṽ(t + t) + α(X(t), V(t)) t + β(X(t), V(t))
√ t ξ , (5.1) whereṼ(t + t) is the post-collision velocity of the heavy particle at time t + t, which only takes into account collisions with the heat bath particles from Ω D . It is either equal to V(t) or computed by (3.1) if a collision with a heat bath particle occurred in Ω D . Note that we dropped the subscript μ in (3.1) and (3.2) to simplify our notation. Equation (5.1) is a generalization of (4.1) to threedimensional simulations where α(X(t), V(t)) t is the drift vector, β(X(t), V(t)) √ t ξ is the noise term and ξ = [ξ 1 , ξ 2 , ξ 3 ] is the vector of three normally distributed random numbers with zero mean and unit variances. Passing t → 0, we observe that the contributions of the collisions from Ω D are given by the Itō stochastic differential equation

Then the average change of the j-th component of the velocity of the heavy molecule caused by collisions with heat bath particles in the time interval (t, t + t) at the surface area
Proof. Let us consider that a heat bath particle which was at point x at time t collided with the heavy molecule at time t + τ ∈ (t, t + t) at the surface point which had coordinate y at time t. Then the coordinate of the surface point at the collision time t + τ was y + τ V(t) and the precollision velocity of the heat bath molecule was v = V(t) + (y − x)/τ . Using equation (3.1), we can write the change of the velocity of the heavy molecule during the collision as The position x of the heat bath particle must be in the half space which lies above the plane tangent to the heavy molecule at the collision point y + τ V(t). It can be parametrized by where c 1 > 0, c 2 ∈ R, c 3 ∈ R, and (y − X(t))/R, η 2 , η 3 is the orthornormal basis in R 3 . Then (5.5) reads as follows: Thus, we have Integrating over c 1 , we deduce (5.4). Using lemma 5.1, we can compute the drift coefficient α(X(t), V(t)) in equation (5.2) as follows: where S(X(t)) is the part of the surface of the heavy molecule which intersects the BD subdomain Ω C , i.e.
Substituting (5.4) into (5.7), we have Using (3.3) and (3.4) and evaluating the surface integrals, we obtain where we dropped the dependence on time t to shorten the resulting formulae. The noise matrix β(X(t), V(t)) will be estimated using β(X(t), 0), i.e. we will only use the first term in the Taylor expansion in V. Using similar arguments as in the proof of (5.6) and (5.7), we have where the last equation can be verified using the same argument as equation (5.10). Note that, by substituting X 1 = R into (5.8), (5.9) and (5.11), we verify the limiting result in lemma 3.1.

(b) Molecular dynamics model [C]
Equations ( . We obtain and α 2 (X, V) and α 3 (X, V) are again given by (5.9). Substituting (3.5) and (3.6) in (5.10) and integrating, we obtain that noise matrix β(X, V) satisfies (5.11). We again note that the special choice X 1 = R in (5.12) can be used to verify the limiting result in lemma 3.2. We consider a three-dimensional generalization of the illustrative problem from figure 3. One heavy particle which starts at position X(0) = [0, 0, 0] with velocity V(0) = [0, 0, 0] is simulated using a three-dimensional generalization of the algorithm [M1]-[M8]. We use the MD model [C] in Ω D = (−∞, 0) × R 2 and the BD model (1.6) and (1.7) in Ω C = (0, ∞) × R 2 . In step [M6], we replace (4.1) with its three-dimensional analogue (5.1), where drift α and diffusion coefficient β are given by (5.12), (5.9) and (5.11). The distribution of X 1 positions of the heavy particle at time t = 1, computed using 10 5 realizations of the multi-scale algorithm, is plotted in figure 4a. The limiting BD result is again given by (4.9). In figure 4b, we plot the time evolution of the mean squared displacement. The mean squared displacement corresponding to the limiting BD model (1.6) and (1.7) is given by (2.13) multiplied by √ 3 because we have three spatial dimensions (dashed line). As expected, the models compare well. The mean squared displacement obtained for the BD model (1.1) is plotted as the dot-dashed line.

Application to protein binding to receptors
In this section, we apply our results to a simplified model of protein binding to receptors on the cell membrane. We consider simple geometry, which is schematically shown in where h > 0 and the interface I is at x 1 = h. Diffusing proteins are modelled as spheres of radius R. We consider that a protein which hits the boundary ∂Ω M will bind to a receptor with probability P, and otherwise it is reflected. This type of a reactive boundary condition is common for BD simulations [20]. In the case of MD, more detailed models of protein binding could be introduced in Ω D [21,22]. However, the main goal of this section is to show how an MD model in Ω D can be coupled with BD simulators which have been developed for simulations of intracellular processes. Therefore, we keep the MD model in Ω D as simple as possible.
If we used BD model (1.6)-(1.7) in Ω C , then the situation would be more or less the same as in §5. However, modern BD simulators of intracellular processes work with the high-friction limit (1.1) rather than (1.6)-(1.7). For example, the software package Smoldyn discretizes (1.1) with a fixed time step and uses (1.2) to update positions of diffusing proteins. In particular, it uses larger values of time step than we used in §5. Then the problem can be formulated as follows: we would like to use the MD model with time step t in Ω D and couple it with the BD model (1.2) with larger time step t, namely if the diffusing molecule is far away from Ω D . We couple these models using the intermediate BD model (1.6)-(1.7). We introduce two additional interfaces where h < h 2 < h 3 < L 1 , as shown in figure 5. We denote i.e. Ω C1 and Ω C2 are two overlapping subdomains of Ω C . We simulate the time evolution of the position X(t) of one protein molecule. If X(t) ∈ Ω C2 , then the BD model (6.1) will be used in Ω C2 until the molecule leaves Ω C2 . Then we switch to the shorter time step t and use the BD model (1.6)-(1.7) in Ω C1 . The protein molecule can leave Ω C1 in two possible ways: (i) The protein molecule crosses the interface I 3 .
Then we revert to the BD model (6.1) which is used in Ω C2 . (ii) The protein molecule crosses the interface I.
Then we use the method from §5 for coupling the MD model in Ω D with the BD model (1.6)-(1.7) in Ω C1 .
As the subdomains Ω C1 and Ω C2 overlap, we can use the limiting result (4.9), which implies that, for times t ≥ γ −1 , the BD model (1.6)-(1.7) is given by the BD model (6.1) shifted by time t * = 3/(2γ ). In particular, we will also add or subtract t * from the time variable whenever we switch between BD models.  The solution of (6.2) and (6.3) with uniform initial condition (x 1 , 0) ≡ const. is given by the black solid line in figure 6a. As we only visualize the distribution along the X 1 -axis in figure 6a, we can further decrease the computational cost by truncating the simulation domain in the x 2 and x 3 directions to the region close to the protein molecule. That is, we only simulate small particles in the subdomain [0, h] × [X 2 (t ) − h 4 , X 2 (t ) + h 4 ] × [X 3 (t ) − h 4 , X 3 (t ) + h 4 ], where t is the time when the protein molecule enters Ω D ∪ Ω C1 . This subdomain (moving window) is shifted accordingly whenever X 2 (t) or X 3 (t) approach its boundary. The probability that the protein is adsorbed to the surface is given as a function of time in figure 6b. It again compares well with the results obtained by the limiting PDE system (6.2) and (6.3).

Discussion
In this paper, a multi-scale approach that uses MD simulations in a part of the computational domain and BD simulations in the rest of the domain has been presented and analysed. The ultimate goal of this research is to use MD to help parametrize BD models of intracellular processes. One application area is modelling proteins in an aquatic environment, which is useful for understanding protein binding to receptors (surfaces) as shown in §6.
The main idea of the presented coupling of MD and BD models is based on using equations (4.1) and (5.1) and estimating drift and diffusion coefficients for velocities of molecules which cross the interface I. This coupling uses the same time step for the BD model (1. for the MD model. In §6, it was shown that this is not a limiting step of this approach, because the BD model (1.6)-(1.7) is only needed in a small part of the domain next to Ω D . Then the coarser BD model (6.1) with larger time step can be used in the rest of the simulation domain, using a suitable overlap region. Another overlap region could be used to couple BD simulations with mean-field PDE-based models [11]. Then multi-scale models which couple BD (of point particles) with coarser reaction-diffusion approaches would be capable of further increasing time scales and space scales of simulations [10,11].
MD models considered in this paper are relatively simple and analytically tractable, describing water molecules as point particles. An important generalization is to consider more complicated MD models of water molecules [23]. For example, Rahman & Stillinger [24] model water molecules as rigid asymmetric rotors. That is, each water molecule is described by six coordinates: the position of its centre of mass and three angles describing molecule orientation. The energy of water solution is given as the sum of kinetic energies (for translation and rotation) and the intermolecular potential which is assumed to be pairwise additive and can be given in several different ways, i.e. the heat bath is given by its Hamiltonian [23][24][25]. I am currently investigating MD models based on Hamiltonian dynamics, with the aim of designing and analysing multi-scale algorithms similar to the algorithm [M1]-[M8] from this paper. The ultimate goal of this research is to design BD models of intracellular processes which make use of modern MD simulations [26,27] in small regions of intracellular space where a BD description is not available or applicable (and an MD level of detail is required). Results will be reported in a future publication.