text:
# What drives the translocation of stiff chains?

## I Introduction

The transfer of genetic material through the membrane surrounding a cell nucleus is fundamental to the understanding of basic cell processes from gene therapy to viral infection. The motion of linear polymers through pores into confined volumes also arises in many other biological contexts , perhaps the most common examples of which include the translocation of proteins from the cytosol into the endoplasmic reticulum, or into mitochondria or chloroplasts. The export of m-RNA through nuclear pore complexes is still another example of great importance. Despite the longstanding and widespread interest in this process, however, our knowledge about it is still rudimentary.

The process of translocation under influence of an external field or a chemical potential gradient has recently been studied extensively . There have also been several theoretical studies which have specifically investigated chain translocation in the presence of binding particles . The major role of these binding particles has been recognized as a ”Brownian ratchet”, a mechanism which was introduced over ten years ago in pioneering work by Simon, Peskin and Oster . They happen to have treated the case of protein translocation, but their arguments apply equally well to nucleic acids. To account for translocation rates fast compared to simple diffusion under a wide variety of conditions and circumstances, they proposed that non-specific binding by globular proteins results in a ”biased” — or ”ratcheted” — motion of the chain. Many experiments confirm that efficient translocation can indeed take place without the involvement of motor proteins. For example, the entire length (about 40 microns) of the DNA which comprises the genome of T5 phage is observed to enter its bacterial cell host without requiring metabolic energy . The experiment of Salman *et al.* on phage $`\lambda `$ DNA shows similarly that translocation of a comparable length of chain can occur without the help of active processes; simple diffusion would require significantly longer times .

Each of the mechanisms mentioned above for chain translocation, namely diffusion and ratcheting, corresponds to a different time scale and to different physics. Simple diffusion requires a characteristic time $`t_d=L^2/2D`$, where $`L`$ is the total length of the polymer and $`D`$ its diffusion coefficient. In the ”ratcheting” scenario, as soon as a specific length — $`\delta `$ — of the polymer enters, a protein binds to it and the chain is no longer able to diffuse backwards because the pore size is too small for the DNA/protein complex to pass through. In this case, the chain simply diffuses from one binding site to the next, and the translocation time is equal to the product of the time it takes for the chain to diffuse the distance $`\delta `$ times the number of ratcheting sites $`M=L/\delta `$, *i.e.*, $`t_{ratchet}=M\delta ^2/2D=L^2/(2MD)=t_d/M`$, corresponding to a speed-up of the translocation time by a factor of $`M`$ over simple diffusion. As pointed out by Simon, Peskin and Oster, this time represents an idealized limit in which the ratcheting mechanism functions ”perfectly”, i.e., as each successive binding site passes into the cell it is bound irreversibly by a protein that prohibits the chain from diffusing backwards. In actuality, however, an entering site is not necessarily bound immediately by a protein, and/or the protein does not stay adsorbed long enough to act as a ratchet at that site; accordingly, the translocation time is increased beyond $`L\delta /2D`$ by a factor that depends on the ratio of on- and off- rates for binding.

In this paper we consider explicitly the effect of binding particles and show that — via a new mechanism — translocation can occur at rates significantly faster than that provided by the ”perfect” ratcheting scenario described above. More specifically, we argue that the particles which bind reversibly to the chain give rise to a net force on the chain that pulls it into the cell. Furthermore, this force accounts fully for the translocation process and embraces the different mechanisms mentioned above, *e.g.*, in a special limit the Brownian ratcheting appears as a particular idealization (and appealing simplification) of the effect of such a force. The magnitude of the force depends in a delicate way on the concentration of the binding particles and on their diffusion coefficient relative to that of the chain. In the overdamped limit the translocation time in the presence of a force $`F`$ is expected to take the form $`t_F=L/v=L\zeta /F`$, where $`\zeta `$ is the friction coefficient of the chain, related to its diffusion constant through the Einstein relation $`D=k_BT/\zeta `$.

Using Brownian Molecular Dynamics (BMD) simulations, we calculate separately the average force and the translocation time in the presence of binding particles and find that they do indeed obey the relation $`t_F=L\zeta /F`$. Consequently, the translocation process is force-driven and the translocation time turns out to be longer or shorter than the ideal ratcheting time, depending on the concentration and diffusion coefficient of the binding particles. These results can be understood in the context of a generalized Fokker-Planck equation for the probability that at time $`t`$ the chain will have entered the cell to a distance $`x`$ and have $`n`$ particles bound to it. The drift terms corresponding to chain entry and particle binding are obtained from derivatives with respect to entry distance and binding number of a Langmuir adsorption free energy for the overall system. ***In the simplest case of the fully equilibrated, fixed chemical potential Langmuir adsorption problem, the 1D pressure of the system corresponds to the force on the chain: $`P_{1D}=f=(k_BT/\delta )\mathrm{ln}(1+\varphi \mathrm{exp}(\beta ϵ))`$ where $`\varphi `$ is the volume fraction of binding particles, $`ϵ`$ their binding energy, and $`\delta `$ the binding site size. In the limit of large binding energies $`(\beta ϵ1)`$ one recovers the $`fϵ/\delta `$ result discussed in Section IV; otherwise, entropic effects are important as well. We thank Pierre-Gilles de Gennes for these observations. We find that the binding process involves important non-equilibrium effects in general, on which depends the actual value of the force pulling the chain in, and hence the ratio $`t_F/t_{ratchet}`$.

## II Simulation

Brownian Molecular Dynamics simulations of a stiff polymer translocating through a pore and into a cell are performed using a coarse-grained model in which the chain is represented by a rigid rod of beads. The beads are rigidly linked to their nearest neighbors along the chain and do not interact with each other. In this way we model DNA as a perfectly straight — rather than the usual semi-flexible — chain, because the focus of our work is on the entering segment of chain which is within a persistence length from the pore. The link between adjacent beads is rigid in order to avoid chain contraction and extension resulting from the binding of particles; we also neglect changes in the shape or twist of the chain due to binding. The distance between monomers along the chain, $`\sigma `$, corresponds to the “footprint” of binding particles in that the center of each bead is considered as an absorbing site. The binding particles are modeled as spherical, interacting with each other through the repulsive part of a Lennard Jones (LJ) potential with diameter $`\sigma `$; the interaction between binding particles and chain monomers is treated by a *full* (12-6) LJ potential. Since the distance between the absorbing sites is equal to the diameter of the binding particles, each site can be surrounded by a maximum of six particles. The cell-particle interactions are taken to vanish for particles within the radius $`R_s`$ of cell and to increase as $`(RR_s)^4`$ for particles at distances greater than $`R_s`$ from the sphere center. This potential is simply a convenient way to describe an interior wall. Finally, we treat the pore itself as being completely ”inert”, having no effect on the chain except to allow it to enter or leave the cell.

