Multi-resolution polymer Brownian dynamics with hydrodynamic interactions

A polymer model given in terms of beads, interacting through Hookean springs and hydrodynamic forces, is studied. Brownian dynamics description of this bead-spring polymer model is extended to multiple resolutions. Using this multiscale approach, a modeller can efficiently look at different regions of the polymer in different spatial and temporal resolutions with scalings given for the number of beads, statistical segment length and bead radius in order to maintain macro-scale properties of the polymer filament. The Boltzmann distribution of a Gaussian chain for differing statistical segment lengths gives a Langevin equation for the multi-resolution model with a mobility tensor for different bead sizes. Using the pre-averaging approximation, the translational diffusion coefficient is obtained as a function of the inverse of a matrix and then in closed form in the long-chain limit. This is then confirmed with numerical experiments.


I. INTRODUCTION
There have been many studies which use a Brownian dynamics (BD) model for polymers with hydrodynamic interactions as a method to coarse grain complex interactions at the atomic level to study macromolecules in biology and materials science, for example stepping kinetics of kinesin 1 , dynamics for λ-phage DNA [2][3][4] and for polystyrene 5,6 .BD models describe polymers as beads connected by Hookean springs 7 .In this work, we extend the BD modelling framework to allow a polymer molecule to be considered on multiple resolutions, so that the statistical segment length and bead size vary along the polymer as well as the timestep used for simulating BD, i.e. we consider multiple resolutions of both spatial and temporal scales.
The use of a multi-resolution model allows us to look at areas of interest on the polymer in additional detail, while modelling less important areas in less detail and maintaining global properties of the polymer.The main benefit of this form of hybrid modelling is in computational savings compared to modelling the entire domain in high resolution, and has become an increasingly popular technique in recent years as a method to look at complex models on larger spatial scales with more biologically relevant timescales [8][9][10][11] .This is useful, for example, in modelling the interaction between a DNA-binding protein and a DNA filament, where only the area near the protein needs to be modelled in high resolution and other areas can be modelled in a lower resolution.A similar approach to the one taken in this paper has been considered in our previous work in the case where hydrodynamic interactions can be neglected 12 .
a) Electronic mail: edward.rolls@pmb.ox.ac.uk b) Electronic mail: erban@maths.ox.ac.uk; http://people.maths.ox.ac.uk/erban/This work acts as an extension to analytic results for properties of the single-scale polymer dynamical simulations with hydrodynamic interactions.Much of the initial analysis of the single-scale model was done by Kirkwood and Riseman 13 , who introduced elements of the BD model and gave an approximation to the translational diffusion coefficient which uses equilibrium averages for internal configurations.Following this, Zimm 14 found an approximation to the translational diffusion coefficient by pre-averaging the inter-particle distances in the hydrodynamic interaction tensor.Öttinger 15 found a more accurate approximation equivalent to the previous work of Fixman 16 by considering the centre of hydrodynamic resistance and manipulating the Langevin equation, which is the work we will build on to obtain equations for the diffusion coefficient, partly because its formulation has a natural extension to multi-resolution modelling.
In this paper, we start by formulating the model using a statistical physics description and using the Boltzmann distribution to form a mechanical model, similar to how the single-scale model is built by Doi and Edwards 7 .We also define the mobility matrix, which comes from extending the Rotne-Prager-Yamakawa tensor 17,18 to allow for different bead sizes by considering fluid dynamics properties 19 .This provides a Langevin equation for the new model, as well as formulations for the distribution of inter-bead distances, including the root mean squared (rms) end-to-end distance.From this Langevin description, we form two approximations for the translational diffusion coefficient: one as a solution to a matrix inversion problem and another in closed form in the long chain limit.This form of the diffusion is used to give us scaling laws for the statistical segment length, bead radius and the number of beads in each region of the polymer.These scaling laws maintain the global properties of the rms end-to-end distance and diffusion.
A number of algorithms have been proposed for efficient BD simulations with hydrodynamic interactions for the single-scale polymer model in the literature [20][21][22][23] .In Section IV, we adapt the Ermak-McCammon algorithm 24 by varying the timestep along the polymer in order for BD simulations of the Langevin equation to take place.We conclude with illustrative computational results confirming the presented theory in Section V.

