Stretching and heating cells with light—nonlinear photothermal cell rheology

Stretching and heating are everyday experiences for skin and tissue cells. They are also standard procedures to reduce the risk for injuries in physical exercise and to relieve muscle spasms in physiotherapy. Here, we ask which immediate and long-term mechanical effects of such treatments are quantitatively detectable on the level of individual living cells. Combining versatile optical stretcher techniques with a well-tested mathematical model for viscoelastic polymer networks, we investigate the thermomechanical properties of suspended cells with a photothermal rheometric protocol that can disentangle fast transient and slow ‘inelastic’ components in the nonlinear mechanical response. We find that a certain minimum strength and duration of combined stretching and heating is required to induce long-lived alterations of the mechanical state of the cells, which then respond qualitatively differently to mechanical tests than after weaker/shorter treatments or merely mechanical preconditioning alone. Our results suggest a viable protocol to search for intracellular biomolecular signatures of the mathematically detected dissimilar mechanical response modes.


Introduction
A common difficulty in medical diagnosis and therapy of intracellular malfunctions is that symptoms may manifest themselves on a much higher level of organization, such as organs or larger parts of the body. Also certain therapeutic practices operate on macroscopic tissue rather than on a cellular level. Physiotherapy, for example, often conjointly employs intense massage and thermotherapy (tissue heating by several degrees over several minutes) to synergistically soothe body stiffness and relieve muscle spasms. A microscopic or even molecular understanding and objective quantification of the effectiveness of such treatment seems elusive, as it would require one to monitor biomolecular and intracellular processes during macroscopic treatment.
Here, we propose that an important first step in this direction can nevertheless be taken, by applying prototypical physiotherapeutic protocols to individual isolated cells and by quantitatively comparing and interpreting their effects onto the mechanical state of the cells. To this end, we employ two different optical stretcher (OS) setups that allow stretching forces with and without concomitant heating to be applied to a large number of identically prepared individual cells. The cells are trapped by two divergent laser beams in a micro-fluidic device, under well-controlled and reproducible conditions. An advantage of the setup is that the heating is not applied slowly or permanently, as is most commonly the case in other methods [1][2][3][4][5][6], where the cells have time to slowly adapt their internal state to the altered conditions. Rather, the heating occurs fast (within 200 ms) and in perfect synchrony with the mechanical deformation, and both can be adjusted in a range comparable to what cells might experience in everyday life or during conventional physical therapies. The protocols allow us to clearly discern the reversible short-time effect and the long-time impact of heating. Further, we can explore experimental conditions that transcend the framework of (linear) thermorheology previously successfully applied in cell mechanics [7,8].
While our technical setup conveniently puts us into the position to both apply a variety of treatments and quantify the resulting mechanical states of the cells, we (still) lack a reliable quantitative real-time control of the underlying intracellular biomolecular processes. One could imagine to microscopically label and track a selected number of molecular responses. But, apart from obvious limitations to its precision, this approach would clearly require one to specifically know in advance what to look for. Here, we take a different route in order to precisely monitor and quantitatively characterize the mechanical state of a cell. Namely, we employ, as a second crucial ingredient apart from the cell stretcher, a schematic theoretical model that i) can quantitatively fit the experimental data over a wide range of conditions and protocols, and ii) allows us to discriminate two qualitatively distinct but experimentally intermingled types of mechanical response, called viscoelastic and inelastic, respectively.
In this context, 'viscoelastic' refers to a viscously damped, but otherwise reversible, elastic response that is amenable to viscoelastic (spring-dashpot-type) modeling. As 'inelastic' we characterize a component of the response that involves a substantially longer memory that is most plausibly due to some form of structural damping [9,10]. Inelastic deformations typically involve the breaking of weak, transient bonds between major load-bearing elements in the cytoskeleton or other major structural changes. As our data indicate, the corresponding relaxation or 'structural healing' proceeds indeed very slowly and takes much longer than the viscoelastic relaxation, so that it cannot be modeled in the same way. It shows almost no tendency to saturate on the times (of about a minute) probed in our experiments. For living cells, which undergo permanent active remodeling, such slow processes are usually not amenable to any purely passive modeling and we do not attempt this.
As the reader may anticipate from the above characterization, our theoretical model is deliberately schematic and phenomenological, and indeed meant to be the simplest conceivable model that enables us to represent and quantitatively analyze the general physical processes that we suspect to underlie the dominant cell mechanical response. Nevertheless, we want to argue that this simplistic level of modeling is most suitable for our purpose. Any more detailed model would, first of all, require stronger a priori assumptions about the relevant molecular mechanisms, similar as in the mentioned molecular labeling scenario. Secondly, it would almost surely involve additional unknown parameters. While some practitioners may see this as an advantage in terms of the model's versatility and capability to produce accurate fits to data, it would in fact lower its explanatory power. With growing complexity, phenomenological models commonly tend to run into the problem of fit redundancies and overfitting, and they do so much more quickly than widely recognized [11][12][13]. In contrast, our coarse minimalistic model allows us to, first of all, consistently and accurately parameterize the cell-mechanical response for a wide range of system parameters and thermorheological protocols; secondly, to disentangle (and make robust qualitative statements about) two qualitatively distinct underlying microscopic mechanisms. Finally, we provide new insights concerning the ongoing debate whether heating stiffens and solidifies [14] or fluidizes and softens cells [1][2][3][4][5][6]. As we detail below, our analysis moreover allows us to draw non-trivial conclusions as to the kind of stretching and heating protocols that are most effective in firmly altering the state of a cell, and the thresholds that therefore need to be overcome. In particular, we can confirm, on the level of individual cells, the common folklore that physical therapy ought to be applied with a certain vigor and persistence to be most effective. Which suggests that the notion of 'light physiotherapy for single cells' might provide a useful qualitative analogy to guide the reader through the remainder of the paper.

