We here introduced the form of bio-energy in living system and eluciduted again the new theory of bio-energy transport along protein molecules in living systems based on the changes of structure and conformation of molecules arising from the energy, which is released by hydrolysis of adenosine triphosphate (ATP). In this theory, the Davydov’s Hamiltonian and wave function of the systems are simultaneously improved and extended. A new interaction have been added into the original Hamiltonian. The original wave function of the excitation state of single particles have been replaced by a new wave function of two-quanta quasicoherent state. In such a case, bio-energy is carried and transported by the new soliton along protein molecular chains. The soliton is formed through self- trapping of two excitons interacting amino acid residues. The exciton is generated by vibrations of amide-I (C=O stretching) arising from the energy of hydrolysis of ATP. The properties of the soliton are extensively studied by analytical method and its lifetime for a wide ranges of parameter values relevant to protein molecules is calculated using the nonlinear quantum perturbation theory. The lifetime of the new soliton at the biological temperature 300 K is enough large and belongs to the order of 10- 10 second orτ/τ0≥700. The different properties of the new soliton are further studied. The results show that the new soliton in the new model is a better carrier of bio-energy transport and it can play an important role in biological processes. This model is a candidate of the bio-energy transport mechanism in protein molecules.

Keywords: Form; Living system; Bio-energy; protein; Biological energy; Soliton; ATP hydrolysis; Amide; Exciton; Life time; Amino acid; Quasi-coherent state.

As it is known, Kal’kar first proposed the idea of aerobic phosphorylation, which is carried out by the phosphorylation coupled to the respiration. Belitser studied in detail the stoichometric ratios between the conjugated bound phosphate and the absorption of oxygen and gave further the ratio of the number of ionorganic phosphate molecules to the number of oxygen atoms absorbed during the respiration, which is not less than two. He thought also that the transfer of electrons from the substrate to the oxygen is a possible source of energy for the formation of two or more ATP molecules per atom of absorbed oxygen. Therefore Belitser and Kal’kar’s research results are foundations establishing modern theory of oxidative phosphorylation of ATP molecules in the cell [1-3].

In such a case we must know clearly the mechanism and properties of the oxidation process, which involves the transfer of hydrogen atoms from the oxidized molecule to another molecule, in while there are always protons present in water and in the aqueous medium of the cell, thus we may only consider the transfer of electrons in this process. The necessary number of protons to form hydrogen atoms is taken from the aqueous medium. The oxidation reaction is usually preceded inside the cell under the action of special enzymes, in which two electrons are transferred from the food substance to some kind of initial acceptor, another enzymes transfer them further along the electron transfer chain to the second acceptor etc. Thus a water molecule is formed in which each oxygen atom requires two electrons and two protons.

The main initial acceptors of electrons in cells are the oxidized forms NAD+ and NADP+ of NAD (nicotine amide adenine dinucleotide or pyridine nucleotide with two phosphate groups) molecules and NADP(nicotine amide adenine nucleotide phosphate or pyridine nucleotide with three phosphate groups) as well as FAD (flavin adenine dinucleotide or flavoquinone) and FMN (flavin mononucleotide).The above oxidized forms of these molecules serve for primary acceptors of electrons and hydrogen atoms through attaching two hydrogen atoms [3], which is expressed by $$NAD{P}^{+}+2H^{+}+2e^{-}\to NADP\cdot H+H^{+}$$ Where NADP molecule becomes the reduced molecule NADPH. The NAD molecule has also the same active center as the NADP molecule; it can be converted to the reduced molecule NAD. H under combining with two atoms of hydrogen according to the reaction [3]: $$NA{D}^{+}+2H^{+}+2e^{-}\to NAD{P}^{+}+H^{+}$$ The NAD⁺ and NADP⁺ are the enzymes, which can perform the reaction of dehydrogenation on compounds containing the group of atoms through removing two hydrogen atoms.

In the presence of enzymes, such as pyridine-dependent hydrogenases, and with the participation NAD⁺ and NADP⁺ molecules two hydrogen atoms, including two protons and two electrons, are removed from this group of atoms. One proton and two electrons combine with the NAD⁺ and NADP⁺ molecule converting them to the reduced forms NADP . H or NAD . H and the second proton is released. This mechanism can be also used to oxidize lactic acid (lactate) with the formation of pyruvic acid (pyruvate) and NAD H, in which the reduced molecules NADP . H and NAD . H serve as electron donors (reducing agents) in other reactions. They are involved in a large number of biosynthetic processes, such as in the synthesis of fatty acids and cholesterol.

Therefore, the molecule NAD • H can serve as an electron donor in the process of oxidative phosphorylation, then the phosphorylation reaction is of [3] $$H^{+}+NAD\cdot H+3H_{3}P{O}_4+3ADP+1/2O_2\to NA{D}^{+}+4H_{2}O+3ATP$$ Where ADP is called the adenosine diphosphate. The abbreviated form of this reaction can be written as $$ADP+P_i\to ATP+H_{2}0$$ Thus three ATP molecules are formed in the reaction, in which the synthesis of ATP molecule are carried out through the transfer of two electrons from the NAD • H molecule along the electron transport chain to the oxygen molecule in the mitochondria. In this way the energy of each electron is reduced by 1.14 eV. The reaction is called the phosphorylation of ADP molecules.

However, an ATP molecule can reacts with water in an aqueous medium, which results in the energy release of about 0.43eV under normal physiological conditions by virtue of some special enzymes. The reaction can be represented by $$AT{P}^{4-}+H_{2}O\to AD{P}^{3-}+HP{O}_4^{2-}+H^{+}+0.43eV$$ Its abbreviated form is of $$ATP+H_{2}O\to ADP+P_i$$ In this process ATP molecules are transformed as ADP molecules and the bio-energy of about 0.43eV is also released. Then it is referred to as de-phosphorylation reaction of ATP molecules.

We know from the above representations that an increase in free energy $\Delta \text{G}$ in reaction and its decrease in reaction depend on their temperatures, concentrations of the ions Mg²⁺ and Ca²⁺ and on the pH value of the medium. Under the standard conditions $\Delta G_0$ =0.32 eV ( ~7.3 kcal/mole). If the appropriate corrections are made taking into consideration the physiological pH values and the concentration of Mg²⁺ and Ca^2* inside the cell as well as the normal values for the concentrations of ATP and ADP molecules and inorganic phosphate in the cytoplasm we can obtain a value of ~0.54eV ( ~12.5 kcal/mole) for the free energy in the hydrolysis of ATP molecule. Hence the free energy for the hydrolysis of ATP molecules is not constant. But it is impossibly the same at different sites of the same cell if these sites have different concentrations of ATP, ADP, Pi, Mg²⁺, Ca²⁺.

On the other hand, cells contain a number of phosphorylated compounds the hydrolysis of which in the cytoplasm is associated with the release of free energy. Then the values for the standard free energy of hydrolysis for some of these compounds are also different.

The enzymes carrying out the above synthesis of ATP molecules from ADP molecules and inorganic phosphate in the coupling membranes of mitochondria are the same as in the cytoplasmic membranes of bacteria, which are mainly composed of F and F, which are joined to each other by the small proteins F5 and F6 . These proteins form the F - F complex or the enzyme ATPase, in which F is composed of five protein subunits and has the shape of a sphere with a diameter of about 9nm which projects above the surface of the membrane in the form of a protuberance. In the coupling membrane of mitochondria and the cytoplasmic membrane of bacteria the complex F - F is positioned so that the enzyme F is on the inside of the membrane [1-3].

The enzyme F can extend from one side of the membrane to the other and has a channel which lets protons through. When two protons pass through the complex F - F in the coupling mitochondrial membrane one ATP molecule is synthesized inside the matrix from an ADP molecule and inorganic phosphate. This reaction is reversible. Under certain condition the enzyme transports protons from the matrix to the outside using the energy of dissociation of ATP molecules, which may be observed in a solution containing isolated molecules of enzyme F and ATP. The largest two proteins in F, which is composed of five protein molecules, take part in the synthesis and dissociation of ATP molecules, the other three are apparently inhibitors controlling these reactions.

After removing enzyme F molecules from mitochondria the remaining F enzymes increase greatly the permeability of protons in the coupling membranes, which confirms that the enzyme F has really a channel for the passage of protons which is constructed by the enzyme F. However, the complete mechanism for the synthesis of ATP molecules by the enzyme ATP-ase is still not clearly known up to now.

As it is known, many biological processes, such as muscle contraction, DNA reduplication, neuroelectric pulse transfer on the neurolemma and work of calcium pump and sodium pump, and so on, are associated with bioenergy transport through protein molecules, where the energy is released by the hydrolysis of adenosine triphosphate (ATP) in the living systems. Thus there here are always biological processes of energy transport from production place to absorption place in the living systems. In general, the bioenergy transport is carried out by virtue of protein molecules. Therefore, the study of the bioenergy transport along protein molecules is a very interesting subject in biology and has important significance in life science. However, understanding the mechanism of bioenergy transport in biomacromolecular systems has been a long-standing problem that remains of great interest today. As an alternative to electronic mechanisms [1], one can assume that the energy is stored as vibrational energy in the C=0 stretching mode (amide-I) of a protein molecular chain of polypeptide. Following Davydov’s idea [2], ones take into account the coupling between the amide-I vibrarional quantum (exciton ) and the acoustic phonon (molecular displacements) in the amino acid residues; Through the coupling, nonlinear interaction appears in the motion of the vibrartional quanta, which could lead to a self-trapped state of the vibrational quantum. The latter plus the deformational amino acid lattice together can travel over macroscopic distances along the molecular chains, retaining the wave shape, energy, momentum and other properties of the quasiparticle. In this way, the bioenergy can be transported as a localized “wave packet” or soliton. This is just the Davydov’s model of bioenergy transport in proteins, which was proposed in the 1970s [2,3].

Davydov model of bioenergy transport work at α − helical proteins as shown in Figure.1.

