Magnetic microrheometry of tumor-relevant stiffness levels and probabilistic quantification of viscoelasticity differences inside 3D cell culture matrices

The progression of breast cancer involves cancer-cell invasions of extracellular matrices. To investigate the progression, 3D cell cultures are widely used along with different types of matrices. Currently, the matrices are often characterized using parallel-plate rheometry for matrix viscoelasticity, or liquid-like viscous and stiffness-related elastic characteristics. The characterization reveals averaged information and sample-to-sample variation, yet, it neglects internal heterogeneity within matrices, experienced by cancer cells in 3D culture. Techniques using optical tweezers and magnetic microrheometry have measured heterogeneity in viscoelasticity in 3D culture. However, there is a lack of probabilistic heterogeneity quantification and cell-size-relevant, microscale-viscoelasticity measurements at breast-tumor tissue stiffness up to ≃10 kPa in Young’s modulus. Here, we have advanced methods, for the purpose, which use a magnetic microrheometer that applies forces on magnetic spheres within matrices, and detects the spheres displacements. We present probabilistic heterogeneity quantification using microscale-viscoelasticity measurements in 3D culture matrices at breast-tumor-relevant stiffness levels. Bayesian multilevel modeling was employed to distinguish heterogeneity in viscoelasticity from the effects of experimental design and measurement errors. We report about the heterogeneity of breast-tumor-relevant agarose, GrowDex, GrowDex–collagen and fibrin matrices. The degree of heterogeneity differs for stiffness, and phase angle (i.e. ratio between viscous and elastic characteristics). Concerning stiffness, agarose and GrowDex show the lowest and highest heterogeneity, respectively. Concerning phase angle, fibrin and GrowDex–collagen present the lowest and the highest heterogeneity, respectively. While this heterogeneity information involves softer matrices, probed by ≃30 μm magnetic spheres, we employ larger ≃100 μm spheres to increase magnetic forces and acquire a sufficient displacement signal-to-noise ratio in stiffer matrices. Thus, we show pointwise microscale viscoelasticity measurements within agarose matrices up to Young’s moduli of 10 kPa. These results establish methods that combine magnetic microrheometry and Bayesian multilevel modeling for enhanced heterogeneity analysis within 3D culture matrices.

