Thermodynamic analysis of a non linear lattice

Thermodynamic analysis of a non linear lattice

R. A. S. Silva^*; E. Drigo Filho; J. R. Ruggiero

UNESP, Campus de São José do Rio Preto, Departamento de Física Rua Cristóvão Colombo, 2265, CEP 15054-000 São José do Rio Preto - SP, Brazil

ABSTRACT

We investigate an extension of the Perard and Bishop (PB) model. In the studied model we get a harmonic one dimensional lattice chain with an additional Morse potential on site. The rotation and vibration motion of each component of the lattice are considered and the coupling for these motions is introduced by a resonance condition. Thermodynamics and structural properties of the system are explored.

Keywords: Nonlinear lattice; Peyrard-Bishop model; Morse Potential; Radius of gyration

1. INTRODUCTION

One dimensional non linear lattice has a long history beginning with the Fermi-Pasta-Ulam (FPU). After that, several works were dedicated to analyze these kind of system. Usually, the dynamical aspects of the problems are studied (for a recent review, see ref.[1]). However, Peyrard and Bishop (PB) had introduced a model [2] to simulate the DNA macromolecule putting attention on thermo-dynamical aspects of the lattices.

In the FPU system the lattice is composed by masses linked by non-linear potentials (the original potential had the form V(x) = ax² + bx⁴). In the PB model, on the other hand, the DNA macromolecule is mimicked by two spring-masses lines with an additional non linear potential joining the two lines at the each adjacent pair of masses, representing H-bonds between the bases pairs. The mathematical treatment of this model leaves to a harmonic lattice with a non linear potential on site. The most frequently potential function in the context is the Morse potential. Some extensions of the original model in [2] are proposed considering, for example, the helicoidal structure of DNA [3].

An important objective of the PB model is to understand the phenomenon of DNA thermal denaturation and gets some knowledge about others fundamental processes involving local denaturation as, for example, the genetic transcription and drugs intercalation. From this type of model is possible to obtain interesting properties, as the average stretching between base pairs as function of the temperature using the transfer integral operator [4]. It is important to note that the original PB model consider only transversal motion of the lattices components, i.e. only vibrational motions are studied. In this work, we proposed an extension of the original PB model [2] by introducing rotational motions for the nucleotides. In this way, both the vibrational and rotational motion of each mass on the lattice is considered. As in the original model, the stretch of the base pairs is driven by the Morse potential and the coupling between the two kinds of motion is obtained through resonance condition. This condition leave us to a geometric link on the model which introduces the radius of gyration of the lattice components.

We consider that the proposed model is still too simple to simulate the DNA macromolecule. However, it permits a precise mathematical analyzes. Then, it is a good theoretical laboratory to study generalizations of the PB model. Following this line of thinking, all the parameters used in this work for the numerical calculations are related with the phenomenological DNA analyzes.

2. THE MODEL

In the present model each nucleotide is represented by a disk of mass m and radius R. These disks are joined together along the ribbon by spring potentials with coupling constants k. The two ribbons are coupled to each other by a non linear potential simulating the hydrogen bonds. Each disk can have two kinds of displacements, one linear (denoted by u and v for each one of the adjacent ribbons) and an angular one, denoted by θ_n and ø_n, respectively. For simplicity, we assumed that the rotations happen around the axis that are parallel to the helix axis, i.e. perpendicular to the hydrogen bond plane (see figure 1). The linear and rotations displacements are all in the same plane. We couple the two ribbons through the V_n potential, V(u_n, v_n, θ_n, ø_n). The dependence in the vibrational coordinates follows the normal PB model [2,5] and the dependence from θ_n and ø_n follows rotational models as, for example, the coordinates cited in the ref. [6,7]. In this way, we suppose that V(u_n, v_n, θ_n, ø_n) = V(u_n - v_n, θ_n - ø_n).

The Hamiltonian for the system can be written as

where

represents the kinetics energy of system, and

is the potential energy. In Eq. (2), and are the linear moments, and are angular moments. For simplicity, the system is considered homogenous, which means that the mass m and the moment of inertia I are equals for all nucleotides. The same hypothesis is assumed by the strength constants k and the angular strength constants x. The Fig. 1 shows a representation of the mechanical model indicated by Eq. (1). Introducing a new set of coordinates:

and

we decouple the x from y motions and α from β ones in the original Hamiltonian (1). The interaction potential V_n, however, remains dependent on a combination of y_n and α_n.

For our purpose, we admit that the potential depends on a linear combination of the variables y_n and α_n, i.e. V(y_n + sα_n ), where s is a parameter coupling the coordinate y with the angular motion described by α. This kind of dependence in the potential is adopted as a first approximation for the real dependence and the parameter s is related to the radius of gyration as shown below. Introducing a canonical transformation for the angular coordinate (α_n; ), _n = sα_n and = /s, we get a new form to the Hamiltonian,

For decoupled system we can identify the natural frequencies for the vibration motion and rotation, as ω_V = 2π(k/m)^(1/2) and ω_R = 2π(ξ/I)^(1/2), respectively. In order to join the angular and the linear motion and uncouple the coordinates in the Hamiltonian, we impose a resonance condition, i.e. ω_V = ω_R. This condition is supported by our main interest, which is to observe the thermodynamical properties of the system, in particular the rupture of the potential (denaturation). In this context, we remember that the oscillation amplitudes increase drastically near the denaturation temperature. Then, it is important to note that the amplitude of motions becomes larger in the resonance condition and the denaturation should be more probable. An additional argument to use this supposition is that the coupling between rotational and vibrational motion became more evident in the resonance condition once in this condition occurs the greatest energy exchange. Therefore, we obtain ξ = ks² when considering I = ms². These constrains permit to identify s as the gyration radius and to connect the parameters ξ and k. Thus, defining a new set of coordinates

the Hamiltonian (6) becomes

As in the PB model, the interaction between two adjacent bases is given by the Morse potential [2,5]

where D is the energy of dissociation of the pair and a is a parameter with inverse of length dimension. This potential is generally chosen to simulate hydrogen bonds in small molecules as well [8,9].

For a chain with N pairs of bases, the classical partition function is given in terms of the Hamiltonian (8), and can be factored as

The integrals in all moments are simple gaussian integrals and gives

and

where k_B is the Boltzmann constant. The Z_x, Z_w and Z_b are contributions from harmonic chain of the oscillator. It is an easy task to calculate analytically these integrals as simple gaussian integrals and the result is given by

and

The remaining term, in the λ variables, comes from an harmonic chain of oscillator with an on-site non linear potential for each oscillator. The Z_λ, partition function can be written as

where F is given by

The integral in Z_λ can be calculated using eigenfunctions and eigenvalues of a transfer integral operator [4]

In the thermodynamic limit (N → ∞), the partition function Z_λ, reduces to

where ε₀ is the ground state eigenvalue for a Schrödinger-like equation given by [4]:

where

The eigenvalue and the normalized eigenfunction for the ground state are, respectively,

and

with d = (kD)^1/2/(k_BTa) and C the normalization constant. In the limit N → ∞ the average stretching < λ > will be dependent only on the ground state eigenfunction as

A criterion to determine the probability of rupture of the non linear Morse potential as function of temperature is to consider that for a distance smaller than 2Å the potential is intact. This criterion is based in the phenomenological knowledge that the hydrogen bond in the DNA is broken when the distances between the basis pair are bigger than two angstroms [5]. Then, the calculation of the probability is done using the follow equation

3. RESULTS

The qualitative behavior of the average stretching < λ > as function of the temperature is obtained from Eq.(23) and shown in Fig 2 while the fraction of closed base pairs versus temperature is obtained from Eq.(24), see Fig. 3. The results are similar to those obtained for a pure vibrational model [5]. However, in order to obtain a qualitative realistic result it is necessary to fix the five parameters introduced in the model, i.e. k, ξ, a, D and s.

