## Abstract

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.

## 1. 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
*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 spatio-temporal 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
_{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 multi-scale 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 [12,5]. 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
*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 [14,12]. 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. *Ω*_{D}⊂*Ω*. We define
*Ω*_{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 [A], [B] and [C] in what follows. More complicated MD approaches are discussed in §7.

The first MD model [A] is introduced in §2. It is a one-dimensional MD model where all particles move along the real line. In particular, the radius *R* does not have to be considered, because it has no influence on the dynamics of large particles. In one dimension, heat bath particles cannot pass each other, which makes the MD model [A] different from three-dimensional models in §3 where heat bath particles (points) do not interact with each other.

In §3, we introduce two three-dimensional models, denoted [B] and [C], where the non-zero radius *R* is one of the key parameters. To make one- and three-dimensional models comparable, we keep *R* fixed in the three-dimensional model and we study the behaviour of all MD models in the limit *m* fixed and pass *M* fixed and pass *m*→0. In what follows, we define the parameter
*M* and *m* have physical units of mass. However, in this paper, all parameters are considered dimensionless for simplicity. We are interested in the limit

All three models [A], [B] and [C] converge in appropriate limits to the Brownian motion of large particles given by (1.1). One can also show that these models converge to the Langevin description [15–17]
*X*_{1},*X*_{2},*X*_{3}] is the position of a diffusing molecule, [*V* _{1},*V* _{2},*V* _{3}] is its velocity, *D* is the diffusion coefficient and *γ* is the friction coefficient. This description can be further reduced to (1.1) in the overdamped limit

Both (1.1) or (1.6)–(1.7) reduce the dimensionality of the problem, making BD less computationally intensive than the corresponding MD simulations. In §4 and §5, we study how MD models [A], [B] and [C] can be used in one part *Ω*_{D} of the computational domain *Ω* and the BD models (1.1) or (1.6)–(1.7) in the remainder *Ω*_{C}, making use of the notation (1.4). We apply our findings to a simplified model of protein binding to receptors in §6. We conclude with discussing our results in §7.

## 2. 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
*μ* 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*, is equal to
*b*−*a*) is the length of the interval [*a*,*b*]. Initial velocities of heat bath particles are given by the normal distribution
*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

### 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 [A] in a finite domain [−*L*,*L*], i.e. we consider a finite number *n*≡*n*(*t*) of heat bath particles which are at positions *x*^{i}∈[−*L*,*L*] with velocities *i*=1,2,…,*n*. We want to formulate boundary conditions of our problem so that the spatio-temporal statistics in [−*L*,*L*] are equivalent to spatio-temporal statistics of the original unbounded process. The following lemma will be useful for designing appropriate boundary conditions.

### Lemma 2.2

*Let* *and* Δ*t*>0. *Let us assume that heat bath particles are distributed according to the Poisson distribution with density* *(2.1)* *in the interval* *Their initial velocities are given according to* *(2.2)* *and there are no particles in the interval* *at time* *t*=0. *Then the average number of particles in the interval* *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* *are distributed at time* *t*=Δ*t* *according to*
*where*

### Proof.

Particles which are at point *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 *t*=Δ*t* is
*x* in interval *t*=Δ*t*
_{μ} 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 step [A2] by calculating the exact time *t*_{c} of their collision, i.e. *t*_{c}∈[*t*,*t*+Δ*t*]. Then we use their pre-collision velocities in the time interval [*t*,*t*_{c}] and post-collision velocities in the time interval [*t*_{c},*t*+Δ*t*] to compute *X*(*t*+Δ*t*) and *x*^{i}(*t*+Δ*t*). All heat bath particles which left the domain [−*L*,*L*] are removed in step [A3].

New heat bath particles are introduced in steps [A4] and [A5]. 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 *π* *erfc*(*z*), where *π* *erfc*(*z*) is given in table 2. We use it with the constants *a*_{1} and *a*_{2} given by
*a*_{1} and *a*_{2} were numerically estimated to maximize the total acceptance probability *ζ*_{2}.

A particle introduced close to the right boundary in step [A5] will have its position sampled according to the probability distribution
*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 *C*_{2} is a normalization constant. To sample random numbers according to the truncated normal distributions in steps [A4] and [A5], we use an acceptance–rejection algorithm which is derived as proposition 2.3 in [18].

