This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.104.055178v1 89/3/1551    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Xing, J. Articles by Oster, G. Search for Related Content PubMed PubMed Citation Articles by Xing, J. Articles by Oster, G.

From Continuum Fokker-Planck Models to Discrete Kinetic Models

Jianhua Xing ^*, Hongyun Wang  and George Oster ^*

^* Departments of Molecular & Cellular Biology and Environmental Science, Policy, and Management, University of California, Berkeley, California; and Department of Applied Mathematics & Statistics, University of California, Santa Cruz, California

Correspondence: Address reprint requests to George Oster, Tel.: 510-642-5277; E-mail: goster{at}nature.berkeley.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
Two theoretical formalisms are widely used in modeling mechanochemical systems such as protein motors: continuum Fokker-Planck models and discrete kinetic models. Both have advantages and disadvantages. Here we present a "finite volume" procedure to solve Fokker-Planck equations. The procedure relates the continuum equations to a discrete mechanochemical kinetic model while retaining many of the features of the continuum formulation. The resulting numerical algorithm is a generalization of the algorithm developed previously by Fricks, Wang, and Elston through relaxing the local linearization approximation of the potential functions, and a more accurate treatment of chemical transitions. The new algorithm dramatically reduces the number of numerical cells required for a prescribed accuracy. The kinetic models constructed in this fashion retain some features of the continuum potentials, so that the algorithm provides a systematic and consistent treatment of mechanical-chemical responses such as load-velocity relations, which are difficult to capture with a priori kinetic models. Several numerical examples are given to illustrate the performance of the method.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
Dynamical studies of molecular motors fall roughly into three categories. Molecular dynamics (MD) purports to follow the motions of all of the atoms by solving Newton's equations using a variety of semi-empirical potential functions that model the forces between atoms (see, for example, Refs. 1–3). At the other extreme, kinetic models of motor dynamics describe the Markov transitions between a discrete set of states. For example, the status of a catalytic site of an ATP-driven motor is frequently represented by four occupancy states: Empty (E), ATP bound (T), ADP/Pi bound (DP), and ADP bound (D). Transitions between states are given by rate constants that may be force-dependent via an exponential Boltzmann factor. A kinetic model is usually based on the assumption that the configuration space is divided into discrete regions (i.e., potential wells), which are separated by rather high potential barriers (4,5). The consequence of high potential barriers is that the system spends most of its time diffusing within the potential well, and barrier-crossing transitions happen rarely, but instantaneously (6). Quite often the assumption of high potential barriers breaks down for molecular motors. For example, a unique feature of molecular motors is that one mechanical degree of freedom is coupled to the chemical reaction. Under high load, when the motor is performing mechanical work, motion along this mechanical coordinate may be slow compared to other dynamical processes in the system. In a theoretical treatment of the mechanical responses of a molecular motor, such as its load-velocity curve, the mechanical coordinate requires more detailed treatment. Kolomeisky and Fisher (7,8) developed an interesting and important generalization of the kinetic model approach applied to molecular motors by introducing some extra kinetic states along the mechanical degree of freedom. (It will be clear regarding the physical meaning of these states, and the relation between the generalized kinetic models and continuous models later in this work. Please note that, in the following discussions, we call it the "generalized kinetic model" to distinguish it from the "chemical state-only kinetic model".) Another alternative is to treat some degrees of freedom continuously as in the Fokker-Planck models discussed below.

The third approach to model molecular motors is intermediate between all-atom MD simulation and discrete state kinetic models. If one can identify collective coordinates that capture the major conformational motions of the protein, then the mechanical forces driving the system along these coordinates can be captured by a set of potential energy functions defined for each chemical occupancy state. These potential energies can be inferred from the molecular structures, and capture the relevant features of the protein geometry. Then the dynamics is studied by solving the continuous governing equations, which consist of Langevin equations along the geometrical coordinates and kinetic (Markov) jumps between the potentials. Thus the Fokker-Planck formalism replaces the discrete states of kinetic models with continuous potential functions defined on geometrical coordinates that represent the major conformational motions of the protein. We shall refer to these as Markov-Fokker-Planck (MFP) models.

