Impact Factor 2.486

The Frontiers in Neuroscience journal series is the 1st most cited in Neurosciences

Original Research ARTICLE

Front. Neurorobot., 08 May 2015 |

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

  • Forschungszentrum Informatik (Centre of Computer Science), Intelligent System and Production Engineering, Interactive Diagnosis and Service Systems, Karlsruhe, Germany

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.


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 multi-molecular 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 c1:  leg 1 is loosely bound to the top of leg 2,

which is tightly attached to the substrate.

Leg-over-leg state c2:  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.


Figure 1. Representation of the two leg states a and b and the two leg-over-leg states c1 and c2. 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).

The two legs are attached to the substrate in two different connection states:

Connection state 1: leg 1 is tightly attached to position pk of the tubulin substrate and leg 2 is weakly attached to site pk + 1.

Connection state 2: leg 2 is tightly attached to site pk and leg 1 is weakly attached to position pk + 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 (1c1, 2c1), and is maintained by step 2, which produces the leg states (2a, 1b). In step 3, the leg-over-leg state (1c2, 2c2) is achieved, and finally, in step 4, the states (1a, 2b) are again established.


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 yes. The tightly bound states are marked by striped circles yes. The first position is fixed by k = 1; the initial state is (1a, 2b).

The following notation is employed for the creation (or annihilation) operators:

Molecular robot r: rleg1 statej positionk; leg2 statej positionk; e.g., rl1ap1; l2bp2.

Substrate s: sstatejpositionk; e.g., sgp1.

Table 1 subsumes the different individual motion patterns of the two legs of r.


Table 1. Representation of the motion patterns of a “leg-over-leg” walking process of a molecular robot r.

The energy transfer process (sr) 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) Bpk and Bpk are applied in this walking phase. The operator transfers during the four steps are summarized as follows (readout from Table 1):

This transfer list of operators delivers the definition of the interaction Hamiltonian Hint and finally the elaboration of the resulting equations of motion.

Interaction Hamiltonian

The local interaction Hamiltonian Hint at positions to p1 to p4 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, g1 and g2, are introduced, where the parameter g1 describes the uneven steps and g2 denotes the even steps.

Hint= g1rl1a p1; l2 b p2se p1sg p2rl1c1 p2; l2 c1p2sg p1se p2 Bp2            +g1sg p1se p2rl1c1 p2; l2 c1 p2  se p1sg p2rl1a p1 ; l2 b p2Bp2            +g2rl1c1 p2; l2 c1 p2sg p1sg p3rl1b p3; l2 a p2Bp2            +g2rl1b p3; l2 a p2sg p3sg p1rl1c1 p2; l2 c1p2Bp2            +g1rl1b p3; l2 a p2sg p2s e p2se p3s g p3rl1c2 p3; l2 c2 p3  Bp3            +g1rl1c2 p3; l2 c2 p3 se p2sg p2sg p3se p3rl1b p3 ; l2 a p2 Bp3            +g2rl1c2 p3; l2 c2 p3sg p2sg p4rl1a p3; l2 b p4Bp3            +g2rl1a p3; l2 b p4sg p4sg p2rl1c2 p3; l2 c2 p3Bp3.    (3.1)

By the replacements p1pk, p2pk + 1, etc., Hint is transferred into Hint, k and has to be summed up by k = 1, 2, …. All operators that are quoted here are Bose operators. Generally, we should involve both Fermi operators (anti-commutation rule) and Bose operators (commutation rule) for the definition of Hint and the resulting calculations of the Heisenberg equations of motions. Here, we use only Bose operators because in this case aggregations of molecules are allowed since the Pauli Exclusion Principle does not have to be applied. We neglect Fermi operators that depict supplementary interactions of Fermions with Bosons.

Endowed with the Hamiltonian Hint, the corresponding full equations of motion have been calculated and are presented in the next sub-chapter.

Full Heisenberg Equations of Motion

The full set of Heisenberg equations of motion is completed by 7 damping constants, γab, etc., 7 fluctuating forces, Fab, etc., and two coupling constants, g1, g2, and is defined by the following expressions, representing coupled nonlinear, delayed, and complex operator equations.

