All-in-one rheology of multicellular aggregates


 Tissues are generally subjected to external stresses, a potential stimulus for their differentiation or remodelling.
While single-cell rheology has been extensively studied, mechanical tissue behavior under external stress
is still poorly known because of a lack of adapted set-ups. Herein we introduce magnetic techniques designed both to form aggregates of controlled size, shape and content
(magnetic molding) and to deform them under controlled applied stresses over a wide range of timescales and amplitudes (magnetic rheometer).
We explore the rheology of multicellular aggregates (F9 cells) using both standard assays (creep and oscillatory response) and an innovative broad spectrum solicitation coupled with inference analysis.
We find that multicellular aggregates exhibit a power-law response with non-linearities leading to tissue stiffening at high stress. Comparing magnetic measurements on aggregates to isolated F9 cells characterization by parallel-plates rheometry, we reveal the role of cell-cell adhesions in tissue mechanics. Thanks to its versatility, the magnetic rheometer thus stands as an essential tool
to investigate model tissue rheology.


Introduction
Mechanical behavior of biological tissues governs various processes ranging from developmental biology to cell migration or tumor invasion [1,2,3,4,5]. Biological tissues are complex 3D multicellular, heterogeneous and active materials. Their intrinsic complexity makes them challenging to characterize as this necessarily entails to perfectly control and reproduce the tissue microscopic and macroscopic architecture. In that sense, multicellular aggregates are highly relevant models of biological tissues for quantitative exploration [6]. Organized in three dimensions and usually composed of a single cell type, their cell number, cell-cell interactions, and overall size are well-controlled. They allow to decouple mechanics from genetic regulation. All these features make them simple, fully-controlled systems of choice, as testied by their widespread current use in developmental biology [7], organogenesis, tumorigenesis [8,9], drug screening and biophysical research [9,10,11,12].
Tissue mechanics is generally found to be time-dependent: elastic on short times scales (seconds to minutes) [13,14,15,16,17] and viscous on longer time scales (minutes to hours) [18]. For longer time scales (days) cell division and death may also contribute to tissue uidization [19]. Aggregate rheology should integrate these multiscale contributions as both intracellular and intercellular parameters are involved [20].
Many dierent techniques were dedicated to determine living cells mechanical features [21], evidencing a power law rheological behavior over a wide range of time scales [22,23,24,25,26,27]. Such a power-law ousts the spring-dashpot models to impose the image of a cell with no dominating timescale [22] but rather a continuum of relaxation times and processes. Single cells are thus characterized by two parameters only: the prefactor and the exponent of the power law. The prefactor quanties the overall stiness of the cell (its typical modulus), while the power-law exponent indicates the relative contributions of elastic (energy storage) and viscouslike (dissipation) behaviors. For single cells, the typical moduli are hundreds of Pascals while the exponent α is found in the range 0.15 − 0.3, thus closer to a solid (α = 0) than to a viscous liquid behavior (α = 1). Both depend on the level of applied stress, with evidence of stress-stiening [28,29] as well as of stress-softening [26], and attempts to reconcile the two opposite observations [30].
The rheology of multicellular aggregates certainly depends on cell-scale mechanical properties and should be inuenced as well by the mechanical properties of intercellular junctions [31] which exhibit viscoelastic properties [32]. However the link between cellular and junctional contributions is not straightforward. The response of suspended epithelial cell monolayers has been characterized by combining conventional viscoelastic and fractional elements with an exponent close to 0.3 [33,34]. Plus power-law rheology has been observed at the organ scale in a shear rheometer [35,36] with an exponent close to 0.2. Understanding how cell-scale mechanical properties are translated at the tissue scale remains an open question. On the experimental side, one major hurdle stems from the diculty to apply controlled stress sollicitations to controlled tissues and to accurately measure their resulting deformations.
Herein we investigate multicellular aggregates rheology. We developed an innovative magnetic compression rheometer providing an unprecedented range and accuracy of the possible stimuli, as well as a two-fold versatility. First, the rheometer provides a wide range of accessible stimulations both in terms of frequency and amplitude. Second, magnetic molding guarantees the formation of controlled cellular spheroids. Combining both features, we subjected magnetized aggregates to standard rheological tests, as well as to broad-spectrum stimulations and tested the interplay between cell-scale mechanics and tissue-scale rheology. These protocols were cross-validated using a micro-plate experiment implemented for multicellular aggregates. Overall, the results demonstrate that multicellular aggregates exhibit a power-law rheological behaviour presenting a stress-stiening. Interestingly, the stiness of this cohesive cell assembly is lower than that of the corresponding individual cell. While cell-cell adhesion and individual cell stiness do not aect the exponent governing the power law response, both contribute to the aggregate stiness.
A magnetic compression rheometer at the tissue scale In order to measure the rheological parameters of a whole multicellular aggregate, we stimulate this material mechanically with an arbitrary signal. We remotely apply magnetic forces on all cells through internalized superparamagnetic nanoparticles [37]. Previous works have evidenced that incorporating magnetic nanoparticles does not modify the cell biological activity [38]. The tissue can be stimulated at will as the magnetic cellular force is created by an electromagnet whose magnetic eld is varied in a time-dependent manner.
The magnetic eld develops a tissue magnetization M V (B) and exerts a roughly constant pulling force f V = M V (B) dB dz per unit volume of aggregate (Figure 1 and Methods). The maximal exerted force is 2.2 kN.m −3 . For z-independant sections (at surface), no simple shear is induced and deformations are limited to uniaxial squeezing. We use a magnetically inspired molding strategy [16] to design cylindrical (with a typical diameter of 2 mm and a typical height h 0 = 1 mm) as well as cuboid (1 mm long square base and 0.6 mm height) multicellular aggregates, corresponding to 10 6 and 2.5 10 5 cells respectively. Cells in these aggregates display well developed cell-cell adhesions (Figure 1e). The stress is always compressive, and we take a positive notation for simplicity. Its local value at any altitude z above the glass plate is simply given by the sum of all magnetic forces exerted on cells located higher within the aggregate: For a uniform volume force f V , the stress increases linearly with the distance from the upper surface. We dene σ as the spatially averaged value of the stress σ ≡ σ loc (z, t) z in the magnetic rheometer which reduces to σ = 1 2 h(t) f V ≈ 1 2 h 0 f V . When h 0 = 1 mm, it reaches about 1 Pa for the maximal value of the applied magnetic eld.
With this set-up (Figure 1f), the aggregate deformation varies in response to the compressive force stimulation imposed through the electromagnet. The displacement δh(t) of the upper aggregate surface is recorded by fast camera imaging and quantied through image analysis. The threshold of detection in displacement is of the order of 0.5 µm while the overall aggregate strain ε ≡ δh/h 0 (positive with the conventions used) does not exceed 5%(i.e., δh ≤ 50 µm).
Our aim is to measure the rheological response R relating the local strain ε loc along the compression axis to the present and earlier values of the stress σ loc along the same axis at the same location: ε loc (z, t) = R ({σ loc (z, t } t ≤t ). The magnetic rheometer provides the applied stress signal σ(t) and the measured strain ε(t) = δh(t)/h 0 at the scale of the aggregate. The latter is equal to the spatial average of the local deformation, ε(t) = ε loc (z, t) z . Hence, assuming a linear local rheological response R, R relates directly the macroscopic strain to the macroscopic applied stress: ε(t) = R ({σ(t )} t ≤t ). (For the nonlinear regime, see Methods).
Thanks to the electromagnet, this set-up is very versatile and aggregates were subjected to various force signals: steps, ramps, sinusoidal and broad spectrum signals with frequencies ranging from 0.1 to 10 Hz (Figure 1g-h). In this work, measurements were carried out on F9 murine cells [39], but the magnetic rheometer can be used with other cell types, provided that the cells are able to form stable multicellular structures and have a sucient magnetic potency based on endocytosis.

Power-law rheology of multicellular aggregates
We rst use two classic signals: stress steps and single-frequency stress sinusoids.
The response to a force step signal of 1.1 kN.m −3 , corresponding to a σ stepping up to σ 0 = 1 2 h 0 f V = 0.55 Pa is recorded, (Supplemental Movie 1). Here, for t ≥ 0 the applied stress is a constant, σ = σ 0 . Hence, R is simply a function of time, and we dene the function J σ 0 (t) through: ε σ 0 (t) ≡ σ 0 J σ 0 (t). We observe (Figure 2a) that it coincides with a power-law function of time over four orders of magnitude (0.01 to 100 s): Here we introduce the classical creep function J(t), and drop the σ 0 subscript for J, assuming that the response is linear. The dimensional parameter τ is set to 1 s by convention. We repeated the same assay over 40 aggregates to sample the distribution of both the exponent α and the prefactor J 0 . Figures 2b and c show their histograms and cumulative histograms (with J 0 in logarithmic scale). The latter are tted by an error function, yielding the average values as well as the standard deviation. We nd α = 0.21 ± 0.03 (mean±sd) and J 0 = (4.8 + 2.3, −0.9) 10 −2 Pa −1 (mean±sd). Our sampling is consistent with normal and lognormal probability laws for α and J 0 respectively, in agreement with [24,25]. As the dispersion between several measurements on a given aggregate is smaller than that observed between different aggregates, the magnetic rheometer is suitable to estimate the inherent distribution of rheological parameters. Further, by comparing the responses of cylindrical and cuboid aggregates, we demonstrate that the rheological parameters are independent of shape (Supplementary Figure S4).
As long as the aggregate is stimulated in its linear regime, its deformation ε(t) in response to an arbitrary imposed stress history σ(t) can be expressed in terms of the above creep function. Indeed, the stress can be viewed as a series of stress steps dσ(t ) = dσ dt (t ) dt at earlier times t ≤ t, and the deformation is expected to coincide with the sum of the corresponding creep functions: To test the power-law response obtained with one step-like stress signal, we superimposed to it a sinusoidal force with amplitude σ osc at angular frequency ω (corresponding to 3 Hz), starting a few seconds after the force step ( Figure 3b). Eq. (2) predicts that in the linear regime, an oscillatory stress yields a sinusoidal strain response with complex amplitude ε osc = J (ω) σ osc , where: with E (ω) the complex Young modulus and Γ the Euler gamma function. The measurements (Figure 3a,b,d) are consistent with the superposition of the earlier power-law response to a force step given by Eq. (1) and the oscillatory response at 3 Hz. It is thus possible to extract a complex modulus value E (ω) for each tested angular frequency ω. Figure 3d displays the amplitude |E | and phase shift ϕ of the complex modulus for several frequencies ω/(2π) ranging from 0.1 to 10 Hz. The amplitude exhibits a power-law dependence on frequency |E * | ∝ ω α over two orders of magnitude, while the phase shift ϕ appears independent of frequency, consistently with Eq. (3) which predicts ϕ = απ/2. E 0 is dened as the amplitude at The oscillatory experiments yield α = 0.22 ± 0.02, hence (2π) α /Γ(1 + α) ≈ 1.64, and E 0 = 40 ± 6 Pa, which corresponds to J 0 = 0.042 ± 0.010 Pa −1 . The result from the oscillatory regime is thus in reasonable agreement (15% less compliant) with the one from the creep experiment. While agreement between these two rheometry protocols has been observed at the single cell level [25], we extend this result to the multicellular scale. and phase (green symbols) of the complex viscoelastic modulus E * of the aggregates measured in two dierent geometries : cylinder (diamonds) and cuboid (squares) and with two distinct set-ups (magnetic and parallel-plates rheometers). Experiments were repeated for 45 cylindrical aggregates and 7 cuboid aggregates. Experimental data are tted by a power-law for the modulus (purple lines) and by a constant for the phase shift (green lines) giving an exponent α = 0.22 ± 0.01 and a stiness E 0 = 40 ± 6 Pa for the magnetic rheometer. Results obtained with a single-cell microplates rheometer for F9 single cells loaded with nanoparticles (open circles) are included for comparison.
Back-to-back comparison with a parallel-plates rheometer As a cross-validation, we used a set-up initially designed for single-cell rheology [40], the parallelplates rheometer [41], and extended it to the scale of multicellular aggregates. In their principle, measurements on cellular aggregates are similar to measurements on single cells. A cuboid aggregate is compressed between two glass microplates, a rigid one and a exible one whose stiness has been previously calibrated. Knowing the force applied to the cellular aggregate and its deformation, its compression viscoelastic modulus is estimated. The design of the exible glass microplate (the force probe) had to be adapted specically: while aggregates are far larger (around 500 µm radius) than single cells, they are softer ( Figure 3d). Thus, original spatulalike plates with wide tips (1 mm in width) and a thin neck around 5 mm long and 100 µm width) were designed to reach a typical stiness of tens of nN/µm (see Methods for details).
The sinusoidal stimulations exerted with the parallel-plates rheometer closely paralleled those exerted with the magnetic rheometer: aggregates were rst compressed, then an oscillatory deformation in the 1 − 2 % range was applied. Studying the response of 7 cuboid aggregates, we nd the same power-law rheology, characterized by an exponent α = 0.27 ± 0.04 and an elastic modulus E 0 = 45 ± 12 Pa. Quantitatively, cellular aggregates thus exhibit a slightly more prominent viscous-like component (higher α) than when tested with the magnetic rheometer, while estimates of the linear elastic modulus E 0 are consistent between the two instruments. This comparison corroborates the results obtained with the magnetic rheometer and conrms the power-law rheology of multicellular aggregates.  Since multicellular aggregates are essentially made of cells adhering to each other, their mechanical behavior is expected to depend on the rheological properties of the individual cells as well as on cell-cell adhesion. To test these dependencies, we rst investigate isolated F9 cells mechanical behavior.