Following Davydov idea [3], the Hamiltonian describing such system has in the form of $$\begin{array}{l}{H}_D={\displaystyle \sum _n\left[{\epsilon }_0{B}_n^{+}{B}_n-J(B_n^{+}{B}_{n-1}+B_n{B}_{n+1}^{+})\right]+}{\displaystyle \sum _n\left[\frac{P_n^2}{2M}+\frac{1}{2}w{(u_n-u_{n-1})}^2\right]}\\ +{\displaystyle \sum _n\left[{\chi }_1(u_{n+1}-u_{n-1})B_n^{+}{B}_n\right]=H_{ex}+H_{ph}+H_{\mathrm{int}}}\text{ (1)}\end{array}$$ where ${\epsilon }_0=0.205\text{ev}$ is is the amide-I quantum energy, -J is the dipole-dipole interaction energy between neighboring sites, $B_n^{+}(B_n)$ is the creation (annihilation) operator for an amide-I quantum excitation (exciton) in the site n, u_n is the displacement operator of amino acid residues at site n, P_n is its conjugate momentum operator, M is the mass of an amino acid molecule, w is the elasticity constant of the protein molecular chains, and ${\chi }_1$ is an nonlinear coupling parameter and represents the coupling size of the exciton-phonon interaction. The wave function of the systems proposed by Davydov is in the form of $$|D_2(t)>=|{\phi }_D(t)>1\beta (t)>={\displaystyle \sum _n{\phi }_n}(t){{\rm B}}_n^{+}\exp \left(-\frac{i}{\hslash }{\displaystyle \sum _n[{\beta }_n(t){{\rm P}}_n-{\pi }_n(t){\text{u}}_n]}\right)10>(2)$$ $$|D_1(t)>={\displaystyle \sum _n\{{\phi }_n(t)B_n^{+}\exp \left({\displaystyle \sum _q\left[{\alpha }_{nq}(t)a_q^{+}-{\alpha }_{nq}^{*}(t)a_n\right]}\right)\}10}>\text{ (3)}$$ where I0 > =I0 > _ex Io _ph I0 > _ex and I0 > _ph are the ground states of the exciton and phonon, respectively, $a_q(a_q^{+})$ is annihilation (creation) operator of the phonon with wave vector q, ${\varphi }_n(t)\text{ and }{\beta }_{\text{n}}(t)=<\Phi \left|u_n\right|\Phi >$ and ${\pi }_n(t)=<\left|\Phi \right|P_n|\Phi >$ ${\alpha }_{nq}(t)=<D_1(t)\left|a_q\right|D_1(t)>$ are some undetermined functions of time. Obviously, ${\varphi }_D(t)>={\displaystyle \sum _n{\varphi }_n(t)B_n^{+}10}{>}_{ex}$ in Eq.(2) is an eigenstate of the number operator $\widehat{N}={\displaystyle \sum _n{B}_n^{+}{B}_n},$ corresponding to a single excitation, i.e., $\widehat{N}|{\varphi }_D(t)>=|{\varphi }_D(t)>$.

The Davydov soliton obtained from Eqs. (1)-(2) in the semi classical limit and using the continuum approximation has the from $${\varphi }_D(x,t)={\left(\frac{{\mu }_D}{2}\right)}^{1/2}\text{sech}\left[\frac{{\mu }_D}{r_0}(x-x_0-vt)\right]\exp \left\{i\left[\frac{\hslash v}{2J{r}_0^2}(x-x_0)-E_{\nu }t/\hslash )\right]\right\}\text{ (4)}$$ Corresponding to an excitation localized over a scale r_/ D_μ , where ${\mu }_D=\frac{{\chi }_1^2}{(1-s^2)wJ},{\text{ G}}_{\text{D}}=4J{\mu }_D,{\text{ s}}^{\text{2}}=\frac{v^2}{v_0^2},{\text{ v}}_{\text{0}}=r_{\text{0}}{\text{(w/M)}}^{\text{1/2}}$ is the sound speed in the protein molecular chains, v is the velocity of the soliton, r0 is the lattice constant. Evidently, the soliton contains only one exciton, i.e., ${\varphi }_D(t)\left|\widehat{N}\right|{\varphi }_D(t)>=1$ This shows that the Davydov soliton is formed through self-trapping of one exciton with binding energy E_BD, $E_{BD}=\frac{-{\chi }_1^4}{3J{w}^2}$

Davydov’s idea yields a compelling picture for the mechanism of bioenergy transport in protein molecules and consequently has been the subject of a large number of works [3-28]. A lot of issues related to the Davydov model, including the foundation and accuracy of the theory, the quantum and classical properties and the thermal stability and lifetimes of the Davydov soliton have extensively been studied by many scientists [7-26]. However, considerable controversy has arisen concerning whether the Davydov soliton is sufficiently stable in the region of biological temperature to provide a viable explanation for bio-energy transport. It is out of question that the quantum fluctuations and thermal perturbations are expected to cause the Davydov soliton to decay into a delocalized state. Some numerical simulations indicated that the Davydov soliton is not stable at the biological temperature 300K [7-11, 24-26]. Other simulations showed that the Davydov soliton is stable at 300 K [10-24], but they were based on classical equations of motion which are likely to yield unreliable estimates for the stability of the soliton [3]. The simulations based on the ID₂ state in Eq.(2) generally show that the stability of the soliton decreases with increasing temperatures and that the soliton is not sufficiently stable in the region of biological temperature. Since the dynamical equations used in the simulations are not equivalent to the equation, the stability of the soliton obtained by these numerical simulations is unavailable or unreliable. The simulation[9] based on the ID1>state in Eq. (3) with the thermal treatment of Davydov[3], where the equations of motion are derived from a thermally averaged Hamiltonian, yields the confusing result that the stability of the soliton is enhanced with increasing temperature, predicting that ID1>- type soliton is stable in the region of biological temperature. Evidently, the conclusion is doubtful because the Davydov procedure, in which an equation of motion for an average dynamical state from an average Hamiltonian, corresponding to the Hamiltonian averaged over a thermal distribution of phonons, is inconsistent with standard concepts of quantum-statistical mechanics in which a density matrix must be used to describe the system. Therefore, any exact fully quantum- mechanical treatment for the numerical simulation of the Davydov soliton does not exist. However, for the thermal equilibrium properties of the Davydov soliton, there is a quantum Monte Carlo simulation [13]. In the simulation, correlation characteristic of solitonlike quasiparticles occur only at low temperatures, about T< 10k, for widely accepted parameter values. This is consistent at a qualitative level with the result of Cottingham et al. [15,21]. The latter is a straightforward quantum-mechanical perturbation calculation. The lifetime of the Davydov soliton obtained by using this method is too small (about 10^-12 - 10^-13Sec) to be useful in biological processes. This indicates clearly that the Davydov solution is not a true wave function of the systems. A through study in terms of parameter values, different types of disorder, different thermalization schemes, different wave functions, and different associated dynamics leads to a very complicated picture for the Davydov model [10-12]. These results do not completely rule out the Davydov theory, however they do not eliminate the possibility of another wave function and a more sophisticated Hamiltonian of the system having a soliton with longer lifetimes and good thermal stability.

Indeed, the question of the lifetime of the soliton in protein molecules is twofold. In Langevin dynamics, the problem consists of uncontrolled effects arising from the semiclassical approximation. In quantum treatments, the problem has been the lack of an exact wave function for the soliton. The exact wave function of the fully quantum Davydov model has not been known up to now. Different wave functions have been used to describe the states of the fully quantum-mechanical systems [4,5]. Although some of these wave functions lead to exact quantum states and exact quantum dynamics in the J=0 state, they also share a problem with the original Davydov wave function, namely that the degree of approximation included when $J\ne 0$ is not known. Therefore, it is necessary to reform Davydov’s wave function. Scientists had though that the soliton with a multiquantum $\text{n}\ge \text{2}$ , for example, the coherent state of Brown et al. [4], the multiquantum state of Kerr et al. [12] and Schweitzer et.al. [15,21], the two-quantum state of Cruzeiro-Hansson [18] and Forner [22], and so on, would be thermally stable in the region of biological temperature and could provide a realistic mechanism for bioenergy transport in protein molecules. However, the assumption of the standard coherent state is unsuitable or impossible for biological protein molecules because there are innumerable particles in this state and one could not retain conservation of the number of particles of the system. The assumption of a multiquantum state (n>2) along with a coherent state is also inconsistent with the fact that the bioenergy released in ATP hydrolysis can excite only two quanta of amide-I vibration. On the other hand, the numerical result shows that the soliton of two-quantum state is more stable than that with a one-quantum state.

Cruzeiro-Hansson [18] had thought that Forner’s twoquantum state in the semiclassical case was not exact. Therefore, he constructed again a so-called exactly two-quantum state for the semiclassical Davydov system as follows [18]: $$\text{I}\varphi \left(\text{t}\right)>={\displaystyle \sum _{\text{n},\text{m}=\text{1}}^{\text{N}}{\phi }_{\text{nm}}}\left(\left\{{\text{u}}_{\text{1}}\right\},\left\{{\text{P}}_{\text{1}}\right\},\text{t}\right){\text{B}}_{\text{n}}^{+}{\text{B}}_{\text{m}}^{\text{+}}|\text{0}{>}_{\text{ex}},\text{ (5)}$$ where $B_n(B_n^{+})$ is the annihilation (creation) operator for an amide-I vibration quantum (exciton), $u_1$ is the displacement of the lattice molecules, $P_1$ is its conjugate momentum, and ${|0\rangle }_{ex}$ is the ground state of the exciton. He calculate the average probability distribution of the exciton per site, and average displacement difference per site, and the thermodynamics average of the variable, $P=B_1^{+}{B}_1-B_2^{+}{B}_2$ as a measure of localization of the exciton, versus quantity $\nu =JW/{\chi }_1^2$ and $ln\beta (\beta =1/K_{B}T)$ in the so-called twoquantum state. Eq.(5), where ${\chi }_1$ is a nonlinear coupling parameter related to the interaction of the exciton-phonon in the Davydov model. Their energies and stability are compared with those of the one-quantum state. From the results of above thermal averages, he drew the conclusion that the wave function with a two-quantum state can lead to more stable soliton solutions than that with a one-quantum state, and that the usual Langevin dynamics ,whereby the thermal lifetime of the Davydov soliton is estimated, must be viewed as underestimating the soliton lifetime.