r˙l1apk;l2bpk+1= ig1[sepksgpkrl1c1 pk+1;l2 c1pk+1sgpk+1sepk+1 Bpk+1]                            +ig2[sgpk1sgpk+1rl1c2 pk; l2 c2 pkB pk]                            γabrl1apk; l2bpk+1+ Fab.    (4.1)
r˙l1bpk+1; l2apk= ig1[se pksg pkrl1c2 pk+1; l2c2pk+1sgpk+1sepk+1B pk+1]                             +ig2 [sgpk1sg pk+1 rl1c1 pk; l2 c1 pkB pk]                             γbarl1b pk+1; l2apk+ Fab.    (4.2)
r˙l1c1 pk; l2 c1 pk= ig1[sepk1sgpk1rl1apk1; l2 b pksgpkse pkB pk]                              +ig2[sgpk1sgpk+1rl1b pk+1; l2a pkB pk]                              +γc1rl1c1pk; l2c1pk+ Fc1.    (4.3)
r˙l1c2 pk; l2 c2 pk= ig1[sgpk1sepk1rl1b pk; l2 a pk1sepksgpkB pk]                                 +ig2[rl1a pk; l2 b pk+1sgpk1sgpk+1B pk]                                 +γc2rl1c2 pk; l2 c2pk+ Fc2.    (4.4)
s˙gpk= ig1[rl1a pk; l2 b pk+1 se pksg pk+1sepk+1rl1 c1 pk+1; l2c1 pk+1           B pk+1]              +ig1[rl1c1 pk; l2 c1 pk sepk1sg pk1se pkrl1 a pk1; l2b pk        B pk]              +ig2[rl1b pk+2; l2 a pk+1 sgpk+2rl1 c1 pk+1; l2c1 pk+1 B pk+1]              +ig2[rl1 c1pk1; l2 c1 pk1 sgpk2rl1b pk; l2a pk1 B pk1]              +ig1[rl1b pk+1; l2 a pk se pksepk+1 sg pk+1rl1 c2 pk+1; l2c2 pk+1     B pk+1]              +ig1[rl1c2 pk; l2 c2 pk sepk1sg pk1se pkrl1b pk; l2a pk1 B pk]              +ig2[rl1 c2pk1; l2 c2 pk1 sgpk2rl1a pk1; l2b pk B pk1]              +ig2[rl1apk+1; l2bpk+2 sgpk+2rl1c2 pk+1; l2c2 pk+1B pk+1 ]              γg sgpk+Fg.    (4.5)
s˙epk= ig1[rl1a pk1; l2 b pk sgpk1se pk1sgpkrl1 c1 pk; l2c1 pk B pk]             +ig1[rl1c1 pk+1; l2 c1 pk+1 sgpksgpk+1se pk+1rl1a pk; l2b pk+1                        B pk+1]             +ig1[rl1b pk; l2 a pk1sgpk1 se pk1sgpkrl1 c2 pk; l2c2 pk B pk]             +ig1[rl1c2 pk+1; l2 c2 pk+1 sg pksgpk+1se pk+1rl1b pk; l2a pk1                        B pk+1] γe se pk+Fe.    (4.6)
B˙ pk= ig1[{rl1a pk1; l2 b pkrl1c1 pk; l2 c1 pk               +rl1b pk; l2 a pk1rl1c2 pk; l2 c2 pk} se pk1sg pk1sg pkse pk]               +ig2[{rl1b pk+1; l2 a pkrl1c1 pk; l2 c1 pk               +rl1a pk; l2 b pk+1rl1c2 pk; l2 c2 pk} sg pk1sgpk+1]               γB B pk+FB.    (4.7)

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.


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):

|ϕ(t0)=rl1a p1; l2 bp2se p1sg p2|ϕ0,    (5.1)

where | ϕ0 〉 delineates the vacuum state. The initial conditions of all relevant states are:

rl1a p1; l2bp2(t0)=1,rl1c1p2; l2c1p2(t0)=0.    (5.2)
sep1(t0)=1,sgp1(t0)=0,sgp2(t0)=1,sep2(t0)=0. Bp2(t0)=0.    (5.3)

We set rl1c2p1; l2c2 p1 = 0 in Equation (4.1) because this state is not present. The resulting modified Equation (5.1) now reads:

r˙l1a p1; l2 b p2= ig1[sep1sgp1rl1c1p2; l2c1p2sgp2sep2Bp2]                                γabrl1a p1; l2b p2.    (5.4)

In a similar way, going through all remaining equations the following reduced set of equations is obtained:

r˙l1c1p2; l2c1p2= ig1[sep1sgp1rl1a p1; l2bp2sgp2sep2Bp2]                                γc1rl1c1 p2; l2c1p2    (5.5)
s˙gp1=  ig1[rl1ap1;l2b p2sep1sgp2se p2rl1c1p2;l2c1p2Bp2]γgsgp1,    (5.6)
s˙gp2= ig1[rl1c1 p2; l2 c1 p2sep1sg p1se p2rl1 a p1; l2b p2Bp2]γgsgp2,    (5.7)
s˙ep1= ig1[ rl1c1p2;l2c1p2sgp1sgp2 sep2rl1ap1;l2bp2Bp2]γesep1    (5.8)
s˙ep2= ig1[rl1a p1; l2b p2 sgp1sep1sgp2rl1c1 p2; l2c1 p2 Bp2]γesep2.    (5.9)
B˙p2= ig1[rl1ap1; l2 bp2 sg p1sgp2se p1se p2rl1c1 p2; l2 c1p2]γBBp2.    (5.10)

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 = γc1 = γc2, and neglecting for the moment the damping of the B operator: γB = 0), delivers a consistent, periodic solution for a walking biped molecular robot:

rl1ap1;l2bp2=cos(f)eγt,rl1c1p2;l2c1p2=sin(f)eγt.    (5.11)
sep1=cos(g),sgp1=isin(g),sep2=cos(h),sgp2=sin(h).    (5.12)
Bp2=b, b .    (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:

f˙= g14sin(2g)sin(2h)b.    (5.14)
g˙= g14sin(2f)eγtsin(2h)b.    (5.15)
h˙= g14sin(2f)eγtsin(2g)b.    (5.16)
b˙= g12g˙ sin(2g).    (5.17)

The last equation can be directly solved by the technique of separation of variables:

b(t)= g14[cos(2g(t))cos(2g(t0))]+b(t0).    (5.18)

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):


Figure 3. Three-dimensional limit cycle defined by the three variables b, g, and h. The initial point is marked by a cross.

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.


Figure 4. Temporal dependence of the real parts of all variables of step 1, with γ = 0. The algebraic symbols are:  abp1=rl1ap1;l2bp2, c1p2=rl1c1p2;l2c1p2, gp1=sgp1, gp2=sgp2, ep1=sep1, ep2=sep2, Bp2=Bp2. The coupling constants are g1 = g2 = 0.1.

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 solution—it 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).


Figure 5. Damped temporal curves representing the real parts of the following operators of the first step: abp1=rl1ap1;l2bp2,c1p2=rl1c1p2;l2c1p2,gp1=sgp1,gp2=sgp2,ep1=sep1,ep2=sep2,Bp2=Bp2. The value of the damping constant is set to γ = 0.0001. The coupling constants are g1 = g2 = 0.1.

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:

|ϕ(t1))=rl1c1p2;l2c1p2se p2Bp2|ϕ0,    (5.19)

The initial conditions are:

rl1c1p2; l2c1 p2(t1)=1,rl1b p3; l2a p2(t1)=0.    (5.20)
sgp1(t1)=1,sep2(t1)=1,sgp3(t1)=0.Bp2(t1)=1.    (5.21)

The Heisenberg equations of motion of the expectation values of step 2 are:

r˙l1c1 p2; l2 c1 p2= ig2[sgp1sgp3rl1b p3; l2a p2Bp2]γc1rl1c1 p2; l2 c1 p2.    (5.22)
r˙l1b p3; l2 ap2= ig2[sgp1sg p3rl1c1 p2; l2 c1p2 Bp2]γbarl1b p3; l2ap2    (5.23)
s˙gp1=  ig2[rl1bp3; l2ap2sgp3rl1 c1p2; l2c1p2Bp2]γgsgp1.    (5.24)
s˙gp3= ig2[rl1 c1p2; l2 c1 p2 sgp1rl1b p3; l2a p2Bp2]γgsgp3.    (5.25)
B˙p2=  ig2[rl1b p3; l2 ap2sgp1sgp3rl1c1 p2; l2 c1p2]γBBp2.    (5.26)

Here, the damping process is turned on again at the beginning, with γ = 0.0005. Figure 6 shows the expected “ringing” effects now represented in phase portraits.


Figure 6. Phase portraits of the real parts of the second step of the following operators: bap3=rl1ap3;l2bp2, c1p2=rl1c1p2;l2c1p2, gp1=sgp1, gp3=sgp3, Bp2=Bp2. The value of the damping constant is γ = 0.0005. All trajectories vary after their starting points (marked by a cross in the first figure) until they end up at the fixed point 0. The coupling constants are g1 = g2 = 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:

|ϕ(t2))=rl1b p3; l2a p2se p2sg p3|ϕ0,    (5.27)

In this case, the initial conditions of all states are:

rl1b p3; l2 a p2(t2)=1,rl1c2p3; l2c2 p3(t2)=0,    (5.28)
sep2(t2)=1,sgp2(t2)=0,sep3(t2)=0,sgp3(t2)=1.    (5.29)
Bp3(t2)=0.    (5.30)

The Heisenberg equations of motion of the expectation values of step 3 are:

r˙l1b p3; l2ap2= ig1[se p2sgp2rl1c2 p3; l2 c2 p3sgp3sep3 Bp3]                             γbarl1b p3; l2ap2    (5.31)
r˙l1c2 p3; l2c2 p3= ig1[sgp2 sep2rl1b p3; l2 a p2sep3sgp3Bp3]                                  γc2rl1c2 p3; l2 c2 p3    (5.32)
s˙gp2= ig1[rl1bp3;l2 ap2sep2sep3 sgp3rl1 c2 p3;l2c2 p3 Bp3]γg sgp2.    (5.33)
s˙gp3= ig1[rl1c2 p3; l2 c2 p3sep2sgp2sep3rl1b p3; l2a p2 Bp3]  γg sgp3.    (5.34)
s˙ep2= ig1[rl1c2p3; l2c2p3 sgp2sgp3sep3rl1b p3; l2a p2Bp3] γe sep2.    (5.35)
s˙ep3= ig1[rl1b p3; l2 a p2sgp2 sep2sgp3rl1 c2 p3; l2c2 p3 Bp3] γe sep3.    (5.36)
B˙p3= ig1[ rl1b p3; l2 a p2rl1c2 p3; l2 c2 p3 sep2 sgp2 sg p3se p3] γBBp3.    (5.37)

Attention is drawn to the symmetry that exists between all equations of motion of the first and third steps. Further, the initial condition for sgp2 in step 3 is sgp2(t2) = 0, while in step 1 sgp2(t0) = 1 was valid. Similarly, the value sep2(t2) = 1 is assumed in step 3, whereas sep2(t0) = 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:

|ϕ(t3))=rl1c2p3; l2c2p3sep3Bp3|ϕ0,    (5.38)

where | ϕ0 〉 defines the vacuum state. In this case, the initial conditions of the two robot states are:

rl1b p3; l2 a p2(t3)=0,rl1c2p3; l2c2 p3(t3)=1,.    (5.39)

The initial states of the substrate at the three positions p2, p3 and p4 are:

sgp2(t3)=1,sep3(t3)=1,sgp4(t3)=0.    (5.40)

The heat-bath operator satisfies the initial condition Bp3(t3) = 1.

The relevant Heisenberg equations of motion of step 4 are:

r˙l1c2 p3; l2 c2p3= ig2[sgp2sgp4rl1a p3; l2bp4Bp3] γc2rl1c2 p3; l2 c2 p3.    (5.41)
r˙l1a p3; l2b p4= ig2[sgp2sgp4rl1c2 p3; l2 c2 p3Bp3] γabrl1ap3; l2bp4.    (5.42)
s˙gp2=  ig2[rl1a p3; l2 b p4 sgp4rl1c2 p3; l2c2 p3Bp3]  γgsgp2,    (5.43)
s˙gp4= ig2[rl1c2 p3; l2 c2 p3sgp2rl1ap3; l2b p4Bp3] γgsgp4    (5.44)
B˙p3= ig2[rl1a p3; l2 b p4rl1c2 p3; l2 c2 p3sgp2sgp4] γBBp3.    (5.45)

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 rl1c1p2;l2c1p2 slightly more oscillations than for rl1ap1;l2bp2 in step 2.


Figure 7. Damped temporal curves representing the real parts of the four molecular robot operators that are consecutively performed by the steps 1–4. The following algebraic equivalents for the selected operators are introduced: abp1=rl1ap1;l2bp2, c1p2=rl1c1p2;l2c1p2, bap3=rl1bp3;l2ap2, c2p3=rl1c2p3;l2c2p3. The value of the damping constant is set to γ = 0.0005. The coupling constants are g1 = g2 = 0.1.

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 sgp3 is different in step 2, and during step 2 (see Figure 2), the ground state sgp3 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.


Figure 8. Damped temporal curves representing the real parts of three selected substrate operators (ground states) that are consecutively performed by steps 1–4. The following algebraic equivalents for the selected operators are introduced: gp2 = sgp2, gp3 = sgp3, gp4 = sgp4. The value of the damping constant is set to γ = 0.0005. The coupling constants are g1 = g2 = 0.1.

