Viscoelastic parameterization of human skin cells characterize material behavior at multiple timescales

Countless biophysical studies have sought distinct markers in the cellular mechanical response that could be linked to morphogenesis, homeostasis, and disease. Here, an iterative-fitting methodology visualizes the time-dependent viscoelastic behavior of human skin cells under physiologically relevant conditions. Past investigations often involved parameterizing elastic relationships and assuming purely Hertzian contact mechanics, which fails to properly account for the rich temporal information available. We demonstrate the performance superiority of the proposed iterative viscoelastic characterization method over standard open-search approaches. Our viscoelastic measurements revealed that 2D adherent metastatic melanoma cells exhibit reduced elasticity compared to their normal counterparts—melanocytes and fibroblasts, and are significantly less viscous than fibroblasts over timescales spanning three orders of magnitude. The measured loss angle indicates clear differential viscoelastic responses across multiple timescales between the measured cells. This method provides insight into the complex viscoelastic behavior of metastatic melanoma cells relevant to better understanding cancer metastasis and aggression.

Eq. S2 and Eq. S3 allows each order of the variable s to be grouped, and the coefficients for each 49 series can then be calculated using known material parameters. If either the force or indentation is 50 provided, and the Lee and Radok assumptions still hold, calculating the material response simply 51 requires taking discrete derivatives of the AFM observables.

52
Since the material models have been described using and , Eq. S1 is transformed into

75
Intermediate combinations, such as the Standard Linear Solid (SLS) or Burgers models, still contain 76 relatively few parameters for describing samples but can be useful for initial approximations or 77 temporally limited datasets. For further discussion of these models, including the strengths and 78 weaknesses of each, the reader is directed to the available literature 3,4 .

112
The Generalized Kelvin-Voigt and Generalized Maxwell models are illustrated in 131 132

145
This description assumes that the model contains both a Glassy Compliance ( ) and a single 149 150

165
It is worthwhile to note that the existence of an elastic term and steady-state fluidity in the

183
There are several critical points to keep in mind when applying the parameterization methods

184
presented. First and foremost, the methodology does not naturally prevent "overfitting" a dataset; timescale gives a clear indication of when the minimum number of terms necessary for a close estimate occurs, but users must be careful to include only orders of magnitude that can reasonably be justified as being truly present in the dataset. For example, the likelihood that an indentation occurring for 1 millisecond contains mechanical artifacts from timescales on the order of 10 0 is low 191 and including a characteristic time on that order could obscure the effects of lower order terms. To

216
In extrapolating material behavior using parameterized models, it is critical to acknowledge the

248
Lastly, and most fundamentally, the assumptions presented during the framework derivations 249 must remain valid throughout the duration of the experiments. In particular, the linearity 250 assumptions dictate that the indentation depth must be "sufficiently small" such that the force  step-inputs) in stress or strain can require users to utilize a stiffness-or compliance-based model.

271
In such situations, the user must understand that an appropriate parameterization of one model,

272
given that it reproduces the input datasets to a high degree, will supply the parameters for its pair 273 through one of several methods:

275
1. Using the extracted parameters to "simulate" the response to a step function in stress or

282
If expediency is important, option 1 is convenient because it does not require conversion to the

283
Laplace domain for either model (although this step is analytically straightforward). In many cases, 284 option 2 is not feasible-beyond a single element, the relationship between models involves fourth 285 order polynomials in the parameters and is ill-posed. Even for two elements, the mathematical 286 treatment required would easily exceed the time or patience of most analysts.

287
The remaining option is to use the s-domain Collocation Method outlined by Tschoegl 3 , which 288 involves using the Laplace Domain forms of the Relaxance and Retardance (Eq. S8 and Eq. S10)

289
with the understanding that ( ) ( ) = 1. As such, evaluating or for a wide range of s-values  Eq. S16 Eq. S17

336
Providing Eq. S8 and Eq. S10 to a symbolic programming tool can thus dynamically generate the 337 coefficients required to simulate the following well-known differential representation of viscoelastic 338 materials: 339 340 Eq. S18 Here the variables M and N represent the highest order of s associated with the desired 341 polynomial series, and are denoted as separate exclusively for the case where u(s) and q(s) do not 342 contain the same order of s. The quantities and represent the coefficient for the ℎ order 343 derivative of the stress and strain, respectively, which can be obtained by rearranging the 344 polynomial series in the complex variable , as described above. Importantly, when the Generalized

345
Kelvin-Voigt representation utilizes a steady-state fluidity term, the polynomial will also feature an 346 order of s -1 , which in the time domain represents integration from the onset of indentation to the 347 current time. The stress and strain in Eq. S18 are then replaced by functions of force and 348 indentation. The following relationship can be used to generate AFM force curve data for the Lee

349
and Radok (Eq. S19) indentation framework: Eq. S19 Once the coefficients and are available, one can simulate the tip motion in time by

352
programmatically beginning with the AFM probe located at an initial height, and the enforcing an   Table S1.  Eq. S21 Alternatively, for the Generalized Kelvin-Voigt, the storage and loss compliance take the form:

488
Supplementary  where SSE denotes the Sum of Squared Errors. 512 513 parameters using an atomic force microscope and static force spectroscopy. Beilstein