# Analytical tools for solitons and periodic waves corresponding to phonons on Lennard-Jones lattices in helical proteins

Francesco d'Ovidio, Henrik Georg Bohr, Per-Anker Lindgård

► To cite this version:

Francesco d'Ovidio, Henrik Georg Bohr, Per-Anker Lindgård. Analytical tools for solitons and periodic waves corresponding to phonons on Lennard-Jones lattices in helical proteins. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2005, 71, pp.26606. 10.1103/PHYSREVE.71.026606. hal-00770747

HAL Id: hal-00770747

https://hal.science/hal-00770747

Submitted on 17 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Analytical tools for solitons and periodic waves corresponding to phonons on Lennard-Jones lattices in helical proteins

Francesco d'Ovidio*

Mediterranean Institute of Advanced Studies IMEDEA (CSIC-UIB), Campus Universitat Illes Balears, E-07122 Palma de Mallorca, Spain

Henrik Georg Bohr†

Centre for Quantum Protein, Department of Physics, Technical University of Denmark, DK 2800 Lyngby, Denmark

Per-Anker Lindgård‡

Materials Research Department, Risø National Laboratory, 4000 Roskilde, Denmark
and Centre for Quantum Protein,

Department of Physics, Technical University of Denmark, DK 2800 Lyngby, Denmark

(Received 20 May 2004; published 15 February 2005)

We study the propagation of solitons along the hydrogen bonds of an α helix. Modeling the hydrogen and peptide bonds with Lennard-Jones potentials, we show that the solitons can appear spontaneously and have long lifetimes. Remarkably, even if no explicit solution is known for the Lennard-Jones potential, the solitons can be characterized analytically with a good quantitative agreement using formulas for a Toda potential with parameters fitted to the Lennard-Jones potential. We also discuss and show the robustness of the family of periodic solutions called cnoidal waves, corresponding to phonons. The soliton phenomena described in the simulations of α helices may help to explain recent x-ray experiments on long α helices in Rhodopsin where a long lifetime of the vibrational modes has been observed.

DOI: 10.1103/PhysRevE.71.026606

PACS number(s): 05.45.Yv, 87.15.—v

## I. INTRODUCTION

The structure in most proteins consists of α helices and the β sheets connected by loops. The correctly folded structure is vital for the protein function [1]. One interesting, but not yet well understood, function is that of a nanomachine which is capable of transforming chemical and radiation energy into mechanical motion and can transport it directionally and without loss over large distances [2]. Considerable effort has been spent on studying the high-frequency amide I mode (C=O stretch at 1600 cm⁻¹), which was proposed by Davydov [3] as a candidate for soliton transport. In general, simplified one-dimensional (1D) models have been studied from a mathematical point of view [4]. However, it is important to understand the overall dynamical behavior, especially at low frequencies and for more realistic models. We shall focus on that aspect in this paper.

In this regard, recent Raman measurements of long-lived excitations at 118 cm⁻¹ (like pump/probe experiments on Bacteriorhodopsin [5]) appear to be of particular relevance. The fact that there are no local vibrational modes corresponding to such an energy points to the interpretation of the long-lived excitations as collective modes, that is, modes that involve a large number of amino acids. Spectroscopic observations at such energies are thus of particular relevance, since they provide information on large protein domains (i.e.,

on the dynamics of secondary and larger structures) and thus in principle on the functions and the conformational changes of a protein. It was suggested [5] that such modes probably corresponded to large-scale vibrations of the protein—and hence should be damped quickly by the interaction with the surrounding—quite contrary to the observed, particularly long, lifetime of this mode. An alternative interpretation may be based on a comparison with simple models of regular polypeptide systems [6]. It was alternatively suggested [7] that the resonance in Bacteriorhodopsin, which consists chiefly of seven long α helices, was dominated by an excitation along the hydrogen bonded chains along the sides of the α helix parallel to its axis. In this case no major large-scale motion is required, and hence far weaker damping might be expected. The possibilities inherent in pump-probe experiments have opened up a renewed interest in the nonlinear modes in α helices. In an even more recent paper [8], direct evidence is reported for solitonlike behavior (or self-trapping) for poly-gamma-benzyl-L-glutamate. Although related to the NH stretch mode (and thus to Davydov-type solitons), this work demonstrates the renewed interest in the problem of nonlinear excitations in biopolymers, and call for further, more realistic theoretical treatments than the previously studied idealized models.

The phenomena of collective and nonlinear modes in proteins are not fully understood from a theoretical point of view, since both models and numerical simulations face the difficulty of a large number of degrees of freedom with complex interactions. In particular, one may expect nonlinear effects leading to coupled modes and solitonlike modes with long lifetimes.

