PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Phys Biol. Author manuscript; available in PMC Aug 1, 2012.
Published in final edited form as:
PMCID: PMC3148805
NIHMSID: NIHMS308405
“Dicty Dynamics”: Dictyostelium motility as persistent random motion
Liang Li,1 Edward C. Cox,2 and Henrik Flyvbjerg3
1Department of Physics, Princeton University, Princeton, NJ 08544, USA
2Department of Molecular Biology, Princeton University, Princeton, NJ 08544, USA
3Department of Micro- and Nanotechnology, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark
Liang Li: liangl.pu/at/gmail.com; Edward C. Cox: ecox/at/princeton.edu; Henrik Flyvbjerg: henrik.flyvbjerg/at/nanotech.dtu.dk
We model the motility of Dictyostelium cells in a systematic data-driven manner. We deduce the minimal dynamical model that reproduces the statistical features of experimental trajectories. These are trajectories of the centroid of the cell perimeter, which is more sensitive to pseudopod activity than the usual tracking by centroid or nucleus. Our data account for cell individuality and dictate a model that extends the cell-type-specific models recently derived for mammalian cells. Two generalized Langevin equations model stochastic periodic pseudopod motion parallel and orthogonal to the amoeba’s direction of motion. This motion propels the amoeba with a random periodic left-right waddle in a direction that has a long persistence time. The model fully accounts for the statistics of the experimental trajectories, including velocity power spectra and auto-correlations, non-Gaussian velocity distributions, and multiplicative noise. Thus we find neither need nor place in our data for an interpretation in terms of anomalous diffusion. The model faithfully captures cell individuality as different parameter values in the model, and serves as a basis for integrating the local mechanics of cell motion with our observed long-term behavior.
1.1. Bottom-up vs. top-down approaches to motility modeling
Most cells respond at some point in their life cycles to external cues by moving toward or away from the source of the cue. Epidermal cells respond to wounds by moving toward the wound site and dividing as they approach it. Neutrophils sense peptides released by bacteria and move singly and in groups toward the bacterial source. Axonal growth cones migrate many hundreds of cell diameters in the developing nervous system. And free living amoebas are capable of sensing food sources and migrating toward them, often over a thousand or more cell diameters. A great deal about the molecular machinery underlying these processes is known, and there are now good cellular models describing the roles played by receptor-ligand interactions, downstream signaling pathways, and changes in the cytoskeleton as cells sense then respond to graded signals in the environment (recently reviewed in [1]). How forces are generated as cells move has received considerable recent attention as new methods for quantifying forces at the sub-cellular level have been developed and exploited [2, 3, 4]. All of these experiments, which might be characterized as molecular or bottom-up approaches, attempt to understand the macroscopic behavior of single cells in terms of their microscopic properties. Despite these many advances, there remain largely unanswered questions about how cells develop polarity when chemotaxing in a gradient—whether it be by a symmetry breaking event in response to a chemoattractant, guided by an internal vector, or by the selection of more or less randomly extended pseudopods—and the features of the response systems that allow for great sensitivity to small gradients across the cell, of order 1% differences between leading and trailing edges [5].
Many cell types also move in the absence of external cues, and as they make the transition from this state to a directed one, additional cellular machinery comes into play as receptors and cytoskeletal components become mobilized at the leading and trailing edge of the chemotaxing cell. Cell motion in the absence of external signals is by comparison to the studies mentioned above, relatively unstudied, and it would seem obvious that the transition to a highly structured motile cell moving toward or away from a signal source will be built on the undirected “ground state” in which the machinery for signal detection and cell motion is not highly structured. An important part of research in this area is to accurately measure and then characterize cell motion. This is a top-down approach, and it is important because it defines the motion that the bottom-up approach attempts to explain.
1.2. Discrete vs. continuous stochastic events in modeling
Recently Takagi et al. [6] modeled the statistics of Dictyostelium trajectories using the NHDF-model of Selmeczi et al. [7, 8]. Simultaneous with Takagi et al., Li et al. explored the dynamics of isolated Dictyostelium cells by tracking the centroid of the cell perimeter [9]. Cells were separated by many cell lengths in these experiments, and thus there were no chemotactic sources, nor other clues to directions. Also, cells were followed for very much longer times, resulting in a richer data set.
Li et al. modeled their trajectories as zig-zag paths composed of straight runs and sharp turns, reminiscent of the run-and-tumble trajectories of E. coli [10], and demonstrated that the run-and-turn data could be accounted for as a simple Markov process, i.e., with no memory extending back beyond the last turn. Bosgraaf and van Haastert soon after demonstrated that this interpretation of trajectories is well founded [11]. They tracked the entire cell perimeter as pioneered in [12, 13, 14, 15], defined pseudopodia as any convex protruding segment of the perimeter, and found that pseudopodia come in two flavors, split and de novo, according to whether they grow out of an existing pseudopod or not. Series of the first kind give rise to a waddling motion with directional persistence, while the interruption of such a series by a pseudopod of the second kind gives rise to a turn. They found good experimental support for the simplest possible statistics, which decides the flavor of the next pseudopod in a Bernoulli trial [16, Fig. 3]. Except for the waddling motion of the straight “runs,” this is a two-dimensional version of Fürth’s persistent random walk in one dimension [17, 18]. In the continuum limit of infinitesimal steps, the two-dimensional version of Fürth’s persistent random walk becomes the run-and-tumble dynamics of E. coli [10], which has the run-and-turn trajectories of Li et al.
Figure 3
Figure 3
Velocity auto-correlation functions for individual cells of strain Ax2. Velocities were defined as (displacement vector)/(time-lapse) for each 5 s time-lapse of each cell’s trajectory.
The trajectory of the centroid of the perimeter is smooth, however, so it is reasonable to model it with a generalized Langevin equation. We demonstrate here how to do this. This demonstration shows in practise that both zig-zaging and smooth interpretations of trajectories are valid. They result from different emphases in data interpretation, and not from fundamental differences in motility patterns.
1.3. “Programmed” periodic motion vs. Markov processes. Diffusion vs. super-diffusion
A dichotomy between “programmed” periodic motion, on the one hand, and Markov processes on the other was discussed in [19]. Our 100% data-driven modeling turns out to combine the two descriptions in a Markov process that on an intermediate time-scale describes “programmed” periodic motion around a persistent direction of motion, while it displays ordinary diffusion on its longest time scale.
The latter ordinary behavior shows that there is neither need nor place in our data for anomalous diffusion of the kind observed in [20, 21]. Our model’s mean squared displacement as function of time has a smooth transition from deterministic to purely diffusive behavior. This transition can be misinterpreted as super-diffusive behavior, if not followed through to its true long-term behavior. We follow it through for the velocity auto-correlation function, for which it is easier, pending good statistics, which we have.
2.1. Experimental Materials and Methods
Dictyostelium discoideum AX4 and AX2 were grown on lawns of Escherichia coli B/r at 22 °C. Vegetative amoebas were harvested after 1.5 days and bacteria removed by centrifugation. The cells were suspended in PB (20 mM KH2PO4, 20 mM Na2HPO4 · 7H2O) and plated on 1% agar in distilled deionized water at a very low density so that the average distances between individual cells are more than 1000 cell diameters. Cell motions were followed by phase contrast microscopy with 15X magnification. Movies were recorded with a CoolSNAP-Pro Color camera from RS Photometrics with a frame every 10 or 5 s, respectively, for AX4 and AX2 cells for 8–10 hr—see movies in [9]. For more details about experimental methods, see [9].
An algorithm was developed in MATLAB to track amoeba. Cells were recognized in each frame by the program, and cell locations were defined as the centroid of the pixels that cover a cell’s contour; see Fig. 1.
Figure 1
Figure 1
A cell (black pixels), its contour (white pixels), and the centroid of the contour (center of magenta circle). The centroid, or center of gravity, of the white contour-pixels, is the point with which a cell is tracked. The width of a pixel corresponds (more ...)
This centroid is more sensitive to pseudopod extensions than the centroid of the cell—see Fig. 1—so its trajectory reveals more about the pseudopod activity that drives the motion. The resulting trajectory consists of a Δt = 5 or 10 s time-lapse sequence of coordinates in the plane: r.gif" border="0" alt="[r with right arrow above]" title=""/>j = r.gif" border="0" alt="[r with right arrow above]" title=""/>(tj), tj = jΔt, j = 0, 1, 2, …; see Fig. 2. The Ax2-cell trajectories studied here have 5 s intervals. The Ax4-cell trajectories used for some of our Supplemental Information have 10 s intervals and are the trajectories used in [9].
Figure 2
Figure 2
Typical trajectory, produced by time-lapse recording of amoeba every 5 s, i.e., with sampling frequency 0.2 Hz. The standard deviation σ of the Gaussian error (due to pixel round-off errors) on our determination of the centroid positions that (more ...)
Care was taken to make the environment experienced by an amoeba homogeneous, isotropic, and constant in time: The agar surface was the same everywhere, in all directions, and at all times, apart from the track left by the amoeba. Figuratively, a cell found itself in a dense fog with no compass or clock. All directions look the same at all times, except if it encounters its own tracks, which several cells do, with some tendency to follow them, see Fig. 2.
If a cell follows its own tracks for a while, that does not change the statistical properties of its trajectory, since the part of the track that is followed was laid down originally on track-free agar. Thus, a cell following its own track merely over-samples a part of a track that was laid down on virgin agar.
We checked that trajectories obey the statistics expected from a constant and isotropic environment: Local time-averages of a cell’s speed should be constant in time, apart from fluctuations due to finite statistics. The start-to-end direction of a trajectory should point to any corner of the world with equal probability, so the distribution of such directions should be even. Figures S4 and S5 demonstrate that this is the case.
2.2. Theoretical Methods
In Subsect. 2.2.1 below, visual inspection of trajectories reveals the necessary minimum of ingredients in a model for this kind of trajectories. In Subsect. 2.2.2, inspection of the velocity auto-correlation functions of individual cells excludes that they differ simply through individual time-scales. This excludes the simple modeling that such scaling invites. The long-time behavior of the velocity auto-correlation functions, however, is seen to be a simple exponential decrease. So the long-time behavior of trajectories is just persistent random motion.
In Subsect. 2.2.3, velocity power spectra of individual cells reveal the nature of the motility model that will describe these data: It is a first-order stochastic (integro-) differential equation according to Occam’s Razor, because most of the spectrum is just f−2-behavior. A broad bump in the spectrum at low frequencies signals an oscillatory component overlying the persistent random motion.
With the stage thus set, it is changed to the amoeba’s frame of reference in which trajectory data are peeled systematically for successive layers of information in Subsects. 2.2.4 through 2.2.9, each corresponding to a part of a model. In Subsect. 2.2.10 the parts are finally assembled to a model in the lab frame of reference.
2.2.1. Reading tracks
Trajectories appear to consist of persistent motion overlaid with noisy periodic motion (Figs. 2 and S1–S3). The noise is so large, however, that a trajectory is difficult to predict based on its past, even one period into its future.
The apparently random component of the motion that we refer to as noise, may be the result of the cell’s deterministic response to local randomness in the structure of the agar it crawls on. Or it may be the deterministic result of complex processes on time-scales that we do not resolve. Or both. Whatever its origin, if it has the statistical properties of noise, then the averaging done when estimating a trajectory’s velocity auto-correlation function, will reveal the deterministic component of the motion, while Fourier analysis will expose a periodic deterministic component as one or more peaks in the velocity power spectrum. We focus on these statistics based on velocities rather than the positions we actually recorded, because velocities remain bounded, hence can be averaged and Fourier transformed meaningfully, while positions are arbitrary, and displacements are unbounded functions of time. Also, a dynamical model of motility should be expressed in terms of velocities and cannot depend on a cell’s position or its orientation.
2.2.2. The velocity auto-correlation function
Figure 3 shows our experimental results for the velocity auto-correlation function defined as
equation M1
(1)
where left angle bracket · · · right angle bracket denotes averaging over time t for the individual cell. It is estimated for a given trajectory and tlag = [ell]Δt by using v.gif" border="0" alt="[v with right arrow above]" title=""/>j = (r.gif" border="0" alt="[r with right arrow above]" title=""/>jr.gif" border="0" alt="[r with right arrow above]" title=""/>j−1)/Δt as the velocity in
equation M2
(2)
The terms subtracted in numerator and denominator remove bias, which otherwise would occur for finite statistics. These subtractions cause another, small bias, which is removed by division with N[ell] − 1 instead of N[ell] in the numerator and by N − 1 instead of N in the denominator.
Each correlation function in Fig. 3 decreases essentially exponentially at time-lags larger than ~2 min. The exponential decrease is overlaid with correlated noise, which resembles an oscillatory component to the correlation, but is not: Its amplitude does not decrease with increasing tlag, but remains constant for hours, up to the longest time-lags we have reliable statistics for (not shown). This noise is due to finite statistics and only looks like a deterministic oscillation because it is strongly auto-correlated for short times and anti-correlated for slightly longer times because the velocity itself is (anti-)correlated in this manner, as we shall see below.
The slow exponential decrease in the velocity auto-correlation function is a quantitative measure of the persistence observed in “oscillation-averaged” trajectories: A trajectory showing oscillations around a somewhat persistent direction of motion will have velocities that oscillate as well, around an “oscillation-averaged” velocity that changes slowly compared to the period of oscillations. In a steady state of motion, this rate of change is constant, giving rise to an exponential decrease with a characteristic time longer than the period of oscillation. The simplest possible examples of this simple exponential decrease—without oscillations on top—are found for Fürth’s persistent random walker and the Ornstein-Uhlenbeck process [17, 22, 23].
An oscillating velocity is positively auto-correlated for time-lags less than a quarter period, negatively correlated for time-lags around half a period, then positively correlated again for times around one period, etc. If the period of such an oscillating velocity fluctuates around a mean value, then its auto-correlation function decreases to zero with increasing time-lag. This behavior is observed in Fig. 3 in each individual correlation function in the manner these functions at time-lags less than 4 min differ from the slowly decreasing exponential functions observed at time-lags larger than 4 min. This interpretation of Figs. 3 and and22 is borne out in Figs. 6 and and77 below.
Figure 6
Figure 6
(A) Auto-correlation function for v[perpendicular] = u. The difference between the theoretical function (black) and experimental data (red) is explained by finite statistics. Inset: The non-diminishing amplitude of the experimental auto-correlation function (more ...)
Figure 7
Figure 7
(A) Auto-correlation function for w. The difference between the theoretical function (black) and experimental data (blue) is explained by finite statistics. Inset: The non-diminishing amplitude of the experimental auto-correlation function at large time (more ...)
Note that the experimental correlation functions in Fig. 3 cannot be made to fall on top of each other by plotting each function against time-lag/(function’s persistence time), since such an individual rescaling of the time axis cannot bring the “knees” in the graphs at 1.5–2 min time-lags to collapse on top of each other. Dictyostelium cells display more than one internal time-scale in their trajectories. Thus Fig. 3 shows that individuality of cells amounts to more than a simple change of scales (compare [7, Fig. 7], [24]). This individuality will be described by a single model below, with different parameter values distinguishing individual cells.
2.2.3. The velocity power spectrum
The power spectrum of the velocities in a given trajectory is just the Fourier transform of the trajectory’s velocity auto-correlation function (the Wiener-Khintchin theorem), so it is not independent information. It displays short-time behavior better, however, as high-frequency behavior. And its experimental values at different frequencies are statistically independent if—or to the extent that—the velocity-process can be modeled with a linear stochastic differential equation with constant coefficients, driven by white or colored additive noise. If this statistical independence is in effect, it is much easier to fit a theory to the power spectrum than to the auto-correlation function, since the latter’s values are as correlated as the velocities themselves.
Figure 4 shows an attempted data collapse of velocity power spectra for different cells. The “bump” or “shoulder” occurring in the same place in each spectrum can be understood as a peak that has been broadened by noise and is standing on a broad-band noise spectrum, similar to the peak seen, e.g., in the power spectrum of a damped harmonic oscillator at finite temperature [25].
Figure 4
Figure 4
Partial data collapse of velocity power spectra. (A) For individual cells, with each cell’s spectral values measured in units of its time-averaged squared speed. In these units all spectra would be identical, if statistical properties of different (more ...)
At frequencies above the broadened peak, all spectra decrease like 1/f2, which, because of our finite sampling frequency fsample = 12 min−1 = 0.2 Hz, shows aliasing near the Nyquist frequency equation M3. Specifically, aliasing turns 1/f2 into [26, 27],
equation M4
(3)
which is shown as a dashed-dotted curve in Fig. 4, Panel B. It is seen to describe the averaged spectra perfectly between 1 and 2 min−1. Above 2 min−1, however, the average of experimental spectra overshoots the aliased 1/f2-behavior by an amount that is explained by a random error on the centroid’s position in each frame, and due to the finite size of pixels in the camera used; see Supporting Material for details. This error’s contribution to the power spectrum is additive and of the form [23]
equation M5
(4)
Thus, apart from this positional noise-spectrum, the spectrum is an (aliased) 1/f2 for f > 1 min−1. This signature behavior points to a first-order stochastic integro-differential equation driven by an additive white noise as the simplest dynamical model for the velocity. The model-free analysis done up to this point thus has narrowed down the search considerably. Qualitatively, the motion described by a first-order stochastic differential equation is less “dynamical,” more “contemplative” than motion described by higher-order differential equations. It does not bounce, it glides. Qualitatively, first-order dynamics is closer to amoeboid motion than higher-order dynamics.
2.2.4. The amoeba’s frame of reference
The velocity auto-correlation function and corresponding power spectrum considered above are both averages over orientation—scalars in mathematical terms. This makes them model-free and easy-to-evaluate, since no interpretation of trajectories is involved. Details are lost by this averaging, however, since an amoeba moves with persistence, because persistence implies a polar cell. A polar cell is likely to have different dynamics for its longitudinal and transverse velocity components, as it indeed does for NHDF-cells (but not for HaCaT-cells) [7]. Consequently, observations of an amoeba’s gait are best made from a frame of reference that remains aligned with the amoeba’s orientation.
There is no objective way to define such a frame, except we cannot proceed as in [7], where the velocity v.gif" border="0" alt="[v with right arrow above]" title=""/> itself was used to define the orientation of motile cells. With the fine time-resolution of our Dictyostelium data, the experimental noise on positions due to finite pixel size occasionally causes an approximately 180° change in the apparent direction of motion, followed immediately by another such ~180° change in direction. A dynamical model, such as a stochastic differential equation, predicts the immediate future from the current and past state of the system. This is done most reliably with a small random component. Consequently, it should not be done from a reference frame, which itself jumps about at random. We therefore choose the velocity of the persistent component of an amoeba’s trajectory as our definition of that amoeba’s orientation at any given time.
There is no objective way to define this persistent component. Any good approximation will do, however. We “reconstruct” the persistent component from the instantaneous velocity by taking a weighted average over past values. The rationale for this is that the cell’s body must have passed approximately “between its footprints,” so the cell’s velocity vector is approximately the local weighted average of the velocity v.gif" border="0" alt="[v with right arrow above]" title=""/> defined by its “footprints.” Let v.gif" border="0" alt="[v with right arrow above]" title=""/>(t) be the amoeba’s velocity at time t. Let γ be the inverse of a characteristic time that is larger than the period of the oscillatory component of the amoeba’s motion. Then
equation M6
(5)
is a local time-average of v.gif" border="0" alt="[v with right arrow above]" title=""/>, which has filtered away the oscillation. Define R.gif" border="0" alt="R" title=""/>(t) as the time-integral of V.gif" border="0" alt="V" title=""/>(t) (see Supporting Material for mathematical details). Figure 5 compares the trajectory of R.gif" border="0" alt="R" title=""/>(t) with the centroid’s trajectory r.gif" border="0" alt="[r with right arrow above]" title=""/>(t).
Figure 5
Figure 5
The trajectory R.gif" border="0" alt="R" title=""/>(t) (red) obtained by integrating V.gif" border="0" alt="V" title=""/>(t) resulting from γ = 0.4334 min−1. R.gif" border="0" alt="R" title=""/>(t) cuts smoothly through the oscillations of the experimental trajectory r.gif" border="0" alt="[r with right arrow above]" title=""/>(t (more ...)
2.2.5. Modeling strategy
Define two orthogonal unit vectors e.gif" border="0" alt="e" title=""/>||(t) and e.gif" border="0" alt="e" title=""/>[perpendicular](t) in the direction of V.gif" border="0" alt="V" title=""/>(t), respectively orthogonal to it,
equation M7
(6)
equation M8
(7)
They define a Cartesian coordinate system that changes direction with V.gif" border="0" alt="V" title=""/>. Define v[perpendicular](t) and v||(t) as the components of v.gif" border="0" alt="[v with right arrow above]" title=""/>(t) after e.gif" border="0" alt="e" title=""/>[perpendicular](t) and e.gif" border="0" alt="e" title=""/>||(t). Experimental results for auto-covariance and power spectrum of v[perpendicular] are shown in red in Fig. 6 for one amoeba, our Cell #2 of strain AX2.
We first formulate a dynamical theory for v[perpendicular], treating it as independent of v||. Then we formulate a theory for v|| by writing it as a sum of two terms: A deterministic part, defined as the expected value of v|| for a given dynamical state of the variable v[perpendicular], plus a stochastic part w that is not correlated with v[perpendicular], because v||’s dependence on v[perpendicular] is captured in v||’s deterministic part. The dynamics of (v[perpendicular], v||) will describe the dynamics of v.gif" border="0" alt="[v with right arrow above]" title=""/> in the reference frame defined by V.gif" border="0" alt="V" title=""/>, while Eq. 5 gives the dynamics of V.gif" border="0" alt="V" title=""/> in terms of the dynamics of v.gif" border="0" alt="[v with right arrow above]" title=""/>. The solution to these coupled equations will give v.gif" border="0" alt="[v with right arrow above]" title=""/> in the lab-frame.
Since our decomposition of experimental data for v.gif" border="0" alt="[v with right arrow above]" title=""/> into components v[perpendicular] and v|| will depend on the value of γ, our modeling will yield a theory with parameter values that depend on γ. Once found, we fit this theory’s parameters, including γ, to data to achieve an objective parametrization.
2.2.6. Model for v[perpendicular]’s dynamics
Let u denote v[perpendicular] to simplify notation. Our experimental estimate for the auto-covariance [var phi](t) = left angle bracketu(t)u(0)right angle bracket is shown in red in Fig. 6. It closely resembles a harmonic function with exponentially decreasing amplitude,
equation M9
(8)
which is shown in black. Here P is the persistence time and ω0 the cyclic frequency of [var phi](t). Discrepancies between data and theoretical function are fully and consistently explained as correlated noise on the experimental result: Statistically identical correlated noise is observed also at all larger values of t for which we have equally good statistics. The theoretical expression shown in black is not a fit to the experimental data, because the correlated noise on the latter makes correct fitting very complicated. Instead, its Fourier transform, the power spectrum
equation M10
(9)
was fitted to the experimental spectrum shown in Fig. 6B, after we had aliased the theoretical spectrum and added the spectrum of positional noise, Eq. 4, to it. With Eq. 9 we have introduced the parametrization σ2 = 2[var phi](c)/P and equation M11, which acquires meaning below.
Note that at the Nyquist frequency in Fig. 6B, the ratio of signal to positional noise is just below one, so more frequent recordings than those we have done here, would be pointless. The sampling frequency we used was optimal, because it reveals the spectrum’s noise-component so well that its signal-component is also determined well.
Note also in Fig. 6B that the experimental power-spectral values appear to scatter about the fitted spectrum in an uncorrelated manner. They do according to our theory, and this is the reason why we fit to the power spectrum: it makes correct fitting easy, and the quality of the fit can be directly inspected by eye. The distribution of the scatter is known theoretically, which simplified correct fitting [28]. Using the Wiener-Khintchin theorem, we find that the values for the parameters (P, ω0, a, σ), which are returned by the fit of the power spectrum, give the parameter values of the auto-covariance as
equation M12
(10)
The sharp peak at 0.4 min−1 in the power spectrum in Fig. 6B corresponds to a period of 2.5 min for the transverse velocity u. This period can also be observed directly in the movie of the moving cell as the period of left-right alternating pseudopod extensions: the average duration of one left-right extending pair of pseudopods is also 2.5 min [29]. These pseudopods are ~3–5-fold slower than those reported in [11].
We restrict our search for a dynamical theory to the few that give rise to the auto-covariance in Eq. 8. The simplest possible is linear, of the form
equation M13
(11)
with σ the constant amplitude of an additive white noise; see Supporting Material for mathematical details. The next-simplest such theory has the same form, except σ depends on u. All other theories that have the same auto-covariance, must be constructed and fine-tuned to ensure that their more complex features cancel with each other to produce the simple auto-covariance in Eq. 8. So Occam’s Razor points unambiguously to Eq. 11, with or without u-dependent σ.
This integro-differential equation is so simple in form that it can be rewritten as two coupled, ordinary differential equations, du
equation M14
(12)
equation M15
(13)
if one introduces the variable equation M16. Thus, if σ is a constant or σ = σ(u(t), U(t)), the state of the centroid’s motion transversely to the amoeboid’s motion is fully specified at any time t by (u(t), U(t)). These two values determine all future values of u and U for given noise ξ, according to Eqs. 12 and 13. This dynamics results in the power spectrum Eq. 9 with parameters
equation M17
(14)
2.2.7. Model for deterministic part of v||
Extension of a pseudopod will typically cause simultaneous motion of the centroid in the transverse and longitudinal directions. Thus, v|| and v[perpendicular] are correlated by nature. We model their correlation as a functional dependence of v|| on v[perpendicular] because we assume cells are non-chiral—i.e., don’t know left from right. Consequently, v|| can be a function of v[perpendicular], an even function. But v[perpendicular] cannot be a single-valued function of v||.
Since the state of transverse motion at any time t is fully given by (u(t), U(t)), we write
equation M18
(15)
where the longitudinal fluctuation w by definition satisfies
equation M19
(16)
everywhere in the (u, U)-plane. By tiling that plane with rectangles of suitable size and estimating left angle bracketv||right angle bracket within each rectangle from the experimental time series for v||(t), we obtain an experimental result for left angle bracketwright angle bracket|(u,U). We fit that result with a low-degree polynomial in u and U and find that a second-degree polynomial is sufficient. Because cells are non-chiral, that polynomial should remain invariant when the sign of u, hence U, is changed. Thus it is of the form
equation M20
(17)
Figure S7 shows a fit of this four-parameter function to 24 different experimental function values, which are seen to be well described by this polynomial. Figure S8 illustrates how much of v||(t) is explained as a deterministic function of (u(t), U(t)), hence how much is the fluctuation w.
2.2.8. Model for stochastic part of v||
Figure 7A shows the experimental auto-covariance of w with a plot of [var phi](t) in Eq. 8 superimposed. The parameter values in [var phi](t) were taken from the fit of the theoretical power spectrum in Eq. 9 to the experimental one in Panel B, after aliasing of Eq. 9 and adding the contribution from positional noise, Eq. 4.
The stochastic dynamics of w is thus seen to be described by a pair of equations identical to those in Eq. 12 describing u’s dynamics, apart from the parameter values being different. w describes a damped periodic pseudopod extension with a large random component, just as u does. They model the same pseudopod at the same time, but w models only that part of its longitudinal motion, which is uncorrelated with its transverse motion.
2.2.9. Model for noise
The model developed above is maximally simple for σ constant with independent values for u and w, but it is not unique. Velocity-dependent amplitudes of the form σu(u(t), U (t)) and σw(w(t), W (t)) result in the same velocity auto-correlation functions and power spectra, with σ2 replaced by equation M21 and equation M22. A difference is seen in the distribution of u-values and w-values, however. These distributions are Gaussians for σ constant, while experimental data show that this is not the case (Fig. S9 and top panel in Figs. S10–S12).
We determine σu(u, U) and σw(w, W) experimentally as follows: Equation 12 states that the expectation value left angle bracketdu/dtright angle bracket|u,U of du/dt for given (u, U) is a linear combination of u and U, and that its fluctuations around this expectation value is delta-function auto-correlated in time with amplitude σu(u(t), U (t)). We consequently measured left angle bracketdu/dtright angle bracket|u,U experimentally by calculating it from our recorded time series, and then calculated the standard deviation of du/dt as a function of (u, U). Figure S13 shows that left angle bracketdu/dtright angle bracket|u,U indeed is essentially linear in (u, U), and left angle bracketdw/dtright angle bracket|w,W in (w, W). Figure 8 shows the standard deviations of du/dt and dw/dt as a functions of (u, U), respectively (w, W), and that (2,2) Padé approximants model these results well. With these functions describing σu(u(t), U (t)) and σw(w(t), W (t)), the dynamics of u and w are correct by construction using Itō calculus. Figure 9 illustrates this: The perfect agreement seen here between theoretical and experimental velocity distributions is not due to fitting to the data shown, but the result of tailoring the theory, as just described, to other aspects of the time-series that is the source of the data shown here. Figures S10–S12, bottom panels, show similar plots for Ax2:Cell 3, 11, and 14.
Figure 8
Figure 8
State-dependence of noise term amplitudes. Colored: experimental results. Black: theory, fitted simultaneously to these and other data. (A, B) Standard deviation σu(u, U) of noise-term in dynamics of u = v[perpendicular], Eq. 12, as function of state (more ...)
Figure 9
Figure 9
Speed distributions. Experimental distributions: for u (red data) and w (blue data). Three black lines plotted on top of each experimental data set: model’s expectation values ± 1 SD. The 1 SD shown is what the model predicts for results (more ...)
2.2.10. “Dicty Dynamics”
The dynamical model that was derived above, describes the velocity v.gif" border="0" alt="[v with right arrow above]" title=""/> as the combined result of two dynamics: The transverse, “waddling” persistent dynamics of v[perpendicular] = u and the longitudinal, correlated dynamics of v|| = C(u, U) + w. They combine to give
equation M23
(18)
at any time. Here e.gif" border="0" alt="e" title=""/>[perpendicular] and e.gif" border="0" alt="e" title=""/>|| change with V.gif" border="0" alt="V" title=""/> according to their definitions in Eqs. 6 and 7, and the dynamics of V.gif" border="0" alt="V" title=""/> follows from Eq. 5,
equation M24
(19)
This is our model for “Dicty Dynamics.”
Figure 10 shows a simulated trajectory of this dynamics. Figures S14–S19 show more such trajectories. They were simulated with parameter values found from the experimental trajectory of AX2:Cell 2; see tables in Supporting Material. These parameter values were determined by simulating the model and taking and processing data from the simulated trajectory in exactly the same manner as was done experimentally, then fitting the model’s parameters, including γ, to achieve the best possible simultaneous agreement between theoretical and experimental data for the power spectra for u and w (Figs. S20 and S21, bottom panels), for the polynomial C(u, U) (Eq. 17 and Fig. S7), for the rational approximates to σu(u, U) and σw(w, W) (Fig. 8), and, in order also to fit long-time behavior well, the velocity-autocorrelation function (Fig. S22, top panel). The simulated trajectory was chosen so long that the various statistics just mentioned, are essentially noise-free.
Figure 10
Figure 10
Simulated trajectory produced by “Dicty dynamics” with parameter values obtained from fit of model to statistics of Ax2:Cell 2; see tables in Supporting Material.
Since power-spectral values are not perfectly Gaussian distributed, but Gamma-distributed, least-squares fitting could not be used; it would result in systematic errors. Maximum-likelihood fitting to the relevant Gamma-distribution was necessary. Computationally, least-squares algorithms are a good deal faster than maximum-likelihood algorithms, so we used the former to obtain the result of the latter with the help of an exact correction factor [28].
The true stochastic errors on our experimental results are known from our fitted theory and were computed by simulating the fitted theory for as long as it took to record the part of the AX2:Cell 2 trajectory that we analyzed, 12.5 hours. This simulation was repeated a large number of times and standard deviations on results were calculated from this large number of in silico-results.
See bottom panels of Figs. S20 and S21 for agreement between experiment and simulated model regarding u and w. Figure S22 shows agreement between experiment and simulated model regarding auto-covariance and power spectrum of v.gif" border="0" alt="[v with right arrow above]" title=""/> and probability density function of v.
3.1. Dynamics
As demonstrated in Figs. 69 and S19–S22, we have a dynamical model that produces trajectories with the same statistical properties as our experimental trajectories. This model is about as good as our data, so the full information-content of the data is captured by the model. The model is an extension of the Ornstein-Uhlenbeck model to the simplest model with oscillatory and persistent random motion that agrees with the data. The logical minimal extension of the Ornstein-Uhlenbeck model, which is studied in [30] with elegant results, cannot produce oscillatory motion. In [9], oscillations in the cell’s direction of motion is designed into the dynamics of the angles θ(t) and y(t), but the additional degrees of freedom required to model the cell’s motility are not introduced. In [31], our oscillatory dynamics for u and w was suggested for the components of v.gif" border="0" alt="[v with right arrow above]" title=""/> itself. The result is a minimal extension of the Ornstein-Uhlenbeck model to oscillatory random motion. Its persistence lasts only for a fraction of an oscillation, however. So that model’s “oscillation-averaged” trajectory R.gif" border="0" alt="R" title=""/>(t) has no persistence, but is purely diffusive Brownian motion at time-scales longer than the period of oscillations. Experimental velocity auto-correlation functions were not shown in [31], and the experimental and theoretical mean squared displacement that is shown, is shown for durations little longer than one period ([31, Fig. 7B]). Consequently, one is not confronted with this mismatch between the theory suggested there and the persistence of experimental data.
Our model is a minimal model that reproduces the data. It may well be the unique minimal model for our kind of data. Our data type is not unique, however. So better data types might be worth considering. A richer data-set could consist of the centroid r.gif" border="0" alt="[r with right arrow above]" title=""/>(t) of the perimeter supplemented with the centroid R.gif" border="0" alt="R" title=""/>(t) of the cell’s image. With such data, V.gif" border="0" alt="V" title=""/> = dR.gif" border="0" alt="R" title=""/>/dt might provide an experimentally defined reference frame from which we could study the dynamics of r.gif" border="0" alt="[r with right arrow above]" title=""/>(t). That would save us the parameter γ and its fit. On the other hand, when R.gif" border="0" alt="R" title=""/>(t) is an experimental quantity, its dynamics is not as simple as it was here, where it was just a first-order filtered version of r.gif" border="0" alt="[r with right arrow above]" title=""/>(t). The dynamics of the centroid R.gif" border="0" alt="R" title=""/>(t) of Dictyostelium cells was already studied by Takagi et al. [6], who demonstrated that it takes the NHDF-model of Selmeczi et al. to capture it. So only if the dynamics of the centroid of the perimeter, r.gif" border="0" alt="[r with right arrow above]" title=""/>(t), simplifies substantially when observed from the reference frame of the centroid, R.gif" border="0" alt="R" title=""/>(t), will this richer data-set simplify modeling.
Our approach to modeling is general in its data-driven systematics. In its present application it exposed
  • the necessity of modeling in a coordinate system that follows the cell;
  • the impossibility of letting the instantaneous velocity v.gif" border="0" alt="[v with right arrow above]" title=""/> define that coordinate system;
  • the simplest possible definition of a “cell velocity” V.gif" border="0" alt="V" title=""/>, which defines the needed coordinate system;
  • the proper interpretation of the velocity V.gif" border="0" alt="V" title=""/> that popped out of the mathematical problem defined by the motility data analyzed in [7], as well as the need for the extra term added in [31, Eq. (3)].
    Other workers have smoothed cell trajectories in various ad hoc ways in order to obtain a low-pass filtered “instantaneous” v.gif" border="0" alt="[v with right arrow above]" title=""/> that can be used to define the needed coordinate system following the cell. Such smoothing, however, throws the child out with the bath-water, if one is interested in more than low-frequency behavior. We have seen here
  • that positional noise accounts for much of the “roughness” of the trajectory that results from recording the trajectory with high time-resolution;
  • that this positional noise, being uncorrelated with itself in time, is easily modeled with just a single fitting parameter, its amplitude;
  • that this modeling allows us to exploit all other information contained in the trajectory at high frequency, such as
    • the velocity auto-correlation functions in Figs. 6 and and77 at short time-separations (positional noise averages to zero, except in the first two values, which are affected in a known manner [23]);
    • our description of the motility model’s noise term in Fig. 8;
    • the speed distributions in Fig. 9, which are unfiltered except from the low-pass filtering inherent in velocities obtained from time-lapse recorded trajectories [23].
Our Dictyostelium model is also a simple protocol for how to characterize trajectories statistically in an optimal manner. One fits the model to trajectory data and the values returned for model parameters are the characterizing statistics. Despite the model’s many parameters, not one is redundant, and more are hardly possible: Every parameter was forced upon us by our insistence on an exhaustive analysis and modeling of trajectory data, but they form a hierarchy according to the relative errors with which they are determined. Those with the largest errors (C3, SUU, and SWW in Tables S3, S5, and S6 in Supporting Material) suggest that all information has been exhausted: If additional parameters are introduced to capture yet finer structures in the data, those parameters will remain undetermined due to the finite nature of the data set.
3.2. Cell biology
It is not the main purpose of this communication to extrapolate the dynamics described here to the biology of cell motion. A bottom-up approach that has our phenomenological model as target for biological explanations, seems much safer for a system as complex as an amoeba. However, the short- and long-term dynamics we observe, are very suggestive and consistent with the known behavior of pseudopods in amoeboid cells, especially their lifetime and dependence on an active cytoskeleton with the properties of an excitable system [32, 11, 33, 34]. What is not understood, is the biological nature of the internal vector modeled here as V.gif" border="0" alt="V" title=""/>, which enables an amoeba to keep its direction of motion for much longer than the lifetime of a single pseudopod, as most clearly revealed by the experiments of Samadani et al. [35].
We conclude that we have a dynamical model that teaches us how to extract the entire information content of Dictyostelium trajectories produced by tracking the centroid of the cell’s perimeter. The features of our model are generic in nature, so the model will probably also describe the trajectories of other cell types with amoeboid motion. For Dictyostelium, we note that individuality in cells is evident in Fig. 3 and in the fitted parameter values (Tables, Supporting Material). The challenge now is to relate these parameters to biology. This is optimally done in experiments that pursue how parameter values depend on other differences observed in cells, and how they respond to changes in stimuli, environment, and within cells. The model is optimally suited, we believe, for capturing even small changes quantitatively when they affect trajectories. This we base on its ability to capture individuality of cells via their trajectories and on a comparison with two alternatives: The run-and-turn description, which is much simpler, but does force an a priori interpretation on trajectories, and the description in terms of contours, which is more detailed with fine resolution of pseudopod activity, but either limited to short time series—time-series of contours shown as in [36, Fig. 9 and and10]10] and [37] simply must be short—or must be compressed by further processing, as done elegantly in [11].
Two qualitative, conceptual results were already produced by our insistence on quantitative modeling of every detail in the data: (i) We resolved more than one time scale in Dictyostelium dynamics, the longest being the persistence time for the general direction of motion that is left, after the faster time scale of the quasi-periodic left-right waddle has been filtered away. (ii) We found that this longest time scale is so long that many trajectories actually are too short to self-average properly, which explains why some workers thought they saw super-diffusive behavior in Dictyostelium trajectories [20, 21]. Such behavior requires a truly exceptional memory of spatial orientation from the cell, as any reasonable mechanism of memory looses memory at a constant rate, resulting in exponentially decreasing correlations, and hence normal diffusive behavior on time-scales much longer than the correlation time. This result of our top-down phenomenological, quantitative approach—that Dictyostelium motility can be understood without mechanisms of memory that one cannot imagine a physical/chemical/biological basis for—simplifies the target of the biochemical bottom-up approach immensely.
Thus our dynamical model’s precise capture of macroscopic Dictyostelium motility provides a natural benchmark for more detailed understandings of motility. Conversely, input from more detailed understandings must be the future of macroscopic motility modeling, now that the limiting factor no longer is lack of experimental data. Only modeling with a solid biophysical/biochemical basis in the vein of [2, 3, 38, 39, 40, 41] can filter the information-deluge that is possible now that time-resolution and duration of measurements are limited only by what is meaningful for a given cell type, and hundreds or thousands of data can be recorded for a single cell at every time-lapse.
Supplementary Material
supplement
Acknowledgments
We thank Ned Wingreen for a critical of reading the manuscript. The research leading to these results has received funding from the European Union’s Seventh Framework Programme (FP7/2007-2013) under grant agreement No. 214566-2 (Nanoscale). It was supported at Princeton by National Institutes of Health grants GM078591 and GM071508.
1. Swaney KF, Huang CH, Devreotes PN. Eukaryotic chemotaxis: a network of signaling pathways controls motility, directional sensing, and polarity. Annu Rev Biophys. 2010;39:265–289. [PubMed]
2. del Alamo JC, Meili R, Alonso-Latorre B, Rodriguez-Rodriguez J, Aliseda A, Firtel RA, Lasheras JC. Spatio-temporal analysis of eukaryotic cell motility by improved force cytometry. Proc Natl Acad Sci (USA) 2007;104:13343–13348. [PubMed]
3. Meili R, Alonso-Latorre B, del Alamo JC, Firtel RA, Lasheras JC. Myosin ii is essential for the spatiotemporal organization of traction forces during cell motility. Mol Biol Cell. 2010;21(3):405–417. [PMC free article] [PubMed]
4. Iwadate Y, Yumura S. Actin-based propulsive forces and myosin-ii-based contractile forces in migrating dictyostelium cells. J Cell Sci. 2008;121(8):1314–24. [PubMed]
5. Hu Bo, Chen Wen, Rappel Wouter-Jan, Levine Herbert. Physical limits on cellular sensing of spatial gradients. Phys Rev Lett. 2010;105(4):048104. [PMC free article] [PubMed]
6. Takagi H, Sato MJ, Yanagida T, Ueda M. Functional analysis of spontaneous cell movement under different physiological conditions. PLoS ONE. 2008;3(7):e2648. [PMC free article] [PubMed]
7. Selmeczi D, Mosler S, Hagedorn PH, Larsen NB, Flyvbjerg H. Cell motility as persistent random motion: Theories from experiments. Biophys J. 2005;89:912–931. [PubMed]
8. Selmeczi D, Li L, Pedersen LII, Nørrelykke SF, Hagedorn PH, Mosler S, Larsen NB, Cox EC, Flyvbjerg H. Cell motility as random motion: A review. Eur Phys J Special Topics. 2008;157:1–15.
9. Li L, Nørrelykke SF, Cox EC. Persistent cell motion in the absence of external signals: A search strategy for eukaryotic cells. PLoS ONE. 2008;3(5):e2093. [PMC free article] [PubMed]
10. Berg HC, Brown DA. Chemotaxis in Escherichia coli analysed by three-dimensional tracking. Nature. 1972;239:500–504. [PubMed]
11. Bosgraaf L, van Haastert PJM. The ordered extension of pseudopodia by amoeboid cells in the absence of external cues. PLoS ONE. 2009;4(4):e5253. [PMC free article] [PubMed]
12. Soll DR. DMS, a computer-assisted system for quantitating motility and dynamics of cytoplasmic flow and pseudopod formation: its application to dictyostelium chemotaxis. Cell Motil Cytoskeleton. 1988;10(Suppl):91–106. [PubMed]
13. Soll DR, Voss E, Varnum-Finney B, Wessels D. The dynamic morphology system: a method for quantitating changes in shape, pseudopod formation and motion in normal and mutant amebae of Dictyostelium discoideum. J Cell Biochem. 1988;37:177–192. [PubMed]
14. Soll DR. The use of computers in understanding how animal cells crawl. Int Rev Cytol. 1995;163:43–104. [PubMed]
15. Soll DR, Wessels D, editors. Motion Analysis of Living Cells. Wiley-Liss; New York: 1998.
16. van Haastert PJM, Bosgraaf L. Food searching strategy of amoeboid cells by starvation induced run length extension. PLoS ONE. 2009;4(8):e6814. [PMC free article] [PubMed]
17. Fürth R. Die Brownsche Bewegung bei Bercksichtigung einer Persistenz der Bewegungsrichtung. Mit Anvendungen auf die Bewegung lebender Infusorien. Z Physik. 1920;2:244–256.
18. Gail MH, Boone CW. The locomotion of mouse fibroblasts in tissue culture. Biophys J. 1970;10:980–993. [PubMed]
19. Hartman RS, Lau K, Chou W, Coates TD. The Fundamental Motor of the Human Neutrophil is Not Random—Evidence for Local Non-Markov Movement in Neutrophils. Biophys J. 1994;67:2535–2545. [PubMed]
20. Upadhyaya A, Rieu JP, Glazier JA, Sawada Y. Anomalous diffusion and non-gaussian velocity distribution of hydra cells in cellular aggregates. Physica A. 2001;293:549–558.
21. Dieterich P, Klages R, Preuss R, Schwab A. Anomalous dynamics of cell migration. Proc Natl Acad Sci (USA) 2008;105:459–463. [PubMed]
22. Uhlenbeck GE, Ornstein LS. On the theory of Brownian Motion. Phys Rev. 1930;36:823–841.
23. Li L, Gradinaru C, Cox EC, Flyvbjerg H. Time-lapse recorded trajectories of motile micro-organisms: How to connect discrete data with continuous models. 2011. To be submitted.
24. Flyvbjerg H, Jobs E, Leibler S. Kinetics of self-assembling microtubules: An ‘”inverse problem” in biochemistry. Proc Natl Acad Sci USA) 1996;93:5975–5979. [PubMed]
25. Nørrelykke SF, Flyvbjerg H. Harmonic oscillator in heat bath: Exact simulation of time-lapse-recorded data, exact analytical benchmark statistics. Phys Rev E. 2011;83:041103. [PubMed]
26. Berg-Sørensen K, Flyvbjerg H. Power spectrum analysis for optical tweezers. Rev Sci Instr. 2004;75:594–612.
27. Press WH, Flannery BP, Teukolsky SA, Vetterling WT. Numerical Recipes. Cambridge University Press; 1986.
28. Nørrelykke SF, Flyvbjerg H. Power spectrum analysis with least-squares fitting: Built-in bias and its elimination, with application to optical tweezers and AFM cantilevers. Rev Sci Instrum. 2010;81:075103. [PubMed]
29. van Haastert PJM. Private communication. 2010.
30. Peruani F, Morelli LG. Self-propelled particles with fluctuating speed and direction of motion in two dimensions. Phys Rev Lett. 2007;99:010602. [PubMed]
31. Shenderov AD, Sheetz MP. Inversely correlated cycles in speed and turning in an Ameba: An oscillatory model of cell locomotion. Biophys J. 1997;72:2382–2389. [PubMed]
32. Postma M, Roelofs J, Goedhart JJ, Loovers HM, Visser AJWG, van Haastert PJM. Sensitization of dictyostelium chemotaxis by phosphoinositide-3-kinase-mediated self-organizing signalling patches. J Cell Sci. 2004;117(14):2925–2935. [PubMed]
33. Andrew N, Insall RH. Chemotaxis in shallow gradients is mediated independently of ptdins 3-kinase by biased choices between random protrusions. Nature Cell Biol. 2007;9(2):193–200. [PubMed]
34. Xiong Y, Huang C, Iglesias P, Devreotes P. Cells navigate with a local-excitation, global-inhibition-biased excitable network. Proc Natl Acad Sci USA. 2010;107(40):17079–17086. [PubMed]
35. Samadani A, Mettetal J, van Oudenaarden A. Cellular asymmetry and individuality in directional sensing. Proc Natl Acad Sci USA. 2006;103(31):11549–11554. [PubMed]
36. Cox D, Condeelis J, Wessels D, Soll D, Kern H, Knecht DA. Targeted disruption of the abp120 gene leads to cells with altered motility. J Cell Biol. 1992;116:943–955. [PMC free article] [PubMed]
37. Titus MA, Wessels D, Spudich JA, Soll D. The Unconventional Myosin Encoded by the myoA Gene Plays a Role in Dictyostelium Motility. Mol Biol Cell. 1993;4:233–246. [PMC free article] [PubMed]
38. Ridley AJ, et al. Cell migration: Integrating signals from front to back. Science. 2003 December;302:1704–1709. [PubMed]
39. Keren K, Pincus Z, Allen GM, Barnhart EL, Marriott G, Mogilner A, Theriot JA. Mechanism of shape determination in motile cells. Nature. 2008;453:475–480. [PMC free article] [PubMed]
40. Keren K, Yam PT, Kinkhabwala A, Mogilner A, Theriot JA. Intracellular fluid flow in rapidly moving cells. Nature Cell Biology. 2009;11:1219–1224. [PMC free article] [PubMed]
41. Weiger MC, Ahmed S, Welf ES, Haugh JM. Directional persistence of cell migration coincides with stability of asymmetric intracellular signaling. Biophys J. 2010;98 [PubMed]