Multi-oscillation microrheology via acoustic force spectroscopy enables frequency-dependent measurements on endothelial cells at high-throughput

Active microrheology is one of the main methods to determine the mechanical properties of cells and tissue, and the modelling of these viscoelastic properties is under heavy debate with many competing approaches. Most experimental methods of active microrheology such as optical tweezers or atomic force microscopy based approaches rely on single cell measurements, and thus suffer from a low throughput. Here, we present a novel method for frequency-dependent microrheology on cells using acoustic forces which allows multiplexed measurements of several cells in parallel. Acoustic force spectroscopy (AFS) is used to generate multi-oscillatory forces in the range of pN–nN on particles attached to primary human umbilical vein endothelial cells (HUVEC) cultivated inside a microfluidic chip. While the AFS was introduced as a single-molecule technique to measure mechanochemical properties of biomolecules, we exploit the AFS to measure the dynamic viscoelastic properties of cells exposed to different conditions, such as flow shear stresses or drug injections. By controlling the force and measuring the position of the particle, the complex shear modulus G*(ω) can be measured continuously over several hours. The resulting power-law shear moduli are consistent with fractional viscoelastic models. In our experiments we confirm a decrease in shear modulus after perturbing the actin cytoskeleton via cytochalasin B. This effect was reversible after washing out the drug. Additionally, we include critical information for the usage of the new method AFS as a measurement tool showing its capabilities and limitations and we find that for performing viscoelastic measurements with the AFS, a thorough calibration and careful data analysis is crucial, for which we provide protocols and guidelines.

with the particle radius , compressibility and , density and of the particle and surrounding medium, respectively. The time average is denoted by ⟨.⟩ with the acoustic pressure and the acoustic velocity . The dimensionless parameter̃ is given bỹ = with the viscous penetration depth =

Lateral translation of free beads
A lateral translation of free beads can be observed when applying an acoustic pressure using the AFS. This can be seen in Supp.

Force distribution of different fields of view of the same chip
The force values differ inside a field of view depending on the position which is shown in the spatial calibration map (see Fig.  3D). Additionally, the values even further differ between different fields of view in the same chip as shown in Supp.

Data analysis of the SFC
In order to obtain a correct calibration map some beads were filtered out using our software Kitsune (available on request); the filter criteria are presented in the following. The raw data obtained with the modified LabVIEW software is an ASCII-file containing the -position, the applied amplitude value % and the frequency at a time point and a force type value with = 0 at time points when no voltage signals were sent to the piezo and = 7 for time points when a pulling force calibration and = 5 for a multi-oscillation force. For the SFC analysis the time points with = 7 were analyzed. The first time point with ≠ 0 was set to = 0 s. The -data was shifted by an offset such that ( = 0) = 0 m. Firstly, the following criteria had to be fulfilled during the applied constant force before being further analyzed.
• There are no -data points that equal zero which typically indicate tracking errors.
• There are no sudden changes in , i.e. between each sampling point the difference must be lesser than 10 m.
• The maximal -data must differ by at least 5 m above the surface to filter out attached beads or weak forces that could not lift the bead as high during the force application time of 1 s.
The relevant -positional data is from the ground to the -node, i.e. the -position at which the force equals zero. However, the acoustic forces inside the AFS is typically also exerted laterally. Therefore, the bead is also displaced laterally. The force is typically still applied even when the bead has reached the -node, but due to a possible slight tilt of the flow chamber or a position-dependent -node the -position of the bead may gradually change. This is taken into account by only analyzing the start-point up to the time point when the bead first reaches the -node. For the estimation when the -node is reached the following steps were performed.
• The full uncropped range is used and Eq. 16 is used to obtain the optimal fit parameters (see description in the Results section).
• The obtained velocity is numerically differentiated to obtain the acceleration.
• The first -node estimation is the -position at the time point when the absolute of acceleration is minimal.
• If the first -node estimation is < 13 m it is set to 13 m.
• New optimal fit parameters are now obtained with thedata from the start point to the first -node estimation.
• If the fit quality, here, represented by the sum of the squared difference of the numerically obtained num and the measured meas , is greater than 1 m 2 the -node estimation is decreased by 0.1 m until the fit quality is smaller than 1 m 2 or the -node estimation is smaller than 13 m.
With the obtained fit parameters the acoustic force profile can be calculated with Eq. 9. The calibration values are saved to be run through the following filters to create the final calibration map. The values obtained from that amplitude for that bead are filtered out if at least one of the following criteria of Filter-Type 1 is fulfilled: • The fit quality (represented by the sum of the squared difference of the numerically obtained num and the measured meas ) is greater than 2 m 2 .
• The fit quality is smaller than 0.001 m 2 .
• The -position or the -position of the bead is < 0 m.
• The force at 1 m changes more than 15 % of the maximal force.
• The force at 1 m is lesser than 0 pN.
• The force at 1 m and the maximal force is greater than 80 pN. Due to the high amount of beads that were not simultaneously tracked they might overlap. Therefore, the mean values of the conversion factor of the beads that are closer than a merge distance of 10 m are used. The map is then created by creating a rectangular grid with mesh size of 0.5 m starting from the minimal measured -position to the maximal measured -position. The grid interpolation is based on a biharmonic spline and is performed by the in-built MATLAB function griddata.

