Interstitial flows regulate collective cell migration heterogeneity through adhesion

The migration behaviors of cancer cells are known to be heterogeneous. However, the interplay between the adhesion interactions, dynamical shape changes and fluid flows in regulating cell migration heterogeneity and plasticity during cancer metastasis is still elusive. To further quantitative understanding of cell motility and morphology, we develop a theory using stochastic quantization method that describes the role of biophysical cues in regulating diverse cell motility. We show that the cumulative effect of time dependent adhesion interactions that determine the structural rearrangements and self-generated force due to actin remodeling, dictate the super-diffusive motion of mesenchymal phenotype in the absence of flow. Interstitial flows regulate cell motility phenotype and promote the amoeboid over mesenchymal motility through adhesion interactions. Cells exhibit a dynamical slowing down of collective migration, with a decreasing degree of super-diffusion. Our findings, suggest a mechanism of Interstitial flow induced directed motion of cancer cells through adhesion, and provide the much needed insight into a recent experimental observation concerning the diverse motility of breast cancer cells.


INTRODUCTION
The collective migration of cells is relevant for many processes, like tissue remodeling events that underlie embryonic morphogenesis, wound repair, and cancer invasion [1][2][3][4]. Collective cancer cell invasion, followed by local or distant metastasis, is a hallmark of cancer [5]. Cancer cells can employ different cellular and molecular modes of invasion, depending on the specific cell-type. Prominent examples include the autonomous mechanism characterized by stochastic birth-death process, as well as reactive mechanisms induced by the local microenvironment. In such migration, cells move as sheets, strands, clusters or ducts rather than individually, and use similar actin and myosin mediated protrusions, and guidance by extrinsic chemotactic, and mechanical cues as used by single migratory cells [3,[6][7][8][9]. Transition from epithelial to mesenchymal states plays an important role in cancer metastasis, is a subject of contemporary medical importance [6,7,10,11]. Tumor cells reside in an extracellular matrix (ECM) containing interstitial fluid that transports nutrients and signaling molecules. The interstitial flow has been shown to affect the morphology and migration of cells such as fibroblasts, cancer cells, endothelial cells, and mesenchymal stem cells [12]. Recent experiment has demonstrated that interstitial flows promote amoeboid over mesenchymal motility of brest cancer cells using a three dimensional microfluidic model [13].
How cancer cells respond to fluid shear stress (FSS) is less understood. Cancer cells escape from the primary tumor, and travel through blood stream to secondary organ sites [14][15][16]. During this metastatic spread, cancer cells are exposed to a microenvironment of dramatically varying physical forces. Cancer cells in the bloodstream experience hydrodynamic forces, which are alien to cancer cells [17][18][19]. FSS is mainly generated during blood flow in the vascular microenvironment, as well as interstitial flow in the tumor microenvironment [20,21]. FSS could contribute to the metastatic process by enhancing tumor cell invasion, and circulating tumor cell (CTC) adhesion to blood vessels. In this article, we develop a theory to describe how cancer cells respond to fluid shear stress. The cumulative effect arising from the nonequilibrium description of living cells using a stochastic birth-death process, and intercellular hydrodynamic interactions, lead to complex dynamics, which may have far reaching implications in our understanding of tumor growth. One of the major difficulties in the study of collective behavior of the cells far from equilibrium is the breakdown of a fluctuation-dissipation theorem (FDT); hence, independent diagrammatic expansions for the response function and the correlation function. The equilibrium distribution is not known and averages can be computed only for the statistical noise. We study the relevant continuum description of collective behavior of a colony of cells in the long time scale, using stochastic quantization technique, originally proposed by Parisi and Wu [22]. We show that non-linear terms that determine the long range hydrodynamic interactions among living cells dictate the complex collective behavior, when a continuum description of the cellular colony is invoked, in the long time limit. The hydrodynamic interactions lead to directed motions of cells. We find that cells exhibit the mean-square displacement of the cells behaves as t α , with, α = 4/3. In the long time limit the motion becomes ballistic.

II. THEORY
We begin by considering the dynamics of a colony of cells in a dissipative environment where inertial effects are negligible. Each cell experiences mechanical forces, such as adhesion, excluded volume interactions due to neighbors, and a random force characterized by Gaussian white noise. The equation of motion for single cell i is [23] The first term on the r.h.s. of Eq.(1) is the unperturbed velocity profile, where u ∞ is the velocity of fluid far from the cell surface, K P the Darcy permeability of the porous medium, z the height, and e x the unit vector parallel to the x-axis.
The second term on the r.h.s. of Eq.(1) is the term associated with the motility of the cells with local fluid velocity. A recent experiment shows that a fluid force present in the tumor microenvironment modulates gene expression downstream of ROCK-LIMK-cofilin activated YAP1 signaling, resulting in developing fine cytoplasmic extensions or filopodia. Therefore, mortile cells exert force dipole f on the fluid. The fluid exerts equal and opposite force on the cell and the net force is zero at all times. The cell is assumed to be at neutral buoyancy. We assume that at low Reynolds number regime, in the absence of fluid fluctuations, the mortile cell is convected along with the fluid at the local fluid velocity, u(r). The dynamics of the cell is then described by ∂r ∂t = u(r), where r is the position of the cell. The flow velocity of the fluid at r is determined by the solution of the Stokes equation, η∇ 2 u(r) = ∇P + F , where F is the force density by the cell on the fluid, i.e., F = N i f δ(r − r i ). We assume that the fluid is incompressible, i.e., ∇ · u = 0. From these set of equations, the Stokes equation is readily solved to obtain the flow field at the center of the cell as, u(r)=k B T N j=1 µ ij f j , where µ ij is the mobility matrix.
The potential U = U r + U a + U s , determining the configuration-dependent forces experinced by the cells consists of three terms: . U r determines the excluded volume interaction with a range λ, U a determines the favorable attractive interactions between cells with a range σ, and interaction strengths v, and κ respectively. The potential term U s describes the cell-substrate or cell surface-ECM interactions as a function of z 0 , the average distance between cell and wall or ECM. Cells interact with the substrate through receptorligand interactions, described by short range potential. Negative sign associated with the adhesion term, U a , indicates a reduction of hydrodynamic interactions between cells [24][25][26]. To obtain an analytical solution, we use a Gaussian potential to describe adhesion. The third term in Eq.(1) is the effect of force acting on cell j creates a hydrodynamic flow-field in the fluid, thereby entraining cell i. Fourth term in the Langevin equation is due to the spacial variation of the cell's self motilities. It is introduced to compensate the flux caused by the position dependent random velocity contributions η i , which are assumed to be Gaussian random vectors exhibiting hydrodynamic correlations according to fluctuation-dissipation Hydrodynamic effects are incorporated via the mobility matrix µ ij , which is obtained from the Green's function of Stokes equation that satisfies the no-slip condition on a planar boundary [27], with the relative coordinate between cell i and j, r = r i − r j and the relative coordinate between cell i and the image of the cell j, R = r i − r j . The wall Green's function G W contains the bulk Oseen tensor, 8πηG 0 , with solvent viscosity η. The finite radii of the spheres are approximately taken into account via a multipole expansion to the second order, in terms of the cell radius a [28], The self motility tensor µ ii is obtained in the limit r i → r j and regularized far away from the surface by the bulk sphere motility µ 0 = 1/(6πηa).
We consider the evolution of the density function for a single cell can be obtained using standard approach [29]. The time evolution of φ(r, t) is given by, Note that the density equation contains same the information as N-body stochastic Langevin equations. We extend the model phenomenologically by adding the source term that describes both cell division and death, and a noise term that breaks the cell number conservation. Our argument follows from the Doi-Peliti formalism [30][31][32], which was introduced in the context of reaction-diffusion process.
The final Langevin equation, for the time-dependent changes in the density, φ(r, t) is φ denotes the strength of the noise corresponding to number fluctuations, and is a function of concentration, which has been derived using stochastic birth and death processes. This is an out of equilibrium problem characterized by the absence of fluctuation-dissipation theorem due to long range hydrodynamic term. Eq. (4) can be studied an-alytically by treating the non-linear terms using a perturbative approach, based on the stochastic quantization scheme [22,33,34].
We assume that the density fluctuates around a constant value, which simplifies the multiplicative noise term (last term in Eq. (4)). Hence, we can define the density using φ(r, t) = φ 0 + φ 1 (r, t), and write down the equation for density fluctuations in Fourier space as follows: To understand the role of non-linearity, we introduce a change of scale r → sr, φ → s χ φ and t → s z t where χ is the exponent corresponding to cell density fluctuation, and z is the dynamic exponent. The nonlinear term (−iq)µ(q)f φ 1 (q)φ 1 (k − q) involves hydrodynamic interaction, which scales as s 2χ−1 . The other non-linear term (−q · k)µ(q)U (q)φ(q)φ(k − q), which arise due to adhesion to the substrate, and cell-cell interaction, which scales as s 2χ . The non-linear term, (−iq)µ(q)f φ 1 (q)φ 1 (k − q), arises due the force dipole on the fluid, determines the scaling solutions in the long time and long distance limit. The non-linear term (−q · k)µ(q)U (q)φ(q)φ(k − q), may dominate in the intermediate time, if interaction strengths are sufficiently strong. This term determines the scaling in the intermediate time regime. µ(q) scales as q −2 at large scale limit. Therefore in the hydrodynamic limit(i.e. r → ∞ and t → ∞), scaling behavior is determined by the long range hydrodynamic term.

A. Limiting cases
In the absence of hydrodynamic interactions, the nonlinear term −(q·k)U (q)φ(q)φ(k−q) determines the scaling in the intermediate regime, and, the linear terms determine the scaling in the long time limit. We find the caging characterized by sub-diffusive motion of cells in the intermediate time and normal diffusion in the long time limit [35]. In the homeostatic limit, where the death and birth balance, we find normal diffusion.

B. Stochastic quantization approach
To understand the dynamics of collections of cells, we use the stochastic quantization method developed by Parisi-Wu in another context. The collective migration of cells described by Eq. (5) is an out of equilibrium problem characterized by the absence of FDT, which relates the correlation and response function in momentum space as C = 1 w ImG. The usual analytic route to get the scaling solution of this problem, one can introduce a response fieldφ. We need to calculate both the response function (G =< φφ >) and correlation function (C =< φφ >) because of the absence of a fluctuation dissipation relation. The key advantage of the present method is that we do not need to calculate both the correlation, and the response functions. The FDT is constructed in fictitious time introduced in the problem. The FDT relation enables us to obtain the scaling of the correlation function, once the scaling of the response function is known. By taking the infinite limit in fictitious time, one can obtain the correlation function in real time. The scaling solution of the problem can be obtained by power counting analysis instead of doing renormalization group calculation.
We now exploit the Parisi-Wu stochastic quantization scheme [22], and introduce a fictitious time 'τ f ', and consider all variables to be functions of τ f in addition to k and w. A Langevin equation in 'τ f ' space is, This ensures that as τ f → ∞, the distribution function will be given by the action S(k, w), because in the τ f -space a fluctuation dissipation theorem (FDT) is preserved. The action S(k, w) can be obtained by writing down the probability distribution P corresponding to the noise term f φ , and the action S(k, w) in terms of φ 1 (k, w) with the help of Eq.(5). The expression for the S is in the appendix.
The correlation functions, calculated from Eq. (14), lead to the required correlation functions of the original theory, in the τ f → ∞ limit. In order to obtain the scaling laws, it suffices to work at arbitrary τ f . It is obvious from Eq. (14) that in the absence of the nonlinear terms, the Greens function G (0) is given by, where ω τ f is the frequency corresponding to the fictitious time τ f . As is customary, the effect of non-linear terms, can be included perturbatively leading to the Dysons' equation Here, we are concerned with the behavior of Σ(k, ω, ω τ f ), which becomes non-linear when expanded to second order. We note that the contribution comes from two different sources (1) a one-loop contribution from the second order term (containing three φ 1 fields) in Eq. (14) (second term in Fig.2) and (2) a two-loop contribution from the first order term ( containing two φ 1 fields) in Eq. (14) (first term in Fig. 2). The contribution arising coming from the term containing three φ 1 fields, in Eq. (14) can be readily obtained by contracting two of the φ 1 fields. The second order term coming from the one loop contribution in Eq. (14) does not have any new momentum dependance. Hence it is the second-order contribution (first term in Fig.(2)), coming from the two-loop contribution in Eq. (14), which is significant. The correlation function is given by the FDT as C = 1 ωτ f ImG. With these observations, Eq. (15) can be written as, Expanding ν ef f ,D about ν and D 0 , respectively, and noting that the renormalization of ν dominates, we get or, ∆ν = 1 2ν Σ(0, ω, ω τ f )

III. RESULTS
Super-diffusive motion: In the intermediate time, hydrodynamic interactions arising due to adhesion determine the scaling behavior. In a self consistent mode coupling theory, we now replace ν by ∆ν in the self energy term Σ(0, ω, ω τ f ) in Fig.(2), use G as given by Eq. (9), and C, which follows from the FDT. According to scale transformation, we know that ω ∼ k z , ω τ ∼ k 2z , G ∼ k −2z , C ∼ k −4z and the vertex factor V ∼ k z . The self energy term in Fig.(2) can be written as dω τ 2π V V GC. By carrying out the momentum count of Σ(0, ω, ω τ f ), and using ∆ν ∼ k z , we find that Σ(k, ω, ω τ f ) ∼ k d−z . Using Eq. (11), we have k z ∼ k d−z , which leads to z = d 2 . In the context of dynamic renormalization group theoty, the usual approach to avoid infrared divergence in the loop integral, we integrate higher momentum and rescaling the momentum. One can end up with the flow equation for the coupling constant and by finding out the fixed-points; one can calculate the dynamic exponents. We do not consider that to be of any particular significance one way or another. By introducing fictitious time, we are able to bring back the fluctuation-dissipation theorem, albeit in fictitious time, and this allows a power counting analysis for the problem. In order to avoid the infrared divergence, we get the following condition, d − z > 0, implies z < d, therefore, MSD exponent (2/z) > 2/d in 3d. The dynamic exponent we obtain z = 5/2 < 3, consistent with the finiteness condition.
Therefore the single cell mean-square displacement behaves as: In this regime, α = 4/3, implying super-diffusive motion. In the absence of shear flow, the adhesion term determines the caging behavior in the intermediate time, and leads to sub-diffusive motion with MSD exponent, α = 0.8 [35]. measurement of MSD using microfluidic model [13].
In the long time limit, the motility term that arises due to force dipoles on the fluid, determines the scaling behavior. In a self consistent mode coupling theory, similar arguments (as in the intermediate time) shows that z = d−1 2 . The single cell mean-square displacement behaves as t α , with α = 2. Hence, we see that the shear flow driven directed motion in a colony of cells arises due to long range hydrodynamic interactions. The scaling exponents depend on the form of interactions used in the theory. For the analytical solution of the present theory, gaussian potential has been used. Other forms of short range interactions may give different numerical values of the exponents, but the qualitative features of the underlying physics described by the theory, is likely to be preserved.

IV. CONCLUSION
In the present contribution, using a new theoretical framework, we provide insight into the dynamics of a colony of cancer cells driven by shear flow, and birthdeath of cells. The conventional practice in dealing with this out of equilibrium problems is to use a set of fictitious fields called response fields, which provide a field theoretic prescription for the response function. In contrast, we propose the introduction of a fictitious time in which a FDT is valid, and thereby only correlation functions need to be calculated. Our approach greatly simplifies the evaluation of scaling exponents. We see that the nonlinear terms in the density evolution equation, arising from mechanical interactions, determine the scaling behavior in the long time limit. The super-diffusive motion of a colony of cells in 3D is induced by hydrodynamic interactions. A recent experiment shows that the fluid force present in the tumor microenvironment modulates gene expression downstream of ROCK-LIMK-cofilin-activated YAP1 signaling [36]. Hence, the fluid force can trigger cell motility by developing fine cytoplasmic extensions or filopodia, in the highly metastatic PC3 human prostate cancer cell line. The filopodia formation modifies the cellcell interaction and cell-substrate interaction, and hence modifies the hydrodynamic interactions in the present theory. The new theoretical framework introduced here provides evidence of shear flow directed collective motion and could explain the invasion of cancer cells under shear flow, observed in a recent experiment [36]. The theory introduced here could help us understand how cancer cells spread by invading adjacent tissues involved in metastasis [37].
We introduce a fictitious time 'τ f ', and consider all variables to be functions of τ f in addition to k and w. A Langevin equation in 'τ f ' space is, with < f φ1 f φ1 >= 2δ(k + k )δ(w + w )δ(τ f − τ f ). This ensures that as τ → ∞, the distribution function will be given by the action S(k, w), because in the τ -space a fluctuation dissipation theorem (FDT) is preserved. The action S(k, w) can be obtained by writing down the probability distribution P , corresponding to the noise term f φ , and the action S(k, w) in term of φ 1 (k, w) with the help of Eq. (13).