| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
* Department of Physiology and Biophysics, State University of New York, Buffalo, New York; and Physics Department and Biophysics Center, University of Illinois at Urbana-Champaign, Urbana, Illinois
Correspondence: Address reprint requests to Frederick Sachs, E-mail: sachs{at}buffalo.edu.
 |
ABSTRACT |
|---|
|
INTRODUCTION |
|---|
). Motor proteins, such as kinesin (2,3), myosin (4), and dynein (5,6), are now studied at the single-molecule level. The motor protein converts chemical energy into mechanical energy by hydrolyzing ATP. The mechanical energy is used for transporting cargo in discrete steps along cytoskeletal filaments, such as microtubules for kinesin and dynein, and actin for myosin. Since motor function is independent of the length of the substrate, the process can be simplified to an infinite chain of identical reaction units (7), where a unit is the set of conformations assumed by the motor protein while it moves one step along the substrate. The location of individual motors as a function of time can be measured with nanometer precision using fluorescence microscopy (6,8,9). Typically, a fluorescent probe is attached to the motor, and a charge-coupled device camera on an optical microscope records the position of the probe. The motors—driven by ATP hydrolysis—track along cytoskeleton filaments immobilized on a glass substrate (8), and their location is traced from frame to frame. The frame rate determines the time resolution of the measurements. The data consist of a time series of the probe position projected along the filament axis. Since the motor proteins generally move forward in discrete steps, the data look like a staircase, although backward steps may occasionally occur, as predicted for all reversible reactions.
A staircase step, i.e., a segment in the data where the measured position of the probe is constant, is called a "dwell". The duration of each step (the dwell time) is stochastic, with exponential distribution. There may be multiple conformational states associated with a single position of the probe—analogously to the multiple "closed" states of ion channels—and transitions between them are not distinguishable as individual events. Nevertheless, the existence of these unobservable states can be inferred statistically from the distribution of the observed events (10). We denote by "state" the combination of biochemical and positional states. Due to the finite sampling time and exponential distribution of event durations, there will always be missed events. The effect of missed events upon the interpretation of reaction kinetics has been extensively studied in the context of single ion channel kinetics (11–14). In the case of molecular motor data, the same considerations apply. This means that not only are the estimated rate constants in error, but the motor may take one or more steps during a single frame and some of the observed movements will reflect multiple steps.
In this article, we present a maximum likelihood algorithm that extracts the most likely (highest probability) dwell time sequence from the data and returns the global parameters that describe the step size distribution and the rate constants (10). The method relies upon the large knowledge base developed over the last 30 years for the analysis of single ion channel data (11–30). Our main objective is to realistically model the mechanochemistry of the motor in terms of the kinetics and step size distributions, and to take into account instrumental limitations, including various noise sources. The noise components can be described by a continuum of states. These include wideband noise arising from random photon statistics, digitizing errors, and the Brownian motion of the motor. There are also low frequency noise sources arising from drift of the microscope, and errors in correcting for the bending of the substrate filament. We have modified and combined algorithms for hidden Markov models (31) to describe the discrete states of the motor, and Kalman filters (32) for the continuous states of the noise. Thus, we model the mechanochemistry with a periodic Markov chain, as required by the identical chemistry of each step, whereas the Kalman filter models the variability of the motor step and of the baseline. A modified Baum-Welch algorithm (31) recursively extracts (i.e., "idealizes") the most likely sequence of motor positions (the dwells), and subsequently calculates maximum likelihood estimates for the state transition probabilities, for the mean and standard deviation of the step size distribution, and for the measurement noise. We found that including the Kalman filter also made the algorithm less sensitive to initial parameter estimates.
The algorithm is general, capable of dealing with arbitrary kinetic models, with a wide set of available a priori constraints and step size distributions. The most common kinetic models in the literature have one or two states per reaction unit, irreversible kinetics and uniform steps. Our method is applicable to any staircase-style data, such as that generated by F1-ATPase, where the rotational angle of the rotary motor is traced from frame to frame (33,34), tracking RNA polymerase along the DNA template (35), or tracking topoisomerase activity (36). A distinct advantage of our method is that it fully utilizes the correlation between sequential dwell times that is a consequence of the topology of connections between states (37,38). Since the algorithm is fast, it allows the user to easily compare different models—quantitatively—by computing the likelihoods. Our method is global, in that the model is fitted to all the data at once, whether in a single or in a collection of data files. In comparison, the t-test (39) is a local extraction method, with the generally implicit assumption of a model with only two states and with uniform steps.
We take advantage of the periodic structure of the Markov model (10) to make the algorithm efficient and fast. Staircase data can be fitted to a traditional finite state Markov model (40), but—with the assumption that the minimum number of states is the number of discrete positions of the motor—large data sets would generate huge transition matrices. For example, if there were 100 steps in the data set, one would have to use a matrix of at least 104 elements. Our algorithm is implemented in Windows, and freely downloadable along with other tools, such as single molecule simulators (41). The user graphically inputs a state model of three consecutive reaction units, and a set of initial parameter values. The model can include constraints such as detailed balance or fixed rates, etc. The program internally creates a list of the most likely state at each data point, and uses it to create a sequential list of the dwell times at each position of the motor. This idealized list is an abbreviated data set (there are fewer dwells than there are data points) that can be used to further explore different kinetic models (10).
Our results show that, as with all fitting algorithms, confidence in the output depends upon the signal/noise ratio (SNR) and the length of the data set. We define the SNR as the ratio of the mean step size (allowing for multimodal distributions) to the standard deviation of the background noise. Ideally, the SNR should be >2. The noise level of current experiments is 
1–3 nm (5,6,8,9), so that steps as small as 8 nm (e.g., kinesin or dynein) can be successfully analyzed. As shown further, the SNR has a strong effect upon the estimated kinetic parameters, and in many cases it is worth sacrificing the time resolution by reducing the bandwidth of the data with appropriate resampling to maintain the validity of the first order Markov assumption. Our studies show that, as expected, the variance of the step size (whether intrinsic to the mechanochemistry or caused by instrumental noise) has a significant effect on the step resolution, and in general we recommend selecting a unit model with only a few steps of fixed amplitude (a multimodal distribution of the step sizes), without intrinsic variance. This multimodal distribution may reflect, for example, the availability of multiple binding sites for the same motor step, or an asymmetrical positioning of the fluorescent probe on the motor.
|
MODEL AND ALGORITHM |
|---|
,43). Myosin V walks with a hand-over-hand mechanism (8) along the actin filament, with the stalk taking 37 nm steps per ATP hydrolyzed. To walk, the motor alternately rotates its two heads about the stalk: the rear head moves 74 nm forward (twice the stalk movement) to become the leading head, and so on (Fig. 1 A). Kinesin has a similar mechanism (9).
|
|
STAIRCASE DATA |
|---|
Pi+2 in Fig. 1 B is "second order". Inferring the jump order from the observed jump amplitude can only be done in a probabilistic way, as discussed below. We denote by µJ and 
J the mean and standard deviation of the jump size corresponding to a displacement between discrete positions of the motor.
A "dwell" is a horizontal segment in the staircase sequence. The observed position of the motor is constant during a dwell. We say observed because the motor may in fact undergo reversible transitions during this time that are too short to be observed. The dwell times have exponential distribution. The "amplitude" of a dwell is the mean vertical position (the y coordinate) of those data points within the dwell time. We denote by µk the mean amplitude of the kth dwell in the staircase sequence. By definition, two consecutive dwells are separated by a jump. However, two consecutive dwells do not necessarily correspond to consecutive positions because of the potential for missed events. We approximate the measurement noise of the probe position (primarily determined by the number of photons per pixel, per frame) as a 
-correlated Gaussian with zero mean and standard deviation M.
|
STEP SIZE DISTRIBUTION |
|---|
The observed jump amplitude—noted "J" in Fig. 1 B—follows the distribution of the step size, but it is also affected by the location of the fluorescent probe within the motor molecule (9). Different jump patterns occur as follows:
, where N refers to the Gaussian distribution function. 
. For the physical models under consideration, the two unequal steps must alternate, hence the weighting factors P1 and P2 are in fact conditional probabilities that alternate between 0 and 1. The sum of two such unequal jumps is equal to twice the jump of the stalk:
. Regardless of fluorophore position, the molecular structure remains the same, and thus it seems reasonable to assume that the combined variance of a pair of unequal jumps is twice the variance of the uniform jump. Hence, we also assume that the relative variance of each jump in the pair is proportional to the relative jump amplitude: 
. 
. Note that, since the amplitude of first order jumps is modeled as a random Gaussian variable, the amplitude of a higher order jump is also a random Gaussian variable. As with any fitting program, the most rapid convergence and most reliable parameters come from the best starting guesses. If one knows a priori the location of the probe within the motor, the appropriately constrained model should be used in the analysis. |
KINETIC MODEL |
|---|
), where each unit corresponds to a position on the substrate. A minimum kinetic model must include at least two distinct states per unit (Fig. 1 C): an ATP binding step and a position translocation step. Although the motor appears stationary, it actually undergoes conformational transitions between these two (or more) states. If NS is the number of states per reaction unit and NP is the number of all the substrate positions occupied by the motor within a given data set, then the process can be described by a Markov model with NS x NP states (10). The frequency of transition between kinetic states is quantified by rate constants, conventionally grouped into a rate matrix Q. For our idealization algorithm, we focus on the probability of observing a transition between any two states, within a sampling interval dt. These probabilities are grouped into a matrix noted A, which can be calculated as 
. The properties of the Q and A matrices—notably the periodicity constraints—and their computational details are given in Milescu et al. (10). Transitions between states within the same reaction unit cannot be directly observed, but their statistical properties can be inferred from the distribution and cross correlations of dwell times. |
LIKELIHOOD FUNCTION |
|---|
 | (1) |