Modeling
As a starting point for the theoretical description of the experimental data, we use the phenomenological glassy wormlike chain model (GWLC) [15][16][17]. The GWLC is a viscoelastic, meanfield-type model which has been introduced as an extension of a standard polymer model, namely the wormlike chain (WLC). It links the well-understood dynamics of the polymeric constituents of the cytoskeleton to the rheology of cytoskeletal networks and whole cells in a minimalistic but quantitatively accurate way. To this end, it represents the interactions of a test chain with its complex surroundings in a schematic but effective way.

WLC and GWLC dynamics
The dynamics of a semiflexible polymer of contour length L is represented via an over-damped Langevin equation for a wormlike chain in the weakly bending limit, where a viscous drag force is balanced by two elastic contributions and a random force representing thermal noise: Here, the vector filed r ⊥ (s, t) is the time and arclength dependent transverse excursion from the straight ground state, ζ ⊥ is a transverse friction coefficient per length, κ is the bending rigidity and f is an optional uniform backbone tension. The random force ξ ⊥ is modeled as Gaussian white noise characterized by its mean and variance: Equation (1) can be solved via an eigenmode decomposition ansatz: For hinged boundary conditions the mode functions W n (s) are sine functions: W n (s) = 2/L sin(s · nπ/L), the Langevin equation for the nth mode amplitude function A n then reads: Here the vector notation has been dropped, because the components of r ⊥ are independent of each other, and the random force ξ ⊥ has likewise been separated into its Cartesian components that obey . For the construction of the GWLC, we reformulate equation (4) by introducing the WLC mode relaxation times τ WLC n : To include the interactions of a test chain with its slowly evolving surroundings, an exponential stretching of the WLC mode relaxation spectrum τ WLC n is employed in the GWLC model. In particular it is assumed that: (a) A test polymer feels its surroundings only at local interaction points, which are homogeneously distributed along the polymer's contour and separated by a distance Λ called the interaction length (b) Every interaction with a contour element of length λ n slows down the WLC mode relaxation times τ WLC n by an Arrhenius factor exp(ε) Thus, the standard WLC mode relaxation times are modified according to: Here λ n = L/(πn) is the (half-) wavelength of the nth mode, N n = λ n /Λ − 1 is the number of interactions associated to the nth mode, and ε is the stretching parameter that controls how much the modes of the test polymer are slowed down by its surroundings. It can loosely be interpreted as the characteristic height of typical free energy barriers in units of k B T in a rough free energy landscape, or simply as a phenomenological measure of the strength of the mutual caging of the polymers.

Connecting the mesoscopic and continuum description
To relate the mesoscopic, meanfield-type description of such a caged test polymer to the overall stress σ applied externally to a polymer network or even to a whole cell and to the observed overall strain γ, we assume affine deformations, i.e., we take a shortcut through the complex geometry by linking the applied stress σ and the backbone tension f felt by the test polymer via an effective mesh size Ξ, which is a convenient and intuitive way to express the density of load-bearing strands: Accordingly, we relate the relative change of the force-dependent end-to-end distance: of the test polymer to the observed strain: where d 0 is the initial end-to-end distance. The response of d(f, t) to a constant step force f can be calculated from the governing Langevin equation. The result reads: where we have introduced the GWLC relaxation rates ω n = 1/τ GWLC n and the initial values A 0 n of the amplitude functions, which, for an initially equilibrated chain are given by equipartition:

Thermorheological softening and structural remodeling
There are several pathways to include thermorheological effects into the GWLC. One pathway, which is suitable for short-time dynamical effects, is to allow characteristic model parameters such as ε, Λ, or Ξ to become (adiabatic) functions of the external fields [15,18,19]. This approach captures internal structural changes that follow the external fields and immediately and fully relax upon reversing or switching off the external field. In this sense, these remodeling effects still fall within a conventional non-linear viscoelastic paradigm.
To include them into our description, we allow for an adiabatic dependence of the GWLC relaxation spectrum τ GWLC n and the effective mesh size Ξ on the imposed stress and temperature protocol. To reduce the number of parameters, we further require a large number of interactions per test polymer N n 1 such that the slowdown of the mode relaxations is characterized by one single number S that we call 'solidity', defined via In summary, viscoelastic remodeling effects are represented in our description by allowing that: (a) The effective mesh size Ξ (b) The solidity S may take different values at different externally applied stress σ and temperature T. The two parameters represent different physics. The length Ξ predominantly characterizes structural and geometrical aspects. Its reduction in response to an external field has a rather straightforward interpretation, namely as a force-or temperature-induced node-dilution in the relevant cytoskeletal network, which thereupon softens 'geometrically'. The dimensionless solidity S comprises information about the strength of the caging of the test polymer. To discriminate between these two structural changes that fall into the viscoelastic category, we propose to use the term fluidization for a reduction of S and reserve the notion of geometric softening for the above-mentioned network-dilution effects controlled by Ξ.
The solidity is affected both geometrically, namely via variations of the average distance Λ between interaction points along the test chain, and energetically, by changes of their strength, quantified by e ε . If ε is viewed as an activation energy, it might be directly temperature dependent [19][20][21] and moreover affected by a tilting of the free energy landscape in the spirit of a generalized Kramers escape rate model [22]: Here, δ is the characteristic width of the free energy wells and barriers. The minus sign corresponds to the intuitive case that an external force lowers the energy barrier, while the plus sign can account for cases in which the bonds can become reinforced in response to an external force (so-called catch bond behavior). Beyond the above adiabatic or viscoelastic processes, the GWLC can accommodate more severe and longer-lived structural remodeling effects described by additional dynamic processes, in which case the model is referred to as the inelastic GWLC [9,18,23]. It then falls into a more general category of material modeling, which also includes (visco-)plastic and pseudoplastic effects, typically associated with the breaking and reforming of weak reversible bonds, cross-links or entanglements. Their mathematical description introduces additional slow dynamic variables into the discussion that may represent both truly (visco-)plastic effects as well as pseudo-plastic effects that may ultimately be reversible due to some self-healing processes. In the context of biomechanics, both effects are often impossible to distinguish on the investigated timescales, or the distinction is of no practical relevance, because the systems actively remodel.
We refrain from entering the level of inelastic modeling in what follows and employ a different strategy instead. We try to describe the data with the discussed reversible, viscoelastic GWLC. A failure of this attempt then provides strong evidence for the occurrence of long-lived effects that fall into the wider inelastic category. We chose this strategy as it has several advantages over strategies including more complex models from the outset: in the present contribution, we are more interested in detecting evidence for long-lived remodeling and not so much in any detailed mechanistic modeling of it. This saves us from the temptation to include more fit parameters and allows us to keep the analysis minimalistic. Moreover, current inelastic models like the inelastic GWCL are meant to describe force-induced bond breaking, while the very pronounced (apparently anomalous) remodeling effects observed in actin networks and living cells [1-3, 7, 8, 14] are not easily included.

Methods and materials
While the OS has been, in the past, used quite extensively to study more applied biological questions, and to lay the basis for potential future medical applications of cell mechanics, in this study we investigate the thermomechanical properties of suspended cells with a photothermal rheometric protocol that can disentangle fast transient and slow 'inelastic' components in the nonlinear mechanical response.