However, by checking carefully Eq.(5), we can find that the Cruzeiro-Hansson wave function[18,24-26] does not represent exactly the two-quantum state. To find out how many quanta the state Eq.(1), indeed contains, the expectation value of the exciton number operator has to be computed. $N={{\displaystyle {\sum }_{n}B}}_n^{+}{B}_n$ , in this state Eq.(5), and sum over the sites, i.e., the exciton numbers N are $$\begin{array}{l}N=<\varphi |{\displaystyle \sum _n{B}_n^{+}}{B}_n|\varphi >={\displaystyle \sum _{ijlmn}{\phi }_{im}^{\ast }}{\phi }_{jl}ex<0|B_i{B}_m{B}_n^{+}{B}_n{B}_j^{+}{B}_l^{+}|0{>}_{ex}(6)\\ ={\displaystyle \sum _{nj}\left({\phi }_{nj}^{\ast }{\phi }_{jn}+{\phi }_{jn}^{\ast }{\phi }_{jn}\right)}+{\displaystyle \sum _{nl}\left({\phi }_{nl}^{\ast }{\phi }_{nl}+{\phi }_{\ln }^{\ast }{\phi }_{nl}\right)}=4\end{array}$$ where we use the relations $$[B_n.B_j^{+}]={\sigma }_{nj},{\displaystyle {\sum }_{nl}{\left|{\phi }_{nl}\right|}^2}=1\text{ (7)}$$ $${}_{ex}\langle 0|B_n^{+}{|0\rangle }_{ex}{=}_{ex}\langle 0|B_n^{+}{B}_n{|0\rangle }_{ex}{=}_{ex}\langle 0|B_n^{+}{B}_n{B}_l{|0\rangle }_{ex}=\mathrm{....}=0(8)$$ Therefore, the state Eq.(5), as it is put forward in Ref. [10], deals with four excitons (quanta), instead of two excitons, in contradiction to the author’s statements. Obviously, it is impossible to create the four excitons by the energy released in the ATP hydrolysis (about 0.43 eV). Thus the author’s wave function is still not relevant to protein molecules, and his discussion and conclusion are all unreliable and implausible in that paper [10].

It is believed that the physical significance of the wave function, Eq. (5), is also unclear, or at least is very difficult to understand. As far as the physical meaning of Eq.(5) is concerned, it represents only a combinational state of single-particle excitation with two quanta created at sites n and m;[18,26] is the probability amplitude of particles occurring at the sites n and m simultaneously. In general, n=≠m and in accordance with the author’s idea. In such a case it is very difficult to imagine the form of the soliton by the mechanism of self- trapping of the two quanta under the action of the nonlinear exciton-phonon interaction, especially when the difference between n and m is very large. Hansson has also not explained the physical and biological reasons and the meaning for the proposed trial state. Therefore, we think that the Cruzeiro-Hansson representation is still not an exact wave function suitable for protein molecules. Thus, the wave function of the systems is still an open problem today.

On the basis of the work of Cruzeio-Hansson [12,13,18,22], Schweitzer [21] and Pang [24-26] proposed a new model of the bioenergy transport in the protein molecules, in which both the Hamiltonian and the wave function of the Dovydov model [24] have been improved. A new coupling interaction between the acoustic and amide-I vibrational modes was added to the original Davydov’s Hamiltonian which takes into account relative displacement of the neighbouring amino acids resulting from dipole-dipole interaction of the neighbouring amide-1 vibrational quanta. Davydov’s wave function has been also replaced with a quasi-coherent two-quanta state to exhibit the coherent behaviors of collective excitations of the excitons and phonons [25-26] which are a feature of the energy released in ATP hydrolysis in the systems. The equation of motion and the properties of the new soliton in the new model are different from those in the Davydov model and as a result the soliton lifetime and stability are greatly enhanced. It is suggested that this model can resolve the controversy on the thermal stability and lifetime of the soliton excited in the protein molecules. The quantum properties of the new soliton will be studied here, but here attention is paid also to the problem of its lifetime and thermal stability at biological temperature 300 K and the lifetime of the new soliton at 300K is calculated in detail by using the generally accepted values of the parameters appropriate to -helical protein molecules in terms of the quantum perturbation theory developed by Cottingham et al. [15], which can take simultaneously into account the quantum and thermal effects. It can be seen that the lifetime of the new soliton at 300 K is long enough to provide a viable explanation of the bio-energy transport in the proteins. The plan of this paper is as follows. In Section 2, the new model, including the extended Hamiltonian and the wave function, is presented. The equations of motion and the new soliton solution in this model are given in Section 3. In Section 4, the properties and thermal stability of the new soliton are discussed, and the possibility of the soliton being a suitable candidate for the mechanism of bioenergy transport in protein molecules is predicted on the basis of results obtained in this paper. In Section 5, the properties of the new soliton are described and its lifetime is calculated by using quantummechanical perturbution method. The detailed discussion of the properties and changes of the lifetimes of the soliton for a large range of parameter values is presented. The conclusions of this investigation are also given in this section.

Results obtained by many scientists over the years indicate that the Davydov model, whether it is the wave function or the Hamiltonian, is indeed too simple, i.e.., it does not denote the elementary properties of the collective excitations occurring in protein molecules, and many improvements of it have been unsuccessful, as mentioned above. What is the source of this problem? It is well known that the Davydov theory on bioenergy transport was introduced into protein molecules from an excitonsoliton model in generally one-dimensional molecular chains [24]. Although the molecular structure of the alpha-helix protein is analogous to some molecular crystals, for example acetanilide (ACN) (in fact, both are polypeptides; the alpha-helix protein molecule is the structure of three peptide channels, ACN is the structure of two peptide channels. If comparing the structure of alpha helix protein with ACN, we find that the hydrogen-boned peptide channels with the atomic structure along the longitudinal direction are the same except for the side group), a lot of properties and functions of the protein molecules are completely different from that of the latter. The protein molecules are both a kinds of soft condensed matter and bio-self-organization with action functions, for instance, self-assembling and self-renovating. The physical concepts of coherence, order, collective effects, and mutual correlation are very important in bio-self- organization, including the protein molecules, when compared with generally molecular systems [25.26]. Therefore, it is worth studying how we can physically describe these properties. It is noted that Davydov operation is not strictly correct. Therefore, it is believed that a basic reason for the failure of the Davydov model is just that it ignores completely the above important properties of the protein molecules.

Let us consider the Davydov model with the present viewpoint. First, as far as the Davydov wave functions, both $|D_1>and|D_2>$ are concerned [3], they are not true solutions of the protein molecules. On the one hand, there is obviously asymmetry in the Davydov wave function since the phononic parts is a coherent state, while the excitonic part is only an excitation state of a single particle. It is not reasonable that the same nonlinear interaction generated by the coupling between the excitons and phonons produces different states for the phonon and exciton. Thus, Davydov’s wave function should be modified [24-26], i.e., the excitonic part in it should also be coherent or quasicoherent to represent the coherent feature of collective excitation in protein molecules. However, the standard coherent [4] and large-n excitation states [12,22] are not appropriate for the protein molecules due to the reasons mentioned above. Similarly, Forner’s and Cruzeiro-Hansson’s two-quantum states do not fulfill the above request. In view of the above discussion, we proposed the following wave function of the protein molecular systems: $$\begin{array}{l}\Phi (t)={|\phi }_P\left(t\right)>|\beta \left(t\right)>=\frac{1}{\lambda }\left[I+{\displaystyle \sum _n{\phi }_n\left(t\right)B_n^{+}+\frac{1}{2!}}{\left({\displaystyle \sum _n{\phi }_n\left(t\right)B_n^{+}}\right)}^2\right]|0{>}_{ex}\times \text{ (9)}\\ \exp \left(-\frac{i}{\hslash }{\displaystyle \sum _n[{\beta }_n(t)P_n-{\pi }_n(t)u_n]}\right){|0\rangle }_{ph}\end{array}$$ where $B_n^{+}$ and $B_n$ are boson creation and annihilation operators for the exciton, $|0{>}_{ex}and|0{>}_{ph}$ are the ground states of the exciton and phonon, respectively $u_n$ and $P_n$ are the displacement and momentum operators of the amino acid residue at site n respectively. The ${\varphi }_n\left(t\right).{\beta }_n\left(t\right)=<\Phi \left(t\right)|u_n|\Phi \left(t\right)>and{\pi }_n\left(t\right)=<\Phi \left(t\right)\left|P_n\right|\Phi \left(T\right)>$ are there sets of unknown functions, $\lambda$ is a normalization constant. It is assumed hereafter that $\lambda =1$ for convenience of calculation, except when explicitly mentioned.