. "Topological" refers to the number of states and their connectivity, and the parameters describe the state transitions and the step size distribution. Note that the absolute value of the likelihood is only used to compare models. The likelihood function includes all possible states occupied by the motor at each time t. Its calculation must take into account all possible state sequences and their intrinsic probabilities
 | (2) |
Of special interest is the sequence of states that maximizes the likelihood
 | (3) |
Due to the memory-less character of a first order Markov process and the assumption of -correlated measurement noise, the argument in the above expression can be simplified to
 | (4) |
Note that for simplicity we have excluded M(), but the dependency on model and parameter values is implicit. The probabilities in the above expression have been presented many times before, for ion channels (30,44,45). Briefly, the state conditional probabilities are the elements of the transition matrix A:
 | (5) |
The "hidden" character of the Markov process comes from the fact that one has to observe the state sequence in the presence of noise, so that neither the state nor the amplitude can be known unequivocally at any data point. The probability of picking the state of correct amplitude at any given data point generally comes assuming a Gaussian distribution of errors. Thus, the conditional probability of making a measurement yt given a state x is the following Gaussian:
 | (6) |
We cannot directly measure the size of the motor steps, due to finite temporal resolution and background noise. Hence, the step sizes must be inferred from the mean dwell amplitudes µk, discussed above. The result is that the model contains some parameters that are stochastic quantities—the µk values—and the likelihood function must be modified to include the probability of a particular set of mean dwell amplitudes, as follows:
 | (7) |