We shall not discuss MD simulation here. The Fokker-Planck equations can be formally obtained from the complete dynamical equations of the system (as in MD simulations) by selecting some primary degrees of freedom, projecting out all the remaining degrees of freedom, and introducing some physically well-justified approximations (9). As studied in the field of chemical dynamics, kinetic models are obtained by approximating the underlying continuum dynamics of the system using a set of discrete states (6,10). As a classical example first studied by Kramers, the dynamics of a double-well system embedded in a heat bath can be described by a single Fokker-Planck equation. By choosing a dividing surface that separates the system into two regions, transitions between these two regions (states) can be approximated by a rate process, and the rate constants can be obtained from the Fokker-Planck solutions. From the viewpoint of numerical computation, kinetic and MFP models are not completely distinct. The master equation for kinetic models consists of systems of ordinary differential equations, whereas MFP models are systems of partial differential equations. When discretized in numerical simulations, both kinetic and MFP models can be reduced to discrete Markov chains, and then the distinction resides in the number of Markov states one assigns to the geometrical coordinates. However, most kinetic models are constructed by selecting the kinetic states a priori, based on biochemical observations, or on an intuitive picture of the protein's motions. MFP models force one to deal more explicitly with the conformational motions, and thus make closer contact with the actual protein geometry. This greater fidelity comes at the cost of having to deal with a continuous geometrical coordinate which, in simulations, is usually discretized into many Markov states (11).

Here we present a new algorithm that generalizes the numerical algorithm developed by Wang et al. (11) for solving the MFP equations. The new algorithm reduces the continuum MFP equations into much simpler discrete jump models, yet retains many of the advantages of the former. This work is motivated by the requirements: first, it is computationally less demanding than the old algorithm of Wang et al. (11); second, it provides sufficient treatment for the mechanical degree(s) of freedom. The latter is particularly important when the slowest dynamics is mechanical rather than chemical transitions—a situation likely confronted in single molecule experiments, where the motor operates under a large viscous load, or a large load force is applied to the motor. However, even in the case where the chemical transitions are rate-limiting, phenomena may arise that are counterintuitive and difficult to capture in a simplistic kinetic model. For example, load velocity curves may not be monotonic, and increasing the load in a range may actually increase the motor velocity (12). An example of nonmonotonic load velocity relation is given in Numerical Examples.

	KINETIC AND MARKOV-FOKKER-PLANCK MODELS

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
Kinetic models represent a system by discrete states, with the dynamics governed by ordinary differential master equations of the form

(1)

Here p is a normalized vector containing state occupation probabilities, K is the transition matrix with its off-diagonal elements k_' giving the chemical transition rate from state ' to state , and the diagonal elements given by Assignment of the kinetic states is usually based on chemical considerations. For example, an ion binding site can be in either an empty or occupied state (13), and a catalytic site of an ATPase can have different nucleotide binding states (4,14). Although Eq. 1 is an evolution equation for the probability vector, the stochastic evolution of an individual system (jumping between the discrete states) is also governed by the transition matrix and can be simulated numerically. Although it is natural to model the occupancy of a catalytic site using a set of discrete states, the mechanical motion is continuous and it is not clear whether a large conformational change can be modeled simply as a chemical transition. To describe the mechanical nature of molecular motors, introduction of some intermediate states becomes necessary. For example, an ATP-binding power-stroke can be modeled by transition from a weakly bound state to a tightly bound state. Fisher and Kolomeisky go further along this line by introducing more kinetic states along the mechanical coordinate to describe more subtle mechanical responses of the system (7,8). They have applied this method successfully to describe the statistics of kinesin and myosin dynamics (15,16).

Next we turn to continuum descriptions. Proteins live in the world of low Reynolds numbers where inertia can be neglected. The timescale of inertia is the time it takes for the motor to forget its current velocity due to friction. For example, the inertial timescale of a 1-µm bead in water is 56 ns (17). The stochastic dynamics of a protein motor generally takes place on timescales much longer than this, and is well described by overdamped Langevin equations of the form (9,18),

,

(2)

where j is the current chemical occupancy state of the system. Here x denotes a mechanical (geometric) coordinate and is a drag coefficient, related to the diffusion coefficient D by the Einstein relation = k_BT/D (k_B is the Boltzmann constant and T the absolute temperature). The value V_j(x) is the potential of mean force as a function of the geometrical coordinate, x, whereas, in chemical occupancy state j, F_Load is the external load force on the motor and f(t) is white-noise (the derivative of a Weiner process). Chemical transitions can accompany motions along the mechanical degree of freedom. The Langevin equation (Eq. 2) is not closed. It governs the stochastic evolution of the mechanical coordinate given the current occupancy state. The dynamics along the chemical coordinates of occupancy states is governed by a discrete Markov model of the same form as Eq. 1, with transition rates that generally depend on the system configuration, x (19). Let K(x) be the matrix of transition rates along the chemical coordinates at position x. The off-diagonal elements k_ji(x) are the transition rates from chemical occupancy state i to state j; the diagonal elements are k_jj(x) = –_{i j}k_ij(x). The Langevin equation (Eq. 2) coupled with a discrete Markov process with transition matrix K(x) describes the stochastic evolution of the motor system.

