format: text

text:
# A Simple Kinetic Model Describes the Processivity of Myosin-V

## INTRODUCTION

Various classes of enzymes, usually termed motor proteins, play important roles in biological processes such as cellular transport, cell division, muscle function and genetic transcription (see, e.g., Lodish et al., 1995). What we may term translocatory motor proteins (in contrast to rotary motor proteins) are epitomized by kinesins, dyneins, myosins, and DNA and RNA polymerases that move under loads along polar linear tracks such as microtubules, actin filaments, and double-stranded DNA, the motion being fueled by the hydrolysis of ATP or related reactions.

Motor proteins may work collectively in large groups, like myosin in muscles, or they may operate individually as do most microtubule-based kinesin and dynein molecules (Leibler and Huse, 1993; Howard, 2001). Those motor proteins that function collectively are typically nonprocessive, i.e., they make at most one mechanical step along their tracks during a catalytic cycle before detaching from the track. On the other hand, individual motors that move vesicles over long distances (up to several microns) need to stay bound to their tracks over many catalytic cycles: such motors are processive. For example, conventional kinesin motors can walk along microtubules taking a hundred or more 8.2 nm steps before dissociating (Howard et al., 1989; Block et al., 1990; Vale et al., 1996).

Recently, single-molecule experiments by Mehta et al. (1999; Mehta, 2001), Rief et al. (2000), Sakamoto et al. (2000), Rock et al. (2001), Veigel et al. (2002) and Nishikawa et al. (2002) have demonstrated that myosin-V and myosin-VI, in contrast to the behavior of other members of the myosin superfamily (Howard, 2001), are also efficient processive molecular motors. Here we will focus on the dynamics of myosin-V.

Myosin-V is a dimeric, two-headed molecule that in the presence of actin readily hydrolyzes ATP to produce ADP and P$`_\text{i}`$ (Mehta, 2001). Kinetic experiments in bulk solution (De La Cruz et al., 1999, 2000; Mehta, 2001) have demonstrated that release of ADP is the rate-limiting step in the actin-activated ATPase cycle. Under conditions of limiting ATP the kinetically prevalent state appears to have both head domains bound to the actin filament as captured in electron micrographs by Walker et al., (2000); but more generally, see the discussions in Mehta (2001) and De La Cruz et al. (2001).

Optical traps equipped with electronic feedback mechanisms have provided valuable information regarding the dynamics of individual myosin-V molecules under low load (see Mehta, 2001). The experiments allow one to monitor the displacement, $`x(t)`$, of a single molecule as a function of the time $`t`$ under different concentrations of ATP, etc., while maintaining a steady external load, $`F`$, which opposes the directed motion of the motor.

The principle experimental findings can be summarized as follows: (i) Myosin-V moves along actin filaments towards the plus or barbed end taking large steps of size averaging 35-38 nm (Mehta et al., 1999) approximating the 37 nm pseudo-repeat of the actin filament (Bray, 2001); (ii) The stepping dynamics depends strongly on the ATP concentration: thus, the mean dwell time $`\tau (F,\text{[ATP]})`$ observed between successive steps (preceding a forward step) at low \[ATP\] (=1 $`\mu `$M) hardly varies with the external load, while under saturating conditions (\[ATP\] $``$ 2 mM) the mean dwell time grows rapidly as $`F`$ approaches the stall force, $`F_S=3.0\pm 0.3`$ pN (at which, on average, the motor just fails to progress); (iii) The overall stepping rate or mean velocity

$$V(F,\text{[ATP]})dx(t)/dt$$
(1)
follows a Michaelis-Menten form in that it is proportional to \[ATP\] at low concentrations but becomes independent of \[ATP\] under saturating conditions; (iv) Tight coupling between chemical and mechanical cycles is valid, with one ATP molecule hydrolyzed per individual myosin-V forward step along an actin filament; but (v) in contrast to the dynamics of conventional kinesin (Coppin et al., 1997; Visscher et al., 1999), myosin-V under load not infrequently exhibits sequences of two or three reverse or backward steps; and, finally, (vi) the addition of ADP to the in vitro solution significantly reduces the turnover rate of ATP (as is to be expected); moreover, the inhibitory effect of ADP scales with the concentration of ATP (Rief et al., 2000) (and even when the mean stepping rate is reduced two-fold the distribution of dwell periods is unaltered).

The growing quantity of information concerning myosin-V has naturally stimulated theoretical discussions of the dynamics. Several models have been proposed and are reviewed by Mehta (2001). In particular, in order to provide an explanation of the observed load-dependence of the processivity, the mean dwell time at temperature $`T`$ has been modeled phenomenologically \[following a proposal of Wang et al. (1998)\] as a sum of two terms, namely

$$\tau (F)=\tau _1+\tau _2\mathrm{exp}(Fd^{}/k_BT),$$
(2)
corresponding, respectively, to putative force-independent and force dependent transitions. It is natural to expect here that $`d^{}`$ corresponds to the observed step size $`d36`$ nm (Fisher and Kolomeisky, 1999; Fisher and Kolomeisky, 2001; Hille, 2001). However, fitting the experimental data of Mehta et al. (1999), which is displayed in Fig. 2 below, necessitates an (effective) step size $`d^{}`$ of 10-15 nm, which is only 30-40 $`\%`$ of the actual step size. This discrepancy is rationalized by asserting that $`d^{}`$ is some “characteristic distance over which load affects the catalysis rate.” Furthermore, this approach fails to account clearly for the observed stalling of the motors at $`F_S3.0`$ pN. Clearly, a more soundly based quantitative theory for processivity of myosin-V seems called for in order to satisfactorily describe the currently available data and to provide testable predictions. The present article aims to meet these requirements.

We present a theoretical analysis of the dynamics of myosin-V using simple, discrete-state stochastic models which have recently been developed and analyzed in detail by Kolomeisky and Widom (1998), Fisher and Kolomeisky (1999a, b; 2001; 2002) and Kolomeisky and Fisher (2000a, b); for brevity these articles will be referenced below as FK’01, KF’00a, etc. This approach has been used successfully in FK’01 to analyze the extensive experimental data on the dynamics of single conventional kinesin molecules moving in vitro along microtubules obtained by Visscher et al. (1999) and Schnitzer et al. (2000). We will demonstrate that most of the currently available experimental data on the processivity of myosin-V can be well accounted for by the simplest ($`N=2`$)-state model embodying a theoretical picture in satisfactory accord with other kinetic and structural experiments. Our treatment also provides specific predictions for as yet unexplored features of myosin-V dynamics that can be tested experimentally and should uncover further details of the stepping mechanism.

## Theoretical Approach

For completeness we first outline briefly the class of stochastic models used in our analysis and the explicit analytical results available for them. In the simplest periodic sequential kinetic model, illustrated schematically in Fig. 1, the protein motor is viewed as moving along a linear periodic track and binding at specific sites located at $`x=ld`$ ($`l=0,\pm 1,\pm 2,\mathrm{}`$) where $`d`$ is a fixed step distance. In a first treatment of myosin-V on actin filaments we may adopt the observed mean value, $`\overline{d}36`$ nm, corresponding to to the helix repeat distance (Bray, 2001). However, as discussed further below, the analysis can be extended to take account of the variations in the individual step sizes seen in the data for myosin-V (Mehta et al., 1999; Rief et al., 2000; Walker et al., 2000; Veigel et al., 2002): the variation seems primarily to result from binding on actin monomers (at spacing 5.5 nm) adjacent to the main 6.5-monomer helix repeat distance (Steffen et al., 2001).

The basic model then supposes that in a catalytic cycle which translocates a motor from binding site $`l`$ to $`l+1`$ the protein undergoes $`N`$ intermediate biochemical transitions from states $`j_l=0_l`$ to $`1_l`$ to $`2_l`$ $`\mathrm{}`$ to $`(N1)_l`$ to $`N_l0_{l+1}`$. Kinetic rates $`u_j`$ and $`w_j`$ are associated with the transitions from state $`j_l`$ forwards to state $`(j+1)_l`$ and backwards to state $`(j1)_l`$, respectively. The state $`0_l`$ represents the motor tightly bound at site $`l`$ in the absence of fuel molecules — ATP in the case of myosin-V. Binding of a fuel molecule is represented by the transition $`0_l1_l`$, unbinding by $`1_l0_l`$. Subsequent hydrolysis and release of products occur in the forward transitions $`1_l2_l\mathrm{}`$. But it is important to note that backwards intermediate transitions and whole steps (possibly associated with reverse hydrolysis) are allowed and observed experimentally.

For this model the mean velocity, $`V(\{u_j,w_j\})`$, (see Eq. 1) may be expressed exactly in a closed analytic form in terms of the rate constants $`\{u_j,w_j\}`$ for any value of $`N`$ (FK’99). Furthermore, similar explicit formulae are available for the dispersion (or effective diffusion constant) of the motion, defined by

