## Abstract

Molecular dynamics (MD) simulations of ions (K^{+}, Na^{+}, Ca^{2+} 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.

## 1. 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
**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
**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, *d***W**≡(*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.2)–(1.5) and show how parameters *η*_{j}, *j*=1,2,3,4, can be connected with diffusion constant *D* used in the BD model (1.1).

The coarse-grained model (1.2)–(1.5) is used as an intermediate model between the all-atom MD model and BD description (1.1). In §5, we show how it can be coupled with the BD model which uses a much larger timestep than the MD model. In §6, the coarse-grained model (1.2)–(1.5) is coupled with all-atom MD simulations. We then show that all-atom MD models of ions can be coupled with BD description (1.1) using the intermediate coarse-grained model (1.2)–(1.5). In §7, we introduce a hierarchy of stochastic coarse-grained models which generalize the 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–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.

## 2. 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 [23,1]. 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],
*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 long-range (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.

## 3. 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**≡(*V* _{1},*V* _{2},*V* _{3}) which evolve according to
*M* is the mass of the ion (given in table 1) and **F**≡(*F*_{1},*F*_{2},*F*_{3}) is the force acting on the ion. We use all-atom MD simulations as described in §2 to estimate diffusion coefficient *D* and second moments of *V* _{i} and *U*_{i}=*F*_{i}/*M*, *i*=1,2,3. They are given for four ions considered in table 2. To estimate *i*th direction *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:
*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
*U*_{i}*V* _{i}〉=0 and 〈*U*_{i}*Z*_{i}〉=0, we obtain that 〈*V* _{i}*Z*_{i}〉=0 at steady state and
*η*_{1} using the MD averages *η*_{1}, we can also estimate the value of *Z*_{i} and equation (1.5) by *U*_{i}, we obtain
*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:
*X*_{i}*V* _{i}〉=*D*, *X*_{i}*Z*_{i}〉=*η*_{1}*D* and
*Z*_{i}, we obtain

## 4. 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 *D* which is the sole parameter of the BD model (1.1). In this section, we explain why the coarse-grained 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–31]:
**X**≡(*X*_{1},*X*_{2},*X*_{3}) is the position of a diffusing molecule, **V**≡(*V* _{1},*V* _{2},*V* _{3}) 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 *t*. To illustrate this, let us parametrize Langevin dynamics (4.1)–(4.2) using diffusion constant *D* and the second velocity moment *ξ*_{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
*D* and ^{+} 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, ^{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 *J*(*v*,*u*) from all-atom MD simulations, we calculate the rate of change of acceleration during each MD timestep
*U*_{i}(*t*+Δ*t*)−*U*_{i}(*t*))/Δ*t* during every timestep and record their average in two-variable array *J*(*v*,*u*) indexed by binned values of *V* _{i}(*t*)=*v* and *U*_{i}(*t*)=*u*. 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
*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.3)–(1.5), i.e. *p*(*v*,*u*,*z*) *dv* *du* *dz* is the probability that *V* _{i}(*t*)∈[*v*,*v*+*dv*), *U*_{i}(*t*)∈[*u*,*u*+*du*) and *Z*_{i}(*t*)∈[*z*,*z*+*dz*). Then the jerk matrix (4.6) of the coarse-grained model (1.2)–(1.5) is

Using (1.4), we rewrite it as
*p*(*v*,*u*,*z*) of (1.3)–(1.5) is Gaussian with mean (0,0,0) and stationary covariance matrix:
*t*=10^{−3} ps) for both the coarse-grained model (1.2)–(1.5) and all-atom MD simulations. The coarse-grained model (1.2)–(1.5) can also be coupled with BD description (1.1), which uses much larger timesteps, as we show in the next section.

## 5. 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 **y**_{i}=(*V* _{i},*U*_{i},*Z*_{i})^{T}, where T stands for transpose, equations (1.3)–(1.5) can be written in vector notation as follows:

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
*η*_{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 table 3, we present the values of eigenvalues *λ*_{j}, *j*=1,2,3, in table 4. The eigenvalues *λ*_{j}, *j*=1,2,3, are distinct. The general solution of the SDE system (5.1) can be written as follows [32]:
*d***y**_{i}=*B***y**_{i} *dt*. Considering deterministic initial conditions, equation (5.3) implies that the process is Gaussian at any time *t*>0. Equations for means, variances and covariances then uniquely determine the distribution of **y**_{i}(*t*) for *t*>0. Equations for means can be written in the vector form as *d*〈**y**_{i}〉=*B*〈**y**_{i}〉 *dt*. Equations for variances and covariances are given in §3 as equations (3.3)–(3.6), (3.9) and (3.16).

There are two important conclusions of the above analysis. First of all, eigenvalues *λ*_{j}, *j*=1,2,3, given in table 4 satisfy
*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
*a*. We solve the system of 10 ODEs for variances and covariances given as equations (3.3)–(3.6), (3.9), (3.11)–(3.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 *Dt*. We observe that there is an approximately constant shift, denoted *t*>0.2 ps. We illustrate this further by plotting *a*. The values of shift

Next, we show how the BD model (1.1) and the coarse-grained model (1.2)–(1.5) can be used in different parts of the computational domain. This coupling will form one component of multiscale methodology developed in §6. BD algorithms based on equation (1.1) have been implemented in a number of methods designed for spatio-temporal modelling of intracellular processes, including Smoldyn [33], MCell [34] and Green's-function reaction dynamics [35]. Smoldyn discretizes (1.1) using a fixed BD timestep Δ*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. 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 *Ω*_{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 2*b*, we present an illustrative simulation where *h*=1 Å. We report averages over 10^{6} simulations of ions, half of them are initiated at **X**(0)=(*h*,0,0), i.e. they initially follow the coarse-grained model (1.2)–(1.5) with zero initial condition for other variables (*V* _{i}(0)=*U*_{i}(0)=*Z*_{i}(0)=0). The second half of ions are initiated at **X**(0)=(−*h*,0,0), i.e. they initially follow BD description (5.5). We plot the (marginal) distribution of ions along the first coordinate (*X*_{1}) at time 10^{3} ps in figure 2*b*. The computed histogram is plotted using bins of length 2 Å, 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
*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 *X*_{i}(0)=0), we can again compute the mean square displacement. As in figure 2*a*, it can be shifted in time to better match with the BD result, 2*Dt*. We denote this time shift as *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 2*b*. Thus, we can use the above approach in *Ω*_{3}∪*Ω*_{4}∪*Ω*_{5} without introducing observable errors in the multiscale method developed in the next section.

## 6. Coupling all-atom MD and BD

Let us consider all-atom MD in domain *Ω*_{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
*Ω*_{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 coarse-grained model (1.2)–(1.5). We still simulate water molecules in *Ω*_{1} and we allow them to 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:
*ω*≥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
*Ω*_{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 ion at its real position **X** in *Ω*_{1}, but also by its copies at periodic locations **X**+(*iL*,*jL*,*kL*), where *Ω*_{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 coarse-grained model (1.2)–(1.5). In figure 4*a*, 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 4*b* for three different values of *ω*. The results compare well with (2*Dt*)^{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 *Ω*_{1} and *Ω*_{2} are given by (6.3), and
*ω*=10, *h*_{1}=*L*/20 and *h*_{2}=*L*/10. Then the BD domain is *Ω*_{1}), i.e. **X**(0)=(0,0,0), and we simulate each trajectory until it reaches the distance 4*L*=99.32 Åfrom the origin. Let *r* from the origin. In figure 5, we plot escape time *r*. We plot the value of *r*=4*L*) 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 plot average *Ω*_{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}.

## 7. 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–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 *j*th fictitious particle, *N*+1) equations (for *i*=1,2,3 and **W**_{j}≡(d*W*_{j,1},d*W*_{j,2},d*W*_{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 *j*th 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 *N*=1 and denoting *η*_{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
^{+} 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 *λ*_{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 *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*. 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
*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)):
*α*_{j,k} are chosen so that the eigenvalues of *D* and moments *α*_{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*, *α*_{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.

## 8. 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–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].

## Competing interests

I have no competing interests.

## Funding

I thank the Royal Society for a University Research Fellowship and the Leverhulme Trust for a Philip Leverhulme Prize.

- Received August 12, 2015.
- Accepted January 11, 2016.

- © 2016 The Authors.

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