Brownian fluctuations dominate the dynamics of molecular motors, and so its trajectory is stochastic. However, experiments generally measure only the average quantities, such as mean positions, velocities, and reaction cycle rates. Average quantities can be studied more efficiently by following the evolution of probability densities that are governed by Fokker-Planck equations of the form

.

(3)

Here _j(x, t) is the probability density of the motor being at position x at time t in chemical occupancy state j. The Fokker-Planck equations and the Langevin equations are equivalent (9). This framework has been applied in many theoretical studies of real and generic molecular motors; we refer the readers to Howard (4), Zwanzig (9), and Reimann (20), and references therein.

Generally speaking, kinetic models are simpler than continuum models. In some cases, analytical solutions can be obtained (e.g., Ref. 8), which can provide valuable physical insight. In principle, the two frameworks—i.e., kinetic models (especially the generalized kinetic models) and MFP models—achieve equivalent descriptions if a large number of kinetic states are included to emulate the continuous geometric coordinate. However, the Fokker-Planck (MFP) models have many advantages over a priori kinetic models despite their greater computational complexity.

1. There are clear connections between the spatial potentials in a Fokker-Planck model and the molecular structure so that structural information can be inferred. We believe that this is an important aspect in modeling protein motors, which is not easy to incorporate into an a priori kinetic model. However, the potential-based kinetic models that we will construct below from MFP models can be related to structural information.
   
2. The force-velocity relation is one of the most important characteristics of a molecular motor. Below we will show by an example that the force-velocity relation may be nonmonotonic. Although such a nonmonotonic force velocity relation is naturally accommodated within the framework of Fokker-Planck models, it is not easy to accommodate it in an a priori kinetic model without referring back to the potentials from which the kinetic model was constructed. In some treatments, one assumes a simple form for the underlying potentials so that the effects of the external load can be added to the rate-constant expression analytically (13).
   
3. Motors like myosin may function in groups where there can be cooperative effects so that the system dynamics is not a simple sum of single motor dynamics (21,22). An a priori kinetic model for a single motor is just a phenomenological model for the behavior of the motor when it is not coupled to other motors. When two or more motors are coupled, the multimotor system can be described by another kinetic model. If the coupling is weak, the overall kinetic states can be treated as a combination of individual motor states. However, when the coupling is strong (e.g., motors are tightly coupled by a rigid filament), chemical states of individual motors become less well-defined, as illustrated in Fig. 1. By contrast, in the MFP framework, treatment of the coupling is straightforward.
   
4. In constructing an a priori kinetic model, one usually makes implicit assumptions, which are not easy to discern when revising the model in light of new data. For example, all the possible reaction pathways may not be treated. Consequently, one generally assumes the existence of a dominant reaction pathway (but see, e.g., Ref. 23); however, for multisubunit motors, this notion breaks down. For example, experiments on helicases reveal many broadly-distributed, concentration-dependent pathways (24). On the other hand, in MFP models all reaction pathways are accommodated and none need be excluded without justification, and the existence of concentration-dependent pathways emerges naturally.
   

View larger version (18K): [in this window] [in a new window]   FIGURE 1  (a) Two rigidly coupled motors are driven by a set of two-state potentials, V1 and V2. (b) The coupled motor system has four configurations. Vij refers to the potential with the two motors in states i and j, respectively. It is clear that if the dynamics of the coupled motor system are described by a kinetic model, there is no simple relation between the rate constants of the compound system and those of each individual motor.

	NUMERICAL METHODS

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
Wang et al. (11) developed an efficient and robust numerical algorithm for solving Fokker-Planck equations, hereafter referred to as the WPE algorithm (11,26). This method has the desired property of ensuring detailed balance while computing the correct mean velocities and the variances in a mean-square sense. In the algorithm, a set of continuous Fokker-Planck equations is approximated by the master equation for a jump process on discrete grid points. Each grid point can exchange population/probability with its nearest neighbors along both the spatial and reaction coordinates. The jump rates along the spatial coordinates are calculated based on local steady-state solutions. The jump rates in the reaction coordinates are taken as the rates at the grid points. The algorithm has been proved to be second-order-accurate and robust. Recently, Fricks et al. (27) extended the algorithm to study how motor dynamics is affected by a viscous load elastically linked to the motor. However, a good numerical solution requires that the distance between two neighboring grid points be small enough so that