$$D=D(\{u_j,w_j\})=\frac{1}{2}\underset{t\mathrm{}}{lim}\frac{d}{dt}\left[x^2(t)x(t)^2\right].$$
(3)
This measures the statistical deviation of the motor trajectories from uniform motion at constant velocity. The knowledge of both the velocity $`V`$ and the dispersion $`D`$, conveniently combined in terms of randomness (Svoboda et al., 1994)

$$r=2D/Vd,$$
(4)
serves to set bounds on $`N`$ via a determination of the number of rate-limiting kinetic biochemical transitions and thus yields valuable information regarding the mechanism of processivity (Visscher et al., 1999; KF’00a; FK’01; FK’02; Koza, 2002).

To account properly for the externally imposed force, $`F`$, it is essential (FK’99; FK’01) to introduce load distribution factors, $`\theta _j^+`$ and $`\theta _j^{}`$ (for $`j=0,1,\mathrm{},N1`$). Then the transition rates may be taken to vary as

$`u_ju_j(F)`$ $`=`$ $`u_j^0\mathrm{exp}(\theta _j^+Fd/k_BT),`$
$`w_jw_j(F)`$ $`=`$ $`w_j^0\mathrm{exp}(\theta _j^{}Fd/k_BT),`$ (5)
where the most reasonable requirement (FK’99; FK’01; Hille, 2001) is

$$\underset{j=0}{\overset{N1}{}}(\theta _j^++\theta _j^{})=1,$$
(6)
which implies that the condition of stall reflects stochastic quasiequilibrium amongst the (on-pathway) intermediate mechanochemical states. Indeed, these expressions embody a picture of load-dependent activation barriers for forward and reverse rates between intermediate states $`j`$ that lie on a multidimensional reaction pathway. The load distribution factors $`\theta _j^\pm `$ provide significant mechanochemical information since they embody a projection of the valleys and cols (or passes) of the reaction pathway onto the force axis, which we suppose is parallel to the motor track. Thus one may identify substeps of magnitude

$$d_j=(\theta _j^++\theta _{j+1}^{})d,$$
(7)
between motor states $`j_l`$ and $`(j+1)_l`$. If the spatial fluctuations of the center of force of the motor in the intermediate states $`j_l`$ and $`(j+1)_l`$ are sufficiently small relative to $`d_j`$ one may hope to identify this substep in suitably averaged traces $`x(t)`$ of individual motor motions (FK’02).

For the present purpose we note that the explicit expressions for the mean velocity, $`V`$, for general $`N`$ lead to a simple relation for the stalling force as defined by $`V(FF_S)0`$, namely,

$$F_S=\frac{k_BT}{d}\mathrm{ln}\left[\underset{j=0}{\overset{N1}{}}(u_j^0/w_j^0)\right];$$
(8)
see FK’99.

The $`N`$-state periodic kinetic model presented in Fig. 1 is, mathematically, an example of the general one-dimensional nearest-neighbor random hopping model for which first-passage questions have been much studied: see van Kampen (1992). Of particular interest here are the so-called splitting probabilities and mean first-passage times. Specifically, in order to analyze observations of motor-protein dwell times, we need the “single-step forward splitting probability,” $`\pi _+(\{u_j,w_j\})`$, defined as the probability that a motor starting at site $`l`$ will arrive at site $`l+1`$ without having undergone sufficiently many intermediate reverse transitions to complete a full backwards step from $`l`$ to site $`l1`$. The corresponding conditional mean single-step first-passage time, $`\tau _+(\{u_j,w_j\})`$, then represents the average time a motor spends at site $`l`$ before leaving and making a forward step to site $`l+1`$. Because of the periodic structure of the $`N`$-state model of Fig. 1 the (rather elaborate) expressions developed by van Kampen (1992) can be simplified considerably even for general $`N`$ (Kolomeisky and Fisher, to be published). Here we quote the simplest $`N=2`$ results which will suffice for our present purposes, namely, for the mean forward-step dwell time,

$$\tau _+=(u_0+u_1+w_0+w_1)/(u_0u_1+w_0w_1),$$
(9)
while the fraction of backward (or reverse) steps is

$$\pi _{}=1\pi _+=w_0w_1/(u_0u_1+w_0w_1).$$
(10)
Finally we mention that the basic model exhibited in Fig. 1 can be extended in various ways while still retaining explicit expressions for $`V`$, $`D`$, etc. In particular, one may allow for detachments or “death” rates, $`\delta _j`$, from the various motor states and for branching (KF’00a), for parallel site-to-site “jumping” (KF’00a), for parallel biochemical processes (Kolomeisky, 2001), and for waiting time distributions and the associated degrees of mechanicity, $`M_j^\pm `$, of the various intermediate processes (FK’99; KF’00b; FK’01). However, the range of observational data so far obtained for myosin-V (unlike that known for kinesin) does not yet warrant consideration of these extensions.

## Analysis of Myosin-V Data

The bulk-solution kinetic data on myosin-V ATPase activity indicate that at least two processes, namely, ATP binding and ADP release, should be taken into account in analyzing the motility (De La Cruz, 1999; Mehta, 2001). While recognizing that a more complete description may require further intermediate states, it is appropriate, therefore to consider first the simplest ($`N=2`$)-state model. Then, as indicated above, the states $`j=0`$ correspond to a myosin-V molecule bound to the actin filament in the absence of ATP — presumably with both heads attached, one behind the other (Walker et al., 2000) — while $`j=1`$ labels myosin-actin complexes with bound ADP. Thus, in the scheme advanced in Fig. 6 of Mehta (2001), the first and last configurations correspond to $`j=0`$ while the four intermediate states are gathered into $`j=1`$; in Fig. 4 of Rief et al. (2000) the $`j=0`$ state corresponds to that labelled V; in Fig. 9 of De La Cruz et al. (2001) the second configuration corresponds to $`j=0`$, the remaining three to $`j=1`$.

It now follows that the forward ATP-binding rate should take the form $`u_0^0=k_0^0`$\[ATP\], where the superscripts 0 denote the limit of zero load: see Eq. 5. On the other hand, the reverse unbinding rate, $`w_1`$, and the forward, ADP release rate, $`u_1`$, should be independent of \[ATP\], but, of course, may depend on $`F`$.

According to standard chemical kinetic arguments, the backward rate $`w_0`$ should, in principle, be proportional to \[ADP\]; and, indeed, the concentration of P$`_\text{i}`$ should also play a role. Note, particularly, in this connection the high affinity of ADP for actomyosin which, as discussed by Mehta (2001) and De La Cruz et al. (2000), had led to significant discrepancies in estimates of steady-state cycling rates. The detailed measurements (Mehta et al., 1999; Rief et al., 2000) have, therefore, been performed with the aid of an ATP regeneration system \[as previously adopted in the kinesin experiments of Visscher et al. (1999)\]. In such a set-up neither the concentration of ADP, nor the that of P$`_\text{i}`$, is monitored. While experiments that do control \[ADP\] and \[ P$`_\text{i}`$\] separately are much to be desired, in their absence we are forced (as in FK’01) to model the ATP regeneration scheme more or less phenomenologically. Thus if (a) we suppose $`w_0^0=k_0^{}\text{[ATP]}^\alpha `$ (which amounts to \[ADP\] $``$ \[ATP\]α), (b) recall that the stall force, $`F_S`$, is given by Eq. 8, and (c) note that the current experimental observations reveal no significant dependence of $`F_S`$ on \[ATP\] (Mehta, 2001), we are led to adopt $`\alpha =1`$. Indeed, in light of the use of ATP-regeneration in the experiments, the proportionality of \[ADP\], and hence of $`w_0`$, to \[ATP\] at low concentrations is to be expected: see also FK’01. It should be remarked, however, that the details of our description of the ATP regeneration scheme play only a minor role in fitting the myosin-V processivity data.

Now in many previous experimental studies of processive motor proteins the mean velocities, $`V(\text{[ATP]},F)`$, have been measured and reported. Such observations must, at least in principle, include some fraction of backward or reverse steps, especially at large loads approaching stall. However, in their experiments on myosin-V (Mehta et al., 1999; Rief et al., 2000) the authors opted to measure only dwell times, separating adjacent steps of mean size $`d36`$ nm (Mehta, 2001), preceding forward steps. Thus their reported dwell times, $`\tau (\text{[ATP]},F)`$, as plotted in Fig. 2A, do not precisely correspond to an “overall mean step time,” say $`\overline{\tau }`$, related to the mean velocity simply via $`\overline{\tau }=d/V`$ — although at low loads, where the fraction of reverse steps is small, $`\overline{\tau }`$ should provide a good approximation; but under near stall conditions, when $`V0`$, the overall mean step time, $`\overline{\tau }`$, diverges to infinity whereas the dwell times $`\tau (F`$$``$$`F_S)`$ remain bounded. Rather, we identify the observed dwell times with the conditional single-step mean first-passage times, $`\tau _+`$, identified above: see Eq. 9. Accordingly, our analysis of the myosin-V data is based upon the expression