Contribution of single cell mechanics
We probe individual F9 cells in a single-cell microplate rheometer [40] (Figure 3d). Briey, the applied deformation on individual F9 cells is in the 10−50 % range, and the elastic modulus is measured at dierent frequencies up to 0.5 Hz. The exponent found for F9 single cells (N = 11), α SC = 0.27 ± 0.03, is consistent with the one measured in the microplate rheometer for F9 cell aggregates. However, F9 single cells are signicantly stier than F9 cell aggregates, with an elastic modulus E SC 0 = 240 ± 80 Pa, leading to J SC 0 = (8 ± 2) 10 −3 Pa −1 . At the single-cell level, α is robust to most inhibitors of cytoskeletal activity [26]. We test it at the scale of multicellular aggregates by investigating the role of the actin cytoskeleton. Aggregates are submitted to latrunculin A, an inhibitor of actin polymerization. On single cells, latrunculin A softens the cells by 2 to 5-fold, depending on cell type and inhibitor concentration [26,22,42]. On aggregates, latrunculin A disrupts the dense cortical actin meshwork visible on normal conditions (small patches appear, Figure 4d). Quantitatively, the exponent α does not change signicantly (Figure 4e). By contrast, latrunculin A-treated aggregates exhibit a twice higher compliance (Figure 4f), as also found in single cells. Changes of individual cell mechanics thus have a direct impact on aggregate stiness.

