Collective durotaxis of cohesive cell clusters on a stiffness gradient

Abstract Many types of motile cells perform durotaxis, namely directed migration following gradients of substrate stiffness. Recent experiments have revealed that cell monolayers can migrate toward stiffer regions even when individual cells do not—a phenomenon known as collective durotaxis. Here, we address the spontaneous motion of finite cohesive cell monolayers on a stiffness gradient. We theoretically analyze a continuum active polar fluid model that has been tested in recent wetting assays of epithelial tissues and includes two types of active forces (cell–substrate traction and cell–cell contractility). The competition between the two active forces determines whether a cell monolayer spreads or contracts. Here, we show that this model generically predicts collective durotaxis, and that it features a variety of dynamical regimes as a result of the interplay between the spreading state and the global propagation, including sequential contraction and spreading of the monolayer as it moves toward higher stiffness. We solve the model exactly in some relevant cases, which provides both physical insights into the mechanisms of tissue durotaxis and spreading as well as a variety of predictions that could guide the design of future experiments. Graphic Abstract


Introduction
The organized motion of cohesive groups of cells, usually referred to as collective cell migration, plays a key role in many instances of morphogenesis, tissue regeneration, and cancer invasion [1][2][3][4][5][6]. The mechanisms by which cells coordinate their motion are diverse and often not fully understood. Recent work has shown that groups of cells may respond to external stimuli as a whole, that is, in the form of collectively organized directed motion, in ways similar to what single cells do. Such collective migration can arise in response to a variety of external stimuli such as gradients in either chemical concentrations or in the stiffness of the environment, which, respectively, lead to collective chemotaxis [7] and durotaxis.
We are interested in the phenomenon of durotaxis, which refers to the directed motion of cells along stiffness gradients of the extracellular matrix, typically toward stiffer regions. This is a well-known phenomenon for single-cell migration [8], which is rather common in many types of cells and has important implications for cancer invasion. More recently, durotaxis has been reported also for collective cell migration [9,10]. a e-mail: jaume.casademunt@ub.edu (corresponding author) Remarkably, large cell monolayers can perform durotaxis collectively even when their constituent cells do not [9], and in some cases, there is an optimal intermediate stiffness for tissue spreading [11,12]. Collective durotaxis has been theoretically described both via hybrid computational models [13][14][15][16] and via a continuum active polar fluid model [17] that generalized previous work on tissue wetting [18]. This continuum model was solved numerically to reveal two possible mechanisms of collective durotaxis [17].
Here, we extend the work in Ref. [17] to provide a more comprehensive classification of the dynamical regimes of the model in terms of physical parameters. Remarkably, we solve the model analytically in some simple but relevant situations, allowing for a better grasp of the physical mechanisms at play. As shown in Ref. [18], the model predictions can be fitted to experimental data to infer physical parameters that are often elusive to direct measurement.
The model describes cell monolayers moving on a substrate as a quasi-two-dimensional viscous fluid with two types of active forces: cell-substrate traction and cell-cell contractility. The competition between both active forces was shown to give rise to the socalled active wetting transition, whereby a tissue either spreads or retracts depending on its size [18]. The same model also predicted a fingering instability of the lead-7 Page 2 of 15 Eur. Phys. J. E (2022) 45 :7 ing edge of the tissue [19]. In addition to the active forces, the model also features two passive forces: an effective viscosity, which arises from cell-cell adhesion, and a friction force due to cell-substrate interactions. All these forces are treated in a coarse-grained way at the supracellular scale. The rationale of the approach is to identify the dynamical behaviors of cell monolayers that are of mechanical origin, explicitly excluding any signaling effects that cannot be encoded in the mechanical parameters of the model. To what extent such purely mechanical approach may succeed as a first step to account for the observed phenomenology is an interesting open question that might be settled by future experiments.