We focus on the dynamics of the chain once one end has been inserted. Let $`x`$ denote the length of chain inside the cell and $`𝐫_𝐢`$ the position of the $`i`$ th binding particle (see Fig. 1). The time evolution of these coordinates is described by the overdamped Langevin equationsWe ignore the force in the $`y`$ and $`z`$ directions on the chain, and let the rod move only in the $`x`$ direction.

$`{\displaystyle \frac{d𝐫_i}{dt}}`$ $`=`$ $`𝐟_iD_0/(k_BT)+𝐛_i`$ (1)
$`{\displaystyle \frac{dx}{dt}}`$ $`=`$ $`FD_{rod}/(k_BT)+B.`$ (2)
Here $`𝐟_𝐢`$ and $`F`$ are the deterministic forces acting on each particle and the rod, respectively, and $`𝐛_𝐢`$ and $`B`$ are the corresponding random — “Brownian” — forces satisfying $`<𝐛_𝐢(t)>=0`$, $`<𝐛_𝐢(t)𝐛_𝐣(t^{})>=6D_0\delta (tt^{})\delta _{ij}`$, $`<B(t)>=0`$ and $`<B(t)B(t^{})>=2D_{rod}\delta (tt^{})>`$. $`D_0`$ is the diffusion coefficient of an individual binding particle, related to its friction coefficient, $`\zeta _0`$, through the Einstein relation $`D_0=k_BT/\zeta _0`$. As for the chain, we introduce an effective diffusion coefficient, $`D_{rod}=k_BT/\zeta _{rod}`$ which in principle may include all the pore-DNA interactions. Since so little is known about these complicated interactions, we have simply taken $`D_{rod}=D_0/(L/\sigma )`$, consistent with the translational diffusion coefficient of a stiff chain being inversely proportional to its length. In all that follows we use $`\sigma `$, $`k_BT`$, and $`\sigma ^2/D_0`$ as the units of length, energy, and time. $`ϵ`$, the LJ binding energy between the particles and the monomers, is set equal to $`5k_BT`$; the diameter of the spherical cell is $`2R_s=24\sigma `$ and the total length of the chain is $`L=16\sigma `$.

Figure 1 is a schematic snapshot of the simulated system. The particles bind predominantly at the tip and then move back along the chain, with particles occasionally adding in empty spaces along the rod; at the same time, particles leave from other parts of the chain and allow for new particles to bind. Geometrically, up to six particles can sit around each monomer, but just as in the familiar Langmuir adsorption problem there are always empty sites on the chain due to entropic factors. The fact that particles are mostly added at the tip is completely a dynamical effect. Under the influence of particle binding, the chain is in general moving too fast to allow for saturation of adsorption along its length. This results in less than the equilibrium number of occupied sites on the chain. In general, particles bind to the tip of the chain by pushing aside some already attached particles. However, the force pulling in the chain is exerted mostly at the entering positions where empty binding sites first appear.

## III Results

Figure 2 shows the Mean First Passage Time (MFPT) versus the length of the chain inside the cell, $`x`$. MFPT is the average, over a large number of trajectories, of the time it takes for the front tip of the chain to first arrive to the position $`x`$. Each curve in the figure corresponds to a different value of $`N`$, the number of binding particles inside the cell.

As mentioned in the INTRODUCTION, there exist different mechanisms for the translocation of the rod under the influence of the binding particles. If the chain simply diffused into the spherical cell, the MFPT would be equal to $`t_d=L^2/(2D_{rod})=2048\sigma ^2/D_0`$ corresponding to the length $`L=16\sigma `$; see the quadratic function depicted by the heavy solid curve in Fig. 2. As shown in the figure, $`t_d`$ lies significantly above the translocation times we find in our simulations; even in the presence of only five binding particles, the MFPT is about three times shorter.

The dotted line in Fig. 2 represents the time that it would take for the entire chain to enter if there were ratchets functioning perfectly at every absorbing site. In this case, the chain simply undergoes successive and independent diffusions between neighboring sites, completing each in a time $`\sigma ^2/(2D_{rod})`$; the MFPT is equal to the product of this time and the number of steps $`(x/\sigma )`$ associated with the entry distance $`x`$. According to the ratcheting mechanism, then, the slope of time versus displacement is simply $`\sigma /(2D_{rod})`$ =$`8\sigma /D_0`$, for $`D_{rod}=\frac{D_0}{16}`$.

The numerical results from our simulations, as well as the theory outlined in the next section, confirm the presence of a quite different translocation mechanism, namely, drift due to a net force exerted on the chain by binding particles. Entry into the cell of successively longer portions of chain “feeds” new binding sites into the system; as each additional particle binds to the chain the free energy of the chain drops, and this reduction gives rise to a force pulling the chain into the cell. In the presence of a constant force, and in the overdamped limit, the corresponding translocation time is $`t_F=L/v`$ where $`v=FD_{rod}/k_BT`$ is the velocity of the chain.

With competing mechanisms operative, the MFPT will reflect predominantly the one with the smallest translocation time. As already remarked, the diffusion time is always much longer than those arising from the other two mechanisms (see Fig. 2). However, the ratcheting time could be longer or shorter than that of the force-driven process according to whether the deterministic force $`F`$ is larger or smaller than the “effective” Brownian ratcheting force. Comparison of the ratcheting velocity of a chain, $`2D_{rod}/\sigma `$, with the usual expression for the velocity of a chain under a constant force, $`v=F_{ratchet}D_{rod}/k_BT`$, shows that the ”effective” Brownian ratcheting force is $`2k_BT/\sigma `$. Consequently, the ratcheting mechanism is expected to be dominant in the presence of weak enough driving forces. Nevertheless, even in this limit, as we shall see below, the attractive force is operative and completely determines the translocation velocity.

