Multi-resolution dimer models in heat baths with short-range and long-range interactions

1. Introduction

Molecular dynamics (MD) approaches, based on the rules of classical mechanics, have been used to study the behaviour of complex biomolecules in biological applications [1,2]. They are written in terms of the positions and velocities of particles, representing either individual atoms or groups of atoms, describing parts of a biomolecule [3–6]. Inter-particle forces in MD models include combinations of short-range and long-range interactions [7,8]. In all-atom MD models, a common example of short-range forces are interactions described by the Lennard–Jones potential [9,10], while Coulomb forces provide an example of long-range forces [7]. Considering coarse-grained or caricature MD models, short-range interaction models include systems when particles only interact through direct collisions [11–14], while long-range interactions also include models where particles interact through harmonic springs [15,16]. Once the inter-particle interactions are specified, MD describes the time evolution of the model as a system of ordinary or stochastic differential equations for the positions of particles, which can also be subject to algebraic constraints, representing bonds between atoms or fixed internal structures of a biomolecule [2,17,18].

Biologically relevant simulations have to be done in aqueous solutions. A number of water models have been developed in the literature to use in all-atom MD simulations, including commonly used three-site (SPC/E, TIP3P) models [19,20]. In coarse-grained MD models, water is often treated with the same level of coarse-graining as other molecules in the system. For example, four water molecules are combined into a single coarse-grained water bead in the Martini model [3], while the Wat Four water model [6] uses four linked beads placed at the corners of a tetrahedron to collectively represent 11 water molecules. In this paper, we consider two theoretical heat baths which enable more analytical progress than solvent models based on all-atom or coarse-grained water models. In both cases, the convergence to the Langevin description of the solute particle can be established in a certain limit [11–16]. Our solute particle will also be treated with the same level of simplicity and described as a simple dimer molecule consisting of two monomers (beads) connected by a spring.

Multi-resolution (hybrid) methods use detailed and coarse-grained simulations in different parts of the simulation domain during the same dynamic simulation [21–24]. Such methods have been developed in different application areas and at different spatial and temporal scales in the literature, including dual-resolution approaches AdResS and H-AdResS for all-atom MD simulations [25–29], methods for coupling Brownian dynamics approaches with lattice-based stochastic reaction-diffusion models [30–32] or methods which make use of continuum mean-field equations for the macroscopic component of the simulation [33–35].

In some multi-resolution MD approaches, the region of high resolution moves together with the large microscopic structure of interest so that the high resolution model is always used for the whole considered structure, which can range in size from a single biomolecule (a protein or a DNA in solution [27,28]) to virus-like particles [36,37]. The structure of interest is placed in the centre of the simulation domain, and it is solvated using a detailed atomistic MD water in its immediate neighbourhood, which is coupled with a coarse-grained water description in the rest of the computational domain.

Another type of multi-resolution modelling is used for modelling of macromolecules where a detailed model of an important part of a macromolecule is coupled with a coarser model of the rest of the macromolecule. For example, atomistic detail of the active part of an enzyme has been coupled with a coarser model of the rest of the protein [38], different resolutions have been used in bead-spring modelling of DNA [39,40] or for modelling of polymer melts [41,42].

In this paper, we study both multi-resolution approaches using a simple dimer model consisting of two monomers (beads) connected by a spring. Similar models, where a macromolecule is described as several beads, representing parts of the simulated biomolecule, connected by springs, have been obtained in the literature using the method of ultra-coarse-graining [43]. Thus, our dimer model can be considered as a caricature of an ultra-coarse-grained model of a macromolecule. We study its behaviour in two theoretical heat baths. Our investigation focuses on multi-resolution (multiscale) descriptions of the solvent which can be described at the microscopic level of individual solvent molecules or at the macroscopic (dimer) level with the introduction of extrinsic random thermal forces on the monomers. We present models of the same dimer with various multi-resolution descriptions for the solvent and highlight the conditions and reasons, when and why, different model approximations of the solvent may be made in simulations.

Our paper is organized as follows. In §2, we introduce the macroscopic dimer model with a macroscopic description for solvent forces. This macroscopic model is fully described by Langevin equations. The Langevin macroscopic model is commonly used in simulation due to ease of implementation and analysis. We discuss in §2 the properties of this description with the intent to use these properties as benchmarks against which to compare microscopic and multi-resolution solvent models for the same dimer. Two theoretical microscopic approaches to model the solvent are introduced and studied through multi-resolution (simultaneous microscopic and macroscopic coupled) modelling in §§3 and 4. One of them is based on (very) short-range interactions, as heat bath particles only interact with the dimer on contact. The other one is at the opposite extreme, as it is based on (very) long-range interactions, where the heat bath is modelled as a system of many harmonic oscillators.

2. The dimer model

In this section, we will talk exclusively about the construction of the model for the dimer which will be used throughout this manuscript. In doing so, we describe the solvent at the macroscopic level as an extrinsically added random force. The result will be a set of Langevin equations. Throughout the manuscript, we will modify the treatment of the solvent forces at various scales and hybrid resolutions but the underlying dimer model will be the same.

We consider a model of a dimer which is described by positions of its two monomers, denoted by X₁ = [X_1;1, X_1;2, X_1;3] and X₂ = [X_2;1, X_2;2, X_2;3], respectively. Each monomer has the same mass, M. We denote by R the vector describing the separation between the monomers, i.e. R = X₂ − X₁, and by R its magnitude R = |R|. The interaction between monomers is given in terms of the potential $\mathit{\Phi }\equiv \mathit{\Phi }(R):[0,\mathrm{\infty })\to \mathbb{R}$, which generates a force on each of the monomers with magnitude Φ′(R).