1. The potential between two neighboring grid points can be well approximated by a linear function (the algorithm constructs the local steady-state solution by assuming a linear potential between two neighboring grid points).
   
2. The chemical transition rate in the cell around a given grid point can be well approximated by the chemical transition rate at that grid point (the algorithm simply uses the rate at the grid point).
   

In studies of molecular motors, system size and lack of precise structural information prevent obtaining reliable potentials directly from first-principle calculations (e.g., molecular dynamics simulations). In a recent article, Xing et al. (28) used the F_o motor of ATP synthase as an example to demonstrate a method of constructing empirical potential energy surfaces from qualitative and quantitative experimental data. Potentials with tunable parameters were first constructed based on experimental observations. Then the parameters were determined by fitting to quantitative experimental data. Finally some of our model predictions were confirmed by new dynamic experiments, and new structural predictions were also made from the potentials. This procedure is analogous to the method of constructing an empirical potential energy surface that is widely used in chemical dynamics studies: potentials determine the forces that govern dynamics, and potentials can be related to structures.

Construction of the potential surfaces itself may be a complicated procedure, which we do not address here. However, the issues we do address are: 1), how to construct a faithful and computationally efficient algorithm, for a given set of potentials; and 2), how to formulate a natural method for constructing simple kinetic models based on the underlying continuum Fokker-Planck models. We will see that both of these two goals are achieved by a generalization of the WPE algorithm. The generalized algorithm retains the use of local steady-state solutions, but relaxes the two assumptions listed above. The remainder of the article is organized in the following order. First we derive the generalized algorithm in one dimension and higher dimensions. Next some numerical tests are presented. Finally, limitations of the method are discussed. The consistency, stability, and convergence of the generalized algorithm are analyzed in Appendix B.

	THE GENERALIZED ALGORITHM

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
We begin by observing that many systems involve phenomena that operate on widely disparate timescales. Thus one can use adiabatic approximations (i.e., singular perturbation) to treat separately degrees of freedom with fast timescale motions that adjust quickly to the slow timescale motions. A well-known example of this is the Born-Oppenheimer approximation in quantum mechanics. Another closer related example is the quasi-steady-state approximation in chemical dynamics studies (29). If one is only interested in long timescale dynamics, then it is usually a good approximation to assume that a steady state is established for the high-frequency degrees of freedom without appreciable evolution of the low-frequency degrees of freedom. Below we describe the generalized algorithm. Derivation details are given in the Appendices.

Algorithm for one-dimensional equation with no reaction
Consider first a one-dimensional Fokker-Planck equation with no chemical transitions:

(4)

This equation describes the stochastic evolution of a particle driven by potential V(x). We divide the computational region (0, L) into N sub-intervals. Let x = L/N and x_i = ix. We call the sub-interval (x_i–1, x_i) the i^th cell in the spatial direction. We use a jump process to approximate the continuous motion of the particle as shown in Fig. 2 (notice that cells, and not grid points, are used in the current algorithm). The particle can jump from a cell to the neighboring cells. For a one-dimensional problem, a cell has two neighboring cells in the spatial dimension. Let r_ii+1 denote the jump rate from (x_i–1, x_i) to (x_i, x_i+1), and r_i+1i the jump rate from (x_i, x_i+1) to (x_i–1, x_i). In each cell, we define the probability of the particle being in that cell as

Then Eq. 4 is approximated by the master equation of the jump process:

(5)

View larger version (9K): [in this window] [in a new window]   FIGURE 2  Schematic illustration of the algorithm. The jump process defined in Eq. 8 approximates the continuum model from Eq. 4. Note that the Markov states are now intervals corresponding to the integration domains.

so that

Consistent with normal statistical mechanics definition, the free energy is related to the partition function defined in the cell. Our definition of the free energy term is also consistent with what used in chemical dynamics (6). The right-hand side of Eq. 4 involves second-order partial differentiation; thus, it takes two conditions to uniquely specify a steady-state solution. The rates r_ii+1 and r_i+1i are determined by equating the numerical flux with the flux of the local steady-state solution specified by

This local steady-state solution can be written as a linear combination of the Boltzmann distribution and the quantity

