A generalised Davydov-Scott model for polarons in linear peptide chains

We present a one-parameter family of mathematical models describing the dynamics of polarons in periodic structures, such as linear polypeptides, which, by tuning the model parameter, can be reduced to the Davydov or the Scott model. We describe the physical significance of this parameter and, in the continuum limit, we derive analytical solutions which represent stationary polarons. On a discrete lattice, we compute stationary polaron solutions numerically. We investigate polaron propagation induced by several external forcing mechanisms. We show that an electric field consisting of a constant and a periodic component can induce polaron motion with minimal energy loss. We also show that thermal fluctuations can facilitate the onset of polaron motion. Finally, we discuss the bio-physical implications of our results.


Introduction
The polaron, a quasi-particle formed by the coupling of an electron to a vibrating lattice, was first theorised by Landau in 1933 [1]. In essence, polaron formation is a process of electron self-trapping. The presence of the electron causes localised distortions in the natural vibrational mode of the lattice, a.k.a. the lattice phonon. In return, if the electromagnetic interaction between the electron and lattice is appropriate, then the phonon distortions can lower the potential well for the electron, thus trapping the electron.
Some twenty years after the inception of the polaron concept, a mathematical description of it was formalised by Fröhlich [2] and subsequently Holstein [3,4]. Since then, properties of the Fröhlich-Holstein polaron have been well studied, with some authors hypothesising an application of dynamical polarons as electron transporters in conductive material [5][6][7]. In 1970s, A.S. Davydov used the basis of polaron theory to explain some biological processes [8]. Specifically he proposed that, in an α-helical protein, a certain intramolecular oscillator can interact with the peptide chain in a way similar to an electron interacting with a crystal lattice. Davydov suggested that this interaction could lead to the localisation and propagation of vibrational energy in the α-helix. Later, A.C. Scott modified Davydov's theory, taking into account the internal geometry of peptide units [9]. Some authors argued that, given the polarisability of peptide units, electron self-trapping is also possible in proteins, and it can be described by the same mathematical model that Davydov and Scott used [10,11]. It should therefore be possible that polaronic a e-mail: b.m.a.g.piette@durham.ac.uk transport of electrons may take place in proteins, too. Recently, Brizhik et al. reported on the properties of static and dynamical polarons in simple molecular chains, and adverted to the applicability of their results to electron transport in biomolecules such as proteins [12][13][14]. Their studies were based on the Davydov-Scott model.
In the current study, we propose a generalisation to the Davydov-Scott model, and use it to explore the properties of polarons in a linear peptide chain. In the generalised model, there is an extra parameter which represents the extent to which the electron-polypeptide interaction is spatially symmetric. In Section 2, we describe our model and explain why the extra parameter is necessary. We also give physical interpretations of all other parameters in the model, justifying the choices of their values where possible. Then, we derive a set of coupled dynamical equations which govern the electron and phonon parts of the polaron, as well as how they interact. In Section 3 we look at solutions to our equations which are stationary, and thus deduce properties of static polarons admissible by our model, such as the polaron's binding energy. The process of solving the equations is carried out analytically as well as numerically. By the former approach, a closed-form expression for the solution is found, but its use is limited, because the solution process involves a few approximations and simplifying assumptions. By the numerical approach, no convenient expression for the solution is possible, but the method solves the equations directly without simplifications. We compare the results produced by the two different methods. Section 4 concerns dynamical polarons. We discover that it is possible to use a suitable external forcing to displace the stationary polaron, and to sustain its motion in such a way that its energy remains highly stable.
We investigate how the polaron's motion depends upon our forcing parameters. We use only numerical methods to obtain our results in Section 4, as well as those in Section 5, where we consider how the polaron's motion is affected by temperature of the environment. For this part, the external forcing from Section 4 remains in place, but we also utilise a parameter which controls the magnitude of the thermal effect. To account for the random nature of thermal fluctuations, we repeat each numerical simulation many times over, taking the average of the results. Finally, we conclude by discussing the physical realisabililty of our mathematical model, particularly how the external forcing which we study in Section 4 may be realised. We also briefly discuss the generalisability of our model to studying electron transport by polarons in α-helices.

The model and dynamical equations
In both Davydov's and Scott's models, the Hamiltonian for a system of excitons interacting with one-dimensional lattice phonons is written in Fröhlich-Holstein formĤ = H e +Ĥ p +Ĥ int , whereĤ e ,Ĥ p andĤ int represent energy contributions from the exciton, phonon and interaction parts, respectively [2,3,9,15,16]. We adopt this Hamiltonian for our model, and following [12][13][14] we consider an additional external Hamiltonian,Ĥ ext , so that our Hamiltonian takes the form whereĤ e describes a tight-binding electron, the stretching and compressing of hydrogen bonds in the peptide chain are phonon oscillations described byĤ p ,Ĥ int accounts for the electron-phonon interaction, andĤ ext represents the effect of an external electric field. We assume that the peptide chain consists of N + 1 identical units and N identical hydrogen bonds. In the tight-binding approximation, we havê The subscript n in equation (2) labels peptide units, which are the unit cells of our lattice.Â † n andÂ n are local electron creation and annihilation operators, respectively. J 0 is the potential energy of a localised electron. Modelling each unit as a point-dipole, we assume the nearestneighbour dipole interaction energy is a constant and write it as −J 1 [17][18][19]. The external Hamiltonian, models the effect of an electric field with strength E(t) on the potential energy of a localised electron with charge −q. The potential energy due to E(t) is set to zero at some arbitrary n 0 , and R = 4.5Å is the equilibrium lattice spacing. Since the electron mass is several orders smaller than the mass of a peptide unit, we take a semi-classical approach where the phonon Hamiltonian,Ĥ p , is a classical one. In the harmonic approximation, the hydrogen bonds are modelled as Hookean springs with force constant K, and thereforeĤ p takes the form where M = 1.774 × 10 −25 kg is the average mass of a peptide unit in a membrane α-helix [20], and we have defined Ω := K/M . U n and P n are, respectively, the displacement and conjugate momentum of the nth unit. Thus, the first and second sums in the expression for H p represent, respectively, the kinetic and potential energies of the lattice. We take the value of Ω to be the natural angular frequency of slow phonons in an α-helix, Ω = 5.5 × 10 12 s −1 [21][22][23]. To derive the interaction Hamiltonian,Ĥ int , Davydov and Scott assumed that the energy of an on-site excitation depends on lattice deformations in its vicinity. For us, the local deformations are namely the amount by which the lengths of hydrogen bonds deviate from equilibrium. If we write the electron energy at site n in a Taylor expansion, the first two terms are J 0 + χG n (S n , S n−1 ), where χ is a constant, G n is a bilinear function, and |χG n /J 0 | 1. Then the interaction Hamiltonian isĤ int = N n=0 χG nÂ † nÂ n . Davydov assumed that S n and S n−1 have equal influence on local excitation energies [15], so G n = (S n + S n−1 ) /2, and H int iŝ Davydov's model is therefore spatially symmetric, sincê A † nÂn is coupled equally to U n+1 and to U n−1 . Scott modified Davydov's model by opting for the antisymmetric G n (S n , S n−1 ) = S n instead [9]. The reason is, while both authors let an intra-peptide C=O oscillator take the role of the exciton, Davydov neglected the internal geometry of the peptide units, but Scott pointed out that every unit has its C=O pair immediately adjacent to the next hydrogen bond in the chain. This leads Scott to assume that A † nÂn is coupled to, without loss of generality, S n , and not S n−1 . Scott therefore had Since we are modelling an electron as opposed to an intrapeptide oscillator, we cannot use Scott's argument to justify assuming thatÂ † nÂ n is coupled to U n+1 and not to U n−1 . Nor should we assume that on-site energies are affected equally by deformations on both sides, as Davydov did. We therefore propose G(S n , S n−1 ) = χ r S n + χ l S n−1 , taking without loss of generality χ r > 0 and 0 ≤ χ l ≤ χ r . Then, By defining so that 0 ≤ β ≤ 1, we can writê We treat χ as an adjustable parameter. Setting β = 0 (χ l = χ r ) gives us the symmetric model of Davydov as per equation (6), whilst setting β = 1 (χ l = 0) produces the antisymmetric model of Scott as per equation (7). The larger β is, the less spatial symmetry our model possesses. Indeed, for β ∈ [0, 1), the ratio of n-(n + 1) coupling strength to n-(n − 1) coupling strength is given by χ r /χ l = (1+β)/(1−β), and this ratio is strictly increasing with β. We write the electronic state of the system as a linear superposition of local excitations [12], where |vac is the vacuum state, and α n ∈ C is the probability amplitude for an electron localised at the nth site, subject to the normalisation condition, N n=0 |α n | 2 = 1.
We proceed to derive dynamical equations for α n and U n . By equating coefficients ofÂ † n |vac on both sides of the Schrödinger equation, i d |Ψ /dt = (Ĥ e +Ĥ int +Ĥ ext ) |Ψ , we obtain We have defined so that equation (13) holds for all n including the boundary terms (n = 0, N). Equations for U n are derived from classical Hamilton equations, dU n /dt = ∂H cla /∂P n and dP n /dt = −∂H cla /∂U n , where H cla := Ψ |(Ĥ p +Ĥ int )|Ψ . These equations are In order that equation (15) holds at the boundaries, we have set This boundary condition is justified because we expect the probability distribution |α n | 2 to be highly localised with half-width of O(1), and because we will be working with large lattices with N 1. Next, we introduce the gauge transformation, which sets J 0 = 2J 1 in equation (13). Physically this represents a shift in the arbitrary reference value from which energy is measured. We then have the discrete Laplacian, −J 1 (α n+1 + α n−1 − 2α n ), on the r.h.s. of equation (13). Meanwhile, to account for the interaction between the peptide chain and its environment, we need to add Langevin terms to the r.h.s. of equation (15) [13,14,24,25]. They are, a damping term describing energy dissipation due to friction, −Γ dU n /dt, where Γ is the viscous damping coefficient; and a stochastic term F n (t), describing random forces due to thermal fluctuations. Specifically, F n (t) is normally-distributed with zero mean and correlation function where k B is the Boltzmann constant and Θ is the temperature of the environment. Scaling time by Ω −1 and length by R gives us the following non-dimensionalised dynamical equations for ψ n and u n := U n /R, for n = 0, 1, . . . , N. where we have defined (20) and where the overdot denotes differentiation with respect to dimensionless time τ , and Equation (19) holds subject to the boundary conditions (14) and (16), as well as the normalisation condition, It is easily verifiable that setting β = 0 and β = 1 in equation (19) produces Davydov's and Scott's dynamical equations, respectively [9,15]. ρ is known as an adiabaticity parameter, as it is the ratio of the characteristic time scale of phonon vibrations to that of electronic phase variations [26]. We fix J 1 , following [12][13][14] (which used a different scaling), at ρ = 2.1. Moreover, since M, R and Ω are fixed, the ratio σ/δ = 1880 is constant. The range of δ which we consider throughout this study correspond to χ ∼ O (10 −11 ) N, agreeing with [9]. Finally, we take γ = 0.05, agreeing with [12][13][14] up to different scaling factors.

