Varying the resolution of the Rouse model on temporal and spatial scales: application to multiscale modelling of DNA dynamics

A multi-resolution bead-spring model for polymer dynamics is developed as a generalization of the Rouse model. A polymer chain is described using beads of variable sizes connected by springs with variable spring constants. A numerical scheme which can use different timesteps to advance the positions of different beads is presented and analyzed. The position of a particular bead is only updated at integer multiples of the timesteps associated with its connecting springs. This approach extends the Rouse model to a multiscale model on both spatial and temporal scales, allowing simulations of localized regions of a polymer chain with high spatial and temporal resolution, while using a coarser modelling approach to describe the rest of the polymer chain. A method for changing the model resolution on-the-fly is developed using the Metropolis-Hastings algorithm. It is shown that this approach maintains key statistics of the end-to-end distance and diffusion of the polymer filament and makes computational savings when applied to a model for the binding of a protein to the DNA filament.

1. Introduction. Over the past 70 years, there have been multiple attempts to dynamically model the movement of polymer chains with Brownian dynamics [24,37,38,47], which have more recently been used as a model for DNA filament dynamics [2]. One of the first and simplest descriptions was given as the Rouse model [38], which is a bead-spring model [2], where the continuous filament is modelled at a mesoscopic scale with beads connected by springs. The only forces exerted on beads are spring forces from adjacent springs, as well as Gaussian noise. Hydrodynamic forces between beads and excluded volume effects are neglected in the model in favour of simplicity and computational speed, but the model manages to agree with several properties of polymer chains from experiments [24,34]. Other models exist, for example the Zimm model introduces hydrodynamic forces [47] between beads, or bending potentials can be introduced to form a wormlike chain and give a notion of persistence length [1], see, for example, review article [2] or books [6,7] on this subject.
Most of the aforementioned models consider the filament on only a single scale. In some applications, a modeller is interested in a relatively small region of a complex system. Then it is often possible to use a hybrid model which is more accurate in the region of interest, and couple this with a model which is more computationally efficient in the rest of the simulated domain [8,9,13]. An application area for hybrid models of polymer chains is binding of a protein to the DNA filament, which we study in this paper. The model which we have created uses Rouse dynamics for a chain of DNA, along with a freely diffusing particle to represent a binding protein. As the protein approaches the DNA, we increase the resolution in the nearby DNA filament to increase accuracy of our simulations, whilst keeping them computationally efficient.
In this paper we use the Rouse model for analysis due to its mathematical tractability and small computational load. Such a model is applicable to modelling DNA dynamics when we consider relatively low resolutions, when hydrodynamic forces are negligible and persistence length is significantly shorter than the Kuhn length between each bead [2]. The situation becomes more complicated when we consider DNA modelling at higher spatial resolutions.
Inside the cell nucleus, genetic information is stored within strands of long and thin DNA fibres, which are separated into chromosomes. These DNA fibres are folded into structures related to their function. Different genes can be enhanced or inhibited depending upon this structure [12]. Folding also minimises space taken up in the cell by DNA [21], and can be unfolded when required by the cell for different stages in the cell cycle or to alter gene expression. The folding of DNA occurs on multiple scales. On a microscopic scale, DNA is wrapped around histone proteins to form the nucleosome structure [45]. This in turn gets folded into a chromatin fibre which gets packaged into progressively higher order structures until we reach the level of the entire chromosome [12]. The finer points of how the nucleosome packing occurs on the chromatin fibre and how these are then packaged into higher-order structures is still a subject of much debate, with long-held views regarding mesoscopic helical fibres becoming less fashionable in favour of more irregular structures in vivo [28].
In the most compact form of chromatin, many areas of DNA are not reachable for vital reactions such as transcription [12]. One potential explanation to how this is overcome by the cell is to position target DNA segments at the surface of condensed domains when it is needed [5,27], so that transcription factors can find expressed genes without having to fit into these tightly-packed structures. This complexity is not captured by the multiscale model of protein binding presented in this paper. However, if one uses the developed refinement of the Rouse model together with a more detailed modelling approach in a small region of DNA next to the binding protein, then such a hybrid model can be used to study the effects of microscopic details on processes over system-level spatial and temporal scales. When taking this multiscale approach, it is necessary to understand the error from including the less accurate model in the hybrid model and how the accuracy of the method depends on its parameters. These are the main questions studied in this paper.
The rest of the paper is organized as follows. In Section 2, we introduce a multiresolution bead-spring model which generalizes the Rouse model. We also introduce a discretized version of this model which enables the use of different timesteps in different spatial regions. In Section 3, we analyze the main properties of the multiresolution bead-spring model. We prove two main lemmas giving formulas for the diffusion constant and the end-to-end distance. We also study the appropriate choice of timesteps for numerical simulations of the model and support our analysis by the results of illustrative computer simulations. Our main application area is studied in Section 4 where we present and analyze a DNA binding model. We develop a method to increase the resolution in existing segments on-the-fly using the Metropolis-Hastings algorithm. In Section 5, we conclude our paper by discussing possible extensions of the presented multiscale approach (by including more detailed models of DNA dynamics) and other multiscale methods developed in the literature. 2. Multi-resolution bead-spring model. We generalize the classical Rouse bead-spring polymer model [38] to include beads of variable sizes and springs with variable spring constants. In Definition 2.1, we formulate the evolution equation for this model as a system of stochastic differential equations (SDEs). We will also introduce a discretized version of this model in Algorithm 1, which will be useful in Sections 3 and 4 where we use the multi-resolution bead-spring model to develop and analyze multiscale models for DNA dynamics.
Definition 2.1. Let N > 1 be a positive integer. A multi-resolution bead-spring polymer of size N consists of a chain of N beads of radius σ n , for n = 1, 2, . . . , N , connected by N − 1 springs which are characterized by their spring constants k n , for n = 1, 2, . . . , N − 1. The positions r n ≡ r n (t) = [r n,1 (t), r n,2 (t), r n,3 (t)] of beads evolve according to the system of SDEs (for n = 1, 2, . . . , N ) ζ n dr n = k n−1 r n−1 (t) − (k n−1 + k n )r n (t) + k n r n+1 (t) dt + 2k B T ζ n dW n , (2.1) where ζ n = 6πησ n is the frictional drag coefficient of the n-th bead given by Stokes' Theorem, η is the solvent viscosity, dW n ≡ [dW n,1 , dW n,2 , dW n,3 ] is a Wiener process, T is absolute temperature, k B is Boltzmann's constant and we assume that each spring constant k n can be equivalently expressed in terms of the corresponding Kuhn length b n by We assume that the behaviour of boundary beads (for n = 1 and n = N ) is also given by equation (2.1) simplified by postulating r 0 (t) = r 1 (t) and r N +1 (t) = r N (t).
In Figure 2.1, we schematically illustrate a multi-resolution bead-spring polymer for N = 21. The region between the 8-th and the 14-th bead is described with the highest resolution by considering smaller beads and springs with larger spring constants (or equivalently with smaller Kuhn lengths). The scalings of different parameters in Definition 2.1 are chosen so that we recover the classical Rouse model [38] if we assume σ 1 = σ 2 = · · · = σ N = σ and b 1 = b 2 = · · · = b N −1 = b. Then equation (2.1) simplifies to where ζ = 6πησ, k = 3k B T /b 2 and we again define r 0 (t) = r 1 (t) and r N +1 (t) = r N (t) in equations for boundary beads. In the polymer physics literature [7], the Rouse model (2.2) is equivalently written as where random thermal noises f n (t) exerted on the beads from Brownian motion are characterized by the moments [7] f n (t) = 0, where n, m = 1, 2, . . . , N and i, j ∈ {1, 2, 3}. For the remainder of this paper, we will use the SDE notation as given in (2.2), because we will often study numerical schemes for simulating polymer dynamics models. The simplest discretization of (2.2) is given by the Euler-Maruyama method [20], which uses the finite timestep ∆t and calculates the position vector r t n of the n-th bead, n = 1, 2, . . . , N , at discretised time t by for ξ n = (ξ n,1 , ξ n,2 , ξ n,3 ), where ξ n,i is normally distributed random variable with zero mean and unit variance (i.e. ξ n,i ∼ N (0, 1)) for i = 1, 2, 3. In order to discretize the multi-resolution bead-spring model, we allow for variable timesteps.
Definition 2.2. Let ∆t > 0 and let j n , n = 1, 2, . . . , N − 1, be positive integers such that j n−1 | j n or j n | j n−1 for n = 2, 3, . . . , N − 1. Let us assume that at least one of the values of j n is equal to 1. We define ∆t n = j n ∆t for n = 1, 2, . . . , N − 1 and we call ∆t n a timestep associated with the n-th spring.
Definition 2.2 specifies that all timesteps must be integer multiples of the smallest timestep ∆t. The timesteps associated with two adjacent springs are also multiples of each other. The time evolution of the multi-resolution bead-spring model is computed at integer multiples of ∆t. One iteration of the algorithm is shown in Algorithm 1. The position of the n-th bead is updated at integer multiples of min{j n−1 , j n }∆t = min{∆t n−1 , ∆t n } by calculating the random displacement due to Brownian motion, with displacement caused by springs attached to the bead also updated at integer multiples of the timesteps associated with each spring, i.e. ∆t n−1 or ∆t n . Considering the situation that all beads, springs and timesteps are the same, then one can easily deduce the following result. Lemma 2.3. Let σ > 0, ζ > 0, k > 0 and ∆t > 0 be positive constants and N > 1 be an integer. Consider a multi-resolution bead-spring polymer of size N with σ n = σ, ζ n = ζ, for n = 1, 2, . . . , N , and k n = k, for n = 1, 2, . . . , N − 1. Let the timesteps associated with each spring be equal to ∆t, i.e. j 1 = j 2 = · · · = j n = 1 in Definition 2.2. Then Algorithm 1 is equivalent to the Euler-Maruyama discretization of the Rouse model given as equation (2.5).
Put r i∆t N := r (i−1)∆t N . end Algorithm 1: One iteration of the numerical algorithm for simulating the multi-resolution bead-spring model introduced in Definition 2.1. We update positions of beads at time t = i∆t where i is a positive integer and each spring is simulated using its associated timestep given in Definition 2.2. Lemma 2.3 shows that the multi-resolution bead-spring model is a generalization of the Rouse model. In the next section, we will study properties of this model which will help us to select the appropriate parameter values for this model and use it in multiscale simulations of DNA dynamics.
3. Macroscopic properties and parameterizations of multi-resolution bead-spring models. We have formulated a multiscale Rouse model which varies the Kuhn lengths throughout the filament, but we would like to keep properties of the overall filament constant regardless of the resolution regime being considered for the filament. We consider a global statistic for the system to be consistent if the expected value of the statistic is invariant to the resolution regime being considered for the filament. We consider the self diffusion constant and root mean squared (rms) end-to-end distance as two statistics we wish to be consistent in our system, which can be ensured by varying the bead radius and the number of beads respectively. The precise way to vary these properties will be explored in this section.
3.1. Self diffusion constant. The self diffusion constant is defined as where r t G is the centre of mass of the polymer chain at time t, which is defined by is an extension to the definition given by Doi and Edwards [7] for the centre of mass of a continuous chain on only one scale. If all beads have the same radius σ (i.e. if σ n = σ for n = 1, 2, . . . , N ), then equation (3.2) simplifies to the centre of mass definition for the classical Rouse model. In this case, the self diffusion constant is given by [7] where N is the number of beads. This result explains the, on the face of it, counterintuitive scaling of equation (3.2) with σ n . If we suppose that each bead had the same density, then the mass of each bead would be proportional to its volume, i.e. to σ 3 n . However, in definition (3.2), we have used weights σ n instead of σ 3 n , because beads do not represent physical bead objects like nucleosomes, but representations of the filament around it, so the bead radius scales with the amount of surrounding filament, which is linear in bead radius in this formulation. If we consider DNA applications, we could imagine each bead as a tracker for individual base pairs at intervals of, say, thousands of base pairs away from each other along the DNA filament. The filament in the model is then drawn between adjacent beads.
This linear scaling with σ n can also be confirmed using equation (3.3) for the classical Rouse model. If we describe the same polymer using a more detailed model consisting of twice as many beads (i.e. if we change N to 2N ), then we have to halve the bead radius (i.e. change σ to σ/2) to get a polymer model with the same diffusion constant (3.3). In particular, the mass of a bead scales with σ (and not with σ 3 ). In the next lemma, we extend result (3.3) to a general multi-resolution bead-spring model. Lemma 3.1. Let us consider a multi-resolution bead-spring polymer of size N and a set of timesteps associated with each spring satisfying the assumptions of Definitions 2.1 and 2.2. Then the self diffusion constant of the polymer evolution described by Algorithm 1 is given by Proof. Algorithm 1 describes one iteration of our numerical scheme. Multiplying the steps corresponding to the n-th bead by σ n and summing over all beads, we obtain how Ω r i∆t G changes during one timestep ∆t. Since ζ n = 6πησ n , tension terms cancel after summation and the evolution rule for Ω r i∆t G simplifies to where ξ i n ∼ [N (0, 1), N (0, 1), N (0, 1)] and function Q(j, i) is defined for positive integers j and i by Let us denote by H the least common multiple of {j 1 , j 2 , . . . , j n }. Every bead is updated in Algorithm 1 at integer multiples of H∆t. We can eliminate function Q from equation (3.5) if we consider the evolution of Ω r t G when time t is evaluated at integer multiples of H∆t. We obtain the evolution rule where we used the fact that the sum of normally distributed random variables is again normally distributed. Dividing equation (3.6) by Ω, we obtain Using definition (3.1), we obtain (3.4).