Setup and experiments
A microfluidic version of the OS was used to investigate the mechanical deformation of cells upon optical stress. In brief, the OS is a dual-beam laser trap (DBLT) capable of trapping and deforming cells through optically induced stress acting on the cell surface. Two optical fibres placed against each other are aligned perpendicular to a square glass capillary (figures 1(A) and (B)). The cells in suspension are delivered into the trapping region through the glass capillary. The flow into the glass capillary is adjusted by the relative difference in heights of an inlet and outlet reservoir connected to the capillary (figure 1(B)). Detailed descriptions of the OS working principle and setup can be found in [24], respectively [25][26][27].
To investigate the influence of mechanical stress and transient heat on the cells, two OS setups with different laser wavelengths (1064 nm and 780 nm) were employed to stretch the cells ( figure 1(A)). The first OS was mounted on an inverted microscope (Zeiss Axiovert 200M) equipped with a LD Plan-NEOFLUAR Ph2 40x/0.60 NA objective. The laser used was a single-mode, continuous-wave fibre laser at a wavelength of λ = 1064 nm (YLM-5-1070-LP; IPG Photonics, Oxford, MA) which in addition to the applied optical stress causes a concomitant temperature increase of 13KW −1 [28]. The second OS device was mounted on an inverted microscope (IX71; Olympus, Melville, NY) equipped with a Plan Fluor 40x/NA 0.75 objective (Olympus). The laser used was a single-mode, continuous-wave fibre laser at a wavelength of λ = 780 nm (Eylsa 780; Quantel, Les Ulis, France) which in addition to the applied optical stress causes an estimated [29] concomitant temperature increase of only 1KW −1 . Protocols carried out with the 1064 nm laser are called heating protocols and protocols carried out with the 780 nm are call non-heating protocols. CCD cameras (AVT MARLIN F-146B; Allied Vision, CA) were attached to the microscopes for image acquisition. For data acquisition and analysis, custom-built LABVIEW software (National Instruments, Austin, TX) was used to track the cell shape during the stretching. For every measured cell, the deformation along the 'major axis' r(t) was stored for every time-frame. 'Major axis' and 'minor axis' are the dimensions of the cell along and perpendicular to the laser optical axis, respectively. The initial radial 'major axis' r 0 was measured during the trapping period (figure 1). The time-varying axial strain: was then evaluated accordingly.
During the experiments, cells were exposed to optical stress pulses of different duration according to the employed protocol (figures 1(E) and (F)). We conducted single-and multiple-stress-pulse experiments. For all protocols and experiments, cells were initially trapped by applying a small amount of optical stress for 1 s. For single-pulse experiments, one stress pulse of 4 s duration (figure 1(E)) and varying laser power was applied after the trapping phases. For our multiple-pulse experiments, two types of pulse protocols were applied: a short-pulse protocol of 10 cycles, each consisting of 0.5 s-long stretching phase and 3.5 s-long relaxation phase, and a long-pulse protocol of four cycles, each consisting of a 4 s-long stretching phase and 6 s-long relaxation phase (figure 1(F)). The laser stretching power was fixed to 1 W in total (0.5 W per fiber) for all multiple-pulse protocols. We investigated two different cell lines: HL60 cells cultured in suspension and HeLa cells cultured adherent to a substrate. To test the contribution of the actin cortex to the overall mechanical properties of the cells, we treated cells with cytochalasin D (CytoD; Sigma-Aldrich). CytoD is an inhibitor of actin polymerization and leads to a destabilization of the cell actin cortex. For the multiple-pulse experiments, all protocols are listed together with the number of measured cells for each

Verification of the amount of heating by the OS
To verify previous results [28,29] on the temperature changes in the OS, we estimated the temperature induced by the optical trap. To do so, we employed poly N-isopropylacrylamide (pNIPAM) microgel beads as thermoresponsive [31] probes. The pNIPAM beads were trapped at low power, P = 0.1 W (total power; 0.05 W per fibre), and a higher power single-pulse was applied identically as in the experiments with cells.
The highest available total power of P = 1.5 W was used in the case of the 780 nm lasers. The increase in laser power induced a noticeable, but minimal, size reduction of the bead (first line of the figure 1(D)). For the case of the 1064 nm laser a total power of P = 0.6 W was enough to induce total collapse of the pNIPAM bead (second line of figure 1(F)). These experiments were carried out at a room temperature of T = 23 • C. A full characterization of our pNIPAM beads is the subject of a different publication, but we found that the critical temperature (temperature at which the beads collapse) is T = 31 • C. Furthermore, the ΔT = 8 K induced by the 0.6 W power of the laser would mean that the temperature increase in our 1064 nm OS setups is 13.3KW −1 . This result is in good agreement with [28].

Sample preparation
HeLa cells (cultured adhered to a tissue culture substrate) were detached from the tissue culture flask and transferred from culture medium to PBS. The cells, of both cell lines, were kept for about 30 min in PBS before the start of the experiment to allow them to adapt to the new conditions. For CytoD treatment, cells were suspended in PBS and incubated with 0.1 μM of CytoD for 30 min at 37 • C prior to the experiment. Experiments were performed at a room temperature of T = 22 • C.

Data selection and noise level estimation
Most cells deform and relax maintaining their volume and their elliptical shape when exposed to periods of optical and thermal stress (figure 2(A), first row). However, some of the studied cells reacted to the stress by spontaneous development of shape irregularities (figure 2(A), second and third row). For cells that showed very pronounced shape irregularities, the cell axis parallel to the optical axis became smaller than the cell axis which is perpendicular to the optical axis. Further, they became unstable in the trap, leading to rotation (figure 2(A), third line). Since, for these cells, the acquisition of an unique major axis length over the entire experiment was impossible, they were not considered in our analysis. The frequencies of cell behaviors observed in the various protocols is presented in figure 2(B). Also, for some experimental protocols (HeLa cells probed by the non-heating protocols) no detectable deformation in response to the applied stress could be found by visual inspection of the deformation curves. As a basis for a quantitative discussion of this behavior, we estimated the noise level of the measurements. For the estimation we used two different methods. On the one hand, deformation data from the initial trapping phase was used to generate deformation histograms. For the histogram for HeLa cells, we used 6734 relative deformation values from 637 cells and for HL60 cells we used 10 292 data points from 857 cells. The resulting histograms are approximately fitted by Gaussian distributions ( figure 3). From the fits, we obtained a mean deformation close to zero and a standard deviation of 0.16% for both HeLa and HL60 cells.
For the second method, we acquired an additional data set where 10 cells were each trapped for 7 s and again no additional optical stress was applied. For each cell, the measured deformation in pixel was converted into relative deformations γ. In total N = 700 deformation values were acquired. As an estimate of the noise level we then computed the mean absolute error, MAE, of the first N − 1 deformation values and the last N − 1 deformation values: which yields a noise level of 0.3%. Both methods are in good agreement when we use two standard deviations as a conservative noise level estimate in the first method. We, thus, chose 0.32% as an upper bound for the noise level.