<footnote>*Email address: dovidio@imedea.uib.es</footnote>

<footnote>†Email address: hbohr@fysik.dtu.dk</footnote>

<footnote>‡Email address: p.a.lindgard@risoe.dk</footnote>
In order to shed some light on this direction, we have, in a previous work [9], modeled an α helix in three dimensions by three coupled Toda [10] lattices. We have shown that such a system has stable soliton solutions that can spontaneously appear if a perturbation parallel to the helix axis is applied. The soliton solution is not destroyed by the 3D coupling. The importance of the Toda theory is that a large number of exact, analytical results have been derived for the Toda potential, which is of a very special exponential form.

It would be highly useful if these results could be used for estimations of the behavior in a realistic situation. In this paper, we show that soliton solutions exist in the more realistic Lennard-Jones (LJ) potential as well, and we demonstrate that the Toda description is accurate in a surprisingly large interval of energies. It has been proven theoretically [11] that the Toda solutions are robust under small perturbations away from exact integrability. However, the predicted lower bounds indicate that it only applies for perturbations of almost infinitesimally small energies. New results in this direction were also obtained very recently by Kosevich et al. [12], who showed that propagating supersonic kink solitons in Fermi-Pasta-Ulam and LJ lattices can be analytically described within the same displacement patterns.

Although a lattice with Lennard-Jones bonds is not directly solvable, in the following we shall provide analytical relations for describing soliton propagation, providing good quantitative agreement with results from direct numerical simulations. We shall actually show from the simulations that in a quite wide range of physical parameters, we can maintain the soliton picture, hence indicating that the soliton phenomenon is inherent to the nature of α-helical dynamics. One might expect that one effect of a coupling between three 1D Toda strands in a 3D structure would be to introduce devastating noise on the single strand—an effect we have not found at low temperatures. However, we have found that the 3D solitons decay in the presence of a static distortion similar to a noise level corresponding to room temperatures. The effect of thermal noise on a 1D Toda lattice was previously systematically studied [13] and it was shown that, in a certain temperature interval, solitons do exist, but their interaction gives rise to delocalization effects and a decay to broadband 1/f noise at low frequencies (corresponding to the transformation into low-frequency phonons).

More realistic 3D models of an α helix have been investigated [14–16] with respect to the Davydov soliton mechanism. We shall not focus on that aspect here. The models differ in detail from ours. Hennig [16] used a confining cage and initiated soliton solutions by nonlinear map methods. Christensen et al. [15] studied an unsuspended 3D helix with Morse interactions. They found both the phonon and the soliton solutions numerically, but did not attempt to compare that with Toda solutions, nor did they attempt a comparison with experiments. Previously, Tuszynski et al. discussed in the context of the Davydov model the emergence of Jacobi elliptic waves, including the cnoidal waves, and calculated numerically the corresponding free energy [17].

## II. COLLECTIVE MODES OF α HELICES AS SOLITONS ON THE HYDROGEN BONDS

The involvement of the hydrogen bonds for the observed excitation at 118 cm⁻¹ can be supported also by direct cal-

culations [7] as follows. Consider a linear chain with a realistic potential. The linear modes correspond to phonons, and the frequency ν of a phonon (at maximum density of states) is given by

$$2\pi\nu = 2\sqrt{k/m}, \qquad (1)$$

where ν is the frequency, k the force constant, and m the average mass of the amino acids. In a chain of amino acids connected by hydrogen bonds H—O of length r₀ = 1.9 Å, k ≈ 1.41 × 10⁴ dyn/cm and m = 1.7 × 10⁻²² g is the average mass of the residues. An estimation can be obtained using 1/λ = ν/c with λ the wavelength and c the speed of light and gives 97 cm⁻¹. If an infinite poly(L-alanine) α helix is considered, a complete normal mode calculation [19] gives a peak at 118 cm⁻¹. Thus an excitation on the H bonds is certainly to be expected near the observed frequency. Both phonons and solitons could exist on the chain in a nonharmonic potential. However, solitons need nonlinearities in the potential in order to exist. Such nonlinearities become relevant at high energies, which can be quantified using a Lennard-Jones potential for the H—O bond. The nonlinearities are hence expected to play a role at the high level of pumping in the pump-probe experiments [5].

## III. REALISTIC POTENTIALS AND ANALYTICAL TOOLS

For a molecular solid, the interaction between neighboring atoms with equilibrium distance r₀ is to a good approximation given by the well-known standard 6-12 Lennard-Jones potential,

$$V_{\text{LJ}}(r) = A \left[ \left(\frac{r_0}{r}\right)^{12} - 2\left(\frac{r_0}{r}\right)^6 \right] \qquad (2)$$

