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

This work investigates multi-resolution methodologies for simulating dimer models. The solvent particles which make up the heat bath interact with the monomers of the dimer either through direct collisions (short-range) or through harmonic springs (long-range). Two types of multi-resolution methodologies are considered in detail: (a) describing parts of the solvent far away from the dimer by a coarser approach; (b) describing each monomer of the dimer by using a model with different level of resolution. These methodologies are then used to investigate the effect of a shared heat bath versus two uncoupled heat baths, one for each monomer. Furthermore, the validity of the multi-resolution methods is discussed by comparison to dynamics of macroscopic Langevin equations.


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][4][5][6]. Inter-particle forces in MD models include combinations of short-range and long-range interactions [7,8]. In allatom 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][12][13][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][12][13][14][15][16]. License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.
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][22][23][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][26][27][28][29], methods for coupling Brownian dynamics approaches with latticebased stochastic reaction-diffusion models [30 -32] or methods which make use of continuum mean-field equations for the macroscopic component of the simulation [33][34][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-coarsegraining [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.

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 1 ¼ [X 1;1 , X 1;2 , X 1;3 ] and X 2 ¼ [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 2 2 X 1 , and by R its magnitude R ¼ jRj. The interaction between monomers is given in terms of the potential F ; F(R):[0, 1) ! R, which generates a force on each of the monomers with magnitude F 0 (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: is the centre of mass of the dimer and V ¼ (V 1 þ 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 V in the following form: Therefore, we have Integrating over t, we deduce Taking the difference of equation (2.4) minus equation (2.2), implementing the over-damped assumption (where g 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 [2F(R)/(MDg)]. 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 ' 0 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 1 ¼ MDg=(k ' 2 0 ) ( 1. Then, L d can be calculated as which is valid up to the first order in 1. 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.

Short-range interaction heat bath
We describe the two monomers as balls with radius r 0 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 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 2 2 X 1 .
Our short-range interaction heat bath is described in terms of positions x j i and velocities v j i , 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 m by We are interested in the parameter regime where m ) 1.
Our MD model is based on elastic collisions of heavy monomers (balls with mass M and radius r 0 ) 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]: where v j i 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 k denote tangential components.
Heat bath models based on elastic collisions (3.4) and (3.5) have been studied by a number of authors [11][12][13][14]. Consider a single monomer in infinite domain R 3 , and let the heat bath consist of an infinite number of particles with royalsocietypublishing.org/journal/rsfs Interface Focus 9: 20180070 positions distributed according to the spatial Poisson process with density (3: 6) This means that the number of points in a subset V of R 3 has its probability mass function given by the Poisson distribution with mean l m jVj, where jVj is the volume of V. Let the velocities of the heat bath particles be distributed according to the Maxwell -Boltzmann distribution where 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 m ! 1.
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 longrange (in §4) interactions, we study whether there is a significant difference between the one-bath and two-bath models as we vary ' 0 , the separation distance, introduced in equation (2.8). In order to study this problem, we make use of multi-resolution modelling.

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 multiresolution approach where we only model the solvent that is within the close vicinity of the dimer. We consider a comoving cubic frame of length L that is centred at X 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 l m given in (3.6) and the velocities are distributed according to f m (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 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, R 3 , we need to introduce new solvent particles at the boundary of the frame. Consider that time is discretized using small time step Dt, i.e. if the current time is t, we want to calculate the state of the system at time t þ Dt. In our simulations of the multiresolution model, we need the probability of introducing a particle at a boundary of the frame (3.10) in a timestep of length Dt 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 comoving frame which over an interval of length Dt has velocity The frame is always translated to occupy the region [0, L] 3 . Thus, the velocities for the solvent particles in the new reference frame are given by w j ¼ v j À V f . We first calculate the density of particles that enter the frame via a particular boundary within a timestep of length Dt. Take, as an illustrative example, the boundary face corresponding to fx 1 ¼ 0g. Consider particles which are in half-space (À1, 0) Â 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, 1) Â R 2 at time t þ Dt. Their density, h(x 1 ), only depends on their first coordinate x 1 [ (0, 1). We can calculate h(x 1 ) by integrating density (3.6) and (3.7) over solvent particles which are at Figure 1. A diagrammatic representation of multi-resolution approaches for a dimer in a heat bath with short-range interactions. (a) Simulation of the whole dimer in a co-moving frame. The green box depicts the co-moving frame that is centred about the dimer. The blue dots correspond to solvent molecules that are explicitly modelled. Solvent molecules are not explicitly modelled in the external grey regions. (b) Simulation of one monomer in a co-moving frame. (c) Simulation with a fixed region of space where an MD model is explicitly used. A dimer molecule can move to the grey region where it is simulated using the Langevin description.
royalsocietypublishing.org/journal/rsfs Interface Focus 9: 20180070 at time t and have the appropriate velocity to reach x 1 [ (0, 1) at time t þ Dt, namely as [14] h( where V f;1 is the first component of the frame velocity and gives us the average number of particles that have entered the frame from the fx 1 ¼ 0g boundary in a time interval of length Dt as In our simulations, we choose a timestep small enough that p in ( 1, we can therefore use p in as the probability of introducing a new solvent particle. Let z ¼ [z 1 ; z 2 ; z 3 ] be the position of the new solvent particle in the coordinate system of the comoving frame. Then coordinates z 2 and z 3 are uniformly distributed in (0, L) and the first coordinate can be sampled from the error function distribution where C 1 is a normalizing constant. Then the position of the new solvent particle in the original coordinates is The velocity, w, of the new particle in the co-moving frame must have a first coordinate exceeding z 1 /Dt in order to reach z 1 in a time interval of length Dt. Noting that w ¼ v new À V f we write down the distribution of the velocity as the following truncated Gaussian distribution 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 ffiffiffi ffi p p 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 b [ R is a constant and C 3 (b) is the normalizing constant given by The algorithm in table 1 does this by generating an exponentially distributed random number h 3 with mean a 1 (b), where To maximize the acceptance probability of this algorithm, we choose its second parameter, a 2 (b), as Then its acceptance probability is depending on b as : We plot the acceptance probability (3.20) in figure 2 for our choices (3.18) and (3.19) of a 1 (b) and a 2 (b) as the solid line. We observe that the acceptance probability (3.20) for b ¼ 0 is equal to 2=p % 63:7%. This value can be improved [14] in the case of b ¼ 0 to 86.3% provided that we choose a 1 ¼ 0.532 and a 2 ¼ 0.814. To obtain a similar improvement for all values of b, we could choose both a 1 (b) and a 2 (b) to maximize the acceptance probability (3.20), rather than postulating that a 1 (b)  Table 1. Acceptance -rejection algorithm for sampling random numbers according to the probability distribution p(z; b) given by (3.16).
-Compute an exponentially distributed random number h 3 by , then choose h 3 as a sample from the probability distribution (3.16). Otherwise, repeat the algorithm.
royalsocietypublishing.org/journal/rsfs Interface Focus 9: 20180070 is given by the piecewise defined function (3.18) and optimizing a 2 (b) only. The acceptance probability (3.20) of the resulting algorithm (which would have a 1 (b) and a 2 (b) 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 1 (b) and a 2 (b) 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 f;1 Dt) (using the acceptance-rejection algorithm in table 1 for b ¼ V f;1 Dt)) and multiplying them by the factor s m Dt ffiffi ffi 2 p . One iteration (i.e. an update of the state of the system from time t to time t þ Dt) 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 Dt is chosen small enough so that at most one collision happens per iteration.
We initialize the two monomers with a separation distance ' 0 and generate a Poisson number (with mean l m L 3 ) 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 þ Dt] using the 'free-flight' positions for each monomer and solvent particle, namely we usê where i ¼ 1, 2 is the monomer number and j ¼ 1, 2, . . . , N(t), is the number of the heat bath particle. Since Dt is chosen so small that only one collision happens during the time interval [t, t þ Dt), most of the 'free-flight' positions of solvent particles are accepted in step [S2] as their updated positions x j i (t þ Dt) 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 þ Dt]. We discretize (3.1) and (3.2) using the forward Euler method as follows: and In step [S6], we use probability p in , given by (3.13), to check whether any solvent particle entered the simulation domain during the time interval (t, t þ Dt]. Since p in is the probability of entering the domain through one of its six sides, we can, for time step Dt small enough that 6p in ( 1, introduce at most one solvent particle through a randomly chosen side with probability 6p 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 in . To derive equation (3.13), we integrated over the half-space (À1, 0) Â 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 0 when Dt is not sufficiently small as boundary effects become more pronounced. To compensate for this effect, we consider the sampled position, x new and velocity v new of the new incoming particle at time t þ Dt 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. Table 2. One iteration of the multi-resolution simulation algorithm of the dimer in a co-moving frame. [S1] Update [S4] Calculate the new centre of the co-moving frame, X f (t þ Dt), by (3.9). Update N(t) by removing solvent particles which now lie outside of the frame (3.10) from the simulation.
[S5] Calculate the velocity of the frame, V f , over the interval [t, t þ Dt] by equation (3.11).
[S6] Generate two random number r 1 and r 2 uniformly distributed in interval (0, 1). If r , 6p in , then choose a side of the cube at random and generate proposed position x new and velocity v new of the new solvent particle according to distributions (3.14) and (3.15).
If r 2 , h acc (x new , v new ), then increase N(t) ¼ 1 and initialize the new solvent particle at position x new with velocity v new .
Continue with step [S1] using time t ¼ t þ Dt.
royalsocietypublishing.org/journal/rsfs Interface Focus 9: 20180070 Namely, we use the acceptance probability in step [S6] given by where Y j , R 3 is the region of the space which consists of points which have exactly j of their coordinates outside of the interval [2L/2, L/2]. For example, in our two-dimensional diagrammatic representation in figure 1a, the lighter grey shading corresponds to region Y 1 while the darker grey shading corresponds to region Y 2 .
In our illustrative simulations, we use algorithm [S1]-[S7] from table 2 together with parameter values r 0 ¼ 0.08, g ¼ 10, D ¼ 1, m ¼ 10 3 , k ¼ 10 6 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 X f (t) ¼ 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 l m given by (3.6). The velocities are distributed according to the Maxwell-Boltzmann distribution f m (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 m ! 1. 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.
In figure 3, we set the separation distance to be ' 0 ¼ ar 0 where a ! 2, such that at this distance apart the monomers are not overlapping. The plot shows the two-sided 99% confidence intervals for (L d 2 ' 0 )/' 0 for a [ f2.25, 2.5, 2.75, 3, 3.5, 4, 4.5, 5g. Firstly, we note that L d . ' 0 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 comoving 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.

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][39][40][41][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. X f (t) ¼ X 1 (t).
One iteration of the algorithm is presented as algorithm where j 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 ( [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-  royalsocietypublishing.org/journal/rsfs Interface Focus 9: 20180070 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 (t), 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).
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 g can be estimated simultaneously. We found that D % 1.0714, which is higher than our parameter value D ¼ 1, and g % 9.6064, which is lower than g ¼ 10 used in our simulations. This could suggest that the value of l m is too low or that of s m 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.

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 coarsegrained 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]: Table 3. One iteration of the multi-resolution simulation algorithm of the dimer in the heat bath with short-range interactions, where the second monomer is simulated by the Langevin dynamics. Update the velocity of the first monomer by (3.23) and the velocity of the second monomer by (3.25).
[M4] Calculate the new centre of the co-moving frame as X f (t þ Dt) ¼ X 1 (t). Update N(t) by removing solvent particles which now lie outside of the frame (3.10) from the simulation. Use steps [S5]-[S6] from the algorithm in table 2 to introduce new solvent particles into the comoving frame (3.10).
where x j i 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 a i,j , j ¼ 1, 2, . . . , N i , i ¼ 1, 2. Equations (4.1) and (4.2) are coupled with the evolution equations for solvent particles. We assume that v j i 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 m dv for j ¼ 1, 2, . . . , N i and i ¼ 1, 2. Unlike in some fictitious particle models [50][51][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 1 and N 2 are considered to satisfy N 1 ) 1 and N 2 ) 1. We are therefore working 'close' to the limit N 1 ! 1 and N 2 ! 1, 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] x is the initial position of the jth heat bath particle corresponding to the ith monomer, v j i (0) is its initial velocity and v i,j ¼ (k i,j /m) 1/2 is its frequency. Substituting for x j 1 and x j 2 in dimer's equations of motion (4.1) and (4.2), we obtain the following coupled system of generalized Langevin equations: where the friction kernel k i (t) and noise term j i ; j i (t) ¼ [j i;1 , j i;2 , j i;3 ] are given by and for i ¼ 1, 2. We assume that initial positions and velocities of solvent oscillators, x hj i;j (t)j i;n (t À t)i ¼ 2k B Td j,n k i (t), where k B is the Boltzmann constant and T is the absolute temperature. Next, we assume that the frequencies v i,j are sampled from a (continuous) exponential distribution with mean v and we set our interaction constants equal to where g . 0 is the friction constant used in equations (2.2) and (2.4). Then friction kernel (4.7) becomes cos (v i,j t): Passing to the limit N i ! 1 allows us to consider the above summation as a continuous integral over the distribution of oscillator frequencies, with both friction kernels k 1 (t) and k 2 (t) converging to the same friction kernel [54] (4:8) Then we can define the limiting friction kernel by which, for our choice of oscillators' frequencies and interaction terms (4.7), satisfies k 1 (t) ¼ 0 for t . 0 and k 1 (0) ¼ 1. 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 v ! 1, 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 B T/(gM).
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 shortrange interaction MD model in §3.1 with the help of multiresolution 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][56][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 ' 0 as in figure 3, expressed as royalsocietypublishing.org/journal/rsfs Interface Focus 9: 20180070 a-multiples of r 0 , although our long-range interaction model does not make use of parameter r 0 . The value of L d is for each value of ' 0 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 ' 0 and sample oscillators' frequencies, v i,j , according to the exponential distribution with mean v. 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.
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 where the heat bath evolution equation (4.4) includes terms corresponding to both monomers m dv j dt ¼ Àk 1,j x j À a 1,j X 1 À k 2,j x j À a 2,j X 2 , for j ¼ 1, 2, . . . , N. Our results will then depend how we choose parameters k i,j and a i,j . For example, if we choose k i,j and a i,j to be the same for both monomeres, i.e. k 1,j ¼ k 2,j ¼ k j and a 1,j ¼ a 2,j ¼ a j , ( 4 :9) for j ¼ 1, 2, . . . , N, then the oscillating frequency of the jth heat bath oscillator is v j ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi 2k j =m p 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 a j , i.e. using a 2 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 ' 0 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, ' 0 . 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 1 and N 2 oscillators, i.e. N ¼ N 1 þ N 2 , and choose our parameters k i,j and a 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 a i,j , the single heat bath case can both increase or decrease the average length of the dimer.

Multi-resolution modelling of dimer
In figure 1b, we use our dimer example to illustrate a multiresolution approach which models a part of a molecule  Table 4. One iteration of the multi-resolution simulation algorithm of the dimer in the heat bath with long-range interactions, where the second monomer is simulated by the Langevin dynamics. [L1] Update velocities of the dimer and solvent particles for a half time step using (4.11).

[L2]
Update positions of the dimer and solvent particles using (4.12).
[L4] Update velocities of the dimer and solvent particles for a half time step using (4.13).
[L5] Continue with step [L1] using time t ¼ t þ Dt. and where A i , for i ¼ 1, 2, is the acceleration of the corresponding monomer. For the first monomer, its acceleration A 1 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 j 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 1 , V 1 and A 1 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 v ¼ 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 1 ¼ 10 5 , 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.
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).

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][39][40][41][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][51][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 royalsocietypublishing.org/journal/rsfs Interface Focus 9: 20180070 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, 1) Â 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, R 3 . However, we now want to replace solvent particles which are in {x 1 , b} ¼ (À1, b) Â 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 fx 1 . bg and fx 1 , bg. In figure 1c, we present a two-dimensional illustration of a monomer when it intersects the interface, fx 1 ¼ bg. Such a monomer is subject to the collisions with heat bath particles on the part of its surface which lies in fx 1 . bg. This has to be compensated by using a suitable random force from fx 1 , bg, 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][68][69] can be generalized to the analysis of behaviour of molecules close to hybrid interfaces in multiresolution 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.