The formula (3.4) is a generalization of equation (3.3) obtained for the Rouse model.
It is invariant to the resolutions provided that the mass of the filament Ω remains constant through selection of the number of beads and bead radius, therefore the self diffusion constant is consistent.

3.2.
End-to-end distance. We define the end-to-end vector R = r N − r 1 from one end of the filament to the other [7]. An important statistic to consider related to this is the root mean squared (rms) end-to-end distance of the filament µ = R 2 1/2 . The expected value of the long-time limit of the rms end-to-end distance, denoted µ ∞ , for the classical Rouse model is given by [7] µ ∞ = lim We generalize this result in the following lemma.
for n = 1, 2, 3, . . . , N − 1, (3.8) and the long-time limit of the rms end-to-end distance is given by Proof. Equations (2.1) describe a system of 3N linear SDEs. However, the SDEs corresponding to different spatial dimensions are not coupled. We therefore restrict our investigation to the behaviour of the first coordinates of each vector in (3.8). Let us arrange the differences of the first coordinates of subsequent beads into the (N − 1)-dimensional vector Then SDEs (2.1) can be rewritten to the system of SDEs for y(t) in the matrix form The stationary covariance matrix, defined by is the solution of Lyapunov equation [18] AC + CA T + BB T = 0. It can be easily verified that the unique solution of this equation Multiplying this result by 3 (the number of coordinates), we obtain (3.8). The end-to-end distance can be rewritten as R = r N − r 1 = N n=2 (r n − r n−1 ). Substituting into (3.9), using (3.8) and the fact that the stationary covariance matrix C is diagonal, we obtain (3.9).
3.3. Optimal model refinement in time and space. Lemmas 3.1 and 3.2 describe theoretical results which have been derived under slightly different assumptions. Lemma 3.1 is formulated as a property of Algorithm 1, but the same result, equation (3.4), also holds when we calculate the self-diffusion coefficient of the SDE formulation of the multi-resolution bead-spring model (2.1). Algorithm 1 is designed in such a way that all force terms corresponding to springs cancel when the evolution equation for r i∆t G is derived (see equation (3.5)). In particular, Lemma 3.1 holds for any choices of the lengths of timesteps associated with different springs.
On the other hand, Lemma 3.2 describes the property of the SDE formulation (2.1). If we use a discretized version of (2.1), then we introduce a discretization error. This error can be made smaller by choosing smaller timesteps. In this section, we show that the smallest timesteps are only required in the regions with the highest spatial resolution. We define a family of optimal multi-resolution (OMR) models designed to have macroscopic properties invariant to resolution regime. (3.10) The first region contains springs indexed by {1, 2, . . . , N 1 } and the j-th region, j = 2, 3, . . . , R, contains springs indexed by {1 + Let us associate with each region an integer resolution s j , where s j |s j−1 or s j−1 |s j , for j = 2, 3, . . . , R, and s 2 j |N j for j = 1, 2, 3, . . . , R, with at least one region in resolution 1 which is the region with the finest detail. Larger values of s j represent coarser representations of the filament. We define the OMR model as the multi-resolution bead-spring model which consists of R regions of consecutive beads and springs. In the j-th region, we have N j springs with Kuhn length b j and associated time steps ∆t j given by

11)
where σ j is the radius of beads which are connected to two springs which have the same Kuhn length b j . We assume that the bead radius of beads on region boundaries sharing springs with Kuhn lengths b j−1 and b j is ( σ j−1 + σ j )/2, for j = 2, 3, . . . , R. Moreover, we assume that the bead radius of the first and last bead of the polymer chain are equal to (σ + σ 1 )/2 and (σ + σ R )/2, respectively.
Substituting scalings (3.11) into (3.2), and using (3.10), we obtain that the OMR model satisfies Substituting into (3.4), we deduce that the OMR model has the same self-diffusion constant as the original detailed model (given by (3.3)). Considering the limit ∆t → 0, we can use Lemma 3.2 and scalings (3.11) to derive the expected rms end-to-end distance for the filament: which is again independent of the choice of resolutions s j , j = 1, 2, . . . , R. As the Kuhn length and bead radius vary across resolutions, it is important to consider the numerical stability of the model [40]. We choose timesteps to be sufficiently small so that solutions do not grow exponentially large. In discretized equations of Algorithm 1, drift terms appear in the form (k n /ζ n )(r n+1 − r n )∆t n , which is proportional in the j-th region of the OMR model to Using scalings (3.11) and assuming that (r n+1 − r n ) is of the same order as the Kuhn length b j , we obtain that the size of (3.12) scales with s j . Assuming that ∆t is chosen in the original fine scale model so that (k/ζ)(r n+1 − r n )∆t is small compared to the Kuhn length b, then the drift term of the OMR model, given by (3.12) is also small compared to b j , the characteristic lengthscale of the OMR model in the j-th region, j = 1, 2, . . . , R.
Next, we compare the number of calculations made by the original detailed singlescale Rouse model with the OMR model. The j-th region has N j springs simulated with timestep ∆t j . Using scalings (3.11), we obtain that we use s 6 j -times fewer calculations in the j-th region by advancing fewer beads over larger timesteps. Assuming that the computational intensity of the simulation of the detailed model in each region is proportional to the size of the region, N j /(N − 1), we can quantify the fraction of computational time which is spent by the OMR model (as compared to the detailed model) by (3.13) For example, if we coarse-grained the detailed model everywhere using the integer resolution s 1 = 2, then (3.10) and (3.13) implies that we speed up our simulations by the factor of 64.

