text:
# Robust control of decoherence in realistic one-qubit quantum gates

## I Introduction

Quantum computing holds an extremely appealing practical promise and exerts a no less remarkable intellectual fascination. Indeed, the quantum algorithms proposed until now could simplify considerably the computational complexity of classically hard problems. While their actual implementation in robust and realistically sized quantum computers remains still a promise, there is general consensus in expecting further significant advances both in developing new quantum algorithms and the hardware they could be implemented in Pitt ; NC .

In the presently accepted paradigm, a quantum computation is a sequence of unitary transformations effected by the quantum computer upon a state in a finite dimensional Hilbert (i.e. a complex Euclidian) space. At the end of the computation, the system is left in a state which, upon measurement (read-out), would yield with high probability the desired result. While the read-out is non-invertible, quantum computation per se is deterministic and reversible. Thus, the realization of any quantum computer relies on maintaining quantum coherence in the system, for a period of time spanning at least the duration of the desired computation. Unfortunately, due to coupling to the environment which is not explicitly accounted for in quantum computation, the evolution within the quantum computer loses its unitary character; in other words, the system decoheres Pitt ; NC ; Gardiner91 ; Giulini ; BEZ ; Weiss99 .

Any quantum computation, no matter how general or convoluted, can be approximated to arbitrary precision using a universal set of quantum gates e.g., the discrete set of gates made up of the CNOT, Hadamard, phase, and $`\pi /8`$ gates Pitt ; NC . Thus, a successful approach to the decoherence problem in any quantum computer requires to address and solve the problem at the elementary gate level. To be specific, in this paper we shall focus on the NOT gate, which will serve as an illustration for any one-qubit gate. The analysis of the two-qubit CNOT gate will be considered elsewhere.

The NOT gate acts on a single qubit, which is the most general superposition state of a two-level system. Since the set containing the identity and the Pauli matrices, $`\{I,\sigma _x,\sigma _y,\sigma _z\}`$, forms a basis in the space of $`2\times 2`$ matrices, the coupling of a two-level system to the outside world can be essentially realized through the matrices $`\sigma _z`$ and $`\sigma _x`$. These two couplings lead to different effects, known as adiabatic and thermal decoherence, respectively.

In general, quantum decoherence is a very fast process, whose speed depends primarily on the size and temperature of the computer, but may depend significantly on other factors as well, such as imperfect state preparation, undesired excitations, inaccuracies in manipulations of the logic gates, explicit stochasticity and noise, etc. For example, decoherence is induced by the spontaneous emission of a photon from a two-level atom, where the environment is represented by the weakly-coupled electromagnetic field. Whereas the evolution of the two-level system and its environment considered together as a large (closed) quantum system is in principle unitary, the subdynamics of the two-level system generally is not. Clearly, this behavior represents a major stumbling block on the path to realizing a usable quantum computer at room temperature.

Several approaches have been proposed to eliminate or mitigate the undesirable effects of decoherence in open quantum systems, including open loop (quantum bang-bang) control Tombesi\_ol ; Viola98 ; Viola99a ; Viola99b , decoherence free subspaces (DFS) Zanardi97 ; Lidar98 , quantum error correction Preskill98 ; Knill00  and quantum feedback Tombesi\_qf ; Wang01 .

The open loop (quantum bang-bang) control has been pioneered by the seminal work of Lloyd and Viola (see Refs. Viola98 ; Viola99a ; Viola99b  and references therein). Vitali and Tombesi Tombesi\_ol  have also considered applying a sequence of frequent parity kicks, as well as an appropriate stochastic modulation, in order to control decoherence. Essentially this type of approach relies on applying rapid control pulses to the system in order to dynamically cancel the system-environment interaction. It has been shown that decoherence can be effectively suppressed if the pulse rate is much higher than the frequency corresponding to the correlation time of the environment. Thus, this open-loop control technique ”freezes” decoherence in a manner analogous to the quantum Zeno effect.

DFS are special subspaces of the Hilbert space, which are inherently immune to decoherence, due to specific symmetries in the system-environment interaction. This concept was first introduced by Zanardi and Rasetti Zanardi97 , and has been recently generalized to noiseless subsystems Knill00 . It has also been shown Viola00 ; WL  that bang-bang control can be used to create the conditions required for large DFS. The successful use of decoherence free subspaces (DFS) to passively stabilize a quantum state has been reported for photon singlets Kwiat00 , for a pair of trapped ions Kielpinski01 , and more recently for a noiseless subsystem via liquid state NMR Viola01 .

Quantum error correction is historically the first method proposed to counter decoherence, and also the only one to date for which an accuracy threshold theorem has been proved Knill96 . It is essentially an active form of quantum feedback control implemented on a system encoding redundant information. Valid states of the system are restricted to special Hilbert subspaces (codes) so that any loss of information due to decoherence (errors) can be diagnosed and reversed. The system is closely monitored to observe errors, which can be corrected later by conditionally performing appropriate feedback operations. It has been shown that fault-tolerant recovery which allows for errors in the correcting feedback, if the error rate is below a critical value.

The efficient implementation of both DFS and quantum error correction techniques requires a large amount of ancillary quantum memory. Thus quantum feedback (closed-loop) and open-loop control strategies have the advantage that in general they do not require any ancillary quantum resources.

Quantum feedback can also be used to stabilize a quantum state against decoherence, as first shown by Tombesi et al. in a series of papers Tombesi\_qf . For example, in a scheme recently proposed by Wang et al. Wang01  for a two-level atom, photocurrent feedback from homodyne detection of spontaneous emission is used to alter the atom-environment dynamics. The quantum state of the atom is driven to and then remains in a pure target state, even though the atom is still spontaneously emitting. Of course, the feedback delay time has to be shorter than the time scale corresponding to the spontaneous emission rate, in order to successfully maintain the atom in the target state. We note that this quantum feedback scheme cannot be used to maintain an unknown state such as that obtained from a quantum computation, it can only drive the system to a known target state.

In this paper, an open-loop control scheme is proposed and simulated, to eliminate the effects of decoherence in a two-level system. While the general philosophy is following in the steps of Lloyd and Viola’s seminal results Viola98 ; Viola99a ; Viola99b , our implementation has several new aspects. More precisely, we present the first implementation of a quantum control scheme in which: (i) decoherence and control are taken to act simultaneously within a realistic model, which allows one to deal with thermal decoherence; (ii) the required control is directly related to and calculated from the decoherence effects, which presents the practical advantage of maintaining the frequency and amplitude of the required controls at manageable levels; and (iii) the effect of imperfect control pulses is assessed. We are aware that determining the decoherence behavior of a quantum system can be difficult in practice, but if we know it then our control scheme can take advantage of this, to tailor more efficient controls.

The remaining of the paper is organized as follows: In section 2 we present the spin-boson model for a two-level system in contact with the environment and an external control. In Section 3 we calculate the reduced density matrices for the two level system after tracing out the environment modes, in the case of adiabatic and thermal decoherence. The control strategy is implemented in Section 4 and the results are discussed in Section 5.

## II Spin-Boson Model

The exact quantum mechanical description of the system starts from the complete Hamiltonian that accounts for the two-dimensional system and the environment which, in principle, is the rest of the Universe. In different situations the environment predominantly manifests itself in different ways, e.g., electromagnetic modes, acoustic modes, etc. Thus, in many instances, a reasonable and generic model is represented as a system of noninteracting boson modes. The two-level system interacting with these modes is known as the spin-boson model Weiss99 . We write the complete Hamiltonian in the form:

$$H=H_s+H_e+H_i+H_c$$
(1)
where

$$H_s=\underset{i=1}{\overset{2}{}}P_{ii}E_i$$
(2)
and

$$H_e=\mathrm{}\underset{q=1}{\overset{\mathrm{}}{}}\omega _{0q}a_q^{}a_q$$
(3)
represent the two-level system and the boson bath. Here $`P_{ii}`$ denote the projection operators $`P_{ii}=|ii|`$, and $`E_i`$ are the corresponding energies. The terms

$$H_i=\mathrm{}ϵ(\alpha \sigma _x+\beta \sigma _z)\underset{q=1}{\overset{\mathrm{}}{}}(\mathrm{\Omega }_q^{}a_q+\mathrm{\Omega }_qa_q^{})$$
(4)
and

$$H_c=\mathrm{}\mathrm{\Omega }_FV(t)\sigma _x\mathrm{cos}(\omega _0t)$$
(5)
describe the interaction of the two level system with the environment and the external control, $`V(t)`$, respectively. The interaction with the environment is parametrized by $`ϵ`$ (whose magnitude indicates the strength of the coupling) and can result in a phase flip or a bit flip (a bit-phase flip can also be obtained either by combining the two effects or by including a $`\sigma _y`$ term). However, from a formal viewpoint, the effects of $`\sigma _x`$ and $`\sigma _y`$ are essentially similar and we shall not consider the latter. On the other hand, the control which is applied as a fast laser pulse, is assumed to act strongly only upon the spin, without significantly affecting the environment, either directly or indirectly.

While the spin-boson model in various guises has been thoroughly studied Leggett87 ; Kampen95 ; Palma96 ; Mozyrsky98 , the consistent inclusion of the control term, $`H_c`$, in the Hamiltonian, and the general derivation of the corresponding reduced description are new.

The dynamics of the two-level system considered in this paper is expressed in terms of parameters normalized to the Rabi frequency, $`\mathrm{\Omega }_F`$, namely the time evolution is described in terms of dimensionless Rabi time units, $`\tau =\mathrm{\Omega }_Ft`$, and all the frequencies, namely $`\omega _q=\omega _{oq}/\mathrm{\Omega }_F`$, $`\omega =\omega _o/\mathrm{\Omega }_F`$, $`g_q=\mathrm{\Omega }_q/\mathrm{\Omega }_F`$ and $`\omega _{12}=\mathrm{\Omega }_{12}/\mathrm{\Omega }_F`$ are also renormalized with respect to the Rabi frequency. The frequency $`\mathrm{\Omega }_{12}`$ is the transition frequency between the energy levels of the two-level system.

The state of the system at $`\tau =0`$ is described by the density matrix:

$$\rho (s,e,0)=\rho (s,0)\rho (e,0)$$
(6)
where

$$\rho (e,0)=\underset{q}{}\rho _q(e,0)$$
(7)
$$\rho _q(e,0)=[1e^{\frac{\mathrm{}\omega _{oq}}{kT}}]^1\underset{n_q}{}e^{\frac{\mathrm{}\omega _{oq}}{kT}n_q}|n_qn_q|$$
(8)
and