, is the goal of the idealization algorithm presented next. |
IDEALIZATION ALGORITHM |
|---|
), or the Viterbi algorithm (46,47) is used to optimally idealize Markov data in the presence of noise when the parameters are known. The result is the most likely state sequence XML, out of all possible sequences. Given that the true parameters are unknown (except for simulations), the strategy is to run an optimizing algorithm based on the Expectation-Maximization framework (48), such as Baum-Welch (31), or segmental k-means (30,49). There are several reasons why we cannot directly apply these algorithms to molecular motor staircase data:
Fortunately, the periodic structure of the Markov model (10) permits a simple solution to the unknown number of states. We advance hierarchically in complexity, starting with a small number of unitary steps, and add more states as it becomes necessary. To gain efficiency, the state space is truncated to eliminate those transitions that are not likely to happen (10). To deal with the complications caused by step size variability and baseline drift, we take an approach based in part on our previous work idealizing single-channel data with nonstationary baseline (50). Thus, we use the hidden Markov model to describe the state transitions of the motor, and a continuous Gaussian model (51) to describe the baseline and—more importantly—the variability of the motor step size. Note that only the relatively small and continuous component of step variability is handled this way, whereas the larger and discrete step variability, such as that due to alternating small and large steps, is explicitly included in the Markov model. Thus, we regard all motor steps of the same kind (e.g., small or large, etc.) as having constant amplitude throughout the data, and account for variability in the observed jumps through the baseline noise distribution. In its simplest form (i.e., an unbiased random walk), the baseline model is given by the equation
 | (8) |