$$\tau (F,\text{[}ATP])=\frac{k_0^0\text{[ATP]}e^{\theta _0^+Fd/k_BT}+u_1^0e^{\theta _1^+Fd/k_BT}+k_0^{}\text{[ATP]}e^{\theta _0^{}Fd/k_BT}+w_1^0e^{\theta _1^{}Fd/k_BT}}{k_0^0\text{[ATP]}e^{\theta _0^+Fd/k_BT}u_1^0e^{\theta _1^+Fd/k_BT}+k_0^{}\text{[ATP]}e^{\theta _0^{}Fd/k_BT}w_1^0e^{\theta _1^{}Fd/k_BT}},$$
(11)
following from Eqs. 5, 6, and 9, with $`d=\overline{d}=36`$ nm.

Then, by systematically exploring the full seven-dimensional parameter space specified by ($`k_0^0,\mathrm{},\theta _0^{}`$) we find that the observed stall force, $`F_S`$, and the dynamics of myosin-V as a function of \[ATP\] and of the load, $`F`$, up to $`F_S`$, are well described by the rates

$`k_0^0=0.70\pm 0.10\mu \text{M}^1\text{s}^1,`$ $`u_1^0=12.0\pm 1.0\text{ s}^1,`$ (12)
$`k_0^{}=(5.0\pm 0.5)\times 10^6\mu \text{M}^1\text{s}^1,`$ $`w_1^0=(6.0\pm 0.5)\times 10^6\text{ s}^1,`$
and the load-distribution factors

$`\theta _0^+=0.010\pm 0.010,`$ $`\theta _1^+=0.045\pm 0.010,`$ (13)
$`\theta _0^{}=0.580\pm 0.010,`$ $`\theta _1^{}=0.385\pm 0.010.`$
It should be noted that consideration of the limits of low and high \[ATP\] and low and high loads, confirm a fair degree of independence of the various fitting parameters. The uncertainties indicated in Eqs. 12 and 13 correspond to the ranges of acceptable fits to the processivity data while constraining the other parameters appropriately. The central values yield the fits presented in Fig. 2 as solid curves.

In respect to our fits for $`k_0^0`$ and $`u_1^0`$ note that the bulk solution kinetic experiments yield an ATP binding rate constant (corresponding to $`k_0^0`$) between 0.7 and 1.6 $`\mu `$M-1s-1, while the ADP-release rate (corresponding to $`u_1^0`$) is about 12-16 s-1 (Mehta, 2001; De La Cruz et al., 1999). The agreement is clearly most satisfactory.

## Discussion

### Mean Velocity and Load Dependence

The quality of the fits in Fig. 2 ensures that the observed (approximate) Michaelis-Menten behavior is respected. Indeed, using the rate and load-distribution parameters in Eqs. 12 and 13 and previous theory (e.g., FK’01) enables us to predict the variation of the mean velocity, $`V`$, with $`F`$ and \[ATP\]: see the solid curves Fig. 3. Evidently, the stall force of about 3 pN seen in the experiments is reproduced. Note also, from the dwell-time data imposed on the predictions in Fig. 3 using $`Vd/\tau `$, that, as anticipated in the discussion before Eq. 11, the approximation $`\tau \overline{\tau }d/V`$ is valid for small loads (up to $`F2.5`$ pN). Indeed, from Eq. 10 (with Eqs. 5, 6, 12, and 13) one finds that the fraction of reverse steps is negligible until $`F2.5`$ pN.

### Load Dependence of Rates

It is notable from Eq. 13 that within the fitting uncertainties there is essentially no load-dependence to the binding of ATP to the myosin-V-actin complex, i.e., $`\theta _0^+0`$; see also Mehta (2001). This contrasts strongly with the properties of conventional kinesin moving on a microtubule where $`\theta _0^+0.13`$ was found in FK’01 for both $`N=2`$ and $`N=4`$ fits. This lack of load-dependence on binding ATP to actin-myosin accounts for the fact that the dwell time remains constant at saturating ATP conditions up to $`F2.3`$ pN: see Fig. 2A.

Nevertheless the other transitions are load dependent with ADP release bearing a modest ($`5\%`$) fraction of the dependence. In parallel to kinesin, however — see FK’01, the reverse transitions carry most of the load-dependence. Indeed, the load distribution pattern (FK’01) for myosin-V is close to a featureless descending ramp. Note that this result is in striking contrast to the implications of the phenomenological expression Eq. 1 which suggests that only forward (i.e., binding and/or hydrolysis) processes need be considered and could exhibit significant load-dependence. Indeed, our analysis indicates that at least three biochemical transitions in the actin-myosin-V ATPase cycle are load-dependent whereas Eq. 1 entails only a single load-dependent process. It seems that this difference is the main reason why fits for the “characteristic distance” $`d^{}`$ in Eq. 1 differ so markedly from the true mean step size $`\overline{d}36`$ nm. Since our analysis recognizes reverse transitions, which, by the fits, occur at a non-vanishing rate that is enhanced under load \[see, again, Eq. 10\], an explanation is provided for the observation of more frequent backward steps in myosin-V at high loads (Rief et al., 2000). Our treatment also provides a basis for a quantitative discussion of the ADP inhibition effect which it would be instructive to explore further experimentally.

### Substeps

A striking feature of the data of Rief et al. (2000) is the observation of “half steps” under high loads ($`2`$ pN). From the published traces the steps appear to correspond to an intermediate state with a mean center of force lying a distance, say $`d_{1/2}`$ forward from the bound-state ($`j=0`$) center with $`d_{1/2}/d0.48\pm 0.04`$. On the other hand, Eq. 7 and the load distribution factors in Eq. 13 indicate a substep with $`d_0/d0.38\pm 0.03`$ (corresponding to $`d_01314`$ nm). Rief et al. (2000) suggest that these “half steps” (always followed by a complementary forward or backward step to complete a movement with $`\mathrm{\Delta }x=d`$ or 0) reflect an “off-pathway state” because they remain rare even under the high loads that uncover their presence. While this suggestion seems most reasonable on the available evidence, our analysis suggests that the half steps might possibly represent genuine substeps (lying on or close to the main reaction pathway), which appear stochastically under high loads when the forward rates, $`u_1(F)`$, have been slowed down while the reverse rates $`w_1(F)`$ are significantly enhanced.

In other experiments Veigel et al. (2002) observed attachments of single myosin-V molecules to an actin filament (stretched between two optically trapped beads) at \[ATP\]=100 $`\mu `$M. After some of the attachment events, “staircases” of from two or three to a dozen forward steps were seen of mean size 36 nm; the staircases typically terminated in an effective stall (signaled by interspersed forward and backward steps) before detachment from the filament: see Fig. 3 of Veigel et al. However, the authors concluded that the first step in each staircase had a mean size of only $`d_1=26.2\pm 2.3`$ nm (similar to the amplitude of isolated attachment events lacking any subsequent steps). A similar initial unitary step of $`20`$ nm was seen in experiments by Moore et al. (2001) on heavy meromyosin-like fragments of myosin-V. These displacements were identified (in both articles) as a “working stroke”; and Veigel et al. saw a comparable step of $`21`$ nm in attachment events of a single-headed recombinant myosin-V. Furthermore, Veigel et al. in their Fig. 5A, report stiffness measurements (using a sinusoidal driving force) which revealed low-stiffness intervals of variable durations (longer at higher loads): the mid-positions of these intervals was about 20 nm further along the actin filament than the preceding higher-stiffness intervals, a displacement similar to the initial “working stroke.”

In our formulation and fits using a single intermediate mechanochemical state prior to completion of a full ($`d=36`$ nm) step, such a $`d_1`$ should as the notation chosen suggests, correspond to a $`dd_022`$ nm “substep.” The agreement of these various findings (within the combined experimental and fitting uncertainties) appears to lend support to our values for the load-distribution factors $`\theta _j^\pm `$. However, corresponding substeps have not been identified at low loads by Mehta, Rief and coworkers. Nevertheless, a detailed examination of the sample stepping records for \[ATP\] = 2 mM and $`F`$ = 1 pN presented in Fig. 2A of Rief *et al.* (2000) reveals plausible indications of substeps in 13 to 16 of a total of around 32 “full” steps of $`36`$ nm, some of the substeps appearing to have dwell times as long as 0.1 - 0.2 s. More favorable conditions for detecting the predicted substeps and checking their dwell times should be realized at low loads and \[ATP\] $`10\mu `$M (which corresponds roughly to the effective Michaelis-Menten concentration, $`K_M`$: see Rief *et al.* (2000) and Fig. 3). Such data not consistent with the present predictions might require the introduction of waiting-time distributions (KF’00b): see also the remarks below concerning randomness.

### Variability of Step Sizes