Contribution of cell-cell adhesion
We next examine the role of cell-cell adhesion in the mechanical properties of cell aggregates, rst by varying the maturation time from 16 h to 48 h, second adding EGTA to inhibit the formation of intercellular adhesions. Indeed, increasing the maturation time leads to more cohesive aggregates : cell-cell adhesions develop with time from the periphery to the center of the aggregate (Figure 4a-b). Conversely, under the action of EGTA, a calcium ion chelator decreasing intercellular adhesion forces without aecting single cell mechanical features [26], the aggregates appear looser, with a lower expression level of cadherins at the membrane (Figure 4c). In these conditions, while the exponent α is not signicantly modied by the degree of intercellular adhesion (Figure 4e), the compliance factor J 0 depends on the cell-cell contacts (Figure 4f): when the aggregate cohesiveness is reduced, either by shortening the maturation time or under the action of EGTA, the aggregate becomes softer (larger value of J 0 ). Conversely, the maturation of cell-cell adhesions stiens the aggregate and at 48 h of maturation, the compliance of the aggregate becomes remarkably close to the one of isolated F9 cells.

Non-linear and broad spectrum rheology
After determining the linear rheological features of the multicellular aggregates , we next probe the non-linearity of the power-law response.
A stress ramp σ(t) =σ 0 t (whereσ 0 is a constant stress rate) is expected to generate a power-law response in the linear regime. Indeed, Eq. (2) predicts: The actual deformation (Figure 5a) shows an apparent deviation from linearity for volume forces above 1.2 kN.m −3 (corresponding to 0.6 Pa stress when h 0 = 1 mm): the observed deformation is smaller than expected from Eq. (4). This stress stiening is conrmed by recording the creep responses to dierent force steps corresponding to stress amplitudes ranging from 0.3 to 1.2 Pa. For each value of the stress, the creep deformation still follows a power-law with almost unchanged exponent (Figure 5c), while the compliance prefactor J 0 is reduced for higher stress levels ( Figure 5d).
To systematically quantify the nonlinearity of the aggregate rheology in a single experiment, we introduce an original protocol in which we impose a large amplitude, broad spectrum stress signal (example on Figure 6a (Figure 6c,d). These signals and their time derivatives form a trajectory in the system phase space. Using a broad spectrum signal allows to explore the phase space eciently and to thus sketch the manifold that corresponds to the constitutive equation.
Various, possibly non-linear, rheological models can be tested and compared quantitatively (Methods). Inspired by the power-law creep linear response and by the litterature [43,44], we tested a linear fractional model, and obtained a dierentiation exponent α = 0.25 ± 0.02 and a compliance prefactor J 0 = 0.045 ± 0.010 Pa −1 (N = 8). It provides a t with reasonable, yet far from perfect, quality (Figure 6b, blue line). Therefore, more general, nonlinear equations with one fractional element were investigated. Using the Akaike Information Criterion (Table 1), we identify as the best model the following nonlinear fractional stress-strain relationship (Eq. (8) in Methods): where I τ α (f (t)) is the α-order fractional integral of the signal f (t) [44]. Maximal likelihood optimization yields the following values of model parameters: the dynamical exponent α = 0.22±0.01, the compliance J 1 = 0.13±0.02 Pa −1 and the critical stress for non-linear eects σ c = 0.5 ± 0.1 Pa (N = 8). Note the good agreement between the value of the order of dierentiation α and previous estimates obtained with classical rheological stimulations. Including either strain or stress thresholds in the model did not improve agreement with experimental data, which ruled out the relevance of plastic behaviour [6] in the F9 multicellular aggregates in that range of deformations.
Eq. (5) also predicts strain histories when the aggregate is subjected to simpler stress signals. The creep response to a step stress sollicitation is a single power-law function, as in the experimental observations, with a constant exponent α (Figure 5c) while the stress dependence of the estimated eective compliance reproduces the one established by step forcing (Figure 5d). Interestingly, an eective exponent that decreases with the applied force is obtained by complementing the model with a second fractional integral acting on the strain ( Figure S5, see Eq. (14) in Methods for additional details).
For each sample, a single broad spectrum experiment allowed us to recover the rheological features identied from stimulations based on both creep functions and harmonic signals and to generalize them in the form of a dierential equation (Eq.5). Taken together, the analysis of all collected data conrms that the rheology of multicellular aggregates is best described by a fractional rheological law with a force-dependent prefactor (stress stiening).

Discussion
Our work shows that magnetic cells may act as actuators to apply controlled internal forces on the aggregate they form. This set-up is thus able to test the rheology of model tissue. Given that estimates of single cell rheological parameters often dier by orders of magnitude when measured by dierent techniques [21], the agreement between magnetic and microplate rheometry is remarkable and validates this internal magnetic rheometry method. The magnetic rheometer allows to exert arbitrary compressive stress signals in the range from 0.1 to 10 Hz. It oers a unique possibility to implement broad spectrum rheology, an approach which aims at providing the constitutive equation of the material using a single experiment per sample. This is reminiscent of the optimally windowed chirp protocol introduced in [45,46], a protocol carried out in Fourier space. In contrast, the present broad spectrum protocol is performed in real space, thereby allowing to fully characterize nonlinearities with only one signal.
The direct measurements of model tissue rheology reported here show that multicellular aggregates behave as viscoelastic materials well described by a nonlinear fractional model with stress stiening. The nonlinear form given in (5) is reminiscent of the dependence found experimentally between stiness and stress for single cells [27,47]. The creep function of the aggregate is a power law function. This behavior is the hallmark of the presence of a continuous distribution of relaxation times, already identied at the individual cell level [25]. The value of α close to 0.2 puts multicellular aggregates in the same range as individual cells in terms of dynamics. This dynamic exponent also proves robust to most perturbations. By contrast, stiness is remarkably sensitive not only to intracellular cytoskeleton (reecting the single cell rheology change in the case of latrunculin-A actin disruption) but also to cell-cell adhesions. The compliance J 0 decreases (stiening) as the aggregate matures ( Figure 4). Indeed as cellular adhesions develop, individual cells become more sharply faceted and expulse interstitial uid from the aggregate [48] and reduce their sliding possibilities in the interstitial regions. This change in geometry, as in liquid foams, corresponds to a stiening of the entire structure since it increases the level of local deformations required to match a given macroscopic deformation. Alongside, we observed a softening (Figure 4) in the presence of EGTA, which is known to weaken cell-cell adhesion and results in less facetting cells.
Single cells appear stier than the multicellular aggregates ( Figure 3), even after maturation up to 48 hours, (Figure 4f). That may result from boundary conditions being very demanding in the single-cell experiment, where a signicant fraction of the cell surface is connected to rigid plates and the cell has to largely deform. By contrast, in measurements on the entire aggregate, the overall deformation imposed can be reached through a combination of cell deformations a b c d modes since the surroundings of each cell is then softer. Note that this argument also applies to 2D situations, making 2D epithelia softer than single epithelial cells, as measured in a suspended geometry [49]. Finally, both intercellular and intracellular components contribute to the multicellular aggregates. Stress stiening, observed here with cell aggregates, has previously been reported for complex polymeric materials as well as living cells [50,51,52]. For polymers, alignment along the direction of the applied force leads to an increased resistance, which contributes to the stiening of the network. This contribution, present at the scale of the cytoskeleton, may extend to the cellular, as well as to the multicellular scale. In tissues, cell-cell arrangements may also play a role as stress can disrupt adhesion and facilitate in turn cell-to-cell movements leading to a stress softening. However, since maximal deformations in our experiments are in the 5% range, the disruption of cell junctions seems unlikely to be relevant but cell deformations and sliding may contribute. Observed non-linearities may therefore be mainly associated with the cellular contribution. One caveat concerns the threshold stress values above which stress stiening is observed, of the order of 1Pa here, vs. 10 − 1000 Pa in single cells [47,53] and several hundreds Pa in monolayers [33,34]. However, the magnetic rheometer allows to exert comparatively small stresses, while most other techniques involve larger stresses. Thus a multi-step stress stiening at dierent force levels may be implied.

Conclusion
The rheology of multicellular aggregates presents many similarities with single cell and cell monolayer rheology in spite of marked dierences in terms of scales and structures. We nd that their responses to a step-compression both follow a power-law in the range from tenths of seconds to minutes, with a comparable dynamical exponent. The similarity in behavior extends in the non-linear regime as multicellular aggregates exhibit stress stiening above stresses of the order of 1Pa. However, the behaviour of multicellular aggregates diers from that of single cells. Their eective compliance is higher than that of single cells, and depends on the maturation of cell junctions. The magnetic rheometer gives the possibility to exert small forces and study precisely the small deformations of millimeter scale multicellular aggregates. Together with model inference and model selection, broad spectrum analysis allows to decipher in a single acquisition the rheological model that best describes a sample. Thanks to its versatility both with respect to the nature of applied signals (from step to multi frequency force signal) and to the nature of biological samples (cell types, drugs and inhibitors), this all-in-one magnetic rheometer either on its own, or coupled with inference analysis could become a key player for the mechanical exploration of biological model tissues.

Nanoparticles and reagents
The nanoparticles used to magnetically label the F9 cells are citrated iron oxide maghemite (Fe 2 O 3 ) superparamagnetic nanoparticles with a characteristic diameter δ 0 7 − 8 nm (polydispersity 28 %, standard deviation s 2ñm). The magnetization curve ( Figure S1) is described by a Langevin function weighted by the nanoparticles size δ: where M np is the magnetization at saturation of a maghemite nanoparticle, ξ(δ) = µ 0 HM np (δ)/(k B T ) and L(x) = coth(x) − 1/x. The size distribution P is approximated by a log-normal law with the characteristic diameter δ 0 and the standard deviation s. As particles are superparamagnetic, M is independent of the direction of the magnetic eld and applied forces are pulling ones.
To inhibit actin lament formation, Latrunculin A (1 µM) purchased from Sigma is used in some experiments during aggregate formation and stimulation. EGTA (0.4 mM) may be added at the same stage to act as a calcium chelator and decrease cell-cell cadherin mediated interactions as well as individual calcium inux.
Immunouorescence was performed using SiR-actin (2 µM for 1 hour) to stain actin laments and DAPI to label nuclei. E-cadherin antibody (ECCD-2 from Thermosher used at 1 : 50 for 2 hours) coupled with a secondary anti-rat GFP antibody is used to enhance this cell-cell adhesion molecule. Magnetic cell labeling and magnetic formation of multicellular aggregates F9 murine cells are grown in DMEM supplemented with 5 % fetal bovine serum (FBS) and 1 % antibiotics (penicillin/streptomycin) at 37 • C and 5 % CO 2 . When cells are at conuence, nanoparticles solution ([Fe] = 2 mM in RPMI) is added for 2 h. Nanoparticles are internalized in cells through endocytosis pathways and are nally localized within intracellular vesicles identied as endosomes [54]. After incubation, cells are washed with fresh RPMI and detached with trypsin to be seeded in 2 % agarose gel cylindrical molds (diameter 2 mm, height 1.5 mm) or cubic molds 1 mm square) obtained by forming a network of small magnetic cylinders (or cubes) within an agarose mold ( Figure S2). Indeed by removing the magnets below, and taking o the magnetic cylinders or cubes within the gel, we achieve to have wells to seed cells within. The seeding of cells (1 million) in the wells is then magnetically driven by the magnets below the dish ( Figure S1A-B). Cell-cell contacts develop overnight to produce a cohesive cylindrical tissue that can be removed from the mold (Figure 1a-d). The resulting aggregates are removed from the well after