with the potential strength A, which is also the depth of the energy minimum at r=r₀. For excitation energies higher than A, the bond is broken. It is nonlinear, and can be expanded around r₀ as

$$V_{\text{LJ}}(r-r_0) = A \left[ -1 + 36 \left(\frac{x}{r_0}\right)^2 - 252 \left(\frac{x}{r_0}\right)^3 + o[x^4] \right], \quad (3)$$

where x=(r-r₀). In one dimension and around a minimum of the same depth at r₀, the potential introduced by Toda [10] has the following form:

$$V_{\text{Toda}}(r-r_0) = -A - \frac{a}{b} + \frac{a}{b}e^{-b(r-r_0)} + a(r-r_0), \quad (4)$$

which likewise can be expanded as

$$V_{\text{Toda}}(r-r_0) = -A + \frac{1}{2}abx^2 - \frac{1}{6}ab^2x^3 + o[x^4]. \qquad (5)$$

Hence the potentials are identical in form up to o[x⁴] and both include a third-order nonlinear term. The harmonic term giving the result Eq. (1) is usually written V(r-r₀) = 1/2kx². By equating the coefficients in the expansion, we immediately get the relation between the parameters,
<figure><img src="image_2.png" /><figcaption>FIG. 1. Comparison of Lennard-Jones and Toda potential.</figcaption></figure>

$$k = ab = 72A/r_0^2, \quad a = kr_0/21, \quad b = 21/r_0, \quad (6)$$

where $r_0$ is the equilibrium position and a, b, and A are parameters (in the following, we often set $r_0$ equal to 1 for simplicity).

As we can see from Fig. 1, the Toda potential is similar to a Lennard-Jones potential in a substantial region around the equilibrium. This is very important since there exist a large number of exact results for the Toda model in one dimension. However, an important difference between the Lennard-Jones and the Toda potentials is that the latter is attractive for all interatomic distances beyond $r_0$, hence no bond breaking is possible. The subject of this paper is to explore how far the Toda analytic, exact results can be used for describing the α helices having realistic (Lennard-Jones-like) potentials and three-dimensional architecture.

The measurement units have been chosen normalizing the length and the mass to the equilibrium distance of the hydrogen bond ($r_0$) and to the average mass (m) of an amino acid, respectively. Time is measured in picoseconds. Calling $x_j$ and $p_j$ the space coordinates and momenta of the jth amino acid, and numbering the amino acids as they appear along the helix, the Hamiltonian is given by

$$H = E_{\text{kin}} + V_{\text{peptide}} + V_H = \frac{1}{2m} \sum_{j=1}^{N} p_j^2 + \sum_{j=1}^{N} V_{a,b}(x_j, x_{j+1}, x_{j-1}) \\ + \sum_{j=1}^{N} V_{c,d}(x_j, x_{j+3}, x_{j-3}), \quad (7)$$

where V are Toda or Lennard-Jones potentials of parameters a, b and c, d (or alternatively A and $A_p$), and the sums have obviously to include only the amino acids with $j-1 > 0, j-3 > 0, j+1 < N$, and $j+3 < N$. The equations of motion are obtained straightforwardly and can be solved numerically with a standard Runge-Kutta integrator. In order to follow the propagation of a wave, it is useful to introduce the energy density per site,

$$E_j = E_{\text{kin},j} + \frac{1}{2}V_{\text{cov},j} + \frac{1}{2}V_{\text{hydr},j}. \quad (8)$$

<figure><img src="image_3.png" /><figcaption>FIG. 2. The cnoidal wave form for different wavelengths and amplitudes. (a,b,c) for $\lambda=(2,4,10)$. The curves on which the atom positions are indicated by a (small,medium,large) size dot correspond to the modulus $k^2=(0.3,0.9,0.99)$. Obviously even for $k^2=0.3$ the waveform is close to a harmonic wave form for all wavelengths, whereas for larger k there are significant differences.</figcaption></figure>

We check the existence of solitary waves by performing a simulation of the system at rest (initially) after a perturbation at time $t=0$ on one of the amino acids at the terminal points of the helix. This is done by setting all $x_j$ variables at the equilibrium values and by choosing $p_j=0$ for all the amino acids except one, which has a momentum corresponding to the energy of the perturbation.

We briefly remind the reader of the exact results for the 1D Toda lattice [10]. The dynamics is described by solitons with the energy

$$\mathcal{E} = \frac{2a}{b}(\sinh \kappa \cosh \kappa - \kappa). \quad (9)$$