Hydrodynamic model
Our model stems from a hydrodynamic approach to cell tissues, a strategy that has proven useful when tissues are organized at a supracellular scale, such that information at the cellular scale is not relevant [20][21][22][23][24]. This is the case in many examples of collective cell migration, where coarse-grained fields such as velocity, cell density, and polarization are treated as smooth fields varying on scales larger than the cell size. Continuum field theories based on linear irreversible thermodynamics, often called active gels theories, were first devised to account for active matter at the cellular scale, such as the cytoskeleton [25][26][27][28], but have more recently been extended to multicellular scales [29].
The basic idea is that tissues can be modeled to some extent as continuous active materials, in such a way that the biological properties are encoded in a series of physical parameters, including passive ones such as viscosity or friction, and active ones such as contractility or traction. These parameters will in general be time and space dependent to account for the biological regulation of the cell properties and interactions. For instance, in a simple model for the spreading of epithelial monolayers [30], it was shown that their effective viscosity increases with time as they become thinner due to the spreading. This type of approach is useful to identify activity-driven hydrodynamic instabilities that can either be avoided or exploited by the biological regulation of parameters [19,31].

Assumptions and model equations
In this paper, we take the simplest possible model of an active fluid that combines active cell-substrate traction and cell-cell contractile forces. This model was introduced in Ref. [18] and was extended in Ref. [17] to account for substrates with non-uniform stiffness. The model is for a two-dimensional active fluid, which describes the quasi-two-dimensional cell monolayer with two continuous fields: the velocity v α and the polarity p α . The polarity is the orientational degree of freedom of the cells which arises from the polarization of its internal cytoskeletal structure and defines the direction along which traction forces are exerted. The tendency of cells to align with their neighbors is accounted for by an effective free energy of the form where K is an effective Frank constant that quantifies the energetic cost of polarity gradients [32]. The constant a is a restoring coefficient that is taken positive such that the unpolarized state (p = 0) is energetically favored in the bulk. We assume that the polarity follows a purely relaxational dynamics ∂ t p α ∝ −δF/δp α which is much faster than the temporal variations of the rest fields [18], such that we can take a quasistatic evolution (∂ t p α = 0) and hence δF/δp α = 0. Then, we have where L c ≡ K/a is the nematic length that characterizes spatial variations of the polarity field [17][18][19]. Since epithelial cells migrate toward free space, we enforce a boundary condition of maximum polarity |p| = 1 directed normally to the tissue edge. Then, L c defines the thickness of a polarization boundary layer near the tissue edge, such that polarity decays from p = 1 at the edge to p = 0 deep into the tissue. Neglecting inertia, the force balance equation is where σ αβ is the stress tensor of the monolayer and f α is the external force density due to the contact with the substrate. These quantities are directly related to the experimentally measured monolayer tension, σ αβ h, and traction stress, T α ≡ −f α h, with h the height of the monolayer [18]. We now take the constitutive equations for a compressible active polar fluid of the form [18,33] where η is the viscosity, ξ is the cell-substrate friction, ζ < 0 is the contractility, and ζ i > 0 is the contact active force (hereinafter referred to as the traction parameter), which accounts for the maximal traction stress T 0 ≡ hζ i exerted by polarized cells on the substrate. A summary of the symbols for the variables and parameters, together with their units and estimates can be found in Table 1.
We assume that, in our 2d description, the cell monolayer is compressible, with ∂ α v α = 0, because the inplane compression and expansion of the cell monolayer can be accommodated by changes in the monolayer height h. We assume that in-plane deformations do not amount to significant changes in pressure as the layer can deform in 3d and, hence, pressure gradients are neglected in front of the rest of the contributions in the . This approximation has been used and discussed for instance in [17][18][19]21,22,30,31]. Tissue growth driven by cell proliferation is also neglected.
In the explicit form of the stress tensor Eq. (4), we have also assumed, following Ref. [18], and in order to reduce the number of parameters and define the simplest possible equations, that the bulk viscosity isη = η, and that the isotropic contractility is given by ζ = ζ/2. Similarly, active stresses not associated with polarization are also neglected, that is,ζ ζ as defined in Eq. (S12) in Ref. [18]).
Furthermore, the profile of polarity p α is dictated directly by the boundary shape, so flow alignment and other elastic effects, which have been addressed in more general models such as in Ref. [31], are here neglected. In the simplest formulation, we assume stress-free boundary conditions, but we also generalize the model to other cases.