A second problem arises for the Davydov Hamiltonian [24-26,28]. The Davydov Hamiltonian takes into account the resonant or dipole-dipole interaction of the neighboring amide-I vibrational quanta in neighboring amino acid residues with an electrical moment of about 3.5D, but why do we not consider the changes of relative displacement of the neighboring amino acid residues arising from this interaction ? It is reasonable to add the new interaction term ${\chi }_2\left(u_{n+1}-u_n\right)\left(B_{n+1}^{+}{B}_n+B_m^{+}{B}_{n+1}\right)$ into the Davydov’s Hamiltonian to represent correlations of the collective excitations and collective motions in the protein molecules, as mentioned above [24-26]. Although the dipole- dipole interaction is small as compared with the energy of the amide-I vibrational quantum, the change of relative displacement of neighboring peptide groups resulting from this interaction cannot be ignored due to the sensitive dependence of dipole-dipole interaction on the distance between amino acids in the protein molecules, which is a kind of soft condensed matter and bio-self-organization. Thus, the Davydov Hamiltonian is replaced by $$\begin{array}{l}H=H_{ex}+H_{ph}+H_{\mathrm{int}}={\displaystyle \sum _n\left[{\epsilon }_0{B}_n^{+}{B}_n-J\left(B_n^{+}{B}_{n+1}+B_n^{}{B}_{n+1}^{+}\right)\right]}+{\displaystyle {\sum }_n[\frac{P_n^2}{2M}+\frac{1}{2}}w{(}^{u_n}]\\ +{\displaystyle \sum _n\left[{\chi }_1\left(u_{n+1}-u_{n-1}\right)\right]}{B}_n^{+}{B}_n+{\chi }_2\left(u_{n+1}-u_n\right)\times \left(B_{n+1}^{+}{B}_n+B_n^{+}{B}_{n+1}\right)\text{ (10)}\end{array}$$ Where ${\epsilon }_0=0.205\text{ev}$ is the energy of the exciton (the C=0 strechiong mode). The present nonlinear coupling constants are ${\chi }_{1}and{\chi }_2$ They represent the modulations of the on-site energy and resonant (or dipole-dipole) interaction energy of excitons caused by the molecules displacements, respectively .M is the mass of a amino acid molcule and w is the elasticity constant of the protein molecular chains. J is the dipole-dipole interaction energy between neighboring sites. The physical meaning of the other quantities in Eq.(6)are the same as those in the above explanations.

The Hamiltonian and wave function shown in Eqs. (9)- (10) are different from Davydov’s. We add a new interaction term ${\displaystyle {\sum }_n{\chi }_2}\left(u_{n+1}-u_n\right)\left(B_{n+1}^{+}{B}_n+B_n^{+}{B}_{n+1}\right),$ into the original Davydov Hamiltonian. Thus the Hamiltonian now has better correspondence between the interactions and can also represent the features of mutual correlations of the collective excitations and of collective motions in the protein molecules. We should point out here that the different coupling between the relevant modes was also considered by Pang [24-26] and others [27-28] in the Hamiltonian of the vibron-soliton model for one-dimensional oscillator-lattice and protein systems, respectively, but the wave functions of the systems they used are different from Eqs. (9)- (10).

Evidently , the present wave function of the exciton in Eq.(9) is not an excitation state of a single particle, but rather a coherent state, more precisely, a quasicoherent state, because it retain only for three terms of the expansion of a standard coherent state, which can be viewed as an effective truncation of a standard coherent state. When ${\phi }_n\left(t\right)$ is small, i.e., $\left|{\phi }_n\left(t\right)\right|<<1$ Pang represented the wave function of the excitons, $|{\phi }_P\left(t\right)>$ , in Eq.(9) as $$\begin{array}{l}|{\phi }_P\left(t\right)>=\frac{1}{\lambda }\left[1+{\displaystyle \sum _n{\phi }_n\left(t\right)B_n^{+}+\frac{1}{2!}{\left({\displaystyle \sum _n{\phi }_n\left(t\right)B_n^{+}}\right)}^2}\right]|0{>}_{ex}~\frac{1}{\lambda }\exp [-\frac{1}{2}{\displaystyle \sum _n{\left|{\phi }_n(t)\right|}^2}]\times \\ \exp \left\{{\displaystyle \sum _n{\phi }_n\left(t\right)B_n^{+}}\right\}|0{>}_{ex}=\frac{1}{\lambda }\exp \left\{{\displaystyle \sum _n\left[{\phi }_n\left(t\right)B_n^{+}-{\phi }_n^{\ast }\left(t\right)B_n\right]}\right\}|0{>}_{ex}\text{ (11)}\end{array}$$ The last representation in Eq.(11) is a standard coherent state. Therefore, the state of exciton denoted by the new wave function $|\phi \left(t\right)>$ has a coherent feature, thus the wave function inEq.11) is normalized at $\lambda =1$ Since ${{\displaystyle {\sum }_n\left|{\phi }_n\left(t\right)\right|}}^2=1$ required in thecalculation, then this condition of $\left|{\phi }_n\left(t\right)\right|<<1$ naturally satisfied for the roteins consisting of several hundreds of amino acids. Just so, the wave function denoted in Eq.(9) represents exactly the coherent features of collective excitations of both the excitons and phonons caused by the nonlinear exciton - phonon interaction, this indicates that the wave function in Eq.(9) is justified for the proteins. However, it is not an eigenstate of the number operator, $\widehat{N}={\displaystyle {\sum }_n{B}_n^{+}{B}_n}$ because of $$\widehat{N}|{\phi }_P>={\displaystyle \sum _n{B}_n^{+}}{B}_n|{\phi }_P\rangle =\left\{{\displaystyle \sum _n{\phi }_n\left(t\right)B_n^{+}}+{\left({\displaystyle \sum _n{\phi }_n\left(t\right)B_n^{+}}\right)}^2\right\}|0{>}_{ex}=2|{\phi }_P\rangle -\left(2+{\displaystyle \sum _n{\phi }_n\left(t\right)B_n^{+}}\right)|0{>}_{ex}(12)$$ Therefore, the $|{\phi }_P\rangle$ represents a coherent superposition of the excitonic state with two quanta and the ground state of the exciton, but it has a determinate numbers of quanta. From the expectation value of number operator $\widehat{N}$ we find that this state contains the number of exciton is $$\begin{array}{l}N=<{\phi }_P\left|\widehat{N}\right|{\phi }_P>={\displaystyle \sum _n\langle {\phi }_P|B_n^{+}{B}_n|{\phi }_P\rangle }=\left\{{\displaystyle \sum _n{\left|{\phi }_n\left(t\right)\right|}^2+\left({\displaystyle \sum _n{\left|{\phi }_n\left(t\right)\right|}^2}\right)\left({\displaystyle \sum _m{\left|{\phi }_m\left(t\right)\right|}^2}\right)}\right\}\\ =\left({\displaystyle \sum _n{\left|{\phi }_n\left(t\right)\right|}^2}\right)\left(1+{\displaystyle \sum _m{\left|{\phi }_m\left(t\right)\right|}^2}\right)=2\text{ (13)}\end{array}$$ namely, it contains only two quanta. Where we utilize Eq.(8) and the following relation [24] is: $${\displaystyle \sum _n{\left|{\varphi }_n\left(t\right)\right|}^2}=1,{\displaystyle \sum _m{\left|{\varphi }_m\left(t\right)\right|}^2}=1,[B_n.B_m^{+}]={\delta }_{nm}\text{ (14)}$$ $$\begin{array}{l}{}_{ex}<0\left|B_n^{+}\right|0{>}_{ex}{=}_{ex}<0\left|B_n^{+}{B}_n\right|0{>}_{ex}{=}_{ex}<0\left|B_n^{+}{B}_m\right|0{>}_{ex}\\ {=}_{ex}<0\left|B_n^{+}{B}_m{B}_l\right|0{>}_{ex}{=}_{ex}<0\left|B_n^{+}{B}_m{B}_l^{+}{B}_n\right|0{>}_{ex}\\ {=}_{ex}<0\left|B_n^{+}{B}_m{B}_i^{+}{B}_l{B}_j\right|0{>}_{ex}{=}_{ex}<0\left|B_n^{+}{B}_m{B}_l^{+}{B}_i{B}_j{B}_n\right|0{>}_{ex}\mathrm{....}=0\end{array}$$ Therefore, the new wave function is completely different from Davydov’s. The latter is an excitation state of a single particle with one quantum and an eigenstate of the number operator, but the former is not. The new state is a quasicoherent state. It contains only two excitons, which come from the second and third terms in Eq.(9), in which each term contributes only an exciton, but it is not an excitation state of two single parties. Hence, as far as the form of new wave function in Eq.(9) is concerned, it is either two-quanta states proposed by Forner [21] and Cruzeiro- Hansson [10,18] or a standard coherent state proposed by Brown et al. [4,2] and Kerr et al’s [13] and Schweitzer et al’s 15,[21] multiquanta states. Therefore, the wave function, Eq. (9), is new for the protein molecular systems. It not only exhibits the coherent feature of the collective excitation of excitons and phonons caused by the nonlinear interaction generated by the exciton-phonon interaction, which , thus, also makes the wave function of the states of the system symmetrical, but it also agrees with the fact that the energy released in the ATP hydrolysis (about 0.43 eV) may only create two amide-I vibrational quanta which, thus, can also make the numbers of excitons maintain conservation in the Hamiltonian, Eq.(10). Meanwhile, the new wave function has another advantage, i.e., the equation of motion of the soliton can also be obtained from the Heisenberg equations of the creation and annihilation operators in quantum mechanics by using Eqs. (9) and (10), but the wave function of the states of the system in other models could not, including the one-quanta state [3] and the two-quanta state [12,22]. Therefore, the above Hamitonian and wave function, Eqs. (9) and (10), are reasonable and appropriate to the protein molecules.

We now derive the equations of motion from Pang’s model. First of all, we give the interpretation of ${\beta }_n(t)$ and ${\pi }_n(t)$ in Eq.(9).We know that the phonon part of the new wave function in Eq.(9) depending on the displacement and momentum operators is a coherent state of the normal model of creation and annihilation operators. A coherent state for the mode with wave vector $\text{q}$ is [3,12,24-26] $$|\alpha (t)\rangle =\exp \left({\displaystyle \sum _q[{\alpha }_q(t)a_q^{+}-{\alpha }_q^{*}(t)a_q]}\right){|0\rangle }_{ph}\text{ (15)}$$ Utilizing the standard transformations $$u_n={\displaystyle \sum _q{\left[\frac{\hslash }{2NM{\omega }_q}\right]}^{1/2}{e}^{iqn{r}_0}(a_{-q}^{+}+a_q)},P_n=i{\displaystyle \sum _q{\left[\frac{M\hslash {\omega }_q}{2N}\right]}^{1/2}{e}^{iqn{r}_0}(a_{-q}^{+}-a_q)}\text{ (16)}$$ we can get [12,23] $|\alpha (t)\rangle =|\beta (t)\rangle$ ,where $|\beta (t)\rangle$ is in Eq.(9), and ${\omega }_q=2{(w/M)}^{1/2}\sin (r_{0}q/2)$ is the distance between neighboring amino acid molecules, and ${\text{a}}_{\text{q}}({\text{a}}_{\text{q}}^{+})$ is the annihilation (creation) operator of the phonon with wave vector $\text{q}$ ,where $$\langle \alpha (\text{t})\left|{\text{a}}_{\text{q}}\right|\alpha (\text{t})\rangle ={\text{a}}_{\text{q}}(\text{t})={\left(\frac{\text{M}{\omega }_{\text{q}}}{\text{2h}}\right)}^{\text{1}/\text{2}}{\beta }_{\text{q}}(\text{t})+\text{i}{\left(\frac{\text{1}}{\text{2M}\hslash {\omega }_{\text{q}}}\right)}^{\text{1}/\text{2}}{\pi }_{\text{q}}(\text{t})$$ $${\beta }_{\text{q}}(\text{t})=\frac{\text{1}}{\sqrt{\text{N}}}{\displaystyle \sum _{\text{n}}{\text{e}}^{-{\text{iqnr}}_{\text{0}}}{\pi }_{\text{n}}(\text{t})},{\pi }_{\text{q}}(\text{t})=\frac{\text{1}}{\sqrt{\text{N}}}{\displaystyle \sum _{\text{n}}{\text{e}}^{-{\text{iqnr}}_{\text{0}}}{\pi }_{\text{n}}(\text{t})},\langle \Phi (\text{t})\left|{\text{P}}_{\text{n}}\right|\Phi (\text{t})\rangle ={\pi }_{\text{n}}(\text{t})\text{ (17)}$$ Utilizing again the above results and the formulas of the expectation values of the Heisenberg equations of operators $u_n$ and $P_n$ , in the state $|\Phi (t)\rangle$. $$i\hslash \frac{\partial }{\partial t}\langle \Phi (\text{t})\left|{\text{u}}_{\text{n}}\right|\Phi (\text{t})\rangle =\langle \Phi (\text{t})\left|\left[{\text{u}}_{\text{n}},H\right]\right|\Phi (\text{t})\rangle ,i\hslash \frac{\partial }{\partial t}\langle \Phi (\text{t})\left|{\text{P}}_{\text{n}}\right|\Phi (\text{t})\rangle =\langle \Phi (\text{t})\left|\left[{\text{u}}_{\text{n}},H\right]\right|\Phi (\text{t})\rangle \text{ (18)}$$ We can obtain the equation of motion for the ${\beta }_n(t)$ as $$M{\ddot{\beta }}_n(t)=w\left[{\beta }_{n+1}(t)-2{\beta }_n(t)+{\beta }_{n-1}(t)\right]+2{\chi }_1[{\left|{\varphi }_{n+1}(t)\right|}^2-{\left|{\varphi }_{n-1}(t)\right|}^2]+$$ $$2{\chi }_2\{{\varphi }_n^{*}(t)\left[{\varphi }_{n+1}(t)-{\varphi }_{n-1}(t)\right]+{\varphi }_n(t)\left[{\varphi }_{n+1}^{*}(t)-{\varphi }_{n-1}^{*}(t)\right]\}\text{ (19)}$$ From Eq.(19) we see that the presence of two quanta for the oscillators increases the driving force on the phonon field by that factor, when compared with the Davydov theory.

