Characteristics of the polymer transport in ratchet systems

Molecules with complex internal structure in time-dependent periodic potentials are studied by using short Rubinstein-Duke model polymers as an example. We extend our earlier work on transport in stochastically varying potentials to cover also deterministic potential switching mechanisms, energetic efficiency and non-uniform charge distributions. We also use currents in the non-equilibrium steady state to identify the dominating mechanisms that lead to polymer transportation and analyze the evolution of the macroscopic state (e.g., total and head-to-head lengths) of the polymers. Several numerical methods are used to solve the master equations and nonlinear optimization problems. The dominating transport mechanisms are found via graph optimization methods. The results show that small changes in the molecule structure and the environment variables can lead to large increases of the drift. The drift and the coherence can be amplified by using deterministic flashing potentials and customized polymer charge distributions. Identifying the dominating transport mechanism by graph analysis tools is found to give insight in how the molecule is transported by the ratchet effect.

Due to the increased complexity of the models with internal structure, numerical methods play a more important part. This is due to the fact that the ratchet effect occurs in a far from equilibrium environment and only simple model systems can be analyzed exactly (see, e.g., [16,17]).
The internal structure is an important aspect for many real-life molecular motors (e.g., the well-studied kinesin [18]). For systems with non-homogeneous potentials, internal states usually play a more important part than in "traditional" transport driven by biased external forces (such as a constant electric field). Even single particle systems based on the ratchet effect have been shown to display many phenomena, of which the current inversion phenomenon is one of the most interesting. Current inversions are found to be rather common and can usually be generated by tuning of variables (e.g., diffusion constant, friction, potential shape and/or period) [1, [19][20][21][22][23][24][25][26]. In view of this it is reasonable to assume that systems with internal dynamics possess even more surprising properties, and knowledge of the correlation between internal states and transport would enable artificial engineering of the molecules and to boost wanted properties such as the velocity. An intriguing possibility considered in this Article is the control of electrophoresis [27] by modifying the internal charge distribution of the molecule. Because of the large number of parameters and different models, it is hard to compare results from different works and form any universal rules about the current or energetic properties for the ratchet effect. Things get even more complicated for complex molecules, for which the results are even more model dependent.
Therefore we think that it is necessary to at least develop some general methodology for how to systematically study and monitor the behavior of these systems. This is indeed one of the key themes of this work.
We have recently studied polymers using the Rubinstein-Duke (RD) model in timedependent periodic potentials [28]. The RD model [29] is a good prototype of a complex molecule since the size of a linear polymer can be easily varied, it is strongly correlated, and the model has been actively studied for two decades [30][31][32][33][34][35][36]. There has also been in-terest towards polymers as Brownian motors recently [5,6,37]. In Ref. [28] we presented a general "toolbox" based on the numerical solutions for master equations and found current inversions for the RD model in the flashing ratchet and traveling potentials. In this Article we extend our work and methodology by considering the efficiency, different potential time-dependency schemes, non-homogeneous charge distributions and the dynamics of the internal states leading to the macroscopic transport properties. We formulate the operators and master equations that are then solved with suitable numerical tools that fall into areas of linear algebra, integration, optimization and graph analysis. Due to the nature of the ratchet effect, most observables that we are interested in (such as drifts and conformational changes) are very small. Therefore we find that a discrete space model that allows numerically exact solutions provides a very useful framework in this context.
It is found that, like in many other studies on the ratchet effect before, varying certain model parameters has a large effect on the velocity, coherence and energetic efficiency. We take this aspect a step further by doing multiple parameter optimization for the RD model in order to maximize the steady state drift. If the internal states and the movement of the polymer are tightly correlated (such as in the RD model), changing the parameters increases the importance of some molecule conformations over the others. We demonstrate this by comparing the expected values for certain characteristic macroscopic properties for polymers, such as head-to-head and total length. We also identify and compare the most important microscopic conformations of the polymer that are responsible for the currents in different situations.
This Article is organized as follows. In Section II a mathematical framework and notations are defined and in Section III we go through the numerical methodology. In Sections IV and V we present our results for transport properties and examine their microscopic origin.
In the Appendix, the operator algebra involved is discussed in detail. Our conclusions are given in Section VI.

II. MODEL
We study the transport of the RD polymer [29] and its modification, the free motion (FM) polymer [28], in temporally and spatially changing driving potentials. Essentially the RD model consists of connected Markovian random walkers (reptons) in continuous time (see Fig. 1). Each repton carries a charge that interacts with the potential. The model was originally developed to study the reptation process of the polymer in a restrictive medium (gel). However, in the context of this study, the model is primarily used as a good prototype of a molecule with a large number of internal states. To study the importance of the bulk motion, the assumption of the reptation can be relaxed, which results in the FM model. The complexity of the polymers can be increased by considering arbitrary charges of the reptons.
In the following Section, formal definitions of the model are given for the implementation for numerical computations. Readers not interested in the technical details may skip this part and proceed to Section B.

A. Stochastic generator and operators
Consider a one-dimensional discrete Markovian random process in continuous time [38].
After the transition rates between all the allowed states (i, j) in the system are given (elements H i,j ), the stochastic matrix H can be defined. For molecular motors, this matrix includes all the internal conformations and spatial positions of the molecule in the potential [39][40][41]. In the case of the Markovian stochastically driven potential, it also includes the states of the external potential. We consider systems with stochastic (type 1 ) and deterministic potential switching schemes with sudden (type 2 ) and smooth (type 3 ) switching.
The potential V (x, t) is assumed to be L and T periodic in space and time (for stochastic switching T is the expectation value). The stochastic matrix for the polymer dynamics is for the type 1 and for the deterministic case, where T s is the expected lifetime of the potential V s , and q ∈ ℜ N the repton charges. The switching of the potential is assumed to be cyclic, i.e., V 1 → V 2 → · · · → V s → V 1 . The operators n s and h s create transitions between the potential states, also Ref. [28]). In Fig. 1 we have an illustration of the six repton polymer in one of its configurations. We fix the direction of the motion such that up arrows indicate the positive direction and vice versa.
The type 2 operator now becomes and for the type 3 we choose V (x, t) = V max (x) sin 2 (πt/T ). The type 2 and 3 potentials are more reasonable for artificial molecular motors that have external driving mechanisms (e.g., electric potential), whereas the type 1 occurs most likely in nature (e.g., ATP driven motors). After the generator is defined, the dynamics is given by the master equation where the elements of the probability vector P (t) include all the individual states y of the system. The stationary state P stat for the type 1 generator means that HP stat = 0 and for types 2 and 3 that P stat (t) = P stat (t + T where ω k is the corresponding value of the macrostate (e.g., the polymer length), n y a microstate operator, and F k a (large) collection of microstates. For the RD-type model, there are 3 N −1 microstates, for which the operators have the form where the function g(y, i) defines the state (A, B or ∅) of each bond i between the reptons i and i + 1. We define the following four macrostate operators for the RD-type model: the zero-bond count (number of ∅-bonds), the kink count (number of AB or BA bond pairs), the head-to-head length (distance between first and last repton) and the total length (maximum distance between two reptons). The head-to-head and total lengths are calculated in the potential direction (the only spatial direction for the one-dimensional model) and for the fully accumulated polymer they both are zero. The corresponding operator definitions of these observables are found in the Appendix. Separating the head-to-head and the total length is important since the polymer can take a U-shape. For example, for the configuration in Fig. 1 the values for ω k of these operators would be 2 for zero-bonds (formed by reptons 1-3), 1 for kinks (reptons 3-5), 1 for the head-to-head length and 2 for the total length.

B. Selection of the rates
Despite the large number of studies with discrete state Brownian motors, the importance of choosing the rates H i,j has not got much attention. By demanding the local detailed balance (no net currents in equilibrium), the usual choices for the rates are [42,43] 1 where Γ sets the time-scale and 1/k B T is the Boltzmann factor. Both of these constants and the lattice constant, are set to 1 in this paper. All three definitions lead to the required P i = exp(−E i )/Z distribution in equilibrium, but generate the different kinds of dynamics when applied to ratchet systems (far from equilibrium) such that the microstate energy E i contains the potential. To demonstrate this, we have plotted in Fig. 2  Although all three curves for flashing and traveling potentials share a similar shape, the scales are different and large differences can be seen in the limit where the temporal period T → 0. Being fast and simple, the Metropolis form is usually the favorite choice for the rates. But especially with ratchet systems it can be a poor choice since it does not take into account the slope of the downhill moves (rate being limited to 1) that is important for the dynamics. This is also true for the Kawasaki form, since it is basically just a smoothened Metropolis function. Since there is no single correct choice for the rates (based on theory), the selection must be made on experimental or model specific grounds. Only exponential (in flashing ratchet) and Metropolis (in traveling potential) dynamics lead to zero drift in this limit, which is a physically more realistic situation and is also consistent with the single Brownian particle model [44]. Therefore we choose these rates in this study.

C. Non-uniform charge distributions
The usual assumption in the studies concerning polymer transport is that all monomers are identical, i.e., they carry identical charge and mass. We relax this assumption and study the effect of the non-uniform charge distributions along reptons. Previous works on the RD model have considered some aspects of this. In Ref. [34], a magnetophoresis model (i.e., one charged head repton) was considered and in Ref. [45] it was shown that when it comes to the drift velocity all charge distributions are equivalent in small fields (i.e., linear response regime) [65]. In Ref. [46] it was noticed that the drift in constant field depends strongly on the position of the charged repton within the polymer and in Ref. [35] nonhomogeneously charged RD polymers in large fields were studied. Recent study of the dimer in the periodic potential show that if the connected particles are non-identical, directed drift can be generated even in the symmetric potential [7].
We want to find the best possible charge distributions q for the RD and FM polymers by finding the largest possible drifts. This leads to a multi-dimensional, nonlinear constrained optimization problem with constrains coming from the charges q i . We choose i q i = Q and q i ≥ 0, where the first constraint simply sets the total charge corresponding to an uniformly charged polymer and the second one fixes the sign of the charges. The optimal charge distribution gives some (indirect) information about the polymer conformations and reptons that dominate the transport (i.e., have the largest impact on the drift). Lastly we note that optimization has been carried out for some single particle systems [47][48][49].

A. Network analysis
The stochastic matrix H can be also treated as a graph with vertices (states) and edges (transitions) that can be analyzed to gain more detailed information of the transport process, as described in this section. Graphs and statistical physics have a long history due to the close similarities between stochastic systems and electric circuits, and in the seminal work of Schnakenberg in 70's many important results between these two were presented [50] (see [51] for some recent developments). Most of the works on this subject deal with the relations between steady state, rates, probability fluxes and entropy. We are however interested in finding the optimal paths within a current graph, which has not gained interest within previous works. Such ideas have however risen in other disciplines such as microbiology [52]. In the following, only basic knowledge of the graph theory is expected (see, e.g., [53]).
For simplicity, we consider only type 1 scheme where the time-dependency of the stationary state does not need to be explicitly dealt with thus making the numerical computations easier.
After the stationary state P stat of H is found, the net currents (edge weights) between the states can be computed. In addition to the stationary state and stochastic generator, for the type 1 system). For additional details, see Ref. [28]. The graphs G and G sign are then formed as follows. Let a directed edge i → j in G and G sign with weights w i,j and w sign i,j = ±H sign i,j w i,j . With the sign in front of the weights w sign i,j , one chooses the direction of interest of the transport (see below). The weights w i,j are probability flows and the weights w sign i,j are mean displacement flows in the stationary state.
Let γ k be a path i 1 → i 2 → · · · → i k in the graph with i x = i y ∀x = y (the path is non-intersecting) and γ k (j) = i j . We then look for the path(s) with X being w i,j or w sign i,j . The resulting path computed with w i,j contains transitions that lead to the largest mean probability flow and we denote it by γ k . Similarly with w sign i,j , one gets the path with the largest mean probability flow and we denote it by γ sign k . We call these paths the dominating processes. The function f is known as the target function. If the system is closed (periodic), the process must eventually return to its starting state and a cycle is formed, in which case γ k = γ 1 . Since the potentials we study are indeed periodic, we concentrate on closed systems from now on. For the cycle γ sign k , the target function defines the mean cycle velocity, i.e., v Whether there is a difference between γ and γ sign depends on the details of the system.
It may turn out that γ only includes transitions that are not responsible for the directed molecule transportation, but instead results from the non-transporting diffusive motion.
Formally this means that k−1 i=1 H sign γ k (i+1), γ k (i) = 0, which we call a stationary process, as the net transport for the cycle is zero. This is indeed typical for the ratchet transport, since the molecule spends most of its time near the minima of the potential, being unable to move until the suitable state of the potential and molecule conformation is reached. Therefore γ sign carries more interesting information as it takes into account the directions and magnitude of the moves. If the path has a property k−1 i=1 H sign γ k (i+1),γ k (i) = 0, we call it a transporting process. It is not guaranteed that the dominating process is a transporting process in either case.
In the literature, the problem in Eq. (3) for cycles is known as the optimum cycle ratio problem (see, e.g., [54]). The graphs G and G sign may include all states of the system or a fraction of them with the rest summed over, hence the level of the coarse graining can be chosen. For example, if one is interested only on the molecule internal dynamics, summing over all states of the potential may turn out useful. For a RD-type model this would mean that the dimension of the graph is reduced by a factor of 1/SL, which also makes the numerical optimization easier.
The dominating processes simply give a collection of the most probable transitions that the molecule can go through successively, thus giving information about the types of processes that are important. The probability for the (complex) molecule to precisely follow such fixed paths is of course very small. Because of this, it would be hard and time consuming to try to identify dominating processes from the simulation or experimental data. Our proposed graph analysis is simple and can in principle be done for all finite discrete stochastic non-equilibrium systems which have non-zero currents. Whether this analysis is worth the effort (i.e., if γ does contain interesting information), depends on the complexity of the system and the importance of the molecule internal dynamics to the transport process.