Reduction to a 1d solvable model
The problem at hand is formally a free-boundary problem, since the boundary of the cell cluster is free to deform and move, as its normal velocity coincides with that of the adjacent fluid. The evolution of the shape and position of the boundary is thus part of the solution of the problem. An example of how a spontaneous symmetry-breaking of the morphology of the boundary can couple to the overall motility of the domain was discussed in the context of cell fragments [42]. In the present study, we are interested in cases where the symmetry is broken by the existence of an external gradient. Since this is the dominant effect causing motion of the domain, here we ignore boundary deformations. We may also assume that the effective surface tension of the tissue is strong enough, and the monolayers small enough, to suppress the active fingering instability that is inherent to this model, as reported in Ref. [18,19]. Our interest is thus to describe the motion of circular monolayers of radius R on a substrate with a stiffness gradient by tracking the position of the center of mass and the monolayer size.
For simplicity, and in order to obtain exact solutions and physical insights, we formulate the problem in a 1d setup, in which the monolayers are strips that are finite in the spreading direction, and infinite in the transverse direction ( Fig. 1, bottom). This setup corresponds to the experiments on collective durotaxis of Ref. [9] and was used also in the numerical study of the present model in Ref. [17]. The present work extends that previous study with a more comprehensive discussion of the wealth of dynamical behaviors allowed by this model and their physical interpretation, in particular taking advantage of explicit analytical solutions.
The basic physics of this 1d formulation (rectangular geometry with translational invariance in the transverse direction) is the same as that in the 2d case with circular monolayers (circular geometry with rotational invariance). The results are equivalent up to geometrical factors, but much simpler in the rectangular geometry, as already illustrated in the preceding studies in both geometries [17][18][19]. Furthermore, the availability of an exactly solvable model with sufficiently simple analytical results is of great theoretical value to gain insights into the physical mechanisms at play, in particular when a relatively large number of parameters are present. Moreover, we will show that some of the limitations of the 1d formulation, such as the lack of the Young-Laplace pressure drop due to tissue surface tension, can be effectively introduced in a simple way into the 1d reduction of the problem.
In the 1d setup, Eqs. (3)(4)(5) reduce to The polarity profile is given by the solution of Eq. (2) satisfying p = ±1 at the two edges x = x + and x = x − < x + , respectively. In terms of the center-of-mass position X ≡ (x + +x − )/2 and the monolayer half-width L ≡ (x + − x − )/2, it reads There are several length scales whose ratios determine different physical scenarios in the model. The scale L c is typically the smallest one, as the polarized boundary layer of the tissue is often thin compared to the system size L and the other length scales [9,18,30]. The so-called screening length λ ≡ 2η/ξ is a measure of the range of hydrodynamic interactions [17,19,27,31], and it defines two important limits: In the so-called wet limit, when λ L, long-ranged hydrodynamic interactions produce non-local effects and the system behaves globally as a whole; in the so-called dry limit, when λ L, the spreading dynamics is governed by local forces, i.e., the two edges behave independently from each other. Another relevant length scale, that we call the active polar length L p [18], arises as the ratio of contractility to traction forces: L p ≡ |ζ|/(2ζ i ). In the wet case, this length defines the critical tissue size for the wetting-dewetting transition, as reported in Ref. [18].
Equation (6) will be solved typically with stressfree boundary conditions. If a normal stress component is required to mimic the effect of an effective surface tension, as if L would be the monolayer radius, we will impose The solution of Eq. (6) provides the spatial velocity profile v(x), from which we obtain v ± as well as the velocity of the center of mass U ≡Ẋ and the spreading velocity V ≡L.
For tissues on a substrate with variable stiffness, the parameters of the passive and active forces on the surface, that is friction and active traction, will be space dependent, ξ(x) and ζ i (x). The relationship between these spatial variations and that of the substrate stiffness must be determined independently of the hydrodynamic model. An explicit derivation requires a detailed knowledge of the molecular mechanisms at play, and a discussion based on empirical data was made for instance in Ref. [17]. Both friction and traction parameters increase with and eventually saturate with increasing substrate stiffness [43][44][45][46][47][48]. However, to avoid introducing more parameters and to make the interpretation of the results more transparent, we mostly consider cases where those parameters are either space independent or have a uniform gradient, hence introducing only two new parameters associated with the stiffness variation, namely ξ ≡ ∂ x ξ(x) and ζ i ≡ ∂ x ζ i (x). This restriction is relaxed in Sect. 3.2. In most cases we take ξ = 0 (uniform friction) and focus on the effect of a uniform traction gradient ζ i > 0 on the net displacement of the monolayer. We then find the velocity profile at any given time by solving Eq. (6) with the initial conditions L 0 ≡ L(0) and X 0 ≡ X(0) for the set of parameters is the initial traction offset.  Table 1, except for Lc = 5 µm (smaller to see better convergence). Only for the largest λ, the critical size L * ≈ 200 µm = Lp approaches the wet limit prediction; for the other two values of λ, the dry approximation is better for smaller sizes