where q_i(x) is the steady-state solution of Eq. 4 specified by q_i(x_i–1) = 0 and To calculate the jump rates r_ii+1 and r_i+1i, we assume that the solution of Eq. 4 in the interval (x_i–1, x_i+1) is approximately in a steady state, which has the general form of

(6)

The probability flux based on the local steady-state solution is given by

(7)

The two constants c₁ and c₂ are determined from the conditions

Solving for c₂ and substituting into Eq. 7, we obtain

Comparing this with the numerical flux leads immediately to the expression for r_ii+1 and r_i+1i,

(8)

where

and

Functions and are the jump rates (normalized by the diffusion coefficient) derived using a local steady-state solution. In the jump process with this set of rates, the local steady state (without reaction term) is preserved exactly. As stated in Wang et al. (11), there are two motivations for using a local steady state without reaction. First, the exact solution can be written out analytically. Second, the local steady-state solution is a good approximation of the real time evolving solution. Consider the local system consisting of two adjacent cells. This local system has a spatial size of 2x. Assume this local system is isolated from other cells. The timescale for the diffusion to drive it to a steady state is of the order (x)². The time for the convection to change this local system appreciably is of the order (x), whereas the timescale for the chemical reaction to change this local system is of the order O(1) (the coefficient may be large). Thus, although this local system is not isolated from other cells, it is in a pseudo-steady state (relaxation to the steady state is much faster than the evolution of the steady state). Therefore, this approach is justified at least in the limit of x converging to zero. Notice that the local steady-state solution is actually time-evolving because the constraints on the local steady-state solution vary with time. Here we have explicitly included the potential, V, in the notations of and As we will see, this notation will be very convenient in discussing two-dimensional problems.

The rates r_ii+1 and r_i+1i automatically satisfy detailed balance

This is an important constraint on the solutions, as discussed in Wang et al. (11) and Elston and Doering (30).

In the absence of chemical transitions, the steady-state solution obtained with this algorithm is exact. Notice that the algorithm does not assume a linear potential in (x_i–1, x_i+1). This is one of the reasons that it can represent the continuous Fokker-Planck equation with a small number of cells. In contrast, the WPE algorithm approximates the potential within (x_i–1, x_i+1) by a linear interpolation. Because we assumed that a steady state is established in each interval (x_i–1, x_i+1) for a time-dependent solution, short-time resolution is lost if only a small number of cells are used in the simulations; this will be illustrated in the numerical examples below.

The jump process illustrated in Fig. 2 is a discrete kinetic model. Thus, the new algorithm provides a natural way of building a discrete kinetic model from the underlying continuum MFP model.

Algorithm for one-dimensional equation with reactions
Consider a variant of Eq. 3 where the external force is absorbed into the potentials:

.

(9)

The discrete reaction coordinate can be visualized as orthogonal to the continuous geometric coordinate. First, we introduce several notations for the numerical discretization:

1. is the probability that the system is in cell (x_i–1, x_i) and in state j.
   
2. is the transition rate in the spatial dimension from cell (x_i–1, x_i) to cell (x_i, x_i+1) in state j
   
3. is the transition rate in the reaction dimension from state j to state j' in cell (x_i–1, x_i).
   
4. is the free energy of cell i in chemical state j.
   

As in the previous section, the continuous motion along the geometric coordinate is approximated by a jump process. Then Eq. 9 is approximated by the master equation of the jump process,

(10)

where the jump rates along the geometrical coordinate, and are calculated using Eq. 9, which was derived for equations with no reaction. Using Eq. 8 here implicitly assumes that the presence of a chemical reaction does not significantly affect the local steady-state solution. This assumption is valid when the cell size is small enough such that the diffusion within the cell is faster than the reaction, and can be considered homogeneous. The jump rates along the chemical coordinate are calculated by averaging the chemical reaction rates in cells with Boltzmann weights:

(11)

The rates and automatically satisfy detailed balance:

The derivation of this result and Eq. 11 is discussed in Appendix A.

Algorithm for two-dimensional equations
Consider a Fokker-Planck model with two geometric coordinates, but with no chemical transitions:

(12)

where D^(x) and D^(y) are the diffusion coefficients in the x- and y-directions, respectively.

We divide the computational region (0, L^(x)) x (0, L^(y)) into N^(x) x N^(y) subregions. Let x_i = ix, and y_j = jy. We call the subregion (x_i–1, x_i) x (y_j–1, y_j) the cell (i, j). We introduce the following notation for the two-dimensional discretization:

1. p_(i,j) is the probability that the system is in cell (i, j).
   
2. r_(i,j)(i+1,j) is the transition rate in the x-dimension from cell (i, j) to cell (i+1, j).
   