Several important observations can be extracted from Fig. 2. The first is that translocation times depend on the number (concentration) of binding particles. This is also true, of course, for the ratcheting mechanism, since higher concentrations of binding particles imply faster “on-rates” ($`k_+`$), hence higher translocation velocities. But, as Simon, Peskin and Oster have themselves emphasized, there is a maximum translocation rate corresponding to the limit of large $`k_+/k_{}`$. This maximum rate — or minimum time, $`t_{ratchet}=L\delta /2D_{rod}`$ — also corresponds to the smallest distance ($`\sigma `$ in the present model) between ratcheting (binding) sites. From Fig. 2 we see that the translocation occurs even faster than the limiting ratcheting prediction when the number of binding particles exceeds $`N=100`$. This result suggests that an additional mechanism is operative, which we show is associated with a net force acting on the chain along its direction of motion.

More explicitly, the slope of MFPT *vs* $`x`$ plots reveals that the average velocity of the chain, $`v`$, remains almost constant throughout the translocation process (except right at the beginning and towards the end). From the relation $`F=\zeta _{rod}v`$ one expects that the average force on the chain also stays constant during this process, in which case the slope of time $`t`$ versus $`<x>`$ will be almost the same as that of MFPT versus $`x`$ iff $`vx/D_{rod}1`$The translocation time of a polymer can be defined in terms of the mean first passage time $`\tau (x)`$, which satisfies the backward Fokker-Planck equation with a reflecting boundary condition at the hole and an absorbing boundary condition at $`x`$. The solution to the backward Fokker-Planck equation in the presence of a constant force $`F=\zeta v`$ is $`\tau (x)=\frac{x}{v}+\frac{D}{v^2}(e^{vx/D}1)`$. If $`v/D1`$, then $`\tau (x)=\frac{x}{v}`$. In the other limiting case, *i.e.*, when $`v/D1`$, then $`\tau (x)=\frac{x^2}{2D}`$. . Using $`F=\zeta _{rod}v`$, we calculated the effective force on the chain from the velocity obtained in our MFPT *vs* length curves (Fig. 2). Specifically, we calculated the average velocity of the chain over the range $`x=4\sigma `$ to $`x=11\sigma `$, *i.e.*, in the region where the MFPT versus $`x`$ curves are almost linear (see Fig. 2).

Alternatively, we can also determine the force acting on the chain directly from the simulations. The squares in Fig. 3 represent the average force exerted on the chain by the binding particles (coarse-grained over $`\sigma `$), as a function of the length inside. This force is seen to be nearly constant; furthermore, it is found to agree within numerical uncertainty, as shown in Fig. 4, with the force calculated (as described above) from the average velocity (*i.e.*, the inverse slope in Fig.2) for all the values of $`N`$ (up to 300) that we treated.

For $`N<50`$, where the forces are of the order of unity or smaller (i.e., smaller than the “effective” Brownian ratchet force $`2k_BT/\sigma `$), Fig. 2 shows that the translocation times begin to be significantly smaller than the perfect ratcheting limit. This scenario is in principle embodied in the basic result of Ref.  in which the translocation time is written as $`t_{ratchet}(1+2K^1)`$ where $`t_{ratchet}=L\delta /2D_{rod}`$ is the ideal ratcheting time and $`K=k_+/k_{}`$ is the effective strength of binding, expressed as a ratio of ”on” and “off” rates. $`K1`$ corresponds to saturated binding and to ideal ratcheting; otherwise, the ratchet mechanism becomes less efficient and the translocation times are longer than $`t_{ratchet}`$. It is notable that even in this weak force situation, where one expects the ratcheting mechanism to dominate, we still observe a force-controlled translocation (*i.e.,* the translocation time is still determined by the force pulling the chain in, as shown in Fig. 4). The reason is that these weaker forces correspond to smaller numbers of less strongly bound particles, and hence also to less efficient ratcheting.

## IV Dynamical Theory

The translocation process simulated above can be described theoretically using a simplified dynamical model. The two relevant variables are the length of the chain $`x`$ inside the cell and the number of particles $`n`$ attached to the chain. As the chain enters a length $`x`$, the number of available absorbing sites increases. With $`ϵ`$ the binding energy of a single particle we can write the (Langmuir adsorption) free energy, $`A`$, of the system as

$`A(x,n)=nϵk_BT\mathrm{log}{\displaystyle \frac{(6x/\sigma )!}{n!(6x/\sigma n)!}}`$ (3)
$`k_BT(Nn)\mathrm{log}{\displaystyle \frac{V}{(Nn)v_0}}.`$ (4)
Here $`N`$ is the total number of particles, $`V`$ is the volume of the spherical cell, and $`v_0`$ is the volume of a single particle. The coefficient 6 in the second term of Eq. (3) is the number of particles that can interact attractively with a single binding site (chain monomer); accordingly, $`6x/\sigma `$ is the total number of available sites on a chain of length $`x`$. The first term in Eq. (3) is the energy gain due to binding; the second is the entropic contribution associated with the total number of ways in which $`6x/\sigma `$ sites can be occupied by $`n`$ particles; and the last term is the (ideal gas) contribution associated with the “free” particles, numbering $`Nn`$.

We consider the dynamics of translocation as a coupled diffusion process involving both the $`x`$ and $`n`$ degrees of freedom. Using Mesoscopic Nonequilibrium Thermodynamics (MNET) , one can derive the Fokker-Planck equation governing the time-dependent probability density $`\rho (x,n,t)`$ that at time $`t`$ a length $`x`$ of the chain has passed through the hole and has $`n`$ particles attached to it:

$`{\displaystyle \frac{\rho (x,n,t)}{t}}={\displaystyle \frac{}{x}}D_{rod}\left({\displaystyle \frac{1}{k_BT}}{\displaystyle \frac{A(x,n)}{x}}\rho +{\displaystyle \frac{\rho }{x}}\right)`$ (5)
$`+{\displaystyle \frac{}{n}}D_n\left({\displaystyle \frac{1}{k_BT}}{\displaystyle \frac{A(x,n)}{n}}\rho +{\displaystyle \frac{\rho }{n}}\right).`$ (6)
Note that there is a drift and a diffusion term for each of the $`x`$ and $`n`$ variables. The force driving the translocation is $`F(x,n)=\frac{A(x,n)}{x}`$, and we see from Eq. (3) that its origin is entropic, arising from the second (“Langmuir”) term in Eq. (3); the binding of the particles gives rise to a force that pulls the chain in. Similarly, the factor $`\frac{A(x,n)}{n}`$ is a “thermodynamic force” driving the particle binding. $`D_{rod}`$ is the spatial diffusion coefficient of the rod (with the usual dimensions of $`length^2time^1`$), while $`D_n`$ (which in general may depend on $`x`$ and $`n`$) is the *kinetic* rate constant (with dimensions of $`time^1`$) for the process of particle binding and unbinding. A crude, but time-honored and physically reasonable, expression for $`D_n`$ comes from the Smoluchowski theory of aggregation dynamics  for diffusing particles:

$$D_n=a\frac{N}{V_{cell}}D_0,$$
(7)
where $`a`$ is a length of order the particle size. $`D_n`$ here is simply proportional to the concentration of binding particles $`N/V`$, and their spatial diffusion coefficient $`D_0`$.

The Fokker-Planck Eq. (5) provides a complete description of the kinetics of both chain entry and particle binding. However, a simpler description can be achieved by considering the possibility of time scale separation. The characteristic times for the entry and the binding processes scale as $`\tau _x1/D_{rod}`$ and $`\tau _n1/D_n`$, respectively. If the binding process is very fast compared to the chain entry (i.e. $`\tau _n/\tau _xD_{rod}/D_n1`$), it is reasonable to assume that the “fast” variable — here the number of attached particles, $`n`$ — will decay very rapidly to its equilibrium distribution. In this case the process can be described by the evolution of the slow variable, the position $`x`$ of the chain. Suppression of the fast variable can be carried out using the standard adiabatic elimination technique , which is essentially equivalent to integration of the Fokker-Planck equation over the equilibrium distribution of the fast variable. The resulting one dimensional Fokker-Planck equation is

$$\frac{\stackrel{~}{\rho }(x,t)}{t}=\frac{}{x}D_{rod}\left(\frac{\stackrel{~}{F}(x)}{k_BT}\stackrel{~}{\rho }+\frac{\stackrel{~}{\rho }}{x}\right),$$
(8)
where

$$\stackrel{~}{F}(x)=\frac{A(x,n)}{x}f_{eq}(n;x)𝑑n,$$
(9)
is the average driving force. Here we have assumed that the spatial diffusion coefficient of the rod is independent of $`n`$, and defined $`\rho (n,x,t)=\stackrel{~}{\rho }(x,t)f_{eq}(n;x)`$, with $`f_{eq}(n;x)`$ the “local” equilibrium distribution,

$$f_{eq}(n;x)=\frac{e^{\frac{A(x,n)}{k_BT}}}{e^{\frac{A(x,n)}{k_BT}}𝑑n}.$$
(10)
In the particular case where $`f_{eq}(n;x)=\delta (nn_{eq}(x))`$, the force driving the translocation process becomes $`F_{eq}(x)=\frac{A(x,n)}{x}|_{n_{eq}(x)}`$, and the number of attached particles is equal to the equilibrium one, $`n_{eq}(x)`$, given by the solution of $`\frac{A(x,n)}{n}|_{n_{eq}}=0.`$

The above quasi- (or “local”) equilibrium approach basically assumes that, as soon as the chain’s advance makes available new sites, particles bind to them immediately. Figure 3 compares, for $`N=100`$, the average force calculated in our simulation (squares) with the “quasi-equilibrium” force given by Eq. (9) (dashed curve). The comparison shows that the actual force that the binding particles exert on the chain is significantly smaller than the one which follows from the assumption that binding equilibration can keep up with chain entry (and this is true for all other $`N`$ values considered). To check the limit in which the simple quasi-equilibrium 1D Fokker-Planck description, (8), provides an accurate description of the translocation process, we slowed down the entry of the rod by decreasing $`D_{rod}`$. The circles in Fig. 3 shows the average force calculated via BMD simulation when $`D_{rod}=1/960`$ (*vs* 1/16). Agreement with the ”quasi-equilibrium” force evaluated from (9) is now excellent, confirming the validity of our model for the thermodynamic free energy $`A(x,n)`$. Comparably good agreement is found between the number of attached particles calculated in the simulation and that predicted by the quasi-equilibrium theory. (Alternatively, we can increase the diffusion coefficient $`D_0`$ of the binding particles, in which case the calculated force also approaches its quasi-equilibrium value, as we checked in the simulation). The free energy given in Eq. (3) was written in the continuum limit under the assumption that the system is dilute. Considering its simplicity, the accuracy of the model in describing the system is surprisingly good.

In general, to compare simulation results with the predictions of our theory, we need to solve the full Fokker-Planck equation for the coupled $`x`$ and $`n`$ degrees of freedom. The numerical solution is obtained by converting Eq. (5) into its equivalent set of Langevin equations , and then solving these equations using standard stochastic algorithms. From these numerical solutions the average force, the mean first passage time, and the average number of adsorbed monomers can be determined. Using the approximate expression for $`D_n`$, given by Eq. (7), with $`a=10.8\sigma `$ (which in fact is very close to the value $`a=4\pi \sigma `$ from the Smoluchowski theory) we have solved Eq. (5) for different values of $`N`$. The resulting MFPT values *vs* $`x`$, and averaged pulling forces, compare well with those from the simulations. As an example, the resulting average force for $`N=100`$ is represented by the solid line in Fig. 3, showing an excellent agreement with the results of the simulation.

## V Conclusion

The present work has attempted to provide a basic theoretical framework for treating the translocation dynamics of a stiff chain as it moves into a cell containing particles which interact attractively with it and bind to it. We conclude that this process is in general *force*-controlled. Obviously, *pure* diffusion occurs only in the absence of any binding particles; and the *rectification* of diffusion —which is the essence of the Brownian ratcheting mechanism— appears as a manifestation of the binding force in a special limit of interaction potentials and particle concentration.

