Quantifying Cytoskeleton Dynamics Using Differential Dynamic Microscopy

Here, we show examples of the analysis done with PyDDM from two different sets of experiments. In one set of experiments, sub-micron tracer beads were embedded in networks consisting of the intermediate filament protein vimentin and imaged using a 100x objective lens in brightfield mode at 100 frames/s (Figure 3A). Vimentin is expressed in mesenchymal cells and is a key determinant of the mechanical properties of the cytoplasm⁶⁵ and the mechanical stability of the nucleus in cells performing confined migration^66,67. So far, reconstituted vimentin networks have been studied primarily by macroscopic rheology^64,68,69, whereas the dynamics have received comparatively little attention^13,70,71. Additional details of these experiments can be found in Supplementary File 2. In the other set of experiments, active cytoskeleton networks were prepared with actin, microtubules, and myosin. Spectrally distinct fluorescent labels allowed the actin and microtubule filaments to be imaged with a two-color laser-scanning confocal microscope using a 60x objective lens at 2.78 frames/s (Figure 3B,C). Actin and microtubule filaments are both important drivers of dynamic cell shape changes, with their actions coordinated by mechanical and biochemical interactions⁷². Additional details of these experiments can be found in¹¹. Individual frames from image sequences taken in these experiments are shown in Figure 3.

Figure 3: Images from the time series analyzed. (A) Brightfield image of 0.6 µm beads in a vimentin network. (B,C) Image of the (B) microtubules and (C) actin in an active actin-microtubule composite taken with a 60x objective on a laser-scanning confocal microscope, using 561 nm excitation light for the microtubule imaging and 488 nm excitation light for the actin imaging. Please click here to view a larger version of this figure.

For images of tracer beads in vimentin networks, movies of 5000 frames with a size of 512 x 512 pixels at 100 frames/s were recorded. From these, the DDM matrix was computed at 60 logarithmically spaced lag times between 1 and 1000 frames, or 0.01 s and 10 s. To estimate the background, B, the mean of the squared Fourier-transformed images, , was computed and set equal to ^55,73. An assumption was made that, over the largest 10% of q-values, this quantity equals B/2 and that B is independent of q. This is the package's default method for estimating B, but other methods are possible by setting the background_method parameter to a different value.

With the parameters A(q) and B determined from , one can extract the intermediate scattering function (ISF) from the DDM matrix. Example ISFs are shown in Figure 4. In Figure 4A, the ISF from images of 0.6 µm diameter beads embedded in a network with a vimentin concentration of 19 µM is shown. In Figure 4B, the ISF for the same type of beads in a network with a vimentin concentration of 34 µM is shown. Interestingly, in neither case did the ISF decay to zero. At large lag times, the ISF should approach zero for ergodic systems. That is, in such systems, density fluctuations should completely decorrelate over large lag times. The fact that the ISF here did not decay to zero could have resulted from inaccurate estimates of A(q) and B, which were used to find the ISF from the computed DDM matrix. Notably, the method employed here can overestimate B in certain scenarios⁶². However, it is more likely that the dynamics of the tracer beads are truly nonergodic as the beads have a comparable size to the network mesh size and may, therefore, become caged. Other data corroborated the finding of nonergodicity. Namely, the bead size, 0.6 µm, was larger than the calculated average value for the mesh sizes of 0.4 µm for the 19 µM concentration and 0.3 µm for the 34 µM concentration. Additionally, the results from single particle tracking of these tracer beads, which are shown later, also showed confined motion.

Figure 4: Intermediate scattering functions at several wavenumbers for vimentin networks. The ISF is plotted as a function of lag time for q values from about 1 to 9 µm^-1. (A) The ISF from images of 0.6 µm beads in a vimentin network with vimentin concentration of 19 µM. (B) The ISF from images of 0.6 µm beads in a vimentin network with vimentin concentration of 34 µM. The long lag time plateau of the ISF at a value well above zero indicates nonergodicity. Please click here to view a larger version of this figure.

