Time dependent stress relaxation and recovery in mechanically strained 3D microtissues

Characterizing the time-dependent mechanical properties of cells is not only necessary to determine how they deform but also to understand how external forces trigger biochemical-signaling cascades to govern their behavior. At present, mechanical properties are largely assessed by applying local shear or compressive forces on single cells grown in isolation on non-physiological 2D surfaces. In comparison, we developed the microfabricated vacuum actuated stretcher to measure tensile loading of 3D multicellular “microtissue” cultures. Using this approach, we here assessed the time-dependent stress relaxation and recovery responses of microtissues and quantified the spatial viscoelastic deformation following step length changes. Unlike previous results, stress relaxation and recovery in microtissues measured over a range of step amplitudes and pharmacological treatments followed an augmented stretched exponential behavior describing a broad distribution of inter-related timescales. Furthermore, despite the variety of experimental conditions, all responses led to a single linear relationship between the residual elastic stress and the degree of stress relaxation, suggesting that these mechanical properties are coupled through interactions between structural elements and the association of cells with their matrix. Finally, although stress relaxation could be quantitatively and spatially linked to recovery, they differed greatly in their dynamics; while stress recovery acted as a linear process, relaxation time constants changed with an inverse power law with the step size. This assessment of microtissues offers insights into how the collective behavior of cells in a 3D collagen matrix generates the dynamic mechanical properties of tissues, which is necessary to understand how cells deform and sense mechanical forces in vivo.


INTRODUCTION
The requirements of crawling, dividing, and contracting require the cell's cytoskeletal network of structural and motor proteins to be tremendously dynamic. This behavior is unique for other soft materials and gives cells and tissues their distinct elastic and dissipative properties. Defining these properties is necessary for understanding not only how cells deform but also how cells sense and transduce external mechanical forces into biochemical signals that direct their behavior in vivo. In that regard, when cells are stretched or come into contact with a stiffer matrix, there are time scale-dependent conformational changes to adhesion and cytoskeletal protein networks, which, in turn, alter ligand-receptor binding affinities to trigger biochemical signaling cascades. 1,2 Through this regulation of biochemical signaling, mechanical forces have been linked to normal development and function, as well as disease progression, including bone, muscle, heart, and lung disorders and cancer. 3,4 In regard to defining the mechanical behavior of living matter, it has long been recognized that cells and tissues exhibit both solid-like elastic and fluid-like viscous properties. 5,6 Traditionally, this behavior, called viscoelasticity, has been described using a network of elastic springs and viscous dashpots. In particular, when these elements are connected in series, solving the constitutive equations for a step change in the length gives an exponential decay in stress with a characteristic time constant(s), which depends upon the elastic modulus of the spring(s) and viscosity of the dashpot(s). 7 With these spring-dashpot models in mind, many studies set out to characterize the viscoelastic behavior of single cells and to link them to specific occurring processes. Although early experimental data could be fit with a single time constant, [8][9][10][11] as the resolution of techniques improved, a power law behavior emerged. 12 In particular, the frequency, creep, and stress relaxation responses of isolated cells have all been shown to be accurately captured by a single power exponent describing a continuous, featureless distribution of timescales. [12][13][14][15][16] Further universal observations of power law rheology across different cell types and techniques and following a range of cytoskeletal drugs has since given traction to the hypothesis that cells belong to a class materials called soft glasses. 12,13 One notable exception to this featureless rheological behavior seems to occur at short timescales (<1 s) and under large volumetric deformations Under these conditions, a characteristic behavior has been reported arising from poroelastic effects caused by the redistribution of cytosolic fluid. 17 In the case of tissue-level mechanics, while a featureless relaxation behavior has also been reported, 6,18 the field has not reached a consensus on the use of power laws in describing their viscoelastic response. Rather, spring-dashpot models with characteristic timescales remain prominently reported in soft tissue mechanics. 19,20 For instance, the viscoelastic behavior of muscle tissue, particularly when the mechanical response is dominated by actin and myosin kinetics, has been shown to deviate from a power-law behavior and instead followed a broad distribution of timescales around a characteristic time constant set by acto-myosin activity. 21 Furthermore, growing cells on a 2D substrate, as largely required for assessing individual cell mechanics, forces an un-natural apical-basal polarity of adhesion complexes. This, in turn, is known to cause vast differences in the distribution and structure of the cytoskeleton. 22 Although it remains unclear, it is not unreasonable to suspect these fundamental changes to the cytoskeleton, caused by the dimensionality of the cell's environment, may alter the mechanical behavior of cells in 3D tissues compared to when studied on 2D substrates. Finally, characteristic detachment rates of cellular adhesions to matrix proteins may further differentiate tissue-level viscoelastic behavior from the power-law rheological model widely seen in isolated cells. 23 Therefore, although our knowledge of the mechanical behavior of isolated cells has greatly advanced, we lack a complete understanding of how a heterogeneous population of cells within a 3D extracellular matrix establishes the time scale varying and non-linear viscoelastic properties of tissues. To answer this research question, 3D cell culture methods that enable the assessment of tensile forces have become of keen interest in the field of cell mechanics. 24 In that regard, both the frequency and the stress relaxation response of cells within bulk reconstituted collagen gels have previously displayed a characteristic timescale behavior following standard linear spring-dashpot models. [25][26][27] Due to their centimeter-scale, however, these cultures tend to have poor cellular organization and low cell density, suffer from slow experimental throughput, are hard to image, and possess a high diffusive barrier for nutrients. These limitations of bulk 3D cell cultures can be largely overcome by shrinking the cell culture size by adopting a lab-on-chip approach, as in the microtissue model. 28 In particular, the high-throughput array of sub-millimeter-scale structures that form around pairs of flexible cantilevers in the microtissue model possesses comparable cell alignment and density to human dense connective and muscle tissues, as well as enables high-resolution live-cell imaging and assessment of acto-myosin dynamics. [28][29][30][31] Although quasi-static mechanics (i.e., contractility and stiffness) of microtissues have already received much attention, 28,29,32,33 the time dependencies in these 3D cell cultures remain unclear.
Accordingly, in this article, we investigated the viscoelastic stress relaxation of microtissue cultures using a microfabricated device that we have developed called the Microtissue Vacuum Actuated Stretcher (MVAS). 30,31 The MVAS allows large deformations, simultaneous measurements of tension, and live imaging of microtissue cultures. Importantly, the description of the dynamic behavior of 3D microtissue cultures that follows in this article qualitatively differs from previous measurements on isolated cells in 2D culture and, thus, raises important questions about our understanding of how an aggregate population of cells produces the time-dependent and nonlinear mechanical properties of tissues.

