Inﬂuence of boundary conditions on stochastic movement of biological molecular motors

The aim of this work is to show the use of statistical physics in biomechanics, especially for modelling a muscle on molecular level. To approximate stochastic (Markovian) evolution in the spatial variable of muscle conﬁguration, we solve the Fokker-Planck equation with the use of the Wang-Peskin-Elston algorithm. The approach is used for numerical simulation in MATLAB of the system time development with different boundary conditions.


Introduction
In nanoscale world, thermal fluctuations (Brownian motion) has a great influence on the movement of all particles which are involved in a microsystem description.This eternal motion cannot be neglected when modelling molecular processes.Nevertheless, the motion is stochastic and cannot be predicted by solving a dynamic equation of motion.That is why we use a statistical prediction for a description of a myosin II molecular motor, which is responsible for muscle contractions.
There are many statistical approaches to describe stochastic systems out of equilibrium.One of them is based on the Langevin equation [15].It can simulate a particle movement in the stochastic environment.The Langevin equation can be generalised to the Fokker-Planck equation [15].We express the stochasticity of the system dynamics by this equation.The equation is modified to the Master equation via the WPE algorithm [19].We solve the Master equation numerically for various boundary conditions.We use MATLAB for the computations and its solver ode15s.Then we study their influence to a control parameter, which is defined as a sum of probabilities at a given time.
The paper is organised as follows.A brief description of myosin II dynamics, especially from the biological point of view, is in Section 2. In following sections, we describe a mathematical approach (Section 3) which is used for a numerical model (Section 4) and numerical results (Section 5) devoted to the dependence of boundary conditions.

Molecular motors -A biological view
In an organic body, it is necessary to ensure the transport of organelles, cell division, muscle contraction, and other biological processes.Three superfamilies of molecular motors are described for executing such processes [11].
Motors from kinesin and dynein superfamilies are mostly involved to mediate transport along microtubules.Kinesin moves primarily from the minus (−) end to the plus (+) end of the microtubule.Two of 14 kinesin classes are exceptions from these rules.Kinesin-13 can enhance the depolymerization of microtubule ends due to the absence of motor activity.Kinesin-14, on the other hand, has motor activity but it moves in opposite direction to other kinesins -from the (+) end to the (−) end [11].
Motors from dynein superfamily provide transport in the same way as kinesin-14 on microtubules.So, they are called retrograde molecular motors.Dyneins are very large proteins, consisting of two large (> 500 kDa), two intermediate, and two small subunits.The family of dynein-related proteins is not very diverse [11].
The last molecular motors group is myosin superfamily.These motors convert an energy stored in adenosine triphosphate molecules (ATP) into work.Every motor from this family contends from six polypeptide subunits -two identical high-molecular-weight polypeptides (Myosin heavy chains) and four light chains.Myosin heavy chains are created by a globular head, a flexible neck and a long tail.
Also in this superfamily, some very different classes exist.For example, myosin I is involved in endocytosis (transport organelles through cell membrane [17]).Another different representative of the myosin superfamily is myosin V which is quite similar in its functionality to kinesins.It also transports a cargo along a filament.In this case, it is an actin filament.
The most common member of the myosin superfamily is myosin II, which was also discovered as the first one, so it is known shortly as myosin [11].Myosin II molecular motor is responsible for a muscle contraction.It is caused by a movement of the myosin II on an actin filament and their connection.The myosin head is attached to the actin filament during consumption energy from ATP hydrolysis.It was shown, one molecule of ATP provides enough energy for several myosin head states [7,8].Myosin II produces 3-4 pN per its cycle [5].
Myosins and actins create together myofilaments.They connect themselves into a myofibril.Packs of myofibrils formulate a group of muscles fibres then.These groups clump into fasciculus and several these fasciculus create a whole muscle [14].Even in one cell, there are hundreds of myosin II motors.Thus, there is a myriad of myosins in every muscle.
The complex myosin movement is usually simplified to a 6-state kinetic model [6].The cycle is created as follows: A molecule of ATP is bounded on the myosin head (state 1).Then the ATP is hydrolyzed to ADP (adenosine diphosphate) and P (phosphorus) -state 2. The hydrolysis provides an energy to bound to actin (state 3), where the atom of phosphorus is released (state 4).The ADP molecule is relaxed from myosin (state 5) and it is replaced by a new ATP molecule (state 6).Then myosin unbinds itself from actin and the cycle is finished.
A little different description to the above mentioned is introduced in [11] where the 5-state myosin cycle is presented.It assumes that during the first state, the myosin head is released from the actin filament and the head binds ATP.In the second state, the ATP is hydrolyzed to ADP and P. It causes the myosin head rotation.Then, the myosin head is bent to the actin filament.In the fourth state, there is the "power-stroke", when the phosphorus atom is released, elastic energy straightens myosin and moves actin filament to the left.In the last state, the ADP molecule is released, new ATP molecule is bounded and the head is released from the actin.The fifth state ends at the same configuration as the first one.So, it can be considered as a 4-state myosin cycle.
The movement of myosin can be more simplified to 3-state cycle [2,3,16] or even to 2-state cycle [1].These models neglected several states of the previous cycle.For example, the 3-state model uses only the second, the third and the fourth states from the 5-state model, which is mentioned above.The 3-state cycle is described in more details in Section 5.