When the dimer is placed into a heat bath, there are additional forces on the two monomers caused by interactions with solvent molecules. The solvent forces can be modelled in a number of different ways and at various scales. In this manuscript, we consider two classes of models to describe the solvent–dimer interactions. The first, presented in §3, models the solvent as a bath of point particles which collide with the monomers and elastic collisions (short-range interactions) contribute to the generation of the forces. In the second case, described in §4, solvent molecules are point particles which oscillate around and interact at a distance (through long-range interactions) with the monomers. The solvent-dimer interactions are the sum of harmonic oscillatory forces acting on each of the monomers. Importantly, both descriptions under suitable assumptions lead to a macroscopic description of the dimer given by the following set of Langevin equations:

where V₁ = [V_1;1, V_1;2, V_1;3] and V₂ = [V_2;1, V_2;2, V_2;3] are velocities of the first and second monomer, respectively, W₁ and W₂ are three-dimensional vectors of independent Wiener processes, D is a diffusion coefficient and γ is a friction coefficient, with dimension ${[\gamma ]=[\text{time}]}^{-1}$.

System (2.1)–(2.4) provides a macroscopic model of the dimer, which we compare with microscopic (or multi-resolution) MD simulations which explicitly model the solvent. Its validity for different MD models can be tested by comparing values of different dimer’s statistics at equilibrium, including its expected length L_d, dimer velocity autocorrelation function C_d(τ) and dimer diffusion constant D_d, defined by

where $\overline{\mathbf{X}}=({\mathbf{X}}_1+{\mathbf{X}}_2)/2$ is the centre of mass of the dimer and $\overline{\mathbf{V}}=({\mathbf{V}}_1+{\mathbf{V}}_2)/2$ is its velocity. These quantities can be obtained analytically for our macroscopic model (2.1)–(2.4) as follows. Adding equation (2.2) and equation (2.4) and noting that the sum of two independent Wiener processes is another Wiener process W with its variance equal to the sum of the variances of the original two processes, we obtain an Ornstein–Uhlenbeck process for $\overline{\mathbf{V}}$ in the following form: Therefore, we have Integrating over τ, we deduce Taking the difference of equation (2.4) minus equation (2.2), implementing the over-damped assumption (where γ is large) and combining the independent Weiner processes into a single Weiner process W gives The stationary distribution corresponding to this process is proportional to exp [−Φ(R)/(MDγ)]. Normalizing, we find the distribution of dimer lengths equal to In the simulations that follow in this manuscript, we shall be assuming the dimer potential acts like a linear spring with a rest length of ℓ₀ and a spring constant of k between the two monomers. That is, we shall assume Each monomer within the dimer is representing a half of a molecule of interest and the value of the spring constant indicates the flexibility in which the molecule can change its shape. In this paper, we consider the parameter regime where the spring constant k is sufficiently large so that the dimer has a well-defined structure. In the limit of large k, we have $\epsilon =\hspace{0.17em}MD\gamma /(k\hspace{0.17em}{\ell }_0^2)\ll 1$. Then, L_d can be calculated as which is valid up to the first order in ɛ. In particular, the presence of heat baths extends the dimer from its rest length on average. In the following two sections, we study two theoretical MD models, where we use equations (2.6), (2.7) and (2.9) to compare the macroscopic theory with the results obtained by MD simulations.

3. Short-range interaction heat bath

We describe the two monomers as balls with radius r₀ and mass M which interact with point solvent particles when they collide with them. In particular, this is a theoretical model of a (very) short-range interaction heat bath. Between collisions, monomers follow Newton’s Second Law of Motion in the form

and where, following our notation introduced in §2, positions and velocities of the monomers are denoted by X_i and V_i, respectively, and R = X₂ − X₁.

Our short-range interaction heat bath is described in terms of positions ${\mathbf{x}}_i^{\hspace{0.17em}j}$ and velocities ${\mathbf{v}}_i^{\hspace{0.17em}j}$, of heat bath particles, where i = 1, 2 is the monomer number and j = 1, 2, 3, …, is the number of the heat bath particle. Notice that this formulation allows us to consider two important cases: (a) each monomer has its own heat bath; (b) a single heat bath is shared by both monomers. By comparing our simulation results in cases (a) and (b), we can explicitly investigate whether there are any significant hydrodynamic interactions between the monomers. In the case (b), we simplify our notation by describing particles of the single heat bath by

In both cases (a) and (b), we assume that all heat bath particles have the same mass, m, and define (dimensionless) parameter μ by We are interested in the parameter regime where μ ≫ 1. Our MD model is based on elastic collisions of heavy monomers (balls with mass M and radius r₀) with point heat bath particles with masses m. We assume that the collisions are without friction, then conservation of momentum and energy yields the following formulae for post-collision velocities [12]: and where ${\mathbf{v}}_i^{\hspace{0.17em}j}$ is the velocity of the heat bath particle which collided with the ith monomer, tildes denote post-collision velocities, superscripts ⊥ denote projections of velocities on the line through the centre of the monomer and the collision point on its surface, and superscripts $\parallel$ denote tangential components.

Heat bath models based on elastic collisions (3.4) and (3.5) have been studied by a number of authors [11–14]. Consider a single monomer in infinite domain ${\mathbb{R}}^3$, and let the heat bath consist of an infinite number of particles with positions distributed according to the spatial Poisson process with density

This means that the number of points in a subset Ω of ${\mathbb{R}}^3$ has its probability mass function given by the Poisson distribution with mean ${\lambda }_{\mu }|\mathit{\Omega }|$, where |Ω| is the volume of Ω. Let the velocities of the heat bath particles be distributed according to the Maxwell–Boltzmann distribution where v = [v₁, v₂, v₃] and Then the monomer’s behaviour is known to converge to the Langevin dynamics [12,14]. In particular, if we consider that each monomer has its own heat bath, we can show that the position and velocity of the monomers, X_i and V_i, converge (in the sense of distributions) to the solution of (2.1)–(2.4) in the limit μ → ∞.