3. r_(i,j)(i,j+1) is the transition rate in the y-dimension from cell (i, j) to cell (i, j+1).
   
4. G_(i,j) = –k_BT is the free energy of cell (i, j).
   

Equation 12 is discretized as

where the numerical fluxes in the two spatial dimensions are

and

(13)

where

Details of the derivation are given in Appendix C. The transition rates in the y-dimension are obtained in a similar way.

For two-dimensional equations with chemical reactions, the jump rates in the reaction dimension are calculated by averaging the chemical reaction rates in cells with Boltzmann weights as before. Equation 13 can be easily generalized for problems with more than two dimensions.

Next we illustrate the algorithm with some simple models.

	NUMERICAL EXAMPLES

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
In the following calculations, the linear equations governing the steady-state solutions were solved using a sparse matrix solver in MatLab (The MathWorks, Cambridge, MA). Integrals used in the jump rates of the XWO algorithm were calculated using the fourth-order Simpson method. MatLab codes are available on request.

A two-state one-dimensional model
Here we consider a simple two-state system designed to capture the main physics of a typical molecular motor where the system switches between a pair of potentials. Each potential is a periodic function with period L. The two potentials are given by (see Fig. 3)

where

and G₀ = 20 k_BT, G = 30 k_BT, and k_BT = 4.1 pN/nm. After one cycle the motor has returned to its initial chemical state, but the motor position has advanced by L, whereas the free energy of the system (motor plus environment) decreases by G (e.g., ions are transported from high concentration region to low concentration region, or ATP molecules are hydrolyzed).

View larger version (16K): [in this window] [in a new window]   FIGURE 3  The potentials used in the one-dimensional numerical example.

where k₀ = 2000 s^–1, and within the window (x_a2, x_b2), with the transition rates

Therefore, the governing equations are

(14)

with If the geometric coordinate, x, is periodic with period L = 2/3 (i.e., a rotary motor), then the diffusion constant is chosen as D = 10⁴ radian²/s; the mathematical formulation is the same for a translational motor.

Case 1
The chemical transition regions are localized around the potential minima, x_a1 = 0.55, x_b1 = 0.65, x_a2 = 1.55, and x_b2 = 1.65. Results are shown in Fig. 4 a. Since the chemical transitions are rather localized, with the WPE algorithm many numerical cells are needed to cover the transition region sufficiently. On the other hand, performance of the XWO algorithm is remarkable.

View larger version (12K): [in this window] [in a new window]   FIGURE 4  Numerical convergence tests for case 1 of the one-dimensional problem. (Left). The force-velocity curves calculated with the WPE algorithm. The numerical cell sizes are L/N, where L = 2/3 is the periodicity. (Middle). The force-velocity curves calculated with the XWO algorithm. (Right). The averaged relative errors for one force-velocity curve estimated by where M is the number of data points.

View larger version (13K): [in this window] [in a new window]   FIGURE 5  Similar to those in Fig. 4, but for the case 2 of the one-dimensional problem. The results with N = 5 in b are obtained with five uneven numerical cells per period (see text for details).

Two-dimensional example
Here we connect two motors with an elastic linkage (see Fig. 6). One motor (with its position denoted by x) drives the other motor (with its position denoted by y). Parameters describing motor x are the same as in case 1, D_x = 10⁴ radian²/s, L_x = 2/3. For the load motor y, the potential forms are similar to those used in the one-dimensional model, except they are reflected at x = 0,

and D_y = 10⁴ radian²/s, L_y = 2/10, G₀ = 12 k_BT, and G = 6 k_BT. The chemical transition regions are located at y_a1 = 0.35, y_b1 = 0.45, y_a2 = 1.35, and y_b2 = 1.45. The elastic linkage is given by

where k 9.6 pN/nm per rad². The motion of motor x exerts a torque on motor y, forcing the latter to move up its free energy gradient. The potential barriers of motor y are not high, so there is a certain probability that motor y slides over a potential peak without making a chemical transition. In other words, motor rotation and chemical transitions are not perfectly tightly coupled.