Statistical description of a muscle at the molecular level
Myosins, due to their size, are influenced by their own Brownian motion and also by bombardment of surrounding fluids molecules [18].For this reason, the most suitable description of muscle at the molecular level is a statistical one.
A molecular motor has many degrees of freedom, but there is only one connected to undirectional motion of the motor.The others can be approximated by a mean field potential because their timescales are much shorter than the time set by the average motor velocity.Another reason, to create a mathematical model only in one spatial dimension, is to make the best connection with an experiment, where the motor is observed along one spatial dimension [18].
The simplest way to mathematically describe the molecular motor through its probability density is by the Smoluchowski equation where ρ is the normalised concentration of the motor's mental copies, therefore can be interpreted as the probability density of the motor presence, D is the diffusion coefficient, ϕ is defined as the potential of an external force F ext [4] by the classic relation The external force does not include only collisions caused by the Brownian motion.It is made by charged surrounding particles in an aqueous solution.
We use one diffusion coefficient for simplicity.It is connected with the Brownian motion by the Einstein relation where ξ is the average friction (drag) coefficient [10].Because molecular motors work under an isothermal condition, the temperature T is constant [9].Thus, the diffusion has a meaning of a (Gaussian) random walk.If a molecular motor is driven by alternating its exposure to multiple energy profiles, the driving equation needs to be coupled with a term describing transitions between these potentials.In the case of the Smoluchowski equation, the equation ( 1) is modified into a system of equations where k ij are transition rates between states.They have to fulfil the condition [18] M j=1 k ij (x) = 0 for ∀i. (5) A little more complex equation (used, for example, in [1]) is the Fokker-Planck one which is changed with respect to the Smoluchowski equation in the term with potential.The term , where V i (x) is the active potential energy produced by a chemical reaction.The prime denotes spatial derivation with respect to x, which changes the energy to the corresponding force.The parameter F Load is external load force of the myosin motor.In case of a transition between potentials, the equation has the form of The approach via the Smoluchowski equation or the Fokker-Planck equation leads to a ratchet model.A large number of ratchet models can be found in [13] for different problems.
In [1], there is a 2-state myosin model connected with potentials which correspond to an on-off ratchet.The model is very simple.The myosin head is detached (the corresponding potential is zero), or attached (the potential is modelled by a Fourier series that guarantees the smoothness of the solution and the periodicity of the potential).There is nothing between these states.In this paper, the ratchet model is enhanced to a 3-state ratchet, i.e. a 3-state model of myosin head (see Section 5).
Other molecular motors approximation are through a direct description of the Brownian force f B .It is a very complex approach and that requires a random number generator [4] to simulate a statistical behavior of the thermal motion (the white noise [13]).The simplest governing equation is called the Langevin equation.It can be written as follows In fact, it is an overdamped Newton law of motion in 1D, where the mass m is typically about 10 −21 kg, and the drag coefficient ξ is about 10 −7 pNs/nm [4].