In reality, all beads representing a macromolecule exist within a single heat bath. Thus, we ask whether the correlations introduced by a bath of solvent which interacts with both monomers has a non-negligible affect on the equilibrium statistics of the dimer. Introducing such coupled heat baths for both short-range (in this section) and long-range (in §4) interactions, we study whether there is a significant difference between the one-bath and two-bath models as we vary ℓ₀, the separation distance, introduced in equation (2.8). In order to study this problem, we make use of multi-resolution modelling.

3.1. Multi-resolution model using a co-moving frame

The solvent in the short-range heat bath interacts with the monomers of the dimer through direct contact. In order to simulate the model for long times, i.e. where the dimer has undergone a large excursion, the simulated domain must be vast as will be the number of solvent particles that must be modelled. We present a multi-resolution approach where we only model the solvent that is within the close vicinity of the dimer. We consider a co-moving cubic frame of length L that is centred at ${\mathbf{X}}_{\mathrm{f}}(t)$, which we here identify with the centre of mass of the dimer at time t, i.e.

Within this frame, we explicitly model the heat bath with solvent particles, i.e. they are simulated in the cubic box Externally, we model the heat bath as a continuum, where the particles are distributed according to the spatial Poisson process with density λ_μ given in (3.6) and the velocities are distributed according to f_μ(v) given in (3.7), see figure 1a for a diagrammatic representation of the multi-resolution framework (drawn for clarity in two spatial dimensions, while all our simulations are three-dimensional). As the dimer moves around in ${\mathbb{R}}^3$, the frame will move with it. In order for the multi-resolution model to capture, the full model where solvent particles are distributed in the entire domain, ${\mathbb{R}}^3$, we need to introduce new solvent particles at the boundary of the frame.

• Download figure
• Open in new tab
• Download PowerPoint

Consider that time is discretized using small time step Δt, i.e. if the current time is t, we want to calculate the state of the system at time t + Δt. In our simulations of the multi-resolution model, we need the probability of introducing a particle at a boundary of the frame (3.10) in a timestep of length Δt and subsequently the distribution of the position x_new and velocity v_new of the new solvent particle. For simplicity, we transform into the coordinate system of the co-moving frame which over an interval of length Δt has velocity

The frame is always translated to occupy the region [0, L]³. Thus, the velocities for the solvent particles in the new reference frame are given by ${\mathbf{w}}^{\hspace{0.17em}j}={\mathbf{v}}^{\hspace{0.17em}j}-{\mathbf{V}}_{\mathrm{f}}$. We first calculate the density of particles that enter the frame via a particular boundary within a timestep of length Δt. Take, as an illustrative example, the boundary face corresponding to {x₁ = 0}. Consider particles which are in half-space $(-\mathrm{\infty },0)\times {\mathbb{R}}^2$ at time t. These particles have not yet been explicitly included in the simulation. Some of them will be in half-space $(0,\mathrm{\infty })\times {\mathbb{R}}^2$ at time t + Δt. Their density, h(x₁), only depends on their first coordinate x₁ ∈ (0, ∞). We can calculate h(x₁) by integrating density (3.6) and (3.7) over solvent particles which are at $x_1^{\mathrm{\prime }}\in (-\mathrm{\infty },0)$ at time t and have the appropriate velocity to reach x₁ ∈ (0, ∞) at time t + Δt, namely as [14] where $V_{\mathrm{f};1}$ is the first component of the frame velocity and $\text{erfc}(z)=2/\sqrt{\pi }{\int }_z^{\mathrm{\infty }}\exp (-s^2)\hspace{0.17em}\text{d}s$ is the complementary error function. Integrating (3.12) over the domain (0, ∞) × [0, L] × [0, L] gives us the average number of particles that have entered the frame from the {x₁ = 0} boundary in a time interval of length Δt as In our simulations, we choose a timestep small enough that $p_{\mathrm{in}}\ll 1$, we can therefore use $p_{\mathrm{in}}$ as the probability of introducing a new solvent particle. Let z = [z₁; z₂; z₃] be the position of the new solvent particle in the coordinate system of the co-moving frame. Then coordinates z₂ and z₃ are uniformly distributed in (0, L) and the first coordinate can be sampled from the error function distribution where C₁ is a normalizing constant. Then the position of the new solvent particle in the original coordinates is ${\mathbf{x}}_{\mathrm{new}}=\mathbf{z}+{\mathbf{X}}_{\mathrm{f}}(t+\mathrm{\Delta }t)-[L/2,L/2,L/2]$. The velocity, w, of the new particle in the co-moving frame must have a first coordinate exceeding z₁/Δt in order to reach z₁ in a time interval of length Δt. Noting that $\mathbf{w}={\mathbf{v}}_{\mathrm{new}}-{\mathbf{V}}_{\mathrm{f}}$ we write down the distribution of the velocity as the following truncated Gaussian distribution where C₂ is a normalizing constant and $H(\cdot )$ is the Heaviside step function, satisfying H(y) = 1 for y ∈ [0, ∞) and H(y) = 0 otherwise. The position and velocity of solvent particles introduced at the other five faces can be done by symmetric modifications of the above distributions.

Random numbers from distributions (3.14) and (3.15) can be efficiently sampled using acceptance–rejection algorithms. We use an acceptance–rejection method for the truncated normal distribution (3.15) presented in the literature [44], while we sample random numbers from the distribution (3.14) using the acceptance–rejection algorithm presented in table 1. This is a generalization of the acceptance–rejection algorithm for sampling random numbers according to the distribution $\sqrt{\pi }\hspace{0.17em}\text{erfc}(z)$ previously used in simulations in the stationary frame [14]. In the case of the distribution (3.14), we need to sample random numbers according to the probability distribution

where $\beta \in \mathbb{R}$ is a constant and C₃(β) is the normalizing constant given by The algorithm in table 1 does this by generating an exponentially distributed random number η₃ with mean a₁(β), where To maximize the acceptance probability of this algorithm, we choose its second parameter, a₂(β), as Then its acceptance probability is depending on β as We plot the acceptance probability (3.20) in figure 2 for our choices (3.18) and (3.19) of a₁(β) and a₂(β) as the solid line.