View larger version (13K): [in this window] [in a new window]   FIGURE 6  Convergence test for the two-dimensional system with the XWO algorithm. (a) Schematic illustration of the two motors coupled by an elastic linkage. (b) The rotation rates (upper) and the relative errors (bottom) as a function of the numerical cell numbers Nx = Ny = N. The numerical cell sizes are (Lx/Nx, Ly/Ny).

	CONCLUDING REMARKS

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
We have presented an efficient numerical algorithm for solving the kinds of Fokker-Planck equations used in molecular motor studies. The new algorithm has been tested on several systems, and is found to be numerically accurate even for fairly large spatial step sizes. The algorithm also provides a natural procedure for constructing simpler kinetic models starting from continuum Fokker-Planck models. The algorithm links continuum Fokker-Planck equations naturally with mechanochemical Markov chains while retaining many features of the original modeling framework. Detailed balance is automatically preserved and, because potentials are not assumed to be linear between two adjacent grid points, satisfactory numerical accuracy can be achieved even for fairly large spatial step sizes. The effective chemical transition rates are calculated by averaging over each numerical cell. With the proper distributions, numerical accuracy is significantly improved over the WPE algorithm. Kinetic models constructed using the algorithm reproduce the force-dependence of transition rates consistent with that in the original continuum Fokker-Planck model, such as the puzzling observations of nonmonotonic load-velocity behavior (12).

We believe that this algorithm provides a new tool for theoretical modeling of molecular motors and other mechanochemical systems. For complex systems such as multimeric motors like the portal protein, helicases, and PilT, the total number of chemical states is quite large. Consequently, the computational cost of using a small spatial step is prohibitively high. For such motors, the new algorithm may be the only viable way to explore a large parameter space. In practice it is easier to identify (and exclude) those sparsely populated high energy chemical states with a potential-based model. Since the computational effort to solve linear equations grows as a power law with the number of cells, the advantage of the XWO algorithm over the WPE algorithm is appealing for complex systems. For example, the new algorithm was also tested on a combined F₁-F_o system (to be addressed in a future article). The F₁ model has 64 chemical states, with eight chemical occupancy states to describe ion binding sites of the F_o motor. The combined F₁-F_o model is a two-dimensional system with 64 x 8 = 512 chemical states. The necessary number of cells for the WPE and XWO algorithms are 10⁹ and 10⁶, respectively. In simulations of the combined F₁-F_o system, the new algorithm yields results similar to those of the WPE algorithm, but with a computational cost several orders-of-magnitude lower.

In treating multidimensional problems, one implicitly assumes that the cross-potential terms can be approximated by linear relationships. This is a serious limitation of the current algorithm. The new algorithm shows greater advantages over the old WPE algorithm for systems where the coupling potential terms are slowly varying compared to the direct (diagonal) terms. The current algorithm can be further improved by combining it with more elaborate finite element treatments (32).

	APPENDIX A: DERIVATION OF EQ. 11

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
To calculate the reaction transition rates and used in discretization Eq. 10, we compare the numerical and the theoretical probability flux along the reaction coordinate between state j and state j' in cell (x_i–1, x_i),

where Since the exact solution _j(x,t) is unknown, we approximate it by a local Boltzmann distribution,

(15)

Substituting this into the theoretical flux and comparing with the numerical flux, we obtain the rates given in Eq. 11. The Boltzmann distribution approximation is justified for two extremes of motion-reaction couplings:

1. For motors in which the chemical reactions are well coordinated with the mechanical motion, the mechanical motion cannot continue until the reaction switches the system to another chemical state. For example, in the F₁ ATPase, each chemical transition occurs at a specific rotational location, and rotation cannot continue until this transition is completed. In this case, the mechanical degree of freedom is thermally equilibrated in the local reaction region, and the probability distribution in the local reaction region is well approximated by the Boltzmann distribution.
   
2. For motors in which the chemical reactions are not affected by the motion, k_jj'(x) is independent of x. In this case, we have independent of i.
   

Of course, as the cell size goes to zero, the rates obtained with the Boltzmann distribution weighting converges to the rates obtained with the exact solution weighting. In this sense the Boltzmann distribution approximation is always justified for small cell size. Numerical examples show that the Boltzmann distribution approximation also works well for moderate cell size.

One advantage of using Boltzmann distributions approximation is that detailed balance is exactly preserved. The exact transition rate functions k_jj'(x) and k_j'j(x) satisfy detailed balance

which implies

Detailed balance for and follows immediately.

	APPENDIX B: CONSISTENCY AND STABILITY

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
In this Appendix, we show that the new method is second-order accurate for smooth potentials. Specifically, we will first show that the method is consistent with the Fokker-Planck equation. Then we will show that the method is stable with respect to the L² norm. Once we have consistency and stability, the Lax equivalence theorem implies convergence (33). The approach of the analysis used here is similar to that used in Wang et al. (11).

