Viscoelastic Networks: Forming Cells and Tissues

Spatiotemporal changes in viscoelasticity are a key component of the morphogenesis of living systems. Experimental and theoretical findings suggest that cellular- and tissue-scale viscoelasticity can be understood as a collective property emerging from macromolecular and cellular interactions, respectively. Linking the changes in the structural or material properties of cells and tissues, such as material phase transitions, to the microscopic interactions of their constituents, is still a challenge both at the experimental and theoretical level. In this review, we summarize work on the viscoelastic nature of cytoskeletal, extracellular and cellular networks. We then conceptualize viscoelasticity as a network theory problem and discuss its applications in several biological contexts. We propose that the statistical mechanics of networks can be used in the future as a powerful framework to uncover quantitatively the biomechanical basis of viscoelasticity across scales.


INTRODUCTION
The viscoelastic or material properties of cells and tissues are key regulators of cell and tissue growth, motion, and homeostasis [1][2][3][4][5]. Viscoelasticity allows living systems to preserve a basic architecture due to their solid-like characteristics, but also at the same time to dynamically reorganize in different shapes and patterns due to their viscous-like characteristics [4,[6][7][8]. Cellular-scale viscoelasticity influences several single-cell functions such as shape, division, and motility, and it is predominantly determined by the physical properties of the underlying cytoskeletal networks [8]. Tissue-scale viscoelasticity was shown to be important in collective morphogenetic processes such as tissue folding, spreading, wound healing and migration, and it is mainly determined by the interplay of cellcell and/or cell-extracellular space interactions [2,4,5]. Advances in biophysical tools measuring viscoelasticity [4,9,10] have revealed an essential and physiologically relevant link between material properties and morphogenesis [11][12][13], opening the challenge to now understand how emergent viscoelasticity is regulated by, and in turn, regulates the mechanochemistry of living systems.
A material is viscoelastic if it displays both viscous and elastic behavior [14]. Our knowledge of viscoelasticity mainly comes from material sciences, where certain physical parameters are welldefined for non-living materials such as glasses, rubbers, metals and polymers [14,15]. Viscoelasticity of such materials is evaluated from the degree of deformation upon constant force application and release, an experimental procedure called creep and recovery test ( Figure 1A). Solid-like objects deform shortly and reversibly under constant force, whereas fluid-like objects irreversibly increase their deformation as long as a force is exerted. Viscoelastic materials exhibit characteristics of both solids and fluids: at short time scales they deform elastically and, at long timescales, they behave as viscous fluids. On the theoretical level, viscoelasticity was independently modelled by Maxwell and Kelvin in the 19th century [14]. Both models abstract a viscoelastic system as a composite structure containing an elastic spring connected to a dashpot containing a viscous fluid ( Figure 1B). The difference between Maxwell's and Kelvin's approaches comes from the disposition of the system: In Maxwell's model, the dashpot and the spring are connected in series, whereas in Kelvin's model, the two components are connected in parallel. As a result, when force is applied on the system, in the former, the deformation of the spring will drag the dashpot, describing the deformation of a viscoelastic fluid; whereas when force is applied in the latter, the dashpot response is restricted by the deformation of the spring, describing thus the deformation of a viscoelastic solid [16]. Although biological materials display characteristics of viscoelastic materials, due to their heterogeneous composition and dynamic nature, fitting such models and frameworks from material science, is conceptually and technically still difficult [16,17]. For instance, during experimental measurements of cell and tissue viscoelasticity, typically an external force is applied to the system such as via a micropipette, or a magnetic field [13,18,19] for a certain time window, during which its deformation is monitored ( Figure 1C'). From such measurements, parameters such as elastic modulus, viscosity, and yield stress can be extrapolated (Supplementary Box S1) [16,20,21]. However, whether such measurements are relevant to the time window of morphogenesis still needs to be addressed, since such measurements show how a tissue deforms when applying an exogenous force for a certain time. How these features compare to the magnitude and duration of the endogenously exerted forces is still unclear [9,22,23]. In addition, during such experimental measurements, the applied force typically triggers a large deformation of the material (nonlinear regime). In contrast however, most theoretical frameworks predicting viscoelasticity upon deformation, are based on forces triggering infinitesimal deformations (linear regime) [24]. Last, similar to the viscoelasticity of non-living materials, viscoelasticity of cells and tissues is a property arising from their underlying structure, defined by the way macromolecules and cells interact [1,8,15,[25][26][27][28][29]. Several theoretical frameworks have for long been used in non-living materials to link the microscopic structure to the macroscopic viscoelastic properties (Supplementary Box S1), showing that viscoelasticity may behave as an emergent property [27][28][29][30][31][32], in the same way, e.g., temperature arises non-trivially from microscopic particle motion in statistical mechanics [14]. How macroscopic viscoelasticity can be predicted by the interactions of the microscopic constituents of living cells and tissues is an open question at the interface of statistical and soft matter physics with molecular and cell biology, and the main topic of this review.
An intriguing empirical observation is that the material properties of the microscale components of cell and tissue viscoelasticity, such as the cytoskeletal elements and the cells, respectively, usually do not match the macroscale material properties of cells and tissues [15,29,33]. Macroscopic viscoelasticity frequently exhibits nonlinear changes that are not observed at the microscopic level. Such examples have been experimentally detected such as the strain-stiffening response of the cytoskeleton networks [26,30,34], phase transitions in the energetic costs of cell movements [12] or abrupt changes in tissue viscosity [13] (Supplementary Box S1, Figures 1C",D"). In the above cases, the mechanical resilience of the individual microscopic components to forces falls short in explaining the macroscopic viscoelastic changes, and thus probing the pattern of interactions between the components instead, is key.
Such nonlinear phenomena set a number of challenges to the theoretical understanding of cell and tissue dynamics. Theoretical analyses of cell and tissue material properties are typically addressed from the biopolymer or cell level, respectively. In the first theoretical approach, the microscopic basis is the mechanical properties of the biopolymer filaments building the cytoskeleton [35] and the macroscopic viscoelastic features are derived from the network geometry and local topology of the filaments [26]. Nonlinear phenomena such as strain-stiffening of cytoskeletal networks have been probed with such models. In the second theoretical approach, the microscopic basis is the tiling patterns of the constituents, such as the cells forming a tissue. In such modeling frameworks, mainly represented by the vertex models (Supplementay Box S1) [25,36,37], rheological properties such as rigidity and fluidity are inferred from the energetic costs for cells to independently move through the tissue [27,28,[38][39][40][41]. This viewpoint comes from the fundamental observation that material deformation can only take place through cell-cell rearrangements. Nonlinear phenomena such as jamming transitions (Supplementary Box S1) have been predicted with such frameworks [27,28,[40][41][42]. A third theoretical approach that is not as frequently applied in active viscoelastic systems, but has been used so far to probe material properties across scales, is network theory [29,31,43,44]. In this framework, the starting point is the topology of the network, e.g., how the system's constituents are connected between them. Of particular relevance for these approaches is the concept of percolation (Supplementary Box S1). Percolation refers to a wide range of phenomena where a sudden shift in the macroscopic properties of a system made of microscopic, interacting units, is observed when a certain threshold of connectivity at the microscopic scale is overcome [45]. In material sciences, percolation transitions underlie many sudden, qualitative changes in the behavior/response of the material, including, among other, shifts in rigidity and force transmission properties [46][47][48][49][50][51][52][53][54]. A paradigmatic example of the role of percolation theory in explaining material properties is found in the exploration of the emergence of cracks when the material is under stress. In this context, the length and width of cracks emerging in the material increase dramatically when the system approaches the critical point of rigidity percolation [55,56] (Supplementary Box S1). Several forms of percolation theory have been applied to cytoskeletal networks, fiber networks and, recently, to "cellular" networks -tissues -to map material properties [29,31,43,44]. As is the case with inanimate materials, the structure of interactions at the microscopic level--and the potential nonlinear shifts arising from small changes in them--are supposed to underlie and, ultimately, explain, the emergence of macroscopic material properties like viscoelasticity.
In this review we summarize and discuss experimental and theoretical work probing cell and tissue viscoelasticity as an emergent property. We will first introduce experimental findings on how viscoelasticity emerges in cytoskeletal networks, extracellular matrix fiber networks and tissues. We will then summarize and classify theoretical frameworks supporting such experimental findings that address cell and tissue material properties. Finally, we will discuss the potential of applying network theory to predict viscoelasticity, and speculate how such an approach could impact our biophysical understanding of the material properties of living systems and their morphogenesis.

VISCOELASTICITY AS AN EMERGENT PROPERTY: EXPERIMENTAL OBSERVATIONS Cellular Viscoelasticity
Changes in cellular-scale viscoelasticity are key for cell physiology [57]. Both the cytosol and cytoskeleton contribute to cellularscale viscoelasticity, with cytoplasmic viscosity dominating processes of macromolecular movement [58] and cytoskeletal viscoelasticity influencing cell morphology, motility and division [8,57]. Given that experimental work suggests that the cytoskeleton is the major determinant of cellular viscoelasticity [59,60], we focus here on its viscoelastic properties. The cytoskeleton is the underlying biopolymer scaffold of living cells, and its viscoelasticity offers a balance between dynamic reorganization and maintenance of the cell body. Cytoskeletal viscoelasticity is a complex phenomenon, since it arises from the mechanical properties and interactions of at least three different biopolymers: actin, microtubules and intermediate filaments.
These fiber structures further self-organize into filaments and heterogeneous networks through mechanisms of entanglement, branching, crosslinking and bundling ( Figure 1C) [61][62][63][64]. Such mechanisms involve several types of linkers, such as crosslinking/bundling proteins a-actinin and fascin, and motor proteins like myosin and kinesin [65]. In general, the mechanical properties of the cytoskeleton depend on the physical properties of the individual filaments, the pattern of linkages between the filaments, and the geometry of the filament arrangement [66,67]. Although cell viscoelasticity is a result of the heterogeneous mechanical properties of all cytoskeletal elements together, our understanding of this process comes mostly from studying each type of element separately (for extensive reviews see [1,3,8]). At the single-filament level, each cytoskeletal biopolymer displays different physical properties, assessed by the ratio of its persistence length (l p ) over its contour length (l c ) (Supplementary Box S1). When l c > l p , the polymers are flexible, and this is the case for the intermediate filaments, which display the shortest l p (from 200 nm to ∼1 μm) and are the softest among the cytoskeletal elements [68,69]. Actin filaments have a higher l p (from 3 to 17 μm) and are semiflexible and microtubules exhibit the highest l p (>1 mm) and are stiff polymers [70]. Already at the single filament level, semiflexible polymers such as actin and vimentin, display a nonlinear increase of their shear modulus at different strain amplitudes [71] (Supplementary Box S1). This nonlinear force-extension relationship becomes more apparent at the network level, where new material properties emerge that are absent at the single filament level. Cytoskeletal elements build networks via various forms of filament interactions that influence the viscoelastic behavior of the whole network. For example, transient and non-covalent interactions with crosslinks turn the network into a viscoelastic material, whereas covalent interactions turn the network into an elastic material [3]. Experimental measurements of reconstituted cytoskeletal networks ( Figure 1C') revealed nonlinear force-extension relationships such as stress-stiffening in the presence of tensile load ( Figure 1C") but also stress-softening in the presence of compressive load [5,30,33,62,71,72] (Supplementary Box S1). Actin and intermediate filament networks are highly strainsensitive, with 10-100 times stiffening appearing at very low strains [71,73]. Experimental work suggests that depending on the density and interactions in the network e.g. dense vs sparse, the macroscopic material properties change following a welldefined phase diagram [30]. Gardel and colleagues have shown that the addition of crosslinkers results in the formation of rigid solid networks, with elastic modulus several orders of magnitude greater [30]. Even when using flexible crosslinkers like filamin, the network displays nonlinear increases of the elastic modulus when increasing strain, reaching values that match the elastic moduli of cells [34]. Similarly, addition of molecular motors to reconstituted actin networks, such as myosin II, leads to a sharp increase of the elastic modulus [74,75]. Microtubules on the other hand, and also weakly cross linked actin networks, decrease their modulus as the applied stress is increasing displaying a stress-softening response [30,[76][77][78]. In conclusion, experimental work shows a very rich collection of nonlinear macroscopic viscoelastic behaviors of the cytoskeleton that represent a challenge for theoretical understanding at the microscopic level.

Extracellular Matrix Viscoelasticity
Besides the emergence of material properties in intracellular cytoskeletal networks, similar behaviors are observed in networks of the extracellular matrix (ECM), the non-cellular material backbone spanning cells and tissues. ECM viscoelasticity is fundamental in cell migration, tissue morphogenesis, organ development, and cancer progression [5,79,80]. The ECM is a heterogeneous network composed by several biopolymer filaments, such as fibronectin, laminin, and collagen, that exhibit various persistence lengths. This can range from the 4-8 nm persistence length l p of flexible hyaluronan biopolymers to a few millimeters' persistence length l p of stiff collagen fibers [81,82]. Given that the ECM is composed of various proteins, enzymes and polysaccharides, probing its viscoelasticity becomes highly complex. Due to its covalent nature of crosslinking, in contrast to the cytoskeletal networks, the ECM is considered an elastic-like network [83]. Collagen networks also display a strain-stiffening response that is in this case emerging from the network level and specifically its connectivity [73,84,85] (Figures 1C,C'''). In this case, the nonlinear behaviour emerges for strain of only 10% increase where stiffness increases by 100x before network rupture [44] ( Figure 1C'').

Tissue Viscoelasticity
Similarly, tissue-scale viscoelasticity has been recently experimentally measured to undergo nonlinear changes [11][12][13][86][87][88][89] resembling phase transitions [90][91][92]. Direct measurements of viscoelastic features such as yield stress and viscosity have been performed in embryonic tissues and spatial and/or temporal drastic changes have been observed [12,13]. In the case of the early zebrafish blastoderm, for example, tissue viscosity was found to abruptly drop by more than an order of magnitude within a few minutes at the onset of morphogenesis [13] (Figures 1D',D'''). In addition, comparison of the yield stress between two neighboring tissues along the zebrafish body axis, the presomitic mesoderm and the progenitor zone, has revealed the presence of solid-like and fluid-like tissues, respectively [12]. However, does cell viscoelasticity scale in such cases with tissue viscoelasticity? Although--to our knowledge--no simultaneous analysis has been performed yet on measuring both cell and tissue scale viscoelasticity under different conditions to quantitatively assess their relationship, several lines of evidence point at the hypothesis that tissue-scale viscoelasticity critically depends on the interaction patterns between cells. Along these lines, it was reported that inhibition of myosin cytoskeletal motors in zebrafish, that is expected to decrease cell-scale viscoelasticity, had no effect on tissue viscosity [29]. Similarly, pharmacological treatments of the zebrafish tailbud with blebbistatin (a pharmacological myosin II inhibitor), which lowers cytoskeletal elasticity, had surprisingly inverse results on the tissue-scale material properties, where the tissue yield stress in treated embryos it is almost double than the control embryos [12]. In contrast, changes in cell rearrangements, cell-cell adhesion, contractility and cell division rates were shown to trigger changes in the tissue scale material properties ( Figure 1D) [11-13, 86-88, 93]. This agrees with extensive theoretical work inferring tissue material phase transitions based on certain cell parameters such as cell motility and adhesion [27,28,32,41], and quantitatively linking tissue rigidity to cell-cell connectivity and adhesion [29], but not directly to the rigidity of individual cells.
Overall, experimental measurements of several viscoelastic characteristics in cells and tissues have revealed nonlinear changes at the macroscale, e.g., viscosity, yield-stress, elastic modulus that do not trivially match similar changes at the microscale, e.g., individual cell and filaments material properties, strongly supporting that viscoelasticity of living systems is an emergent property.

VISCOELASTICITY AS AN EMERGENT PROPERTY: THEORETICAL MODELS
Several theoretical frameworks have been developed to establish a micro-macro link that can explain cell and tissue viscoelasticity. Numerous models exist in describing complex viscoelastic behaviors in chemical polymers where we recommend to the reader for a more specialized relevant literature [94][95][96]. Here however, we will summarize experimentally-based models belonging to three categories, based on the abstraction used to represent the biological system: In cytoskeleton networks, the building blocks are the biopolymer filaments and their mechanical properties, in vertex models, the geometrical properties of individual cells tilling the tissue and finally, in the topological models, the local topological arrangements of e.g., cell-cell contacts.

Modeling Viscoelasticity in Cytoskeleton Networks
The theoretical modeling of filament networks must account for the particular rheological phenomena these systems show, such as the strain stiffening response, stiffening tunability, and recoverable network fluidization [8,15,26,30,33,71,77,[97][98][99]. These models consider, at the microscopic level, filament properties like stiffness or length, and the local geometric and topological patterns of cross-linking, which project to the macroscopic level as material properties. Given these parameters, qualitative shifts in the response of the network are expected while increasing the density of filaments: First, beyond a certain threshold of density, the phenomenon of geometric percolation is observed--usually referred to simply as percolation, purely based on the network topology (Supplementary Box S1). At higher filament density, another qualitative shift in the properties of the network is observed, when the stiffness percolation threshold is overcome (Figure 2A). Beyond this filament density threshold, any stress applied at any point of the network will propagate throughout the whole system, meaning that the Young modulus, E, transits from E 0 to E > 0 [26,100]. Stiffness percolation is related to rigidity percolation (Supplementary Box S1), where stresses are no longer absorbed locally, but rather globally. Rigidity percolation, however, only refers to the topological structure [48,101], whereas, in the case of filament networks, further parameters are considered, as discussed below.
A well accepted model describing the phase space of filament viscoelasticity defined by the parameters of filament density and individual filament properties is the Worm-like Chain (WLC) model [15,102]. In the WCL model, biopolymers are considered as elastic rods or fibers with finite resistance to bending. Geometrically, such fibers are depicted by an inextensible curve with an energy penalty for bending. Let r(s, t)denote the path the curve takes in space (and time), parameterized by s, the arc length along the curve and χthe bending modulus. The functional accounting for the energy costs of fiber deformation will be given by: The second derivative accounts for the local curvature of the filament. The bending modulus χ has units of energy times length. The above energy function is penalizing any increase of the filament curvature and χ gives us the scale of such energy penalties. Thermal fluctuations play an important role here, due to the microscopic size of the filaments. As it is standard in statistical mechanics [103], non-zero temperature regimes imply the presence of stochastic, brownian fluctuations whose scale is k B T, where T is the room temperature and k B is the Boltzmann constant. Knowing the scale of the local stochasticity due to temperature enables the definition of a characteristic length for the system, as: which corresponds to the persistence length l p (Supplementary Box S1). Together with the contour length l c (Supplementary Box S1), the filaments are classified according to their stiffness as flexible, semi-flexible and stiff. If l p /l c ≪ 1, the filament is flexible and thermal agitation can induce traverse fluctuations. In this regime, entropic elasticity (Supplementary Box S1) dominates the dynamics, such as during cytoskeleton stretching [30,71,75,104]. In the case where l p /l c ≫ 1 , filaments are considered stiff and no fluctuations induced by thermal agitation are allowed. In this regime, enthalpic elasticity (Supplementary Box S1) dominates the deformations of the filament, such as during cytoskeleton bending or buckling observed in branched actin networks in the cell lamellipodia [77,105]. Finally, when l p /l c ≈ 1, filaments are considered semi-flexible and transverse undulations due to thermal fluctuations are possible, although attenuated. Whereas both the flexible and rigid regime show linear response to strain, the response in the semiflexible regime is nonlinear. To characterize the phase space, we observe that the system can transit between different regimes by changing the filament density ρ (proportional to 1/l c ) and filament length L playing the role of l p (proportional to the chain's molecular weight). Taking these two parameters as the coordinates of a phase diagram, one can identify four regimes: Affine entropic, affine enthalpic, non-affine, solution [15,26,33] ( Figure 2B, Supplementary Box S1). In this phase diagram, the phase transition from a fluid regime to a rigid regime is found at the border between the solution and the nonaffine regime ( Figure 2B), whose functional shape can be approximated as: Experimental work on reconstituted F-actin networks, revealed that linear and nonlinear strain-stress relationships can be explained by entropic and enthalpic models, respectively. In the absence of crosslinkers, actin networks generally form weak elastic gels mimicking the elastic nature of the filaments [106,107]. At the entropic regime, elasticity comes from the resistance of each polymer/filament against stretching [33,107]. At the enthalpic regime, while increasing strain, filaments first bend, reorganize along the direction of shear strain and the network deformation arises from the enthalpic stretching of the aligned filaments [108][109][110]. Further work showed that by decreasing the concentration of cross-linkers the network transits from affine to non-affine [74,75]. In conclusion, the modeling of biopolymer networks within this framework has proven powerful enough to account for the special viscoelastic properties of these systems.

Cell-Based and Energy Minimization in Vertex models
When modelling tissue viscoelasticity, the fundamental units, cells, are considered to exhibit certain properties [25,[111][112][113][114][115][116][117] (Supplementary Box S1) arising from the cell cytoskeleton [118]. Parameters, such as departure from an ideal cell shape and active fluctuations condense the material response of the cell and its effect to the tissue architecture under stress. There are several abstractions, depending on where the emphasis is placed concerning the energetic cost of the tissue deformations or configurations. Here, we will briefly mention Cellular Potts models and Centroid models, and focus more on the Vertex models (Supplementary Box S1) which provide a widely applicable framework to the understanding of biological tissue properties. Cellular Potts models idealize the tissue architecture as a mesh in which each point can be in several states and, accordingly, can represent a part of the cell, a contact point, or a free space. Each state of the mesh point has a particular contribution to the overall energy of the system, and may depend critically on the state of its neighbors [116]. Cellular Potts models may be considered within the much broader family of network models. In turn, centroid models base the analysis on the assumption that cells can be represented by their centroid position within the geometry of the tissue. Energetic contributions are based on geometrical considerations between centroids, such as distance between them [117,119,120]. The source of energetic contributions can be considered somehow complementary or even opposed [120] to the one considered in vertex models and the energetic costs associated with different configurations can be associated with material properties of the modeled tissue.
Vertex models have recently received a lot of attention due to their potential in describing a wide range of tissue properties [25, 27, 36, 111-115, 117-119, 121-123]. Moreover, the accurate study of their mathematical properties revealed a plethora of interesting physical properties [27,28], such as second order rigidity [121], that could play a relevant role in the modeling of biological tissues. Importantly, these properties may imply an interesting departure from the framework proposed by inanimate material science, and proposes a theoretical framework that extends to biological materials. In principle, vertex models have been postulated to model confluent tissues, with recent approaches extending the framework to non-confluent tissues [32]. The essential problem vertex models try to answer is: what are the energetic costs for cell migration within a tissue? In that framework, fluid states are those by which the movement of cells can happen at almost no cost, and rigid states will correspond to those states by which moving a cell -in particular, performing a T1 transition -implies a positive energy penalty, to be paid either in the form of external work or by the cells themselves ( Figure 2C). Arguably, the response to an external stress will be, at least partially, driven by the possibility of cell rearrangements which, in turn, depend on the ease of movements of cells within the tissue. The energetic contributions that configure the overall energy of the tissue come from the resistance against compression and the departure from some preferred shape in cells which, in most cases, is introduced as the preferred relation perimeter/area in 2D projections [27,113,121,122]. Having N cells, the functional accounting for the energy of the (2D) system reads: where A j and A 0j are the actual and preferred areas of cell j; P j and P 0j are the actual and preferred perimeter of cell j. K Aj and K Pj are the area and perimeter moduli, respectively. The first term models volume incompressibility and the second term models the active contractility of the actomyosin subcellular cortex. As shown, the underlying complex properties of the filament network which are themselves the outcome of a multidimensional problem have been absorbed by the scalar parameters K Aj and K Pj . The above equation can be non-dimensionalized in length if we divide it by A 0 √ , resulting in an effective shape index s 0 : s 0 P 0 A 0 √ s 0 will act as the control parameter in the different phase transitions between the material regimes of the tissue ( Figure 2C). The study of the energetic costs of motility within the tissue is performed through the analysis of the energy barriers arising when a T1 transition occurs (Supplementary Box S1). The energy barrier corresponding to a T1 transition can, in consequence, be computed as the energy difference, based on Eq. 1, associated to the reduction of the length of an edge to 0 [124,125]: If a T1 transition can be performed at no energy cost, the tissue is in a fluid configuration. Otherwise, the tissue is considered rigid. Remarkably, one can observe a well-defined phase transition from ΔH 0to ΔH > 0 as a function of s 0 . In particular, for values s 0 < 3.81, the energy imbalance is no longer zero, meaning that T1 transitions happen at a finite energy cost. Interestingly, in disordered systems, the shear modulus vanishes at s 0 < 3.81 [25,27]. Given the composition of tissues, any deformation must involve rearrangement of cells and, therefore T1 transitions. Consistently, one would expect that the shear modulus goes hand in hand with the energy barriers. Remarkably, this is not the case, and for values 3.71 < s 0 < 3.81, the tissue has zero linear response but non-zero high-order response [121]. In consequence, this phenomenon has been called "second order rigidity" and the material regime has been named "soft solid" regime [126]. In spite of the numerical evidence that the emergence in the linear regime of the shear modulus coincides with the emergence of theoretical energy barriers in disordered systems, it is currently still unclear what is the theoretical background that could explain the emergence of viscoelastic properties from the theoretical framework of vertex models. Finally, an extension of this vertex model has been developed for non-confluent tissues where stochastic fluctuations in cell surface tensions, density and cell rearrangements control rigidto-fluid transitions [32]. This contribution detaches the framework of vertex models from the structural constraints that confluency imposes, hence broadening their application in biological tissues.

Topological Models Based on Cell Contact Networks
The presence of interstitial fluid (Supplementary Box S1) in nonconfluent tissues, such as embryonic tissues and tumors [12,13,127], opens the possibility to apply even simpler theoretical frameworks to study their viscoelasticity [29]. The reason is that abandoning confluence liberates the system from a lot of implicit constraints at the structural level. Within nonconfluent tissues, the range of potential structural patterns increases enormously and topological models can exploit this potential heterogeneity. Here we discuss how topological models can be connected to the rheological properties of non-confluent tissues.
Topological models consider only the structural pattern of connections as the source of the material properties of the tissues. When the system is abstracted at the topological level, its structure is represented by a network defined by nodes and connections among them (Figure 2A). It is important to stress that, at first approximation, no other component, such as link properties or geometric embedding, is considered. The topological analysis, therefore, distills the structure of the system at the level of microscopic minimal components and combinatorial relations among them. In that sense, the basic observables of these models are, for example, the number of links connecting a given node to other nodes of the network, or the existence of paths, within the graph, between a given pair of them. In general, the approximation of a random network, in which the number of connections per node fluctuates stochastically according to some general constraint, quite accurately describes the behavior of real systems -given a suitable choice of constraints [92,[128][129][130][131]. For example, in the case that a network is representing the contacts of cells in epithelial (2D) tissues, it must belong to the class of planar networks, namely, those networks that can be extended in a 2D surface without displaying any link overlap [131]. In spite of the apparent simplicity of the approach, the study of networks at the topological level displays a wide range of non-linear phenomena, such as phase transitions or self-organized criticality [92,[129][130][131]. Particularly relevant is the phenomenon of percolation [91,92,128,129], briefly mentioned in previous sections (Supplementary Box S1). We will focus on the emergence of the so-called rigid cluster percolation (Supplementary Box S1), due to its important implications in cells and tissue material properties. Rigid cluster percolation is based on generic rigidity theory. Given a graph of N nodes and N L links, a graph is rigid if none of its nodes can be moved independently without constraining or stretching a link ( Figure 3A). Despite that the informal definition provided above for rigidity percolation appeals to material deformations e.g., stretching, compressing links, it turns out that the identification of rigid regions in a graph is a purely topological problem: it relies on the identification of actual degrees of freedom remaining in the network through application of the pebble game algorithm [47] based on a theorem considering only the topology of the network demonstrated by [101] (see Figure 3A). A natural question arises: what are the conditions leading to rigidity in a network? The answer is based on Maxwell's constraint counting [48,132], where very large networks made of nodes N and N L links are considered, the links acting as pairwise constraints, limiting the possibilities of independent motion of the nodes connected by the link. For simplicity, a random triangular lattice is considered (e.g., networks in Figure 2A) where the probability that a link exists is p e . If p e 1, all the links are present and the average connectivity 〈c〉, the average number of links per node, is 〈c M 〉 6 --recall we are considering very large networks where the boundary effects are negligible. In consequence, the number of links of this network will be fairly approximated by: A network will be rigid if the constraints absorb all the degrees of freedom for the motion of nodes ( Figure 3A). Considering a triangular lattice embedded in a 2D plane, each node has, a priori, d 2 degrees of freedom (position x and y). The number of degrees of freedom remaining in the network or Floppy modes, F, will be approximated by: Note that here we consider that all links are responsible for an independent constraint. We must notice that this is an important simplification [51,132] of the problem, but has proven useful in terms of both simplicity and predictive power. Starting from p e 0, the number of floppy modes is expected to decrease as long as p e grows. The key question is to identify the p e by which F 0 and, thus, the whole network is expected to be rigid. This is known as the isostatic point (Supplementary Box S1) and, using the expression for E found in Eq. 2, and setting F 0, we obtain: N〈c M 〉 that lead us to: which, in terms of average network connectivity, implies 〈c〉 4. That is, if the probability of link existence is p e > 2 / 3 , one expects the network to be rigid, and no independent movements could, in principle, be performed without imposing work over the system. Nevertheless, the emergence of rigidity is a much more complex phenomenon. In the case that p e < 2 / 3 , for example, there are constraints already acting in the system, so one would expect to see rigid regions within the network. At the same time, the probabilistic nature of the reasoning for finding the isostatic point may induce one to think that some regions of the network may remain floppy even for p e > 2 / 3 The answer is that, in very large systems, the relative size of rigid regions at p e < 2 / 3 , is negligible and that, for p e > 2 / 3 one observes the emergence of the giant rigid cluster (GRC) spanning almost all the network (Figure 2A). The emergence of the GRC is abrupt and has all the features of high order phase transitions. Therefore, the isostatic point p e 2 / 3 is the critical point of a phase transition called generic rigidity percolation, or simply, rigidity percolation [47,51].
In order to connect rigidity and topology with the material properties, the mechanics of the links, considered as springs, should also be considered. In that context, p e < 2 / 3 (subcritical regime) implies that one can perform a differential deformation over the network at no cost. On the contrary, if p e > 2 / 3 , any deformation implies an energetic cost in the form of external work performed over the network, as some spring will have to be unavoidably stretched or compressed. Therefore, in the case of very large systems, the first natural consequence of the emergence of the rigid cluster is that the Young modulus E will be E 0 in the subcritical regime and E > 0 in the supercritical regime. In consequence, topological phase transition that results into the emergence of the GRC projects into the material properties, implying a qualitative shift in the material response of the system under deformations. How does it project specifically into the viscoelastic behavior of the network? To understand that, a minimal ingredient of viscoelastic behavior should be introduced within the springs. This is performed by considering that each spring may update its rest length at random [29,133]. Specifically, in the simplest approach, the spring updates its rest length at random at every time unit step with probability p 1 τ . Formally, at time t the rest length is l(t) and the actual length, due to some external stress, is l(t) + Δl(t), p(l(t + 1) → l(t) + Δl(t)) 1 τ In Supplementary Box S2 we show how the microscopic dynamics of energy dissipation gives rise to macroscopic viscosity. In particular, if E is the Young modulus of the network, each link has the same spring constant k, and all such springs update their rest length at random with probability p 1 τ at every time step, one is led to: where γ, is the strain and σ, the applied stress ( Figure 1B). Eq. 4 is a constitutive equation for a Maxwellian viscoelastic material [14] ( Figure 1B). The above equation enables us to identify: as the viscosity of the system. What is relevant here is that we have a direct relation between the Young's modulus of the system, E, and the viscosity of the material, η, up to a constant that is the average lifetime of springs: The faster the update of the spring rest length is, the more fluid the behavior of the material is. On the contrary, in the limit of no updating, the material is only elastic. Moreover, it is established above that the emergence of a finite Young's modulus depends critically on the rigidity regime of the network. In conclusion, non-zero viscosity will emerge as the consequence of the emergence of the GRC, at least at the linear level. The phase transition observed for the emergence of the rigid cluster must therefore leave a footprint in the viscoelastic behavior of the network at the critical point. The above framework is able to bridge, in spite of its simplicity, microscopic topological patterns to macro-structural properties [47,51], without a priori reference to the mechanical properties of the constituents. These patterns, in turn, can be formally mapped to macroscopic material observables, such as viscoelasticity, hence demonstrating a potentially widely applicable framework to quantitatively link microscopic interactions to macroscopic viscoelasticity across scales.

RIGIDITY PERCOLATION PROBING VISCOELASTICITY ACROSS SCALES
Recent experimental work indicates that analysis at the purely topological level has the potential to indeed probe viscoelastic properties across scales. These models have been applied to several forms of biological networks, such as the cytoskeleton, ECM and cellular networks-tissues, to probe their material properties [29,31,44].
In the case of the cytoskeleton, both experimental and theoretical studies have shown that network connectivity is an essential parameter for cytoskeletal network mechanics [26]. If the actomyosin cytoskeleton is considered as a passive system, then rigidity percolation can predict the elastic modulus of the system based solely on its connectivity. In fact, a phase transition is proposed to occur in such networks, where mechanical rigidity emerges at the isostatic point -see Eq. 3 [31] (Figures 3B,B'). Such a framework can be expanded to active systems, where network connectivity together with motor activity can be further used as parameters to predict the contractile behavior of the actin cytoskeleton. Experiments on actin gels, where connectivity is regulated by the density of fascin crosslinkers and the motor activity by the density of myosin, showed characteristics of a rigidity percolation transition [43]. Briefly, weakly crosslinked networks (low connectivity) showed local contractions, medially crosslinked networks (higher connectivity) formed distinct contractile clusters within the network with a certain rigid cluster size distribution and, strongly crosslinked networks (highest connectivity) exhibited global network contraction associated with network fracture [43] ( Figure 3B). The authors further propose that the motors have the ability to reduce connectivity via forcing the crosslinkers to unbind, in order to avoid network fracture and thus the interplay of motor activity and crosslinking drives active gels to a critically connected state that can balance between local and global contractions [43].
Rigidity properties of ECM networks, and in particular type I collagen fiber networks, have also been well-described by rigidity percolation theory [134]. Studies combining experiments and theory suggest that the shear modulus of collagen fibers shows a strong correlation with the collagen volume fraction, and that these networks display connectivity near the percolation threshold [134][135][136][137][138]. Further experimental work however, has revealed that collagen networks with connectivity slightly below the isostatic threshold, can also become rigid in the presence of large deformation instead, thus in such cases passive rigidity percolation may not be sufficient to explain ECM rigidity [44,85]. In particular, increasing shear deformation in sub-isostatic networks leads to nonlinear increase of the elastic modulus of such networks along different connectivity values, an observation highlighting the possibility of incorporating the active nature of such systems when applying rigidity percolation theory [139].
Recently, the concept of rigidity percolation has been applied in non-confluent embryonic tissues to map tissue rigidity/ viscosity (Supplementary Box S2). Although tissues do not form physically crosslinked networks as the cytoskeleton or ECM, they can be approached as "cellular networks", where the nodes are the cells and the connecting links the adherensjunctions ( Figure 3C) [140]. This theory was applied to the zebrafish blastoderm which undergoes an abrupt and dramatic loss in viscosity at the onset of morphogenesis [13,29]. These changes in blastoderm viscosity were probed via rigidity percolation analysis over cell contact networks of the blastoderm. The size of the GRC was analyzed as a function of connectivity, and it was found that the GRC size correlates with the experimentally observed changes in tissue viscosity [29] ( Figure 3C). The emergence and disappearance of the GRC around the critical point matched the empirical observations where embryos whose cell contact network displays an average connectivity below the critical point, display a small GRC and are fluidized, whereas embryos whose network displays an average connectivity above the critical point display a big GRC and are rigid ( Figure 3C). This work further traced hallmarks of phase transitions, such as the diverge of macroscopic observables and its critical exponents at criticality, showing that rigidity percolation theory can be applied in embryonic tissues in vivo to link macroscopic tissue rigidity to the microscopic cell connectivity of these tissues [13,29].

DISCUSSION AND OUTLOOK
Living cells and tissues behave like viscoelastic materials [20,141], a long-standing observation that has only recently been linked to cell and tissue physiology [2][3][4][5]. Spatiotemporal regulation of viscoelasticity has been shown to influence essential biological processes, such as cell motility, proliferation, wound healing and the morphogenetic processes of body axis elongation and tissue spreading during embryonic development [11-13, 86, 88, 142-144]. Tracing the microscopic regulators of viscoelasticity is, however, a challenging task: typically, the mechanical properties of these microscopic building blocks do not match trivially the emerging viscoelastic behavior of cells and tissues at the macroscopic scale. Among many examples, we find the nonlinear increase of the elastic modulus of cytoskeletal and fiber networks in response to strain, or abrupt drastic changes in tissue viscosity without associated mechanical changes at the cellular level [12,29,30,33]. Given that the mechanical properties of tissues are regulated at the microscopic level, e.g., from the properties of the microscopic building blocks and their interactions, quantitatively bridging the microscale to macroscale is fundamental in order to understand the emergence of viscoelasticity [145]. Several theoretical approaches shed light on the biological mechanisms by which viscoelasticity can emerge in a system. However, the application of such theories is still far from comprehensive, given several challenges -such as the active, non-equilibrium nature of living systems.
At the cellular scale, the viscoelasticity of networked biopolymer filaments forming the cytoskeleton and ECM is most frequently modeled based on the mechanics, geometric alignment and local topology of the biopolymer fibers. It is worth mentioning that most of the research activity was performed in networks composed of actin or microtubules or intermediate filaments. In fact, biopolymer networks show a much wider heterogeneity, since the cytoskeleton is a dynamic mixture of all the cytoskeletal elements, interacting with each other, and these interactions were shown to affect their mechanical properties. Percolating networks of actin and vimentin for example, display synergistic effects in the elastic modulus, which becomes much greater than the sum of the elastic moduli of the two networks alone [146]. Another synergistic effect on actin and IFs networks is the recently reported phenomenon of super-elasticity, observed during epithelial morphogenesis [147]. Similarly, actin and MTs composite networks were shown to exhibit reinforcement against compression and mechanical enhancement [148][149][150], and a similar phenomenon was also observed when compressing composite ECM networks of collagen and hyaluronan [151,152]. Future experimental and theoretical work on composite networks is expected to provide a more complete understanding of cell and ECM viscoelasticity.
At the tissue scale, the vertex models are widely used to represent tissue viscoelasticity and have the special advantage of collapsing the emergence of macroscopic properties on a single cellular parameter, such as cell shape, and its associated scalar parameters, such as compressibility or departure from preferred area/perimeter [27,28]. Since, in a tilling of cells, deformation can only arise through cell-cell rearrangements, vertex models mostly focus the analysis on the study of the energy penalties associated with these cell-cell rearrangements. Such rearrangements can for example be induced by differences in cell-cell adhesion and cortical tension or active tension fluctuations [27,28,32,125]. It is implicitly postulated that the results at that level of abstraction will project into the material properties of the system. Even though numerical approaches are coherent with the predictions of the models, a direct bridge between micro and macro scales in tissue viscoelasticity has yet to be clearly described. The recent application of rigidity percolation in tissues however, provides a quantitative link between the topological patterns of cell-cell contacts and tissue rigidity regime (quantified by the size of the GRC). In the case of the zebrafish blastoderm, the topological approach of rigidity percolation was sufficient to capture the floppy and rigid regimes of the tissue by one single (microscopic) measurable parameter, the average number of cell-cell contacts in different cell types and, as a result, match the observed (macroscopic) changes in viscoelasticity [29]. Cell connectivity was further shown to be defined by the biophysics of cell-cell contacts and specifically the cell-cell interfacial tension at the contact. In this biological context, experimental testing of the phase transition parameters revealed that changes in connectivity and cell-cell adhesion were driven by cell division, and not by cell rearrangements, cell shape or active tension fluctuations [29]. Since cell-cell adhesion is key in tissue rigidity theoretical frameworks so far, we speculate that some parameters may be common, such as cell-cell adhesion strength [29,32,118,123], and some others may be context-specific such as cell division and rearrangements. Future experimental work has the potential to disentangle the physiological role of several cellular parameters in rigidity transitions.
In all of the above models, incorporating detailed dynamic analyses of the microscopic parameters that can account for local heterogeneities, such as in adhesion strength (between the cells and with their environment), shape differences or heterogeneous motility patterns, will increase their potential to model absolute viscosity values and provide a more accurate and representative image of tissue viscoelasticity. A key challenge is the choice of the viscoelastic model relevant to the biological system. Here, we have extensively described how Maxwell viscoelasticity has the potential to be linked to maxwell rigidity. However, other biological systems may be better described with different models. For example, the Kelvin model was recently used to describe the phenomenon of arrested coalescence in multicellular aggregates from the adherent and contractile protrusion interactions between the cells [153]. In addition, both Maxwell and Kelvin viscoelastic models can describe different viscoelastic regimes during density/packing dependent collective cell migration [154]. Another important challenge is to understand if and how the timescale of a biological process taking place at the microscale is relevant to the timescale of a biological process taking place at the macroscale when bridging scales in viscoelasticity. For example, how macromolecular motion within the cytoskeleton network (milliseconds to seconds) influences cell shape changes driven by the mechanics of the cytoskeleton (seconds to minutes), or how cellular motion within a tissue (minutes) influences tissuescale fluidization (minutes) and spreading (hours)? Several theoretical and experimental frameworks should be developed do bridge scales in time and space [142,155].
Disentangling the connection between microscale behavior and emerging, macroscopic properties is not a novel goal: it has been the long-sought target of statistical mechanics in order to connect thermodynamics to a solid microscopic basis [156]. Moving towards a broader conception of statistical mechanics encompassing the living phenomena requires the introduction of the microscopic role of the biological building blocks -which are far more complex than gas particles, for example. Nevertheless, establishing the critical point in the microscopic dynamics of the building blocks that would trigger a macroscopic phase transition would create a rich toolbox for biology, regardless of the theoretical approach being used. Phase diagrams or morphospaces, accounting for what is possible in the relations between tissue organization and material properties, can be defined that will allow further exploration for the different regimes the system can occupy and their grounding. In the case of the zebrafish blastoderm, for example, it was experimentally possible to position the system in the vicinity of criticality. Hallmarks of criticality, such as divergence of macroscopic observables with associated power laws, were determined experimentally in the living embryo, showing that tissue morphogenesis in vertebrates may start close to a critical point of a rigidity transition. This indicates that embryonic tissues may be at optimal fitness [157] since they are able to easily switch between rigid and floppy regimes by slightly changing their connectivity at almost zero energetic cost [29]. Through this lens, one can explore fundamental questions concerning morphogenesis. How do local heterogeneities in the microscopic parameters influence the emergence of macroscopic viscoelasticity? How do noisy biological systems [158][159][160] guarantee stability during development when they are poised at criticality?
Networked systems and criticality have been for long used to understand an extremely rich palette of macroscopic phenomena occurring in the natural world [92,131,161,162], that now include the viscoelastic characteristics of biological systems, from the nanometer scale of the cytoskeleton to the micrometer scale of tissues and embryos. Such observations indicate that an efficient organizing strategy of complex biological systems may be to behave as networked systems close to criticality.

AUTHOR CONTRIBUTIONS
BC-M and NP designed, drafted, revised and finalised the work.