Molecular quantum robotics: particle and wave solutions, illustrated by “leg-over-leg” walking along microtubules

Remarkable biological examples of molecular robots are the proteins kinesin-1 and dynein, which move and transport cargo down microtubule “highways,” e.g., of the axon, to final nerve nodes or along dendrites. They convert the energy of ATP hydrolysis into mechanical forces and can thereby push them forwards or backwards step by step. Such mechano-chemical cycles that generate conformal changes are essential for transport on all different types of substrate lanes. The step length of an individual molecular robot is a matter of nanometers but the dynamics of each individual step cannot be predicted with certainty (as it is a random process). Hence, our proposal is to involve the methods of quantum field theory (QFT) to describe an overall reliable, multi–robot system that is composed of a huge set of unreliable, local elements. The methods of QFT deliver techniques that are also computationally demanding to synchronize the motion of these molecular robots on one substrate lane as well as across lanes. Three different challenging types of solutions are elaborated. The impact solution reflects the particle point of view; the two remaining solutions are wave based. The second solution outlines coherent robot motions on different lanes. The third solution describes running waves. Experimental investigations are needed to clarify under which biological conditions such different solutions occur. Moreover, such a nano-chemical system can be stimulated by external signals, and this opens a new, hybrid approach to analyze and control the combined system of robots and microtubules externally. Such a method offers the chance to detect mal-functions of the biological system.


Introduction
Molecular robotics, which operates on a nano scale, has in the last decade witnessed impressive growth, (e.g., Murata et al., 2013). The current topics in this field are molecular machines (Balzani et al., 2008;Roux, 2011;Fukuda et al., 2012;Seeman, 2014) and, even more important in the context of this contribution, biped DNA walkers (Sherman and Seeman, 2004;Shin and Pierce, 2004;Omabegho et al., 2009;Lund et al., 2010). In a previous contribution, the biological process of muscle contraction by the "one-legged" motion of myosin II along actin-filaments has been described using QFT methods (Haken and Levi, 2012). The results of different types of synchronization of multimolecular systems have been conclusive in the sense that beside the more classical impact solution, two quantum mechanical solutions exist that can describe the synchronization of billions of unreliable molecules. This result can usually not be modeled by pure classical particle solutions. Despite the theoretical existence of these innovative solutions, knowledge of a great set of experimentally confirmed data is unfortunately still now very limited. Nevertheless, these theoretical outcomes encourage us to continue with our quantum theoretical approach, since we are convinced that in the near future confirming experimental data will be available.
This paper focuses on a description of the "leg-over-leg" motion of molecular robots along tubulin strands. These biological nano-robots (motor proteins) walk along such lanes in a four-step process, and are fueled by ATP consumption (in physics, energy exchange with a heat bath) after the first and third steps. From a system perspective, this denotes that we are modeling an open system that is not in a thermal equilibrium. The C-terminal domain, e.g., of kinesin-1, act as a gripper and is attached to cargo, e.g., vesicles that are transported along axons and dendrites, which are both connected to neuron cell bodies. Very similar processes occur for the much larger dynein, with two legs (in biological terms two "heads"), which also transports vesicles along axons and dendrites. Due to this similarity, this paper focuses on a general description of the motion processes of both molecular robots. We abandon the modeling of cargo transport in this paper for the sake of understandability.
In this paper, a set of such molecules is regarded as a swarm of molecular robots that must be configured, synchronized, and in real time supported by energy to perform the next step (Levi and Haken, 2010). Due to the utilization of a quantum theoretical approach, such molecules can be described as particles or as waves (according to matter/field dualism). The following three different solutions are presented: (1) Impact (stroke) solution of one molecular robot on a single substrate lane, (2) Coherent motion of many robots on and between parallel substrate lanes, (3) Running wave that synchronizes the motions of many robots on one substrate lane.
In all three solutions, a local B-field is activated within a lane and across different lanes. Such a field not only controls the process of energy consumption but also, even more importantly, acts as a signal field that synchronizes the steps of the molecular robots. This field is generated by the molecular robots themselves and is not externally injected. The impact solution pertains to the motion of a single molecular robot on one lane. Damping effects transform the typical wave characteristics of a QFT approach to particle behavior. The appearance of a sequence of impacts pushes the molecule (particle) forwards or backwards step by step.
The next two solutions plainly reveal the wave features of our approach. The coherent motion solution distinguishes between the individual robots walking on a lane and across the lanes. In contrast to the second solution, the running wave solution presents a result that can be obtained in special restrictive conditions.