From a purely mathematical point of view, we have to compare the two Equations (5.25) and (5.34) for sgp3 in the second and third steps. In the second step, there is a strong mutual dependence between sgp3 and sgp1. In Equation (5.34), we mainly have to consider the dependence of sgp3 and sep3. The operator product sep2sgp2 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 Bp2 and Bp3. 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, Bp2 starts with Bp2(t0) = 0 and ends with Bp2(t1) = 1. In step 2, it starts with this initial value, Bp2(t2) = 1, and ends permanently with Bp2(t2) = 0. Such interplay is valid in the same manner for Bp3(t3). 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.


Figure 9. Damped temporal curves showing the overlay of the real (gray) and imaginary (red) parts of the operators Bp2 and Bp3, which are consecutively outlined through steps 1–4. The algebraic equivalents of the selected operators are: Bp2 = Bp2, Bp3 = Bp3. The damping constant is set to γB = 0.001. The coupling constants are g1 = g2 = 0.1.

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 Bp2 and Bp3. The contours of each step can be considered as an envelope around the oscillating solutions.


Figure 10. Damped temporal curves showing the overlay of the real (gray) and imaginary (red) parts of the operators Bp2 and Bp3, which are consecutively outlined through steps 1–4. The algebraic equivalents of the selected operators are Bp2 = Bp2 and Bp3 = Bp3. The value of the damping constant is γB = 0.0001. The coupling constants are g1 = g2 = 0.1.

The next two sections focus on wave solutions. Firstly, we combine modes that synchronize within a lane and between lanes (similar to light modes or laser modes, e.g., Haken, 1970). Secondly, we present a solution that constitutes a running wave.

Synchronized Wave-based Motion

Coherent Motion on Parallel Substrate Lanes

It is assumed that there are M different lanes, and on each lane up to L positions (lane length) are available. The different walking lanes of robots are distinguished by the index m = 1, …, M, and the position of a molecular robot on lane m is denoted by pm = 1, …, L. The B fields are also coupled in a lane and across different lanes. Therefore, we introduce a new two-dimensional vector l by combining both parameters l = (m, pm). This combination modifies the Equations (4.1–4.7) by an additional index, as e.g., demonstrated by the reformulation of Equation (4.1):

r˙l1al; l2bl+(0,1)= ig1[sel sglrl1c1l+(0,1); l2 c1l+(0,1)sgl+(0,1) sel+(0,1)                                Bl+(0,1)]+ig2[sgl(0,1) sgl+(0,1) rl1c2 l; l2 c2lBl]                                γab,l rl1al; l2bl+(0,1)+ Fab, l    (5.46)

Further, we define another two-dimensional vector k=(kj=2π njL, kj+1=2π nj+1L), nj, nj+1=,1, 0, 1,  which defines a two-dimensional wave number vector. This second vector has been introduced since the operator Bl will be developed by a plane wave approximation:

Bl=1L kB˜keik·l.    (5.47)

The synchronization activities of the two heat-bath operators Bl, Bl on the molecular robots are described as external signal field modes that are coupled by the constant Jll.

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 pm of l. The use of the coupling constant Jll for  B˜k requests a summation along a lane and between lanes:

Bl= 1LkB˜˙keik·l      = ig1[m{r(l1apm1; l2bpm), m r(l1c1pm; l2 c1pm),m            + r(l1bpm; l2 apm1), m r(l1c2 pm; l2 c2pm),m}            × sepm1,msgpm1,m sgpm,msepm,m]            + ig2[m{r(l1b pm+1; l2 apm),m r(l1c1 pm; l2 c1pm),m            + r(l1a pm; l2bpm+1),m r(l1c2 pm; l2c2pm),m} sgpm1sgpm+1,m]             γB, l1L kB˜keik·l+FB,l+1LlJll kB˜keik·l.    (5.48)

We multiply both sides of (5.48) with 1Lleik·l and use the orthogonal relation 1Llei(k·lk·l)=δkk:

B˜˙k= ig1[l,m1Leik·l{r(l1apm1; l2bpm),mr(l1c1pm;l2c1pm),m             +  r(l1bpm; l2 apm1), m r(l1c2 pm; l2 c2pm),m}             × sepm1,msgpm1,m sgpm,msepm,m]             +ig2[l,m1Leik·l{r(l1b pm+1;l2apm),mr(l1c1 pm;l2c1pm),m             +r(l1a pm;l2 bpm+1),mr(l1c2 pm;l2c2pm),m}sgpm1, msgpm+1,m]              γB, l1Lleik·l kB˜keik·l+1Lleik·lFB,l             + 1Lleik·llJll kB˜keik·l.    (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 = ll′. Further, it is supposed that the coupling constant J depends only on the distance Jll = Jd and the lane length L is defined by = LxLd. It follows that:

llJd1Leik·lik·l= δkk1Lddei(k+k)·d2Jd                                           = 1Lddeik·dJd=ωk    (5.50)

where ωk has the dimension of a frequency. We continue with the simplification process and set

Γkk=1Ll ei(kk)·l γB, l=δkkγB,    (5.51)

if all damping constants are set equal to

γB,l= γB.    (5.52)
F˜k= 1Lleik·lFl.    (5.53)

Finally, after all these simplifications the resulting formula is obtained:

B˜˙k= ig11L[l,m eik·l{r(l1apm1; l2 bpm),  m r(l1c1 pm; l2 c1pm),m            +  r(l1bl; l2 apm1), m r(l1c2 pm; l2 c2pm),m}            × sepm1,msgpm1,m sgpm,msepm,m]            +ig21L[l,meik·l{r(l1bpm+1;l2apm),mr(l1c1 pm;l2 c1pm),m            +r(l1a pm;l2bpm+1),mr(l1c2pm;l2c2pm),m}sgpm1, m sgpm+1,m]             γ B˜k+F˜k+ωkB˜k.    (5.54)

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.

B˜˙0= ig11L[m {r(l1apm1; l2 bpm),  m r(l1c1 pm; l2 c1pm),m             +  r(l1bpm; l2 apm1), m r(l1c2 pm; l2 c2pm),m}             ×sepm1,msgpm1,m sgpm,msepm,m]             + ig21L[ m{r(l1bpm+1; l2apm),m r(l1c1 pm; l2 c1pm),m             +r(l1a pm;l2 bpm+1),mr(l1c2 pm;l2 c2pm),m}sgpm1,m sgpm+1,m]              γB˜0+F˜0+ω0B˜0.    (5.55)

The last term, ω0B˜0, can be deleted in the interaction representation and the expectation value of F˜0 can be ignored. These two modifications lead to the following equation:

B˜˙0= ig11L[m {r(l1apm1; l2 bpm),  m r(l1c1 pm; l2 c1pm),m           +  r(l1bpm; l2 apm1), m r(l1c2 pm; l2 c2pm),m}           × sepm1,msgpm1,m sgpm,msepm,m]           +ig21L[ m{r(l1b pm+1; l2apm),m r(l1c1 pm; l2 c1pm),m           +r(l1a l; l2 bpm+1),m r(l1c2 pm; l2 c2pm),m}sgpm1,m sgpm+1,m]            γB˜0.    (5.56)

Comparing the result of Equation (5.56) with the original expression (4.7), two main differences are immediately observable. At first, the two coupling constants, g1 and g2, have to be replaced by g11L and g21L. Second, there is a summation over all lanes m = 1, …, M. The elaborate calculations of this subchapter show that the synchronization of motion can easily be accomplished (at least for B˜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 pk = k, pk + 1 = k + 1, … Under these assumptions, the approach that guides us to a running wave solution is formulated as follows:

rl1a pk; l2b pk+1= RabeiK/(k+1),    (5.57)
rl1c1 pk+1; l2c1pk+1= Rc1eiK/(k+1),    (5.58)
sgpk= SgeiK/k,    (5.59)
sepk= SeeiK/k,    (5.60)
sgpk+1= S˜geiK/(k+1),    (5.61)
sepk+1= S˜eeiK/(k+1),    (5.62)
Bpk+1= BbeiK/(k+1),    (5.63)

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:

R˙ab= ig1SeSg S˜gS˜eRc1BbeiK/(k+1)γRab        = ig1D Rc1BbeiK/(k+1) γRab,    (5.64)
R˙c1= ig1SeSgS˜gS˜eRabBbeiK/(k+1) γRc1       = ig1DRabBbeiK/(k+1) γRc1.    (5.65)

In doing so, we introduced the following abbreviations:

D1=SeSg=SgSe,D2=S˜gS˜e=S˜eS˜g,D=D1D2= D2D1.

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.


Figure 11. Undamped temporal curves, showing the real parts of the operator products D1, D1 (abridged to D1k), D2, D2 (abridged to D2k), D, D (abridged to Dk) and Bp2Bp2, consecutively outlined for step 1. The damping constant is γ = 0, the coupling constants are constant: g1 = g2 = 0.1.

In the next step, we reformulate the Equations (5.6–5.10) by including the D–terms:

s˙g= ig1SeD2RabRc1BbeiK/(k+1)γSg,    (5.66)
s˙e= ig1SgD2Rc1Rab Bb eiK/(k+1)γSe,    (5.67)
S˜˙g=  ig1S˜eD1Rc1RabBbeiK/(k+1)γS˜g    (5.68)
S˜˙e=  ig1S˜gD1RabRc1BbeiK/(k+1)γS˜g,    (5.69)
B˙b= ig1DRabRc1eiK/(k+1)γBBb.    (5.70)

The calculations are continued by substituting the three operators Rab, Rc1 and Bb by

Rab=ρabeγt,Rc1=ρc1eγt,Bb=ρbeγBt.    (5.71)

The derivatives of the first two expressions read as follows:

R˙ab=ρ˙abeγtγRab,R˙c1=ρ˙c1eγtγRc1.    (5.72)

Insertion of r˙ab in (5.64) provides the expression

ρ˙ab=ig1DeγBt(ρc1ρb)eiK/(k+1)ig1GeiφeγBt(ρc1ρb),    (5.73)

where φ˙ = 0 and both φ and G are real.

By inserting R˙c1 in (5.65), the following expression is obtained:

ρ˙c1=ig1DeγBt(ρabρb)eiK/(k+1)ig1G2eiφeγBt((ρabρb).    (5.74)

It follows, by insertion of Bb, given by Equation (5.71) into expression (5.70), that

ρ˙b= ig1De(2γγB)t(ρabρc1)eiK/(k+1)       ig1G2eiφe(2γγB)t(ρabρc1).    (5.75)

To reproduce the r.h.s. of these three Equations (5.73–5.75), the following assumption is pursued:

ρab=ρ0eF(t),ρc1=ρ0eF(t)2eiφ,ρb=ρ0eF(t)2e2iφ,    (5.76)

where F(t)=0tG(τ)dτ, ρ0=const..

With this ansatz the derivatives of the three operators ρab, ρc1, ρb can be cast in the following form:

ρ˙ab=iGρab=iG1ρ0(ρc1ρb)eiφ,    (5.77)
ρ˙c1=iG2ρc1=iG2 1ρ0(ρabρb)eiφ,    (5.78)
ρ˙b=iG2 ρb=iG21ρ0(ρabρc1)eiφ.    (5.79)

From expression (5.83), the value of G is concluded as:

ρ0Re(D)eiK/(k+1)eiφ=G.    (5.80)

Both sides of this expression should be real, therefore the exponents must be zero: φ = − K/(k + 1) and ρ0Re(D) = G. A second solution is formed by the same φ but with ρ0Re(D)=G2. This last result comes from the two Equations (5.74) and (5.75).

Due to the determination of the two variables φ and K, running wave solutions of the molecular “leg-over-leg” walking can be expressed for the three operators rl1a p1; l2b p2, rl1c1 p2; l2c1p2, Bp2 in the conventional form (k = 1, 2,..; G = ρ0Re(D) = const.):

rl1a pk; l2b pk+1= ei(Gtφ(k+1)) ρ0eγt,    (5.81)
rl1c1 pk+1; l2c1pk+1= ei(Gtφ(k+1)) ρ0eγt,    (5.82)
Bpk+1=ei(Gtφ(k+1)) ρ0eγBt.    (5.83)

Such a solution is not surprising because it is well known that a composition of different modes (more than only the ground mode k = 0) can generate running waves (see Section Coherent Motion on Parallel Substrate Lanes).


The central goal of this contribution was to support the hypothesis that quantum mechanical effects in molecular biology—especially in human brains—cannot be neglected. This assumption is mainly justified by the size of the interacting objects (nano size) and by the experimental results showing that even greater molecules up to about 100 atoms can demonstrate quantum behavior, e.g., in double-slit experiments (Haken and Levi, 2012). The nano size argument is further supported by the established experimental technology to construct molecular machines and walkers from DNA.

The most striking results to emerge from the produced solutions is that neural processes like the motion of molecular robots (kinesin and dynein) in axons and dendrites can be modeled by an approach that is particle-based (impact solution) or purely wave-oriented (synchronized modes solution and running wave solution). In addition, wave-based solutions also seem to be very relevant for the interactions between different neural layers.

We are aware that our predictions suffer from limitations and have to be considered as a first, basic approach due to the following three reasons. Firstly, the number of the many parameters (in total 16) was reduced to four (two coupling constants and two damping constants). Secondly, there are still not enough experimental data to fix the values of the last mentioned four parameters or even the great set of all relevant available biological data. Thirdly, a greater number of experimentally approved parameters could lead to a higher generalization of the achieved results.

A first extension of our approach could be the integration of tunneling effects, the inclusion of fluctuation forces, and the consideration of the cargo transport. All these additional considerations would greatly increase the predictive power of the outlined model. In this way, it would be feasible to describe e.g., the walk of a molecular robot along the DNA and the transcription of a nucleotide into RNA, whereby for each step the effective potential barrier must be tunneled.

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


The research leading to these results has received funding from the European Union Seventh Frame Program (FP7/2007-2013) under grant agreement no. 604102 (Human Brain Project, Neurorobotics Platform, SP10). We express also our gratitude to Dr. B. Schenke for his numerical calculation support of this contribution.


Alberts, B., Johnson, A., Lewis, J., Raff, M., Roberts, K., Walter, P., et al. (2008). Molecular Biology of the Cell, 5th Edn. New York, NY: Garland Science.

Google Scholar

Balzani, V., Venturi, M., and Credi, A. (2008). Molecular Devices and Machines: Concepts and Perspectives for the Nanoworld, 2nd Edn. Weinheim: Wiley–VCH.

Google Scholar

Fukuda, Y. T., Hasegewa, Y., Sekiyama, K., and Aoyama, T. (2012). Multi-locomotion Robotic Systems, a New Concept of Bio-inspired Robotics. Berlin: Springer.

Google Scholar

Haken, H. (1970). Laser Theory. Encyclopedia of Physics. Berlin: Springer.

Haken, H., and Levi, P. (2012). Synergetic Agents—Classical and Quantum. From Multi–robot Systems to Molecular Robots. Weinheim: Wiley–VCH.

Holmes, P. H., Lumley, J. L., Berkooz, G., and Rowley, C. W. (2012). Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge Monograph on Mechanics. Cambridge: Cambridge University Press.

Google Scholar

Levi, P., and Haken, H. (2010). “Towards a synergetic quantum field theory for evolutionary symbiotic multi-robotics,” in Symbiotic Multi–robot Organisms, eds P. Levi and S. Kernbach (Berlin: Springer), 25–54.

Google Scholar

Lund, K., Manzo, A. J., Dabby, N., Michelotti, N., Johnson-Buck, A., Nangreave, J., et al. (2010). Molecular robots guided by prescriptive landscapes. Nature 465, 206–210. doi: 10.1038/nature09012

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Murata, S., Konagaya, A., Kobayashi, S., Saito, H., and Hagiya, M. (2013). Molecular Robotics: A New Paradigm of Artifacts, New Generation Computing, Vol. 31. Berlin: Ohmsha Ltd. and Springer.

Google Scholar

Omabegho, T., Sha, R., and Seeman, N. C. (2009). A bipedal DNA Brownian motor with coordinated legs. Science 324, 67–71. doi: 10.1126/science.1170336

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Roux, B. (ed.). (2011). Molecular Machines. River Edge, NJ: World Scientific.

Google Scholar

Seeman, N. C. (2014). “Molecular machines made from DNA, to be published,” in Proceedings of the 20th International Conference on DNA Computing and Molecular Programming. (Kyoto: Kyoto University).

Sherman, W. B., and Seeman, N. C. (2004). A precisely controlled DNA biped walking device. Nano Lett. 4, 1203–1207. doi: 10.1021/nl049527q

CrossRef Full Text | Google Scholar

Shin, J. S., and Pierce, N. A. (2004). A synthetic DNA walker for molecular transport. J. Am. Chem. Soc. 126, 10834–10835. doi: 10.1021/ja047543j

PubMed Abstract | Full Text | CrossRef Full Text | Google Scholar

Keywords: molecular robotics, biped walking on microtubules, motion of kinesin-1 and dynein alog axons and dendrites, particles and wave solutions, neuro-robotics

Citation: Levi P (2015) Molecular quantum robotics: particle and wave solutions, illustrated by “leg-over-leg” walking along microtubules. Front. Neurorobot. 9:2. doi: 10.3389/fnbot.2015.00002

Received: 19 September 2014; Accepted: 23 April 2015;
Published: 08 May 2015.

Edited by:

Marc-Oliver Gewaltig, Ecole Polytechnique Federale de Lausanne, Switzerland

Reviewed by:

Hiroaki Wagatsuma, Kyushu Institute of Technology, Japan
Valery E. Karpov, National Research University Higher School of Economics, Russia

Copyright © 2015 Levi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Paul Levi, Forschungszentrum Informatik (Centre of Computer Science), Intelligent System and Production Engineering, Interactive Diagnosis and Service Systems, Haid-und-Neu-Street 10-14, 76131 Karlsruhe, Baden-Württemberg, Germany,