Simulation results.
In this section we show that simulations of the OMR method match the original single-scale Rouse model. We also compare this to analytic results predicted from equations (3.4) and (3.7) for the rms end-to-end distance and the self diffusion constant of a filament in an equilibrium state. For the detailed model, we choose the parameters: ∆t = 0.8 µs, b = 60 nm, σ = 1.2 nm, (3.14) where the Kuhn length is chosen to be longer than the persistence length of DNA [16], and the other parameters are chosen arbitrarily. For the remainder of this paper, we shall use k B = 1.4 × 10 −23 J K −1 , T = 300 K and η = 1 cP, the viscosity of water.  3.4.1. Comparison at equilibrium. We compare two resolution regimes for the same system, with the single scale model considering the full system in high resolution and a multiscale model considering the middle 10% of the filament in high resolution and the remainder in low resolution. The corresponding parameters of the OMR model are given in Table 3.1. The OMR model contains 69 beads connected by 68 springs, while the original detailed model is given by 501 beads connected by 500 springs.
We generate the initial configuration of the polymer filament using a multiscale generalisation of the Freely Jointed Chain (FJC) model [14,25]. The chain is generated iteratively, with the (n + 1) th bead in the chain placed uniformly at random on the surface of the sphere with radius b n centred on the n-th bead, so that a chain is produced with unconstrained random angles [2]. We run the model for one second and estimate both the rms end-to-end distance and the self diffusion constant of a filament in an equilibrium state. The results given in Table 3.2. From the results we can see that the OMR model accurately maintains macroscopic properties at a fraction of number of calculations of the detailed model, where 99.94% of the calculations used are used to update beads in the fine resolution region.
The length of the simulated DNA was chosen to be short enough for easy computation, but long enough to be comparable to real simulations of the DNA behaviour [36]. Those simulations showed the diffusion of DNA to be a thousand times larger than our simulations, this is due in part to neglecting hydrodynamic forces, which gives a different form of the expected self diffusion constant [7] compared to (3.4).