Solutions for a uniform substrate
We first consider as a reference the case with no stiffness gradient, so that ζ i = 0, ξ = 0, and consequently there is no net monolayer displacement: U = 0. This case was studied in the wet limit, ξ → 0 in Ref. [18] in circular geometry, and in the wet-dry crossover and in rectangular geometry in Ref. [19]. The exact solution of this case is given in Appendix A. Taking γ = 0 and assuming L c L, the expression in the wet limit λ L for the spreading velocity takes the simple form which recovers Eqs. (5) and (7) from Ref. [19]. In the dry limit L c λ L, we obtain In the wet limit, there is a critical tissue size L * ≈ L p that defines the so-called active wetting transition of Ref. [18]. This transition distinguishes whether the cluster is expanding (positive spreading velocity V > 0 for L > L * ) or contracting (negative spreading velocity V < 0 for L < L * ). The wet limit is represented by the dashed line in Fig. 2. In this limit, the spreading velocity V does not depend on λ, and the transition from contraction to expansion takes place at L = L p . The condition V = 0 thus defines the wetting-dewetting transition reported in Ref. [18]. However, here we refer to spreading and avoid the term 'wetting,' which usually refers to the local motion of a fluid front on a substrate. This precision is meant to avoid confusion in  Table 1 except for λ, which takes values λ = 40, 100, 200, 300, 450 µm cases where the center of mass of the tissue is moving. In those cases, the soft tissue edge may recede with respect to the substrate while the tissue globally expands. We discuss such examples in Sect. 3.1.
In the dry limit, the spreading transition is controlled by the screening length λ. Equation (9) shows that there is a critical λ * ≈ L p such that for λ < λ * the cluster contracts (V < 0), regardless of its size L, and for λ > λ * the cluster always expands (V > 0). This result in the dry limit is represented by the dotted lines in Fig. 2, which do not depend on the tissue size L and exhibit the spreading transition at V = 0 for λ = L p .
The full velocity and stress profiles (Fig. 3) are not used here, but they allow the model predictions to be tested against experimental data. These profiles provide a simple visualization of where the system stands in the wet-dry axis, and of the forces and flows in the cell monolayer. For example, the stress plateau in the bulk (darkest curves in Fig. 3b) is a signature of the wet limit (large λ). Respectively, two peaks of width L c near the edges (lighter curves in Fig. 3b) are indicative of the dry limit (small λ). In this case, the velocity profiles feature a plateau of null velocity in the bulk (lightest curve in Fig. 3a). In this situation with a uniform substrate stiffness, the profiles of stress and velocity are, respectively, even and odd with respect to the center of the monolayer.