Microtissues are viscoelastic
After three days of static culture in the MVAS device, fibroblasts cells had compacted the collagen matrix to form dense 3D microtissues suspended between the cantilevers (Fig. 1). Based upon the deflection of the force-sensing cantilever, microtissues had developed a resting tension of 8.7 6 0.4 lN (N ¼ 79). The cells were mostly aligned with the direction of tension development shown by both the longitudinally oriented actin cytoskeleton and nuclei [ Fig. 1(b)]. As previously reported, 28,29 these observations indicated that microtissues formed under tensile constraints, as in the MVAS device, are tightly compacted highly organized structures and, in these respects, broadly resemble dense connective and muscle tissues.
In addition to static tensile measurements, the MVAS enables the assessment of mechanical behavior in response to changes in strain through a vacuum-driven planar stretch [ Fig. 1(a)]. In such conditions, following a step change in the length, microtissue tension rose sharply and then relaxed to a new equilibrium point as commonly seen in viscoelastic solids [ Fig. 1(c)]. Upon returning the microtissues to their initial length, the tension dropped past and then slowly recovered to the original resting value (P > 0.05, repeated measures t-tests). The stress relaxation and recovery responses were well conserved among microtissues [average standard deviation (SD): 0.3 lN] and were highly repeatable in a given microtissue as subsequent loading cycles were largely superimposable (supplementary material II).
Stress relaxation and recovery responses together indicated that microtissues modulate their tension in response to step changes in the length toward their resting values. Importantly, however, microtissues did not completely relax upon a step strain, but rather seemed to be capable of maintaining residual stresses. Consequently, at each strain amplitude, microtissues reached a unique tension equilibrium that was independent of the loading history.
In order to further characterize the tension (T) relaxation and recovery responses, we searched for an appropriate mathematical model to describe these viscoelastic behaviors [Figs. 2(a) and 2(b), respectively]. Classically, engineers have modeled viscoelastic effects with a series of springs, described by their stiffness constants (k) and dashpots, described by their viscosity. 7 In particular, the standard linear solid (SLS) model consists of a spring and a dashpot in series (i.e., a Maxwell body) in parallel to another spring. Following a step change in the length, this model predicts that tension will exponentially decay with a time constant, s, toward a residual value [T SLS t ð Þ ¼ e k 1 e Àt=s þ k 2 À Á ]. The amplitude of this exponential function is given by the spring constant k 1 , while the residual stress is described by the spring constant k 2 . As such an SLS model has been previously applied to early experiments on isolated cells, [9][10][11] 3D cell cultures, 25 and ex vivo tissue strips, 7 we begun by considering this model for describing microtissue mechanics. In that regard, Figs. 2(a) and 2(b) show that a SLS model could describe microtissue mechanics at short time-scales but failed to capture the long time tail in both relaxation and recovery responses, respectively.
In contrast to an exponential model, slow relaxation and recovery at long timescales are exemplified by a power law rheological model. This second model predicts that stress relaxation/recovery should follow a power law with a dimensionless constant (a), which may describe a specific continuous distribution of discrete viscoelastic behaviors [T PL ðtÞ ¼ e k 1 t a þ k 2 ð Þ ]. Since a long time tail behavior has been widely reported in single cells [12][13][14][15][16] and tissues, 6,18 we next considered a power law model for describing microtissue mechanics. We found that such a power law could capture the slow relaxation behavior of microtissues at long timescales but predicted faster than experimental dynamics at shorter times.
From considering these first two models, it was clear that in order to capture both the exponential-like stress decay at shorttime scales and the slow power-law relaxation at long time-scales, we needed a model that first appeared as an exponential but crossed over to a more slowly relaxing tail behavior at longer elapsed times. This behavior is characteristic of the so-called stretched exponential function [Eq. (1)]. In this model, one may FIG. 1. Microtissue tension dynamically relaxed and recovered with changes in the length. Microtissues were grown in our MVAS device (a). In the MVAS, changes in the microtissue length are driven by applying a regulated vacuum pressure to a chamber that borders one side of each microtissue well. The cantilever closest to the vacuum chamber (shown here on the right) is actuated to stretch the microtissue in plane, while changes in tension can then be measured by tracking the deflection of the opposing cantilever (on the left). Microtissues are organized 3D cell cultures freely suspended between the cantilevers. Max projections of confocal stacks, orthogonal views and high magnification images, are shown in (b). Both the actin cytoskeleton (green) and the nuclei (blue) are predominately aligned between the cantilevers. Scale bars in (a) and (b) represent 200 and 100 lm, respectively. Following a step in strain (inset), microtissue tension (N ¼ 79) sharply increased and then relaxed to a new set point as expected in a viscoelastic solid (c). Then, upon returning the microtissue to its original length, the tension recovered.
interpret the power law constant, b, as a dimensionless descriptor of a specific distribution of discrete timescales, which broadens from a single timescale behavior as the value of b decreases from 1 to 0. In other words, a stretched exponential function mimics a sum of Maxwell bodies with a specific distribution of time constants. However, in the absence of a physical model for this distribution, this method is an ineffective use of variables compared to a stretched exponential. A further discussion of stretch exponential FIG. 2. Stress relaxation and recovery in microtissues followed stretched exponential trajectories. Viscoelastic modeling of stress relaxation and recovery responses is shown in (a) and (b), respectively. In both panels, representative data were normalized to the peak tension. Stretched exponentials capture relaxation and recovery behaviors over three decades of time. To assess the dynamics of the relaxation and recovery, responses were normalized to their amplitudes using tension measurements immediately following step length changes, T 0 , and after 160 s, T Ã . By normalizing responses in this manner, it is clear that relaxation occurred much quicker than stress recovery (c). Yet, stress relaxation and recovery appeared to share the same spatial locations within a given microtissue in terms of both viscoelastic deformation in the longitudinal (e xx ) and transverse (e yy ) directions immediately following the change in the length (d). The scale bar represents 100 lm.
APL Bioengineering ARTICLE scitation.org/journal/apb model parameters and simulated stress relaxation responses is presented in supplementary material III, Although we have not come across previous reports that have modeled the viscoelastic behavior of cells or tissues with stretched exponentials, they have been widely used to describe relaxation in glassy, disordered systems that display a broad distribution of timescales. 35,36 In our microtissue experiments, however, a stretched exponential captured both stress relaxation and recovery data over three decades of time (R 2 > 0.99) with average fitting constants summarized in Table I. Importantly, the creep that inherently accompanied our method of measuring microtissue tension relaxation was insufficient to significantly change either the model or the fitting constants (supplementary material IV).
The absolute amplitudes of the stretch exponential stress relaxation and recovery responses, given in Table I, agreed well with each other (0.55 6 0.03 vs 0.50 6 0.02 lN/%); P > 0.05, repeated measure t-tests). Moreover, the spatial distributions of strain immediately following the step change in the length, describing the locations of stress relaxation and recovery, were well correlated [ Fig. 2(d)]. Because the regions that underwent stress relaxation were spatially and quantitatively linked to stress recovery, the viscoelastic response of microtissues likely occurred through a broad distribution of reversible viscoelastic cell deformations.
The spring constant k 2 in our mathematical model [Eq. (1)] characterizes the recoverable, residual stress that microtissues elastically store, as described above. On average, microtissues had a residual stiffness of 0.42 6 0.02 lN/%. The source of this elasticity is unclear, but it significantly contributes to the mechanics of microtissues. It may arise from the contribution of the matrix or perhaps is indicative of the existence of a stress threshold up to which cytoskeletal elements can behave elastically but beyond which they yield. In recovery curves, k 2 is approximately equal to zero, suggesting that microtissues fully recovered to their original prestressed state upon being returned to their initial length.
The dynamics of the relaxation and recovery responses are shown in Fig. 2(c). In this figure, microtissue tension, T, has been normalized to the response amplitude measured across the initial 160 s. Normalizing the curves in this manner revealed that relaxation occurred much quicker than stress recovery (s ¼ 14 6 1 vs 30 6 2 s, repeated t-test, P < 0.001). Furthermore, relaxation responses had a significantly, albeit slightly, higher power law constant (b ¼ 0.69 6 0.02 vs 0.62 6 0.02, repeated t-test, P < 0.01). Therefore, although stress relaxation and recovery are spatially and quantitatively linked, the dynamics of these responses varied considerably; recovery occurred more slowly and over a slightly more broadly distributed set of time scales.
Taken together, in this section, we have described the viscoelasticity of microtissues and found that this behavior can be captured through a relatively efficient stretched exponential model with a residual elasticity. In the remainder of this article, we aimed to gain further insight into how the mechanics of microtissues are controlled by underlying molecular and structural mechanisms in the cells by first examining how this behavior is affected by strain amplitude, and second, how microtissue viscoelasticity changes in a range of pharmacological treatments targeting specific cytoskeletal elements.