The Model
The "leg-over-leg" walking of a kinesin molecule (dynein) is modeled as a bipedal molecular robot r, where the two heads (light chains) are considered as two legs, and the "coiled-coil tail" as an effector (Alberts et al., 2008). The track is a microtubule surface (substrate lane), along which r moves in discrete steps. During such a walk, a robot r can take one of four leg states and one of two walking states: Leg state a: leg 1 or leg 2 points backwards. Leg state b: leg 1 or leg 2 points forwards. Leg-over-leg state c 1 : leg 1 is loosely bound to the top of leg 2, which is tightly attached to the substrate. Leg-over-leg state c 2 : leg 2 is loosely bound to the top of leg 1, which is tightly attached to the substrate. Walking state 1: r is moving forward with leg 1. Walking state 2: r is moving forward with leg 2. Figure 1 shows the four different leg states.
The two legs are attached to the substrate in two different connection states: Connection state 1: leg 1 is tightly attached to position p k of the tubulin substrate and leg 2 is weakly attached to site p k + 1 . Connection state 2: leg 2 is tightly attached to site p k and leg 1 is weakly attached to position p k + 1 of the substrate.
The level of attachment to the substrate in these two connection modes is defined by the two possible states of the substrate: Ground state g: represents a weak attachment of a leg to the substrate molecule. Excited state e: represents a strong attachment of a leg to the substrate molecule.
The discrete periodic movement of the molecular robot r is characterized by a four-step cycle (Figure 2) that starts with the two leg states (1a, 2b), continues with the transference of these states in step 1 into the leg-over-leg state (1c 1 , 2c 1 ), and is maintained by step 2, which produces the leg states (2a, 1b). In step 3, the leg-over-leg state (1c 2 , 2c 2 ) is achieved, and finally, in step 4, the states (1a, 2b) are again established. The following notation is employed for the creation (or annihilation) operators: Molecular robot r: r † leg 1 state j position k ; leg 2 state j ′ position k ′ ; e.g., r † l 1 ap 1 ; l 2 bp 2 . Substrate s: s † state j position k ; e.g., s † gp 1 .

State a State b
State c 1 State c 2 1 , 2 1 , 2 = 1 2 = 1 2 FIGURE 1 | Representation of the two leg states a and b and the two leg-over-leg states c 1 and c 2 . The walking direction is from left to right. For simplicity, we present in this diagram each leg as tightly connected to the substrate (crossed circle). The loosely bound states connected to the substrate are omitted (see Figure 2). FIGURE 2 | Demonstration of one complete "leg-over-leg" cycle of a molecular robot r which walks along a substrate lane. The loosely bound states are marked by crossed circles . The tightly bound states are marked by striped circles . The first position is fixed by k = 1; the initial state is (1a, 2b).
TABLE 1 | Representation of the motion patterns of a "leg-over-leg" walking process of a molecular robot r.
: tightly bound; : looseley bound; : leg up. The initial state is (1a, 2b). Table 1 subsumes the different individual motion patterns of the two legs of r.
The energy transfer process (s −→ r) occurs in steps 1 and 3 and is combined with the excited status of the substrate. In addition, the heat-bath operators ("fueling" and synchronization operators) B † p k and B † p k are applied in this walking phase. The operator transfers during the four steps are summarized as follows (readout from Table 1): initial state step 1 "leg-over-leg" state r † l 1 a p 1 ; l 2 b p 2 s † e p 1 s † g p 2 −→ r † l 1 c 1 p 2 ; l 2 c 1 p 2 s † g p 1 s † e p 2 B † p 2 . "leg-over-leg" state step 2 reversed initial state . reversed initial state step 3 "leg-over-leg" state This transfer list of operators delivers the definition of the interaction Hamiltonian H int and finally the elaboration of the resulting equations of motion.