• Download figure
• Open in new tab
• Download PowerPoint

We observe that the acceptance probability (3.20) for β = 0 is equal to $2/\pi \approx 63.7\text{\%}$. This value can be improved [14] in the case of β = 0 to 86.3% provided that we choose a₁ = 0.532 and a₂ = 0.814. To obtain a similar improvement for all values of β, we could choose both a₁(β) and a₂(β) to maximize the acceptance probability (3.20), rather than postulating that a₁(β) is given by the piecewise defined function (3.18) and optimizing a₂(β) only. The acceptance probability (3.20) of the resulting algorithm (which would have a₁(β) and a₂(β) given by a lookup table, rather than by using formulae (3.18) and (3.19)) is plotted in figure 2 as the dashed line for comparison. However, in our illustrative simulations, we use the acceptance–rejection algorithm in table 1 with the values of a₁(β) and a₂(β) given by (3.18) and (3.19).

Comparing equations (3.16) and (3.14), we observe that we can sample random numbers from the distribution (3.14) by sampling random numbers from the distribution $p(z;V_{\mathrm{f};1}\mathrm{\Delta }t)$ (using the acceptance-rejection algorithm in table 1 for $\beta =V_{\mathrm{f};1}\mathrm{\Delta }t)$) and multiplying them by the factor ${\sigma }_{\mu }\mathrm{\Delta }t\sqrt{2}$.

One iteration (i.e. an update of the state of the system from time t to time t + Δt) of the multi-resolution simulation algorithm in the co-moving frame is given as algorithm [S1]–[S7] in table 2. It evolves the positions and velocities of both monomers together with the positions and velocities of N(t) solvent particles, where N(t) does depend on time t. To formulate algorithm [S1]–[S7], we assume that the timestep Δt is chosen small enough so that at most one collision happens per iteration.

We initialize the two monomers with a separation distance ℓ₀ and generate a Poisson number (with mean λ_μL³) of solvent particles in our simulation domain, the cubic frame (3.10). The solvent particles are initially placed uniformly in the frame (3.10), where we remove particles overlapping with monomers (before we begin our simulation) to get the initial number, N(0), of simulated solvent particles. Their initial velocities are drawn from the Maxwell–Boltzmann distribution (3.7).

In step [S1], we update the system over the time interval (t, t + Δt] using the ‘free-flight’ positions for each monomer and solvent particle, namely we use

and where i = 1, 2 is the monomer number and $j=1,2,\hspace{0.17em}\dots ,N(t)$, is the number of the heat bath particle. Since Δt is chosen so small that only one collision happens during the time interval [t, t + Δt), most of the ‘free-flight’ positions of solvent particles are accepted in step [S2] as their updated positions ${\mathbf{x}}_i^{\hspace{0.17em}j}(t+\mathrm{\Delta }t)$ and only the solvent particle colliding with a monomer is further updated.

In step [S3], we update the velocities of the monomers by solving (3.1) and (3.2) over one time step [t, t + Δt]. We discretize (3.1) and (3.2) using the forward Euler method as follows:

and where ${\tilde{\mathbf{V}}}_i$, for i = 1, 2, is either the post-collision velocity (if a collision happened in step [S2]) or is equal to V_i(t). In steps [S4]–[S5], we update the position and velocity of the frame. We remove solvent particles which are outside of the simulation domain and update N(t) accordingly.

In step [S6], we use probability $p_{\mathrm{in}}$, given by (3.13), to check whether any solvent particle entered the simulation domain during the time interval (t, t + Δt]. Since $p_{\mathrm{in}}$ is the probability of entering the domain through one of its six sides, we can, for time step Δt small enough that $6p_{\mathrm{in}}\ll 1$, introduce at most one solvent particle through a randomly chosen side with probability $6p_{\mathrm{in}}$. The initial position and velocity of the introduced solvent particle are sampled according to distributions (3.14) and (3.15) or their symmetric modifications, taking into account through which side of the cubic frame (3.10) the particle entered the frame.

There is one little caveat in our derivation of $p_{\mathrm{in}}$. To derive equation (3.13), we integrated over the half-space $(-\mathrm{\infty },0)\times {\mathbb{R}}^2$, meaning that once we consider all six faces of the cubic frame (3.10) we have over-counted twice at the edges and three times at the corners (as it is highlighted with darker grey shading in our illustrative diagram in figure 1a). This will have negligible effect if we choose L sufficiently large. However, it can bias our simulation for values of L comparable with the monomer size r₀ when Δt is not sufficiently small as boundary effects become more pronounced. To compensate for this effect, we consider the sampled position, ${\mathbf{x}}_{\mathrm{new}}$ and velocity ${\mathbf{v}}_{\mathrm{new}}$ of the new incoming particle at time t + Δt and calculate its previous position at time t by

If y is in the regions which were counted twice or three times in our derivation, we reject the proposed introduction of the new solvent particle with the corresponding probability. Namely, we use the acceptance probability in step [S6] given by where ${\mathcal{Y}}_{\hspace{0.17em}j}\subset {\mathbb{R}}^3$ is the region of the space which consists of points which have exactly j of their coordinates outside of the interval [−L/2, L/2]. For example, in our two-dimensional diagrammatic representation in figure 1a, the lighter grey shading corresponds to region ${\mathcal{Y}}_1$ while the darker grey shading corresponds to region ${\mathcal{Y}}_2$.