Microtissue viscoelasticity is nonlinear
To assess whether there were nonlinearities in the viscoelastic behavior of microtissues, stress relaxation and recovery measurements were performed at two different step strain amplitudes on the same microtissues. Shown in Fig. 3(a) are average response curves for both 3.3 6 0.2% and 10.9 6 0.8% step strains (N ¼ 17). The changes in fitting constants for stress relaxation and recovery between these strain amplitudes are summarized in Tables II and III, respectively. A larger step strain amplitude decreased both k 1 and k 2 spring constants for relaxation responses by 26 6 4 and 11 6 5%, respectively (repeated measures T-tests; for k 1 : P< 0.001 and for k 2 : P < 0.05). Likewise, the absolute k 1 spring constant fitted to recovery responses decreased by 28 6 5% (repeated measures T-test; P < 0.01). Altogether, these responses were indicative of a strain amplitude softening behavior that followed from abrupt loading.
The relaxation responses were normalized as shown in Fig. 3(b) to show how the dynamics varied with the strain amplitude. It was found that a larger step strain decreased s by 56 6 4% and b by 17 6 4% (repeated t-tests; for s: P < 0.001, for b: P < 0.01). Therefore, by increasing the step size, microtissues relaxed more quickly and over a more broadly distributed set of time scales.
In contrast to stress relaxation responses, the time and power law constants of stress recovery were step strain amplitude invariant [ Fig. 3(c)] (repeated t-tests; both P > 0.05). Moreover, they did not depend upon the strain at which the microtissue recovered (supplementary material V). Therefore, whereas the rate of stress relaxation in microtissues was nonlinear with strain amplitude, stress recovery was always seen to be a linear process.
To further define the relationships between stretch exponential fitting constants and the strain amplitude, 160 responses with step strains ranging from 1.75 6 0.06% to 14.3 6 0.2% were discretized into strain bins and analyzed with regression models. In accordance with repeated measures observations, Fig. 3(d) shows that the relaxation time constant decreased with the strain amplitude and surprisingly followed a power law dependence (non-linear regression: R 2 ¼ 0.96 P < 0.001). Further, in agreement, the recovery time constant was seen to be strain-invariant (linear regression: P > 0.05). The remaining fitting constants to the relaxation and recovery responses are summarized in supplementary material VI and VII, respectively.
Thus far in this section, we have shown that microtissue mechanics varied considerably with the imposed strain amplitude. To link these mechanical behaviors to changes in the local spatial viscoelastic deformations, the local strain distribution was characterized over time immediately starting after the change in the length. Figures 4(a) and 4(b) show the integrated raw and normalized local viscoelastic deformations in the transverse direction in representative microtissues over 100 s of stress relaxation and recovery, respectively. As with the fitting  constants, the spatial distribution of local viscoelastic deformations that occur following a step change in the length varied considerably with strain. In that regard, viscoelastic deformations during stress relaxation and recovery in the transverse direction increased with the step size and became transversely focused to the center of the tissue but longitudinally diffuse, spanning the entire tissue length.
To compute a mean metric of the integrated relaxation response through time, transverse strain fields were averaged across the region of interest [ Fig. 4(c)]. Following a small step strain (1.7 6 0.2%), microtissues showed comparably little relaxation in the transverse direction. In comparison with an intermediate step strain (3.4 6 0.2%), microtissues relaxed to reach a new FIG. 3. The rate of stress relaxation was nonlinear, while recovery was strain invariant. To assess viscoelastic non-linearity, the step size was varied (inset). Stress relaxations for a 3.3 and 10.9% step are shown in (a) (N ¼ 17). After normalizing for response amplitude over the initial 80 s, (b) shows that the rate of stress relaxation increased with the strain amplitude. The rate of recovery was, however, invariant upon the step size (c). To further understand the nonlinearities in the viscoelastic behavior, 160 step strain experiments were completed across a range of strain amplitudes. Panel (d) shows that, as the strain amplitude increased, the time constant for stress relaxation decreased with an inverse power law with a near unity exponent, while in contrast, the rate of recovery was again invariant upon the step size.
À266 4 (ÃÃÃ) À11 6 5 (Ã) À56 6 4 (ÃÃÃ) À17 6 4 (ÃÃ)   Fig. 4(c)], the mean viscoelastic deformation during recovery showed slower dynamics as expected   Fig. 4(d)]. Yet, following small and intermediate step strains, both tensile relaxation and compressive recovery occurred with similar spatial distributions and amplitudes, indicating that the gross structural viscoelastic deformation that occurred after a step strain was a reversible process. In contrast, however, transverse viscoelastic deformation recovery following a large step strain was quantifiably greater than the final equilibrium in relaxation curves because of the delayed contraction.
In comparison with the viscoelastic deformation in the transverse direction, the deformation in the longitudinal direction was less extensive being constrained by the cantilevers. Longitudinal viscoelastic deformation was also highly tissue-specific and mainly located near the fronts of the cantilevers rather than the body of the tissue. Since final relaxation and recovery longitudinal strain fields shared similar amplitudes and microtissue tension fully recovered, this behavior likely did not indicate tissue detachment from the cantilevers (supplementary material VIII). Rather, the longitudinal viscoelastic deformation response more likely reflected stress concentrations formed around the front of the cantilevers created through a reduced fraction of longitudinally oriented actin filaments [ Fig. 1(b)].