To examine the extent to which the simple Brownian ratchet mechanism can account for chain translocation due to binding particles, we have performed several simulations using different values for $`\delta `$, the distance between binding sites. We varied $`\delta `$ from $`1\sigma `$ to $`4\sigma `$, *i.e.*, from $`\delta `$ smaller than to larger than the range of the Lennard-Jones attractive interaction (about $`2\sigma `$) . When $`\delta `$ is sufficiently larger than the range of interaction — say for $`\delta =4\sigma `$ — we find that the chain performs a cycle of drift and diffusion movements. As an absorbing site enters the cell the chain is first pulled by the attractive force acting on the binding site, and then diffuses between adsorbing sites, during which time no net force is acting on the chain. In this case, the effect of the binding force is twofold: it pulls the chain in the region where it is acting, and also impedes its backwards diffusion. In the limit where the range of the force is small compared to the distance between sites, the contribution of the drift to the translocation is negligible, but nevertheless the force is still rectifying the diffusion. In these circumstances, the effect of the force can be most simply described by the Brownian ratchet idea. Specifically, the rectification arises from the free energy penalty for moving an absorbing site out of the cell (as provided approximately by our Langmuir adsorption model). The efficiency of the rectification depends on the ratio of free energy penalty to the thermal energy $`k_BT`$.

We have also found that the effective force can be significantly (several times) smaller than the value one would estimate from a quasi-equilibrium treatment of the binding and entry dynamics, i.e., from assuming that the time scale for particles to diffuse and bind is much shorter than that for chain entry. More explicitly, only if we increase sufficiently the friction coefficient of the chain do we find agreement between the force calculated from the simulations and that obtained from the quasi-equilibrium solution to the coupled entry/binding Fokker-Planck equation. This result indicates that there can be significant differences between the *stalling* force measured in single-molecule experiments and the actual *value* of the force during the process of translocation.

In summary, as established directly from our BMD simulations, we find that chain translocation can be understood in terms of the pulling force arising from particle binding. Moreover, this mechanism leads to translocation times that can become distinctly smaller than ratcheting. By decreasing the diffusion coefficient of the chain relative to that of the binding particles, the simulations give a still larger effective force pulling the chain along its length into the cell. The maximum value of this force, attainable in the limit $`N1`$ and $`x\mathrm{}`$ is precisely $`6ϵ/\sigma `$, *i.e.* the drop in free energy per unit length when all binding sites are occupied. These results are nicely confirmed by solving directly a coupled Fokker-Planck equation for the chain entry and particle binding dynamics. We conclude that physically realistic situations are in general more complicated than a ratcheting mechanism in which it is assumed that the only effect of particle binding at a chain site is to prohibit its diffusing back through the pore.

The authors would like to acknowledge helpful discussions with Y. Kantor, M. Kardar, P. G. de Gennes, J. Widom and M. Deserno. This research was supported by the National Science Foundation under Grants No. CHE99-88651 and CHE00-76384.


format: text

text:
# Gigahertz-Peaked Spectrum Radio Sources

## 1. Initial discovery

The history of peaked-spectrum radio sources can be traced back to Ken Kellermann’s graduate student days at the Owens Valley Radio Observatory in the early 1960s. At that time, the revised 3C catalog had just been published (Bennett 1962), and there was great interest in obtaining accurate broad-band radio spectra of these bright sources. For his Ph.D. thesis, Kellermann used the variable spacing interferometer at Owens Valley to obtain flux densities at four frequencies between 750 and 3200 MHz for 160 objects from the 3C list and other radio catalogs available at that time. Approximately half of the sample had straight, steep spectra consistent with incoherent synchrotron emission from a population of electrons having a power-law energy distribution.

The remaining radio sources had somewhat curved spectra, and two unusual objects: CTA 21 (0316+161) and CTA 102 (2230+114), had convex profiles that peaked around 900 MHz (Kellermann et al. 1962), hence the origin of the term “gigahertz-peaked spectrum” (GPS) sources. These two CTA objects had been discovered serendipitously by Dan Harris at Owens Valley several years earlier, and subsequently generated a great deal of excitement in the Soviet Union as possible beacons from extraterrestrial civilizations. Kardashev (1964) had argued that considering the effects of quantum noise and competing background emission, the most efficient spectral profile for transmitting radio signals across the galaxy was remarkably similar to that of CTA 21 and CTA 102. This coincidence fueled intense speculation as to their origin. The situation grew more interesting when Sholomitskii (1965) reported highly unusual flux density variations in CTA 102 at 32 cm. This motivated Kellermann (1966) to look for “artificial” features as predicted by Kardashev (1964) in the spectrum of PKS 1934-638, but none were present.

Kellermann et al. (1962) had realized early on that the 3C sources with distinctly curved spectra had the highest brightness temperatures, but incorrectly attributed this to spectral ageing processes. Slysh (1963) was subsequently able to fit the curved spectra of 3C 48 (0134+329) and CTA 21 using formulae he derived for a synchrotron self-absorbed source, and calculated their angular sizes to be exceedingly small (140 and 10 milliarcsec, respectively). This led Shklovsky (1965) to reason, using his spherically expanding supernova model, that peaked-spectrum sources had to be quite young and should display substantial flux density variations over time.

At around the same time, Dent (1965) reported significant flux density changes in 3C 279, and Kellermann & Pauliny-Toth (1968) began to regularly monitor a sample of bright AGN at wavelengths ranging from 2 to 40 cm. They soon discovered that flux variability was quite common in flat-spectrum AGN, but virtually absent in the peaked-spectrum sources. This posed serious problems for Shklovsky’s simple spherical expansion model for AGN variability, which was later dropped in favor of more complex shock-in-jet models (e.g., Aller et al. 1985). Several groups, including Maltby & Moffet (1965), failed to confirm Sholomitskii’s discovery of variability in CTA 102, although later low-frequency observations of this source (e.g., Hunstead 1972) showed flux density variations that were probably associated with interplanetary scintillation.

### 1.1. Early VLBI observations

By 1966 there was overwhelming indirect evidence that many of the brightest radio-loud AGN were exceedingly compact. Both the high turnover frequencies of the peaked-spectrum sources and light travel time arguments applied to source variability suggested that some AGN could only be studied using very long interferometric baselines exceeding several million wavelengths. Kellermann (1964) had pointed out two years earlier that that the peaked-spectrum sources were all intrinsically very luminous, and therefore could provide a number of good targets for very long baseline interferometry (VLBI).

