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.
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 , the MAPK pathway  and signal trasduction in Escherichia coli chemotaxis . 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=[X1,X2,X3] and its diffusion constant by D, a simple model of Brownian motion is given by the (overdamped) Langevin equation 1.1 where Wi, i=1,2,3, are three independent Wiener processes . 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 , MCell , Green's function reaction dynamics  and the first-passage kinetic Monte Carlo method . 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 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  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  or actin dynamics in filopodia ), 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 , a bimolecular reaction occurs whenever the distance of reactants is less than the reaction radius 1.3 where DA (resp. DB) 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 . 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. or . Our main goal is to design and analyse multi-scale 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 [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 This limit can be achieved in many different ways. For example, we can keep m fixed and pass , or we can keep M fixed and pass m→0. In what follows, we define the parameter 1.5 This parameter is dimensionless, even if we assume that 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] 1.6 1.7 where [X1,X2,X3] 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 . We overview the results which relate MD models [A], [B] and [C] with Brownian motion in §2 and §3.
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 xi and velocities vi, i=1,2,3,…, of heat bath particles, and positions Xi 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 . 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 2.1 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, is equal to where (b−a) is the length of the interval [a,b]. Initial velocities of heat bath particles are given by the normal distribution 2.2 Let us consider a model with a single heavy particle, i.e. N=1. Then its location X1 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 vi, their velocities are updated using the conservation of mass and momentum by 2.3 and 2.4 where tildes denote post-collision velocities. Using (1.5), the equations (2.3) and (2.4) can be rewritten as 2.5 The following result can be shown for the above MD model [A].
Let γ>0 and D>0. Let us consider the heavy particle of mass M with initial position Xμ(0)=X0 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 2.6 where X(0)=X0 and V (0)=V 0. That is, X and V solve the one-dimensional version of equations (1.6) and (1.7).
This lemma can be proved using the main theorem in Holley , 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 , 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 xi∈[−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.
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 2.7 The positions x and velocities v of these particles are distributed according to 2.8 where H(⋅) is the Heaviside step function. In particular, the positions of the particles at point are distributed at time t=Δt according to 2.9 where is the complementary error function.
Particles which are at point 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 at time t=Δt is Thus, we proved (2.9). Integrating this formula over x in interval , we obtain the average number of particles which are in the interval at time t=Δt 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 step [A2] by calculating the exact time tc of their collision, i.e. tc∈[t,t+Δt]. Then we use their pre-collision velocities in the time interval [t,tc] and post-collision velocities in the time interval [tc,t+Δt] to compute X(t+Δt) and xi(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 . An acceptance–rejection algorithm for sampling random numbers from π erfc(z) is given in table 2. We use it with the constants a1 and a2 given by 2.10 The values of constants a1 and a2 were numerically estimated to maximize the total acceptance probability 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 2.11 where C1 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 a1 and a2 given by (2.10). In step [A5], we also sample the velocity of the new particle using the truncated Gaussian distribution 2.12 where C2 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 .
In figure 1, we present illustrative results computed by the algorithm [A1]–[A6]. We use μ=103, γ=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 1a, 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 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 2.13 where 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. , as the dot-dashed line in figure 1b. If we neglect the exponential terms in (2.13), we obtain 2.14 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.
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 xi and velocities vi, i=1,2,3,…, of heat bath particles, and positions and velocities , i=1,2,…,N, of heavy particles of mass M≫m, where m is the mass of a heat bath particle . 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  3.1 3.2 where vj 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.
(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].
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 3.3 Let the velocities of heat bath particles be distributed according to 3.4 and v=[v1,v2,v3]. 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
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. . 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 . They were interested in the limit m→0, which is equivalent to . Considering the scaling m3/2f(vm1/2) of the velocity distribution of heat bath particles (with density scaled as λ/m1/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 γ>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 3.5 Let the velocities of heat bath particles be distributed according to 3.6 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  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.
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 3.7 and the average number of particles in the semi-infinite cuboid at time t=Δt is 3.8
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 at time t=Δt is equal to which is formula (3.8).
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 at time t=Δt, 3.9 where (⋅)+ denotes a positive part. Integrating this formula over x in the semi-infinite cuboid , we obtain Substituting (3.5) for λμ 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 2d, 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 2d 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 2d, 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 2a–c, where we consider only one large (heavy) particle, i.e. N=1. The large particle can either be in ΩC (figure 2a), or be in ΩD (figure 2c) or crossing the boundary as is shown in figure 2b. Our geometry is given by (1.4), where 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 4.1 where 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 , 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 4.2 If we explicitly modelled heat bath particles in ΩC, then they would be distributed according to the Poisson distribution with density λμ in the interval . 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 4.3 where λμ and σμ are given by (2.1) and (2.2). In the limit , we have . Thus, we use the Taylor expansion in (4.3) to get 4.4 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 4.5
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 4.6 where C1 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 a1 and a2 given by (2.10). In step [M5], we also sample the velocity of the new particle using the truncated Gaussian distribution 4.7 where C2 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 . 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 4.8
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 105 realizations of the algorithm [M1]–[M8], 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  4.9 In figure 3b, we plot the time evolution of the mean squared displacement. This figure can be directly compared with figure 1b, 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 1b 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 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: 4.11 Substituting (4.4) and (4.5), we derive (1.7). In particular, we have verified the limiting result in lemma 2.1.
5. From three-dimensional molecular dynamics models [B] and [C] to Brownian dynamics
We use a simple multi-scale geometry where domain is divided into two half spaces. Heavy molecules are simulated in both half spaces. In , we use the MD model [B] or [C]. It is coupled with the BD model given by (1.6) and (1.7) in This set-up is a three-dimensional version of multi-scale problems which are schematically drawn in figure 2. Boundary conditions for heat bath particles at the interface 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 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=[X1,X2,X3] with velocity V=[V 1,V 2,V 3] and radius R. It intersects the interface I if X1(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 5.1 where 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 three-dimensional simulations where α(X(t),V(t))Δt is the drift vector, 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 5.2 To estimate drift α and diffusion coefficient β, we separately consider MD models [B] and [C] in §5a and §5b, respectively.
(a) Molecular dynamics model [B]
The following lemma will be useful to estimate the drift coefficient α.
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)=[X1(t),X2(t),X3(t)] and with velocity V(t)=[V 1(t),V 2(t),V 3(t)]. Let y=(y1,y2,y3) be a given point on the surface of the heavy molecule at time t, i.e. 5.3 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+dy) is ψj(y) dy, where 5.4
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 5.5 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 c1>0, and (y−X(t))/R, η2, η3 is the orthornormal basis in Then (5.5) reads as follows: Thus, we have 5.6 Substituting (3.4) for fμ and integrating over τ, c2 and c3, we have Integrating over c1, we deduce (5.4). □
Using lemma 5.1, we can compute the drift coefficient α(X(t),V(t)) in equation (5.2) as follows: 5.7 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 5.8 and 5.9 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 5.10 for i=1,2,3. Substituting (3.4) for fμ, (3.3) for λμ and using β(X,V)=β(X,0), we obtain 5.11 where the last equation can be verified using the same argument as equation (5.10). Note that, by substituting X1=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 §5a, 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 5.12 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 X1=R in (5.12) can be used to verify the limiting result in lemma 3.2.
(c) Illustrative numerical results
In §5a,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 the BD model (1.6) and (1.7) in 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 X1 positions of the heavy particle at time t=1, computed using 105 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 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.
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,L1]×[0,L2]×[0,L2], where L1>0 and L2>0. The cell membrane is modelled by one side of the cuboid, namely which is shaded grey in figure 5. Our goal is to model the binding of diffusing proteins to receptors on the cell membrane with an MD level of detail. Therefore, we define ΩD as a part of the intracellular space which is close to the cell membrane ∂ΩM, i.e. where h>0 and the interface I is at x1=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 . 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 , namely 6.1 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<h2<h3<L1, 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 I3.
Then we revert to the BD model (6.1) which is used in ΩC2.
(ii) The protein molecule crosses the interface I.
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 105 realizations. Initial positions of the protein molecule are uniformly distributed along the X1-axis. The histogram of positions (along the X1-axis) at time t=1 is plotted in figure 6a. Interfaces I, I2 and I3 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 6.2 with boundary conditions  6.3 The solution of (6.2) and (6.3) with uniform initial condition ϱ(x1,0)≡const. is given by the black solid line in figure 6a. As we only visualize the distribution along the X1-axis in figure 6a, we can further decrease the computational cost by truncating the simulation domain in the x2 and x3 directions to the region close to the protein molecule. That is, we only simulate small particles in the subdomain [0,h]×[X2(t′)−h4,X2(t′)+h4]×[X3(t′)−h4,X3(t′)+h4], where t′ is the time when the protein molecule enters ΩD∪ΩC1. This subdomain (moving window) is shifted accordingly whenever X2(t) or X3(t) approach its boundary.
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 . 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 . For example, Rahman & Stillinger  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.
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.
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.