In our illustrative simulations, we use algorithm [S1]–[S7] from table 2 together with parameter values r₀ = 0.08, γ = 10, D = 1, μ = 10³, k = 10⁶ and L = 0.72 for the one-bath case. In figure 3, we compare simulation results of the average length of the dimer at equilibrium, L_d, for the one-bath and two-bath models. Since the two-bath case uses uncoupled heat baths, we can further improve the efficiency of our algorithm by centring the co-moving frame corresponding to each heat bath on the corresponding monomer, i.e. we use ${\mathbf{X}}_{\mathrm{f}}(t)={\mathbf{X}}_i(t)$ for the heat bath corresponding to the ith monomer in step [S4] (instead of the centre of mass (3.9)) and choose smaller value of L in the two-bath case, namely L = 0.32. In both one-bath and two-bath models, the solvent particles are distributed according to the spatial Poisson process with density ${\lambda }_{\mu }$ given by (3.6). The velocities are distributed according to the Maxwell–Boltzmann distribution f_μ(v) given by (3.7). We note that in the two-bath case, our model converges to the Langevin dynamics (2.1)–(2.4) as μ → ∞. This allows us to attribute any changes between the one-bath case and the Langevin model to the correlations induced by sharing a heat bath. The asymptotic analytic result obtained for the Langevin model, equation (2.9), is plotted as the black solid line for comparison.

• Download figure
• Open in new tab
• Download PowerPoint

In figure 3, we set the separation distance to be ℓ₀ = αr₀ where α ≥ 2, such that at this distance apart the monomers are not overlapping. The plot shows the two-sided $99\text{\%}$ confidence intervals for (L_d − ℓ₀)/ℓ₀ for α ∈ {2.25, 2.5, 2.75, 3, 3.5, 4, 4.5, 5}. Firstly, we note that L_d > ℓ₀ in each of the models as predicted in (2.9). There seems to be reasonable correspondence between the one- and two-bath models, with the confidence intervals overlapping. This suggests that the correlations we lose by approximating a larger co-moving frame around both monomers with two smaller dedicated frames around each monomer are negligible, allowing us to increase efficiency without biasing our overall results. In the next section, we build on this observation and present a multi-resolution framework which replaces one of the smaller dedicated frames by a coarser model of the heat bath, written in terms of the Langevin dynamics.

3.2. Monomers with different resolution

As the length of a polymer (i.e. numbers of monomers) increases, a model incorporating solvent particles around each of the monomers becomes increasingly computationally expensive. However, a fully coarse-grained Langevin model of a polymer such as the Rouse model [39] can lack the required level of detail. Thus, some multi-resolution approaches for simulating macromolecules only model an important (small) part of a macromolecule using a detailed modelling approach [38–42]. In our case, we can mimic such methodologies by modelling the first monomer with explicit solvent with a heat bath of physical molecules, while the second monomer is modelled using the Langevin equations (2.3) and (2.4). Such a multi-resolution approach is schematically shown in figure 1b. To simulate this model, we use a co-moving frame, given by equation (3.10), which is centred around the first monomer, i.e. ${\mathbf{X}}_{\mathrm{f}}(t)={\mathbf{X}}_1(t)$.

One iteration of the algorithm is presented as algorithm [M1]–[M5] in table 3. To begin, we initialize the particle positions and velocities in the similar way as in the case of algorithm [S1]–[S7], with the only difference that the cubic frame (3.10) is now centred around the first monomer. steps [M1] and [M2] are directly equivalent to steps [S1] and [S2]. In step [M3], we update the position and velocity of the second monomer by

where ξ is sampled from the normal distribution with zero mean and unit variance. That is, we have replaced the heat bath of the second monomer by solving the corresponding Langevin equation (2.1)–(2.4) using the standard Euler–Maruyama integrator. There have been other schemes developed in the literature for discretizing the Langevin equation such as van Gunsteren & Berendsen [45] and the Langevin Impulse integrators, which capture the Langevin dynamics more accurately especially in the presence of forces, such as the spring force between the monomers [46]. Another option would be to consider the BBK integrator [47], which we use in §4.1, where we present a multi-resolution algorithm for the long-range interaction heat bath model and discretize the Langevin equation using a combination of the velocity Verlet and Euler–Maruyama integrators, see equations (4.11)–(4.15). An additional approach is the Verlet scheme [48] that approximates the velocity using a central difference discretization rather than the forward difference approach used in the Euler–Maruyama method, or Runge–Kutta methods [49], which could further reduce the error of the multi-resolution simulations.

In order to compare simulations of the multi-resolution model with simulations of the Langevin model (2.1)–(2.4), we use the velocity autocorrelation function of the dimer, C_d(τ), given by equation (2.5). It has been analytically calculated for the Langevin description in equation (2.6). In figure 4, we present numerical estimates of the velocity autocorrelation function of the multi-resolution model from long time simulation data, using definition (2.5).

• Download figure
• Open in new tab
• Download PowerPoint

Our results compare well with the theoretical result for the Langevin model, though it seems like there is a slightly raised value for C_d(0). Using (2.7), we can estimate the diffusion constant of the dimer D_d by numerically integrating the velocity auto-correlation function in interval [0, 1]. We obtain D_d ≈ 0.529, while its theoretical value for the dimer model is given in equation (2.7) as D/2 = 0.5. Another approach is to fit the exponential function, in the form equation (2.6), to the computational result presented in figure 4. In this way, the values of both D and γ can be estimated simultaneously. We found that D ≈ 1.0714, which is higher than our parameter value D = 1, and γ ≈ 9.6064, which is lower than γ = 10 used in our simulations. This could suggest that the value of ${\lambda }_{\mu }$ is too low or that of σ_μ is too high in our simulations. However, when these quantities are measured during the simulations we do not observe any deviation. This suggests that, rather than our sampling methods, there are small errors introduced by our implementation of the moving frame, or more profound boundary effects introduced by the small size of the frame. A potential problem in the implementation of the co-moving frame, is that solvent particles that leave the frame never return. For a stationary frame, this is valid as the monomer cannot interact with a particle that leaves. However, for a co-moving small frame centred about the monomer, a solvent particle could leave the frame and return at a later time in the simulation. This is not taken into account in the presented algorithms.

4. Long-range interaction heat bath

