On the relationship between velocities, tractions, and intercellular stresses in the migrating epithelial monolayer
11 On the relationship between velocities, tractions, and intercellular stresses in the migrating epithelial monolayer
Yoav Green , Jeffrey J. Fredberg , and James P. Butler Harvard T.H. Chan School of Public Health, Boston, MA 02115, USA. Department of Mechanical Engineering, Ben-Gurion University of the Negev, Beer-Sheva 8410501, Israel Department of Medicine, Brigham and Women’s Hospital, Harvard Medical School, Boston, MA 02115, USA. * Correspondence: [email protected] The relationship between velocities, tractions, and intercellular stresses in the migrating epithelial monolayer are currently unknown. Ten years ago, a method known as Monolayer Stress Microscopy (MSM) was suggested from which the intercellular stresses could be computed given a traction field. The core assumption of MSM is that the intercellular stresses within the monolayer behave similarly to passive systems like a Hookean solid (an elastic sheet) or a Newtonian fluid (thin fluid film), implying a relation between the displacements/velocities and tractions. Due to the lack of independently measured intercellular stresses, validation of MSM is difficult. An alternative approach, which we give here, is based on simultaneous measurements of the monolayer velocity field and the cell/substrate tractions. With limited assumptions, the velocity field suffices to compute tractions, which we can then compare directly with those measured by traction force microscopy. We find that the calculated tractions and measured tractions are uncorrelated. Since both classical MSM and a purely viscous description of the relation between displacements or velocities and tractions depends on a linear constitutive law, it follows that some modification of these approaches is needed. One possible resolution is the inclusion of an active force. To this end, we give a new relationship between the active force density and the measured velocity (or displacement) field and tractions, which by Newton’s laws, must be obeyed. I. INTRODUCTION.
During wound healing, cancer metastasis, development, and asthmatic airway remodeling, cells comprising a confluent epithelial layer migrate collectively [1–5]. Within the epithelial monolayer, each constituent cell exerts intercellular stresses on neighboring cells, and exerts traction forces [6] on its substrate. While traction forces exerted by a monolayer have been measured for two decades [7–9] , their relationship to measured cellular velocities remains unknown, and the relationship between intercellular stresses and tractions remains unresolved [10].
Ten years ago [11,12], our group suggested a method to recover the induced stresses, , within a monolayer given a measured traction vector field, T . The method was termed Monolayer Stress Microscopy (MSM) [11,12] and we denote the recovered stresses as ( ) MSM T . The principle assumption was that the stresses could be described by the same equations describing simple passive systems such as Hookean solids [13] or Newtonian fluids [14,15], where there is a linear relationship between stress and strain [13], or stress and strain rate, respectively [15]. Since our first work, [11], this approach has been adopted by many, and various derivative formulations have been suggested (Hookean solid [10–12,16–28] or Newtonian fluid [10,18,20,29–31], see recent review [10] for more details). It is clear that this approach is exact in 1D. Because where the divergence of the stress is given by the traction, no assumptions are required regarding the constitutive law in 1D. In 2D, however, the situation is different, and a constitutive law is required. Additional difficulties can be attributed to the effects of the boundary conditions (described below). Tambe et al. [12] addressed the Poisson ratio issue by computing intercellular stresses with a wide variety of Poisson ratios. They found the dependence upon the Poisson ratio to be small, thus lending credence to the MSM approach. On the other hand, it has also been suggested that under certain conditions, including a Poisson ratio of zero, the 2D problem decouples into two 1D (exact) problems. This statement is incorrect. Nevertheless, MSM remains to be validated in the following sense: the tractions are the independent variable, along with boundary conditions at the edge of the monolayer or field of view (see Section II below) and then ( )
MSM T is calculated. But we currently have no direct measures of intercellular stress by which to compare these. An equivalent test validating the use of a passive constitutive law is to note that a measured local velocity or displacement field, together with an assumed viscosity or shear modulus, suffices to compute the tractions. And these in turn can then be compared with the experimentally measured tractions using traction force microscopy (TFM) [7–9]. This comparison is the focus of this work. We find that tractions computed from the measured velocity field are uncorrelated with tractions measured by TFM. This paper is structured as follows. In Sec. II we review MSM used for stress recovery within a monolayer [11]. Since its conception ten years ago [11], many derivative models have been suggested that use either a Hookean solid [10–12,16–28] or a Newtonian fluid [10,18,20,29–31] constitutive equation. We review both formulations and most embedded assumptions of MSM. In Sec. III, we present our alternative approach and show that MSM is an incomplete theory. Sec. IV discusses the implications of our result, which is that the constitutive equation needs modification. One of the suggested modifications is the inclusion of active stresses. A variety of such active stresses formulations have been put forward on an ad hoc basis [10,17,19–21,26–31]. By contrast with ad hoc formulations, here we put forward a fundamental relationship that links active force density to substrate tractions and velocity fields, and is derived with a minimal amount of assumptions. Because it follows directly from Newton’s laws and does not depend on a specific biophysical rheology or molecular mechanism, this relationship is general and robust. Finally, we conclude in Sec. V with a few short remarks. II.
MONOLAYER STRESS MICROSCOPY
In this section, we review MSM as suggested by us and its implementation by others with various derivative models. Sec. II.A presents the general three dimensional (3D) governing equations. Sec. II.B discusses constitutive equations while Sec. II.C discusses incompressibility. In Sec. II.D, we reduce the governing equation to two dimension (2D). Sec. II.E discusses the various boundary conditions (BCs) needed at the edge of the monolayer. A.
3D model
In a 3D bulk where inertia is neglected and there are no body forces [Figure 1(a)], Newtons’s second law is D D , (1) where D are the 3D stresses. In 3D, the 2D tractions, T , are accounted for as BCs [Figure 1(b)] while BCs at the edge of the monolayer are discussed in Sec. II.E. In the remainder of this work, we will refer to T as those measured by traction force microscopy – TFM [8,9] (i.e. the forces the cells apply on the substrate). Eq. (1) is an indeterminate vector of three equations comprising six different terms ( , , , , ) ij i j x y z . Hence, solution of this equation requires additional equations in the form of either a constitutive equation or compatibility equation. Figure 1. (a) 3D schematic of a thin film of fluid (or epithelial monolayer) on a rigid substrate. Not drawn to scale – the length and width of the film is substantially larger than the height. (b) x z cross-section. The top surface is stress free while at the bottom, tractions are balanced by viscous shear stresses. In this work, tractions ( x T and y T ) refer to the measured tractions by the cells on the substrate, cell on substrate T . The dashed lines refer to boundary conditions at the edges discussed in Sec. II.E. B. Elastic and viscous constitutive laws
The original MSM model [11,12] was formulated in 2D. There, the 2D equations comprised three terms ( , , ) xx yy xy . To close these set of equations, the Beltrami-Mitchell compatibility equation was supplemented [11,12]. However, using this equations is tantamount to assuming that the monolayer behaves as a simple Hookean constitutive law [13] characterized by two material constants: shear modulus, G , and Poisson ratio, . In this work, for the sake of simplicity, we consider the case of , which is the limiting case of an incompressible elastic material (i.e. volume conservation – see Sec. II.C). For an incompressible elastic medium the constitutive law [13]
13 3 3 3 3 ( ) ( )
TD D D D tr G
I u u , (2) where D u is the 3D displacement vector field. However, without time-dependent terms, in Eq. (1) the equations for an incompressible Hookean solid and incompressible Newtonian fluid are identical, except for a time derivative in the latter case. This is easily seen by inspection of the constitutive equation for an incompressible fluid [14,15]
13 3 3 3 3 ( ) ( )
TD D D D tr I v v , (3) where D v and are respectively the 3D velocity vector field and viscosity. It is immediately clear that Eqs. (2) and (3) are identical with the identification of displacement with velocity, and with shear modulus and viscosity. How can these two distinctly different physical models give the same prediction? From the experimental perspective, what is typically measured is the displacement of cells, u , between two consecutive frames, say at t and t t . Importantly, if the stress field at t is assumed to be zero, then whether one uses the displacement field, Eq. (2), or the velocity field, Eq. (3), (from / t v u ) the stresses from both approaches are necessarily the same. This addresses the issue of whether a Hookean description can be valid when cells are clearly moving. The answer is yes insofar as for any pair of images, the assumption a stress-free state in the first image yields the exact same stress field as the Newtonian fluid description. Importantly, the fluid description is more general in that does not require an assumption of a zero stress state in any of the image frames, nor does it require an assumption of small strains. Thus both models are essentially identical, here we adopt the viscous formulation as it requires fewer assumptions from the experimental perspective. In what follows, we will use the more conventional notation for the fluid’s constitutive equation ( ) TD D D D p I v v , (4) whereby the pressure, p , takes the place of the trace of the stress tensor, ( ) p tr . Regardless of the chosen experimental interpretation, since the models are identical, both formulations are equivalent and share the discrepancy that will be discussed in Sec.III. Inserting any of these constitutive equations [Eq. (2)-(4)] into Eq. (1) leads to D p v . (5) C. Incompressibility
Cells are comprised mainly of water and as such are essentially incompressible. In either the Hookean solid formulation [13] or Newtonian fluid formulation [14,15], an incompressible material must also satisfy volume conservation. In 3D, this reads as D D v . (6) Nevertheless, cells are not isovolumic insofar water can be transported into or out of the cell via osmotic stress or water channels. Intercellular fluid flow within a monolayer occurs over timescales of hours [32–34]. These water fluxes comprise a source term in the continuity equation [Eq. (6)]. Indeed, it has been shown that osmotic pressure results in cell area and volume oscillations with a period of four hours and amplitude of 20% [32,33]. Further, proliferation and apoptosis, and cell extrusion form the layer represent additional source terms. [26,35]. Because these phenomena are slow (~hours) compared with the time scales of interest here ( t ~minutes), we take the monolayer to be isovolumic, as have others [10,19,20,24,36–40]. D. Reduction from 3D to 2D
In the case of plane stress, the 3D equations for either elastic bodies or fluids can be reduced to thin fluid film in 2D [13]. The reduction is based on the assumption that the height of the monolayer, h , is much smaller than characteristic in-plane lengths and is taken to be uniform [9,10,18]. Here we give a brief derivation in the incompressible case (see chapter 67 of Sokolnikoff [13]). In plane stress, the shear stresses are zero xz yz , (7) and the z component of the normal stress is zero ( zz ). The constitutive law [Eq.(4)] yields , , zz z z z z p v v p . (8) The comma subscript denotes partial differentiation. From 3D incompressibility [Eq.(6)] we have D D z z x x y y D v v v v v , (9) where ( , ) x y is a 2D operator and ( , ) TD x y v v v . Equating Eqs. (8)-(9) yields the 2D incompressibility condition / (2 ) D p v . (10) In contrast to a 3D formulation, where the tractions are imposed as boundary conditions [41], in the plane stress reduction the tractions assume the role of a body force [10–12,16–24,29]. It is common to write the 2D governing equation of force balance as / D D h T . (11) However, it is beneficial to use the constitutive equation [Eq.(4)]. Then, this reads as , , , , ( ) / ij j i i jj j ji i p v v T h , (12) Using Eq. (10), Eq. (12) becomes
32, , , / ij j i i jj i p v T h , (13) or, alternatively, in vector notation
232 2 2 2 / D D D p h v T , (14) Typically, MSM is presented by Eq. (11) along with the incompressible Beltrami-Mitchell equation ( ) /
D xx yy h T . (15) We call these two equations [Eqs. (11) and (15)] the classical MSM formulation, in which the tractions are the independent variable and the stresses are the dependent variable. By contrast, in Eqs. (10) and (14), the tractions are again the independent variable, but now the velocities are the dependent variable, from which the stresses can be calculated. We call these two equations the Stokes formulation. As noted in Sec. II.B it is important to recognize that both formulations are identical. If all the BCs at the edge of the monolayer be given in term of the stresses, then Eqs. (11) and (15) can be solved. However, as we will discuss in the next sub-section, typically the BCs are given in term of the velocities, and as such, Eqs. (11) and (15) cannot be solved without the introduction of an intermediate variable of displacement or velocity through the constitutive law. In strong contrast, Eqs. (10) and (14) can be solved directly and locally for tractions with velocities taken as the independent variable without the need for specific BC assumptions. This will be used below in Sec. III, where the tractions obtained from Eqs. (10) and (14) are used to directly compare with experimentally measured tractions. In the remainder of this work, where we consider the 2D MSM problem and for the sake of brevity, we drop the 2D subscript. E. Edge boundary conditions
Originally, MSM [11,12] estimated intercellular stresses directly from the traction field. Solution of the governing equations requires assumptions about the boundary conditions at the free edges of the cellular domain or at the boundary of the field of view. This latter boundary was termed the “optical edge” [12]. In what follows we will refer to the conditions at both the 0 free and optical edges as “edge conditions”, to distinguish these from the 3D BCs at the top and bottom surfaces discussed in Figure 1. Here we revisit the question of what are the most appropriate edge conditions necessary for a computation of intercellular stresses. We begin with a short review of our previous discussion [12]. Figure 2 shows a schematic representation of an island of cells expanding into free space where four scenarios are considered. The island’s boundary is given by a solid black line. The red circle represents an obstacle. Each square represents a possible field of view, and is color coded for four possible cases. In all cases, we assume velocities (or displacements) and tractions are accessible through measurement.
Figure 2. Schematic of a cell island approaching an obstacle. The cell island is given by a solid black line. Different scenarios, where the cell island is imaged fully or partially, are color coded. See text for case 4 edge conditions.
Case 1 (dashed purple line). The entire island is imaged. In this situation, the edge conditions at the free edges are stress free:
0, 0 n t . (16) Case 2 (dashed-dotted green line). Only part of the island is imaged. For the free edges, the edge conditions are given by Eq. (16), as in Case 1. However, for the optical edges (at the 1 boundaries of the field of view) an assumption must be made. Previously [12], we argued that at the “optical” boundaries the condition is zero normal flux and zero shear stress
0, 0 v n t . (17) Here we suggest that, instead of this ad hoc assumption, it be replaced by using additional measurements of velocities or displacements at the optical edge (e.g. with PIV) such that at the optical edge, meas v v . (18) Other approaches are possible, including requiring zero tractions at the cell free boundary [23,24], with an apparent stress jump at the same boundary, the origins of which remain open. Case 3 (dotted blue line) is solely determined by optical edges, for which Eq. (18) is the appropriate edge condition. Case 4 (long-dashed dotted red line in Figure 2), when cells have contacted an obstacle. Here we must make a distinction between two kinds of obstacles, each of which as two limiting subtypes. The first is a physical obstacle. In this case, the first subtype is no-slip and no flux, for which, at the interface, v . (19) The second subtype is a pure slip (i.e. zero shear stress) with no flux, for which, again at the interface, the edge condition is
0, 0 t v n . (20) The other kind of obstacle is chemical, for example, the boundary between a region coated with collagen and region without. Here, the first subtype is where there are no cellular interactions outside the observed boundary of the island, for which the edge condition is stress-free [Eq. (16)]. The second subtype is where there may be cellular interactions and non-zero tractions 2 outside the visible island. For example, the presence of protruding lamellipodia may exert tractions but may not be visible on a phase image. In this case, it is most appropriate to use zero normal velocity and zero stress [Eq.(20)]. Which condition is most befitting depends on the nature of the obstacle and demands attention to the biology of cell/ECM interactions. Finally, we note the approach of two recent works [30,31]. These works differ from this one, whose starting point is a purely passive constitutive equation, in that their initial assumption was that the stress tensor included an active term [similar what will be discussed in Sec. IV]. The inclusion of their stress term added an additional free parameter that could be determined by applying a tangential and two normal BCs at the free edge
0, 0, 0 t n v n . (21) It should be noted that this approach worked in their quasi 1D geometry of a circular island. However, it remains to be seen if a single fitting parameter is suitable for arbitrary 2D geometries.
III.
MSM VALIDATION
This section is divided as follows. Sec. III.A discusses the inherent difficulties in validating MSM. Sec. III.B presents our new approach, which bypasses all these difficulties and allows for validation. Sec. III.C gives a short review of the experimental results used in this work. In in Sec. III.D we compare theory and experiments, and show that MSM is incomplete and requires further evaluation. The implications of which are given in Sec. III.E. A. Difficulties in validating MSM
A direct validation of MSM would require an independent experimental measurement of stresses, which is extremely difficult. In the Stokes formulation, an independent measurement of stress is no longer required, because the stresses follow directly from Eq. (4). Rather, the validation requires that both the tractions, T , and velocities, v (these will be denoted by meas T
3 and meas v , respectively), be simultaneously measured. Fortunately, both are easily measured: meas T by TFM [7–9] and meas v by PIV, optical flow or other experimental methods. For the Stokes approach, one end goal would be an analytical solution for the velocity given a traction field, ( ) calc v T , in 2D or 3D, and comparing this with meas v . However, this is not without difficulties. First, the governing equation [Eqs. (10) and (14)] are a set of coupled linear partial differential equations. Second, these equations are subject to varying types of BCs (see Sec. II.E). Third, unless the geometry is extremely simple, an analytical solution is typically not possible, thus, necessitating numerical evaluation. In this work, we suggest an alternative approach that negates the need for solving for ( ) v T . However, we did find two ( ) v T solutions: Appendix A provides a derivation of a 3D solution, ( ) D v T , for a monolayer where we assume that the monolayer is infinite in the , x y plane and solve the governing equation via a Fourier Transform. However, due to the complexity of the expressions in k -space, the final calculation requires numerical evaluation. Appendix B provides a derivation of the simpler 2D plane-stress problem, which is typically the most relevant for MSM, we calculated these inverse functions and give an exact solution for ( ) D v T . B. Alternative approach to MSM
Once more, the governing equations in 2D are / 2 p v , (22) / p h v T . (23) Inserting Eqs. (22) into Eq. (23) leads to [13] h v v T t . (24) 4 In classical MSM [11,12,16–25,29], ( ) MSM T , is calculated. In the Stokes formulation, Eq. (24), is solved for ( ) v T [and then ( ) Stokes T is calculated via Eq. (4)]. As previously stated solving for ( ) v T is non-trivial. These difficulties can be entirely circumvented. As T and v are both measured, an alternative approach is to solve for T given v ; denoted as ( ) calc T v . Here, we compute the calculated tractions, denoted ( ) / calc calc meas h t T v [Eq. (24)]. In contrast to ( ) v T , which is difficult and BC dependent, calculating ( ) calc t v is simple; it is algebraic, and in-plane edge conditions are not needed. Also, ( ) Stokes v stresses can be calculated directly (up to a normalization of ) from the velocities from Eqs. (4) and (22). We compared calc t with meas T based upon measurements from a previous work wherein meas v and meas T were obtained independently and concurrently in the same experiment [16]. C. Experimental setup
Briefly, cells were cultured until they were confluent in a confined region bounded by a barrier [16]. Subsequently, the barrier was lifted so that the monolayer could migrate into a cell-free space containing a circular obstacle into which the cells cannot migrate. The velocity field was measured by particle image velocimetry (PIV) [9] and the traction field was measured by TFM [8,9]. To smooth the spontaneous spatial fluctuations that are known to characterize collective cellular migration, velocity and traction fields were averaged over six identical experiments, and these fields were smoothed further using a Gaussian filter with a 2 pixel standard deviation [42,43]. The advancing monolayer encountered and then encompassed the circular obstacle (Figure 3a) [16]. The corresponding velocity components ( , ) x y v v around the circular obstacle 5 resembles that of flow around a cylinder [14], in which case the flow field divides at an upstream stagnation point [Figure 3(b), (c)].
Figure 3. (a) Phase image of the migrating cells encompassing the obstacle. The field of view is 725x725 m . Measured velocities (b) x v and (c) y v (color bar is in units of
10 [ / ] nm s ). Raw data from Ref. [16]. D. Comparison of meas T and calc h t After applying the same Gaussian smoothing filter [42,43] to the measured x -component of traction measx T , the heterogeneity and punctate nature of this field became apparent [Figure 4a)]. The corresponding x -component of traction, calcx ht , is also heterogeneous, but more highly punctate due to its origin in the derivatives of the velocity field [Figure 4(b)]. The values chosen for the height and viscosity are discussed below (see Sec. IV.A. and Appendix C). A pixel-by-pixel scatter plot of these two traction fields, reveals that the peak of measx T is centered around 4 Pa while that of calcx ht is centered around zero (Figure 5). This difference can also be seen in Figure 4(d) and (e) which show the heat maps of | | meas T and | | calc h t , respectively, as well as their vector fields. It is evident that meas T [Figure 4(d)] appears to be dominated by the x -component, measx T , and is pulling away from the free edge [9,16], calc h t [Figure 4(e)] appears to lack a distinct direction and is diffuse. To quantify the degree of correlation, we computed the Pearson correlation coefficient between measx T and calcx t , given by 6 cov( , ) [ ( ) ( )] meas calc meas calcx x x x x T t std T std t . This is the covariance of the measured and calculated tractions divided by their respective standard deviations, and is therefore dimensionless and independent of h . The calculated correlation coefficients are x and y indicating little correlation (the y -component data are given in the Supplementary Material [41]). To ensure that the lack of correlation is not due to the ensemble average or the effect of the interaction with the obstacle we also compared the measured and calculated tractions from a single (non-averaged) dataset pre- and post-interaction with the obstacle [41]. Tractions for these datasets were also uncorrelated. Figure 4.(a) Measured tractions [ ] measx
T Pa used in Ref. [16]. (b) Calculated tractions [ ] calcx ht Pa with h m , and kPa s . (c) The active force density, active meas measx x x hF T ht [see Eq. (28) below]. (d)-(f) A zoomed view of the marked green box from panel (c), showing (d) [ ] meas Pa T , (e) [ ] calc h Pa t and (f) [ ] active h Pa F . The green arrows show the direction of each respective vector of (d)-(f) where the lengths of the arrows are proportional to the magnitude. Figure 5. Scatter plot of [ ] measx
T Pa versus [ ] calcx ht Pa . The color bar is proportional to the count density. E. Implications for MSM
Until now, the foundations of ( )
MSM T have not been verified experimentally because there are no independent measurements of stress. On the other hand, MSM can be tested for internal consistency by a specific comparison of meas T and ( ) calc meas T v . Here we have done that. Our finding of a lack of correlation between meas T and ( ) calc meas T v [or ( ) calc t v ] necessarily implies ( ( )) ( ) calc meas measStokes MSM
T v T . In words, intercellular stresses computed from measured tractions are not the same as intercellular stresses computed from the measured velocity field. This implies that any passive linear formulation, including MSM, is an incomplete formulation. 8 We point out, that by contrast to the 2D situation, stress recovery using MSM in 1D is straightforward and unambiguous, and a constitutive law is not required. For example, a simple balance of forces requires that stresses in the x direction, xx x h T dx . This approach has been successfully used to calculate 1D stresses [9,17,19,20,26–28,30,31] in cases where the symmetry of the geometry (planar or circular) allows, through averaging, a reduction to 1D. More complex geometries require solutions in 2D and the need for compatibility equations arising from a constitutive law. For example, Zimmerman et al [21] recently suggested that perhaps stress recovery was method-independent. In that work [21], they compared 2D stresses recovered from 2D tractions using MSM versus the stresses recovered using the molecular dynamic simulation method of Hardy [44]. Heat maps of recovered stresses in that report [21] show close agreement between the two methods at longer length scales, supporting their claim that stress recovery is method-independent, but at shorter length scales indicate discrepancies that approach 400%. At this scale intercellular stress recovery is strongly method-dependent, and as such, the issue of intercellular stress recovery remains open. IV.
RESOLVING THE MSM PARADOX
The main finding of this work is that tractions calculated from the measured velocity field, ( ) calc meas
T v , do not correspond to tractions that are measured directly, meas T . Hence, the underlying assumptions of MSM require reevaluation. Namely, the constitutive equation must be modified. There are three natural candidates: a linear visco-elastic rheology, a non-linear rheology (visco-elastic or otherwise), or an active rheology. Linear visco-elastic approaches have been adopted in various biological systems (cell aggregate under aspiration [38], cell aggregate spreading [45], single cell twisting via magnetic tweezers [46,47]. However, to the best of our knowledge, such modeling has not been attempted on monolayers. This approach should be considered in future works on monolayer 9 dynamics. This statement, of future consideration, also holds for the even more difficult non-linear rheologies. Here, we consider the last of these candidates – active rheology. We do this because currently this is the most favored approach in investigating monolayers dynamics [10,17,19–21,26–31] and there is some experimental evidence supporting this approach. Below, we show a new relation for the behavior of active stresses that can be derived from Newton’s law and requires a minimal amount of assumptions. A. Active stresses as a resolution
Previous works [10,17,19–21,26–31] considered the possibility that the stress tensor additively comprises both a passive, p , and active, a , term, such that p a . (25) It is then assumed that [10,17,19–21,26–31], the passive term is either an elastic solid [Eq. (2)] or viscous fluid [Eq. (4)]. To the extent that the velocities are representative of the motion of the composite monolayer, we insert p viscous and Eq. (25) into Newton’s law [Eq.(11)] ( ) / a p h T . (26) The active force density, denoted active F , is defined as active a F . (27) From Eqs. (4),(26) and (27) we have active F explicitly in terms of measured quantities: . [ 3 ( )] active meas calcmeas meas meas h hh F T tT v v .. (28) active F includes not only the tractions and velocities, but also the height h and the viscosity . It follows from Eq. (28) that in the limit of low viscosity active F is dominated by the tractions meas T , whereas in the limit of high viscosity active F is dominated by passive component meas ht . 0 Note that while h may be estimated accurately [9,10,18], the same is not true for . Appendix C shows a number of estimates using different methods [41]. Here we illustrate activex F for the case where h m and kPa s , the latter being our best estimate resting on the fewest assumptions, and also falling in an intermediate range such that the contributions of both meas T and meas ht can both be seen [Figure 4(c)]. For these values, and for the data set described above, we used Eq. (28) to compute activex hF [Figure 4(c)], which shows clear contributions from both the traction term [Figure 4(a)] and the passive velocity terms [Figure 4(b)]. For example, it can be observed that activex hF is similar to measx T in that the tractions at the free edge are largely pointing away from the obstacle [Figure 4(a)], as was previously shown [9]. Similarly, there is fine scale structure in activex hF which reflects its origin in the passive component [Figure 4(b)]. Figure 4(f) is a vector map showing the directions active F where the length of the arrows are proportional to active F . Note that the active forces are largely pointing away from the obstacle, similar to the tractions, as was shown in [9,16]. B. Origin of active stresses
Recent experiments have confirmed the existence of active forces in epithelial monolayers. For example, in Ref. [26] it was found that stresses are linearly associated with strains but with a nonzero offset, expressed as a E . This is interpreted as reflecting an active component in both the modulus and the offset [26] and supports the active rheology suggested by Eq. (25). Two recent works also used Eq. (25) in their derivations to explain wetting [31] and fingering instabilities [30] observed in epithelial monolayers. However, they didn’t consider the question of whether or not that the passive formulation is incomplete. With respect to the origin of active forces, velocities, and tractions, several previous works have suggested possible phenomenological forms for the active stress [10]. These 1 potential origins of active F differ substantially in the literature. Some works suggest the origin is in the concentration of active contractile units [19,20] while other suggests in-plane cell polarization [19,20,29–31,48], non-zero tension [26], cell division [26,35]. Importantly, without a compatibility equation for a , integrating Eq. (28) to calculate a , and thus the total stress, , is not possible. The current state of knowledge regarding active F is inconclusive as to which of the above models is most appropriate [10]. Nevertheless, Eq. (28) describing the active stress is derived from Newton’s laws, with a minimal number of embedded assumptions, and must be satisfied. V. CONCLUSIONS
We have shown that, based on a linear passive constitutive law (either elastic or viscous), tractions computed from the measured velocity field are uncorrelated with measured tractions. This implies that MSM as a method is incomplete and requires further evaluation. We suggest that intercellular active stresses may resolve this inconsistency, but that if the active stresses are not derivable from a constitutive law, then it may not be possible to explicitly compute intercellular stresses. However, we do give an explicit expression [Eq.(28)] for the active force density as a function of the measured tractions and velocities. Regardless of the proposed origin of the active stresses, this condition, derived from Newton’s 2 nd law, must be satisfied. ACKNOWLEDGMENTS
We thank Dr. Jae Hun Kim for generously providing us with the raw experimental data used in Ref. [16], and used in this work. We also thank Dr. Bo Lan for helpful discussions. This work was supported by National Institutes of Health (NIH) Grants U01CA202123, PO1HL120839 and T32HL007118. 2
APPENDIX A: 3D SOLUTION FOR ( ) v T A. Governing equations in real space
In 3D, the governing Stokes equations for an incompressible fluid of uniform density are v , (29) p v , (30) where ( , , ) Tx y z v v v v is the 3D velocity vector in Cartesian coordinates and p is the pressure[Figure 1 (a)]. Equations (29) and (30) express conservation of mass and force balance, respectively. For boundary conditions, we take the top (apical) surface, z h , to be stress free ( 0) n [Figure 1(b)]: , , ( ) 0 x z z x z h v v , (31) , , ( ) 0 y z z y z h v v , (32) , (2 ) 0 z z z h v p . (33) The comma subscript denotes partial differentiation. Note that we assume the film height to change only slowly ( h ) such that curvatures terms are negligible [49–51]. At the bottom (basal) surface that is in contact with the substrate, z , the shear stress must be balanced by the tractions, and we take zero normal flux: , 0 / x z xz v T , (34) , 0 / y z yz v T , (35) z z v . (36) 3 B. Governing equations in k space The form of Eqs. (29)-(36) suggests using 2 dimensional Fourier transforms, defined for any function ( ) f r by ( ) ( ) i f d rf e k r k r , integrated over the plane. Here ( , ) T x y r , the wave vector ( , ) T k , k k , | | r r and i is the imaginary unit. With this notation, Eqs. (29)-(36) transform to , x y z z i v i v v , (37) ( ) 0 x zz x i p v k v , (38) ( ) 0 y zz y i p v k v , (39)
2, , ( ) 0 z z zz z p v k v . (40) Note that this approach to solving the field equations is exactly the same as that used in Ref. [9], albeit with different top and bottom boundary conditions. In the current application, the transformed stress-free boundary conditions are given by, for the top surface , ( ) 0 x z z z h v i v , (41) , ( ) 0 y z z z h v i v , (42) , (2 ) 0 z z z h v p , (43) and for the tractions and homogeneous flux at the bottom , 0 / x z xz v T , (44) , 0 / y z yz v T , (45) z z v . (46) C.
3D Solution in k space There are 6 solutions to Eqs. (37)-(40); these are given in Ref. [9]. The two linear combinations of which satisfy the four homogeneous conditions [Eqs. (41)-(43) and (46)] can be written in vector notation cosh[ ( )]( ) cosh 0 k h zz kh , (47) cosh sinh cosh1 2 cosh 2( ) cosh sinh coshsinh cosh sinhsinh 02(cosh ) sinh 0cosh (cosh 2 (cosh ) s kz kz kzk h hkz kz kz kz kzik kz ik kz ik kzkzhk kz kzik kz ik kz hk inh ) kz , (48) where hk hk . The corresponding solutions for the pressure are ik hkp p ik kz kz . (49) In these terms, write the linear combination as , xyz vv A B p A p B pv v , (50) where , A B are chosen to satisfy the inhomogeneous traction conditions [Eqs. (44)-(45)]. Explicitly, this requires // xy TAM TB , (51) where tanh 2tanh 2 khM k kh . (52) 5 This leads to // xy TA M TB . (53) Specifically, coth , 2 y x x y T T T TA kh Bk k . (54) Note that in experiments with transverse averaging (such that y v ), and where the height is spatially uniform, there is a simple volume conservation relation between ( ) z v z h and the height averaged x v at the two edges, say x a and x b . Explicitly, [ ( ) ( )] ( '; ) ' 0 bx x za h v x b v x a v x z h dx . (55) In 2D expanding islands, the equivalent expression of volume conservation is ( , ; ) 0 zC D h ds dxdyv x y z h n v , (56) where C and D are the boundary and domain, respectively. If proliferation is present, these expressions give the volume averaged mean proliferation rate, given measurements of both edge and height velocities. APPENDIX B: 2D SOLUTION FOR ( ) v T A.
2D solution in k space The 2D governing equations are / 2 p v , (57) / p h v T . (58) Solution of Eqs. (57)-(58) for the case of infinite and periodic edge BCs (similar to Appendix A) yields ( ) p p G k T , (59) 6 ( ) 2 p ihk kG k , (60) )43( hk hk k T k Tv . (61) With tractions as the only input, inserting Eq. (59)-(61) into the constitutive law [Eq. (4)] confirms the scaling argument in Refs. [11,12] that the stresses are independent of the viscosity (or shear modulus). B.
2D solution in real space
We can now calculate Green’s function in real space, i.e. the impulse response corresponding to a delta function of tractions at the origin. This is given by the inverse Fourier transform ( ) (4 ) ( ) i f d kf e rk r k . We start with the pressure [Eq. (59)], ii ip p ikr i ed k e d k e d khk hkrdk d e dk J kr dkJ krh k h k h hr kk rr rk kG r G rk . (62) Thus, the pressure is given by the convolution, integrated over all space, ( ) ( ' )' 4 | |' ' p d r h r r rr rT . (63) For the velocity, we identify that we need to calculate three different terms , , k k k [ ( , ) T k ]. The first term is i ii i FT d k e d k ei d k e i ix d k ex xk k kk k k kk r rr rk , (64) 7 To find i kd k e k r , consider the impulse response to the 2D Laplace equation, G r , for which the solution is (2 ) ln r . Its 2D Fourier transform is G k . This implies that (2 ) ln i kd k e r rk from which it is immediate that xFT x r rx rk , (65) k yFT rr , (66)
41 3 2 xyFT rrk . (67) APPENDIX C: VISCOSITY ESTIMATES
Despite the apparent inconsistencies of using only a passive stress tensor, we use the pure fluid dynamical description to estimate the monolayers viscosity. We employ two different independent methods. First, consider the scaling argument that viscous stresses in a thin film balance the tractions, ~ /
Th u . For typical values of ~ 10[ ] T Pa , ~ 10[ / ] u nm s , and ~ 5[ ] h m [9,11,17] one gets ~ 5[ ] ~ 5 10 cell water kPa s with
10 [ ] water
Pa s . Next, we considered a simple coarse grained calculation meas calcrms rms h T t where rms is the root-mean-square average; this gives the result / (2.3 10 ) water O . These values, together with others from the literature are shown in Table 1. Our values are on the lower end of the reported range, but we note that all these estimates are extremely large compared with what we consider physically reasonable. This is another indication that a purely passive viscous description is inadequate, and points to the necessity of active stresses to account for the relationship between velocities and tractions. Table 1. Viscosity estimates from various experimental setups. The first two rows are from this work. Method [References]
10 ( / ) water
Scaling ( / Th u ) Coarse graining ( meas calcrms rms h T t ) Low frequency '' G [46] Micropipette aspiration [38]
Expanding island assays [29] -10 Expanding island – active stresses [30,31] (0.5-5) REFERENCES [1] P. Friedl and D. Gilmour, Collective cell migration in morphogenesis, regeneration and cancer, Nat. Rev. Mol. Cell Biol. , 445 (2009). [2] P. Friedl, E. Sahai, S. Weiss, and K. M. Yamada, New dimensions in cell migration, Nat. Rev. Mol. Cell Biol. , 743 (2012). [3] L. Atia, D. Bi, Y. Sharma, J. A. Mitchel, B. Gweon, S. A. Koehler, S. J. DeCamp, B. Lan, J. H. Kim, R. Hirsch, A. F. Pegoraro, K. H. Lee, J. R. Starr, D. A. Weitz, A. C. Martin, J.-A. Park, J. P. Butler, and J. J. Fredberg, Geometric constraints during epithelial jamming, Nat. Phys. , 613 (2018). [4] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim, B. Gweon, J. Notbohm, R. Steward Jr, S. Burger, S. H. Randell, A. T. Kho, D. T. Tambe, C. Hardin, S. A. Shore, E. Israel, D. A. Weitz, D. J. Tschumperlin, E. P. Henske, S. T. Weiss, M. L. Manning, J. P. Butler, J. M. Drazen, and J. J. Fredberg, Unjamming and cell shape in the asthmatic airway epithelium, Nat. Mater. , 1040 (2015). 9 [5] E. Bazellières, V. Conte, A. Elosegui-Artola, X. Serra-Picamal, M. Bintanel-Morcillo, P. Roca-Cusachs, J. J. Muñoz, M. Sales-Pardo, R. Guimerà, and X. Trepat, Control of cell–cell forces and collective cell dynamics by the intercellular adhesome, Nat. Cell Biol. , 409 (2015). [6] P. Roca-Cusachs, V. Conte, and X. Trepat, Quantifying forces in cell biology, Nat. Cell Biol. , 742 (2017). [7] M. Dembo and Y.-L. Wang, Stresses at the Cell-to-Substrate Interface during Locomotion of Fibroblasts, Biophys. J. , 2307 (1999). [8] J. P. Butler, I. M. Tolić-Nørrelykke, B. Fabry, and J. J. Fredberg, Traction fields, moments, and strain energy that cells exert on their surroundings, Am. J. Physiol. - Cell Physiol. , C595 (2002). [9] X. Trepat, M. R. Wasserman, T. E. Angelini, E. Millet, D. A. Weitz, J. P. Butler, and J. J. Fredberg, Physical forces during collective cell migration, Nat. Phys. , 426 (2009). [10] S. Banerjee and M. C. Marchetti, Continuum models of collective cell migration, in Adv. Exp. Med. Biol. (Springer International Publishing, Cham, 2019), pp. 45–66. [11] D. T. Tambe, C. Corey Hardin, T. E. Angelini, K. Rajendran, C. Y. Park, X. Serra-Picamal, E. H. Zhou, M. H. Zaman, J. P. Butler, D. A. Weitz, J. J. Fredberg, and X. Trepat, Collective cell guidance by cooperative intercellular forces, Nat. Mater. , 469 (2011). [12] D. T. Tambe, U. Croutelle, X. Trepat, C. Y. Park, J. H. Kim, E. Millet, J. P. Butler, and J. J. Fredberg, Monolayer Stress Microscopy: Limitations, Artifacts, and Accuracy of Recovered Intercellular Stresses, PLoS ONE , e55172 (2013). [13] I. S. Sokolnikoff, Mathematical Theory of Elasticity, Mathematical Theory of Elasticity (Krieger Pub Co, Malabar, Fla, 1983). [14] I. G. Currie, Fundamental mechanics of fluids,
Fundamental Mechanics of Fluids (McGraw-Hill, 1974). 0 [15] G. K. Batchelor, An Introduction to Fluid Dynamics,
An Introduction to Fluid Dynamics (Cambridge University Press, 2000). [16] J. H. Kim, X. Serra-Picamal, D. T. Tambe, E. H. Zhou, C. Y. Park, M. Sadati, J.-A. Park, R. Krishnan, B. Gweon, E. Millet, J. P. Butler, X. Trepat, and J. J. Fredberg, Propulsion and navigation within the advancing monolayer sheet, Nat. Mater. , 856 (2013). [17] X. Serra-Picamal, V. Conte, R. Vincent, E. Anon, D. T. Tambe, E. Bazellieres, J. P. Butler, J. J. Fredberg, and X. Trepat, Mechanical waves during tissue expansion, Nat. Phys. , 628 (2012). [18] Y. Kondo, K. Aoki, and S. Ishii, Inverse tissue mechanics of cell monolayer expansion, PLOS Comput. Biol. , e1006029 (2018). [19] S. Banerjee, K. J. C. Utuje, and M. C. Marchetti, Propagating Stress Waves During Epithelial Expansion, Phys. Rev. Lett. , 228101 (2015). [20] J. Notbohm, S. Banerjee, K. J. C. Utuje, B. Gweon, H. Jang, Y. Park, J. Shin, J. P. Butler, J. J. Fredberg, and M. C. Marchetti, Cellular Contraction and Polarization Drive Collective Cellular Motion, Biophys. J. , 2729 (2016). [21] J. Zimmermann, R. L. Hayes, M. Basan, J. N. Onuchic, W.-J. Rappel, and H. Levine, Intercellular Stress Reconstitution from Traction Force Data, Biophys. J. , 548 (2014). [22] V. Nier, S. Jain, C. T. Lim, S. Ishihara, B. Ladoux, and P. Marcq, Inference of Internal Stress in a Cell Monolayer, Biophys. J. , 1625 (2016). [23] T. Das, K. Safferling, S. Rausch, N. Grabe, H. Boehm, and J. P. Spatz, A molecular mechanotransduction pathway regulates collective migration of epithelial cells, Nat. Cell Biol. , 276 (2015). [24] M. Vishwakarma, J. D. Russo, D. Probst, U. S. Schwarz, T. Das, and J. P. Spatz, Mechanical interactions among followers determine the emergence of leaders in migrating epithelial cell collectives, Nat. Commun. , 3469 (2018). 1 [25] H. Jang, J. Notbohm, B. Gweon, Y. Cho, C. Y. Park, S.-H. Kee, J. J. Fredberg, J. H. Shin, and Y. Park, Homogenizing cellular tension by hepatocyte growth factor in expanding epithelial monolayer, Sci. Rep. , 45844 (2017). [26] R. Vincent, E. Bazellières, C. Pérez-González, M. Uroz, X. Serra-Picamal, and X. Trepat, Active Tensile Modulus of an Epithelial Monolayer, Phys. Rev. Lett. , 248103 (2015). [27] M. Uroz, S. Wistorf, X. Serra-Picamal, V. Conte, M. Sales-Pardo, P. Roca-Cusachs, R. Guimerà, and X. Trepat, Regulation of cell cycle progression by cell–cell and cell–matrix forces, Nat. Cell Biol. , 646 (2018). [28] R. Sunyer, V. Conte, J. Escribano, A. Elosegui-Artola, A. Labernadie, L. Valon, D. Navajas, J. M. García-Aznar, J. J. Muñoz, P. Roca-Cusachs, and X. Trepat, Collective cell durotaxis emerges from long-range intercellular force transmission, Science , 1157 (2016). [29] C. Blanch-Mercader, R. Vincent, E. Bazellières, X. Serra-Picamal, X. Trepat, and J. Casademunt, Effective viscosity and dynamics of spreading epithelia: a solvable model, Soft Matter , 1235 (2017). [30] R. Alert, C. Blanch-Mercader, and J. Casademunt, Active Fingering Instability in Tissue Spreading, Phys. Rev. Lett. , 088104 (2019). [31] C. Pérez-González, R. Alert, C. Blanch-Mercader, M. Gómez-González, T. Kolodziej, E. Bazellieres, J. Casademunt, and X. Trepat, Active wetting of epithelial tissues, Nat. Phys. , 79 (2019). [32] S. M. Zehnder, M. K. Wiatt, J. M. Uruena, A. C. Dunn, W. G. Sawyer, and T. E. Angelini, Multicellular density fluctuations in epithelial monolayers, Phys. Rev. E , 032729 (2015). [33] S. M. Zehnder, M. Suaris, M. M. Bellaire, and T. E. Angelini, Cell Volume Fluctuations in MDCK Monolayers, Biophys. J. , 247 (2015). 2 [34] M. Guo, A. F. Pegoraro, A. Mao, E. H. Zhou, P. R. Arany, Y. Han, D. T. Burnette, M. H. Jensen, K. E. Kasza, J. R. Moore, F. C. Mackintosh, J. J. Fredberg, D. J. Mooney, J. Lippincott-Schwartz, and D. A. Weitz, Cell volume change through water efflux impacts cell stiffness and stem cell fate, Proc. Natl. Acad. Sci. , E8618 (2017). [35] J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost, and F. Jülicher, Fluidization of tissues by cell division and apoptosis, Proc. Natl. Acad. Sci. , 20863 (2010). [36] M. H. Köpf and L. M. Pismen, A continuum model of epithelial spreading, Soft Matter , 3727 (2013). [37] S. Pawlizak, A. W. Fritsch, S. Grosser, D. Ahrens, T. Thalheim, S. Riedel, T. R. Kießling, L. Oswald, M. Zink, M. L. Manning, and J. A. Käs, Testing the differential adhesion hypothesis across the epithelial−mesenchymal transition, New J. Phys. , 083049 (2015). [38] K. Guevorkian, M.-J. Colbert, M. Durth, S. Dufour, and F. Brochard-Wyart, Aspiration of Biological Viscoelastic Drops, Phys. Rev. Lett. , 218101 (2010). [39] M. Basan, J. Elgeti, E. Hannezo, W.-J. Rappel, and H. Levine, Alignment of cellular motility forces with tissue flow as a mechanism for efficient wound healing, Proc. Natl. Acad. Sci. , 2452 (2013). [40] M. Basan, J.-F. Joanny, J. Prost, and T. Risler, Undulation Instability of Epithelial Tissues, Phys. Rev. Lett. , 158101 (2011). [41] See Supplemental Material at http://link.aps.org/supplemental/for supplementary information., (n.d.). [42] MATLAB, MATLAB (The MathWorks, Natick, Massachusetts, United States., 2016). [43] To further reduce the noise we applied a Gaussian smoothing kernel with a standard deviation of 2 pixels (‘imgaussfilt’ in Matlab [42]). The derivatives of the velocity were calculated in Matlab using the function ‘gradient’ [42]. 3 [44] R. J. Hardy, Formulas for determining local properties in molecular‐dynamics simulations: Shock waves, J. Chem. Phys. , 622 (1982). [45] S. Douezan, K. Guevorkian, R. Naouar, S. Dufour, D. Cuvelier, and F. Brochard-Wyart, Spreading dynamics and wetting transition of cellular aggregates, Proc. Natl. Acad. Sci. , 7315 (2011). [46] B. Fabry, G. N. Maksym, J. P. Butler, M. Glogauer, D. Navajas, and J. J. Fredberg, Scaling the Microrheology of Living Cells, Phys. Rev. Lett. , 148102 (2001). [47] N. Bonakdar, R. Gerum, M. Kuhn, M. Spörrer, A. Lippert, W. Schneider, K. E. Aifantis, and B. Fabry, Mechanical plasticity of cells, Nat. Mater. , 1090 (2016). [48] S. Banerjee and M. C. Marchetti, Substrate rigidity deforms and polarizes active gels, EPL Europhys. Lett. , 28003 (2011). [49] A. Oron, S. H. Davis, and S. G. Bankoff, Long-scale evolution of thin liquid films, Rev. Mod. Phys. , 931 (1997). [50] R. V. Craster and O. K. Matar, Dynamics and stability of thin liquid films, Rev. Mod. Phys. , 1131 (2009). [51] S. Howison, Practical Applied Mathematics: Modelling, Analysis, Approximation, Practical Applied Mathematics: Modelling, Analysis, Approximation (Cambridge University Press, 2005).(Cambridge University Press, 2005).