Temperature-dependency of the resonance frequency
A short version of the SFC was performed inside a field of view with more than 14 beads at different locations. Only one amplitude was used, therefore the force is represented instead of the conversion factor. The measured data were in Δ = 1 • C temperature steps and Δ = 0.01 MHz frequency steps. The mean force of more than 14 beads at different positions in one field of view is shown in Supp.

Resonance frequency shift inside different medium solutions
At a fixed frequency and temperature, the force depends on the medium as shown in For every condition the conversion factors are randomly split into three subgroups for statistical analysis. The significance stars * * * represents ≤ 0.001 using a two-sample -test with three random subgroups per condition after determining a statistically significant difference between the respective condition groups using a one-way ANOVA with < 0.001.

Temperature change during force application
It is expected that the temperature changes upon force application. The in-built temperature controller inside the chip holder cannot measure the temperature inside the fluid chamber, but measures the chip's exterior instead. Therefore, an external temperature sensor (disassembled from an AFS G1 chip holder) is used to measure the temperature of the bottom glass underneath the fluid chamber. For this, thermal grease is spread onto the bottom glass. However, care has to be taken of that the thermal grease or the sensor may leave scratches at the bottom of the chip which renders the affected FoVs useless. In addition, if thermal grease is used, the interface glass to air is then changed to glassgrease which typically increases the acoustic energy transmission and thus results in a weaker force. This can lead to a lower temperature increase compared to the situation without the thermal grease.
During the SFC the temperature hardly changes (Supp. Fig.  S.7A †). The exact time of the force application during the measurement of 90 s was not known, however, no relevant temperature change could be measured.
When applying higher amplitudes, such as during the mOsc experiment, the temperature drastically increases (Supp.

HUVEC culture inside an AFS chip
To demonstrate the possibility of culturing a HUVEC monolayer inside the AFS chip, the cell development is recorded about three minutes after seeding with a sampling time of 1 min. For the recording, the chip is placed on the AFS equipment with the temperature controller instead of the incubator. Due to the closed flow system, the CO 2 environment is not of importance. Figure  S.8 shows selected time frames and Movie 2 † shows the culture inside the AFS every minute. Inside the AFS chip HUVECs started to spread out onto the surface after about 2 h and were fully spread out a day after seeding. In this case, after 93 h a monolayer of HUVEC had formed. During the incubation, the medium has been renewed with a flow rate of 1.66 l min . The medium renewal was essential to culture HUVEC monolayers inside the AFS (see Supp.

Microrheological Experiment Theory of microrheology
The mechanical properties of complex materials can be measured using microrheological experiments. To interpret the obtained data the microrheological analysis undergoes a few assumptions as follows. The first assumption is that the cell is behaving as a linear viscoelastic material. This can be considered valid for a sufficiently small deformation of the material. For a linear material the superposition principle holds, thus the relation of the strain ( ) and the applied stress ( ) has the form of the hereditary integral with the creep compliance ( ). The Laplace-transformation of Eq. S.11 † yields̄ with the complex frequency in the Laplace-domain, or in the Fourier-domain identify = , hencẽ Considering an oscillatory applied stress ( ) = 0 of frequency = 2 the resulting steady state strain is also oscillatory ( ) = 0 ( − ) with the same frequency , but with an additional phase shift . For a purely linear elastic material represented by a single spring element, the phase shift equals zero ( = 0 • ) and the strain immediately follows the applied stress; for a purely linear viscous material, a single dash-pot element, the strain lags behind the stress by = 90 • . Here, the material of interest is viscoelastic and the phase shift is therefore ∈ (0 • , 90 • ).
The viscoelastic complex modulus * can be defined as the ratio of the stress to strain * ( ) = ( ) ( ) = 0 0 . (S.14) The relation between the viscoelastic complex modulus and the creep compliance with Eq. S.13 is then The second assumption is the choice of the model that describes the viscoelastic material. In this case, a model based on fractional calculus is chosen. The following calculations are mostly based on 5 . The underlying assumed element is that the creep compliance increases as The fractional integration and differentiation can be obtained, shown by performing fractional integration followed by ordinary integer differentiation, see Schiessel et al. 5 , and Weyl's fractional integral can be notated as −∞ ≡ d − d − . Re-substituting = + 1 yields the rheological constitutive equation (RCE) and with = 2 and comparing with Eq. S.14 the previously described phase shift is Thus, for a purely linear elastic material the exponent equals zero, for a purely viscous material, = 1, and for a viscoelastic material ∈ (0, 1). Due to the scaling invariance of Eq. S.20 † the time constant is set to 0 = 1 s and 0 ≡ ′ 0 0 . Note that the RCE and the resulting complex shear modulus for the SFE (Eqs. S.19 †, S.20 †) represent a power-law behavior as often 6,7 whereas the real part represents the elastic storage modulus and the imaginary part the viscous loss modulus. It is noted that 0 can be regarded as an apparent shear modulus for the time scale of 0 = 1 s. Similar to known mechanical models, such as the Maxwell or the Kelvin-Voigt model, in which springs and dash-pots are arranged in series or parallel, the same arrangements can be applied using single fractional elements instead of simple elements. These arrangements guarantee that the equations provide mechanical and thermodynamical stability and thus are physically meaning-ful. 5 Two fractional elements can be arranged in series (Generalized Maxwell model, GM) or in parallel (Generalized Kelvin-Voigt model, GKV). The calculations for the GM is shown later in the Supplemental Information (Eq. S.36 †). For the GKV, the elements are arranged in parallel, meaning the stresses add and the RCE and the expression of the complex shear modulus for the GKV is The simplification yields expressions (Eqs. S.28 †, S.29 †) with three independent parameters ( 0 , , ). While for the SFE (Eq. S.20 †) there is only one power exponent , the GKV shows two power exponents ( , ). Therefore, the GKV is able to model crossover events in which at a certain and only frequency the real part of the complex modulus equals the imaginary part. This crossover frequency for a set 0 = 1 s for the GKV is for > without loss of generality since and are interchangeable and ∈ (0.5, 1) and ∈ (0, 0.5). For the special case where + = 1 the crossover frequency is = 1 0 (1 rad s for 0 = 1 s). Later, in the ESI † (Fig. S.10 †) certain selected cases for different , are shown for the complex shear modulus obtained with the SFE, GM and GKV. Briefly, for frequencies ≪ the complex modulus behaves similar to for > 0 and for frequencies ≫ , ∼ . The special case of = 0 shows a plateau in the real part of the complex modulus for frequencies ≪ while the imaginary part is ∼ ∀ . The interpretation of the GKV parameters ( 0 , , ) are the following. 0 can be considered as an apparent shear modulus at the time scale 0 = 1 s similar to the interpretation of 0 for the SFE. However, note that GKV,0 = 2 ⋅ SFE,0 and is to be considered in a comparison. The power parameters and define the crossover frequency for a set 0 (see Eq. S.30 †). Therefore, the single information of one power parameter merely shows the complex shear modulus of the material over a certain frequency range, while the tuple ( , ) yields information of the complex shear modulus over the whole frequency range for a material with one crossover event.
The model using single fractional elements can be further extended by arranging multiple SFEs, such as the known mechanical models Zener model or Poynting-Thomson model with simple elements. However, this increases the amount of parameters and thus, it will not be used in this framework.
The third assumption for the microrheological analysis is dependent on the experiment. For experiments where the force is transmitted to the sample by an external bead neither the stress nor the strain are often directly accessible, but the force acting on the bead and the bead displacement instead. In the paper by Kollmannsberger et al. 6 , for a magnetic tweezers experiment, the strain is estimated as the bead displacement divided by the bead radius and the stress is estimated as the applied force divided by the bead cross section, i.e. the complex mod-ulus̃ ( ) =̃ ( ) ( ) 1 . While in the paper by Balland et al. 7 , for an optical tweezers experiment, the resulting complex modulus is calculated depending on the immersion of a bead inside the material characterized by an immersion half-angle , i.e. .

Generalized Maxwell Model
As mentioned, the underlying assumption is that the creep compliance behaves as a power law (Eq. S.16 †). In the following the results and calculations mostly based on 5  Note that for < the imaginary part is higher than the real part for the GM and vice versa for the GKV and accordingly, for > the real part greater than the imaginary part for the GM and vice versa for the GKV. The complex shear modulus cannot be described with one power exponent close to the crossover frequency. Figure  S.10B shows a case with = 0 and + ≠ 1. Here, the crossover frequency is not the same for the GM and the GKV. Also, for fre- Supplementary Figure S.11 † shows the fit to a real data set with the different models. Due to the existence of a crossover frequency, the SFE model was used separately for the real and imaginary part of the complex shear modulus.

Bead immersion half-angle correction
The complex shear modulus is defined as the ratio of the stress to strain. However, during the experiment only the force and the displacement of the bead with radius is measured. Here, to estimate the stress and strain the complex modulus is calculated with a correction ( ) based on the bead immersion half-angle into the cell̃  with the constant coefficients shown in Table S.11. The correction factor is also shown in Supp. Fig. S.12 † for both corrections. Due to the fact that the cell height at the location of the bead was unknown during the experiment, the correction 1 ( ) was used for the analysis. However, the height of HUVECs usually does not exceed 3 m at the cell body and is ≥ 1 m at cell periphery. 10 The difference between the other correction ( 2 , ℎ 2 ) can be considered as a possible error range due to the estimation of the stress and strain.

Comparison to optical tweezers
Phagocytosed 1 m polystyrene beads in HUVECs were used to measure the viscoelastic properties of HUVECs with optical tweezers. The measurement time for one bead was about 100 s. Searching the next bead and setting up the required conditions for the measurement required on average 6 min. Thus, measuring 15 beads would take about 115 min and these results would correspond to different times as the measurements are shifted. Using the AFS, however, we could on average measure 15 beads simultaneously, drastically increasing the throughput while also obtaining property values every 100 s for each bead at the same time.
As used in our results with the AFS, the results of optical tweezers measured on HUVECs can also be described using the generalized Kelvin-Voigt (GKV) model (Eq. S.27). Since we could only measure phagocytosed beads inside HUVECs, it is expected that the obtained results differ from the AFS where we had measured the cortex. However, the GKV model can well describe the measured data, see Supp.

Dynamics of the shear modulus with the GKV model
Images of the cells with beads during the control measurement are shown in Supp. Fig. S.14 †. Figure 7E shows the obtained parameter 0 using the GKV model. To test for a significant difference the obtained values are grouped into the time intervals in which the conditions had changed, e.g. a higher flow with cyto B insertion or the washing out process. The significance tests are shown in Supp.

Data post-processing and analysis
The obtained microrheological data has been post-processed using our software Kitsune to obtain correct values for the viscoelasticity. The software filters for inter-particle errors in the tracking, corrects for drift in the long time measurements and includes the local distribution of the conversion factor. Details about the analysis are shown in the next section. Briefly, the software corrects for the following three main problems: a) During the constant medium flow, other particles can interfere with the tracking of the beads that will cause erroneous values in the ( , , )-position data. Typically, these erroneous values only appear for a short time, i.e. one or two data points, due to the relatively high flow rate compared to the sampling frequency. This problem is addressed by replacing the erroneous value with the median value of data points temporally close to it.
b) The long measurement time may introduce mechanical drifts or allows drifts caused by active, but slow, cell movements over time. The drift is corrected by obtaining the drift profile with segmented, continuous fits of polynomial functions of second order and subtracting them from the data. This correction can be per-formed because the cell movement and activity themselves are not of interest in this case.
c) The lateral inhomogeneity of the force distribution inside the field of view will cause a wrong assignment of the actual force acting on the beads. This can be corrected by assigning a conversion factor of the amplitude to force for each bead according to their current lateral position with the obtained spatial map from the SFC, as described in the next section.