The fits to the data so far described have utilized a fixed step size, $`d`$, taken equal to the observed mean step size $`\overline{d}36`$ nm that corresponds closely, as mentioned above, to the known (half) repeat distance of the actin filament double helix (Bray, 2001). But separate single-molecule experiments by Steffen et al. (2001) using myosin-S1 motor domains indicate “target zones” for binding to the filament consisting of three adjacent accessible actin monomers at spacings $`\mathrm{\Delta }d5.5`$ nm, the active zones repeating along the filament helices at $``$36 nm intervals. Furthermore, the processivity data for myosin-V reveal significant variations in individual step sizes about the mean, $`\overline{d}`$. The observations \[see: Mehta et al. (1999) Table 1; Rief et al. (2000) Fig. 2B; Walker et al. (2000) Fig. 2; Veigel et al. (2002) Fig. 3b\] are consistent with about 60$`\%`$ of the steps being of size $`d_{(0)}=36`$ nm while 20$`\%`$ each are of sizes $`d_{(\pm )}=d_{(0)}\pm \mathrm{\Delta }d=41.5`$ and 30.5 nm; only a few percent of longer or shorter steps appear. The fact that 40-45$`\%`$ of the observed steps deviate from $`d_{(0)}=36`$ nm raises the possibility that our fits using a unitary step size might be misleading or especially sensitive to the spread in sizes.

To address this issue note, first, that steps of distinct sizes, say $`d_{(k)}`$, should be expected to have different mean dwell times: an ideal set of experimental observations would, then, report the corresponding $`\tau _{(k)}(F,\text{[ATP]})`$ and their probabilities, say $`p_k`$. An analysis using Eq. 11 with $`d`$ replaced by $`d_{(k)}`$, etc., could subsequently be performed for each set and might possibly prove revealing. To a leading approximation one may suppose the various dwell times will be independent: in that case, the overall mean dwell time should be given by

$$\tau =\underset{k}{}p_k\tau _{(k)}.$$
(14)
More realistically, however, if the target-zone picture is valid, there will be correlations between successive steps: thus on average a short step, say of size $`d_{()}`$, must be followed immediately by a longer step, of size $`d_{(0)}`$ or $`d_{(+)}`$, and vice-versa. In principle, such correlations are open to observation and one might, indeed, expect the dwell times to depend on the size of the previous step, say $`d_{(k)}^{}`$, as well as on the step to be made. Theoretically the situation could clearly be modeled by a Markov process. \[See, e.g., Steffen et al. (2001).\]

In the absence of such more detailed observations, however, we may test the sensitivity of our fits by further exploratory calculations. As an extreme case, suppose 50$`\%`$ of the steps are of magnitude $`d_{(+)}=41.5`$ nm and 50$`\%`$ of size $`d_{()}=30.5`$ nm. How would the fits change from those assuming a unitary step $`d_{(0)}=36`$ nm? An answer is displayed by the dashed curves plotted in Fig. 2. These have been obtained by using Eq. 14 with $`p_+=p_{}=1/2`$ and computing $`\tau _{(+)}`$ and $`\tau _{()}`$ from Eq. 11 using $`d_{(+)}`$ and $`d_{()}`$ together with the same zero-load rates and load distributions factors given in Eqs. 12 and 13. As evident from Fig. 2, there is no significant change in the quality of the fits — even though it would not be unreasonable to suppose that the rates and load factors might have some dependence on the $`\pm 15\%`$ changes in step-size. One might say that “the averages win out” — a not unexpected conclusion.

In fact we may go further and study the effects of correlated step sizes by utilizing the expressions for $`N`$-state periodic models (KF’03) with $`N`$ an integral multiple of $`N_0`$, the number of intermediate states in the basic catalytic cycle. In our analysis we have $`N_0=2`$ and so can utilize an $`N=2+2=4`$ periodic system to describe alternating long and short steps of sizes $`d_{(+)}`$ and $`d_{()}`$ (with, of course, the same previous average step size $`d_{(0)}`$). If we again use the zero-load rates and distribution factors in Eqs. 12 and 13, and compute the mean velocity as a function of load, we obtain the dashed curves presented in Fig. 3. Once more the deviations from the $`d=\overline{d}`$ results are negligible at loads $`F<2`$ pN, while at higher loads sufficiently precise data might reveal discrepancies.

We conclude, therefore, that the consequences of replacing a distribution of step sizes by the mean $`\overline{d}`$ are not significant at current levels of experimental precision. Conversely, unless fairly precise experimental data can be obtained that are categorized by step length, there may be little more that can be reliably determined by fitting such observations.

### Randomness

As mentioned previously, the fluctuation statistics of motor motion are effectively captured in the randomness parameter, $`r`$, as defined in Eq. 4. The fits presented in Eqs. 12 and 13 suffice to predict the variation of $`r`$ with \[ATP\] under various loads (or vice versa) assuming that all the rate processes may be adequately represented as standard kinetic transitions: see FK’99 and KF’00. The corresponding predictions for $`r(\text{[ATP]})`$ are presented in Fig. 4 for loads $`F=0.4`$ and 2.5 pN. At low \[ATP\] the randomness is close to unity indicating that only one rate-limiting process is effective in this concentration range. However, under a low load a marked dip to $`r0.5`$ occurs around $`\text{[ATP]}=1020\mu `$M: this, in turn, is indicative of two competing rate processes that both play a role in this “crossover” regime. On the other hand, at high loads that approach stall, $`r`$ rises rapidly above unity; however, this is primarily a consequence of the vanishing of the velocity $`V`$ when $`FF_S`$ since $`r`$ must then diverge: see also Fig. 5(B), below.

It must be noted, however, that the analogous predictions, on the basis of an ($`N=2`$)-state kinetic model, for the the randomness of kinesin are not supported by the data of Visscher et al. (1999). Rather, for low loads and $`\text{[ATP]}30\mu `$M, the randomness falls rapidly and remains below 0.5 up to saturation concentrations: because of the bound $`r1/N`$ (FK’99, Koza, 2002), this is inconsistent with a kinetic description. Thus the data for conventional kinesin demand $`N=4`$ (or more) states (in accord with the usual biochemical picture of ATP hydrolysis). Alternatively, and, in light of certain experiments (Nishiyama et al., 2001), possibly more realistically, one may invoke a waiting-time distribution to describe the process of hydrolysis and ADP release with a mechanicity $`M_10.6`$ (KF’00, FK’01). Thus measurements of $`r(F,\text{[ATP]})`$ for myosin-V might well prove equally revealing of the mechanism by failing to verify the behavior predicted by Fig. 4!

### Reverse or Assisting Loads

Another interesting and potentially instructive set of predictions can be advanced for the behavior under negative or assisting loads, $`F<0`$. Such experiments have been performed for kinesin by Coppin et al. (1997). Although their data posed certain problems (in particular, a significantly lower overall processivity under low loads) the same load distribution factors (and similar rates) provided a not unreasonable ($`N=2`$) fit (FK’01) simply by extending the analog of Eq. 11 to negative values of $`F`$. The corresponding predictions for the dwell time and for the randomness as a function of $`F`$, extending down to $`3`$ pN, are displayed in Fig. 5.

A caveat must, however, again be raised in light of subsequent experiments on kinesin by Block (2001) and coworkers. The validity of the extension of Eq. 11 to negative $`F`$ clearly rests on a mechanistic/geometric assumption, namely, that changing abruptly the direction at which the coiled-coil myosin tail leaves the junction with two heads (or motor domains), i.e., from trailing upwards and backward ($`F>0`$) to pulling upwards and forward ($`F<0`$) does not result in a corresponding abrupt change in the mechanics of ATP binding, unbinding, or hydrolysis. If the junction were a perfect universal swivel joint, then as $`F`$, which is just the component of the total load force, say $`\stackrel{}{F}`$, parallel to the track, passes through zero the stresses and strains within motor should, indeed, vary smoothly. However, the junction cannot be totally torsion free and if, for example, the tail were to rest against part of the head in one configuration but become dissociated in the other, then the smoothness assumption embodied in Eqs. 5 would fail. Indeed, just such an abrupt change of behavior has since been found by Block and coworkers (2001) for kinesin. Clearly, comparable experiments on myosin-V are desirable and should prove informative.

## Conclusions

In summary, we have presented a simple two-state stochastic model, with allowance for fluctuating step sizes, which describes well essentially all the available experimental data on single-molecule myosin-V processivity. It reveals that ATP binding is load-independent, while ADP release is weakly load-dependent; but (as for kinesin) the loading forces strongly affect the reverse transition rates. Our analysis is consistent with the observation of tight coupling between catalytic cycles and mechanical steps, i.e., one ATP molecule is consumed per individual step, and with ATP binding and ADP release rates measured in bulk solution. It also indicates that an intermediate myosin-ADP-actin complex has its center of force advanced by 13-14 nm forward from the position prior to ATP binding, in reasonable agreement with various observations indicating a subsequent “working stroke” of around 22 nm. We have discussed specific predictions for the dwell times, mean velocity, and randomness of myosin-V motors in various experimental regimes including the imposition of assisting loads. Further experiments are needed in order to investigate the validity of our theoretical description and to uncover other mechanochemical features of myosin-V.

## Acknowledgments

Discussions with Jonathon Howard and Steven Block concerning kinesin and insightful comments from John Sleep and the referees concerning myosin-V have been appreciated. A.B.K. acknowledges the financial support of the Camille and Henry Dreyfus New Faculty Awards Program (under Grant NF-00-056). M.E.F. is grateful to the National Science Foundation for support (under Grant No. CHE-99-81772).