We now derive the equation of motion for the${\phi }_n$ . A basic assumption in the derivation is that $|\Phi (t)\rangle$ in Eq. (9) is a solution of the time-dependent equation [24-26]: $$\text{i}\hslash \frac{\partial }{\partial \text{t}}|\Phi (\text{t})\rangle =\text{H}|\Phi (\text{t})\rangle \text{ (20)}$$ The left-hand side of Eq. (16) has [12,23] $$i\hslash \frac{\partial }{\partial t}|\Phi (t)\rangle =\left\{i\hslash \left({\displaystyle \sum _n{\dot{\varphi }}_n(t)B_n^{+}+{\displaystyle \sum _n{\dot{\varphi }}_n(t){\varphi }_n(t)B_n^{+}{B}_n^{+}{|0\rangle }_{ex}}}\right)\right\}|\beta (t)\rangle$$ $$+|{\phi }_P(t)\rangle \left\{{\displaystyle \sum _n\left\{{\beta }_n(t)P_n-{\pi }_n(t)u_n+\frac{1}{2}\left[{\beta }_n(t){\dot{\pi }}_n(t)-{\dot{\beta }}_n(t){\pi }_n(t)\right]\right\}|\beta (t)\rangle }\right\}\text{ (21)}$$ Now left-multiplying the both sides of Eq.(21) by $\langle \Phi (\text{t})|$ ,the lefthand side of Eq.(21) can be $$\text{i}\hslash \langle \Phi (\text{t})\left|{\text{u}}_{\text{n}}\right|\Phi (\text{t})\rangle \text{i}\hslash {\displaystyle \sum _{\text{n}}{\varphi }_{\text{n}}^{*}(\text{t}){\varphi }_{\text{n}}(\text{t})}\left({\displaystyle \sum _{\text{m}}{\varphi }_{\text{m}}^{*}(\text{t}){\varphi }_{\text{m}}(\text{t})+\text{1}}\right)$$ $$=\frac{\text{5}}{\text{4}}{\displaystyle \sum _{\text{n}}\left[{\dot{\beta }}_n(t){\pi }_n(t)-{\dot{\pi }}_n(t){\beta }_n(t)\right]}{\displaystyle \sum _{\text{n}}{\left|{\varphi }_{\text{n}}(\text{t})\right|}^{\text{2}}}\text{ (22)}$$ Similarly, for the right-hand side of Eq. (20) we can have [12,23] $$\langle \Phi (\text{t})\left|\left({\text{H}}_{\text{ex}}+{\text{H}}_{\text{ph}}+{\text{H}}_{\text{int}}\right)\right|\Phi (\text{t})\rangle =\{{\displaystyle \sum _{\text{n}}\left\{{\epsilon }_{\text{0}}\left|{\varphi }_{\text{n}}{(\text{t})}^{\text{2}}\right|-\text{J}{\varphi }_{\text{n}}^{*}(\text{t})\left[{\varphi }_{\text{n}+\text{1}}(\text{t})-{\varphi }_{\text{n}-\text{1}}(\text{t})\right]\right\}}\times \left(1+{\displaystyle \sum _m{\left|{\varphi }_m(t)\right|}^2}\right)$$ $$+\{{\displaystyle \sum _n\{{\chi }_1\left[{\beta }_{n+1}(t)-{\beta }_{n-1}(t)\right]{\left|{\varphi }_n(t)\right|}^2+{\chi }_2\left[{\beta }_{n+1}(t)-{\beta }_{n-1}(t)\right]}$$ $$\times {\phi }_n^{*}(t)\left[{\phi }_{n+1}(t)-{\phi }_{n-1}(t)\right]\}\}\left(1+{\displaystyle \sum _m{\left|{\phi }_m(t)\right|}^2}\right)+\frac{5}{2}W(t){\displaystyle \sum _n{\left|{\phi }_n(t)\right|}^2}\text{ (23)}$$ where $$W(t)=\langle \beta (\text{t})\left|H_{ph}\right|\beta (\text{t})\rangle ={\displaystyle \sum _n\left(\frac{1}{2M}{\pi }_n^2(t)+\frac{1}{2}w{\left[{\beta }_n(t)-{\beta }_{n-1}(t)\right]}^2\right)}+{\displaystyle \sum _q\frac{1}{2}\hslash {\omega }_q}\text{ (24)}$$ and utilizing Eqs.(8) and (12)-(14) and the relationships can be obtained: $${\displaystyle \sum _{\text{n}}\left[{\beta }_{\text{m}+1}(t)-2{\beta }_{\text{m}}(t)+{\beta }_{\text{m}-1}(t)\right]{\beta }_m(t)}=-{\displaystyle \sum _{\text{n}}{\left[{\beta }_{\text{m}+1}(t)-{\beta }_{\text{m}-1}(t)\right]}^2}$$ $$\langle \Phi (\text{t})\left|{\displaystyle \sum _n\left(B_n^{+}{B}_{n-1}+B_n{B}_{n-1}^{+}\right)}\right|\Phi (\text{t})\rangle ={\displaystyle \sum _n\left[{\varphi }_n^{*}(t){\varphi }_{n+1}(t)+{\varphi }_{n-1}^{*}(t)\varphi (t)\right]\left(1+{\displaystyle \sum _m{\left|{\varphi }_m(t)\right|}^2}\right)}$$ $$\langle \Phi (\text{t})\left|{\displaystyle \sum _n\left(u_{n+1}-u_{n-1}\right)\left(B_n^{+}{B}_n\right)}\right|\Phi (\text{t})\rangle ={\displaystyle \sum _n\left\{\left[{\beta }_{\text{m}+1}(t)-{\beta }_{\text{m}-1}(t)\right]{\left|{\varphi }_n(t)\right|}^2\right\}\left(1+{\displaystyle \sum _m{\left|{\varphi }_m(t)\right|}^2}\right)}$$ $$\begin{array}{l}\langle \Phi (\text{t})\left|{\displaystyle \sum _n\left(u_{n-1}-u_n\right)\left(B_n^{+}{B}_{n-1}+B_n{B}_{n-1}^{+}\right)}\right|\Phi (\text{t})\rangle =\\ {\displaystyle \sum _n\left\{\left[{\beta }_{\text{m}+1}(t)-{\beta }_{\text{m}-1}(t)\right]\left[{\varphi }_n^{*}(t){\varphi }_{n+1}(t)+{\varphi }_{n-1}^{*}(t)\varphi (t)\right]\right\}}\times \left(1+{\displaystyle \sum _m{\left|{\varphi }_m(t)\right|}^2}\right)\text{ (25)}\end{array}$$ From Eqs.(20)-(23) we can obtain $$i\hslash \frac{\partial }{\partial t}{\varphi }_n(t)={\epsilon }_0{\varphi }_n(t)-J\left[{\varphi }_{n+1}(t)+{\varphi }_{n-1}(t)\right]+{\chi }_1\left[{\beta }_{n+1}(t)+{\beta }_{n-1}(t)\right]{\varphi }_n(t)$$ $$-{\chi }_2\left[{\beta }_{n+1}(t)+{\beta }_n(t)\right]\times \left[{\varphi }_{n+1}(t)+{\varphi }_{n-1}(t)\right]$$ $$+\frac{5}{2}\left(W(t)-\frac{1}{2}{\displaystyle \sum _m\left[{\dot{\beta }}_m(t){\pi }_m(t)-{\dot{\pi }}_m(t)\beta (t)\right]}\right){\varphi }_n(t)\text{ (26)}$$ In the continuum approximation we get from Eqs.(19) and (26) $$i\hslash \frac{\partial }{\partial t}\varphi (x,t)=R(t)\varphi (x,t)-J{r}_0^2\frac{{\partial }^2}{\partial x^2}\varphi (x,t)-G_p{\left|\varphi (x,t)\right|}^2\varphi (x,t)\text{ (27)}$$ and $$\frac{\partial \beta (x,t)}{\partial \xi }=\frac{\partial \beta (x,t)}{\partial x}=-\frac{4\left({\chi }_1+{\chi }_1\right)}{w\left(1-s^2\right)r_0}{\left|\varphi (x,t)\right|}^2\text{ (28)}$$ here $\xi =x-v_t$ $R(t)={\epsilon }_0-2J+\frac{5}{2}\left\{W(t)-\frac{1}{2}{\displaystyle \sum _m\left[{\dot{\beta }}_m(t){\pi }_m(t)-{\dot{\pi }}_m(t)\beta (t)\right]}\right\}$ and $s=v/v{}_0$ . The soliton solution of Eq.(27) is thus $$\varphi (x,t)={\left(\frac{{\mu }_p}{2}\right)}^{1/2}\sec h\left[\left({\mu }_p/r_0\right)\left(x-x_0-V_t\right)\right]\times \exp \left\{i\left[\frac{\hslash v}{2J{r}_0^2}\left(x-x_0\right)-E_v\frac{t}{\hslash }\right]\right\}\text{ (29)}$$ with $${\mu }_P=\frac{2{\left({\chi }_1+{\chi }_2\right)}^2}{w\left(1-s^2\right)J},G_P=\frac{8{\left({\chi }_1+{\chi }_2\right)}^2}{w\left(1-s^2\right)}\text{ (30)}$$ Although forms of the above equations of motion and the corresponding solution, Eqs.(27)- (30), are quite similar to that of the Davydov soliton, the properties of new soliton have very large differences from the latter because the parameter values in the equation of motion and the solution Eqs.(27) and (29), including $R(t),G_P,$ and ${\mu }_P$ have obvious distinctions from that in the Davydov model. A straightforward result of Pang’s model is to increase the nonlinear interaction energy $G_P\left(G_P=2G_D\left[1+2\left({\chi }_2/{\chi }_1\right)+{\left({\chi }_2/{\chi }_1\right)}^2\right]\right)$ and the amplitude of the new soliton and decrease its width due to an increase of [24-26] when compared with Davydov soliton [3], where ${\mu }_D=x_1^2/w(1-s^2)J$ , and $G_D=4x_1^2/w(1-s^2)$ are the corresponding values in the Davydov mode [3-8]. Thus the localized feature of the new soliton is enhanced. Therefore its stability against the quantum fluctuation and thermal perturbations increased considerably as compared with the Davydov soliton.

The energy of soliton in Pang’s model becomes [24-26] $$E=<\Phi (t)\left|H\right|\Phi (t)=\frac{1}{r_0}{\displaystyle {\int }_{-\infty }^{\infty }2[J{r}_0^2}{\left(\frac{\partial \varphi }{\partial x}\right)}^2+R{\left|\varphi (x,t)\right|}^2-G_p{\left|\varphi (x,t)\right|}^{4}dx+$$ $$\frac{1}{r_0}{\displaystyle {\int }_{-\infty }^{\infty }\frac{1}{2}\left[M{\left(\frac{\partial \beta (x,t)}{\partial t}\right)}^2+w{r}_0^{}{\left(\frac{\partial \beta (x,t)}{\partial x}\right)}^2\right]}dx=E_0+\frac{1}{2}{M}_{sol}{v}^2\text{ (31)}$$ The rest energy of the new soliton is $$E_0=2({\epsilon }_0-2J)-\frac{8{(x_1+x_2)}^4}{3w^{2}J}=E_s^0-W\text{ (32)}$$ where $W=[2{(x_1+x_2)}^4]/3w^{2}J$ is the energy of deformation of the lattice. The effective mass of the new soliton is $$M_{sol}=2m_{ex}+\frac{8{\left(x_1+x_2\right)}^4\left(9s^2+2-3s^4\right)}{3w^{2}J{\left(1-s^2\right)}^3{v}_0^2}\text{ (33)}$$ We utilize Eqs.(8) and (12)-(14) in the above calculations. In such a case, the binding energy of the new soliton is $$E_{BP}=\frac{-8{\left(x_1+x_2\right)}^4}{3J{w}^2}\text{ (34)}$$ $E_{BP}$ is larger than that of the Davydov soliton. The latter is $E_{BD}=-x_1^4/3J{w}^2$ .They have the following relation: $$E_{BP}=8E_{BD}\left[1+4\left(\frac{x_2}{x_1}\right)+6{\left(\frac{x_2}{x_1}\right)}^2+4{\left(\frac{x_2}{x_1}\right)}^3+{\left(\frac{x_2}{x_1}\right)}^4\right]\text{ (35)}$$ We can estimate that the binding energy of the new soliton is about several decades larger than that of the Davydov soliton .This is a very interesting result. It is helpful to enhance thermal stability of the new soliton. Obviously, the increase of the binding energy of the new soliton comes from its two-quanta nature and the added interaction ${\displaystyle {\sum }_i{\chi }_2\left(u_{n+1}-u_n\right)\left(B_{n+1}^{+}{B}_n+B_n^{+}{B}_{n+1}^{}\right)}$ in the Hamiltonian of the systems, Eq.(10). However, we see from Eq.(35) that the former plays the main role in the increase of the binding energy and the enhancement of thermal stability for the new soliton relative to the latter due to ${\chi }_2<{\chi }_1$ . The increase of the binding energy results in significant changes of properties of the new soliton, which are discussed as follows.

In comparing various correlations to this model, it is helpful to consider them as a function of a composite coupling parameter like that of Pouthier and Spatchek et al. [32-35] and Scott [6] again, it is convenient to define another composite parameter [3,24-26] that can be written as $$4\pi {\alpha }_P={\left({\chi }_1+{\chi }_2\right)}^2/2w\hslash {\omega }_D\text{ (36)}$$ where ${\omega }_D={\left(w/M\right)}^{1/2}$ is the band edge for acoustic phonons (Debye frequency). If, $4\pi {\alpha }_P>1$ it is said to be weak. Using widely accepted values for the physical parameters for the alpha helix protein molecule [2-23],$$J=1.55\times {10}^{-22}J.w=(13-19.5)N/m.M=(1.17-1.91)\times {10}^{-25}kg$$ $${\chi }_1=62\times {10}^{-12}N.{\chi }_2=(10-18)\times {10}^{-12}N.r_0=4.5\times {10}^{-10}m.\text{ (37)}$$ We can estimate that the coupled constant lies in the region of $4\pi {\alpha }_P=0.11-0.273$ , which is not a weakly coupled theory, the coupling strength is enhanced as compared with the Davydov model, the latter is $4\pi {\alpha }_D=0.036-0.045$ Using the notation of Bullough et al. [29,30], Teki et al. [31,32], and Pouthier et al. [33- 35] $$\gamma =J/2\hslash w_D\text{ (38)}$$ In terms of the two composite parameters, $4\pi {\alpha }_P$ and $\gamma$ ,the soliton binding energy for Pang’s model can be written by $$E_{BP}/J=8{\left(4\pi {\alpha }_P/\gamma \right)}^2/3,M_{sol}=2m_{ex}\left[1+32{\left(4\pi {\alpha }_P\right)}^2/3\right]\text{ (39)}$$ From the above parameter values, we find $\gamma =0.08$ Utilizing this value, the $E_{BP}/J$ versus $4\pi \alpha$ relations in Eq.(39) are plotted in Fig.1. However, $E_{BP}/J={\left(4\pi {\alpha }_P/\gamma \right)}^2/3$ for the Davydov model (here $M_{sol}^{\text{'}}=m_{ex}\left[1+2{\left(4\pi {\alpha }_P/\gamma \right)}^2/3\right],4\pi {\alpha }_D={\chi }_1^2/2w\hslash {\omega }_D),$ then the $E_{BD}/J$ versus $4\pi {\alpha }_D$ relation is also plotted in figure 2 from this figure we see that the difference of soliton binding energies between two models becomes larger with increasing $4\pi \alpha$ [24-26].

Also, we see clearly from Eqs. (28)-(32) and (35) that the localized feature of the new soliton is enhanced due to increases of the nonlinear interaction and of the binding energy of the new soliton resulting from the increases of exciton-phonon interaction in Pang’s model. Thus, the stability of the soliton against quantum and thermal fluctuations is also enhanced considerately [24-26].

As a matter of fact, the nonlinear interaction energy forming the new soliton in Pang’s model is ${\text{G}}_{\text{P}}=\text{8}{\left({\chi }_{\text{1}}+{\chi }_{\text{2}}\right)}^{\text{2}}/\left(\text{1}-{\text{s}}^{\text{2}}\right)w=3.8\times {10}^{-32}J$ ,and it is larger than the linear dispersion energy $\text{J}=\text{1}.\text{55}\times {\text{10}}^{-\text{32}}\text{J}$ ,i.e., the nonlinear interaction in Pang’s model is so large that it can actually cancel or suppress the linear dispersion effect in the equation of motion ,thus the new soliton is stable in such a case according the soliton theory [2,33-35]. On the other hand, the nonlinear interaction energy in the Davydov model is only ${\text{G}}_{\text{D}}=\text{4}{\chi }_{\text{1}}{}^{\text{2}}/\left(\text{1}-{\text{S}}^{\text{2}}\right)w=1.8\times {10}^{-21}J$ and it is about three to four times smaller than ${\text{G}}_{\text{P}}$ .Therefore, the stability of the Davydov soliton is weaker as compared with the new soliton. Moreover, the binding energy of the new soliton in Pang’s model is $E_{BP}=\left(4.16-4.3\right)\times {10}^{-21}J$ in Eq.(31), which is somewhat larger than the thermal perturbation energy, $k_{B}T=4.13\times {10}^{-21}J$ at $300K$ and about four times larger than the Debye energy $k\Theta =\hslash {\omega }_D=1.2\times {10}^{-21}J$ (there ${\omega }_D$ is the Debye frequency). This shows that transition of the new soliton to a delocalized state can be suppressed by the large energy difference between the initial (solitonic) state and final (delocalized) state, which is very difficult to compensate with the energy of the absorbed phonon. Thus ,the new soliton is robust against quantum fluctuations and thermal perturbation, therefore it has a large lifetime and good thermal stability in the region of biological temperature .In practice, according to Schweitzer et al. s studies (i.e the lifetime of the soliton increases as ${\mu }_p$ and $T_0=\hslash V_0{\mu }_p/K_B\pi$ increase at a given temperature)[15] and the above obtained results, an inference could roughly be drawn that the lifetime of the new soliton will increase considerably as compared with that of the Davydov soliton due to the increase of ${\mu }_p$ and $T_0$ because the latter are about three times larger than that of the Davydov model. On the other hand, the binding energy of the Davydov soliton ${\text{E}}_{\text{BD}}={\chi }_{{}_{\text{1}}}^{\text{4}}/{\text{3w}}^{\text{2}}\text{J}=0.188\times {10}^{-21}J$ and it is about 23 times smaller than that of the new soliton, about 22 times smaller than $K_{B}T$ , and about 6 times smaller than $K_B\Theta$ , respectively. Therefore, the Davydov soliton is easily destructed by the thermal perturbation energy and quantum transition effects. Thus it indicates that the Davdov soliton has a very small lifetime, and it is unstable at the biological temperature $300K$ .This conclusion is consistent at a qualitative level with the result s of Wang et al [13,14] and Cottingham et al.[15,21].

One can sum up the differences between Pang’s model and Davydov’s model, Eqs.(1)-(4), as follows. First, the parameter μp is increased (μp =2μD $\left[1+2(\frac{{\chi }_2}{{\chi }_1})+{(\frac{{\chi }_2}{{\chi }_1})}^2\right],).$ Secondly the nonlinear coupling energy becomes Gp = ${\frac{8({\chi }_1+{\chi }_2)}{w(1-s^2)}}^2$ (G_p=2G_D $\left[1+2(\frac{{\chi }_2}{{\chi }_1})+{(\frac{{\chi }_2}{{\chi }_1})}^2\right]$ , where $G_D=\frac{4{\chi }_1^2}{w(1-s^2)}$ is the nonlinear interaction in the Davydov model, resulting from the two-quanta nature and the enhancement of the coupling the coefficient $({\chi }_1+{\chi }_2)$ . For $\alpha$ helical protein molecules, and using the parameter values listed in Eq. (37) the values of the main parameters in Pang’s model can be calculated. These values and the corresponding values in the Davydov model are simultaneously listed in table 1.

From table 1 we can see clearly that the new model produces considerable changes in the properties of the new soliton, such as large increase of the nonlinear interaction, binding energy and amplitude of the soliton, and decrease of its width as compared to that of the Davydov soliton. This indicates that the soliton in Pang’s model is more localized and more robust against quantum and the thermal stability has been enhanced [2,27,28] which implies an increase in lifetime for the new soliton. From Eq. (19) it can also be found that the effect of the two-quanta nature is larger than that of the added interaction. We can therefore refer to the new soliton as quasi-coherent.

In the above studies, the influences of quantum and thermal effects on soliton state, which are expected to cause the soliton to decay into delocalized states, we postulate that the model Hamiltonian and the wave function in Pang’s model together give a complete and realistic picture of the interaction properties and allowed states of the protein molecules. The additional interaction term in the Hamiltonian gives better symmetry of interactions. The new wave function is a reasonable choice for the protein molecules because it not only can exhibit the coherent features of collective excitations arising from the nonlinear interaction between the excitons and phonons, but also retain the conservation of number of particles and fulfill the fact that the energy released by the ATP hydrolysis can only excite two quanta. In such a case, using a standard calculating method [2,26]and widely accepted parameters we can calculate the region encompassed of the excitation or the linear extent of the new soliton, $\Delta X=2\pi r_0/{\mu }_p$ to be greater than the lattice constant r0, i.e., $\Delta X>r_0$ as shown in table 1. Conversely, we can explicitly calculate the amplitude squared of the new soliton using Eq. (29) in its rest frame as $|\phi (x){|}^2=\frac{{\mu }_p}{2}\sec h^2(\frac{{\mu }_{p}x}{r_0})$ Thus the probability to find the new soliton outside a range of width r0 is about 0.10. This number can be compatible with the continuous approximation since the quasi-coherent soliton can spread over more than one lattice spacing in the system in such a case. This proves that assuming the continuous approximation used in the calculation is valid. Therefore we should believe that the above calculated results obtained from Pang’s model is correct.

The lifetime of the soliton in the protein molecules is an centre problem in the process of bioenergy transport because the soliton possess certain biological meanings and can play an important role in the biological process, only if it has enough long lifetimes. Therefore, to calculate the lifetime of the new soliton in Pang’s model has important significance.

For convenlence of calculation, we here represent the wave function of the system in Eq. (9) by [24-26] $$|\Phi (t)>=|\varphi (t)>|\beta (t)>={\text{U}}_{\text{1}}\text{10}{＞}_{\text{ex}}{\text{U}}_{\text{2}}\text{10}{＞}_{\text{ph}},\text{ (40)}$$ where $${\text{U}}_{\text{1}}=\frac{1}{\lambda }[1+{\displaystyle \sum _n{\varphi }_n(t)B_n^{+}+}\frac{1}{2!}({\displaystyle \sum _n{\varphi }_n(t)B_n^{+}{)}^2}]\text{ (40a)}$$ $${\text{U}}_{\text{2}}=\exp \left\{-\frac{i}{\hslash }{\displaystyle \sum _n[{\beta }_n(t)P_n-{\pi }_n(t)u_n}]\right\}\text{ (40b)}$$ $$=\exp \left\{\frac{1}{\sqrt{N}}{\displaystyle \sum _q{\alpha }_q(t)a_q^{+}-{\alpha }_q^{*}(t)a_q}\right\}\text{ (40c)}$$ where we assume ${{\displaystyle \sum _i\left|{\varphi }_i\right|}}^2=n$ ,where n is an integer, denotes the number of particle. The wave function, Eq.(40), does not only exhibit coherent properties, but also agrees with the fact that the energy released in the ATP hydrolysis (about 0.43eV) excites only two amide-I vibrational quanta, instead of multiquanta (n>2)[24- 26]. Therefore, the Hamitonian and wave function of the systems, Eqs. (9)-(10), or (40) are reasonable and appropriate to the protein molecules. Using the standard transformation in Eq.(16), where ${\omega }_q=2{(w/M)}^{1/2}\sin (\frac{r_{0}q}{2})$ ,Eq.(10) becomes $$\begin{array}{l}H={\displaystyle \sum _n\left[{\epsilon }_0{B}_n^{+}{B}_n^{}-J(B_n^{+}{B}_{n+1}^{}+B_{n+1}^{+}{B}_n^{})\right]}+{\displaystyle \sum _q\hslash {\omega }_q(a_q^{+}{a}_q+\frac{1}{2})}\\ +\frac{1}{\sqrt{N}}{\displaystyle \sum _{q.n}[g_1(q)B_n^{+}{B}_n^{}+g_2(q)(B_n^{+}{B}_{n+1}^{}+B_n^{+}{B}_{n+1}^{})](a_q+a_{-q}^{+})e^{in{r}_0^{}q}}\text{ (41)}\end{array}$$ where $$g_1(q)=2{\chi }_{1}i{\left[\frac{\hslash }{2M{\omega }_q}\right]}^{1/2}\sin r_{0}q;g_2(q)={\chi }_2{\left[\frac{\hslash }{2M{\omega }_q}\right]}^{1/2}(e^{i{r}_{0}q}-1)\text{ (42)}$$ In a semi classical and continuum approximations, from Eq.(41) we can obtain the envelope soliton solution Eq.(29) in Pang’s model, we now represent Eq.(29) by the following form[24-26] $$\varphi (x,t)={\left(\frac{{\mu }_p}{2}\right)}^{1/2}Sech[\frac{{\mu }_p}{r_0}(x-vt)]\exp \left[\frac{i}{\hslash }\left(\frac{{\hslash }^{2}vx}{2J{r}_0^2}-E_{sol}t\right)\right]\text{ (43)}$$ where $${\mu }_p=\frac{2{({\chi }_1+{\chi }_2)}^2}{w(1-s^2)J}\text{ (44)}$$ The energy of the new soliton is $${\alpha }_q(t)=\frac{i\pi ({\chi }_1+{\chi }_2)}{w{\mu }_p(1-v^2/v_0^2)}{\left[\frac{M}{2\hslash {\omega }_q}\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}({\omega }_q+qv)\csc h(\pi q{r}_0/2{\mu }_p)e^{iqvt}={\alpha }_q{e}^{iqvt}\text{ (46)}$$ This treatment yields a localized coherent structure of the excitons with size of order 2 r0/μp that propagates with velocity v and can transfer energy ES01< . Unlike bare excitons that are scattered by the interactions with the phonons, this soliton state describes a quasi-particle consisting of the two excitons plus a lattice deformation and hence a priori includes interaction with the acoustic phonons. So the soliton is not scattered and spread by this interaction of the vibration of amino acids (lattices), and can maintain its form, energy, momentum and other quasiparticle properties moving over a macroscopic distance. The bell-shaped form of the soliton Eq. (43) does not depend on the excitation method. It is self-consistent. Since the soliton always move with velocity less than that of longitudinal sound in the chain they do not emit phonons, i.e., their kinetic energy is not transformed into thermal energy. This is one important reason for the high stability of the new soliton. In addition the energy of the soliton state is below the bottom of the bare exciton bands, the energy gap being for small velocity of propagation. Hence there is an energy penalty associated with the destruction with transformation from the soliton state to a bare exciton state, i.e, the destruction of the soliton state requires simultaneous removal of the lattice distortion. We know in general that the transition probability to a lattice state without distortion is very small, in general, being negligible for a long chain. Considering this it is reasonable to assume that such a soliton is stable enough to propagate through the length of a typical protein structure. However, the thermal stability of the soliton state must be calculated quantitatively. The following calculation addresses this point explicitly [24- 26].We now diagonalize partially the model Hamiltonian in order to calculate the lifetime of the soliton, Eq. (43), using the quantum perturbation method [14] Since one is interested in investigating the case where there is initially a soliton moving with a velocity v on the chains, it is convenlent to do the analysis in a frame of reference where the soliton is at rest. We should then consider the Hamiltonian in this rest frame of the soliton, $\tilde{H}-\text{vP},$ where P is the total momentum, and $\text{P}={\displaystyle \sum _q\hslash q(a_q^{+}{a}_q-B_q^{+}{B}_q)}$ ,where $B_q^{+}=\frac{1}{\sqrt{N}}{\displaystyle \sum _n{e}^{iqn{r}_0}}{B}_n^{+}$ . Also, in order to have simple analytical expressions we make the usual continuum approximation. This gives $$\begin{array}{c}\tilde{H}={\displaystyle {\int }_0^{L}dx2\left[({\epsilon }_0-2J){\varphi }^{+}(x)\varphi (x)+J{r}_0^2\frac{\partial {\varphi }^{+}}{\partial x}\frac{\partial \varphi }{\partial x}-\frac{i\hslash v}{2}\left(\frac{\partial {\varphi }^{+}}{\partial x}\varphi (x)-{\varphi }^{+}(x)\frac{\partial \varphi }{\partial x}\right)\right]}\\ +{\displaystyle \sum _q\hslash ({\omega }_q-qv)a_q^{+}{a}_q+\frac{1}{\sqrt{N}}{\displaystyle \sum _q}2[g_1(q)+2g_2(q)](a_{-q}^{+}+a_q){\displaystyle {\int }_0^{L}dx{e}^{ikx}{\varphi }^{+}(x){\varphi }^{}(x)}}\text{ (47)}\end{array}$$ where $\varphi (x)$ represents now the field operator corresponding to Bn in the continuum limit (whereas before it only indicated a numerical value). Here L=Nr0, -π< kr0< π, and ${\omega }_q$ ≈(w/M)^1/2 r0•|q| ,x=nr₀. Since the soliton excitation is connected with the deformation of intermolecular spacing, it is necessary to pass in Eq.(47) to new phonons taking this deformation into account. Such a transformation can be realized by means of the following transformation of phonon operators [29] $${\text{b}}_{\text{q}}=a_p-\frac{1}{\sqrt{N}}{\alpha }_q,b_q^{+}=a_q^{+}-\frac{1}{\sqrt{N}}{\alpha }_q^{*},\text{ (48)}$$ which describe phonons relative to a chain with a particular deformation, where b_q (b_q⁺) is the annihilation (creation) operator of new phonon. The vacuum state for the new phonons is $$|\tilde{0}{\rangle }_{ph}=\exp \left[\frac{1}{\sqrt{N}}{\displaystyle \sum _q({\alpha }_q(t)a_q^{+}-{\alpha }_q^{*}(t))a_q}\right]|0>{}_{ph}\text{ (49)}$$ which is a coherent phonon state[30], i.e., $|\tilde{0}{\rangle }_{ph}=0$ . The Hamiltonian $\tilde{H}$ can now be rewritten[ 24-26] as $$\begin{array}{c}\tilde{H}={\displaystyle {\int }_0^{L}2dx\varphi (x)[{\epsilon }_0-2J+V(x)-J{r}_0^2\frac{{\partial }^2}{\partial x^2}+i\hslash \frac{\partial }{\partial x}]}\varphi (x)+\\ {\displaystyle \sum _q\hslash ({\omega }_q-qv)[b_q^{+}{b}_q+\frac{1}{\sqrt{N}}({\alpha }_q{b}_q^{+}+{\alpha }_q^{*}{b}_q^{+})]}+W\text{'}+\\ \frac{1}{\sqrt{N}}{\displaystyle \sum 2[g_1(q)+2g_2(q)](b_{-q}^{+}+b_q){\displaystyle {\int }_0^{L}dx{e}^{iqx}{\varphi }_{}^{+}(x)}}\varphi (x)(50)\end{array}$$ where $$W\text{'}=\frac{1}{N}{\displaystyle \sum _q\hslash }({\omega }_q-qv)|{\alpha }_q^{}{|}^2,V(x)=\frac{1}{N}{\displaystyle \sum _q[g_1(q)+2g_2(q)]({\alpha }_{-q}^{*}+{\alpha }_{-q}^{})e^{iqx}}\text{ (51)}$$ To describe the deformation corresponding to a soliton in the subspace where there is $${\displaystyle {\int }_0^{L}dx{\varphi }_{}^{+}(x){\varphi }_{}^{}(x)}=1$$ from Eq(45) in such a case. From the above formulae we can obtain $$\text{V}\left(\text{x}\right)=-2J{\mu }_p{}^2\sec h^2({\mu }_{p}x/r_0)\text{ (52)}$$ In order to partially diagonalize the Hamiltonian Eq.(50) we introduce the following canonical transformation[14,23] $$\varphi (x)={\displaystyle \sum _j{A}_j}{C}_j(x),{\varphi }^{+}(x)={\displaystyle \sum _j{C}_j^{*}(x)A_j^{+}}\text{ (53)}$$ where $${\displaystyle \int C_1^{*}(x)C_j(x)dx={\delta }_{lj},{\displaystyle \sum _j{C}_j^{*}(x^{\prime })C_j(x)=\delta (x-x^{\prime }),{\displaystyle \int dx|C_j(x){|}^2=1}}}\text{ (54)}$$ The operators A_s⁺ and A_k⁺ the creation operators for the bound states C_s(x) and delocalized state C_k(x), respectively. The detailed calculation of the partial diagonalization and of corresponding C_s(x) and C_k(x) are described in Appenix A. The partially diagonalized Hamiltonian obtained is as follows $$\begin{array}{c}\tilde{H}=W\text{'}+E_s{A}_s^{+}{A}_s+{\displaystyle \sum _k{E}_k{A}_k^{+}{A}_k+}{\displaystyle \sum _q\hslash ({\omega }_q-qv)b_q^{+}{b}_q}+\\ \frac{1}{\sqrt{N}}{\displaystyle \sum _q\hslash ({\omega }_q-qv)(b_q^{+}{\alpha }_q+{\alpha }_q^{\ast }{b}_q})(1-A_s^{+}{A}_s)+\frac{1}{\sqrt{N}}{\displaystyle \sum _{k{k}^{\prime }q}F(k,k^{\prime },q)(b_{-q}^{+}+b_q)}{A}_{k^{\prime }}^{+}{A}_k\\ -\frac{1}{\sqrt{N}}{\displaystyle \sum _{kq}\tilde{F}(k,q)(b_{-q}^{+}+b_q)}(A_s^{+}{A}_{-k}-A_k^{+}{A}_s)\text{ (55) }\end{array}$$ and $$C_s(x)={(\frac{{\mu }_p}{2r_0})}^{1/2}\sec h({\mu }_{p}x/r_0)\exp [i\hslash xv/2J{r}_0^2],$$ $$\text{with }{E}_s=2\left[{\epsilon }_0-2J-\frac{{\hslash }^2{v}^2}{2J{r}_0^2}-{\mu }_p^{2}J\right]\text{ (56a)}$$ $$C_k(x)=\frac{{\mu }_p\tanh ({\mu }_{p}x/r_0)-ik{r}_0}{\sqrt{N{r}_0}[{\mu }_p-ik{r}_0]}\exp [ikx+\frac{i\hslash vx}{2J{r}_0^2}],\text{ (56b)}$$ with $$E_k=2\left[{\epsilon }_0-2J-\frac{{\hslash }^2{v}^2}{2J{r}_0^2}-J{(k{r}_0)}^2\right]$$ where $$F(k,k^{\prime },q)=2[g_1(q)+2g_2(q)]{\displaystyle {\int }_0^{L}dx{e}^{iqx}{C}_{k^{\prime }}^{\ast }}(x)C_k(x)$$ $$\approx 2[g_1(q)+2g_2(q)]\left\{1-\frac{i{\mu }_{p}q{r}_0}{[{\mu }_p+i(k+q)r_0][{\mu }_p-ik{r}_0]}\right\}\approx F[k,(k+q),q]{\delta }_{k^{\prime }k+q}\text{ (57)}$$ $$\tilde{F}(k,q)=2[g_1(q)+2g_2(q)]{\displaystyle {\int }_0^{L}dx{e}^{iqx}{C}_{k^{\prime }}^{\ast }}(x)C_s(x)$$ $$=\frac{2\pi }{{\sqrt{2\mu }}_p}[g_1(q)+2g_2(q)]\left\{\frac{iq{r}_0}{[{\mu }_p+ik{r}_0]}\right\}\sec h[\pi (k-q)r_0/2{\mu }_p]\text{ (58)}$$ where ${\alpha }_q$ is determined by V(x) and the condition $({\omega }_q-\text{vq)}{\alpha }_q=({\omega }_q+\text{qv)}{\alpha }_q^{\ast }$ which is required to get the factor, $(1-A_s^{+}{A}_s)$ in the $\tilde{H}$ in Eq.(55). Thus we find $${\alpha }_q=\frac{i\pi ({\chi }_1+{\chi }_2)}{w{\mu }_p(1-s^2)}{\left[\frac{M}{2\hslash {\omega }_q}\right]}^{1/2}({\omega }_q+qv)\csc h(\pi q{r}_0/2{\mu }_p)\text{ (59)}$$ and $W\text{'}=\frac{2}{3}{\mu }_p^{2}J$ For this ${\alpha }_q$ the $|\tilde{0}{>}_{\text{ph}}$ in Eq.(49) is just the coherent phonon state introduced by Davydov. However, the bound state C_s(x) in Eq.(56a), unlike the unbounded state C_k(x) in Eq.(56b), is selfconsistent with the deformation. Such a self-consistent state of the intramolecular excitation and deformation forms a soliton which in the intrinsic reference frame is stationary. For the new soliton described by the state vector $|\psi >=\frac{1}{\sqrt{2!}}{(A_s^{+})}^2|0{>}_{ex}|\tilde{0}{>}_{ph}$ he average energy of $\tilde{H}$ in Eq.(55) is $$<\psi \left|\tilde{H}\right|\psi >=2({\epsilon }_0-2J-\frac{{\hslash }^2{v}^2}{4J{r}_0^2})-\frac{4}{3}J{\mu }_p^2\text{ (60)}$$ Evidently, the average energy of $\tilde{H}$ in the soliton state $|\psi >$ Eq.(60), is just equal to the above soliton energy E_sol, or the sum of the energy of the bound state in Eq.(56a), Es, and the deformation energy of the lattice,W’, i.e., $<\psi \left|\tilde{H}\right|\psi >=E_{sol}=E_s+W’$ This is an interesting result, which shows clearly that the quasi-coherent soliton formed by this mechanism is just a self-trapping state of the two excitons plus the corresponding deformation of the amino acid lattice. However, it should be noted that $|\psi >$ is not an exact eigenstate of $\tilde{H}$ owing to the presence of the terms in $\tilde{H}$ with $A_k^{+}{A}_s$ and $A_s^{+}{A}_{-k}$