Coarse-grained models of molecular systems can be written in terms of beads interacting through coarse-grained force fields. Each bead represents a collection of atoms and a coarse-grained potential energy can be constructed from detailed all-atom MD. Such an approach can usually provide a good description of equilibrium properties of molecular systems, but it does not necessarily lead to correct dynamics if the time evolution of the system is solely based on the Hamiltonian dynamics corresponding to the coarse-grained potential energy surface [50]. Dynamical behaviour can be corrected by introducing additional degrees for freedom (fictitious particles) interacting with each coarse-grained bead [50–52]. Fictitious particles can then be subject to suitable friction and noise terms to correct the dynamics.

Considering our dimer molecule model as an example of a coarse-grained molecule, written in terms of two coarse-grained beads (monomers) interacting through coarse-grained potential energy (2.8), then each monomer could be coupled with one or several fictitious particles interacting with the monomer through a suitable harmonic spring term [50,51]. Our long-range interaction heat bath is based on this approach, by assuming that the ith monomer, i = 1, 2, is coupled with N_i harmonic oscillators, in a manner similar to well-known theoretical heat bath models [15,16]. Then equations (3.1) and (3.2), expressing Newton’s Second Law of Motion, include additional terms as follows [15]:

and where ${\mathbf{x}}_i^{\hspace{0.17em}j}$ is the position of the jth solvent particle which interacts with the ith monomer through a harmonic spring with spring constant k_i,j and interaction constants α_i,j, $j=1,\hspace{0.17em}2,\hspace{0.17em}\dots ,N_i$, i = 1, 2. Equations (4.1) and (4.2) are coupled with the evolution equations for solvent particles. We assume that ${\mathbf{v}}_i^{\hspace{0.17em}j}$ is the velocity of the jth solvent particle interacting with the ith monomer. Moreover, we assume that all oscillators have the same mass, m. Using Newton’s Second Law of Motion, we get the following evolution equations for the heat bath oscillators: and for $j=1,2,\hspace{0.17em}\dots ,\hspace{0.17em}{N}_i$ and i = 1, 2. Unlike in some fictitious particle models [50–52], we do not include friction and random forces into equation (4.4) for solvent, because we assume that we explicitly model all solvent particles, i.e. N₁ and N₂ are considered to satisfy N₁ ≫ 1 and N₂ ≫ 1. We are therefore working ‘close’ to the limit N₁ → ∞ and N₂ → ∞, in which we can get the convergence of our long-range interaction heat bath to the Langevin dynamics as discussed below. In practice, it is impossible to include all solvent molecules in simulations and friction and noise terms are still included to control temperature of the simulated system [2,53]. We can solve the solvent equations of motion (4.3) and (4.4) to give [2,54] where ${\mathbf{x}}_i^{\hspace{0.17em}j}(0)$ is the initial position of the jth heat bath particle corresponding to the ith monomer, ${\mathbf{v}}_i^{\hspace{0.17em}j}(0)$ is its initial velocity and ω_i,j = (k_i,j/m)^1/2 is its frequency. Substituting for ${\mathbf{x}}_1^{\hspace{0.17em}j}$ and ${\mathbf{x}}_2^{\hspace{0.17em}j}$ in dimer’s equations of motion (4.1) and (4.2), we obtain the following coupled system of generalized Langevin equations: and where the friction kernel κ_i(τ) and noise term ξ_i ≡ ξ_i(t) = [ξ_i;1, ξ_i;2, ξ_i;3] are given by for i = 1, 2. We assume that initial positions and velocities of solvent oscillators, ${\mathbf{x}}_i^{\hspace{0.17em}j}(0)$ and ${\mathbf{v}}_i^{\hspace{0.17em}j}(0)$, are both independently sampled according to their equilibrium distributions. Then noise autocorrelation function is given by the generalized fluctuation–dissipation theorem where k_B is the Boltzmann constant and T is the absolute temperature. Next, we assume that the frequencies ω_i,j are sampled from a (continuous) exponential distribution with mean $\overline{\omega }$ and we set our interaction constants equal to where γ > 0 is the friction constant used in equations (2.2) and (2.4). Then friction kernel (4.7) becomes Passing to the limit N_i → ∞ allows us to consider the above summation as a continuous integral over the distribution of oscillator frequencies, with both friction kernels κ₁(τ) and κ₂(τ) converging to the same friction kernel [54] Then ${\int }_0^{\mathrm{\infty }}\kappa (\tau )\hspace{0.17em}\text{d}\tau =\gamma .$ Moreover, we can define the limiting friction kernel by which, for our choice of oscillators’ frequencies and interaction terms (4.7), satisfies κ_∞(τ) = 0 for τ > 0 and κ_∞(0) = ∞. Thus, the limiting kernel is a multiple of the Dirac delta function. Therefore, the position and velocity of the monomers, X_i and V_i, converge to the solution of (2.1)–(2.4) in the limit $\overline{\omega }\to \mathrm{\infty }$, provided that each monomer has its own separate heat bath. Moreover, we obtain the Einstein–Smoluchowski relation for the diffusion constant of the monomer as D = k_BT/(γM).

As in §3, we have explained our MD model of the dimer using the case where each monomer has its own heat bath. We now turn our attention to the case when monomers share their heat bath. This has been studied in the case of the short-range interaction MD model in §3.1 with the help of multi-resolution modelling in a co-moving frame, as schematically shown in figure 1a. In the case of long-range interactions, a co-moving frame is less straightforward to implement because we need to take into account that particles outside of the simulated box do exert (long-range) forces on particles in our simulation domain. Some multi-resolution techniques in the literature solve this problem by introducing suitable overlap (bridging, blending) regions [51,55–57], where molecules which are near the simulation domain exert some partial forces on the simulated molecules.