FIGURE LEGENDS

FIGURE 1. Specification of the simplest $`N`$-state periodic stochastic model. A motor in state $`j_l`$ can undertake a forward transition at rate $`u_j`$ or it can make a backward transition at rate $`w_j`$. The bound state $`N_l`$ is identified with $`0_{l+1}`$.

FIGURE 2. Fits to the data of Mehta et al. (1999) for the mean dwell times of myosin-V: (A) as a function of external load, $`F`$, at different ATP concentrations; (B) as a function of \[ATP\] under an external load $`F=0.4`$ pN and a prediction for $`F=2.3`$ pN. The solid curves represent Eq. 11 with the central parameter values in Eqs. 12 and 13; the dashed curves represent the mean dwell times predicted for a 50:50 mixture of short, $`d_{()}=30.5`$ nm, and long, $`d_{(+)}=41.5`$ nm steps using the same values for the other parameters: see the subsection Variability of Step Sizes, below. (Note that in part (B) the dashed curve for $`F=0.4`$ pN cannot be distinguished from the solid curve.)

FIGURE 3. The force-velocity or ($`F`$, $`V`$) dependence of myosin-V at various concentrations of ATP as predicted using the parameter values in Eqs. 12 and 13: solid curves. The corresponding dashed curves follow from a model with alternating long and short steps ($`d_{(+)}=41.5`$ nm and $`d_{()}=30.5`$ nm) but otherwise the same zero-load rate constants and load distribution factors, $`\theta _j^\pm `$. The superimposed data bars (for \[ATP\]$`=1\mu `$M and 2 mM) derive from the observed dwell times by using the approximate relation $`Vd/\tau `$ (with $`d=36`$ nm); they track the predictions for $`V(F)`$ fairly well because of the paucity of reverse or backward steps under loads $`F2.5`$ pN.

FIGURE 4. Predictions for the variation of the randomness, $`r`$, of myosin-V as a function of \[ATP\] at low ($`F=0.4`$ pN) and high external load ($`F=2.3`$ pN).

FIGURE 5. Predicted behavior of myosin-V under assisting (i.e., negative) and resisting (positive) external loads, $`F`$, for two ATP concentrations: (A) mean dwell time; (B) randomness. See the text for appropriate caveats.


text:
# Images of the Dark Soliton in a Depleted Condensate

## I  Introduction

Let us start with a trivial but important remark on measurement in quantum mechanics: in general a single measurement on an $`N`$-particle system is neither equivalent to many measurements on the $`N`$-particle system, nor to $`N`$ measurements on $`N`$ one-particle systems. The only exception is the case of a product $`N`$-particle state where all particles are described by the same single-particle wavefunction. Not surprisingly the particles must be bosons forming a Bose-Einstein condensate, and a measurement on this $`N`$-particle state can be replaced with $`N`$ measurements on $`N`$ one-particle systems in the same single-particle wavefunction.

The importance of this ”nonequivalence” principle can be seen in the explanation of an interference experiment Ref.interf  in an excellent paper by Javanainen and Yoo JY , and then elaborated in Refs.localization . In Ref.JY  the authors consider a Fock state

$$|\frac{N}{2},\frac{N}{2}$$
(1)
with $`N/2`$ atoms in a condensate of wavefunction $`\varphi _0(x)`$ and $`N/2`$ in a condensate with an orthogonal wavefunction $`\varphi _1(x)`$. A single particle density matrix of the state (1) predicts that on average a density measurement will give a distribution of atoms

$$\frac{N}{2}\left(|\varphi _0(x)|^2+|\varphi _1(x)|^2\right)$$
(2)
without any interference fringes between the two condensates. This single particle distribution is an average over many density measurements on the same Fock state (1). Simulation of a destructive density measurement JY  based on the full $`N`$-particle wavefunction (1), and the actual experiment interf , give distributions with an interference term

$`{\displaystyle \frac{N}{2}}`$ $`(`$ $`|\varphi _0(x)|^2+|\varphi _1(x)|^2+`$ (3)
$`+\varphi _0(x)\varphi _1^{}(x)e^{+i\phi }+\mathrm{c}.\mathrm{c}.).`$
For a given experimental realization the phase $`\phi `$ is chosen randomly from the interval $`[0,2\pi )`$. The outcome of a density measurement is the same as if before the measurement all $`N`$ atoms were not in the state (1) but in an absurd single condensate with a wavefunction

$$\frac{1}{\sqrt{2}}\left(\varphi _0+\varphi _1e^{i\phi }\right)$$
(4)
with a relative phase $`\phi `$ that is not known before we actually measure the density distribution. It turns out localization  that it is enough to measure only a small fraction of the large total $`N`$ to prepare the remaining atoms in one of the condensate wavefunctions (4), or, in other words, to establish a definite phase $`\phi `$. Destructive measurement of atomic positions, which effectively annihilates measured atoms from the trap, drives the remaining atoms into a condensate state with a randomly picked condensate wavefunction. These interesting conclusions cannot be obtained without full knowledge of the $`N`$-particle state (1): the single particle density distribution (2) is a (misleading) average over many experiments with different outcomes $`\phi `$. To predict possible outcomes of individual measurements on an $`N`$-particle state it is essential to know its $`N`$-particle wavefunction.

A perfect condensate is a state where all $`N`$ atoms are in the same single particle state described by a single particle wavefunction $`\varphi _0(\stackrel{}{x})`$,

$$\left(\widehat{a}_0^{}\right)^N|0.$$
(5)
When we include interactions between atoms but neglect depletion of atoms from the condensate wavefunction, then $`\varphi _0`$ is a solution of the celebrated Gross-Pitaevskii equation GPE . Including small quantum depletion in the framework of the number-conserving Bogoliubov theory BT  leads to pair- correlated eigenstates Z . For example, the state without any quasiparticles (the Bogoliubov vacuum) has a pair-correlated form

$$(\widehat{d}^{})^{\frac{N}{2}}|0=\left(\widehat{a}_0^{}\widehat{a}_0^{}+\underset{k=1}{\overset{\mathrm{}}{}}\lambda _k\widehat{a}_k^{}\widehat{a}_k^{}\right)^{\frac{N}{2}}|0.$$
(6)
Operators $`\widehat{a}_k`$ annihilate in a basis of functions $`\varphi _k`$ orthogonal to the condensate wavefunction $`\varphi _0`$. This general pair-correlated ansatz (6) has been known for a while, see e.g. the review by Leggett Leggett , but a general solution for the coefficients $`\lambda _k`$ and corresponding wavefunctions $`\varphi _k`$ has been found only very recently in Ref.Z  (the solution in the special case of a uniform condensate was given in Ref.Leggett ). The two-particle creation operator $`d^{}`$ has to commute with all quasiparticle annihilation operators (see Appendix A)

$$[\widehat{b}_m,d^{}]=0,\text{for all}m.$$
(7)
The $`N`$-particle state (6) is a foundation on which one can build theory of BEC entirely in the language of $`N`$-particle wavefunctions.

In our two recent papers greyI ; greyII  we studied depletion from a condensate with a dark soliton. The dark soliton state is a collectively excited condensate with a notch in the condensate wavefunction. In the Thomas-Fermi limit of strong interactions the notch behaves like a soliton with a wavefunction close to the notch proportional to

$$\varphi _0(x)\mathrm{tanh}(xX),$$
(8)
where $`X`$ is a position of the soliton measured in units of the healing length. The anomalous Bogoliubov mode with a wavefunction

$$\varphi _1(x)_X\varphi _0(x)=\mathrm{cosh}^2(xX)$$
(9)
dominates quantum depletion from the condensate within the soliton. In Ref.greyII  we found an average number of atoms $`dN`$ depleted into the non-condensate mode $`\varphi _1`$ and a single particle density distribution which close to the soliton at $`X=0`$ looks like

$`p_1(x)`$ $`=`$ $`(NdN)|\varphi _0(x)|^2+dN|\varphi _1(x)|^2.`$ (10)
For $`dN=0`$ the distribution has a hole at $`X=0`$, but with increasing quantum depletion $`dN`$ the hole is filling with atoms. This single particle result means that on average over many experiments the notch of a depleted dark soliton appears filled with atoms (it appears grey, not dark), but it gives no clue what to expect as an outcome of a single destructive density measurement. We address this issue in the following sections.

The paper is organized as follows. In Sections II,III,IV we apply the number conserving Bogoliubov theory (shortly summarized in Appendix A) to the dark soliton, and simulate density measurement on different $`N`$-particle quantum states that all may pass under the same label: “depleted dark soliton”. We find that images of the “depleted dark soliton” depend qualitatively on the actual $`N`$-particle state. To simulate the measurements we adapted the numerical algorithm of Ref.JY  which is described in Appendix B. We conclude in Section V.

## II  Stationary Dark Soliton

The dark soliton state DS  of a condensate in a 1D harmonic trap can be defined as an antisymmetric solution of the stationary Gross-Pitaevskii equation