In figure 1, we present illustrative results computed by the algorithm [A1]–[A6]. We use *μ*=10^{3}, *γ*=10 and *D*=1. We initialize the position and velocity of the heavy particle as *X*(0)=0 and *V* (0)=0 and we use the algorithm [A1]–[A6] with time step Δ*t*=10^{−7} in the interval [−*L*,*L*] where *L*=20. In figure 1*a*, we present 30 illustrative trajectories of the heavy particle *X*(*t*) computed for *t*∈[0,10]. The mean squared displacement given by the MD model [A] is plotted in figure 1*b* 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
*b*. We also plot the mean squared displacement corresponding to the overdamped limit (1.1), i.e. *b*. If we neglect the exponential terms in (2.13), we obtain
*b* as the dotted line. We will use (2.14) later in §6 to couple the overdamped BD model (1.1) with MD simulations.

## 3. Three-dimensional molecular dynamics models [B] and [C]

MD models [B] and [C] are three-dimensional generalizations of the MD model [A]. They are described in terms of positions **x**^{i} and velocities **v**^{i}, *i*=1,2,3,…, of heat bath particles, and positions *i*=1,2,…,*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.

MD models [B] and [C] are both based on elastic collisions of heavy molecules (balls with mass *M* and radius *R*) with point bath particles with masses *m*. As the collisions are without friction, conservation of momentum and energy then yields the following generalization of formulae (2.5) for post-collision velocities [16]
**v**^{j} is the velocity of the heat bath molecule which collided with the *i*th 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.

### (a) Molecular dynamics model [B]

MD model [B] will use the normal distribution for velocities of heat bath particles. The following lemma generalizes lemma 2.1 to the three-dimensional MD model [B].

### Lemma 3.1

*Let* *γ*>0, *D*>0 *and* *R*>0. *Let us consider the MD model [B] where heat bath particles are distributed according to the Poisson distribution with density*
*Let the velocities of heat bath particles be distributed according to*
*and* **v**=[*v*_{1},*v*_{2},*v*_{3}]. *We will consider one heavy molecule in such a heat bath, i.e*. *N*=1. *Then the position and velocity of the heavy molecule*, **X**_{μ} *and* **V**_{μ}, *converge (in the sense of distribution) to the solution of* *(1.6)* *and* *(1.7)* *in the limit*

### Proof.