The profile of the soliton, in terms of the displacements $x_n$ from the equilibrium distance, is given by
<figure><img src="image_4.png" /><figcaption>FIG. 3. The reduced dispersion curve $\omega=2\pi\nu/\sqrt{ab/m}$ as a function of $q=2/\lambda$ for the Toda lattice for the modulus $k^2=(0.3,0.9,0.99)$. For $k^2=0.3$, the dispersion curve $\omega(q)$ is very close to the harmonic one, which peaks at $\omega(1)=2$.</figcaption></figure>

$$x_n = -\frac{1}{b}\log\left[1 + \frac{m}{ab}\beta^2\mathrm{sech}^2(\kappa n \pm \beta t)\right]. \quad (10)$$

The speed $v$ is

$$v = \frac{\beta}{\kappa} = \sqrt{\frac{ab}{m} \frac{\sinh \kappa}{\kappa}}. \quad (11)$$

The parameter $\kappa$ completely determines the energy $\mathcal{E}$ and the dynamics of the solitons ($1/\kappa$ is proportional to the width of the soliton).

Another important family of solutions to the Toda lattice are the cnoidal waves, which are periodic and equivalent to phonons. An interesting properties of cnoidal waves is that they are made of equally spaced solitons, with the distance between them being the wavelength, $\lambda$. Thus a train of solitons can be connected to phonons, and much can be learned about the cnoidal waves by studying a single soliton.

In general, the periodic solution to the Toda lattice is written in terms of the Jacobi elliptic functions as

$$x_n = -\frac{1}{b}\log\left[1 + \frac{dn^2(2K[n/\lambda - vt] - E/K)}{sn^{-2}(2K/\lambda) - 1 - E/K}\right], \quad (12)$$

where $K=K(k)=\int_0^1 [(1-k^2x^2)(1-x^2)]^{-1/2}dx$ and $E=E(k)=\int_0^1 [(1-k^2x^2)/(1-x^2)]^{1/2}$ are elliptic integrals of the first and second kind.

The dispersion relation between the wavelength $\lambda$ and the frequency $v$ of the periodic wave is

$$2\pi\nu = \sqrt{\frac{ab}{m}\frac{\pi}{K}} \left(\frac{\mathrm{sn}^{-2}(2K/\lambda) - 1}{K}\right)^{-1/2} = \left|_{k \to 0} \sqrt{\frac{ab}{m}} 2|\sin(\pi/\lambda)|\right|. \quad (13)$$

Hence it is tending to the harmonic one when the amplitude of the wave described by Eq. (12) is small, corresponding to the limit where the modulus $k$ of the elliptic integrals is small and $K \sim E \sim \pi/2$ [notice the dispersion curve Eq. (13) is periodic with a period $2K/r_0$ which is different from $\pi/r_0$ for finite $k$ values]. It is interesting that the functional form of the wave profile and the dispersion only depends on $k$ and the wavelength $\lambda$, whereas it does not depend on the

<figure><img src="image_5.png" /><figcaption>FIG. 4. Space-time plot of a one-dimensional lattice of elements connected by the Lennard-Jones potential. The system is initialized with a wave corresponding to a perfect Toda soliton [Eq. (10)] of low energy. The parameters of the soliton are $\kappa=0.5$, $E=9.39 \text{ cm}^{-1}$. In this energy range, the Toda and Lennard-Jones potentials are quite similar and thus the Toda soliton survives for a long time even if the potential is not Toda.</figcaption></figure>

constants $a, b$ describing the potential. Hence it is the same for the set Eq. (6) corresponding to the 6-12 LJ potential as for the more anharmonic one corresponding to the 10-12 LJ potential [18].

It is instructive to discuss these results in physical terms. In Fig. 2, we show the cnoidal wave form for different wavelengths and amplitudes. For small wavelengths, the largest excursions are towards the weak side of the potential, hence the energy is lower than for the harmonic potential, i.e., a renormalized dispersion is expected, whereas, surprisingly, for longer wavelength the largest excursions are in the direction of the steep side of the potential, hence a higher energy. These are compression modes corresponding to a train of well-separated pure solitons. In Fig. 3, we plot the corresponding dispersion curves versus the wave vector $q=2/\lambda$. The sequence is the same as in Fig. 2. We notice the dispersion is reduced for short-wavelength, large-amplitude modes. However notice, for the Toda potential, that the solutions are exact and hence there is no damping of the modes.

## IV. RESULTS FOR 1D LATTICES