In what follows, we do not truncate the simulated domain, but we consider a different multi-resolution approach in §4.1. Before then we discuss results comparable to figure 3, i.e. we compare simulations with a single heat bath and two heat baths for the case of our long-range interaction MD model. The results are presented in figure 5, where we use the same values of ℓ₀ as in figure 3, expressed as α-multiples of r₀, although our long-range interaction model does not make use of parameter r₀. The value of L_d is for each value of ℓ₀ calculated from a long simulation over 200 dimensionless time units, where the first 100 time units are used to equilibrate the system, while the second half of each simulation is used to compute L_d. To initialize this model, we start with monomers separated by the rest length ℓ₀ and sample oscillators’ frequencies, ω_i,j, according to the exponential distribution with mean $\overline{\omega }$. Their positions and velocities are sampled from the Maxwell–Boltzmann distribution. For the two-bath model, each dimer particle is separately initialized with its own set of oscillators around their respective positions in space.

• Download figure
• Open in new tab
• Download PowerPoint

In figure 5, we observe that in the case of the two-bath model we obtain results which match well with equation (2.9) for our parameter values. These results are also directly comparable with the results obtained for the two-bath case in figure 3. The situation is more complicated in the case of simulations with a single heat bath with N oscillators. Then, using notation (3.3), we can rewrite (4.1) and (4.2) as

for i = 1, 2, where the heat bath evolution equation (4.4) includes terms corresponding to both monomers for $j=1,2,\hspace{0.17em}\dots ,\hspace{0.17em}N$. Our results will then depend how we choose parameters k_i,j and α_i,j. For example, if we choose k_i,j and α_i,j to be the same for both monomeres, i.e. for $j=1,2,\hspace{0.17em}\dots ,\hspace{0.17em}N,$ then the oscillating frequency of the jth heat bath oscillator is ${\omega }_{\hspace{0.17em}j}=\sqrt{2k_{\hspace{0.17em}j}/m}$ and we can subtract the evolution equations for monomers to obtain This equation does not contain any heat bath variables. Using (4.7) to select α_j, i.e. using ${\alpha }_{\hspace{0.17em}j}^2=2\gamma \overline{\omega }/(Nm\pi {\omega }_{\hspace{0.17em}j}^2)=\gamma \hspace{0.17em}\overline{\omega }/(N\pi k_{\hspace{0.17em}j})$, we get Using potential (2.8), we conclude that we effectively obtain a shorter rest length of the spring which gives the following approximation: This result is plotted in figure 5 together with results obtained by illustrative simulations. We use a long-time simulation of length 200 dimensionless time units, with monomers initially placed at separation ℓ₀ and averaging over the second half of the simulation (of length 100 dimensionless time units) to obtain the presented values of dimer’s expected length L_d.

In figure 5, we observe that the average dimer length, L_d, during our single heat bath simulations is smaller than the natural length of the spring, ℓ₀. However, this conclusion is only a consequence of our choice of parameters (4.9). An opposite phenomenon can be observed in simulations for other parameter regimes. For example, if we divide our oscillators into two groups consisting of N₁ and N₂ oscillators, i.e. N = N₁ + N₂, and choose our parameters k_i,j and α_i,j such that

then our ‘one-bath’ case is effectively equal to the two-bath case for which we have the result given in equation (2.9) presented in figure 5. In particular, depending on our choices of k_i,j and α_i,j, the single heat bath case can both increase or decrease the average length of the dimer.

4.1. Multi-resolution modelling of dimer

In figure 1b, we use our dimer example to illustrate a multi-resolution approach which models a part of a molecule using a detailed MD approach, while using a coarser description of the rest of the molecule. Here, in the same manner as carried out for our short-range model in §3.2, we illustrate such a multi-resolution approach using our long-range interaction MD model. We use the Langevin model (2.1)–(2.4) to coarse-grain one of the monomers, while the other monomer is modelled in detail using the MD model with its heat bath described by harmonic oscillators (4.3) and (4.4). As in figure 4, we again calculate numerical estimates for the velocity autocorrelation function, C_d(τ) in equation (2.5), from long time simulations of the dimer after equilibrium has been reached.

The pseudo-code of one iteration our multi-resolution algorithm is presented as algorithm [L1]–[L5] in table 4. algorithm [L1]–[L5] is based on the velocity Verlet integrator, where both monomers are updated by

where A_i, for i = 1, 2, is the acceleration of the corresponding monomer. For the first monomer, its acceleration A₁ is defined as the right-hand side of equation (4.1) divided by M, i.e. For the second monomer, we use the BBK integrator [47], i.e. we define its acceleration as where ξ is sampled from the normal distribution with zero mean and unit variance. The corresponding solvent oscillator integrator is identical to the scheme (4.11)–(4.13), with X₁, V₁ and A₁ replaced by x^j, v^j and a^j, respectively, where acceleration a^j is defined as the right-hand side of equation (4.4) divided by m, i.e. The results obtained by algorithm [L1]–[L5] are compared with analytic results given by equation (2.6) for the Langevin model in figure 6. We see that there is a good correspondence between these, suggesting that the value $\overline{\omega }=100$ is large enough to create an accurate Dirac delta approximation from the kernel function (4.8), along with having a large enough number of oscillators, N₁ = 10⁵, in our heat bath for our other approximations to hold. If these conditions did not hold, we would see that our kernel function has a different form (for example, decaying at a slower rate), and in this case we would have to use a generalized Langevin model as our coarse-graining approach in order to capture the dynamics of the dimer with sufficient accuracy.

• Download figure
• Open in new tab
• Download PowerPoint

The diffusion constant of the dimer, D_d, can again be estimated by numerically integrating the velocity auto-correlation function. Integrating our results from figure 6 over interval [0, 1], we obtain D_d ≈ 0.510, which compares well with the theoretical value, D/2 = 0.5, given by equation (2.7).

5. Discussion and conclusion