$$\frac{1}{2}_x^2\varphi _0+\frac{1}{2}x^2\varphi _0+g|\varphi _0|^2\varphi _0=\mu \varphi _0$$
(11)
In this paper we use the dimensionless oscillator units. In the Thomas-Fermi (TF) limit of $`g1`$ the dark soliton state is

$$\varphi _0(x)F(x)\mathrm{tanh}\left(x/l_0\right),$$
(12)
with

$$F(x)=\sqrt{\frac{2\mu x^2}{2g}},$$
(13)
where the chemical potential $`\mu (3g/2)^{2/3}/2`$ and the healing length $`l_0\sqrt{2}\left(2/3g\right)^{1/3}`$.

A perfect condensate with all $`N`$ atoms in the dark soliton mode $`\varphi _0`$ is not an energy eigenstate because interactions between atoms continuously deplete atoms from the condensate wavefunction $`\varphi _0`$. The Bogoliubov theory is a systematic way to describe small quantum fluctuations around the Gross-Pitaevskii solution, and to find stationary states of a depleted condensate.

Solution of the Bogoliubov-de Gennes (BdG) equations (see Appendix A) reveals many different modes of the quasiparticle spectrum. However, only the negative energy (anomalous) mode (9) has a wavefunction localized in the soliton notch Law ; greyII ; Z . ¿From all modes the anomalous mode will contribute the most to the density of incoherent atoms in the soliton notch. In the following, similarly as in Ref. greyII , we truncate the Bogoliubov spectrum to the anomalous mode alone and look for the Bogoliubov vacuum state in the form

$$|0_b:N(\widehat{a}_0^{}\widehat{a}_0^{}+\lambda \widehat{a}_1^{}\widehat{a}_1^{})^{\frac{N}{2}}|0$$
(14)
with $`\widehat{a}_1^{}`$ creating in the mode $`\varphi _1`$. The anomalous mode solution of the BdG equations (35) with $`\omega =\frac{1}{\sqrt{2}}`$ is

$$u_1=\frac{1}{2}\left(f_++f_{}\right),v_1=\frac{1}{2}\left(f_+f_{}\right),$$
(15)
with functions

$`f_+(x)={\displaystyle \frac{\sqrt{3g}}{\sqrt{8}\mathrm{cosh}^2\left(\frac{x}{l_0}\right)}},f_{}(x)=\sqrt{{\displaystyle \frac{2}{3}}}F(x),`$ (16)
see Ref.greyII . For $`g1`$ equations (15) are dominated by $`f_+`$, $`f_+|f_+f_{}|f_{}`$, and the two functions become the same, $`u_1v_1\frac{1}{2}f_+`$. As a result in the TF limit the anomalous mode does not mix with other modes and the dominant eigenfunction is

$$\varphi _1(x)=\frac{f_+(x)}{\sqrt{f_+|f_+}},$$
(17)
as has been verified numerically in Ref.Z . To find the asymptotic behavior of $`\lambda `$ for $`g1`$ we have to be more careful and still keep $`f_{}`$ for a moment. The eigenvalue can be obtained from the general prescription presented in Appendix A and it is

$$\lambda \frac{f_+|f_+f_+|f_{}}{f_+|f_++f_+|f_{}}1\frac{2^{7/6}3^{1/3}}{g^{2/3}}.$$
(18)
Knowing $`\lambda `$ we can find an average number of depleted atoms $`dN`$ in the Bogoliubov vacuum state (14). In the particle Fock representation the unnormalized vacuum state (14) is a sum

$`|0_b:N`$ $``$ $`{\displaystyle \underset{k=0}{\overset{\frac{N}{2}}{}}}\lambda ^k{\displaystyle \frac{\sqrt{(N2k)!(2k)!}}{(N/2k)!k!}}|N2k,2k.`$ (19)
For depletion greater than unity but small with respect to the total number of atoms $`N`$, i.e. $`1dNN`$, we may use the Stirling formula, and then replace the sum by an integral. As a result we find expressions for $`dN`$ and its dispersion $`D(dN)`$

$`dN{\displaystyle \frac{1}{\mathrm{ln}\lambda }}{\displaystyle \frac{g^{2/3}}{42^{1/6}3^{1/3}}},`$ (20)
$`D(dN)=dN.`$ (21)
The single particle density distribution in the vacuum state (14,19) is

$$p_1(x)=(NdN)|\varphi _0(x)|^2+dN|\varphi _1(x)|^2.$$
(22)
The first term on the RHS is the dark soliton density profile with a hole at $`x=0`$, while the second one stands for the depleted atoms localized in the soliton notch. A good measure of depletion is a ratio of atomic density in the notch at $`x=0`$ to the density near the notch

$$\frac{dN|\varphi _1(0)|^2}{NF(0)^2}N^{1/3},$$
(23)
suggesting that for sufficiently large number of atoms, the notch will not be visible. For the parameters of the Hannover experiment Hannover , where $`g7500`$ and $`N=1.5\times 10^5`$, $`dN60`$ atoms and the ratio is $`15\%`$.

However, Eq. (22) predicts an average over many experimental realizations. For a single experiment one has to choose randomly positions of $`N`$ atoms according to the $`N`$-particle probability distribution. In general this is rather a formidable task, but, as shown by Javanainen and Yoo JY  (see Appendix B), for a state spanned on only two modes, and with a sequential algorithm that chooses a position of $`n+1`$-st atom after positions of $`n`$ atoms have already been “measured”, this problem is polynomial in $`N`$.

In general we expect that a single measurement of a depleted dark soliton will show the same atomic density distribution as if all $`N`$ atoms were initially prepared in a condensate with one of wavefunctions

$$\varphi _0(x)\varphi _1(x)qe^{i\phi }.$$
(24)
The real parameters $`(q,\phi )`$ are chosen at random and fluctuate from experiment to experiment. In the special case of $`\phi =0,\pi `$ the condensate wavefunction (24) is simply a displaced dark soliton

$$\varphi _0(x)\varphi _1(x)qF(x)\mathrm{tanh}\left(\frac{xq\sqrt{\frac{3l_0}{4}}}{l_0}\right).$$
(25)
In general, when $`\phi 0,\pi `$, the condensate wavefunction (24) gives a density distribution

$`N|\varphi _0\varphi _1qe^{i\phi }|^2=`$
$`N\left[\varphi _0\varphi _1q\mathrm{cos}\phi \right]^2+N\varphi _1^2q^2\mathrm{sin}^2\phi .`$ (26)
For small $`q`$ this distribution can be interpreted as a soliton shifted to $`q\sqrt{\frac{3l_0}{4}}\mathrm{cos}\phi `$, but with a nonzero atomic density $`N\varphi _1^2q^2\mathrm{sin}^2\phi `$ in the notch. For $`\phi 0,\pi `$ the soliton appears as a displaced grey soliton.

Adopting the scheme of Javanainen and Yoo we generated results of individual density measurements of the Bogoliubov vacuum state (14,19). A generic outcome is a displaced dark soliton shown in Fig. 1. In Fig. 2 we show how definite values of $`\mathrm{cos}(\phi )`$ and $`q`$ are established in the course of measurement of subsequent atomic positions. It is enough to measure only a small fraction of atoms to suppress fluctuations of $`\phi `$ and $`q`$, see Fig. 2.

A density measurement of the Bogoliubov vacuum state (14,19) always gives a displaced dark soliton with $`\phi =0`$ in Eq.(24). At first sight this may seem obvious because the wavefunctions $`\varphi _0`$ and $`\varphi _1`$ are real, and all amplitudes in the state (19) are also real. However, this simple argument does not take into account the possibility that the measurement of a small fraction of atoms may localize the state of remaining atoms not in a single condensate (24), but in a superposition of several condensates localization  with different parameters $`(q,\phi )`$. Indeed, as we can see in Eq.(26), density distributions for $`(q,+\phi )`$ and $`(q,\phi )`$ are the same so the density measurement cannot distinguish $`+\phi `$ from $`\phi `$. A real superposition of two condensates with opposite phases has real amplitudes, just like the state (19). We checked that a measurement of a Fock state $`|Nn,n`$, instead of the Bogoliubov vacuum state (14,19), indeed results in a real superposition of two condensates with nonzero $`+\phi `$ and $`\phi `$, which appears as a displaced grey soliton. What makes a difference is the fact the Bogoliubov vacuum (19) is a superposition over many Fock states. It is the phase coherence between different Fock states, all with real amplitudes, that enforces the observed $`\phi =0`$.

## III  Soliton Condensate without Initial Depletion

In current experiments Hannover ; DS2  the dark solitons are generated with the help of the phase imprinting PIM . As we argued in Ref. greyII , when the depletion in the condensate before the imprinting is negligible, then the soliton state right after imprinting can be idealized by the state without any atoms depleted from the solitonic condensate $`\varphi _0`$ to the anomalous mode $`\varphi _1`$,

$$|\mathrm{\Psi }(t=0)=|N,0.$$
(27)
This state is not a stationary state. To obtain time evolution of this state we consequently truncate the Bogoliubov Hamiltonian to the anomalous mode where

