Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

Kengo
Nishi
^{a},
Maria L.
Kilfoil
*^{bc},
Christoph F.
Schmidt
*^{ad} and
F. C.
MacKintosh
*^{efghi}
^{a}Third Institute of Physics-Biophysics, University of Göttingen, 37077 Göttingen, Germany
^{b}Alentic Microscience Inc., Halifax, NS B3H 0A8, Canada
^{c}Department of Physics & Atmospheric Science, Dalhousie University, Halifax, NS B3H 4R2, Canada. E-mail: maria.kilfoil@dal.ca
^{d}Department of Physics, Duke University, Durham, NC 27708, USA. E-mail: cfschmidt@phy.duke.edu
^{e}Department of Chemical & Biomolecular Engineering, Rice University, Houston, TX 77005, USA. E-mail: fcmack@gmail.com
^{f}Center for Theoretical Biological Physics, Rice University, Houston, TX 77030, USA
^{g}Department of Chemistry, Rice University, Houston, TX 77005, USA
^{h}Department Physics & Astronomy, Rice University, Houston, TX 77005, USA
^{i}Department of Physics and Astronomy, Vrije Universiteit, 1081HV Amsterdam, The Netherlands

Received
20th December 2017
, Accepted 21st March 2018

First published on 21st March 2018

Passive microrheology typically deduces shear elastic loss and storage moduli from displacement time series or mean-squared displacements (MSD) of thermally fluctuating probe particles in equilibrium materials. Common data analysis methods use either Kramers–Kronig (KK) transformation or functional fitting to calculate frequency-dependent loss and storage moduli. We propose a new analysis method for passive microrheology that avoids the limitations of both of these approaches. In this method, we determine both real and imaginary components of the complex, frequency-dependent response function χ(ω) = χ′(ω) + iχ′′(ω) as direct integral transforms of the MSD of thermal particle motion. This procedure significantly improves the high-frequency fidelity of χ(ω) relative to the use of KK transformation, which has been shown to lead to artifacts in χ′(ω). We test our method on both model and experimental data. Experiments were performed on solutions of worm-like micelles and dilute collagen solutions. While the present method agrees well with established KK-based methods at low frequencies, we demonstrate significant improvement at high frequencies using our symmetric analysis method, up to almost the fundamental Nyquist limit.

The practical implementation of passive microrheology is subject to at least three limitations: (1) the temporal bandwidth and spatial resolution limits of the method used to measure the fluctuations, (2) the artifacts introduced by the analysis of these fluctuations to derive the micromechanical probe-particle response, and (3) the accuracy and appropriateness of models such as generalizations of the Stokes formula that are employed to relate the particle response to bulk material properties. Here, we address primarily the second of these limitations by introducing an improved method of analysis to obtain the local, micromechanical compliance or response function χ of the probe particle from displacement time-series data. Our analysis presupposes that either a direct displacement time series x(t) or the mean-squared displacement (MSD) 〈(x(t′ + t) − x(t′))^{2}〉 of particle motion is measured over about three decades or more in t. Various experimental techniques can be used for such measurements. The MSD of an ensemble of probe beads can, for instance, be obtained using light scattering methods, including dynamic light scattering (DLS) and diffusing wave spectroscopy (DWS), where the intensity fluctuations of the coherently scattered light can be related to the MSD.^{15} These methods offer high spatial resolution and a wide bandwidth (∼1 nm and 10–10^{5} Hz for DWS). Moreover, the averaging extends over hundreds or more probe particles, resulting in good statistics. This averaging, however, can become problematic if particles reside in inhomogeneous micromechanical environments.

Alternatively, laser interferometry can be used to track thermal motions of single beads.^{8,9,16,17} This method delivers ∼1 nm spatial resolution and a frequency range of 0.1 Hz up to 100's of kHz, and averaging over multiple probe particles must be done sequentially. Finally, individual bead motions can be imaged directly over time using video microscopy, employing standard optical microscopes and digital cameras. The high-frequency limit is determined by the camera, and is typically less than 100 Hz, but can be extended up to 10s or even 100s of kHz with the use of specialized high-speed cameras.^{18} In this approach, 10s to 100s of beads can be tracked simultaneously, resulting in good statistics of the ensemble-averaged MSD. Importantly, both of the latter tracking methods retain single-bead trajectories for further evaluation. Among other things, this can be used to directly identify inhomogeneities in the micromechanical environments of the probe particles. From the displacement time series x(t), one calculates the power spectral density (PSD), typically by fast Fourier transformation (step (a) in Fig. 1), or the mean-squared displacement (MSD) (step (b) in Fig. 1).