Compacting the filament.
We have shown that the dynamics of some multiscale systems match the analytic results at equilibrium predicted by equations (3.4) and (3.7). It is also important to show that non-equilibrium dynamics of the detailed system can be replicated by the multiscale systems. We consider a filament shrunk so that all beads begin at the origin and compare how the models expand towards equilibrium. We compare four systems (C1, C2, C3 and C4) given in Table 3.3 which, by construction, have the same rms end-to-end distance, but which vary in their level of detail and the number of resolutions considered.
Model C2 in Table 3.3 is the most detailed model which uses the same parameters as in (3.14). We compare the time evolution of rms end-to-end distance of the considered polymer chain models. The results (averaged over 10 3 realizations) are shown in Figure 3.1. 4. DNA binding model. To show the usefulness of the OMR model, we have devised an illustrative model for the relationship between a DNA binding protein and a segment of DNA. There is a large variety of different binding proteins, for example transcription factors or polymerases. In general, binding protein dynamics are not fully understood and vary greatly depending on the function of the protein [17]. In Brownian dynamics models, a binding protein binds to DNA with a given probability if it comes within a certain distance of the binding site on the DNA filament [10]. We will use a simple version of this approach and assume that a binding protein binds to DNA if it comes within a 'binding radius' of a bead in the bead-spring polymer chain.