Pharmacological behaviors
To assess the role of several key cellular proteins in governing microtissue viscoelasticity, stress relaxation and recovery were measured following various pharmacological treatments meant to either disrupt the cytoskeleton or alter actomyosin activity. Normalized responses and fitting constants are shown in Fig. 5(a) and Table IV, respectively. Each fitting constant was normalized and compared with paired t-tests to their own pre-treatment control. DMSO was used as a loading control and did not affect either relaxation or recovery curves. Furthermore, tension measurements taken prior to and following step strain experiments did not vary with any treatment, indicating that stress relaxations were fully recoverable under each condition and that incubation periods were sufficient to reach an equilibrium state.
First, consistent with our understanding that microtubules act as compressive struts that oppose actomyosin activity, 37 microtubule depolymerization with nocodazole increased the resting tension in microtissues by 3.7 6 0.4 lN (P < 0.001), as well as spring constants for stress relaxation (k 1 and k 2 both P < 0.001) and recovery (k 1 P < 0.001). While the time or power law constants for stress relaxation (s and b both P > 0.05) were not affected, microtubule depolymerization increased the stress recovery time constant by 87 6 25% and decreased the power law constant by À7 6 4% (P < 0.01 and P < 0.05, respectively). This suggests that although the relaxation rate is not affected, the recovery rate and, more modestly, the broadness of the distribution of internal timescales in the response may be altered through the regulatory action of microtubules on acto-myosin activity.
Second, as expected, myosin inhibition with blebbistatin decreased resting microtissue tension by À2.8 6 0.5 lN (P < 0.001) and reduced relaxation and recovery spring constants (various P values). Similar to microtubule depolymerization, blebbistatin treatment did not affect either the relaxation time or power law constants (s and b P > 0.05). Blebbistatin treatment did decrease the recovery time constant by À27 6 8% (P < 0.01), which indicates that myosin motor activity, also in part, determines the recovery speed. However, unlike microtubule depolymerization, myosin inhibition did not affect the power law constant (P < 0.0) and, thus, does not determine the broadness of the distribution of internal time scales.
Finally, actin depolymerization (cytochalasin D) and decellularization (triton-X) had similar effects on the mechanical behavior of microtissues. Both treatments decreased resting microtissue tension (À5.0 6 0.4 lN and À8.8 6 0.8 lN, respectively; both P < 0.001) and relaxation and recovery spring constants (various P values). Neither treatment displayed observable stress relaxation or recovery responses, and for this reason, it was not possible to accurately fit time and power law constants. Nevertheless, these responses revealed that deformation in actin microfilaments is chiefly reasonable for the viscoelastic behavior of microtissues.
To examine how microtubules, actin microfilaments, myosin motor activity, and the extracellular matrix contributed to the local viscoelastic deformation of microtissue during relaxation and recovery, strain fields were assessed following a small step strain (1.7 6 0.2%) with the pharmacological treatments. Microtubule depolymerization and myosin inhibition behaved similar to control microtissues in their transverse (Figs. 5(b)-5(d)] and longitudinal (supplementary material IX) deformations. In contrast, actin depolymerization and decellularization lead to viscoelastic deformation behaviors that were more consistent with a larger step strain in terms of the spatial distribution, amplitude, and dynamics of the strain fields, apart from lacking an intermediary contraction response. This finding seems to implicate the depolymerization of actin and the contribution from matrix proteins in the gross deformation behavior and nonlinear dynamics observed with large strains.