A specific method for the numerical solution of transport equations
The Smoluchowski equation ( 4) and the Fokker-Planck equation ( 6) are partial differential equations of the second order in space and of the first order in time.So, it is necessary to discretize them for their numerical solution.For time derivation, it is possible to use the wellknown method of Runge-Kutta (or similar).In the space dimension, it is possible to use one of the standard methods (for example, the Euler method), but it does not preserve the detail balance [18].It leads to an additional numerical error in a stationary state.An algorithm, which preserves the detail balance, is well described in [19].The algorithm is called WPE according to its authors.Basically, the problem for simulation of a probability density is replaced by a problem of a probability of jumps between nodes.Then, the solved equation has the form of the master equation where p n is the probability of finding the motor at the node x n at time t and it satisfies the relationship for the space step Δx.
Parameters F n+1/2 and B n+1/2 are the forward and the backward transition rates between nodes x n and x n+1 , respectively (n = 1, . . ., N).These parameters are determined by a stationary solution of the original equation ( 1) or (6).For the Fokker-Planck equation ( 6), it can be written as follows and where Δφ n+1/2 is defined as and φ(x) = V (x) − xF Load is called the effective potential [13].In case of absorbing boundary conditions (see Fig. 1), every molecular motor particle vanishes at the boundary.Thus, the transition rates have to change.On the left end of the domain (x = 0) the backward rate B abs has the form of On the right end (x = L), the forward transition rate F abs is given by the equation The governing equations are and at the ends of domain, respectively.The algorithm allows also reflecting boundary conditions, which is drawn in Fig. 2.

Three-state model
In 1999, a three-state model of myosin II was introduced in the article [2] and further developed in articles [3,16].The last mentioned model is based on an ensemble approach.Authors call these mechanochemical states as the unbound state (motor is loaded with ADP and P), the weakly-bound state (myosin is bent to the actin filament) and the post-power-stroke state, see Figs. 3 and 4. The transition rates between every two states were taken from the article [3] to enable a comparison with the earlier works in the future research.The transitions between the unbound and the weakly-bound states is k 12 40 s −1 and k 21 2 s −1 , respectively.The Fig. 4. Scheme of mechanochemical states of myosin II molecular motor as springs.The unbound state is denoted as 1 , the weakly-bound state as 2 , and the post-power-stroke state as 3 [3] transitions between the weakly-bound and the post-power-stroke states have similar values for both directions, i.e. k 23 k 32 10 3 s −1 .The last transition from the post-power-stroke state to the unbound one is irreversible and can occur with the rate where k m 2.5 × 10 −3 Nm −1 is the spring constant of myosin neck linkers.These linkers are schematically shown in Fig. 4. The parameter d 8 nm corresponds to the length of the myosin power-stroke.The force F 0 = 12.6 pN corresponds to the internal motor load.
The aim of the presented model is to statistically predict the position of the myosin molecular motor using the transition rates values from [3].For the simplification, the rate k 31 is taken as the median value of k 31 (x).The WPE algorithm is applied to the Fokker-Planck equation.
The myosin head movement is controlled by the action potential V (x) and the applied load force F Load (see above) in the Fokker-Planck framework.This is adjusted according to the meaning of each state considered in the model.The unbound state is represented by a constant potential, which equals zero in our model.In fact, however, due to the usage of the WPE algorithm, only the difference in space is important, see (12).The other potentials are given by the Fourier series and respectively (see Fig. 5).The Fourier series are often used for ratchet models in the Fokker-Planck framework [1,18].The amplitude ΔG represents the energy derived from one ATP hydrolysis [1].
Fig. 5.The potential for the unbound state (V 1 ), the weakly-bound state (V 2 ) and the post-power-stroke state (V 3 ) In that article, ΔG equals 25 k B T. The experiments [7,8], however, show that one ATP molecule provides more energy than the treshold necessary for a single myosin state.That is why we choose a half value of ΔG approximately, i.e.ΔG = 12 k B T. The diffusion coefficient (and the next parameter) is chosen according to [1] as D = 5.47 × 10 7 nm 2 s −1 .The length L is the characteristic distance between two actin binding sites, i.e.L = 36 nm [6].
The temperature is considered as a constant parameter and its influence is hidden in the chosen values of parameters D and k ij .We do not study an explicit dependency on the temperature of the surroundings.The external load F Load is set to zero in our study.
The initial conditions are chosen as equilibrated with canonical distribution of energy, i.e.
where φ i is the effective potential φ i = V i − xF Load for a given state i.