These findings provided the scientific motivation for several groups around the globe to independently pursue the development of VLBI techniques. Early VLBI observations between Green Bank and Haystack in 1967 at a wavelength of 18 cm (Clark et al. 1968) showed that many AGN were still unresolved, which underscored the need for even longer baselines at shorter wavelengths. The following year, Kellermann et al. (1968) obtained the first inter-continental fringes of several AGN at 6 cm, between Green Bank and Onsala, Sweden. These fringe detections gave the first direct measurement of sizes 1 mas or smaller for several AGN, including the peaked-spectrum source 2134+004.

The direct confirmation of small source sizes using VLBI renewed a debate launched two years earlier by Kellermann (1966) in his analysis of the spectrum of PKS 1934-638: was the low-frequency spectral turnover in peaked-spectrum sources due to free-free absorption (FFA) by thermal plasma, or synchrotron self-absorption (SSA) from relativistic electrons? Or was it related to the presence of cold plasma, as described by the Razin (1960) effect?

Surprisingly, this debate still continues nearly forty years later, with recent papers discussing either FFA (e.g., Bicknell, Dopita, & O’Dea 1997) or SSA models (e.g., Snellen et al. 1999). The main difficulty is that both models contain enough free parameters to fit virtually any convex spectrum, albeit with somewhat different constraints (see O’Dea 1998). Multi-frequency VLBI absorption studies (e.g., Marr, Taylor & Crawford 2001; Shaffer, Kellermann, & Cornwell 1999) may eventually settle this issue, but at present it remains plausible that both mechanisms play a role in creating the spectral turnovers in these sources.

## 2. Classification schemes

The standard view of peaked-spectrum sources today is that they are compact, powerful, young radio sources that reside in gas-rich environments at the centers of active galaxies (O’Dea 1998). Other typically quoted properties include low fractional polarization, low apparent jet speeds, and low flux density variability. However, in reality most individual peaked-spectrum sources display only some of these characteristics, and exceptions to the “classic” definition are unfortunately the norm.

Part of the problem is that the current definition is based only on the overall radio spectrum of an AGN, which is often the sum of many complex, time-variable emitting regions in a relativistic jet. The spectral shape can also be heavily influenced by the external environment, as in the case of FFA. It could be argued, therefore, that a peaked spectrum by itself tells us very little about an AGN, except for the obvious fact that it is dominated by an emitting region of small angular size. This suggests that the term “GPS” should either be abandoned entirely, or qualified in more detail so that it can be used as a true physical classification, and not solely a phenomenological one. Below I suggest some simple ways that the GPS classification could be modified so that it favors intrinsic rather than observer-biased properties.

### 2.1. “Masquerading” blazars

A first step toward a more definitive classification scheme would be to eliminate flat-spectrum, variable AGN (i.e., blazars) from the GPS class. What distinguishes the blazars from GPS sources is their bright, flat-spectrum cores, but in some blazars the signature of the core can be occasionally hidden in its integrated radio spectrum. Kovalev et al. (2002) have shown that the radio spectra of many blazars change dramatically with time, often showing remarkably GPS-like convex profiles during flux density outbursts. This can naturally be explained by the ejection of a new feature in the jet, whose individual self-absorbed spectrum briefly dominates the total radio flux. Many of the “high-frequency peakers” found by Torkinoski et al. (2001) appear to fall into this category.

A small fraction of blazars also display persistent features in their jets that often outshine the flat-spectrum core. These are likely associated with either stationary shocks or bends in the jet, and can give rise to a relatively stable, convex spectrum. In most cases, however, the signature of the flat-spectrum core can be found in another region of the spectrum, well away from the spectral peak of the stationary component. A well-known example is the blazar 4C 39.25 (0923+392), which has a convex spectrum above 1 GHz, but a distinctively flat spectrum over the region $`0.1<\nu <1`$ GHz. Ironically, CTA 102 also shows the distinct flat-spectrum and variability of a blazar at frequencies above 3 GHz, and should be excluded from the GPS category. As it turns out, the only reason CTA 102 was originally called a GPS source is that it was discovered during a relatively inactive state when its spectrum happened to be convex.

In order to eliminate these “masquerading” blazars from the peaked-spectrum class, two additional restrictions should be imposed. First, a bona fide peaked-spectrum source should display a canonical convex spectrum that a) peaks at any radio frequency, and b) maintains this spectrum consistently at all epochs. A canonical peaked spectrum has been constructed by de Vries et al. (1997), and has $`\alpha >0.5`$ for $`\nu <\nu _m`$, where $`S_\nu \nu ^\alpha `$ and $`\nu _m`$ represents the spectral turnover frequency. In the high frequency part of the spectrum well past the turnover where the source is optically thin, the spectral index is steeper than $`0.5`$.

### 2.2. Compact-steep spectrum sources

The current classification scheme is also complicated by the nebulous division between the GPS and compact steep-spectrum (CSS) sources. These are often grouped together based on their many shared characteristics (see the review by O’Dea 1998). For example, when plotted together on the intrinsic size vs. $`\nu _m`$ plane, they form a continuous trend having a single slope ($`\mathrm{size}\nu _m^{0.65}`$; O’Dea & Baum 1997). The CSS sources originally got their name since they were unresolved with connected interferometers ($`\theta <1^{\prime \prime }`$), and had overall steep spectra ($`\alpha 0.7`$). However, subsequent low-frequency observations of CSS revealed that many in fact had spectral turnovers at a few hundred MHz. The term GPS then became synonymous with sources with $`\nu _m1`$ GHz, and those with lower frequency turnovers remained CSS sources.

This completely unphysical, arbitrary division generates additional redshift biases in “complete” samples of CSS and GPS sources, since it depends on the observed turnover frequency. Some authors (e.g., Fanti et al. 1985) have used an intrinsic size criterion of 1 kpc to divide CSS and GPS sources, but again this does not correspond to any physical quantity. Specifically, there is no evidence of a bi-modality in the observed size distribution. Instead, 1 kpc merely corresponds to the typical size of a source with a turnover frequency of 1 GHz. Furthermore, not all GPS/CSS have measured redshifts, making some size determinations impossible. To make matters worse, many intrinsically small GPS sources appear to have associated, large-scale emission, that is perhaps the result of past, episodic jet activity in the source (see review by Lara et al. 2002). It is unclear how to assign a characteristic size to these sources. Finally, there is again the issue of “masquerading” blazars, many of which tend to have small projected jet sizes due to foreshortening. These sources can display steep radio spectra if their flat-spectrum core is under-luminous compared to the jet. This appears to be the case for the superluminal CSS sources 3C 138 (0518+165), 3C 147 (0538+498), and 3C 216 (0906+430).