The MSD and PSD are entirely equivalent, as they are simply related by Fourier transformation. Thus, either of these can be obtained (indicated by steps (a) and (b) in Fig. 1) to a level of accuracy and over a time/frequency range limited by noise and detection bandwidth. From the PSD, one can directly determine the imaginary part χ′′ of the response function via the FDT (Fig. 1 step (c)). A significant advantage of this approach is its rigor up to this point: the accuracy of the result is limited only by experimental noise and detection band-width. A significant practical limitation, however, of prior implementations of this method enters with the calculation of the real part of the response function χ′, which is needed to determine the full complex shear modulus. This has led to a loss of accuracy at high frequencies over as much as a decade below the Nyquist limit, which imposes a cut-off in the required Kramers–Kronig integral transformation.

A common alternative approach starts with the calculation of the MSD of particle motion (Fig. 1 step (b)). One can determine the complex shear modulus in an approximate manner by performing a Laplace transform on the MSD (Fig. 1(e)) and fitting with an assumed functional form. One can then transform this fitted function to obtain an estimate of the frequency-dependent shear modulus using the generalized Stokes formula. Even assuming that the generalized Stokes formula is valid, the result can only be as good as the assumed functional form for the fitting. The choice of that function is either empirical or based on an expected form of the shear modulus. It thus represents an uncontrolled approximation: functions that look very similar in the Laplace domain can represent very different functions when continued to the Fourier (frequency) domain. In practice, this indirect approach often depends on some knowledge or expectation of the rheological properties of the medium.

Our aim here is to develop a more direct and still rigorous and objective approach to the calculation of the local micromechanical compliance χ(ω) directly from the MSD (Fig. 1(d)). From this, the complex shear modulus can also be determined directly (Fig. 1(f)), provided the response can be accurately modeled in terms of the macroscopic rheology, e.g., via a Stokes-like formula. In Section 2, we review the relevant linear response theory and present our analysis method, applying it to simulated data, constructed for a viscoelastic medium. In Sections 3 and 4 we present experimental results, comparing our new method with other approaches. Finally, in Section 5 we discuss the limits of our approach and compare with alternative approaches.

(1) |

For a spherical particle of radius a moving in a Newtonian liquid of viscosity η, the natural observable A of interest is the particle velocity v(t), which responds to a mechanical force f(t) that is conjugate to the particle position x. On time scales long enough to neglect inertial effects, eqn (1) reduces to an instantaneous response v = μf, where μ = 1/(6πηa) is the particle mobility or inverse drag coefficient. Here, we assume that the Stokes formula describes the drag force on the particle. This limit corresponds to above. Similarly, in the limit of a purely elastic response with no inertia, eqn (1) reduces to an instantaneous response in the particle position x(t) ∝ f(t). For a particle in an incompressible elastic medium with shear modulus G, the Stokes formula generalizes to above. More generally, by Fourier transforming this linear response relation in eqn (1), we obtain (ω) = (ω)(ω), where

(2) |

(3) |

The fluctuation–dissipation theorem (FDT) states that the correlation function 〈x(t)x(0)〉 and response function χ are related by^{19}

(4) |

In view of the strictly even symmetry of 〈x(t)x(0)〉 as functions on the full real line, it is convenient to define the even and odd functions

(5) |

(6) |

2kT_{O}(ω) = iωC(ω). | (7) |

(8) |

2kTχ′′(ω) = ωC(ω), | (9) |

From the power spectrum C(ω) and from χ′′(ω), one can determine χ′(ω) = _{E}(ω), and therefore the full using a Kramers–Kronig transformation^{8–10}

(10) |

M(t) = 〈(x(t) − x(0))^{2}〉 = 2[〈x^{2}〉 − 〈x(t)x(0)〉]. | (11) |

(12) |

