Competing elastic and viscous gradients determine directional cell migration

Mathematical Biosciences 380 (2025) 109362 Available online 17 December 2024 0025-5564/© 2025 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs Original research article Competing elastic and viscous gradients determine directional cell migration Pablo Saez a,b,∗, Pallavi U. Shirke c, Jyoti R. Seth c, Jorge Alegre-Cebollada d, Abhijit Majumder c aLaCàN, Universitat Politècnica de Catalunya-BarcelonaTech, 08034 Barcelona, Spain bInstitute of Mathematics of UPC-BarcelonaTech.–IMTech, Barcelona, Spain cDepartment of Chemical Engineering, Indian Institute of Technology Bombay (IITB), 400076 Mumbai, India dCentro Nacional de Investigaciones Cardiovasculares (CNIC), 28029 Madrid, Spain ARTICLE INFO Keywords: Clutch model Active gel models Viscotaxis Durotaxis Cell adhesion Mechanotransduction ABSTRACT Cell migration regulates central life processes including embryonic development, tissue regeneration, and tumor invasion. To establish the direction of migration, cells follow exogenous cues. Durotaxis, the directed cell migration towards elastic stiffness gradients, is the classical example of mechanical taxis. However, whether gradients of the relaxation properties in the extracellular matrix may also induce tactic responses (viscotaxis) is not well understood. Moreover, whether and how durotaxis and viscotaxis interact with each other has never been investigated. Here, we integrate clutch models for cell adhesions with an active gel theory of cell migration to reveal the mechanisms that govern viscotaxis. We show that viscotaxis is enabled by an asymmetric expression of cell adhesions that further polarize the intracellular motility forces to establish the cell front, similar to durotaxis. More importantly, when both relaxation and elastic gradients coexist, durotaxis appears more efficient in controlling directed cell migration, which we confirm with experimental results. However, the presence of opposing relaxation gradients to an elastic one can arrest or shift the migration direction. Our model rationalizes for the first time the mechanisms that govern viscotaxis and its competition with durotaxis through a mathematical model. 1. Introduction Cell migration determines key biological processes such as embryonic development, tissue regeneration, wound healing, and tumor invasion [1–3]. For decades, there has been increasing interest in understanding how cells organize themselves to migrate. There is also growing concern about controlling cell migration, as it may allow us to arrest tumor metastasis [4–6], enhance tissue regeneration, and design better biomimetic materials [7,8]. Three forces cooperate to enable cell migration on two-dimensional (2D) substrates (see Fig. 1). First, the F-actin network, along with myosin motors that pull on actin filaments, creates an active viscoelastic flow, or retrograde flow [9, 10]. Second, the F-actin network continuously polymerizes, pushing the cell membrane outward [11,12], creating a differentiated actin network compared to the contractile actomyosin network. These two competing actin networks, the first moving inwards and the latter outwards, determine the protruding velocity of the cell contour and, therefore, cell migration. Third, cells create adhesions that link the intracellular actomyosin network with the extracellular matrix (ECM). Cell adhesion is mediated by molecules that sense and respond to the mechanical properties of the ECM [13–15]. Through the formation and disassembling of adhesion complexes, cells tune the friction between the intracellular network and the ECM. ∗Corresponding author at: LaCàN, Universitat Politècnica de Catalunya-BarcelonaTech, 08034 Barcelona, Spain. E-mail address: [email protected] (P. Saez). But then, how does a cell organize itself to determine the direction of migration? Most cells require stimuli to break the symmetry in the intracellular force generation or extracellular force transmission to enable directed cell migration. Otherwise, cells would just spread with symmetric distributions of F-actin retrograde flow, polymerization, and adhesion formation [16,17]. Mechanical [3,18–20] and chemical signals [21,22] have been well studied. In terms of mechanical cues for cell migration, which we analyze here, most studies have focused on durotaxis, the motion of cells towards elastic stiffness gradients of the ECM [18,19,23,24]. Durotaxis is expressed across different cell types, including fibroblasts [25], smooth muscle cells [26], neurons [24], and cancer cells [27]. When cells are placed on ECM with an elastic stiffness gradient, they establish an asymmetric adhesive state [28,29], with one side of the cell more adhesive than the other. In response to this friction gradient between the cell and the substrate, the retrograde flow also polarizes [29] and establishes a prospective cell front and rear. Moreover, the direction of the adhesion gradient determines whether the cell front orients towards the positive or negative elastic stiffness gradient [29]. However, extracellular matrices in vivo are not purely elastic but also present viscous properties [30]. The viscoelasticity of the ECM https://doi.org/10.1016/j.mbs.2024.109362 Received 16 June 2024; Received in revised form 3 December 2024; Accepted 4 December 2024 Mathematical Biosciences 380 (2025) 109362 2 P. Saez et al. Fig. 1. (a) Sketch of the main forces acting in cell migration velocity (v). A contractile actomyosin flow, (𝑣𝐹black arrows), competes with protrusive forces at the front (𝑙𝑓(𝑡)) and rear (𝑙𝑟(𝑡)) of the cell (𝑣𝑝, blue arrows). At the cell-ECM interface, intracellular forces are transmitted extracellularly through a friction force that appears due to adhesion mechanisms. (b) A standard linear model describes the viscoelastic behavior of the ECM, determined by 𝑘𝑒,𝑘𝑣, and 𝜇. On the right are the force and displacement curves of a relaxation test. (c) Detail of an adhesion complex, made of adaptor proteins (top) and results of a clutch model (bottom). Results show the probability of bound binders, 𝑃𝑏, showing cycles of formation and rupture of the adhesion complex, the velocity of the actin filaments, 𝑣, and traction forces on the substrate, 𝑃. mediates several cell functions, including the formation of adhesion complexes [31–34], cell differentiation [34,35], proliferation [34,36], lamellipodial and filopodial protrusion [37], and cell spreading [31– 33,38]. Most studies have analyzed the effect of homogeneous viscous properties of the ECM on cell function. How cells migrate along gradients in the relaxation properties of the ECM has only been preliminarily analyzed experimentally, without exploring the underlying physical processes that determine viscotaxis [39] or through physical models that neglect the heterogeneous internal retrograde flow and actomyosin distribution [40], which is cause and effect of directed cell migration. Here, based on the same principles that drive durotaxis, we reveal through physical models mechanisms of cell viscotactic migration in 2D. Moreover, because elastic and relaxation gradients may appear simultaneously in vivo and in vitro, we analyze how durotaxis and viscotaxis compete with each other, which has not been addressed by any previous work. A clear understanding of how elastic and viscous gradients affect each other can allow us to better engineer biomaterials for specific biological functions. In this work, we first describe the mathematical models used for cell adhesion (the clutch model) and for the cell motility (an active gel model). Next, we discuss how we prepared the gel and tested mechanically. We finish the methods section by describing cell culture, time lapse microscopy and cell migration analysis of our work. The, in our results, we first show the cell adhesion behavior as a function of the elasticity and viscoelasticity of the ECM, which we later use to propose a mathematical model of single cell viscotaxis. Then, we explain theoretically how viscotaxis appear and show that viscotaxis can compete with durotaxis to control cell migration. Finally, we conclude discussing the main results of our work. 2. Materials and methods 2.1. Clutch model for cell adhesions based on Montecarlo simulations We use previous clutch models to understand the response of cell adhesion to the viscoelasticity of the ECM [29,41–43] (Fig. 1(b-c)). It considers several molecular clutches that attach to the ECM on one side and to the adaptor talin rod on the other. At its other end, talin is linked to actin filaments pulled by myosin motors. The computation of the clutch model relies on repeated random sampling, which we solve by a Monte Carlo (MC) simulation, during which many events of engagement and disengagement occur, resulting in cycles of formation and rupture of adhesion complexes (see Fig. 1c). We use a constant time step, 𝛥𝑡 =0.005 s, and a final time of 𝑡𝑓= 100 s, which was shown to provide accurate results [43]. Then, the MC simulations are averaged over time to obtain the results at specific mechanical properties of the ECM. During the cycling process of cell adhesion, the actomyosin network pulls on the adhesion molecules thanks to the contractile forces of the myosin motors. The actin velocity, 𝑣, is given by the force-velocity relation 𝑣 =𝑣𝑢(1 −𝐹∕𝐹𝑠𝑡𝑎𝑙 𝑙),(1) where 𝑣𝑢is the free velocity of the actin flow. 𝐹𝑠𝑡𝑎𝑙 𝑙is the stall force of the myosin motors, being 𝐹𝑠𝑡𝑎𝑙 𝑙=𝑛𝑚𝐹𝑚, where 𝐹𝑚is the force required to stall the activity of one myosin motor and 𝑛𝑚the number of motors in the system. 𝐹is the reaction force in the substrate, which balances the force on all the molecular clutches, 𝐹=∑𝑛𝑒𝑛𝑔 𝑖=1 𝐹𝑐 ,𝑖, being 𝑛𝑒𝑛𝑔 the number of engaged clutches. The displacement of the engaged clutches is computed as 𝛥𝑧 =𝑣𝛥𝑡, and the displacement of the substrate at the current time step, 𝑧𝑛+1, is obtained by solving the balance of forces between the 𝑛𝑒𝑛𝑔 engaged molecular clutches and the substrate, from where we can later obtain 𝐹. The force at each 𝑖th clutch, 𝐹𝑐 ,𝑖, is given by 𝐹𝑐 ,𝑖 =𝜅𝑐(𝑧𝑐 ,𝑖 −𝑧𝑠𝑢𝑏), where 𝜅𝑐is the stiffness of each molecular clutch (see [41–43] for details) and 𝑧𝑐 ,𝑖 is the displacement of the 𝑖th molecular clutch, 𝑧𝑠𝑢𝑏 is the displacement of the substrate. Once 𝐹is computed, the cell traction can be obtained, 𝑃=𝐹∕𝜋 𝑎2, where (𝜋 𝑎2)is the area occupied by a circular adhesion complex. Moreover, once the force at each clutch, 𝐹𝑐 ,𝑖, is computed, we can calculate the binding and unbinding events at each molecular clutch. The attachment of the 𝑛𝑐molecular clutches to the ECM, which represents the integrin-ECM link, follows a catch-bond behavior [42,44]. Several molecular clutches bind to the ECM with a constant rate 𝑘𝑜𝑛, and unbind with a dissociation rate 𝑘∗ 𝑜𝑓 𝑓, that depends on the force on each molecular binder, 𝐹𝑐, as 𝑘∗ 𝑜𝑓 𝑓=𝑘𝑜𝑓 𝑓 𝑠𝑒(𝐹𝑐∕𝐹𝑠𝑙 𝑖𝑝 𝑏)+𝑘𝑜𝑓 𝑓 𝑐𝑒(−𝐹𝑐∕𝐹𝑐 𝑎𝑡𝑐 ℎ 𝑏),(2) where the first and second terms represent a slip and catch bond, respectively. 𝑘𝑜𝑓 𝑓 𝑠and 𝑘𝑜𝑓 𝑓 𝑐and 𝐹𝑠𝑙 𝑖𝑝 𝑏and 𝐹𝑐 𝑎𝑡𝑐 ℎ 𝑏are the unloaded dissociation rates and the characteristic bond rupture forces in the slip and catch pathways, respectively. Mechanosensitivity is promoted by talin unfolding, which induces vinculin recruitment and subsequent integrin clustering [42]. This process of integrin crowding is usually termed adhesion reinforcement. We introduce talin reinforcement by considering talin unfolding according to Bell’s model at a rate 𝑘∗ 𝑢𝑛𝑓 [42]. Then, it can either refold again or vinculin can bind to it at a force-independent rate 𝑘𝑜𝑛𝑣. If talin unfolds, Mathematical Biosciences 380 (2025) 109362 3 P. Saez et al. Table 1 Left: Parameters adopted for the clutch model of cell adhesion. The model parameters for the catch bonds with talin reinforcement are obtained from [42,43]. The unbinding rate is considered to follow a cycle of mechanical reinforcement, following 𝑘∗ 𝑜𝑓 𝑓= 8.10−4 𝑒(𝐹𝑐∕8.16) + 10.14 𝑒(−𝐹𝑐∕6.24) + 900 𝑒(−𝐹𝑐∕0.01) [42]. Right: Parameters adopted for the cell migration model and related references. * indicates that the value was adjusted in this work. Clutchs parameter Migration parameters 𝑎(nm) 1700 𝜇𝐹(kPa s) 10 [48,49] 𝑛𝑐1200 𝜂0(kPa s/μm2) 0.05 [48,49] 𝜅𝑐(pN∕nm) 1000 𝜁(kPa) 0.05 [47,48] 𝑘𝑜𝑛𝑡(μm2∕s) 2.11 × 10−4 𝑘0.1 [49] 𝑑𝑖𝑛𝑡 (int /μm2) 300 𝐷𝐹(μ m∕s) 0.3 [48,50] 𝐹𝑚(pN) 2 𝐷(μ m∕s) 0.4 [48,50] 𝑛𝑚800 𝑘𝑑0.1 [51,52] 𝑣𝑢(nm∕s) 110 𝑘𝑝0.1 [51,52] 𝑘𝑜𝑛𝑣 (s−1)1 × 108𝑣𝑝 0(μ m∕s) 0.55 [53] 𝑖𝑛𝑡𝑎𝑑 𝑑(int /μm2) 24 𝐿𝑏(μm) 7 * 𝑚𝑟(int /μm2) 15 000 𝜏𝑠𝑡𝑎𝑙 𝑙(nN∕μm) 0.4 * 𝛥𝑡 (s) 0.005 𝑡𝑓(s) 100 it can either refold again or vinculin can bind to it. Then, integrin density increases by 𝑖𝑛𝑡𝑎𝑑 𝑑= 24 integrins/𝜇 𝑚2and the binding rate of integrins results in 𝑘𝑜𝑛 =𝑘𝑜𝑛𝑡𝑑𝑖𝑛𝑡, where 𝑘𝑜𝑛𝑡 is the true binding rate characterizing each integrin-fibronectin bond, and 𝑑𝑖𝑛𝑡 is the density of integrins in the adhesion complex. We do not allow integrin density to go above a maximum integrin density, 𝑚𝑟, determined by the separation between integrins (see [43] for details). If the clutch unbinds before vinculin binds to talin, integrin density is decreased by 𝑖𝑛𝑡𝑎𝑑 𝑑, reflecting the fact that adhesions lose integrins if force application is decreased [45,46]. To model the viscoelastic behavior of the ECM, we use the Maxwell representation of a Standard Linear Solid (SLS) model such that the stress–strain relation is given as 𝐹+𝜇∕𝑘𝑣 𝐹=𝑘𝑒𝑧𝑠𝑢𝑏 +𝜇(𝑘𝑣+𝑘𝑒)𝑧𝑠𝑢𝑏∕𝑘𝑣.(3) 𝑘𝑒is the long-term stiffness, 𝑘𝑣is the relaxation stiffness, which represents the drop to the long-term stiffness as the material relaxes, 𝑧𝑠𝑢𝑏 =𝑣, and 𝜇is the viscosity of the ECM (see Fig. 1(b)). The time derivatives are discretized using a backward Euler method, which results in an equation for 𝑧𝑛+1 𝑠𝑢𝑏 as a function of 𝐹𝑛+1. Here, we assume that the friction between the cell and the ECM is 𝜂=𝑃∕𝑣 [29,47], so friction changes as a function of the ECM stiffness similar to how actin velocity and tractions change. This relation assumes that strong and stable adhesions, which are represented by large traction forces, induce large friction in comparison to small and weak adhesion complexes. All model parameters related to the clutch model are summarized in Table 1. 2.2. Active gel model and numerical solution of the system of equations Recently, we and others have shown that durotaxis emerges from an asymmetry in the cell adhesion strength [28,29]. Hence, our results on these friction gradients due to gradients in the relaxation properties of the ECM strongly suggest that single cells may also express viscotaxis. To analyze the exposure of single motile cells to relaxation gradients of the ECM, we use an active gel model of cell migration [29, 47,49,54]. We consider a 1D domain, 𝛺, with moving coordinates 𝑥(𝑡) ∈ [𝑙𝑟(𝑡), 𝑙𝑓(𝑡)], where 𝑙𝑟(𝑡)and 𝑙𝑓(𝑡)represents the rear and front boundaries of the cell and, therefore, the cell length is determined as 𝐿(𝑡) =𝑙𝑓(𝑡) −𝑙𝑟(𝑡)(Fig. 1). The retrograde flow is modeled with a contractile viscous gel in contact with the ECM. The constitutive relation for the stress in the active fluid is given as 𝜎𝐹=𝜇𝐹𝜕𝑥𝑣𝐹+𝜁 𝜌𝐹𝜌𝑀.(4) 𝑣𝐹is the velocity of the retrograde flow, 𝜇𝐹is the shear viscosity, 𝜁the active contraction exerted by the myosin motors and 𝜌𝑀and 𝜌𝐹are the densities of myosin motors and F-actin, respectively (see below). The balance of linear momentum for the active gel is then 𝜕𝑥𝜎𝐹=𝜂𝐹𝑣𝐹.(5) We assume zero stresses on the cell boundaries. The right-hand side represents the friction between the intracellular network and the ECM, where the clutch model for cell adhesion described above couples with the active gel model for cell migration. The friction parameter is 𝜂𝐹=𝜂0+𝜂, where 𝜂is described above and 𝜂0=0.05 kPa⋅s∕μm2is added to the effective friction 𝜂as baseline friction of the retrograde flow with the surrounding cytoskeletal structures. To compute the velocity of actin network polymerization against the membrane, which we consider as a separated F-actin structure, we use previous models (see [53,55,56] for details) where actin polymerization decreases when the membrane tension increases as 𝑣𝑝= 𝑣𝑝 0[1 −𝜏(𝐿(𝑡))∕𝜏𝑠𝑡𝑎𝑙 𝑙]8from a free-growth velocity 𝑣𝑝 0.𝜏𝑠𝑡𝑎𝑙 𝑙is a tension required to stall the polymerization of the actin network. The membrane tension is calculated as a linear Hookean spring, 𝜏(𝐿(𝑡)) =𝑘(𝐿(𝑡) − 𝐿𝑏), 𝐿(𝑡)< 𝐿𝑏, where 𝑘is the spring constant and 𝐿𝑏accounts for the resting length and the buffer membrane length of the cell membrane. The model assumes zero compressive stresses. The outward polymerization velocities, 𝑣𝑝 𝑟,𝑓 , compete with the inward retrograde flow velocity, 𝑣𝐹|𝑟,𝑓 , at the cell rear and front. This competition makes the cell ends expand or retract with velocity  𝑙𝑟,𝑓 (𝑡) =𝑣𝑝 𝑟,𝑓 −𝑣𝐹|𝑟,𝑓 , respectively (Fig. 1) and the migration velocity is defined as 𝑣= ( 𝑙𝑟(𝑡) − 𝑙𝑓(𝑡))∕2, the mean value of the front and rear velocity of the cell. The distribution of the F-actin density, 𝜌𝐹(𝑥, 𝑡)is modeled through the following transport equation, 𝜕𝑡𝜌𝐹+𝜕𝑥(𝑤𝜌𝐹−𝐷𝐹𝜕𝑥𝜌𝐹)=𝑘𝑝−𝑘𝑑𝜌𝐹,(6) where the right-hand side includes the polymerization and depolymerization terms with rates 𝑘𝑝and 𝑘𝑑, respectively [52], 𝐷𝐹is the diffusive parameter of the F-actin. Similarly, the myosin motors bound to the F-actin network [57], 𝜌𝑀, as 𝜕𝑡𝜌𝑀+𝜕𝑥(𝑤𝜌𝑀−𝐷 𝜕𝑥𝜌𝑀) = 0,(7) where 𝐷is the effective diffusion parameter [49,50]. F-actin and bound myosin motors are also drifted in the cell frame with velocity 𝑤=𝑣𝐹−𝑣. We impose zero fluxes on the cell boundaries to reflect that no F-actin or myosin motors can enter or leave the cell domain. To solve the system of PDEs described above, we use a staggered finite element method to discretize the system in space, and an implicit second-order Crank–Nicholson method to discretize the equations in time [58,59]. To avoid undesired oscillations from the numerical solution of the parabolic equations if the problem becomes convective dominant, i.e. 𝑃 𝑒 >1, we include the Stream-Upwind Pretrov Galerkin (SUPG) stabilization term. This strategy allows us to use finite elements of constant size ℎ=𝑙𝑓(𝑡)∕𝑁, where 𝑁is the number of elements and, consequently, keep the number of elements of our domain constant. The complete finite element and time discretization have been presented previously [57]. All model parameters related to the migration model are summarized in Table 1. 2.3. Gel preparation and mechanical testing All the experimental procedures and experimental data described in this work for the gel mechanical properties are taken from a previous publication [39]. Viscoelastic properties, represented by the storage modulus (G’) and loss modulus (G’’), were measured for each gel. Creep measurements were also performed with both gels where constant shear stress, 𝜎0, was applied and strain 𝛾was measured for an hour. Mathematical Biosciences 380 (2025) 109362 4 P. Saez et al. Fig. 2. (a) Friction of the cell-ECM adhesion, 𝜂, for 𝑘𝑒=1 pN/nm and a range of values of the relaxation stiffness, 𝑘𝑣, and viscosity, 𝜇. (b) Friction as a function of the viscosity, 𝜇, for 𝑘𝑣=0.2 pN/nm (dash line) and 𝑘𝑣=20 pN/nm (solid line). (c) Friction as a function of the relaxation stiffness, 𝑘𝑣, for 𝜇=2 pN⋅s/nm (dash line) and 𝜇=20 pN⋅s/nm (solid line). Here, we fit the experimental compliance (J= 𝛾(t)/𝜎0) with the compliance expression of a SLS model, 𝐽(𝑡) =1 𝑘𝑒 (1 − [𝑘𝑣 𝑘𝑒+𝑘𝑣 ]𝑒𝑥𝑝(− 𝑘𝑒𝑘𝑣 𝜂(𝑘𝑒+𝑘𝑣)𝑡)),(8) using the non-linear least squares fitting algorithm in MATLAB. Since the initial compliance 𝐽(0) = 1∕(𝑘𝑒+𝑘𝑣)is known and is equal for both high-loss and low-loss gels (as seen from the experimental data), 𝑘𝑒+𝑘𝑣 was constrained to remain constant for the two fits. 2.4. Cell culture, time lapse microscopy and cell migration analysis Again, the experimental data described in this work for the cell culture are taken from a previous publication [39]. In short, umbilical cord human Mesenchymal stem cells (UC-hMSCs) between passages 4 and 7 were used for the experiments. Then, the paths of the cells were tracked using the ImageJ software (National Institutes of Health, Bethesda, USA) plugin with manual tracking developed by Fabrice Cordelières at Institut Curie, Orsay, France, and further migration analysis was done using the ‘‘Chemotaxis and Migration Tool’’ plugin developed by IBIDI Software. were used for the experiments. We refer to our previous work for details [39]. 3. Results 3.1. Adhesion behavior as a function of the elasticity and viscoelasticity of the ECM We first analyze the cell adhesion behavior for substrates of different long-term stiffnesses, relaxation stiffnesses, and viscosities (Fig. 2). Changes in long-term stiffness are related to durotaxis. In agreement with previous models [29,42], the friction increases as we increase 𝑘𝑒 (see Fig. 2(a)). To analyze the viscoelastic effects, we fix the long-term stiffness, 𝑘𝑒=1pN/nm, and analyze the effect of the relaxation stiffness 𝑘𝑣and viscosity 𝜇. We show an increase in the friction between the cell and the substrate as we increase 𝑘𝑣and 𝜇(Fig. 2(b)). For a relaxation stiffness lower than ≈1 pN/nm, friction vanishes for all values of ECM viscosities, 𝜇 <0.5pN⋅s/nm (Fig. 2(c)). For values of 𝑘𝑣bigger than ≈1 pN/nm, friction increases as the viscosity of the ECM increases (Fig. 2(c)). For constant values of viscosity, we show a transition for friction values. For viscosities larger than ≈10 pN⋅s/nm, we observe monotonically increasing values of friction, that reach ≈2 kPa⋅s∕μm2. For viscosities smaller than ≈4 pN⋅s/nm, we see a biphasic behavior. Friction increases up to ≈0.7 kPa⋅s∕μm2at 𝑘𝑣≈10 pN/nm, and then it decreases, describing a negative gradient for the relaxation stiffness (Fig. 2(c)). Our results confirm previous data on how cell adhesion responds to the elastic stiffness of the ECM [29,41,42], showing that cells migrate towards positive gradients of friction created by the asymmetric strength in cell adhesions [28,29]. When cells express talin reinforcement, as in our model, this migration occurs towards the positive gradient of the elastic stiffness [29,42]. Further, we show that an increase in the viscosity for constant values of relaxation stiffnesses behaves equivalently to a stiffness increase in pure elastic substrates with cells expressing talin reinforcement [42]. However, an increase in the relaxation stiffness induces, for increasing values of viscosity, a transition from an adhesion behavior reminiscent of non-reinforced adhesion complexes to reinforced ones (Fig. 2(c)). 3.2. Theory predicts viscotaxis in single cells We simulate cells placed on substrates of length 𝐿𝑠equal to 200 and 1000 μm (Fig. 3). To obtain values of friction in the moving coordinate system that represents the migrating cell, 𝑘𝑒,𝑘𝑣, and 𝜇become functions of space that follow a linear relation with positive (e.g. 𝑥𝑘𝑣∕𝐿𝑠) or negative slope (e.g. 𝑘𝑣(𝐿𝑠−𝑥)∕𝐿𝑠). Because the absolute viscoelastic properties in which cells are initially located modify the friction they experience, we also analyze three different initial locations. In these three samples, we fix the long-term stiffness to 𝑘𝑒= 1 pN/nm and we vary independently 𝑘𝑣and 𝜇from 10−2to 102(Fig. 2b-c). We consider that cells undergo viscotaxis if the migration velocity is larger than 𝑣= 1nm/s, otherwise we assume that they spread and remain stationary. Cells on the 200 μm sample express the strongest viscotactic response compared to the 1000 μm sample, which did not express viscotaxis (Fig. 3). In all cases, the strongest viscotactic response is obtained for cells located at the largest gradient of the sample. Therefore, to analyze the differences between the viscous properties of these matrices in detail, we focus on the 200 μm sample Fig. 3(a). Cells initially located at the center of the sample show a persistent positive viscotaxis with velocities of ∼5 nm/s for samples of varying 𝑘𝑣and fixed 𝜇=2 and 20 pN⋅s/nm. Same behavior is observed for samples of varying 𝜇and 𝑘𝑣≈20 pN/nm. On the other hand, cells initially located at the center of the sample but of a constant 𝜇= 0.02 pN⋅s/nm remained stationary because there was no friction gradient for the cells to polarize. The same behavior is observed for samples of varying 𝜇and 𝑘𝑣≈0.02 pN/nm (see Fig. 2b-c). However, when cells are initially located at that same middle spot but for samples of varying 𝑘𝑣and constant 𝜇= 2 pN⋅s/nm, they express negative viscotaxis as they move towards negative gradients of the viscosity, which correspond to regions where the friction gradient is negative (Fig. 2c). Our results show that in most cases cells migrate towards the positive gradients of the relaxation parameters 𝑘𝑣and 𝜇of the ECM (Fig. 3), i.e. expressed positive viscotaxis. However, there are specific combinations of viscous properties that induce negative viscotaxis. Our results suggest, in agreement with previous results on durotaxis [28, 29], that gradients in the relaxation stiffness of the ECM can express both positive and negative viscotaxis. However, gradients in 𝜇can only induce positive viscotaxis. All these results are explained by the adhesion strength shown in the previous section (Fig. 2). In other words, our results indicate that cells will always migrate with the cell front oriented with the strongest adhesion or, equivalently, towards positive gradients of friction [29]. The asymmetric adhesion forces induce asymmetric retrograde flow which, given the constant polymerization velocity at both cell edges, establishes the cell front. 3.3. Viscotaxis is mediated by intracellular asymmetric motile forces To get a deeper understanding of how viscotaxis emerges and progresses, we analyze the intracellular motile forces that coordinate Mathematical Biosciences 380 (2025) 109362 5 P. Saez et al. Fig. 3. Migration velocity over time for substrates of 200 μm (a) and of 1000 μm (b) in length. Left and right columns are for cells located at points ii and iii, respectively. Cells located at point i did not migrate (data not shown). For each substrate, results in the first row are for fixed values of 𝜇to 0.02 (blue), 2 (red), and 20 (black) pN⋅𝑠/nm while cells migrate on a substrate of varying 𝑘𝑣. In the second raw, we fixed values of 𝑘𝑣to 0.02 (blue), 0.2 (red), and 20 (black) pN/nm while cells migrate on a substrate of varying 𝜇. Fig. 4. For the black line case (i in Fig. 3a, bottom-left) we show the evolution of intracellular variables in space and time. At steady state, (a) actin (black) and myosin (blue) densities and (b) retrograde flow, (c) polymerization velocity (solid) and the retrograde velocity at the cell membrane (dash) in time, and (d) kymograph of the friction 𝜂. cell migration in space and time. We focus on the case that showed the strongest viscotactic response (sample with a viscous gradient of 200 μm in length, fixed 𝑘𝑒=0.1 pN/nm and 𝑘𝑣= 20 pN/nm, and cells initially located at 𝜇= 20 pN⋅s/nm). As cells migrate, they travel through regions of different viscosities and, therefore, the friction computed through the clutch model changes in space and time (Fig. 4d). This heterogeneous friction along the cell creates an asymmetric retrograde flow (Fig. 4b, which competes with the constant polymerization velocity 𝑣𝑝at the cell ends (Fig. 4c). The right end protrudes faster and, therefore, represents the cell front. The asymmetric retrograde flow also polarizes the actomyosin network (Fig. 4a). All other cases of viscotaxis described above exhibit similar results, with the model variables polarized towards the positive or negative direction, and at a smaller or larger extent, depending on the strength of the cell friction. Because viscotaxis and durotaxis emerge from the same force principles, the strength of the cell adhesion and, therefore, the strength of the asymmetric friction values, our results are in agreement with the behavior observed during durotaxis [28,29]. 3.4. Competition between viscotaxis and durotaxis Then, we wonder how viscotaxis competes with durotaxis. Gradients of the elastic stiffness may coexist with gradients of the relaxation properties of the ECM in vivo. Therefore, it is important to understand how these two different but coexisting mechanical signals may control cell migration. To do so, we take our previous model of durotaxis [29] and put it in competition with our current model of viscotaxis. We name this cooperative response as viscoelastotaxis. We use again the strongest viscotactic response described above. However, we let the long-term stiffness 𝑘𝑒change to enable durotaxis. We also use the strongest durotactic response for a sample of 200 μm in length where the long-term stiffness changes between 10−2to 102 pN/nm. We simulate matrices where the long-term and viscosity stiffness change along the same direction or opposite to each other. For aligned viscoelastic and elastic gradients, the friction profile is similar to pure viscotaxis (Fig. 5a), although the migratory ability is enhanced and cells reach a maximum velocity of 6 nm/s. However, when the two mechanical gradients of the ECM are in opposite directions, cells may express positive, negative, or null viscoelastotaxis. We show a minimum in friction at ≈1 pN/nm relaxation stiffness (Fig. 5(a). We term this spot diverging stiffness. Cells placed initially on the left of the diverging stiffness express negative viscoelastotaxis while cells on the right express positive durotaxis. Cells located at the diverging stiffness would tend to migrate away from it contrary to durotaxis, where an optimum rigidity marks a goal location for the cells to migrate [29]. Finally, we use previous experimental results [39] to validate our computational model. Cells were located in substrates of different viscous and elastic properties. In summary, two drops of pre-polymer solutions, which we also tested mechanically, are placed in contact and diffuse with each other to form a relaxation gradient (see [39] and Materials & Methods for details). The frequency sweep test showed that these materials have the same storage modulus (G’) but different loss moduli (G’’). From the small amplitude oscillation regime, we show that the elastic moduli were nearly the same for both gels, which also translates to equal short-term for both gels (see Fig. 6a). These two polymers were also tested in a creep test for applied stress of 𝜎0=10 Pa. High loss and low loss gels deform equally at short time scales (t<2 s). Fitting of the material parameters to the SLS model Mathematical Biosciences 380 (2025) 109362 6 P. Saez et al. Fig. 5. (a) Effective friction for the aligned elastic and relaxation gradients (solid line) and opposing elastic and relaxation gradients (dash line). Gradient direction is specified by a triangle. Friction increases monotonically for increasing values of the elastic stiffness, while it presents a concave shape with a minimum at ≈ 1pN/nm, where the diverging rigidity localizes. (b) Migration velocities of single cells in a sample of 200 μm in length exposed to aligned and opposing durotactic and viscotactic signals. Black line corresponds to the case of aligned gradients. For opposing cues, cells are initially located at the red cross (𝑘𝑒=0.1 pN/nm and 𝑘𝑣=10 pN/nm), blue cross (𝑘𝑒=1 pN/nm and 𝑘𝑣1=pN/nm) and yellow cross (𝑘𝑒=10 pN/nm and 𝑘𝑣=0.1 pN/nm) in panel a. Colors in velocities in panel b correspond to these colors in panel a. Fig. 6. (a) Experimental frequency sweep G’ modulus (diamonds) and G’’ modulus (circles) for Low Loss (blue) and High Loss (red) gels. (b) Compliance (J) of the Low Loss (blue) and High Loss (red) gels. The experimental compliance (circles) was fitted following the SLS model (lines). (c) A collection of tracks for many cells obtained from multiple independent experiments of hMSCs migrating on the gels with gradients in the viscoelastic properties. The material varies from 𝑘𝑒= 1.06 pN/nm, 𝑘𝑣= 0.28 pN/nm and 𝜇=38.5 pN⋅s/nm on the left to 𝑘𝑒= 0.22 pN/nm, 𝑘𝑣= 1.11 pN/nm and 𝜇=46.2 pN⋅s/nm on the right (scale bar 50 μm). (d) Normalized velocities with respect to the control case in (c) . Colors represent the control case in black, fixed 𝜇in blue, fixed 𝑘𝑒in red, fixed 𝑘𝑒and 𝑘𝑣in green, and fixed 𝑘𝑣in magenta. (Fig. 6b). showed that the first material has 𝑘𝑒= 1.06 pN/nm, 𝑘𝑣= 0.28 pN/nm and 𝜇=38.5 pN⋅s/nm and the second has 𝑘𝑒= 0.22 pN/nm, 𝑘𝑣= 1.11 pN/nm and 𝜇=46.2 pN⋅s/nm. This represents a complex viscotactic state where gradients of long-term and viscous modulus and viscosity are established in opposite directions and the prediction of what direction cells would move towards is not clear. Then, cells were seeded on the substrate and the migration direction was measured. 23 cells out of 33 analyzed migrated towards the positive gradient of the long-term modulus, in opposite directions to the viscous gradients (Fig. 6c). Then, we simulate the same viscoelastotaxis case. Our results reproduce the migration towards the positive gradient of the long-term modulus (see control case, Fig. 6d). To further explore the model results, we performed several simulations to predict how visco-elastotaxis would behave if the material properties changed. We set each of the viscoelastic parameters constant and leave the others as the control case (Fig. 6d). Our results show that if we keep 𝑘𝑒or 𝑘𝑒and 𝑘𝑣constant migration direction is reversed. A closer look into the migration velocities indicates that the opposing gradients in 𝑘𝑒and 𝑘𝑣act as competing mechanical stimuli. When one of these parameters is kept constant, the migration direction increases by one order of magnitude, indicating that in the control case, they were acting against each other. When 𝜇was fixed, we observed a mild change in migration velocity. These results indicate that the relaxation time scale is far from the timescale of the adhesion cycle determined by the clutch model. Indeed, when 𝑘𝑒and 𝑘𝑣are fixed, migration velocity decreases by one order of magnitude. 4. Discussion A deep understanding and precise control of cell migration are key for medicine and bioengineering as it may allow us to not only understand how tumor cells invade and metastasize [4,5] or how cells orchestrate tissue regeneration [1,60], but also to engineer cell migration for better biomimetic tissue designs [61,62]. In this work, we have shown through physical arguments how single cells may express viscotaxis, as previously shown by an experimental study. Most living tissues are viscoelastic; therefore, cells migrating across these materials feel a combination of viscoelastic properties. Understanding how gradients in the relaxation properties of the ECM may impact cell migration is a fundamental question in biology, medicine, and bioengineering. Here, we show that the same physical mechanisms that enable durotaxis dictate viscotaxis. Our results show that gradients of the relaxation properties of the ECM may activate an oriented cell migration thanks to an asymmetric adhesion strength, or cell-ECM friction. We also show that specific combinations of viscoelastic properties may induce positive and negative viscotaxis. Specifically, gradients in the viscosity of the ECM always induce positive viscotaxis when the relaxation stiffness is large enough compared to the long-term stiffness, or no migration at all if the relaxation stiffness is negligible because no relaxation of the material occurs at the time scales of the cellular adhesion. However, gradients in the relaxation stiffness show both positive and negative viscotaxis. If the relaxation stiffness is smaller or similar compared to the long-term stiffness, cells will migrate with a weak positive viscotaxis, because the material’s relaxation will also be weak. As the values of the relaxation stiffness increase, positive viscotaxis is enhanced when a large viscosity of the material allows for material relaxation. However, if the viscosity is low compared to the time scale of cell adhesion, the effect of high relaxation stiffness is hindered. We also confronted viscotaxis and durotaxis and found that a diverging minimum in the cell friction makes cells migrate away from it, contrary to what happens with the optimum stiffness in durotaxis [28, 29]. Other mathematical models of durotaxis have also been used to explain how durotaxis emerged [63–66]. This diverging point appears due to the sum of the opposing gradients in the long term and relaxation stiffnesses, or friction. The point at which the minimum of these added contributions appears corresponds to the diverging minimum. Our results also indicated that durotaxis seems to be more efficient in controlling oriented cell migration than viscotaxis. However, a viscous gradient could also shift the direction of durotaxis, or inhibit it. When the relaxation timescale of the material is much larger than the adhesion timescale, or when the relaxation stiffness is low compared to the long-term stiffness, cells mostly sense the nearly elastic long-term Mathematical Biosciences 380 (2025) 109362 7 P. Saez et al. stiffness of the substrate. In these cases, the viscotactic response is weak and cell migration would be due to durotaxis or would not occur. These opposing responses of viscoelastotaxis can be used for designing relaxation gradients of the ECM to either arrest or shift an initially oriented durotactic response (see Fig. 5). Or, the other way around, elastic stiffness gradients can be manipulated to change pure viscotactic responses, as shown in Fig. 5. There is also limitation in our model that should be addressed in the future. The first one is due to the 1D approach we chose, which limits the analysis of shape changes in the cell or the study of migration in more complex environments. Well-established numerical methods could be used in the future to address migration in higher dimensions [67,68]. We also used a simple SLS model. Perhaps, certain ECMs are better reproduced with other viscoelastic models. Moreover, we focused exclusively on mechanical signals and, if they are not present, the cell would just spread symmetrically. However, there are always other signals that may also coexist with the mechanical gradients, and that can be responsible for the natural polarization of the cell. These signals can trigger downstream polarization of GTPases [69] that may self-polarize due to noisy stimuli, which would eventually polarize the downstream signals that, eventually, induce cell migration. Such interaction between mechanical and chemical signals should be addressed in the future. In summary, we have added the effect of a viscous matrix to the notion of durotaxis and rationalized how viscotaxis works. Our model describes for the first time the competition between elastic and viscous gradients of the ECM and establishes a new, or more complex mechanical signal to control cell migration. CRediT authorship contribution statement Pablo Saez: Writing – review & editing, Writing – original draft, Software, Methodology, Investigation, Formal analysis, Conceptualization. Pallavi U. Shirke: Data curation. Jyoti R. Seth: Formal analysis. Jorge Alegre-Cebollada: Writing – review & editing, Conceptualization. Abhijit Majumder: Writing – review & editing, Supervision. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments P.S acknowledges support from the Ministerio de Ciencia, Innovación y Universidades (MCIU), Grant PID2019-11094GB-100 funded by MICIU/ AEI /10.13039/501100011033. The CNIC is supported by the Instituto de Salud Carlos III (ISCIII), MCIU, and the Pro CNIC Foundation and is a Severo Ochoa Center of Excellence (grant CEX2020001041-S funded by MCIU). References [1] P. Friedl, D. Gilmour, Collective cell migration in morphogenesis, regeneration and cancer, Nat. Rev. Mol. Cell Biol. 10 (7) (2009) 445–457. [2] R. Mayor, S. Etienne-Manneville, The front and rear of collective cell migration, Nat. Publ. Gr. 17 (2) (2016) 97–109. [3] S. Van Helvert, C. Storm, P. Friedl, Mechanoreciprocity in cell migration, Nat. Cell Biol. 20 (1) (2018) 8–20. [4] P. Friedl, K. Wolf, Tumour-cell invasion and migration: diversity and escape mechanisms, Nat. Rev. Cancer 3 (5) (2003) 362–374. [5] H. Yamaguchi, J. Wyckoff, J. Condeelis, Cell migration in tumors, Curr. Opin. Cell Biol. 17 (5 SPEC. ISS.) (2005) 559–564. [6] D.T. Butcher, T. Alliston, V.M. Weaver, A tense situation: forcing tumour progression, Nat Rev Cancer 9 (2) (2009) 108–122. [7] R. Langer, D.A. Tirrell, Designing materials for biology and medicine, Nature 428 (6982) (2004) 487–492. [8] U.G.K. Wegst, H. Bai, E. Saiz, A.P. Tomsia, R.O. Ritchie, Bioinspired structural materials, Nature Mater. 14 (1) (2014) 23–36. [9] D. Pantaloni, C. Le Clainche, M.F. Carlier, Mechanism of actin-based motility, Science (80-. ) 292 (May) (2001) 1502–1506. [10] T.D. Pollard, G.G. Borisy, Cellular motility driven by assembly and disassembly of actin filaments, 2003. [11] M.J. Footer, J.W.J. Kerssemakers, J.A. Theriot, M. Dogterom, Direct measurement of force generation by actin filament polymerization using an optical trap, Proc. Natl. Acad. Sci. USA 104 (7) (2007) 2181–2186. [12] C.H. Schreiber, M. Stewart, T. Duke, Simulation of cell motility that reproduces the force-velocity relationship, Proc. Natl. Acad. Sci. USA 107 (20) (2010) 9141–9146. [13] B. Geiger, J.P. Spatz, A.D. Bershadsky, Environmental sensing through focal adhesions, Nature Rev. Mol. Cell Biol. 10 (1) (2009) 21–33. [14] J.T. Parsons, A.R. Horwitz, M.A. Schwartz, Cell adhesion: integrating cytoskeletal dynamics and cellular tension, Nature Rev. Mol. Cell Biol. 11 (9) (2010) 633–643. [15] P. Kanchanawong, G. Shtengel, A.M. Pasapera, E.B. Ramko, M.W. Davidson, et al., Nanoscale architecture of integrin-based cell adhesions, Nature 468 (7323) (2010) 580–584. [16] H. Wolfenson, T. Iskratsch, M.P. Sheetz, Early events in cell spreading as a model for quantitative analysis of biomechanical events, Biophys. J. 107 (11) (2014) 2508–2514. [17] J. Fouchard, C. Bimbard, N. Bufi, P. Durand-Smet, A. Proag, A. Richert, O. Cardoso, A. Asnacios, Three-dimensional cell body shape dictates the onset of traction force generation and growth of focal adhesions, Proc. Natl. Acad. Sci. USA 111 (36) (2014) 13075–13080. [18] A. Shellard, R. Mayor, Durotaxis: The hard path from in vitro to in vivo, Dev. Cell 56 (2) (2021) 227–239. [19] R. Sunyer, X. Trepat, Durotaxis, Curr. Biol. 30 (9) (2020) R383–R387. [20] G. Charras, E. Sahai, Physical influences of the extracellular environment on cell migration, Nat. Rev. Mol. Cell Biol. 15 (12) (2014) 813–824. [21] R.R. Kay, P. Langridge, D. Traynor, O. Hoeller, Changing directions in the study of chemotaxis, Nature Rev. Mol. Cell Biol. 9 (6) (2008) 455–463. [22] P.J.M. Van Haastert, P.N. Devreotes, Chemotaxis: signalling the way forward, Nat. Rev. Mol. Cell Biol. 5 (8) (2004) 626–634. [23] J.A. Espina, C.L. Marchant, E.H. Barriga, Durotaxis: the mechanical control of directed cell migration, FEBS J. (2021). [24] D.E. Koser, A.J. Thompson, S.K. Foster, A. Dwivedy, E.K. Pillai, et al., Mechanosensing is critical for axon growth in the developing brain, Nat. Neurosci. 19 (12) (2016) 1592–1598. [25] C.M. Lo, H.B. Wang, M. Dembo, Y.L. Wang, Cell movement is guided by the rigidity of the substrate, Biophys. J. 79 (1) (2000) 144–152. [26] J.Y. Wong, A. Velasco, P. Rajagopalan, Q. Pham, Directed Movement of Vascular Smooth Muscle Cells on Gradient-Compliant Hydrogels, Langmuir 19 (5) (2003) 1908–1913. [27] S.P. Singh, M.P. Schwartz, J.Y. Lee, B.D. Fairbanks, K.S. Anseth, A peptide functionalized poly(ethylene glycol) (PEG) hydrogel for investigating the influence of biochemical and biophysical matrix properties on tumor cell migration, Biomater. Sci. 2 (7) (2014) 1024–1034. [28] A. Isomursu, K.Y. Park, J. Hou, B. Cheng, M. Mathieu, et al., Directed cell migration towards softer environments, Nature Mater. (2022). [29] P. Sáez, C. Venturini, Positive, negative and controlled durotaxis, Soft Matter (2023). [30] A. Elosegui-Artola, The extracellular matrix viscoelasticity as a regulator of cell and tissue dynamics, Curr. Opin. Cell Biol. 72 (2021) 10–18. [31] E.E. Charrier, K. Pogoda, R.G. Wells, P.A. Janmey, Control of cell morphology and differentiation by substrates with independently tunable elasticity and viscous dissipation, Nature Commun. 9 (1) (2018) 1–13. [32] Z. Gong, S.E. Szczesny, S.R. Caliari, E.E. Charrier, O. Chaudhuri, et al., Matching material and cellular timescales maximizes cell spreading on viscoelastic substrates, Proc. Natl. Acad. Sci. USA 115 (12) (2018) E2686–E2695. [33] O. Chaudhuri, L. Gu, M. Darnell, D. Klumpers, S.A. Bencherif, et al., Substrate stress relaxation regulates cell spreading, Nature Commun. 6 (2015) 1–7. [34] A.R. Cameron, J.E. Frith, J.J. Cooper-White, The influence of substrate creep on mesenchymal stem cell behaviour and phenotype, Biomaterials 32 (26) (2011) 5979–5993. [35] O. Chaudhuri, L. Gu, D. Klumpers, M. Darnell, S.A. Bencherif, et al., Hydrogels with tunable stress relaxation regulate stem cell fate and activity, Nature Mater. 15 (3) (2016) 326–334. [36] S. Nam, O. Chaudhuri, Mitotic cells generate protrusive extracellular forces to divide in three-dimensional microenvironments, Nat. Phys. 14 (6) (2018) 621–628. [37] K. Adebowale, Z. Gong, J.C. Hou, K.M. Wisdom, D. Garbett, et al., Enhanced substrate stress relaxation promotes filopodia-mediated cell migration, Nature Mater. 20 (9) (2021) 1290–1299. [38] O. Chaudhuri, S.T. Koshy, C.B. Da Cunha, J.W. Shin, C.S. Verbeke, et al., Extracellular matrix stiffness and composition jointly regulate the induction of malignant phenotypes in mammary epithelium, Nature Mater. 13 (10) (2014) 970–978. Mathematical Biosciences 380 (2025) 109362 8 P. Saez et al. [39] P.U. Shirke, H. Goswami, V. Kumar, D. Shah, S. Beri, et al., Viscotaxisdirected migration of mesenchymal stem cells in response to loss modulus gradient, Acta Biomater. 135 (2021) 356–367. [40] W. Shu, C.N. Kaplan, A multiscale whole-cell theory for mechanosensitive migration on viscoelastic substrates, Biophys. J. 122 (1) (2023) 114–129. [41] C.E. Chan, D.J. Odde, Traction dynamics of filopodia on compliant substrates, Science (80-. ) 322 (5908) (2008) 1687–1691. [42] A. Elosegui-Artola, R. Oria, Y. Chen, A. Kosmalska, C. Pérez-gonzález, et al., Mechanical regulation of a molecular clutch defines force transmission and transduction in response to matrix rigidity, Nat. Cell Biol. 18 (2016) 540. [43] C. Venturini, P. Sáez, A multi-scale clutch model for adhesion complex mechanics, PLOS Comput. Biol. 19 (7) (2023) 1–27. [44] F. Kong, A.J. García, A.P. Mould, M.J. Humphries, C. Zhu, Demonstration of catch bonds between an integrin and its ligand, J. Cell. Biol. 185 (7) (2009) 1275–1284. [45] N.Q. Balaban, U.S. Schwarz, D. Riveline, P. Goichberg, G. Tzur, et al., Force and focal adhesion assembly: a close relationship studied using elastic micropatterned substrates, Nature Cell Biol. 3 (5) (2001) 466–472. [46] D. Riveline, E. Zamir, N.Q. Balaban, U.S. Schwarz, T. Ishizaki, S. Narumiya, Z. Kam, B. Geiger, A.D. Bershadsky, Focal contacts as mechanosensors: Externally applied local mechanical force induces growth of focal contacts by an Mdia1dependent and rock-independent mechanism, J. Cell. Biol. 153 (6) (2001) 1175–1186. [47] B. Rubinstein, M.F. Fournier, K. Jacobson, A.B. Verkhovsky, A. Mogilner, Actin– myosin viscoelastic flow in the keratocyte lamellipod, Biophys. J. 97 (7) (2009) 1853–1863. [48] K. Larripa, A. Mogilner, Transport of a 1D viscoelastic actin–myosin strip of gel as a model of a crawling cell, Physica A 372 (1 SPEC. ISS.) (2006) 113–123. [49] T. Putelat, P. Recho, L. Truskinovsky, Mechanical stress as a regulator of cell motility, Phys. Rev. E 97 (1) (2018) 1–11. [50] E.L. Barnhart, K.C. Lee, K. Keren, A. Mogilner, J.A. Theriot, An adhesiondependent switch between mechanisms that determine motile cell shape, PLoS Biol. 9 (5) (2011). [51] C.A. Wilson, M.A. Tsuchida, G.M. Allen, E.L. Barnhart, K.T. Applegate, et al., Myosin II contributes to cell-scale actin network treadmilling through network disassembly, Nature 465 (7296) (2010) 373. [52] D. Raz-Ben Aroush, N. Ofer, E. Abu-Shah, J. Allard, O. Krichevsky, et al., Actin turnover in lamellipodial fragments, Curr. Biol. 27 (19) (2017) 2963–2973.e14. [53] A. Mogilner, G. Oster, Force generation by actin polymerization II: The elastic ratchet and tethered filaments, Biophys. J. 84 (3) (2003) 1591–1605. [54] J. Prost, F. Jülicher, J.F. Joanny, Active gel physics, Nat. Phys. 11 (2) (2015) 111–117. [55] M.M. Kozlov, A. Mogilner, Model of polarization and bistability of cell fragments, Biophys. J. 93 (11) (2007) 3811–3819. [56] K. Keren, Z. Pincus, G.M. Allen, E.L. Barnhart, G. Marriott, et al., Mechanism of shape determination in motile cells, Nature 453 (7194) (2008) 475–480. [57] J. Betorz, G.R. Bokil, S.M. Deshpande, S. Kulkarni, D. Rolando, et al., A computational model for early cell spreading, migration, and competing taxis, J. Mech. Phys. Solids (2023) 105390. [58] O.C. Zienkiewicz, R.L. Taylor, The finite element method (Fluid Dynamics), Methods 3 (2000) 347. [59] J. Donea, A. Huerta, Finite Element Methods for Flow Problems, Wiley, 2003. [60] F. Qu, F. Guilak, R.L. Mauck, Cell migration: implications for repair and regeneration in joint disease, Nat. Rev. Rheumatol. 15 (3) (2019) 167–179. [61] D.E. Ingber, V.C. Mow, D. Butler, L. Niklason, J. Huard, et al., Tissue engineering and developmental biology: Going biomimetic, Tissue Eng. 12 (12) (2006) 3265–3283. [62] B. Zhang, A. Korolj, B.F.L. Lai, M. Radisic, Advances in organ-on-a-chip engineering, Nat. Rev. Mater. 3 (8) (2018) 257–278. [63] E.A. Novikova, M. Raab, D.E. Discher, C. Storm, Persistence-driven durotaxis: Generic, directed motility in rigidity gradients, Phys. Rev. Lett. 118 (2017) 1–5. [64] A.R. Hassan, T. Biel, T. Kim, Mechanical model for durotactic cell migration, ACS Biomater. Sci. Eng. 5 (8) (2022) 3954–3963. [65] A.A. Malik, B. Wennberg, P. Gerlee, The impact of elastic deformations of the extracellular matrix on cell migration, Bull. Math. Biol. 82 (2019) 49. [66] A.A. Malik, P. Gerlee, Mathematical modelling of cell migration: stiffness dependent jump rates result in durotaxis, J. Math. Biol. 78 (7) (2020) 2289–2315. [67] W.K. Liu, Y. Liu, D. Farrell, L. Zhang, X.S. Wang, Y. Fukui, et al., Immersed finite element method and its applications to biological systems, Comp. Meth. App. Mech. Eng. 195 (2006) 1722–1749. [68] A. Moure, H. Gomez, Phase-field modeling of individual and collective cell migration, Arch. Comp. Meth. Eng. 28 (2) (2021) 311–344. [69] A.J. Ridley, M.A. Schwartz, K. Burridge, R.A. Firtel, M.H. Ginsberg, G. Borisy, J.T. Parsons, A.R. Horwitz, Cell migration: Integrating signals from front to back, Science 302 (5651) (2003) 1704–1709.