$`H`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\widehat{b}_1^{}\widehat{b}_1`$ (28)
with

$`\widehat{b}_1`$ $`=`$ $`{\displaystyle \frac{\widehat{a}_0^{}\widehat{a}_1\lambda \widehat{a}_0\widehat{a}_1^{}}{\sqrt{N(1\lambda ^2)}}}.`$ (29)
With this Hamiltonian the initial perfect condensate $`|\mathrm{\Psi }(t=0)`$ evolves

into a depleted state $`|\mathrm{\Psi }(t)`$. In Fig. 3 we show how a single particle density evolves in time. After about 15 ms (for the parameters of the Hannover experiment Hannover ) the hole in the single particle density distribution fills up with depleted atoms.

We simulated single realizations of destructive density measurements of the state $`|\mathrm{\Psi }(t)`$. A generic density distribution after evolution for 10 ms is shown in Fig. 4. We see that the soliton is dark and it is displaced with respect to the trap center.

Inspection of the state after 10 ms shows that $`Nk,k|\mathrm{\Psi }(10\mathrm{ms})e^{ik\phi }`$ with a tiny $`\phi 0.01`$. In principle the time evolution might generate any phase $`\phi [0,2\pi )`$, but the parameters of the experiment Hannover  are such that there is not enough time to get any more significant $`\phi `$ from the idealized initial state (27) before the linearized Bogoliubov theory breaks down.

## IV  Generic Non-Stationary Soliton States

In a realistic experiment the dark soliton is created in neither the stationary nor the zero depletion state discussed in the previous two Sections. The two states may appear “natural” to a theorist’s eye accustomed to such highly idealized examples. However, these states are not realistic from experimental point of view. The soliton is phase imprinted on a condensate with a $`10\%`$ thermal cloud i.e. in a highly excited non-stationary state. We do not know what the actual state, or ensemble of states, is, but we can still test how robust is the picture developed in the previous Sections with respect to imperfections of the quantum state.

To get an idea what are images of a generic non-stationary state we model the state by a state with random phases. The Hamiltonian of the system is diagonal in the Fock basis of quasiparticles so we expect a generic non-stationary state in to have apparently random phases in particle Fock representation. To be more specific, we take as an example the Bogoliubov vacuum (19) with randomized phases $`\alpha _k`$,

$$|\mathrm{rand}\underset{k=0}{\overset{N/2}{}}e^{i\alpha _k}\lambda ^k\frac{\sqrt{(N2k)!(2k)!}}{(N/2k)!k!}|N2k,2k.$$
(30)
Single particle density in this state is exactly the same as for the Bogoliubov vacuum, see Eq. (22), but possible outcomes of density measurements can be very different in the two states. Indeed, a generic outcome in the state $`|\mathrm{rand}`$ shown in Fig. 5 reveals a displaced grey soliton with $`\mathrm{cos}\phi =0.45`$. The phase disorder present in the state $`|\mathrm{rand}`$ translates into random $`\phi `$’s.

As a second example, we take the state $`|\mathrm{rand}`$ as an initial state instead of the state (27), evolve it for 0.7 ms, and then measure atomic positions. In Fig.6 we show a single particle density after 0.7 ms and a corresponding generic outcome of density measurement. Not surprisingly, we again see a displaced grey soliton, this time with $`\mathrm{cos}\phi =0.48`$.

These two examples demonstrate that phase disorder in particle Fock representation can change qualitatively possible results of density measurements.

## V  Conclusion

In Refs.greyI ; greyII  we have considered influence of quantum fluctuations on a density profile of a Bose-Einstein condensate excited to a solitonic state. We have shown that the soliton notch may be filled with atoms that are depleted from the condensate due to interactions between particles. The analysis has been performed by calculating a single particle density matrix. However, such a matrix provides informations about the density profile that is a result of average over many experimental realizations. To obtain predictions for a single experimental realization one has to choose randomly positions of atoms according to the multiparticle probability density.

In this paper we concentrate on the analysis of outcome of the destructive density measurement in a single experiment. Applying the number conserving Bogoliubov theory of a BEC we calculate $`N`$-particle quantum states of a depleted condensate with a dark soliton. The $`N`$-particle wavefunctions are used to see what happens after the depleted condensate is subject to an ideal destructive (density) measurement of atomic positions. We have found that:

* For a BEC excited to a stationary state, i.e. the Bogoliubov vacuum state without any quasiparticles, a destructive measurement reveals density distribution which appears like a displaced dark soliton with an empty notch. Average displacement is proportional to the average number of incoherent atoms in the notch.
* For the non-stationary state that evolves from an initial state without any incoherent atoms in the soliton notch (i.e. a perfect condensate with a dark soliton) a destructive measurement reveals, similarly like in the previous case, a displaced dark soliton. The longer the evolution lasts the larger the depletion of the condensate becomes (and thus also average displacement of the soliton).
* The previous two cases are related to highly idealized states of a BEC. To analyze influence of imperfections (in the preparation of the soliton condensate) on the result of density measurement we introduced random phases in amplitudes of different Fock states corresponding to the Bogoliubov vacuum state. In such a non-stationary state the measurement reveals a displaced grey soliton. Both the average displacement and the average density of atoms in the notch of this grey soliton are proportional to the number of atoms depleted from the condensate.

Very different $`N`$-particle quantum states may correspond to the same single particle density matrix. However, different states result in qualitatively different results of destructive density measurements. In view of the strong dependence of measurement results on a quantum state, it is an open question what $`N`$-particle state is relevant for the present day dark soliton experiments Hannover . The phase imprinting method PIM  prepares neither an energy eigenstate nor the initial state without any depleted atoms. Certainly, further experiments would be helpful to settle this problem.

Quantum fluctuations of a condensate are an intriguing result of interatomic interactions that introduce two-particle and higher correlations into the state of the system. BEC as a correlated system might be a very interesting toy model for condensed matter physics. Depletion to the notch of the dark soliton gives a unique opportunity to detect quantum fluctuations simply by density measurement. A dark soliton in a strictly one-dimensional trap is well within the reach of present technology.

## APPENDIX

### V.1 The N-particle State of BEC

Here we briefly summarize main conclusions of the number-conserving Bogoliubov theory of BEC BT ; Z . A field operator $`\widehat{\psi }(\stackrel{}{x})`$ is split into a condensate part $`\widehat{a}_0\varphi _0(\stackrel{}{x})`$ and a non-condensate part $`\delta \widehat{\psi }(\stackrel{}{x})`$,

$$\widehat{\psi }(\stackrel{}{x})=\varphi _0(\stackrel{}{x})\widehat{a}_0+\delta \widehat{\psi }(\stackrel{}{x}).$$
(31)
To zero order in the “small” fluctuation operator $`\delta \widehat{\psi }(\stackrel{}{x})`$ the $`N`$-particle state is a perfect condensate

$$\left(\widehat{a}_0^{}\right)^N|0$$
(32)
with all $`N`$ atoms in the condensate wavefunction $`\varphi _0(\stackrel{}{x})`$. The $`\varphi _0`$ solves a stationary Gross-Pitaevskii equation (GPE)

$$\mu \varphi _0=\frac{1}{2}^2\varphi _0+V(\stackrel{}{x})\varphi _0+g|\varphi _0^2|\varphi _0.$$
(33)
To zero order the interaction affects only the shape of the condensate wavefunction $`\varphi _0(\stackrel{}{x})`$ through the nonlinear term $`g|\varphi _0^2|\varphi _0`$ in the GPE.

Quantum depletion from this condensate shows up in the second order of the perturbation theory where the $`N`$-particle state becomes a pair-correlated Bogoliubov vacuum state

$$|0_b:N\left(d^{}\right)^{\frac{N}{2}}|0=\left(\widehat{a}_0^{}\widehat{a}_0^{}+\underset{k=1}{\overset{\mathrm{}}{}}\lambda _k\widehat{a}_k^{}\widehat{a}_k^{}\right)^{\frac{N}{2}}|0.$$
(34)
The sum runs over an orthonormal basis of non-condensate modes $`\varphi _k(\stackrel{}{x})`$ orthogonal to $`\varphi _0(\stackrel{}{x})`$. The eigenvalues $`\lambda _k`$ and the eigenmodes $`\varphi _k(\stackrel{}{x})`$ have been calculated in full generality only very recently Z . Their construction runs in the following steps. Its full justification can be found in Refs.BT ; Z .

* Solve the stationary GPE (33) to get $`\varphi _0(\stackrel{}{x})`$ and $`\mu `$. Note that $`\varphi _0`$ does not need to be the ground state, it can be an excited stationary state with a dark soliton or vortex.
* Solve the linear Bogoliubov-de Gennes equations for the Bogoliubov normal modes of small fluctuations around $`\varphi _0`$,

$`{\displaystyle \frac{1}{2}}^2U_m+V(\stackrel{}{x})U_m+2g|\varphi _0|^2U_m+g\varphi _0^2V_m=`$
$`\mu U_m+\omega _mU_m,`$
$`{\displaystyle \frac{1}{2}}^2V_m+V(\stackrel{}{x})V_m+2g|\varphi _0|^2V_m+g(\varphi _0^{})^2U_m=`$
$`\mu V_m\omega _mV_m`$ (35)
to get the eigenvalues $`\omega _m`$ and “raw” Bogoliubov eigenmodes $`U_m(\stackrel{}{x}),V_m(\stackrel{}{x})`$.
* Project the raw $`U_m,V_m`$ on the subspace orthogonal to $`\varphi _0`$,

$`u_m`$ $`=`$ $`U_m\varphi _0\varphi _0|U_m,`$
$`v_m`$ $`=`$ $`V_m\varphi _0\varphi _0|V_m,`$ (36)
to get the Bogoliubov eigenmodes $`(u_m,v_m)`$.
* Normalize the modes $`(u_m,v_m)`$ so that a norm

$$u_m|u_mv_m|v_m=+1.$$
(37)
Ignore modes with negative norm.

The normalized Bogoliubov modes define bosonic quasiparticle annihilation operators

$$\widehat{b}_m=\frac{\widehat{a}_0^{}}{\sqrt{N}}u_m|\delta \widehat{\psi }\frac{\widehat{a}_0}{\sqrt{N}}\left(v_m^{}|\delta \widehat{\psi }\right)^{}$$
(38)
which transfer atoms between the condensate and the non-condensate modes, but conserve the total number of atoms $`N`$. For small depletion the Hamiltonian can be approximated by a sum of harmonic oscillators

$$\widehat{H}=\underset{m}{}\omega _m\widehat{b}_m^{}\widehat{b}_m.$$
(39)
The $`N`$-particle Bogoliubov vacuum (34) is the eigenstate of $`\widehat{H}`$ without any quasiparticles,

$$\widehat{b}_m|0_b:N=0.$$
(40)
The $`N`$-particle state (34) is annihilated by all $`\widehat{b}`$’s when the $`\widehat{d}^{}`$ defined in Eq.(34) commutes with all $`\widehat{b}`$’s,

$$[\widehat{b}_m,\widehat{d}^{}]=0.$$
(41)
In general construction of such a $`\widehat{d}^{}`$ runs along the following steps.

* Choose an orthonormal basis of states $`\stackrel{~}{\varphi }_k(\stackrel{}{x})`$ orthogonal to $`\varphi _0`$, and then calculate matrices $`\stackrel{~}{U}`$ and $`\stackrel{~}{V}`$ with elements

$`\stackrel{~}{U}_{mk}=u_m|\stackrel{~}{\varphi }_k,`$ (42)
$`\stackrel{~}{V}_{mk}=\stackrel{~}{\varphi }_k|v_m.`$ (43)
* Invert $`\stackrel{~}{U}`$ and then calculate a matrix

$$\stackrel{~}{Z}=\stackrel{~}{U}^1\stackrel{~}{V}.$$
(44)
* Diagonalize $`\stackrel{~}{Z}`$ to get its eigenvalues $`\lambda _k`$ and orthonormal eigenmodes $`\varphi _k`$. In the basis of $`\varphi _k`$ the $`\stackrel{~}{Z}`$ becomes a $`Z=\mathrm{diag}\{\lambda _1,\lambda _2,\mathrm{}\}`$.

Having calculated $`\lambda _k`$ and $`\varphi _k`$ we can construct the operator $`\widehat{d}^{}`$ and the Bogoliubov vacuum state (34). Excited eigenstates can be obtained by repeated action of the quasiparticle creation operators $`\widehat{b}_m^{}`$ on the vacuum state.

### V.2  Density Measurement

In an ideal density measurement positions of all $`N`$ atoms are measured. The probability distribution for the $`N`$ positions in a $`N`$-particle quantum state $`|\mathrm{\Psi }`$ is given by

$`p_N(x_1,\mathrm{},x_N)`$
$`\mathrm{\Psi }|\widehat{\psi }^{}(x_1)\mathrm{}\widehat{\psi }^{}(x_N)\widehat{\psi }(x_N)\mathrm{}\widehat{\psi }(x_1)|\mathrm{\Psi }`$ (45)
We want to simulate density measurements by generating from this distribution typical outcomes $`(x_1,\mathrm{},x_N)`$. In general this is a formidable task.

The problem becomes tractable when, for some physical reasons, it is possible to truncate the single particle Hilbert space to just two modes, say, $`\varphi _0(x)`$ and $`\varphi _1(x)`$ (generalization to $`2,3\mathrm{}`$ modes is straightforward). With only two modes we can adopt the algorithm of Javanainen and Yoo JY . In their algorithm positions $`(x_1,\mathrm{},x_N)`$ are not generated all at once but one after another. To begin with, $`x_1`$ is chosen randomly with a reduced single particle probability distribution

$$p_1(x_1)\mathrm{\Psi }|\widehat{\psi }^{}(x_1)\widehat{\psi }(x_1)|\mathrm{\Psi }.$$
(46)
Once the actual $`x_1`$ is chosen we calculate a state vector

$$|\mathrm{\Psi }_1=\widehat{\psi }(x_1)|\mathrm{\Psi }=\left[\widehat{a}_0\varphi _0(x_1)+\widehat{a}_1\varphi _1(x_1)\right]|\mathrm{\Psi },$$
(47)
and store it in the memory. After $`(k1)`$ steps we know $`(x_1,\mathrm{},x_{k1})`$ and a state vector

$$|\mathrm{\Psi }_{k1}=\widehat{\psi }(x_{k1})|\mathrm{\Psi }_{k2},$$
(48)
calculated for the actual $`x_{k1}`$. To generate the next coordinate $`x_k`$ we use a conditional probability distribution

$`p(x_k`$ $`|`$ $`x_{k1},\mathrm{},x_1)={\displaystyle \frac{p_k(x_k,\mathrm{},x_1)}{p_{k1}(x_{k1},\mathrm{},x_1)}}`$ (49)
$`{\displaystyle \frac{\mathrm{\Psi }_{k1}|\widehat{\psi }^{}(x_k)\widehat{\psi }(x_k)|\mathrm{\Psi }_{k1}}{\mathrm{\Psi }_{k1}|\mathrm{\Psi }_{k1}}}.`$
Note that this is a function of $`x_k`$ only, all coordinates $`x_{k1},\mathrm{},x_1`$ have already been fixed. Once we get the actual outcome $`x_k`$ we calculate $`|\mathrm{\Psi }_k=\widehat{\psi }(x_k)|\mathrm{\Psi }_{k1}`$ and store it in the memory.

The key simplification discovered in Ref.JY  is that all single particle probability distributions that we encounter, like Eqs.(46,49), have the same functional form

$`p(x)`$ $`=`$ $`A|\varphi _0(x)|^2+B|\varphi _1(x)|^2+`$ (51)
$`C\varphi _0^{}(x)\varphi _1(x)+C^{}\varphi _0(x)\varphi _1^{}(x).`$
To find the real $`A,B`$ and the complex $`C`$ we need to sample $`p(x)`$ at only four values of $`x`$. When $`\varphi _0,\varphi _1`$ are real, like for the dark soliton, then $`C`$ is also real, and just $`3`$ sampling points are enough. The total algorithm to generate $`N`$ positions $`(x_1,\mathrm{},x_N)`$ scales like $`N^2`$, a slight improvement over the $`N^3`$ in Ref.JY .

To summarize, we provide an input $`N`$-particle state $`|\mathrm{\Psi }`$ (in a Fock representation build on the two single particle modes $`\varphi _0,\varphi _1`$) and the algorithm generates a random string of coordinates $`(x_1,\mathrm{},x_N)`$ according to the $`N`$-particle probability distribution (45). Each application of the algorithm returns a random outcome $`(x_1,\mathrm{},x_k)`$ of a single realization of density measurement on the state $`|\mathrm{\Psi }`$.

### Acknowledgments

We acknowledge helpful discussions with James Anglin, Yvan Castin, and Jakub Zakrzewski. This research was supported in part by KBN grants 2 P03B 092 23 (JD) and 5 P03B 088 21 (ZPK and KS).


text:
# 1 Summary of the quality of the fits for different scenarios considered in this report, for √𝑠≥5 GeV: DB02Z – the 2002 Review of Particle Properties database with new ZEUS data, DB02Z-CDF – with the CDF point removed; DB02Z-E710/E811 – with E710/E811 points removed. The first line gives the overall 𝜒²/𝑑𝑜𝑓 for the global fits, the other lines give the 𝜒²/𝑛𝑜𝑝 for data sub-samples, the last line gives in each case the parameter controlling the asymptotic form of cross sections.

## Acknowledgments

The COMPAS group was supported in part by the Russian Foundation for Basic Research grants RFBR-98-07-90381 and RFBR-01-07-90392, K.K. is in part supported by the U.S. D.o.E. Contract DE-FG-02-91ER40688-Task A, E.M. is a visiting fellow of the Belgian FNRS. We thank Professor Jean-Eudes Augustin for the hospitality at LPNHE-Université Paris 6, where part of this work was done.