Linear traction profile
The presence of a stiffness gradient should in general affect the interactions between the cells and the substrate, thus altering both traction and friction forces. The dependence of the traction and friction parameters on substrate stiffness must be determined independently of the model, either empirically or from a microscopic model of cell-substrate interactions. In this section, to obtain analytical solutions, we take the simplest possible spatial dependence of these parameters: a linear profile of traction , and a uniform friction coefficient, with ξ = 0. The corresponding results will be applicable locally to more general traction profiles as long as ζ i L/ζ i 1. The exact results for a linear traction profile, together with some approximate expressions, are given in Appendix B. An important exact result for this case of uniform traction gradient ζ i is that the spreading velocity V is that on a uniform substrate V u with the traction evaluated at the monolayer center, that is V (ζ 0 i , ζ i ) = V u (ζ 0 i ). Therefore, the spreading behavior is independent of the existence of a traction gradient. More generally, in cases where the traction gradient is not quite uniform, the spreading velocity will be relatively insensitive to that gradient. The velocity of the center of mass, however, is sensitive to the existence of a traction gradient, which gives rise to the phenomenon of durotaxis.
Next, we discuss the results in the dry and wet limits. Taking γ = 0, the expression of the edge velocities in the dry limit where ζ ± i are the local values of the traction at the edges. The corresponding center-of-mass velocity, neglecting 2L 2 c in front of Lλ, reads and the spreading velocity is V dry = V u,dry (ζ 0 i ), with . Although L appears in Eq. (11), giving a linear increase in U with L (dotted lines in Fig. 4), U can be rewritten in terms of the traction difference emphasizing that the spreading dynamics is local in the sense that the two edges behave independently from the other. The traction difference then directly drives tissue durotaxis.
On the contrary, in the wet limit 1 L c L λ, the two edges are coupled through hydrodynamic interactions, and we get which yields a center-of-mass velocity and a spreading velocity V wet = V u,wet (ζ 0 i ). Both v ± and U depend on the system size L and the traction 1 The strict wet limit λ → ∞ (ξ → 0) with finite Lc is illdefined for a the case of a nonzero ζ i , because force balance cannot be globally satisfied unless there is friction.  Table 1, except for Lc = 5 µm gradient ζ i , which illustrates that the two edges are hydrodynamically coupled. We provide a summary of results for tissue durotaxis U and spreading V in the wet and dry limits in Table 2.
Two main conclusions emerge which are general in the whole wet/dry range for this case of uniform traction gradient ζ i and uniform friction (ξ = 0). On the one hand, the center-of-mass velocity U is proportional to the traction gradient ζ i and independent of the traction offset ζ 0 i . U has the same sign as ζ i , and there is durotaxis to stiffer regions as long as the traction is a monotonically increasing function of the stiffness. On the other hand, the spreading velocity depends on the traction offset and not on the traction gradient. Accordingly, Fig. 2 still applies in the present case, and durotaxis is independent of whether the monolayer is spreading or contracting (Fig. 4).
In fact, the following situations are possible. First, the monolayer can contract either with the two edges moving in opposite directions (v − > 0 and v + < 0) or in the same direction (0 < v + < v − ). In the former case, both edges are retracting, or dewetting. In the latter case, the + edge is wetting and the − edge is dewetting. Second, the monolayer can expand, or spread, if both edges move away from each other (v − < 0 and v + > 0), both wetting the substrate, but also if both edges move in the same direction (0 < v − < v + ), with the + edge wetting and the − one dewetting.
It is thus clear that the condition of spreading or contraction, which is a property of the cell monolayer as a whole, and the condition of wetting or dewetting, which refers to the direction of motion of each tissue edge, are two distinct conditions that only coincide when the center of mass does not move (U = 0), as in Ref. [18]. We show examples of these distinct situations in Figs. 5 and 14 in Appendix C.
The repertoire of dynamical behaviors contained in the model as a function of parameters is quite rich. Table 2 Summary of the results. In the wet limit, the critical length is such that if L < L * ≈ Lp there is contraction (V < 0) whereas if L > L * there is expansion (V > 0). In the dry limit, both regimes are defined by λ * ≈ Lp

Uniform
Stiffness gradient (ζi, ξ unif.) (lin. ζi(x), unif. ξ) Spreading and center-of-mass velocities can be plotted against monolayer size (Fig. 6) and traction offset (Fig. 7), which are two quantities that can be easily varied and controlled in experiments [9]. Importantly, in addition to being independent of traction offset, the durotactic velocity U does not depend on the contractility either, which is a parameter that is more difficult to infer from experiments, and is assumed to be uniform throughout the system. An increase in either monolayer size L or traction gradient ζ i implies an increase in the difference of local tractions at the edges, ζ + i − ζ − i , and thus an increase in durotactic velocity U . The spreading velocity V , which is independent of the traction gradient ζ i , increases with the monolayer size L, the screening length λ, and the traction offset ζ 0 i , and it decreases with the contractility |ζ|.
The velocity and stress profiles, plotted in Fig. 8 for a range of λ, are qualitatively similar to those of the uniform-stiffness substrate (Fig. 3), except that they become asymmetric due to the stiffness gradient.