$$\rho (s,0)=\underset{i,j=1}{\overset{2}{}}\rho _{ij}(0)P_{ij}.$$
(9)
In general, the decoherence experienced by a two-level system has two distinct timescales, corresponding to the adiabatic and thermal regimes separately.

Adiabatic decoherence is responsible for the decay of transverse $`(x,y)`$ polarization (the off-diagonal density matrix elements). It acts on a relatively short timescale, such that a superposition state will typically decay to a mixture of states. For larger, macroscopic sub-systems, adiabatic decoherence ensures that quantum superpositions of distinct states will not be observed. The mixture of states generated by adiabatic decoherence will then continue to relax to a stationary value over a longer thermal timescale.

Thermal decoherence generally causes a change of all density matrix elements, and leads to the exponential decrease of the excited state population of the two-level system, due to energy exchanges between the two-level system and its environment. It is commonly neglected in discussing decoherence of quantum systems since its timescale is longer than the timescale for adiabatic decoherence, which is deemed responsible for eliminating most quantum effects ZHP . We note that often thermal decoherence is refered to as (quantum) dissipation, and adiabatic decoherence is simply refered to as decoherence.

To simplify the calculations and render them more transparent, it is convenient to write the evolution of the system in the interaction representation, and consider separately the adiabatic ($`\alpha =0`$, $`\beta =1`$) and thermal ($`\alpha =1`$, $`\beta =0`$) decoherence regimes. A combined treatment of these regimes is straightforward, though computationally cumbersome. The unitary evolution operator in the interaction representation, $`U_I`$, is related to the unitary evolution operator in the Schrödinger representation, $`U`$, by

$$U_I(s,e,\tau )=e^{iH_0\tau /\mathrm{}}U(s,e,\tau )$$
(10)
where $`H_0=H_s+H_e`$.

Thus the evolution of the density matrix of the system and environment is

$$\frac{\rho }{\tau }=\frac{i}{\mathrm{}}[H_I,\rho ]$$
(11)
where the interaction Hamiltonian contains both the interaction with the environment and with the applied control:

$$H_I=H_{Ic}+H_{Ii},$$
(12)
$$H_{Ic}=\frac{\mathrm{}}{2}V(\tau )\sigma _x,$$
(13)
$$H_{Ii}=\{\begin{array}{cc}ϵ\mathrm{}\sigma _z_q(g_q^{}a_qe^{i\omega _q\tau }+g_qa_q^{}e^{i\omega _q\tau })\hfill & \alpha =0\text{}\beta =1\hfill \\ ϵ\mathrm{}[\frac{\sigma _x+i\sigma _y}{2}_qg_q^{}a_qe^{i(\omega _{12}\omega _q)\tau }+\frac{\sigma _xi\sigma _y}{2}_qg_qa_q^{}e^{i(\omega _{12}\omega _q)\tau })]\hfill & \alpha =1\text{}\beta =0\hfill \end{array}$$
(14)
These interaction-picture Hamiltonians have been calculated by using the rotating wave approximation, and assuming a zero detuning $`\delta =\omega _{12}\omega =0`$.

The evolution equation for $`\rho `$ admits the formal solution

$$\rho (s,e,\tau )=U_I(s,e,\tau )\rho (s,e,0)U_I^{}(s,e,\tau )$$
(15)
where the evolution operator satisfies the equation

$$\frac{dU_I(s,e,\tau )}{d\tau }=\frac{i}{\mathrm{}}[H_{Ic}+H_{Ii}]U_I(s,e,\tau )$$
(16)
given the initial condition $`U_I(s,e,0)=1`$. The formal solution of this equation is Scully\_qo ,

$$U_I(s,e,\tau )=𝒯[\mathrm{exp}\{\frac{i}{\mathrm{}}_0^\tau 𝑑\tau ^{}(H_{Ic}(\tau ^{})+H_{Ii}(\tau ^{}))\}],$$
(17)
where $`𝒯[]`$ represents a time-ordering operator.

Since the various terms in the interaction Hamiltonian $`H_I`$ do not commute, the complete evolution operator, $`U_I(s,e,\tau )`$, cannot be calculated exactly, except in special cases. For instance, Viola and Lloyd Viola98  have shown that if the interaction with the environment $`H_{Ii}`$ contains only $`\sigma _z`$, the evolution operator $`U_I(s,e,\tau )`$ can be substantially simplified by using $`\pi `$-pulses for the control. However this type of control would not allow the same simplification if the interaction with the environment $`H_{Ii}`$ contains $`\sigma _x`$ or $`\sigma _y`$.

The general Baker-Hausdorff theorem Gardiner\_qn  can be used to expand the evolution operator into an infinite product of exponentials,

$$U_I(s,e,\tau )=e^{\frac{i}{\mathrm{}}_0^\tau 𝑑tH_{Ic}(t)}\times e^{\frac{i}{\mathrm{}}_0^\tau 𝑑tH_{Ii}(t)}\times e^{(\frac{i}{\mathrm{}})^2_0^\tau 𝑑t_0^t𝑑t^{}[H_{Ic}(t),H_{Ii}(t^{})]}\times \mathrm{}$$
(18)
We assume that the effect of the control and interaction Hamiltonians is relatively small at all times. Neglecting the commutators, the evolution operator can be written as a first order approximation in the magnitude of the control pulses, $`V(\tau )`$, and the coupling strength parameter, $`ϵ`$, of the system-environment interaction,