DISCUSSION
In this article, we assessed the dynamic mechanical behavior of physiologically relevant 3D microtissue cell cultures. As expected in living active matter, microtissues displayed timescale-varying mechanics. In response to a step increase in the length, the tension in microtissue cultures quickly rose and then relaxed to a new equilibrium stress, resembling a viscoelastic solid. Then, when microtissues were returned to their initial lengths by a recovery strain, their tension quickly dropped below initial measurements and then slowly, yet fully, recovered. As these responses differed from the widely observed weak power laws that describe the relaxation behavior of cells in 2D culture, [12][13][14][15][16] our findings in 3D microtissues may reflect additional considerations necessary to completely understand the viscoelastic response of living tissues. In that regard, while our field has largely been focused on characterizing the mechanical response of individual cells, tissue mechanics is perhaps instead generated by an aggregate of active cytoskeletal dynamics in a heterogeneous population of interconnected cells and through interacting with a soft predominantly elastic 3D matrix.
In particular, unlike previous work in 2D cell culture, we found that tensile stress relaxation and recovery in microtissues followed stretched exponential behavior. To the best of our knowledge, stretched exponentials have not been previously used to describe the viscoelastic behavior of cells or cell cultures. Still, however, our findings are in close agreement with previous reports from bulk fibroblast populated collagen gels. 25,27 In that regard, we reported characteristic time constants of similar magnitudes with time dependencies that are more broadly distributed than what can be explained by a single exponential function (i.e., a SLS model). On the contrary, there does exist a strong precedent for the use of stretched exponentials to describe relaxation in disordered condensed matter systems. In fact, both polymer networks 38 and glasses near their transition temperature 35,36 relax according to stretch exponentials. Notably, such a glassy polymer response has already been linked to the mechanical behavior of the cytoskeleton. 39 Furthermore, our average values of b (0.64 6 0.01 and 0.61 6 0.01 for relaxation and recovery, respectively) are surprisingly close to a universal experimental value of 0.6 in glasses dominated by Brownian motion, 40 which is theoretically justified by the so-called "trapping model." 41 Mathematically, this model is beyond the scope of this work, but it is relatively easy to conceptualize. In a trapping model, the material consists of randomly distributed traps. During relaxation, these traps capture molecules that diffuse through the material, and as the traps are filled, the remaining molecules must travel longer to find unoccupied traps. Thereby, the relaxation rate falls with a power law given by the dimensionality of the model: in 3D environments, b ¼ 0. 6. Although a trapping model can be a great tool to conceptualize stretched exponentials, the question remains whether we can draw out any physiological meaning or link any mechanism(s) to account for the stretched exponential relaxation and recovery in microtissue behavior. First, similar to power laws, stretched exponentials can be expanded into a superposition of exponential functions with a nontrivial distribution of relaxation times. 42 In this respect, the power law constant describes the broadness of the distribution. It should not be overly surprising that relaxation processes exist across a broad spectrum of timescales in something as complex and heterogeneous as the network of cytoskeletal proteins inside single cells and, in a more collective sense, when considering an aggregate of cells and extracellular matrix. Thus, our fitting constants may be interpreted as a general description of microtissue behavior reflecting the broad viscoelastic heterogeneity in the cytoskeleton and among cells rather than capturing a specific process. What may be more surprising was that the same model, describing an inter-related set of exponential functions, could be applied over a range of pharmacological treatments and strain amplitudes with only small changes to the power law constant, while changes in the other fitting constants were large. At the moment, an explanation for this behavior remains unclear. Perhaps future work in generating a more mechanistic modeling approach may provide further insight into the stretch exponential behavior of microtissues and their power law constant.
In addition to displaying stretch exponential dynamics, we found that the relaxation spring constants (k 1 and k 2 ) of microtissues remained tightly coupled throughout an assortment of pharmacological treatments targeting specific proteins as well as in measurements at different step strain amplitudes. In that regard, all treatments followed a single linear relationship [ Fig. 6(a)]. Although assigning specific meaning to these spring constants would be speculative, our results indicated they both largely depended upon the actin cytoskeleton and, to a lesser extent, myosin activity, the presence of microtubules, and the elasticity of the matrix. Previously, this coupling behavior between energy loss and residual elasticity in tissues has been interpreted as evidence suggesting that residual elastic and dissipative stresses are borne from the same origin(s). 43,44 Even so, it is surprising that targeting different cytoskeletal elements can produce responses that follow a single universal relationship, unless an understanding of mechanisms at a hierarchy beyond the single protein level is required. In other words, rather than following responses in single molecules, as our field has been studying in a traditional reductionist approach, the viscoelasticity of cells and tissues is perhaps more influenced by the association and interaction of proteins with each other. 45 For example, through bundling collagen filaments, cells may control the residual elasticity of the matrix, and thereby, any change in acto-myosin dynamics will also be seen as a change in the stiffness of the matrix. Further evidence in support of this view has also been reported through various universalities in single cell mechanics. 13,37,46,47 The degree of stress relaxation was also strongly coupled to the stress regained during recovery [ Fig. 6(b)]. Furthermore, for the most part, local microtissue deformation during relaxation and recovery occurred with similar spatial distributions and comparable amplitudes. Altogether, these results suggest that stress relaxation and recovery depended upon reversible remodeling in the same structural elements within the microtissue. That being said, a notable exception was observed in the transverse deformation response to a large amplitude (10%) strain, where a delayed contraction followed stress relaxation. This contractile behavior following a large step strain is consistent with stretch activated contraction responses previously reported in isolated fibroblast cells 48 and 3D cultures, 49 and it seems to act separately from stress relaxation and recovery, causing permanent structural remodeling.
In contrast to the amplitudes, the dynamics of the relaxation and recovery responses differed greatly. Relaxation rates were nonlinear with step strain amplitude, decreasing with what appeared to be a power law, but were relatively invariant with microtubule depolymerization or myosin inhibition. In contrast, recovery rates were invariant with step strain amplitude but increased with microtubule TABLE IV. The pharmacologically induced changes in stress relaxation and recovery fitting constants. Each treatment was compared with its pretreatment control using paired t-tests ( Ã P < 0.05, ÃÃ P < 0.01, and ÃÃÃ P < 0.001).