Given that the dynamics are likely nonergodic, the ISFs are fit to the form , where C is the nonergodicity factor ³². This form of the ISF has been used in previous studies of non-ergodic dynamics, such as that of colloidal gels^32,74 or tracer particles in actin-microtubule networks ¹⁰. The dotted black lines in Figure 4 show the fits along with the data. From these fits, one can now look at the q-dependence of the decay time, τ, and of the nonergodicity parameter, C.

Figure 5: Decay time vs. wavenumber for vimentin networks. From fits to the ISF, the decay time τ is determined for a range of q values. For clarity, we are not showing the value of τ for every q, but just a logarithmically spaced set. In blue (tan) is the data from images of 0.6 µm beads within vimentin networks with a vimentin concentration of 19 µM (34 µM). The error bars represent the standard deviations in τ across multiple movies (four movies for the data with the 19 µM network [blue] and five movies for the data with the 34 µM network [tan]). Red dash-dotted lines mark estimated bounds for our temporal and spatial resolution, as described in the results. The solid black line shows  scaling, which would indicate diffusive motion. Neither data set follows this scaling. Rather, beads in the 19 µM network show subdiffusive motion ( with β > 2), and beads in the 34 µM network show confined or caged motion. Please click here to view a larger version of this figure.

The decay times showed a large amount of uncertainty, both at the low q and high q extremes, as seen in Figure 5. The error bars on this plot show the standard deviation among four videos analyzed for the lower vimentin concentration case or five videos analyzed for the higher concentration. To understand the source of the large uncertainty at these extremes, consider both the temporal and spatial resolution. Approximate limits of the resolution are shown with three red dash-dotted lines. The two horizontal lines correspond to the minimum and maximum lag times probed. Given the frame rate of 100 frames/s and the maximum lag time corresponding to 1000 frames (20% of the total video duration), accuracy was lost when measuring dynamics occurring faster than 0.01 s or slower than 10 s. At the lower q-values, the fitted values for τ were greater than 10 s. Therefore, large uncertainties should be expected in decay times that are larger than the maximum lag time. At the higher end of the q-range, the decay time approached the minimum lag time of 0.01 s but remained above it. Rather than being limited by the temporal resolution, at these higher q values, the spatial resolution may be the limiting factor. Given the pixel size of 0.13 µm, the largest value for q was about 24 µm^-1. However, the diffraction-limited resolution does not necessarily allow accurate measurements of the dynamics at these high spatial frequencies. Approximating the optical resolution as  leads to an upper wavenumber limit of about 16 µm^-1, given the objective lens's numerical aperture, NA, of 1.4 and wavelength of light, . This is demarked by the vertical red dash-dotted line in Figure 5. Indeed, the data were noisy at large values of q. Even before this approximate upper limit of q, increased uncertainty in τ was seen, and this could be from overestimating q_max. Poorer optical resolution than predicted may be because an oil immersion lens was used to image beyond the coverslip into an aqueous sample or because the condenser lens was imperfectly aligned.

