Chaotic Dynamics of Polyatomic Systems with an Emphasis on DNA Models
CChaotic Dynamics of Polyatomic Systems with anEmphasis on DNA Models
Malcolm Hillebrand
A thesis presented for the degree of Doctor of Philosophy in theDepartment of Mathematics and Applied Mathematics
University of Cape TownDecember 2020 eclaration
I, Malcolm Hillebrand, hereby declare that the work on which this thesis is based is myoriginal work (except where acknowledgements indicate otherwise) and that neither thewhole work nor any part of it has been, is being, or is to be submitted for another degreein this or any other university. I authorise the University to reproduce for the purposeof research either the whole or any portion of the contents in any manner whatsoever.Signature: Date: 23 rd of December 2020.2 bstract In this work we investigate the chaotic behaviour of multiparticle systems, and in par-ticular DNA and graphene models, by applying various numerical methods of nonlin-ear dynamics. Through the use of symplectic integration techniques—ef fi cient routinesfor the numerical integration of Hamiltonian systems—we present an extensive analy-sis of the chaotic behaviour of the Peyrard-Bishop-Dauxois (PBD) model of DNA. Thechaoticity of the system is quanti fi ed by computing the maximum Lyapunov exponent(mLE) across a spectrum of temperatures, and the effect of base pair disorder on the dy-namics is investigated. In addition to the inherent heterogeneity due to the proportionof adenine-thymine (AT) and guanine-cytosine (GC) base pairs, the distribution of thesebase pairs in the sequence is analysed through the introduction of the alternation index.An exact probability distribution for arrangements of base pairs and their alternationindex is derived through the use of Pólya counting theory. We fi nd that the value of themLE depends on both the base pair composition of the DNA strand and the arrange-ment of base pairs, with a changing behaviour depending on the temperature. Regionsof strong chaoticity are probed using the deviation vector distribution, and links be-tween strongly nonlinear behaviour and the formation of bubbles (thermally inducedopenings) in the DNA strand are studied. Investigations are performed for a wide vari-ety of randomly generated sequences as well as biological promoters. Furthermore, theproperties of these thermally induced bubbles are studied through large-scale moleculardynamics simulations. The distributions of bubble lifetimes and lengths in DNA areobtained and discussed in detail, fi tted with simple analytical expressions, and a physi-cally justi fi ed threshold distance for considering a base pair to be open is proposed andsuccessfully implemented. In addition to DNA, we present an analysis of the dynami-cal stability of a planar model of graphene, studying the behaviour of the mLE in bulkgraphene sheets as well as in fi nite width graphene nanoribbons (GNRs). The well-attested stability of the material manifests in a very small mLE, with chaos being a slowprocess in graphene. For both possible kinds of GNR, armchair and zigzag edges, themLE decreases with increasing width, asymptotically reaching the bulk behaviour. Thisdependence of the mLE on both energy density and ribbon width is fi tted accuratelywith empirical expressions. 4 ublications and Presentations ofResults From this PhD research, results have been presented through publications in interna-tional peer reviewed journals:1.
M. Hillebrand , G. Paterson-Jones, G. Kalosakas and Ch. Skokos,
Distributionof Base Pair Alternations in a Periodic DNA Chain: Application of Pólya Count-ing to a Physical System , Regular and Chaotic Dynamics, No. 2, 135 (2018),doi.org / / S1560354718020016.2.
M. Hillebrand , A. Schwellnus, G. Kalosakas and Ch. Skokos,
Heterogeneity andChaos in the Peyrard-Bishop-Dauxois Model of DNA , Physical Review E, , 022213(2019), doi.org / / PhysRevE.99.022213.3.
M. Hillebrand , B. Many Manda, G. Kalosakas, E. Gerlach and Ch. Skokos,
ChaoticDynamics of Graphene and Graphene Nanoribbons , Chaos , 063150 (2020),doi.org / / M. Hillebrand , G. Kalosakas, Ch. Skokos and A.R. Bishop,
Distributions of Bub-ble Lifetimes and Bubble Lengths in DNA , Physical Review E, , 062114 (2020),doi.org / / PhysRevE.102.062114.Parts of the results of this thesis have been presented, at conferences and seminars, asfollows:1. Seminar Presentation:
Chaotic Dynamics in DNA Chains , Centre for Theoreticaland Mathematical Physics, University of Cape Town (Cape Town, South Africa,March 2018).2. The International Conference on Nonlinear Localization in Lattices – NLL 2018,Poster Presentation:
Chaotic behaviour of the Peyrard-Bishop-Dauxois model of DNA (Spetses, Greece, June 2018).3. Seminar Presentation:
Chaotic Dynamics of DNA and Graphene , Professur für As-tronomie und Lohrmann-Observatorium, Technische Universität Dresden (Dres-den, Germany, February 2019).4. Postgraduate Science Symposium, Stellenbosch University:
Applying Dynamicsto Biology: Understanding the Chaos of our DNA (Stellenbosch, South Africa,September 2019).5. The 62nd Annual Congress of the South African Mathematical Society, Oral Pre-sentation:
Chaotic Dynamics of DNA (Cape Town, South Africa, December 2019).5 ontents fi cance of Open Base Pairs and Denaturation . . . . . . . . . 262.2 The Peyrard-Bishop-Dauxois Model of DNA and its Precursors . . . . . 282.2.1 The Poland-Scheraga Model . . . . . . . . . . . . . . . . . . . . . . . . 292.2.2 The First Peyrard-Bishop Model . . . . . . . . . . . . . . . . . . . . . 292.2.3 Peyrard-Bishop-Dauxois: Nonlinear Coupling . . . . . . . . . . . . 312.2.4 The Extended Peyrard-Bishop-Dauxois Model . . . . . . . . . . . . 35 / GC composition . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.2.1 Computing an average mLE . . . . . . . . . . . . . . . . . . . . . . . . 696haotic Dynamics of Polyatomic Models4.2.2 Average mLE Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.3 Heterogeneity: The Effect of Alternations . . . . . . . . . . . . . . . . . . . 744.3.1 Biological Sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.4 Comparison with the Extended PBD Model . . . . . . . . . . . . . . . . . . 804.5 Deviation Vector Distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.5.1 DVD and Displacement at Low Energies . . . . . . . . . . . . . . . 824.5.2 DVD and Displacement at Higher Energies . . . . . . . . . . . . . 854.5.3 DVD of Promoter Sequences . . . . . . . . . . . . . . . . . . . . . . . 864.5.4 Effect of Base Pair Sequence on the DVD . . . . . . . . . . . . . . . 894.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
CONTENTS Malcolm Hillebrand 7 hapter 1Introduction and Numerical Methods
This thesis is a foray into the wonderful world of nonlinear dynamics, being appliedto perhaps surprising places. We will explore some of the possibilities of implementingnumerical methods of nonlinear dynamics, especially chaotic dynamics, to polyatomicsystems. The two systems explored are both fascinating in their own right: DNA isentrenched in the public eye as a cornerstone of modern biology (as a cursory glancethrough any textbook will reveal, e.g. [ ] ), and graphene has made waves after its intro-duction as a supermaterial for the future [ ] . There have been many ways of studyingthese systems, as will be discussed in detail in the relevant chapters, but it is not par-ticularly common to approach this problem from a mathematical physics standpoint,let alone from a chaotic dynamics point of view. This gives us the opportunity here topush forward the understanding of the physical systems through applying these tech-niques, but also to develop the understanding of how we can use chaotic dynamics ina variety of situations. With regard to the fi rst aim, we hope to bring out fresh in-formation about the systems that can contribute to further studies. The signi fi canceof DNA and graphene (as well as the physical robustness of the models used to studythem) provides ample motivation for us to use any tools we have at our disposal to tryand better understand the functioning and potential uses of these substances. The par-ticular well-documented biophysical importance of openings in the DNA double helix(a fundamentally dynamic process) in the context of transcriptional activity adds furtherphysical motives to the study. In general, it is advantageous to implement different ap-proaches to studying a problem, as they are likely to provide some new understandingand insight. In an increasingly interdisciplinary age, the second aim is an important one– exploring the nonlinear dynamics side of physical systems, and identifying what thepossibilities are for the application of chaos theory to these real substances. This ap-proach could give answers to some important questions, such as: How can we use chaosindicators to understand the phenomenal structural strength of graphene? Is there a wayto probe the transcriptional properties of DNA through nonlinear dynamics signatures?Furthering the state of knowledge in this direction is a promising avenue to pursue.So here we will give outlines of the structure of DNA and graphene, as well as theirbiochemical and physical signi fi cance, and also of the various numerical and theoreticalmethods we have used in studying these polyatomic systems. These methods includechaos detection techniques, mathematical aspects (as Chapter 3 will reveal, sometimes very mathematical), high performance computing ideas – including multi-thread par-allelism on central processing units (CPUs) and graphical processing units (GPUs), as8haotic Dynamics of Polyatomic Modelswell as adjustments of these numerical techniques and tricks speci fi c to our applicationof chaotic dynamics.In what follows in this chapter, we will brie fl y introduce the notions of Hamiltoniandynamics and chaos, as well as the main players in our study: The maximum Lyaponuvexponent (mLE) as the primary chaos indicator we have used, the deviation vector dis-tribution (DVD) as a means of identifying regions of strong chaoticity, and the essentialnumerical techniques of symplectic integration and parallel computing for integratingsystems of differential equations.The remainder of this thesis is structured as follows. In Chapter 2 we present somebackground to the dynamical study of DNA, with a particular emphasis on the modelsused in our study. Chapter 3 describes the mathematical basis for our quanti fi cation ofheterogeneity in DNA sequences, using an extension of Pólya’s enumeration theorem.Results from the investigation of the chaotic dynamics of DNA are given in Chapter 4,with emphasis on the effect of heterogeneity on the system’s chaoticity. Taking a differ-ent approach to studying DNA, Chapter 5 presents results from an investigation intothe physical properties of DNA, with a particular focus on thermally induced open-ings in the double strand called bubbles. Chapter 6 then presents a dynamical model ofgraphene, some computational details for the implementation of parallelisation, and astudy of the chaotic dynamics of 2D graphene models. Finally, in Chapter 7 we concludethe study, summarise our fi ndings and discuss future outlooks. Within the broad scope of dynamical systems, we are particularly focussing here onHamiltonian systems. In any Hamiltonian system, the dynamics are governed by aHamiltonian function H . This function generally serves as the energy of the systemas well as a generator of the time evolution [ ] . If we have generalised positions q , andconjugate momenta p which de fi ne the state of our dynamical system, then we can writethe Hamiltonian expressed in terms of these variables, and possibly also time as H ( q , p , t ) .In a closed or isolated system, the Hamiltonian is simply the sum of the system’s kinetic ( T ) and potential ( V ) energies, giving H ( q , p , t ) = T + V , (1.1)where it is quite straightforward to understand the Hamiltonian as the total energy ofthe system. If in addition to this the Hamiltonian is not explicitly time dependent, sothat ∂ H ∂ t =
0, (1.2)then we call the dynamical system autonomous. The Hamiltonian will usually of coursestill depend on time implicitly, through the time-dependent coordinates and momenta.So we would have an autonomous Hamiltonian in the form H ( q ( t ) , p ( t )) ,where the explicit t -dependence has fallen away.Chapter 1 Malcolm Hillebrand 9haotic Dynamics of Polyatomic ModelsSigni fi cantly, time independent Hamiltonians conserve the total energy of the sys-tem (the value of the Hamiltonian function itself) [ ] . This is extremely useful not onlyfor simplifying analytical calculations, but also makes a lot of computational analysismuch easier. For instance, it immediately gives us a conserved quantity to check ournumerical integration schemes with – if our computation does not conserve the initialvalue of the Hamiltonian (up to numerical accuracy), then we know we have a prob-lem in the machinery somewhere. We will typically measure this energy conservationthrough the relative energy error, E Re l = | H ( t ) − H ( ) || H ( ) | , (1.3)which tells us how close we are staying to the initial energy value.Quite often in autonomous systems, the Hamiltonian is “separable”, i.e. (1.1) can beexpressed as H ( q , p ) = T ( p ) + V ( q ) , (1.4)where T ( p ) is the kinetic energy depending only on the momentum of the system, and V ( q ) is the potential energy, a function of only the generalised coordinates.A further useful aspect of Hamiltonian mechanics is the existence of the Poissonbracket. The Poisson bracket of two functions f and g of the generalised positions andmomenta in an N -dimensional system with coordinates ( q , p ) is de fi ned as [ ] { f , g } = N (cid:1) i = (cid:2) ∂ f ∂ q i ∂ g ∂ p i − ∂ f ∂ p i ∂ g ∂ q i (cid:3) . (1.5)As our main use case for the Poisson bracket, the time derivative of a function with noexplicit time dependence in a Hamiltonian system is given by d fd t = { f , H } , (1.6)which provides a route to Hamilton’s equations of motion, among other things.In general, much of the signi fi cance of Hamiltonian systems consists in the analyticaland numerical methods that it is possible to construct for them. Especially for conser-vative Hamiltonian systems, the beautiful properties of the conserved volume of thephase space – the space spanned by the set of coordinates q and p – lend themselves tosome very neat computations, which in our case enable the ef fi cient study of the chaoticdynamics of these systems through symplectic integration (see Section 1.2.1), and thevariational equations (Section 1.3.2).Relevant details of chaotic behaviour will be given in the discussion of the maximumLyapunov exponent (Section 1.3), but as a brief introduction we will present the basicframework for chaotic motion. An excellent introduction to chaotic dynamics is givenin [ ] , for a range of dynamical systems. For now, we are concerned with chaotic dy-namics of Hamiltonian systems, but of course the basic characteristics of fully chaoticsystems hold here as well [ ] :1. Sensitivity to initial conditions.2. Topological mixing.10 Malcolm Hillebrand Chapter 1haotic Dynamics of Polyatomic Models3. Density of periodic orbits.The sensitivity to initial conditions is by and large well understood. For any giventrajectory in a chaotic region of the phase space an initially nearby orbit will eventuallydiverge signi fi cantly from it. This will be our crucial point in identifying and quantify-ing chaos later. The topological mixing means that a chaotic orbit will eventually visitevery region in the phase space, to an arbitrarily small closeness. The result of this is thatthe chaotic phase space cannot be separated into different subsystems, but rather has tobe treated as a single whole. Finally, the density of unstable periodic orbits means that inthe phase space each point will have an unstable periodic orbit in a small neighbourhoodaround it, yielding some sense of regularity in the system.Consequently, much of the study of chaotic or mixed (partly chaotic and partly reg-ular phase space) systems involves trying to identify chaotic regions or trajectories basedon these criteria. For low-dimensional systems it is sometimes possible to directly in-tegrate a collection of initial conditions and identify chaotic regions through observingthe dynamical behaviour, such as in the archetypical Lorenz system for atmospheric fl ow [ ] . The phase space of these low dimensional systems can also be studied throughtaking surfaces of section, such as the Poincaré surface of section (PSS) [ ] , which canreveal chaotic and regular regions. This has for example successfully been applied toastronomical models, such as the well known Hénon-Heiles system [ ] .Higher dimensional systems (for instance of hundreds of degrees of freedom) makethis direct approach challenging to say the least. While methods exist for visualising fourdimensional phase spaces such as can arise in the study of three dimensional Hamilto-nian models [
9, 10 ] , visualising fi ve hundred dimensions is a task beyond most humanminds. Fortunately, there are direct methods available for identifying chaotic motion,and distinguishing it from regular behaviour. The method of chaos identi fi cation andquanti fi cation used in our work is the mLE, described in detail in Section 1.3. Furtheroptions include the smaller (SALI) and generalised (GALI) alignment index methods [ ] , which are very ef fi cient ways of identifying chaotic trajectories. Recent reviewsof various modern chaos detection techniques can be found in [ ] . In any research project based largely on computational simulations, it is paramount touse effective and ef fi cient numerical methods. In light of this, there was a signi fi canteffort in the course of this study to continuously improve the state of the computer codeused for integration and analysis, and throughout the process new ideas were consideredfor optimisation. The two major tools used in the optimisation process were symplecticintegrators, and code parallelisation. Here we brie fl y discuss these two aspects, with afocus on the practical side of how they were actually implemented in the project. For the numerical integration of the equations of motion of Hamiltonian systems, aspeci fi c class of methods known as symplectic integrators have been introduced [ ] withgreat success [ ] . These are integration routines that were derived bearing in mind thespeci fi c properties of Hamiltonian phase space, and particularly of conservative Hamil-tonian systems, in order to provide accurate long-time simulations. A great resource forChapter 1 Malcolm Hillebrand 11haotic Dynamics of Polyatomic Modelsa thorough overview of the topic is [ ] , but here we will present a practical outline ofwhat these methods are and how they are implemented in our case.The basic idea of a symplectic integration method is instead of approximately inte-grating the exact Hamiltonian (as would happen with a traditional Runge-Kutta inte-grator for instance), to rather exactly integrate an approximate Hamiltonian, by usinga sequence of exact symplectic mappings [ ] . The bene fi t of this approach is that byconserving the symplectic structure of the phase space, we in fact conserve the energyprecisely. This is a result of the symplectic preservation, which corresponds to conserv-ing the volume of the phase space and consequently the energy of the system. The magicof this particular approach lies in the fact that it is possible to design these integrationmethods based on a desired Hamiltonian, and end up with a routine that integrates avery similar Hamiltonian with complete accuracy, allowing us to integrate this slightlyperturbed Hamiltonian for arbitrarily long times, and with arbitrary closeness to theoriginal Hamiltonian set by the time step.As outlined very well in [ ] , a basic approach to constructing these algorithms is toapproximate the exact canonical transformation ( q , p ) → ( q ( t ) , p ( t )) , which is usuallyimpossible to write down, with a series of simpler canonical transformations. So let ustake a separable Hamiltonian, of the form (1.4) H = T ( p ) + V ( q ) .With this Hamiltonian there is an associated time evolution operator of the system’sstate exp ( ˜ H t ) , where ˜ H represents an operator that takes the Poisson bracket of itsargument with the Hamiltonian: ˜ H ( X ) = { X , H } . Thanks to the separability of theHamiltonian, this time evolution operator can be written in terms of two independentoperators A ( X ) = { X , T } and B ( X ) = { X , V } corresponding to the kinetic and potentialterms respectively as [ ] exp ( ˜ H t ) = exp [( A + B ) t ] . (1.7)The trick now is to try to approximate this exponential operator as accurately aspossible, by only using simple operations of exp ( At ) and exp ( B t ) which can be obtainedanalytically in the case of Hamiltonians in the form of (1.4) [ ] . Of course for mostcases the operators do not commute, and the exponential does not simply split out intoexp [( A + B ) t ] = exp ( At ) exp ( B t ) , but rather a sequence of exponentials with carefullychosen coef fi cients c i , d i is necessary to approximate the true operator. This sequence ofexponentials is necessarily fi nite, and for an approximation with n applications of eachoperator we have exp [( A + B ) t ] ≈ n (cid:4) i = exp ( c i At ) exp ( d i B t ) . (1.8)How to choose these coef fi cients, and in general the art of constructing symplectic in-tegrators is a subject that has attracted much attention (see e.g. [ ] ).The accuracy of a given symplectic integration scheme is de fi ned by the order of theerror in the approximation of the true time evolution operator. In the approximation of(1.8), the error term will be of some fi xed order in the time step. So with an integrationtime step of τ , we would have more preciselyexp [( A + B ) τ ] = n (cid:4) i = exp ( τ c i A ) exp ( τ d i B ) + (cid:4) ( τ m + ) , (1.9)12 Malcolm Hillebrand Chapter 1haotic Dynamics of Polyatomic Modelswhere the integration scheme is said to be of order m ∈ (cid:1) . Exactly what m is for a givenmethod is dependant on the number of terms in the product and the particular choice ofcoef fi cients c i and d i . The aim with designing these methods is of course to reduce theorder of the error as far as possible, while having to perform as few operations at eachtime step as possible.We will now discuss some practical implementation details of these methods. Interms of the computations, a single application of an A or B operator is the same asupdating the positions q or momenta p respectively. As described in [ ] , this meansthat for a single time step we have a sequence of updates to the positions and momentausing the usual Hamiltonian equations of motion: p i + = p i − c i τ ∂ V ( q i ) ∂ q , (1.10) q i + = q i + d i τ ∂ T ( p i + ) ∂ p . (1.11)The index i refers to the substeps within one integration time step of the process; so thesystem has only been evolved by τ time units after all these substeps have been com-puted.The simplest form of symplectic integration is to simply take a fi rst order approxima-tion and set c = d =
1, and use only a single substep. This is commonly referred to as thesymplectic Euler method for numerical integration [ ] . For more involved multistepintegrators, these substeps are applied consecutively with the appropriate coef fi cients,and produce accurate evolutions of the system with a varying degree of precision de-pending on the method and the system under study.As a speci fi c example of symplectic integration in action, we can consider the fourth-order, six-stage Symplectic Runge-Kutta-Nyström integration scheme SRKN b [ ] . Thisis the integrator that is used in obtaining the results of Chapters 4 and 5, to integrate theequations for the DNA model. The integrator takes the formexp [( A + B ) t ] ≈ exp ( c B t ) exp ( d At ) . . . exp ( c B t ) exp ( d At ) exp ( c B t ) . (1.12)Using the iterative procedure of Eqs. (1.10) and (1.11), this routine requires six updatestages for the A operator, and seven stages for the B operator. The coef fi cients for thisintegration scheme are given in Table 3 of [ ] , showing the optimised values found to tominimise errors. With these values, we are able to apply the iterative procedure directly,and implement this integrator as a computational algorithm.As an illustration of the bene fi t of symplectic integration, let us consider the simula-tion of the DNA model studied in this thesis. The details of the model and its equationsare discussed in Chapter 2, but here we will simply look at one aspect of the system – theconservation of energy, as this is the primary advantage of SIs in our application. Fig-ure 1.1 shows the relative energy error (1.3) throughout the integration process for thesymplectic scheme for a particular initial condition, compared to that of the commonly-used fourth order Runge-Kutta integrator (RK4). Here the issue with the non-symplecticscheme quickly becomes apparent – with the same integration time step, the RK4 inte-grator and SRKN b have a very similar total CPU time (as seen in the inset), but evenfrom the beginning the accuracy of the RK4 method is much worse than the symplecticscheme. Even using a time step ten times smaller, resulting in very slow integration, theChapter 1 Malcolm Hillebrand 13haotic Dynamics of Polyatomic Models GNU Parallel
The tool we have used for the simple running in parallel of multiple independent pro-grams is GNU Parallel [ ] . The process here is quite straightforward. Say that we havea compute node with 24 CPUs, as was our case on the CHPC cluster, and plenty ofmemory for our application. We note that in most of our numerical integrations thecomputer code is strongly compute bound, i.e. never coming close to running out ofmemory even when parallelised. If we also have a large number of individual programsto run, each using only one CPU core, then it would be convenient to be able to runthese simultaneously (24 at a time) rather than having to run them inef fi ciently in series.This is precisely what GNU parallel enables us to do – arrange our multiple individual fi les to be fed through the processors in batches.In the simplest case, we will take serial programs, where each program uses only asingle core. This situation can arise in several ways; the main way which it was used inour application was running many cases of the DNA model (see Section 4.1.1), wherethe same code needs to be run with a large number of different parameters and initialcondition. These cases can be split into separate executable fi les, each running the codewith the required conditions, completely independently of the other programs, and runusing a single CPU core. Whenever this splitting is possible, it is inherently the mostef fi cient form of parallelisation possible, guaranteeing 100% ef fi ciency as the numberof cores is increased (assuming no memory-related slowdown, since there is no needfor communication between the cores). The other situation where this parallelisationoccurs in our work is the analysis of large numbers of data fi les, which can be processedin parallel independently. Then again it is highly ef fi cient to analyse the data fi les indifferent blocks, each block being processed by a separate program.This is the fi rst choice for all the parallelisation we have used. If individual runs donot take extremely long, as is the case in all DNA results (see Chapters 4 and 5), then wehave used only this form of independent parallelism.GNU Parallel can also be used when each individual program is not purely serial.So for instance if each individual run requires 4 cores, and we still have access to the 24core node, then we can make use of GNU Parallel to run 6 instances of the code on eachnode. This allows for a balancing between the speedup of an individual program dueto internal parallelisation, and the ef fi ciency that comes from being able to run manyprograms at once. OpenMP and CUDA
Turning our attention to this internal parallelisation, we will roughly outline the meth-ods used in our computations. Primarily, our codes have been parallelised using CPUmultithreading, and particularly the OpenMP framework (see e.g. [ ] ). This multi-threading has been used in the relatively complex two dimensional graphene model dis-cussed in Chapter 6.Brie fl y, in applying multithreading instead of having one CPU core perform all thecomputations one after the other, we will divide the work between multiple cores andhave them perform the operations simultaneously. In the context of lattice integra-tion, this means when the integration step occurs we split the integration of the sitesbetween multiple cores. So if we assign 4 cores to the integration, that means we willgive each core one quarter of the atoms to integrate, and these cores will execute simulta-Chapter 1 Malcolm Hillebrand 15haotic Dynamics of Polyatomic Modelsneously. During each application of the integration operators exp ( τ c i A ) and exp ( τ d i B ) ,the equations of motion for each atom are computed and the position / momentum up-dated. Since these equations of motion depend strictly on the positions and momenta atthe previous time step, there is no problem with computing them simultaneously. Eachatom can be updated seperately, independently of the changes to the other atoms. Atthe end of each step, the code ensures that all threads (or cores) have fi nished all theirassigned sites, in case some operations take slightly longer or shorter. This then avoidscomplications of trying to continue the integration with the next step while some siteshave not in fact been updated yet.When parallelising code it is always necessary to take precautions to ensure that thereare no race conditions, such as problems arising from different threads attempting to ac-cess data before it has been updated or trying to modify memory that is being used byanother thread. This demands frequent synchronisation checkpoints to be sure that allthreads have completed their assigned tasks. Because of this, as well as other inef fi cienciesin the parallelisation process (e.g. cache usage, necessarily serial parts, slower individualprocessors), there are limits to how much speedup can be achieved by this kind of mul-tithreading [ ] . In Section 6.2, we will present some results discussing this imperfectspeedup, and how we identi fi ed a reasonably optimal number of threads to use for eachindividual integration.A second possibility for multithreading is using GPU parallelisation. For this, wewrote the integration code in CUDA (for an excellent introduction, see [ ] ), as themost ef fi cient available way of making use of the NVIDIA GPU architecture, the mostcommon today. The primary distinction between GPUs and CPUs is that GPUs tendto have many more processing cores (for instance, the Tesla V100 GPUs deployed in theCHPC cluster have an impressive 5120 CUDA cores, compared to the 24 CPU proces-sors available on a single node), but with lower processing speeds (the V100 has a peakof 1455MHz, while the CHPC Intel Xeon processors have a peak of 2600MHz). Thisresults in GPUs being much more powerful for computations involving signi fi cant par-allelisation, and the integration of very large numbers of atoms. The design of CUDAprogramming also lends itself to very ef fi cient ways of handling enormous lattices (tensor hundreds of thousands of atoms), through various techniques [ ] .In CUDA programming, the same ideas as for CPU parallelisation are followed: Weassign to each core a number of sites, which it then goes through and integrates. Dueto the complexities of interfacing between the GPU and normal CPU processing, alongwith necessarily very careful memory management, integrating the lattice using CUDAis more dif fi cult than using CPU parallelisation, and only provides meaningful advan-tages with very large lattices. Up until every GPU core is in use (for our application,that corresponds to at least 5120 atoms in the lattice since we are using a GPU with 5120cores), all computations take the same amount of time – with the lower clock speed, thismeans that we are unlikely to be doing better than in the CPU parallelisation regime.However, once this point is reached, there will be the expected linear relationship be-tween compute time and lattice size.As with the CPU performance, scaling results for an implementation of the dynam-ical model of graphene in CUDA are presented in Section 6.2.16 Malcolm Hillebrand Chapter 1haotic Dynamics of Polyatomic Models General Optimisation and Implementation Details
Apart from the parallelisation of the integration code, here some of the other compu-tational details are presented. As previously mentioned, numerical integration of equa-tions of motion is an intensive task largely limited by the number of arithmetic opera-tions required at each iteration. This means that a vital step in the optimisation of ourcode is to fi rst ensure that the equations are written in a computationally ef fi cient form –shortening equations as much as possible, checking that all cancellations are performed,and minimising the number of complex operations (e.g. trigonometric or exponential)that are performed at each step.Beyond this, several optimisation steps were taken in the course of the developmentof the integration code for the DNA and graphene models. They are listed here, as theyare in general very effective steps to take when optimising code of this sort.1. Focussing on the integration steps, which are performed millions of times in thecourse of a run. Making sure that only absolutely necessary computations wererepeated each time; anything that can be taken out of this step should be.2. Choosing an ef fi cient integrator for each model. This involved running compar-ative tests to identify the best option given the speci fi cs of the model, the size ofthe time step required, the accuracy required, and so on.3. Identifying quantities that were calculated repeatedly, and pulling them out at thebeginning of each step as pre-calculated constants.4. Once a basic version of the code is working, rework the equations to fi nd anyways of saving time; perhaps certain elements of the integration can be computedtogether ef fi ciently (for instance sines and cosines in several programming lan-guages).5. Code pro fi ling to identify bottlenecks, and inform the decisions of where to focusattention. Tools such as gprof and NVprof are useful in this regard.6. Becoming familiar with the compiler-speci fi c optimisations available, in terms ofvectorisation of operations, interprocedural optimisation, and optimisations forcertain processors.7. For parallelised code, running scaling tests to have a clear idea of an effective num-ber of cores to use for a given system size. This is always a trade-off between theoptimal ef fi ciency, which is running in serial, and how long the code takes to run.For these projects, a combination of codes written in Fortran90, C ++ , Python,Mathematica and Octave were used. The simulations were performed using Fortran90(DNA models) and C ++ (graphene model), compiled using IFort and ICC / ICPC re-spectively. Data analysis used primarily Python, with Fortran90, Octave and Mathe-matica where useful. In particular, the Python libraries numpy , scipy and matplotlib were used extensively.Linux shell utilities such as bash, AWK, sed, and grep were invaluable in manipulat-ing fi les, arranging data, and generally saving large amounts of time on repetitive tasks.Chapter 1 Malcolm Hillebrand 17haotic Dynamics of Polyatomic Modelsby Oseledec [ ] , where his Multiplicative Ergodic Theorem opened the way for thecomputation of the LCEs, and in particular the largest of the LCEs, the maximal Lya-punov characteristic exponent (mLCE), or just the maximal Lyapunov exponent (mLE).This led to the development of techniques to estimate the mLE, based on the relation-ship between the LCEs and the exponentially fast divergence of perturbed chaotic orbits.The fi rst estimation of the mLE, χ , was presented as the asymptotic value when t → ∞ of a fi nite time mLE (ftmLE), X ( t ) [ ] , which can be computed directly based on themean exponential rate of growth of the deviation, the details of which were outlined inthe seminal work of Benettin et al [
32, 33 ] . The “standard method” developed in thiswork has remained the most commonly used technique for computing mLEs, provingto be a very effective tool in a diverse array of dynamical contexts [ ] .Before further discussion of the mLE, let us actually de fi ne it. Given some referenceorbit, e.g. x ( t ) in Fig. 1.2, and a deviation from this trajectory called w ( t ) de fi ning anearby orbit y ( t ) (Fig. 1.2), which evolves through time, then the mLE is de fi ned to bethe quantity χ = lim t →∞ t ln (cid:2) || w ( t ) |||| w ( ) || (cid:3) , (1.13)where || · || denotes the Euclidean norm (although any choice of norm produces the cor-rect result [ ] ). Based on this formulation, we can draw some conclusions, many ofwhich were discussed originally in [ ] . Firstly, the fact that there is a prefactor of 1 / t when we are taking the in fi nite time limit means that for any non-growing or slowly-growing deviation vector, the mLE value is going to be zero. This corresponds to thestatement that for regular orbits the mLE is zero, since the implication of a slowly grow-ing deviation is that the orbit is not chaotic. For a chaotic orbit we have an exponentiallygrowing deviation vector [ ] , and then despite taking the logarithm of || w ( t ) || this willcontinue to grow in time, and will counteract the 1 / t to produce a constant value. So wehave the second statement that chaotic orbits produce positive mLE values. Immediatelywe now have our desired chaos identi fi er; provided that we can somehow determine themLE of a given orbit, we know whether or not it is chaotic. This alone is extremelyuseful, especially in “mixed” systems, where regions of the phase space exhibit chaoticmotion and other regions are regular. By calculating mLE values in different regions, wecan then identify islands of stability in a chaotic sea, or regions of chaos in a generallyregular phase space.The second powerful aspect of the mLE is that it not only identi fi es chaos, but itquanti fi es the chaoticity. We can see that if a deviation of orbit x A , w A ( t ) , grows fasterthan the deviation from an orbit x B , w B ( t ) , that the resultant mLE of orbit x A is go-ing to be larger, simply by the de fi nition. This means that by calculating the mLE fordifferent orbits in a system, we can identify stable and chaotic regions, but we can alsodistinguish the more and less chaotic behaviours. The estimation of the mLE also facili-tates the comparison of the chaoticity of different states of the same system (for physicalsystems such as studied in this work, these states could correspond to the different en-ergy or temperature values), and even the comparison between the chaoticity of differentsystems.As a consequence of these properties, the mLE is a tremendously useful tool forstudying dynamical systems, when it can be calculated. The obvious catch with thisstatement is being able to calculate the mLE. In many cases it is completely impossibleto theoretically calculate the rate of growth of the deviation vector, which in turn rendersChapter 1 Malcolm Hillebrand 19haotic Dynamics of Polyatomic Modelsthe mLE incalculable. In other cases even if some estimates can be made for the growthof the deviation vector, the exact calculation of the in fi nite limit creates problems.As introduced in [ ] and discussed in clear detail in [ ] , there have been very ef-fective algorithms developed to numerically estimate the value of the mLE. If we aregoing to use a numerical approach, the fi rst issue we have to clear away is the in fi nitelimit, as in fi nite compute times are not possible to achieve. To this end, the ftmLE, anestimation of the mLE with the in fi nite limit truncated to a certain time, is used. ThisftmLE is de fi ned according to X = t ln (cid:2) || w ( t ) |||| w ( ) || (cid:3) , (1.14)where from (1.13) we can see that χ = lim t →∞ X .The evolution of the ftmLE itself is useful, even apart from the asymptotic behaviourwhich is its most important aspect. For clearly chaotic orbits, the ftmLE can generallybe used fairly straightforwardly – if the mLE is a positive constant, then it follows thatat some point the ftmLE must converge to that constant value. So then by pursuing thenumerical integration of two nearby initial conditions, we can continue until reachinga point where it is clear that the ftmLE is practically no longer changing with time. Forregular orbits, or effectively regular orbits, it becomes a little less clear. While one couldin principle continue a computation until the ftmLE reaches some sort of “numericalzero”, based on the machine precision of the computer, this could take an extraordinaryamount of time. Fortunately though, it was shown early on [ ] that for regular orbits,the ftmLE should tend to zero with a gradient of 1 / t . This can be seen from the def-inition of the ftmLE, (1.14), since a small growth in the deviation vector size will getcompletely lost in the overwhelming strength of the prefactor of 1 / t , leading to the evo-lution of the ftmLE being dominated by the slope of 1 / t . Additionally, we note thatas long as the deviation vector grows linearly in time, the ftmLE will evolve accordingto a modi fi ed power law of ln ( t ) / t , which will become closer to 1 / t in the long timelimit. This fact is particularly useful in light of the behaviour of so-called “weak chaos”regimes [ ] , where the ftmLE appears to be continuously decreasing, but this is not are fl ection of a regular system. As a result of this phenomenon, it has proved useful bothin studying one dimensional lattices [ ] and two dimensional lattices [ ] to be able tocompare the slope of the ftmLE with the expected 1 / t gradient. These works furthercon fi rm the usefulness of studying the behaviour of the ftmLE, as well as the fi nal mLE.A practical concern with the computation of the ftmLE is the fact that we are tryingto investigate an exponentially growing variable (the norm of w ), over potentially verylong times while we wait for the ftmLE to converge. The combination of exponentialgrowth and long times is a recipe for numerical over fl ow disaster, and we need to fi nda way to circumvent this eventuality. As usual, the mathematical details can be foundin [ ] , but the fi nal outcome of practical importance is as follows. Instead of computingthe ftmLE X according to the norm of the deviation vector only at the fi nal time, itis possible to instead iteratively update the ftmLE at each time step of the numericalintegration, and renormalise the deviation vector (retaining its direction) at each step.With this, it is possible to write the ftmLE after n integration steps of the system witha time step of τ as [ ] X ( t = n τ ) = n τ n (cid:1) i = ln (cid:2) || w ( i τ ) |||| w ( ) || (cid:3) . (1.15)20 Malcolm Hillebrand Chapter 1haotic Dynamics of Polyatomic ModelsSo after each integration time step, the value of the deviation vector norm is used toupdate the ftmLE estimate, and then the deviation vector is renormalised, thus avoidingthe problem of exponentially growing values.As a fi nal comment on the use of the mLE in dynamical systems, let us note thatit is apparent from (1.13) that the units of the mLE are inverse time. By inverting the χ value, we fi nd the so-called Lyapunov time T L , which is an estimation of how longit takes for the system to become chaotic [ ] . This is a useful quantity in several con-texts, providing both an expected horizon for chaotic behaviour to appear in dynamicalsystems such as astronomical motions [ ] , and a characteristic time scale for chaos inphysical models.What has not actually been discussed yet is how exactly this deviation vector is com-puted through time. Here we present the two most common methods of evolving thedeviation vector: The two particle method (2PM) simply using two nearby initial condi-tions and evolving them separately, and the tangent map (TM) method using the so-calledvariational equations. The simplest, most direct approach to fi nding the evolution of the deviation vector is todirectly integrate two nearby initial conditions. Since this approach involves effectivelysimulating two nearby particles, this is known as the two particle method (2PM). Inthis process, we choose an initial condition we would like to investigate, and create thesystem with these coordinates. Then we initialise the deviation vector, typically by ran-domly perturbing each coordinate, and then normalise the vector to the required norm.Finally, we create a second instance of the lattice with initial conditions described by thedeviation from the initial lattice.More clearly, say we have a phase space point x ( ) = x , containing the initial po-sition and momentum coordinates of the orbit we wish to characterise as regular orchaotic. We then choose the deviation vector we wish to evolve, δ x ( t = ) . Using thisdeviation, we can de fi ne a second orbit y ( t ) with initial condition y ( ) = x + δ x ( ) .These two initial conditions can be integrated simply according to the Hamiltonianequations of motion of our system, and in general at some time t we can fi nd the cur-rent state of the deviation vector by fi nding the difference between the two trajectories, δ x ( t ) = y ( t ) − x ( t ) . In the evolution of these two trajectories, we are able to performthe renormalisation procedure discussed above, avoiding the problem of exponentialgrowth of the deviation vector.In an ideal world, this would be the end of it; we can choose our deviation howeverwe wish and then simply evolve the two initial conditions. However, there are someconsiderations that need to be made. Foremost among these is that we have to havea nearby initial condition for the second trajectory. This means we need to somehowchoose a suf fi ciently small deviation vector that our integration will correctly simulatetwo almost equivalent initial conditions. A fi rst thought would be why not make thisdeviation as small as possible? After all, we would really like to simulate an in fi nitesi-mally small deviation. What we are able to do is to minimise the deviation size based onthe available machine precision. Due to the renormalisation procedure, which needs tobe performed accurately, a lower limit on the norm of the deviation vector has been sug-gested to be about 10 − [ ] . Any smaller deviation will start to be affected by machineChapter 1 Malcolm Hillebrand 21haotic Dynamics of Polyatomic Modelsprecision errors since the square of the norm will be smaller than 10 − , which is be-yond the accuracy possible with double precision fl oating point numbers. The trade-offfor a small deviation vector is that the accuracy demanded of the integration is corre-spondingly higher – for a deviation vector norm of 10 − R , the relative energy error of theintegration should be on the order of 10 − R − or smaller to maintain an accurate simula-tion of the two nearby orbits [ ] . This accuracy is only possible through the decreasingof the integration time step (assuming that the integration scheme is optimal for this ap-plication), which means longer integration times. In some applications a larger deviationvector may produce accurate results, which allows for faster integration times. The exactchoice of deviation vector thus has to be decided on a case by case basis.When using this method, it is also advantageous to only update the ftmLE compu-tation and renormalise the deviation vector after a gap of several time steps. By renor-malising the extremely small deviation vector at every time step, there is a substantiallyincreased chance for fl oating point roundoff errors to accrue, yielding an inaccurate es-timation of the mLE. Managing fl oating point errors is particularly important sincechaotic systems are inherently very sensitive to small changes in the coordinates. Assuch a small accrual of inaccuracies in the deviation vector size will very likely becomeapparent in a badly behaved ftmLE. Much like the size of the deviation vector, an optimalrenormalisation time will depend on the system and the time step used. Approximatelyone time unit is recommended in [ ] , but with very long integration times or large timesteps, longer renormalisation times may be necessary. Once again, the aim in choosingrenormalisation times and time steps is to balance long-term accuracy with the desirefor precision.So implementing the 2PM, we have a direct way of estimating the mLE. While it isan old method [ ] , it is generally accurate and still relevant in many applications today(see e.g. [ ] ). A major advance in the road to fast and accurate computations of the mLE for dynamicalsystems, and Hamiltonian systems in particular, was brought about by the introductionof the so-called variational equations to integrate the evolution of the deviation vectorthrough time [ ] . The central idea here is rather than evolving two trajectories in thephase space and fi nding the difference between these two orbits, to directly simulate theevolution of a deviation vector in the tangent space. This is especially useful in the caseof Hamiltonian systems, where the phase space has many useful properties, and we areable to study the tangent space with relative ease.So if we wish to evolve a deviation vector w ( t ) from a reference orbit x ( t ) , then wecan do this using the variational equations. Following [ ] , we can de fi ne the Hamil-tonian equations of motion for our solution x ( t ) of an N -dimensional system (with acorrespondingly 2 N -dimensional phase space) as x ( t ) = g ( x ) = J N · ∂ H ∂ x (1.16)where ∂ H / ∂ x is the 2 N dimensional vector given by ∂ H ∂ x = (cid:2) ∂ H ∂ q , ∂ H ∂ q , . . . , ∂ H ∂ q N , ∂ H ∂ p , ∂ H ∂ p , . . . , ∂ H ∂ p N (cid:3) T , (1.17)22 Malcolm Hillebrand Chapter 1haotic Dynamics of Polyatomic Modelswith the q i and p i the position and momenta coordinates of the vector x . The matrix J N is the antisymmetric matrix J N = (cid:2) N I N − I N N (cid:3) , (1.18)where N is the N × N matrix of zeros and I N is the N -dimensional identity matrix. Wecan see that (1.16) is of course just the vectorised form of the well-known Hamiltonianequations of motion ˙ q i = ∂ H ∂ p i , ˙ p i = − ∂ H ∂ q i . (1.19)Now to study the tangent dyanamics of this orbit x ( t ) we need to take anotherderivative to fi nd the tangent vector, and use this to evolve the deviation w ( t ) forwardin time. To this end, we de fi ne the variational equations (the equations of motion forthe deviation vector) as ˙ w ( t ) = ∂ g ∂ x ( x ( t )) · w = J N · D H ( x ( t )) , (1.20)where D H ( x ( t )) denotes the Hessian matrix. We note that this is simply the lineari-sation of the evolution of the nearby orbit, x + w , arising from Taylor expanding theequations of motion for x + w and taking only linear terms in w . This leaves us with 2 N coupled linear differential equations for the motion of w ( t ) , albeit with coef fi cients thathave to be computed after each time step of the integration, as in general they dependon the current state of the orbit x ( t ) .There are several advantages to this approach, the most signi fi cant of which are re-lated to the main aims of any numerical integration procedure: Speed and accuracy.Since we have more or less exact equations for the evolution of the deviation vector, weare freed from a number of the constraints of the 2PM. Firstly, we no longer have toworry about the size of the deviation. Any deviation norm will be correctly evolvedby the variational equations, and since we are not simulating a second integration inthe phase space, we are at liberty to choose any reasonable deviation vector norm. Itis generally convenient (both mathematically and computationally) to set the norm ofthe deviation vector to unity. The fact that the size of the components of the deviationvector are no longer on the same order as the machine precision of our computers alsomeans that we can renormalise the deviation vector after every time step at no cost tothe precision of the integration. The second main advantage is that the variational equa-tions provide a much more precise estimation of the deviation vector’s evolution, andhence of the mLE itself [ ] . We are consequently able to use results found through thismethod with a great deal of con fi dence (having of course ascertained that the equationsare correct!), rather than having to take much care as with the 2PM that no numericalerrors are creeping in to the process.The major drawback to this method is that it requires the computation of the Hessianmatrix of the Hamiltonian function. In many cases (especially for low dimensional sys-tems and 1D lattices) this is not a problem, and explicitly fi nding the Hessian matrix andconsequently the variational equations is fairly straightforward. In other cases though, fi nding these second derivatives may involve a signi fi cant amount of work. This thenresults in the very pragmatic trade off between mathematical dif fi culty in the derivationChapter 1 Malcolm Hillebrand 23haotic Dynamics of Polyatomic Modelsof the variational equations, and the numerical challenges of the 2PM. For a given sys-tem, it makes sense to consider both options and choose the method that will pose thefewest practical problems in the course of computing the mLE, based on the complexityof the Hamiltonian and the computational power available.We note that there are other methods for computing the mLE from the deviationvector evolution, such as singular value and QR decomposition procedures (see for ex-ample [ ] and references therein) which are nevertheless closely related to the “standardmethod”, but in all the investigations carried out in our work we have used either the2PM or the variational equations. Before we abandon the deviation vector as nothing but a tool for computing the mLE, itis worth taking a closer look at what else we can learn about the chaoticity of the systemfrom it. It is well established that in a chaotic system, the deviation vector eventuallyaligns with a direction de fi ned by the mLE, corresponding to the most chaotic behaviourof the system [ ] . From this it is natural to look more closely at this direction of thedeviation vector, as by considering it we can fi nd useful information about the instabilityand the regions of active chaos [ ] . This concept also fi ts the intuitive understandingof the deviation vector as a measure of chaos – consider for example a one dimensionallattice. If at a particular site the component of the deviation vector remains very small,it implies that for that site the two nearby orbits yield very similar results. Thus thesystem is in some sense “locally” relatively stable at that site. Conversely, if at a givensite the deviation vector grows extremely fast, then we know that there the nearby orbitsare producing very different results, commensurate with chaotic behaviour. So we cantrack these components of the deviation vector, and use them to identify chaotic hotspotsin our systems [
23, 34, 35, 40 ] .We de fi ne the deviation vector distribution (DVD) as the normalised distribution ofthe deviation at each site, ξ i = w i + w i + N (cid:5) Nj = w j + w j + N , (1.21)where the deviation vector w has N components of displacement deviations δ q fol-lowed by N momentum deviations δ p . Thus ξ i measures the (squared) relative total“deviation size” in the i th direction, or corresponding to site i in a lattice. By comput-ing ξ through the course of the deviation vector evolution, we can build up a picture ofhow the instability spreads and moves through the lattice. Wherever the DVD is largestcorresponds to the highest concentration of chaoticity, and the regions where the DVDis small are regions of relative stability.The DVD has been used to study the phenomenon of breaking Anderson locali-sation in disordered one dimensional lattices [ ] , as well as to probe the underlyingmechanisms of chaotic wave packet spreading in both one dimensional [ ] and twodimensional lattices using various models, such as the Klein-Gordon and DNLS mod-els [ ] , as well as Hamiltonian mean fi eld [ ] and Hertzian and Fermi-Pasta-Ulam-Tsingou models [ ] . These studies have provided good evidence that the DVD can bea useful measure in understanding chaotic dynamics in a variety of systems.24 Malcolm Hillebrand Chapter 1 hapter 2Modelling DNA Now we turn to discussing the material of the research, through a brief introduction tothe basics of DNA and a summary of the progress of the most relevant DNA models.This chapter will outline the structure of DNA, as far as the dynamics are concerned,and discuss the DNA model. This will include an overview of the ways that DNA hasbeen modelled in the past, with an eye to how the models used in this investigation farein reproducing experimental results.
The rich biochemical structure of DNA is deserving of several textbooks’ worth of dis-cussion (e.g. [
1, 42, 43 ] ), but fortunately for this speci fi c application many of the biolog-ical details can be omitted. Indeed, the general model used for our studies, the Peyrard-Bishop-Dauxois (PBD) model, makes some drastic simpli fi cations for the purposes ofaccurately capturing the dynamical behaviour [ ] . The model itself will be dis-cussed in the next section, but here we will fi rst look at the basic structure of DNA. Fora thorough review of DNA in the context of nonlinear dynamics, Peyrard’s 2004 reviewpaper is highly recommended [ ] .Double-stranded DNA exists in the fairly well-known double helix structure. Fig-ure 2.1 illustrates this basic formation. The fundamental building blocks of the strandare tripartite nucleotides – comprised of a phosphate group, a sugar ring, and a nitroge-nous base. These monomers form the overall DNA polymer, the double-stranded DNAhelix. The strands (made up of the phosphate and sugar backbones) are bound togetherby hydrogen bonds that form between the bases of the nucleotides. The discovery ofthis base-pair structure of DNA in 1953, and the relative simplicity of their nature, wasgroundbreaking in the fi eld of biochemistry [ ] . There are four bases found in thesenucleotides: adynine, guanine, thymine, and cytosine. These bases fall into two groups;the purines (adenine and guanine) and the pyramidines (thymine and cytosine). Purineshave a double-ringed structure, while the pyrimidines have a single-ringed structure.One of the more surprising aspects of Watson and Crick’s discovery was that these basesform pairs very strictly with one complementary base. Rather than either purine bond-ing with either of the pyrimidines, they discovered that the only two possible naturallyoccurring base pairs are adenine-thymine (AT) and guanine-cytosine (GC). These basepairs turn out to be the essence of the genetic code that governs characteristics of livingorganisms, with the sequence of AT and GC pairs encoding essential information. They25haotic Dynamics of Polyatomic Modelsa particular DNA sequence “instructs” proteins in their development, de fi ning the se-quence of amino acids in the protein by the formation of a particular polypeptide chain,which then undergoes protein folding. In this process of communication between DNAand protein, the DNA strand opens up, under the in fl uence of the RNA polymeraseenzyme, to copy the base sequence of the relevant gene to a strand of RNA [ ] . In par-ticular, when transcription occurs, the RNA polymerase moves along the DNA strand,opening up the double helix to expose the base sequence to the RNA molecule. Speci fi cbase pair sequences are known to denote the transcriptions start and end sites, ratherthan amino acids, and the RNA polymerase moves between these two set sequences, al-lowing a single gene to be transcribed at a time. This RNA molecule will then providethe genetic information for the polypeptide chain, the sequence of amino acids, and thusthe fi nal protein product.What all this means is that the opening of base pairs, indeed base pairs opening wideenough for the RNA copying to occur, is an essential part of the everyday activity ofDNA molecules (even though these openings are created by an enzyme). This activityof base pairs is precisely the opportunity for the physicists, and particularly the dynam-icists, to provide input. The opening and closing of base pairs is a fundamentally dy-namical process – while in the case of DNA transcription it requires an external agentto perform the zipping and unzipping, similar openings occur entirely spontaneuslyin the form of thermal fl uctuations. These thermal fl uctuations are termed “bubbles”,and can take various forms; breathers in a lattice, long-lived stable openings, travellingwaves similar to transcription bubbles, erratic openings, or any number of ways [ ] .Wherever there is energy in the strand, these openings can occur. The tendency of theseopenings to even exist with very large separations, sometimes larger than 7Å, or threetimes the natural length of a hydrogen bond [
47, 54 ] , means that nonlinear aspects ofthe dynamics will certainly be explored. While the dynamics of the stacking force andhydrogen bonds may be near linear at short distances, when probing the large displace-ment regimes, nonlinearity becomes important, as will be seen in the following sectionon the different models used to describe DNA dynamics.A particularly important property of DNA, from a biological as well as dynamicalstandpoint, is the phenomenon of denaturation, the complete separating of the doublehelix [ ] . In some sense this is an extreme case of bubbles forming, and encompass-ing the entire strand. The experimental approach to measuring denaturation is perhapssurprisingly simple for such microscale work. Due to the electronic structure of thenucleotide bases, there is a dramatic increase in the absorption of ultra-violet light withwavelength of 260nm when the base pairs separate, in the region of a 30–40% increase [ ] . This allows the fraction of open base pairs to be quite accurately measured experi-mentally – as the solution is slowly heated, the amount of light absorbed can be tracked,and this corresponds precisely to the fraction of open base pairs. Experiments have beencarried out on a variety of synthesised DNA strands, initially on homopolymers, strandswith a single base pair type, but soon expanding to heterogeneous strands under vari-ous conditions (see [
47, 56, 57 ] and references therein). This thermal denaturation isa remarkable process, since as the temperature of DNA in a solvent is increased, thereis an abrupt increase in the fraction of base pairs that are separated [
56, 57 ] . The typeof solvent, salinity and various other factors can slow the process somewhat, but theoverall picture is very clear, of a sudden change. Figure 2.2 shows an illustration of theultra-violet absorption curve as the temperature increases.Chapter 2 Malcolm Hillebrand 27haotic Dynamics of Polyatomic Modelschecked against is the denaturation temperature – as discussed above, there are experi-ments fi nding fairly precise melting temperatures for different base pair compositions.It is thus possible to model DNA by declaring that the most important aspect is whetheror not a given base pair is open or closed, and then to reproduce some physical measure-ments of critical melting temperatures.There are of course several ways of approaching this task, including a signi fi cantbody of work focussing on solitons as explanations for open base pairs [
64, 65 ] , as wellas the self-consistent phonon method [ ] . But since here we are primarily concernedwith discrete models, and of course particularly the PBD model of DNA, the discussionwill be focussed on the precursors to this model and the model itself. The starting point for this type of model was a relatively simple Ising-like model, thePoland-Scheraga model developed in 1966 [ ] . The Poland-Scheraga model treats eachbase pair as a lattice site that can be either open or closed, in the same way as the ferro-magnetic Ising model has two-state magnetic dipole moments of ±
1. This type of modelhas also been called a “helix-coil” model, as the two states (closed and open) correspondto sections of the DNA molecule that are in a helix form (closed base pairs) and a coil(open base pairs, the helix is separated) [ ] . The essential idea of this model is to pro-vide a rate ratio for the helix to grow or to shrink. If the helix grows, then more basepairs close, and this is often called zipping. Correspondingly, should the helix shrink,more base pairs are opening and the process is known, unsurprisingly, as unzipping. Themain idea of this model is to accurately reproduce melting curves, essentially reportingthe correct numbers of open base pairs. It has been moderately successful in this [ ] , and is still relevant for testing more recent results against [ ] . The model wasalso improved by considering the problem from a random walk point of view, and par-ticularly self-avoiding random walks, to describe the melting transition more accuratelyby Fisher [ ] . However, despite this success it does not consider the dynamics at all,beyond considering the movement of bubbles in the double strand, and so it does nothave the same scope as a more detailed model would have. The Peyrard-Bishop (PB) model arises from developing the Ising-like PS approach, andhas been successfully used to reproduce experimental results of the thermodynamicproperties of DNA [
44, 45 ] . Essentially, this model takes the binary, discrete, variableat each site of a base pair being open or closed, and now allows this quantity to varycontinuously. This moves the model from purely describing the presence of bubblesand denaturation to actually modelling the nonlinear dynamics of the base pair fl uctuta-tions. This of course can be used to study the open / closed nature of base pairs, throughchoosing thresholds for open base pairs (this will be discussed in detail in Chapter 5),but inherently allows a more thorough investigation of the nonlinearity of the base pairinteraction as pointed out by Prohofsky [ ] . A signi fi cant advantage of the model isthat it presents a Hamiltonian dynamical system, which allows for ef fi cient and effectiveanalysis of the dynamics in general, and the chaoticity in particular.In the basic version of this PB model, the effects of the left-right ordering of basepairs on the dynamics is not considered. So an AT base pair is functionally identical toChapter 2 Malcolm Hillebrand 29haotic Dynamics of Polyatomic Modelsa TA base pair, and similarly GC and CG. The reduction from the complex molecularDNA structure to a straightforward lattice is illustrated in Fig. 2.3. The overarchingidea is to reduce the biological complexity, which is tremendously hard to model infaithful detail, to a simpler structure that is amenable to modelling through commonmolecular dynamics techniques. In this regime, the biological complexities will be ex-pressed through carefully chosen and fi tted parameters, optimising an approximationrather than attempting to exactly reproduce details. → → → F IGURE
The simpli fi cation process for the PB model of DNA. The ini-tial double helix is “ fl attened” into a straight line, with the twisting forcesaccounted for in the potential. Next, the precise details of AT / TA andGC / CG base pair orientation are ignored, resulting in a simple one di-mensional lattice of two possible element types – AT and GC base pairs.The fi nal simpli fi cation is to ignore even this distinction, and create a ho-mogeneous lattice of base pairs, where the type is neglected. We will seehow more developed models step these simpli fi cations back, and take moredetail into account. Now that the structure is simpli fi ed to a comprehensible lattice, we can move on todescribe the potential used to model its dynamics. The fi rst, fundamental, Hamiltoniandeals with the bases themselves as the elementary entities, and so the dynamical variables u n and v n are the displacements from equilibrium of the n th base in each strand. Here,we have not yet moved to the base pair as the fundamental object – so these u n and v n displacements describe the positions of the actual bases themselves, in each backbonestrand. In these variables, for an N -base pair double strand the Hamiltonian reads [ ] H = N − (cid:1) m (cid:6) ˙ u n + ˙ v n (cid:7) + D (cid:6) e − a ( u n − v n ) − (cid:7) + k (cid:8)(cid:6) u n − u n − (cid:7) + (cid:6) v n − v n − (cid:7) (cid:9) , (2.1)where we have assumed that the masses of each nucleotide have the same value m . Thetwo potential terms model the on-site potential for each base pair and the stacking in-teraction respectively. The on-site potential is modelled by a Morse potential [ ] . Thishas been used extensively in modelling of chemical potentials [
49, 73 ] , and accuratelyimitates the combination of hydrogen bond attraction and phosphate interactions, pro-viding a simple but effective way of approximating the base pair forces. Crucially, it hasa sharp increase for negative displacements [ see Fig. 2.6(a) ] , which is physically desirable– bases should very rarely come into close contact. The potential has a minimum at zero30 Malcolm Hillebrand Chapter 2haotic Dynamics of Polyatomic Modelsdisplacement, corresponding to the desired behaviour that the base pairs should returnto precisely zero displacement. It also has a plateau at large displacements. This meansthat as the hydrogen bonds stretch signi fi cantly, the force pulling the base pair back toequilibrium diminishes, exactly as occurs when the bonds stretch too far and break. Inthe case of a simple harmonic potential, which has many of the desired characteristics– the minimum at zero, and prohibiting small displacements – the large stretching ofbase pairs would not be possible, as there is no plateau for large displacements. Thismakes the Morse potential a qualitatively appealing option for describing the base pairinteraction. Here, the parameters of the Morse potential (the depth D and characteristicwidth a ) do not depend on the site at all, meaning that the Morse potential representsan average potential for both base pair types. This is an obvious simpli fi cation, but itfacilitates analytical treatment of the system [
44, 47 ] .The other term of the potential simply approximates the stacking interaction be-tween adjacent bases as a straightforward harmonic coupling. As we will see, this is ade fi nite oversimpli fi cation, but as with the average Morse potential this allows for rela-tively simple mathematical analysis.Returning to the discussion of the strand dynamics, a change of variable simpli fi es theanalysis signi fi cantly, as what we are primarily interested in is the displacement betweenbase pairs – thus, we are concerned with out-of-phase motion of the two strands. If wedenote the in-phase motion by the variable x n , and out-of-phase motion by y n , i.e. x n = ( u n + v n ) (cid:8) y n = ( u n − v n ) (cid:8) y n to provide the dynamical infor-mation about the displacements between the base pairs. This then gives the Hamiltonianfor the out-of-phase displacements after rescaling the parameter a as H = N − (cid:1) n = p n m + D ( e − ay n − ) + k ( y n − y n − ) . (2.2)It is important to note that after this change of variables, the variable y n does notexactly represent the physical displacement between base pairs – in order to recoverthe actual displacement, the factor of (cid:8) fi nes adynamical system which can be studied numerically and theoretically. In [ ] , wherethe model was introduced, dynamical quantities were studied using statistical physicstechniques, and particularly the transfer integral method (see [
47, 63 ] and referencestherein). One main result produced was to estimate the mean stretching of the bondsbetween base pairs, depending on the chosen coupling constant. Figure 2.4 shows thisresult, and most signi fi cantly shows the sharp increase in the mean stretching near thephase transition point. While this result does not explicitly relate to the fraction ofopen base pairs and thus denaturation, it certainly gives an indication that the modelapproximately exhibits the desired behaviour at the transition. While this fi rst model of the nonlinear dynamics of DNA is useful for studying variousproperties, and admitted extensive theoretical treatment as a straightforward model withChapter 2 Malcolm Hillebrand 31haotic Dynamics of Polyatomic Models − . . . . . . . y ( ˚ A) . . . . . . . V ( e V ) ATGC (a) y n . . . . . . . y n − . . . . . . . V ( e V ) . . . . . . . (b) F IGURE (a) The Morse potential for an AT base pair (blue) and a GC basepair (green). We see the effect of the coef fi cient D n , dictating the heightof the plateau. The higher plateau of the GC base pair is indicative of thestronger attractive on-site force. The constant a n affects the narrownessof the well. (b) The coupling potential between two base pairs. If the dis-placements are equal (i.e. no relative displacement), then the coupling po-tential contributes nothing, as would make sense. As displacement grows,the coupling potential increases exponentially. variety of DNA sequences. Through allowing for heterogeneity in the Morse potential,allowing the parameters D and a to depend on the base pair type, and fi tting to exper-imentally determined denaturation curves, parameters have been found that producegenerally accurate results for a diverse array of base pair sequences [ ] .The parameters chosen by Campa and Giansanti [ ] are not signi fi cantly differentto the original proposed PBD parameters, but the addition of heterogeneity allows fora much more precise reproduction of experimental results. The PBD potential, nowconsidering all these factors, is V = N (cid:1) n = D n ( e − a n y n − ) + k ( + ρ e − α ( y n + y n − ) )( y n − y n − ) . (2.4)Here we see the combination of the heterogeneous Morse potential (note the subscriptson the parameters D n and a n ), and the nonlinear coupling. The parameters of Campaand Giansanti have become accepted as near optimal [
54, 60, 75 ] , and are the parametersused throughout this investigation. The parameters are given in Table 2.1.Parameter Value k / Å ρ α − D AT D GC a AT − a GC − T ABLE
Optimised parameters for PBD simulations.
34 Malcolm Hillebrand Chapter 2haotic Dynamics of Polyatomic ModelsFigure 2.6 shows the potential energy contributed by each term. Figure 2.6(a) illus-trated the Morse potential, for the two sets of parameters. We see that the larger welldepth D for the GC base pairs naturally results in a higher plateau, corresponding to thegreater amount of energy required to break the triple hydrogen bonds. The characteris-tic width of the GC potential is also narrower, with the well consequently being muchsharper than the AT case. So while it requires more energy to break the hydrogen bonds,the physical displacement required to reach the fl attened region is in fact lower than itis for an AT base pair. The coupling potential is inevitably more dif fi cult to visualise, asit depends on both neighbouring base pairs. In Fig. 2.6(b) we see the coupling energy asa function of the displacement of two base pairs. This plot shows basically the expectedbehaviour – a valley of zero energy when the displacements are equal, and so in relativeequilibrium, and nonlinear increases as soon as the displacements differ. While the heterogeneous PBD model with a nonlinear coupling is an effective modelfor DNA, being used in studies of various properties of DNA (e.g. [ ] ),there are still dynamical aspects of DNA that are not explicitly addressed in this po-tential. One such aspect is the dependence of the stacking interaction on the base pairsequence. It is known that the stacking interaction changes depending on the actual basesinvolved, and the stacking energy of neighbouring base pairs can vary signi fi cantly [ ] . This stacking interaction depends not only on the base pair type (AT or GC), butalso on the exact ordering of the bases. So now whether a base pair is AT or TA be-comes signi fi cant, as this affects the stacking force. This is actually a large part of themotivation behind introducing a sequence-dependent term to the stacking force, as thedifference between CC / GG and CG / GC stacking energies is signi fi cant [ ] . Underthe PBD model’s framework, these two sequences are considered precisely equivalent –the details of CC / GG or CG / GC are completely lost in the fi nal simpli fi cation step.However, if we wish to accurately model these different behaviours, we need to accountfor this sequence dependence. It is noteworthy that the difference in melting tempera-tures between these two sequences – one where one strand consists of purely guaninebases and the other of cytosine, and one where both strands alternate – is around 22K,which is signi fi cant (see Fig. 2.7 for an illustration of these two sequences).In 2009, an extended PBD model (ePBD) was proposed, that took the sequence-dependence of the stacking interaction into account [ ] . Here the stacking potential isadapted very slightly to account for this dependence by allowing the constant k to varyaccording to the sequence. The updated stacking potential is then W = k n (cid:10) + ρ e − α ( y n + y n − ) (cid:11) ( y n − y n − ) , (2.5)where the index n on the constant k allows for the coupling between every stacked basepair to be different. This means that k of course takes on ten possible values, one for eachpossible step from one base pair to another. To fi nd values of k n that accurately modelboth homogeneous and heterogeneous sequences, the authors of [ ] fi tted the meltingtemperatures obtained by the simulations to experimentally observed melting tempera-tures for speci fi c sequences. Use was also made of previously available stacking constantratios, and the known stacking energies. The optimised values of k n are presented inChapter 2 Malcolm Hillebrand 35haotic Dynamics of Polyatomic Models (a) (b) F IGURE (a) A sequence with guanine bases on the left strand (blue), andcytosine on the right (orange). This corresponds to the CC / GG case,where the stacking energy and melting temperature is much lower. (b)A sequence with alternating guanine and cytosine bases. This (CG / GC)is a much more robust arrangement, as the stacking energies are higher,resulting in a signi fi cantly higher melting temperature. C G A TC 0.0192 0.028 0.025 0.0229G 0.0249 0.0192 0.019 0.0226A 0.0226 0.0229 0.0228 0.023T 0.019 0.025 0.0193 0.0228T
ABLE
Numerical values of the sequence-dependent stacking constants k n in eV / Å , following the base sequence of the top strand. The row denotesthe fi rst base, and the column denotes the second base in the step. So k AG ,the stacking constant for the step from A to G, is found in the third row(A), and second column (G). Fig. 2.8 as they are given in [ ] . The numerical values are given in Table 2.2. We seethat there is a noticeable difference between the extremes, and also that the CC / GG –CG / GC discrepancy is illustrated in the vastly different stacking constants. The othernotable fact is that the stacking constants are almost all lower than the value used in thePBD model, suggesting that overall the base pairs are slightly more weakly coupled inthe extended model.The purpose of all this effort in extending the model is to more accurately reproducephysical results, particularly melting temperatures. The ePBD model has been success-fully shown to match experimental melting temperatures for several sequences, as wellas to produce other results that are biologically consistent and useful [
54, 75 ] . In furtherchapters, the behaviour of this extended model will be compared to the simpler PBDmodel, both in chaotic dynamics and in physical properties such as melting tempera-tures and bubble formation.36 Malcolm Hillebrand Chapter 2 hapter 3The Alternation Index and itsProbability Distribution Before actually getting into the dynamical study of DNA, we will take a rather math-ematical side path in order to prepare an important tool for the dynamical study. Inthe course of studying the effects of heterogeneity of DNA, we need to have some un-derstanding of what it is precisely that makes one DNA sequence “heterogeneous”, andanother sequence “homogeneous”. As discussed in some depth in the preceeding chap-ter, we know that there are two kinds of base pair. It follows quite naturally that theheterogeneity, and disorder, in this DNA model comes from the combination of thetwo types of base pair.It is intuitively quite clear that a sequence of only AT base pairs should be consideredto be homogeneous. There is no disorder whatsoever, and no distinction between givensites in the lattice (here the one dimensional lattice is of course the DNA double strand).Similarly for a sequence of purely GC pairs, while it would be dynamically different tothe AT sequence, it is nonetheless very homogeneous. The extension of this natural stepis to introduce heterogeneity based on the mixing of these base pairs in a single sequence.So a sequence comprised of say fi fty GC base pairs and thirty- fi ve AT base pairs wouldbe considered heterogeneous and disordered – it is now signi fi cant which lattice site weconsider, as it will have different dynamical properties depending on what base pair typeit is hosting.However, this still leaves us a little short of a complete description of heterogene-ity. We can separate homogeneous and heterogeneous sequences, but we would verymuch like to be able to compare the heterogeneity of two sequences. In this chapter, wewill introduce the idea of the alternation index as a measure of this heterogeneity, andpresent the mathematical machinery used to derive a completely rigourous probabilitydistribution for it. This chapter is based on the work in [ ] . We will consider some particular DNA sequences to get some idea of what we are look-ing for. But fi rst, in order to make things clearer for this section, in all representationsAT and GC base pairs will be represented by black and white circles. This is moti-vated by the fact that for now we are forgetting all the DNA details – all we want tounderstand is the mixing of base pairs. So a simple and clear representation serves our38haotic Dynamics of Polyatomic Modelspurposes neatly. Filled black circles or beads will represent GC base pairs, and emptycircles or white beads will represent AT base pairs. Figure 3.1 shows one particular ar-rangement of base pairs. In this con fi guration there are ten base pairs, four AT and sixGC. Throughout this investigation we will quantify the ratio of AT to GC by referringto the percentage of AT base pairs in the chain. So in this case, where four of the tenbase pairs are AT, the AT percentage P AT = IGURE
A schematic representation of a DNA sequence. Here we haveten base pairs, four of which are AT (white beads), and six of which areGC (black beads). Thus the AT percentage P AT = We can take the same base pairs and mix them around a bit, and get a large numberof different arrangements (such as in Fig. 3.2). These arrangements are referred to asdisorder realisations, when linking back to the idea of DNA as a one dimensional lattice.So for a given number of base pairs in the sequence, and a given AT percentage, thereare multiple disorder realisations possible.F
IGURE
Here we still have ten base pairs, four of which are AT, and six ofwhich are GC, but now there is some mixing – there are three distinct clus-ters of GC base pairs, and three distinct clusters of AT base pairs. So withthe same AT percentage and number of base pairs, we have a sequence thatwe would expect to be meaningfully different to the previous arrangementin Fig. 3.1. There is a more thorough mixing in this arrangement, whichwe consider to be more homogeneous.
The question now becomes what is different between these realisations? We knowthat the sequence of base pairs is genetically signi fi cant, which is encoded in the factthat the PBD model accounts for base pair types in the Morse potential. It is also nat-ural than when considering heterogeneity the AT percentage is signi fi cant – of coursesequences with different compositions will yield differing dynamical behaviours. Buteven within the same AT percentage, as these two realisations have illustrated, quite dif-ferent arrangements are possible.In order to quantify these differences in arrangements we will introduce the alter-nation index , α . This alternation index simply counts the number of times the DNAsequence changes base pair type, as it is traversed. So for instance, consider Fig. 3.1.Starting from the left, we fi nd that from the third to fourth base pairs, the base pairtype switches from GC to AT. So this is the fi rst alternation. Then from the seventh tothe eighth base pair, once again we have an alternation back from AT to GC. This is thesecond and fi nal alternation. We now also realise that boundary conditions play a signif-icant role in the alternation index – if we were to impose periodic boundary conditionson the DNA molecule then we would end up where we started. This will be discussedin more detail later where we do make use of periodic boundary conditions.Applying the same counting process, we see that the sequence in Fig. 3.2 has six alter-nations. This tells us that the more “well-mixed” a sequence is, the higher the alternationChapter 3 Malcolm Hillebrand 39haotic Dynamics of Polyatomic Modelsindex. We can also see that for this case, of four AT base pairs mixed between six GCbase pairs, there is a maximum possible value of α of eight – this corresponds to every ATbase pair being isolated from each other and dispersed among the GC mass. Of course,whether the majority of base pairs are GC or AT is irrelevant – if we had four GC basepairs and six AT base pairs, the maximum value of α would still be eight. Similarly,we can see that the minimum number of alternations in a heterogeneous sequence (so atleast some of both AT and GC base pairs) is one, or two if we impose periodic boundaryconditions. So we now have a way of quantifying how heterogeneous a sequence is – thesmaller the alternation index, the more “chunky”, and hence heterogeneous, it is.What we would like to do then, is to study the dynamics of DNA based on this mea-sure of the heterogeneity. This is the motivation for introducing the alternation index.So having de fi ned it, we are in principle in a position to study the effects of heterogene-ity – we can vary both the AT percentage and the alternation index, and see what resultswe fi nd. However, we would like to perform a more informed investigation. In partic-ular, we would like to know which values of α are more probable, and which values areextremely improbable. Ideally this would be in the form of a probability distributionfunction (PDF), where the exact most likely value can be found, and the probabilityof extreme cases can be quanti fi ed. As it turns out, this is a surprisingly dif fi cult prob-ability distribution to fi nd. Especially when imposing periodic boundary conditions, fi nding a clean way of calculating the number of possible combinations is extremelytricky. There have been calculations making use of Markov chains to approximate theprobabilities of different arrangements of base pairs in non-periodic sequences [ ] .However, in our periodic case there turns out to be a mathematically beautiful, and mostimportantly exact, way of calculating the number of different arrangements, and hence fi nding the probability of each. We are about to dive in to a mass of mathematics, which we will ultimately use to de-scribe the probability of various values of α , in a sequence of N AT AT base pairs and N GC GC base pairs. But before being subjected to pages of group theory, let us see themotivation for using all this machinery. Our aim is essentially simple: Take a num-ber of base pairs, N , of which N AT are AT pairs and N GC are GC. If we were to pick acompletely random arrangement of these base pairs, we would like to know what theprobability is for some particular number of alternations α . Here we are imposing pe-riodic boundary conditions, so we additionally require that the fi rst and last base pairsare neighbours.This distribution depends on two variables, N AT and N GC (note N = N AT + N GC ). Itwould be natural to consider a binomial distribution, or some kind of general statisticaldistribution. However, as it turns out trying to select only sequences with the samealternation index is a rather dif fi cult task in this approach. We instead consider theproblem from a more fi rst principles standpoint: What we actually need is the numberof possible arrangements of base pairs with a given alternation index, for every possiblealternation index. From this we could simply normalise the distribution and we wouldhave our desired PDF.A crucial fact to recognise in this problem is that our dif fi culty (restricting sequencesto a given alternation index) in fact allows for a simpli fi cation which leads to a new type40 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Modelsof solution. Our periodic DNA sequence is, for the purposes of this study, equivalentto a necklace of black and white beads. So associating black beads with GC base pairsand white beads with AT base pairs, we can address this as a classical combinatoricalquestion, rather than a purely statistical one. With this in mind, we can recognise thatactually a necklace of black and white beads can be written as a necklace of containers ofbeads, with each container holding a certain number of base pairs. Let us consider someexamples to illustrate this idea.Figure 3.3 shows the periodic representation of the sequence shown in Fig. 3.2. Thisstep is fairly straightforward – we are simply connecting the boundary beads to form aloop. F IGURE
The necklace of black and white beads corresponding to the ar-rangement in Fig. 3.2. The sequence has just been looped around so thatthe fi rst and last beads are neighbouring. So having this necklace, we can represent it through a sequence of containers ratherthan beads (Fig. 3.4). This necklace of containers takes each homogeneous chunk ofthe sequence and compresses it to a single element, which gives the colour of the beadscontained in that chunk as well as the number of beads. So the two white beads compressto a single white container labelled with “2”. Note that even the single beads are nowrepresented through containers, even though the number of beads contained is simplyone. F IGURE
The representation of the necklace in Fig. 3.3 as a necklace ofcontainers. The number on each container corresponds to the numberof beads “contained” within it. So each homogeneous sequence of beadscollapses to a single element, which then represents that sequence as thenumber of beads in it. The containers are grouped in pairs to show thealternations more clearly – each pair of containers corresponds to two al-ternations. Note also that we could expand this necklace of containersunambiguously into the original necklace of black and white beads.
The important thing to note here is that the number of containers is equal to thealternation index. This is where we are taking advantage of the fact that the actual num-ber of beads in each container does not affect the alternation index. So as long as theChapter 3 Malcolm Hillebrand 41haotic Dynamics of Polyatomic Modelscorrect number of beads are somehow distributed among these containers (of coursewith at least one bead in each), then we are satisfying the criteria of set numbers of blackand white beads and a set number of alternations. This reformulation of the problem isactually signi fi cantly easier to approach.The approach we use to count the possible number of bead distributions for eachvalue of the alternation index is based in Pólya counting theory [ ] . In the follow-ing sections we will introduce the necessary tools to solve the problem – group theoryand especially cyclic groups, fi nite power series, and the crux of the theory, Burnside’sLemma and Pólya’s Enumeration theorem with an extension for bipartite sets. Let us dive into some mathematics, primarily from the fi eld of group theory, to establishthe requisite tools to address the rather complex problem of counting necklaces subjectto constraints. These topics are covered in most algebra texts; see for instance [
84, 85 ] .The discussion here is based largely on the contents of those works. Mathematical groups are well-known abstract entities, consisting of a set with an op-eration. So formally, given a set A with a binary operation · , they are considered incombination to be a group if the following hold [ ] :1. ∃ i ∈ A | ∀ a ∈ A i · a = a · i = a (Existence of an Identity).2. a · ( b · c ) = ( a · b ) · c ∀ a , b , c ∈ A (Associativity).3. ∀ a ∈ A ∃ b ∈ A | a · b = b · a = i (Existence of Inverses). De fi nition 3.2.1 (Group) . Any set with a binary operation satisfying properties 1-3 is agroup.
A particularly useful, and hence well-studied, branch of group theory is that of per-mutation groups, where the operation on the set of items is a permutation of the groupitems. These permutations are rearrangements of the elements of the set. So the ele-ments of the set C = {
1, 2, 3 } could be written in six different ways: 123, 132, 321, 312,231, 213. Each of these ways of writing the elements is a permutation of the set C .The symmetric group on a set A is the set of permutations of A . De fi nition 3.2.2 (Symmetric Group) . The symmetric group on a set A, S A is de fi ned asS A = { σ : A → A | σ a permutation } This group is important for a number of reasons; not only is it useful in itself, but alsoevery group whose elements interchange set elements in any way is a subgroup of thesymmetric group, as all possible interchanges are contained in the set of permutations.Within these permutations, it is useful to introduce the notion of cycles . A cycle, asthe name suggests, is a permutation that is cyclically periodic under repeated application.We say that a given permutation is an n -cycle if after n applications of the permutationthe set elements have returned to their original positions.42 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Models Example 1.
Consider the set A = {
1, 2, 3, 4 } , and a permutation σ that swaps the fi rst andthird elements and swaps the second and fourth elements, i.e. σ ( ) = , σ ( ) = , σ ( ) = , σ ( ) = . In this example, we see that clearly after two applications of σ , the set returns to itsoriginal state: σ ( σ ( A )) = A . We also see that there are two distinct cycles occurringwithin σ . The fi rst cycle sends 1 → →
1, while the second cycle sends 2 → →
2. We denote these in cycle notation as ( ) and ( ) respectively. This is equivalentto calling them ( ) and ( ) , as these mean the same thing – interchanging 1, 3 and2, 4 respectively. What is also clear is that each of these cycles are 2-cycles. After twoiterations, each element is back where it started. Thus, we can write σ as the compositionof two 2-cycles, σ = ( )( ) . De fi nition 3.2.3 (Cycle) . A permutation σ ∈ S A is called a cycle if there are distinct ele-ments in the set A = { a , a , . . . , a m } , such that1. σ ( a i ) = a i + for i ∈ {
1, 2, . . . , m − } .2. σ ( a m ) = a .3. σ ( a ) = a for any element of A / ∈ { a , a , . . . , a m } . So any n − cycle can be written as ( a a . . . a n ) . In a natural way, we say two cyclesare disjoint if they do not interact at all. So if we have two disjoint sets { a , a , . . . , a m } and { b , b , . . . , b k } , then the cycles σ = ( a a . . . a m ) and σ = ( b b . . . b k ) are alsodisjoint.With these de fi nitions, and using the fact that for this problem we are only concernedwith fi nite groups, we can state a theorem that will be very useful in further considera-tions. Theorem 1 (Cycle Decomposition Theorem) . For a fi nite set A, any permutation σ ∈ S A can be written as a composition of a fi nite number of pairwise disjoint cycles. So σ = ( a a . . . a m ) . . . ( a k a k . . . a km k ) . (3.1) This decomposition is unique, up to ordering of the cycles.
Being able to write permutations as a product of distinct cycles turns out to be crucialfor the machinery in the coming sections, as we will soon see.
In order to understand the next tool in the necklace counting process, we need a basicnotion of a group action . A group action of some group G on a set B is a mapping ofthe elements of G to the symmetry group S B . This mapping is in fact a homomorphismthat identi fi es the group elements of G with permutations in S B , in some sense allowingcommunication between the group and the set. We label the group action of G on B as A G on B .For instance, let us take G as the cyclic group of order four C , which is the groupof rotations for sets of four elements. Then the action of G on a set of four points P = { p , p , p , p } would be to rotate the points through their various positions. So weChapter 3 Malcolm Hillebrand 43haotic Dynamics of Polyatomic Modelshave the group action A G on P , where we abbreviate notation by saying that the actionof a group element g on an element p is g p , rather than A G ( g )( p ) . If g is the elementof G that performs a single rotation, then g p i = p ( i + ) %4 . Similarly, g p i = p ( i + ) %4 andso on. This is the idea of a group action, where the group elements act on the set toproduce permutations of the set.Having this group action, we can de fi ne some important concepts. The fi rst we willaddress is the idea of fi xed points . De fi nition 3.2.4 (Fixed Points) . Given a group G acting with action A G on a set B , theset of fi xed points of an element g ∈ G is denoted by
Fixed g = { b ∈ B | g b = b } .This simply de fi nes the set of all elements that are left invariant by the action of thegiven element. For this reason this set is also referred to as the invariant set of g .In the context of counting distinct necklaces, there is an intuitive link here betweenthe ideas of “elements that are unchanged by a permutation” and “necklaces that are thesame when permuted”. This will of course become clearer in due course.The next subset to consider is the orbit of some element of the set. De fi nition 3.2.5 (Orbit) . The orbit of an element b ∈ B under the action of a group G isthe set of all elements that G maps b to
Orbit b = { g b | g ∈ G } .While this may be less clear at fi rst thought than the idea of fi xed points, it is alsonot an intuitively strange idea. It is somewhat analagous to the range of a conventionalmathematical function.Thirdly, we will consider the subgroup called the stabiliser of an element b . De fi nition 3.2.6 (Stabiliser) . The stabiliser of an element b ∈ B under the action of agroup G is the group of all elements of G that leave b unchanged.
Stab b = { g ∈ G | g b = b } This is of course similar to the set of fi xed points – where the fi xed points are set ele-ments unchanged by a given group element, the stabiliser is comprised of group elementsthat are idempotent on a given set element.The fi nal concept to introduce here is the cycle index [ ] . The cycle index encodescombinatorical information in the form of a polynomial, and in particular the coef fi -cients of each term in the polynomial provide information about the action of groupson sets. When we make use of the cycle index for counting distinctly coloured necklaces,its utility will become more readily apparent. For now though, we settle for a de fi nitionand some illustrative examples. De fi nition 3.2.7 (Cycle Index) . Given a fi nite group G, acting on a fi nite set B with ac-tion A G . If the cardinality of B is n, then the cycle index of the group action is given byZ G ( y , y , ..., y n ) = | G | (cid:1) g ∈ G y c ( g ) y c ( g ) ... y c n ( g ) n .44 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic ModelsIn this de fi nition, the c i refer to the number of i − cycles in the (unique) disjoint cycledecomposition of the permutation A G ( g ) . So for example if there are twelve cycles oflength seven, c = [ ] . • The cyclic group C n , which is used to account for rotational symmetries: Z C n ( y , y , ..., y n ) = n (cid:1) p | n ϕ t ( p ) y n / pp . (3.2)Here ϕ t ( p ) represent Euler’s totient function, counting the number of integersless than p , that are also relatively prime to p . • The dihedral group D n , which adds re fl ectional symmetries to the rotations of thecyclic group: Z D n ( y , y , ..., y n ) = (cid:12) n (cid:5) p | n ϕ t ( p ) y n / pp + y y ( n − ) / , for n odd, n (cid:5) p | n ϕ t ( p ) y n / pp + (cid:13) y y ( n − ) / + y n / (cid:14) , for n even.(3.3)In our bead-counting context, the two cases for odd and even n correspond to theslightly different meaning of a re fl ection for necklaces with odd or even numbersof beads. This is discussed in more detail later, when the dihedral group is used forthe counting of distinct necklaces.We have now introduced, albeit brie fl y, the fundamental building blocks needed inthe process of counting distinct necklaces. There are some other concepts that will bede fi ned as they come up in the natural fl ow of outlining the counting theory, rather thanfront-loading this chapter with an enormous mass of abstract ideas without context. With these tools, we turn back to the discussion of different ways of arranging beads ina necklace. To complete the change from DNA sequences to pure combinatorics, wewill refer to these arrangements of beads as colourings . This is a very natural changeof terminology – after all, we can consider the base necklace to be an empty templateof N blank beads. Then to create some particular necklace we colour N AT beads withwhite, and N GC with black. So by colouring the “blank template”, we have formed ournecklace of beads.Let us recapitulate the aim in this new language: We wish to count the number of distinct colourings of the necklace subject to some restrictions, with colours of black orwhite. We do now need to spend some time thinking about these restrictions. Firstlythere are the natural restrictions mentioned previously, that we have N AT white beadsand N GC black beads, and a total of α alternations in the chain. But there are further,more subtle restrictions we need to impose. For instance, in principle the two neck-laces in Fig. 3.5 are different, although they are in practice the same necklace. We needto account for this symmetry somehow. All necklaces that are equivalent under somerotation should be subtracted from the count of total possible arrangements.Chapter 3 Malcolm Hillebrand 45haotic Dynamics of Polyatomic ModelsF IGURE
Two apparently distinct necklaces. However, it is clear that theyare simply rotated versions of the same necklace. Harking back to theoriginal problem, periodic DNA molecules are clearly unaffected by a ro-tation like this. So for our purposes we would like to consider these twonecklaces to be equivalent.
Furthermore, for the PBD model of DNA, there is a re fl ectional symmetry. Due tothe fact that it does not matter which way we traverse the sequence in this model, thetwo necklaces in Fig. 3.6 are functionally equivalent. This fi nally completes our set ofF IGURE
Two arrangements which are again in principle different, but inactual fact are considered equivalent as they are simply re fl ections of eachother. restrictions on the necklace. The problem now is to count all possible colourings, orpermutations of beads in the necklace, and then to discount all duplicate colourings. The basic lemma underlying the theory of counting distinct colourings is often referredto as Burnside’s Lemma [ ] . But in order to make use of this, we need fi rst to translatethe colouring problem to a problem using sets and groups. So we shall consider a set B , containing all the beads we wish to colour. The cardinality of this set is of course N . Then we wish to colour each of these beads with an element of the set of possiblecolours Ω . In our case this set Ω contains just two colours, black and white, but all thatfollows applies to more general sets of colours. Our colouring is then a mapping from B to Ω , giving each bead in B a colour from Ω . We will denote this map as ϕ : B → Ω .With this description of the colourings, we can label the set of all possible colourings ϕ from B to Ω as Ω B . Explicitly, Ω B : = { ϕ : B → Ω } Next, we will use a group action to account for the symmetries in the necklaces. Ifwe have some symmetry group G with a group action A G , then this group action on theset B describes the symmetry we need to account for. Intuitively, the chosen group G will be associated with certain symmetries. The group action A G on the set B will then,in the course of the algorithm, allow us to eliminate these repeats. What is important to46 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Modelsnote here is that if we have this action A G on B , which maps to permutations of B , thenthese permutations can apply in exactly the same way to the colours associated with thepermuted elements. So if some group element g takes an element b ∈ B to c ∈ B , thenthe same element g would map ϕ ( b ) ∈ Ω to ϕ ( c ) ∈ Ω B . Concisely, the group action A G on B gives rise to an associated group action ˜ A G on Ω B such that g ϕ : b (cid:13)→ ϕ ( g b ) .Now we have the descriptions we need – we have a set of elements, and a set of possi-ble colourings, as well as a group which will account for the relevant symmetries. Notethat the number of distinct possible colourings is equal to the number of distinct orbitsof the group action ˜ A G . Intuitively, the key point is that all colourings that are equiv-alent under symmetry are included in a single orbit. This is the fi rst step to reducingthe number of excess colourings. The second step is of course to eliminate orbits thatthemselves are duplicates. So fi nally counting only these distinct orbits, we have actuallycounted all the truly independent colourings of the necklace.Burnside’s lemma gives us a very ef fi cient way to calculate the number of distinctorbits of a group action on a set [ ] : Lemma 1 (Burnside’s Lemma) . In the group action of a fi nite group G on a set S, thenumber of distinct orbits is given by the average number of fi xed points of the elements in G.N DO = | G | (cid:1) g ∈ G (cid:15)(cid:15)(cid:15) Fixed g (cid:15)(cid:15)(cid:15) (3.4)On the face of it, this may not appear particularly helpful, but in many cases it isactually deceptively simple. To illustrate the process more clearly, consider a simpleexample. Example 2.
Compute the number of cyclically distinct arrangements of black and whitebeads in a fi ve-bead necklace, i.e. colourings that are equivalent up to rotations of the neck-lace. In this simple fi ve bead case, we can actually explicitly construct all the colouringsthat are not rotations of another. All the possible arrangements, without any restric-tions, are shown in Fig. 3.7.Figure 3.8 shows the arrangements of necklaces that are not rotationally equivalent.We can see that there are eight such distinct colourings, and one can convince oneselfof this by comparing every possible colouring in Fig. 3.7 and seeing that the remainingnecklaces are in fact not unique under rotation.Now, let us see what Burnside’s lemma gives us. The formula requires the cardinalityof the symmetry group, and the number of necklaces fi xed or left invariant by each element of the group. Since the symmetry in question is rotation, we will use the cyclicgroup of order fi ve, C . There are of course fi ve elements in this group, rotations of onethrough four places as well as the identity permutation which leaves each bead as it is.So now we have N DO = (cid:1) g ∈ C | Fixed g | . (3.5)For each element in C , i.e. each rotation, we need to count how many necklaces (out ofall the possible necklaces) remain unchanged by a rotation. Since there are fi ve beads oftwo colours, the total number of possible necklaces is 2 =
32. If we now consider eachrotation case by case:Chapter 3 Malcolm Hillebrand 47haotic Dynamics of Polyatomic ModelsF
IGURE
All possible arrangements of the fi ve bead necklace. It is imme-diately clear for instance that all necklaces with four black beads and onewhite bead are in fact rotationally equivalent. • Rotation by 0 places (Identity). As this leaves every necklace invariant, | Fixed R | = • Rotation by 1 place. Now the only necklaces left invariant are the pure blackand pure white necklaces. Rotating every other necklace gives a slightly differentnecklace each time. So then | Fixed R | = • Rotation by 2 places. Once again, because of the odd number of beads in thenecklace, there are only the same two invariant necklaces, giving | Fixed R | = • Rotation by 3 places. As above. So | Fixed R | = • Rotation by 4 places. Again, as above. | Fixed R | = IGURE
Rotationally distinct arrangements of black and white beads ina fi ve bead necklace. Clearly none of these necklaces can be constructedmerely by rotating another. However, if we consider any other arrange-ment of beads, we will fi nd that it is simply a rotation of one of these eightfundamental necklaces.
48 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic ModelsNow, we can compute the number of distinct orbits using (3.5). So we get N DO = ( + + + + ) = = fi ned in Section 3.2 andrestate Burnside’s lemma as N DO = Z G ( k , k , ..., k ) , (3.6)where the symmetry is given by a group G , and there are n beads coloured with k colours. Note that the cycle index takes n arguments, so in (3.6) k enters n times intothe bracket.If we now apply this process to our example, we have a much easier route to thenumber of distinct colourings. Recalling from earlier (3.2) the cycle index of the cyclicgroup, we can calculate the cycle index of the order fi ve group as Z C ( x , x , x , x , x ) = x + x , (3.7)since the only integers in {
1, 2, 3, 4, 5 } that divide 5 are 1 and 5, and ϕ t ( ) = fi nd the number of distinct colourings of a fi ve bead necklace colouredwith two colours from Burnside’s lemma is very straightforward: N DO = Z C (
2, 2, 2, 2, 2 ) = · + · = + =
8, (3.8)just as we had before.Among the obvious advantages of this approach in general is that we can triviallycompute the possibilities for different numbers of colours by changing k . While we doneed to compute the cycle index for each value of n we are interested in, since we havea general formula for most symmetry groups this is not a major hardship, and certainlyfar less work than manually fi nding the fi xed points of every permutation in G .So we now have a robust algorithm to calculate the number of distinct necklacesof n beads with k colours, with various symmetries. This is great progress. If all wewanted was to fi nd the number of possible DNA sequences with given length, this wouldsuf fi ce. However it addresses neither of our chief concerns: Nowhere does this methodaccount for speci fi c numbers of black and white beads, nor does it consider the numberof alternations of these beads. Some more work is required in order to include theserestrictions.Chapter 3 Malcolm Hillebrand 49haotic Dynamics of Polyatomic Models The fi rst topic we will address is that of specifying set numbers of black and white beads.In order to account for numbers of different kinds of beads, we need a mathematicalway of identifying the bead types. The way we do this is by introducing the concept ofa weight for each colour. So for each colour c in our set of colours Ω , we associate withit a weight w c ∈ (cid:1) . The function that assigns these weights is simply w : Ω → (cid:1) , where w ( c ) = w c .The weight of a particular colouring ϕ : B → Ω , we can fi nd the total weight of thiscolouring by taking the sum of the weights of all the colours: (cid:15)(cid:15)(cid:15) W ϕ (cid:15)(cid:15)(cid:15) = (cid:1) ≤ i ≤ | B | w ( ϕ ( b i )) . (3.9)These weights then in some way take note of the colours of each bead in a given colour-ing.What is signi fi cant about the total weight is that it allows us to identify colour-ings with certain numbers of colours. More concretely, take the pertinent case of twocolours, black and white, for a necklace of N beads in a set B . Assign w white = w black =
2, with N white + N black = N in some particular colouring k . Then the totalweight will be given by | W k | = (cid:1) i = N w ( k ( b i )) , b i ∈ B = N white + N black ,and with the restriction that N white + N black = N , there is of course only one total weightthat possible for a speci fi c choice of N white and N black . This means that we can rephrasethe problem of counting arrangements with N white white beads and N black black beadsas a problem of counting the number of distinct colourings with a given weight. Forinstance, in the above case, if we were to set N =
10 and N white = N black =
7, then weare looking for colourings of 10 beads satisfying | W k | = · + · = fi nd some way of counting these distinct colourings.What has not appeared in the above discussion of course is this notion of distinctcolourings. We still want to account for equivalence classes of colourings, and onlycount non-equivalent colourings. Of course, counting non-equivalent colourings is thesame as counting the equivalence classes of the colourings, up to some desired symmetry.So now we are back in the territory of Burnside’s lemma, where we wish to account forthese symmetries with some group action. If we have some symmetry group G actingon our set of beads B , then what is important to note is that the orbits of A G on B allhave the same total weight. This is the crux of the idea – because in a given orbit, wehave the same set of coloured beads being shuf fl ed around (rotated, re fl ected, etc.), thetotal weight of course remains the same. The same beads, however permuted, have thesame total weight.So in order to count the number of distinct colourings, we need to count the numberof distinct orbits under the desired group action. Correspondingly, to count distinct50 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Modelscolourings with given numbers of beads, we need to count the distinct orbits with therequisite total weight.In order to generalise Burnside’s lemma to this case, we will use Pólya’s powerfultheorem for counting distinct colourings with chosen weights [ ] . Before moving tothe actual theorem however, we do need to introduce one more piece of mathematics: Generating functions . These are essentially simple objects. The idea of a generating func-tion is to encode information about a certain problem or process in the coef fi cients ofa power series [ ] . Frequently they are used in counting problems, where these coef fi -cients can represent solutions for different cases. This is the application used here – wewill see that the computation reduces to fi nding the desired coef fi cient in the generatingfunction series. An ordinary generating function takes the form [ ] f ( x ) = (cid:1) n = c n x n , (3.10)which is a power series in x with coef fi cients c n .We are however interested in a more speci fi c class of generating functions. What wewant to investigate are colourings with a particular total weight. So at this point, weneed a generating function for a given assignment of weights. Take an assignment ofweights w , which assigns each colour in Ω a natural number in (cid:1) . If there are a fi nitenumber of colours with each weight, i.e. the inverse function w − ( n ) maps to a fi nitenumber of colours, then we can de fi ne a generating function for this weight assignmentas the following: f w ( x ) = (cid:1) n = | w − ( n ) | x n , (3.11)where | w − ( n ) | is the number of colours with a particular weight n .Now we are ready to understand Pólya’s enumeration theorem. While there arefurther mathematical tricks that will be introduced in due course to facilitate rapid cal-culations, we can present the framework of the theorem. Theorem 2 (Pólya’s Enumeration Theorem) . Take a set of objects B , a set of colours Ω , afunction w : Ω → (cid:1) assigning weights to each colour with an associated generating functionf w , with a group G acting on the set B through the action A G . This action passes to a groupaction ˜ A G on the set Ω B , and we have that the generating function for the number of distinctorbits of ˜ A G by total weight isf ˜ A G DO ( x ) = Z G ( f w ( x ) , f w ( x ) , . . . , f w ( x n )) , (3.12) where Z G is the cycle index of the group G. This theorem allows to obtain information about the number of colourings witha desired total weight from the cycle index of the symmetry group. As a generatingfunction, the coef fi cients of this series give the number of colourings with a certain totalweight . So the number of colourings with a total weight of k is given by the coef fi cientof the k th term in the polynomial.To illustrate this, let us return to the earlier example (Example 2), where we com-puted the total number of possible colourings of cyclically distinct necklaces with twocolours. Now let us restrict further, and say we want only necklaces with three whiteChapter 3 Malcolm Hillebrand 51haotic Dynamics of Polyatomic Modelsbeads and two black beads. Note that this restriction was not available to us in Burnside’slemma, so we do in fact need to use Pólya’s enumeration theorem.We will once again assign weights very simply: w ( white ) = w ( black ) =
2. So wehave that the total weight W = N white + N black , and N = N white + N black . Recalling ourreformulated aim, which is to count colourings with a given total weight, we need toidentify the desired total weight in this case. Since we want N white = N black = W = · + · =
7. The next step is then to produce the generating functionpolynomial, and fi nd the coef fi cient of the term in x , as this coef fi cient will tell us thenumber of colouring with total weight of seven.The ingredients we need are 1) the relevant cycle index, and 2) the generating functionfor these weights. We of course already have the cycle index (3.7), so let us look at thegenerating function. The generating function will be as given in (3.11). In our case, wehave two colours, so the only weights that have colours assigned to them are the numbers1 and 2. So expressing our generating function, we will have f w ( x ) = ∞ (cid:1) n = | w − ( n ) | x n = | w − ( ) | x + | w − ( ) | x . (3.13)Now we can use the fact that each weight corresponds to only one colour – the weight 1corresponds to the colour white, and the weight 2 to the colour black. This means that | w − ( ) | = | w − ( ) = | , and our generating function becomes simply f w ( x ) = x + x . (3.14)So we now have all the pieces to construct the fi nal polynomial. Pólya’s enumerationtheorem tells us that f ˜ A G DO ( x ) = Z G ( f w ( x ) , f w ( x ) , . . . , f w ( x n )) f ˜ A C DO ( x ) = Z C ( f w ( x ) , f w ( x ) , f w ( x ) , f w ( x ) , f w ( x ))= Z C ( x + x , x + x , x + x , x + x , x + x )= ( x + x ) + ( x + x )= x + x + x + x + x + x . (3.15)This is now our generating function polynomial for distinct colourings under rotationalsymmetry, by total weight. So we can read off the coef fi cients and recover the number ofdistinct colourings with any distribution of black and white beads. Our original aim was N white = N black =
2, with total weight of seven. Clearly there are two possible distinctcolourings in this case, as the coef fi cient of the x term is equal to two. Pleasingly, wealso get all the information about other arrangements at the same time – so we can fi ndthat the number of arrangements with two white and three black beads is also two forinstance. We can in fact gain all the information we gained by painstakingly constructingFig. 3.8, just by reading each coef fi cient.So this illustrates the power of the method. We can, in general, use this to fi nd thenumber of distinct colourings with particular numbers of black and white beads. Butof course, this is not the end – we have a further restriction, as well as a dif fi culty onthe horizon. The further restriction is that we want to enforce a certain number of52 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Modelsalternations between differently coloured beads. Thus we need to move from the no-tion of coloured beads to the containers of beads introduced in Section 3.1.1, wherethe number of alternations maps directly to the number of containers. The challengewhich quietly raised its head in the course of the example above is the expansion of thepolynomial. While in this case we had a very small necklace, and could expand the poly-nomial without undue dif fi culty, for large polynomials this becomes computationallyvery expensive.We will address this, and some other computational details, in due course. But fornow, we need to turn to this idea of the containers, and how we can use Pólya’s enumer-ation theorem to help us fi nd what we want. So let us recall the container construction described in Fig. 3.4. The obvious advantageof this description is that if we are working with these containers, then the problemsabout the alternation index are immediately solved – if we have 2 M alternations, thenwe simply have 2 M containers, M containers of each colour. So far so good, there isin principle no problem with using these containers as the fundamental elements ofour necklace. But what else do we require in order to use Pólya’s counting theorem?The symmetry group is fi ne, we can still choose the same group to account for relevantsymmetries. So if we have this, and of course the corresponding cycle index, all thatremains is to fi gure out the weights. Conveniently, there is also a very natural way totreat the weights – instead of weighting white and black beads with 1 and 2 respectively,we will now colour each container with the number of beads it contains . So if a containerholds seven beads, its weight is seven. Notice a critical point from this assignment ofweights: the total weight (the crucial element we are looking at) of the black containerswill give the total number of black beads, and the same for the white containers.What is not clear in this situation is how to deal with the different coloured con-tainers; how does the symmetry group act on the containers? And how does Pólya’stheorem apply if we have to treat the different colours separately? To this end, we needto consider the special case of bipartite sets and group actions.So let us take a set B , which is separated into two partitions, so B = Y ⊔ Z , where ⊔ denotes a disjoint union. In our case, we can see that our set will be the set of containers,and the partitions the sets of black and white containers. What we then need to consideris a colouring that colours elements of the set Y from a set Ω Y , and colours elements of Z from Ω Z . Thus a valid colouring ϕ according to this bipartite viewpoint is any colouring ϕ : B → Ω Y ⊔ Ω Z if it takes only elements of Y to Ω Y and elements of Z to Ω Z , with nocrossover. So concisely, ϕ ( a ) ∈ Ω Y ⇐⇒ a ∈ Y and ϕ ( a ) ∈ Ω Z ⇐⇒ z ∈ Z .Now that we have made sure that our colouring preserves the bipartite nature of theset, all that remains is to discuss the group action, which should do something similar.Our only additional requirement on the group action, in order for us to be able to statea version of Pólya’s theorem for bipartite sets, is that the group action needs to maintainthe separation of the sets. So we will call this kind of group action partition-preserving . De fi nition 1 (Partition-Preserving Group Action) . Given a group G, and a partitionedset B = Y ⊔ Z , we say that the group action A G is partition preserving if for each g ∈ G,g a ∈ Y ⇐⇒ a ∈ Y and g ∈ G , g a ∈ Z ⇐⇒ a ∈ Z .
Chapter 3 Malcolm Hillebrand 53haotic Dynamics of Polyatomic ModelsWith this requirement, the group action in some way splits into two effectively iden-tical actions acting simultaneously on the two partitions. In particular, it means thatthe cycles of the cycle decomposition of elements of the group action are themselvespartition-preserving. So if we write the decomposition of A G ( g ) , for some element g ∈ G , we will have c · c · c · · · · c m , and each c i permutes elements of the same setpartition. Not all cycles have to be in the same partition – c j and c k can permute ele-ments in partition Y and Z respectively, but no cycle can itself move elements betweenthe partitions.We will then denote by c Ym ( g ) the number of m − cycles in the cycle decomposition of A G ( g ) , that are con fi ned to the partition Y . Similarly, c Zm gives the number of m − cyclesin partition Z . From this we have a fairly complete bookkeeping of the group actionon the partitioned set. Particularly, we are able now to de fi ne the bipartite cycle index ,which will allow us to move on to the fi nal statement of a Pólya counting theorem forthese bipartite sets. De fi nition 2 (Bipartite Cycle Index) . Given a set B = Y ⊔ Z of n elements, and a groupG, the bipartite cycle index of a partition-preserving group action A G on B is de fi ned as ˜ Z G ( y , . . . , y n , z , . . . , z n ) = | G | (cid:1) g ∈ G y c Y ( g ) . . . y c Yn ( g ) n z c Z ( g ) . . . z c Zn ( g ) n (3.16)Now, we have all the tools to state our bipartite counting theorem [ ] . Theorem 3 (Bipartite Pólya Enumeration Theorem) . Take a fi nite group G, with a partition-preserving group action A G on a fi nite set B = Y ⊔ Z . Let the set of colours be Ω = Ω Y ⊔ Ω Z ,with weights w Y : Ω Y → (cid:1) + and w Z : Ω Z → (cid:1) + and corresponding generating functionsgiven by f Y and f Z . If the set of valid colourings of B is ϕ , then the group action A G passesnaturally to a group action ˜ A G on ϕ , and the generating function by total weight for thenumber of distinct orbits of ˜ A G isf ˜ A G DO ( a ) = ˜ Z G ( f Y ( a ) , . . . , f Y ( a m ) , f Z ( a ) , . . . , f Z ( a m )) . (3.17)So now we are able to use this theorem to solve our actual problem: We can identifyterms with the desired total weight in each variable, which corresponds to the numberof black and white beads, and the number of alternations is given by the total numberof elements (so the cardinality of the set B will be equal to α ). All that remains is toconstruct the bipartite cycle index for the dihedral group, since that is the symmetrywe are concerned with, and then to fi nd ef fi cient ways of computing the coef fi cients ofrelevant terms of the polynomial. Now, in order to solve our original problem, we need to adapt the relevant symmetrygroup to the bipartite structure of the necklace of containers. The symmetry groupunder consideration is still the dihedral group, as this accounts for the rotational andre fl ective symmetries in our necklace. It is worth noting here that the PBD model of54 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic ModelsF IGURE
An illustration of a necklace of six containers, showing the actionof a rotation on the containers. The black and white containers are rotatedindependently.
DNA does not have a prescribed direction of traversal – whether one moves clockwiseor counter-clockwise along the sequence, the dynamics are identical. This means thatre fl ectionally equivalent sequences are dynamically equivalent. In DNA models with aset direction of traversal (such as the ePBD model) there is no longer this re fl ective sym-metry, as a re fl ection would mean essentially moving along the sequence in the oppositedirection. In these models the symmetry group simpli fi es to the cyclic group. However,in this study we are sticking to the PBD model, and so we use the dihedral group.In order to construct the bipartite cycle index of the dihedral group D M , we look atthe two components separately: The contribution from the cyclic component, and thecontribution of re fl ective symmetry. The action of rotations on the containers is verystraightforward. Figure 3.9 shows the partition-preserving rotations, simply movingeach container around, while maintaining the separation of black and white beads.So the group generated by these rotations, the cyclic group, has the usual cycle indexform, but adapted to the bipartite nature by the inclusion of two sets of variables. Recall(3.2), for the cycle index of the cyclic group in one set of variables. We adapt it to thebipartite case by simply adding in the second set of variables in a natural way:˜ Z C M ( y , y , ..., y M , z , z , ..., z M ) = M (cid:1) p | M ϕ t ( p ) y M / pp z M / pp . (3.18)This accounts for half of the elements of D M , as the dihedral group consists of halfrotations and half re fl ections.Accounting for the re fl ections requires a little more thought however. Consider thetwo cases given in Figs. 3.10 and 3.11, corresponding to re fl ections with M odd or even.Recall that there are a total of 2 M containers in the necklace. We can see that in the oddcase, the re fl ections are fairly straightforward affairs – one white container stays where itis, one black container stays where it is, and every other container switches places withthe container opposite (forming 2-cycles). Thus the total cycle index for the case of M being odd is˜ Z D M ( y , ..., y M , z , ..., z M ) = M (cid:1) p | M ϕ t ( p ) y M / pp z M / pp + y z y ( M − ) / z ( M − ) / . (3.19)Now we need to consider the slightly more complex case of M being even. For this,we see that in half the re fl ections we fi x two white containers and in half the re fl ectionswe fi x two black containers. In all these cases the rest of the containers split into 2-cycles. So to account for this we will have two terms in the cycle index – one to accountChapter 3 Malcolm Hillebrand 55haotic Dynamics of Polyatomic ModelsF IGURE
A re fl ection with M odd. Note that this fi xes one black and onewhite container. F IGURE
A re fl ection with M even. This now fi xes one pair of eitherblack or white containers. for the half fi xing two white containers, and one to account for the half fi xing two blackcontainers. The fi nal form for even M is then˜ Z D M ( y , ..., y M , z , ..., z M ) = M (cid:1) p | M ϕ t ( p ) y M / pp z M / pp + (cid:13) y y ( M − ) / z M / + z z ( M − ) / y M / (cid:14) .(3.20)Now armed with our bipartite cycle index, we are in principle set to create polyno-mials, and fi nd numbers of possible arrangements from the coef fi cients. Particularly, wealso have as a corollary from our bipartite Pólya counting theorem the bivariate gener-ating function [ ] f ˜ A G DO ( a , b ) = ˜ Z G ( f Y ( a ) , . . . , f Y ( a m ) , f Z ( b ) , . . . , f Z ( b m )) , (3.21)which gives a generating function for the number of distinct colourings of the set B bytotal weight in Ω Y and Ω Z .We still have not addressed the issue of actually computing the coef fi cients of theterms in this polynomial, but let us fi rst look at an example to illustrate why we veryquickly reach a point where brute force techniques do not work. Example 3.
Find the number of distinct necklaces with seven white beads, six black beads,and ten total alternations.
So, we have α = M = N W = N B =
6. Thus we have the case where M isodd, and we can immediately identify our generating function using the bipartite Pólya’s56 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Modelstheorem with the cycle index given by (3.19). f ˜ A G DO ( a , b ) = ˜ Z G ( f Y ( a ) , . . . , f Y ( a m ) , f Z ( b ) , . . . , f Z ( b m ))= ˜ Z D (cid:6) f W ( a ) , . . . , f W (cid:6) a (cid:7) , f B ( b ) , . . . , f B (cid:6) b (cid:7)(cid:7) = (cid:1) p | ϕ t ( p ) f / pW ( a p ) f / pB ( b p ) + f W ( a ) f B ( b ) f W (cid:6) a (cid:7) f B (cid:6) b (cid:7) = (cid:16) ( a + a + . . . ) ( b + b + . . . ) + ( a + a + . . . )( b + b + . . . ) (cid:17) + ( a + a + . . . )( a + a + . . . ) ( b + b + . . . )( b + b + . . . ) . (3.22)Here we the generating functions f W and f B are used as introduced earlier, simply powerseries in their argument. Note here that while in previous examples we truncated thepower series to only two terms, since in principle we can have any number of colourshere (recall that the colours are the number of beads in each container, which for longnecklaces could become arbitrarily large), we will work with the full power series. Thisleaves us with some dif fi culties. The fi rst is that we have in fi nite series to multiply –this is obviously impossible. Since we are only interested in the coef fi cient of terms in a b (recall this is the total weight required for the given numbers of beads), we cantruncate the series though, so that we only expand terms that can contribute this term.In this case, we can truncate the fi rst polynomials (which are raised to the power of fi ve)at fourth order since higher powers have no effect on the end polynomial below termsof ninth order. We are not actually concerned with the second pair of polynomials –since the powers only appear in multiples of fi ve, none of these will have an impact onour targetted term. The fi nal four polynomials can also be truncated. Since the highestpower concerned is eight, the two squared polynomials can be truncated to two terms,since all other terms will ultimately enter in order nine or above. Similarly, the non-squared polynomials can be truncated to four terms, as after this the added terms do nothave an impact. So for our purpose, we just need to expand the following polynomial,and extract the coef fi cient of a b :110 (cid:16) ( a + a + a + a ) ( b + b + b + b ) (cid:17) + ( a + a + a + a )( a + a ) ( b + b + b + b )( b + b ) . (3.23)Needless to say, this is best performed using a computer algebra system. Expanding thetwo multiplications using Sage [ ] we fi nd the relevant parts of the expansion to be110 (cid:6) . . . 75 a b . . . (cid:7) + (cid:6) . . . 3 a b . . . (cid:7) . (3.24)But we now have the fi nal result – the coef fi cient of a b , the number of necklaces withseven white and six black beads up to dihedral symmetry, is 75 / + / = f w ( x ) n grows exponentially with n , so it is imperative that we fi nd a better method.To save us this trouble, we will introduce some more mathematical short cuts. The fi rst of these is to consider our generating functions as formal power series. The theoryof formal power series (see e.g. [ ] ) gives a large number of properties that we canpotentially take advantage of, but in particular we will use the fact that there is a formof the binomial theorem that holds. Lemma 2. If ( − x ) − n denotes the formal inverse of ( − x ) n , then we have that ( − x ) − n = ∞ (cid:1) i = (cid:2) n + i − n − (cid:3) x i , (3.25)Just as the binomial theorem allows us to calculate binomial coef fi cients, so we areable to use this theorem, with some massaging, to fi nd the coef fi cients in the powerseries expansions. The key here is that the computation of these binomial coef fi cients isperformed in linear time (computation time increasing as a linear function of n ), so thisis a major boost in performance from the naive expansion of enormous polynomials.Particularly, we have another lemma we can state about the powers of the generatingfunctions f w ( x ) . Lemma 3.
The generating function f w ( x ) n can be written as a formal power series asf w ( x ) n = ∞ (cid:1) i = (cid:2) n + i − n − (cid:3) x n + i , (3.26)This lemma is straightforward to prove from Lemma 2. Since x f w ( x ) = x + x + . . . = f w ( x ) − x , we can see that clearly f w ( x ) = x ( − x ) − , which gives us that f w ( x ) n = x n ( − x ) − n . Thus, we now have the result from Lemma 2.Putting these results into the fi nal form, for fi nding relevant coef fi cients, we havethat the coef fi cient of x p in the polynomial f w ( x a ) b is given by (cid:16) f w ( x c ) d (cid:17) p =
1, if d = c = d = c > d > c ∤ p or p < c d , (cid:6) p / c − d − (cid:7) , otherwise. (3.27)and that the coef fi cient of x p in the polynomial f w ( x a ) b · f w ( x a ) b is given by (cid:16) f w ( x c ) d · f w ( x c ) d (cid:17) p = p (cid:1) i = (cid:16) f w ( x c ) d (cid:17) i (cid:16) f w ( x c ) d (cid:17) p − i (3.28)where the coef fi cients on the right hand side are computed according to (3.27).Now, despite the daunting appearance of these formulas, we have a computationallyef fi cient and fundamentally straightforward way of computing the coef fi cients of thepolynomial terms. The complexity of the general formulas diminishes rapidly whenconfronted with concrete cases, so let us now illustrate their use by redoing example 3.58 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic ModelsThe cycle index is given in (3.22). Once again, let us look at this one polynomial set ata time. Recalling that we want terms in a b , for the fi rst pair, we have ( a + a + . . . ) ( b + b + . . . ) .So the based on the procedures above, we want to fi nd [ f ( a ) ] , the coef fi cient of theterm in a , for the generating function f ( a ) raised to the fi fth power. This means thatin (3.27), we have c = d = p =
7, and hence (cid:16) f ( a ) (cid:17) = (cid:2) − − (cid:3) = (cid:2) (cid:3) = b , we fi nd (cid:16) f ( b ) (cid:17) = (cid:2) − − (cid:3) = (cid:2) (cid:3) = fi cient of a b in the fi rst polynomial pair is 15 · =
75. As before,the second polynomial pair contributes nothing, as according to (3.27) all the coef fi cientsare 0.The fi nal polynomial, from the re fl ective component, requires us to make use of(3.28), as we have a product of generating functions with different powers. Considering fi rst the series in b . What we want is (cid:16) f ( b ) · f ( b ) (cid:17) = (cid:1) i = [ f ( b )] i (cid:16) f ( b ) (cid:17) − i ,and after checking each term in the sum, we realise that the only nonzero term comesfrom the case where i =
2, and this gives us (cid:16) f ( b ) · f ( b ) (cid:17) = (cid:2) − − (cid:3)(cid:2) / − − (cid:3) = (cid:2) (cid:3)(cid:2) (cid:3) = a .Here we now have (cid:16) f ( a ) · f ( a ) (cid:17) = (cid:1) i = [ f ( a )] i (cid:16) f ( a ) (cid:17) − i ,and we need to work through each term in the series to fi nd the coef fi cient. Here wehave nonzero contributions from the terms where i = i =
3, corresponding toterms from a · a and a · a . So the overall contribution is (cid:16) f ( a ) · f ( a ) (cid:17) = (cid:2) (cid:3)(cid:2) (cid:3) + (cid:2) (cid:3)(cid:2) (cid:3) =
2. (3.30)Now adding together all the contributions, we fi nd that in the end we have110 ( ) + ( + ) = fi cient way of calculating the number of possible distinct neck-laces with a given number of white beads, a given number of black beads, and a setnumber of alternations from one colour to another. In order to produce the distribu-tion then, all one needs to do is to select the numbers of black and white beads, and thenstep through every allowed value of α for those beads. Then normalising the numberof possible necklaces for each α by dividing through by the total number of necklacesacross all values of α .This was a fairly indirect route to fi nding what looked like being a relatively simpleprobability distribution. Especially since, as we will see in the next subsection, MonteCarlo simulations suggested that the distribution should be more or less normal. Butthis approach, for all its apparent abstract complexity, does illustrate the beauty andpower of applying branches of pure mathematics to very practical, applied problems.Ultimately, we have successfully created an ef fi cient probability distribution for DNAsequences, as desired. Here we present some plots displaying properties of this distribution, as numbers of ATand GC base pairs are varied, as the ratios are varied, and so forth.We begin by comparing the theoretical results to Monte Carlo simulations. Thishas a twofold purpose – fi rstly it serves as a check for the theory, that we have not gotlost somewhere in the mathematics. Secondly, it gives an idea of how many randomrealisations are necessary in order to more or less cover the spectrum of possible cases.Figure 3.12 shows the comparison between the theoretical calculation and Monte Carlosimulations. We see that the results are very close, which is a useful validation of thetheory. It is also worth noting that since the possible number of DNA sequences with N AT = N GC =
50 and α =
50 is on the order of 10 , which makes it unsurprising thatthe Monte Carlo approximation is imperfect.However, in order to check the convergence of the Monte Carlo method and un-derstand how the number of simulations affects the accuracy, we can compute the totalabsolute difference between the true distribution and the approximate distribution. Wede fi ne this difference d as d ( N MC ) = (cid:1) α | P MC ( N MC , α ) − P ( α ) | . (3.31)Here P MC ( N MC , α ) is the probability of fi nding a sequence with α alternations, usingMonte Carlo simulations with N MC realisations. We see in Fig. 3.12(c) that as N MC isincreased, we see the expected improvement in accuracy. There is however a signi fi cantslowing down in the improvement around N MC = α values fairly well.The other point that arises immediately after looking at Fig. 3.12(a) and (b) is that theshape of the distribution is more or less normal. This suggests that we should be able to fi t the distribution with a Gaussian function, which would give us some useful insightsinto the properties of the distribution as the parameters ( N AT , N GC ) change. Figure 3.1360 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Models (cid:0) (cid:1)(cid:0) (cid:2)(cid:0) (cid:3)(cid:0) (cid:4)(cid:0) (cid:5)(cid:0)(cid:0) (cid:0) (cid:0)(cid:6)(cid:0)(cid:0)(cid:0)(cid:6)(cid:0)(cid:1)(cid:0)(cid:6)(cid:0)(cid:2)(cid:0)(cid:6)(cid:0)(cid:3)(cid:0)(cid:6)(cid:0)(cid:4)(cid:0)(cid:6)(cid:5)(cid:0)(cid:0)(cid:6)(cid:5)(cid:1)(cid:0)(cid:6)(cid:5)(cid:2)(cid:0)(cid:6)(cid:5)(cid:3)(cid:0)(cid:6)(cid:5)(cid:4) (cid:0) (cid:0) (cid:0) (cid:1) (cid:0)(cid:1)(cid:2) (cid:0) (cid:0)(cid:1) (cid:0) (cid:1)(cid:2)(cid:0) (cid:2)(cid:3) (cid:0) (cid:3)(cid:2) (cid:0)(cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:1)(cid:10)(cid:11)(cid:4)(cid:1)(cid:8)(cid:4)(cid:3)(cid:12)(cid:13)(cid:7)(cid:9) (cid:0) (cid:1)(cid:0) (cid:2)(cid:0) (cid:3)(cid:0) (cid:4)(cid:0) (cid:5)(cid:0)(cid:0) (cid:0) (cid:0)(cid:6)(cid:0)(cid:0)(cid:0)(cid:6)(cid:0)(cid:1)(cid:0)(cid:6)(cid:0)(cid:2)(cid:0)(cid:6)(cid:0)(cid:3)(cid:0)(cid:6)(cid:0)(cid:4)(cid:0)(cid:6)(cid:5)(cid:0)(cid:0)(cid:6)(cid:5)(cid:1)(cid:0)(cid:6)(cid:5)(cid:2)(cid:0)(cid:6)(cid:5)(cid:3)(cid:0)(cid:6)(cid:5)(cid:4) (cid:0) (cid:0) (cid:0) (cid:1) (cid:0)(cid:1)(cid:2) (cid:0) (cid:0)(cid:1) (cid:0) (cid:1)(cid:2)(cid:0) (cid:2)(cid:3) (cid:0) (cid:1)(cid:2) (cid:0)(cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:1)(cid:10)(cid:11)(cid:4)(cid:1)(cid:8)(cid:4)(cid:3)(cid:12)(cid:13)(cid:7)(cid:9) (cid:0)(cid:1)(cid:1)(cid:1)(cid:1) (cid:2)(cid:1)(cid:1)(cid:1)(cid:1) (cid:3)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0) (cid:0)(cid:1) (cid:1)(cid:4)(cid:1)(cid:1)(cid:1)(cid:4)(cid:1)(cid:2)(cid:1)(cid:4)(cid:1)(cid:5)(cid:1)(cid:4)(cid:1)(cid:6)(cid:1)(cid:4)(cid:1)(cid:7)(cid:1)(cid:4)(cid:0)(cid:1)(cid:1)(cid:4)(cid:0)(cid:2) (cid:0) (cid:0) (cid:1) (cid:0)(cid:1)(cid:2) F IGURE
Comparison of the theoretically determined probability dis-tribution (empty circles) and Monte Carlo generated results using 20 000realisations ( fi lled stars) for two cases. (a) N AT = N GC =
40 and (b) N AT = N GC =
50. We see that the two methods agree fairly closely, as de-sired. The Monte Carlo results will never agree exactly, due to fi nite statis-tics, but these results do validate the theoretical method and show thatwith a fi nite number of Monte Carlo simulations a good approximationto the true distribution can be made. Panel (c) Shows the total absolute dif-ference d for increasing numbers of statistics. The points are averaged over fi ve Monte Carlo runs with each number of realisations, and the errorbarshows the standard deviation from this averaging. Chapter 3 Malcolm Hillebrand 61haotic Dynamics of Polyatomic Models (cid:0) (cid:1)(cid:0) (cid:2)(cid:0) (cid:3)(cid:0) (cid:4)(cid:0) (cid:5)(cid:0)(cid:0) (cid:0) (cid:0)(cid:6)(cid:0)(cid:0)(cid:0)(cid:6)(cid:0)(cid:1)(cid:0)(cid:6)(cid:0)(cid:2)(cid:0)(cid:6)(cid:0)(cid:3)(cid:0)(cid:6)(cid:0)(cid:4)(cid:0)(cid:6)(cid:5)(cid:0)(cid:0)(cid:6)(cid:5)(cid:1)(cid:0)(cid:6)(cid:5)(cid:2)(cid:0)(cid:6)(cid:5)(cid:3)(cid:0)(cid:6)(cid:5)(cid:4) (cid:0) (cid:0) (cid:0) (cid:1) (cid:0)(cid:1)(cid:2)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:9)(cid:1)(cid:7)(cid:10)(cid:11)(cid:12)(cid:3)(cid:13)(cid:14)(cid:3)(cid:2)(cid:1)(cid:15)(cid:7)(cid:16) F IGURE A fi tting of the probability distribution in the N AT = N GC = fi t wellindeed. For this case, we fi nd that the mean value is at α = σ = shows the fi tting of the data with a Gaussian function of the form P ( α ) = σ (cid:8) π e − ( α − α σ ) . (3.32)We can explore some features of the distribution by fi xing the number of AT basepairs, and then varying the number of GC base pairs. Of course the same results wouldbe found if the number of GC base pairs was held constant. We would expect that in caseswhere there are an overwhelming number of one type of base pair (say there were fi veAT base pairs scattered amongst two thousand GC base pairs) that the alternation indexwould almost always be maximal (in the example, ten). Similarly, it would be naturalfor the distribution to be more spread out when the number of base pairs are similar.Figure 3.14 shows some cases of these distributions. Our expectations are indeed borneout, as we see the extreme cases being substantially narrower than the more balancedcases. This behaviour motivates a slight further investigation; what is the general trendof the parameters of the Gaussian fi tted function? We are expecting that as N GC is varied,at the extreme ends the standard deviation would be small and the maximum probabilityrelatively large. The mean value should increase, until for very large N GC it approachesnear saturation, and should be close to 200.In Fig. 3.15 we see more or less the expected results. The mean value increases rapidlyat fi rst, but then exhibits a very slow asymptotic tendancy to the maximum allowed valueof α = N GC ≈ N AT ≈ N GC . There is a slowing decrease towards large N GC ,suggesting that it will only extremely slowly start to resemble a very sharp distribution.The maximal value follows a similar but inverted trend, giving high probabilities for theunbalanced cases, with a minimum where N AT = N GC .62 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Models (cid:0)(cid:1) (cid:2)(cid:1)(cid:1) (cid:2)(cid:0)(cid:1) (cid:3)(cid:1)(cid:1) (cid:0) (cid:1)(cid:4)(cid:1)(cid:1)(cid:1)(cid:4)(cid:1)(cid:0)(cid:1)(cid:4)(cid:2)(cid:1)(cid:1)(cid:4)(cid:2)(cid:0)(cid:1)(cid:4)(cid:3)(cid:1)(cid:1)(cid:4)(cid:3)(cid:0)(cid:1)(cid:4)(cid:5)(cid:1) (cid:0) (cid:0) (cid:0) (cid:1) (cid:0) (cid:0)(cid:1) (cid:0) (cid:1)(cid:2)(cid:2)(cid:1) (cid:0) (cid:2)(cid:3) (cid:0) (cid:3)(cid:4)(cid:2)(cid:2)(cid:0) (cid:0)(cid:1) (cid:0) (cid:1)(cid:2)(cid:2)(cid:1) (cid:0) (cid:2)(cid:3) (cid:0) (cid:4)(cid:2)(cid:2)(cid:0) (cid:0)(cid:1) (cid:0) (cid:1)(cid:2)(cid:2)(cid:1) (cid:0) (cid:2)(cid:3) (cid:0) (cid:1)(cid:2)(cid:2) (cid:0) (cid:0)(cid:1) (cid:0) (cid:1)(cid:2)(cid:2)(cid:1) (cid:0) (cid:2)(cid:3) (cid:0) (cid:5)(cid:4)(cid:0) (cid:0)(cid:1) (cid:0) (cid:1)(cid:2)(cid:2)(cid:1) (cid:0) (cid:2)(cid:3) (cid:0) (cid:4)(cid:2)(cid:0) (cid:0)(cid:1) (cid:0) (cid:1)(cid:2)(cid:2)(cid:1) (cid:0) (cid:2)(cid:3) (cid:0) (cid:3)(cid:4) F IGURE
Probability distributions for DNA sequences with differentnumbers of GC base pairs, while N AT =
100 is kept constant. As antici-pated, the two extreme cases (very few and very many GC base pairs) cor-respond to the narrowest, tallest distributions, while the more balancedmiddle region exhibits broader, shorter distributions. Plotted points arethe exact distribution, with the curves denoting Gaussian fi ts. (cid:0) (cid:1)(cid:0)(cid:0) (cid:2)(cid:0)(cid:0)(cid:0) (cid:2)(cid:1)(cid:0)(cid:0) (cid:3)(cid:0)(cid:0)(cid:0) (cid:3)(cid:1)(cid:0)(cid:0) (cid:0) (cid:0)(cid:1) (cid:0)(cid:1)(cid:0)(cid:2)(cid:0)(cid:0)(cid:2)(cid:1)(cid:0)(cid:3)(cid:0)(cid:0)(cid:3)(cid:1)(cid:0) (cid:0) (cid:0) (cid:0)(cid:1)(cid:2) (cid:0) (cid:1)(cid:0)(cid:0) (cid:2)(cid:0)(cid:0)(cid:0) (cid:2)(cid:1)(cid:0)(cid:0) (cid:3)(cid:0)(cid:0)(cid:0) (cid:3)(cid:1)(cid:0)(cid:0) (cid:0) (cid:0)(cid:1) (cid:3)(cid:4)(cid:5)(cid:1)(cid:6)(cid:7)(cid:8) (cid:0) (cid:0) (cid:0)(cid:1)(cid:2) (cid:0) (cid:1)(cid:0)(cid:0) (cid:2)(cid:0)(cid:0)(cid:0) (cid:2)(cid:1)(cid:0)(cid:0) (cid:3)(cid:0)(cid:0)(cid:0) (cid:3)(cid:1)(cid:0)(cid:0) (cid:0) (cid:0)(cid:1) (cid:0)(cid:4)(cid:0)(cid:1)(cid:0)(cid:4)(cid:2)(cid:0)(cid:0)(cid:4)(cid:2)(cid:1)(cid:0)(cid:4)(cid:3)(cid:0)(cid:0)(cid:4)(cid:3)(cid:1)(cid:0)(cid:4)(cid:5)(cid:0)(cid:0)(cid:4)(cid:5)(cid:1) (cid:0) (cid:1) (cid:2) (cid:0) (cid:0) (cid:1) (cid:1) (cid:2) (cid:3) (cid:0)(cid:1)(cid:2) F IGURE
The variation with the number of GC base pairs N GC of (a) themean value of the distribution, (b) the standard deviation of the distribu-tion, and (c) the maximum probability of the distribution. Chapter 3 Malcolm Hillebrand 63haotic Dynamics of Polyatomic Models (cid:0)(cid:1)(cid:1) (cid:2)(cid:1)(cid:1) (cid:3)(cid:1)(cid:1) (cid:4)(cid:1)(cid:1) (cid:5)(cid:1)(cid:1) (cid:6)(cid:1)(cid:1) (cid:0) (cid:1)(cid:7)(cid:1)(cid:1)(cid:1)(cid:7)(cid:1)(cid:5)(cid:1)(cid:7)(cid:0)(cid:1)(cid:1)(cid:7)(cid:0)(cid:5)(cid:1)(cid:7)(cid:2)(cid:1) (cid:0) (cid:0) (cid:0) (cid:1) (cid:0)(cid:1)(cid:2) (cid:0) (cid:0)(cid:1) (cid:0)(cid:0) (cid:2)(cid:3) (cid:1)(cid:2)(cid:0)(cid:2) (cid:0) (cid:0) (cid:1)(cid:2)(cid:2)(cid:2)(cid:0) (cid:0) (cid:3)(cid:2)(cid:2)(cid:0) (cid:0) (cid:4)(cid:2)(cid:2) (cid:0)(cid:1)(cid:1) (cid:2)(cid:1)(cid:1) (cid:3)(cid:1)(cid:1) (cid:4)(cid:1)(cid:1) (cid:5)(cid:1)(cid:1) (cid:6)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:2) (cid:0) (cid:0)(cid:1) (cid:0)(cid:0) (cid:2)(cid:3) (cid:1)(cid:2)(cid:0)(cid:3) (cid:0) (cid:0) (cid:1)(cid:2)(cid:2)(cid:0) (cid:0) (cid:3)(cid:4)(cid:2)(cid:0) (cid:0) (cid:5)(cid:4)(cid:2) (cid:0)(cid:1)(cid:1) (cid:2)(cid:1)(cid:1) (cid:3)(cid:1)(cid:1) (cid:4)(cid:1)(cid:1) (cid:5)(cid:1)(cid:1) (cid:6)(cid:1)(cid:1) (cid:0)(cid:0)(cid:1)(cid:2) (cid:0) (cid:0)(cid:1) (cid:0)(cid:0) (cid:2)(cid:3) (cid:1)(cid:2)(cid:0)(cid:3) (cid:0) (cid:0) (cid:1)(cid:2)(cid:3)(cid:2)(cid:0) (cid:0) (cid:4)(cid:2)(cid:2)(cid:0) (cid:0) (cid:5)(cid:3)(cid:2) F IGURE
Probability distributions for three cases of total number of basepairs, for the ratio of N GC : N AT equal to (a) 1 : 1, (b) 2 : 1, (c) 6 : 1. Plottedpoints are the exact distribution, with the curves denoting Gaussian fi ts. (cid:0) (cid:1)(cid:0)(cid:0) (cid:2)(cid:0)(cid:0) (cid:3)(cid:0)(cid:0) (cid:4)(cid:0)(cid:0) (cid:5)(cid:0)(cid:0)(cid:0) (cid:0) (cid:0)(cid:5)(cid:0)(cid:0)(cid:1)(cid:0)(cid:0)(cid:6)(cid:0)(cid:0)(cid:2)(cid:0)(cid:0)(cid:7)(cid:0)(cid:0) (cid:0) (cid:0) (cid:0)(cid:1)(cid:2) (cid:8)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:3)(cid:14)(cid:5)(cid:8)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:1)(cid:14)(cid:5)(cid:8)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:5)(cid:14)(cid:5) (cid:0) (cid:1)(cid:0)(cid:0) (cid:2)(cid:0)(cid:0) (cid:3)(cid:0)(cid:0) (cid:4)(cid:0)(cid:0) (cid:5)(cid:0)(cid:0)(cid:0) (cid:0) (cid:0)(cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:0)(cid:5)(cid:1)(cid:5)(cid:2)(cid:5)(cid:3)(cid:5)(cid:4) (cid:0) (cid:0) (cid:0)(cid:1)(cid:2) (cid:6)(cid:7)(cid:8)(cid:9)(cid:10)(cid:11)(cid:3)(cid:12)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:10)(cid:11)(cid:1)(cid:12)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:10)(cid:11)(cid:5)(cid:12)(cid:5) (cid:0) (cid:1)(cid:0)(cid:0) (cid:2)(cid:0)(cid:0) (cid:3)(cid:0)(cid:0) (cid:4)(cid:0)(cid:0) (cid:5)(cid:0)(cid:0)(cid:0) (cid:0) (cid:0)(cid:6)(cid:0)(cid:0)(cid:6)(cid:5)(cid:0)(cid:6)(cid:1)(cid:0)(cid:6)(cid:7)(cid:0)(cid:6)(cid:2)(cid:0)(cid:6)(cid:8)(cid:0)(cid:6)(cid:3)(cid:0)(cid:6)(cid:9)(cid:0)(cid:6)(cid:4)(cid:0)(cid:6)(cid:10) (cid:0) (cid:1) (cid:2) (cid:0) (cid:0) (cid:0) (cid:1) (cid:1) (cid:2) (cid:3) (cid:0)(cid:1)(cid:2) (cid:11)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:3)(cid:17)(cid:5)(cid:11)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:1)(cid:17)(cid:5)(cid:11)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:5)(cid:17)(cid:5) F IGURE
The variation with total number of base pairs N = N AT + N GC of different parameters of the Gaussian fi t (a) The mean value, (b) the stan-dard deviation, and (c) the maximum probability. The fi nal property we investigate is keeping the ratio of N GC : N AT constant, andincreasing the total number of base pairs. Here the expectation is a little less obvious;the mean value should increase as N = N AT + N GC increases, but the other parametersare less clear. We see some of the distributions plotted in Fig. 3.16. It is immediatelyclear that as we saw previously, the more unbalanced ratios give narrower distributions.The larger ratios also appear to show a slower increase in the mean value.Once again the changes in parameters α , σ α and maximum probability have beenplotted, in Fig. 3.17. Here we see that in fact the mean value of the fi tted Gaussianincreases linearly in all cases. The slope does however depend on the ratio; for the ratioof 1 : 1, the slope m = m = m = N , but then slowingdown for longer sequences. As before the maximum probability more or less displaysan inverse behaviour, with a sharp decrease for small N and then decreasing very slowlyfor larger values of N . Using methods of Pólya counting theory, and extending Pólya’s Enumeration Theoremto bipartite sets using adapted cycle indices along with generating functions as formalpower series, we have constructed an extremely ef fi cient algorithm for computing theprobability distribution of the alternation index α as the numbers of AT and GC basepairs in a DNA sequence are varied. Monte Carlo simulations have validated our theoret-ical results, and show that around 5000 random realisations are suf fi cient to approximate64 Malcolm Hillebrand Chapter 3haotic Dynamics of Polyatomic Modelsthe distribution of the alternation index. This distribution has various properties whichcan be studied, and signi fi cantly can be fi tted with a Gaussian function which enables alarge amount of investigation. Most importantly for the primary task of this thesis how-ever, we now have a valid way of answering the question “what are the most and leastlikely values of α to occur in a random DNA sequence”, and we are able to quanti fi ablymeasure the impact of the alternation index and hence heterogeneity on the dynamicsof DNA molecules.Chapter 3 Malcolm Hillebrand 65 hapter 4Chaotic Dynamics of DNA Models In this chapter, we will fi nally move on to some nonlinear dynamics, and study theeffects of heterogeneity on the chaoticity of DNA molecules in the PBD model. ThemLE is computed to quantify the chaoticity, and the variation of this mLE with ATcontent is studied across a spectrum of energy densities. The in fl uence of heterogeneityis probed more deeply through studying the effect of the alternation index α that we haveintroduced in Chapter 3 on the chaoticity. Further, we study the behaviour of the mLEas the temperature increases past the melting temperature, comparing this regime to verylow temperatures. This is followed by a discussion of the distribution of deviation vector(DVD), and what can be gleaned about the chaotic hotspots in the strands, especiallyfocussing on biologically signi fi cant sequences. Finally, some similar results (mLE aswell as DVD) are produced for the sequence-dependent ePBD model, and compared tothe original model.The work in this chapter is based around the results presented in [ ] . As has been outlined in Chapter 2, the study of DNA through mathematical models isa well-established endeavour. Here we now turn particular focus to the studies of the nonlinear dynamics of DNA. As noted earlier, this is a crucial aspect of the modellingDNA behaviour, responsible for reproducing most of the signi fi cant properties of DNA.These include the biologically important formation of long-lived bubbles at transcrip-tionally active sites [ ] (more about this in Chapter 5). The nonlinearity of the modelalso introduces the localisation of energy, and particularly the formation of nonlinear lo-calised modes [ ] . The oscillations near the minimum of the Morse potential well leadto periodic breather solutions, related to the low-frequency vibrational modes presentin DNA molecules. The existence of these localisations in the absence of disorder canin some sense be attributed to “dynamical disorder”, due to the nonlinear potential.This nonlinearity causes certain sites to become self-trapping high amplitude, low fre-quency breathers, which are unable to disseminate their energy due to the dynamicalrestrictions on the phonon spectrum [ ] . The presence of these nonlinear fl uctuationsimpacts charge transport in DNA molecules where “vibrational hotspots”, regions witha high propensity for bubbles or breathers to form (i.e. localised nonlinear modes), cancause anomalous subdiffusion of the charge [
93, 94 ] . It is again interesting that the non-linearity can fi ll the role typically associated with disorder, as in these studies of charge66haotic Dynamics of Polyatomic Modelstransport only homogeneous DNA polymers have thus far been considered [ ] . Mostcritically of course, the nonlinearity leads to what we noted previously, the abrupt de-naturation curve which is so important to meaningfully modelling the biodynamics ofDNA [ ] .Here we should also discuss the phase transition mentioned in Chapter 2. The exis-tence of phase transitions in one dimensional systems with short-range interactions hasbeen the topic of much discussion in the past. In a large number of cases, these transi-tions are forbidden due either to the fi ndings of van Hove [ ] , or by the argumentationof Landau [ ] . As such, the emergence of the DNA melting transition in the PBDmodel as an apparent fi rst order phase transition is a slight surprise – even in the homo-geneous cases studied early on [
44, 45 ] , the abrupt denaturation process demonstratedproperties of a true thermodynamic transition. However, due to a number of factors,the arguments against phase transitions in one dimensional systems to not actually applyto the PBD model in question. In particular, van Hove’s theory demands that the poten-tials all be functions of ( y n − y n − ) , so the presence of an on-site potential which is only afunction of y n means that we cannot necessarily apply this theory [
97, 98 ] . The secondrequirement for the absence of a phase transition in a one-dimensional system is that theenergy of the domain walls between different regions in the lattice should be fi nite [ ] .In the DNA case, since the energy of domain walls between regions of open and closedbase pairs is in fi nite, the phase transition is not in fact forbidden by this theory [
47, 97 ] .Thus we have a true phase transition [ ] , and further study of the dynamics of the PBDmodel of DNA to understand this phase transition has provided useful results [ ] . Focusing in even more closely from general nonlinear dynamics to the chaotic behaviourof the PBD model, there has been fairly little research into identifying and quantifyingthe chaotic behaviour of the model. However, an initial investigation of the chaoticdynamics has been performed [ ] . In that investigation, the chaoticity is quanti fi edthrough the mLE, which is estimated using a combination of the transfer integral methodand Riemannian geometry. Here the harmonic coupling is used to expedite the analyt-ical calculations, but numerical results are presented for both the harmonic and anhar-monic couplings. The result of that work was the fi nding that the mLE could serveas a dynamical order parameter, showing all the requisite characteristics. The mLE in-creases with energy density (or equivalently temperature), and demonstrates a changein behaviour near a critical value corresponding to the phase transition. This result isshown in Fig. 4.1, reproduced from [ ] . The harmonic coupling was found to lead to asecond-order phase transition, while the anharmonic coupling shows a more fi rst-order-like transition. This additional thermodynamical signi fi cance associated with the mLEprovides another point of importance to the study of the chaotic dynamics of DNAmolecules [ ] .It is with this in mind that the work presented here was undertaken. These earlyresults indicate that there is value in a thorough understanding of the chaotic dynamics atplay in DNA dynamics, and the work of Barré and Dauxois [ ] with the homogeneousmodel suggests that there may be more depth to be gained by studying more completeversions of the model with the advantage of modern computational resources.Chapter 4 Malcolm Hillebrand 67haotic Dynamics of Polyatomic Models In order to build a statistically robust picture of the mLE’s behaviour, we need to per-form suf fi ciently many simulations and combine them effectively. The procedure herewas to perform, at each energy density, and for each AT percentage, one hundred totalsimulations. These simulations are further split up into ten realisations (base pair con- fi gurations) which each then have ten runs with different initial conditions. This allowsthe separate investigation into whether the individual realisation give signi fi cantly differ-ing mLE values, which led to the investigation of the alternation index, and the resultspresented in Section 4.3. For each of these simulations, the fi nite time mLE (ftmLE) χ ( t ) was computed and its evolution tracked. Figure 4.2 shows the evolution of onesuch run. In this plot we see that χ ( t ) displays the behaviour expected for a chaotic tra-jectory, initially decreasing before saturating to a positive value. This fi nal value is theestimate for the mLE of this trajectory.The ultimate article of interest in this is of course essentially only this fi nal value ofthe ftmLE. It is however instructive to observe these individual curves, as the choice ofthe fi nal integration time needs to be made with the complete saturation of the mLEin mind. In Fig. 4.3(a) we can see that after 10 ps, there is still a noticeable deviationbetween the various curves. We would like, as far as possible, to eliminate any kind of fi -nite time inaccuracy from our calculations, and be sure that the only source of spreadingof values is actually inherent to the system. To this end, after looking at the curves andchecking their variation and saturation, the fi nal time of 10 ps was chosen as the pointat which the ftmLE has practically ceased to change in time, as here there are very fewnoticeable changes in the curves. Figure 4.3(b) shows the ftmLE averaged over multipleinitial conditions and realisations 〈 χ ( t ) 〉 (shown by the solid dark line), along with theindividual curves that make up the average (shown by lighter grey lines). The error barsin the average ftMLE are simply the standard deviation of the 100 runs making up theaverage. This shows quite clearly that the fi nal values of χ ( t ) are very similar in all cases,giving reason to believe that it is reasonable to assign a single mLE value to each pointwith given E n and P AT .In all cases, the initial conditions are given by the displacements y n = fi nal measurement of the mLE at a given energy den-sity, for a given AT percentage, we average the last points of the ftmLE for each run,giving an average value χ . The standard deviation of this average gives the uncertainty,completing the measurement.Chapter 4 Malcolm Hillebrand 69haotic Dynamics of Polyatomic Models t (ps) . . . . . χ ( t )( p s − ) F IGURE
The time evolution of the fi nite time mLE χ ( t ) at E n = P AT = t (ps) . . . . . χ ( t )( p s − ) (a) t (ps) . . . . . χ ( t )( p s − ) (b) F IGURE
The time evolution of the fi nite time mLE at E n = P AT = P AT fairly well. The error bars of the average valueare the standard deviation of the 100 runs making up the averaged curve. Having this machinery to estimate the mLE χ for any chosen parameters, we can nowlook at some of the results. The fi rst thing to consider is just a simple computation of χ for a large number of values of E n , for different AT percentages. Given the dynam-ical differences between the base pair types encoded into the model, we would expectdifferent behaviours between percentages.Figure 4.4 shows this result. We see that there is a different behaviour dependingon P AT , as well as an increase in the chaoticity with increased energy density. Thereare however two distinct regimes at play here: The fi rst where the mLE increases withenergy density, and then the second where the mLE fl attens out before rapidly droppingto zero. We will split the discussion to analyse these sections. This is expedited by the factthat there is a physical justi fi cation for separating the regions – the fi rst is the physicallysigni fi cant region before the melting transition, the complete separation of the DNA70 Malcolm Hillebrand Chapter 4haotic Dynamics of Polyatomic Models .
00 0 .
02 0 .
04 0 .
06 0 . E n (eV) . . . . . . χ ( p s − ) P AT = 0% P AT = 10% P AT = 30% P AT = 50% P AT = 70% P AT = 90% P AT = 100% F IGURE
The average mLE χ for several AT percentages, across a spec-trum of energy densities E n . Line connections guide the eye only. Theerror bars correspond to the standard deviation in the average calculationof χ . double helix, and the second is what happens after this complete denaturation. Pre-Melting mLE Behaviour
The fi rst point here is to establish where in fact the melting transition takes place. Themelting temperature (in Kelvin) of a general DNA molecule in the PBD model is givenby [ ] T m = − P AT , (4.1)following the experimentally veri fi ed linear variation, and a melting temperature of365K for a pure GC sequence, and 325K for a pure AT sequence, again in agreementwith experiments [ ] .So now we know that we should cut the results at the melting temperature appro-priate to each P AT . But we do need to estimate the temperature of the system, since weare performing simulations in the microcanonical ensemble where energy rather thantemperature is kept fi xed. Thus we use the usual representation of temperature as theaverage kinetic energy in a one dimensional system scaled according to the Boltzmannconstant T = 〈 K n 〉 k B , (4.2)where the average kinetic energy 〈 K n 〉 is calculated as 〈 K n 〉 = mn n (cid:1) j = p j , (4.3)with n particles of mass m .Thus for each simulation, we calculate the temperature, and then take an ensembleaverage for all (100) runs at that energy density and arrive at an estimation of the averagetemperature corresponding to that energy density. With that, it is possible to identifyChapter 4 Malcolm Hillebrand 71haotic Dynamics of Polyatomic Models .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 . E n (eV) . . . . . χ ( p s − ) P AT = 0% P AT = 10% P AT = 30% P AT = 50% P AT = 70% P AT = 90% P AT = 100% F IGURE χ for each AT percentage, as a function of E n , as in Fig. 4.4, butstopping when melting occurs. the point where the temperature crosses the melting point, and then to stop the curvein the plot of the mLE before that point. The plot of χ as a function of energy densitybefore the melting point is given in Fig. 4.5.Here we see that indeed the whole section of the plot where χ drops to zero, as wellas some of the preceding decrease, is omitted. This now leaves us with the physicallyrelevant region of interest. In this region there are apparently three distinct behavioursvisible in the mLE. While there is still an overall trend of increasing chaoticity with en-ergy density, we see that at lower energies, generally sequences with higher AT contentare more chaotic. Below E n ≈ (cid:2) E n (cid:2) P AT corresponds to a larger χ value.This is followed by a small region where the effect of the AT / GC composition is min-imal; all cases more or less give similar values of χ between 0.025 (cid:2) E n (cid:2) P AT value. Thereafter comes a region where the rate of increase in χ slows down, as the curves fl atten out towards the melting point. There is also a fairlyclear trend here that the low energy behaviour has switched around, and now higherGC content corresponds to higher χ .The effect of the base pair composition on the chaoticity is consistent, with the moreweakly bonded AT base pairs tending to exhibit more chaotic behaviour at low energies,and the GC base pairs taking longer to break apart and thus showing more chaoticity athigher energies.Figure 4.6 shows the variation of the mLE as a function of the temperature, where thetemperature is found from (4.2). The shape of the curves is slightly different, as the rela-tionship between energy density and temperature is nonlinear (see Section 5.2.2). Duealso to the fact that the AT percentages has an effect on the energy-temperature relation-ship, points in Fig. 4.6 at the same energy are not necessarily at the same temperature,especially as energy increases.It is also good to note that the overall behaviour exhibited here is in agreement withthe results of [ ] , showing the same shape of curves, but of course different values dueto the heterogeneity and use of different parameters. In the work of Barré and Dauxois,72 Malcolm Hillebrand Chapter 4haotic Dynamics of Polyatomic Models T (K) . . . . . χ ( p s − ) P AT = 0% P AT = 10% P AT = 30% P AT = 50% P AT = 70% P AT = 90% P AT = 100% F IGURE χ for each AT percentage, as a function of the temperature T ,stopping at the melting temperature. The temperature is calculated accord-ing to (4.2) for each energy density. only a homogeneous version of the PBD model was considered, thus omitting the dif-ferent behaviours evident in the more general model. The other major difference comesfrom the use of different units, and particularly the resulting time units, which meansthat our results differ from those in Fig. 4.1 by two orders of magnitude. Post-Melting mLE Behaviour
Some time after the melting point we see that the average mLE decreases suddenly tozero in Fig. 4.4. This drop to zero of course signi fi es a transition to a completely regularsystem. We see that there are also some intermediate points where the average mLEdrops, but not all the way to zero. These can be explained by seeing that the system isat this point mixed – some initial conditions lead to chaotic trajectories, while otherscorrespond to regular trajectories.We do fi rst have to note that this post-melting regime is physically irrelevant – themodel makes no claim to accurately model the dynamics past the melting point. In fact,given that the two strands can completely separate and decouple, this situation is in somesense impossible to model. We can brie fl y consider the model from a purely dynamicalperspective though, and look at what happens at this transition.The regularisation of the motion can be explained from the potential. Recallingthat the nonlinearity of both the coupling and the on-site Morse potential arises fromexponential terms of the form e − y n and e − ( y n + y n − ) , we can see that as these displacementsincrease beyond a certain point, the exponential terms rapidly decrease to near zero,and the remaining terms are simply a constant and a harmonic coupling. This systemis now completely linear, and we would expect completely regular behaviour. Fromthis observation we can also understand the reason for the χ to increase more slowlyat higher energy densities. As more of the base pairs remain in widely separated states(see Chapter 5 for a detailed discussion of the fraction of open base pairs), these basepairs are in some sense contributing linearly to the overall dynamics. So while the basepairs that remain dynamically bound are increasingly active, and thus more unstable andChapter 4 Malcolm Hillebrand 73haotic Dynamics of Polyatomic Models t (ps) − χ ( t )( p s − ) F IGURE
The evolution of ten fi nite time mLEs χ ( t ) for P AT =
70% at E n = chaotic, the base pairs that reach the breaking point and separate form a stable dynamicalbehaviour in wide-open bubbles. They can of course fl uctuate between open and closed,and different base pairs will contribute differently to the chaoticity at different times, butthe average picture is more or less the same.Figure 4.7 shows the evolution of a few χ ( t ) values for different initial conditionsat E n = P AT =
70% case, which is above the melting point. Here wecan see the hybrid behaviour of the system. Some of the ftmLE estimations remainpositive, indicating that chaotic behaviour remains prevalent and the double strand isnot completely denatured. Others however show the steady decrease towards zero, at aslope of − Having seen that the base pair makeup of the sequence is important for the dynamicsof the system, we can now turn to a more detailed investigation of the heterogeneity.74 Malcolm Hillebrand Chapter 4haotic Dynamics of Polyatomic ModelsOne question that arises from the previous section is “where does the spread of mLEvalues come from?” Is it just a fact that slightly different arrangements of base pairs andslightly different initial conditions give different values for χ ? Is this variation inherentin any molecules, even very similar ones? Or, as we will investigate here, can we breakdown this spread of values into different classes of physically similar molecules, by con-sidering the “clumping” of base pairs? To study this “clumping”, we will make use ofthe alternation index α introduced in Chapter 3.So now at last we can reap bene fi t of ploughing through the mathematics required toproduce a probability distribution for α . What we really want to know is how the het-erogeneity affects the chaotic dynamics. For instance, if we had a DNA double strandconsisting of alternating AT and GC base pairs, we might expect the dynamics to berather distinct from that of a molecule with two solid chunks of AT and GC pairs.The alternation index allows us to quantify this heterogeneity, and so we can comparethe different extreme cases of very well mixed, homogeneous sequences, and extremelychunky, heterogeneous, sequences. The probability distribution enables the identi fi ca-tion of these extremes, as well as the most likely value of α . In principle, one couldproduce a plot showing the variation of χ with α for all possible values of α , but it ismuch more interesting to look at the extremes and see if there is any visible differencethere before continuing with a series of long computations. So here results are presentedfor a collection of AT percentages, comparing the cases of1. The most likely value of α .2. Extremely small values of α (more heterogeneous sequences).3. Extremely large values of α (more homogeneous sequences).Figure 4.8 show the value of χ for fi ve different AT percentages. In each plot, thevariation of χ with energy is shown for fi ve values of α . In order to cover the spec-trum of possibilities we use the most likely value, the lowest possible value, the largestpossible value, and two extra points near the extremes: the largest possible value lessfour, and the lowest possible value plus four. This means we will have a set of fi ve cases, { α min , α min +
4, ¯ α , α max − α max } . These give an idea of the behaviour in the threeregimes (unmixed, most probable, and well-mixed). We expect that the most probablecase will re fl ect similar behaviours to the overall picture presented in Section 4.2, sincerandomly chosen sequences will of course be drawn more often than not from the centreof the distribution.It is apparent from Figs. 4.8(a) and 4.8(e) that where there is a signi fi cant proportionof a single base pair type, the alternation index does not signi fi cantly affect the chaoticity.Thus we can see that if there are only ten AT (GC) base pairs in a sequence of ninety GC(AT) pairs, it does not make a difference to the chaoticity where the ten remaining basepairs are scattered. In the intermediate cases however [ Figs. 4.8(b),4.8(c) and 4.8(d) ] , wedo see some clear differences between the regimes. It is of course most notable in the50% case, where the balanced number of base pairs lends extra signi fi cance to each oftheir position, we will fi rst consider Fig. 4.8(c). Here the two well-mixed, high α cases(red and purple lines) clearly exhibit a smaller mLE than the other, more heterogeneous,cases. This is most pronounced at low energies, where below E n = α =
50 case is more than double that of the α =
96 and α =
100 cases. This continuesat higher energies, although less extremely. The green α =
50 threads its way between theChapter 4 Malcolm Hillebrand 75haotic Dynamics of Polyatomic Models .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . E n (eV) . . . . . . χ ( p s − ) (a) α = 2 α = 6 α = 16 α = 18 α = 20 0 .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . E n (eV) . . . . . . χ ( p s − ) (b) α = 2 α = 6 α = 42 α = 56 α = 600 .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . E n (eV) . . . . . . χ ( p s − ) (c) α = 2 α = 6 α = 50 α = 96 α = 100 0 .
00 0 .
01 0 .
02 0 .
03 0 . E n (eV) . . . . . . χ ( p s − ) (d) α = 2 α = 6 α = 42 α = 56 α = 600 .
00 0 .
01 0 .
02 0 .
03 0 . E n (eV) . . . . . . χ ( p s − ) (e) α = 2 α = 6 α = 16 α = 18 α = 20 F IGURE χ for fi ve different values of α , for fi ve AT percentages. (a) P AT = P AT = P AT = P AT = P AT = α values are shown in the legend, corresponding to twoextremely low values, two extremely high values, and the most probablevalue.
76 Malcolm Hillebrand Chapter 4haotic Dynamics of Polyatomic Models .
00 0 .
01 0 .
02 0 .
03 0 . E n (eV) . . . . . . χ ( p s − ) P AT = 0% P AT = 50% , α = 100 P AT = 50% , α = 50 F IGURE χ across the energy density spectrum, comparing the extremelywell-mixed 50% AT sequences with α =
100 (orange) with the pure GCsequence (blue). The α =
50 curve for the 50% AT case is also shown(green) to show the most common behaviour of this percentage. high α lines and the low α lines, which are in turn substantially more chaotic. The blue α = α = α =
50 case, again with an emphasis on the lower energies. Herethe difference lasts all the way up to close to E n = α ) lead to distinctly morechaotic behaviour. Figs. 4.8(b) and 4.8(d) show for the P AT =
30% and P AT =
70% casesrespectively that the same trend holds true that the green α =
42 case (the most probable)stays between the two extremes, while the α = α = α is relatively large is also apparent.For instance, where in Fig. 4.8(c) the green line is, in the beginning at least, similarly farfrom both extremes, in Figs. 4.8(b) and 4.8(d) the most probable case shows a de fi niteaf fi nity for the higher α cases, with a much larger gap between the orange α = α =
42 line. This is most signi fi cant for lower energies, as we see a regionabove E n ≈ χ , the lines do not really separate again, and we see moreor less the same behaviour for all values of α up to the melting point.Overall, these results are fairly consistent – the more well-mixed a sequence is, the lesschaotic it is. This sits well with the intuition that a lattice with the disorder evenly spreadshould be more stable than a lattice with extended regions of weaker on-site potentials,where larger displacements are prone to occur.What is also notable is to look at the particular behaviours of the extreme cases. Letus focus in on the 50% AT case, as the effects of α on the mLE are strongest here. Whatwe would like to see is how the extreme cases of α compare to other AT percentages.Figure 4.9 shows the most stable case for a sequence of 50% AT base pairs, α = .
00 0 .
01 0 .
02 0 .
03 0 . E n (eV) . . . . . χ ( p s − ) P AT = 100% P AT = 50% , α = 2 F IGURE χ across the energy density spectrum, comparing the com-pletely unmixed 50% AT case with α = (orange curve), compared to the pure GC sequences (blue curve). This emphasises howextremely stable the very homogeneous molecules are at low energies – here the P AT = α =
100 curve is well below the pure GC curve up to large energy densities. We didsee however, in Fig. 4.5 that at very low energy densities, the pure GC case is not the moststable. In that regime we should be comparing to the 50% case itself, where the mostprobable case is shown in Fig. 4.9 as the green curve. So we can combine this comparisonwith the result from Fig. 4.8(c) showing the large α cases to be noticeably more stablethan the most probable case, and say that the most stable DNA con fi gurations at lowenergy densities (or equivalently, low energies) are relatively homogeneous sequenceswith extremely well-mixed base pairs.Additionally, we can study the opposite extreme – what happens when we have acomplete separation of AT and GC base pairs ( α = χ values between the P AT = α = P AT = fi cantly tothe chaoticity. As the energy increases, and the melting point of the pure AT sequenceis being reached, the GC cluster reluctantly contributes to the dynamics, and we see themLE increases past the value for the pure AT case, as the GC base pairs are excited intononlinear behaviour.This ties into the point that the length of the DNA sequence does not affect themolecule’s chaoticity – so here the dynamics are approximately representative of a fi ftybase pair pure AT sequence, which is more or less the same as the one hundred base paircase. Intuitively this is justi fi able, as we are saying that the chaotic dynamics are beingdominated by the fi fty AT base pairs, which are more easily excited than the GC cluster.So here we see that in just the single heterogeneous case of DNA sequences made up ofhalf AT and half GC base pairs we can more or less see both the most and least stablebehaviour exhibited by all DNA sequences, solely by varying the alternation index.78 Malcolm Hillebrand Chapter 4haotic Dynamics of Polyatomic Models .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 . E n (eV) . . . . . χ ( p s − ) AdMLPRandom (a) .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . E n (eV) . . . . . χ ( p s − ) Lac OperonRandom (b) F IGURE (a) χ values across the energy density spectrum for the AdMLP(blue), and random sequences with the same AT / GC composition, P AT = χ values for the Lac Operon (blue) and random sequences with P AT = This con fi rms the early impression that the heterogeneous nature of DNA moleculeshas a strong in fl uence on its dynamics and stability. Apart from the general effects ofthe AT / GC composition noted earlier in Section 4.2, we see here that the particularordering of these base pairs is important in determining the chaoticity of the sequence.The region in the energy / temperature spectrum where all AT / GC compositionsyield the same mLE is noteworthy, especially as this trend survives into the more de-tailed study through the alternation index, with all cases there also showing this samebehaviour. It is also very interesting to note the effect of the alternation index on thelower temperature region, where there is more scope for variation in the mLE, suggest-ing that the chaoticity of DNA sequences is more sensitive to the base pair ordering atlow temperatures. By varying α in sequences of equal parts AT and GC content we cansee more or less the entire spectrum of chaotic behaviour, from behaviour similar topure AT sequences to behaviour showing properties of pure GC sequences. This doessuggest that there is relatively signi fi cant variation possible in the chaotic dynamics duepurely to the particular arrangement of base pairs in the molecule. To complete the analysis of the effect of heterogeneity, and generally sequence-dependence,on the chaotic dynamics of the PBD model, we look at two biologically signi fi cant DNAsequences and compare their mLE values with those of random sequences with the samecomposition, i.e. exactly the same base pairs, with exactly the same alternation index,but otherwise arranged randomly within these constraints. The sequences chosen are an86-base pair segment of the adenovirus-associated major late promoter (AdMLP), and a129-base pair sequence from the lac operon, in both cases taken from around the tran-scription start site (TSS), the site along the strand where transcription begins [ ] . Thesesequences are well-known transcriptional promoters, and have been studied in the con-text of bubble dynamics and DNA transcription activity [
50, 100 ] . The precise sequenceof the upper strand of each promoter is given below.Chapter 4 Malcolm Hillebrand 79haotic Dynamics of Polyatomic Models .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . E n (eV) . . . . . χ ( p s − ) P AT = 0% P AT = 10% P AT = 30% P AT = 50% P AT = 70% P AT = 90% P AT = 100% F IGURE
The average mLE 〈 χ 〉 as a function of energy density, for theePBD model. Results are shown up to the melting point. This fi gure is thecorresponding result to Fig. 4.5, but now for the ePBD model. AdMLP:
GCCACGTGACCAGGGGTCCCCGCCGGGGGGGTATAAAAGGGGGCGGACCTCTGTTCGTCCTCACTGTCTTCCGGATCGCTGTCCAG lac operon:
GAAAGCGGGCAGTGAGCGCAACGCAATTAATGTGAGTTAGCTCACTCATTAGGCACCCCAGGCTTTACACTTTATGCTTCCGGCTCGTATGTTGTGTGGAATTGTGAGCGGATAACAATTTCACACAGG
We note that the AT percentage of the AdMLP is P AT = lac operonis P AT = fi cant difference in the χ values between the promoters and randomsequences. All the variations that are present are within uncertainty of each other, andfor the most part the lines are more or less identical. This suggests that there is nospecial chaotic behaviour inherent in the biological sequences. The chaotic dynamicsare controlled by the heterogeneity, and particularly the alternation index, and whilethere is still some variation of the chaoticity of the sequences (as evidenced by the errorbars), the biological signi fi cance of a sequence does not impact the value of the mLE. We will present a very short comparison between the results of mLE computations forthe PBD model and the ePBD model. As we will see, the outcomes are very much thesame for both models, which limits the usefulness of an exhaustive comparison betweenthem.Figure 4.12 shows the same results we saw in Fig. 4.5, but recomputed for the ex-tended model. The overall behaviour is identical; we have the same three regions of80 Malcolm Hillebrand Chapter 4haotic Dynamics of Polyatomic Models .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . . . . . . χ ( p s − ) (a) PBDePBD0 .
01 0 .
02 0 .
03 0 . . . . . . χ ( p s − ) (b) .
01 0 .
02 0 .
03 0 . E n (eV) . . . . χ ( p s − ) (c) F IGURE χ as a function of energy density E n for the PBD (blue)and ePBD (orange) models, with (a) P AT = P AT = P AT = χ for the two models are withinuncertainty of each other. Chapter 4 Malcolm Hillebrand 81haotic Dynamics of Polyatomic Modelsdistinct behaviour, and the same increase with energy density.Further comparison is easier to do on a case by case basis, in order to investigatewhether or not there is any systematic difference. Figure 4.13(a) shows the estimation of χ for the PBD and ePBD models, across the energy spectrum for P AT = k , the ePBD model predicts very slightly more stablebehaviour than the PBD model. Since the two results are always within uncertainty ofeach other however, we cannot make any strong conclusions from this, with the generaltakeaway still being that the two models give more or less identical results.The heterogeneous cases all exhibit similar behaviours; the P AT =
50% result isshown as a representative case. Figure 4.13(b) shows that once again the ePBD modelgives a marginally less chaotic behaviour than the PBD model, but the lines once againfollow the same pattern, and remain within uncertainty of each other.The pure GC case is shown in Fig. 4.13(c), once again displaying the same tenden-cies. Here however the difference, particularly at larger energies, is slightly greater thanin the pure AT case, suggesting that the sequence-speci fi c coupling constant has a moresigni fi cant variation from the PBD value. This is the only case where some points areactually different even up to uncertainty; thus we could say that near the melting tran-sition, the ePBD model gives a lower estimation of the chaoticity than the PBD model,for the pure GC case.Overall, while the impression is clear that the exact nature of the sequence-dependencedoes not have a strong effect on the chaoticity of the DNA molecule, it is also apparentfrom these results that the ePBD model yields slightly less chaotic behaviour than thePBD model, in all cases showing smaller mLE values. The fi nal element we will mention in this study of the chaotic dynamics of DNA are thedeviation vector distributions (DVDs). Here we want to investigate any link betweenthe physically observable dynamics – displacements, and particularly the formation ofbubbles – and the presence of strong nonlinearity. For instance, would regions of largedisplacement correspond to regions of signi fi cant nonlinearity? Or would they perhapsbe devoid of chaotic behaviour, since the distant base pairs are only interacting relativelylinearly? So our investigation of the DVD follows along these lines: Observing thedisplacements and the DVD in different contexts, investigating how different energiesaffect the dynamics, and how base pair sequences affect the dynamics. To investigate the general behaviour of the DVD in the low-energy regime, it is instruc-tive to consider both individual cases, and averages over multiple runs. Figure 4.14 showsthe time evolution of the DVD spectrum along with the displacement y n , for a particularcase with P AT =
30% at E n = fi cally to understand the effect of the particular base pair sequenceon the DVD, we consider only averages over many initial conditions, so as not to beoverly distracted by the effect of a particular initial condition on the localisation of theDVD.Figure 4.19 shows the average DVD (upper panel) and displacement (lower panel)for the AdMLP, over a period of time. In general, this shows the characteristics foundin Section 4.5.2, as these results are in the same energy regime. However, the dynamicalbehaviour of the AdMLP has been studied previously, and in particular the regions ofhigh and low bubble probability [ ] , which gives us additional information that allowsus to place the DVD results in context. Apart from the clear large displacements in theAT base pairs, particularly the TATA box around base pair −
29 (a sequence of base pairsindicating the position of the TSS and the transcription direction [ ] ), we see that there isa region of large displacements centred around the TSS. Here the base pairs are labelledaccording to the biological convention, where the base pairs are numbered relative tothe TSS, which is conventionally given the label 1. So for instance in the AdMLP, the64 th site is the TSS, and so this is marked with 1. The downstream sites are labelledincreasingly in the negative direction, relative to this site, so site 0 would be labelled −
63. Upstream sites are labelled increasing in the positive direction so that site 84 islabelled with 21. The DVD generally avoids the TATA box, where we have consistentlarge displacements, which is in agreement with our earlier discussion of the DVD athigher temperatures. In addition to this however, we see that the DVD also avoids theregion downstream of the TSS, which is known to be a low-probability region for largebubbles [ ] . This again leads us to the conclusion that the DVD does not localise nearlong-lived, large, consistent bubbles, or in regions of little activity, but rather preferringthe less predictable fl uctuations in the DNA strand.The same results for the lac operon, shown in Fig. 4.20, are still instructive. Onceagain, we know the bubble probability pro fi le for this sequence, allowing us to comparethis with the DVD [ ] . Here we observe that the main areas avoided by the DVDappear to be around the TSS, and near sites −
59 and below −
39. These once again cor-respond to fairly high bubble probabilities (in the TSS case) and relatively low bubbleprobabilities (the downstream segment). Broadly these results support the understand-ing that the DVD avoids the extreme cases, with the DVD not concentrating at regionsof large amplitude bubbles, or at regions of very small displacements.It is worth noting at this point that the interpretation of the DVD is that the nonlin-earity of the system is strongest at sites where the DVD is concentrated [
23, 34 ] . Thismeans that the DVD results here imply that the system is least chaotic in regions wherethe openings are systematically large, or where they are systematically small. Intuitivelythis matches with our previous analysis of the two near-linear borders of the system –the case of large displacements, where the potential effectively linearises, and the caseof very little movement at low temperatures. The regions of large long-lived bubbleswill exhibit approximately linear behaviour, and the regions of small displacements willbe relatively less chaotic than other sites where the nonlinearity is being more stronglydrawn out by the larger displacements.88 Malcolm Hillebrand Chapter 4haotic Dynamics of Polyatomic Modelscentre with a 10-base pair island of AT base pairs (Fig. 4.22), and GC base pairs (Fig. 4.23)respectively, we can fi nd if there is something more to be seen about the effect of AT / GCarrangement on the DVD. In order to make the comparison as fair as possible, the exactsame initial conditions were used in all three of these cases shown here (Figs.4.21, 4.22and 4.23), but with the different base pair realisation.In Fig. 4.22 it is apparent that the deviation vector has no overwhelming preferencefor its location. The displacement once again shows what we expect, a tendency to focusat the AT base pairs, resulting in a channel of high displacement in the centre. The DVDnow does follow this excited region, showing a moderate tendency to focus in this centralisland. Around 14% of the DVD is focussed in this 10-base pair region, demonstratinga slightly above-average concentration, but nothing extraordinary. What is interestingis that once again the striping of the DVD reveals that despite the 45–55 base pair splitin favour of AT base pairs, overall 48% of the DVD is focussed at GC base pairs. Thismeans that outside the central region, we more or less maintain the previous fraction of55% of the DVD at GC base pairs (48 / ≈ fi nity or aversion to particular base pair types. In all thecases examined in this section, it is clearly dependent not only on the particular base pairrealisation, but also on the initial condition. By averaging over initial conditions we canreveal the dependence on the realisation alone, since the effect of the initial conditionsthemselves is minimised by the averaging, showing the tendency of the DVD to avoidconsistently highly displaced regions. In this chapter, we investigated the chaotic dynamics of the PBD model of DNA, as wellas the sequence-dependent ePBD model. Using numerical methods, we studied the effectof AT and GC composition on the chaoticity of DNA sequences through the estimationof the mLE. We computed the mLE for DNA strands with several AT percentages, aswell as fi xing the alternation index and studying how the mLE can be affected by theordering of base pairs even if the total number of AT and GC base pairs are the same.Additionally, we computed the DVD for different DNA sequences, and using singlecases as well as averaged data studied the localisation of chaoticity in various cases.Simply investigating the chaoticity of the model as a function of energy, the compu-tation of the mLE informs us that the chaoticity increases with energy, with a slowingin the increase near the melting temperature. The effect of the percentage of AT basepairs in the sequence on the chaoticity is signi fi cant, with the low-energy regime yield-ing more chaoticity from AT-rich sequences, before a mid-energy region between energyChapter 4 Malcolm Hillebrand 91haotic Dynamics of Polyatomic Modelsdensities of 0.025eV and 0.035eV shows little effect from changing the AT content, and fi nally at higher energies the GC-rich sequences become more chaotic.The heterogeneity can be probed further, by use of the alternation index α . Com-puting the mLE at different values of α we fi nd that the more well-mixed a sequenceis, the higher the value of α , the less chaotic it is. Less heterogeneous sequences, withfew alternations, provide more chaotic behaviour than the most likely values of α . Forsequences consisting primarily of one base pair type (for instance 90% AT base pairs),the alternation index has very little effect.The addition of more sensitive sequence-dependence through the use of the extendedPBD (ePBD) model yields no signi fi cant changes to the chaotic dynamics.Investigation of the deviation vector distribution (DVD) yields that the nonlinearityin the system focusses near regions of large displacement, but not at those regions. Inparticular, the DVD avoids regions of consistently large bubbles, such as near the TATAbox of biological promoters, as well as regions of very little fl uctuational activity, wherethe displacements are consistently small. The DVD concentrates rather in regions ofunpredictable and changing displacements. The DVD jumps around between these re-gions of changing displacement, often without apparent cause. However, at very lowtemperatures, the DVD concentrates in homogeneous islands in the sequence. This ap-pears to be independent of the base pair that the island is comprised of, but the DVDwill concentrate in islands of both AT and GC base pairs rather than spread through themore heterogeneous segments of the strand.The effect of heterogeneity in DNA is clearly visible in the chaotic dynamics, andthere are consequently several possibilities for future investigations. For instance, thestudy of different models, such as modi fi cations of the PBD model to account for thetwisting of the DNA strand explicitly [ ] , or other adaptations. There is also scopeto study the chaotic dynamics through the computation of covariant Lyaponov vec-tors (CLVs) [ ] , which have the potential to provide valuable information about thenonlinearity of the models. Further investigation of the DVDs, possibly together withCLVs, is another direction that could prove fruitful. More insight into the relationshipbetween bubbles and concentrated nonlinearity holds promise for the better understand-ing of DNA dynamics.92 Malcolm Hillebrand Chapter 4 hapter 5Melting and Bubble Properties ofDNA The notion of thermally induced openings in DNA, or bubbles, has been introduced inboth the introductory DNA chapter (Chapter 2) and mentioned in the context of theDVD in Chapter 4. We will now move into a thorough investigation of these bubbles,and try to build an understanding of bubble occurrences probabilities, as well as bubblelifetimes in typical DNA sequences. The biochemical signi fi cance of bubbles has beenexplored extensively in the literature, with results being found indicating a relationshipbetween transcriptionally signi fi cant sites – the promoter TSS, as well as notable reg-ulatory sites – and the spontaneous formation of bubbles [
50, 54, 104, 105 ] . Againstthis background of biological interest, in this chapter we will present results pertainingto the formation of bubbles in DNA strands. This includes a relationship between thetemperature and energy density, a distance threshold for considering base pairs open,and then distributions for bubble lengths in DNA and bubble lifetime probability dis-tributions. A comparison of bubble lifetimes for a particular biological promoter withresults obtained for general sequences is also presented, adding context to our fi ndings.The results of in this chapter are based on the research presented in [ ] . As an introduction to the study of bubbles, we can look at a number of essential as-pects of bubbles that have warranted attention in the past. The biological importance oftranscriptional sites has meant that many projects have been devoted to studying the dis-tribution of bubble occurrence across the DNA sequence [
50, 52–54, 75, 100, 101, 105,107–111 ] . These investigations look at different features of thermally induced openings(i.e. bubbles) in DNA molecules, and compare dynamically active sites with experimen-tally determined transcription sites. One of the most signi fi cant suggestions to arisefrom studying the probability of formation and lifetimes of bubbles at the TSS of knownpromoters, and comparing this to known nonpromoter regions, is that DNA dynam-ically directs its own transcription [
50, 54 ] . By performing numerical simulations atphysiological temperatures, the probabilities and lifetimes of bubbles can be estimated,and from these distributions particular conclusions can be drawn – for instance, that“soft” regions in promoters, i.e. regions with a high propensity for thermal openings,correspond to transcription binding sites [ ] . While the majority of studies (includ-93haotic Dynamics of Polyatomic Modelsing the works cited here) perform simulations and track the appearance and lifetimes ofbubbles directly, there are also works using methods such as autocorrelation functions [ ] , which have revealed further properties of base pair openings. As an importantnote with regard to computations with the PBD or similar models, simulations can beperformed either at precise temperatures through Langevin dynamics, e.g. [ ] , or atapproximate temperatures through constant-energy simulations, e.g. [ ] . The use ofLangevin dynamics has the added in fl uence of simulating the behaviour of DNA in a so-lution, imitating the stochastic interaction with neighbouring particles and the effects ofthe larger system on the DNA molecule (see for instance [ ] and references therein).As a result of these works, the relationship between bubble formation and DNAtranscription is strong enough to suggest that the sequence of base pairs alone is notenough to provide the relevant information for transcription. The (randomly occur-ring) thermal openings are important in the transcription process as they provide pointsalong the molecule that are much more likely to attract the RNA polymerase whichinitiates the transcription [ ] . Consequently, there is strong motivation purelyfrom this transcriptional point of view to pursue investigations into the properties ofDNA bubbles, in a variety of transcriptionally signi fi cant sequences.There have also been studies of the distribution of bubble lengths, the probabili-ties of longer or shorter bubbles forming in different sequences [
60, 113 ] . These havefound that in general the likelihood of long bubbles forming in a DNA molecule areextremely small, with the probability of bubble occurrence decreasing with increasingbubble length according to a power law modi fi ed exponential.The signi fi cance of openings in DNA also goes beyond the transcriptional aspect.For instance, in studies of charge propagation along DNA strands, the presence of bub-bles has been found to signi fi cantly impact the transport of the charge [
93, 114–118 ] . Thelocalisation of charge in DNA molecules, also called charge trapping, has been linked tothese fl uctuations in the strand, and in particular regions with a high propensity for bub-ble occurrence have a strong impact on charge trapping [
94, 119 ] . Controlling chargetransport has applications in different branches of bioelectronics, as the use of DNA-likewires for practical engineering considerations is a promising avenue for new develop-ments [ ] .Research has not been restricted to stationary openings in DNA either, with mobilediscrete breathers providing a number of interesting aspects regarding the dynamics andfunctioning of DNA [ ] . These mobile discrete breathers have also been linked to thetrapping of charges in DNA, with recent results suggesting that the presence of thesebreathers can cause targetted movement of charge along the strand [ ] .Additionally, results from dynamical openings in DNA have been used to studymore properties of DNA such as free energy landscapes through modifying the PBDmodel [ ] , and investigating the process of transcription by identifying protein bind-ing sites through their free energy characteristics [ ] . The correspondence of theseresults with experimental fi ndings provides promising signs for the possibility of usefulpredictions about promoters that are currently not well studied. While the study of bubbles in DNA is clearly worthy of attention, when examiningthem from a molecular dynamics perspective there is the obvious problem of choosing94 Malcolm Hillebrand Chapter 5haotic Dynamics of Polyatomic Modelsa threshold for considering a base pair to be in a bubble. Since there is not a particularstandard for this value, there are often wildly varying criteria applied, even within thesame project (for some examples of results using different thresholds see [
50, 107, 113 ] and references therein). Frequently the choice of a particular threshold value is madeto emphasise a certain point – for instance the presence of uncharacteristically large am-plitude bubbles at signi fi cant sites [ ] . In other cases the implemented threshold valuedoes not play a major role in the actual outcomes, affecting the exact numerical detailsbut not the qualitative result [ ] .Thus the proposition of a threshold which can be justi fi ed as a general measure ofwhether a base pair is open or closed would be useful. While in some cases the value willstill be varied to investigate a speci fi c property, a general criterion would enable morestandardised results for cases where the exact threshold value is somewhat arbitrary.Returning to the signi fi cance of the denaturation curves and melting temperatures,a logical criterion is to consider the de fi nition of the melting temperature for DNAstrands, which provides a reference point for open base pairs. Since melting is said tohave occurred when precisely half the base pairs are open, this gives us a relationshipbetween what we want (criterion for an open base pair) and something well established(the melting temperature). The approach now seems fairly straightforward – performsimulations, and retrieve data for the displacements of each base pair, and investigatethis data near the melting temperature. Then by some process, settle on a thresholdwhich gives precisely 50% of base pairs open at the melting temperature. These meltingtemperatures in Kelvin are taken from Kalosakas and Ares [ ] , where they have beencalculated for DNA polymers of various AT and GC constituents and are given in (4.1),as T m = + P GC .In principle the same arbitrary-sequence computations could be applied to particularsequences where the melting temperature is known, if the boundary conditions can becorrectly managed. However, if the aim is to provide a general threshold, there maybe some variance in individual cases, especially in the ePBD model where the sequencedependent stacking interaction provides a signi fi cant difference in melting temperaturesbetween similar molecules. We also note that for consistency and comparability withpreviously published results [ ] , in this chapter the AT / GC content of DNA sequencesis quanti fi ed by the GC percentage P GC rather than the AT percentage that has been usedin prior chapters.For our primary concern of general sequences, let us start with trialling a thresholdvalue of 1Å, which has been used in previous studies as a useful value [
60, 113 ] . In allthese simulations, we have used 100 setups with random initial conditions, and in thecase of heterogeneous sequences, random base pair realisations, with 1000 base pairs inthe sequence. Periodic boundary conditions are used to avoid edge effects. We run thesimulations for 10ns to thermalise, and then record data (displacement of each base pair)every 1ns for the next 10ns. This data can then be analysed using the desired threshold,in this case 1Å, to count the number of base pairs open beyond the threshold. The resultsshown in Fig. 5.1 give the fraction of open base pairs f o , dividing the total number ofrecorded open base pairs by the total number of recorded base pairs.Figure 5.1(a) shows the fraction of open base pairs f o in the DNA molecule as thetemperature T changes for fi ve GC percentages, with a threshold of 1Å for consider-ing the base pair to be open. Note that the temperature is once again computed as theChapter 5 Malcolm Hillebrand 95haotic Dynamics of Polyatomic Models
240 260 280 300 320 340 360 T (K) . . . . . . f o (a) P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0%
240 260 280 300 320 340 360 T (K) . . . . . . f o (b) P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% F IGURE
The fraction of open base pairs f o as a function of temperature T for the PBD model given for fi ve GC percentages: P GC =
0% (blue circles), P GC =
25% (red diamonds), P GC =
50% (green triangles), P GC =
75% (yel-low stars) and P GC = fi cantly as GC content increases, giving f o well below50%. rescaled mean kinetic energy, T = 〈 K n 〉 / k B . We can immediately see that this thresh-old value is far too high to predict half the strand to be separated, as in all cases onlyabout 20% of the base pairs are marked open by the time the melting temperature isreached. So we know that this is a fairly high threshold, and for our current purpose itis not a good choice. It does allow us to understand results that use this threshold bet-ter though – for instance, analysis performed with this (or even larger thresholds) willre fl ect properties of very widely separated bubbles, not considering the contributionsfrom “slightly” open bubbles.Moving past this threshold, in Fig. 5.1(b) we consider a threshold value of 0.2Å.This is a much better choice, giving us approximately 50% of all base pairs open at themelting temperature. It is not perfect though, and displays a clear trend of markingtoo many AT base pairs open and not enough GC base pairs. This is consistent withthe weaker bonding between AT base pairs, resulting in a higher propensity to sepa-rate. This phenomenon has been documented in the study of bubble lengths of differ-ent DNA molecules [ ] , where AT base pairs have been seen to be more likely to openthan GC base pairs. We thus need to consider some sort of heterogeneous approach tothe threshold, in order to account for the different behaviours.Having taken note of the region that this threshold falls in, we can now considerthe properties of our model to see if there is a justi fi cation for any particular choiceof threshold, and particularly any reasonable way of choosing one threshold for GCbase pairs and one for AT base pairs. The natural place to look is at the on-site Morsepotential, which primarily governs the distances between the bases. Referring all the96 Malcolm Hillebrand Chapter 5haotic Dynamics of Polyatomic Models
240 260 280 300 320 340 360 T (K) . . . . f o (a) P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0%
240 260 280 300 320 340 360 T (K) . . . . f o (b) P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% F IGURE
The fraction of open base pairs f o as a function of temperature T , for fi ve GC percentages: P GC =
0% (blue circles), P GC =
25% (reddiamonds), P GC =
50% (green triangles), P GC =
75% (yellow stars) and P GC = y thrAT = y thrGC = a from the Morse potential gives avery accurate representation of 50% of the base pairs open at the meltingtemperature. (b) Using an average threshold in between the two values, y = P GC = way back to (2.4), we recall that the Morse potential has the form V n = D n ( e − a n y n − ) ,where the constants D n and a n depend on the base pair type at the particular site. Thissuggests a natural length scale based on the parameter a n – this provides the characteristicwidth of the Morse potential well, in units of inverse length. Consequently, there is goodreason to investigate the choice of threshold y thrAT = / a AT and y thrGC = / a GC . Calculatingthese quantities, we get y thrAT ≈ y thrGC ≈ fi nd an almost exact presentation of 50% of the base pairs being open at the meltingtemperature. The only small blemish is the pure GC case (purple crosses), where around51% of the base pairs are open, while for all other GC percentages almost exactly 50% ofbase pairs are open. This closeness very strongly motivates for this choice of thresholdas a general measure of when a base pair is open or not.For cases where only a rough measure is needed, a homogeneous averaged thresholdvalue of y thrAT and y thrGC , y = P GC = We are able to immediately make use of our new threshold in the investigation of theproperties of the ePBD model of DNA. Recall that the only difference between the PBDChapter 5 Malcolm Hillebrand 97haotic Dynamics of Polyatomic Models
240 260 280 300 320 340 360 T (K) . . . . . . . f o P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% F IGURE
The fraction f o of open base pairs at a given temperature T forthe ePBD model, when the distinct threshold values y thrGC = / ≈ y thrAT = / ≈ fi ne the opening of a GC and anAT base pair respectively. The solid line indicates 50% open base pairs.We see that using these thresholds, we can obtain approximate meltingtemperatures of T m = ( + P GC ) K, which retains the experimentallyobserved linear relationship. The dotted lines show the proposed meltingtemperatures for each GC percentage. and ePBD models is the sequence dependence built into the ePBD model (given by thevalues shown in Table 2.1), by varying the coupling strength depending on the exactordering of the bases in the strand [ ] . In light of this, and particularly due to the factthat the Morse potential [ from (2.2) ] remains unchanged between the two models, withthe parameters chosen to reproduce physical results, the threshold we have found for thePBD model can be used for the ePBD model in exactly the same way.In order to fi nd a general formula for melting temperatures of DNA molecules withthe PBD model, a fairly complex procedure had to be followed involving fi ttings of bub-ble pro fi les [ ] . The advantage of being able to carry the threshold values over fromthe PBD model circumvents this lengthy process in the ePBD model, as we can in factdirectly fi nd the melting temperatures by reversing the process applied above – insteadof fi nding the threshold to match known melting temperatures, we can use the knownthreshold to locate the correct melting temperature. Figure 5.3 shows the result of fol-lowing this approach by plotting the fraction of open base pairs f o against the tempera-ture T for the ePBD model. For all fi ve GC percentages shown, we see that the thresh-olds predict very clear melting temperatures which follow the known linear relationshipwith the GC content P GC [ ] . The obtained melting temperatures T M are shown bythe dashed vertical lines in Fig. 5.3, with the horizontal solid black line marking 50%open base pairs. The melting temperature T M of DNA sequences as a function of theGC content P GC in the ePBD model is thus given by T M = ( + P GC ) K. (5.1)The precise matching of this relation to the numerically observed melting pointsin Fig. 5.3 is convincingly close, and we are able to use this equation to estimate themelting temperatures of DNA molecules of various GC percentages. Since the ePBDmodel is designed to be more sensitive to the exact base pair sequencing, there is a lot98 Malcolm Hillebrand Chapter 5haotic Dynamics of Polyatomic Models T (K) . . . . . . . . E n ( e V ) (a) P GC = 100% P GC = 50% P GC = 0% 0 50 100 150 200 250 300 350 T (K) . . . . . . . . E n ( e V ) (b) P GC = 100% P GC = 50% P GC = 0% P GC (%) . . . . . γ ( e V / K ) × − PBD: γ = [4 . . P GC ] × − ePBD: γ = [5 . . P GC ] × − (c) PBDePBD F IGURE
The energy density–temperature relationship for (a) the PBDmodel, and (b) the ePBD model for three GC percentages, P GC =
0% (bluecircles), P GC =
50% (green triangles), and P GC = fi tting to (5.2) foreach GC percentage shown. (c) The variation of the parameter γ of (5.2)with GC content, for the PBD model (blue circles) and the ePBD model(empty orange squares). The fi tted linear equations are marked by solid(PBD) and dashed (ePBD) lines, with the equations displayed in the panel.The error of the last signi fi cant digit of the obtained parameters is shown inparentheses after each fi tted value, e.g. 4.97 ( ) corresponds to 4.97 ± of variability for every GC percentage. For instance, pure GC sequences are known tohave melting temperatures varying by more than 25K, and this variance is captured bythe ePBD model [ ] . So we should be careful when applying this relation to speci fi ccases, but for all of our investigations here, we are concerned with average behavioursof many realisations, meaning that this relation is certainly applicable. With the melting temperature (5.1) for the ePBD model, we can brie fl y interrogate therelationship between energy density E n and temperature T for the PBD and ePBD mod-els, now that we know the relevant bounds for the ePBD model. To this end, we per-formed simulations at several energy density values and computed the average temper-ature across these runs. For each GC percentage, at each energy density, 100 differentrandom initial conditions were run, also with random base pair realisations in hetero-geneous cases, using periodic sequences of 1000 base pairs. Each simulation was ther-malised for 10ns, then data recorded for the next 1ns every 0.01ns.The results of these simulations are shown in Fig. 5.4(a), for the PBD model and inFig. 5.4(b) for the ePBD model. Only three GC percentages are shown, so as not toobscure the results, P GC =
0% (blue circles), P GC =
50% (green triangles), and P GC = Chapter 5 Malcolm Hillebrand 99haotic Dynamics of Polyatomic Models100% (purple crosses). GC-rich sequences [ such as the pure GC case, purple crosses inFig. 5.4(a) ] tend to have a higher energy density at a given temperature, particularly athigher temperatures. In all cases we see that at low energies, there is an evident linearrelationship between energy density and temperature, with the GC content making verylittle difference to the results. As the temperature increases however, a nonlinear trendemerges, and for temperatures above about 150K, the linear approximation is certainlyinsuf fi cient. In order to model this change, we have fi tted the data with a linear termproportional to the Boltzmann constant k B , which is accurate for low energies, but alsoadded a cubic term to more accurately follow the higher temperature behaviour, E n = k B T + γ T . (5.2)The constant γ is free to be fi tted for the two models.The fi tted curves are shown in Fig. 5.4(a) and (b) by solid and dashed lines respec-tively. We see that in all cases (5.2) provides a very good fi t to the data, despite the sim-plicity of the functional form. The values of the fi tted constant γ are shown in Fig. 5.4(c)as the GC content is varied, for both models. The linear increase of γ with the GC per-centage allows the fi tting of a straight line to the data, resulting in a simple relation whichis included in the fi gure. There is an overall tendency for γ to be larger in the ePBDmodel, although the changes with GC content are very similar for the two models.Equation (5.2) then gives us a very useful way of referring back and forth between E n and T for the two models, combining (5.2) with the values of γ given in Fig. 5.4(c)even providing us with a direct way of computing E n for a given T and GC percentage.For investigations in the microcanonical ensemble, working with constant energy ratherthan temperature, it can frequently be useful to have a link between the two quantities,as often T is a more physically useful quantity to consider. Of particular importance toour work here, where we will look largely at physiological temperatures, is that we canclearly tell from the outset that for different GC percentages we need to use different E n values in order to achieve the same T values. So in order to compare results at the sametemperatures, we can make use of this relationship (5.2) to fi nd the appropriate energyvalue corresponding to T = With these initial results all cleared, we can investigate some properties of bubbles inDNA sequences using our heterogeneous threshold. The fi rst aspect we will look at isthe probability distribution for bubbles of various lengths to occur in a DNA strand at T = [ ] ,as well as for temperatures in the range of 270-350K [ ] . These investigation made useof Monte-Carlo methods, and used a single threshold value of 1Å (which, as we saw inFig. 5.1, does not necessarily capture small bubbles correctly), rather than our use ofmolecular dynamics and heterogeneous thresholds here.For this study, we performed 8 000 simulations of DNA sequences with N = − − − − − P L ( l ) (a) P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% l (bp) − − − − − P L ( l ) (b) P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% . . β l (c) PBD: β l = 0 . . P GC ePBD: β l = 0 . . P GC PBDePBD0 . . κ ( bp ) (d) PBD: κ = 0 . . P GC ePBD: κ = 0 . . P GC P GC (%) C (e) PBD: C = 0 . − . P GC ) + 1 . C = 2 . − . P GC ) + 1 . F IGURE
The distribution P L of bubble lengths in the PBD (a) and ePBD(b) models for fi ve GC percentages, namely P GC =
0% (blue circles), P GC =
25% (red diamonds), P GC =
50% (green triangles), P GC = P GC = fi ttings to (5.3). (c)–(e):The fi tting parameters of (5.3) to the results of (a) and (b), the stretchingparameter β l , the characteristic length κ , and the prefactor C respectivelyas a function of the GC content P GC . The PBD parameters are shown bysolid blue circles, and the ePBD parameters by empty orange squares. For β l and κ , the variation of the parameters with GC content has been fi ttedwith a straight line, and the equation of the lines for the two models aredisplayed in each panel. The constant C is fi tted with a simple exponentialfunction, and this equation is also presented in the fi gure. recorded for a further 10ns, taken every 0.1ns. At each recording step, for each base pairthe information was stored as to whether it was open or closed, based on the thresh-old related to its base pair type. Thus, for each run data was recorded for a total of 1000 base pairs at 100 steps, resulting in 100 000 data points (open / closed base pair) perrun, which were analysed to fi nd the lengths of all bubbles occurring. Combining allthese runs together, we had a total of 800 000 000 recorded data points, providing strongstatistics even up to the extreme ends of the probability distributions. Analysis of thefull data set was compared with using only sections of the data to check that our resultsare not affected by a lack of data, and we found that our results were virtually identicaleven when only using 10% of the total data. This gives us a good deal of con fi dence inthe data, and allows us to trust the tails of the distribution more strongly.The results of this analysis are shown in Fig. 5.5(a) and (b), showing the probabilitydistribution P L ( l ) of bubble lengths for the PBD and ePBD models respectively. Someinitial points are immediately apparent – for instance, the heterogeneity of the sequencesappears to be counterbalanced by the base-pair dependent threshold for short bubbles,as bubbles with l < fi cant variation in probabilitywith GC content. However, as the lengths increase, the base pair composition becomesChapter 5 Malcolm Hillebrand 101haotic Dynamics of Polyatomic Modelssigni fi cant once again, as sequences with fewer GC base pairs demonstrate a markedlyhigher probability for long bubbles ( l > [ purple crosses in Figs. 5.5(a) and (b) ] and pure AT (blue dia-monds) cases are the least and most likely to exhibit long bubbles, with the intermediatepercentages showing a continuous trend of more long bubbles being present when theGC content is lower.These overall behaviours are observed in both the PBD model [ Fig. 5.5(a) ] and theePBD model [ Fig. 5.5(b) ] , but there are some small differences between the two. Pri-marily, the ePBD model exhibits a larger probability for bubbles to appear at almost alllengths, and particularly for long bubbles. This effect is strongest in the AT-dominantstrands (e.g. the P GC =
25% case, red circles in Figs. 5.5(a) and (b), and the pure AT case,blue diamonds), where the tail of the distribution stretches noticeably longer than inthe PBD model. The pure GC case [ purple crosses in Figs. 5.5(a) and (b) ] appears to beaffected very little by the choice of model, with no discernible differences between thetwo plots. This behaviour results in a much more signi fi cant spread of results for theePBD model, with the distributions distinctly more separated than in the PBD case.We also note that these general trends have been seen previously, using a uniformthreshold value [ ] , where particularly for large bubbles the behaviour is fairly similar.One distinct point of difference however is the behaviour for short bubbles, where in ourcase the base pair dependent threshold results in the probabilities of bubble occurrencebeing almost identical for all GC percentages.These obtained P L ( l ) distributions can be fi tted fairly accurately by a stretched ex-ponential, of the form P L ( l ) = C exp (cid:8) − ( l / κ ) β l (cid:9) , (5.3)with C , κ and β l parameters free to be fi tted. This follows the data accurately even upto the tails in Figs. 5.5(a) and (b). These fi ts are shown for the PBD model in Fig. 5.5(a)as solid curves, and in Fig. 5.5(b) by dashed curves for the ePBD model, con fi rming thatthey provide a very good approximation of the numerical data.The dependence of the fi tting parameters of (5.3) on the GC content P GC is illustratedin Fig. 5.5(c)–(e). In Fig. 5.5(c) we see a linear increase in the stretching exponent β l as theGC content is increased. This holds for both the PBD model [ blue circles in Fig. 5.5(c) ] and the ePBD model [ empty orange squares in Fig. 5.5(c) ] , albeit with different relationsfor the two models. The linear increase has in turn been fi tted with a straight line, andthe equations of these lines are given in the panel. For the characteristic length κ , wealso see a linear increase with the GC percentage, in Fig. 5.5(d). This linear functionhas also been fi tted, and the equations are shown in the panel. For the prefactor C , thechange with GC content is no longer linear, but rather exponential [ Fig. 5.5(e) ] . Theexponential decrease can be fi tted accurately with a simple exponential function for bothmodels, as shown in Fig. 5.5.In general, the values of the fi tting parameters β l , κ and C con fi rm the impressionfrom Figs. 5.5(a) and (b) that the difference between the PBD and ePBD models is morepronounced for lower GC content. For instance, the stretching exponent β l [ Fig. 5.5(c) ] and the prefactor C [ Fig. 5.5(e) ] both are almost identical for pure GC sequences, while κ [ Fig. 5.5(d) ] is very close. At the other end of the spectrum though, pure AT sequences( P GC = P GC . . . . . . . h l i ( bp ) PBD: h l i = 1 . − . P GC ) + 0 . h l i = 1 . − . P GC ) + 1 . F IGURE
The average bubble length 〈 l 〉 (5.4) as a function of GC percent-age P GC , for the PBD model (blue circles) and the ePBD model (emptyorange squares). The data are fi tted with an exponential function, whoseequations is shown in the fi gure. The solid blue and dashed orange curvesshow these fi ts for the PBD and ePBD models respectively. An additional quantity we can calculate from the bubble length distribution P L ( l ) is the average bubble length 〈 l 〉 , which can be calculated according to the equation (af-ter [ ] ) 〈 l 〉 = (cid:5) l l P L ( l ) (cid:5) l P L ( l ) . (5.4)Here P L ( l ) is the probability of a bubble of length l occurring [ Figs. 5.5(a) and (b) ] .Figure 5.6 shows the computed 〈 l 〉 values for both models. We see a clear decrease inthe average bubble length 〈 l 〉 as the GC content increases (also previously observed fora constant threshold in [ ] ), for both models. The ePBD model (orange squares inFig. 5.6) shows generally longer average bubble lengths, although once again the twomodels more or less converge for pure GC sequences. A large range of 〈 l 〉 values is alsoseen with the ePBD model, as a result of the greater sensitivity inherent in the model.The variation of the average bubble length with GC content can be fi tted accuratelywith an exponential function, where the solid blue (dashed orange) curves show this fi tfor the PBD (ePBD) model in Fig. 5.6. The equations of the curves are embedded in the fi gure. Another, less studied, aspect of denaturation bubbles is their lifetime. In a typical DNAdouble strand, how long would one expect a bubble to persist? In order to address thisquestion, we can once again perform MD simulations, and by using detailed records ofthe base pair displacements, fi nd distributions of the bubble lifetimes. The aim here is to fi nd (at physiological temperature) the probability of a particular bubble lifetime givena) the GC percentage P GC of the sequence, and b) the length l of the bubble. We alreadyknow from Section 5.3 that long bubbles form extremely rarely, but here we will go onestep further, and see what we can learn about the lifetimes of these bubbles.Firstly, let us outline the algorithm used to track and record bubbles, as this is afairly non-trivial task. We note that this is much more complex than the analysis toChapter 5 Malcolm Hillebrand 103haotic Dynamics of Polyatomic Modelsproduce the bubble length distributions, which does not require the tracking of bubblesthrough time. For a single simulation (one run of a single sequence with random initialcondition), the following steps are followed to create the distribution of bubble lifetimes1. For the given initial condition, perform the MD simulation to create, at each timestep, a record of whether the base pair is open or closed, for each base pair in thesequence. In our simulations, the open / closed data is de fi ned using the thresholdsde fi ned in Section 5.2.2. For every time step, move through the sequence and record, at the start site ofeach bubble, the length of each bubble that occurs.3. Check every recorded bubble (site and length) against the record from the previoustime step. • If there is a bubble starting where previously there was no bubble recordedat that site, begin a record of that bubble – a tuple of (length, lifetime). • If a bubble starting at this site survives with the same length, then incrementits lifetime by one time step. • If a bubble starting at this site survives but changes length, then end therecord of that bubble, and begin a new tuple at that site with the new length. • If there is no bubble present where one was at the previous time step, thenend the record and do not create a new one. • If there is no bubble present, and no bubble was previously present, then donothing.4. Finally, record the list of tuples of (length, lifetime) at each site.5. Create the distribution of these (length, lifetime) data points, for each length.This procedure allows us to study the bubble lifetime distribution for bubbles ofvarying lengths, and effectively combine many different runs to create a coherent pic-ture of the general behaviour. As in all our investigations, we perform this analysis basedon the GC content P GC of the sequences, by using many runs of random base pair real-isations with the same numbers of GC and AT base pairs.So we can now begin the investigation by creating bubble lifetime distributions forbubbles of varying lengths. For each GC percentage P GC studied, we have used 1 000different realisations with random initial conditions, integrating a sequence of 100 basepairs with periodic boundary conditions. A thermalisation time of 10ns was again usedto ensure equilibration of the system, and then data were recorded every 0.01ps for thenext nanosecond. This data recording whether each base pair is open or closed wasanalysed using the algorithm described above, and we created normalised distributionsof the observed bubble lifetimes, for lengths of up to l =
10 base pairs. We note that in allcases we have more than 200 000 recorded bubbles (with shorter bubbles having severalmillion records), which gives a minimal likelihood for statistical anomalies to occur.Results could be presented for longer bubbles, but as noted in Section 5.3, these longbubbles are exponentially rare, and consequently the amount of data for these bubblesis signi fi cantly lower than for smaller bubbles.104 Malcolm Hillebrand Chapter 5haotic Dynamics of Polyatomic Models P l ( t ) l = 1 l = 6 P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% P l ( t ) l = 2 l = 7 P l ( t ) l = 3 l = 8 P l ( t ) l = 4 l = 9 . . . . . . t (ps) P l ( t ) l = 5 . . . . . . t (ps) l = 10 F IGURE
Bubble lifetime distributions P l ( t ) in the PBD model for fi veGC percentages, P GC =
0% (blue circles), P GC =
25% (red diamonds), P GC =
50% (green triangles), P GC =
75% (yellow stars) and P GC = l = fi ttings to the stretched exponential (5.5). Chapter 5 Malcolm Hillebrand 105haotic Dynamics of Polyatomic Models P l ( t ) l = 1 l = 6 P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% P l ( t ) l = 2 l = 7 P l ( t ) l = 3 l = 8 P l ( t ) l = 4 l = 9 . . . . . . t (ps) P l ( t ) l = 5 . . . . . . t (ps) l = 10 F IGURE
Similar to Fig. 5.7, but showing bubble lifetime distributions P l ( t ) for the ePBD model, for the same fi ve GC percentages, P GC = P GC =
25% (red diamonds), P GC =
50% (green triangles), P GC =
75% (yellow stars) and P GC = fi tted, shown by dashed curves.
106 Malcolm Hillebrand Chapter 5haotic Dynamics of Polyatomic ModelsThe resulting bubble lifetime distributions P l ( t ) are shown in Figs. 5.7 and 5.8, forthe PBD and ePBD model respectively. In each panel of the fi gure the distributions areshown for a particular bubble length l (as denoted by the label on the panel), for fi ve GCpercentages, P GC =
0% (blue circles), P GC =
25% (red diamonds), P GC =
50% (greentriangles), P GC =
75% (yellow stars) and P GC = l = fi rst panel of Figs. 5.7 and 5.8) –there appears to be a kind of double-peaked structure, with an initial peak around 0.1psand a second following later at either 0.025ps (for GC-rich sequences) or 0.45ps (forAT-rich sequences). The peaks are also much broader in the pure AT sequences [ bluediamonds in Fig. 5.7 ] . This structure may be the result of some interplay between thecharacteristic times of l = l increases,this double peak feature does remain, although less clearly. For l = fi rst peak is dramatically less pronounced, while the second peak is still visible but alsomuch smaller in height. Past l =
3, we have much smoother distributions, with thepeaks becoming less evident. The distributions also demonstrate a fairly clear tendencyto exhibit fewer and fewer long-lived bubbles, although the last few cases l =
7, 8, 9, 10in Fig. 5.7 are extremely similar visually.The differences between the PBD results in Fig. 5.7 and the ePBD results displayedin Fig. 5.8 are not particularly strong. We fi nd the same behaviour for single base pairopenings ( l = P l ( t ) = A exp (cid:13) − ( t / τ ) β (cid:14) , (5.5)having as fi tting parameters A , τ , and β . This provides a very good approximation ofthe true distribution in all cases apart from the anomalous l = fi tting ofthe distributions to this function are shown in Figs. 5.7 and 5.8 by solid and dashed linesrespectively. In both these fi gures we see that for the l = l = fi ts very closely, and gives us a veryuseful tool to understand the general likelihood of bubbles of particular lengths andlifetimes occurring.The fi tting parameters for the stretched exponential functions (5.5) are shown inFig. 5.9 for the PBD and ePBD models. First considering the stretching exponent β [ Fig. 5.9(a) ] , we see that as bubble lengths increase β rapidly approaches 1, correspond-ing to a standard exponential function in (5.5). This con fi rms the early impression fromthe distributions themselves [ Figs. 5.7 and 5.8 ] that for long bubbles we had a consistent,exponential, behaviour. The very large value of β for l = P l ( t ) distribution for l = β beyond the fi rst two points, i.e. for l ≥
3, and that the two modelshave extremely close values for β in all cases.Chapter 5 Malcolm Hillebrand 107haotic Dynamics of Polyatomic Models l (bp) . . . . . . β (a) P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% l (bp) . . . . . . . τ ( p s ) (b) P GC = 100% P GC = 75% P GC = 50% P GC = 25% P GC = 0% . . . c β PBDePBD P GC (%) . . λ β ( bp ) . . c τ ( p s ) P GC (%) . . τ ( p s ) F IGURE
The parameters of (5.5) fi tted to the bubble lifetime distributionsof Figs. 5.7 and 5.8. (a) The stretching factor β as a function of the bubblelength l for fi ve GC percentages, P GC =
0% (blue circles), P GC = P GC =
50% (green triangles), P GC =
75% (yellow stars)and P GC = β with increasing bubble length is fi tted with a simple exponential, (5.6),which is shown by solid and dashed curves for the PBD and ePBD modelsrespectively. The corresponding parameters c β and λ β of the fi tting areshown by the insets as functions of the P GC values, with a quadratic fi tting(solid blue curve for PBD, dashed orange for ePBD). (b) The characteristictime τ , as a function of the bubble length l . Once again an exponentialdecrease is evident for both models, fi tted with (5.7). The insets show thevariation of the parameters c τ and τ of (5.5) with P GC , this time fi tted bystraight lines. The characteristic time τ [ Fig. 5.9(b) ] also exhibits an exponential decrease with bub-ble length l , although nothing so drastic as for β . Here we see a more steady decreasetowards a constant asymptotic value, which is dependent on both the model and the GCpercentage P GC . The ePBD model (depicted by empty symbols in Fig. 5.9) has a con-sistent tendency to produce a larger τ value for all bubble lengths than the PBD model.The effect of the model is apparently stronger for AT-rich sequences, with the differencebetween the PBD and ePBD value of τ being most signi fi cant for the pure AT case [ bluediamonds in Fig. 5.9(b) ] .Since both these quantities vary exponentially with l , we are able to perform a fi ttingwith simple exponential functions of the form β = c β exp ( − l / λ β ) +
1, (5.6) τ = c τ exp ( − l / λ τ ) + τ . (5.7)These fi ts are shown in Figs. 5.9(a) and (b) by solid curves for the PBD model and dashedcurves for the ePBD model, and the insets of these fi gures show the dependence of theparameters c β , λ β , c τ , and τ on the GC percentage. First addressing β , we see fromthe inset of Fig. 5.9(a) that the values of c β and λ β are very similar for the two mod-els, with no consistent difference between the two. Furthermore, we can perform the108 Malcolm Hillebrand Chapter 5haotic Dynamics of Polyatomic Modelssimplest fi tting that captures the correct behaviour of these values, which in this caseis a quadratic function. These fi ttings are shown by the solid blue and dashed orangecurves in the inset of Fig. 5.9(a). The fact that the fi tted curves are almost identicalalso suggests that the effect of the model on the value of β is minimal, and we canuse the same values of c β and λ β for both cases while preserving accuracy. The com-puted equations for the parameters are c β = ( )( P GC ) − ( ) P GC + ( ) , and λ β = − ( )( P GC ) + ( ) P GC + ( ) .For the fi ttings to the values of τ [ (5.7) and Fig. 5.9(b) ] , we found that a value of λ τ = c τ and τ of (5.5) both decrease linearly withincreasing GC percentage P GC . Fitting these values with straight lines, we see that for c τ the fi tted lines are effectively identical for both models, while there is a distinct differencebetween the models in the values of τ [ see inset of Fig. 5.9(b) ] . The computed equationsfor the parameters are given by c τ = ( ) − ( ) P GC , and then for the PBD model τ = ( ) − ( ) P GC and for the ePBD model τ = ( ) − ( ) P GC .Finally, the prefactor A in (5.5) is considered to be free in our fi tting procedure, butit is fundamentally de fi ned by the normalisation of the distribution, which it maintainsto an accuracy of around 5%. Note also that the lines of best fi t shown in Figs. 5.9(a)and (b) are plotted using the parameter values estimated using the quadratic and linear fi ttings (shown in the insets of the fi gure), con fi rming that these fi tted equations stillproduce an accurate representation of the true situation.With these parameters for (5.5), we are now able to predict the probability of a bub-ble occurring with a particular length and lifetime, given only the GC content of theDNA sequence (and of course choosing a model). One point that is worth emphasisinghere is how extremely rare long-lived large bubbles are in the general DNA sequences,regardless of the AT / GC composition. The combination of the near-exponential rar-ity of long bubbles in general [ see Figs. 5.5(a) and (b) ] with the exponential decrease inprobability of long-lived bubbles [ Figs. 5.7 and 5.8 ] means that the chances of fi ndingsuch a bubble by pure chance are “doubly exponentially small”. This provides a usefulbackdrop to studies performed on bubble pro fi les in biological promoters, highlightingthe fact that the presence of large long-lived bubbles at transcriptionally active sites isdistinctly anomalous (see e.g. [
50, 54 ] )A further property of the bubble lifetime distributions that we can investigate is theaverage lifetime; exploring how long the average bubble lives at different lengths and instrands of varying GC content. We can estimate the mean bubble lifetime 〈 t 〉 l directlyfrom the distributions, by using the formula 〈 t 〉 l = N b (cid:1) j = t j P l ( t j ) δ t , (5.8)where here the sum is over all the N b =
500 bins of the histogrammed lifetime data, withwidth of δ t = t j is the timeat the centre of the j th bin, and P l ( t j ) is the numerical probability of a bubble lifetimeoccurring within that bin from Figs. 5.7 and 5.8.The result of performing this calculation for each bubble length is shown in Fig. 5.10for a range of GC percentages, with the average bubble lifetime for the PBD modelChapter 5 Malcolm Hillebrand 109haotic Dynamics of Polyatomic Models . . . . . h t i l ( p s ) (a) P GC = 100% P GC = 88% P GC = 75% P GC = 63% P GC = 50% P GC = 37% P GC = 25% P GC = 12% P GC = 0% l (bp) . . . . . h t i l ( p s ) (b) P GC = 100% P GC = 88% P GC = 75% P GC = 63% P GC = 50% P GC = 37% P GC = 25% P GC = 12% P GC = 0% . . . . δ ( p s ) (c) PBD: δ = 0 . − . P GC ePBD: δ = 0 . − . P GC PBDePBD1 . . . . α ( bp ) (d) PBD: α = 1 . − . P GC ePBD: α = 1 . − . P GC P GC (%) . . . . B ( p s ) (e) PBD: B = 0 . . P GC ePBD: B = 0 . . P GC F IGURE
The numerical mean bubble lifetimes 〈 t 〉 l for various GC per-centage values P GC (5.8) as a function of the bubble length l , for (a) thePBD model, and (b) the ePBD model. Larger GC percentages correspondto shorter mean lifetimes, and shorter bubbles exhibit longer lifetimes.The decrease in mean lifetime with increasing bubble length is fi tted by asimple exponential (5.9). These functions are shown by solid curves in (a)and dashed curves in (b). The parameters of this fi t are shown in panels(c)–(e) as a function of the P GC value, the asymptotic limiting value δ , thecharacteristic length α , and the prefactor B respectively. Results for thePBD model are shown by solid blue circles, and for the ePBD model byempty orange squares. All of these quantities vary more or less linearlywith the GC percentage, and the fi ts are shown by solid blue lines (PBD)and dashed orange lines (ePBD), with the equations of the lines shownwith respect to the bubble length l in each panel. shown in Fig. 5.10(a), while the ePBD results are seen in Fig. 5.10(b). Similarly to thebehaviour of τ [ (5.7) and Fig. 5.9 ] , we see an exponential decrease in the mean bub-ble lifetime 〈 t 〉 l with increasing bubble length l . Sequences with more AT base pairsclearly exhibit a longer average lifetime, with a monotonic increase in 〈 t 〉 l as the GCcontent P GC decreases at every bubble length. In addition, comparing the two models [ Figs. 5.10(a) and (b) ] we see that in general the ePBD model predicts longer averagelifetimes, with a difference of around 0.02ps for each GC percentage. The differencebetween the two models is largest for the pure AT case [ blue diamonds in Figs. 5.10(a)and (b) ] , and smallest in the pure GC case (purple crosses).The exponential decrease in mean lifetime 〈 t 〉 l can be fi tted accurately by the simplefunctional form 〈 t 〉 l = B exp ( − l / α ) + δ , (5.9)with B , α and δ free fi tted parameters. The fi ts of the data to this equation are shownby solid curves in Fig. 5.10(a) and dashed curves in Fig. 5.10(b), verifying that this formis an appropriate function to model the behaviour.110 Malcolm Hillebrand Chapter 5haotic Dynamics of Polyatomic ModelsFigures 5.10(c)–(e) give the variation of the fi tting parameters with GC percentage,the limiting constant δ in Fig. 5.10(c), the characteristic length α in Fig. 5.10(d), and thenthe prefactor B in Fig. 5.10(e). In all cases we see a more or less linear change with GCpercentage. For δ [ Fig. 5.10(c) ] , we have a very clear linear decrease as the GC percent-age increases, for both the PBD model (blue circles) and the ePBD model (empty orangesquares). The ePBD model gives a consistently larger limiting value, which con fi rms ourobservation of the differences between the results in Figs. 5.10(a) and (b). The equationof the line of best fi t is given in the panel. The characteristic length α [ Fig. 5.10(d) ] showsa less clearly linear behaviour, but nevertheless using a straight line fi t gives a reasonableapproximation. This time the PBD model exhibits higher values, with the slope onceagain fairly similar for the two models. Finally the prefactor B [ Fig. 5.10(e) ] is also fi ttedwith a linear function (the equation is also given in the panel), with the ePBD modelgiving higher values, as well as a steeper slope. The values of B for small values of P GC are fairly similar, but the discrepancy is much larger for GC rich sequences. A natural step after fi nding these general bubble lifetime distributions and properties isto compare them with results from particular biological promoters, in the same way aswe explored whether the mLE is different for these promoters in Section 4.3.1. Here wewill focus on the AdMLP, since it has been the focus of many previous studies (startingfrom [ ] ), and there is an established theme of large, long-lived bubbles forming attranscriptionally active sites [ ] . Consequently, it will be instructive to see if there aredifferences in the bubble lifetime distributions inherent to the promoter sequence, aswell as the average lifetimes.A direct comparison of the lifetime distributions P l ( t ) are presented for the ePBDmodel in Fig. 5.11, showing the distributions for the AdMLP (empty circles in Fig. 5.11),as well as for sequences with exactly the same base pair composition (corresponding to P GC = / ≈ l = l >
7, the green triangles of the AdMLP are consistently abovethe red diamonds of the random sequences. These are small differences, but as the un-certainties of the mean value computation are very small (less than 1%), this adds moresigni fi cance to the slightly different behaviours of the bubble lifetime distributions inFig. 5.11. It is evident that while the observed trend was minor, the consistency is suf fi -Chapter 5 Malcolm Hillebrand 111haotic Dynamics of Polyatomic Models P l ( t ) l = 1 AdMLP P GC = 66% l = 3 AdMLP P GC = 66% . . . . . . t (ps) P l ( t ) l = 2 AdMLP P GC = 66% . . . . . . t (ps) l = 8 AdMLP P GC = 66% F IGURE
A comparison between the bubble lifetime distributions P l ( t ) of the AdMLP (shown by blue empty circles) and randomly ordered se-quences comprised of the same base pairs (shown by the red solid line).Results are shown for four bubble lengths, l = l = l = l = cient to affect the mean lifetime noticeably. It is also consistent with the previous studiesthat the AdMLP exhibits longer-lived bubbles, even if our results are not focussing onany particular regions of the promoter, where anomalously long-lived large bubbles havebeen found [ ] .These results are however very encouraging, as they suggest that the average distribu-tions we have found can be meaningfully applied to speci fi c sequences, without fear thatin particular cases our results are not applicable. This allows the general results of thissection to be used as a benchmark for the expected bubble lifetime behaviour of speci fi csequences, meaning that investigations into the transcriptionally signi fi cant regions ofthe sequence can be con fi dently compared to the general distribution in order to quan-tify how unusual the observed behaviours are. Using these bubble lifetime distributionsis a promising future line of research to understand anomalous behaviour of functionalDNA regions. In this chapter we have studied properties of bubbles in DNA strands by introducingbase pair dependent threshold values based on the characteristic length of the Morsepotential, that satis fi es the de fi nition of the melting transition occurring when 50% ofbase pairs in the strand are open. Using these thresholds, we have investigated the energy-temperature relationship in DNA, distributions of bubble lengths and distributions ofbubble lifetimes through extensive numerical simulations. In all aspects, both the PBDmodel and the sequence-dependent ePBD model were studied.112 Malcolm Hillebrand Chapter 5haotic Dynamics of Polyatomic Models l (bp) . . . . . . . h t i l ( p s ) PBD AdMLPPBD P GC = 66%ePBD AdMLPePBD P GC = 66% F IGURE
The mean bubble lifetime 〈 t 〉 l (5.8) for the AdMLP and ran-domly generated sequences with the same base pairs in random arrange-ments. The results are shown for the PBD model (AdMLP: purple crosses,random sequences: yellow stars) and the ePBD model (AdMLP: green tri-angles, random sequences: red diamonds). Using the derived threshold values, an estimation for the average melting tempera-ture of the ePBD model was presented, which fi ts well with experimental understanding.A nonlinear relationship (5.2) was found between the temperature T (calculated throughthe averaged kinetic energy) and the energy density E n of the system, where a cubic cor-rection to the low-temperature relation E n = K B T produces an accurate approximation.Similar results are found for the PBD and ePBD models.The distributions of the lengths of bubbles in DNA was found (Fig. 5.5), which areapproximated well by stretched exponential functions (5.3). Due to the considered het-erogeneous thresholds, the probability for short bubbles to occur is similar for any basepair composition, but for long bubbles ( l > fi ttedby stretched exponentials (5.5), although for long bubbles ( l (cid:3) [ ] could be probed by detailed examination of bubbleprobabilities and lifetimes at lower temperatures.114 Malcolm Hillebrand Chapter 5 hapter 6Dynamics of Graphene This chapter represents a slight change of tack, moving away from the DNA models stud-ied in the previous chapters, and presenting an investigation into the chaotic dynamicsof a planar model of graphene. First we address the necessary background of graphene asa material as well as the computational machinery required to study a complex 2D sys-tem numerically, including scaling ef fi ciencies for CPU and GPU parallel integration.Thereafter we look at results from the computation of the mLE for an “in fi nite” sheetof graphene, modelled as a graphene sheet with periodic boundary conditions to imitatethe bulk behaviour of the material. Further, we will compute the mLE for fi nite widthgraphene nanoribbons (GNRs) with both zigzag and armchair edges. A short discus-sion of various simpli fi cations of the model is also presented, brie fl y investigating thein fl uence of the 2D hexagonal geometry on the chaoticity of the system.The results in this chapter are based around the research presented in [ ] . Since bursting onto the scene in 2004 [ ] , graphene has attracted a spectacular amountof attention due to its remarkable properties. Graphene in itself is a surprisingly simplematerial – just a single 2D layer of carbon atoms connected in a hexagonal lattice, co-valently bonded together. As far as the bonding details go, each carbon atom forms aninteratomic σ bond (the strongest covalent bond, formed by the end-to-end overlappingof atomic orbitals) with each of its three neighbouring atoms, as well as a single π bond (aweaker covalent bond, formed by lateral overlap of atomic orbitals) perpendicular to theplane [ ] . The σ bonds are the primary cause of the exceptional structural strength ofgraphene (and in fact the strength of this honeycomb structure in general) [ ] , whilethe π bonds are the primary contributors to the electronic properties of graphene thatpresent it as a candidate for so many technological innovations [ ] .Among these applications are high energy supercapacitors [ ] , ef fi cient transis-tors [ ] , electromechanical resonators [ ] as well as advanced sensors [ ] . Theseengineering applications are possible in large part due to the extraordinarily high elec-tron mobility [ ] and the similarly exceptional thermal conductivity of graphene [ ] . Further applications of graphene to structural engineering are highly promisingdue to its extremely high intrinsic strength, with a tensile strength measured at 130 GPa115haotic Dynamics of Polyatomic Modelsand a Young’s modulus of 1TPa [ ] . This strength translates into very high criticalbuckling strains under tension and compression [ ] , meaning that graphene has po-tential for extensive uses in structure development and reinforcement.The dynamical aspects of graphene are of particular relevance to this work, andindeed MD simulations have been used extensively to investigate its properties. Forinstance, MD simulations have been used to study the interface interactions betweengraphene and different metals, studying the cohesive energy and bonding strength fromdynamical considerations, which provides valuable information for the development ofengineering applications [ ] . MD simulations have also been useful in the further anal-ysis of the thermal conductivity of graphene [ ] , and in particular to study the effectsof the size of the sheet considered. A signi fi cant advantage of these numerical techniquesis the ability to consider differently sized GNRs, thin strips of graphene, without hav-ing to painstakingly create an accurate sample in a laboratory. These investigations haverevealed that the ribbon length plays a major role in the thermal conductivity, and alsothat applying strain to the ribbon can decrease the conductivity [ ] . In addition to theconventional 2D forms of graphene, the material has further scope for physical applica-tions in various conformations through controllable defect engineering, such as carbonnanotubes, extremely strong cylindrical forms of graphene [ ] , nanocages formed bycareful chemical folding of the graphene sheet [ ] and indeed a host of possibilities [ ] .All these properties and avenues for further development present graphene as a sig-ni fi cant part of future technologies. Understanding its stability and dynamical proper-ties is consequently a very important part of the process of identifying effective use casesand searching the depths of its great potential. Compared to the problem of modelling DNA, graphene presents a number of new dif fi -culties. Foremost is that graphene is not only a two dimensional material, but the hexag-onal lattice structure requires some detailed consideration to model accurately. Thus thepotential governing the interatomic interactions needs to account for both the covalentbonding between atoms (not dissimilar to the bonding between the bases in DNA basepairs), but it also needs to fi rmly preserve the angles forming the hexagons. One simpli- fi cation over DNA modelling is that in graphene all atoms have the same equilibriumdistance, and disorder is only introduced through changing the masses of atoms (exclud-ing the considerably more complex task of modelling defects at an MD scale, which isnot considered here).To achieve accurate modelling, a number of approaches have been used, such as MDsimulations using empirical force fi elds [ ] and using Monte Carlo methods with theseforce fi elds [ ] . Tight-binding models have also been used with some success, partic-ularly in studying the nonlinear stress response and elasticity of graphene [ ] . Othermethods used to study the stability and strength of graphene include continuum me-chanics [ ] and ab initio calculations [ ] . All these methods have met with reason-able success, accurately reproducing many of the experimentally determined propertiesof graphene.In this work, we consider an in-plane model derived from fi rst principle density func-tional theory (DFT) calculations, where the empirically determined force fi elds are fi tted116 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Models i i + 1 jj + 1 ( i, j ) F IGURE
The hexagonal structure of graphene. Each of the atoms (repre-sented by dots) is at equilibrium separated by r = j ), while the zigzag columns are as shown by the red line (markingcolumn i ). The intersection of these is the point ( i , j ) , marked in blue. with analytical expressions to produce a relatively simple atomistic Hamiltonian modelof graphene [ ] . This model splits the molecular potential into two parts: An inter-atomic potential to model the covalent bonds, and an angle potential to account for theeffects of bending in the lattice, which maintains the hexagonal shape. Among the ad-vantages of this approach is that accurate estimations of the empirical potentials can befound by careful exploration of the DFT calculations, and consequently the analyticalexpressions for the potential can be closely fi tted, resulting in a very effective model.Before progressing further with the discussion of the model, it is useful to fi rst de- fi ne the lattice we use. Figure 6.1 illustrates the honeycomb structure of graphene, andintroduces the labelling convention for rows and columns in the model. One can thinkof an N x × N y graphene lattice as being constructed by taking N x zigzag columns with alength of N y atoms and placing them together into the honeycomb structure. This leadsnaturally to the x (horizontal) index being given by the number of the zigzag column,counting from the left – as shown in Fig. 6.1, marked as i and coloured in red. The y (vertical) index is then given by the number of the atom along the zigzag column count-ing from the bottom, marked with j , coloured green. The origin of the indexing systemis taken at the bottom left corner, with increasing indices up and to the right.With this structure, we notice that there are two distinct cases for the possible neigh-bours of each atom. Some atoms, such as the blue atom ( i , j ) in Fig. 6.1, have a neigh-bouring atom to the right as well as above and below. All atoms have a top and bottomneighbour, but depending on where along the zigzag chain a particular atom is, it canhave either a left or a right neighbour. These two possibilities are illustrated in Fig. 6.2,where the case with a left neighbour is shown in Fig. 6.2(a) and the right neighbour inFig. 6.2(b). So the blue atom in Fig. 6.1 is an example of case (b) in Fig. 6.2. Now re-calling that the left bottom corner atom in Fig. 6.1 is indexed as (
0, 0 ) , we see that thisatom is an example of case (a), noting that we are using periodic boundary conditionshere. Moving either up one atom or to the right by one atom (remember that the x axisis indexed by zigzag columns), we move from case (a) to case (b). In general, by lookingChapter 6 Malcolm Hillebrand 117haotic Dynamics of Polyatomic Models ( a ) ( i, j )( i − , j ) ( i, j + 1)( i, j − φ φ ( b )( i, j ) ( i + 1 , j )( i, j + 1)( i, j − φ φ F IGURE
A schematic representation showing the two possible cases foratoms in the lattice – either having a neighbour to the left or to the right.The fi rst case (a) occurs when the sum of the indices i + j is even, and case(b) when i + j is odd. Using our angle labelling convention, the angle φ in case (a) would be written as i , j φ i − ji , j + . at the points in Fig. 6.1, we fi nd that every point that has i + j even has a left neighbour [ case (a) ] , and every point with i + j odd has a right neighbour [ case (b) ] . We will usethis even / odd notation throughout to refer to the two cases.The equilibrium interatomic distance in graphene is known to be r = [ ] , and the equilibrium angles necessary to create the hexagonal lattice are all φ = π / fi culty being accounting forthe odd / even disparity. Adjusting for these two cases, we can write the position vectorsfor each atom in the lattice as r i , j = (cid:22) ( i − ) r (cid:8) cos (cid:13) φ (cid:14) + (cid:9) ( j − ) r sin (cid:13) φ (cid:14) (cid:23) , i + j even (cid:22) ( i − ) r (cid:8) cos (cid:13) φ (cid:14) + (cid:9) + r cos (cid:13) φ (cid:14) ( j − ) r sin (cid:13) φ (cid:14) (cid:23) , i + j odd. (6.1)As one more required piece of labelling, we need a way of referring to each angle inthe lattice. Each atom has three associated angles, formed between itself and its threenearest neighbours. Consider the angle labelled φ in Fig. 6.2(a). This angle is formedbetween the points ( i , j + ) , ( i − j ) and ( i , j ) , and centred at ( i , j ) . So we wouldlike to unambiguously refer to this angle, while clearly maintaining the fact that it iscentred at ( i , j ) . To this end we will de fi ne the slightly cumbersome but concise anglenotation i , j φ m , nk , l for the angle formed between the points A = ( m , n ) , B = ( i , j ) and C = ( k , l ) , where moving around the central point in a counterclockwise direction thepoints are ordered AB C . Thus the angle i , j φ m , nk , l = A ˆ B C , in the conventional geometriclabelling. With this method, we label the angle φ from Fig. 6.2(a) with i , j φ i − ji , j + . Notethat our counterclockwise convention is satis fi ed here; starting from ( i , j ) and moving118 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Modelsthe lattice, is given by a quadratic potential well for the angle, with a cubic correction [ ] : V b (cid:13) i , j φ m , nk , l (cid:14) = d (cid:8) i , j φ m , nk , l − φ (cid:9) − d ′ (cid:8) i , j φ m , nk , l − φ (cid:9) , (6.3)This expression has the distinct advantage of being a simple polynomial, but comes atthe cost of being expressed in terms of the angles themselves, rather than a convenientcartesian form. Nevertheless, the parameters d and d ′ can be found through the fi tting ofDFT calculations, and as seen in Fig. 6.3(b) this form of the potential accurately capturesthe behaviour of the bending response. The optimised values are d = / rad and d ′ = / rad .So we now have expressions for each contribution to the overall potential, and needonly to combine all these together to form the complete Hamiltonian. First consideringthe Morse potential components, we can see that summing over every point and addingup the contribution of all three bonds connected to the point will result in an overcount-ing by a factor of two – every bond will have its energy counted by the atoms on eitherend. The second consideration we need to make with this sum is to decide for each atomwhether the Morse potential needs to be added for the neighbouring atom on the rightor on the left. So for all points we will have contributions from the bonds going up anddown, but we have to discriminate on a case by case basis depending on the parity of i + j for the left or right bonds.For the bending potential, we have to make the same even / odd distinction, since thepotential of (6.3) needs to be evaluated with the correct three angles for each point. Assuch, the fi nal potential will be split into two sums: One over all atoms with i + j even,and one for i + j odd.With all this in mind, we can now write the Hamiltonian by adding in the kineticenergy, and then the two sums over the potential terms, ensuring to pass the correct an-gles as arguments to the bending potential V b . Thus the expression of the Hamiltonianis H = (cid:1) i , j m i , j (cid:24)(cid:13) ˙ x i , j (cid:14) + (cid:13) ˙ y i , j (cid:14) (cid:25) + (cid:1) i + j even (cid:26) (cid:8) V s (cid:13) r i , j + i , j (cid:14) + V s (cid:13) r i , j − i , j (cid:14) + V s (cid:13) r i − ji , j (cid:14)(cid:9) + V b (cid:13) i , j φ i − ji , j + (cid:14) + V b (cid:13) i , j φ i , j − i − j (cid:14) + V b (cid:13) i , j φ i , j + i , j − (cid:14) (cid:27) + (cid:1) i + j odd (cid:26) (cid:8) V s (cid:13) r i , j + i , j (cid:14) + V s (cid:13) r i , j − i , j (cid:14) + V s (cid:13) r i + ji , j (cid:14)(cid:9) + V b (cid:13) i , j φ i , j + i + j (cid:14) + V b (cid:13) i , j φ i + ji , j − (cid:14) + V b (cid:13) i , j φ i , j − i , j + (cid:14) (cid:27) , (6.4)where we have written in terms of the velocities ˙ x and ˙ y , with p x = m ˙ x and p y = m ˙ y .Note that the mass m is indexed as well, to allow for doping with othe carbon isotopessuch as C atoms. Unless otherwise stated, we use C atoms, with the mass taken as12u, but where we use C atoms the mass is 13u. In the two sums, we have included theappropriate stretching terms – up and down neighbours in both cases, and then the left120 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Modelsneighbour when i + j is even, and the right neighbour when i + j is odd. The bendingpotential terms are computed for the three angles that occur in each case, and here thereis no double counting.Having set up the Hamiltonian, we are now in principle ready to compute the equa-tions of motion, and potentially even the variational equations, and proceed with nu-merical simulations. While H (6.4) has a fairly involved form, it is at least possible toexpress it concisely, and the equations of motion take the usual form m i , j ¨ x i , j = − ∂ H ∂ x i , j , m i , j ¨ y i , j = − ∂ H ∂ y i , j . (6.5)The differentiation of the Morse potential V s (6.2) is straightforward enough, and theseterms can be computed easily. While the bending potential V b (6.3) poses more dif fi cul-ties due to the explicit angular dependence, the equations of motion can nevertheless becomputed exactly, and these have been presented in [ ] . Simulating the graphene lattice dynamics is a rather computationally challenging taskfor two main reasons:1. The complexity of the equations to integrate. The equations of motion themselvesinvolve performing a large number of operations, and either using the TM or the2PM approaches will exacerbate the integration time requirements.2. The relatively large number of atoms to be integrated in a 2D lattice (at least severalthousand in most cases for reasonable results).The integration is carried out using the fourth-order ABA864 symplectic integrationmethod [ ] , which has been found to be one of the most ef fi cient integrators for latticeintegration [ ] .In order to tackle the complexity of the task, we fi rst consider the choice betweenthe 2PM and the TM method for computing the mLE. For most models, the variationalequations are an ef fi cient way of estimating the mLE (see discussion in Section 1.3). Inthis case though, the variational equations are extremely complex (see [ ] ), and in factthe 2PM method is both simpler to implement and more ef fi cient. As such we made useof the 2PM for all mLE computations, following the method outlined in Section 1.3.1.This leaves us only integrating the equations of motion (albeit twice), but does requirea smaller time step in order to maintain the requisite accuracy for the small deviation.Since in general for the 2PM, we require a deviation that is somehow small, but doesnot go beyond machine precision, as discussed in Section 1.3.1 a lower limit on the normof the deviation vector has been suggested to be about 10 − [ ] . In our application, thelimit of 10 − is in fact somewhat smaller than can be accurately computed due to thelarge number of sites in the 2D lattice. For lattices of several thousand atoms such as wewill want to use, the deviation vector will have on the order of ten thousand components(since it has components in the x , y , p x and p y directions). So more realistically, for thecomponents to remain larger than 10 − in magnitude, a deviation vector norm of 10 − is necessary, although for small lattices a smaller deviation vector can be used.Chapter 6 Malcolm Hillebrand 121haotic Dynamics of Polyatomic Models No. of CPUs . . . . . . Sp ee dup F a c t o r × × × F IGURE
The total speedup factor as a function of number of CPU cores,simulating the graphene lattice up to a time of 10 ps, integrating twonearby initial conditions and computing the ftmLE. Three lattice sizes areshown here, N x = N y = ( N = N x · N y = ) (blue circles), N x = N y = ( N = ) (orange squares), and N x = N y = ( N = ) (green triangles). The dashed straight line shows the perfect speedup rate.We see that the speedup drops off quite noticeably after 4 CPU cores forthe 36 ×
24 cases, and although for the larger lattice even up to 10 cores stillprovides a distinct improvement, after 4 cores the speedup becomes less ef- fi cient. The speedup past 12 CPUs is quite minimal in all cases, althoughfor the 80 ×
60 lattice there is still slow improvement.
It is also signi fi cant that even at the lowest energies considered in this investigation,the average fl uctuations in atom positions and velocities are on the order of 10 − (wherepositions are measured in Å and velocities in Å / ps). So taking the deviation norm of10 − , if we have a lattice of more than a thousand atoms, as is typically the case, thenthe individual components of the deviation vector will be on the order of 10 − , which isvery clearly smaller than the typical fl uctuations, providing further assurance that thisdeviation is indeed small on the scale of the natural dynamics of the lattice.The issue of a large number of atoms with complex equations of motion to be sim-ulated is perfectly suited to parallel computing, since the lattice can be broken up intoblocks of atoms that are integrated in parallel. This is possible due to the fact that changesin the position and momentum of a particular atom only affect its neighbours in the nextiteration of the integration, allowing each integration step to be performed for each atomindependently (see Section 1.2.2). By parallelising the integration, we can substantiallyspeed up the integration process, especially as the lattice size increases.In order to fi nd a near-optimal number of CPU cores to use for the integration pro-cedure, we performed scaling tests with varying numbers of CPU cores comparing thetotal wall time (real time elapsed, rather than CPU hours) taken to integrate the latticeup to a common fi nal time. These tests were performed with three different lattice sizes,to get a better grasp of what is optimal for different sizes, especially since we expect largerlattices to be more amenable to parallelisation speedup. Figure 6.4 shows the speedupobtained (where a speedup factor of 2 corresponds to halving the compute time, and afactor of 3 to a reduction to one third, and so on) by using multiple cores to parallelise theintegration, with lattices of N x × N y = × =
864 atoms (blue circles), 60 × = × = fi rst a signi fi cant speedup observed in all cases. Theaddition of a second core is almost perfectly ef fi cient, with a speedup factor of near 2.This keeps the point near the line of optimal speedup, which is when we have a speedupfactor equal to the number of cores used. Going to 4 cores is not quite as ef fi cient, but isnevertheless a strong improvement, with a speedup of close to perfect for the two largerlattices, with the 36 ×
24 lattice (blue points) starting to show a decrease in ef fi ciency.For this lattice, using more than 4 cores seems to be fairly wasteful, as the speedup dropsoff dramatically after this point. For the larger lattices however, there is still merit in go-ing beyond 4 cores, even if the ef fi ciency decreases. Up to around 10 cores, the speedupfrom each extra core is noticeable, and a speedup by a factor of 7 or 8 is de fi nitely useful.Past this point though, while there is some speedup it is highly inef fi cient.These results suggest that using 4 cores for most simulations will yield fairly ef fi cientspeedup, dramatically reducing the time needed to perform the computations, but notimposing a large overhead cost on the parallelisation. We are also then able to run 6simulations simultaneously on each node.Another aspect of the scaling to consider is a closer look at how the compute timescales with size, for a fi xed number of CPU cores, and to compare this with using aGPU. As a fi rst step, we will present some comparative results using CPU and GPUresults, primarily to con fi rm that the GPU implementation is in fact producing correctresults. Two main results are presented in Fig. 6.5(a) and (b), the relative energy error,and the estimation of the ftmLE. In both cases we see that the two computation methodsproduce near-identical results, with some small variations inevitable in the simulationof a chaotic system. Importantly, the energy remains conserved up to a small constanterror. This is a strong indication that for both the CPU and GPU implementation theequations of motion are being accurately integrated. The estimation of the ftmLE isalso identical up to the fi nal time of 10 ps considered here, where the value has saturatedto the fi nal mLE. In both these measurements there are slight differences inevitable inthe simulation of a chaotic system, but are effectively identical. These results allow us tocontinue with the GPU code in con fi dence that the results being produced are consistentwith the CPU results which have previously been carefully tested.In Fig. 6.5(c), we see the effect on the total simulation time of increasing the size of thelattice. Results are shown for integration with 4 and with 24 CPU cores (green trianglesand orange squares respectively), as well as the GPU implementation (blue circles). Fromthe results of Fig. 6.4, we have already seen that larger lattices bene fi t more from the useof many CPU cores. In the CPU results of Fig. 6.5(c), this is borne out by the muchfaster increase in compute time required when using only 4 CPUs (green triangle) whencompared to the 24 CPU integration (orange squares). We also see that very quickly the24 CPU parallelisation becomes faster than the 4 CPU case, requiring around a third ofthe time to integrate a lattice of 4096 atoms. The GPU parallelisation [ blue circles inFig. 6.5(c) ] is clearly substantially faster when it comes to integrating very large lattices,and is even faster to integrate the 4096 atom lattice. For very small lattices it is hard tocompare the times, but since the GPU takes the same time to integrate all lattices smallerthan 5120, it is not at all ef fi cient for integrating lattices of only a few thousand atoms.It is however unfair to just compare the net compute time, for a number of reasons.Most obviously, the comparison between the 4 and 24 CPU cases ignores the fact thatwe can run 6 instances of the simulation at once when using only 4 CPUs, while weChapter 6 Malcolm Hillebrand 123haotic Dynamics of Polyatomic Models t (ps) . . . . . | H ( t ) − H ( ) | / H ( ) × − (a) GPUCPU t (ps) − − χ ( t )( p s − ) (b) GPUCPU N . . . . . . W a ll t i m e ( h o u r s ) (c) GPU24 CPUs4 CPUs 0 10000 20000 30000 N C P U t i m e ( h o u r s ) (d) GPU24 CPUs4 CPUs F IGURE (a) The comparison of the time evolution of the relative en-ergy error | H ( t ) − H ( ) | / H ( ) for identical initial conditions in a latticeof N x × N y = × = − , as required. (b) The time evolution of the ftmLE χ ( t ) forboth computation methods, for the same initial condition. Again we seethat they both produce the same results, con fi rming that the GPU mLEcomputation algorithm is also correctly implemented. (c) The total com-pute time required to integrate lattices of size N atoms up to a fi nal timeof t = ps, using a GPU (blue circles), 24 CPU cores (orange squares)and 4 CPU cores (green triangles). The lines connect the points to guidethe eye. From this plot it is clear that for lattices of more than 4000 atomsthe GPU requires the least time to integrate the system, and the 24 CPUsare faster than 4 CPUs. (d) The rescaled CPU time (see text) required tointegrate the system, for the same three cases as (c). Now it is clear that the4 CPU integration is more ef fi cient for small lattices, of fewer than 8000atoms. After this point the GPU is more ef fi cient despite each GPU hourbeing equivalent in cost to 40 CPU hours (see text). As the lattice sizeincreases, the GPU has by far the smallest gradient, indicating its strongperformance for very large lattices.
124 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Modelscan only run a single instance with 24 CPUs. Similarly, the GPU is a much more inten-sive processor than a single CPU, and comparing even 24 CPUs in parallel with a singlehigh-powered GPU is not necessarily a reasonable comparison. Relating CPU hours toGPU hours is thus a slightly tricky issue. We would then like to answer the questionof “what is the most ef fi cient way to integrate the lattice at different sizes?” by somede fi nite measure. One logical such measure is to use the CPU hour - GPU hour rela-tionship implemented by the computing cluster, which equates 1 GPU hour to 40 CPUhours [ ] . Thus we can normalise the compute times to be in terms of CPU hours,which yields a much fairer estimation of the compute cost for performing many simula-tions. Of course if only one or two simulations need to be performed, then the resultsof Fig. 6.5(c) will apply more pertinently.These rescaled results are shown in Fig. 6.5(d). With this more practical approach,the 4 CPU case demonstrates much more ef fi cient performance for lattices of fewer than8000 atoms. The 24 CPU parallelisation is never more ef fi cient than either of the alter-natives, illustrating the problem with simply scaling up CPU multithreading. Despitethe increase by a factor of 40 though, the GPU still massively outperforms the CPUcode for atoms of more than 10000 atoms. This is consistent with the intended designof the GPU, which is to be ef fi cient in performing highly parallelisable computationsfor very large systems of equations, which is exactly what we have with lattices of tensof thousands of atoms. As such, for any very large lattice simulations the GPU codeis undoubtably the best choice, but for small lattice computations which are generallysuf fi cient for our purposes, using around 4 CPUs is the most ef fi cient option.From this section, we can take away several conclusions. For most lattice sizes, using4 or 6 CPU cores yields an ef fi cient speedup, providing a good balance between reduc-ing run times and not wasting CPU time on inef fi cient parallelisation. Simulating thesystem using GPU parallelisation is very effective for large lattices, and would be anef fi cient choice for investigations using lattices of more than around 5000 atoms. Fur-ther optimisation of the GPU code would also allow for even better ef fi ciency in theintegration of large lattices; for instance careful pro fi ling to fi nd inef fi cient computa-tion blocks, and further minimisation of GPU-CPU data transfers. For our applicationhowever, where we do not need to consider such very large lattices, the best choice ofparallelisation method is to combine OpenMP multithreading with GNU Parallel torun several parallelised CPU integrations alongside each other. Having set up the graphene model and outlined the numerical processes, we can proceedwith some fi rst results of simulating the system. While the focus of this investigation isthe chaotic dynamics of the lattice, in order to con fi rm the validity of the used modeland our implementation of the equations of motion, we will reproduce some knownresults for the dynamical behaviour of graphene.In all these results, we impose periodic boundary conditions on all edges to elimi-nate edge effects. The simulations are run using N x =
60 atoms in the x direction and N y =
48 atoms in the y direction, but changing the lattice size has no effect on the con-sidered quantities, when using periodic boundaries. For each energy density E n (totalenergy divided by the number of atoms), we use ten different setups with random ini-tial conditions (random excitations of the momentum in both the x and y directions,Chapter 6 Malcolm Hillebrand 125haotic Dynamics of Polyatomic Models t (ps) T ( K ) (a) . . . . E n (eV) T ( K ) (b) C C F IGURE (a) A single case of the evolution of the lattice temperature intime, for energy density E n = T and energy density E n for the graphene model. The straightline shows the linear function of (6.6), con fi rming that the relationshipholds not only at low temperatures, but even extending far beyond com-monly experienced temperatures. The blue circles show results for thelattice using a mass of 12u ( C), while the red squares show the resultsfor a graphene sheet comprised entirely of C atoms. We see that forboth masses the same energy-temperature relation is found, as the pointsoverlap and are indistinguishable. rescaled to give the required total energy value), run for 5ns to thermalise, and then datacollected every 10ps for the next 5ns.The fi rst result we consider is the energy density-temperature relationship. Since weare still working in the microcanonical ensemble with constant energy, we estimate thetemperature T of the 2D system according to T = Kk B N , (6.6)with K being the total kinetic energy, the system having a total of N atoms and k B beingthe Boltzmann constant.Figure 6.6(a) shows the estimated temperature through time for a single run, at anenergy density of E n = fi fty data points, we can fi nd an average temperature measurementfrom this run.By computing this for several initial conditions (ten in this case), and averaging overall the data points, we can numerically fi nd the relationship between the energy density E n and temperature T . These results are shown by the blue circles in Fig. 6.6(b). Therewe see a purely linear relationship, extending far beyond the low-temperature region (aswas seen in the DNA in Fig. 5.4). This straight line is in fact given by the low-temperaturegradient for a 2D system of E n = k B T , (6.7)126 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Modelswhere k B is the Boltzmann constant. The fact that no correction is needed to thislow-temperature relation even when the temperatures reach extreme levels above 3000Kbears testament to the stability of the graphene lattice. It is also very likely that the re-striction of the model to in-plane motion adds to the stability of these high-temperatureresults, as the magnitude of out-of-plane motion is typically much greater at high tem-peratures.In order to test the model for doping with C atoms, following the same procedurewe also compute the temperature for a graphene lattice comprised entirely of C [ redsquares in Fig. 6.6(b) ] as the most extreme case of this doping. We see that for this too thelinear relation is maintained across the full energy density spectrum, with no discernibledifference between the energy-temperature relation of the two isotopes.Apart from the energy-temperature, we can look at the distributions of interatomicdistances and atom speeds. For both of these quantities, we have a theoretical expecta-tion for their distributions, which will allow us to support both the physicality and theaccuracy of our computations. First, the distribution of the interatomic distances in thelattice should be normal [ ] . So after equilibration of the lattice, we should fi nd thatthe distances r between neighbouring atoms are distributed according to P ( r ) = σ (cid:8) π e − ( r − µ σ ) , (6.8)where µ is the mean value (expected to be equal to r = σ is the standarddeviation.Figure 6.7(a) shows the distribution of distances r , using fi ve different initial con fi g-urations and 500 snapshots from each run. The distributions are shown for fi ve energydensities, both to check that the Gaussian form is maintained, and also to see how theshape of the curve changes. In each case, the symbols show the distributions from thesimulations, and the solid curves show a fi tting with (6.8). The normal distribution fi tsthe data exactly for all energies, and correctly fi nds a mean value of µ = r = E n = fi t yieldsa correspondingly small standard deviation. As the energy increases, the distribution fl attens, but maintaining the mean value of µ = fi rm theoretical prediction, in the form of theMaxwell-Boltzmann distribution [ ] . We expect that the speeds should be distributedaccording to the distribution P ( v ) = mvk B T e − mv kB T , (6.9)since we have a 2D lattice. This is of course fi xed for a given temperature (or in lightof the linear energy-temperature relation [ (6.7) and Fig. 6.6(b) ] , a given energy density),meaning that we should fi nd the atom speed distributions to follow (6.9) without anyneed for fi tting. In Fig. 6.7(b) we see the theoretical prediction (solid curves) along withthe numerical results obtained from our simulations (symbols) for several energy den-sities. In all cases, the speed distributions are exactly as they should be, following theMaxwell-Boltzmann distribution precisely. This serves as another con fi rmation of ourtemperature results, as fi tting (6.9) to the distributions would be an alternative methodof computing the temperature for a given energy density.Chapter 6 Malcolm Hillebrand 127haotic Dynamics of Polyatomic Models . . . . . r ( ˚ A) P D F (a) E n = 0 . E n = 0 . E n = 0 . E n = 0 . E n = 0 . . . . . . . v ( ˚ A/ps) P D F (b) E n = 0 . E n = 0 . E n = 0 . E n = 0 . E n = 0 . F IGURE
The distribution of (a) the distances r between atoms, and (b)the particle speeds v in the graphene lattice. For each case the results areshown at fi ve energy densities, E n = E n = E n = E n = E n = r = fi t-ting well in each case. As the temperature increases the distances spreadmore widely, consistent with the overall more active dynamics. The ve-locities also show the expected distribution, fi tting closely to the Maxwell-Boltzmann distribution, (6.9). Like the distances, the velocity distributionis broader for higher temperatures. With these computational methods and validated results, we are able to estimate the mLEof graphene lattices of different varieties, and investigate graphene’s stability and chaoticdynamics. Here we will study graphene in two forms: A periodic sheet to approximatethe dynamics of an in fi nite lattice, and fi nite width GNRs, with one pair of edges leftfree while the other two are joined with periodic boundary conditions. As our fi rst task, we would like to understand the overall chaotic dynamics of graphene,without considering boundary effects or fi nite size issues. To this end, we will simulatean in fi nite graphene lattice by taking a reasonably large graphene lattice and applyingperiodic boundary conditions to all the sides. The resultant “torus” will allow for evenspreading of energy, and should provide an indication of the behaviour of an in fi nitelattice.The immediate question we need to answer is what constitutes a reasonably largelattice. In order to quantify this, we performed several simulations with varying latticesizes, and computed the mLE at each size. Once we reach a point where the mLE nolonger changes with lattice size, this will give a strong suggestion that we have reacheda suf fi ciently large lattice size that the effects due to the fact that the sheet is not trulyin fi nite are negligible.Figure 6.8(a) shows the evolution of the ftmLE at an energy density of E n = N x × N y = × =
384 (solid blue line), 36 × =
864 (dashed128 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Models t (ps) − − χ ( t )( p s − ) (a) × × × ×
60 80000 90000 1000009 . . × − t (ps) − − − χ ( t )( p s − ) (b) × . . . × − F IGURE (a) The evolution of the ftmLE χ ( t ) for a periodic graphenelattice, for four lattice sizes, N x × N y = × =
384 (solid blue curve),36 × =
864 (dashed orange curve), 60 × = × = E n = fi ed. Wesee that for all lattices, the fi nal mLE value is the same up to a precisionof around 5 × − , but the 24 ×
16 lattice systemically yields a slightlyhigher value. The 60 ×
48 and 80 ×
60 lattices however give extremelysimilar results. These results also illustrate that for E n = fi nal value by a time of 10 ps. (b) The time evolutionof the ftmLE χ ( t ) for a periodic graphene lattice, at an energy density of E n = N x × N y = × = ps the value has saturated to aconstant, as con fi rmed by the inset which shows a magni fi cation of thelast stage of the evolution. orange line), 60 × = × = fi nal value of the ftmLE saturates atthe same value up to a spreading of 1 × − ps − , as seen in the inset. The blue curveof the smallest lattice size, 24 ×
16, does exhibit a slightly larger fi nal saturated value forthe mLE, which we have also observed in other cases, suggesting that this size is likelytoo small to be reliable as a good estimator of the dynamics of the in fi nite lattice. The36 ×
24 lattice, given by the dashed orange curve, gives an mLE value very slightly belowthe two larger lattices, but within a range of 2 × − ps − in a measurement on the scaleof 10 − ps − , which is negligible in the fi nal analysis. However, the two largest lattices,60 ×
48 (green dotted curve) and 80 ×
60 (red dash-dotted curve) are almost identicaleven in the magni fi ed image presented in the inset. This inspires con fi dence that theeffects of the fi nite lattice size have been neutralised by the time we have a lattice of N x × N y = ×
48 atoms, and most likely even before this. Since the mLE values nolonger change even at a tiny scale past this point, we will use this lattice size of 60 × fi nal points of the ftmLE as an estimate of the actual mLE forthat initial condition, we can average over multiple simulations to construct an averageChapter 6 Malcolm Hillebrand 129haotic Dynamics of Polyatomic ModelsmLE for a given energy density.To be sure that all simulations have completely converged to the fi nal mLE value,we have integrated most cases up to a fi nal time of 10 ps [ as in Fig. 6.8(a) ] , as this gives aclear saturation. For particularly small energies ( E n = E n = ps, as these cases take longer to thermalise andfor the ftmLE to saturate. In Fig. 6.8(b), we see the evolution of the ftmLE up to a timeof 10 ps for E n = fi nal time.The inset in Fig. 6.8(b) shows the last stage of the ftmLE evolution for all the consideredcases, clearly displaying that the ftmLE has ceased to decrease and has stabilised at the fi nal value of the mLE. We also note that the choice of initial conditions for the latticehas no actual bearing on the mLE value, with different types of initial excitation trialledresulting in identical fi nal ftmLE values.To study the overall dynamics of the system now, we can take exactly the same ap-proach as for the DNA (see Section 4.2), and collate average mLE values from multipleruns to build up a picture of the energy dependence of the system’s chaoticity. To thisend, we will perform the same computations shown in Fig. 6.8(b) for multiple energydensities E n , and average over the ten runs to fi nd an average mLE value for each E n value. We will use the fi nal value of the ftmLE, averaged over the ten runs, with anerrorbar given by the standard deviation of these data as an estimation of the system’smLE X .Figure 6.9 presents the average mLE X across the energy spectrum for the periodicgraphene lattice. The computation of the mLE was performed for both a lattice of purely C atoms, shown by red squares in Fig. 6.9 and C atoms, shown by blue circles. Wesee that in both cases the mLE increases monotonically with the energy density, withthe C lattice showing slightly more stable behaviour, arising from the higher mass ofeach atom. Unlike the DNA, here we see no changing in the behaviour of the mLE, theincrease with energy (or equivalently temperature, since the relationship between thetwo quantities is linear) is smoothly monotonic. This implies that there is no changein the structural stability in the graphene shell in the range of energy densities consid-ered here, which cover temperatures up to around 4000K. The extraordinary stability ofgraphene is once again illustrated here, with even extreme conditions unable to provokeany irregularity in its dynamics (although investigating the effects of out-of-plane atomicmotion would be interesting here).At low energies the mLE is more or less linearly dependent on the energy density E n ,curving slightly upwards as it moves past E n ≈ fi tted accurately with a simple quadratic function, X ( E n ) = β E n + γ E n . (6.10)where β and γ are free constants to be determined. Performing this fi tting for C we fi nd β = ± − eV − , and γ = ± − eV − . Similarly, for C we fi nd that β = ± − eV − , and γ = ± − eV − .Due to the similarity of the mLE values, these fi tted values are fairly close together, onceagain demonstrating the minor effect that C doping has on the chaotic dynamics andstability of the lattice. The small variance in each coef fi cient also re fl ects the closenessof the quadratic fi t.Considering the Lyapunov time T L = / X of the system (see Section 1.3), also pro-130 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Modelsvides some useful insight into the process of chaotisation in graphene. Estimating T L byinverting the mLE values in Fig. 6.9(a) will give us the time (in picoseconds) by whichthe system has reached a chaotic state. What we are really interested in however is howthis time relates to the natural time scales of the lattice. This will allow us to understandthe chaotisation time in the speci fi c context of the graphene sheet.In the lattice, one clear time scale is speci fi ed by the frequency of the normal modes,and in particular the highest frequency of the normal modes. The normal modes of themodel used here have been computed through phonon dispersion relations [ ] , andwe can use the highest frequency thus obtained, the transverse optical (TO) mode for thetime scale of the system. Measured in units of cm − , using wavenumber notation, thefrequency of this mode is around ω W ≈ − . In order to convert this to an actualfrequency in Hz, we rescale ω = π c ω W , where c is the speed of light. With thisrescaling, we fi nd a frequency of ω ≈ τ ≈ × − psfor each oscillation of the system from inverting this frequency.So now, instead of just presenting the Lyapunov time T L , we are able to present therescaled Lyapunov time (cid:28) T L = T L / τ , (6.11)which gives us the number of oscillations of the highest frequency mode before chaossets in. Figure 6.9(b) shows (cid:28) T L as a function of the energy density E n , demonstrating thevery large number of normal mode oscillations needed before chaotisation occurs. Atthe lowest energies, more than 10 oscillations are required before the system becomesfully chaotic, with the room temperature case of E n = × oscillations. As the energy increases, chaotisation happens more and more quickly, buteven at the extreme energies of E n ≈ oscillations are needed. Fromthe results of Fig. 6.9(b) we can say that in graphene, the onset of chaos is very slow,which is consistent with the material’s stable behaviour. Linearisation of the Model
As a short aside, we will brie fl y consider the possibilities for simplifying the model, inorder to better understand the source of chaos in graphene. Considering the potentialfunctions of eqs. (6.2) and (6.3), there is evidently nonlinearity inherent in both of them.The question is then what the effects are of the 2D hexagonal geometry itself, as well asthe actual functional form of the potentials. For instance, if we kill the bending potential V b (6.3) completely and linearise the Morse potential to a harmonic coupling, how doesthe mLE behave? And what effect does the quadratic part of the bending potential V b have; since it is quadratic in the angles, it is certainly nonlinear in the coordinates x and y , but it is not clear how this will interact with the interatomic coupling.In order to investigate these effects, we performed simulations with three variationsof the full nonlinear potential. These three cases are as follows:1. Only the harmonic interatomic coupling: V s = k ( r − r ) . (6.12)2. Only a harmonic interatomic coupling with a quadratic angular potential: V s = k ( r − r ) , V b = k ′ ( φ − φ ) . (6.13)Chapter 6 Malcolm Hillebrand 131haotic Dynamics of Polyatomic Models . . . . . E n (eV) . . . . . . . . X ( p s − ) (a) C C . . . . . E n (eV) (cid:1) T L (b) C C F IGURE (a) The average mLE X as a function of the energy density E n ,for the periodic graphene lattice. Results are shown for a lattice composedsolely of C atoms (red squares) and C atoms (blue circles). The error-bars are given by the standard deviation of the mLE from each run, butsince the values are extremely close the errors are not visble on the scale ofthe plot. The solid red and dashed blue curves show a fi t to (6.10) for Cand C respectively. (b) The rescaled Lyapunov time (cid:28) T L (6.11) shown ina semilogarithmic plot on the y axis, counting the number of oscillationsof the fastest normal mode of the lattice before chaotisation occurs. The fi tted curves are given by (6.10) rescaled like the data.
3. A full linearisation of the equations of motion, by fi rst order approximation witha Taylor expansion. To obtain the true linearisation, we consider the equations ofmotion of the Hamiltonian (6.4) in vector form. If the positions of the atoms inthe lattice are given by a vector u = (cid:13) x , x , . . . , x N x ,1 . . . , x N x , N y , y , y , . . . , y N x ,1 . . . , y N x , N y (cid:14) ,then the equations of motion are given by the accelerations of the atoms, m ¨ u i = − ∂ H ∂ u i ,which can be compressed into a vector form as m ¨ u = − ∂ H ∂ u . (6.14)We note that this expression corresponds to the equations of motion in (6.5). Inthis form, we are able to perform a Taylor expansion around the equilibrium po-sition u . Keeping only the term in fi rst order of u − u , we have m ¨ u = − ∂ H ∂ u ( u ) − ∂ H ∂ u ( u )( u − u ) . (6.15)Evaluated at the equilibrium positions, all the values in the vector ∂ H / ∂ u i arezero, and the matrix of elastic coef fi cients is given by Q i , j = ∂ H ∂ u i ∂ u j . (6.16)132 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Models t (ps) − − − − ¯ χ ( t )( p s − ) Full Nonlinear Model V s = k ( r − r ) , V b = 0 (Case 1) V s = k ( r − r ) , V b = k ′ ( φ − φ ) (Case 2)Linearised Model (Case 3) F IGURE
The ftmLE averaged over ten initial conditions, for four differ-ent variations of the graphene model at E n = ( t ) / t . We thus have the fi nal linearised equations of motion as¨ u = − Q m ( u − u ) . (6.17)For each of these cases, we can compute the system’s mLE, using k = D for the har-monic potential and k ′ = d / fi ed versions of the model, as well asthe full model. Each curve corresponds to the average over ten simulations with randominitial conditions, at an energy density of E n = fi nd that the addition of aquadratic angular potential (case 2 – green dotted curve) stabilises the system dramati-cally, with the mLE of this potential being almost two orders of magnitude below thefull system. Finally, the result of the linearised equations of motion (case 3) is displayedby the dash-dotted red line, demonstrating the expected continuous decrease towardszero as in this case we have no chaos at all.From these results, a number of interesting points arise. The simplest form of thepotential, the simple harmonic coupling (case 1), actually results in the most chaoticbehaviour. This is at fi rst remarkable, as we expect that harmonic oscillators shouldproduce linear behaviour. However, due to the 2D geometry of the system, and possiblyexacerbated by the hexagonal structure, there is a large amount of motion transverse tothe direction of the oscillators. Including the full 2D range of motion in the potential,without ignoring this transverse motion, results in a strongly nonlinear interaction.The addition of the quadratic angular potential (case 2) results in a signi fi cantly morestable system. At fi rst this is counterintuitive as well, since this is the part of the poten-tial we expected to be nonlinear. However, this potential V b acts as a potential wellChapter 6 Malcolm Hillebrand 133haotic Dynamics of Polyatomic Modelsfor each of the angles, which maintains the hexagonal shape much more strongly thansimply allowing the atoms to move around constrained only by the springs (atom-atominteractions). The effect of this added rigidity is that the motion is restricted, and thustransverse displacements are a lot smaller and the system closer to a linear behaviour.The fact that case 2 is still nonlinear is due to both the small transverse motions as wellas the inherently nonlinear behaviour of the angular potential.Finally, the ftmLE of the completely linearised potential decreases towards zero aswe expect, since in this case the system is actually linear. The slope of ln ( t ) / t is not assteep as the 1 / t slope that is expected for a non-growing deviation vector, but this is justdue to a linear growth of the deviation vector [ ] , which still meets the requirementfor regular motion of a subexponential deviation growth (see Section 1.3). Apart from the in fi nite sheet or bulk behaviour of graphene, we can also study the dy-namics of a very physically relevant graphene structure – the GNR. A GNR is sim-ply a narrow strip of graphene, with the transverse direction containing relatively fewatoms, and the longitudinal direction containing a large number of atoms. In our simu-lations, we will approximate the very long ribbon using periodic boundary conditionsonce again in one direction, and the leave the ribbon edges free to move in the otherdirection. This of course leaves us both options as to which edges we leave free – thearmchair or the zigzag edges. The choice of which edge is free describes the type of rib-bon, and thus the labels “armchair” and “zigzag” GNRs refer to the free edge. Harkingback to Fig. 6.1, we can see that by choosing the free edges on top and bottom, we wouldhave an armchair GNR. If we instead had the left and right edges free, this would be azigzag GNR.In the same vein as with the periodic sheet, we can build up a plot of the averagemLE of the GNRs by performing multiple runs at each energy density, and averagingthe mLE values from each run to provide a single data point. Here however we have anadditional parameter that we are concerned with, beyond the energy density. For theGNRs, we would like to study how the width affects the stability of the lattice, rangingfrom very narrow ribbons of only a few atoms to larger ribbons where the dynamicsshould start to approximate the in fi nite lattice. Based on the geometry of the hexagonallattice (see Fig. 6.1), we can calculate the width of the nanoribbons using some trigonom-etry, yielding W A = (cid:8) (cid:6) N y − (cid:7) r /
2, (6.18)for the width of an armchair ribbon, while W Z = ( N x / − ) r , (6.19)for their zigzag counterparts. Here again N x refers to the number of atoms in the x direction (the number of zigzag columns), and N y to the number of atoms in the vertical y direction.So now we can perform a series of computations, estimating the mLE for a) differentenergy densities, and b) different ribbon widths, for both GNR types. For these runs,we used periodic boundary conditions on the longitudinal direction, with a length of 60atoms, as increasing the length beyond this was found to have no effect on the dynam-ics. The results of these runs are presented in Fig. 6.11, showing the mLE as a function134 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Models W (nm) . . . . . . . . X ( p s − ) (a) E n = 0 . E n = 0 . E n = 0 . E n = 0 . E n = 0 . W (nm) . . . . . . . . X ( p s − ) (b) E n = 0 . E n = 0 . E n = 0 . E n = 0 . E n = 0 . F IGURE
The average mLE of GNRs as a function of ribbon width W in nm, for several energy densities. (a) Armchair ribbons, and (b) zigzagribbons. The horizontal dashed lines indicate the value of the mLE of theperiodic sheet at that energy density, which are shown in Fig. 6.9(a). Thedata for each E n value in both plots are fi tted with (6.20), shown as solidcurves. of ribbon width for several energy densities: E n = E n = E n = E n = E n = E n = E n = .
05 0 .
10 0 .
15 0 . E n (eV) . . . . A ( p s − ) (a) ArmchairZigzag 0 .
00 0 .
05 0 .
10 0 .
15 0 . E n (eV) . . . n (b) ArmchairZigzag F IGURE
The parameters of (6.20) fi tted to the data of Fig. 6.11, for arm-chair (green squares) and zigzag (blue circles) nanoribbons across severalenergy densities. (a) The coef fi cient A , and (b) the exponent n . The dashedlines are to guide the eye. For both GNR types, the mLE values can be fi tted quite accurately with a decreasingHill function, with a constant limiting value added, X ( W ) = A + W n + X b . (6.20)This constant is simply given by the periodic sheet mLE value for the particular energydensity [ see Fig. 6.9(a) ] . As seen by the solid fi tting curves in Fig. 6.11, this functionalform captures the rate of decrease effectively, allowing us to understand how the increas-ing ribbon width affects the chaoticity of the lattice at different energies, and to comparethe behaviour of the two GNR types.The computed value of the two parameters, A and n , of the fi tting are given inFig. 6.12. The scale coef fi cient A is shown in Fig. 6.12(a), exhibiting a distinct increasewith energy density commensurate with the overall trend of increasing mLE values withenergy. This parameter A is related to the range of values exhibited by the mLE, whichcan be seen more clearly by rewriting (6.20) as X ( W ) − X b = A / ( + W n ) . Because ofthis, the larger values of A seen for the armchair GNRs are consistent with what we see inFig. 6.11, where the range of mLE values is larger for the armchair case [ Fig. 6.12(a) ] thanthe zigzag ribbons [ Fig.6.12(b) ] . For instance, let us take the extreme case of E n = [ purple crosses in Figs. 6.11(a) and (b) ] . Because of the anomalous decrease in the mLEfor the smallest width zigzag ribbons [ W ≈ ] , the largest mLEvalue is just below 0.005ps − , while the limiting value is approximately 0.0032ps − . Thisdifference of ∼ − is substantially smaller than the corresponding difference ofaround ∼ − for the armchair case [ Fig. 6.11(a) ] .The exponent n decreases with increasing energy density, showing that for higherenergies there is in fact a relatively slower decrease towards the limiting value. Onceagain, the armchair ribbons have consistently larger values of n compared to zigzagGNRs, which demonstrates that the mLE of the armchair ribbons tends to decreasemore quickly towards the bulk value. Although this rate of decrease is systemicallyfaster in armchair GNRs, for both cases there is an eventual decrease to the asymptoticvalue corresponding to the mLE of the in fi nite sheet.136 Malcolm Hillebrand Chapter 6haotic Dynamics of Polyatomic Models In this chapter, in order to study the chaoticity of graphene we have used a Hamilto-nian model which very effectively simulates the dynamics of the material while remain-ing relatively easy to use and computationally manageable. The unavoidable computa-tional dif fi culties inherent in studying large two dimensional lattices, compounded bythe distinctly non-trivial equations of motion arising from the geometrical complexity ofthe graphene model, were met using CPU multithreading and optimisation techniques,which allowed us to integrate a large number of atoms (typically on the order of 3000)in a reasonable time frame. A GPU implementation of our code was also written usingCUDA, which drastically outperforms the CPU code for very large lattices, consistingof tens of thousands of atoms. Future work that involves integrating such large systemswould de fi nitely bene fi t from using the GPU implementation of the code rather thanCPU parallelisation.We have computed the mLE of various graphene structures using the 2PM, due tothe extreme complexity of the variational equations. The mLE of an in fi nite graphenelattice was estimated by approximating this system through the use of periodic bound-ary conditions, and the mLE was found to increase quadratically with energy density [ Fig.6.9(a) ] . Replacing the C atoms in the lattice with C atoms results in a slightlymore stable structure due to the higher masses, with a marginally smaller mLE valueestimated. By inverting the mLE, Lyapunov times were also explored, and compared tothe natural time scales of graphene lattices given by the highest frequency normal modes.While chaos sets in faster as energy density increases, the number of oscillations of thenormal modes with the highest frequency required before the system is fully chaotic isover 10 even up to temperatures of 4000K [ Fig. 6.9(b) ] .Some simple modi fi cations of the model have been investigated, by approximatingthe nonlinear Morse potential for interatomic interactions with a simple harmonic cou-pling, as well as by using only a quadratic angular potential for the bending potential.We found that a hexagonal lattice with only harmonic interatomic coupling was morechaotic than even the full nonlinear system, while the addition of a quadratic angularpotential to the harmonic coupling dramatically stabilises the system by discouragingmotion transverse to the oscillators (Fig. 6.10). We have also con fi rmed that the fulllinearisation of the equations of motion using a Taylor approximation leads to regularbehaviour.The stability of GNRs was also studied using the mLE, with both edge types – arm-chair and zigzag – investigated. We found that the mLE decreases as the GNR widthincreases for both cases (Fig. 6.11), with both the range of mLE values and the rate of de-crease being higher for armchair ribbons than zigzag, quanti fi ed by fi tting to a decreasingHill function (6.20). In both cases the mLE converges to the value of the in fi nite sheetasymptotically as the width increases.There is plenty of scope to expand this work, and to use the model (6.4) for furtherinvestigations. One example would be studying the effect of stretching or compressionof the lattice, as computing the mLE under these conditions could lead to insight aboutthe stability of graphene. The addition of a term in the potential to account for out-of-plane deformations, as outlined in [ ] , would allow for a more robust analysis of hightemperature regions, as the current model is likely overstating the stability of graphenein these extreme conditions.Chapter 6 Malcolm Hillebrand 137haotic Dynamics of Polyatomic ModelsFurther investigations based on the results of this chapter would be to look at theeffects of geometric anharmonism on the chaoticity of 2D (and potentially 3D) lattices.Our fi nding that the harmonic interatomic coupling yields strongly chaotic behaviour,and that the addition of a simple angular potential is capable of stabilising the systemso much, raises questions about whether this is a particular property of the honeycomblattice, and thus a contributing factor to graphene’s extraordinary strength, or if it holdstrue in different 2D geometries. So performing estimations of the mLE and carryingout other stability computations across triangular, rectangular and hexagonal latticescould provide useful insights. Results obtained by using both harmonic and anharmoniccouplings could also be compared, to give an understanding of when, counterintuitively,a nonlinear potential produces less chaotic results.138 Malcolm Hillebrand Chapter 6 hapter 7Conclusions and Perspectives This thesis has been an exploration of the use of methods of chaotic dynamics in biolog-ical and chemical systems. The aim was to provide both new insights into the systemsthemselves and a deeper understanding of how chaotic dynamics, and particularly themaximum Lyapunov Exponent (mLE), can be used to study physical models. We haveexplored the nonlinear behaviour of the Peyrard-Bishop-Dauxois (PBD) model of DNA,investigating the effects of heterogeneity and temperature on the chaoticity, quanti fi edby the mLE. We have further probed the biological signi fi cance of DNA by studying theproperties of bubbles in the double strand in detail. Beyond DNA, we have performedan extensive study of the dynamical stability of graphene, again using the mLE as a quan-ti fi er for both large graphene sheets and fi nite width graphene nanoribbons (GNRs). Incombination, these studies have in turn led to a better grasp of the value that the mLEprovides in the study of biochemical models. Here we summarise the salient fi ndingsfrom these research efforts, as well as pointing towards a number of exciting possibleresearch avenues to be followed in the future.We made use of symplectic integration methods as well as parallel computing tech-niques to numerically simulate the dynamics of the Hamiltonian systems, with an em-phasis on computing the mLE as accurately and quickly as possible. The model ofDNA which forms the core focus of this thesis, the PBD model, uses a relatively sim-ple mesoscale potential approach to accurately approximate key dynamical features ofDNA, particularly melting curves. This allows the ef fi cient computation of macroscopicdynamical quantities, which is perfect for our application.The investigation of DNA is presented in three distinct parts, each corresponding toa chapter of this thesis. In Chapter 3, we introduced the notion of the alternation index α , which counts the number of times that the base pair type [ adenine-thymine (AT) orguanine-cytosine (GC) ] changes in the sequence. This value α quanti fi es the heterogene-ity of a given sequence by describing how well-mixed the base pairs are, whether thereis a homogeneous clumping of base pairs (small α ) or a heterogeneous mixture of basepairs in the strand (large α ). In order to understand which arrangements of base pairsare more or less likely, which would facilitate studying the effect of heterogeneity on thesystem’s dynamics, we found a probability distribution for α given a particular numberof AT and GC base pairs in a periodic DNA sequence. Although it was a challengingmathematical problem, we used Pólya counting theory and some tools from group the-ory to fi nd this distribution, which enabled us to construct an algorithm for computingthe number of possible arrangements of base pairs with a given number of alternations,139haotic Dynamics of Polyatomic Modelsand hence a complete distribution for all possible values of α . While Monte Carlo sim-ulations using around 20 000 runs yield a good approximation of the true distribution,the Pólya algorithm is more ef fi cient even for short sequences, and is substantially fasterfor long sequences. These distributions now allow the investigation of high and lowprobability regions, as well as understanding how probable the extreme cases are. Withan eye to future developments, while there is no obvious generalisation of this quantityto 2D lattices one could certainly fi nd an analagous measurement to quantify the het-erogeneity of a 2D lattice. The use of α and a 2D or 3D generalisation could be veryinformative for the study of heterogeneity in lattices.The next phase of our investigation was to compute the mLE of the considered DNAmodels and study their chaotic dynamics (Chapter 4). Particularly, our interest lay inthe heterogeneity; fi rstly the inherent heterogeneity arising from DNA sequences com-prising of both AT and GC base pairs in various ratios (which can be quanti fi ed byconsidering the percentage of AT or GC base pairs in the strand), and secondly the het-erogeneity within a set group of base pairs given by the clumping or mixing of basepairs (quanti fi ed by α ). The mLE was computed for sequences with varying AT and GCcontent, as well as across a spectrum of energy density (or equivalently, temperature) val-ues. In all cases, the mLE increases with energy, as the DNA molecule destabilises. Atlow energies, a higher concentration of AT base pairs results in more chaotic behaviour,while GC-rich sequences are more chaotic at high energies. Moving past the meltingpoint (and the physical validity of the model), we see that the mLE decreases to zero, be-having like a dynamical order parameter for the system. Looking at the effect of α , wefound that the more well-mixed heterogeneous sequences produce consistently more sta-ble (i.e. less chaotic) dynamics. Particularly for sequences with a roughly equal numberof AT and GC base pairs, more homogeneous molecules are substantially more chaoticthan the typical sequences. Beyond the global information gleaned from the mLE, thelocal chaoticity was investigated by computing the deviation vector distribution (DVD)for a number of sequences. For all cases, the DVD tends to concentrate near regions oflarge displacement, avoiding areas of very little motion, but also not focussing at highlydisplaced sites. In particular, this means that generally the DVD avoids regions with avery high or very low probability of bubble occurrence. Results for the PBD model andthe sequence-dependent ePBD model are found to be very similar with regard to theirchaoticity.The fi nal piece of our tripartite study of DNA was a thorough investigation in Chap-ter 5 of the properties of thermally induced openings, called bubbles. As a starting point,we introduced physically motivated threshold values for how widely separated a basepair needs to be in order to be considered open, or part of a bubble. A DNA strandis de fi ned to be melted when 50% of the base pairs are open, and using this criterionwe could fi nd very effective base pair-dependent thresholds, one for AT base pairs andone for GC. Using these thresholds, we found a remarkably precise melting tempera-ture for the ePBD model, and established a relationship between the temperature andenergy density for both the PBD and ePBD models. This relationship was fi tted witha simple cubic polynomial. Following this, and by using a large number of moleculardynamics simulations to create reliable statistics, we found distributions for the lengthsof bubbles in DNA sequences of varying base pair content at physiological tempera-tures. Using the base-pair-discrimating thresholds we introduced in our work meansthat probabilities for short bubbles (4 or fewer base pairs open) are similar for all se-140 Malcolm Hillebrand Chapter 7haotic Dynamics of Polyatomic Modelsquences regardless of AT / GC composition, but longer bubbles are more likely to occurin AT-rich sequences. The bubble length distribution was accurately fi tted by a stretchedexponential function, and the parameters provided in Section 5.3 (and also in [ ] ) al-low researchers to estimate the probability of a bubble of a given length occurring in anydesired sequence. Additionally, the ePBD model consistently predicts that more bubblesoccur than in the PBD model, particularly for large bubbles in AT-rich sequences. The fi nal, and most signi fi cant, result of our study was the determination of bubble lifetimedistributions for the two models, again at physiological temperatures. For a variety ofbase pair compositions and bubble lengths, the lifetime distributions have been accu-rately computed. These distributions were also fi tted with stretched exponential func-tions, which provide an excellent approximation for bubbles of two or more base pairs.Computing the mean bubble lifetime for each length reveals an exponential decrease inaverage lifetime as the bubble length increases, with longer lifetimes for sequences withmore AT content. The same trends were observed for the PBD and ePBD models, withgenerally longer lifetimes predicted using the sequence-dependent ePBD version of themodel. The general results found from simulating random sequences were comparedto simulations of the well-studied AdMLP sequence. While the ePBD model predictedslightly longer lifetimes in the AdMLP, consistent with its known propensity for largelong-lived bubbles, the distributions for random sequences and the AdMLP are very sim-ilar using both models, con fi rming that our results are applicable to speci fi c promoters,which holds great promise for aiding future investigations.In Chapter 6 we then considered a Hamiltonian model of graphene and studied thechaoticity of several graphene structures. Firstly, an in fi nite graphene sheet was mod-elled by applying periodic boundary conditions to all edges. There, in the same veinas the DNA, the mLE was computed for a spectrum of energy density values. In thegraphene however, there is no phase transition, and the mLE continues to grow quadrat-ically until the model loses meaning at extreme temperatures. Mild disorder was intro-duced in the graphene sheet by doping with C atoms, but even considering a sheetcomprised of purely C atoms the mLE was only very slightly affected. The mLE wasfound to be very small (orders of magnitude smaller than the mLEs computed for theDNA) and consequently the related Lyapunov time (the inverted mLE) revealed thatchaos takes a very long time to manifest in graphene when considered on the naturaltime scale of the system’s fastest normal modes. Several simpli fi cations of the nonlinearsystem were investigated, with the angular potential found to be a strong stabilising in- fl uence on the dynamics. Apart from this bulk behaviour, fi nite width GNRs were alsostudied, with the mLE computed as a function of ribbon width for several energies. Forboth zigzag and armchair ribbons the mLE was found to decrease as the width grows ac-cording to a decreasing Hill function towards the bulk value, with the armchair ribbonsbeing slightly more chaotic at very narrow widths.In all these investigations, and particularly the graphene simulations, computationaloptimisation played a fundamental role. The use of multithreading methods such asGNU Parallel and OpenMP, as well as optimisation techniques like code pro fi ling, haveenabled the pushing of numerical boundaries that have been previously impassable. Thehundredfold improvement in computational time from the fi rst implementation of thegraphene equations of motion to the fi nal version speaks to the value of devoting timeto writing ef fi cient code and optimising bottleneck routines. While this is not a majorresult of this thesis per se , computational optimisation is an increasingly important aspectChapter 7 Malcolm Hillebrand 141haotic Dynamics of Polyatomic Modelsof numerically studying dynamical systems. Particular aspects of the optimisation andparallelisation process were discussed in Chapter 1.Looking over the results for both systems, and the behaviour of the mLE in eachcase, we can make some more general observations. Firstly, it is gratifying that the mLEof the very stable graphene model is so much smaller than that of the soft DNA. Thispurely dynamical quanti fi cation of the systems’ stability matches neatly with graphene’sknown mechanical stability. Further work to consider a number of molecular com-pounds of varying mechanical strength (such as other proteins, or crystal lattices) withappropriate dynamical models and compute the mLE in each system could provide use-ful insight about this relationship. Studying the effects of different force models, suchas the Lennard-Jones interatomic potential [ ] or the anisotropic Gay-Berne potential [ ] , could also prove interesting. We also saw that in the DNA as the energy increases,openings in the strand cause local linearisation of the system due to the form of the inter-molecular Morse potential. These linearisations lead to a slowing down of the increasein the mLE with energy, giving a fl attening behaviour as melting is approached. Dueto the complete lack of this behaviour in the 2D graphene sheet, the mLE continues togrow with no changing behaviour. In this way the mLE describes the overall dynam-ics of the systems effectively across the energy spectrum. As energy increases, therewill be an inevitable decrease in stability due to stronger atomic fl uctuations, but in theDNA system the mLE captures the process of denaturation which does not occur in thegraphene lattice. Further investigation of this idea would be interesting – for instance,if a graphene sheet is placed under strain would there be some characteristic change inthe behaviour of the mLE as the breaking point is approached? Or in other systemsthat exhibit local linearisations of the potentials, do we see behaviours in the mLE likethat of the DNA? The capturing of this global dynamical behaviour through the mLEcould be a very promising avenue of research, which could lead to better understandingof structural stability as well as dynamical stability in a variety of compounds.Another strong possibility for future work is to look more closely at the localisa-tion of strong nonlinearity. This is more relevant in the heterogeneous DNA system,where particular regions of the sequence are of biological interest, but could also applyto graphene and related models. The behaviour of the DVD in the DNA suggests that itcould provide information about the formation of bubbles or other physical phenom-ena. Further study of the DVD, as well as local stretching of deviations and covariantLyapunov vectors would be very interesting based on the results we have found so far.The contribution to the overall stability of graphene from the angular potential is an-other potentially signi fi cant point, particularly linked to the appeareance of chaos in thepurely harmonically coupled lattice. Since any 2D lattice with harmonic coupling thatallows atoms full 2D range of motion will exhibit chaotic motion, it could be instruc-tive to further investigate how this chaoticity depends on the geometry (e.g. hexagonal,triangular, rectangular or Kagome lattices), as well as how a potential well for the anglesbetween atoms affects the stability in each case. Since there are many physical materialswith varying lattice structures, this could in turn inform solid state experimentation.In all numerical studies, there is always a limitation that the results can only be as ac-curate as the model studied allows by its faithful representation of the physical system.In our case, the models used have been rigorously tested and used extensively, givingstrong evidence that our numerical simulations are applicable to the true systems. How-ever, the exploration of other models, which may capture different dynamical aspects of142 Malcolm Hillebrand Chapter 7haotic Dynamics of Polyatomic Modelsthe real systems more accurately, could always provide deeper understanding of thesesystems. For instance, DNA models that include twisting force terms [ ] , differ-ent statistical mechanics approaches [ ] , or other ideas (see e.g. [ ] and referencestherein) would be worth studying in more detail.The logical next step regarding the complexity of the graphene model is to considerout of plane torsional effects, rather than con fi ning all motion to a 2D plane. The useof a model including such terms (such as in [ ] ) would allow a more realistic investi-gation of higher-energy regions, as well as provide insight to the signi fi cance of this outof plane motion. While this addition would make the equations of motion much morecomplex, and the computational expense greater, it would be a very useful step towardsa dynamical understanding of graphene.In conclusion, we note that our research has provided new and signi fi cant knowl-edge about the dynamics and stability of DNA and graphene, along with additionalevidence of the usefulness of nonlinear dynamics methods (like the mLE and the DVD)for studying complex physical systems. This work has laid the foundation for a numberof exciting new prospects, which could be in fl uential in a broad array of applied sciences.Chapter 7 Malcolm Hillebrand 143 ibliography [ ] C. Calladine and H. Drew,
Understanding DNA (Academic Press, 2004). [ ] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos,I. V. Grigorieva, and A. A. Firsov, Science , 666 (2004). [ ] H. Goldstein, C. P. J. Poole, and J. L. Safko,
Classical mechanics (3rd ed.) (AddisonWesley, 2002). [ ] S. Strogatz,
Nonlinear dynamics and chaos (Boca Raton: CRC Press, 2015). [ ] R. L. Devaney,
An introduction to chaotic dynamical systems (2nd ed.) (WestviewPress, 2003). [ ] E. Lorenz, J. Atmos. Sci. (1963). [ ] S. N. Rasband,
Chaotic dynamics of nonlinear systems (Wiley, 1990). [ ] M. Henon, Physica D , 412 (1982). [ ] P. A. Patsis and L. Zachilas, Int. J. Bifurcat. Chaos , 1399 (1994). [ ] M. Richter, S. Lange, A. Bäcker, and R. Ketzmerick, Phys. Rev. E , 022902(2014). [ ] C. Skokos, Journ. Phys. A , 10029 (2001). [ ] C. Skokos, T. C. Bountis, and C. Antonopoulos, Physica D , 30 (2007). [ ] C. Skokos and T. Manos, Lect. Notes Phys. , 129 (2016). [ ] C. Skokos, G. A. Gottwald, and J. Laskar, Lect. Notes Phys. (2016). [ ] R. De Vogelaere, Technical Report, University of Notre Dame
No. 4 (1956). [ ] E. Hairer, C. Lubich, and G. Wanner,
Geometric numerical integration (Springer-Verlag Berlin Heidelberg, 2006). [ ] E. Forest and R. D. Ruth, Physica D (1990). [ ] H. Yoshida, Phys. Lett. A (1990). [ ] R. I. McLachlan and G. R. W. Quispel, Acta Numerica , 341 (2002). [ ] E. Forest, J. Math. Phys. , 1133 (1990). [ ] S. Blanes, F. Casas, and A. Murua, Bol. Soc. Esp. Mat. Apl. , 89 (2008). [ ] S. Blanes, F. Casas, A. Farrés, J. Laskar, J. Makazaga, and A. Murua, App. Num. Math. , 58 (2013). [ ] B. Senyange, B. M. Manda, and C. Skokos, Phys. Rev. E , 052229 (2018). [ ] S. Blanes and P. Moan, Journ. Comp. App. Math. , 313 (2002).144haotic Dynamics of Polyatomic Models [ ] Centre for High Performance Computing, (2020) https://chpc.ac.za . [ ] O. Tange,
GNU parallel 2018 (Ole Tange), p. 1. [ ] B. Chapman, R. Van De Pas, and G. Jost,
Using OpenMP: portable shared memoryparallel programming (MIT Press, 2007). [ ] J. Sanders and E. Kandrot,
CUDA by example: an introduction to general-purposeGPU programming (Pearson, 2010). [ ] C. Skokos, “The Lyapunov characteristic exponents and their computation”, in
Dynamics of small solar system bodies and exoplanets (Springer Berlin Heidelberg,Berlin, Heidelberg, 2010), pp. 63–135. [ ] A. M. Lyapunov, Int. J. Cont. , 531 (1992). [ ] V. I. Oseledec, Trans. Moscow Math. Soc. , 197 (1968). [ ] G. Benettin, L. Galgani, and J.-M. Strelcyn, Phys. Rev. A , 2338 (1976). [ ] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Meccanica (1980). [ ] C. Skokos, I. Gkolias, and S. Flach, Phys. Rev. Lett. , 064101 (2013). [ ] B. Many Manda, B. Senyange, and C. Skokos, Phys. Rev. E , 032206 (2020). [ ] G. Contopoulos,
Order and chaos in dynamical astronomy (Springer, 2002). [ ] L. Mei and L. Huang, Comp. Phys. Comm. , 108 (2018). [ ] A. Bountis, K. Kaloudis, T. Oikonomou, B. Many Manda, and C. Skokos, Int. J. Bi-furc. Chaos, In Press (2020). [ ] G. Contopoulos, L. Galgani, and A. Giorgilli, Phys. Rev. A (1978). [ ] L. Miranda Filho, M. Amato, Y. Elskens, and T. Rocha Filho, Comm. NonlinearSci. , 236 (2019). [ ] A. Ngapasare, G. Theocharis, O. Richoux, C. Skokos, and V. Achilleos, Phys.Rev. E , 032211 (2019). [ ] J. Watson, T. Baker, S. Bell, A. Gann, M. Levine, and R. Losick,
Molecular biologyof the gene, fi fth edition (Cold Spring Harbor Laboratory Press, 2003). [ ] W. Saenger,
Principles of nucleic acid structure (Springer-Verlag, 1984). [ ] M. Peyrard and A. R. Bishop, Phys. Rev. Lett. , 2755 (1989). [ ] T. Dauxois, M. Peyrard, and A. R. Bishop, Phys. Rev. E , 684 (1993). [ ] T. Dauxois, M. Peyrard, and A. R. Bishop, Phys. Rev. E , R44 (1993). [ ] M. Peyrard, Nonlinearity , R1 (2004). [ ] J. D. Watson and F. H. C. Crick, Nature , 737 (1953). [ ] T. Dauxois, Phys. Lett. A , 390 (1991). [ ] C. H. Choi, G. Kalosakas, K. Ø. Rasmussen, M. Hiromura, A. R. Bishop, andA. Usheva, Nucl. Ac. Res. , 1584 (2004). [ ] S. Ares, N. K. Voulgarakis, K. Ø. Rasmussen, and A. R. Bishop, Phys. Rev. Lett. , 035504 (2005). [ ] G. Kalosakas, K. Ø. Rasmussen, and A. R. Bishop, Chem. Phys. Lett. , 291(2006).Chapter 7 Malcolm Hillebrand 145haotic Dynamics of Polyatomic Models [ ] T. S. van Erp, S. Cuesta-López, and M. Peyrard, Eur. Phys. J. E , 421 (2006). [ ] B. S. Alexandrov, V. Gelev, S. W. Yoo, L. B. Alexandrov, Y. Fukuyo, A. R.Bishop, K. Ø. Rasmussen, and A. Usheva, Nucl. Ac. Res. , 1790 (2009). [ ] R. Thomas, Biochim. Biophys. Acta , 231 (1954). [ ] R. M. Wartell and A. S. Benight, Phys. Rep. , 67 (1985). [ ] R. B. Inman and R. L. Baldwin, J. Mol. Biol. , 172 (1962). [ ] D. W. Ussery, in
Encyclopedia of genetics , edited by S. Brenner and J. H. Miller(Academic Press, New York, 2001), pp. 550–553. [ ] N. C. Harris and C.-H. Kiang, J. Chem. Phys. B , 16393 (2006). [ ] G. Kalosakas and S. Ares, J. Chem. Phys. , 235104 (2009). [ ] M. Peyrard and T. Dauxois, Math. Comput. Simul. , 305 (1996). [ ] C. W. Hsu, M. Fyta, G. Lakatos, S. Melchionna, and E. Kaxiras, J. Chem. Phys. , 105102 (2012). [ ] N. Theodorakopoulos, Proc. Int. Conf. on Localization and Energy Transfer inNonlinear Systems (Escorial,Spain) (2002). [ ] S. W. Englander, N. R. Kallenbach, A. J. Heeger, J. A. Krumhansl, and S. Litwin,P. Natl. Acad. Sci. , 7222 (1980). [ ] S. Yomosa, Phys. Rev. A , 2120 (1983). [ ] Y. Gao and E. W. Prohofsky, J. Chem. Phys. , 2242 (1984). [ ] D. Poland and H. A. Scheraga, J. Chem. Phys. (1966). [ ] R. D. Blake, J. W. Bizzaro, J. D. Blake, G. R. Day, S. G. Delcourt, J. Knowles,K. A. Marx, and J. SantaLucia J, Bioinformatics , 370 (1999). [ ] R. Blossey and E. Carlon, Phys. Rev. E , 061911 (2003). [ ] M. E. Fisher, J. Chem. Phys. , 1469 (1966). [ ] E. W. Prohofsky, Phys. Rev. A , 1538 (1988). [ ] P. M. Morse, Phys. Rev. , 57 (1929). [ ] L. V. Yakushevich,
Nonlinear physics of DNA (John Wiley and Sons, Ltd, 1998). [ ] A. Campa and A. Giansanti, Phys. Rev. E , 3585 (1998). [ ] B. S. Alexandrov, V. Gelev, Y. Monisova, L. B. Alexandrov, A. R. Bishop, K.Rasmussen, and A. Usheva, Nucl. Ac. Res. , 2405 (2009). [ ] T. Dauxois and M. Peyrard, Phys. Rev. E , 4027 (1995). [ ] S. G. Delcourt and R. D. Blake, J. Biol. Chem. , 15160 (1991). [ ] J. Šponer, P. Jure ˇ cka, I. Marchan, F. J. Luque, M. Orozco, and P. Hobza, Chem.-Eur. J. , 2854 (2006). [ ] M. Hillebrand, G. Paterson-Jones, G. Kalosakas, and C. Skokos, Regul. Chaot. Dy-nam. , 135 (2018). [ ] M. Régnier, Discrete App. Math. , 259 (2000). [ ] S. Robin and J. J. Daudin, J. Appl. Probab. , 179 (1999).146 Malcolm Hillebrand Chapter 7haotic Dynamics of Polyatomic Models [ ] S. Robin and S. Schbath, J. Comp. Biol. , 349 (2001). [ ] G. Pólya and R. C. Read,
Combinatorial enumeration of groups, graphs, and chem-ical compounds (New York: Springer, 1987). [ ] N. Jacobson,
Basic Algebra I, second edition (Dover Publications, 2009). [ ] I. N. Herstein,
Abstract algebra, third edition (New York: Wiley, 1999). [ ] R. A. Brualdi, “Pólya counting”, in (Prentice Hall, 2010), pp. 541–581. [ ] J. H. van Lint and R. M. Wilson,
A course in combinatorics (Cambridge Univer-sity Press, 2001). [ ] W. Burnside,
Theory of groups of fi nite order (Cambridge University Press, 1897). [ ] The Sage Developers,
Sagemath, the sage mathematics software system (version 9.0) , (2020). [ ] O. Zariski and P. Samuel,
Commutative algebra, vol. II (Springer-Verlag BerlinHeidelberg, 1960). [ ] M. Hillebrand, G. Kalosakas, A. Schwellnus, and C. Skokos, Phys. Rev. E ,022213 (2019). [ ] M. Peyrard and J. Farago, Physica A , 199 (2000). [ ] G. Kalosakas, K. Ø. Rasmussen, and A. R. Bishop, J. Chem. Phys. , 3731(2003). [ ] G. Kalosakas, Phys. Rev. E , 051905 (2011). [ ] L. van Hove, Physica , 137 (1950). [ ] L. D. Landau and E. M. Lifshitz,
Statistical physics, third edition (Pergamon Press,Oxford, 1980). [ ] N. Theodorakopoulos, T. Dauxois, and M. Peyrard, Phys. Rev. Lett. , 6 (2000). [ ] J. A. Cuesta and A. Sánchez, J. Stat. Phys. , 869 (2004). [ ] J. Barre and T. Dauxois, Europhys. Lett. , 164 (2001). [ ] A. Apostolaki and G. Kalosakas, Physical Biology , 026006 (2011). [ ] C. H. Choi, Z. Rapti, V. Gelev, M. R. Hacker, B. Alexandrov, E. J. Park, J. S.Park, N. Horikoshi, A. Smerzi, K. Ø. Rasmussen, A. R. Bishop, and A. Usheva,Biophys. J. , 597 (2008). [ ] M. Zoli, J. Chem. Phys. , 115101 (2011). [ ] F. Ginelli, H. Chaté, R. Livi, and A. Politi, J. Phys. A.-Math.-Theor. , 254005(2013). [ ] H. M. Sobell, P. Natl Acad. Sci , 5328 (1985). [ ] G. Kalosakas, K. Ø. Rasmussen, A. R. Bishop, C. H. Choi, and A. Usheva, Eu-rophys. Lett. , 127 (2004). [ ] M. Hillebrand, G. Kalosakas, C. Skokos, and A. Bishop, Phys. Rev. E , 062114(2020). [ ] B. S. Alexandrov, V. Gelev, S. W. Yoo, A. R. Bishop, K. Ø. Rasmussen, and A.Usheva, PLoS Comput. Biol. , 1 (2009).Chapter 7 Malcolm Hillebrand 147haotic Dynamics of Polyatomic Models [ ] R. Tapia-Rojo, D. Prada-Gracia, J. Mazo, and F. Falo, Phys. Rev. E , 021908(2012). [ ] H.-H. Huang and P. Lindblad, J. Biol. Eng. , 10 (2013). [ ] L. B. Nowak-Lovato K. and. Alexandrov, A. Banisadr, A. Bauer, A. Bishop, A.Usheva, F. Mu, E. Hong-Geller, K. Rasmussen, W. Hlavacek, and B. Alexandrov,PLoS Comput. Biol. , e1002881 (2013). [ ] R. Tapia-Rojo, J. Mazo, J. Hernandez, M. Peleato, M. Fillat, and F. Falo, PLoSComput. Biol. , e1003835 (2014). [ ] M. Allen and D. Tildesley,
Computer simulation of liquids: Second edition (OxfordScience, Nov. 2017), pp. 1–626. [ ] S. Ares and G. Kalosakas, Nano Letters , 307 (2007). [ ] D. Hennig, E. Starikov, J. Archilla, and F. Palmero, J. Biol. Phys , 227 (2004). [ ] P. Maniadis, G. Kalosakas, R. K.Ø., and A. Bishop, Phys. Rev. E , 021912(2005). [ ] G. Kalosakas, K. Ngai, and S. Flach, Phys. Rev. E , 061901 (2005). [ ] A. Chetverikov, W. Ebeling, V. Lakhno, A. Shigaev, and M. Velarde, Eur. Phys. J. B , 101 (2016). [ ] L. Gu and H.-H. Fu, New. J. Phys. , 053032 (2016). [ ] E. Diaz, R. Lima, and F. Dominguez-Adame, Phys. Rev. B , 134303. [ ] A. Chetverikov, W. Ebeling, V. Lakhno, and M. Velarde, Phys. Rev. E , 0522203(2019). [ ] D. Chen, S. Aubry, and G. Tsironis, Phys. Rev. Lett , 4776 (1996). [ ] J. Marmur and P. Doty, J. Mol. Biol. , 109 (1962). [ ] M. Hillebrand, B. Many Manda, G. Kalosakas, E. Gerlach, and C. Skokos, Chaos , 063150 (2020). [ ] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim,Rev. Mod. Phys. , 109 (2009). [ ] D. R. Cooper and et al, ISRN Condensed Matter Physics , 501686 (2012). [ ] C. Liu, Z. Yu, D. Neff, A. Zhamu, and B. Jang, Nano Letters , 4863 (2010). [ ] Y.-M. Lin, A. Valdes-Garcia, S.-J. Han, D. B. Farmer, I. Meric, Y. Sun, Y. Wu,C. Dimitrakopoulos, A. Grill, P. Avouris, and K. A. Jenkins, Science , 1294(2011). [ ] J. S. Bunch, A. M. van der Zande, S. S. Verbridge, I. W. Frank, D. M. Tanenbaum,J. M. Parpia, H. G. Craighead, and P. L. McEuen, Science , 490 (2007). [ ] F. Schedin, A. K. Gaim, S. V. Morozov, E. W. Hill, B. P, M. I. Katselnslson, andK. S. Novoselov, Nat Mater. , 652 (2007). [ ] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V.Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature , 197 (2005). [ ] K. Bolotin, K. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, and H.Stormer, Solid State Commun. , 351 (2008).148 Malcolm Hillebrand Chapter 7haotic Dynamics of Polyatomic Models [ ] S. Ghosh, I. Calizo, D. Teweldebrhan, E. P. Pokatilov, D. L. Nika, A. A. Ba-landin, W. Bao, F. Miao, and C. N. Lau, Appl. Phys. Lett. , 151911 (2008). [ ] C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science , 385 (2008). [ ] G. Tsoukleri, J. Parthenios, K. Papagelis, R. Jalil, A. C. Ferrari, A. K. Geim, K. S.Novoselov, and C. Galiotis, Small , 2397 (2009). [ ] A. Xu and M. J. Buehler, J. Phys.-Condens. Mat. , 485301 (2010). [ ] Z. Guo, D. Zhang, and X.-G. Gong, Appl. Phys. Lett. , 163103 (2009). [ ] A. Sgouros, M. M. Sigalas, K. Papagelis, and G. Kalosakas, J. Phys.: Condens. Mat-ter , 125301 (2014). [ ] L. Zhang, X. Zeng, and X. Wang, Sci. Rep. , 3162 (2013). [ ] A. P. Sgouros, G. Kalosakas, M. M. Sigalas, and K. Papagelis, RSC Adv. , 39930(2015). [ ] Y. Huang, J. Wu, and K. C. Hwang, Phys. Rev. B , 245413 (2006). [ ] K. V. Zakharchenko, M. I. Katsnelson, and A. Fasolino, Phys. Rev. Lett. ,046808 (2009). [ ] E. Cadelano, P. Palla, S. Giordano, and L. Colombo, Phys. Rev. Lett. , 235502(2009). [ ] C. Reddy, S. Rajendran, and K. Liew, Nanotechnology , 864 (2006). [ ] F. Liu, P. Ming, and J. Li, Phys. Rev. B , 064120 (2007). [ ] G. Kalosakas, N. N. Lathiotakis, C. Galiotis, and K. Papagelis, J. Appl. Phys. , 134307 (2013). [ ] A. Ngapasare, “Dynamical behavior of graphene models”, MSc Thesis (Univer-sity of Cape Town, South Africa, 2017). [ ] B. Senyange and C. Skokos, Eur. Phys. J. Spec .Top. , 625 (2018). [ ] D. Frenkel and B. Smit, in
Understanding molecular simulation (second edition) ,edited by D. Frenkel and B. Smit, Second Edition (Academic Press, San Diego,2002), pp. 63–107. [ ] Z. G. Fthenakis, G. Kalosakas, G. D. Chatzidakis, C. Galiotis, K. Papagelis, andN. N. Lathiotakis, Phys. Chem. Chem. Phys. , 30925 (2017). [ ] J. Jones and S. Chapman, P. R. Soc. Lond. A , 441 (1924). [ ] J. G. Gay and B. J. Berne, J. Chem. Phys. , 3316 (1981). [ ] M. Barbi, S. Cocco, and M. Peyrard, Phys. Lett. A , 358 (1999). [ ] G. Weber, N. Haslam, N. Whiteford, A. Prugel-Bennett, J. Essex, and C. Neylon,Nature Phys. , 55 (2006). [ ] M. Manghi and N. Destainville, Phys. Rep.631