Model
We apply the multiscale aspect of the model by adjusting the resolution regime of the filament dynamically, so that as the binding protein moves close to the filament, nearby filament sections increase in resolution. To increase the resolution on-the-fly, we develop a Markov Chain Monte Carlo (MCMC) scheme. The model is presented as a two-resolution model for simplicity, but it can easily be extended to more resolution levels. We denote by s the difference in resolution between high and low regions. Using the notation introduced in Definition 3.3, we consider filaments which include regions with higher resolution (where s j = 1 in Definition 3.3) and regions with lower resolution (where s j = s). The rms end-to-end distance for a compacted filament to re-extending towards the analytic (equilibrium) rms end-to-end distance in magenta for systems C1-C4 given in Table 3.3.

Resolution increase with MCMC.
We first present a framework for increasing the resolution between adjacent beads. If we wish to increase the resolution in a region of low resolution between beads r A and r B which lie a distance d apart, then by Definition 3.3 we introduce s 2 − 1 new springs with Kuhn length b. We place the new beads between r A and r B to lie a distance of d ′ apart, which is in general not equal to b, because d is in general not equal to sb. We select d ′ so that the rms end-toend distance for the new chain formed is to equal d at the time of its creation. Using equation (3.7), we require d ′ = d/s to be the distance between each bead. It should be noted that we are generating beads by the FJC model, introduced in Section 3.4.1 to be such that the distance between each new bead is the same. This is for simplicity as when we apply the dynamics there is a fast transition for these bond lengths to revert to the equilibrium distribution given by the Gaussian chain model [7].
In order to apply MCMC, we seek the probability function ψ of a chain F = {r i | i = 1, 2, . . . , M } of M = s 2 − 1 beads a distance d ′ apart between r A and r B . We first consider the probability function for the first M − 1 beads {r i<M } = {r i | i = 1, 2, . . . , M − 1}, which is proportional to the circumference of the circle of points on which the last bead r M can be placed, which we call the circle of allowable points Γ. This is given as the set or the set of all points a distance d ′ from both r M −1 and r B . In some cases Γ = ∅ when |r M −1 − r B | > 2d ′ . Note that if ψ({r i<M }) is proportional to the circumference of circle Γ it is also proportional to its radius.
Provided Γ = ∅, to find the radius of Γ we consider the triangle with points r M −1 , u which is an arbitrary point on Γ, and v which is the midpoint between r M −1 and r B . By construction there is a right angle at v. The radius of Γ is given by the distance between u and v, which is l = d ′2 − (λ/2) 2 , where λ is the distance between r M −1 Generate a proposal chainF = {r i | i = 1, 2, . . . , M − 1} with each spring length d ′ starting at r A from the FJC model.
Repeat and generate a new proposal chain. end Algorithm 2: The algorithm for picking the next chain F i+1 in the Metropolis-Hastings algorithm given the previous chain F i . and r B . We therefore find the probability density function for the chain: with normalisation constant K, r 0 = r A and δ the Dirac delta function. To select from the probability distribution ψ, we apply the Metropolis-Hastings algorithm [4] as outlined in Algorithm 2. Given the first M − 1 beads, and provided Γ = ∅, the density function for the final bead is given by Algorithm 2 is a Metropolis-Hastings algorithm with a candidate-generating distribution which uses rejection sampling [31] to keep selecting chainsF = {r i | i = 1, 2, . . . , M − 1} with distance d ′ between each bead until we find a chain where the second last beadr M −1 is within a distance 2d ′ from the endpoint r B , which is then used as the candidate. Based on 10 4 simulations we see approximately a 77.5% acceptance of the candidate for M ≥ 3 beads generated between beads a distance b apart, invariant to the number of beads for M ≥ 3. To ensure the Markov chain is sufficiently well mixed and that we have 'lost memory' of the initial chain, we choose the 10 th filament in the sequence generated from Algorithm 2; in 10 3 simulations with a resolution difference of 3, we see only one case where the initial distribution is the same as the 10 th .