In the previous sections, we have seen that a Toda potential fits a Lennard-Jones potential in a neighborhood of a minimum, up to the third order. By continuity, this means that solutions valid for a Toda lattice can be used for representing solutions for Lennard-Jones lattices, for a time length that grows when $E \to 0$. In order to test this idea, we study the behavior of a 1D Lennard-Jones lattice initialized with a condition that gives rise to an integrable soliton solution in a Toda lattice. Since we are especially interested in long-lived solutions, we use a one-dimensional lattice with periodic boundary conditions. This geometry ensures that a soliton can travel over the lattice indefinitely, even if the number of sites is finite, with the wave going out from one end and
<figure><img src="image_6.png" /><figcaption>FIG. 5. Same as Fig. 4 but with a soliton of higher energy ($\kappa = 2, E = 1031$). In this region, the Toda and Lennard-Jones potentials are different, but nevertheless the Toda soliton persists for a long time on the Lennard-Jones lattice.</figcaption></figure>

coming back at the opposite side. As a first simulation, we initialize the lattice with $N=40$ points connected by Lennard-Jones potentials with a soliton of low energy, taken from Eq. (10) and setting $\kappa=0.5$.

This choice corresponds to a total energy for the lattice of 9.39 cm⁻¹, and to a maximum energy per site of 3.14 cm⁻¹. The soliton can be qualitatively seen by making a space-time plot of the lattice. In Fig. 4, the first 20 ps are plotted and show that the soliton persists without any apparent deformation. A similar behavior is observed for all the time of the simulation (8000 ps). The persistence of the soliton at such low energy is indeed to be expected, due to the similarity of the Toda and Lennard-Jones potentials. However, the same behavior appears to exist for much higher energies: Fig. 5 reports a soliton for the case of $\kappa=2$, corresponding to a total energy of 1294 cm⁻¹ and to a maximum energy per site of 1031 cm⁻¹. This result is surprising since the Toda and Lennard-Jones potential greatly differ at this energy (see Fig. 1).

In order to study the persistence of the soliton solutions more quantitatively, we plot a projection of the trajectories

<figure><img src="image_7.png" /><figcaption>FIG. 6. Projection of the trajectories of Toda-like solitons of different energy (30–1300 cm⁻¹, increasing patterns) evolving on a 1D Lennard-Jones lattice. The trajectories stay near to close orbits for 8000 ps.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Projection of the trajectories of Toda-like cnoidal waves of different energy (30–1300 cm⁻¹, increasing patterns) evolving on a 1D Lennard-Jones lattice ($k^2=0.1, 0.3, 0.6, 0.98, \lambda=4r_0$). The trajectories stay even nearer to the close orbits for the 8000 ps than the solitonlike solutions.</figcaption></figure>

on the two-dimensional space given by position (bond length $\Delta_1=x_1-x_0$) and speed of the first lattice element. This method is a good way for studying the persistence of a soliton on a periodic lattice. In fact, a soliton that travels without changing its shape comes back to the same position after a loop over the lattice. The overall dynamics is thus a cycle.

Repeating the simulation for $\kappa=0.875, 1.25, 1.625$, and $2$ (corresponding to the energy range 30–1300 cm⁻¹) and integrating the system for 400 ps, we plot all the trajectories in Fig. 6. The picture shows that the energies of the first two cases evolve on closed orbits up to deviations of about 1%. That confirms the observation obtained from the space-time plots that Toda soliton solutions persist on a Lennard-Jones lattice for a long time (8000 ps, in the simulation), up to energies of more than 1000 cm⁻¹. Figures 7 and 8 show that a similar agreement is also obtained with cnoidal waves ($k^2=0.1, 0.3, 0.6, 0.98$) for both $\lambda=4r_0$ and $\lambda=2r_0$. Hence, our "phonons," the cnoidal waves, on a 1D lattice, have lifetimes well above 8000 ps—even for large-amplitude modes. We remark that long-lived protein modes near 100 cm⁻¹ oscillate of the order of 1500 cycles, which corresponds to having a lifetime of the order of 500 ps.

For high energies, one expects that the differences between the two potentials become relevant and that a Toda solution can represent a Lennard-Jones solution for a very short time only. The fact that on the contrary the two solu-

<figure><img src="image_9.png" /><figcaption>FIG. 8. Same as Fig. 7 for $\lambda=2r_0$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 9. Difference between Lennard-Jones and Toda potentials, as a function of the (Lennard-Jones) potential depth. The upper branch corresponds to $x<r_0$, while the lower branch to $x>r_0$. Notice that on the lower branch, the two potentials differ only a few $cm^{-1}$ over an energy range of several hundreds of $cm^{-1}$.</figcaption></figure>