: Comsol simulation of magnetic-field gradients. A-B: Simulation shows the degree of homogeneity for the magnetic-field gradients within the experimental workspace over the duration of the sinusoidal forces exertion (i.e. 20 s). Horizontal red lines indicate the boundaries of the magnetic hotspot, which is the experimental workspace. The simulation is based on [1] and uses rotational symmetry in cylindrical coordinates.
Materials and methods: data-processing pipeline. To measure microscale viscoelasticity, videos on magnetic and reference sphere displacements are captured. The videos needs to be transformed into the viscoelasticity data. The data-processing pipeline for the purpose is written in Python and operates mostly automatically. Still, user inputs are always required in the step 1 and in the steps 3-4 if the quality of the data is lower. The pipeline works as follows: (a) Manually initialize bounding boxes for magnetic and reference spheres (b) Tracker binarizes the bounding boxes with Otsu-Binarization (c) Tracker finds the biggest "blob" for each bounding box and computes the center point (d) The bounding box is centered based on the computed new center point (e) Repeat (b-d) over the entire video 2. Track association, involving with matching the IDs over repeated measurements. The same magnetic sphere should have the same ID over the repeats. This step is only necessary if the spheres are not tracked in the same order over the repeats (i.e bounding boxes are not drawn in the same order for each repeat).
(a) DBSCAN is used to cluster tracked the xy-coordinates (b) New IDs for each magnetic sphere based on the cluster assignment 3. Magnetic sphere radius estimation, including an estimate of the magnetic sphere radius.
(a) The first repeat of a measurement set should have a sweep video through the z-stack (b) Estimate the optimal binarization over the video (c) Binarize (d) Find the best focus by calculating the Laplacian of a bounding box and calculating the variance (for each magnetic sphere separately). The best focus should have the highest variance. i. If reference sphere is too far or right next to the magnetic sphere (i.e. we define that the distance range between a magnetic and a reference sphere is from 50 µm to 250 µm) ii. If amplitude is zero iii. If the shape of the sinusoidal signal is not symmetric enough (i.e. detrended displacement signal should have similar area for the two half periods) Materials and methods: the |G100| |G30| ratio. By using the given formula for absolute shear modulus (|G|) in Eq. 4, we can derive an approximation for the |G100| |G30| ratio, as followed: This information is plugged into the calculation of the ratio: • The ratio between the volumetric force constants:F volumetric,100 • We want to explore the |G| values at the resolution limit range of the microrheometer:p sphere,100 ≈p sphere,30 < 100 nm Therefore, we conclude that the 100 µm magnetic spheres are expected to measure approximately 8 times larger stiffness levels than the 30 µm magnetic spheres.
Materials and methods: model comparison. As described in the methods and materials section, the heterogeneity model (Model 1, Fig. 2) was compared against four other variants of the model (Models 2-5, Figs S1-4). The results are summarized in Table S1, where a higher expected log pointwise predictive density (ELPD) value denotes a better predictive performance. The Model 1 has been chosen as the primary model due to its ability to capture interpretable experimental design and has the secondhighest ELPD for |G| estimation. The Model 2 has the highest ELDP, which is however within the standard errors for the Models 1-2, and the heterogeneity quantification by the Model 2 is uninterpretable as described in Materials and Methods. The choice of the Models 1-2 for estimating |G| has a slight improving impact on the model performance based on the ELPD, while the most impacting aspect on the performance is to have a unique noise for each magnetic sphere. The unique noise term is needed, because particle tracking and radius estimation are done separately for each magnetic sphere producing different levels of noise. For estimating ϕ, adding hierarchy to the model has no noticeable effect on the performance based on this ELPD analysis (i.e. Model 1 has the lowest ELPD, yet, all the models have a close ELPD that are within standard errors of the Model 3). To summarize, the choice of the Model 1 was mainly done based on the performance to estimate the |G| and the interpretability (i.e. Models 2-5 are partially uninterpretable for the purpose).      Figure S6: Parallel-plate rheometry of viscoelastic properties of 3D-culture matrices with the spheres, and without the spheres. A-B: No significant differences in shear modulus (|G|) and phase angle (ϕ) are found between the matrices without the spheres, and with the magnetic/reference spheres. The data is for agarose matrices (a 0.9% concentration) and GrowDex matrices (a 1.25% concentration). Comparable viscoelasticity measurements have been previously established for accurate, cell-scale bulk properties of collagen matrices [2]. Because magnetic spheres can measure both the GrowDex and collagen matrices properties, it is reasonable to assume that GrowDex-collagen matrices are similarly measurable.
Materials and methods: calibration. For calibration, the initial estimation method is based on calculating the volumetric force constant as a mean value of each individual magnetic sphere's volumetric force, showing a method-related scatter of individual sphere results (Table S2). Here, instead, we use a more accurate, enhanced estimation method. Our approach is directly based on the definition of the force from Stokes' law and a linear regression model (Eq. 6-9). The data in Table S2 Table S2, this enhanced estimation improves reliability.  Figure S7: Graph of calibration values based on the initial and enhanced estimates ofF volumetric . The graph compares theF volumetric constants by the initial estimate (calculated separately for each sphere) and the enhanced estimate (calculated using Stokes' law). Both of the estimates provide similar mean values. Importantly, the graph illustrates how the enhanced estimate has consistently lower uncertainty than the initial estimate. The estimates show that theF volumetric values are increasing linearly as a function of the current i grad for the 30 µm spheres. Each dataset is either for the 30 µm or the 100 µm spheres.   Figure S8: Viscoelasticity of the matrix types based on microrheometry and parallel-plate rheometry. A-B: Absolute shear modulus (|G| ) and phase angle (ϕ ) results of the microrheometry measurements for the used matrices types are shown. The compositions of the matrices are as followed: GrowDex has a concentration of 1.25 %; GrowDex-collagen has a GrowDex's concentration of 0.45 % and a collagen concentration of 2.0 mg/mL; fibrin has a concentration of 30 mg/mL; and agarose has a concentration of 0.5 %. For both plots, the dark grey spheres show values measured by individual magnetic spheres, while the red squares are the mean values of the spheres. Standardly, the plots' box and whiskers indicate the lower/upper extremes, the 25 % and the 75 % percentiles, and the median. C-D: Absolute shear modulus (|G|) and phase angle (ϕ) results of the parallel-plate rheometry measurements for the used matrices types are shown. Black vertical dashed line indicates the timepoint when the microrheological measurements have been started for all matrix types apart from the fibrin matrix. Fibrin-matrix measurements were started after 30 min. The curves with darker colors indicate the mean values of the data, while the lighter shaded-color regions indicate the lower 5 % and the upper 95 % percentiles for the data. Figure S9: Absolute shear modulus in relation to the displacement signal amplitude. Absolute shear modulus (|G|) has a multiplicative inverse dependence on the amplitude of the displacement signal (p sphere ). An increased uncertainty of |G| at low displacement amplitude values is shown (i.e. smaller differences in the displacement amplitude cause larger differences in |G| at lower displacement amplitudes compared to larger displacement amplitudes). Figure S10: SEM micrographs of agarose matrices with the 30µm and the 100µm spheres. SEM micrographs illustrate the occasional deviation of (A:) the 30µm and (B:) the 100µm magnetic spheres from a perfect spherical shape, which is expected to be an error source in our microrheometry methods and detectable in the calibration.