Interaction Hamiltonian
The local interaction Hamiltonian H int at positions to p 1 to p 4 is defined by expression 3.1, where the first step starts at the initial site k marked in Figure 2. Further, two real coupling constants, g 1 and g 2 , are introduced, where the parameter g 1 describes the uneven steps and g 2 denotes the even steps.
H int =hg 1 r † l 1 a p 1 ; l 2 b p 2 s † e p 1 s † g p 2 r l 1 c 1 p 2 ; l 2 c 1 p 2 s g p 1 s e p 2 B p 2 +hg 1 s † g p 1 s † e p 2 r † l 1 c 1 p 2 ; l 2 c 1 p 2 s e p 1 s g p 2 r l 1 a p 1 ; l 2 b p 2 B † p 2 +hg 2 r † l 1 c 1 p 2 ; l 2 c 1 p 2 s † g p 1 s g p 3 r l 1 b p 3 ; l 2 a p 2 B † p 2 +hg 2 r † l 1 b p 3 ; l 2 a p 2 s † g p 3 s g p 1 r l 1 c 1 p 2 ; l 2 c 1 p 2 B p 2 +hg 1 r † l 1 b p 3 ; l 2 a p 2 s g p 2 s † e p 2 s e p 3 s † g p 3 r l 1 c 2 p 3 ; l 2 c 2 p 3 B p 3 +hg 1 r † l 1 c 2 p 3 ; l 2 c 2 p 3 s e p 2 s † g p 2 s g p 3 s † e p 3 r l 1 b p 3 ; l 2 a p 2 B † p 3 +hg 2 r † l 1 c 2 p 3 ; l 2 c 2 p 3 s † g p 2 s g p 4 r l 1 a p 3 ; l 2 b p 4 B † p 3 +hg 2 r † l 1 a p 3 ; l 2 b p 4 s † g p 4 s g p 2 r l 1 c 2 p 3 ; l 2 c 2 p 3 B p 3 . (3.1) By the replacements p 1 → p k , p 2 → p k + 1 , etc., H int is transferred into H int, k and has to be summed up by k =

Full Heisenberg Equations of Motion
The full set of Heisenberg equations of motion is completed by 7 damping constants, γ ab , etc., 7 fluctuating forces, F † ab , etc., and two coupling constants, g 1 , g 2 , and is defined by the following expressions, representing coupled nonlinear, delayed, and complex operator equations. (4.6) To solve the Equations (4.1-4.7), we utilize the semi classical approach and replace the operators by their expectation values in a coherent state representation. In doing so, the operators become complex numbers and the expectation values of the fluctuating forces are zero.

Solutions
We start the description of the three a fore-mentioned solution with the impact solution and are firstly searching for solutions for each step separately. Later, one overall solution will be outlined by combining all four separate solutions. We are hereby concentrating on one robot walking on one lane. The walking processes of several robots on different lanes will be explained in the succeeding sections.

Impact Solution: First Step
As mentioned above, the procedure begins with the first step (Figure 2), at position k. The initial state is defined as (starting with position k = 1): where | φ 0 delineates the vacuum state. The initial conditions of all relevant states are: We set r † l 1 c 2 p 1 ; l 2 c 2 p 1 = 0 in Equation (4.1) because this state is not present. The resulting modified Equation (5.1) now reads: In a similar way, going through all remaining equations the following reduced set of equations is obtained: Before presenting the complete set of numerical solutions for all Equations (5.4-5.10), we prefer to assure ourselves that the right solution can also be generated by analytical methods and not only by numerical calculations. This should convince us that we are on the right track.
The following approach (setting equally the robot-relevant damping constants γ = γ ab = γ ba = γ c 1 = γ c 2 , and neglecting for the moment the damping of the B operator: γ B = 0), delivers a consistent, periodic solution for a walking biped molecular robot: (5.13) By inserting these expressions into the Equations (5.4-5.10), four coupled equations for the four variables f, g, h, b are finally obtained:ḟ = g 1 4 sin 2g sin 2h b. (5.14) The last equation can be directly solved by the technique of separation of variables: The calculation is started with the assumption that γ = 0 is valid in order to start with "perfect" symmetry. Afterwards, we will switch to the damping process in order to observe how this symmetry will be broken. Incidentally, it is not surprising if all solutions of the Equations (5.14-5.17) are periodic (wave-like) if the damping process is switched off. Figure 3 demonstrates this prediction in a phase portrait (Holmes et al., 2012) that also includes the expected periodicity of the variable b (B-field): Even if the general complex solutions of the Equations (5.4-5.10) are calculated, symmetrical patterns are visible. Figure 4 shows the real parts of all operators with the expected symmetry. The imaginary parts show a very similar behavior; therefore they are not presented here or in the next sections.
It is well known that processes without any damping are artifacts. Therefore, we will now include damping processes. Even if "perfect symmetry" is not realistic, it is worth taking it into consideration since it serves as a kind of "roadmap" to the new individual trajectories if γ deviates from zero and becomes positive.
Even if the damping constant is assigned the small value γ = 0.0005 the symmetry (periodicity) will be broken. The limit cycle shown in Figure 3 is dissolved and the trajectory goes directly to zero (fixed point). The same happens with the general solutionit goes directly to zero. Thus, in both cases a large γ value should be chosen.
If the value of γ is continuously decreased, the ringing effect increases along the trajectory to γ = 0. Therefore, we can achieve all possible trajectory patterns, from a direct path to zero to complete oscillation (no damping).
Even if the value of the damping constant slightly diminishes, to γ = 0.0001, the first effect can be observed; after the first peak the trajectory converges to zero (Figure 5).
The results of all remaining steps, 2-4, will not be presented in the same detail as in step 1. Only the actual initial conditions and the corresponding equations of motion are given. At the end of this subchapter, we present all four steps in sequence.