Model description.
On top of our multiscale Rouse model, we introduce a binding protein, modelled as a diffusing particle with diffusion constant D p . The position of the binding protein at time t is denoted p(t). We assume that it evolves according to the discretized Brownian motion p(t + ∆t) = p(t) + 2D p ∆t ξ, Snapshots of the DNA binding model, taken from an illustrative OMR simulation. In (a) the filament is represented by blue in its lowest resolution with the binding protein in green. This is a snapshot from just before the protein comes within resolution increase radius r I . In (b) the protein has moved closer than this radius, so we increase the resolution in this region. The higher resolution region is shown in purple. We continue the simulation until either the protein drifts away from the filament or binds to it, as we see in (c).
where ξ ∼ [N (0, 1), N (0, 1), N (0, 1)] and ∆t is the time step of the highest resolution. Associated with the binding protein is a binding radius d b , such that if any bead along the filament comes within a distance d b from the binding protein then it 'binds' to the filament and the simulation finishes. We also consider the protein to have drifted to infinity and fail to bind if it gets a distance d ∞ from all beads which exist in the lowest resolution of the OMR model.
The binding protein also has its resolution increase radius r I (where r I > d b ), so that if the protein is within a distance r I from a bead with an adjacent spring in low resolution then we change the spring to be high resolution and introduce new beads as described in Section 4.1, and reformulate bead radius, timestep and Kuhn length according to the OMR model in Definition 3.3. See Figure 4.1 for an example of the resolution increasing as the filament comes within range of the protein.
As well as zooming in to specific areas of the filament, we include a provision that if the protein moves away far enough from the zoomed in filament section that it will not be interacting with the filament, then we zoom out. This is implemented as a zoom out factor Z, so that if the protein moves Z r I away from both beads on the boundary of a region in high resolution, then we 'zoom out' by changing the resolution of the region to the lower resolution and removing all beads inside.
To initialize the simulation, we generate the filament using the FJC model, with R regions each containing one spring and the beads on the region boundaries denoted q j (0), j = 1, 2, . . . , R + 1. Since all regions are in low resolution at time t = 0, scalings (3.11) are given by for all regions j = 1, 2, . . . , R. We use R = 100 in our illustrative simulations. We then place the protein uniformly at random on the sphere with radius d 0 centred at the middle bead of the chain (i.e. at the (1+R/2)-th bead). Once the initial configuration has been generated, we compute the time evolution of the filament and the protein iteratively using Algorithm 3, which describes one timestep of the method. At each time step, the multi-resolution bead-spring polymer is described as in Definition 2.1 by the positions of N ≡ N (t) beads at r n ≡ r n (t) = [r n,1 (t), r n,2 (t), r n,3 (t)], n = 1, 2, . . . , N (t). We initially have N (0) = R + 1 beads at positions r n (0) = q n (0), for Calculate the position of the chain at time t + ∆t by Algorithm 1. Update both r n (t + ∆t) for n = 1, 2, . . . , N (t) and q j (t + ∆t) for j = 1, 2, . . . , R + 1. Calculate the protein position at time t + ∆t by equation (4.1).
The protein binds to the filament. STOP the simulation.
The protein drifts to infinity. STOP the simulation. else for j = 1, 2, . . . , R do if (the j-th region is in the low resolution) then if |q j − p| < r I OR |q j+1 − p| < r I then Introduce new beads according to Algorithm 2. Adjust bead radius, timestep and Kuhn length in the j-th region using Definition 3.3, i.e. use σ, ∆t and b, respectively. end else (the j-th region is in the high resolution) if |q j − p| > Z r I AND |q j+1 − p| > Z r I then Remove all beads inside the j-th region. Adjust bead radii σ j , timestep ∆t j and Kuhn length b j in the j-th region by (4.2). end end end end Algorithm 3: One iteration of the algorithm for the DNA binding model. n = 1, 2, . . . , N (0). As time progresses, some regions are refined, so the value of N (t) changes and indices of some beads are relabelled. To simplify the presentation of Algorithm 3, we denote by q j ≡ q j (t), j = 1, 2, . . . , R + 1, the positions of boundary beads of each region at time t, independently of the fact whether the chain was refined or not in the corresponding region.
We have confirmed the consistency of the model where resolutions are changed dynamically against analytic results. It performs well with a rms end-to-end distance of 2.997 µm and diffusion of 7.49·10 −5 µm 2 s −1 compared to expected results from (3.4) and (3.9) of 3 µm and 7.34 · 10 −5 µm 2 s −1 , respectively, after 10 3 simulations, with parameters given in (3.14) and Table 4.1.