Alternatively, one can simply perform one-sided cosine and sine transformations of eqn (12) to obtain

(13) |

(14) |

The combination of eqn (13) and (14) represent a symmetric derivation of χ′ and χ′′, in contrast to the use of a Kramers–Kronig transformation. This should improve the high-frequency results for G′(ω) and G′′(ω) for viscoelastic materials. As noted, our approach assumes that Ṁ(t) vanishes at long times. Assuming this, the approach in eqn (13) and (14) provides a mathematically direct implementation of the FDT, since it involves a single time derivative. Recently, Evans et al. introduced a data analysis method for microrheology based on Fourier transformation of the second derivative of M(t) with respect to time.^{22–24} Yanagishima et al. have also demonstrated a similar method for determining the rheology by analysis of the velocity autocorrelation function.^{25} These approaches are particularly effective for systems with a long-time fluid-like response, for which neither M(t) nor its first derivative Ṁ(t) can be directly transformed. The use of the second derivative is, however, more complicated with discretely sampled experimental data with noise. The resulting analysis is necessarily less direct than the present approach based on Ṁ(t) (see ESI†). In practice, many microrheology experiments employ optical traps that result in bounded M(t) and vanishing Ṁ(t). Moreover, the main challenge we aim to address here concerns the high-frequency limits of microrheology that arise from finite sampling rates. Thus, our approach based on Fourier transformation Ṁ(t) should provide a practical approach in most cases.

(15) |

(16) |

Fig. 2 (a) Calculation of χ′(ω) (in units of 1/(2kT)) from simulated mean-squared displacement data: exact χ′(ω) for idealized semiflexible polymer network, together with χ′(ω) obtained from sampled exact expressions for PSD (KK FFT and KK Simpson) or from the numerical derivative and exact integration of M(t) (symmetric method). The corresponding Nyquist frequency is shown as the vertical dashed line. (b) The real and imaginary parts of the micromechanical stiffness K = 1/χ (in units of 2kT) are shown. Exact results are shown for comparison with three different methods, based on a Kramers Kronig integral (KK FFT), the symmetric method described here, as well as the method described by Evans et al.^{22,24} Again, vertical dashed line indicates the Nyquist frequency. |

Here, the derivative Ṁ was computed numerically using the five-point stencil method,^{27} with which the derivative can be obtained to order δ^{4} in the spacing δ between consecutive data points. This is more accurate than a simple (order δ^{2}) local slope of M(t) obtained from pairs of consecutive data points. Specifically, we use

(17) |

(18) |

(19) |

In order to compare with prior analysis methods, we used the exact χ′(ω) from eqn (16), again truncating the summation at n = 11. We also numerically sampled the exact χ′′(ω) at the Nyquist frequency (f = 100 in our units). From this, the real part of the response function χ′(ω) was then calculated using a Kramers Kronig integral. The KK integral was evaluated over a finite frequency range.^{9,21} This can be done in either of two ways, as indicated in Fig. 2a. In the first of these (KK FFT), the KK integral was calculated as the convolution of the functions 1/ω and χ′′(ω) as shown in the first line of eqn (10). To speed up this calculation, we performed additional Fourier and inverse Fourier transformations using fast Fourier transformations (FFT) in the last line of eqn (10). In the second method (KK Simpson), χ′ is obtained with an improved integration algorithm for eqn (10) using Simpson's rule. This provides, at the cost of speed, a substantial improvement over simple FFT, particularly at high frequencies. Both KK-based methods still exhibit high-frequency artifacts due to finite bandwidth. In contrast, we observe improvement in the high-frequency region when avoiding the Kramers–Kronig integral altogether, using instead the new symmetric method described above.

In Fig. 2b, we show the corresponding results for the complex stiffness K(ω) = 1/χ(ω), which is proportional to the complex shear modulus, e.g., for a response governed by the Stokes formula. Here, in addition to the exact results, we show the results using a Kramers Kronig integral (KK FFT), as well as the symmetric method described above. We find significant improvement on both K′ and K′′ using the latter method. Interestingly, we still observe a larger error in K′ than K′′ over the last approximately 1/3 decade below the Nyquist frequency. This is likely due to the greater sensitivity to small values of t in the cosine transform in eqn (13) than the sine transform in eqn (14). We have also included a comparison to the method of Evans et al.,^{22,24} based on the second derivative of the MSD, M(t). Here, too, we found a greater error in the evaluation of K′ than in that of K′′.