B. Motor efficiency
The efficiency of the molecular motor is an important aspect, especially for non-artificial molecular motors that have limited energy available. In the literature, there are several definitions of the efficiency for Brownian motors, see, e.g., Refs. [1, [55][56][57][58][59]. Here we adopt the basic thermodynamic definition that relies on the constant load force F on the polymer, which means that the output power of the motor is vF . The input power W in comes from turning the potential on, thus forcing the polymer periodically in a higher energy state depending on its location. This approach is different from the model where the molecule gains constant amount of energy by, e.g., ATP hydrolysis. We assume that the energy is dissipated, when the polymer goes back to lower energy state, i.e., this energy is not taken into account by reducing it from the input energy. By assuming that transitions between potentials of type 1 system are cyclic (i.e., V 1 → V 2 → . . . V S → V 1 ), the input power for stochastic and deterministic potential schemes can be written as where ǫ and every ǫ s include L3 N −1 states. Since the type 2 potential has discontinuities . The efficiency is defined by η = vF W in . Although the efficiency of the flashing ratchet model is very low for single particles (see, e.g., Ref. [56]), it can be greatly increased for some many particle systems as shown in the recent work [12,13]. Besides the efficiency, we are also interested in the stopping force F stop which, when applied, causes the average drift go to zero. It is expected that the stopping force gets larger as N increases, as seen in Ref. [5].

C. Algorithms
When dealing with large linear systems (of the order of 10 5 states and beyond), one must really pay attention to the convergence properties and therefore the choice of the numerical methods are important. In this Paper we have three types of numerical problems to solve P stat . For the fully stochastic system (type 1) we used the Arnoldi and bi-gradient stabilized (BiGradStab) methods (drift and diffusion), for on/off deterministic system (type 2) adaptive Runge-Kutta 4-5 method and for smooth continuously deterministic system (type 3) quasi-minimal residual (QMR) method. The solution of the type 1 problem is a straightforward eigenstate computation, the other two are more involved integration problems. All computations were performed in Matlab with a modern desktop computer. Solving stationary states for the type 2 and 3 potentials were the most time consuming parts of the computations.
When solving the stationary state for type 1, a random initial vector is good enough choice, but for types 2 and 3 this is not the case. A better initial guess is needed to reduce the computation time. We found that the stationary state of the mean-field operator is easy to compute and a good one to begin with. In many cases, previous solutions can be also used (e.g., when varying T ). A random initial state however serves as a good check of the numerics, since the results must not depend on the choice of the initial state.
The stationary solution for the type 3 can be found with the same manner as for the type 2 (RK45), which however requires that the operator is available for all t ∈ [0, T ] and are either re-build every step or loaded from the memory. The other way (which we used) is to solve the larger linear equation problem as a first order discretization in time, where ∆t = T /M, M being the number of discretization steps. We found that M = 30...60 is accurate enough. In the matrix form this leads to the problem H P = A, where H includes H(t) for all M time-steps and the discretization operator, and the normalization is preserved with A i = 1∀i = LY, 2LY, ..., MLY otherwise zero. As before, the time-dependent diffusion coefficient is found by solving another linear problem. For these linear systems the QMR method turned out to be well converging (LSQR is also a fool-proof method, but very slow).
To maximize or minimize the velocity v(q, T ) for charges and the temporal period, nonlinear optimization can be carried out with the standard sequential quadratic programming method. To find the velocity, the generator H(q, T ) must be constructed several hundred/thousand times because of changes in the transition rates. Efficient implementation presumes that this process is fast, which is achieved for example by manipulating the required matrix elements directly in the memory instead of re-building the whole matrix. The choice of the initial state is crucial (as usual for optimization problems) and a random state is used with several repetitions to confirm the global optimal point. A symmetric initial charge distribution easily leads to a local optimal point with a symmetric charge distribution (as seen in Section IV C). If q is fixed, optimization can be replaced by interpolation, since function v(T ) is very smooth.
The best known exact algorithms to find the optimal cycle ratio have the complexity O(nm) [60], where n and m are number of vertices and edges, but in practice these algorithms are not the fastest ones [54]. We applied an improved version of the Howard's method [61] implemented in the Boost C++ library. There also exist brute-force methods to efficiently find (enumerate) all cycles in graphs [62], but this approach is limited to very small networks and/or cycle lengths. We also tested a simple greedy algorithm where we begin from a single edge with the largest weight and start to grow the path by always choosing the edge with the largest weight available at the moment, until the path form a cycle (i.e., crosses itself). This method however works poorly and an optimal solution is found only for very simple cases (e.g., a polymer in strong static field), where the results are also easy to guess beforehand. In general situations, the optimal path contains transitions that cannot be chosen by a simple greedy algorithm.

IV. RESULTS FOR THE DIFFERENT POTENTIAL AND POLYMER TYPES
Since both RD and FM models include a large number of parameters, some of them must be fixed, primarily those that have a minimal qualitative impact on the results. In addition to N (reptons), other parameters in the models have the following interpretations: • The environment ↔ the potential V (x, t) = V (x + L, t + T ) • The medium ↔ tube deformation Ω (0 for RD, 1 for FM) • The polymer internal fine-structure ↔ charges in q The single most important parameter is the period T of the potential, which is also one of the easiest one to control in experimental set-ups. The parameter Ω models the porosity and viscosity of the medium by either restricting polymer strictly into the reptation tube (Ω = 0) or not (Ω = 1). As before in Ref. [28], we set S = 2 and L = 3 to achieve a both maximal N/L ratio and keep feasible matrix sizes. The flashing ratchet is V 1 (1) = V max , V 1 (2) = V max /2, V 1 (3) = 0 and V 2 (x) = 0 ∀x, and the traveling potential V 1 (1) = V 2 (2) = V max and zero for V 1 (2), V 1 (3), V 2 (1) and V 2 (3). In Fig. 2 of Ref. [28] there is an illustration of these potentials. Time symmetry parameter x = T 1 /T is fixed to 1/2 for the flashing ratchet potential and 1/4 for the traveling one. The maximum potential strength V max has only a small effect on the results and is set to unity (with one exception in Fig. 7) [66]. The direction of the potentials is set up in such way that the expected "main drift" is always positive and the inverse drift (if present) is negative.
With the definitions in Section II, we study the following three types of time dependent potentials • Type 1: stochastic on/off switching • Type 2: deterministic on/off switching • Type 3: deterministic smooth cosine-type modulation.

A. Comparison of time-dependency schemes
First we compare the differences of the potential time-dependency schemes in the flashing ratchet potential, for which the differences are more distinct. In Fig. 3  Some clear differences between the schemes can be seen. The maxima for the drift and the Peclet numbers are reached for smaller T for type 1 than for types 2 and 3. The type 2 scheme has the largest v and type 3 the smallest, and the same goes for Pe. However, this order changes for the inverse drifts, where types 2 and 3 are equally good. The timedependency scheme turns out to have an effect on the current inversion phenomena, since the type 3 scheme is able to invert all RD polymers with N > 2, whereas types 1 and 2 only those with N > 5. Despite this, the differences between types 2 and 3 are small (type 2 being slightly "better") and we now concentrate only on types 1 and 2. with F * = F/F stop and η * max = η max (F * )/ max η max (F * ) for each polymer size, which reveal the shapes of the curves.
We notice that for the FM polymers the efficiency is generally larger and they can maintain their drift in an opposing field better than the RD polymers in the ratchet. When plotted as a function of E, there is a constant stopping field for all N > 3 FM polymers in both potentials with values around −0.0026/−0.0016 for type 1 and −0.0038/−0.0043 for type 2 ratchet/traveling potentials. This results from the fact that the reptons of the FM polymer are less correlated than those of the RD polymer and the FM polymer thus behave more independently . For the ratchet, the type 2 scheme is found to be 2-4 times more efficient and can withstand almost double load force when compared with the type 1. The stopping force is larger for FM polymers. For the traveling potential, differences are more drastic, as for the type 2 scheme the stopping force is about two times and the efficiency almost one order of magnitude larger when compared to the type 1 scheme. Rescaled curves reveal that despite the large differences in scales, shape of the curves are almost identical for all polymer lengths and both types.
The numerical values of the efficiency are very small. This is a generally known trait   In each case, the rightmost curve is for N = 9 and the bold lines (the less interesting special case N = 2) are shared for both RD and FM polymers. Insets: Rescaled data η * max as a function of F * , with black triangles for RD polymers and blue squares for FM polymers.
especially for flashing ratchet models [56], but it also results from the choice of the rates, since the velocity plays dominating role for the efficiency. By the use of the optimized parameters (e.g., V max , x, q), efficiency could be increased by couple orders of magnitude.
Results show that F stop increases as a function of N, which is in agreement with some previous work [5,12]. The efficiency η max however decreases as the polymer gets longer for all other but the type 1 traveling ratchet, which is surprising.

C. Non-uniform charge distributions
Extensive computations were carried out to find the charge configurations with the largest possible v in forward and backward transport and Pe for various polymers and parameters.
It was found that changes in the drift are so large that one can safely limit to maximizing v alone, since in this case Pe is dominated by the drift. In the following, some of the optimization results are presented for the 8-repton polymers in the type 1 potentials. The basic model with an uniform charge distribution (q i = 1 ∀ i) is also shown for comparison.
In In conclusion, the charge distribution has a large effect on the polymer transport velocity and coherence on the flashing and traveling potentials. Since the drifts generated by the ratchet effect are generally very small and difficult to observe, this could be of interest from the point of view of applications. In the next Section we show that different distributions also lead to different kinds of transport mechanisms.    For the ratchet, the maximum positive current (black lines) is a result of small changes in the polymer average shape, which is caused by the fact that only a single near-head repton is charged and the rest of the polymer is in pure random motion. The maximum negative current (blue (light gray) lines) however is a result of more complex processes, which cause much more variation in the average shape, even more than for a polymer with uniformly For the traveling potential, the curves are more distinguished from each other and are  kink dynamics, large differences are shown in zero-bond dynamics. Note that for positively optimized polymer, values remain unchanged during "on → off" switching and are therefore not shown in the figure. This is because, in the steady state, the potential has no effect on the conformations of the polymer, which would require more than one charged reptons. For the traveling potential, the time-evolution of the observables is more complex.

B. Network analysis
To further understand the formation of the net drift, we now turn to the network analysis of the steady state currents. We concentrate on the RD polymer of the type 1 in the flashing ratchet and the traveling potential with uniform and optimized charge distributions. The temporal periods T are chosen such that they result in the maximum current (4 values of T for both potential types). The graphs G sign containing the steady-state net currents between the states are then computed. We have summed over all the potential states (SL degrees of freedom) so that only the internal states of the polymer remain. After these steps we have eight different graphs with 5832 non-zero directed edges in each of them.
Let us first analyze these G sign graphs by defining the arrays S with elements S i (i = We now turn to the dominating transport cycles of the polymer motion by analyzing the paths in G sign . This results in cycles with lengths of the order of 10. It is found that the common transportation type is such that we call "s 1 -s 2 -scheme" consisting of cyclically accumulated (lengths s 1 and s 2 with |s 1 − s 2 | = 1) and elongated parts of the polymer.
Corresponding to the direction of moves, this scheme can be either positive (up) or negative (down). To illustrate the scheme, we have sketched the positive 4-5 scheme in Fig. 15. The numbered arrows indicate the order and direction of the corresponding repton moves. After all marked moves are done, the initial state is recovered and the cycle is repeated. In the five situations out of eight studied here, the dominating cycle is the s 1 -s 2 -scheme.
In Fig. 16 we show the remaining three situations that are not of the type above. Note that for negative transport in the ratchet with the uniform charge distribution, the mechanism is almost the negative 4-5-scheme.
In Table I   is not much difference between the leading mechanisms for forward or backward motion, and for uniformly charged polymer in ratchet it is actually the same. One can therefore conclude that the current inversion for the RD model is not caused by some abrupt 'phase transition', but gradual changes in the probability distribution along internal states.
We carried out a similar analysis also for the full system without summing over S and L, in which case cycles have up to 30 states and there are some modifications to the pure s 1 − s 2 schemes. However, these cycles are too lengthy to be reported here. It was found that sometimes summing over the potential states is necessary to find a non-stationary cycle and sometimes the summing leads to stationary cycle.

VI. DISCUSSION
We have analyzed the properties of Rubinstein-Duke polymers with some modifications, including tube breaking and non-uniform charge distributions, in time-dependent potentials.
The aim of this work was to further study the properties of complex molecules in out-ofequilibrium conditions and especially the ratchet effect.  In the first part of the study, we extended the previous work reported in Ref. [28] by considering deterministic ratcheting mechanisms, the energetic efficiency and optimized charge distributions of the polymers. It was found that the deterministically flashing potential is superior when compared to a smoothly varying and stochastic potential for velocity, coherence and efficiency. However, despite "scaling differences" in drift and diffusion, the time-dependency scheme seems to have a minor effect on the qualitative results. By using the stochastic scheme, we computed the optimal charge distributions to maximize the steady-state velocity in flashing ratchets and traveling potentials. The differences between these and the uniformly charged polymers were found to be drastic. Changing the charge distribution also changes the mechanism of how the polymer reshapes itself with respect to the potential.
In the second part, the current inversion phenomenon was investigated in detail by using the optimal charge distributions. The expected values of certain macroscopic observables (e.g., length and zero-bond count) were computed and large differences between differently charged polymers were found. To find how the polymer actually moves in the non-equilibrium steady state, we proposed a simple graph analysis method to find most probable series of state transitions (=path) based on the probability currents. For a pe-riodic system such a path is found as a solution of the optimal cycle ratio problem. This method is suitable in situations where a huge network is generated by some automated fashion or measurements and cannot be analyzed "manually" (e.g., Kinesin network in Ref. [63]).
This method was then used to identify the dominating processes of the polymer transport and was found to be very useful to piece together polymer motion. However, the general usefulness of this analysis depends on the model and it would be of interest to test it for other complex out-of-equilibrium systems and also with non-periodic boundary conditions. M i,y,l (q) = R(q, l + f (i, y))(n A,i,y,l n ∅,i+1,y,l + n ∅,i,y,l n B,i+1,y,l − a i,y,l a † i+1,y,l − b † i,y,l b i+1,y,l ) + L(q, l + f (i, y))(n ∅,i,y,l n A,i+1,y,l + n B,i,y,l n ∅,i+1,y,l − a † i,y,l a i+1,y,l − b i,y,l b † i+1,y,l ) + ΩR(q, l + f (i, y))(n A,i,y,l n B,i+1,y,l + n ∅,i,y,l n ∅,i+1,y,l − a i,y,l b i+1,y,l − b † i,y,l a † i+1,y,l ) + ΩL(q, l + f (i, y))(n B,i,y,l n A,i+1,y,l + n ∅,i,y,l n ∅,i+1,y,l − b i,y,l a i+1,y,l − a † i,y,l b † i+1,y,l ), where Ω = 0 for RD polymers and Ω = 1 for FM polymers, and where n z , n h , n g = 0, 1, ..., N − 1 and n k = 0, 1, ..., N − 2. One can verify that #F G i ≥ #F H i holds for all i. By using above sets F and equation (2), measure operands can be constructed and expected values computed. The practical procedure to form all the required operators, especially the previous observables, is explained below.

Operator construction
Since the stochastic generator and measurement operands used in this work are slightly more complex than in the previous works regarding the RD model, we show in some details how the idea of the recursive operator construction work in the current case. Whereas small operands can always be build directly, recursive construction is practically a must for large systems and nowadays widely used in DMRG computations [34,64]. For simplicity, we concentrate only on (discrete) state measure operators, which in the natural basis are diagonal matrices.
Let O i 1 , ..., O i y i be a set of macrostate operators for the system with i sites, which includes all the necessary operators that are required when adding a new site. Here site is a general term, which for example could mean single particle states for classical systems and spin states for quantum systems. By using the usual product state formalism, assume that the new sites are added on the right such that |new state = |old state ⊗ |new site .
The basic algorithm to add new sites (until N) goes as follows The number of required operands is therefore ∝ m 2 for total length and ∝ m for others.
We now consider a concrete example for a total length operator, which is the most complex operator used in this paper. When one enlarges the size of this operator with new particles, one must keep track of the maximum distances of the rightmost repton from all the other reptons. For example, in Fig. 1 these distances would be 2 (from repton nr. 4) and 0 (no any reptons below). We define these as up (u) and down (d) distances. Total distance is then d + u.