Linear friction and saturation profiles
In this section, we relax the restriction of a uniform friction coefficient. This is a more realistic situation since both active traction and passive friction rely on the dynamics of cell-substrate adhesion molecules [33], and hence they both depend on substrate stiffness. Previous works indeed support that cell-substrate friction increases with substrate stiffness [43][44][45][46]. To illustrate the role of this effect on tissue durotaxis, and for the sake of simplicity, we consider a linear friction increase ξ(x) = ξ 0 + ξ x, with ξ = 0. The problem with space dependent ξ can no longer be solved analytically. Solving Eq. (6) numerically with a finite-difference method, we find that the center-of-mass velocity now decreases with the traction offset (Fig. 9). This is because now larger traction correlates with larger friction, which leads to smaller velocities. Accordingly, the spreading velocity grows more slowly with traction offset than in the uniform-friction case.
As already mentioned, the use of linear profiles is particularly convenient from a theoretical point of view since it avoids introducing too many parameters. However, to obtain a more realistic description and to compare with experimental data, other profiles may be more appropriate. A simple feature that can be implemented in the model is the fact that the increase in traction and friction with substrate stiffness must eventually saturate. Then, following Ref. [17] and the references therein, we consider profiles of the form in terms of the spatially varying Young modulus E(x) of the substrate. For in vitro experiments such as those of Ref. [9], a simple choice is to prepare the substrate with a linear stiffness profile E = E 0 + E (x − X). Numerical results for this more general case are qualitatively similar to those in Fig. 9. However, at high stiffness, the saturation of traction and friction makes the tissue dynamics approach those of the uniform-stiffness case, with vanishing durotactic velocity U . The approach to this durotaxis-free regime at high stiffness is controlled by the new parameters ζ ∞ i , ξ ∞ , E * and E introduced in Eq. (14).

General case
For a given set of parameters, η, ζ, L c , initial conditions X 0 , L 0 , boundary conditions, and imposed profiles ζ i (x) and ξ(x), our model supplies a velocity profile v(x) as the solution of the equation where p(x) = p(x; X, L, L c ) is given by Eq. (7). The position of the center of mass X(t) and the cluster size L(t) satisfy the differential equationṡ  Table 1 Table 1 Table 1 and changing λ = 40, 100, 200, 300, 450 µm. Note that with the values of vx at the tissue edges we would obtain the spreading and durotactic velocities in Fig. 7b, for ζ 0 i = 0.05 kPa/µm , the traction offset changes with time according to ζ 0 i (t) = ζ 0 i (0) + ζ i (X(t) − X 0 ). In the rest of this section, we focus on this case with and take no friction gradient (ξ = 0). Fig. 9 Durotactic velocity U (a) and spreading velocity V (b) when there is a positive gradient of the friction coefficient, for parameters in Table 1, varying the friction gradient ξ = 0, 10 −4 , 5 · 10 −4 , 10 −3 , 3 · 10 −3 kPa·s/µm 3 , and taking a stiffness offset ξ0 = 2 kPa·s/µm 2

Uniform traction gradient
Since the durotactic motion is toward increasing traction (U > 0 for ζ i > 0), the local traction offset ζ 0 i increases with time. In a uniform traction gradient, the durotactic velocity U is insensitive to the local traction offset. Therefore, the increasing traction offset does not lead to an increasing durotactic speed. However, U depends also on the monolayer size L, which may grow or decay according to the sign of the spreading velocity V . In general, U increases monotonically with L. Therefore, if L is large enough, the monolayer spreads and then the center-of-mass velocity U increases during the evolution as L increases. As a conclusion, monolayer spreading produces increasingly faster durotaxis.
For small enough L, the monolayer initially contracts (V < 0). However, as the tissue moves toward stiffer regions, it may reach values of traction that are large enough to change the sign of V and produce a transition to spreading. In this case, the evolution of L is nonmonotonic in time, corresponding to initial contraction followed by spreading. Finally, if L is even smaller, the durotactic velocity U may not be sufficient to reach sufficiently large values of traction to reverse the sign of V to produce spreading. In this case, the monolayer contracts (dewets) completely into a three-dimensional spheroid.
The asymptotic behavior of the system at long times is thus either indefinite expansion or the collapse. In both situations, the model is no longer adequate as additional physics will take over at long times. In the case of asymptotic spreading, even if the traction forces do not saturate, other forces such as elastic forces may eventually slow down and even suppress the spreading as discussed below. In the case of monolayer retraction and collapse with L → 0, the quasi-two-dimensional description breaks down and a more elaborate treatment of the three-dimensional structure of the tissue becomes necessary.
For a given set of parameters, the tissue width L controls the transitions between the three possible spreading dynamics. For L ≥ L * , the monolayer expands for all times V (t) > 0. Here, L * is the critical size for spreading on uniform substrates (V (L * ) = 0), which we discussed before. An explicit and exact expression  Table 1. Here, L * = 276.35 µm and L * c ≈ 213 µm. The tissue contracts when the normalized L and U decrease and V < 0, whereas the tissue expands when L and U increase and V > 0. The regime with initial contraction and later expansion presents an almost constant durotactic velocity U and tissue width L. The corresponding edge velocities together with U and V in each case are shown in Fig. 15 for L * is given by equating Eqs. (A.3) and (B.10) to zero. For intermediate values of L, L * c < L < L * , the monolayer initially contracts (V < 0) and later expands (V > 0) indefinitely. Finally, for L ≤ L * c , the monolayer contracts for all times (V (t) < 0). The three regimes are illustrated in Fig. 10.