Parameter values
All model parameters are summarized in table 2. We used only 3 free fit parameters, namely the 'native' solidity S 0 in absence of external force and temperature fields, its reduction ΔS during the stretching phase, and the force f acting on the test polymer. All other parameters were kept constant. While it is intuitively clear that S 0 and ΔS should be system-and protocol-dependent, this is actually also the case for the force f, which depends on the applied optical stress and on the effective mesh size Ξ. Because of the coupling of Ξ and f via the externally controlled stress in equation (7), a detected change of f is indicative of a corresponding reciprocal change of Ξ (i.e., increase in f → decrease in Ξ).
Besides these fit parameters, the persistence length l p of the test polymer is of particular relevance, because it allows to theoretically discriminate between viscoelastic contributions arising from different cytoskeletal polymers. For cells that were not treated with CytoD, we used a persistence length of l p = 10 μm, i.e. we assumed the external stress to be predominantly borne by actin filaments. In contrast, for cells treated with CytoD, which is known to impair the formation of a proper actin network, we used a smaller value l p = 3 μm to acknowledge that part of the stress is now transmitted via intermediate filaments with a typical persistence length of l p = 1 μm [32]. Thus, the effective test polymer is thought of as an interpolation between an actin filament and an intermediate filament. A quantitative change of the l p value (like l p = 2 μm instead of l p = 3 μm) did not affect the fit results qualitatively, but only rescaled the values of the fit parameters. However, employing the persistence length of pure actin or pure intermediate filaments, did not allow to describe the data consistently.

Fitting procedure
Although we have reduced the parameter space to only 3 fit parameters, there is still some parameter indeterminacy and sloppiness left in the model. This means that the effect of changing one parameter can to some extent be emulated by an appropriate change of the other two parameters [12,13]. Thus, fitting the model to the data is a nontrivial task, as there is a considerable risk of being stuck in a local maximum of fit fidelity. In addition, we are not aiming for fits of individual data sets, only, but search for a consistent description of groups of related data sets. Consistency in this context means that the mathematical description of each of the investigated systems obeys the following consistency conditions, even when probed using different experimental protocols (i.e. different pulse duration or different laser wavelength).
(a) S 0 does neither depend on the used laser wavelength nor on the pulse duration (b) ΔS, f do not depend on the pulse duration Condition (a) implies that during the relaxation phase, when the laser is switched off, the structural state quantified by S returns to the same value independently of whether the system has been heated or not in the previous stretching phase. Condition (b) means that the imposed experimental time scale does not influence the amount of thermorheological softening. Such an effect would be considered indicative of inelastic processes.
To cope with the described problems of sloppiness and consistency, we used a three-step fitting procedure. First, the model was evaluated on a parameter grid and the sum of squared residuals were calculated for each data set at each point of the parameter grid. With this step, we tried to minimize the risk of getting stuck in a local maximum of fit fidelity. Secondly, we searched the combined parameter space of several data sets and listed the results with the lowest sum of squared residuals that also obeyed the consistency conditions. In the third step, we compared fits of similar quality and chose the one providing the most natural and straight-forward interpretation. We preferred fits for which one parameter did not have to be changed for the description or took identical values in two or more data sets.

Single-pulse experiments
We started our study by probing HL60 cells with a standard OS protocol [7,8]. The creep response to a single stress-pulse followed by a relaxation period was observed. A typical curve is shown in figure 1(E). We varied the laser wavelength and stretching power and compared the peak compliance defined as the ratio of strain at the end of the stretching phase and global geometrical factor (GGF) [34][35][36], which is a measure of stress which also takes into account the several parameters characterizing the experimental conditions. The results are shown in figure 4. When the cells were stretched by the 780 nm laser, the compliance did not depend on the stretching power. We thus conclude that we operated in the regime of linear rheology and that the mild temperature increases (1KW −1 [29] ) induced by the 780 nm laser were insignificant. For the 1064 nm laser, which induced a stronger temperature increase (13KW −1 [28]), we found a clear increase in compliance with laser power. The fact that there was no temperature dependency of the compliance in the presence of low heating and the same amount of applied stress, suggests that we observed temperature-induced softening. We note that the temperature effect was far too pronounced to be captured by the explicit T-dependence contained in the model.
In physiotherapeutic terms, the above could be summarized by saying that we applied treatment of different stretching intensities combined with or without subsidiary heat application, and confirmed that subsidiary heat application helps to reduce cellular stiffness [7,8].