Impact Solution: Second Step
The description continues with a presentation of the second step. For reasons of clarity, we set k + 1 = 2; the initial state is defined by: . The value of the damping constant is set to γ = 0.0001. The coupling constants are g 1 = g 2 = 0.1.

Impact Solution: Third Step
The presentation continues with a description of the third step. For reasons of comprehensibility better readability, we set k+2 = 3. The initial state is fixed by: In this case, the initial conditions of all states are: The Heisenberg equations of motion of the expectation values of step 3 are:ṙ † l 1 b p 3 ; l 2 ap 2 = ig 1 s e p 2 s † gp 2 r † l 1 c 2 p 3 ; l 2 c 2 p 3 Attention is drawn to the symmetry that exists between all equations of motion of the first and third steps. Further, the initial condition for s † gp 2 in step 3 is s † gp 2 (t 2 ) = 0 , while in step 1 s † gp 2 (t 0 ) = 1 was valid. Similarly, the value s † ep 2 (t 2 ) = 1 is assumed in step 3, whereas s † ep 2 (t 0 ) = 0 was correct in step 1.

Impact Solution: Fourth Step
The presentation continues with a description of the fourth step.
Here, we set +2 = 4, and the initial state is given by: where | φ 0 defines the vacuum state. In this case, the initial conditions of the two robot states are: The initial states of the substrate at the three positions p 2 , p 3 and p 4 are: The heat-bath operator satisfies the initial condition The relevant Heisenberg equations of motion of step 4 are:

Impact Solution: All Four Steps in Succession
Now all four partial solutions are collected and put together into one sequence. The first consecutive view describes the changes in the four robot states during the four separate walking steps. Figure 7 illustrates the corresponding real parts (imaginary parts are again very similar, therefore their representation is skipped) of the transitions of these walking states on a reduced scale of 5000 steps for better visibility. Clearly, two typical effects are observable. The similar steps 1 and 3 are regularly performed and go directly from their initial values to zero. Steps 2 and 4 also conform to this pattern, but one on a broader scale, because step 4 shows for r † l 1 c 1 p 2 ;l 2 c 1 p 2 slightly more oscillations than for r † l 1 ap 1 ;l 2 bp 2 in step 2. Figure 8 summarizes the sequence of the real parts of the ground states of the substrate for all four steps. Here another effect comes to light. The behavior of s † gp 3 is different in step 2, and during step 2 (see Figure 2), the ground state s † gp 3 does not exist and has to be generated, whereas during step 3 this ground state already exists and must be annihilated and transformed into an exited state.
From a purely mathematical point of view, we have to compare the two Equations (5.25) and (5.34) for s † gp 3 in the second and third steps. In the second step, there is a strong mutual dependence between s † gp 3 and s † gp 1 . In Equation (5.34), we mainly have to consider the dependence of s † gp 3 and s † ep 3 . The operator product s ep 2 s † gp 2 is time-independent and self-consistent, while substrate states g and e produced themselves reciprocally (see also Figure 10). The three remaining bordering operators (2 r-operators, one B-operator) are equivalent in both expressions.
The next focus of attention is the behavior of the synchronizing heat-bath operators B † p 2 and B † p 3 . Here, we consider at first the effect of the damping value γ B = 0.001 and fix all other damping constants to zero. Figure 9 reveals the behavior of the overlay of the real parts (gray) and imaginary parts (red) of these two operators. Concerning the real parts, both operators are "strong" during steps 2 and 4 (the fueling process) and "weak" in the two other steps. During step 1, B † p 2 starts with B † p 2 (t 0 ) = 0 and ends with B † p 2 (t 1 ) = 1. In step 2, it starts with this initial value, B † p 2 (t 2 ) = 1, and ends permanently with B † p 2 (t 2 ) = 0. Such interplay is valid in the same manner for B † p 3 (t 3 ). The trajectories of the imaginary parts of these two operators are distinct in that way that the real parts of these two operators in steps 2 and 4 are no longer dominant.
The situation greatly changes if we decrease the corresponding damping constant. Figure 10 portrays the temporal diagram of the overlay of the real and imaginary parts of B † p 2 and B † p 3 . The contours of each step can be considered as an envelope around the oscillating solutions.