In this paper, we have used two theoretical heat baths. Although these heat baths are based on qualitatively different descriptions of solvent–dimer interactions, they both lead to the Langevin description, given in equations (2.1)–(2.4), in a certain limit. In particular, we can use this limiting process to coarse-grain a part of the simulated dimer molecule, while using a detailed MD model to describe the rest of the molecule. Such a multi-resolution approach has potential to significantly speed up computer simulations of dynamics of macromolecules [38–42], provided that it is combined with additional multiscale and multi-resolution methodologies, discussed below.

Our long-range interaction model leads to the system of generalized Langevin equations, given by equations (4.5) and (4.6). Although we have worked in the parameter regime where the generalized Langevin equations can be well approximated by the system of Langevin equations given by (2.1)–(2.4), this will not be the case in other parameter regimes and for more realistic solvent descriptions, especially when the memory kernel is estimated from MD simulations [58,59]. One possible strategy in this case is to couple a detailed MD model with a stochastic coarse-grained model which is written with the help of additional variables [50–52]. To improve the efficiency of simulations further, one can then coarse-grain such a generalized Langevin description using a Brownian dynamics approach [14,60]. Brownian dynamics modelling can be further coupled with stochastic reaction–diffusion modelling based on lattice-based (compartment-based) methods [22]. Lattice-based models are very attractive for simulations of intracellular processes, because they enable modelling of spatio-temporal processes in the whole cell or its significant parts [61]. Coupling Brownian dynamics with compartment-based approaches has been used in a number of applications, including multi-resolution modelling of actin dynamics in filopodia [62,63] or for modelling intracellular calcium dynamics [64].

In this paper, we have investigated multi-resolution approaches, schematically described in figure 1a,b. Another class of multi-resolution approaches in the literature considers a fixed subdomain of the computational domain where a detailed modelling approach is used, which is coupled with a coarser model in the rest of the simulation domain [21,22]. Such an approach is useful, for example, when modelling intracellular ion dynamics. Ions pass through an ion channel in single file and an MD model has to be used to accurately compute the discrete, stochastic, current in the channel [65,66], while the details of the behaviour of individual ions are less important away from the channel where copy numbers may be very large. Thus, we can improve efficiency of our simulations if we allow ions to pass between regions with an explicitly modelled heat bath and a region where their trajectories are described by coarser stochastic models [51].

A similar multi-resolution approach can also be designed for our illustrative dimer model. It is schematically shown in figure 1c, where we identify the region with explicitly simulated heat bath as $\{x_1>b\}=(b,\mathrm{\infty })\times {\mathbb{R}}^2$, where b is the fixed position of the boundary. We are again interested in the behaviour of the dimer in the MD model which would be considered in the full space, ${\mathbb{R}}^3$. However, we now want to replace solvent particles which are in $\{x_1<b\}=(-\mathrm{\infty },b)\times {\mathbb{R}}^2$ by a coarser Langevin description (2.1)–(2.4). To do that, we have to carefully consider how we handle the transfer of monomers between {x₁ > b} and {x₁ < b}. In figure 1c, we present a two-dimensional illustration of a monomer when it intersects the interface, {x₁ = b}. Such a monomer is subject to the collisions with heat bath particles on the part of its surface which lies in {x₁ > b}. This has to be compensated by using a suitable random force from {x₁ < b}, so that the overall model is equivalent to (2.1)–(2.4) in the Langevin limit. Such correction terms can be derived analytically for the case of a spherical monomer in our short-range interaction heat bath and are presented in [14,54]. They can be used to couple the MD model with its corresponding Langevin description, which can be further coupled with Brownian dynamics, simulated using a much larger time step [14].

Mathematical analysis of multi-resolution methodologies can make use of the analysis of the model behaviour close to the boundaries of the computational domain. For example, derivations of reactive (Robin) boundary conditions of macroscopic models from their corresponding microscopic descriptions [67–69] can be generalized to the analysis of behaviour of molecules close to hybrid interfaces in multi-resolution schemes [21,30,31]. Analysis of open boundaries of MD schemes (i.e. boundaries which can transfer mass, momentum and energy) can lead to further understanding of multi-resolution schemes such as AdResS and hybrid continuum-particle dynamics [70], which enable efficient simulation of biomolecules at realistic physiological conditions [71].

Equations for coupled detailed/coarse-grained models can be systematically derived using Zwanzig’s projection method, which has been used to address co-existence of atoms and beads (larger coarse-grained units) in the same dynamic simulations [72,73]. The equations of motion take the form of dissipative particle dynamics, which have been coupled with atomistic water simulations to design multi-resolution schemes in the literature [74]. Other multi-resolution methods couple atomistic water with specially designed coarse-grained water models [75] or with a continuum approach [35]. Coupling discrete and continuum approaches can also be done for different molecular species present in the system and our choice of a modelling approach for each species can be based on its relative abundance [76–78].

One of several important points which have been left out from our discussion is the discretization of time. Although our illustrative simulations use the same time step for both the MD model and the Langevin description, this is not the most efficient or desirable strategy, because the MD model requires much smaller time step than the corresponding Langevin equation. There is potential to design more efficient schemes by updating the coarser description only at certain multiples of the time step which is used in the most detailed model [39]. This is also the case when a modeller further coarse-grains the Langevin description into a Brownian dynamics model which uses even large timesteps [14].

Data accessibility

The computer codes used to compute illustrative results in figures 3–6 have been uploaded as part of the electronic supplementary material. Figures 1 and 2 contain no data.

Authors' contributions

R.G. and D.W. wrote computer codes and performed computer simulations to obtain illustrative results presented in figure 3 (D.W. and R.G.), figure 4 (D.W. and R.G.), figure 5 (R.G.) and figure 6 (R.G.). All authors worked on the analysis of both (short-range and long-range interaction) MD models, wrote the paper and gave final approval for publication.

Competing interests

We declare we have no competing interests.

Funding

This work was supported by funding from the Engineering and Physical Sciences Research Council (EPSRC) (grant no. EP/G03706X/1) and from the Royal Society.

Footnotes

One contribution of 15 to a theme issue ‘Multi-resolution simulations of intracellular processes’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4444454.