Collagen type I (rat tail) was purchased as a stock solution with a concentration of 11.32 mg ml^{−1} in 20 mM acetic acid (BD Biosciences, Heidelberg, Germany). For the experiments, type I collagen was diluted to the desired final concentration of 0.57 mg ml^{−1} in water. A small quantity (<0.001% by volume) of COOH surface-modified polystyrene beads with a diameter of 1 or 2 μm (Kisker Biotech GmbH & Co. KG, Steinfurt, Germany) was added to each sample before measurement.

As shown in Fig. 3, the KK FFT method becomes unreliable at high frequencies, roughly a decade below the Nyquist frequency. The present symmetric method proves to be more robust in the high frequency region, where the KK methods fail (see ESI† for algorithm). Experimental noise in the MSD is a potential issue when performing the discrete derivative to evaluate Ṁ(t). Random experimental noise, however, is strongly suppressed at short times or high frequencies because there are good statistics for the first several points of the MSD. The improvement obtained by the symmetric method is significant because high bandwidth in the measurement of viscoelastic response is one of the main motivations for microrheology. We have also applied the method of ref. 22 to analyze our sampled MSD data (see ESI†). The results are shown in Fig. 3. Both the present symmetrical method, as well as that of ref. 22, show an unphysical down-turn in G′ that begins approximately a factor of four below the Nyquist frequency, which represents an improvement over the loss of accuracy over somewhat more than a full decade when using the KK integral approach.

Microrheology data were also recorded for collagen solutions with a concentration of 0.57 mg ml^{−1}. The bead diameter was 2 μm, the sampling rate was 10 kHz, and the measuring time was 100 s. The filaments in the solution were fibrillar collagen bundles with a diameter of ∼200 nm and contour lengths of tens of μm, as described in the literature,^{31} rather than single collagen triple helices (diameter ∼ 2 nm). In such a network with inhomogeneity on large scales, the Stokes formula for the particle response is unlikely to be applicable: on the scale of the probe particles we use, we cannot assume that the material can be treated accurately as a viscoelastic continuum. Nevertheless, these reconstituted collagen solutions lack active elements, such as molecular motors or contractile cells, and they are in thermal equilibrium. In this equilibrium state, we are justified in assuming linear response theory to relate the position fluctuations of a single probe particle to the local micromechanical compliance χ(ω), even if the Stokes formula is invalid. Hence, instead of the complex shear modulus G(ω), we report our results in terms of a complex spring constant K(ω), defined as K(ω) = K′(ω) − iK′′(ω) = 1/χ(ω).

Fig. 4 shows the results for a dilute collagen solution. Although the micromechanical compliance is not necessarily directly related to the macroscopic shear modulus, its frequency dependence is similar to that of the complex shear modulus measured in semiflexible filament networks that can be approximated as viscoelastic continua. The new symmetric method again is clearly less affected by the finite bandwidth of recording than the KK FFT method, across the entire frequency range. In particular, the high-frequency regime of K′(ω) calculated by the symmetric method shows power law scaling for more than a decade in frequency with a slope of ∼0.75. This scaling corresponds to theoretical predictions for both semiflexible polymer fluctuations and the rheological response of semiflexible networks.^{26,32,33}

The method presented here is most readily implemented for viscoelastic or solid-like media, for which the MSD grows sub-linearly or is bounded in time. For media exhibiting long-time liquid-like or viscous behavior, the integrals in eqn (13) and (14) can become ill-defined. In practice, this is not a problem for optical trap-based microrheology, since trapping can be used to regularize these integrals at long times. Moreover, the main point of the method introduced here is to address problems in resolving short-time or high-frequency behavior. Other methods can be used to resolve the low-frequency response.