Stationary polaron solutions
We derive stationary polaron solutions to equation (19), subject to zero electric field ( = 0) and zero temperature (f n = 0). We consider analytical and numerical methods separately, and compare the results.

Analytical results
When f n = 0 andu n =ü n = 0, equation (19b) becomes which holds if Putting equation (24) into equation (19a) and requiring = 0 gives us Defining and the composite parameters we can rewrite equation (25) as We note that, since M, Ω and J 1 are all fixed, the parameter λ inherits the adjustability of χ. Now, in a stationary state, the time dependence of ψ n can be at most a variation of the phase factor. Following [15], we consider the ansatz where ξ is a real, continuous variable with −N R/2 ≤ ξ ≤ N R/2, φ is a real, twice-differentiable function, and H 0 and k are constants. In particular, H 0 is an energy eigenvalue, in the sense that In the limit N 1, R becomes small compared to the domain size, which enables us to invoke the continuum approximation, implying Putting equations (30)-(33) into equation (29), then dividing by exp (iρH 0 τ + ikξ) and retaining terms up to O(R 2 ), we obtain The last term in equation (34) is equivalent to Equating imaginary parts of equation (34) gives us k = 0. After the scaling x := ξ/R, the real part of equation (34) becomes The subscript x denotes differentiation with respect to x. We seek φ(x) which satisfies equation (36) for all x, not just when x = n − N/2. Then, from such a φ(x) we will be able to recover the discrete solution ψ n (τ ) via ξ = xR and equation (30). Further to being globally defined, we require that φ(x) has vanishing derivatives at infinity, and satisfies the normalisation condition, If η = 0 (i.e. β = 1), then equation (36) reduces to the nonlinear Schrödinger equation with a cubic nonlinearity, which has a well-known solution satisfying all the above constraints [9], we rewrite equation (36) as We dedicate the remainder of this section to analysing equation (39). It is an autonomous equation for φ(x), which allows us to define h(φ) := φ x , and write We define y(φ) := h 2 , so that y φ = 2hh φ , and multiply equation (39) by 2 to obtain The l.h.s. of equation (41) is the total derivative of (1 + 2ηφ 2 )y with respect to φ. We therefore have The integration constant C is determined by considering the limit x → ∞, in which φ 2 → 0 and y ≡ (φ x ) 2 → 0. We therefore have C = 0. Now we note that, if H 0 ≤ 0, then the r.h.s. of equation ( We then define Φ := φ 2 , and equation (43) becomes If equation (44) has a solution which is globally nonnegative and twice-differentiable, has vanishing derivatives at infinity, and satisfies then we claim that Φ(x) must attain its global upper bound of 2H 0 /λ at some finite x, and that every local maximum of Φ must also be a global maximum. The proof of this claim is as follows. Since Φ cannot be identically zero, and since lim x→±∞ Φ(x) = 0, Φ(x) must have at least one turning point, at some finite x and non-zero Φ. But we observe from equation (44) that Φ x vanishes if and only if Φ = 0 or 2H 0 /λ. Therefore, Φ(x) must attain its global upper bound of 2H 0 /λ at least once, and no other local maximum value is possible. This concludes the proof. We further propose that wherever Φ(x) attains its maximum value, say at x = x max , the second derivative Φ xx does not vanish there. The proof is as follows. On the one hand, we have On the other hand, equation (36) is equivalent not only to equation (39) but also to , as required. A corollary of the above proposition is that there must exist some neighbourhood of x = x max containing no maxima of Φ(x) other than x max itself. Without loss of generality, let x max = 0. Suppose the corollary is false, so that every neighbourhood of x = 0 contains some non-zero x at which Φ(x) is maximal. Then, there must exist some sequence x n approaching 0 such that Φ(x) is maximal at every x n , with Φ(x n ) = Φ(0). But this leads to a contradiction. Indeed, for every x n we have the Taylor expansion where the first derivative is absent because Φ(x) has a maximum at 0. It then follows that Φ xx (0) = lim n→∞ 2(Φ(x n ) − Φ(0))/x 2 n = 0, which contradicts the previous proposition. Therefore the corollary is proven. Since 0 has a maxima-free neighbourhood, we say that Φ(x) has an isolated maximum at 0.
We note that we can indeed require that x = 0 is a maximum of Φ(x), because equation (44) is translationally invariant: if Φ(x) is a solution then so is Φ(x − c) for any constant c. We exploit this invariance, requiring that Φ(x) satisfies Now we claim that there exists b > 0, which may be infinite, such that lim x→b Φ(x) = 0, and Φ(x) = 0 for all x ∈ (0, b). The proof of this claim is as follows. If Φ(x) = 0 for all x ∈ (0, ∞), then we are done. If Φ(x) = 0 for some x ∈ (0, ∞), then the set {x ∈ (0, ∞) : Φ(x) = 0} must have a minimum. If it does not, then there would be a sequence x n > 0 such that, as n → ∞, x n → 0 and Φ(x n ) → 0; but this would contradict the continuity of Φ(x) at x = 0. Thus, letting b equal the least positive zero of Φ(x), then we are done. Next, we propose that no other solution on [0, b) exists, and the proof is as follows. If Φ(x) has any maxima in (0, ∞), then the set must have a minimum, because otherwise we would have a contradiction to the fact that x = 0 is an isolated maximum of Φ(x). Let x 1 = min M, and suppose Since Φ x vanishes only at maxima and minima, it follows that Φ x is non-vanishing on (0, b), and therefore Φ(x) is strictly decreasing on [0, b). That is, any solution to equation (44) on [0, b) satisfying all the aforementioned constraints must also satisfy where On is continuous and non-zero, so the reciprocal function 1/g(Φ) is continuous and bounded, and therefore Riemann integrable. But g(Φ) approaches 0 as Φ → Φ 0 , meaning 1/g(Φ) becomes unbounded.
, and the uniqueness of Φ(x) follows.
In summary, we have so far established the following. If equation (44) has a solution which is globally non-negative and twice-differentiable, has vanishing derivatives at infinity, and has the property that its integral over R is 1, then equation (44) has a solution, say Φ(x), which satisfies all the above constraints as well as the condition (46), and there exists some b > 0 which may be infinite such that Φ(x) is strictly decreasing on [0, b), and lim x→b − Φ(x) = 0, and Φ(x) is the unique solution on [0, b). Moreover, using exactly the same arguments as above, it can also be shown that there exists some a < 0 which may be infinite such that Φ(x) is strictly increasing on (a, 0], and lim x→a + Φ(x) = 0, and Φ(x) is the unique solution on (a, 0]. On the interval of uniqueness, (a, b), Φ x is given by where g is defined by equation (48), and sgn is the sign function. Now we describe a method which, given λ and η, determines the unique Φ(x) on (a, b), and also determines a, b, H 0 in the process. Indeed we will show that for any λ and η, the interval of uniqueness for Φ(x) must be (a, b) = R. The fact that (a, b) = R shall have the following subtle consequence. Note that the derivation of equation (44) involved a multiplication by Φ ≡ φ 2 . Thus, the deduction from equation (44) back to equation (39) holds on the condition that Φ = 0. Since a and b are the smallest (in magnitude) zeros of Φ(x), we see that the equivalence between equations (44) and (39) breaks down outside the interval (a, b). That is to say, equations (44) and (39)  The method is as follows. For x ∈ [0, b), consider the coordinate transformation, where is a bijection from (0, Φ 0 ] to (0, 1], and the inverse sech function, arsech, is a bijection from (0, 1] to [0, ∞). Therefore, all the coordinate transformations are invertible. For x ∈ (a, 0], we consider exactly the same transformations as equation (52). Differentiating Z with respect to x we find, for all x ∈ (a, b), where we have used definition (48) of g(Φ). Moreover, by definition we have Y = sechZ, and it follows that 2ηΦ = 2ηΦ 0 Y 2 = 2η(2H 0 /λ)sech 2 Z. Defining we rewrite equation (54) as Due to equation (52), we have Z(x = 0) = 0. We can therefore solve equation (56) as follows. implying Now we can determine the values of a and b. In the limit Z → +∞, the definition of the coordinate transformations, as per equation (52), dictates that we must have either x → a or x → b. At the same time, equation (58) dictates that we must have x → ±∞, because the arctan function on the r.h.s. of equation (58) is bounded, whilst the arsinh function diverges to +∞. It therefore follows that (a, b) = R. The next step is to rewrite equation (58) as an expression for x in terms of Φ, so that we can invert the expression to find Φ(x) for x ∈ R. By definition (52) we have cosh 2 Z = 1/Y 2 = Φ 0 /Φ, and it follows that Since Z is by definition non-negative, we must take the positive square root, We claim that, given Φ 0 > 0 and x ∈ R, equation (59) uniquely determines a value of Φ > 0. The proof is as follows. If x = 0, then immediately from equation (59) we have Φ = Φ 0 , and we are done. If x = 0, consider the function where x and Φ 0 are parameters. Differentiating equation (60) with respect to Φ, we find Thus, G is strictly decreasing for Φ > 0. Since G(Φ) → ∞ in the limit Φ → 0, and G(Φ 0 ) = −sgn(x) √ H 0 x < 0, and G is continuous, we must have G vanishing at exactly one value of Φ ∈ (0, Φ 0 ). This concludes the proof. Moreover, we observe that in equation (59) In practice, given any Φ 0 > 0 and x ∈ R, we can compute Φ(x) by locating the zero of G(Φ). However, the value of Φ 0 cannot be freely chosen. Instead, it is determined by the normalisation condition (45) which, since Φ(x) is an even function on R, now reads where the Z + x is the positive-x branch of Z x , as per equation (56). It then follows that Multiplying equation (64b) by 2 √ ν, replacing H 0 by λΦ 0 /2, and replacing ν by 4ηH 0 /λ = 2ηΦ 0 , we obtain the following transcendental equation for Φ 0 .
To show that exactly one solution to equation (65) exists, we consider the function Differentiating F (Φ 0 ) with respect to Φ 0 , we find and F is continuous, we must have F vanishing at exactly one value of Φ 0 > 0. In practice, given parameters λ and η, we can compute Φ 0 by locating the zero of F (Φ 0 ), and Φ 0 uniquely determines the energy eigenvalue, H 0 = λΦ 0 /2. We can then feed the value Φ 0 into equation (59), and then for every x ∈ R we can find Φ(x) by means we have already described. In summary, given λ and η, equations (46), (59) and (65) together constitute an analytical solution to equation (44); and as we have already proven, it must be the unique global solution to equation (44) which satisfies all the constraints we have imposed.
The eigenvalue H 0 provides a link between Φ(x) and the binding energy of the stationary polaron. By equations (11), (17) and (30), where k = 0 and ξ = xR, the polaron state is written in terms of local excitations as |Ψ = N n=0 α nÂ † n |vac , and in the limit N 1, we have So the stationary By definition, the polaron's binding energy, E b , is its total internal energy measured with respect to J 0 . In units of J 1 , we have An expression for Ĥ p in terms of Φ(x) can be found by using equations (4) and (24). Since the polaron is stationary, the kinetic part of Ĥ p is zero, so we have It then follows that We have made use of definition (28) of λ, η, the fact that |ψ n | = |φ n | for all n, as well as the fact that |φ n | 2 is approximated by Φ(n − N/2). Figure 1 shows how various aspects of the stationary polaron depends upon the symmetry parameter β, and the effective coupling parameter λ. Recall that the former is a measure of the spatial symmetry of the electron-phonon interaction, and the latter measures the strength of this interaction. These are the only two parameters that affect the stationary polaron's physical properties (as η is merely a convenient combination of β and λ). Figure 1 shows how Φ 0 and the half-width of the polaron varies with β and λ. We define the half-width as the distance between the two x-values at which Φ(x) = Φ 0 /2, and it is a measure of how localised the polaron is. As one would expect, the half-width is negatively correlated with Φ 0 , which is the maximum height of Φ(x). The figure shows Φ 0 increasing with λ, and half-width decreasing with λ, and the rate of change of each quantity is greater given larger values of β. That is to say, the more spatially asymmetric the electron-lattice interaction is, the more influential λ is. The figure also has the following implication on the accuracy of Φ(x) as an approximation to the discrete stationary solution to equation (19). In a discrete solution, ψ n = exp(iρH 0 τ )φ n , the physical interpretation of |ψ n | 2 is the probability of the electron being localised at the nth lattice site. Therefore, the normalisation condition is defined in terms of a sum, N n=0 |ψ n | 2 = 1, and consequently we must have |ψ n | 2 ≤ 1 for all n. When a continuum solution Φ(x) is used to approximate the discrete one, we have the relation Φ 0 ≡ max |ψ n | 2 . Thus, any continuum solution with Φ 0 > 1 cannot be reliable as an approximant. When β = 1, Φ 0 exceeds 1 if λ is greater than 8, since Φ 0 = λ/8. On the other hand, when β = 0, we computed Φ 0 for λ up to 100, and Φ 0 remains less than 0.6.
In Figure 1 we see that H 0 increases with λ whilst the polaron's binding energy gains magnitude, meaning the larger λ is the more energy is required to break up the polaron. Once again, the larger β is, the more rapidly these quantities vary with λ. We note that the thick (black) curve for H 0 , corresponding to β = 1, is exactly the graph of H 0 = λ 2 /16, as per equation (38). In Figure 1, we see that a polaron which is more strongly bound has a larger Φ 0 and a smaller half-width, i.e., it is more localised.
The gradient of curves in Figure 1 vary with β, and the variation is more pronounced when β is close to 1. This suggests that the system is highly sensitive to variations in β when β is large, but not so when β is small. Moreover, as λ decreases, curves corresponding to different values of β begin to converge; specifically this happens when λ ≈ 1. This suggests that when λ is small, the extent to which the electron-phonon interaction is spatially symmetric has little bearing on the physical properties of stationary polarons.

Numerical solutions
In this section we solve equations (29) and (31) directly, using a numerical scheme, but not without the help of analytical results from Section 3.1. We then compare the resulting stationary polaron states with the ones we obtained via continuum approximation. Expanding equation (29) using the definitions of Δψ n and Δ |ψ n | 2 , we have Any solution ψ n to equation (74) is an attractor of the following map [27]. where We take the approximate solution from Section 3.1 as initial guess, and repeatedly apply equation (75) until convergence. When converged, ψ n is the stationary solution to equation (74), and H(ψ n ) is equal to H 0 . In practice, on a grid with N = 200, convergence is typically reached within O(10 5 ) iterations, which amounts to O(10 1 ) seconds of computing time. We have computed stationary solutions for various β and λ, and some results are presented in Figure 2. Figure 2a contains information about two key aspects of the stationary polaron state: the electron probability distribution, and the polaron binding energy. Qualitatively speaking, it is in agreement with predictions of the continuum approximation, as per Figure 1: as λ increases, the polaron becomes more localised, and more strongly bound. Moreover, the effect of increasing λ is more profound given larger values of β. However, further comparison between Figures 1 and 2a reveals a noteworthy difference. When β = 1, Φ 0 is a linear function of λ, whereas Figure 2a suggests that max |ψ n | 2 , which is approximated by Φ 0 , is not linearly dependent on λ. In fact, given any β, max |ψ n | 2 grows significantly faster with λ than Figure 1 predicts. Despite that, the growth of max |ψ n | 2 in Figure 2a eventually stalls, when λ becomes sufficiently large. This is a manifestation of a fundamental difference between the continuum and discrete equations, which we discussed in Section 3.1: the continuum equations place no limit on how large Φ 0 can be, whereas the discrete system limits max |ψ n | 2 to 1. Figure 2b shows a selection of |ψ n | 2 solutions. Comparing all the dotted (red) lines, which correspond to λ = 1.0 at various values of β, we see that they are essentially identical. This confirms the belief that when λ is close to 1, systems with different β-values unify. The figure also shows some stationary solutions to the other half of equation (19), namely u n , which is expressed in terms of the stationary |ψ n | 2 solution as per equation (24). Recall that physically u n represents the displacement of the nth molecule from its equilibrium position. In order that the point-dipole model for lattice units is valid, the lattice distortion must satisfy the condition |u n+1 − u n | 1 [19]. This condition is indeed fulfilled in the stationary polaron state, since according to Figure 2b we have max |u n+1 − u n | ∼ O(10 −2 ).
Comparing all the dashed (blue) lines in Figure 2b, which correspond to λ = 4.4 at various values of β, enables us to make the following observation. When β = 0, the u n solution is centred at the location of max |ψ n | 2 , in the sense that its graph is rotationally symmetric about n = 80. This agrees with our intuition that when β = 0, i.e. when the electron-phonon interaction is spatially symmetric, the electron in the stationary state causes equal lattice distortion to its left and right. As β increases, the maximum magnitude of u n remains the same, but the centre of u n shifts away from the location of max |ψ n | 2 , in response to the decrease in spatial symmetry. When β = 1, the molecule at the location of max |ψ n | 2 (n = 120 in this case) is barely displaced, whereas molecules to the right of this point are displaced considerably. Now, the potential energy in the lattice is a sum over terms of the form (u n+1 − u n ) 2 , which is the square of the gradient of the u n graph at site n. In the steady state, this gradient is zero except at a few sites around the location of max |ψ n | 2 , and it is clear that solutions corresponding to larger values of β have steeper gradients there. We therefore conclude that, in the stationary state, systems with greater spatial asymmetry store more potential energy in the lattice.
In Figure 2b we also see a comparison between some |ψ n | 2 solutions and their counterpart continuum approximations, Φ(x). In particular, we look at the thick solid (black) lines and their accompanying thin solid (black) lines. The comparison reveals that, fixing λ, in this case λ = 2.6, Φ(x) is a more accurate approximant to |ψ n | 2 when β is smaller. As β approaches 1, it becomes apparent that Φ(x) under-estimates the height of the |ψ n | 2 profile. If λ is sufficiently large, however, Φ(x) becomes an over-estimate of the profile height. This all comes down once again to the fact that the continuum equations do not limit the height of the Φ(x) solution.

Dynamical polarons at zero temperature
In this part of the study we explore properties of polarons which propagate along the peptide chain, under an external forcing (τ ), and zero temperature (f n (τ ) = 0). Physically, (τ ) may represent the strength of a time-dependent electric field. We solve equation (19) as an initial value problem, using the stationary ψ n and u n solutions which we computed in Section 3.2 as the initial configuration of the system. We prescribe a suitable (τ ), setting n 0 to the location where the stationary |ψ n | 2 attains its maximum. Then we integrate the system forward in time using the 4th-order Runge-Kutta method. To ensure numerical stability, the integration time-step is set at Δτ = 0.01.
As we forward-integrate the system, we keep track of several scalar quantities associated with the polaron, such as its half-width and binding energy. Most importantly, we keep track of the polaron's position, defined as follows. If |ψ n | 2 attains its maximum at lattice siten, then polaron position is the vertex location of the parabola extrapolated from three points: (n, |ψn| 2 ), (n − 1, |ψn −1 | 2 ), (n+1, |ψn +1 | 2 ). We note that, if the polaron is dynamical, then the stationary solution given by equation (70) is no longer valid, and therefore we cannot take equation (73) as the expression for the binding energy. Instead, the binding energy E b as per equation (71) will be computed directly from the numerical solutions.

Constant or periodic electric fields
The most obvious choice of (τ ) is a constant, Using moderately-localised stationary states (max |ψ n | 2 = 0.64) as initial conditions, we computed polaron trajectories under various values of¯ . Our results show that, given β = 1 and λ = 3.0, a constant forcing of anȳ ∼ O(10 −2 ) induces nothing but small oscillations of the polaron around its initial position. An example of trajectory is presented in Figure 3.
As¯ is increased beyond 0.1, we find that eventually the forcing does become strong enough to dislodge the electron from its potential well, and propel the polaron along the peptide chain. However, as the polaron propagates, the magnitude of its binding energy decreases rapidly, and the polaron "delocalises", i.e. breaks up into unbound components, within several hundred time units. A direct manifestation of the polaron's energy loss and eventual delocalisation is that the |ψ n | 2 profile loses height and gains local peaks at lattice sites far away from the global maximum. For the sake of consistency, throughout the remainder of this study we shall say that a polaron has delocalised if its maximum height drops to below 0.1, as it must then be the case that other local peaks have magnitudes comparable to the global maximum. Figure 4 shows an example of a constant forcing large enough to cause polaron displacement, and it illustrates the resultant rapid delocalisation of the polaron. Figure 4 shows the trajectory of a polaron which, within roughly 600 time units, is displaced by just over 300 lattice sites. Its binding energy steadily decreases in magnitude, until the polaron delocalises at τ ≈ 600. Meanwhile, Figure 4 shows the electron probability distribution, |ψ n | 2 , at the time of delocalisation. This |ψ n | 2 profile has evolved from an initial configuration possessing a maximum height of 0.64, and no local peaks apart from the global maximum.
If the polaron's binding energy decreases in magnitude, then the polaron's ability to transport energy is diminished. Beyond the example shown in Figure 4, all our results are consistent with the hypothesis that, regardless of β and λ, a constant causes the polaron to undergo either small periodic oscillations, or rapid losses in energy. We would like to find ways to displace the polaron without significant energy loss. Therefore, we must look for forms of other than constants. The next most natural choice of is periodic, where A is the amplitude and T is the period. Physically this may represent an electromagnetic plane wave which is monochromatic, i.e. coherent. Under periodic (τ ) with A up to 0.2, regardless of β and λ we find that the polaron simply oscillates about its initial position. The polaron's oscillatory motion has a period which coincides with T , and an amplitude which is positively correlated with A. An example of such trajectories is shown in Figure 3. While the polaron remains highly stable over time, its position averaged over its periodic remains constant. Thus, if we want polarons which transport energy from one lattice site to another, we must again look for an alternative form of (τ ).

Periodic electric fields with non-zero mean
Having studied the effects of constant forcing and periodic forcing in Section 4.1, and discovered that neither serves to displace the polaron with minimal energy loss, in this section we consider forcing of the form Equation (79) represents the combination of the two types of forcing considered previously, with a constant component and a sinusoidal one. One may also think of (t) as a mean-shifted periodic forcing (MSPF). In particular, the mean¯ is chosen to be lower than the constant forcing which is required to displace the polaron, in the manner of Figure 4. Therefore,¯ on its own would not give the electron enough energy to escape its potential well. But we hope that the component A can periodically push the electron energy over the threshold, resulting in polaron motion. Another possible advantage of this setup is that A may periodically lower the electron energy, slowing it down and giving the lattice time to "catch up", thus making the polaron motion more sustainable than it would be under a constant forcing. Mathematically, (τ ) depends on three independent parameters,¯ , A and T . Before investigating the effect of each of these parameters, we present Figure 5, which is a direct comparison with Figure 4.
We have replaced the constant forcing = 0.15, which resulted in Figure 4, with an MSPF which has the same maximum amplitude as before. The difference is that now this maximum amplitude is reached once every period T . Figure 5 shows that, within roughly 10 periods, the polaron is displaced by nearly 400 lattice sites. Contrary to the uniform manner in which the polaron moves in Figure 4, now the polaron moves towards one end of the peptide chain and then the other, within each period of (τ ). The overall displacement of the polaron is due to the fact that each movement to one end of the chain is larger than the subsequent swing back the other way. We note that while the polaron moves slightly further compared to where the numerator is the displacement of the polaron, which we denote by D. Our results show that, of the pa-rameters¯ , A and T , the dominant factor which deter-mines the polaron's velocity is the constant component¯ . We will discuss this in more depth in relation to Figure 6. Figure 5 also shows how the polaron's binding energy, E b , varies in time. Following an initial drop in magnitude, E b mostly oscillates between −0.75 and −1.5, until another sharp decrease in magnitude leading up to delocalisation at τ ≈ 4900. There is an important observation to be made here. In Figure 4, we see that when the polaron reaches lattice site n = 300, E b is about −0.6. In Figure 5, the polaron's average position over the 6th period is roughly 300, and the average binding energy over this period is −1.2. That is to say, under the MSPF, the polaron is carrying twice as much energy when it reaches n = 300, compared to when it reaches n = 300 under the constant forcing. Even though the constant forcing gets the polaron to n = 300 in less time, we consider the MSPF a better mechanism for polaron transport, because the polaron binding energy is more stable. Indeed, the same can be said when the destination n is anything larger than 30. If the destination is n < 30, then the constant forcing takes the polaron to n in such a small amount of time that it causes no more variation in binding energy than the MSPF does. In general, all our results are consistent with the hypothesis that, by splitting a constant forcing into constant and sinusoidal components, we lower the polaron's velocity but increase its stability and lifetime. We say that, compared to the constant forcing, the MSPF is a better long-distance transport mechanism, where speed can be sacrificed for energy efficiency.
In Figure 5 we see another aspect of the polaron's motion, namely, how the height and half-width of the |ψ n | 2 profile vary with time. Following an initial decrease, the profile height, max |ψ n | 2 , mostly oscillates between 0.2 and 0.4, until a sudden drop to 0.1, leading to delocalisation. Meanwhile, the half-width mostly oscillates between 1.5 and 4, following an initial growth. The peaks in the half-width, as well as the troughs in max |ψ n | 2 , occur precisely when the polaron turns from moving in one direction to moving in the other. This suggests that when the polaron accelerates, it "spreads out", and so the half-width widens and max |ψ n | 2 drops. We observe this phenomenon in all our results. Based on our observations, we theorise that a polaron's directed motion may be explained physically as follows.
Since the forcing (τ ) is the effect of an electric field, it first-and-foremost provides the electron with extra energy. This is evident in the dramatic energy variation during the first period of (τ ) (see Fig. 5). Following this, it becomes much easier for the electron to overcome the significantly diminished polaron binding energy, E b . This is why the onset of polaron motion always follows a drastic drop in magnitude of E b . Whenever | (τ )| becomes large enough to give the electron sufficient energy to overcome E b , the electron is dislodged from its potential well and propelled along the lattice. If the electron-lattice coupling is strong enough, then the lattice distortion can keep up with the electron, and so the polaron can remain intact. Whenever | (τ )| drops below the binding threshold, the electron-lattice interaction slows down the electron and causes its probability distribution to spread out. This is why the half-width of |ψ n | 2 always peaks at times when the polaron's instantaneous velocity is zero. If | (τ )| remains below the binding threshold for long enough, then the polaron's position can plateau, as is seen in Figure 6, particularly in the lowermost solid (black) lines in Figure 6. If (τ ) has a large enough periodic component A, then it is possible for | (τ )| to overcome the threshold twice per period: once with > 0, once with < 0. If the electron moves towards large n in the > 0 case, then it will move towards small n in the < 0 case. This explains the backwards swing exhibited by some polaron trajectories during each period of motion. The fact that¯ = 0 ensures that the electron always spends more time moving one way than the other, hence the overall directedness of the polaron trajectories. This point is most clearly demonstrated by the trajectories in Figure 6 where the period T = 2000.
Within each subfigure of Figure 6, we can compare trajectories with the same A (represented by line type) but different β (represented by starting position). This reveals the effect of varying the spatial symmetry of electron-lattice interaction. A greater spatial asymmetry, i.e. a larger β, causes the polaron to be more susceptible to displacement. Also within each of subfigure in Figure 6, we can compare trajectories with the same β but different A. The larger A is, the more the polaron oscillates back and forth during each period of motion. We can also compare a trajectory in Figure 6a with a certain A and β to that in Figure 6b with the same A and β. This suggests that the overall velocity of the polaron is determined by¯ , in the sense that the larger¯ is, the more the polaron moves per unit time. Finally, we can compare a trajectory in Figure 6b with a certain A and β to that in Figure 6c with the same A and β, hoping to see the effect of varying T . However, this comparison is not particularly enlightening. We therefore present Figure 7, which not only provides more insight into the effect of T , but also helps to quantify our observations, and reinforce our theories, about the effects of¯ and A.
For several combinations of¯ and T , we examine how the polaron's lifetime, displacement and velocity vary with A. The results are displayed together, in Figure 7, so that the effects of¯ and T can be gauged also. Firstly we consider the lifetime, τ d . We computed all our numerical solutions up to τ = 50 000, which is several times larger than the typical lifetime of a polaron that moves under the MSPF. If the polaron is not displaced by the MSPF, then it is effectively permanent, in the sense that its energy oscillates instead of dissipating over time, and it would have a lifetime far exceeding 50 000. Thus, in Figure 7  Next, we look at the polaron's displacement, D. When A is small, the polaron does not move barring small oscillations, the types of which we saw in Figure 3. As A reaches critical value A c , the polaron turns from being quasi-stationary to moving by several hundred lattice sites during its lifetime. Evidently, the value of A c is independent of T . Note that we only consider the displacement of polarons whose lifetimes are at least 2T , and we set the displacement of polarons with shorter lifetimes to zerosee for instance the dashed (blue) lines in the centre and bottom-middle subfigures.
Whilst the value of A c does not depend on T , the amount of displacement caused by A c does. However, it is unclear from our results what their correlation is. As A increases beyond A c , the qualitative behaviour of D is that it decreases. This is due to the fact that increasing A causes the polaron to delocalise more quickly, and therefore the polaron has less time to move. To understand how A affects the amount of polaron displacement per unit time, we examine the polaron's (average) velocity, V , as per definition (80). When A is small, V is zero. As A reaches critical value A c , the velocity becomes typically O(10 −2 ). Exactly what value this critical velocity V c takes depends on¯the larger¯ is, the larger V c is. As A increases beyond A c , sometimes V simply decays away -see for instance the topright subfigure, where¯ = 0.005. Sometimes, however, V grows before its decay -see for instance the middle-right and bottom-right subfigures, where¯ = 0.03 and 0.1 respectively. Such behaviour is possible when the polaron lifetime decays with A more quickly than the displacement does. When this happens, there may exist some optimal amplitude, A = A m , at which the polaron attains maximum velocity, V m . A m may coincide with A c -see for instance the top-right subfigure. Meanwhile, the middleright and bottom-right illustrate clearly that, for different values of T , the critical A c remains the same, whereas the optimal A m changes. Qualitatively speaking, the larger T is, the smaller A m is. Whilst the value of A c does not depend on T , it does depend on¯ -we see this by comparing any row of subfigures in Figure 7 to any other row. But how do A c and correlate? Our results show that, as¯ grows, A c drops, but crucially the combined magnitude comb =¯ + A c remains roughly constant. Specifically, in the top row we see¯ = 0.005 and A c = 0.167, in the middle row we havē = 0.030 and A c = 0.142, and the bottom row shows = 0.100 and A c = 0.072, each case giving comb = 0.172 when A reaches critical. Recall that, when using a straightforward constant forcing =¯ , there is also a threshold value for¯ , below which the polaron simply exhibits small oscillations, and above which the polaron moves at high speed but delocalises very quickly. It is noteworthy that this threshold is = 0.154 (given β = 0.6), which is significantly lower than the critical combined amplitude of 0.172. In other words, = 0.154 causes polaron displacement, = 0.153 does not; and if one wishes to add on a periodic component A sin(2πτ /T ) in order to move the polaron, one needs A ≥ 0.019, making¯ + A far exceed what¯ is required on its own to move the polaron. This phenomenon is observed across all values of β.
In practice, then, what would make a good combination of forcing parameters, which propel the polaron with decent speed but does not cause large energy dissipation too quickly? First of all, a large¯ results in a large velocity but an energetically unstable polaron which delocalises rapidly -so rapidly that it may move the polaron less far in its lifetime than a small¯ does. The middle column of Figure 7 precisely illustrates this point. Meanwhile, a small¯ results in long-living polarons which can move very far, because of how stable they are, but they would take more time to reach the same destination, compared to polarons under a large¯ . On balance, a moderate value of¯ such as 0.03 is preferable. Secondly, once a¯ is chosen, it remains to choose A and T , and it is obvious that the ideal choice of A is the optimal amplitude, A = A m . Meanwhile, if T is small, such as T = 100, A m would be large. On the other hand, if T is large, such as T = 2000, the value of the maximum velocity would be small. We observe both of these extremes very clearly in the middle-right subfigure of Figure 7. Once again, these observations are not specific to β = 0.6, but universal for all values of β. Overall, we believe that the best MSPF parameters which we have tested are such combinations where¯ ≈ 0.030, T ≈ 500, and A ≈ A m which, given β = 0.6, is A m = 0.157. We will discuss the relationship between A m and β in Section 4.3.
We note one anomaly which we observe in Figure 7 but did not expect. When¯ = 0.005 (top row), if T = 2000 and A is large enough, then the displacement (and therefore velocity) can take large negative values, meaning the polaron moves in the opposite direction to what we expected, and with large speeds. Whilst we are uncertain as to what causes this counter displacement, it is certainly another reason to reject small¯ and large T when choosing forcing parameters.

The relevance of β
To produce Figure 7, we fixed β = 0.6. How would the figure have looked if β had been different? Our results show that qualitatively it would exhibit the same behaviour, characterised by critical amplitudes A c , and optimal amplitudes A m . Quantitatively, the values of A c and A m would change. It is therefore natural to investigate how they change with β. After all, our generalisation to the Davydov-Scott model is manifest in the extra parameter β. Firstly we establish the following preliminary result.
Recall that the stationary polaron, upon which we impose the MSPF, is characterised by two quantities: its probability distribution, specifically its maximum localisation probability, max |ψ n | 2 , and its binding energy. These are in turn determined by the symmetry parameter β and effective coupling parameter λ. As β varies, so does the value of λ required to keep max |ψ n | 2 constant. This correlation is shown in Figure 8. It is clear that λ(β) is a decreasing function. We have made sure that whenever we altered β we also took λ = λ(β), so that all of our moving polarons begin as stationary states which share the same probability distribution. An alternative would have been to take whatever λ is required to keep the binding energy constant. Our results show that if we had decided to keep the binding energy constant at, say, −2.5, then max |ψ n | 2 would have been 0.53 at β = 0, or 0.81 at β = 1. It is a central feature of our model that two stationary polarons with the same maximum localisation probability need not have the same binding energy, and vice versa.
Having established the relation λ(β), we can study how A c and A m depend on β. We do so by fixing¯ and T , and working out what A c and A m are for various {β, λ(β)}. For instance, in Figure 7 we saw that if β = 0.6, λ = λ(0.6) = 5.0,¯ = 0.03 and T = 500, then A c = 0.142 and A m = 0.157. What if we fix¯ and T , and vary β? The result is displayed in Figure 9. When β = 0, we have A c = 0.12, and when β = 1 we have A c = 0.1. In fact,  A c is minimal when β = 1, which suggests that a system with antisymmetric electron-phonon interaction is most conducive to polaron displacement by MSPF. One might expect that a system with symmetric interaction would be least conducive to polaron displacement, and therefore A c should be maximal when β = 0. This is not the case. We observe that A c is maximal when β = 0.6, which models a system with moderately asymmetric electron-phonon interaction. Meanwhile, polaron velocity produced by critical forcing, V c , is maximal when β = 0.87.
The optimal amplitude, A m , varies little when β is less than 0.7, but decays sharply when β increases beyond 0.7, to such an extent that it almost equals the critical amplitude A c . The optimal velocity, V m , is typically of the same order of magnitude as V c . When β = 0.87, V m and V c almost coincide.
Earlier, based on Figure 7, we asserted that the period T of the MSPF has little effect on the value of A c . Indeed, our results show that, if we had produced Figure 9 with T fixed at either 100 or 2000, the A c curve would have been virtually unaffected. We also conjectured that the critical amplitude A c is negatively and linearly correlated with¯ , so that the combined magnitude comb =¯ + A c remains constant as¯ varies. This is supported by our results. Indeed, Figure 9 is produced with¯ fixed at 0.03; but if we had produced Figure 9 with¯ fixed at either 0.005 or 0.1, the A c curve would simply have been shifted along the vertical axis, by an amount equal to the difference between the new¯ and 0.03.

Dynamical polarons at non-zero temperature
We study the effect of random fluctuations which result from non-zero temperatures in the environment surrounding the lattice. The randomised forcing on the lattice is represented by the normally-distributed f n (τ ), which by definition (21) must have, for τ ≥ 0, the following first and second moments.
where θ is the dimensionless temperature, and Θ is the temperature. Compared to the dimensional F n (t) that we introduced in Section 2, the τ which now appears in f n (τ ) is a discrete index. We have followed the standard procedure of replacing δ(t − t ) by 1/Δτ , up to non-dimensionalisation constants. Beginning with a stationary polaron, which we computed numerically in Section 3.2, we integrate the system of equation (19) forward in time from τ = 0. Using a random number generation algorithm, we generate a new vector f n before each integration step. If θ is large, we find that it can cause large distortions in the lattice and rapid delocalisation of the polaron, due to excessive energy input to the system. For appropriate values of θ, we see that the polaron's binding energy undergoes small fluctuations, but on the average it tends to shift towards zero. After some time, the binding energy stabilises. For example, given β = 0.6 and λ = 5.0, the stationary polaron has binding energy −2.47.
Integrating from τ = 0 with θ = 0.0003, we find that after τ ∼ O(10 4 ) the binding energy settles, on average, around −2.15. The period of time required for a polaron to reach such a thermal equilibrium is the thermalisation phase of the polaron dynamics. During this phase, the forcing (τ ) on the electron is kept at zero, but the electron nevertheless undergoes small fluctuations around its initial position, due to its coupling to the thermalised lattice. Our results show that, irrespective of β, we are unable to raise θ > 0.001, because such a large θ induces excessive lattice distortions which cause the polaron to delocalise before reaching thermal equilibrium.
In each simulation, we integrate the system, with (τ ) = 0, until the polaron reaches thermal equilibrium. Then we reset τ = 0, and "turn on" the forcing (τ ) for τ ≥ 0. We examine how the polaron subsequently moves, under combinations of (τ ) and f n (τ ). Since the thermalisation phase raises the polaron energy, we expect that a thermalised polaron would be easier to displace, in the sense that it would require a smaller (τ ) to displace it, compared to the zero temperature case. Indeed, our results confirm this. To obtain our results in this section, every dynamical simulation, with a set of chosen parameters {β, λ(β),¯ , A, T, θ}, is run 100 times, and averages of quantities such as polaron lifetime and displacement are then taken. Figure 10 is to be compared directly with Figure 7, which contained results for β = 0.6 and θ = 0. Specifically, Figure 10 is to be compared with the solid (black) lines in the middle row of subfigures in Figure 7, for which two of the parameters in =¯ + A sin(2πτ /T ) were fixed: = 0.03, and T = 500. We saw that, given said parameter values, the critical amplitude was A c = 0.142. When we have a non-zero θ in the system, we define A c to be the smallest A for which the average polaron displacement (over 100 simulations) exceeds 10 lattice sites. According to this definition, when β = 0.6,¯ = 0.03, T = 500 and θ = 0.0001, we see in Figure 10 that A c = 0.121, which is significantly lower than the case of θ = 0. Fixing said values of¯ , T and θ, we find that the value of A c depends on β in a manner shown in Figure 11. That is, A c is minimal when β = 1, suggesting that an antisymmetric electron-phonon interaction makes it easiest to displace the polaron. Meanwhile, A c is maximal when β ≈ 0.5, suggesting that, counter-intuitively, what makes displacing the polaron most difficult is not a symmetric electron-phonon interaction, but a moderately asymmetric one. Indeed, this A c (β) function is very similar to the one in Figure 9, where we also had¯ = 0.03, T = 500 fixed, but θ = 0. Now with θ = 0.0001, the A c (β) curve in Figure 11 is significantly lower. This means that, regardless of β, a non-zero temperature makes it easier to displace the polaron by the forcing =¯ + A sin(2πτ /T ), in the sense that a smaller combined magnitude¯ +A is required. It is also noteworthy that, under a non-zero temperature, the onset of polaron motion is more gradual, in the sense that a critical amplitude results in a very small velocity, V c . Indeed, comparing V c in Figure 9 with V c in Figure 11, we see that the latter is 2 orders of magnitude smaller.
More can be said about Figure 10. When we raise the temperature to θ = 0.0005, we find that the polaron is displaced (on average) by hundreds of sites even if A = 0. This suggests that, given β = 0.6 and¯ = 0.03, there exists some critical temperature θ = θ c between 0.0001 and 0.0005, for which just the combination of (τ ) =¯ and f n (τ ) is sufficient to displace the polaron, and no periodic component in (τ ) is needed. θ c is critical in the sense that,  if θ is any lower than θ c , then the combination of (τ ) =¯ and f n (τ ) does not energise the polaron enough to move it, and a non-zero A is required. Indeed our results show that, given β = 0.6 and¯ = 0.03, the critical temperature is θ c = 0.00032. Furthermore, we have investigated how θ c changes as we vary β and¯ , and the results are shown in Figure 12. We observe the qualitative trend that, the larger¯ is, the less thermal energy is required to make up for the extra energy that the polaron needs in order to move. We also observe that, in general, the larger β is, the less thermal energy is required to displace the polaron. This fits in nicely with our understanding that, when β is close to 1, we have an electron-phonon interaction which is biased towards one end of the lattice, making the electron more susceptible to displacement. In Figure 13 we present a typical polaron trajectory when θ = θ c . Under this critical temperature, some simulations would produce no polaron displacement at all, but most trajectories would be similar to that in Figure 13, clearly showing a directed movement. We propose to explain the shape of these trajectories as follows. First of all, the combined magnitude of the MSPF would be much lower than what is required to move the polaron under zero temperature. In Figure 13 for example, we havē + A = 0.03, whereas the critical value under zero temperature, as we discovered in Section 4.2, is¯ + A = 0.172. Even when there is a non-zero temperature, the polaron still spends the majority of its lifetime oscillating around its localisation site by small amounts. However, occasionally the random forces on the lattice sites in the vicinity of the electron causes a large distortion, such that the effective potential barrier for the electron is significantly lowered, and the electron can escape the well. Once it does that, it is propelled towards one end of the lattice by¯ . But before the electron has time to move far, the random forces may have further distorted the local lattice sites in such a way that a high potential barrier is restored. This then traps the electron again, giving the polaron time to recover its integrity, before the next random time at which the electron jumps out of its potential well. This explains why a trajectory under the critical temperature appears jagged, showing the polaron "hopping" one or two sites at a time, in stark contrast with the smooth and regular polaron motion exhibited in Figure 6.

Discussions and conclusions
In this study we have presented a new mathematical model describing polaron dynamics in linear peptide chains. The model is dependent on a symmetry parameter, β, which measures the extent to which the interaction between the polaron's electron and phonon components is spatially symmetric. We have shown that when β takes its extreme values, 0 and 1, the model reduces to existing ones for which it was assumed that the electron-phonon interaction was, respectively, symmetric and antisymmetric. We have justified the physical neccessity of including β in the model, in that one should not simply assume the electron to be coupled equally strongly to lattice points on either side, or to be coupled only to the lattice point on one side. Instead, the spatial symmetry should be determined by the adjustable parameter β.
Apart from β, we have also identified two composite parameters which are most vital to the intrinsic properties of the polaron. Firstly there is the adiabaticity parameter, ρ, measuring the characteristic time scale separation between the electron and phonon, which we justifiably fixed throughout the study. Then there is the effective coupling parameter, λ, measuring the strength of the electron-phonon interaction. The combination of β and λ determines the two aspects of the stationary polaron: its maximum localisation probability, and its binding energy. We have computed both of these quantities as functions of β and λ. Moreover, in the infinite lattice limit, we have obtained stationary polaron solutions by analytically integrating the system, and the results are in good agreement with our numerical solutions on a finite lattice.
Our main results relate to using an external forcing to displace the polaron, in a manner which causes minimal energy loss and which, crucially, is directed. Such polaron dynamics could be achieved only if the electron is dislodged from its self-trapping potential well, and the local lattice distortions propagate coherently with the electron, and some mechanism exists which ensures the electron always moves towards one end of the lattice. If the second condition is not met, then over time the electron probability density function would become broader, leading to delocalisation of the polaron. We have found that a constant external force,¯ , on the electron is insufficient for displacing the polaron, unless¯ is larger than some threshold value, but then the forcing causes rapid energy loss and delocalisation. We have also found that a sinusoidal force, A sin(2πτ /T ), on the electron is never sufficient for displacing the polaron, throughout the range of A that we tested. We then combined the constant and sinusoidal forces, resulting in the mean-shifted periodic forcing (MSPF), =¯ + A sin(2πτ /T ). We have discovered that, for each¯ which is insufficient on its own to displace the polaron, there is some critical value A c , such that the polaron is displaced if and only if A ≥ A c . There is also an optimal value A m , such that the polaron attains maximum velocity at A = A m . The value of A c is irrespective of the period T , whilst A m is negatively correlated with T . As¯ is decreased, A c increases, in such a way that the combined magnitude¯ + A c remains constant. This suggests that there is a certain amount of extra energy that the electron needs in order to overcome the polaron binding, and how much of it comes from the constant or sinusoidal part is inconsequential, as long as the two parts combine to a large enough overall amplitude. Nevertheless, the split between¯ and A does determine the manner in which the polaron propagates, specifically its velocity and stability. The velocity is predominantly determined by¯ , and positively correlated with it; but the stability of the polaron is negatively correlated with¯ . By comparing three sets of {¯ , A} with the same combined + A, namely {0.005, 0.167}, {0.03, 0.142}, {0.1, 0.072} (while keeping all other parameters fixed), we found that the combination {0.03, 0.142} produces optimal balance between polaron velocity and stability.
We have examined how the aforementioned phenomena depends upon β. To do so, we needed a way of isolating the effect of varying β. This posed a difficulty, because if we were to fix λ and vary β then both the maximum localisation probability and binding energy of the stationary polaron would change. We would then be comparing dynamical behaviours of dissimilar polarons. We therefore decided to vary λ with β, in a way that allowed us to generate a set of stationary polarons, one for each combination of {β, λ(β)}, such that they all had the same maximum localisation probability. Then we launched these polarons using the same external forcing and compared the results. We have found that β = 1, representing a spatially antisymmetric electron-phonon interaction, produces a polaron which is easiest to move, in the sense that the least amount of forcing is required. We have also found that the symmetric model, β = 0, does not make the polaron most difficult to move (β ≈ 0.6 does that). This hints at the existence of some intrinsic mechanism in the β = 0 model which pushes the electron towards one end of the lattice, despite it being coupled to the other end equally strongly.
We have also studied the MSPF under non-zero temperatures, θ > 0. The manifestation of thermal effects is random forces on the lattice points. We have found that a non-zero θ facilitates polaron propagation, in the sense that it lowers the critical amplitude A c , for any given¯ . Moreover, a non-zero θ results in a gradual onset of polaron motion, meaning the rate of change of polaron velocity with respect to A near A = A c is small, compared to the onset under θ = 0. Our results have also shown that, whenever there is polaron propagation, whether θ = 0 or θ > 0, the relative displacements between neighbouring lattice points remain under O(10 −2 ). This is a necessary condition which allows us to model the lattice points as point dipoles.
Some of the choices of parameters in the MSPF may be justified physically as follows. It is well known that across the plasma membrane of a living cell, a resting membrane potential is maintained by intercellular chemical processes [28]. It is also well known that within the plasma membrane there exist highly stable transmembrane regions of proteins, for instance the human prolactin receptor 2N7I [29], and the rat monoamine oxidase A 1O5W [30], both of which are α-helical structures spanning the entire membrane width. Given a constant potential difference of ΔV across a linear, homogeneous, isotropic dielectric medium with constant width d, the effective electric field inside the medium is given by E 0 = ΔV /(κd), where κ is the dielectric constant (a.k.a. relative permittivity) of the medium [31]. Assuming the plasma membrane is such a medium, we can then attribute the physical origin of¯ to the resting membrane potential, and calculate¯ using equation (21), namelȳ = qE 0 R/( Ω). For many cells the values of the E 0 and d are well established. As an example, one may look at human red blood cells (erythrocytes), one of the most widely studied cells in nature, and point to [32] for the value E 0 = −8.4 mV, as well as [33,34] for d = 78Å. However, the value of κ for a membrane is highly contentious, due to the fact that it depends sensitively upon a large variety of biophysical attributes of the membrane, such as hydration [35], pH value [36], and structural stability [37]. In a recent review, it was reported that the value of κ in literature ranges from 1 to 40 [38]. Feeding these values of E 0 , d and κ into the equation for¯ , we find that¯ ranges from 0.0033 to 0.13. We have taken care to ensure that in this study the values of¯ falls strictly within this range. For a physical origin of the periodic term, A sin(2πτ /T ), one could look to common electromagnetic radiations which fill the environment around us in the modern age, such as the radiation from telecommunication transmitters. In particular, the values of T which we have considered, 100, 500 and 2000, respectively match the frequencies of the IEEE 802.11ad protocal Wi-Fi band, the K u band frequencies for satellite communications and broadcasting, and the UHF band frequencies for cellular communications [39,40]. However, the amplitudes of the aforementioned radiations are much smaller than the values of A for which we have observed noteworthy results. For instance, treating the mobile telephone transmitter as an omni-directional dipole with peak power P , we can estimate the amplitudeÃ of its output waves at operational distance d, by using the well-known formula P/(4πd 2 ) = 0 cÃ 2 /κ, where 0 , c, κ are the vacuum permittivity, speed of light, and relative permittivity of the medium, respectively. Feeding P = 1 W [41] and κ ≤ 40 [38] into the formula, we obtain dimensionless A ≤ 6 × 10 −6 /(d/m). This means that in order to obtain A = 0.1, one needs the operational distance d to be O(10 −5 ) metres, which is unrealistic. It is therefore clear that, in a real cell environment, the effect on a polaron due to a combination of resting membrane potential and random thermal forces are dominant over any external electromagnetic radiation that may commonly be present. This highlights the importance of our observation that, given a constant electric field¯ , there exists some critical temperature θ c , such that the polaron undergoes directed drift if and only if θ ≥ θ c . In other words, a combination of¯ and θ can be sufficient for displacing the polaron, in a manner which causes minimal energy loss and which is directed, just as a combination of¯ and A sin(2πτ /T ) can. It has been reported that stationary polarons formed on 3-dimensional lattices, such as an α-helix, are more strongly bound and therefore can be stabler during propagation [42]. It is our hope that our model can be adapted to study such 3-D systems, and that the stabilising effect of the helical geometry will enable us to raise θ, from our current values of O(10 1 )K to physiological temperatures.
We would both like to thank Larissa Brizhik for kindly answering some questions.