Synchronized Wave-based Motion
Frontiers in Neurorobotics | www.frontiersin.org which defines a two-dimensional wave number vector. This second vector has been introduced since the operator B † l will be developed by a plane wave approximation: The synchronization activities of the two heat-bath operators B † l , B l on the molecular robots are described as external signal field modes that are coupled by the constant J ll ′ .
We insert this hypothesis in the "parallel" extension of Equation (4.7), indicating that our approach now applies to parallel M lanes. This therefore implies that l is fixed and over all paths we sum up m at the various positions with respect to the fixed position p m of l. The use of the coupling constant J ll ′ forB † k requests a summation along a lane and between lanes: (5.48) We multiply both sides of (5.48) with 1 √ L l e −ik ′ ·l and use the orthogonal relation 1 L l e i(k·l−k ′ ·l) = δ kk ′ : (5.49) The last, rather complicated term can be simplified if we introduce the center of gravity x = l + l ′ /2 and the distance d = l − l ′ . Further, it is supposed that the coupling constant J depends only on the distance J ll ′ = J d and the lane length L is defined by = L x L d . It follows that: where ω k has the dimension of a frequency. We continue with the simplification process and set if all damping constants are set equal to Finally, after all these simplifications the resulting formula is obtained: . The value of the damping constant is set to γ = 0.0005. The coupling constants are g 1 = g 2 = 0.1. . The damping constant is set to γ B = 0.001. The coupling constants are g 1 = g 2 = 0.1. + r † (l 1 a p m ; l 2 bp m + 1 ),m r (l 1 c 2 p m ; l 2 c 2 p m ),m s gp m − 1 , m s † The basic solution to Equation (3.67) can be achieved if only the term k ′ = 0 that corresponds to the mode with an infinite wavelength is kept. We will now solve this "reduced" equation. The last term, ω 0B † 0 , can be deleted in the interaction representation and the expectation value ofF † 0 can be ignored. These two modifications lead to the following equation: calculations of this subchapter show that the synchronization of motion can easily be accomplished (at least forB † 0 ) by replacing the coupling constants and performing an addition of all lanes at the same positions.

Running Wave Solutions
We again concentrate on a single robot that moves on a lane and restrict ourselves for simplicity to step 1 since the remaining steps can be handled very similarly. Furthermore, the damping constants with respect to the molecular robot are equaled: γ ab = γ ba = γ , the damping constants of the substrate are inoperative: γ g = γ e = 0, and the damping constant of the B-field is assumed to be effectively γ B = 0. The positions on the lane are p k = k, p k + 1 = k+1,... Under these assumptions, the approach that guides us to a running wave solution is formulated as follows: r † l 1 a p k ; l 2 b p k + 1 = R † ab e iK/(k+1) , (5.57) r † l 1 c 1 p k + 1 ; l 2 c 1 p k + 1 = R † c 1 e iK/(k+1) , where K is a real variable that will be later defined. Insertion of these expressions into the Equations (5.4) and (5.5) provides the following formulas: In doing so, we introduced the following abbreviations: D 1 = S e S † g = S † g S e , D 2 =S gS † e =S † eS g , D = D 1 D 2 = D 2 D 1 .
FIGURE 11 | Undamped temporal curves, showing the real parts of the operator products D 1 , D † 1 (abridged to D 1 k), D 2 , D † 2 (abridged to D 2 k), D, D † (abridged to Dk) and B † p 2 B p 2 , consecutively outlined for step 1. The damping constant is γ = 0, the coupling constants are constant: g 1 = g 2 = 0.1.
The corresponding bosonic operators (and Hermitean adjoint operators) commutate, thus they can be interchanged. Figure 11 shows the -like peaks of these operator products at the same time slots. At the same time, e.g., the ground state g will be annihilated and the exited state e created. The periodicity of this self-consistent process is due to its calculation without damping (revealing again the original mathematical symmetry). In a more realistic case, damping effects are activated and the curves converge after the first peak to zero. Therefore, for a fixed time t, both the real and imaginary parts of D and/or D † are constant expressions.
In the next step, we reformulate the Equations (5.6-5.10) by including the D-terms: The calculations are continued by substituting the three operators R † ab , R † c 1 and B † b by (5.71) The derivatives of the first two expressions read as follows: (5.72)