- J. D. Ferry, Viscoelastic Properties of Polymers, Wiley, 1980 Search PubMed .
- R. B. Bird, Dynamics of Polymeric Liquids, Wiley, 1987 Search PubMed .
- R. G. Larson, Constitutive Equations for Polymer Melts and Solutions, Butterworths, 1988 Search PubMed .
- M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Clarendon Press, 1988 Search PubMed .
- W. W. Graessley, The Entanglement Concept in Polymer Rheology, Springer-Verlag, 1974 Search PubMed .
- T. G. Mason and D. A. Weitz, Phys. Rev. Lett., 1995, 74, 1250 CrossRef CAS PubMed .
- T. Mason, H. Gang and D. Weitz, J. Opt. Soc. Am., 1997, 14, 139 CrossRef CAS .
- F. Gittes, B. Schnurr, P. D. Olmsted, F. C. MacKintosh and C. F. Schmidt, Phys. Rev. Lett., 1997, 79, 3286 CrossRef CAS .
- B. Schnurr, F. Gittes, F. C. MacKintosh and C. F. Schmidt, Macromolecules, 1997, 30, 7781 CrossRef CAS .
- F. C. MacKintosh and C. F. Schmidt, Curr. Opin. Colloid Interface Sci., 1999, 4, 300 CrossRef CAS .
- E. M. Furst and T. M. Squires, Microrheology, Oxford, 2017 Search PubMed .
- H. Callen and T. Welton, Phys. Rev., 1951, 83, 34 CrossRef .
- R. Kubo, Rep. Prog. Phys., 1966, 29, 255 CAS .
- A. Levine and T. Lubensky, Phys. Rev. Lett., 2000, 85, 1774 CrossRef CAS PubMed .
- M. Gardel, M. Valentine and D. Weitz, Microscale Diagnostic Techniques, Springer Verlag, 2005, ch. Microrheology Search PubMed .
- T. G. Mason, K. Ganesan, J. H. van Zanten, D. Wirtz and S. C. Kuo, Phys. Rev. Lett., 1997, 79, 3282 CrossRef CAS .
- F. Gittes and C. F. Schmidt, Opt. Lett., 1998, 23, 7 CrossRef CAS PubMed .
- J. Liu, M. Gardel, K. Kroy, E. Frey, B. Hoffman, J. Crocker, A. Bausch and D. Weitz, Phys. Rev. Lett., 2006, 96, 118104 CrossRef CAS PubMed .
- D. Chandler, Introduction to Modern Statistical Mechanics, Oxford University Press, 1987 Search PubMed .
- P. Chaikin and T. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press, 2000 Search PubMed .
- L. D. Landau and E. M. Lifshitz, Statistical Physics, Pergamon Press, Oxford, 1980 Search PubMed .
- R. M. L. Evans, M. Tassieri, D. Auhl and T. A. Waigh, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 012501 CrossRef CAS PubMed .
- R. M. L. Evans, Br. Soc. Rheol. Bull., 2009, 50, 76 Search PubMed .
- M. Tassieri, R. M. L. Evans, R. L. Warren, N. J. Bailey and J. Cooper, New J. Phys., 2012, 14, 115032 CrossRef .
- T. Yanagishima, D. Frenkel, J. Kotar and E. Eiser, J. Phys.: Condens. Matter, 2011, 23, 194118 CrossRef PubMed .
- F. Gittes and F. C. MacKintosh, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 58, R1241 CrossRef CAS .
- B. Fornberg, Math. Comput., 1988, 51, 699 CrossRef .
- M. Buchanan, M. Atakhorrami, J. F. Palierne, F. C. MacKintosh and C. F. Schmidt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 011504 CrossRef CAS PubMed .
- F. Gittes and C. F. Schmidt, Opt. Lett., 1998, 23, 7 CrossRef CAS PubMed .
- M. Atakhorrami, J. I. Sulkowska, K. M. Addas, G. H. Koenderink, J. X. Tang, A. J. Levine, F. C. MacKintosh and C. F. Schmidt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 061501 CrossRef CAS PubMed .
- M. Shayegan and N. R. Forde, PLoS One, 2013, 8, e70590 CAS .
- F. Amblard, A. C. Maggs, B. Yurke, A. N. Pargellis and S. Leibler, Phys. Rev. Lett., 1996, 77, 4470 CrossRef CAS PubMed .
- D. C. Morse, Macromolecules, 1998, 31, 7044 CrossRef CAS .

## Footnote |

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c7sm02499a |

This journal is © The Royal Society of Chemistry 2018 |