Results.
We compare the results between the OMR model with the classical Rouse model, with the entire filament in the highest resolution. We run Algorithm 3 with the parameters given in (3.14) and Table 4.1. At different initial starting distances the model runs until the protein is either bound or has escaped the filament. In Figure 4.2, we present the probability of binding of the protein to the filament, P b (d 0 ), as a function of the initial distance d 0 ∈ [d b , 10 d ∞ ] of the protein from the middle bead of the filament. We estimate P b (d 0 ) as a fraction of simulations which end up with the protein bound to DNA. Each data point in Figure 4  the value of P b (d 0 ) estimated from 10 3 independent realizations of the process. If d 0 < d b , then the protein is immediately bound to DNA, i.e. P b (d 0 ) = 1 for d 0 < d b .
If d 0 = d ∞ , then the probability of binding is nonzero, because the initial placement, d 0 , is the distance of the protein from the centre of the filament. In particular, the minimum distance from protein to filament is less than or equal to the initial placement distance, d 0 , and the simulations (with the possibility of binding) take place even if d 0 = d ∞ . Due to computational constraints of the single-scale model we consider a selection of initial distances at points d 0 = 10 −2+ℓ/3 µm, ℓ = 0, 1, . . . , 9 (black points), where error bars give a 95% confidence interval based on the Wilson score interval for binomial distributions [44]. We run simulations for more initial distances, d 0 = 10 −2+ℓ/9 µm, ℓ = 0, 1, . . . , 27 (blue line), using the computationally efficient OMR model and present our results as the blue line in Figure 4.2. We see that P b (d 0 ) is very similar between the single-scale and OMR models. The model also succeeds in reducing computational time. For 10 3 simulations with the protein starting 1 µm from the middle bead, with parameters given in Table 4.1, the OMR model represented a 3.2-times speedup compared to the detailed model, with only a 3-times resolution difference. We expect for larger resolution differences to see greater improvements in speed.
5. Discussion. In this paper we have extended basic filament modelling techniques to multiple scales by developing OMR methods. We have presented an MCMC approach for increasing the resolution along a static filament segment, as well as an extension to the Rouse model to dynamically model a filament which considers multiple scales. The bead radius, as well as the number of beads associated with each resolution, is altered to maintain consistency with the end-to-end distance and diffusion of a filament across multiple scales, as well as the timestep to ensure numerical convergence.
We have then illustrated the OMR methodology using a simple model of protein binding to a DNA filament, in which the OMR model gave similar results to the single- scale model. We have also observed a 3.2-times speed-up in computational time on a model which considers only a 3-times increase in resolution, which illustrates the use of the OMR approach as a method to speed up simulations whilst maintaining the same degree of accuracy as the more computationally intensive single-scale model. The speed-up in computational time could be further increased by replacing Brownian dynamics based on time-discretization (4.1) by event-based algorithms such as the FPKMC (First passage kinetic Monte Carlo) and GFRD (Green's function reaction dynamics) methods [33,43].
When considering the zooming out of the DNA binding model, note that it is generally possible to zoom in and out repetitively, as long as the dynamics are such that we can generate a high resolution structure independent from the previous one (i.e., once we zoom out, the microscopic structure is completely forgotten). However, particularly in the case of chromatin, histone modification and some DNA-binding proteins may act as long-term memory at a microscopic scale below the scales currently considered. To reflect the effect of the memory, some properties of the microscopic structure should be maintained even after zooming out. Fractal dimension may serve as a candidate of indices [29], which can be also estimated in living cells by singlemolecule tracking experiments [41].
The OMR method could be applied to modern simulations of DNA and other biological polymers which use the Rouse model [19] in situations where certain regions of the polymer require higher resolutions than other regions. The model considered in this report uses Rouse dynamics, which is moderately accurate given its simplicity, but as we zoom in further towards a binding site, then we will need to start to consider hydrodynamic forces and excluded volume effects acting between beads. Models which include hydrodynamic interactions such as the Zimm model [47] have previously been used to look at filament dynamics [1,11]. Therefore it is of interest to have a hybrid model which uses the Rouse model in low resolutions and the Zimm model in high resolutions. The combination of different dynamical models might give interesting results regarding hierarchical structures forming as we move between resolutions.
As we go into higher resolutions, strands of DNA can be modelled as smooth [2], unlike the FJC model where angles between beads are unconstrained. The wormlike chain model of Kratky and Porod [23], implemented via algorithm by Hagermann and Zimm [15], gives a non-uniform probability distribution for the angles between each bead. Allison [1] then implements the Zimm model dynamics on top of the static formulation to give bending as well as stretching forces. Another interesting open multiscale problem is to implement this at higher resolutions, with the Rouse model at lower resolutions, in order to design a hybrid model.
To introduce even more realism, we would see individual histones and consider forces between these as in the model of Rosa and Everaers [35] which includes Lennard-Jones and FENE forces between beads. As we approach an atomistic level, it may be interesting to consider a molecular dynamics approach to modelling the DNA filament. Coarser Brownian dynamics models can be estimated from molecular dynamics models either analytically [8] or numerically [9], depending on the complexity of the molecular dynamics model. A variety of structure-based coarse-grained models have been used for chromatin (e.g. [30]), also with transcription factors [42]. Multiscale modelling techniques (e.g. [22] with iterative coarse-graining), as well as adaptive resolution models (e.g. [46] for solvent molecules), have been developed. We expect these studies will connect with polymer-like models at a certain appropriate length and time scale. On top of this, models for the target searching process by proteins such as transcription factors could be improved (for example, by incorporating facilitated diffusion under crowded environment [3]).
The need for developing and analyzing multiscale models of DNA which use one of the above detailed simulation approaches for small parts of the DNA filament is further stimulated by recent experimental results. Chromosome conformation capture (3C)-related techniques, particularly at a genome-wide level using high-throughput sequencing (Hi-C [26]), provide the three-dimensional structure of the chromosomes in an averaged manner. Moreover, recent imaging techniques have enabled us to observe simultaneously the motion and transcription of designated gene loci in living cells [32]. Simulated processes could be compared with such experimental results. Recent Hi-C experiments also revealed fine structures such as loops induced by DNAbinding proteins [39]. To develop more realistic models, information about the binding sites for these proteins may be utilized when we increase the resolution in our scheme.