Magnetic eld applicator
The magnetic eld used to deform the magnetized aggregate was generated with an electromagnet composed of a custom-made water-cooled coil (6 Ω, 14 mH) wrapped around a soft iron polar piece (diameter 55 mm, height 85 mm) allowing the passage of a current up to 10A. In its upper part, the diameter of the polar piece goes down from 55 mm to 5 mm. The top of the polar piece was 3 mm below the center of the aggregate. The electromagnet tension was controlled with a TGA 1242 signal generator (Aim&TTi) and amplier (Kebco BOP).
The shape of the polar piece was optimized to maximize the magnitude of the magnetic eld gradient dB/dz along the z coordinate (30T.m −1 when the magnetic eld reaches its maximum value 120 mT) as well as its homogeneity in the aggregate volume. Nanoparticles magnetize in the presence of the magnetic eld B. As a result, the aggregate displays a magnetization M V = M sat V φ(B/B sat ) per unit volume (for φ, M sat V and B sat , see above). The magnetic eld generated at the center of the aggregate linearly depends on the applied tension (see Figure S3a) and the magnetic eld gradient was almost homogeneous (less than 20 % variation over a size exceeding one millimeter), whatever the applied tension ( Figure S3b) ensuring a homogenous force on cells.
The set up was calibrated and characterized in order to determine the bulk attening magnetic force generated as a function of the applied tension V (t) using the set-up as a tensiometer [16]. Ferrouid drops (aqueous solution of the nanoparticles used for cell labelling) of various size deposited on a superhydrophobic surface in a silicon bath, were attened. Fitting each experimental prole with the Laplace equation for capillarity in non-wetting conditions provided the surface tension. This measurement was validated by comparison with classic sessile drop experiment ( Figure S3c). The maximum value of the magnetic force density was measured to be close to 2.5 kN.m −3 .