tions are close to each other for higher energies as well can
possibly be explained by noticing that the agreement be-
tween a Toda and Lennard-Jones potential is not symmetric
around the minimum. Small differences of energies corre-
sponding to distances larger than the equilibrium correspond
to large deviation between the two potentials, since the
Lennard-Jones potential has an asymptotically flat behavior
and the Toda potential has an asymptotically linear growth.
On the contrary, for distances smaller than the equilibrium
position, the two potentials keep a fairly good agreement
even for comparatively larger energy differences (see Fig. 9).
The persistence of Toda solitons of high energy on a
Lennard-Jones lattice can be explained from the fact that a
soliton is a compression wave and thus the potentials are
explored for distances below the equilibrium value, in the
region where the two potentials best match. In the last sec-
tion, we will come back to this issue and argue that the
persistence of the soliton may actually have a deeper and
general meaning, being related to the persistence of inte-
grable solutions in the Hamiltonian system.

On this background we can also discuss the cnoidal
waves. For long wavelengths, the above picture holds and
the excitations for the LJ potential are well described by the
Toda solutions. For small wavelength, the difference in the
weak side of the potentials will play a role for large ampli-
tudes, and tends to give larger renormalizations of the dis-
persion curves (actually depending on *a*, *b*). To estimate the
importance of this, let us calculate the amplitude *X* of the
zone boundary mode, which is at the highest density of
states—and hence corresponds to the largest optical re-
sponse. It has an energy of *hv*<sub>max</sub>=165 K *k*<sub>B</sub>. For a harmonic
oscillator, the energy depends on the amplitude as *E*
= 1/2*m*(2*π*v*)<sup>2</sup>*X*<sup>2</sup>. We have to use quantum mechanics to
accurately calculate the population of the mode: *n*
= 1/[exp(*hv*/*k*<sub>B</sub>*T*) − 1] ~ 1.36 at room temperature, *T*
= 300 K. This, on the other hand, yields the energy *E*(*T*)
=h*v*(n + 1/2). Hence *X*=(1/2*π)√(2*n*+1)*h*/*m*v*=0.05*r*<sub>0</sub> at room

<figure><img src="image_11.png" /><figcaption>FIG. 10. Three-dimensional model of an α helix. Amino acids are mass points (dots) connected by Lennard-Jones potentials giving hydrogen bonds (continuous line) and peptide bonds (dashed line). The helix here has open boundary conditions.</figcaption></figure>

temperature. Comparing with Fig. 1, we notice the Toda and
LJ potentials are almost equal for such amplitudes, and ac-
cordingly the Toda results should be valid to a good approxi-
mation for the 115 cm⁻¹ mode. By comparing with Fig. 2(a),
we find that this amplitude corresponds approximately to that
for the k²=0.9 curve. For this, according to Fig. 3 the renor-
malization is approximately 10%.

V. 3D LATTICES WITH α-HELIX GEOMETRY

Let us now move to a lattice geometry closer to an actual
α helix. The system we are considering is depicted in Fig.
10. Each amino acid is represented by a point, and Lennard-
Jones potentials of the form of Eq. (2) are used to model
hydrogen bonds as well as peptide bonds. The parameters
used are A=355.93 cm⁻¹ for the hydrogen bonds and A
=1281.5 cm⁻¹ for the peptide bond. In a previous work [9],
we considered a similar geometry of points connected by
Toda potentials. We showed that such a system admits a
solitary wave propagating along the whole helix when an
instantaneous perturbation parallel to the hydrogen bonds is
given to the first amino acid of the chain. Due to the analogy
between Lennard-Jones and Toda potentials, as discussed in
the previous section, it is no surprise that the same phenom-
enon arises in the case of Lennard-Jones bonds. However,
the following two scenarios are both possible: either a single
soliton is propagating along the peptide bonds (the back
<figure><img src="image_12.png" /><figcaption>FIG. 11. Energy vs speed for different strength of the peptide bond, for solitons on a 3D $\alpha$ helix with Lennard-Jones potential. The curve does not change.</figcaption></figure>

bone) or three entrained solitons are propagating along the three lattices of hydrogen bonds (as we argued for the case of a helix made of Toda potentials).

The two situations point to very different roles of the bonds. In the first case (one soliton over the peptide bond chain), the hydrogen bonds confine the soliton while the peptide bonds provide a (helical) lattice where the soliton travels. In the second case (three solitons over the hydrogen bond chains), the weaker hydrogen bonds form three parallel, one-dimensional lattices that support solitons and the stiffer peptide bonds act as a strong coupling. Such a coupling has different effects on the soliton dynamics: after the initial perturbation (that is applied on one chain only), the peptide bonds mirror and entrain identical waves on the other two lattices, and after a short transient, the effect of the coupling becomes negligible, since the solitons are already entrained and propagate on identical lattices (in a more realistic simulation, the peptide bond would correct small propagation differences arising from lattice and steric dishomogeneities). This possible scenario has important consequences from a modeling viewpoint. In fact, after the formation of the triple-soliton solution, one can expect that each entrained solitary wave can be described to a good approximation as a *free* soliton on a one-dimensional lattice completely characterized by the hydrogen bond. This allows us to use the result of the previous section and model quantitatively Lennard-Jones solitons on a 3D helix with 1D Toda solitons.

However, we have still to show that such a scenario corresponds to three solitons on each H-bond chain and not one soliton on the peptide bond chain. In order to show that the solitons arise on the H bonds and that the peptide bonds are negligible, after the introduction of the triple solitons, we present numerical simulations where the behavior of the Lennard-Jones solitons is studied at different strength of the peptide and hydrogen bonds. We shall see that the soliton dynamics is greatly affected by changes in the hydrogen bond, being only slightly influenced by changes in the peptide bond.

<figure><img src="image_13.png" /><figcaption>FIG. 12. As in Fig. 11, but changing the strength of the hydrogen bond. The curve changes.</figcaption></figure>

We perform simulations over a three-dimensional $\alpha$-helix modeled as described in Fig. 10. We induce solitons of different energies using perturbation of variable energy (in the range 200–2000 cm$^{-1}$, inducing solitons between 100 cm$^{-1}$ and 500 cm$^{-1}$). The soliton speed is computed by finding the slope of the line that connects the soliton tips in the space-time plot. The soliton energy is computed by summing the energy of the elements around the soliton tips in a window of five elements and averaging in time. In this way, we construct the curve that gives energy versus speed. We construct curves first for a range of peptide strengths (in a range of 10% around the parameters given in the previous section), then repeating the simulations but varying the hydrogen bonds instead.

The outcome of the simulations (Figs. 11 and 12) is quantitatively consistent with our hypothesis: the characteristics of the solitons of the 3D LJ helix are strongly determined by the hydrogen bonds only and the energy-speed curve is in quantitative agreement with the analytical prediction of a soliton traveling over a one-dimensional Toda lattice (Fig. 13) on the H bonds.

## VI. DISCUSSION

We have studied a string of masses interacting with a Lennard-Jones potential and shown that the excitation spectrum up to surprisingly high energies can be accounted for by the theory for the Toda lattice (with parameters fitted to the Lennard-Jones potential). An explanation of the high-energy range lies in the fact that the solitons are compression waves, thus the region of the potential that is involved is the one where Toda and Lennard-Jones potentials are closer. Likewise we have found that the family of periodic solutions, the cnoidal waves, are very stable—if not more—also in the Lennard-Jones potential systems. This is because these are the exact solutions of the Toda system with periodic boundary conditions, whereas a single soliton is only approximately so—unless the tails from all periodic solitons are properly included.
<figure><img src="image_14.png" /><figcaption>FIG. 13. Energy vs speed normalized by the sound speed of the hydrogen bond for solitons on a 3D $\alpha$ helix with Lennard-Jones potential. The continuous line is the prediction using Toda theory on the one-dimensional lattice.</figcaption></figure>

Next, we have studied an unsuspended $\alpha$ helix in three dimensions with different strength of LJ interactions along the helix backbone and along the three hydrogen bonds along the sides of the helix. It is again found that soliton solutions exist in a surprisingly large energy interval. A perturbation inflicted on one of the chains is rapidly transmitted to the two neighboring chains and a triple soliton is propagated lossless along the helix. No overall bending of the helix is observed due to the soliton motion and the three solitons stay phase-locked close to each other. This is contrary to the finding in a somewhat similar study of a 3D helix in a cage [16] in which a bending is found together with an asymmetric propagation along the tree chains. In that model, the radius of the helix decreased in the soliton region, whereas in our case we find an expansion, which is also found by Christensen et al. [15]. The reason for the difference might be either the differences in the models or perhaps rather the quite different methods of solution. In this paper, we have used a minimalistic model and used the most straightforward theoretical and numerical methods. Our soliton corresponds to a compression wave, which physically would be expected to be accompanied by an expanding radius, as we have found. If excited by a pull instead of by a push along the helix axis direction, a compression soliton is also found (due

to the recoil of the returning mass point). We have demonstrated the accuracy with which the speed of the simulated LJ soliton can be described analytically by the corresponding exact Toda soliton.

Summarizing, we have demonstrated that in a realistic 3D $\alpha$ helix described by LJ interactions, a soliton(like) solution is highly stable and can transport large energy quanta without loss and without making overall movements of the helix, even when not confined in a cage. The solitons run on hydrogen bonds whose natural linear excitations are phonons with energies in the region of $100 \text{ cm}^{-1}$ and can be described analytically using Toda theory.

An interesting development of this work is to understand the mechanism by which the soliton decays in the lattice phonons that are then detected in spectroscopic experiments. In fact, in the case of a periodic Toda potential, the soliton persists forever and no transfer of energy to phonons is possible. Since the lattice is fully integrable, even perturbations of the soliton lead to periodic or quasiperiodic solutions. Such perturbed solutions can be expected as tori nested around the periodic soliton solution. Due to integrability, such tori do not interact with the soliton, but in the case of the 3D $\alpha$-helix lattice, we can expect the boundary conditions of the helix as well as the deviation of the Lennard-Jones potential from the Toda shape to break integrability and heavily affect the phase portrait. The soliton can be reflected at the helix boundaries and still correspond to a periodic solution. But in the 3D helix, according to KAM theory, tori that are resonant with the soliton solution are broken into nonintegrable motion, allowing energy flow between phonons and solitons. However, as we have shown in Sec. IV, the Lennard-Jones potential explored during the soliton motion is very similar to the Toda potential, and thus the energy transfer will take place at a slow rate. Such slow decay of the solitons into the phonon bath will thus provide an energy source for phonons of that energy. This might be an explanation for the unexpectedly long lifetime of excitations in the $100 \text{ cm}^{-1}$ region observed by pump-probe free-electron laser experiments [5].

## ACKNOWLEDGMENT

The authors are thankful to Emilio Hernández-García for useful discussion.

[1] C. Branden and J. Tooze, *Introduction to Protein Structure* (Garland, New York, 1991).

[2] H. Wang and G. Oster, Nature (London) **396**, 279 (1998).

[3] A. S. Davydov, J. Theor. Biol. **38**, 559 (1973).

[4] A. C. Scott, Phys. Rep. **217**, 1 (1992); A. V. Zolotaryuk, K. H. Spatschek, and A. V. Savin, Phys. Rev. B **54**, 266 (1996); K. Kundu, Phys. Rev. E **61**, 5839 (2000).

[5] A. Xie, A. F. G. van der Meer, and R. H. Austin, Phys. Rev. Lett. **88**, 018102 (2002).

[6] S. E. M. Colaianni and O. F. Nielsen, J. Mol. Struct. **347**, 267 (1995).

[7] O. F. Nielsen, P.-A. Lindgård, and H. G. Bohr (unpublished).

[8] J. Edler, R. Pfister, V. Pouthier, C. Falvo, and P. Hamm, Phys. Rev. Lett. **93**, 106405 (2004).

[9] F. d'Ovidio, H. G. Bohr, and P.-A. Lindgård, J. Phys.: Condens. Matter **15**, s1699 (2003).

[10] M. Toda, *Nonlinear Waves and Solitons* (Kluwer Academic Publications. Dordrecht, 1983).

[11] E. Fermi, J. Pasta, S. Ulam, and M. Tsingou, in *The Many-Body Problems*, edited by D. C. Mattis (World Scientific, Sin-
gapore, 1993) (reprinted).

[12] Yu. A. Kosevich, R. Khomeriki, and S. Ruffo, Europhys. Lett. 66, 21 (2004).

[13] M. Jensen and W. Eberling, Physica D 141, 117 (2000).

[14] O. H. Olsen, M. R. Samuelsen, S. B. Petersen, and L. Nørskov, Phys. Rev. A 38, 5856 (1988); A. La Magna, R. Pucci, G. Piccitto, and F. Siringo, Phys. Rev. B 52, 15273 (1995); A. V. Zolotaryuk, P. L. Christiansen, and A. V. Savin, Phys. Rev. E 54, 3881 (1996).

[15] P. L. Christiansen, A. V. Zolotaryuk, and A. V. Savin, Phys. Rev. E 56, 877 (1997).

[16] D. Hennig, Phys. Rev. B 65, 174302 (2002).

[17] J. A. Tuszynski, M. L. A. Nip, P. L. Christiansen, M. Rose, and O. Bang, Phys. Scr. 51, 423 (1995).

[18] One may for the hydrogen bonds use a 10-12 Lennard-Jones potential [20] and include an angular dependence. Such a potential has a weaker attraction and more anharmonicity. For this potential, which can be expanded as

$$V_{LJ}^{10-12} \sim -A \left[ -\frac{1}{5} + \frac{12}{10} \left(\frac{x}{r_0}\right)^2 - 100 \left(\frac{x}{r_0}\right)^3 + o[x^4] \right],$$

the parameter relation to the Toda parameters corresponding to Eq. (6) is

$$k = ab = 24A/r_0^2, \quad a = kr_0/25, \quad b = 25/r_0.$$

We shall not investigate this model here. We expect the qualitative results would be the same, although the energy range of validity of the Toda description will be restricted.

[19] S. H. Lee and S. Krimm, Biopolymers 46, 283 (1998).

[20] A. Irbäck, J. Phys.: Condens. Matter 15, s1797 (2003).