Influence of boundary conditions on the total probability
Obviously, the total probability of finding a myosin head over the whole space at each time moment has to be equal to 1. So, this is a good condition to check the boundary condition stability.The time integration is done by MATLAB solver ode15s.
The absorbing boundary condition cannot keep the total probability equal to 1 by its definition.On the other hand, the potential can be modified to catch the myosin head in the area by barriers at the boundary nodes.It can represent, for example, a resistance of the surroundings or a strength in the myosin neck to prevent its rupture.
The reflecting boundary conditions should keep the probability condition by its definition and addition of some barriers at boundaries should not disturb the probability condition.
We have used these barriers having their maximum value 28 k B T, which is high enough to prevent diffusion to exceed it.
The change in the potentials induces also a change in the initial condition via equation ( 20), but the probability axiom is valid for all cases.Thus, the change in the initial distribution of probability density could not cause the change in the control parameter.
According to this result, it is possible to conclude that the absorbing condition with barriers is the best way to obtain a stable total probability value.The system with reflecting boundaries without barriers may be considered as a good approximation.The results in Table 1 show that the reflecting boundary condition with barriers is not suitable for the numerical solution of the problem.It seems that we obtained a strongly overdamped system, which the used numerical integrator cannot deal with.Other numerical integrators in MATLAB such as ode45 provided even worse results -they diverged.Lets us describe the time evolution of individual components of the probability distribution ρ i , i = 1, 2, 3.These components may be interpreted as probabilities of obtaining the myosin head at individual configurations (see Fig. 4 for the unbound, weakly-bound, and post-power-stroke states).
The time evolution of the unbound state ρ 1 presents no significant change.Notice that its potential is constant, i.e.V 1 (x) = 0, except the barriers at the ends of the domain.The only change is its small decrease caused by the transitions to the weakly-bound state, see Fig. 6.
In the weakly-bound state ρ 2 , a non-constant potential is present.There is an evident change during the time evolution, see Fig. 7, where a new local maximum at the beginning of the domain is created.It is caused by a superposition of its own potential V 2 and the boundary barriers.Fig. 6.Probability densities at the beginning and at the end of simulation for the unbound state Fig. 7. Probability densities at the beginning and at the end of simulation for the weakly-bound state Fig. 8. Probability densities at the beginning and at the end of simulation for the post-power-stroke Fig. 9. Probability densities at the end of simulation for the post-power-state in more suitable scale In the post-power-stroke state, big qualitative change happened, see Fig. 8.The maximum caused by the initial condition seems to disappear in the time evolution, but it is only much smaller than the initial condition and thus, not visible, see Fig. 9.There is a similar situation as in the weakly-bound state, just like a mirror twin.The decrease is caused by the transition to the other states.

Conclusions
We tested a particular numerical method for solving the system of partial differential equations (the Fokker-Planck ones) describing the time evolution of probability distribution of a myosin head.The motivation for this problem comes from the modelling of the myosin II as a representative of a molecular motor in a biological muscle cell.We solved the Fokker-Planck equation with four different boundary conditions: absorbing and reflecting boundary in combination with and without potential barriers at the borders of the domain.
As the first conclusion, we found the instability of the conditions which describe reflexion with border barriers.Instead, the omission of the border barriers makes the reflecting border conditions more stable.The expression of this fact is a rather good conservation of the total probability.Nevertheless, the most suitable boundary conditions for numerical calculations seem to be absorbing boundaries with barriers.Therefore, we used them in our modelling.
The obtained numerical results describe well the time evolution of the spatial probability distribution.It is evident how transitions between individual myosin states (unbound, weaklybound, post-power-stroke, see Fig. 3) occur during the time evolution at the characteristic times 10 ms.The values of model parameters used in our numerical simulation were taken from the literature [1,3].The spatial scale is tens of nanometers (36 nm).
The presented results provide a good starting point for more detailed simulations of the problem.Especially, a comparison with experimental results seems to be very important as well as a good interpretation of the model parameters.It will be also interesting to use our calculations in finding some new details concerning the information approach to the statistical physics of mesoscopic systems [12].It is a very fascinating problem, for instance, whether the myosin head might be considered as Maxwell's demon [9].Specifically, we see as useful to contribute to this question, to do numerical simulations of the problem of myosin head localisation (a measurement) and subsequent feedback control of the system.

Fig. 3 .
Fig. 3. Mechanochemical states of myosin II molecular motor used in the three-state model with possible transitions between them [3]

Table 1 .
The obtained total probability at different time steps with different boundary conditions