is Gaussian noise modeled as 
that adds random changes to bt. The standard deviation 
is assumed small relative to the step size, as the baseline is expected to deviate little between two consecutive data points. To allow for step variability, 
is augmented at jump points by a quantity proportional to the step variance and the jump order. Note that, in addition to baseline drift and step variability, the Gaussian process implicitly includes the uncertainty in the initial step size estimate.
The Kalman filter (32,51) is commonly used to estimate the continuous state sequence generated by a noisy Gaussian process—the baseline position in our case. For each data point, the filter provides the most likely baseline from a Gaussian probability distribution, with mean bt and variance Vt. Obviously, we cannot apply the Kalman filter directly to the staircase data, as the Markov data are summed with the baseline, and there is some cross talk between the Markov and Kalman estimators, since a variation in probe position could be attributed to either process. It is the job of the optimizing algorithm to best separate the two processes. Note that, even though the baseline drift is essentially deterministic, its direction is initially unknown. Hence, from an algorithmic point of view, the baseline can be conveniently modeled as a random process with small variance. Our tests showed that adding a deterministic bias is an unnecessary complication.
We implemented the idealization as an Expectation-Maximization (EM) optimizer, based on a modified Baum-Welch algorithm. The EM optimizer alternates two computational steps:
These two steps are iterated until satisfying a convergence criterion (e.g., only a small change in likelihood). The final Expectation step finds the most likely sequence of Markov states and calculates the likelihood, whereas the Maximization step obtains the maximum likelihood parameters. Next, we explain the two steps of the algorithm.
Expectation—state inference
The Expectation is divided into Forward and Backward steps (31). The Forward step recursively calculates the probability of observing the data points {y0,...yt} and ending up in state i. This probability is denoted 
. The Backward procedure does the complementary calculation: 
is the probability of observing the data {yt+1,...yT}, having begun in state i. Additional quantities are also calculated: 
is the probability that the state is i at time t, given the entire data sequence {y0,...yT}; 
is the probability that the state is i at time t and j at time t + 1, given the entire data sequence. All these probabilities are calculated from the current set of parameters, as presented in detail in the Appendix. The 
, ß, 
, and 
calculated here are the dependent variables used in the Maximization step. Idealization is the specification of the most likely state for each data point, i.e., the state index i that maximizes t:
 | (9) |
Knowing the state, we also know the position, and from that we can calculate the expected amplitude µk of each dwell. Note that the sequence of most likely states 
so obtained is not necessarily equal to the most likely sequence of states XML (Eq. 3). Strictly speaking, to obtain XML, we should run Viterbi (46,47) as the final Expectation step, but our tests showed that for all practical purposes, the two state sequences are identical. Our tests also showed that replacing the Forward-Backward procedure with Viterbi throughout the computation, i.e., using a segmental k-means algorithm (30,49), results in significantly poorer performance.
Maximization—parameter reestimation
We want to estimate several parameters: the state transition probabilities stored in the A matrix (10) and the parameters controlling the amplitude distribution, i.e., the measurement noise M, the mean µJ and the standard deviation J of the step size distribution, and the mean amplitude of each dwell µk. The transition probabilities aij are reestimated (updated) with the standard formula 
, but adjusted to satisfy the periodicity constraints of the model (10). Reestimation of M, µJ, J, and µk is more complicated. For a given dwell sequence, there are two independent sources of variance in the data: the variability of the step size, and the measurement noise. Thus, any parameter estimator should not only minimize the total data variance, in a sum of squares sense, but should also correctly partition this variance (SS) between the intrinsic step size variability (SSJ) and the measurement noise (SSM). If we were to consider a model with zero step variability, then SSJ would be zero, and only SSM would have to be minimized. Unfortunately, in the general case, SSJ and SSM cannot be minimized independently, because they are both functions of µk, and because the amplitudes of consecutive dwells are correlated through the step size distribution. Minimizing SS with respect to M, µJ, J, and µk simultaneously requires a nonlinear optimization with a large number of free parameters. As a compromise, we made a linear approximation: we first estimate µk and µJ given the previous M and J, and then estimate M and J, given the new µk and µJ.
The expression of SSJ may take different forms for different models. In the general case, the following issues must be considered:
All these difficulties are handled by our algorithm, but to avoid confusion, we will illustrate the calculations only for two simple cases that apply to many types of experimental data. For irreversible models with uniform step size (i.e., unimodal step size distribution), SSJ takes the following form:
 | (10) |
 | (11) |