Some generalizations of the model
More general profiles of traction and friction can also be used for studying the time evolution. As long as the profiles are both monotonically increasing, the qualitative behavior is similar. The three spreading regimes discussed above, separated by the critical lengths L * and L * c , still exist, but their expressions and values change. For the case of uniform traction gradient ζ i and no friction gradient (ξ = 0), given the initial monolayer size L 0 , we predict the critical lengths as a function of the traction offset, as shown in Fig. 11.
As mentioned above, the two possible asymptotic behaviors of the monolayer dynamics are not particularly interesting. This is because the traction profile cannot grow indefinitely, and other physical effects will either stop the extreme stretching of the cells in the case of spreading or enable the formation of a  Table 1 three-dimensional cell aggregate in the case of contraction. The latter will not be considered here because it requires essential modifications of the model that are deferred to future work. However, different effects may be easily introduced in our current model either to slow down the indefinite spreading of the cluster size or even to stop it.
The first possibility is to introduce of an effective surface tension γ at the tissue edge, as already mentioned in Sect. 2.1. The introduction of this surface tension can be understood if we interpret our 1d model as an approximation for a circular monolayer of radius L. This surface tension slows down the spreading process, less effectively for larger monolayers. For contracting monolayers, surface tension accelerates the contraction.
A less trivial but more determinant modification is to introduce an effective elastic force that prevents excessive cell stretching. This type of force has been introduced at a phenomenological level for single cells to favor a characteristic cell size. It was used for instance in Refs. [49,50] in effective 1d models for single-cell motility. Such an elastic force, together with the Young-Laplace pressure drop due to surface tension γ, can be implemented in the following boundary condition for the stress: Here, k is an elastic constant, and L r is a characteristic size of the cell monolayer, which is proportional to the number of cells if the cell size is somehow regulated. For a uniform traction gradient ζ i and no friction gradient (ξ = 0), the center-of-mass velocity U turns out to be independent of both surface tension and elasticity (see Appendix B). The spreading velocity V , however, is affected, respectively, giving (b) (a) Fig. 12 Spreading velocity V as a function of tissue size L for the parameter values in Table 1 but with (a) surface tension γ = 0, 20, 40, 80 kPa·sµm (k = 0), and (b) elastic constant k = 0, 0.05, 0.1, 0.2 kPa (L r = 150 µm and γ = 0). Here, to showcase its effects, we take values of γ larger than what is measured experimentally for cell aggregates (Table  1). Respectively, we take k comparable to ζiLc ≈ σ when either surface tension or elasticity are added. The surface tension γ always decreases the spreading velocity, while the elastic term contributes with a different sign depending on whether the monolayer size is larger or smaller than L r , always in the direction of approaching the reference value L r (Fig. 12). Both effects affect the spreading dynamics, changing for instance the critical lengths, but the phenomenology and qualitative evolution of the monolayer typically remains unchanged. However, for large k and L > L r (Fig. 13a), a monolayer that starts spreading can change to contraction. In this case, similar to surface tension (Fig. 13c), elasticity slows down expansion and accelerates contraction. On the other hand, if L < L r (Fig. 13b), elasticity accelerates expansion and slows down contraction, although only very large values of k (k T L c ∼ σ), presumably not biologically possible, enable a transition from contraction to expansion.