For simplicity, we present the proof of the consistency and stability for the Fokker-Planck equation:

(16)

The approach can be extended to Fokker-Planck systems with an arbitrary number of states.

The numerical method for this Fokker-Planck equation is

(17)

where in the time dimension it is discretized using the Crank-Nicholson method.

Consistency
We want to show that the truncation error of the method is second-order in both the time and spatial dimensions. The local truncation error is the residual of the method applied to an exact solution. Let (x,t) be an exact solution of Eq. 16. When we substituting into Eq. 17, the residual term is the local truncation error. To find the order of the local truncation error of Eq. 17, we expand every term in r_ii+1 and r_i+1i, around x_i.

Consider two functions:

and

Expanding yields

and

where and are some smooth functions evaluated at x_i.

Substituting and into r_ii+1 and r_i+1iwe obtain

(18)

where is some smooth function evaluated at x_i.

Expanding around x_i, we have

where and c_i^(k) is some smooth function evaluated at x_i.

Substituting and into the numerical flux yields

(19)

where is some smooth function evaluated at x_i. Using we arrive at

(20)

From this result, it is straightforward to show that the local truncation error of Eq. 17 is second-order in both the time and spatial dimensions.

Stability
Next we prove that the method Eq. 17 is stable with respect to the two-norm. We first introduce some notation:

Multiplying both sides of Eq. 17 by and summing over i, we obtain

Using summation by parts and using periodic boundary conditions yields

Using the identity

we get

(21)

Since r_ii+1 and r_i+1i are both positive, the first term on the right side of Eq. 21 is nonpositive.

Using Eq. 18, we have

Since V''(x) is bounded, when x is small enough there exists a constant C such that

(22)

Using summation by parts and Eq. 22, we can write the second term on the right side of Eq. 21 as

Thus, Eq. 21 becomes

which immediately leads to

Therefore, Eq. 17 is stable with respect to the two-norm.

	APPENDIX C: ALGORITHM FOR TWO-DIMENSIONAL EQUATIONS

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS THE GENERALIZED ALGORITHM NUMERICAL EXAMPLES CONCLUDING REMARKS APPENDIX A: DERIVATION OF... APPENDIX B: CONSISTENCY AND... APPENDIX C: ALGORITHM FOR... ACKNOWLEDGEMENTS REFERENCES

 
We will use approximate local steady-state solutions to determine the numerical transition rates. In the one-dimensional case discussed above, the general steady-state solution has only two undetermined coefficients. However, unlike the one-dimensional case, the exact steady-state solution of the two-dimensional equation cannot be determined from the probabilities of two neighboring cells. Therefore, we use a local mean-field approximation to separate the dependence on the two spatial dimensions. Note that a global mean-field treatment is generally not a good approximation, since it omits correlations between different degrees of freedom.

Separable potentials
First, we consider a simple case where the potential can be decomposed into two single variable functions:

This allows us to separate the dependence on the two spatial dimensions exactly. We seek a steady-state solution in (x_i–1, x_i+1) x (y_j–1, y_j) satisfying

Using a method similar to the one we used in the one-dimensional case, we obtain that the probability flux through the right boundary of (x_i–1, x_i) x (y_j–1, y_j) is

Comparing the numerical flux with we immediately obtain

(23)

From the definition of and we can express the transition rates in terms of V(x, y) instead of ₁(x) and ₂(y):

(24)

where

This is Eq. 13 in the text.

The transition rates in the y-dimension can be obtained in a similar way.

Nonseparable potentials
Next we consider the general case V(x,y) = ₁(x) + ₂(y) + ₁₂(x,y). For the general case, the decomposition can still be done locally and approximately by Taylor expanding only the cross-term ₁₂(x, y).

When the cell size is not small, the approximate decomposition may still be valid as long as the coupling term is slowly varying. If we use Eq. 23 to calculate the transition rates, we have to write out the approximate decomposition explicitly. It is not obvious how to find an approximate decomposition that preserves detailed balance exactly. Fortunately, Eq. 13 (which can be obtained by an inverse procedure of the decomposition) allows us to calculate the transition rates directly from V(x, y). All we need to know is that the decomposition can be done locally and approximately. Another advantage of Eq. 13 is that it preserves detailed balance. The free energy of cell (i, j) is

Thus, the transition rates given by Eq. 13 satisfy

	ACKNOWLEDGEMENTS

TOP ABSTRACT INTRODUCTION KINETIC AND MARKOV-FOKKER-PLANCK... NUMERICAL METHODS