The decision about the best values for the parameters is based in experimental results for DNA presented in the literature as, for instance, scattering of neutrons [7], single-molecule experiments [9] and low frequencies modes of vibration [10]. In this way, we obtain the limit Da/2 < 75 pN [7] and ξ = 2 eV [7]. An important parameter introduced in this work is the radius of gyration s. It appears linked to the moment of inertia of each disc. A good approximation to the radius of gyration is s = 4 Å obtained from crystallographic results of the nucleotides, contained in the Protein Data Bank (www.rcsb.org/pdb/home/home.do). The nucleotide is understood as a nitrogen base (Thymine, Adenine, Guanine and Cytosine), the pentose sugar and the phosphate group, as usual, see, for example, ref. [11]. Following the above criterions we obtain an optimal set of parameters as: ξ = 2 eV, k = ξ/s² = 0.125 eV Å^-2, a = 2.22 Å^-1 and D = 0.018 eV . From these values, we note that Da/2 ≈ 32 pN < 75 pN [9]. Using this set of parameters the model became completely defined. In this case, the resulting denaturation temperature is T_D ≈ 350 K (solid curve, in Figs. 2 and 3).

The denaturation temperature obtained is compatible with the experimental results for DNA observed, for example, in absorbance UV light at 260 nm [5]. In Fig 2 and 3 we observe the behavior of denaturation curves for different values for the radius of gyration with fixed values of ξ, a and D. We observe that the melting temperature decrease when s increase.

4. CONCLUSION

The extension of the PB model discussed in this work permits a complete mathematical treatment in analytical way. It also permits to introduce structural factor for the analyzed lattice, i.e. the radius of gyration. The model leave us to determine the melting temperature for the used parameters, ξ = 2 eV , k = 0.125 eV Å^-2, a =2.22Å^-1, D = 0.018 eV and s = 4 Å. Some variations in the value of the radius of gyration are also shown in Fig. 2 and 3.

The results are quite sensitive to this value and we observe that the denaturation temperature increase if the gyration radius decrease. However, it is possible an adjustment of the temperature of denaturation by changing the parameters a and D inside the range 0 < Da/2 < 75 pN. The result present in the figures 2 and 3 show a strong correlation between vibration and rotation motion and they indicate the necessity of more accuracy experimental results in order to fix the parameters, particularly, the Morse potential parameters a and D.

The introduction of rotation motion makes the model closer to the DNA real system. However, the presented model is still very simple to simulate this macromolecule. In particular, the assumption that there is a linear dependence between the vibration variable (y) and the rotational one (α) looks like artificial for a complex system as the real DNA. The model can be improved by using a more realistic connection between these two kinds of variable.

Acknowledgments

This work was partially supported by FAPESP, CNPq, FUNDUNESP and CAPES (Brazilian agencies). EDF wishes to thank the Isaac Newton Institute for Mathematical Sciences, University of Cambridge, for support during the programme on Combinatorics and Statistical Mechanics (January-June 2008), where this work was discussed and completed.

Received on 23 July, 2008

• [1] D. K. Cambell, P. Rosenau, and G. M. Zaslavsky, Chaos 15, 015101 (2005). In the same volume there are several works about this subject.
• [2] M. Peyrard and A. R. Bishop, Phys. Rev. Lett. 62, 2755 (1989).
• [3] S. Zdravkovic, M. Satarić and J. Tuszyński, Biophysical Implications of the Peyrard-Bishop-Dauxois Model of DNA Dynamics. J. Comp. Theo. Nanosc. 1, 171 (2004).
• [4] D. J. Scalapino, M. Sears, and R. A. Ferrell, Phys. Rev. B 6, 3409 (1972).
• [5] M. Peyrard, Nonlinearity 17, r1-r40 (2004).
• [6] L. V. Yakushevich, Phys. Lett. A 136, n.7, 413 (1989).
• [7] V. K. Fedyanin and L. V. Yakushevich, Stud. Biophys. 103, 171 (1984).
• [8] Y. Gao, K. V. Devi-Prosad, and E. W. Prohofsky, J. Chem. Phys. 80, 6291 (1984).
• [9] S. Zdravković and M. V. Satarić, Single-molecule unzippering experiments on DNA and Peyrard-Bishop-Dauxois model. Phys. Rev. E 73, 021905(1-11) (2006).
• [10] E. Drigo Filho and J. R. Ruggiero, Phys. Rev. A 44, n.12, 43 (1991).
• [11] W. Saenger, Principles of Nucleic Acid Structure, Springer Verlag New York (1984).

Publication Dates

• Publication in this collection
  22 Sept 2008
• Date of issue
  Sept 2008

History

• Received
  23 July 2008