model with dynamic Monte-Carlo
simulations
THESIS
submitted in partial fulfillment of the requirements for the degree of
MASTER OFSCIENCE
in PHYSICS
Author : Yahel Houston
Student ID : s2547651
Supervisor : Helmut Schiessel
Second corrector : John van Noort
Building up a nucleosome sliding
model with dynamic Monte-Carlo
simulations
Yahel Houston
Instituut-Lorentz, Universiteit Leiden P.O. Box 9500, 2300 RA Leiden, The Netherlands
March 25, 2021
Abstract
The properties of DNA can be segmented into two main categories; chem-ical and mechanchem-ical. Whilst the chemchem-ical nature of DNA is important to its functionality, we here focus on the mechanical nature of DNA, specif-ically relating to its geometrical properties. The mechanical properties of DNA are strongly linked to the underlying sequence of the DNA which in turn affects its ability to wrap around the fundamental packaging unit, the nucleosome. In this thesis, we explore the dynamics of a simplified model of DNA that is wrapped around a nuceleosome. Using dynamic Monte-Carlo simulations, we concentrate on calculating two main physical quan-tites; the diffusion constant and the mobility. This thesis also represents building up a model from rather basic principals and adding more com-plexity after each step, this allows us to develop a deeper understanding of the dynamics of the system.
1 Introduction 1 2 Model 3 3 Methods 9 3.1 Monte-Carlo 9 3.2 Measurable Quantities 10 4 Results 13 4.1 Initial Setup 13
4.2 Diffusion and drift in the absence of an external potential 14
4.3 Diffusion in a periodic potential 15
Chapter
1
Introduction
Genetic information used by all living organisms is carried within a molecule called Deoxyribonucleic acid (DNA) and is made up of 4 base pairs: Ade-nine (A), GuaAde-nine (G), Thymine (T), Cytosine (C). These pair up as A-T and G-C, forming nucleotide pairs which bind together to form a DNA sequence. DNA as a chain forms a double helix from two right-handed chains of nucleotides. The length of a DNA molecule is much larger than the cells in which they live in and as such must be wrapped up into multi-layered structures. The most fundamental packaging unit is the nucleo-some, where ∼ 1 & 3/4 turns of DNA are wrapped around an octamer of histone proteins [1], this can be imagined like a string wrapped around a spool. The intrinsic mechanical structure of the DNA is encoded into the DNA sequence itself [2], meaning some DNA sequences have a higher preference to be wrapped around the nucleosome which reflects on the sequences ability to be bent sharply. This therefore means, the energy re-quired for bending the DNA around the nucleosome is dependent on how it is bent and which sequence is bent.
Previously, Monte-Carlo simulations were used to investigate whether mechanical properties of DNA can be used to predict nucleosome position rules; where nucleosomes are located with respect to a DNA sequence [3]. The model used is accurate but rather complex, thus, understanding the dynamics of the model can be hard. Here, the idea of ”mutation Monte-Carlo” algorithms are introduced, where one can manipulate the sequence position, but also the sequence itself. A previous post doc of Prof. Schies-sel also used Monte Carlo simulations to predict the bending energy of DNA given its sequence and its position. The bending and unbending of DNA is important in processes such as DNA transcription, which is the initial stages into copying a DNA onto mRNA
In this thesis, computer simulations are used to investigate the basic physics of a mock DNA strand and we slowly build up a more compli-cated model representing the dynamics of DNA purely represented by its intrinsic geometry. We begin by explaining the model and its origins and then build more complexity into the system, along the way confirm-ing certain physical behaviours such as the diffusion of the nucleosome along DNA. This thesis therefore provides a solid basis for further work into investigating the dynamics of DNA based upon its intrinsic geom-etry. The end goal is to investigate what happens when one end of the DNA is pulled whilst allowing for sequence mutation, DNA is forced to slide in these situations. Due to the intrinsic geometry of DNA, this allows us to use a Frenkel-Kontorova like model from which we can calculate the friction constant from the Einstein relation.
In the next chapter, we introduce the model used and the simplifica-tions made from previous work. Next, we introduce the use of Monte-Carlo simulations and discuss what physical quantities we aim to calcu-late through the simulations. Finally, results from the simulations are pre-sented with concluding remarks and a short discussion on future work.
Chapter
2
Model
As stated earlier, we wish to start from a very basic model of DNA wrapped around a nucleosome, however we must first discuss how we arrive to thinking about this basic model. The model we use is first based off of the work by Eslami-Mossallam et al [3] and the idea of ”nucelosome po-sitioning codes”. This is the idea that desired mechanical properties of DNA is written into their sequence. Here, Monte-Carlo simulations are used to force constraints onto the DNA-Nucleosome crystal, representing how DNA is constrained on a real nucleosome. We stress here, that these calculations are dependent on the underlying sequence and the sequences mechanical properties determine how it behaves. In this paper, DNA is modelled using rigid base pairs which describes the DNA double helix as a combination of the position and orientation of its base pairs. The double helix is 147 base pairs long and it wrapped lefthandedly ∼1 & 3/4 times around a nucleosome core, which consists of eight core histone (histone octamer). Here, DNA attaches itself to the histone octamer through its mi-nor groves, once every 10 basepairs, both of these can be seen in fig.2.1 a,b. To simplify the complexity of DNA, the model treats base pairs as rigid plates. This model has six degrees of freedom, 3 translational (shift, slide, rise) and 3 rotational (tilt, roll, twist) per base pair step. For a given base pair, these parameters have intrinsic values which are based off crystallog-raphy of DNA protein complexes and atomistic simulations. The shape of the DNA is accounted by modelling binding sites between the histone and the DNA and allowing for some structural relaxation of the DNA. All the above degrees of freedom contributed some energy to the system. In this model, it is assumed that the energy deformation between the actual state and the intrinsic shape is quadratic in nature (Hookes Law):
E = 1
2(q−qeq)
TQ(q−q
eq) (2.1)
Here qeq is the intrinsic state vector of the six degrees of freedom and Q is a stiffness matrix of a given base pair, both of which are based off of chemistry and whose values can be found in literature [4] [5]. Whilst the model used is very realistic and is able to accurately predict nucleosome positioning rules, the model is too complex to get a full understanding of the underlying reasons as to why these rules come about. A simpler model was proposed by Zuiddam et al [6], where the DNA is forced into an ideal superhelix, which is representative of how DNA is bent on top of a nucleosome, a schematical representation of this is shown in fig. 2.1a. The actual configuration of DNA can be defined by creating a superhelix using super-helical coordinates. The energetic costs are quadratic in the deviations of the ideal superhelix from the intrinsically prefered geometry. The superhelix can be parameterised as follows; we can consider that the DNA that is wrapped around the nucleosome is L = 147 base pairs long and is wrapped left handedly α = 1.84 times. The radius is 4.19nm with a pitch P of 2.59nm which means that the superhelix has an effective radius of Re f f = p R2+ (P/2π)2. We can define the position along the superhelix as
¯r(s) = [Rcos(s/Re f f), Rsin(s/Re f f),−(P/2πRe f f)s] (2.2) The rotational orientation of a base pair plate with respect to the origin can be described by 3 orthonormal vectors (see fig. 2.1c). The double helix shape can be defined from the orthonormal vectors so the double helix twists around the super helix:
[ ˆx(p), ˆy(p), ˆz(p)] = [ ˆn(p)cos(θ p + φ)−ˆb(p)sin(θp + φ),
ˆn(p)sin(θ p + φ)−ˆb(p)cos(θp + φ),
ˆt(p)]
(2.3)
where ˆn, ˆb, ˆt are the set of frenet serret vectors with p(s) = s(L−1)/2πRe f fα. The value of p can act as a base pair position index, with p = 1, ..., n being the nthrest position of the base pair. Here the values of θ, φ are set to repre-sent how much the double helix is twisted. We set θ = 2π/10 to represent that the minor groove faces towards the octamer every 10 basepairs, this can be seen in fig. 2.1. The phase φ = −147π/10 is set such that the cen-tral position of the dinucleotide steps (73,74) correspond to the position
5
(a) (b)
(c)
Figure 2.1:(a)Representation of DNA as a chain of rigid plates. Here rigid plates are fixed, whereas in our model we allow them to move (figure from [6]). (b) A simplified view of DNA wrapped around a histone core, with minor (green) and major (red) grooves shown. (c) Orthogonal vectors of the rigid base pairs.
of maximal roll, so that the major grove faces the histone octamer. Using the coordinate system (eq. 2.3) and the superhelical coordinate p, we can define the roll, tilt and twist angles along the superhelical axis:
[qrollp , qtiltp , qtwistp ] = [γcos(2π p/10−147π/10),
γsin(2π p/10−147π/10), qtwist]
(2.4) where γ ≈ 0.0796rad which comes from the constant curvature of the su-perhelix and qtwist ≈2π/10.17rad which is slightly smaller than the angle of θ = 2π/10, which reflects the fact that the path is superhelical.
The energy of a dinucleotide base pair at a step a, b is given as a sum of all the energies of these parameters (ignoring twist since its energetic contributions are much less important than roll and tilt):
where:
Erollp = 1 2Q
roll(a, b)[q
roll p− ¯qroll(a, b)]2 (2.6) Q are basepair step dependent stiffness’s and q’s are the intrinsic values. A similar term is held for the tilt energy. In addition to ignoring the twist, we also assume that the DNA is unshearable, which helps make the model more analytically solvable.
This thesis begins to work with the ideas from this paper, namely us-ing the same parameterisation as above but also usus-ing its intrinsic shape to act as an energy landscape to define the DNA’s motion. We start by using the following two parameters: the roll and rise. We use the rise to reflect the position of a base pair but also the fact that in DNA, base pairs do in-fact move. This is one of the key differences to the aforementioned paper where base pairs are fixed in position, whereas we allow them to move. Here we use a 1D version of eq.(2.1) and assume that the force between base pairs act like springs, where q acts as the base pairs position (old and new) and k acts as a spring constant, which has a value for every dinu-cleotide. The motivation of allowing base pairs to move has a biological nature since nucleosomes can in fact slide, in a phenomenon referred to as nucleosome sliding [7]. As explained earlier, DNA has an intrinsic shape, which is sequence dependent, therefore when we allow a nucleosome to move, we expect this sequence dependence to become large and appar-ent. This could be apparent if we see a resistance to movement since we are moving the DNA into an un-preferred state and so there is an energy penalty.
How can we start to think about this motion? We can first imagine choosing a minor and a major groove where the DNA is attached to the nucleosome core. If we force the DNA to move, we have something akin to corkscrew motion and the sequence of DNA is pushed through. Here, the base pair originally located at the minor grove may be more energet-ically favourable as the intrinsic shape fits the required geometry. If we push this base pair 5bp further, this base pair step is at a position where there is a major grove, which may be a geometry where it is energetically more costly. We have therefore moved the base pair out of its preferential geometry, whose energetic cost is sinusoidal in nature. This comes from eq.(2.6) by inserting the expression for qrollp from eq.(2.4). We can therefore represent this as the energy landscape that the base pair feels if we force it to move.
Considering both the rise and the roll, a system like this can be mod-elled using the Frenkel-Kontorova (F-K) model. This model was first
in-7
troduced to describe the structure and dynamics of a crystal lattice near a defect, but have also been used in numerous areas such as crystal sur-faces, Josephson junction and biomolecules [8][9]. It describes a chain of particles with nearest neighbour interactions subject to a periodic poten-tial. Although being rather simple, the F-K model can exhibit rich and complex dynamics, such as the formation of kinks [10].
In a very simple way, we can begin using the Hamiltonian as:
H= Usub+ Uint (2.7) Uint = 1 2k( ¯x−x) 2 (2.8) Usub = A 2(1−cos(2πx/λ)) (2.9)
Where A is the periodic potential height, λ is the wavelength, x reflects the position and ¯x is the intrinsic position. Here is how we can link this to the nucleosome model; we consider each particle in the F-K model as a base pair step. Usub corresponds to the periodic potential that is felt by the base pair as we move the DNA around the nucleosome in a corkscrew motion coming from eq.(2.6). Uint corresoinds to the elastic energy coming from the rise in the DNA, which is the distance between two base pairs. One can also easily implement a force into this system by another external potential: U = −f ×x, which has the effect that the chain of springs and particles are stretched out.
Now, in this thesis we work with a very simplified version, meaning we do not consider parameters extracted from real DNA. Instead, we work with some standard parameters that easily allow us to check whether the physics of our system matches theoretical expectations whose numbers that are easier to work with. Here, we concern ourselves mainly with the following: Thermal Energy, Diffusion constant and the Einstein relation. These quantities offer nice indicators of the dynamics of our simulations and whether we are missing something in our understanding of the sys-tem or if there is something wrong with the code itself. For example, when implementing a periodic potential, we know in the lower limit (A →0) the dynamics should replicate regular diffusion.
The end goal is twofold: 1. Produce a model where we start from some-thing simple and when adding more complexity, we take a step back and fully understand all of the dynamics. This allows us to possibly reach a point where we build to a model as complex as [6], [3] but we have a bet-ter understanding of the underlying dynamics. 2. As briefly mentioned
earlier, we are also interested in the sequence dependence on nucleosome sliding. Using this model, we can force nucleosome sliding by adding a pulling force. Using Monte-Carlo simulations, we can also allow the se-quence to simultaneously mutate, investigating different rates of mutation and their affects on the dynamics. This could mean that we get different types of apparent motion of the DNA sliding, sometimes rough, ragged and other times smoother. This can be characterised through the mobility µor the friction constant. We choose Monte-Carlo simulations rather than molecular dynamic simulations for this reason; implementing simultane-ous sequence mutation and position moves is much easier.
Chapter
3
Methods
3.1
Monte-Carlo
Monte-Carlo (MC) simulations represent a large class of algorithms that use random numbers to solve deterministic problems. To investigate the properties of our model we use the Metropolis-Hastings algorithm which works as follows;
• Initialise the system as all particles at equilibrium distance
• We select a random particle and propose a change in its position. • Accept the proposed change if dE < 0, or if dE > 0 accept with
probability given by the Boltzmann factor, p = exp(−dE/kBT), where dE is the change in energy of the system
The change in energy dE is a change in energy of the Hamiltonian (eq.2.7). The energy change is calculated by proposing a change in po-sition (a move) by adding or subtracting a distance labelled ”step size”. From this, we can calculate the energy of the current position and the new position. The move is affected by the energy cause by Hookes law and the position of the particle along the periodic potential - is the spring moving towards or away from the equilibrium position of the spring and is the particle moving up or down the periodic potential. Initial settings of the simulations are described in a section below.
To correctly account for the dynamics of the system we need to care-fully consider one small detail; time within the simulations. When run-ning a regular MC simulation, we define one MC step as N times one iteration step. This means that the simulation time is therefore N ×”MC
time”. However, when calculating a dynamical quantity such as the dif-fusion constant D, we need to re-scale this time. This is because, if one takes in face value this time quantity, it contains both the number of ac-cepted and rejected moves. This is important when we have to consider that particles can be connected by springs and a move by the particle has an associated energy change. We need to retain the fact that our particles only ”move” when a move is accepted, thus, we should re-scale the time value by a quantity: 1/(1−R), where R is the rejection ratio, defined as (steps rejected)/(total number of steps). This means that, if all moves are accepted then we do not re-scale the time, or if half the moves are rejected then we re-scale the time by 2. Algorithms like this are referred to as ”dy-namic Monte Carlo”, however methods including time re-scaling are not commonly found in literature. Equivalence between Brownian dynamics and the dynamic Monte Carlo with the re-scaling by the rejection ratio has been shown [11]. We aim to have a rejection ratio of around 10-60%, whilst this is a rough number, it is often stated as the rejection range for MC simulations.
3.2
Measurable Quantities
Here, we briefly describe some of the quantities that we are interested in and how we calculate them.
When first treating the system as a chain of particles and springs, one can use the thermal energy as a first sanity check where the thermal en-ergy of the system is given by the equipartition theorem: 1/2kbT(N − 1) = 1/2kx2 where N is the number of particles (so N-1 is the number of springs). This can be achieved by summing up the potential energy of spring system over all the time steps. To correct for the timing here, we also re-scale such that in the simulations, the thermal energy is calculated by∑ Epot×timestep/(time×N) where the sum is the potential energy over all particles over all time. In other limits presented (such as no spring con-stant and no periodic potential) this offers a nice check about the energy of the spring system.
The next quantity we concern ourselves with is the diffusion constant for a single particle which can move in a random direction with increments of size a within a time interval δ:
D = a 2
2δ (3.1)
3.2 Measurable Quantities 11
displacement vs time, with the gradient being 2D. This is a quantity that allows us to assess the motion of our simulation in the absence of any forces, i.e similar to its thermal/ free motion. We now briefly deal with the impact of the spring constant k on a chain of particles and measuring the diffusion constant from its centre of mass. Consider a chain of 2 particles with spring constant k and first consider the limit k → 0. Here, we es-sentially have two independent particles diffusing and during a MC step, both particles have the chance to move. On average, half of the moves are the particles moving together (centre of mass does not change) and half of the moves the particles move apart (centre of mass is changed by a). Therefore the diffusion constant is half of that for a single particle. This can be generalised as follows: DN = D/N where N is the number of par-ticles. However, in the limit of k 6= 0, we need to consider the fact that moves can be rejected since the change in energy dE > 0 is only accepted with the probability of the Boltzmann factor. Therefore, as described in section 3.1 we introduce the time scaling of 1/(1−R).
We are also interested in the diffusion of a particle in the presence of a periodic potential. In this, we expect different regimes of particle dif-fusion, dependent on 2 factors: periodic potential height and spring stiff-ness, both parameters of these can be easily adjusted and tested. The fol-lowing expression is introduced in the context of Chromatin Dynamics by Schiessel [12] to study the sliding of nucleosomes along DNA via ”Twist Defects”. The full derivation comes from Dietrich [13], where the mobility is calculated for a (cosine) periodic potential, this can then be related to the diffusion constant by the Einstein Relation (see below):
D0 = DIo−2 A 2kBT (3.2) where A is the barrier height and Iois the modified Bessel function. In the limit of a lower potential barrier height, we expect to recover the normal relation (the Bessel function tends to 1). However, for different spring constants with N particles we must think a little carefully and consider to cases k → 0 and k >> 1. For k → 0, we essentially have the case where there are N independent particles. Thus for this case, we can consider that D00 = D0/N. For the case of k >> 1, we have to consider that the particles are strongly bound, where you can consider them acting as one ”rod” together. If we imagine now an attempt for this rod to move over one of the barriers, essentially all of the individual beads move over at the same time, thus we have an effective barrier height of N×A. We perform simulations within both of these limits.
particles ability to move (for example with a charge, the ability to move through a solid). This can be calculated via two methods: 1. Einstein relation D = µkBT, and in response to a force 2. µ = Vd/F. Here, we can use the diffusion constant for N particles and the Einstein Relation as a theoretical expectation. This can be compared to simulations of pulling on a string and calculating µ from Force vs Velocity curves.
Chapter
4
Results
4.1
Initial Setup
First, we describe our initial experimental setup. As stated earlier since we are using the work in this thesis as a build-up, we use rather simple parameters to start with. If we first consider a simple 2 particle system connected by one spring we consider some of the following parameters: the value of the spring constant, equilibrium distance between particles, the temperature, the step size of the Monte-Carlo simulation. These are shown in the table below.
Quantity Value
Spring Constant 10 pN/nm Equilibrium Distance 1nm
Temperature 300K
Step Size 0.02nm
Using these values, we expect the thermal energy to be 1/2kBT ∼ 2.07×10−21J which we can multiply by N-1 to account for the number of springs we have. We choose the step size (as in when we propose a move how far can it move) in order to have an acceptance ratio in the range of 10-60%. This is due to the fact that when calculating the change in energy from one position to another, it is dependant on the step size. There is indeed a trade-off here, especially when one has to consider an external potential. If the step size is too large, the dE can be too large and all moves are rejected and on the contrary, if the step size is too small, then all moves are accepted. We choose the equilibrium distance to be on the scale of nanometres since the distance between base pairs is on the order of
0.4nm. These parameters can be easily adjusted to be based on parameters of DNA wrapped around a nucleosome, however we note that in order to achieve efficient MC simulations, we must adjust the step size to be on the correct order and to achieve an acceptable acceptance ratio. It is also im-portant to emphasize that in most of the results shown below, we mostly care about the accuracy of the results, in that our simulations can accu-rately replicate certain physical quantities, meaning we are certain that we understand what our simulations are doing.
4.2
Diffusion and drift in the absence of an
ex-ternal potential
Diffusion is simulated via the Monte-Carlo (MC) algorithm outlined in section 3.1, the results for the diffusion of N=4 and N=4 are shown in fig. 4.1. Here we highlight the difference between the ”regular” MC algorithm and a dynamic MC algorithm. This can be studied by considering a chain of N = 4 particles with k = 0 pN/nm, fig. 4.1a and k = 10 pN/nm, fig. 4.1b. For k = 0 pN/nm, all moves are accepted in the MC algorithm since the particles are free to move without any energy change in its proposed move, as such the re-scaling has no effect which can be seen as the gra-dients having the same value. For k = 10 pN/nm, proposed moves can be rejected, so we must introduce the time re-scaling. The difference be-tween the with & without re-scaling can be seen clearly in fig. 4.1b where the gradient calculated without the re-scaling is out of bounds of the error. With the re-scaling of 1/(1−R), the calculated diffusion constants agree well and are within the bounds of error.
We now investigate the mobility of our physical system, which as de-scribed earlier give a description of the ease of an object to move. Here we perform simulations simultaneously. In one simulation, we let a chain of N particles (we choose N = 2) to diffuse freely, from which we use the Einstein relation to calculate the mobility from the diffusion constant -µ = D/(kBT). This serves as our value to compare against, since we know our accuracy of the diffusion constant is good. We also then perform a pulling experiment on the chain of N = 2 particles, where we pull on one end of a chain of particles with the other end free to move. From this, the drift velocity can be calculated by the average position against time, given from a straight line fit. The drift velocity is then plotted against the force, where the mobility is then given as the gradient of this plot. An example of this in fig. 4.2 for a chain of 2 particles. The calculated % difference
4.3 Diffusion in a periodic potential 15
(a) (b)
Figure 4.1: (a)N=4 particles are free to diffuse independently (k=0). Calculated Diffusion constants are D = (5.1±0.1)×10−5nm2/ts for both with and without
re-scaling. The theoretical diffusion constant (eq.3.1) is D = 5×10−5nm2/ts so is within the bounds of error. (b) Diffusion for a chain of N=4 particles with k = 10. The calculated diffusion constants are: D = (4.9±0.2)×10−5 nm2/ts with re-scaling and D = (2.73±0.08)×10−5 nm2/ts without re-scaling. The theoretical diffusion constant is D = 5×10−5 nm2/ts, so with the re-scaling the calculated value is within the bounds of error as the theoretical value. In the units ts stands for time step
between the expected value and that of the force pulling lie within the bounds of error.
Unfortunately, pulling simulations with the inclusion of a periodic po-tential were not obtained due to a lack of time. However, this is where the next steps of the thesis should leave off from since this is where we would start to replicate nucleosome sliding with the inclusion of real parameters.
4.3
Diffusion in a periodic potential
Next we look at adding the periodic potential, which is defined from the quantity Usub, eq.(2.9). Whilst this is not the exact form that we expect (eq.2.6), it can be easily adjusted to take any form or shape we wish and for the moment tried to keep it as simple as possible since, motion and dynamics in Frenkel-Kontorova models can be rather complex. First, we check the diffusion constant against barrier height A. As stated earlier, in the low limit A →0 we expect D = a2/2δ , in the larger limit we expect D to take the form of eq.(3.2).
We note here that there are in fact two different time scales at play for the periodic potential. There is a short scale diffusion, where the particle is
Figure 4.2: Forced pulling on a chain of 2 particles connected by a spring. Sim-ulations were ran for 500 time steps, in steps of 5 with 400 repeats per time step. Forces run from 0pN to 40pN in steps of 2pN. The calculated mobility µ from the simulations is µ = 0.0095±0.0001 nm2/N ts, with a theoretical value of µ = 0.0096 nm2/N ts.
4.3 Diffusion in a periodic potential 17
Figure 4.3: Diffusion constant simulated whilst varying the periodic potential amplitude with 1 particle. Simulations were ran for 800 time steps in steps of 2 with 200 repeats. These also include the standard parameters shown in table 1.
trapped within the well of periodic potential and a longer scale diffusion (hopping diffusion) when a particle has overcome a barrier. These can be clearly seen in fig. A1 of the appendix and are indicated by the two different regions of MSD vs time. These effects become more apparent when the potential barrier is larger. We take the diffusion constant from the hopping diffusion, which is from the last half of the data of the position against time.
Fig. 4.3 shows the theoretical diffusion constant (eq.3.2) against the simulated diffusion over the periodic potential with varying the periodic potential amplitude. Our results show a good agreement between the ex-pected diffusion constant and the simulated ones. There are slight devia-tions, mainly in the middle region of the plot, this could be due to errors arising from the goodness of fit from the MSD vs time.
Variations in the wavelength of the periodic potential were also inves-tigated. Since there is no explicit expression for the wavelength in eq(3.2), we expect no variation in the simulated diffusion constant. This can be seen in fig. 4.4, where although the values fluctuate quite a lot are still within the bounds of error. The reason for this variation is the lack of re-peats and means there is greater statistical error. The reason there are less repeats is because of the increased number of MC steps used, meaning
Figure 4.4: Varying the periodic potential wavelength. These simulations were ran for 75000 timesteps in steps of 200 with 100 repeats.
that the real time running of the simulations was also larger.
Next, we investigate how the diffusion constant changes with the spring constant k whilst they diffuse over a periodic potential, with a potential barrier with energy equivalent to roughly the thermal energy and wave-length of the equilibrium distance. We choose the step wave-length of the sim-ulations to be much smaller than the wavelength (SL = 0.02nm compared to λ = 1nm). This means that the particles feel the potential and the change in energy is not too large that every move is rejected. Here, we believe there again are two regimes for the diffusion constant that dif-fer slightly. For very small values of k, we can imagine that a chain of beads almost move independently, in that they are very weakly coupled. If we look back to eq.(3.1), we can consider how diffusion constant is af-fected by the number of particles N. If we have N independent particles, and calculate the diffusion constant (from their centre of mass), that we can simply divide the diffusion constant eq.(3.2) by N such that we have D0 = D/N× Io−2(A/(2kBT)). However, in the regime where k is large (k >> 1) we can estimate the N particles as acting almost as one. There-fore, if they need to overcome a potential barrier, we can imagine that they do this as one ”rod” rather than going over independently. If, say one par-ticle of the chain moved over a period independently, then we get what is called a kink, which for a large value of k will result in a highly frustrated
4.3 Diffusion in a periodic potential 19
system (this is unlikely to be accepted moves in a MC simulation). This can be expressed explicitly with N included:
D0 = D NI −2 o N A 2kBT (4.1) A brief investigation of the variation of the Bessel Functions shown above as a function of the number of particles and the barrier height are shown in fig. A2 & A3 of the appendix. This allows us to see how these quanti-ties diverge and the differences between the regimes of diffusion that we simulate (without periodic potential and periodic potential for a chain of particles with k <<1, k >>1). This is important as they smooth out small discrepancies that could be missed.
For both high and low values of k, simulations match theoretical pre-dictions well, shown in fig. 4.5. However, towards very large values of k, the simulations begin to fall off quickly. This is due in part to the MC simulations themselves; the values of k become so large that any incre-mental changes in the position and thus the change in energy dE also be-come very large. This means that in the MC simulation, most, if not all moves are rejected rather than accepted, thus the particles remain stuck. The solution to this could be to change the step size to much smaller val-ues, so that the energy change dE remains somewhat reasonable in the Boltzmann distribution. However, the trade-off for this would be that the simulations must be ran for longer in order for the particles to move a meaningful distance (say 1 equilibrium distance). These simulations are already rather computationally expensive (due mostly to calculating ex-ponentials) and thus the real time running of the simulations would also become very long for a home laptop. A work around to this could be to store an access pre-calculated values so that the computer does not have to calculate an exponential each iteration.
(a)N = 2 (b)N = 4
(c)N = 2,3,4
Figure 4.5: Varying the spring constant for a constant periodic potential ampli-tude. (a), (b) are figures for N =2,4 particles for larger k values with A = 1×
Chapter
5
Conclusion
In the work presented, we have shown that Monte-Carlo simulations of a simplified rigid-base pair model match theoretical expectations of quan-tities such as the diffusion constant and the mobility within a degree of error. We start with a relatively simple 1D chain of springs and add more complexity to the system by pulling on it and adding a periodic potential. Adding both of these means we can replicate the full model as presented by Zuiddam and allows us to have a better understanding of the system dynamics. By breaking down the model and building it up in smaller steps, it allows you to realise certain behaviours and characteristics of the model, for example, the diffusion over a periodic potential gives you two distinct diffusion regimes.
For further work, the starting point would be to calculate the mobility for a chain of connected springs over a periodic potential, which can be calculated by pulling on the chain with a varying force. This is where even more complex and intricate behaviour could start to appear, such as the presence of kinks. Here it would also be interesting to look at the effect of the size of the periodic potential and the value of the spring constant. The next step would be to incorporate the real parameters found in DNA such as the values for rise, spring constants between base pairs and periodic potential barrier heights.
Appendix
Figure A1: Mean Squared Displacement for a single particle in a periodic poten-tial with A = 5×10−20J (where the thermal energy is∼2×10−21J). We can see that a straight line fit (which is normally used to calculate the Diffusion constant) does not work here since there are two distinct diffusion regimes. Instead, the second half of the data is used for the straight line fit. However, we see that the variation in the mean square distance (MSD) varies greatly for the larger time values which leads to an increased error (compare this to fig.4.1).
23
(a) (b)
(c)
Figure A2: Diffusion constant as a function of the number of particles. In these figures B1 is eq.4.1 and B2 is D0 = D/N Io−2(A/(2kBT)). Both of these refer to
dif-fusion in a periodic potential in the two regimes of k <<1 and k>>1, whereas D/N refers to diffusion of N particles in the abscence of a periodic potential. (a) compares B1 to D/N, we see for this value of A that B1 falls very sharply, for a lower value of A this would move closer to D/N. (b) compares B2 to D/N. There is a small difference between the two, If you were to increase A there would be greater differences.(c) compares B1 and B2. Here A is set at 7.5×10−21, with the diffusion constant set by D = a2/2δ with a = 0.02nm and delta = 2. These plots
al-low us to see both the difference between diffusion in the abscence and prescence of a periodic potential, these subtle differences could be confusing for analysis.
(a) (b)
(c)
Figure A3: Diffusion constant as a function of the amplitude of the periodic potential A. In these figures D/N = 1/N× a2
2δ B1 is eq.4.1 and B2 is D0 =
D/N×I−2
o (A/(2kBT)). Here, A is varied from 1×10−22J to 1×10−20J, which
is not a large variation but allows us to see how quick the values of the diffusion constant diverge. The number of particles is set at N = 10 with the same other parameters as set in fig.A2.
Bibliography
[1] K. Luger, R. Richmond, and A. Mader. “Crystal structure of the nucleosome core particle at 2.8A resolution”. In: Nature 389 (1997), pp. 251–260.DOI: 10.1038/38444.
[2] E. Segal and J Widom. “A genomic code for nucleosome position-ing”. In: Nature 442 (2006), pp. 772–778.DOI: https://doi.org/10.
1038/nature04979.
[3] B. Eslami-Mossallam et al. “Multiplexing Genetic and Nucleosome Positioning Codes: A Computational Approach”. In: PLos ONE 11 (6) (2016).DOI: https://doi.org/10.1371/journal.pone.0156905.
[4] W. K. Olson, A Gorin, and X Lu. “DNA sequence-dependent de-formability deduced from proteinˆaDNA crystal complexes”. In: Pro-ceedings of the National Academy of Sciences 95, 11163 (1998).
[5] F. Lankas, J. Sponer, and T Cheatham. “DNA basepair step deforma-bility inferred from molecular dynamics simulations”. In: Biophysical journal 85, 2872 (2003).
[6] M. Zuiddam, R. Everaers, and H. Schiessel. “Physics behind the me-chanical nucleosome positioning code”. In: Physical Review E 96(5) (2017).DOI: https://doi.org/10.1103/PhysRevE.96.052412. [7] F. Mueller-Planitz, H. Klinker, and P. Becker. “Nucleosome sliding
mechanisms: new twists in a looped history”. In: Nature Structural Molecular Biology 20 (9) (2013), pp. 1026–1032. DOI: 10.1038/nsmb.
2648.
[8] S. Yomosa. “Soliton excitations in deoxyribonucleic acid (DNA) dou-ble helices”. In: Physical Review A, 27 (4) (1983), pp. 2120–2125.DOI: 10.1103/physreva.27.2120.
[9] I.F Lyuksyutov. “Two-Dimensional Crystals”. In: English Translation: Academic Press, Boston (1988).
[10] Braun and Y. Kivshar. “Two-Dimensional Crystals”. In: Nonlinear dy-namics of the Frenkel-Kontorova model 306 (1998), pp. 1–108. DOI: 10. 1016/s0370-1573(98)00029-5.
[11] E. Sanz and D. Marenduzzo. “Dynamic Monte Carlo versus Brow-nian dynamics: A comparison for self-diffusion and crystallization in colloidal fluids”. In: The Journal of Chemical Physics 132 (19) (2010), p. 194102.DOI: 10.1063/1.3414827.
[12] H. Schiessel and I. Kulic. “Chromatin Dynamics: Nucleosomes go Mobile through Twist Defects”. In: Physical Review Letters 91 (14) (2003).DOI: 10.1103/physrevlett.91.148103.
[13] W. Dieterich and I. Peschel. “Diffusion in Periodic Potentials”. In: Z. Physik B 27 (1977), pp. 177–187.