Relaxation
Recovery À8.8 6 0.8 (ÃÃÃ) À86 6 6 (ÃÃÃ) À65 6 7 (ÃÃ) … … À109 6 10 (ÃÃÃ) … … À0.1 6 0.2 depolymerization and decreased with myosin inhibition. This suggests that myosin contributed only as a passive cross-linker during relaxation, while active myosin crossbridge cycling was partially responsible for stress recovery. In either case, the actin cytoskeleton was the principal source of the viscoelastic behavior of microtissues since depolymerization of actin filaments had the biggest effect on step strain responses outside of decellularization. There also appeared to be some passive component to the recovery as there was a small degree of viscoelastic deformation upon returning samples with actindepolymerized cells and as well as in decellularized microtissues to their initial lengths. Interestingly, this deformation response was consistent with larger step strains, suggesting that the loss of the actin cytoskeleton and how the cells associate with the matrix may be in part be responsible for the non-linear relaxation rates. It is not uncommon for relaxation in soft biological tissues to be nonlinear. For example, nonlinearities have been previously reported in lungs, 50 heart valve, 51 ligament, 52 and muscle. 53 Recently, nonlinear relaxation rates in cells were argued to arise from a poroelastic effect occurring through the redistribution of the viscous cytosolic fluid between local regions within the network of cytoskeletal filaments. 17 It is, however, debatable that the same poroelasticity effect is the main determinant in our measurements, as our time constants differed more than a decade from those previously reported. Moreover, poroelasticity does not offer an explanation as to why the recovery rates were linear with strain amplitude and much slower than relaxation rates.
An alternative mechanism conceives stress relaxation in cells as a friction force produced by cytoskeletal filaments sliding smoothly past one another. 54 The main support for this hypothesis is that local shear measurements on cells in 2D culture follow the structural damping law; the ratio of dissipative to elastic stresses is invariant of the timescale at which they are assessed. 12,13 This hypothesis, however, does not offer any explanation as to the nonlinearity observed in relaxation rates and the imbalance with recovery dynamics. Furthermore, it is not compatible with either standard linear solid (SLS) viscoelasticity or a stretched exponential version of the SLS model utilized in this paper, as both models produce characteristic rate dependencies that are inconsistent with structural damping. Interestingly, in contrast to local shear measurements, characteristic rate dependencies similar to our observations in microtissues have been previously reported in individual cells 11 and 3D cultures 25 when probed in tension at a length scale of the entire cell/culture. It would still be of interest, however, in future work to check if the frequency response of microtissues is consistent with their stretched exponential relaxation in the time domain.
A last possible mechanism to explain the viscoelastic response in microtissues is that stress relaxation and recovery reflect the rupturing and reforming of bonds within the cytoskeleton (i.e., within actin filaments and between actin and its crosslinkers) and between the cells and the extracellular matrix. In support of this hypothesis, there is a strong precedent that mechanical stretch depolymerizes actin filaments 46,55-57 and perturbs myosin binding. [58][59][60][61][62] Indeed, we have shown here that microtissues strain-soften in an amplitude-dependent manner. Further, in regard to this hypothesis, the imbalance in microtissue relaxation and recovery rates can simply be explained by differences in bond destruction and formation rates; it is not unreasonable to suspect that it takes less time to pull bonds apart from that for proteins in a cell to diffuse through Brownian motion, correctly reorient themselves, and then finally re-form chemical bonds, especially if FIG. 6. Underlying relationships in microtissue mechanics. The k 1 and k 2 stress relaxation spring constants were linearly coupled throughout an assortment of pharmacological treatments and step lengths (a) (linear regression; P < 0.001; R 2 ¼ 0.94). Likewise, the amplitude of stress relaxation was coupled to the amplitude of stress recovery (b) (linear regression; P < 0.001; R 2 ¼ 0.92).
additional enzymes (i.e., profilin) and molecules (i.e., ATP) are needed. Furthermore, nonlinearities in relaxation behavior have previously been captured by models in which stress is relieved through sequential rupturing of Maxwell bodies, where each micro-yield event passes the stress on to other regions of the tissue. 63,64 The results from these models give an exponential decay with a long-tail power law. However, as chemical bonds have a finite yield strength, in reality, curves may flatten, behaving as stretched exponentials. Admittedly, however, this model has yet been used to explain nonlinearities beyond quasi-linearity, such as the inverse power law between the relaxation rate and strain amplitude that we observed.