Detailed data post-processing and analysis
The relevant, measured data are the positional data in , , , the amplitude corresponding to the exerted acoustic force and the sampling time . A complete data set consists of the measured data in the full duration of acquisition, e.g. a duration of 3700 s. The following shows the automatized steps of the analysis using the self-written MATLAB software Kitsune (version 0.98).
• Find the start and end time of the mOsc corresponding to the data with force type = 5 in the modified LabVIEW software. (Full force span) • Segment the full force span in sections of 500 s evaluation time intervals in 100 s time steps.
• Analyze each of the sections separately with the following steps.
• Get the mean , -position during the analyzed section and calculate the force value according to the conversion factor from the calibration map.
• Subtract the force values of the section with its mean to avoid the errors from the zero-padding in a later step. This is valid because the oscillation is of interest and not the initial force step.
• Clean the -data from tracking errors, e.g. due to particles interfering with the tracking's region of interest (RoI) during flow. -Create the Vandermonde Matrix for a polynomial of 2nd order ( = 2).
-Constraint the fit: The starting point has to be the ending point of the previous Drift-Section; for the first Drift-Section the starting point of the data is also its starting point.
-Fit the data using lsqlin and polyval which are in-built functions of MATLAB.
-Append all Drift-Sections to a data set. (Drift-Correction) -Subtract the section with the Drift-Correction.
• Filter out sections (after the correction from above) that have high errors, e.g. a particle interfered with the RoI for too long or the tracked particle is no longer being tracked.
-Calculate the standard deviation (std) of the -position in the section.
-Set the results to NaN and stop the analysis if at least one of the following criteria are fulfilled:

SFC Protocol
This protocol is used to perform the Stokes Force Calibration (SFC) on Acoustic Force Spectroscopy (AFS) chips. It also serves as a guide. It is part of the Supplemental Information of the publication Multi-oscillation microrheology via Acoustic Force Spectroscopy enables frequency-dependent measurements on endothelial cells at high-throughput. The calibration method is described in the main text of the publication and the recommended software Kitsune to analyze the SFC data is available on github https://github.com/A141-User/Acoustic-Kitsune. This protocol is for the G2 AFS chips; other generations may differ.

Conditions of the main experiment
The conditions of the main experiment has to be set beforehand. Firstly, to ensure that the AFS can perform under the conditions and, secondly, to perform the SFC under the same conditions due to the dependency of the force on the respective conditions. Conditions: • Temperature ( < 40 • ) • Medium (liquid during experiment) • Bead size (preferably 0.5 m < < 20 m) • Bead type (e.g. polystyrene, silica, etc.) • Desired force range (pN-nN) • (optional) to think about: force type, e.g. constant force, force ramp, etc.
length of force application

Preparation
Prior to a successful calibration there are still a few aspects to be taken care of.
• This calibration method is based on video-capturing, meaning that a camera with a sampling rate of > 30 Hz is recommended.
• The surface of the AFS fluid chamber is made of glass. The beads in combination of the medium may bind to the surface, e.g. due to charge. Therefore, the glass or the beads can be coated accordingly to prevent specific binding on the surface.
• The bottom of the chip should be cleaned to prevent dirt that can disturb the tracking and may also alter the force values.
• The inside of the fluid chamber should be cleaned, e.g. with ethanol or standard bleach (e.g. A1727 Sigma-Aldrich) and thoroughly washed out afterwards.
• The position of the AFS chip is fixed and can only be moved in a controlled manner to ensure to measure at the same field of view (FoV). 5. wait until they are on the ground 6. apply a higher amplitude for 1 s 7. repeat steps 5. and 6. until a desired amount of amplitudes are measured (preferably 4 different amplitudes smaller than the highest possible to capture) 8. end measurement 9. repeat steps 1-8 until a desired amount of beads are measured, i.e. the measured beads are well-distributed throughout the whole FoV (which can be quickly checked using Kitsune's "Bead distribution") The experimental procedure is completed and the -positional data of every bead has been measured. For every bead, theposition has to be evaluated during the force application. The lateral location of each bead is known and the conversion factor is assigned to each bead to obtain the distribution of the conversion factor. In order to obtain the calibration map of the conversion factors, the distribution has to be interpolated. The analysis can be performed using the steps shown in the main text of the publication or using Kitsune.

Further steps
After obtaining the calibration map for one FoV at one specific condition, the calibration is basically successful. A few other steps can be performed, as listed in the following.
• find other FoVs and obtain the maps (calibration), make sure to assign the maps to the FoV accordingly • The analysis software Kitsune is available on request and is based on Matlab R2017 (The MathWorks). It is written for data generated by a modified version of the Generation 2 AFS LabVIEW tracking software. However, it can also analyze the standard version.
• The density of the bead has to be greater than the density of the medium under the conditions.
• If a force type other than the constant force or force ramp is needed, the LabVIEW software may have to be modified.
• If a high force amplitude is applied, the change of the temperature inside the fluid chamber during force application should be measured.