II. MULTI-RESOLUTION BEAD-SPRING MODEL
As a model for a polymer, we use an extension to the bead-spring model which has existed for over 60 years 7 .This method of modelling has N beads connected with (N − 1) Hookean springs, neither of which seek to have physical significance as such, but represent a coarsegrained description of the direction the polymer is coiled.Our multi-resolution extension allows different beads to have different sizes and the statistical segment length between adjacent beads to vary.
To form the multi-resolution model, we start with a static description of the chain as an extension of the Gaussian chain model, where much of the analysis follows from the treatment given by Doi and Edwards 7 , but with varying bond lengths.We use this to derive a potential for spring constants using the Boltzmann distribution which in turn is used to form the dynamic model with hydrodynamic interactions given from a multi-resolution extension of the Rotne-Prager-Yamakawa tensor [17][18][19] .

A. Static Description
We consider a bead-spring polymer model with N beads, where the positions of the beads of the chain are given by r n , so that R n = r n+1 − r n , for n = 1, 2, . . .N − 1, are the corresponding bond vectors.In multi-resolution bead-spring model, the distribution ψ n for the n th bond vector R n , n = 1, 2, . . ., N − 1, can vary along the chain: where the statistical segment lengths R 2 n = b 2 n are allowed to vary along the filament.In the special case b = b n for all n = 1, 2, . . ., N − 1, equation (1) reduces to the standard Gaussian chain model 7 .Using (1), the distribution of r mn = r n − r m , for n = m, is given by where we define If we define the rms end-to-end distance µ = R 2 1/2 , where R is the vector from the first to the last bead, then we find Using (1), the conformational distribution function of the chain is given by: Considering a polymer chain at equilibrium, a mechanical model with potential has an identical Boltzmann distribution to equation (4).

B. Dynamic Model Description
Using potential (5), we form a Langevin equation for the dynamic model 25 for n th bead at time t where H mn ∈ R 3×3 is a positive-definite symmetric mobility matrix, dW m ∈ R 3 is a Wiener process and we define B mn ∈ R 3×3 so that where T is the absolute temperature, superscript T (roman font) denotes the transpose of a matrix and k B is the Boltzmann constant.The matrix B mk exists due to the positive-definiteness of H mn , and can be found by performing a decomposition (e.g. the Cholesky decomposition).In choosing the mobility matrix, we would like to vary the bead sizes along the filament, which can be used to ensure that macroscopic properties of the polymer remain constant.Therefore, the bead size σ n becomes a parameter of n th bead.Rolls et al. 12 use the diagonal mobility tensor which works as an extension of the Rouse model: where η is the dynamic viscosity and I ∈ R 3×3 is the identity matrix.One of the purposes of this paper is to extend model (8) to include hydrodynamic interactions.Many models which include hydrodynamic interactions where beads are of equal sizes use the Rotne-Prager-Yamakawa 17,18 tensor.This has been extended by Zuk et al. 19 to allow for beads of different sizes.To formulate it, we denote for beads m and n: where σ m (resp.σ n ) is the radius of of bead m (resp.n) and r mn is the distance between beads.We also denote by r mn the unit vector between beads, i.e.
r mn = r n − r m r mn .
Then the Rotne-Prager-Yamakawa-type mobility tensor is given by 19 which is positive-definite, symmetric and continuous for sufficiently small σ n for all n.It is also incompressible, so that ∇ m • H nm = 0, which simplifies the Langevin equation ( 6) to where B mn ∈ R 3×3 are given by (7) which exists by the positive-definite symmetric property of H mn and interbead force F m is found by differentiating the potential (5)  to get

III. APPROXIMATION OF THE DIFFUSION COEFFICIENT AND IDEAL PARAMETERISATION
We take a similar approach to Öttinger 15 to find an approximation for the diffusion of the multi-resolution model.This is then used in conjunction with knowledge of the distribution of the rms end-to-end vector to inform the scaling for a multiscale simulation, so that properties of interest match up to the 'ground truth' high-resolution model.