CONCLUSIONS
The viscoelastic behavior of microtissues likely does not come down to a single physiological mechanism, but rather it is an amalgamation of many physical remodeling events occurring at interrelated time and length scales. Accordingly, the measured stress relaxation and recovery of microtissues followed generalized stretch exponential behavior, differing from the power law rheology seen in isolated cells. We further found that stress relaxation rates were nonlinear and largely imbalanced with a linear recovery response. With that said, the contributions of specific cytoskeletal elements (actin, myosin, and microtubules) to microtissue mechanics did not qualitatively differ from our previous understanding of their roles in the mechanical behavior of cells grown on 2D substrates. In particular, actin was found to be predominantly responsible for the viscoelasticity of microtissues. Taken together, however, our results collapsed onto a single relationship between residual elastic stresses (k 2 ) and stress relaxation amplitudes (k 1 ), indicating that the mechanics of microtissues may follow a universal behavior set not by the role of individual proteins but rather by the physical architecture of interconnected cells and their surrounding matrix as a whole. To conclude, the assessment of dynamic mechanics of microtissues has yielded further insights into how tissues gain their mechanical properties from the apparent roles of cells and from the distribution and interaction of their cytoskeletal proteins. This knowledge is a crucial consideration for understanding how physical forces are sensed and regulate cell behavior in health and disease.

MVAS device
MVAS devices were fabricated out of polydimethylsiloxane (PDMS) using mold replication and plasma bonding steps outlined previously. 30,31 Briefly, devices consist of 60 microtissue wells, each containing two cantilevers spaced apart by 500 lm and bordered on one side by a vacuum chamber [ Fig. 1(a)]. Devices consist of three layers replicated from masters created by standard photolithographic techniques. First, the top device layer forms the open-top microtissue wells and has enclosed vacuum chambers. Second, in the middle is a thin flexible membrane that is fabricated with the cantilevers. Finally, the bottom layer has enclosed vacuum chambers that match the geometry of the top layer and empty chambers that allow for the middle membrane to stretch in plane. When a vacuum is applied via an external electronic vacuum regulator (SMC ITV0010) controlled through Labview software, the cantilever closest to the vacuum chamber moves to stretch the microtissue. The opposing cantilever is used as a passive force sensor by optically tracking its deflection and utilizing its known spring constant (k cantilever ¼ 0.83 N/m). The dimensions of the force sensing cantilever were empirically chosen to be stiff enough to reduce the amount of creep that inherently occurs during stress relaxation because of our method of force measurement, yet also soft enough in order to have a sufficient signal to noise ratio in tension measurements.