$$U_I(s,e,\tau )\{\begin{array}{cc}exp\{\frac{i}{\mathrm{}}_0^\tau 𝑑tH_{Ic}(t)\}\times exp\{\frac{i}{\mathrm{}}_0^\tau 𝑑tH_{Ii}(t)\}\hfill & \alpha =0\text{}\beta =1\hfill \\ exp\{\frac{i}{\mathrm{}}_0^\tau 𝑑t(H_{Ic}(t)+H_{Ii_x}(t))\}\times exp\{\frac{i}{\mathrm{}}_0^\tau 𝑑tH_{Ii_y}(t)\}\hfill & \alpha =1\text{}\beta =0\hfill \end{array}$$
(19)
where $`H_{Ii_x}`$ and $`H_{Ii_y}`$ correspond to the terms in the system-environment interaction Hamiltonian $`H_{Ii}`$, proportional to $`\sigma _x`$ and $`\sigma _y`$ respectively.

Performing the time integrals in the exponents, we obtain

$$U_I(s,e,\tau )=\{\begin{array}{cc}\mathrm{exp}\{i\sigma _xI(\tau )\}\times \mathrm{exp}\{2\sigma _zQ_{}(\tau )\}\hfill & \alpha =0\text{}\beta =1\hfill \\ \mathrm{exp}\{i\sigma _xI(\tau )\sigma _xQ_{}(\tau )\}\times \mathrm{exp}\{i\sigma _yQ_+(\tau )\}\hfill & \alpha =1\text{}\beta =0\hfill \end{array}$$
(20)
where

$$I(\tau )=\frac{1}{2}_0^\tau 𝑑\tau ^{}V(\tau ^{})$$
(21)
$$Q_{}(\tau )=\frac{ϵ}{2}\underset{q}{}(M_qa_q^{}M_q^{}a_q)$$
(22)
$$Q_+(\tau )=\frac{ϵ}{2}\underset{q}{}(M_qa_q^{}+M_q^{}a_q)$$
(23)
and

$$M_q(\tau )=\{\begin{array}{cc}\frac{g_q}{\omega _q}(1e^{i\omega _q\tau })\hfill & \alpha =0\text{}\beta =1\hfill \\ \frac{g_q}{\omega _{12}\omega _q}(1e^{i(\omega _{12}\omega _q)\tau })\hfill & \alpha =1\text{}\beta =0\hfill \end{array}$$
(24)
## III Evolution of the Reduced Density Matrix

To evaluate the various effects of decoherence we have to calculate the reduced density matrix of the two-level system, by tracing out the environment modes:

$$\rho _{ij}(s,\tau )=Tr_e\{\underset{k,l=1}{\overset{2}{}}i|U_IP_{kl}\rho (e,0)U_I^{}|j\rho _{kl}(0)\}$$
(25)
for $`i,j=1,2`$.

### III.1 Adiabatic case $`(\alpha =0,\beta =1)`$

Upon expansion of the exponentials in Eq. (20), we obtain the following approximate expressions for the evolution operator and its adjoint:

$$U_I(s,e,\tau )=E_0+iE_x\sigma _x+E_y\sigma _y+E_z\sigma _z$$
(26)
and

$$U_I^{}(s,e,\tau )=E_0iE_x\sigma _xE_y\sigma _yE_z\sigma _z$$
(27)
where

$`E_0`$ $`=`$ $`\mathrm{cosh}(2Q_{}(\tau ))\mathrm{cos}(I(\tau ))`$ (28)
$`E_x`$ $`=`$ $`\mathrm{cosh}(2Q_{}(\tau ))\mathrm{sin}(I(\tau ))`$ (29)
$`E_y`$ $`=`$ $`\mathrm{sinh}(2Q_{}(\tau ))\mathrm{sin}(I(\tau ))`$ (30)
$`E_z`$ $`=`$ $`\mathrm{sinh}(2Q_{}(\tau ))\mathrm{cos}(I(\tau )),`$ (31)
with $`I(\tau )`$ and $`Q_{}(\tau )`$ defined by Eqs. (21)-(22).

Using these expressions in Eq. (25), we get the reduced density matrix of the two-level system:

$`\rho _{11}`$ $`=`$ $`\rho _{11}(0)\mathrm{cos}^2I+\rho _{22}(0)\mathrm{sin}^2Ii[\rho _{12}(0)\rho _{21}(0)]e^{g_{ad}}\mathrm{cos}I\mathrm{sin}I,`$ (32)
$`\rho _{22}`$ $`=`$ $`\rho _{22}(0)\mathrm{cos}^2I+\rho _{11}(0)\mathrm{sin}^2I+i[\rho _{12}(0)\rho _{21}(0)]e^{g_{ad}}\mathrm{cos}I\mathrm{sin}I,`$ (33)
$`\rho _{12}`$ $`=`$ $`\rho _{12}(0)e^{g_{ad}}\mathrm{cos}^2I+\rho _{21}(0)e^{g_{ad}}\mathrm{sin}^2I+i(\rho _{22}(0)\rho _{11}(0))\mathrm{cos}I\mathrm{sin}I,`$ (34)
$`\rho _{21}`$ $`=`$ $`\rho _{21}(0)e^{g_{ad}}\mathrm{cos}^2I+\rho _{12}(0)e^{g_{ad}}\mathrm{sin}^2Ii(\rho _{22}(0)\rho _{11}(0))\mathrm{cos}I\mathrm{sin}I,`$ (35)
where $`I=I(\tau )`$, and

$$g_{ad}:=g_{ad}(\tau )=\gamma _0^{\mathrm{}}𝑑\omega G(\omega )(1\mathrm{cos}\omega \tau )\mathrm{coth}\frac{\beta _0\omega }{2}$$
(36)
is the decoherence function obtained by Palma Palma96 . The dimensionless constant $`\gamma `$ depends on the dipole moment of the two-level system and on the Rabi frequency. The strength of the decoherence rate experienced by the two-level system is determined by $`\gamma `$, ie., the decoherence rate is weak for $`\gamma <<1`$.

The function $`G(\omega )`$ is the spectral function

$$G(\omega )=\omega ^{n2}e^{\omega /\omega _c},$$
(37)
where $`n`$ is the dimensionality of the system (usually $`n=3`$), and $`\omega _c`$ is the (usually very large) cutoff frequency, which ensures the convergence of the improper integral. In the case of phonon modes this cutoff frequency is the Debye frequency. We note that the only dependence on the two-level system is contained in $`\gamma `$.

The expressions (32)-(35) for the reduced density matrix elements $`\rho _{ij}`$ $`(i,j=1,2)`$ retain the correct limit behaviors, and the trace of the density matrix is always equal to one. Indeed, in the limit of zero control we recover Palma’s result Palma96 , whereby only the non-diagonal elements are affected by the decoherence process, while the diagonal elements remain unchanged. When both the decoherence and the control are equal to zero, we recover the ideal situation in which the state remains unchanged.

### III.2 Thermal decoherence $`(\alpha =1,\beta =0)`$

Proceeding as before, we use the approximation in Eq. (20), for $`\alpha =1`$ and $`\beta =0`$, to obtain the following expressions for the evolution operator and its adjoint in the thermal decoherence case:

$$U_I(s,e,\tau )=E_0+iE_x\sigma _x+E_y\sigma _y+E_z\sigma _z$$
(38)
and

$$U_I^{}(s,e,\tau )=E_0iE_x\sigma _xiE_y\sigma _y+iE_z\sigma _z$$
(39)
where

$`E_0`$ $`=`$ $`\mathrm{cos}(I(\tau )+iQ_{}(\tau ))\mathrm{cos}(Q_+(\tau ))`$ (40)
$`E_x`$ $`=`$ $`\mathrm{sin}(I(\tau )+iQ_{}(\tau ))\mathrm{cos}(Q_+(\tau ))`$ (41)
$`E_y`$ $`=`$ $`\mathrm{cos}(I(\tau )+iQ_{}(\tau ))\mathrm{sin}(Q_+(\tau ))`$ (42)
$`E_z`$ $`=`$ $`\mathrm{sin}(I(\tau )+iQ_{}(\tau ))\mathrm{sin}(Q_+(\tau )),`$ (43)
with $`I(\tau )`$, $`Q_{}(\tau )`$ and $`Q_+(\tau )`$ defined by Eqs.(21)-(23). To calculate $`U_I^{}(s,e,\tau )`$, we have used the fact that $`Q_{}(\tau )`$ and $`Q_+(\tau )`$ commute in the first order approximation for the coupling strength parameter, $`ϵ`$.

The environment modes are then traced out to obtain the elements of the reduced density matrix:

$`\rho _{11}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[1+(\rho _{11}(0)\rho _{22}(0))e^{2g_{th}}\mathrm{cos}(2I)i(\rho _{12}(0)\rho _{21}(0))e^{g_{th}}\mathrm{sin}(2I)]`$ (44)
$`\rho _{12}`$ $`=`$ $`[Re\{\rho _{12}(0)\}+Im\{\rho _{12}(0)\}\mathrm{cos}(2I)]e^{g_{th}}{\displaystyle \frac{i}{2}}[(\rho _{11}(0)\rho _{22}(0))e^{2g_{th}}\mathrm{sin}(2I)]`$ (45)
$`\rho _{21}`$ $`=`$ $`[Re\{\rho _{21}(0)\}+Im\{\rho _{21}(0)\}\mathrm{cos}(2I)]e^{g_{th}}+{\displaystyle \frac{i}{2}}[(\rho _{11}(0)\rho _{22}(0))e^{2g_{th}}\mathrm{sin}(2I)]`$ (46)
$`\rho _{22}`$ $`=`$ $`{\displaystyle \frac{1}{2}}[1(\rho _{11}(0)\rho _{22}(0))e^{2g_{th}}\mathrm{cos}(2I)+i(\rho _{12}(0)\rho _{21}(0))e^{g_{th}}\mathrm{sin}(2I)]`$ (47)
where the decoherence function, $`g_{th}`$, is given by

$$g_{th}:=g_{th}(\tau )=\gamma _0^{\mathrm{}}𝑑\omega \frac{1cos[(\omega _{12}\omega )\tau ]}{(\omega _{12}\omega )^2}\omega ^3\mathrm{coth}(\beta _0\omega /2)\mathrm{exp}(\omega /\omega _c)$$
(48)
In the limit of zero control, the matrix elements read:

$`\rho _{11}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1+e^{2g_{th}}(\rho _{11}(0)\rho _{22}(0)))`$ (49)
$`\rho _{12}`$ $`=`$ $`\rho _{12}(0)e^{g_{th}}`$ (50)
$`\rho _{21}`$ $`=`$ $`\rho _{21}(0)e^{g_{th}}`$ (51)
$`\rho _{22}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1e^{2g_{th}}(\rho _{11}(0)\rho _{22}(0)))`$ (52)
## IV Control Strategy

Our control strategy is based on the following idea: by equating the elements of the reduced density matrix in the presence of decoherence (adiabatic or thermal) and (unknown) control (Eqs. (32)-(35) or (44)-(47), respectively) with the elements of the density matrix undergoing a unitary evolution, we can, in principle, determine the control that eliminates the effect of the decoherence and momentarily restores unitarity. In our case, this simply means that the matrix elements should be restored to their initial values, i.e. the unitary evolution is the identity operator. Indeed, since the primary goal here is to show robust elimination of decoherence, we consider that the NOT transformation has already been effected and the state is now waiting to be involved in the next operation, as specified by the ongoing quantum algorithm. Since this “waiting time” may be much longer that the time needed to realize the “$`\pi `$-pulse” of the NOT gate, we simply ignore the latter and concentrate on the former. Inclusion of the “$`\pi `$-pulse” due to a NOT gate has been actually implemented and does not lead to any significant change other than an unnecessay complication of the formalism. This separation of the quantum gate operations from the decoherence control operations allows us to treat any one-qubit gate.

To numerically implement the control strategy, the real and imaginary parts of the four complex elements of the density matrix calculated in Section III are used in a prescribed order (see below). While this particular order is not essential, we shall define and use it consistently here, to make it easier to follow our control strategy. Indeed, since a single control cannot realize the desired effect (namely the instantaneous restoration of the ideal behavior) for all eight elements at once, we have to adjust these elements in turn. The control cycle is made up of 8 control steps, since there are 8 equations to solve for the (real and imaginary) components of the $`2\times `$2 complex density matrix. This leads to a sequence of eight real transcendental equations to be solved in turn.

Ideally, the density matrix of a $`2\times 2`$ system is usually completely described by three independent quantities. It may happen that, even in the non-ideal case, for particular Hamiltonians, certain symmetries impose certain relations between the matrix elements. However, to maintain a general and systematic character of the approach, we prefer to treat them all as independent variables, at the expense of a slight lengthening of the control process. Since this process eventually stabilizes into a rather short cycle, we do not consider this to be a serious drawback.

The algorithm is applied identically for either the adiabatic or the thermal case; thus, we describe it for a generic decoherence function, denoted $`g`$. In the first cycle of eight time steps, we start by considering $`\rho _{11R}(1):=\rho _{11R}(g(1),I(1))`$ after one time step. To determine the control pulse required to set $`\rho _{11R}(1)`$ equal to its unitary value after applying a NOT gate, we find a control value $`I(1)`$ which solves the equation $`\rho _{11R}(1)=\rho _{11R}(0)`$.

Finding the appropriate solutions of the above transcendental equation requires a certain care. Indeed, if there are no solutions in a given interval, it is impossible to restore the unitary behavior. On the other hand, if there are multiple solutions, it becomes very difficult to guarantee that the same control is used for every cycle. To obtain a unique and consistent solution, we had to choose carefully a specific interval. Of course, this caveat applies for the solutions of the other equations as well (see below).

Denoting the solution of this equation by $`I(1)`$, we reset $`\rho _{11R}(1)`$ exactly to its original value by using a control pulse $`I(1)`$. After the first time step, all the other elements will have suffered the effect of decoherence, $`g(1)`$, and the effect of the control pulse, $`I(1)`$:

$`\rho _{11R}(1)`$ $`=`$ $`\rho _{11R}(g(1),I(1))=\rho _{11R}(0)`$ (53)
$`\rho _{11I}(1)`$ $`=`$ $`\rho _{11I}(g(1),I(1))`$ (54)
$`\rho _{12R}(1)`$ $`=`$ $`\rho _{12R}(g(1),I(1))`$ (55)
$`\rho _{12I}(1)`$ $`=`$ $`\rho _{12I}(g(1),I(1))`$ (56)
$`\rho _{21R}(1)`$ $`=`$ $`\rho _{21R}(g(1),I(1))`$ (57)
$`\rho _{21I}(1)`$ $`=`$ $`\rho _{21I}(g(1),I(1))`$ (58)
$`\rho _{22R}(1)`$ $`=`$ $`\rho _{22R}(g(1),I(1))`$ (59)
$`\rho _{22I}(1)`$ $`=`$ $`\rho _{22I}(g(1),I(1))`$ (60)
By definition, the restoration of unitarity executed by the control $`I`$ is exact.

The control strategy advances through the next seven time steps in a similar fashion. The required control pulses are determined by solving each of the other seven density matrix equations in the order presented above.

For example, at the second time step, we consider the equation for $`\rho _{11I}(2)=\rho _{11I}(g(2),I(1)+I(2))`$. To determine the control pulse required to set $`\rho _{11I}(2)`$ equal to its desired value, we find the control $`I(2)`$ which solves the equation $`\rho _{11I}(2)=\rho _{11I}(0)`$. After applying this control, the values of the matrix elements after the second time step are given by:

$`\rho _{11R}(2)`$ $`=`$ $`\rho _{11R}(g(1),I(2))`$ (61)
$`\rho _{11I}(2)`$ $`=`$ $`\rho _{11I}(g(2),I(1)+I(2))=\rho _{11I}(0)`$ (62)
$`\rho _{12R}(2)`$ $`=`$ $`\rho _{12R}(g(2),I(1)+I(2))`$ (63)
$`\rho _{12I}(2)`$ $`=`$ $`\rho _{12I}(g(2),I(1)+I(2))`$ (64)
$`\rho _{21R}(2)`$ $`=`$ $`\rho _{21R}(g(2),I(1)+I(2))`$ (65)
$`\rho _{21I}(2)`$ $`=`$ $`\rho _{21I}(g(2),I(1)+I(2))`$ (66)
$`\rho _{22R}(2)`$ $`=`$ $`\rho _{22R}(g(2),I(1)+I(2))`$ (67)
$`\rho _{22I}(2)`$ $`=`$ $`\rho _{22I}(g(2),I(1)+I(2))`$ (68)
Note that the decoherence function used to calculate $`\rho _{11R}(2)`$ is different from those for all the other matrix elements, since $`\rho _{11R}(1)`$ was reset to $`\rho _{11R}(0)`$ after the first time step. In general, resetting any one of the matrix elements to its corresponding original value implies that the cumulative effect of the control pulses applied since the last correction (eight time steps in the past) cancels out the effect of decoherence since the last correction.

To complete the first cycle we proceed to solve the equations for the next six matrix elements, to determine the required control pulses. By the end of the first cycle of eight time steps each of the matrix elements depends on different decoherence and control values:

$`\rho _{11R}(8)`$ $`=`$ $`\rho _{11R}(g(7),I(2)+I(3)+I(4)+I(5)+I(6)+I(7)+I(8))`$ (69)
$`\rho _{11I}(8)`$ $`=`$ $`\rho _{11I}(g(6),I(3)+I(4)+I(5)+I(6)+I(7)+I(8))`$ (70)
$`\rho _{12R}(8)`$ $`=`$ $`\rho _{12R}(g(5),I(4)+I(5)+I(6)+I(7)+I(8))`$ (71)
$`\rho _{12I}(8)`$ $`=`$ $`\rho _{12I}(g(4),I(5)+I(6)+I(7)+I(8))`$ (72)
$`\rho _{21R}(8)`$ $`=`$ $`\rho _{21R}(g(3),I(6)+I(7)+I(8))`$ (73)
$`\rho _{21I}(8)`$ $`=`$ $`\rho _{21I}(g(2),I(7)+I(8))`$ (74)
$`\rho _{22R}(8)`$ $`=`$ $`\rho _{22R}(g(1),I(8))`$ (75)
$`\rho _{22I}(8)`$ $`=`$ $`\rho _{22I}(g(8),I(1)+I(2)+I(3)+I(4)+I(5)+I(6)+I(7)+I(8))=\rho _{22I}(0)`$ (76)
a pattern which will recur in every cycle of eight time steps.

After one time step in the second cycle of eight time steps, we are back at the starting point. To determine the control pulse required to set $`\rho _{11R}(9)`$ equal to its initial value, we find a control value $`I(9)`$ which solves the equation $`\rho _{11R}(g(8),I(2)+I(3)+I(4)+I(5)+I(6)+I(7)+I(8)+I(9)))=\rho _{11R}(g(1),I(1))=\rho _{11R}(0)`$.

In the following time step, set $`\rho _{11I}(10)=\rho _{11I}(g(8),I(3)+I(4)+I(5)+I(6)+I(7)+I(8)+I(9)+I(10))`$, and so on.

This procedure is repeated as long as the quantum state of the two-level system has to be maintained, i.e. the “waiting time” mentioned previously until the next logic gate is applied. We note two important things for the applicability of the scheme. First, the knowledge of the decoherence function is needed only for a finite period of time (in the example above only eight time steps). Second, after initial transients, controls will stabilize and the whole control cycle will repeat itself periodically. In other words, for each sequence the controls can be calculated off-line and applied in the required order. This behavior has indeed been observed: the values of the control pulses were rapidly stabilized after the first cycle of 8 time steps, as shown in Figures 1 and 2. The stabilization time depends on initial conditions, and the magnitude of the decoherence.

## V Discussion of the Results

The graphs shown in Figs. 1 and 2 illustrate the typical evolution of the matrix elements while applying this control strategy in both adiabatic and thermal decoherence situations.

As expected, the relative size of the time step between control pulses determines the amount of deviation of the matrix elements from their unitary values. Also, the size and frequency of the control needed to restore perfectly the ideal situation depends on the strength of the decoherence. For $`\gamma =1`$ (extremely strong decoherence), the frequency of the control can be decreased only to about twice the Rabi frequency. If we decrease the frequency even more, we cannot restore exact unitarity, at least not by using this scheme. Since the control algorithm described here is based on periodic perfect restoration of the ideal situation, we shall not discuss this case any further. However, imperfect restoration of unitarity is likely to happen and is $`\mathrm{𝑛𝑜𝑡}`$ at all hopeless: this situation will be analyzed in future work. As the strength of the decoherence, $`\gamma `$, decreases (strong to medium decoherence), we can also decrease the frequency of the control pulses while keeping the same amount of deviation of the matrix elements from their unitary values. Of course, higher frequencies result in better restorations of unitarity.

The same observations apply for thermal decoherence. Since the latter affects all elements of the density matrix, it is reasonable to expect that the same amount of control will be less efficient here than in the adiabatic situation. This is indeed the case, as illustrated in Fig. 2.

If the controls are perturbed by noise, the control is still very efective. It is important to note that the amplitude of the noise is calibrated with respect to the state and not with respect to the amplitude of the control needed in the noiseless situation, which in some cases is extremely small.

In conclusion, we presented an open loop control scheme for quantum gates which is: (i) realistic - the control and decoherence act simultaneously, and the model is able to deal with adiabatic as well as thermal decoherence; (ii) relatively simple to calculate, (iii) efficient to implement, and (iv) robust to noise both in the decoherence functions as well as in the control pulses.

The results presented in this paper depend, among other factors, on the model Hamiltonian, and the approximation used in the tracing out procedure. However, modifications of either would not alter the main features of the proposed control strategy.

While our control scheme shares several common features with the quantum feedback scheme proposed by Wang et al. Wang01  and the quantum bang-bang control, there are also significant differences between our scheme and these approaches, which we briefly discuss here.

The quantum feedback eliminates specifically the decoherence effects of spontaneous emission, whereas the present scheme is more general. Both the feedback scheme and the approach proposed here aim to restore (at least approximately) the unitarity i.e., maintain the information contained in the quantum state. However our scheme assumes prior knowledge of the initial state in order to maintain it unchanged, whereas the quantum feedback scheme drives the state to a target state regardless of initial conditions. Finally, our scheme does not require any measurements on the two-level system or the environment, and does not involve any feedback.

Similarly to quantum bang-bang control schemes Tombesi\_ol ; Viola98 ; Viola99a ; Viola99b , the control proposed here relies on a series of relatively fast control pulses to dynamically cancel the system-environment interaction. However, in our approach, we tailor the control according to the specific nature of decoherence present in a physical system. This has certain advantages, namely: (i) controls are applied only at the required frequency and magnitude, without overtaxing the available resources; (ii) the limits of applicability are also clear, beyond which a specific open loop scheme would not be able to restore exact unitarity; (iii) since the environment is not supposed to change significantly during the quantum computation, the control sequence stabilizes very quickly to a set of periodic controls, as shown in Figs. 1 and 2 (this implies that in most situations of interest the controls can be computed off-line and then simply and efficiently applied in the prescribed sequence); and (iv) robustness to noise appears quite strong, both in the adiabatic and thermal cases.

The main difference between our scheme and the bang-bang control scheme proposed by Viola and Lloyd Viola98  is that, on the one hand we can handle the more difficult case involving thermal decoherence, and on the other hand the control pulses are tailored to the task instead of being generic. Indeed, a cycle of two predetermined $`\pi /2`$-pulses is used in Viola and Lloyd’s scheme for adiabatic decoherence. The evolution of these pulses is disjoint from the evolution of the system and the environment, in order to ensure that the cycle evolution can be simplified to an expression containing only $`\sigma _z`$. However, the control proposed in Ref. Viola98  cannot be used as such in the case of thermal decoherence, which has a system-bath interaction Hamiltonian proportional to $`\sigma _x`$, even if the control Hamiltonian is taken to be a linear combination of $`\sigma _x`$ and $`\sigma _y`$. On the contrary, in our scheme the evolution of the pulses is not disjoint from the evolution of the system and the environment, which allows us to handle a general combination of thermal and adiabatic decoherence by employing the BCH approximation. Moreover our controls do not have a priori assigned values. In general, the stabilized pulses have a much smaller amplitude than the $`\pi /2`$-pulses, for the same pulse duration, and they are applied at a slower rate. In principle, this improves the efficiency of the control.

The limitation of our scheme is its reliance on knowledge of the target state that we want to maintain, e.g., the initial values of the reduced density matrix of the two-level system. Since the target state changes during the course of the quantum computation, the applicability of this control scheme is presently limited to short calculations. We intend to explore this aspect together with different Hamiltonians, approximations, and control schemes in the future.

## Acknowledgments

This research was supported in part by the U.S. Department of Energy, Office of Basic Energy Sciences. The Oak Ridge National Laboratory is managed for the U.S. DOE by UT-Battelle, LLC, under contract No. DE-AC05-00OR22725. The authors thank P. Tombesi for pointing out several important references. We are particularly indebted to L. Viola for several pertinent and insightful comments, suggestions, and questions that led to the clarification and improvement of this paper.


format: text

text:
# Translocation Properties of Primitive Molecular Machines and Their Relevance to the Structure of the Genetic Code

## 1 Introduction

The origin of the genetic code has been a subject of intense research since its structure was completely elucidated in the early 1970’s. In subsequent years, the scientific community has produced several theories in order to explain why the genetic code has this structure. Among these theories, the most prominent ones are the stereochemical theory, the frozen-accident theory and the coevolutionary theory -. Roughly speaking, these theories try to account for the structure of the genetic code by looking at the interactions between codons and amino acids, the biosynthetic relationships among different amino acids and how the metabolic pathways between them have been selected throughout evolution. Nevertheless, the fact that all the codons are made up of three nucleotides, has mostly been taken for granted and barely brought into question.

One of the most widely used arguments found in the literature to explain the trinucleotide codon structure of the genetic code, was given by Sidney Brenner in 1961 . According to this argument, codons are made up of three nucleotides (or bases, for short) because there are 20 amino acids to be specified by the genetic information expressed by a 4 letter “alphabet” (the four bases A,G,C,U). If codons were composed of only two bases, there would be only 16 different combinations ($`4^2`$), which are not enough to specify for 20 amino acids. If instead, codons were made up of more than three bases, there would be at least 256 combinations ($`4^4`$), and these are too many for only 20 amino acids. Hence, less than three bases per codon are not enough, and more than three would imply an excessive degeneration of the code. The result coming out from this argument is that three bases per codon is the optimal “bit of information” that can be used in order to specify for the 20 different amino acids by means of a 4 letter “alphabet”.

The above argument, however, does not constitute an explanation by itself, mainly because it only moves the question of “why three?” to the questions of “why twenty?” or for that matter “why four?”. There is no reason for the genetic information to codify for only 20 amino acids since living organisms use more than those specified by the genetic code . In addition, this argument assumes that all the codons must have the same length (number of bases), even though more efficient codes can be obtained by allowing the length of the codons to vary . Finally, *given* that 20 amino acids have to be specified by using 4 different bases, Brenner’s argument leads to the simplest code that might be thought of. But even in such a case, simplicity has to be accounted for as a relevant criterium.

In this work we address the question of the origin of the three-base codon structure of the genetic code from a dynamical point of view. We consider a simple molecular machine model which captures some of the principal features of the interaction between primitive realizations of the ribosome and of the mRNA. Our main objective is to present a dynamical scenario, compatible with prebiotic conditions, of how the triplet structure of the genetic code could have arisen. The model we propose is a follow up of the one introduced by Aldana, Cocho and Martínez-Mekler  and is consistent with the current evidence suggesting the “RNA world” hypothesi1s . In this scheme the crucial molecules involved in the prebiotic and protobiotic processes, that eventually led to codification and translation mechanisms of the genetic information, were RNA related.q

In our model, based on the setup depicted in Fig. 1, a short one-dimensional polymer composed of $`M`$ monomers interacts with a much longer one, via electrostatic forces. In order to avoid confusion, from now on we will refer to the short polymer as “the chain”, and to the long polymer simply as “the polymer”. The electrostatic interaction between the chain and the polymer is due to the presence of electric charges, or multipolar moments, in the monomers of both the chain and the polymer.

The charges of the monomers of the chain and of the polymer are assigned at random following a uniform distribution. Therefore, the resulting chain-polymer interaction potential has a random profile. The chain is allowed to move along the polymer, but is constrained to remain at a fixed perpendicular distance $`\sigma `$ from it. Consequently, transport is one-dimensional. One of our main results is to show that under very general conditions, a dynamics is attained in which the chain moves along the polymer in effective “steps” whose mean length is three monomers. We argue that this dynamical feature may be one of the underlying causes of the three base codon structure of the genetic code.

This paper is organized as follows: section 2 describes in detail the model and the assumptions introduced. In section 3 we recall some statistical aspects of our previous analysis  of the random interaction potentials between the chain and the polymer for the simplest case in which the former is composed of just one particle ($`M=1`$). We exhibit numerically that, even in this simple case, the mean distance $`\overline{d}`$ between consecutive minima along the interaction potential is very close to three: $`\overline{d}3`$ (taking the monomer length as spatial unit). After retrieving the analytical expression for this distance , we then look into the probability $`𝐏_m(d)`$ of two neighboring potential minima being separated by a distance $`d`$. Subindex $`m`$ refers to the number of different types of monomers in the polymer and in the chain. This probability function shows that, even though the mean distance $`\overline{d}`$ is close to three, the most probable distance between consecutive minima is $`d^{}=2`$ for $`m>2`$. In section 4 the monomer charges along the polymer are assigned in correspondence with protein-coding regions of the genome of real organisms (e.g. *Drosophila* or *E. coli*) instead of at random. For this case, the probability function $`𝐏_m(d)`$ is modified so that not only the mean distance is $`\overline{d}3`$, but also the most probable one happens to be $`d^{}=3`$.

In section 5 we introduce the more realistic case $`M>1`$, which takes into account the fact that the ribosome is not a point particle, that it has spatial structure and presents several simultaneous contact points between its own rRNA and the mRNA polymer. For small chain lengths ($`M10`$), the probability distribution $`𝐏_m(d)`$ is indicative of wide fluctuations and has a form strongly dependent on the particular assignment of charges in the chain. One of our main findings is that for such chains the most likely configurations are those in which both, the mean distance $`\overline{d}`$ and the most probable one $`d^{}`$ are equal to three, *even when the monomer charges along the polymer and the chain are assigned at random*. In section 6 we analyze the dynamics resulting from the model when an external force is pulling the chain, forcing it to move as a rigid object along the polymer. The power spectrum of the velocity of the chain reveals that, under some very general circumstances, for small chain lengths ($`M10`$), there is a sharp periodicity in the dynamics of the system, with a slowing down of the velocity of the chain every three monomers. Finally, section 7 is devoted to the discussion of the results and their relevance to the origin of the genetic code.

## 2 The model

The model we propose consists of a chain of $`M`$ monomers interacting with a very long polymer composed of $`N`$ monomers, with $`N\mathrm{}`$ (see Fig.1). The chain is constrained to remain at a given distance $`\sigma `$ perpendicular to the polymer and is allowed to move in along the polymer, we shall define $`x`$ as its position in this direction relative to the polymer. We will denote the monomer charges in the chain and in the polymer by $`\{p_i\}`$ and $`\{q_j\}`$, respectively. We should mention that by “charge” we do not necessarily mean Coulomb charge. Both $`\{p_i\}`$ and $`\{q_j\}`$ could be dipolar moments, induced polarizabilities, or similar quantities resulting from electrostatic interactions between chain monomers and polymer monomers with potentials of the form $`1/r^\alpha `$, where $`\alpha `$ characterizes the “charge” type.

We will assume that all the monomers in the chain, and separately, all the monomers in the polymer, are of the same nature, namely, all of them are either Coulomb charges, or dipolar moments, or polarizable molecules, etc. In addition, taking into account that in the origin of life conditions the genetic molecules were not yet likely to convey any structured information, we will consider the charges $`\{p_i\}`$ and $`\{q_j\}`$ to be discrete independent random variables, acquiring one of the $`m`$ different values $`\xi _1,\xi _2,\mathrm{},\xi _m`$ with the same probability. Hence, the probability function $`P(q)`$ for both, the $`\{p_i\}`$ and $`\{q_j\}`$ variables, will be

$$P(q)=\frac{1}{m}\underset{j=1}{\overset{m}{}}\delta (q\xi _j)$$
(1)
where $`\delta (q)`$ is the Dirac delta function. In general, in this work we will take the values $`\xi _1,\xi _2,\mathrm{},\xi _m`$ as integers. Parameter $`m`$ represents the number of different types of monomers from which the polymer and the chain are made of. For the case of real genetic sequences $`m=4`$, but we will not restrict the value of $`m`$ to be 4.

All the monomers will have the same length $`L`$, which we take as the spatial unit of measure: $`L=1`$. We also assume the charge in each of the monomers to be uniformly distributed along the length $`L`$, so that the charge density $`\lambda _j(x)`$ in the jth-monomer of the polymer, for example, is a constant whose value is $`\lambda _j(x)=q_j/L`$. Nevertheless, it is worth mentioning that the dynamics of the model does not depend strongly on the particular shape of the monomer charge density $`\lambda _j(x)`$, as long it is a smooth function of $`x`$ (“smooth” in the sense of differentiability).

With the preceding assumptions, the interaction potential $`V_{ij}(x)`$ between the ith-monomer in the chain and the jth-monomer in the polymer is given by

$$V_{ij}(x)=Kp_iq_j_{x+i1}^{x+i}_{j1}^j\frac{dx^{}dx^{\prime \prime }}{\left[(x^{}x^{\prime \prime })^2+\sigma ^2\right]^{\alpha /2}}$$
(2)
where $`K`$ is a constant whose value depends on the unit system used to measure the physical quantities. In the above expression, $`L`$ has already been set equal to $`1`$. Parameter $`\alpha `$ characterizes the kind of interaction between the chain and the polymer: $`\alpha =1`$ corresponds to an ion-ion interaction, $`\alpha =2`$ represents an ion-dipole interaction, and so forth. Note that this parameter does not depend on the indices $`i`$ and $`j`$, since all the monomers in the chain are of the same nature and those in the polymer are themselves of the same nature, differing from each other only in the value of the charge they contain. The overall interaction potential $`V_\sigma ^\alpha (x)`$ between the whole chain and the entire polymer is given by the superposition of the individual potentials $`V_{ij}(x)`$:

$$V_\sigma ^\alpha (x)=\underset{i=1}{\overset{M}{}}\underset{j=1}{\overset{N}{}}V_{ij}(x)$$
(3)
Equations (2) and (3) establish the type of random potentials we will be considering. Our first aim is to analyze the spatial structure of these potentials, giving their statistical characterization. This will be done in the three following sections. Subsequently, we will consider the dynamics of the chain moving along the polymer interacting with it by means of a random potential, subject to an external driving force and seek under what conditions, if any, transport in “steps” of three monomers can be achieved.

## 3 Random potentials: $`M=1`$

Let us start with the simplest case $`M=1`$, in which the chain consists of just one monomer. We will refer to this case as the “single-monomer-chain” case, and to the chain simply as “the particle”. The reason to consider this simple situation is twofold: on one hand, it is useful in order to introduce the relevant ideas behind the model. On the other hand, it is simple enough as to obtain exact analytical results in a more o less straightforward way. In previous work we have already analyzed some statistical properties of the random potentials given by expressions (2) and (3) for the case $`M=1`$ . After a short review of some of those results we center our attention on the probability distribution $`𝐏_m(d)`$.

The overall particle-polymer interaction potential is given by

$$V_\sigma ^\alpha (x)=\underset{j=1}{\overset{N}{}}q_j_x^{x+1}_{j1}^j\frac{dx^{}dx^{\prime \prime }}{\left[(x^{}x^{\prime \prime })^2+\sigma ^2\right]^{\alpha /2}}$$
(4)
Note that in the previous expression we have set $`Kp_1`$, the only charge in the chain, equal to one. Fig.2 shows three graphs of the potential $`V_\sigma ^\alpha (x)`$ for $`\sigma =0.5`$ and different values of the parameter $`\alpha `$. To generate these graphs, the following probability function for the charges $`\{q_j\}`$ was used:

$$P(q)=\frac{1}{6}\underset{\mu =3}{\overset{3}{}}\delta (q\mu )$$
(5)
namely, each one of the variables $`\{q_j\}`$ acquired one of the six different values $`\{\pm 1,\pm 2,\pm 3\}`$ with probability $`1/6`$ ($`m=6`$). In Fig.2, the random realization of the charges $`\{q_j\}`$ along the polymer was the same for the three graphs. As can be seen from this figure, the distribution of maxima and minima along the potential does not change by varying the value of the parameter $`\alpha `$, in the sense that all the maxima and minima remain essentially at the same positions. What occurs as $`\alpha `$ takes larger and larger values is that the potential becomes a step-like function.

Fig.3 presents an analogous situation, but now keeping $`\alpha `$ constant ($`\alpha =2`$) and varying $`\sigma `$. The behavior of the potential is similar to the previous case: the potential becomes a step-like function as $`\sigma `$ decreases and the positions of the maxima and minima are not appreciably modified.

The above considerations exhibit that for small values of $`\sigma `$, say $`0<\sigma 1`$, the distribution of maxima and minima along the potential is entirely determined by the distribution of charges along the polymer and is independent of the particular values acquired by the parameters $`\alpha `$ and $`\sigma `$. Therefore, in order to find out the distribution of maxima and minima along the interaction potential, it is possible to substitute the continuous random potential given by expressions (2) and (3), by the equivalent step-like potential defined by

$$V(x)=\underset{j=1}{\overset{N}{}}V_j\left[H\left(x(j1)\right)H(xj)\right]$$
(6)
where $`H(x)`$ is the Heaviside function111The Heaviside function $`H(x)`$ is defined as $`H(x)=0`$ if $`x<0`$ and $`H(x)=1`$ if $`x0`$, and $`V_j`$ is a random variable whose value is directly proportional to the charge $`q_j`$ of the jth-monomer in the polymer:

$$V_j=Kp_1q_j$$
(7)
where $`p_1`$ is the charge of the particle (the only monomer in the chain).

Since the random charges $`\{q_j\}`$ are statistically independent, so are the $`\{V_j\}`$. Expression (6), which we will refer to as the step-like limit, is suitable for the analytical determination of the probability function of the distances between consecutive potential minima. This probability function gives important information concerning the dynamics of the system. If some external force is acting on the particle (or the chain), forcing it to move in one direction (right or left) along the polymer, the particle will spend more time in the energy minima than in the maxima. Such a movement may be interpreted by considering the particle as “jumping” from one minimum to the next (see Fig.4). It is worth noticing that the mean distance between consecutive minima in the potentials shown in Fig.2 and Fig.3 is nearly three. In Fig.2 there are 33 minima distributed among 100 monomers, and consequently the mean distance between consecutive minima in this case is $`\overline{d}=100/333.03`$. Analogously, the mean distance between neighboring minima in Fig.3 is $`\overline{d}=100/342.94`$. Therefore, it is expected that in its motion along the polymer, the velocity of the particle will slow down, on average, every three monomers, being momentarily “trapped” in each of the potential minima.

By using the step-like limit, in reference  we have shown that the mean distance $`\overline{d}`$ between consecutive potential minima for a long polymer (the large N limit) is given by

$$\overline{d}=\frac{6m}{2m1}$$
(8)
where $`m`$ is the number of different monomer types. The above equation shows that the mean distance $`\overline{d}`$ is always between 3 and 4, and approaches 3 asymptotically as $`m\mathrm{}`$. In particular, for $`m=4`$ (the biological value) we have $`\overline{d}3.43`$.

In order to characterize the fluctuations around the mean distance $`\overline{d}`$, it is useful to compute the probability distribution function $`𝐏_m(d)`$, which we recall, gives the probability of two consecutive minima being separated by a distance $`d`$ when there are $`m`$ different types of monomers. In the step-like limit, this computation is carried out by counting all the configurations of the step-like variables $`\{V_j\}`$ in which there are two minima, one at $`V_j`$ and the other at $`V_{j+d}`$, with no other minima in between. The situation is illustrated in Fig.5a. Since in the step-like limit the interaction potential is constant along every monomer, we will adopt the convention to measure the distance $`d`$ between two adjacent minima from the mid point of the first minimum to the mid point of the second one, as illustrated in Fig.5b. With the above convention, the resulting distances $`d`$ can only acquire integer or half-integer values. For a finite number $`m`$ of different charges, the explicit calculation of the probability function $`𝐏_m(d)`$ consists mainly on counting configurations, though conceptually straightforward, it involves a considerable amount of algebra. Here we present the final expressions:

$$𝐏_m(d)=\frac{6}{2m1}\underset{k=0}{\overset{d2}{}}\frac{2k+1}{m^{k+d+2}}N_m(dk)$$
(9)
if $`d`$ is integer: $`d=2,3,4,\mathrm{}`$, and

$$𝐏_m(d)=\frac{12}{2m1}\underset{k=0}{\overset{d^{}2}{}}\frac{2k+1}{m^{k+d^{}+3}}N_m(d^{}k)$$
(10)
if $`d`$ is half-integer: $`d=5/2,d=7/2,d=9/2,\mathrm{}`$, where $`d^{}=Int[d]`$. In the above expressions, $`N_m(d)`$ is a polynomial whose degree and coefficients depend on $`m`$. For $`m=2`$, $`m=3`$ and $`m=4`$ the polynomials are given by

$`N_2(d)`$ $`=`$ $`1`$
$`N_3(d)`$ $`=`$ $`1+2d+2d^2`$ (11)
$`N_4(d)`$ $`=`$ $`2+{\displaystyle \frac{1}{8}}\left[30d+73d^2+22d^3+3d^4\right]`$
For the case of $`m=\mathrm{}`$ we have also derived a closed expression , which has a much simpler form:

$$𝐏_{\mathrm{}}(d)=3\frac{2^d}{(d+3)!}\left(d^2+d2\right)$$
(12)
The preceding distributions are plotted in Fig.6. It can be seen from this figure that the most probable distance $`d^{}`$ between consecutive potential minima is $`d^{}=2`$, except for the case $`m=2`$ in which $`d^{}=3`$. Hence, according to the transport mechanism suggested in Fig.4, whenever there are more than two different types of monomers, the particle will move along the polymer in “jumps” whose mean length is close to three, but whose most probable length is actually two. The difference between the mean distance $`\overline{d}`$ and the most probable distance $`d^{}`$ is due to the presence of “tails” in the probability function $`𝐏_m(d)`$. Namely, to the fact that $`𝐏_m(d)`$ has non zero values even for large $`d`$. Nevertheless, in section 5 we will show that these “tails” can be shrunk almost to zero when the chain is made up of more than one particle ($`M>1`$). This is one of the main results of this paper.

To end this section, it is worth mentioning that half-integer distances between two neighboring minima occur when one or both of these minima extend over several monomers (see Fig.5b). In these configurations, the charges of the adjacent monomers constituting the extended minimum have the same value. Configurations in which groups of adjacent equally charged monomers occur, are less likely than configurations in which adjacent monomers have different charges, and the former tend to disappear as $`m`$ increases (see Fig.6).

## 4 Real genetic sequences: $`M=1`$

The charges $`\{q_j\}`$ along the polymer can be assigned in correspondence with the genetic sequence of an organism, rather than in a random way. The purpose of doing so is to find out how the potential minima and maxima along real genetic sequences are distributed, and to compare the resulting distribution with the one corresponding to the random case. Since genetic sequences are made out of four different bases (A, U, C and G) we consider four different possible values $`\xi _1,\xi _2,\xi _3,\xi _4`$ for the charges $`\{q_j\}`$, i.e. $`m=4`$. To proceed further, it is necessary to establish a correspondence between the charge values $`\xi _1,\mathrm{},\xi _4`$ and the four bases A, U, C and G. An arbitrary look up table is the following:

$$\begin{array}{cc}\text{base}& \text{charge value}\\ & \\ \text{A}& 2\\ \text{U}& 1\\ \text{C}& +1\\ \text{G}& +2\end{array}$$
(13)
With the above correspondence, if the jth-base in a given genetic sequence happens to be A, for example, then the charge of the corresponding jth-monomer in the polymer will be $`q_j=2`$.

Figure Fig.7a shows the probability distribution $`𝐏_4(d)`$ computed numerically by using a *Drosophila melanogaster* protein-coding sequence, 45500 bases in length (several genes were concatenated to construct this sequence). The mean distance between consecutive potential minima for this sequence is $`\overline{d}3.15`$ and, as follows from the figure, the most probable distance is $`d^{}=3`$. Therefore, in the “real sequence case” not only is the mean distance $`\overline{d}`$ very close to 3, but also the most probable distance $`d^{}`$ turns out to be 3. A comparison of Fig.6c with Fig.7a, shows that for protein-coding genetic sequences, the potential minima along the polymer are more often separated by three monomers than in the random case. When protein-coding sequences are used, the value of $`𝐏_4(d)`$ increases at $`d=3`$ and decreases at $`d=2`$.

The above behavior does not occur when non-coding sequences of the genome are used for the monomer charge assignment along the polymer. Fig.7b shows the probability function $`𝐏_4(d)`$ for the case in which the monomer charges are in correspondence with an intergenic sequence of the *Drosophila’s* genome. The length of the sequence is again 45500 bases. As can be seen from the figure, in this case the probability function $`𝐏_4(d)`$ looks much more like the one obtained in the random case.

It is important to remark that the behavior of $`𝐏_4(d)`$ exhibited in Fig.7a for real protein-coding sequences does not depend on the particular correspondence (13) between bases and charge values being used, as long as they are of similar order of magnitude and they can allow for an order relation. These conditions hold for the four bases A, T, C and G which have charges of the same order of magnitude 222The real physical charge values will depend on the system of units, which is contained basically in the constant $`K`$ appearing in expression (2)and are fulfilled by our present choice of $`\{\pm 1,\pm 2\}`$. The order relation is necessary for the interaction potential to have maxima and minima. It is worth asking the effect that changes in the order relation have on the probability function $`𝐏_4(d)`$. There are 24 possible order relations among the four bases A, T, C, and G ($`4!`$ permutations). The 24 probability functions $`𝐏_4(d)`$ corresponding to these permutations are plotted in Fig.8a for *Drosophila’s* protein-coding sequence. As can be seen from the figure, the probability functions basically overlap, independently of the particular order relation between the bases, with a peak at $`d^{}=3`$. The invariance of $`𝐏_4(d)`$ under base permutations also holds for non-coding sequences as is shown in Fig. 8.b., where $`d^{}=2`$ for intergenic sequences of Drosophila *Drosophila’s*. This value of $`d^{}`$ suggests that non-coding sequences behave as random structures.

The “peaking at $`d^{}=3`$” of the probability function $`𝐏_4(d)`$ seems to be a general characteristic associated with the protein-coding sequences of living organisms, not only with *Drosophila*. In Fig.9 we show the probability functions obtained from protein-coding sequences of different organisms, and in all the cases the probability functions present their highest value at $`d^{}=3`$ (the mean distance $`\overline{d}`$ is also very close to $`3`$). The fact that the above characteristic is absent in non-coding genetic sequences may be interpreted in evolutionary terms. Genetic sequences directly involved in the protein-translation processes were selected (among other things) as to bring the distance between consecutive potential minima closer to 3, both in mean and frequency of occurrence. This interpretation raises a question: how likely is it to obtain a randomly generated sequence with a structure similar to that of protein-coding sequences? In other words, if we generate a random sequence and compute its probability function $`𝐏_m(d)`$, how likely is it to come up with a probability function peaking at $`d^{}=3`$?.

In order to answer this question, we define the parameters $`p_{3/2}`$ and $`p_{3/4}`$ as

$`p_{3/2}`$ $`=`$ $`{\displaystyle \frac{𝐏_4(3)}{𝐏_4(2)}}`$
$`p_{3/4}`$ $`=`$ $`{\displaystyle \frac{𝐏_4(3)}{𝐏_4(4)}}`$
If the probability function $`𝐏_4(d)`$ associated with a given sequence has its highest value at $`d^{}=3`$, then the corresponding parameters $`p_{3/2}`$ and $`p_{3/4}`$ will both be greater than one. Otherwise, one or both of these parameters will be smaller than one.

Fig.10a is a plot of $`p_{3/2}`$ vs. $`p_{3/4}`$ for 1000 random sequences, each one consisting of 500 bases (which is a typical length of sequences coding functional proteins ). It can be seen from the figure that only a small fraction of the points (about 0.224) fall in region I, for which $`p_{3/2}>1`$ and $`p_{3/4}>1`$. The rest of the points fall in region II, in which $`p_{3/2}<1`$ and $`p_{3/4}>1`$. Therefore, the probability of having a random sequence, 500 bases long, whose consecutive potential minima are more often separated by a distance $`d^{}=3`$, is close to $`0.224`$.

On the other hand, Figs.10b-f show similar graphs, but using protein-coding sequences of real organisms. These graphs were constructed by analyzing short coding sequences 500 bases in length. The fraction of points falling in region I ($`p_{3/2}>1`$ and $`p_{3/4}>1`$) for the different organisms of Fig.10 is summarized in the following table:

| Organism | Frac. region I |
| --- | --- |
| random sequence | 0.224000 |
| *Chlorella vulg.* | 0.788235 |
| *Drosophila* | 0.780220 |
| *Deinococcus rad.* | 0.803525 |
| *Myc. tuberculosis* | 0.802643 |
| *Aeropyrum pernix* | 0.719973 |

These results show that protein-coding sequences are much more likely to have their consecutive potential minima separated by a distance $`d^{}=3`$.

To end this section, we should make a further comment. The fact that the most probable distance $`d^{}`$ between consecutive potential minima is $`d^{}=3`$ for protein-coding sequences, is a consequence of the order in which the bases or codons appear along the sequence and are not related to the relative weight (fraction) of their occurrence. For example, Fig.11a shows the probability function $`𝐏_4(d)`$ corresponding to a protein-coding sequence of *Escherichia coli*, 45000 bases in length. In Fig.11b the codon composition, expressed as fractional occurrence, of this sequence is depicted. To construct this last graph, we labelled each codon in an arbitrary way, by assigning an integer number, from 0 for AAA up to 63 for GGG. On the horizontal axis the number of the codon is plotted, and on the vertical axis its fractional occurrence on the *E. coli* sequence. Finally, Fig.11c shows the probability function $`𝐏_4(d)`$ corresponding to a random sequence, of the same length and codon composition as the one used in Fig.11a. As can be seen, for the “randomized” *E. coli* sequence the most probable distance is no longer $`d^{}=3`$ as in Fig.11a, but $`d^{}=2`$.

## 5 Extended chain: $`M>1`$

In this section we consider the case in which the chain is composed of more than one single monomer. As mentioned before, this reflects the fact that the ribosome is not a point particle after all. At each given time, the ribosome interacts with the mRNA at several points, giving rise to a collective interaction. Electron microscopy has revealed that the mRNA thread passes across the ribosome throughout a *tunnel* about 10nm long and 2nm in diameter . Since nucleotide length is about 0.5nm, around 20 nucleotides may be in a position to interact simultaneously with the ribosome. Taking the above into account, we will work with a small chain, assuming $`M=10`$ as a reasonable length.

In the step-like limit, the interaction potential equation can be written by equation (3) now with the variables $`\{V_j\}`$ given by:

$$V_j=Kq_j\underset{i=1}{\overset{M}{}}p_i$$
(15)
The main modification caused by using equation (15) instead of (7) for the random variables $`\{V_j\}`$ is that they are no longer statistically independent. Due to the collective interaction between the chain and the polymer, these variables are strongly correlated. This makes the problem of finding the probability function $`𝐏_m(d)`$ too difficult for analytical treatment. Therefore, we will present only numerical results.

In our simulations, $`𝐏_m(d)`$ presents an erratic behavior. The shape of this function strongly depends on the particular realization of monomer charges in the chain. For example, Fig.12a shows the probability function $`𝐏_m(d)`$ for a polymer 500 monomers in length and a particular realization of charges in a 10-monomer chain. The charge values we used to generate these graphs were $`\{\pm 1,\pm 2\}`$. Fig.12b shows a similar graph but for a different realization of charges along the chain (the realization of charges in the polymer was the same in both cases). As can be seen, these two graphs differ considerably: the first presents a very sharp maximum at $`d^{}=2`$, whereas the second does so at $`d^{}=3`$. Wide fluctuations in the probability function $`𝐏_m(d)`$ were always present in our numerical simulations.

This fluctuating behavior did not occur in the single-monomer-chain case. It is apparent from Fig.10a that the fluctuations of the probability function $`𝐏_m(d)`$ in that case are rather small: the probability function has its highest value at $`d^{}=2`$ for the majority of the realizations (77.6$`\%`$) in the case $`M=1`$. Also, the fluctuations are concentrated in a small region of the $`(p_{3/2},p_{3/4})`$ plane.

In order to have an idea of the magnitude of the fluctuations of the probability function $`𝐏_m(d)`$ in the extended-chain case, we can make use of the parameters $`p_{3/2}`$ and $`p_{3/4}`$ defined in expression (LABEL:p32). The fluctuations of the probability function in the extenden-chain case are shown in Fig.13. To construct this figure, a polymer 500 monomers in length was used along with a 10-monomer chain. The charge values in the polymer as well as in the chain were again $`\{\pm 1,\pm 2\}`$. The parameters $`p_{3/2}`$ and $`p_{3/4}`$ were calculated for 4500 charge realizations in the chain (the realization of charges in the polymer was the same in all the cases). Fig.13 is the plot of $`p_{3/2}`$ versus $`p_{3/4}`$ for these 4500 realizations. Two remarks are worth mentioning about this figure:

* The points are spread out over a much wider region than in the single-monomer-chain case (Fig.10a).
* The majority of the points fall in the region in which both $`p_{3/2}>1`$ and $`p_{3/4}>1`$.

The following table shows the fractions of the points falling in each of the four regions of the graph:

| Region | Fraction |
| --- | --- |
| $`p_{3/2}>1`$, $`p_{3/4}>1`$ | 0.52978 |
| $`p_{3/2}<1`$, $`p_{3/4}>1`$ | 0.33800 |
| $`p_{3/2}>1`$, $`p_{3/4}<1`$ | 0.07778 |
| $`p_{3/2}<1`$, $`p_{3/4}<1`$ | 0.05444 |

Therefore, when a collective interaction between the polymer and the chain prevails, a remarkable property arises: The probability of having a random interaction potential, whose consecutive minima are more often separated by three monomers, is the largest.

## 6 Dynamics

Recent experimental evidence suggests that the ribosome-mRNA system presents a ratchet-like behavior in the protein synthesis translocation process . In this view, the ribosome is tightly attached to the mRNA thread in the absence of GTP. This is so because the channel in the ribosome through which the mRNA passes, is more or less closed. When a GTP molecule is supplied (and transformed into GDP), this channel opens leaving the mRNA thread free to move one codon. Subsequently the mRNA passage in the ribosome closes again, trapping the mRNA molecule.

In this clamping mechanism several physicochemical factors are involved, which if taken into account in detail would lead to complex dynamical equations hard to handle. In this work our approach is to look into the behavior of oversimplified molecular models which might capture some of the essential dynamical features of the system and may shed some light on how this mechanism could have arisen in the origin of life conditions.

In our modelling the dynamics of the system is governed by the application of an external force $`F_{ex}`$ to the chain in the horizontal direction, i.e. parallel to the polymer. By this means the chain will be forced to move along the polymer. In principle, the force $`F_{ex}`$ may be time dependent, but we will restrict ourselves to a constant term. This force might come from a chemical pump (like GDP) or from any other electromagnetic force present in prebiotic conditions. The only purpose of this force in our model is to drive the chain along the polymer (which is assumed to be fixed, $`N\mathrm{}`$ limit), avoiding it from getting trapped in some of the minima of the polymer-chain interaction potential $`V(x)`$. Therefore, we will also assume that $`F_{ex}`$ satisfies $`|F_{ex}|>`$ max$`|V(x)/x|`$.

Our analysis relies on Newton’s equation of motion in a high friction regime, where inertial effects can be neglected. This regime actually exists in biological molecular ratchets similar to the one we are considering . Under such conditions, the Newton’s equation of motion acquires the form

$$\gamma \frac{dx}{dt}=\frac{V_\sigma ^\alpha (x)}{x}+F_{ex}$$
(16)
where $`\gamma `$ is the friction coefficient. In what follows, we will set $`\gamma =1`$, which is equivalent to setting the measure of the time unit. The above, though a deterministic equation, gives rise to a random dynamics due to the randomness of the interaction potential $`V_\sigma ^\alpha (x)`$. In order to start analyzing this random dynamics, let us consider first the single-monomer-chain case $`M=1`$. In this case, as before, we will refer to the chain simply as “the particle”.

### 6.1 Single-monomer chain: $`M=1`$. Random sequences.

In Fig.14a we show a typical realization of the velocity of the particle $`v(x)`$ as a function of its position $`x`$ along the polymer. This graph was constructed by solving numerically the equation of motion (16), using the fourth order Runge-Kuta method. The parameter values used were $`\alpha =2`$ and $`\sigma =0.5`$, and the monomer charge values were $`\{\pm 1,\pm 2\}`$ (case $`m=4`$). Fig.14b shows the local-transit times of the particle along a short segment of the polymer (40 monomers in length). This transit time is represented in arbitrary units, and was computed by counting how many time steps the particle spent in every spatial interval $`\mathrm{\Delta }x`$ throughout the polymer. In the graph shown, the value of $`\mathrm{\Delta }x`$ was $`\mathrm{\Delta }x=0.25`$. It is apparent from this figure that, in its way along the polymer, the particle spends more time in certain regions than in others, the former being more or less regularly spaced along the polymer.

In order to find out the spatial regularities in the dynamics of the system, it is convenient to take the Fourier transform of the velocity $`v(x)`$ of the particle. Let us call $`\widehat{v}(k)`$ the Fourier transform of $`v(x)`$, $`k`$ being the Fourier variable conjugate to $`x`$.

Fig.15 shows the Fourier power spectrum of the velocity, $`|\widehat{v}(k)|^2`$, for two different realizations of monomer charges in the polymer. The parameter values in Fig.15a and Fig.15b are $`\{\sigma =0.5,\alpha =1\}`$ and $`\{\sigma =0.5,\alpha =4\}`$ respectively. These graphs were computed for the case $`m=4`$, using the charge values $`\{\pm 1,\pm 2\}`$. From the figure, it is evident that there exists a dominant frequency $`k^{}`$ in the power spectrum of the velocity (the highest peak), whose corresponding spatial periodicity is $`\lambda ^{}=2\pi /k^{}3.3`$. The power spectrum reveals a dynamical regularity in the motion of the particle throughout the polymer. This regularity is inherited from the one present in the random potential, in the sense that the particle spends more time in the minima than in the maxima. The consequence is a slowing down of the velocity nearly every three monomers, which is reflected in the power spectrum. Our interpretation is that the peak occurring in the power spectrum of the velocity conveys the information on the average distance $`\overline{d}`$ between consecutive potential minima, which for the case $`m=4`$ is $`\overline{d}3.4`$.

### 6.2 Single-monomer chain: $`M=1`$. Real sequences.

As in section 4, we can assign the charges along the polymer in correspondence with the genetic sequence of real organisms. In order to do that, we will use the same base-charge correspondence as in expression (13). The objective is to find out how the dynamics of the system changes when using real genetic sequences instead of random ones. In what follows, the value of the parameters $`\alpha `$ and $`\sigma `$ will be $`\alpha =1`$, $`\sigma =0.5`$.

The power spectrum of the velocity of the particle throughout the polymer, when using protein-coding sequences of different organisms, is shown in Fig.16. To generate these graphs, short coding sequences of several organisms, each 10000 monomers in length, were used. Two points are worth noticing in this figure. First, the peak in the power spectrum is much higher than in the random case. This indicates that there is a much more well defined periodicity in the dynamics generated by interaction potentials when protein-coding sequences are used. Second, the spatial periodicity reflected in the peak of the power spectrum is much closer to $`3`$ than in the random case.

This dynamical behavior is not present when real but non-coding sequences are used. For example, in Fig.17 we show the power spectrum of the velocity of the particle along the polymer, for two cases in which the monomer charges in the polymer were assigned in correspondence with intergenic regions of two organisms. As can be seen, the structure of such spectra is similar to the one obtained in the random case. In this sense, intergenic regions again seem to have a random structure.

The fact that the power spectrum corresponding to protein-coding sequences exhibits a very sharp periodicity at $`\lambda ^{}3`$, whereas the one corresponding to non-coding sequences does not, has already been reported in the literature . Nonetheless, in these previous works the power spectrum of the “bare” genetic sequences is analyzed, namely, without considering any kind of interaction potential or dynamical behavior. What we have shown here, though, is that this “structural” periodicity around three transforms into a *dynamical periodicity* in the motion of the particle along the polymer.

### 6.3 Extended chain: M=10

The most interesting dynamics occurs when an extended chain is interacting with the polymer. In such a situation, a collective interaction prevails. At every moment there are several contact points between the chain and the polymer. As we have already pointed out, collective interaction between the chain and the polymer gives rise to a widely fluctuating probability function $`𝐏_m(d)`$. The same occurs with the power spectrum of the velocity of the chain along the polymer. However, these fluctuations, far from being annoying, produce a much richer dynamical behavior than in the single-monomer-chain case.

In Fig.18 we show the power spectra of the velocity of the chain along the polymer for two different random realizations of monomer charges in the chain. The charges in the polymer were the same in both cases. These graphs were constructed with a polymer 500 monomers in length and a 10-monomer chain. The parameter values used were $`\sigma =0.5`$ and $`\alpha =1`$. Also, the charge values were, as above $`\{\pm 1,\pm 2\}`$. From the figure, it is apparent that the power spectrum of the velocity exhibits a very well defined dominant frequency, *even though the charges in both the polymer and the chain were assigned at random*. The power spectrum in Fig.18a presents a dominant frequency corresponding to a spatial periodicity $`\lambda ^{}2`$, whereas the corresponding periodicity in the power spectrum shown in Fig.18b is $`\lambda ^{}3`$. We should explain this difference.

The power spectrum appearing in Fig.18a was constructed by using a polymer and a chain whose associated probability function $`𝐏_4(d)`$ has the same shape as the one of Fig.12a. Namely, for this system the probability function has a very high value at $`d^{}=2`$. On the other hand, the power spectrum in Fig.18b corresponds to a system whose probability function has a very sharp peak at $`d^{}=3`$, as the one shown in Fig.12b. From our numerical simulations, we can conclude that *whenever the probability function has a sharp maximum at a distance $`d^{}`$, the power spectrum of the velocity also presents a sharp peak corresponding to a spatial periodicity $`\lambda ^{}=d^{}`$.*

As we have seen, in the collective-interaction case the most probable configurations are those in which the probability function $`𝐏_4(d)`$ has its highest value at $`d^{}=3`$ (see Fig.13). Therefore, if we assign at random the monomer charges in the polymer and in the chain, with high probability we will come up with a dynamics possessing a very well defined periodicity: the chain will move along the polymer in “jumps” whose length is nearly three monomers.

## 7 Summary and discussion

The results presented throughout this work suggest a possible scenario for the origin of the three base codon structure of the genetic code. In this scenario, primitive one dimensional molecular machines, initially with a random structure, exhibited a regular dynamics with a “preference” for a movement in steps of three bases. By “steps” we mean a slowing down of the velocity of the chain along the polymer, nearly every three monomers (see Fig.14b). Even in the simplest case in which the chain consists of only one monomer, the above dynamical regularity is apparent. We can think of the dynamics of primitive molecular machines as being “biased towards three”.

The preceding property is quite robust inasmuch as it hardly depends on the particular kind of interaction between the polymer and the chain. On one hand, the kind of electrostatic potentials we have used is representative of the actual interaction potentials between particles occurring in Nature. These potentials are characterized in our model by the parameter $`\alpha `$. We have also seen that the distribution distances between neighboring maxima and minima along the interaction potential, characterized by the probability function $`𝐏_m(d)`$, does not depend on this parameter (for small values of $`\sigma `$), i.e. it will be the same whether the interaction is coulombian or dipolar or of any other (electrostatic) type.

On the other hand, the spatial distribution of interaction potential minima also does not depend on the particular values of the monomer charges $`\{q_j\}`$ and $`\{p_i\}`$, as long as these values are of the same order of magnitude. An important feature that the charges must comply with is that they take more than two different values. This allows for an order relation to be established among the different types of monomers, leading to a maxima and minima structure of the interaction potential. The probability function $`𝐏_m(d)`$, which gives the probability of two consecutive minima being separated by a distance $`d`$, only depends on $`m`$, namely, on the number of different types of monomers. As $`m`$ increases, the mean distance $`\overline{d}`$ between consecutive potential minima approaches three. Nevertheless, in the particle (single-monomer-chain) case, the most probable distance is $`d^{}=2`$ (for $`m>2`$).

Still in this particle case, considerable changes take place when the charges along the polymer are assigned in correspondence with protein-coding genetic sequences of real organisms. In this case, not only is the mean distance $`\overline{d}`$ between neighboring potential minima nearly three, but also the most probable one, $`d^{}`$, happens to be three. This is a remarkable property of protein-coding sequences, perhaps acquired throughout evolution. Furhtermore the fact that this “refinement” is absent in non-coding sequences of real organisms, strongly suggests that it is a consequence of the dynamical processes involved in the protein synthesis mechanisms.

This interpretation is supported by the results obtained when the *dynamics* of the particle moving along the polymer is considered. In the random sequence case, there are dominant frequencies in the power spectrum of the particle velocity related with the spatial regularities of the interaction potential. Moreover in the protein-coding sequence case, the power spectrum of the velocity shows a very well defined periodicity corresponding almost exactly to a spatial distance $`\lambda ^{}=3`$. Again, this behavior does not occur for non-coding sequences of real organisms, which are not involved in the translation processes.

A richer dynamics emerges when the chain is composed of several monomers. In this, more realistic, collective-interaction case, the probability function $`𝐏_m(d)`$ presents very wide fluctuations, depending on the particular assignment of monomer charges in the chain. Nevertheless, the most probable configurations are those for which the probability function has its highest value at $`d^{}=3`$. For these configurations, the power spectrum of the chain velocity along the polymer exhibits a very well defined spatial periodicity at $`\lambda ^{}3`$.

Our results suggest an origin of life scenario in which primordial molecular machines of chains moving along polymers in quasi one-dimensional geometries, that eventually led to the protein synthesis processes, were biased towards a dynamics favoring the motion in “steps” or “jumps” of three monomers. The higher likelyhood of these primitive“ribosomes” may have led to the present ribosomal dynamics where mRNA moves along rRNA in a channel conformed by the ribosome. Dynamics may have acted in this sense as one of the evolutionary filter favoring the three base codon composition of the genetic code.

Acknowledgements
We would like to thank Leo Kadanoff, Sue Coppersmith, Haim Diamant and Cristian Huepe for very useful discussions and corrections. This work was sponsored by the DGAPA-UNAM project IN103300, the MRSEC Program of the National Science Foundation (NSF) under Award Number DMR 9808595 and by the NSF Program DMR 0094569. M. Aldana also acknowledges CONACyT-México a posdoctoral grant.