A. Diffusion Approximation
To find an approximation to the diffusion for the mobility tensor with hydrodynamic terms included, we use the pre-averaging approximation 7 , introduced by Zimm 14 .Considering near-equilibrium dynamics, we replace the mobility tensor H mn with its equilibrium average H mn eq , using Ψ from equation ( 4): If we assume that σ n + σ n+1 < b n for all n ≤ N − 1, then our equilibrium distribution H mn eq for off-diagonal entries m = n becomes 1 r mn eq I + r mn ⊗ r mn eq where we have used that the distribution of r mn is independent of r mn .Using r mn ⊗ r mn eq = I/3, the second term cancels and we obtain H mn eq = H mn I, where Using (2), we obtain In the single-scale model where b n = b and σ n = σ for all n, equation ( 12) generalises to the equation for the preaveraged tensor in Doi and Edwards 7,14 .Consequently, by pre-averaging equation (10), we find where Following Öttinger 25 , we define the hydrodynamic center of resistance r h by Multiplying equation (13) through by l n and summing over all n, we get Using (11), the first term on the right hand side is zero and the second term is a linear combination of Wiener processes, which is itself a Wiener process with translational diffusion coefficient Using (14) and the definition of l n , we obtain This forms a matrix equation to provide the pre-averaged approximation for the translational diffusion coefficient.

