Coupling all-atom molecular dynamics simulations of ions in water with Brownian dynamics

Molecular dynamics (MD) simulations of ions (K+, Na+, Ca2+ and Cl−) in aqueous solutions are investigated. Water is described using the SPC/E model. A stochastic coarse-grained description for ion behaviour is presented and parametrized using MD simulations. It is given as a system of coupled stochastic and ordinary differential equations, describing the ion position, velocity and acceleration. The stochastic coarse-grained model provides an intermediate description between all-atom MD simulations and Brownian dynamics (BD) models. It is used to develop a multiscale method which uses all-atom MD simulations in parts of the computational domain and (less detailed) BD simulations in the remainder of the domain.


Introduction
Molecular dynamics (MD) simulations of ions in aqueous solutions are limited to modelling processes in relatively small domains containing (only) several thousands of water molecules [1,2]. Ions play important physiological functions in living cells which typically consist of 10 10 to 10 12 water molecules. In particular, processes which include transport of ions between different parts of a cell cannot be simulated using standard all-atom MD approaches. Coarser models are instead used in applications. Examples include Brownian dynamics (BD) simulations [3] and mean-field Poisson-Nernst-Planck equations [4]. In BD methods, individual trajectories of ions are described using where X = (X 1 , X 2 , X 3 ) is the position of the ion, D is its diffusion constant and W i , i = 1, 2, 3, are three independent Wiener processes [5]. BD description (1.1) does not explicitly include solvent molecules in the simulation. Moreover, in applications, equation (1.1) can be discretized using a (nanosecond) timestep which is much larger than the typical timestep of MD simulations (femtosecond) [6]. This makes BD less computationally intensive than the corresponding MD simulations.
Longer timesteps of BD simulations enable efficient simulations of ion transport between different parts of the cell, but they limit the level of detail which can be incorporated into the model. For example, intracellular calcium is regulated by the release of Ca 2+ ions from the endoplasmic reticulum via inositol-4,5-triphosphate receptor (IP 3 R) channels. BD models in the literature use equation (1.1) to describe trajectories of calcium ions [3,7]. The conformational changes between the open and closed states of IP 3 R channels are controlled by the binding of Ca 2+ to activating and inhibitory binding sites. BD models postulate that binding of an ion occurs with some probability whenever the distance between the ion and an empty site is less than the specific distance, the so-called reaction radius [8,9]. Although details of the binding process are known [10,11], they cannot be incorporated into coarse BD models of calcium dynamics, because equation (1.1) does not correctly describe short-time behaviour of ion dynamics.
The calcium-induced calcium release through IP 3 R channels is an example of a multiscale dynamical problem where MD simulations are important only in certain parts of the computational domain (close to an IP 3 R channel), while in the remainder of the domain a coarser, less detailed, BD method could be used (to describe trajectories of ions). Such multiscale problems cannot be simulated using MD methods, but there is potential to design multiscale computational methods which compute the desired information with an MD-level of resolution by using MD and BD models in different parts of the computational domain [12].
In [12], three relatively simple and analytically tractable MD models are studied (describing heat bath molecules as point particles) with the aim of developing and analysing multiscale methods which use MD simulations in parts of the computational domain and less detailed BD simulations in the remainder of the domain. In this follow-up paper, the same question is investigated in all-atom MD simulations which use the SPC/E model of water molecules. In order to couple MD and BD simulations, we need to first show that the MD model is in a suitable limit described by a stochastic model which does not explicitly take into account heat bath (water) molecules. In [12], this coarser description was given in terms of Langevin dynamics. Considering all-atom MD simulations, the coarser stochastic model of an ion is more complicated than Langevin dynamics. In this paper, it will be given by (1.4) and where X ≡ (X 1 , X 2 , X 3 ) is the position of the ion, V ≡ (V 1 , V 2 , V 3 ) is its velocity, U ≡ (U 1 , U 2 , U 3 ) is its acceleration, Z ≡ (Z 1 , Z 2 , Z 3 ) is an auxiliary variable, dW ≡ (dW 1 , dW 2 , dW 3 ) is white noise and η j , j = 1, 2, 3, 4, are parameters. These parameters will be chosen according to all-atom MD simulations as discussed in §3. In §4, we show that (1.2)-(1.5) provides a good approximation of ion behaviour. In §5, we further analyse the system (1. coarse-grained model (1.2)-(1.5) and can be used to fit additional properties of all-atom MD. This generalization is formulated in terms of fictitious particle models [13] which are themselves based on earlier work in non-equilibrium statistical physics [14][15][16][17]. In particular, we show that the dynamics of the intermediate coarse-grained model (1.2)-(1.5) can be equivalently described by an appropriately formulated fictitious particle model [13]. We conclude by discussing related methods developed in the literature in §8.

Molecular dynamics simulations of ions in SPC/E water
There have been several MD models of liquid water developed in the literature. The simplest models (e.g., SPC [18], SPC/E [19] and TIP3P [20]) include three sites in total, two hydrogen atoms and an oxygen atom. More complicated water models include four, five or six sites [21,22]. In this paper, we use the three-site SPC/E model of water which was previously used for MD simulations of ions in aqueous solutions [1,23]. In the SPC/E model, the charges (q h = 0.4238 e) on hydrogen sites are at 1 Å from the Lennard-Jones centre at the oxygen site which has negative charge q o = −0.8476 e. The HOH angle is 109.47 • . We use the RATTLE algorithm [24] to satisfy constraints between atoms of the same water molecule.
We investigate four ions (K + , Na + , Ca 2+ and Cl − ) at 25 • C using MD parameters presented in [23]. Let us consider a water molecule and let us denote by r i0 (respectively, r i1 and r i2 ) the distance between the ion and the oxygen site (respectively, the first and second hydrogen sites). The pair potential between the water molecule and the ion is then given by [1,23], where A io and B io are Lennard-Jones parameters between the oxygen on the water molecule and the ion, k e is Coulomb's constant and q i is the charge on the ion. The values of parameters are given for four ions considered in table 1. We express mass in daltons (Da), length in ångströms (Å) and time in picoseconds (ps), consistently in the whole paper. Using these units, the parameters of the Lennard-Jones potential between the oxygen sites on two SPC/E water molecules are A oo = 2.6334 × 10 8 Da Å 14 ps −2 and B oo = 2.6171 × 10 5 Da Å 8 ps −2 .
We consider a cube of side L = 24.83 Å containing 511 water molecules and 1 ion, i.e. we have 8 3 = 512 molecules in our simulation box. In the following section, we use standard NVT simulations where the temperature is controlled using Nosé-Hoover thermostat [25,26] and the number of particles is kept constant by implementing periodic boundary conditions. In particular, we assume that our simulation box is surrounded by periodic copies of itself. Then the longrange (Coulombic) interactions can be computed using several different approaches, including the Ewald summation or the reaction field method [27,28]. We use the cut-off sphere of radius L/2 and the reaction field correction as implemented in [1]. This approach is more suitable for multiscale methods (studied later in §6) than the Ewald summation technique. The MD timestep is for all MD simulations in this paper chosen as t = 10 −3 ps = 1 fs.

Parametrization of the coarse-grained model of ion
In MD simulations, an ion is described by its position X ≡ (X 1 , X 2 , X 3 ) and velocity V ≡ and where M is the mass of the ion (given in This value is reported in table 2. In the same way, the reported values of V 2 i are computed as averages over all three dimensions. Diffusion constant D can be estimated by calculating mean square displacements or velocity autocorrelation functions. In table 2, we report the values of D which were estimated in [1] by calculating mean square displacements. Let us consider the coarse-grained model (1.2)-(1.5) and let · denotes an average over many realizations of a stochastic process. Multiplying equations (1.3) and (1.4) by V i and U i , respectively, we obtain the following ordinary differential equations (ODEs) for second moments: Consequently, we obtain that U i V i = 0 and U i Z i = 0 at steady state. Multiplying equations (1.3)-(1.5) by V i , U i and Z i , respectively, and taking averages, we obtain Using U i V i = 0 and U i Z i = 0, we obtain that V i Z i = 0 at steady state and This equation is used in  calculating the second moment of This value is reported in the last column of Using U i Z i = 0 and V i Z i = 0, we obtain at steady state Multiplying equation (1.2) by X i , V i , U i and Z i and equations (1.3)-(1.5) by X i and taking averages, we obtain the following system of ODEs for second moments: Consequently, at equilibrium, we obtain Using (3.7) and (3.10), we have Consequently, we obtain at steady state Using (3.15), we get The values calculated by (

Accuracy of the coarse-grained model of ion
The coarse-grained model (1.2)-(1.5) has four parameters η j , j = 1, 2, 3, 4. To parametrize this model, we have used four quantities estimated from detailed MD simulations, diffusion constant D and steady-state values of V 2 i , U 2 i and Z 2 i . In particular, the coarse-grained model (1.2)-(1.5) will give the same values of these four quantities, including the value of diffusion constant D which is the sole parameter of the BD model (1.1). In this section, we explain why the coarsegrained description given by (1.2)-(1.5) can be used as an intermediate model to couple BD and MD models.
We begin by illustrating why Langevin dynamics (which is used in [12] for a similar multiscale problem) is not suitable for all-atom MD simulations studied in this paper. In [12], a few (heavy) particles with mass M and radius R are considered in the heat bath consisting of 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. In this case, it can be shown that the description of heavy particles converges in an appropriate limit to Brownian motion given by equation (1.1). One can also show that the model converges to Langevin dynamics (in the limit m/M → 0) [29][30][31]: and is its velocity, D is the diffusion coefficient and γ is the friction coefficient. In [12], Langevin dynamics (4.1)-(4.2) is used as an intermediate model which enables the implementation of BD description (1.1) and the original detailed model in different parts of the computational domain. Langevin dynamics (4.1)-(4.2) describes a diffusing particle in terms of its position and velocity, i.e. it uses the same independent variables for the description of an ion as the MD model (3.1)-(3.2). Langevin dynamics can be further reduced to BD model (1.1) in the overdamped limit γ → ∞. However, it cannot be used as an intermediate model between BD and all-atom MD simulations considered in this paper, because it does not correctly describe the ion behaviour at times comparable with the MD timestep t. To illustrate this, let us parametrize Langevin dynamics (4.1)-(4.2) using diffusion constant D and the second velocity moment V 2 i estimated from all-atom MD simulations. To get the same second moment of velocity, Langevin dynamics requires that we choose Discretizing equation (4.2), the ion acceleration during one timestep is where (ξ 1 , ξ 2 , ξ 3 ) is a vector of normally distributed random numbers with zero mean and unit variance. Using (4.3), the second moment of the right-hand side of (4.4) is Using the MD values of D and V 2 i for K + which are given in table 2 and using MD timestep t = 10 −3 ps, we obtain that the second moment (4.5) is equal to 4.44 × 10 5 Å 2 ps −4 . On the other hand, U 2 i estimated from all-atom MD simulations and given in table 2 is 4.86 × 10 3 Å 2 ps −4 which is one hundred times smaller. The main reason for this discrepancy is that Langevin dynamics postulates that the random force in equation (4.2) acting on the particle at time t is not correlated to the random force acting on the particle at time t + t. However, this is not true for all-atom MD simulations where random force terms at subsequent timesteps are highly correlated. Since Langevin dynamics is not suitable for coupling MD and BD models, we need to introduce a stochastic model of ion behaviour which is more complicated than Langevin dynamics. The coarse-grained model (1.2)-(1.5) studied in this paper is a relatively simple example of such a model. Its parametrization, discussed in §3, guarantees that the coarse-grained model (1.2)-(1.5) well approximates all-atom MD simulations at equilibrium. They both have the same value of diffusion constant D and steady-state values of V 2 i , U 2 i and Z 2 i . Next, we show that the coarsegrained model (1.2)-(1.5) also compares well with all-atom MD simulations at shorter timescales. We consider the rate of change of acceleration (jerk or the scaled derivative of force). We define the average jerk as a function of current velocity and acceleration of the ion: To estimate J(v, u) from all-atom MD simulations, we calculate the rate of change of acceleration during each MD timestep i.e. we run a long (nanosecond) MD simulation, calculate the values of (U i (t + t) − U i (t))/ t during every timestep and record their average in two-variable array Since the estimated J(v, u) only weakly depends on u, we visualize our results in figure 1 using two functions of one variable, v, namely where p u (u) is the steady-state distribution of U i estimated from the same long-time MD trajectory. As before, we use all three dimensions to calculate the averages J(v, u) and p u (u). Function J 1 (v) (which gives jerk at the most likely value of U i ) is plotted using crosses and function J 2 (v), the average over U i variable, is plotted using circles in figure 1. In order to compare all-atom MD simulations with the coarse-grained model (1.2)-(1.5), we calculate the corresponding jerk matrix J(v, u) for the coarse-grained model. We denote by p(v, u, z) the stationary distribution of the stochastic processes (1. . Then the jerk matrix (4.6) of the coarse-grained Using (1.4), we rewrite it as The stationary distribution p(v, u, z) of (1.3)-(1.5) is Gaussian with mean (0, 0, 0) and stationary covariance matrix: -

From the coarse-grained model to Brownian dynamics
Let us consider the three-variable subsystem (1.3)-(1.5) of the coarse-grained model. Denoting where T stands for transpose, equations (1.3)-(1.5) can be written in vector notation as follows: where matrix B ∈ R 3×3 and vector b ∈ R 3 are given as Let us denote the eigenvalues and eigenvectors of B as λ j and ν j = (ν 1j , ν 2j , ν 3j ) T , j = 1, 2, 3, respectively. The eigenvalues of B are the solutions of the characteristic polynomial Since η 1 , η 2 and η 3 are positive parameters, we conclude that real parts of all three eigenvalues are negative and lie in interval (−η 2 , 0). Using the values of η j , j = 1, 2, 3, given in  general solution of the SDE system (5.1) can be written as follows [32]: where c ∈ R 3 is a constant vector determined by initial conditions and matrix Φ(t) ∈ R 3×3 is given as There are two important conclusions of the above analysis. First of all, eigenvalues λ j , j = 1, 2, 3, given in table 4 satisfy where Re denotes the real part of a complex number. There is a spectral gap between the first eigenvalue and the complex conjugate pair of eigenvalues. If we used this spectral gap, we could reduce the system to two evolution equations for times t 1/|λ 1 |. However, there is no spectral gap to reduce the system to Langevin dynamics (4.1)-(4.2). In particular, we again confirm our conclusion that a coarse-grained approximation of ion behaviour is not given in terms of Langevin dynamics. Our second conclusion is that on a picosecond time scale, we can assume stationarity in (5.1) to get  .14) and (3.16). We consider (deterministic) zero initial conditions, i.e. X i (0) = V i (0) = U i (0) = Z i (0) = 0. All moments are then initially equal to zero. We plot the mean square displacement X 2 i as a function of time. We compare it with the mean square displacement of BD model (1.1) which is given as 2Dt. We observe that there is an approximately constant shift, denoted t * 1 , between both solutions for times t > 0.2 ps. We illustrate this further by plotting X where (ξ 1 , ξ 2 , ξ 3 ) is a vector of normally distributed random numbers with zero mean and unit variance. We use discretization (5.5) of BD model (1.1) in this paper. BD time step T has to be chosen much larger than the MD timestep t. We use T = 0.5 ps, but any larger timestep would also work well. We could also use a variable timestep, as implemented in the Green's-function reaction dynamics [35]. In §6, we consider all-atom MD simulations in domain Ω ⊂ R 3 . Our main goal is to design a multiscale approach which can compute spatio-temporal statistics with the MD-level of detail in relatively small subdomain Ω 1 ⊂ Ω by using BD model (5.5) in the most of the rest of the computational domain. This is achieved by decomposing domain Ω into five subdomains Ω j , j = 1, 2, 3, 4, 5 (see equation (6.1) and discussion in §6). We use MD in Ω 1 , the coarse-grained model (1.2)-(1.5) in Ω 3 and the BD model (5.5) in Ω 5 . The remaining two subdomains, Ω 2 and Ω 4 , are two overlap (hand-shaking) regions where two different simulation approaches can be used at the same time [12,36]. In the rest of this section, we focus on simulations in region Ω 3 ∪ Ω 4 ∪ Ω 5 which concerns coupling the coarse-grained model (1.2)-(1.5) with the BD model (5.5). We use the coarse-grained model in Ω 3 ∪ Ω 4 and the BD model (5.5) in Ω 4 ∪ Ω 5 . In particular, we use both models in the overlap region Ω 4 . Each particle which is initially in Ω 3 is simulated according to (1.2)-(1.5) (discretized using timestep t) until it enters Ω 5 . Then we use (5.5) to evolve the position of a particle (over BD timesteps of length T) until it again enters Ω 3 when we switch the description back from the BD model to the coarse-grained model. In order to do this, we have to initialize variables V i , U i and Z i , i = 1, 2, 3. We use deterministic initial conditions, V i (0) = U i (0) = Z i (0) = 0, discussed above.
In figure 2b, we present an illustrative simulation where Ω 3 ∪ Ω 4 ∪ Ω 5 = R 3 for simplicity. i.e. the overlap region Ω 4 is equal to one bin (visualized as a green bar). Grey (respectively, blue) bars show the density of ions in Ω 3 (respectively, Ω 5 ). We compare our results with the analytical distribution computed for BD description (1.1) at time t = 10 3 ps given by The computed histogram compares well with (5.6), although we can observe a small error: the green bar is slightly taller than the corresponding value of (5.6). If we wanted to further improve the accuracy, we could take into account that there is time shift t * 1 , discussed above, introduced to the multiscale approach by using the deterministic initial conditions, V i (0) = U i (0) = Z i (0) = 0, for ions entering domain Ω 3 . Another possibility is to sample the initial condition for V i , U i and Z i from a suitable distribution. If we use the stationary distribution of subsystem (1.3)-(1.5), then V 2 i does not evolve and is equal to Substituting this constant for V 2 i into (3.12), the system of 10 ODEs for second moments of (1.2)-(1.5) simplifies to four ODEs (3.11)-(3.14). Solving system (3.11)-(3.14) with zero initial conditions (assuming X i (0) = 0), we can again compute the mean square displacement. As in figure 2a, it can be shifted in time to better match with the BD result, 2Dt. We denote this time shift as t * 2 . Its values are given in table 4. We observe that t * 2 is negative and t * 1 is positive for all four ions considered in table 4. Both time shifts t * 1 and t * 2 (together with optimizing size h of the overlap region) could be used to further improve the accuracy of multiscale simulations in Ω 3 ∪ Ω 4 ∪ Ω 5 [12]. However, our main goal is to introduce a multiscale approach which can use all-atom MD simulations in Ω 1 . Since MD simulations are computationally intensive, we will only consider 100 realizations of the multiscale method in §6. In particular, the Monte Carlo error will be larger than the error observed in figure 2b. Thus, we can use the above approach in Ω 3 ∪ Ω 4 ∪ Ω 5 without introducing observable errors in the multiscale method developed in the next section.

Coupling all-atom MD and BD
Let us consider all-atom MD in domain Ω ⊂ R 3 which is so large that direct MD simulations would be too computationally expensive. Let us assume that a modeller only needs to consider the MD-level of detail in a relatively small subdomain Ω 1 ⊂ Ω, while, in the rest of the computational domain, ions are transported by diffusion and BD description (1.1) is applicable. For example, domain Ω 1 could include binding sites for ions or (parts of) ion channels. In this paper, we do not focus on a specific application. Our goal is to show that the coarse-grained model (1.2)-(1.5) is an intermediate model between all-atom MD and BD which enables the use of both methods during the same dynamic simulation. To achieve this, we decompose domain Ω into five subdomains, denoted as Ω j , j = 1, 2, 3, 4, 5 (as it is schematically illustrated in figure 3). These sets are considered pairwise disjoint (i.e. Ω i ∩ Ω j = ∅ for i = j) and In our illustrative simulations, we consider the behaviour of one ion. If the ion is in Ω 1 , then we use all-atom MD simulations as described in §2. In particular, the force between the ion and a water molecule is obtained by differentiating potential (2.1), provided that the distance between the ion and the water molecule is less than the cut-off distance (L/2). Let us denote the force exerted by the ion on the water molecule by F iw (r i0 , r i1 , r i2 ), where r i0 (respectively, r i1 and r i2 ) is the distance between the ion and the oxygen site (respectively, the first and second hydrogen sites) on the water molecule. We use periodic boundary conditions for water molecules in Ω 1 . Whenever the ion leaves Ω 1 , it enters Ω 2 where we simulate its behaviour using the coarsegrained model (1.2)-(1.5). We still simulate water molecules in Ω 1 and we allow them to Ion simulated using all-atom MD in W 1 .
Ion described by BD model (5.5) in W 4 and W 5 .
If an ion is in W 1 or W 2 , then water molecules are simulated using all-atom MD in W 1 .

Figure 3.
Schematic of multiscale set-up. Note that the schematic is drawn in two spatial dimensions to enable better visualization, but all models are formulated and simulated in three spatial dimensions. (Online version in colour.) experience additional forces exerted by the ion which is present in Ω 2 . These forces have the same functional form, F iw , as in MD, but they have modified arguments as follows: where ω ≥ 0 is a parameter and dist(X, Ω 1 ) is the (closest) distance between the ion at position X and subdomain Ω 1 . If the ion is in region Ω 3 ∪ Ω 4 ∪ Ω 5 , then water molecules in Ω 1 are no longer simulated. We use the coarse-grained model (1.2)-(1.5) to simulate the ion behaviour in Ω 3 and the BD model (5.5) in Ω 5 . Overlap region Ω 4 is used to couple these simulation methods as explained in §5. In §5, we have already presented illustrative simulations to validate the multiscale modelling strategy chosen in region Ω 3 ∪ Ω 4 ∪ Ω 5 . Next, we focus on testing and explaining the multiscale approach chosen to couple region Ω 1 with Ω 2 . The key idea is given by force term (6.2) which is used for MD simulations of water molecules in Ω 1 when an ion is in Ω 2 . This force term has two important properties: (i) If an ion is on the boundary of Ω 1 , i.e. X ∈ ∂Ω 1 , then dist(X, Ω 1 ) = 0 and force (6.2) is equal to force term F iw (r i0 , r i1 , r i2 ) used in Ω 1 . (ii) If ω dist(X, Ω 1 ) ≥ L/2, then force (6.2) is equal to zero. Property (i) implies that formula (6.2) continuously extends the force term used in MD. In particular, water molecules do not experience abrupt changes of forces when the ion crosses boundary ∂Ω 1 . Property (ii) is a consequence of the cut-off distance used (together with the reaction field correction [1]) to treat long-range interactions. In our illustrative simulations, we use and Property (ii) implies that extra force (6.2) is equal to zero on boundary ∂Ω 2 \ ∂Ω 1 which is the boundary between regions Ω 2 and Ω 3 . This is consistent with the assumption that ions in region Ω 3 ∪ Ω 4 ∪ Ω 5 do not interact with water molecules in region Ω 1 .
If an ion is in Ω 1 , we use all-atom MD as formulated in §2. Periodic boundary conditions are implemented in MD simulations. Water molecules are subject to forces exerted not only by the The boundary between Ω 1 and Ω 2 is visualized using the black dashed line. We use ω = 1 in (6.2). (b) The mean square displacement in the first coordinate of K + ion simulated in Ω 1 ∪ Ω 2 and computed as the average of 100 realizations for ω = 1 (blue circles), ω = 2 (black crosses) and ω = 10 (green squares).
ion at its real position X in Ω 1 , but also by its copies at periodic locations X + (iL, jL, kL), where i, j, k ∈ Z. When the ion moves to Ω 2 , one of its copies is in Ω 1 . Force term (6.2) is designed in such a way that the strength of interaction decreases (for every copy of the ion) with the distance, dist(X, Ω 1 ), between the real position of the ion and Ω 1 . In particular, force term (6.2) ensures that there are continuous changes of all forces when the ion moves between regions Ω 1 , Ω 2 and Ω 3 .
In figure 4, we present results of simulations of K + ion in region Ω 1 ∪ Ω 2 . We consider 100 realizations of a multiscale simulation with one ion. Its initial position is X(0) = (−L/2, 0, 0) which lies on boundary ∂Ω 1 . We simulate each realization for time 10 ps which is short enough that all trajectories stay inside the ball of radius L/2 centred at X(0). Then X 1 -coordinate of the trajectory determines whether the ion is in Ω 1 or Ω 2 . If X 1 (t) ≥ −L/2, then the ion is in Ω 1 and it is simulated using all-atom MD. If X 2 (t) < −L/2, then the ion is in Ω 2 and evolves according to the coarsegrained model (1.2)-(1.5). In figure 4a, we use (6.2) with ω = 1 and plot X 1 coordinates of all 100 realizations. We observe that the computed trajectories spread on both sides of boundary ∂Ω 1 (dashed line) without any significant bias. The mean square displacement is presented in figure 4b for three different values of ω. The results compare well with (2Dt) 1/2 which is the mean square displacement of one coordinate of the diffusion process.
We conclude with illustrative simulations which are coupling all-atom MD with BD. We use domain Ω ∈ R 3 decomposed into five regions as in equation (6.1), where Ω 1 and Ω 2 are given by (6.3), and , , where ω = 10, h 1 = L/20 and h 2 = L/10. Then the BD domain is Ω 5 = R 3 \ [−7L/10, 7L/10] 3 . We place an ion at the origin (centre of MD domain Ω 1 ), i.e. X(0) = (0, 0, 0), and we simulate each trajectory until it reaches the distance 4L = 99.32 Å from the origin. Let T (r) be the time when a trajectory first reaches distance r from the origin. In figure 5, we plot escape time T (r) as a function of distance r. We plot the value of T (r) for each realization as a blue point. The largest computed escape times (for r = 4L) are 38 506 ps for K + and 47 212 ps for Na + . They are outside the range of panels in figure 5, but the majority of data points are included in this figure. We also  K + Na + Figure 5. Escape time T (r) to reach distance r from the origin computed by the multiscale method. We consider (a)K + ion; and (b) Na + ion. We plot escape times for individual realizations (blue points), the mean escape time estimated from 100 realizations (red solid line) and the theoretical 95% confidence interval (6.4) (green area). We use ω = 10 in (6.2).
plot average T (r) (red solid line) together with 95% confidence intervals. They are compared with theoretical results obtained for the BD model (1.1). The escape time distribution for the BD model (1.1) has mean equal to T (r) = r 2 /(6D) and standard deviation r 2 /(3 √ 10 D). The corresponding theoretical 95% confidence interval (for 100 samples) is This interval is visualized as the green area in figure 5. We note that it would be relatively straightforward to continue the presented multiscale computation and simulate ion diffusion in domains covering the whole cell. The most computationally intensive part is all-atom MD simulation in Ω 1 ∪ Ω 2 . However, once the ion enters Ω 5 , we can compute its trajectory very efficiently. We could further increase the BD timestep in parts of Ω 5 which are far away from Ω 4 , or we could use event-based algorithms, like Green's-function reaction dynamics [35] or First-passage kinetic Monte Carlo method [37], to compute the ion trajectory in region Ω 5 .

A hierarchy of stochastic coarse-grained models
In §6, the coarse-grained model (1.2)-(1.5) has been used in intermediate regions (denoted Ω 2 , Ω 3 and Ω 4 in figure 3) to couple all-atom MD and BD. It is a simple model which has all the necessary properties for this task. The developed multiscale algorithm enables the use of MD together with modern spatio-temporal simulation algorithms for intracellular processes. The coarse-grained model (1.2)-(1.5) (parametrized by four constants η 1 , η 2 , η 3 and η 4 ) provides a better description of ion dynamics than BD (parametrized by one constant D). However, it does not capture all the details of MD. In this section, we introduce a hierarchy of coarse-grained models which generalize (1.2)-(1.5) and include an increasing number of parameters. We illustrate that these models can be used to fit additional characteristics of all-atom MD simulations. There has been a lot of approaches in the literature to develop coarse-grained models at equilibrium by constructing suitable coarse-grained potential energy functions with fewer degrees of freedom [38][39][40]. However, the same coarse-grained potential energy function can correspond to many different dynamic coarse-grained models. In particular, the conservative Hamiltonian dynamics of coarse-grained variables on a coarse-grained potential surface is usually a poor approximation of the real dynamics of all-atom MD models [41,42]. In recent work, Davtyan et al. [13] use fictitious particles with harmonic interactions with coarse-grained degrees of freedom (i.e. they add quadratic terms to the potential function of the system) to improve the fit between an MD model and the dynamics on a coarse-grained potential surface. Each fictitious particle is also subject to a friction force and noise. We will apply the fictitious particle approach to develop a hierarchy of coarse-grained models which generalize (1.2)-(1.5).
Let us consider N fictitious particles interacting with an ion at position X ≡ (X 1 , X 2 , X 3 ) with velocity V ≡ (V 1 , V 2 , V 3 ). LetX j ≡ (X j,1 ,X j,2 ,X j,3 ) be the position of the jth fictitious particle, j = 1, 2, . . . , N. LetṼ j ≡ (Ṽ j,1 ,Ṽ j,2 ,Ṽ j,3 ) be its velocity. Then the fictitious particle model [13] of an ion can be written as the following set of 6(N + 1) equations (for i = 1, 2, 3 and j = 1, 2, . . . , N) and where white noise vectors dW j ≡ (dW j,1 , dW j,2 , dW j,3 ) are mutually independent. Each fictitious particle can be characterized by four constants: three of them correspond to three force terms on the right-hand side of (7.4) and the fourth one is the fictitious particle mass. In order to write the fictitious particle model in a similar way as other models in this paper, we adsorbed the masses of the ion and fictitious particles into the constants on the right-hand sides of equations (7.2) and (7.4). In particular, the jth fictitious particle is characterized by four positive constants α j,1 , α j,2 , α j,3 and α j,4 . To connect the hierarchy of the stochastic fictitious particle models (7.1)- (7.4) with the coarse-grained model (1.2)-(1.5), we introduce new variablesŨ j,i = α j,1 (X j,i − X i ) and Z j,i = α j,1Ṽj,i . Then equations (7.2)-(7.4) read as follows: j,i dt, (7.5) dŨ j,i = (Z j,i − α j,1 V i ) dt (7.6) and dZ j,i = −(α j,2Zj,i + α j,3Ũj,i ) dt + α j,1 α j,4 dW j,i . (7.7) Using N = 1 and denoting U i =Ũ 1,i and Z i =Z 1,i , we obtain the coarse-grained model (1.2)-(1.5) where η 1 = α 1,1 , η 2 = α 1,2 , η 3 = α 1,3 and η 4 = α 1,1 α 1,4 . Thus, the coarse-grained model (1.2)-(1.5) is equivalent to the fictitious particle model (7.1)-(7.4) for N = 1. Moreover, stochastic coarse-grained models in the hierarchy (7.1)-(7.4) provide generalizations of the coarse-grained model (1.2)-(1.5) for N ≥ 2 and can be used to fit additional details of MD simulations. One commonly used MD characteristic is the velocity autocorrelation function In figure 6, we plot the velocity autocorrelation function estimated from MD simulations of K + and Na + ions. Using (5.3), we can find an analytical expression for the velocity autocorrelation function of the coarse-grained model (1.2)-(1.5) as follows:  where matrix Φ(t) ∈ R 3×3 is given as , where λ j , j = 1, 2, 3, are eigenvalues of matrix B given by (5.2) and ν j = (ν 1j , ν 2j , ν 3j ) T are the corresponding eigenvectors. In figure 6, we plot velocity autocorrelation functions (7.9) as green dotted lines for K + and Na + ions. Using t = 0 in (7.9), we obtain C(0) = V 2 i . Since we have parametrized the coarse-grained model (1.2)-(1.5) to give the same value of V 2 i as obtained by the corresponding MD model, the velocity autocorrelation function of the coarse-grained model agrees with MD simulations for t = 0. In figure 6, we observe that the coarse-grained model approximates C(t) well for t between 0 and 30 fs. It is also a good approximation for large values of t, because C(t) → 0 as t → ∞, but we can clearly see a difference for intermediate values of t. Since we have parametrized the coarse-grained model (1.2)-(1.5) by fitting the diffusion constant D, the coarse-grained model also approximates well the integral of C(t), because of the Green-Kubo formula If a modeller wants to improve the approximation of C(t) for intermediate values of t, the fictitious particle model (7.1)-(7.4) can be used for N > 1 as a generalization of the coarse-grained model (1.2)-(1.5). We will illustrate this by using N = 3. Equations (7.5)-(7.7) can then be rewritten in vector notation as follows (cf. (5.1)): Let us denote the eigenvalues and eigenvectors ofB asλ j andν j , j = 1, 2, . . . , 7, respectively. In what follows, we will assume that coefficients α j,k are chosen so that the eigenvalues ofB are distinct and have negative real parts. Then the velocity autocorrelation function of the fictitious particle model (7.10) can be computed as follows (cf. (7.9)) where matrixΦ(t) ∈ R 7×7 is given asΦ(t) = (exp(λ 1 t)ν 1 | exp(λ 2 t)ν 2 | · · · | exp(λ 7 t)ν 7 ), i.e. each column is a solution of the ODE system dỹ i =Bỹ i dt. Vector V iỹi is the first column of the equilibrium covariance matrix ỹ iỹ T i which can be obtained as the solution of the linear system Since we want the fictitious particle model (7.1)-(7.4) to be a generalization of the coarse-grained models (1.2)-(1.5), we select its parameters so that the model (7.1)-(7.4) has the same diffusion constant D and moments V 2 i and U 2 i as calculated by all-atom MD simulations in table 2. This yields three conditions between 12 parameters α j,k , j = 1, 2, 3, k = 1, 2, 3, 4, of the fictitious particle model (7.1)-(7.4). In particular, we have a lot of freedom to choose the 12 model parameters. In figure 6, we plot the velocity autocorrelation function (7.11) for specific choices of parameters α j,k , illustrating that we can select these parameters to approximate velocity autocorrelation functions obtained by all-atom MD. To obtain these results, we use a simple acceptance-rejection algorithm based on equations (7.11)-(7.12). We start with an initial guess of 12 parameters α j,k and calculate the L 1 -error between the velocity autocorrelation function (7.11) and the MD result (blue solid line in figure 6). Then we perturb the parameters α j,k in a way that the resulting model still has the same D, V 2 i and U 2 i and we use the new set of parameters α j,k to recalculate the L 1 -error between the velocity autocorrelation function (7.11) and the MD result. If the error decreases, then we accept the new set of parameters. We iterate these steps until the desired level of error.
In figure 6, we observe that the velocity autocorrelation function computed by (7.11) for N = 3 (red dashed line) approximates well the MD result. In a similar way, we could fit other autocorrelation functions (if required) by coarse-grained models in the hierarchy (7.1)-(7.4) for sufficiently large N. However, as we have seen in §6, coupling all-atom MD and BD can be achieved by the coarse-grained model (1.2)-(1.5), i.e. by using N = 1.

Discussion
In this paper, we have introduced and studied the coarse-grained model (1.2)-(1.5) of an ion in aqueous solution. We have parametrized this model using all-atom MD simulations for four ions (K + , Na + , Ca 2+ and Cl − ) and showed that this model provides an intermediate description between all-atom MD and BD simulations. It can be used both with MD timestep t (to couple it with all-atom MD simulations) and BD time step T (to couple it with BD description (1.1)). In particular, the coarse-grained model enables multiscale simulations which use all-atom MD and BD in different parts of the computational domain.
In §6, we have illustrated this multiscale methodology using a first passage type problem where we have reported the time taken by an ion to reach a specific distance. Possible applications of this multiscale methodogy include problems where a modeller considers all-atom MD in several different parts of the cell (e.g. close to binding sites or ion channels) and wants to use efficient BD simulations to transport ions by diffusion between regions where MD is used. The proposed approach thus enables the inclusion of MD-level of detail in computational domains which are much larger than would be possible to study by direct MD simulations.
Although the illustrative simulations in §6 are reported over distances of the order of 10 2 Å, this is not a restriction of the method. Most of the computational time is spent by considering all-atom MD in Ω 1 ∪ Ω 2 . BD uses a much larger timestep which enables us to further extend the BD region Ω 5 (and consequently, the original domain Ω). Moreover, if we are far away from the MD domain Ω 1 , we can further increase the efficiency of BD simulations by using different BD timesteps in different parts of the BD subdomain Ω 5 [12], or by using event-based BD algorithms [35,37]. The computational intensity of BD simulations can be further decreased by using multiscale methods which efficiently and accurately combine BD models with lattice-based (compartment-based) models [43,44]. Such a strategy has been previously used for modelling intracellular calcium dynamics [3,7] or actin dynamics in filopodia [45], and enables us to extend both temporal and spatial extents of the simulation.
In §7, we have presented a systematic procedure for including additional auxiliary variables and parameters to the coarse-grained model (1.2)-(1.5). We have illustrated that the generalized models can be used to fit additional details of all-atom MD simulations. The generalization of the coarse-grained model (1.2)-(1.5) is based on the fictitious particle approach [13] which introduces fictitious particles with harmonic interactions with coarse-grained degrees of freedom. Such special heat baths have been previously studied in the context of the generalized Langevin equation [14][15][16][17]. In §7, we have shown that an appropriately formulated fictitious particle model which uses one fictitious particle per ion [13] has the same dynamics as the coarse-grained model (1.2)-(1.5).
In the literature, MD methods have been used to estimate parameters of BD simulations of ions [46]. There has also been a lot of progress in systematic coarse-graining of MD simulations [40,47]. The approach presented in this paper not only uses all-atom MD simulations to estimate parameters of a coarser description, but it also designs a multiscale approach where both methods are used during the same simulation. Methods which adaptively change the resolution of MD on demand have been previously reported in [48,49]. They include algorithms which couple all-atom MD with coarse-grained MD. The coarse-grained model studied in this work does not include any water molecules and has different application areas. One of them is modelling of calcium induced calcium release through IP 3 R channels [3] which is discussed as a motivating example in §1. MD simulations in this paper use the three-site SPC/E model of water. An open question is to extend our observations and analysis to other MD models of water, which include both more detailed water models with additional sites [21,22] and coarse-grained MD models of water [50].