In light of these difficulties, the CSS classification should probably be reserved for physically small ($`<20`$ kpc) sources that truly show no spectral turnover down to the lowest observed frequencies. Some CSS might very well be foreshortened, but these should be re-classified as blazars if their viewing angles can be constrained, for example, by measurements of superluminal motion. Furthermore, any CSS that shows a MHz-frequency spectral turnover and has a stable, canonical convex spectrum should be combined with existing GPS sources under the blanket classification “peaked-spectrum radio source”.

### 2.3. Compact symmetric objects: a misnomer?

With the advent of hybrid imaging techniques for VLBI data in the early 1980’s, it became possible to further classify GPS and CSS sources on the basis of their parsec-scale morphologies. Phillips and Mutel (1980) popularized the “compact double” class on the basis of VLBI observations of CTD 93 (1607+268) and 3C 395 (1901+319), which showed bright lobe-like features separated by a few hundred pc. As the sensitivity of VLBI arrays steadily improved, more and more “double” sources were seen to have faint emission that connected the lobes to a weak core component. The term “compact double” was soon abandoned in favor of “compact symmetric object” (CSO)

Although the term CSO is now widely used in the literature, it is a somewhat misleading description of these sources. First, the word “compact” can have different meanings depending on the context. VLBI researchers tend to use this term to describe a source with a large VLBI-to-single dish flux density ratio (i.e. brightness temperature), whereas most other astronomers would take it to mean a physically small object. Indeed, most AGN that are classified as CSOs do satisfy both of these definitions, but a growing number of CSOs (and other peaked-spectrum sources as well) appear to be associated with kpc-sized regions of weak, extended emission (Stanghellini et al. 1998).

The “symmetric” part of the class definition is also problematic. Wilkinson et al. (1994) first introduced the term CSO to describe sources such as 2352+495, which possess a fairly high degree of morphological symmetry between their jet and counter-jet. However, the term CSO is now widely applied to any GPS source that shows two-sided VLBI jets, regardless of apparent symmetry. This is illustrated by many CSO candidates in the COINS sample of Peck and Taylor (2000).

The current definition of a CSO, then, relies solely on two properties: a GPS spectrum, and a two-sided VLBI jet morphology, with source compactness and symmetry playing insignificant roles. For the remainder of this article, therefore, I will drop the term “CSO” in favor of “two-sided peaked-spectrum source”.

## 3. Jet kinematics of two-sided peaked-spectrum sources

Since 1994, Kellermann and the VLBA 2 cm survey collaboration have been regularly monitoring a large sample of AGN in order to learn more about the kinematics and evolution of relativistic jets (Zensus et al. 2002). Although the goals of the program are mainly statistical, the sample also includes several peaked-spectrum sources which we are studying in closer detail. Here we report on two particularly interesting two-sided peaked-spectrum sources: 4C +12.50 (PKS 1345+125) and OQ 208 (1404+286). For the remainder of this paper we assume an accelerating cosmology with $`\mathrm{\Omega }_\lambda =0.7`$, $`\mathrm{\Omega }_\lambda =0.3`$, and $`H_o=70`$ $`\mathrm{k}ms^1Mpc^1`$.

### 3.1. 4C +12.50

The powerful radio source 4C +12.50 lies within an ultraluminous infrared galaxy at $`z=0.12`$ that is rich in both molecular and atomic gas. The host contains two optical nuclei separated by $`3.5`$ kpc that appear to have undergone a recent merger, as evidenced by its tidal tail structure and distorted isophotes (Gilmore & Shaw 1986). Our VLBA observations (Lister et al. 2003) reveal a well-defined jet structure with a total extent of $`100`$ mas = 220 pc (Fig. 1). The brightest feature in our full-track VLBA image (A1) has an inverted spectrum and was identified as the core by Stanghellini et al. (2001). Like most peaked-spectrum sources, the integrated polarization of 4C +12.50 as measured by the VLA is very low ($`m<0.2\%`$), however, our analysis has revealed small regions with extremely high fractional polarization in the southern jet. These range up to $`m=30\%`$ at the major bend (D1) and $`m=60\%`$ at the very tip of the jet (D3). The electric vectors in both regions are nearly perpendicular to the total intensity contours.

We have detected two components in the southern jet (A2 and A3) that display substantial apparent motion over our five year sampling interval. Their mildly superluminal speeds ($``$ 1.2 c) and the presence of the bright northern counterjet suggest that the jet axes of 4C +12.50 lie fairly close to the plane of the sky. We obtained a good fit to the overall bent morphology of the jet and counter-jet using a simple constant-wavelength helix of opening angle 23 projected at 82 to the line of sight. Although some helical jets are the result of material being ejected ballistically from a precessing nozzle (e.g., SS 433; Hjellming and Johnston 1981), this is not the case for 4C +12.50, as there is ample evidence for streaming of material along the curved ridgeline of the jet. These include the alignment of the polarization gradient and inferred magnetic field vectors along the jet ridgeline, and the highly polarized flux at specific points in the jet (i.e., the major bend and southern tip). These features can naturally be explained as shocks in the flow that have ordered the magnetic field of the jet. The helical ridge line is likely the result of growing Kelvin-Helmholtz instabilities (e.g., Hardee 1987) that are driven by a small precession of the jet nozzle. Such a precession might be expected for the black hole of a galaxy such as this one that has undergone a recent merger event.

The relativistic speeds of the inner jets are significantly higher than the outer lobe expansion speeds measured for other two-sided peaked-spectrum sources ($`v/c`$ $`0.3`$ c; Owsianik, Conway & Polatidis 1999). The southern lobe is too diffuse to determine an accurate advance speed, but if it is typical of other sources, this implies that a significant amount of jet kinetic energy must be dissipated over a relatively small distance ($`200`$ pc). This is perhaps being channeled into the large expansion of the jet at D1 or a terminal shock at D3.

### 3.2. OQ 208