Signal generation and synchronization
The tension signal is transferred to the generator using the Waveform Manager Plus software. To synchronize the deformation signal with the force generation, a LED, powered by a pulse generator that is synchronized with the signal generator, was used to overexpose a single image of the movie at the beginning of the deformation sequence.

Image analysis
Deformations were analyzed using the image correlation based tracking algorithm on ImageJ: Tracker Install (developped by Olivier Cardoso https://sites.google.com/site/oliviercardoso/Home/plugins). It extracts the displacement of the aggregate upper interface as a function of time (δh(t) in Figure 1f-h). Vibrations of the setup (mainly due to the water circulation use to thermalize the container) were compensated by tracking a xed point on the container bottom (adsorbed dust for instance) and subtracting its motion from the interface displacement. The precision on the interface displacement is ±0.5 µm. Sampling frequencies are between 10 to 100 Hz depending on the experiments. Shape-independance of cellular aggregates rheology Figure S4: Inuence of the shape on the power law exponent. Comparison of the compliance a and of the α exponent b in two dierents geometries : cuboids and cylinders. No dierences are noticeable.

Microplate rheometer
A microplate rheometer was used both for single cell and for aggregate measurements [41]. In such a device, the sample is compressed between two parallel microplates. One is rigid. The other is exible. Its deection δ is related to the force F = kδ applied to the sample through the calibrated stiness k.
However, since the set-up was initially designed to test single cells (size around 10 µm), the exible microplate had to be adapted to both the size (500 µm) and the stiness of cell aggregates. Thus, new spatula-like plates have been created with wide tips (1 mm in width) and a thin neck (around 5 mm long and 100 µm width) to be able to test large but very soft objects like cell aggregates. The stiness of the spatula was 27.8 nN/µm. For single cell measurements, as F9 cells are very soft, we used a microplate of 0.6 nN/µm. The exible microplates were otained by heating and pulling borosilicate plates of 10 cm x 1 mm x 200 µm [55]. The stiness of the exible microplates was calibrated using a standard microplate itself calibrated following the protocol detailed in [41]. The setup was mounted on a Leica DMIRB inverted microscope (Leica Microsystems, Rueil-Malmaison, France), and samples were visualized using a x10 objective for aggregates and a x63 objective for single cell measurements. The stress σ applied to the sample was dened as the the force F = kδ divided by the contact area A between the sample and the plate: σ = F/A. For single cells, the contact areas between the plates and the cell were estimated by assuming a circular contact, i.e. A = πD 2 /4. The apparent contact diameters D F and D R (respectively for the exible and rigid plates) were measured on bright eld images. Since D F and D R are usually dierent, we retained their mean square value, leading to σ = 8F π(D 2 F +D 2 R ) . As for the cuboid cell aggregates of dimensions L × L × h, they initially lie on their square face. They are seized laterally between both microplates (Figure 3c). Hence, the contact surface is rectangular and its area A = Lh is determined through bright eld imaging. L is measured directly during the experiment. It may slightly dier from the measurement obtained before the microplate handling. h is estimated by using images taken before the microplates grasped the aggregates: a correction is applied to take into account the possible deformation due to handling (it is computed from the corresponding L values). Finally, the uniaxial cell or aggregate strain perpendicular to plates was dened as ε = L−L 0 L 0 , where L 0 (resp. L) is the distance between the parallel microplates before (resp. after) compression. The apparent modulus is then given by E = σ/ε. Dynamical measurements were carried out in the frequency range from 0.1 to 3.2 Hz for cell aggregates, and from 0.1 to 0.4 Hz for single cells. In the latter case, measurements at higher frequencies were dicult to analyze: the very low stiness of the exible plate (at the limit of the technique), required to measure the very low viscoelastic modulus of F9 cells, makes it very sensitive to all viscous ows induced by the plate motion in the complex geometry around the plates and between both plates around the aggregate.

Broad spectrum signal
The stress signal we applied for the broad spectrum method is a superposition of several sinusoidal signals: where σ 0 is used to tune the amplitude and ϕ(t) = tanh(λ t) is an envelope function used to smooth the beginning of the signal over the rst few seconds. The frequency range is chosen as 0.1-10 Hz. The n = 32 frequencies f i are equally spaced logarithmically within this range. The phases φ i are chosen randomly within [0, 2π]. In order to damp high frequencies, the magnitude of the real prefactor ρ i is chosen to decrease with the frequency. More precisely, ρ 2 i is drawn uniformly within [0, f 2 min /f 2 i ], which ensures that the complex amplitude ρ i e i φ i is distributed uniformly in the disk of radius f min /f i , a scaling reminiscent of 1/f noise.

Rheological models
Generally speaking, a rheological model may be dened by a p-dimensional parameter vector θ and by an equation of the form ε loc (t) = R({σ loc (t )} t ≤t , θ), where R is an integro-dierential operator acting on the function σ loc (t).
The rheological models we considered, dened in the following paragraphs by order of increasing complexity, contain fractional integrals [44], dened as: where an arbitrary time scale τ is set to τ = 1 s and where I τ α (f ) has the same dimension as f .
Linear fractional rheology This model has two parameters (both positive): where α is dimensionless and J 1 is in Pa −1 .
Stress-strain relationship The volumic nature of the applied magnetic force implies that the stress is not uniform within the samples. Rather, it increases approximately linearly from the free surface of the aggregate up to its maximum value at the substrate surface: σ loc (z) = σ max (1 − z/h 0 ) = 2σ u, where σ is the spatially averaged applied stress and u = 1 − z/h 0 is the reduced distance from the upper surface of the aggregate. As a result, the local deformation ε loc also depends on altitude via the (nonlinear) dierential equation ε loc (z, t) = R {σ(z, t )} t ≤t , θ .
Correspondingly, the overall aggregate deformation can be obtained by averaging the local deformation eld. For a response R and a homogeneous magnetic force f V , the macroscopic strain is given by: When the rheology R is linear, the spatial average commutes with R, and global variables are also related through the same operator: ε(t) = R {σ(t )} t ≤t , θ . For instance, for the linear fractional constitutive equation (9) one also has ε = J 1 I τ α (σ) and for a step stress of magnitude σ 0 , using Eq. (8) we nd a power law response: Comparing this result with denition (1), we obtain the creep compliance J L 0 = J 1 /Γ(1 + α) for the linear case.
More care needs to be exercised when the rheology R is nonlinear. A useful example is the case of Model (10), for which the expected overall deformation reads: For a step stress, the compliance now depends on the step magnitude σ 0 in this non-linear case (see Figure 5e): Parametric inference In order to assess to what degree the rheological model, dened by R, is suitable to describe the experimental data {ε exp (t)} t , we assume that the experimental noise φ(t) is white and Gaussian (amplitude s) and contributes additively within the measured deformation value: We then determine the values of the parameters θ and s, i.e., p + 1 real numbers, that maximize the likelihood function for N t time points (typically N t 7 10 3 here). Indeed, the likelihood L indicates how likely it is to observe the actual experimental data {ε exp (t)}, given the rheological model (encoded in R), the parameter set ( θ, s) and the imposed stress history {σ(t)}. The larger the likelihood, the better the model and its parameters t the experimental data.
Once the parameter values of each model have been optimized, we need to be able to select the best model. One diculty is that the various models contain dierent numbers of parameters. In order to avoid both too simple models and too complex models (since overtting may reduce reproducibility and/or predictability), we compute for each model the Akaike information criterion [56] dened as AIC = −2 logL + 2(p + 1) .
The AIC provides an objective estimation of the quality of each model, penalizing the number of parameters p and favouring a high value of the likelihood function L [57]: The absolute AIC values are not particularly meaningful, but their relative values provide a comparison between the dierent models, where the best model corresponds to the lowest AIC value. On Figure 6b, the experimentally measured deformation ε(t) (black points) is compared with the predictions provided by the linear fractional model (9) and the non-linear fractional models (10) given the applied stress σ(t) displayed on Figure 6a. The linear fractional model (9) fails to reproduce the time evolution of the measured deformation, while the non-linear fractional model (10) yields satisfactory predictions. Inspection of Table 1 shows that estimated AIC values agree with these qualitative observations. Further, all nonlinear fractional models (10-12) have a lower AIC than that of the linear fractional model (9). Among (10-12) the lowest AIC is obtained by (10). The AIC obtained for (10) is consistent with that of the more complex model (13); however most parameter values obtained for (13) are consistent with zero. For these reasons, we designate (10) as the best fractional model with one fractional integral.
Finally, introducing a second fractional order in model (14) allows to predict almost perfectly the deformation history ( Figure S5a). The relevance of this model is conrmed by comparing phase space plots on Figures 6cd and S5b. However, numerical integration shows that the power-law behavior then predicted for step stress assays is characterized by an exponent that depends on the stress amplitude σ 0 , at variance with experimental observations (see Figure S5c). Further, the model fails to reproduce quantitatively the σ 0 dependence of the creep compliance J 0 ( Figure S5d). Figure S6: Supplemental movie of an aggregate submitted to a force step corresponding to a 2 Pa sollicitation during 10 s Large eld movie of an aggregate while it is attened by a maximal force step for 10 s.