B. Behaviour in the Long Chain Limit
Our analysis in the previous section has used a general multi-resolution model consisting of N beads with sizes σ n , n = 1, 2, . . ., N , connected by N − 1 springs with statistical segment lengths b n , n = 1, 2, . . ., N − 1.In applications to multiscale computations, we are mostly interested in chains which are split into M regions (where M ≪ N ) of constant statistical segment length.In what follows, we will use lower case greek subscripts α (resp.β and γ) to denote regions, while n (resp.m and k) are indices refering to numbers of individual beads and springs along the polymer chain.We assume that the α-th region contains N α springs, with statistical segment length b α , α = 1, 2, . . ., M .Summing over all regions, we have In this section, we simplify equation (15) in the long chain limit, N → ∞, which is taken in such a way that the fraction of springs in each region, N α /(N − 1), remains a constant, i.e. if (N − 1) doubles in size then each individual region also doubles in size.Equation ( 12) defines a function of two integer variables m and n.We will map the discrete function H mn into a continuous function H(x, y) by generalizing the approach of Öttinger 25 and Fixman 16 .Assuming that bead m lies in region α, its continuous approximation in interval [−1, 1] will be defined by where The continuous analogue of a summation of arbitrary function f n over all beads will then be a weighted integral where we define b(x) as a piecewise constant function given by b In addition to (16), we also write for bead n in region β This leeds to a transformation (m, n) → (x, y), which gives the continuous approximation of The definition of the inverse of the matrix H mn (given as Multiplying both sides by b(y), integrating over y and using (18), we obtain where Using the the method of Auer and Gardner 16,26 , we solve equation (19) for φ(z) to obtain To return to the quantity of interest, H −1 nm , we apply the mapping from equation (17) to give where Γ is the gamma function.Substituting in (15), we obtain diffusion constant in the long chain limit

C. Scaling of Parameters
As the use of a bead-spring model is to give a coarsegrained representation of a filament, statistical segment length and bead radius are not physical qualities so we allow these parameters of the model to vary in order to achieve desired statistics of interest for the polymer.In our previous work 12 , the whole-system statistics of interest have been the rms end-to-end distance µ and translational diffusion coefficient D h of a polymer chain.In this paper, we will consider three quantities which multiresolutions simulations should preserve: the end-toend distance µ, diffusion coefficient D h and the strength of hydrodynamic interactions 27 , defined in terms of the parameter h * by where b is the statistical segment length and σ is the bead radius.From a theoretical standpoint, a value for h * ≈ 0.25 minimises the effect of chain length 28 .Similar values can also match experimental results for viscoelastic properties, for example the Flory-Fox parameter can match experimental values 27 for h * ≈ 0.267.In the multi-resolution model, in order to maintain a consistent value for the strength of hydrodynamic interactions, we therefore scale parameters in order to keep h * constant throughout simulations.The multi-resolution polymer simulations will be compared to the 'ground-truth' model, which will be the single-scale model of the polymer in the maximum detail required.In single-scale models, we can modify the whole-system statistics by varying the statistical segment length b, the bead radius σ and the total number of beads N .In single-scale models, we need to select a level of detail for the entire chain as an additional constraint, but by modelling on multiple scales, we instead get to select the resolution of different regions of the polymer, so that only regions of particular interest need to be in the highest level of detail.To parameterize the 'groundtruth' model, we select b to give the desired value for the rms end-to-end distance µ, from equation (3), i.e. b = µ(N −1) −1/2 .Selecting a value for σ is a bit more subtle than for models without hydrodynamic interactions 12 , as the inclusion of the hydrodynamic interactions mean the leading order long chain diffusion approximation in equation ( 20) is independent of σ, i.e. we cannot use D h .We use the strength of hydrodynamic interaction h * to select an appropriate value of σ.For all simulations in this paper, we use σ = b/4.
Once we have defined the 'ground-truth' model, with statistical segment length b, bead radius σ and total beads N , we can seek to define the scalings for the multiresolution model, where different regions coarse-grain the original model to differing extents.We divide the polymer into M regions and assume that the 'ground-truth' chain contains N α consecutive springs in the α th region, for α = 1, 2, . . ., M .Each region of the multi-resolution model has an associated (integer-valued) resolution s α such that s 2 α | N α .Larger values of s α represent coarser regions, and s α = 1 gives the 'ground-truth' model.In the α th region of the multi-resolution model, we have N α springs with statistical segment length b α given by where the definition of the bead radius, σ α , is slightly modified for the boundary beads; scalings (21) apply to beads where both adjacent springs have the same statistical segment length b α .On the boundaries between regions α and α + 1, for α = 1, 2, . . ., M − 1, we take the bead radius to be ((σ 2 α + σ 2 α+1 )/2) 1/2 , and for end beads at the start and end of the polymer we take σ 1 / √ 2 and σ M / √ 2, respectively.By applying scalings (21), equation (3) gives the expected rms end-to-end distance for the filament at equilibrium to be µ = b(N −1) 1/2 , i.e. it is equal to the 'ground-truth' model.The translational diffusion coefficient for the polymer in the long chain limit, equation (20), is also invariant to the number of regions, as well as the size and resolution of each region, and the strength of hydrodynamic interactions is constant along the filament.

IV. SIMULATION METHOD
We solve the Langevin equation for the polymer in equation (10) by using a modified version of the Ermak-McCammon algorithm 24 , for which different regions have different timestep sizes.The key idea for the modified algorithm is to keep track of the behavior of beads modelled with a higher resolution (and with smaller timesteps) to give an average of the hydrodynamic forces exerted on the coarsely modelled beads between the larger timesteps on which they are modelled.
If the 'ground-truth' model uses timestep ∆t, then with the notation of Section III C we define the timestep associated with the α th region as where s α is the (integer-valued) resolution of the α th region, α = 1, 2, . . ., M .A requirement for the resolution value is that for any two regions α 1 and α 2 that either s 3 α1 |s 3 α2 or s 3 α2 |s 3 α1 , to ensure the timesteps of the coarser regions match up to those for the finer regions.We choose this scaling to ensure numerical stability of simulations so that the size of the tension term for a bead is smaller than the statistical segment length with adjacent beads.In the case of a bead lying between two regions, we take the timestep to be the minimum value of the timesteps given by each region.
The 'ground-truth' model updates time at integer multiples of ∆t, i.e. we compute the polymer state at times t = i∆t, where i = 0, 1, 2, 3, . . . .Considering the multiresolution model, we can formally write the update rule (from time i∆t to time (i + 1)∆t) for the n th bead, for n = 1, 2, . . ., N , in region α n as where H nm is the mobility tensor given in equation ( 9), F mn and ρ n are discretized force and noise terms given below, and the function Q is defined for integers i and j by To define discretized force and noise terms, we denote by α n (resp.β n ) the resolution region for the n th bead (resp.spring).Note that for beads with both adjacent springs in the same region we will see α n = β n , however between regions a bead takes the smaller timestep of the adjacent regions, so we may see α n = β n .In the update rule (22), the timestep is incorporated in the force term.
We define the force multiplied by the time step for the m-th bead by This force term is used as a part of a tension term which includes a memory component for larger timesteps, as An explanation of the model running on multiple timesteps with region s1 in twice finer resolution than s2, giving timesteps 8-times smaller.We run a timestep using equation (22) for each bead, so that for every timestep we simulate the beads in s1 and every 8 timesteps for s2.The force terms F m(i∆t) where bead m is in s1 is summed over the small timesteps to give Fmn(i∆t) on the larger timesteps.We use the same concept to find noise terms Xmn(i∆t).
explained in Table I: The random displacement term ρ n has a multivariate Gaussian distribution defined by the moments where we use the maximum as this term is only expressed on the larger of the two timesteps associated with the beads, as laid out below.To calculate the ρ n terms, we an adapted version of the Ermak-McCammon algorithm 24 , so that when we flatten the tensor H mn into a matrix H, we do so by re-ordering the beads so that beads with smaller timesteps have a smaller index than beads with larger timesteps.Having made this adjustment, we can reduce the computational load of the Cholesky decomposition done by the Ermak-McCammon algorithm (an O(N 3 ) calculation 29 ), by only calculating the submatrix made up of the rows and columns of the matrix corresponding to beads which are being calculated on that particular timestep.Therefore, if there are N 0 beads which move on a given timestep, then this gives a 3N 0 × 3N 0 submatrix.
We use the Cholesky decomposition outlined by Ermak and McCammon 24 to get the lower triagonal matrix B, such that H = BB T with elements given by, for n to give noise terms ρ n with the calculation ν nm X mn (i∆t), which are then reordered and formatted to give ρ n in a N 0 × 3 matrix.The random terms Xmn now include a 'memory' similar to the tension terms so that for the terms X m (t) drawn from a Gaussian normal distribution such that X m (t) = 0 and X m (t 1 )X n (t 2 ) = 2k B T δ mn ∆t αn δ(t 1 − t 2 ).

V. SIMULATIONS
We compare the full BD modelling with the dynamics using the pre-averaged tensor H mn eq from equation (12)  in place of H mn in the Ermack-McCammon algorithm given above.Similar to previous papers simulating beadspring models 30 , we choose to simulate with unit paramwhich in our case has k B T = 1 and η = 1.We shall also hold the rms end-to-end distance constant with µ = 1, which we shall maintain by varying the statistical segment length b as a function of the bead number N as appropriate given the scalings explained in Section III C. To ensure numerical stability of simulations, we found ∆t = 10 −2 b 3 η/k B T to be a good value to use.
In order to study the diffusion coefficient in simulations, we need to extend its definition for the simulations of the multi-resolution model.We define the mass of the polymer Ω as where σ n is the radius of the n th bead, n = 1, 2, . . ., N .Using scalings (21), we observe that Ω is invariant to our choice of resolutions in the multi-resolution scheme.This allows us to define the centre of mass of a polymer r G at time t as where r n (t) is the position of the n th bead at time t.With this we retain the definition of the translational diffusion coefficient as Note that for the single scale simulation this reduces to the standard definition for the translational diffusion coefficient.
In this section we compare the translational diffusion coefficient approximations given in the form of an inverse matrix in equation (15) as well as the long chain limit approximation in equation (20) to BD simulations for both the full mobility matrix as well as the pre-averaged approximation, as well as the rms end-to-end distance, The difference between the matrix formulation for the diffusion coefficient (15) and the long chain limit (20) as N gets large, for both the single scale-system (red dotted line) and the multiscale simulation, which has the middle 25% of the polymer in high resolution (blue dashed line).Parameters are given in Section V. Note that the multi-scale system needs N ≥ 32 by construction for the coarse grained particles to be placed.
which has expected value 1 in all simulations by design.We run the BD simulations using both a pre-averaged and a non-pre-averaged mobility tensor for a total of 10 4 timesteps (in the case of multi-resolution simulations, this refers to timesteps associated with the higher resolution beads), where each result is given as an average over 500 runs, and contrast this to the diffusion approximations.
In the results tables we include 95% confidence intervals for the translational diffusion coefficient and the rms end-to-end distances.For the end-to-end distance we calculate the 95% confidence interval for R 2 and take the square root for both lower and upper bounds to give a range (note that this is not symmetric about the rms value).
We consider three illustrative examples.The first one is a single-scale system, so that M = 1 and s 1 = 1.We use it as a control to compare the other simulations to.Our results are presented in Figure 1 and Table II.As we can see in the table, the matrix formulation D MF is covered within the 95% confidence interval all except one of the pre-averaged and non-pre-averaged values of N .The analytic value of µ = 1 fits in the confidence interval for all simulations.
The next system to consider is one in which the middle 25% of the polymer is in high resolution, while the remainder is coarse-grained by a factor of 2. Therefore we define M = 3 with s 1 = 2, s 2 = 1, s 3 = 2, and N 1 = 3N/8, N 2 = N/4, N 3 = 3N/8.The diffusion and the rms end-to-end distance of this polymer is given in Table III.The matrix formulation for the diffusion D MF is contained in the confidence intervals for all values of  II.Results for diffusion and the rms end-to-end distance in single-scale simulations.The subscript MF represents the matrix formulation, PA is for pre-averaged and NPA is non-pre-averaged.The 95% confidence intervals for µPA and µNPA are given by CIµ PA and CIµ NPA , respectively.
N , and the value of µ = 1 is contained in the confidence interval for all values of N .The convergence of matrix formulation to the long chain limit is shown in Figure 1.
The final system considered had an 8-times resolution increase in the middle 10% of the the polymer.This uses M = 3 with s 1 = 8, s 2 = 1, s 3 = 8, and N 1 = 9N/20, N 2 = N/10, N 3 = 9N/20.We perform the simulation for N = 1280, from which we can report D = 0.189, µ = 1.04 for the pre-averaged case and D = 0.180, µ = 1.01 where we do not use pre-averaging.The matrix formulation for the translational diffusion coefficient gives a value D MF = 0.190.The end-to-end distance in the pre-averaged case narrowly falls out of the 95% confidence interval, but the other three stastics lie within this range.
As can be seen from the simulations, there is good agreement both between the pre-averaged and non preaveraged tensors, as well as between the diffusion approximations and the results from the simulations.In total out of 44 observations, we had two fall outside of the 95% confidence interval.The overall goal of doing this coarse-graining is to improve the speed of simulations.In Figure 2, we compared the timings between the single-scale and multi-scale models for identical parameters as were used to produce Table III.There is a pronounced difference between the multi-scale model and single-scale model without pre-averaging, most of which comes from having to use the Cholesky decomposition on smaller matrices, while the smaller difference in the model with pre-averaging comes from updating fewer beads in each timestep.

VI. DISCUSSION
In this paper we have extended the bead-spring model for a polymer including hydrodynamic interactions to a multi-resolution model in order to gain computational efficiency for BD modelling.By considering a multiresolution Gaussian chain model, we have utilised the Boltzmann distribution in order to form a Langevin equation for the multi-resolution model.From this we used a similar approach to Öttinger 15 in order to derive an integral equation for the diffusion of the polymer using the pre-averaging approximation, which was then manipulated to find a closed form equation for the diffusion in the long chain limit.This gave scaling laws for key parameters of the polymer at different scales.The developed multi-resolution approach keeps the rms end-to-end distance and the diffusion of the polymer invariant to the choices of how we split the polymer up into different resolutions.These scaling laws have been then supported by illustrative simulations, which used an adapted version of the Ermak-McCammon algorithm 24 .
This work has been looking at extensions to a polymer model in a theta solvent, which is not the most general state that a polymer can exist in.To extend this model further, it is of interest to include excluded volume forces to allow for the study of a good solvent, of which there has been much analytical work to derive terms for the rms end-to-end distance 31,32 in the single-scale model.Another possible extension is to look at including additional forces between monomers to get more realistic spring forces which are used in many recent studies of polymers 3,33,34 , for example to form a wormlike chain model 35,36 .
There have also been many recent developments in the algorithms which are used to study polymers with hydro- dynamic interactions [20][21][22][23] , which could improve the computational efficiency of the multi-resolution modelling even further than the Ermack-McCammon algorithm 24 , which has been used here as a demonstration of the scalings.
FIG.1.The difference between the matrix formulation for the diffusion coefficient (15) and the long chain limit(20) as N gets large, for both the single scale-system (red dotted line) and the multiscale simulation, which has the middle 25% of the polymer in high resolution (blue dashed line).Parameters are given in Section V. Note that the multi-scale system needs N ≥ 32 by construction for the coarse grained particles to be placed.

FIG. 2 .
FIG. 2.The CPU times to simulate identical systems using four different algorithms: the multi-scale model using preaveraging (red dotted line), the single-scale model with preaveraging (purple dashed line), the multi-scale model without pre-averaging (blue dot-dashed line) and the single-scale model without pre-averaging (yellow solid line).