Conclusions
In this work, we have presented a comprehensive study of collective cell durotaxis based on a continuum model of cell monolayers as a two-dimensional active fluid on a gradient of substrate stiffness. The stiffness gradient affects tissue dynamics through a spatially dependent traction and friction coefficients ζ i (x) and ξ(x). We analytically solve the model in a one-dimensional setup. The tissue dynamics is characterized by two main observables: the velocity of the center of mass and the spreading velocity. For the simple case of a uniform traction gradient ζ i and uniform friction (ξ = 0), the spreading velocity is exactly the same as that for the uniform substrate case, so the spreading behavior is independent of the existence of a traction gradient. The velocity of the center of mass, instead, is finite and pro- Fig. 13 Examples of monolayer spreading dynamics to illustrate the effect of elasticity k (first row) and surface tension γ (second row). In each plot, curves of the same color show the evolution of the position of the edges x±(t), filling the area between them to represent the tissue width. In the first row, γ = 0 and L r = 150 µm. In (a), the initial size is L0 = 215 µm and k = 0, 0.03, 0.05, 0.5 kPa. In (b), L0 = 100 µm and k = 0, 2, 3, 5 kPa. The elastic constant k increases from lighter to darker curves. In (c), k = 0 and only the L0 = 215 µm case is shown with γ = 0, 1, 3, 10 mN/m, which also increases from lighter to darker green.
Other parameter values are those in Table 1 portional to ζ i . Therefore, the cell monolayer performs durotaxis as long as the traction is a monotonically increasing function of the substrate stiffness. These conclusions are locally valid for more general traction profiles provided that the gradient does not change significantly over the monolayer width.
We have analyzed the durotactic dynamics as a function of physical parameters, for example discussing the wet and dry limits that result from comparing the monolayer size L to the hydrodynamic screening length λ. For broad ranges of values and profiles physical parameters, and for different boundary conditions, we have characterized the different regimes that result from combining states of spreading and contraction with states of interface wetting and dewetting. All of them give rise to durotactic motion. The durotactic velocity increases with both the traction gradient ζ i and the monolayer size L, as seen in Ref. [13]. Moreover, for uniform ζ i and uniform friction, the durotactic velocity is independent of the contractility ζ and the traction offset ζ 0 i . Therefore, the same monolayer placed at different positions along the stiffness gradient would have the same durotactic velocity but different spreading dynamics.
However, for non-uniform friction (ξ = 0), as well as for traction and friction profiles that saturate with stiffness, the model predicts lower velocities for larger stiffness offsets, thus recovering the results in Ref. [9,13,14]. At high stiffness, parameter saturation makes the system asymptotically approach the dynamics on uniform substrates, with vanishing durotaxis. The spreading velocity increases with the traction offset ζ 0 i , the monolayer size L, and the hydrodynamic length λ, but it decreases with contractility.
In addition to the predictions for local durotaxis and spreading, we have discussed the temporal evolution of a monolayer along the stiffness gradient as it changes its position and size. We have identified three regimes for the evolution of monolayer size. Large monolayers spread indefinitely, small monolayers contract indefinitely, and monolayers in an intermediate size range display a non-monotonic evolution whereby they switch from contraction to spreading at a finite time. These three regimes are separated by two critical lengths, which we determined analytically in some simple cases, and we illustrated numerically for more general situations. We also discussed the effect of additional physical ingredients such as surface tension and elastic forces that oppose large deformations of the tissue. We have shown that they typically slow down the expansion and accelerate the contraction.
Our model is relatively simple and strongly predictive, so it could be tested in experiments and used to infer parameter values from experimental data. It could also guide the design of further experiments on collective durotaxis. Nevertheless, the model has obvious limitations: It is restricted to cell monolayers, and it is unrealistic in its long-time behavior. Addressing these limitations would require to include additional physics such as effects from the three-dimensional, multiplelayer structure resulting from monolayer contraction, additional forces to prevent indefinite spreading, and other ingredients such as cell proliferation. All these generalizations of the model are deferred to future work. Finally, it is a question of great interest to elucidate to what extent a purely mechanical description, with no need to invoke biochemical signaling, can account for the observed phenomenology in different forms of collective cell migration. any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendices Appendix A: Uniform substrate
The solution to Eq. (6) for a constant traction ζi and friction ξ is obtained assuming a normal component of the stress in the boundaries due to two different effects: an effective surface tension γ (interpreting our 1d model as an approximation for a circular cluster of radius L), and an effective elastic stiffness k, accounting for a mean-field-type linear elastic interaction as in Ref. [49] that prevents the tissue from excessive stretching, being L r the reference length (justification more extended in Sect. 4.3). Thus, and further, in the wet (L λ) and dry (L λ) cases, v wet (A.6) Setting γ = 0 and k = 0 and neglecting Lc in front of λ or L, we obtain Eqs. (8) and (9), respectively.