Microtissue fabrication
Microtissues consisting of 3T3 fibroblasts in a 3D collagen matrix were cultured as previously described. 28,30 Briefly, the MVAS was sterilized with 70% ethanol and treated with 0.2% Pluronic F-127 (P6866, Invitrogen) for two minutes to reduce cell adhesion. 250 000 cells were resuspended in 1.5 mg/ml rat tail collagen type I (354 249, Corning) solution containing 1Â DMEM (SH30003.02, Hyclone), 44 mM NaHCO 3 , 15 mM d-ribose (R9629, Sigma Aldrich), 1% FBS, and 1 M NaOH to achieve a final pH of 7.0-7.4. The cell-collagen solution was pipetted into the MVAS and centrifuged to load $650 cells into each well. The excess collagen was removed, and the device was transferred into the incubator for 15 min to initiate collagen polymerization. Additional $130 cells were centrifuged into each well and allowed to adhere to the top of the tissues. Excess cells were washed off. Cell culture media were added and changed every 24 h.

Imaging
All images were acquired on a TiE A1-R laser scanning confocal microscope (LSCM) (Nikon) with standard LSCM configurations using appropriate laser lines and filter blocks. To assess morphology, microtissues were fixed in situ with paraformaldehyde for 10 min and permeabilized with 0.5% Triton-X for 3 min. The actin cytoskeleton was stained with Alexa Fluor 546 Phalloidin (Fisher, A22283), and the nuclei were stained with DAPI (Fisher, D1306).

Force measurements
After three days of static culture, we measured both stress relaxation of microtissues following a step strain and stress recovery after subsequently returning them to their initial length. All measurements were completed at 37 C with 5% CO 2 . Changes in microtissue tension were deduced from the visible deflection of the force-sensing cantilever, which was calculated by subtracting the difference in top and bottom positions. Bright field images of the tops of the cantilevers were captured at 15 frames per second during both stress relaxation and recovery. From these images, the positions of the top of the cantilevers were tracked using pattern matching with adaptive template learning in Labview. To aid in tracking the bottom positions of the forcesensing cantilevers, they were fluorescently labeled by doping the PDMS mixture with Rhodamine-B prior to curing. Fluorescence images taken from the bottom of the cantilevers were captured before and after stress relaxation and recovery experiments. The bottom position of the force-sensing cantilever was measured from these images using a simple centroid algorithm in Matlab. A validation of our approach with an elastic standard and characterization of the noise floor without an attached load are given in supplementary material I.
To investigate nonlinearities in the viscoelastic behavior of microtissues, stress relaxation and recovery were measured for various step strain values. The strain, e, was defined as the percent change in the length between the innermost edges of the tops of the cantilevers once the microtissue had fully relaxed or recovered (eðtÞ ¼ ðlengthðtÞ= length o À 1Þ Â 100).
To assess the role of individual cytoskeletal proteins in contributing to the viscoelasticity of microtissues, measurements were taken following 20min of incubation with 10lM nocodazole (Noco), a microtubule polymerization inhibitor, 5 lM blebbistatin (Bleb), a myosin-II inhibitor, 10 lM cytochalasin D (CytoD), an actin polymerization inhibitor, or 0.5% Triton-X, to decellularize microtissues. 0.5% DMSO was used as a loading control. As mechanical properties varied between microtissues, each microtissue was compared with its own pre-treatment value where indicated. To prevent crossover in response from multiple drugs, only a single treatment was administrated to any given microtissue.

Quantification of local viscoelastic deformation
To quantify the spatial distribution of viscoelastic deformation that occurs following a change in the length, local strains were estimated across microtissues using a method described previously. 30 Briefly, starting at the frame immediately following the step change in the length, inter-frame displacements were estimated every one second for 100 s in Labview at two-pixel spacing across a region of interest using a four level-pyramid based Lucas and Kanade algorithm 34 with sub-pixel precision and a window size of 17 Â 17 pixels. In Matlab, the displacement field was first smoothed using a LOWESS surfacefitting algorithm, then the inter-frame strain tensor was calculated by finding the gradient of the smoothed displacement field, and finally, the inter-frame strain tensor was integrated to estimate the local total strain field. Importantly, because the strain distribution was integrated starting the frame immediately following the step change in the length, it reflected the viscoelastic behavior of the microtissues rather than how they deform when stretched. Thereby, this metric is a measure of the spatial viscoelastic deformation response.

Data analysis and statistics
All numerical data are presented as mean 6 standard error. Statistical tests as described in the results were performed using Originlab 8.5 (Northampton, MA), with p < 0.05 considered statistically significant. Fits to stress relaxation data were performed using Matlab's curve fitting toolbox.

Ethics approval
Ethic approval was not required for the use of NIH3T3 cell lines or for any experiment conducted within this study.

SUPPLEMENTARY INFORMATION
See the supplementary material for methodological validations and further step strain experimental results. There are no conflicts to declare.

DATA AVAILABILITY
The data generated during the current study are available from the corresponding author upon reasonable request.