This formulation makes the calculation less sensitive to mistakes in classifying a jump as corresponding to either a small or a large step.
The minimization of SSJ with respect to µp values and µJ (or µJ1 and µJ2) is done by imposing the conditions 
and 
, and solving the resulting matrix equation. Once µp and µJ are determined, the step variance is calculated as 
. The measurement error SSM can be calculated simply as 
, where the sum over all dwells (indexed by k) is in fact reduced to only one term, as only one dwell time interval includes the measurement yt. From this, the measurement variance is obtained as 
, where N is the number of data points.
|
MATERIALS AND METHODS |
|---|
). The simulated data were sampled like the experimental data at 0.5 s intervals. For those experiments designed to test the idealization algorithm as a function of SNR and motor step variability, the random number generator was initialized with the same seed, resulting in identical dwell sequences.
Stage-stepping data
These control data were obtained by attaching a fluorescent probe to the substrate and programming the microscope stage to move in a prescribed pattern, thus characterizing the instrumental resolution. A single Cy3 molecule was attached to a 20-mer dsDNA segment immobilized on a glass coverslip. The stage was translated by a 0.7 nm-resolution piezoelectric stage (8), according to stochastic dwell time sequences created by simulation of the model shown in Fig. 2 A that describes stalk-attached probes (Fig. 1 C1). Specifically, the model had one kinetic state per reaction unit and irreversible kinetics with forward rate constant kF = 0.25 s–1. The SNR was varied by changing the step size to 8, 13, and 30 nm. Despite the desired precision of the stage, both visual inspection and the idealization program show somewhat higher variability in the actual jump amplitude, possibly a result of limited digital-to-analog resolution of the driver. No significant baseline drift is visible. The data sets were recorded with a sampling time of 0.5 s, and were between 46 and 176 s long.
|
). Biotinylated actin filaments were immobilized on a glass surface coated with biotin-BSA and streptavidin. Myosin V was added to the sample flow chamber after actin immobilization, excess was washed off, and the surface was excited by objective-type TIR and the emission was imaged with a charge-coupled device camera. Each fluorescent spot had full width at half maximum of 
280 nm. Each spot contained 5,000–10,000 photons, so the two-dimensional Gaussian centroid could be fitted with 3 nm error in the peak position for typical spots, and 1.5 nm for brighter spots. The spots displayed quantal bleaching indicative of a single molecule. In the absence of ATP, the fluorescent spots were immobile. The addition of 300 nM ATP led to discernable steps, with the rate increasing with [ATP]. Only the steps of singly labeled myosins were analyzed.
Kinesin data
Human ubiquitous kinesin was labeled in the head region with a single Cy3 dye as described previously (9). Sea urchin axonemes (a microtubule rich structure) were immobilized onto the glass surface by flowing a suspension through a chamber. The chamber was then rinsed and perfused with kinesin. Positional stability measurements of labeled kinesin bound to axonemes in the absence of ATP showed that the axonemes were well attached to the glass and kinesins were strongly bound to the axonemes; 150 nM ATP caused the kinesins to walk along the axonemes.
|
RESULTS |
|---|
by changing the measurement noise M. As a convenient measure of motor step variability, we used the coefficient of variation 
, and varied it by changing the standard deviation of the motor step size J. The next step was to analyze control stage-stepping data, where there was only instrumental noise, but the transition points were known. Finally, we analyzed the behavior of myosin and kinesin as a full test, although there are no a priori models with which to compare the results. This last analysis confirmed the adequacy of the various assumptions made about the kinetic and noise models of experimental data.
Computer simulated data sets
We generated staircase data with irreversible and reversible models having uniform or alternating small and large steps. We have tested kinetic models of increasing complexity, but since similar results were obtained, we discuss here only the simpler models shown in Fig. 2, having one state per reaction unit. In all cases, we used a forward rate constant kF = 0.25 s–1, and a backward rate constant kB = 0.0 s–1 for the irreversible and kB = 0.05 s–1 for the reversible models. The SNR of the simulated data was varied between 1and 10. Except for data with alternating steps, the jump amplitude was µJ = 10 nm. For each SNR, CVJ was varied between 0% and 20%, by changing J between 0 and 2 nm. Note that µJ is not important in its absolute value, but only as relative to J and M. Idealization examples are presented in Fig. 3. Note that data generated with larger models (i.e., with two or more states per reaction unit) could be well idealized with a single state model, since the differences in kinetic complexity have only second order effects on idealization. The results obtained from irreversible models having uniform steps were further analyzed, since that is the default model used in the literature, and the results are presented in Figs. 5 and 6.
|
|
|
5, and only with a few mistakes at SNR as low as 2 (e.g., Fig. 3 A, trace 7, c and d). Note that since a model is fitted to an entire data set, a few errors in detection have little effect upon the final parameters. This level of performance can be achieved only when the intrinsic step variance is zero (traces 1, 4, and 7). Increasing the step variance (CVJ > 5–10%; traces 2, 3, 5, and 6) causes two kinds of errors. First, small jumps are undetected in the presence of wideband noise. This leads to an overestimation of µJ, since the next detected jump is two or more steps high, and the corresponding rate constants are smaller since transitions have been missed. Secondly, the visible jumps may be misclassified as jumps of order 1 when they were actually jumps of order 2 or higher, and vice versa. The question of jump order occurs regardless of SNR, but it can be treated by including all possible combinations in the reestimation procedure. Examples of this error are marked in trace 3, a and b, where the true jump orders are 3 and 1, but were incorrectly estimated as 2 and 2, respectively.
Note that since we knew these data were obtained with irreversible models, we imposed the constraint that all entries in the A matrix corresponding to backward steps were zero. This condition, when enforced at the start, propagates throughout the iterations. Without this constraint, especially at low SNR, noise artifacts are more likely to be misidentified as backward steps. We want to reiterate the value of a priori constraints: the more constraints, the more precise the results. The penalty is that the accuracy may suffer because the model is wrong. Modeling must be hierarchical and only expand in complexity when the likelihood indicates significant improvement in the fit.
Reversible models with uniform steps
Idealization results for simulated data are presented in Fig. 3 B, obtained using the model in Fig. 2 A. For the same SNR, reversible data are more difficult to idealize than irreversible data as they lack the constraint of no backward transitions. The result is that it is more likely to miss a transition (trace 3, a and b) or to mistake noise for a movement of the motor (trace 3, c). Missing transitions has the primary effect of underestimating rates, whereas mistaking noise spikes for motor transitions results in overestimated rates and possibly overestimated step variance. At SNR 5, reversible data were idealized without mistakes (traces 1 and 2).
Irreversible models with alternating small and large steps
Examples of idealization results for simulated data are presented in Fig. 3 C, obtained using the model in Fig. 2 B. The small steps are intrinsically more difficult to detect than the large steps. However, the constraint that the steps do not occur independently but must alternate allows reliable idealization of both small and large steps. The example traces shown in Fig. 3 C are idealized without mistake, even when the small step has a corresponding SNR = 2 (trace 4). Note that a simpler jump detector, such as the t-test, which uses only local data information (the means of two samples), may not be so successful at resolving small jumps.
Idealization of experimental data
Stage-stepping control data
The stage-stepping data were idealized using the model in Fig. 2 A, with results similar to the computer simulations, as shown in Fig. 4 A. Although these data were generated with zero step variance (J = 0.0 nm), we measured a finite jump variance illustrating the sensitivity of the method and the limitations of the instrument. The data with µJ = 30.0 (trace 1) and 13 nm (trace 2), with SNR 10, were idealized without error. In contrast, the idealization of 8 nm steps data (traces 3–6), with SNR 4, had a few mistakes because of the lower SNR and the more prominent experimental artifacts. Note that true error rates can only be known with computer simulated data, where the properties of the input data are known to arbitrary precision.
|
J = 1.86 nm.
Myosin motor data
Examples of myosin data are shown in Fig. 4 B, for alternating steps (the dye attached near the motor midpoint, traces 2 and 3) and for uniform steps (the dye attached near the head region, trace 4). We idealized with irreversible, single-state models, with alternating steps (Fig. 2 B), or with uniform steps (Fig. 2 A). As myosin's step size is bigger than kinesin's, the SNR in this case is better and the CVJ is smaller. Hence, all these traces were idealized without difficulty. For these particular examples, we estimated 39/35 nm alternating steps, and 72 nm uniform steps, in good agreement with the literature (8).
Effects of SNR and step variance—simulated data
The effects that the SNR and the step variance have on the accuracy of estimates are quantified in Fig. 5 A. Clearly, for zero step variance, the estimates of µJ, J, M and kF are excellent, even at SNR = 2 (Fig. 5 A, red lines). Estimates are still accurate for CVJ = 10% (Fig. 5 A, blue lines) and SNR = 2, with the exception of J. All parameters are significantly less accurate when the step variance is increased to 20% (Fig. 5 A, green lines), and the idealized dwell sequence has more mistakes. These results were obtained without accounting for all possible jump orders between consecutive dwells, as discussed, but even so, the errors were <2–3%. We advise that photon integration time should be reduced, despite an increase in measurement noise, as long as the SNR remains within the acceptable range. The increased time resolution reduces the number of missed events, whereas the global nature of the fit tends to remove the influence of the excess photon noise.
Effect of missed events
Amplitude discrimination improves with SNR. We recommend using an offline digital filter (e.g., the one in the QuB program), rather than increasing the photon integration time, since the software filter is reversible. The data should be resampled at the Nyquist limit after filtering (nominally two data points per cycle at the cutoff frequency) to satisfy the Markov assumption. Our algorithm is not designed to deal with higher order Markov processes, but extension is possible (25,53). The net effect of filtering is to improve the SNR at the expense of time resolution. The penalty is that there are more higher order jumps (missed events). Although higher order jumps are not a problem when the step variance CVJ = 0, they pose a serious challenge when CVJ > 10%.
We tested the effects of missed events on the idealization of staircase data with irreversible and uniform steps. Data simulated as in Fig. 3 A, with SNR 5 and CVJ = 0...20%, were downsampled by factors of 2, 4, 8, and 16, which progressively eliminated shorter dwells, and increased the fraction of higher order jumps, noted fH. We intentionally left the bandwidth constant to separate the effects of SNR—determined by bandwidth—from the effects of missed events. The dwells averaged 8 points in duration for the original data, but only 0.5 points when downsampled by a factor of 16. The results are shown in Fig. 5 B. For zero jump variance, all parameters were estimated accurately, even when the fraction of higher order jumps fH 0.7. For 10% jump variance, the estimates of µJ and kF are still good, but estimates of J and M are incorrect when fH is high. For CVJ = 20%, the estimates are not usable.
Statistical distribution of estimates
The estimates of µJ, J, and M have a Gaussian distribution, as illustrated in Fig. 6 A, for simulated irreversible staircase data with uniform steps, in this case with SNR 5 and 10, and CVJ = 10%. The width of the distribution depends on SNR and CVJ (results not shown for the latter). No outliers were observed for any of the estimated parameters. As previously discussed and shown in Fig. 5 A, the distributions of J and M are biased toward larger and smaller values, respectively, as the SNR decreases. The estimated parameters are well correlated with the true parameters used in simulation, as shown in Fig. 6 B for µJ, J, and kF.
Convergence of the algorithm
The algorithm is always started with an initial guess of the parameters, and it is desirable that different initial guesses lead to the same solution, i.e., there is a global optimum. We tested the influence of the starting values, as shown in Fig. 7. In the first experiment, all parameters were chosen from a random hypercube in parameter space whose range was ±20% of the true values. With SNR 5 and CVJ = 10%, the program converged to the correct answer in 5–10 iterations (Fig. 7 A). The rate of convergence increased with SNR, and decreased with CVJ. The algorithm did not converge at SNR 1. In general, the algorithm converged monotonically.
|
TOP
ABSTRACT
INTRODUCTION
10%.