The radio source OQ 208 is associated with the broad line Seyfert 1 radio galaxy Mrk 668, which like 4C +12.50, shows signs of tidal distortion (Stanghellini et al. 1993). Being one of the closest known peaked-spectrum sources ($`z=0.077`$), we are able to study its jet kinematics with high spatial resolution (1 mas $`1.5`$ pc). Its radio structure is very compact, with $`>99\%`$ of its cm flux originating from a region with a projected diameter of $`12`$ pc (Fig. 1). Aaron (1996) discovered a faint component $`33`$ mas SW of this region at position angle $`107^{}`$. All of the parsec-scale structure lies within a large diffuse halo of diameter $`30`$ kpc (de Bruyn 1990).

We have assembled ten individual VLBA epochs at 2cm from the VLBA 2 cm survey and other programs that span the time period 1995-2001. For the purposes of our proper motion analysis, we have modeled the inner jet structure using eight discrete circular Gaussian and delta-function components labeled C1-C8 in Fig. 1. The precise location of the center of activity in this source is not well known. Our analysis shows the distance between C1 and C5 is shrinking at a rate of 54 $`\mu \mathrm{a}sy^1`$ ($`v/c`$ = 0.28), which immediately rules out the possibility of the core being in the eastern region. Component C6 lies roughly halfway between the eastern and western “lobe” structures, and is a logical candidate for the core. However, the distance between C6 and the eastern structure is increasing at rate of $`180`$ $`\mu \mathrm{a}sy^1`$, while the distance from C6 to the western structure is decreasing by 140 $`\mu \mathrm{a}sy^1`$. If we assume there are no inward motions, the core must lie between C5 and C6, with C6 representing a jet feature that is moving to the SW.

We are unable to determine any individual component velocities without a stationary reference point, however, C3 and C7 are separating at a rate of 60 $`\pm 4\mu \mathrm{a}sy^1`$. If OQ 208 has been steadily expanding at this rate, its parsec-scale emission is only $`1100\pm 73`$ years old. This is typical of kinematic ages derived for other two-sided peaked-spectrum sources (Owsianik et al. 1999).

## 4. Summary

Since Ken Kellermann’s initial discoveries regarding gigahertz-peaked-spectrum sources nearly forty years ago, these powerful radio sources have greatly increased our understanding of active galaxies. They not only provided strong incentive for the development of modern-day VLBI techniques and synchrotron emission theory in compact radio sources, but are now being used to confirm the scenario by which AGN jet activity is triggered by galactic mergers. Unfortunately, the exposure of GPS sources to the larger astronomical community has been somewhat hindered by a complex and awkward classification scheme. However, with the small modifications discussed in this article, it should be possible to implement a scheme that distinguishes GPS sources on the basis of intrinsic, physical properties, and not observer-biased ones.

#### Acknowledgments.

The author would like to thank the members of the VLBA 2 cm survey collaboration for helpful comments on this manuscript.

## References

Aaron, S. E., 1996, Ph.D thesis, Brandeis University.

Aller, H. D., Aller, M. F., & Hughes, P. A. 1985, ApJ, 298, 296

Bennett, A. S. 1962, MmRAS, 68, 163

Bicknell, G. V., Dopita, M. A., & O’Dea, C. P. O. 1997, ApJ, 485, 112

Clark, B. G., et al. 1968, AJ, 153, 705

Dent, W. A. 1965, Science, 148, 1458

de Bruyn, A. G. 1990, in Lecture Notes in Physics 377, Variability of Active Galaxies, eds. W. J. Duschl, S. J. Wagner & M. Camenzind, Heidelberg, Germany, 105

de Vries, W. H., Barthel, P. D., & O’Dea, C. P. 1997, A&A, 321, 105

Fanti, C., et al. 1985, A&A, 143, 292

Gilmore, G. & Shaw, M. A. 1986, Nature, 321, 750

Hardee, P. E. 1987, ApJ, 318, 78

Hjellming, R. M. & Johnston, K. J. 1981, ApJ, 246, L141

Hunstead, R. W. 1972, Astrophysical Letters, 12, 193

Kardashev, N. S. 1964, Soviet Astr., 8, 217

Kellermann, K. I. 1964, ApJ, 140, 969

Kellermann, K. I. 1966, Aus. J. Phys, 19, 195

Kellermann, K. I., et al. 1962, Nature, 195, 692

Kellermann, K. I. & Pauliny-Toth, I. I. K. 1968, ApJ, 152, 639

Kellermann, K. I. et al. 1968, ApJ, 153, L209

Kovalev, Y. Y., et al. 2002, Pub. Astr. Soc. Aus., 19, 83

Lara, L. et al. 2002, New Astronomy Reviews, 46, 89

Lister, M. L., et al. 2003, ApJ, in press

Marr, J. M., Taylor, G. B., & Crawford, F. 2001, ApJ, 550, 160

Maltby, P. & Moffet, A. T. 1965, ApJ, 142, 409

O’Dea, C. P. 1998, PASP, 110, 493

O’Dea, C. P. & Baum, S. A. 1997, AJ, 113, 148

Owsianik, I., Conway, J. E., & Polatidis, A. G. 1999, New Ast. Rev., 43, 669

Peck, A. B. & Taylor, G. B. 2000, ApJ, 534, 90

Phillips, R. B. & Mutel, R. L. 1980, ApJ, 236, 89

Razin, V. A., 1960, Izvestia Vys. Ucheb. Zaved. Roadiofiz., 3, 921

Slysh, V. I. 1963, Nature, 199, 682

Shaffer, D. B., Kellermann, K. I., & Cornwell, T. J. 1999, ApJ, 515, 558

Shklovsky, J. 1965, Nature, 206, 176

Sholomitskii, G. B. 1965, Inform. Bull. on Var. Stars, 83

Snellen, I. A. G., et al. 1999, New Astronomy Review, 43, 675

Stanghellini, C., et al. 1993, ApJ, 88, 1

Stanghellini, C., et al. 1998, A&AS, 131, 303

Stanghellini, C., et al. 2001, A&A, 377, 377

Tornikoski, M., et al. 2001, AJ, 121, 1306

Wilkinson, P. N., et al. 1994, ApJ, 432, L87

Zensus, J. A., et al. 2002, AJ, 124, 662