The MD model [B] and heat bath distributions (3.3) and (3.4) satisfy the assumptions of theorem 2.1 in Dürr *et al.* [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 *m*^{3/2}*f*(**v***m*^{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].

### Lemma 3.2

*Let* *γ*>0, *D*>0 *and* *R*>0. *Let us consider the MD model [C] where heat bath particles are distributed according to the Poisson distribution with density*
*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.

### Lemma 3.3

*Let* *and* Δ*t*>0. *Let us assume that heat bath particles are distributed according to the Poisson distribution with density* λ_{μ} *in the half space* *their initial velocities are distributed according to* *f*_{μ}(**v**) *and there are no particles in the half space* *at time* *t*=0. *Let us assume that* λ_{μ} *and* *f*_{μ}(**v**) *are either given by* *(3.3)* *and* *(3.4)* (*MD model* [*B*]) *or given by* *(3.5)* *and* *(3.6)* (*MD model* [*C*]).

*Then the positions* **x** *and velocities* **v** *of heat bath particles in the half space* *are distributed at time* *t*=Δ*t* *according to*
*and the average number of particles in the semi-infinite cuboid* *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 (3.3) and (3.4). Integrating (3.7) over positions and velocities (see the proof of lemma 2.2), we conclude that the average number of particles in the semi-infinite cuboid *t*=Δ*t* is equal to

Next, let us consider that λ_{μ} and *f*_{μ}(**v**) are given by (3.5) and (3.6). Integrating (3.7) with respect to **v**, we get the density of particles at *t*=Δ*t*,
_{+} denotes a positive part. Integrating this formula over **x** in the semi-infinite cuboid _{μ} and (3.6) for *σ*_{μ}, we obtain (3.8). □

Lemma 3.3 can be used to specify boundary conditions for simulations of the MD models [B] and [C] in finite domains as we did for the one-dimensional case in lemma 2.2. In §5, we will use lemma 3.3 to develop and analyse multi-scale approaches which can efficiently and accurately compute results with an MD level of detail in a (relatively small) subdomain *Ω*_{D}⊂*Ω* by using coarser BD simulations in the remainder. The geometry of the desired multi-scale method is formulated using (1.4) where an MD model is used in *Ω*_{D}, a coarser BD model is used in *Ω*_{C} and these models are coupled across the interface *I*. The situation is schematically shown in figure 2*d*, which presents a two-dimensional version of our multi-scale set-up. Here, point particles describe heat bath molecules which are used in *Ω*_{D}. Large biomolecules of interest are denoted as grey circles. They are simulated using BD in *Ω*_{C}. The solid line denotes interface *I*.

The schematic in figure 2*d* is presented in two spatial dimensions to better visualize the problem geometry. MD models [B] and [C] are formulated in a three-dimensional physical space. In the three-dimensional version of figure 2*d*, the cloud of point particles would cover the grey ball. To get some insights into this multi-scale problem, we start with the one-dimensional MD model [A].

## 4. 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 2*a*–*c*, where we consider only one large (heavy) particle, i.e. *N*=1. The large particle can either be in *Ω*_{C} (figure 2*a*), or be in *Ω*_{D} (figure 2*c*) or crossing the boundary as is shown in figure 2*b*. Our geometry is given by (1.4), where
*X*(*t*)−*R*,*X*(*t*)+*R*). Let us consider that the large particle intersects the interface *I*, as is shown in figure 2*b*. 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
*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 *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
*Ω*_{C}, then they would be distributed according to the Poisson distribution with density λ_{μ} in the interval *H*(⋅) is the Heaviside step function. Using (2.2), we obtain
_{μ} and *σ*_{μ} are given by (2.1) and (2.2). In the limit

Equations (4.4) and (4.5) are used in the multi-scale algorithm in table 3. The first two steps [M1] and [M2] are the same as [A1] and [A2]. As heat bath particles are only simulated in the subdomain *Ω*_{D}=(−*L*,0), we remove all particles which left *Ω*_{D} during the time interval (*t*,*t*+Δ*t*) in step [M3]. 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
*C*_{1} is a normalization constant. The probability distribution (4.6) can be justified using the same argument as lemma 2.2 and equation (2.11). To sample from the probability 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 *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

In figure 3, we present illustrative results computed by the algorithm [M1]–[M8]. We consider one heavy particle which starts at position *X*(0)=0 with velocity *V* (0)=0, as we did in figure 1. The distribution of its position at time *t*=1, computed using 10^{5} realizations of the algorithm [M1]–[M8], is plotted in figure 3*a*. It is compared with the distribution obtained by the limiting BD model (2.6), which is, for *t*≫*γ*^{−1}, given by [19]
*b*, we plot the time evolution of the mean squared displacement. This figure can be directly compared with figure 1*b*, because we use the same parameter values. The results computed by the multi-scale algorithm [M1]–[M8] compare well with the results given by the BD model (2.6). We have already shown in figure 1*b* that the limiting BD model (2.6) also compares well with the MD simulations. In particular, the algorithm [M1]–[M8] is able to compute results with the MD-level precision by using coarser BD models in a part of the computational domain.

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

## 5. From three-dimensional molecular dynamics models [B] and [C] to Brownian dynamics

We use a simple multi-scale geometry where domain

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 three-dimensional 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
*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 three-dimensional simulations where ** α**(

**X**(

*t*),

**V**(

*t*))Δ

*t*is the drift vector,

**ξ**=[ξ

_{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

**and diffusion coefficient**

*α***, we separately consider MD models [B] and [C] in §5**

*β**a*and §5

*b*, respectively.

### (a) Molecular dynamics model [B]

The following lemma will be useful to estimate the drift coefficient ** α**.

### Lemma 5.1

*Let* *γ*>0, *D*>0, *R*>0 *and* Δ*t*>0. *Let us consider the MD model [B] where the positions and velocities of heat bath particles are distributed according to* (3.3) *and* (3.4). *Let us consider one heavy molecule in such a heat bath, i.e.* *N*=1, *with the position of its centre to be at* **X**(*t*)=[*X*_{1}(*t*),*X*_{2}(*t*),*X*_{3}(*t*)] *and with velocity* **V**(*t*)=[*V* _{1}(*t*),*V* _{2}(*t*),*V* _{3}(*t*)]. *Let* **y**=(*y*_{1},*y*_{2},*y*_{3}) *be a given point on the surface of the heavy molecule at time* *t*, *i.e.*
*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* (**y**,**y**+*d***y**) is *ψ*_{j}(**y**) *d***y**, where

### 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 pre-collision 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
**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
*c*_{1}>0, **y**−**X**(*t*))/*R*, ** η_{2}**,

**is the orthornormal basis in**

*η*_{3}*f*

_{μ}and integrating over

*τ*,

*c*

_{2}and

*c*

_{3}, we have

*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:

*S*(

**X**(

*t*)) is the part of the surface of the heavy molecule which intersects the BD subdomain

*Ω*

_{C}, i.e.

*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

*i*=1,2,3. Substituting (3.4) for

*f*

_{μ}, (3.3) for λ

_{μ}and using

**(**

*β***X**,

**V**)=

**(**

*β***X**,

**0**), we obtain

*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 (5.6), (5.7) and (5.10), which are derived in §5*a*, are applicable to both MD models [B] and [C]. To estimate the drift coefficient ** α**(

**X**,

**V**) for the MD model [C], we substitute (3.5) for λ

_{μ}and (3.6) for

*f*

_{μ}in (5.6) and (5.7). We obtain

*α*

_{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.

### (c) Illustrative numerical results

In §5*a*,*b*, we have observed that the only difference between MD models [B] and [C] is a different formula for the coefficient *α*_{1}(**X**,**V**) in (5.2), given by (5.8) and (5.12), respectively. The remaining terms in (5.2) are the same, given by (5.9) and (5.11). In this section, we present an illustrative computation with the MD model [C], but the same results can also be obtained with the MD model [B] (results not shown). An illustrative computation with the MD model [B] is presented later in §6.

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 ** α** 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 4

*a*. The limiting BD result is again given by (4.9). In figure 4

*b*, 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

## 6. 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 figure 5. Our computational domain is a part of the intracellular space next to the cell membrane given as the cuboid *Ω*=[0,*L*_{1}]×[0,*L*_{2}]×[0,*L*_{2}], where *L*_{1}>0 and *L*_{2}>0. The cell membrane is modelled by one side of the cuboid, namely
*Ω*_{D} as a part of the intracellular space which is close to the cell membrane ∂*Ω*_{M}, i.e.
*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 *Ω*_{D}. We couple these models using the intermediate BD model (1.6)–(1.7). We introduce two additional interfaces
*h*<*h*_{2}<*h*_{3}<*L*_{1}, as shown in figure 5. We denote
*Ω*_{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.

In figure 6, we present illustrative results computed by averaging over 10^{5} realizations. Initial positions of the protein molecule are uniformly distributed along the *X*_{1}-axis. The histogram of positions (along the *X*_{1}-axis) at time *t*=1 is plotted in figure 6*a*. Interfaces *I*, *I*_{2} and *I*_{3} are also shown in this plot. As we used a very simple model of the protein binding, we can compare it with the mean-field limit given by the solution of the PDE
*ϱ*(*x*_{1},0)≡const. is given by the black solid line in figure 6*a*. As we only visualize the distribution along the *X*_{1}-axis in figure 6*a*, 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 6*b*. It again compares well with the results obtained by the limiting PDE system (6.2) and (6.3).

## 7. 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.6)–(1.7) as 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–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.

## Funding statement

The research leading to these results has received funding from the European Research Council under the *European Community's* Seventh Framework Programme *(FP7/2007-2013)*/ERC *grant agreement* no. 239870.

## Acknowledgements

I would like to thank the Royal Society for a University Research Fellowship; Brasenose College, University of Oxford, for a Nicholas Kurti Junior Fellowship and the Leverhulme Trust for a Philip Leverhulme Prize.

- Received January 14, 2014.
- Accepted March 28, 2014.

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.