For the 0.6 µm beads embedded in the less concentrated network (19 µM vimentin), it can be observed from the log-log plot of the decay time vs. wavenumber that the decay time decreased with wavenumber in a way consistent with a power law (Figure 5). However, it does not seem to follow what would be expected for normal diffusive motion, where . Rather, τ decreased more steeply with increasing q. This is indicative of subdiffusive motion, which often occurs for beads in crowded environments such as these. Fitting τ(q) over the range of 1.4 µm^-1 to 12.3 µm^-1 to a power law of the form τ = 1/Kq^β yields the transport parameters K = 0.0953 μm^β / s and β = 2.2. For those more accustomed to thinking about normal diffusion vs. subdiffusion in terms of the mean squared displacement (MSD) of tracer particles as a function of lag time (i.e., MSD = K' Δt^α), it is helpful to recognize that the subdiffusive scaling exponent in the MSD equation, α, is equivalent to α = 2 / β. In other words, the value of β = 2.2 is consistent with a subdiffusive scaling exponent in the MSD equation of α = 0.9. One would set PyDDM to fit τ(q) over this range of q-values by specifying the indices of the array of q with either the parameter Good_q_range in the YAML file or by passing the optional argument forced_qs to the function generate_fit_report. The range of q from 1.4 µm^-1 to 12.3 µm^-1 would, for the data here, correspond to indices of the array of q from 15 to 130.

For the 0.6 µm beads in the more concentrated network (34 µM), the decay time showed little dependence on q. This is likely due to the nonergodicity of beads in a network with a smaller mesh size. To probe the nonergodicity in this system, the nonergodicity parameter, C, should be plotted as a function of q, as in Figure 6. For the 0.6 µm beads in the 19 µM vimentin network, C ≈ 0.2 with little dependence on q (not shown). However, for the network with 34 µM vimentin and for a network with an even higher concentration of 49 µM vimentin, the log of C was proportional to q² as shown in Figure 6. This relationship between C and q is expected for confined motion. For beads trapped within pockets of the network, the MSD is expected to plateau at long enough lag times (i.e., , where is the MSD and δ² is the maximum MSD). Since the ISF can be expressed in terms of the MSD as , and since the nonergodic ISF goes to C at long lag times (i.e., ), the relationship  is obtained^32,75. Therefore, one can use C(q) to find δ², and this yielded δ² = 0.017 μm² and 0.0032 μm² for the 34 and 49 µM vimentin networks, respectively (corresponding to δ = 0.13 μm and 0.057 μm).

Figure 6: Nonergodicity parameter vs. wavenumber for vimentin networks. From fits to the ISF, the nonergodicity parameter C is determined for a range of q values. In tan (red) is the data from images of 0.6 µm beads within vimentin networks with a vimentin concentration of 34 µM (49 µM). The error bars represent the standard deviations in τ across multiple movies (five movies for the data with the 34 µM network [tan] and four movies for the data with the 49 µM network [red]). The y-axis has logarithmic scaling. One observes a q-dependence of C that follows , which allows for extracting the maximum mean squared displacement, δ². Fits to  are shown with the solid lines. Please click here to view a larger version of this figure.

One can use other methods to extract the confinement size δ from the data as well as the subdiffusive exponent found from examining τ(q) for beads within the 19 µM vimentin network. Firstly, one can use the method described by Bayles et al.⁷⁶ and Edera et al.⁷⁷ to extract the MSD from the DDM matrix. Notably, this method requires no fitting of the DDM matrix. One only needs to compute the DDM matrix, D(q, Δt), and  (from which A(q) and B can be determined). Then, to find the MSD, one uses the relationship . Note that this method to find the MSD assumes that the distribution of particle displacements is Gaussian, though previous work has shown that, in certain cases, MSDs derived from DDM do agree with MSDs from particle tracking, even when the displacements are non-Gaussian⁷³. For this system, as expected⁷⁸, there is non-Gaussianity in the distribution of large displacements, as seen in Figure S1. In the PyDDM package, the function extract_MSD should be executed, which returns . Secondly, one can use single particle tracking to find the MSD. Though DDM can be used to analyze images where either the high density of particles or the limited optical resolution prohibits accurate particle localization, for the images of 0.6 µm beads in vimentin networks, we were able to localize and track beads using the trackpy software (https://github.com/soft-matter/trackpy)⁷⁹. This particle tracking software package uses the algorithms described by Crocker and Grier⁸⁰.

Figure 7: Mean squared displacement vs. lag time for vimentin networks. The MSD was determined using two methods. First, the MSD was computed from the DDM matrix (shown with solid symbols). Next, the MSD was determined by using single-particle tracking (SPT) to find particle trajectories (open symbols). Error bars are determined in the same way as described in the previous two figure legends. (A) MSDs for 0.6 µm beads in the 19 µM vimentin network indicate subdiffusive motion, with good agreement between the two methods of finding the MSD. (B) MSDs for 0.6 µm beads in the 49 µM vimentin network indicate caged motion, with good agreement between the two methods of finding the MSD and with the maximum MSD found from the nonergodicity parameter. Please click here to view a larger version of this figure.

The MSDs vs. lag time for 0.6 µm beads in the 19 µM vimentin network and in the 49 µM vimentin network are shown in Figure 7. In both cases, the MSD determined from DDM agreed well with the MSD found through single-particle tracking (SPT). Furthermore, for the less concentrated network, the subdiffusive scaling exponent (α in ) was about 0.9. This is consistent with the τ(q) scaling of  found by fitting the ISF to determine τ(q) (that is, 2/2.2 = 0.9). For the more concentrated network, the MSD plateaus at longer lag times. The maximum MSD found by analyzing the q-dependence of the nonergodicity parameter (shown in Figure 7B with the horizontal line at δ² = 0.0032 μm²) was approximately the same value that the MSDs from both SPT and DDM seemed to be plateauing toward. There is a discrepancy between the longest lag time MSDs determined from DDM and SPT in Figure 7A. While this may be due to a limited number of long lag time trajectories, it may also be the case that further optimizing the range of q values for which the DDM matrix is used to estimate  for each lag time (as done by Bayles et al.⁷⁶ and Edera et al.⁷⁷) would improve our results, and such optimization will be the focus of future work.

These experiments where image sequences were recorded of tracer beads embedded in a network of vimentin intermediate filaments allowed for independent analyses: DDM (using the package described here) and SPT (using trackpy). Both analyses can reveal the degree of subdiffusion and confinement length, allowing one to use two independent image analysis techniques to provide complementary metrics. There are additional quantities one can compare from SPT and DDM. For example, heterogeneity in the dynamics of the sample can reveal itself as non-Gaussianity in the distribution of particle displacements (i.e., the van Hove distribution) determined from SPT, as well as in an ISF determined from DDM that fits to a stretched exponential^34,35. Figure S1 shows the van Hove distribution for the 0.6 µm particles in vimentin networks and discusses the stretching exponent found from fitting the ISFs—metrics used in tandem in previous studies to demonstrate the heterogeneous dynamics of particles within biomimetic systems^9,10,47 or other crowded environments ³⁴. As another example, the ISF can be calculated from particle trajectories measured with SPT and compared with the DDM-acquired ISFs. While the mean squared displacements and displacement distributions are the metrics most often pulled from SPT analysis, one can also compute the ISF from particle trajectories, , using  (see Figure S2). This ISF can be compared with DDM-generated ISFs and used to reveal dynamics not apparent in the MSD⁵⁹.

While acquiring images of tracer particles within a network may allow one to use the complementary analysis methods of SPT and DDM, it is important to note that an advantage of DDM over SPT is that it does not require images of beads (or other features) that can be easily localized and tracked. To demonstrate this point, we next highlight the analysis of active networks of actin and microtubule filaments, where fluorescent labeling of actin and tubulin allows for the imaging of both filament types, distinguished from each other via different fluorophores, with a multi-color laser-scanning confocal microscope.

Images were acquired with a laser-scanning confocal microscope of actin-microtubule networks with activity driven by myosin (rabbit skeletal muscle myosin II; Cytoskeleton #MY02). Details of the experiments and results have been previously described¹¹, and the representative results shown here are from the analysis of two movies provided in the supplemental materials (movies S1 and S4) for¹¹. Both image sequences were recorded at 2.78 frames/s for 1000 frames.

To analyze these images, the DDM matrix was calculated for 50 lag times ranging from 0.4 s to 252 s (1 frame to 700 frames). The DDM matrix was then fit to the model , with the intermediate scattering function being . There are, therefore, four fitting parameters: A, τ, s, and B. The results of these fits are shown in Figure 8. It was observed that the DDM matrix for a particular q-value had a plateau at low lag times, increased with lag time, and then plateaued (or showed signs of beginning to plateau) at large lag times. The DDM matrix for the lower values of q did not reach a plateau at long lag times. One should, therefore, expect poor accuracy in the measurement of the decay time for these low q (large length scale) dynamics.

The characteristic decay times, τ, from the fits to the DDM matrix are shown in Figure 9. Results are presented for an active actin-microtubule composite network (similar to movie S1¹¹) and for an active actin network (similar to movie S4¹¹). Both networks were prepared with the same concentrations of actin and myosin, but the actin-only network was created without tubulin, as described in¹¹. For these two types of active networks, the observed power law relationship was . This scaling indicates ballistic motion and that the myosin-driven contraction and flow dominate over the thermal motion of the filaments. From τ = (vq)^-1, a characteristic velocity, v, of about 10 nm/s for the active actin-microtubule network and 75 nm/s for the active actin network could be found. These values are consistent with the particle image velocimetry analysis of the same videos shown in¹¹. The  scaling did not hold at the lower q values for the active actin-microtubule composite network. This is likely because the true decay times for this actin-microtubule composite network at the lower q values are longer than the maximum lag time of the computed DDM matrix. The maximum lag time is indicated with the horizontal red line in Figure 9, and the decay times deviated from the expected  scaling near these longer times.

Figure 8: DDM matrix vs. lag time for an active actin-microtubule composite network. The DDM matrix for several values of q is plotted as a function of lag time from a movie of a composite network composed of 2.9 µM actin monomers, 2.9 µM tubulin dimers, and 0.24 µM myosin. These data show the analysis of just the microtubule channel of a multicolor time series of images. Please click here to view a larger version of this figure.

Figure 9: Decay time vs. wavenumber for active actin-microtubule networks. From fitting the DDM matrix, the decay time, τ, as a function of wavenumber, q, is found. Plotted is τ vs q for images of an active actin-microtubule network (analyzing just the microtubule channel) in brown and for images of an active actin network in green. Both networks have the same concentrations of actin and myosin (2.9 µM and 0.24 µM, respectively); the actin-microtubule composite has 2.9 µM of tubulin dimers. The decay times for the active actin network are much smaller than the decay times for the active actin-microtubule network, which indicates faster motion of the active actin network. In both cases, the dynamics are ballistic as the data follows a  trend. Inset: the plot of the ISFs vs. the lag time scaled by the wavenumber (Δt × q) shows a collapse of the ISFs over a range of q values. This also indicates ballistic motion. The ISFs shown in this inset are from the active actin network. Please click here to view a larger version of this figure.

For this data of active networks, we chose to fit the DDM matrix, . This contrasts with what was done for the data of beads in the vimentin network, where A (q) and B were estimated without any fitting to isolate the ISF, f (q, Δt). In this case, for the active network data, A and B were left as fitting parameters because the methods used to estimate B did not result in good fits. The default method to estimate B is to compute  and to assume that, at large q, this goes to B/2. However, this method overestimated B for this data, which was seen in the fact that, when calculating the ISFs from B estimated in this way (not shown), the ISFs were greater than 1 at early lag times (whereas they should go from a maximum of 1 to either zero or some nonergodicity parameter with increasing lag time). One can select other methods for estimating B using the parameter background_method. One of these other methods is to estimate B to be the minimum of the DDM matrix at early lag times (set with background_method=1). A similar method was used by Bayles et al.⁷⁶, though they did not assume B was constant with q. Another option is to estimate B to be the average value over all lag times of the DDM matrix at the maximum q (set with background_method=2). These different methods for estimating the background, as well as the results for allowing B to be a free fitting parameter, are shown in Figure 10. From those plots, one can see that the amplitude, A, did not reach zero at the largest q values probed, since  did not plateau at large q (Figure 10B), and since D(q_max, Δt) went from a lower lag time plateau to some higher lag time plateau (i.e., at q_max, there was a non-zero A; Figure 10D). Therefore, neither estimating B as nor as  would be appropriate. One should inspect vs. q and D(q_max, Δt) vs. Δt before deciding on how (or if) to estimate B.

Figure 10: Background vs. wavenumber for active actin-microtubule networks. From fitting the DDM matrix, one can find the background, B, as a function of wavenumber, q. Shown is B vs. q for images of an active actin-microtubule network (analyzing just the microtubule channel) determined from these fits with the purple symbols. The three solid lines in (A) show estimates of the background found without any fitting. The top, darkest line in (A) shows the estimated background using , which may be appropriate if plateaus to a constant value at large q. From (B), note that has yet to reach a constant value at the largest q probed. Therefore, using this method overestimates the background. The bottom line in (A) shows the estimated background using . If the DDM matrix shows a low lag time plateau as shown in (C) with the red line, then this method may be appropriate for estimating the background. The middle, lightest line in (A) shows the estimated background from . This method may be appropriate if, at q_max, the amplitude, A, has reached zero. From (D), it is seen that the amplitude is non-zero and, therefore, this method overestimates the background. Please click here to view a larger version of this figure.

Supplementary Figure S1: Probability distributions of particle displacements. Probability distributions of particle displacements show non-Gaussianity for vimentin concentrations of 34 µM and 49 µM. Single-particle tracking of 0.6 µm diameter beads was performed in vimentin networks of different concentrations. Different lag times are shown in the displacement distributions for the three conditions. (A) The distribution of particle displacements in a 19 µM vimentin network are fit with a Gaussian function. The width of the Gaussian increases with increasing lag time. (B) The distribution of particle displacements in a 34 µM vimentin network shows more non-Gaussianity, especially at large displacements, than for the 19 µM case. (C) The distribution of particle displacements in a 49 µM vimentin network also shows non-Gaussianity. Further, the widths of the distributions do not increase with lag time as significantly as in the samples with lower vimentin concentrations, indicating confined motion. Non-Gaussian van Hove distributions (seen for all vimentin samples but most apparent in the higher concentrations) are associated with heterogeneous dynamics as often seen in the transport of particles in crowded and confined environments. Another indicator of heterogeneous transport that is determined from DDM analysis is the stretching exponent used to fit the intermediate scattering function (the parameter s in the equation for the ISF used here:  + ). The average stretching exponents over the q-range of 0.4 µm^-1 to 9.4 µm^-1 are, from highest vimentin concentration to lowest, 0.53 ± 0.07, 0.64 ± 0.02, and 0.86 ± 0.04 (mean ± standard deviation). Please click here to download this File.

Supplementary Figure S2: The intermediate scattering functions from DDM and SPT. The intermediate scattering functions (ISF) for five different wavenumbers are shown. The ISF versus lag time found through DDM is plotted with circular markers, and the ISF computed from single-particle trajectories with open squares. Dotted black lines show the fits to the DDM-acquired ISFs. The ISF is calculated from single-particle trajectories, , using . In (A), the ISF is shown for 0.6 µm particles in the 19 µM vimentin networks. In (B), the ISF is shown for 0.6 µm particles in the 34 µM vimentin networks. The discrepancies in the ISF found from DDM and SPT are likely due to a limited number of long lag time trajectories. Please click here to download this File.

Supplementary File 1: Protocol for using DDM. The input and output of the steps shown in the protocol are presented. Please click here to download this File.

Supplementary File 2: Details of sample preparation and example parameter files for vimentin networks. Detailed steps for sample preparation and image acquisition on vimentin networks are provided. Additionally, an example parameters file for the analysis of data presented in the representative results section on vimentin networks is also provided. Please click here to download this File.