Multiple-pulse experiments
To go beyond the results of the single-pulse experiments and search for evidence of long-lived structural changes, we performed multiple-stress-pulse experiments with several consecutive cycles, each consisting of a stretching and a relaxation phase. We probed HL60 cells cultured in suspension and HeLa cells cultured adhered to a substrate both with and without CytoD treatment, using short-pulse and long-pulse protocols with and without concomitant heating. An overview of the protocols is given in table 1 and figure 1.
The resulting deformation curves are shown in figure 5. Several conclusions can be drawn directly from the deformation curves. Concomitant heating of cells during the stretching phase leads to more than twofold increased deformations, which confirms the temperature-induced softening found in the single-pulse experiments. Similarly increased deformations are observed when long stretching cycles instead of short ones are applied. When comparing the investigated systems, two stable trends can be observed. Firstly, the treatment with CytoD leads to more deformable cells. Secondly, HeLa (cultured adherent to a substrate) cells are much stiffer than HL60 (cultured in suspension) cells. When probed using the heating show the results for cells treated with CytoD (indicated by the gray background). When HeLa cells were probed by the non-heating protocols we could not detect a deformation in response to the applied stress. An example of this behavior given in (H) (blue line shown for the first 7 s). The overall noise floor of 0.32% is shown as a green band throughout. All fits are of high fidelity indicating that data sets that we tired to fit are well described by the reversible viscoelastic GWLC. protocol (1064 nm laser) the observed deformations of HeLa cells are roughly twofold decreased compared to the HL60. When the non-heating protocol (780 nm laser) is used, deformability is even lower, as the deformation curves fall beneath the noise level. For all other measurements, the general phenomenology is always similar. A creep response is followed by an incomplete relaxation such that residual deformation builds up from cycle to cycle.
To draw further conclusions, in particular concerning long-lived structural remodeling, the data was analyzed using the GWLC model. We tried to find a consistent description which includes both fluidization and geometric softening effects of the type introduced above. A failure of the viscoelastic model to describe the experimental data well is then a strong indication for the occurrence of long-lived inelastic remodeling effects.
We divided the data sets into three groups and compared two fitting procedures. Group 3 contains the data sets where we could not detect a clear signal above the noise level (HeLa cells probed using the non-heating protocols). Group 2 contains three sets (HL60 cells both with and without CytoD treatment probed using the long-pulse heating protocol and HeLa cells treated with CytoD and probed using the long-pulse heating protocol), which are the candidates for inelastic effects as they show deviations from the typical shape of all other deformation curves, which we gather in Group 1. In a first fitting attempt, we tried to find a consistent fit of high fidelity for all data sets in Group 1. The corresponding fits are shown in figure 5. The fit parameters are summarized in table 3. All fits are of high fidelity. Thus, we conclude that these systems are dominated by viscoelastic effects and exhibit no signal of long-lived inelastic effects.
In a second fit attempt, we tried to find a consistent fit for all data sets of group 1 and 2, i.e., we also included the data sets corresponding to deformation curves with somewhat irregular shapes (red curves in figures 5(C) and 6(C), figures 5(G) and 6(G), figures 5(H) and 6(H)). The resulting fits are shown in  The overall noise floor of 0.32% is shown as a green band throughout. Overall the fit fidelity is still good, while for HL60 cells treated with CytoD and probed using the short-pulse heating protocol (red curve in (E)) it is only fair, and for untreated HL60 cells probed using the long-pulses heating protocol (red curve in (C)) as well as for HL60 and HeLa cells treated with CytoD and probed using the long-pulse heating protocol (red curves in (G) and (H)) the fit fidelity is at best mediocre, indicating the occurrence of long-lived inelastic effects, not captured by the reversible, viscoelastic GWLC. figure 6. While fits for the data sets in group 1 are of high or fair fidelity, the best fits that we could achieve for the data sets in group 2 were of mediocre fidelity. Together with the results from the first fitting attempt this implies that the model, including only adiabatic fluidization and geometric softening, is not able to describe the data sets in group 2, which in turn is strong evidence for the occurrence of long-lived inelastic structural rearrangement effects. This suggests that, in light-physiotherapy of single cells, more vigorous treatments, such as longer stretching cycles and additional heating, are required in order to induce long-lived structural changes.
The model moreover allows us to quantitatively compare the different protocols and investigated systems in group 1. Across all protocols we find viscoelastic, thermorheological softening effects that can be accounted for by either geometric softening or geometric and fluidization effects, both implemented as adiabatic parameter changes.
When comparing the non-heating protocol (780 nm laser) to the heating protocol (1064 nm laser), we find that the higher deformations induced by additional heating are chiefly caused by larger forces f, which are increased roughly between 2 and 6 times. This can be interpreted as an increased amount of geometric softening, as a larger mesh size Ξ implies an increased force when the external stress is fixed, as in our case. The effect of additional heating on the amount of fluidization is not as pronounced and more subtle: in HL60 cells treated with CytoD, additional heating leads to increased values of ΔS, i.e., to more fluidization. In the model, there is a natural interpretation for this effect: intuitively, an increase in Ξ (geometric softening) in response to additional heating should be associated with an increase in Λ as both represent, in a sense, the network density, while ε should decrease (equation (13)). Thus, one would expect that S decreases with increasing Λ and Ξ, or that a sparser network deviates more from classical rubber elasticity. However, for untreated HL60 cells we find the opposite effect, namely that the reduction in solidity during the stretching phase ΔS is decreased in the presence of additional heating. One could try to explain this effect by the fact that Λ and Ξ could be two distinct characteristic length scales capturing different aspect of the complex network geometry. Another possible interpretation of the diverse effects of concomitant heating is based on temperature-activated actomyosin activity [37][38][39]. This might create additional tension, which may contribute to the increased force values inside the cell and could (partially) compensate the fluidization effects, in analogy to the more solid-like behavior of adherent cells induced by actomyosin activity [14]. This interpretation resonates well with the observation that the effect vanishes upon CytoD treatment and also with a recent study [40] providing evidence that myosin II activity can make suspended cells appear softer. Besides network dilution effects and temperature activated actomyosin activity, there are other possible interpretations. Some crosslinkers have been reported to become stronger at higher temperatures [20] or show catch bond behavior [41] which could increase S via increased values of ε, overcompensating the decrease in Λ. Lastly, we cannot exclude the possibility that the physical mechanisms of the pronounced temperature sensitivity of actin [19,42] might not be captured accurately by the model. On the present level of our understanding, we can neither rule out one of the interpretations nor can we exclude that we observe a superposition of different effects. Such redundancies, allowing for diverse interpretations of cell rheological data, are not entirely unexpected, since they already afflict the rheology of highly purified in-vitro cytoskeletal assays [42].
Concerning the effect of CytoD we find two trends. On the one hand, CytoD makes cells conform better with the rubber-plateau response predicted by classical polymer models as indicated by a substantial increase in solidity (roughly between 2.5 and 5 fold). On the other hand, CytoD reduces the temperature sensitivity of the cells, since the differences in ΔS and f between the heating and the non-heating protocol become smaller. These results resonate well with the assumption that, for cells treated with CytoD, a major part of the stress is transmitted via intermediate filaments, if intermediate filaments are less temperature sensitive than actin filaments and more in line with classical polymer models. This interpretation is consistent with studies on reconstituted biopolymer networks which indeed indicate a high temperature sensitivity of actin [19,42] and a more rubber-like behavior at least of some intermediate filament networks compared to actin networks [43,44]. Moreover, the effects of CytoD treatment can be interpreted consistently as a reduction in actomyosin activity which can lead to more rubber-like behavior in suspended cells [40].
When comparing the fit results for the two cell lines HeLa (cells cultured adherent to a substrate) and HL60 (cells cultured in suspension), the higher stiffness of HeLa cells is mostly attributed to the at least 80% higher solidity, while the force f acting on the individual test polymer is decreased slightly (differences are below 20%) suggesting that the cytoskeletal network in cells that are cultured adherent to a substrate (HeLa) are only moderately denser but more tightly interconnected compared to cells that are cultured in suspension (HL60).
To investigate why we did not observe deformations above the noise level for HeLa cells in group 3 (an example of this behavior is given in figure 5(H)), we estimated a noise corridor and used the model and the fit results to extrapolate the theoretical deformation curves into that noise corridor. The results are shown in figure 7. For the first extrapolation (red curves in figure 7) we assumed that S and ΔS take the same value as when the cells are probed using the heating protocol. Further, we assumed that the relative amount of temperature induced increase in force f is the same for HeLa and HL60 cells e.g.: This extrapolation led to deformation curves within the noise corridor for untreated cells probed using the short-pulse protocol, but not for the other conditions. We, thus, tested the assumption that undisturbed HeLa cells are stiffer than HL60 cells but also more affected by thermorheological softening effects. Since geometric softening induced by heating dominated the thermorheological softening in HL60 cells, we tested whether increased amounts of geometric softening could explain the behavior of HeLa cells. In particular, we employed larger decreases in forces f when the HeLa cells were probed by the non-heating protocol instead of the heating protocol compared to the decrease in force found in HL60 cells. We found that a four times larger decrease in force (purple curve in figure 7) was sufficient to make the deformation curves vanish inside the noise corridor, while an only two times larger decrease in force (orange curve in figure 7) was not. This suggests altogether that the denser and more tightly connected networks of cells that are cultured adherent to a substrate (HeLa) are also more susceptible to thermorheological softening effects possibly related to temperature-activated actomyosin activity. In physiotherapeutic terms, this suggests that for more severe cases of body stiffness and muscle spasms concomitant thermotherapy is required to induce even a short-term relief.
Besides our approach to thermorheological softening and long-lived structural remodeling there are two other approaches to which we can relate our analysis. First, the practice of thermorheology and time-temperature superposition, well known in polymer physics [45,46] which was more recently extended to living cells [7,8]. It was found that the creep responses measured at different temperatures can be superimposed either by time-temperature superposition alone, i.e. by a rescaling using a time-shift factor, or by an additional rescaling of the modulus. The cells were classified as thermorheologically simple or thermorheologically complex, accordingly. As this analysis allows cells to be classified with respect to their temperature-dependent mechanics without any theoretical modeling, one might be tempted to use this approach to identify long-lived structural remodeling. A possible classification could be that thermorheologically simple and complex behavior indicate different short-lived effects while even more complex behavior could indicate inelasticity and long-lived structural changes. However, this approach has several drawbacks. Most importantly, thermorheology rests on the assumption of linear viscoelastic material behavior. Thus, it is currently unable to distinguish inelastic effects from non-linear viscoelastic effects, which makes the detection of long-lived structural rearrangement effects problematic. In addition, there is no straightforward way to apply time-temperature superposition to our multiple-pulse protocols, which are vital for detecting inelastic effects or long-lived structural changes.
These problems notwithstanding, there is also an interesting similarity between the framework of thermorheology and our modeling approach: similar to time-temperature superposition, fluidization, i.e. a change in solidity S in our model also means a recalling of time scales: However, we do not rescale a single time scale but all the mode relaxation times with factors depending on their mode number, which is therefore more appropriately addressed as 'time-stretching'. There exists also a loose analogy between thermorheologically complex behavior and the combined occurrence of fluidization and geometric softening. Both lead to increased deformations that could be absorbed into some additional rescaling of the deformation-axis. Yet, geometric softening does so in a non-linear fashion and also affects the mode relaxation times, since already the underlying WLC relaxation times τ WLC n are explicitly force-dependent (5).
Another intuitive approach to non-linear softening effects in cells has been proposed in reference [47]. There, the total deformation was split into a reversible viscoelastic part obeying the Boltzmann-superposition-principle, which follows a power-law and an irregular (visco-)plastic part described by another power-law. This perspective, although it might be simple and powerful in certain situations (in particular in the context where it was first proposed for), is not refined enough for our present purpose. For many of our datasets, we can describe apparently plastic effects, i.e., the build-up of residual deformation from cycle to cycle, by the reversible but non-linear viscoelastic version of the GWLC introduced above.

Conclusions
With the aim of uncovering the cell-physiological effects of separate and conjoint mechanical and thermal treatment we scrutinized the nonlinear mechanical response of single cells to stretching and heating pulses applied selectively and reproducible by a versatile optical stretcher. Just as in ordinary physiotherapy a single intervention is seldom enough to introduce a long-lived change, we found it suitable for single cells to perform a series of stretching and heating pulses. In cell mechanics, long-lived mechanical effects of stretching and heating are hard to detect without a suitable rheological model. This is why we made a major effort to consistently fit a large data set from different cells and a range of rheological protocols with the viscoelastic GWLC model. Marked deviations from the fits only occurred for a small subset of the data, namely for the more vigorous heating and stretching protocols, and were identified as indicative of long-lived alterations of the mechanical state of the cell. The necessity of using intense thermorheological protocols to induce long-live structural changes or even to deform the cells at all illustrates the robustness of cells as a mechanical system. Our results also provide insight into the ongoing debate concerning the effect of temperature on cell mechanics: while many studies report that cells become softer or fluidize with increasing temperature [1-3, 7, 8, 39, 48], there is also evidence [14] for the opposite scenario, namely that cells become stiffer and more solid-like with increasing temperature. One may argue that results from different studies, where different cell types are probed and analyzed by different methods and different measures for stiffness are used, are not easily comparable. Our discussion refines this argument, as we can disentangle different possible thermorheological effects and illuminates different underlying mechanisms using the framework of the GWLC. In this context, softening corresponds to node-dilution in the relevant cytoskeletal network while fluidization corresponds to reduction of caging effects. Since we find that softening can be accompanied by more elastic (decreased ΔS) or less elastic (increased ΔS) behavior, we conclude that, depending on the experimental protocol and investigated system, thermomechanical effects might lead to apparent softening or stiffening, which both can be accompanied by more solid-like or more fluid-like behavior in the sense of a decreased or increased loss angle measured in a linear rheology experiment.
Based on our mathematical model, we were able to identify likely candidates for the structural components of the cytoskeleton that are predominantly responsible for the measured cell response to stretching and heating. These expectations could in the future be combined with other (e.g. microscopic, biochemical, . . . ) evidence in order to pinpoint the specific molecular origins of the demonstrated antagonistic effects of heating and the revealed dissimilar cell mechanical response patterns, and possibly even their mechanistic correlation with macroscopic treatments and symptoms.