Enhanced persistence and collective migration in cooperatively aligning cell clusters
EEnhanced persistence and collective migrationin cooperatively aligning cell clusters
Vincent E. Debets , , Liesbeth M.C. Janssen , ∗ , Cornelis Storm , ∗ Department of Applied Physics,Eindhoven University of Technology, P.O. Box 513,5600 MB Eindhoven, The Netherlands Institute for Complex Molecular Systems,Eindhoven University of Technology, P.O. Box 513,5600 MB Eindhoven, The Netherlands ∗ Most cells possess the capacity to locomote. Alone or collectively, this allows them to adapt,to rearrange, and to explore their surroundings. The biophysical characterization of such motileprocesses, in health and disease, has so far focused mostly on two limiting cases: single-cell motilityon the one hand, and the dynamics of confluent tissues such as the epithelium on the other. The in-between regime of clusters, composed of relatively few cells, moving as a coherent unit has receivedless attention. Such small clusters are, however, deeply relevant in development but also in cancermetastasis. In this work, we use cellular Potts models and analytical active matter theory tounderstand how the motility of small cell clusters changes with N , the number of cells in the cluster.Modeling and theory reveal our two main findings: Cluster persistence time increases with N whilethe intrinsic diffusivity decreases with N . We discuss a number of settings in which the motileproperties of more complex clusters can be analytically understood, revealing that the focusingeffects of small-scale cooperation and cell-cell alignment can overcome the increased bulkiness andinternal disorder of multicellular clusters to enhance overall migrational efficacy. We demonstratethis enhancement for small-cluster collective durotaxis, which is shown to proceed more effectivelythan for single cells. Our results may provide some novel insights into the connection betweensingle-cell and large-scale collective motion and may point the way to the biophysical origins of theenhanced metastatic potential of small tumor cell clusters. Introduction
Many cell types—even those that otherwise are largelystationary— possess an innate capacity to migrate, in-dividually and autonomously, on two-dimensional (2D)substrates or in three-dimensional (3D) matrices. Prop-erly regulated, cell migration contributes crucially to or-ganismal functioning, as it drives vital processes such asmorphogenesis, tissue formation, wound healing and theinflammatory response. In pathology, cell migration like-wise features prominently, and nowhere more so than incancer metastasis. Cancer remains one of the leadingcauses of death in the developed world [1] and the vastmajority of deaths due to cancer (approximately 90%)are a consequence of metastasis [2–5]. In metastasis, cellsdetach from a primary tumor and invade the surroundingextracellular matrix (ECM)—i.e. the three-dimensionalcellularized material that provides structural support totissue— migrating towards blood or lymphatic vessels.Once there, migratory cancer cells traverse the vesselwall (intravasation) and pass into the circulation systemas circulating tumor cells (CTCs). Eventually, some ofthese CTCs may once again pass the vessel wall and nav-igate the local ECM to seed a secondary tumor [3, 4, 6].The elimination of malignant tumors through earlydetection and timely resection, possibly combined with ∗ [email protected], [email protected] chemo-, radiation-, and immune therapy, is the principaldirective in treatment. Despite the seemingly straight-forward sequence of metastatic events (often referred toas the metastatic cascade ), this process remains poorlyunderstood; effective countermeasures that directly inter-fere with metastasis itself are scarce [3]. The process con-tinues to hold surprises too: Where it was long generallyheld that distant metastases were mostly seeded by sin-gle tumor cells [7], recent experimental studies reveal sig-nificant contributions to metastasis from so called CTC clusters : heterogeneous cell clusters consisting of approx-imately 2 to 20 cells that have collectively detached froma single primary tumor and are collectively undertakingthe entire metastatic cascade; invading, intravasating,and circulating as one conserved unit [2–4, 6–10]. Theseclusters are dangerously potent: A study of spontaneousbreast cancer in mice revealed that over 97% of all ob-served metastases originated from CTC clusters ratherthan single CTCs. Other work highlights the importanceof the cell-cell adhesion mediators such as E-cadherins inmetastasis, and likewise suggest that CTC clusters maypossess a metastatic potential that is at least 50 times(and possibly over a 100 times) greater than for individ-ual CTCs [2, 6, 8, 10, 11]. CTC clusters are associatedto lowered overall survival and lowered progression-freesurvival in a range of cancer types [3]. It was convinc-ingly shown that CTC clusters indeed remain a singleunit throughout the journey from primary tumor to dis-tant site; the pathway in which polyclonal CTC clusters a r X i v : . [ phy s i c s . b i o - ph ] S e p would assemble from single CTCs at some point duringmetastasis is highly improbable [6]. Finally, while collec-tive metastasis is our main motivation, we note that col-lectively moving clusters also play a crucial role in manydevelopmental processes: Refs. [12] and [13] emphasizetheir importance in e.g. the neural crest, in mesoblasts ingastrulation, and in the extension of chick somites form-ing the sclerotome, among a number of other appearancesin biology.Overall, experimental findings have opened up a com-pletely new field of study focusing on relatively small cellclusters in biology [12]. In order to explain in particu-lar the enhanced metastatic potential of CTC clusters anumber of hypotheses have been brought forward, includ-ing the cooperation of heterogeneous cell types within theCTC cluster, shielding from attacks by immune cells, adifferential capacity for sensing and responding to chem-ical gradients [14, 15], and the protection from pressuresand shear forces while in the bloodstream [2, 7, 16]. Yet,much remains unknown about the genesis, transit andthe settlement of CTC clusters during metastasis [2].The purpose of this article is to examine three physical-mechanical aspects of cluster motility. First: clusters areobviously larger than single cells. How do multiple er-ratic individual motile tendencies, with varying degreesof coordinated organization, add up to the collective mo-tion of a small cluster of identical cells? Second: Howdoes in-cluster heterogeneity (in intrinsic motility) affectmotility at the cluster level? And third: How do thesealtered properties affect the ability of a cluster to per-form durotaxis [17–20]; that is - to move directedly inthe presence of a rigidity gradient? The latter has beenshown to improve in large aggregates [17]; how does itplay out in smaller clusters?To address these research questions, we combinecoarse-grained simulations with analytical active mat-ter theory. Specifically, we use the cellular Potts model(CPM) to simulate cell (cluster) motion. This model isaugmented to capture two important features of collec-tive motility: directional persistence and cell-cell align-ment. Directional persistence captures the tendency ofindividual cells to persist directionally for some amountof time [20–23]. It is quantified by a persistence time,which corresponds to the average time it takes a cellto deviate significantly from an initial course. Cell-cellalignment refers to the tendency of densely packed motilecells to mutually inform the direction of their motion[13, 24, 25], and has been invoked to explain collectivemotility in dense systems [26]. This may happen eitherby direct physical interactions such as volume exclusionand traction forces where cell-cell adhesions drag neigh-bors along, or in a more indirect fashion through contactinhibition of locomotion (CIL) [27]. Although the lat-ter tends to cause cells to move away from each other,in dense systems this effect also suppresses convergentrelative motion and is thus generally manifested as aparallel-aligning field. In this work, we model persis-tence and alignment using a Langevin and Vicsek-type [28] approach, respectively. More specifically, we im-plement persistent migration in two and three dimen-sions using a Langevin description for the stochastic ro-tational diffusion of the cells’ instantaneous direction ofmotion. Alignment is implemented in a Vicsek-like feed-back mechanism, and quantified by the relative weightassigned to neighbour velocities when updating the veloc-ity of a given cell. To rationalize our CPM simulation re-sults, we also develop an analytical model for finite-sizedclusters composed of (aligning) active Brownian particles(ABPs), providing more theoretical insight into the clus-ter migration efficacy as a function of cluster size andcluster heterogeneity.The paper is organised as follows. We start with in-troducing the CPM and demonstrate how a persistentrandom walk and a Vicsek cell-cell alignment term areimplemented in the CPM. We validate the implementa-tion of persistence and alignment by analyzing the tra-jectories of single cells and cell clusters exploring homo-geneous environments. We then discuss the theory ofABPs, and use it to provide an analytical underpinningof the numerically observed behaviors. Finally, we relatethe enhanced persistence of clusters to cell transport ina more complex, durotactic environment. We concludeby summarising the main findings and provide some fu-ture directions and topics where our results may have animpact. Cellular Potts Model
To simulate the motion of a CTC (cluster) throughthe ECM we employ the so-called cellular Potts model(CPM) [29, 30]. This model, credited for explicitly rep-resenting the cell shape, has been successfully appliedto a wide variety of biological phenomena involving forinstance blood vessel network formation, cancer cell in-vasion, and collective cell motion [25, 31–35]. The CPMis a variation on the classic Potts model [36] and consistsof integer spins σ ( x ) ≥ a ), whose sites are char-acterised by their position in space x . Biological cells arerepresented as (simply connected) domains of equal spin σ ( x ) >
0, while the medium or ECM is assumed to behomogeneous (to exclusively focus on cell-cell alignment)and depicted by σ = 0 (see fig. 1).Cell movement can then be imposed on the system viaa modified Metropolis Monte Carlo algorithm [37–39]: acandidate lattice site i is randomly chosen and its spinvalue σ ( x i ) is attempted to be changed to a randomlypicked adjacent spin value σ ( x j ). The attempt is (pro-vided all cells remain simply connected) accepted witha Boltzmann probability, i.e. min(1 , e − ∆ H /T ), where T parameterizes the energy associated with membrane fluc-tuations [29]. The parameter is suggestively called T toemphasize the temperature-like role it plays in tuningdynamics from quiescent to actively disordered, but we A Experiment
Multicolored collective invasion A µ m mTomato mGFP DAPI CD E
Multicolored disseminated tumor cell cluster mTomato mGFP B Dissemination Model Outcome
Tumors invade collectivelybut obligately disseminateas single cellsTumors both invade anddisseminate collectively mTomato mGFP F-Actin
Mono-clonalMulti-clonal
Multicolored intravascular emboli xyxzyz mTomato CFP DAPI VE-Cadherin + CD31
LumenEmbolus20 µ m F Circulating tumor cell clusters
EcadCollective Invasion IntravascularEmbolusDisseminatedTumor Cell Cluster mTomatomGFPDAPI
E-cadherin+ multicellular clusters
Multicolored Red Only Cyan Only mTomatoCFP
DAPIK14 G Cluster size (
10 60 o f E ve n t s Cluster size ( K + ce ll s ( % o f t o t a l ) H µ m µ m 20 µ m10 µ m10 µ m Fig. 2.
Direct observation of polyclonal collective invasion, polyclonal disseminated tumor emboli, and polyclonal CTC clusters. ( A ) Representative micro-graph of polyclonal collective invasion arising from mosaic ROSA mT/mG ;MMTV-PyMT transplanted tumors stained with phalloidin (F-actin) and DAPI ( n =
75 units, 5 tumors). ( B ) Schema of two potential outcomes for disseminated tumor cell clusters at the tumor stromal interface. ( C ) Representative micrographsof a polyclonal disseminated tumor cluster (yellow arrow) in the x - y plane with successive images along the z axis ( Left panels) and reconstructed 3D image(
Right ) ( n =
25 units, 5 tumors). ( D ) Representative micrograph of a polyclonal disseminated tumor embolus contained within a vessel. The transplanted tumoris composed of mTomato + and CFP + tumor cells. Injection with VE-Cadherin and CD31 fluorescently labeled antibodies marked functional vasculature.( E ) Representative micrographs demonstrating E-cadherin + polyclonal collective invasion, dissemination, intravascular embolus (from left to right). Yellowhash marks: vessel lumen. ( F ) Representative micrographs of CTC clusters composed of mTomato + and CFP + tumor cells and stained for K14 and DAPI ( n = n =
13 mTomato + clusters, n = + cluster). ( G ) The number of events for each CTC cluster size is presented as a histogram ( n =
134 events, 3 transplanted mice). ( H ) The median percentage of cells that are K14 + in CTCs of different cluster sizes are presented as a boxplot ( n =
17 clusters).[Scale bars, 2 mm ( A , Left ), 40 μ m ( A , Right ), 20 μ m ( C – E ), and 10 μ m ( F ).] Cheung et al. PNAS | Published online February 1, 2016 | E857 C E LL B I O L O G Y P N A S P L U S B CPM C Trajectories
FIG. 1. (A) Experimental observation of a circulating tumor cell cluster (CTC cluster) taken from [10]. (B) Visualisation in2D of CPM cells representing the experimentally observed CTC cluster. Different colors represent different spins σ with themedium ( σ = 0) shown in white. (C) Example trajectories of an aligned four-cell CPM cluster with the black dots denotingthe end point of each trajectory. Inset shows the respective cluster with the arrows denoting the polarity vector p σ of each cell. stress that it is not an actual temperature. ∆ H is the re-sulting change in a phenomenological Hamiltonian H ; thelatter accounts for all physically relevant terms, which inthe original model are (approximately) constant cell vol-ume and a finite cell-cell interfacial tension yielding aHamiltonian [29, 30] H = H volume + H adhesion = λ (cid:88) σ ( V σ − V σ, ) + (cid:88) i,j J σ ( x i ) ,σ ( x j ) (1 − δ σ ( x i ) ,σ ( x j ) ) | x i − x j | . (1)Here V σ denotes the volume (in 3D) or area (in 2D) ofcell σ , i.e. the number of lattice sites with σ ( x ) = σ ; V σ, is the preferred volume or area of the corresponding cell.The parameter λ represents the strength of the volumeconstraint and the first sum is taken over all cell spins σ >
0. For the adhesion term the sum is taken over near-est and next-nearest neighbours i, j with J σ,σ (cid:48) (= J σ (cid:48) ,σ )denoting the adhesion coefficient between cell σ and cell σ (cid:48) (or the medium σ (cid:48) = 0), | x i − x j | the distance be-tween the neighboring sites, and δ σ,σ (cid:48) the Kronecker deltawhich ensures that only lattice site pairs of different cellscontribute to the surface energy. Note that generally J σ,σ (cid:48) >
0, and that by choosing different values for thecoefficient J between two cells and between a cell andthe medium we may implement preferential cell-cell ad-hesion. To quantify the evolution of the system withinthe CPM we introduce the Monte Carlo step (MCS) as atime measure [29, 30, 38]. The MCS is defined as N l at-tempts to change a spin value, with N l the total numberof sites in the lattice; it ensures that on average each lat-tice site is updated once every MCS, thereby decouplingthe time step from the actual system size [38]. Activity & Persistence
In its original formulation, i.e. eq. (1), cell dynamics inthe CPM arises solely from fluctuations in the cell volumeand interfacial area (or cell area and interfacial length in2D systems). As a result, the cells do not experience anydirectional bias. In real life, however, cells migrate ac-tively, and may exhibit biased, directional motion. Thismay be because of external guiding cues such as the localorganization of the extracellular matrix, and more gener-ally in response to gradients of some kind. In such cases,the motion is called a taxis . The most well-known ofthese tactic motions is chemotaxis, in which cells moveupstream in gradients of beneficial compounds such asnutrients or oxygen. To implement such directed mo-tion, which we assume in one form or other to feature inCTCs migrating through the ECM, an additional energybias ∆ H a is incorporated in the change of the Hamilto-nian ∆ H . This bias promotes attempts that move thecell along a preferred direction which we shall call the polarization , and is given by [25, 39, 40]∆ H a = − (cid:88) σ = σ ( x i ) ,σ ( x j ) κ σ ∆ R σ ( σ ( x i ) → σ ( x j )) · p σ . (2)Here, p σ denotes the (unit) polarization vector of cell σ ,i.e. the direction in which the cell is currently moving.∆ R σ ( σ ( x i ) → σ ( x j )) is the center-of-mass displacementof cell σ that would result if the proposed move wereaccepted, and κ σ > σ . Isolated cells in experiments generally exhibit persis-tent motion. That is, the direction of motion drifts onsome characteristic timescale. Indeed, it has been shownthat single cell motility in 2D can be accurately de-scribed by a persistent random walk (PRW) [21, 22], andthat consequently its mean square displacement (MSD)is given by [23, 41–43] (cid:10) ( r ( t ) − r (0)) (cid:11) = 4 D a τ ( e − t/τ + t/τ −
1) + 4 D t t, (3)where r ( t ) is the position of a cell at time t . The dis-placement of the cell is comprised of two parts; a purelydiffusive part characterized by a passive diffusion coef-ficient D t , and a persistent contribution quantified bya persistence time τ and an active diffusion coefficient D a ≡ v τ / v the active cell speed). At very shorttimes t (cid:28) τ , the resultant motion is diffusive (MSD ∝ t )with diffusion coefficient D t , ballistic (MSD ∝ t ) at in-termediate time scales t ≈ τ , and diffusive again withan enhanced diffusion coefficient D t + D a at long times t (cid:29) τ .Although the PRW accurately describes cell motility in2D, the correct description in 3D involves an anisotropicpersistent random walk model where two persistent ran-dom walks for a primary and nonprimary direction ofmotion of the environment are combined [22]. To makeour general point, and to facilitate comparison with ananalytical model we will present later on in this paper, wemostly restrict our simulations to 2D cell (cluster) mo-tion. Nonetheless, we note that an extension to 3D leadsto similar results (see appendix A). Interpreting our po-larity vector as the instantaneous direction of motion ofour active cell, we impose a PRW by letting p σ undergorotational diffusion. This is implemented by expressingit in terms of its polar angle φ σ : p σ = [cos( φ σ ) , sin( φ σ )],and having it evolve in time according to a discretisedangular Langevin Dynamics process [23, 41–44]: φ σ ( t + ∆ t ) = φ σ ( t ) + (cid:114) τ σ Γ( t ) . (4)Here ∆ t is the time step of the update which we set to1 MCS, τ σ is the implemented persistence time of cell σ (given in units of MCS), and Γ( t ) is a stochastic whitenoise term with zero mean, (cid:104) Γ( t ) (cid:105) = 0, and a varianceequal to ∆ t ; (cid:104) Γ( t )Γ( t (cid:48) ) (cid:105) = ∆ t δ ( t − t (cid:48) ). Vicsek Alignment
With the update scheme given by eq. (4), we have in-corporated the persistent random walk into the CPMthrough reorientations of the polarity vector p σ . Thisvector represents the currently preferred direction of mo-tion, and may be interpreted as an internal polarizationof the motile machinery, i.e. the instantaneous polariza-tion direction of cytoskeletal stress fibers or, in a morepragmatic sense, as the orientation of the leading edgeof the cell [25] (even though the direction of movementdoes not always line up perfectly with either of these twodirections). For now we treat p σ as a proxy for some in-ternal or external bias direction that guides the motion.This brings us to a principal feature of this work: theeffect of cell-cell alignment. As detailed in the introduc-tion, it is reasonable to assume that when cells are in 𝒑 " 𝒑 " 𝒑 " 𝒑 " FIG. 2. Visualisation of the polarity vector p σ of a CPM cell σ (colored in yellow) and its direct neighbor polarities p σ (cid:48) with which cell σ aligns according to a Vicsek-type model.Dashed arrows denote the polarity vectors of the cells in thecluster that are not in direct contact with cell σ . contact with each other, they may influence each oth-ers’s direction of motion and thereby alter the polar-ity vector p σ of nearby cells. Inspired by the capacityfor neighbor-induced migration alignment in the contextof, for instance, wound healing [26], we let this interac-tion between cells manifest itself as a tendency to aligntheir respective polarities in parallel fashion. This is im-plemented numerically using an adaptation of the well-known Vicsek model [28, 45], extending the update rulepresented in eq. (4) to φ σ ( t +∆ t ) = arg (cid:32) γ p σ ( t ) + (cid:88) σ (cid:48) p σ (cid:48) ( t ) (cid:33) + (cid:114) τ σ Γ( t ) , (5)where γ is a weight factor that controls the degree ofalignment, arg( a ) denotes the angle of a vector a in po-lar coordinates and the sum is taken over all cells σ (cid:48) thatare in direct contact with cell σ . We shall call two cellsin direct contact when there is at least one site with spin σ on the one cell that shares a boundary with the latticesite of a nearby cell σ (cid:48) (see fig. 2). Note that this defini-tion of neighbourhood slightly deviates from the originalformulation of the Vicsek model, which aligns all par-ticles within an interaction radius. When γ → ∞ thealignment disappears (all direct contacts carry zero rel-ative weight in the update scheme, recovering eq. (4)),whereas for γ = 1 the polarity vector p σ instantaneouslytakes on the average direction of itself and its neighborsafter each update and we have perfect local alignment.We will call this the fast-aligning regime, and it is alsohow the alignment in the original Vicsek model is imple-mented [28, 45]. Simulation Details
Each simulation starts with initiating a model CTC clus-ter by placing N cells square cells—domains of equal sizeand each with a unique spin σ > a . The systemis then equilibrated by running the CPM simulation for500 MCS including only the original Hamiltonian eq. (1),without the active energy bias eq. (2). This is done toallow the cells to relax to a natural, smoothly convexshape. After this equilibration stage we assign polarityvectors drawn from a uniform distribution to each cell,set the cluster center of mass to R c ( t ) ≡ t = 0. We then runthe the actual simulation using the Hamiltonian eq. (1)including the active energy bias eq. (2). We proceed totrack the cluster (or single cell) center of mass R c ( t n ) atfixed time intervals ∆ t = t n +1 − t n = 1 MCS to generatethe motile trajectory of the cluster (see fig. 1 for exampletrajectories).For now we will assume all cell parameters to be spa-tially independent, and each cell to be identical. That is,we set κ σ = κ , τ σ = τ , V σ, = V , J σ,σ (cid:48) = J cell − cell for σ, σ (cid:48) >
0, and J σ,σ (cid:48) = J cell − substrate for σ ∨ σ (cid:48) = 0. Addi-tionally, for all 2D simulations in this work we have fixedthe simulation pseudo-temperature T = 1, the targetarea of the cells V = 64 a , the area constraint strength λ = 1, the cell-cell line tension J cell − cell = 0 .
5, and thecell-substrate line tension J cell − substrate = 1. The posi-tive difference between J cell − substrate and J cell − cell impliescells prefer boundaries with other cells over boundarieswith the substrate. Furthermore, the fact that both val-ues are individually positive implies that all boundariesexperience a positive (contractile) line tension. Thus,effectively, this choice of J ’s encodes both cell-cell ad-hesion and cortical tension. Combined with an activeenergy bias and cell persistence time of typically κ = 5and τ = 500 MCS, respectively, these parameters en-sure that cells tend to stick together and, by mapping a ∼ µ m and MCS ∼ . ∼ µ m, a speed of ∼ µ m / h, and a persistence time of ∼ H r . This bias term forces the cellsto prefer a circular shape, thus penalizing e.g. fingering-type structures. In principle, this would be taken care ofby the positive cortical tension, but for small systems wefind that lattice effects on the shape are non-negligible.The ∆ H r term may be physically interpreted as a bend-ing rigidity of the cell cortex/perimeter, and we have ver-ified that its precise value does not influence our mainfindings (see appendix B) for more details). Migration in Uniform Environments: Phe-nomenology
Single Cell Motion
Before proceeding to the effects of cell-cell alignment inclusters, we first validate our implementation of the PRWinto the CPM. For this we study the MSD of isolatedcells for different values of τ and κ . To analyze the dif-fusive process we plot the instant diffusion coefficients D in = MSD / t that follow from the calculated MSDsand fit the results with a PRW [eq. (3)]. This is demon-strated in fig. 3. The accurate fit confirms that indeedthe MSD follows a PRW. By plotting the instant diffusioncoefficient we clearly recognize the transition from an ini-tial ’slow’ diffusive process (constant D in = D t ) via anintermediate ballistic regime (manifesting as D in ∼ t ),to again a diffusive process with an increased diffusioncoefficient D in = D t + D a in the long time limit.We further characterise the cell motion from each fit-ted MSD by extracting the persistence time τ p (we addthe subscript to distinguish the fitted persistence timefrom the implemented one; we will do this throughout),the active diffusion coefficient D a (or an average activecell speed v ), and the ’thermal’ diffusion coefficient D t .The resulting values of these parameters are shown as afunction of both τ and κ in fig. 3. These results confirmthat the persistence time τ which we implement in theCPM is also the time observed in the simulation, i.e. by τ p , and that it is independent of the polarity strength κ . Also, as anticipated, the active diffusion coefficient D a scales linearly with τ and quadratically with κ . Tak-ing into account that for a PRW D a ∝ v τ , we concludethat τ p = τ and v ∝ κ . This implies that indeed wedial in the persistence time and active speed of individ-ual cells directly with τ and κ . The observed values for D t remain constant upon changing both τ and κ . This,too, is as expected: The passive motion originates fromthe pseudo-thermal fluctuations in cell area and shapeeffected by the parameter T . This intrinsic randomnessis completely independent of all other parameters.Thus, consistent with earlier work in e.g. [25] and [39],we have demonstrated that an active PRW can effectivelybe mapped onto the CPM. In a broader context than cellmotility we note that this implementation also providesa good framework for the study of a larger range of ac-tive (soft) materials, and is not necessarily limited to adescription of biological cells. Aligned Cell Cluster Motion
We now turn to the collective motion of CTC clusters,focusing specifically on the role of cell-cell alignment andcluster size. Let us first consider the case of fast align-ment ( γ = 1) for a cluster of N cells identical cells. Theresults for D in , extracted from the corresponding MSDs, Time (MCS) -2 -1 D i n = M S D / t( a / M C S ) τ = 200 τ = 500 τ = 1000 τ = 2000 τ = 5000 Time (MCS) -2 -1 D i n = M S D / t( a / M C S ) = 2= 4= 5= 6= 8 τ (MCS) -2 D ( a / M C S ) , τ p ( M C S ) D a Simulation D t Simulation τ p SimulationTheory / Fit -2 D ( a / M C S ) , τ p ( M C S ) D a Simulation D t Simulation τ p SimulationTheory / Fit
A BC D
FIG. 3. (A-B) Plots of the instant diffusion coefficients D in = MSD / t (markers) which have been calculated from 2D singlecell CPM simulations and are fitted with a PRW using eq. (3) (lines). Results are obtained for (A) different implementedpersistence times τ with a fixed value of κ = 5 and (B) different polarity strengths κ with a fixed value of τ = 500 MCS.(C-D) Persistence time τ p , active diffusion coefficient D a , and thermal diffusion coefficient D t obtained from the fits shown in(A-B) respectively. Red dotted lines denote the predicted or fitted values corresponding to τ p = τ and κ ∝ v . Data has beenobtained by time-ensemble averaging over 50 trajectories each consisting of 50000 MCS. are shown in fig. 4. The MSDs are well-defined, indicat-ing sufficiently large sample sizes, and are still accuratelyfitted with a PRW; we conclude that an aligning cluster,too, moves according to a PRW. The question, then, ishow its parameters depend on cluster size and alignment.The resulting fit parameters, plotted as a function of N cells in fig. 4, allow us to extract these dependencies. Weobserve that the ’thermal’ diffusion coefficient decreaseswith the number of cells. This can be attributed to theincreased size of the cluster that results in weaker relativefluctuations in shape and size. A power law fit yields D t ∼ N − . , showing that this decrease is roughly linearwith the number of cells. However, note that these fits(and the subsequent ones) are only made over one decadein N cells and so the significance of the fitted powers islimited.The persistence time of the cluster, on the other hand,is seen to increase linearly with the number of cells: τ p ∼ N cells τ . An intuitive explanation for this may befound in the fact that when the cells are strongly align-ing (sufficiently small γ ), (almost) all cells must simul-taneously reorient towards the same direction in orderto permit the entire cluster to change its course. This suggests that larger clusters of cells prone to alignmentgenerally continue to move along the same direction forlonger times, which corresponds to an increasing persis-tence time.Similar to the persistence time, the active diffusion co-efficient D a also increases linearly with N cells . We mayunderstand this scaling by noting that for a PRW, theactive diffusion coefficient is expected to scale linearlywith the cluster persistence time D a ∝ τ p and thus, byextension, also with N cells . This holds exactly when thecluster speed is independent of the cluster size, and re-mains the same as that of a single isolated cell. We do see,however, that cell clusters actually have a slightly largeractive speed than single cells. This is demonstrated in theinset of fig. 4, which shows the average absolute velocity v as function of N cells . However, a power law fit yields v ∝ N . ; hence there is only a weak dependence of thecell cluster speed on the number of cells that does notstrongly influence the persistent motion of the cluster.Interestingly, comparable results are obtained for 3Dsimulations (see appendix A). This suggests that the ef-fect of the Vicsek alignment of cell polarities on the cellcluster motion is the same in 3D as it is in 2D. In par- Time (MCS) -4 -3 -2 -1 D i n = M S D / t( a / M C S ) N cells = 1 N cells = 2 N cells = 4 N cells = 9 N cells -4 -2 D ( a / M C S ) , τ p ( M C S ) D a Simulation D t Simulation τ p SimulationTheory / Fit N cells -2 -1 v ( a / M C S ) SimulationFit
A B
FIG. 4. (A) Plots of the instant diffusion coefficient D in = MSD / t for fast-aligning ( γ = 1) 2D cell clusters consisting ofa variable amount of N cells identical cells (markers). The results have been fitted with a PRW fitted using eq. (3) (lines).Inset: average cell (cluster) velocity v as a function of N cells including a power law fit. (B) Cluster persistence time τ p , activediffusion coefficient D a , and thermal diffusion coefficient D t obtained from the fits shown in (A). Red dotted lines show a powerlaw fit for D t and a comparison and fit of respectively τ p and D a to the derived results of the fast-aligning active Brownianparticle theory, i.e. eq. (13). Simulation parameters used: κ = 5 and τ = 500 MCS. Data has been obtained by time-ensembleaveraging over 30 trajectories each consisting of 50000 MCS. ticular, it shows that within a more extended 3D modelsetup, the effect of alignment still results in a linear de-pendence of the persistence time on the number of cellsin a cluster.Up until this point, we have imposed fast neighboralignment within the cluster by setting γ = 1. In or-der to assess the influence of the relative weight of theneighboring polarizations, we have calculated MSDs forcell clusters experiencing weaker alignment (by setting γ = 50). For this γ , the MSDs and resulting fit pa-rameters do not change noticeably; we find roughly thesame results as for γ = 1 [47]. Thus, even moderatedegrees of alignment still produce highly cooperative mi-gration in the cluster. There is a finite bound on this ef-fect, however—by increasing γ further, the system entersa regime where the alignment is not sufficiently stronganymore and the cluster may even disintegrate. Indeed,our simulations show that when γ passes a critical value(roughly, of the numerical order of the persistence time τ ), the cluster of cells quickly falls apart into single cellsand an analysis of the center-of-mass trajectory becomesmeaningless. We can understand this by realising that inthe case of no or weak alignment cells often want to travelin different directions (opposite polarity vectors) for longtimes and can then actively pull themselves loose fromthe other adjacent cells. Furthermore, the disintegrationof the cluster suggests that there exists a critical degree ofalignment γ c beyond which alignment is not sufficientlystrong to keep the cluster together.Thus, we have demonstrated that fast alignment of thecells in the CPM will increase the persistence of the clus-ter, allowing it to move more directionally. This happensat the cost of a decrease in the translational diffusion coefficient. In the case we consider, and which we as-sume to most closely represent actual cellular behavior,the overall motion is dominated by its active contribu-tion; ( D a (cid:29) D t ). As a result, the decrease in the ’ther-mal’ diffusion coefficient will hardly influence the overallmotion. This implies that aligned clusters can, on aver-age, cover more distance than single cells within a giventimeframe, provided that v is sufficiently large. In thecontext of CTC clusters, this might allow them to reachtargets such as blood vessels more easily. It also suggeststhat when clusters experience an externally imposed po-larity (through e.g. tracks or anisotropy in the ECM),they are generally better able to follow such tracks col-lectively compared to single cells, enhancing the effectsof contact guidance.Finally, we consider the effect of fast alignment for aheterogeneous cluster. In particular, we investigate howone less persistent cell influences the motion of an oth-erwise more persistent cluster. Indeed, individual cellstypically show a variety of persistence times in experi-ments [22], and thus CTC clusters will consist of a het-erogeneous mixture of cells. To study this effect, we sim-ulate the motion of a fast-aligning cluster consisting of N cells = 4 cells, 3 of which have a ’large’ persistence time(denoted τ large ) of 1300 MCS, and one has a variable’small’ persistence time (denoted τ small ). As before, weretrieve the cluster persistence time from fitting the col-lective MSD. The resulting values are plotted as a func-tion of τ small in fig. 5. It shows that the collective ben-efit of alignment can become much smaller or even non-existent ( τ p < τ large ) by adding a single cell with a smallpersistence to an existing cell cluster. We can understandthis by realising that a small persistence time correspondsto a rapid reorientation of the cell’s polarity. If the reori-entation becomes too fast, the cell will drag along othercells towards this polarity as well, which results in a de-crease of the cluster persistence. This demonstrates thatto exhibit the additional directionality of aligning clustermotion, the spread in individual persistence times shouldnot be too large. τ small (MCS) τ p ( M C S ) CPM SimulationFast - aligning ABP theory FIG. 5. Plot of the cluster persistence time τ p as a functionof the implemented small persistence time τ small . Circles de-note the results obtained from PRW fits of the MSD which inturn has been calculated from fast-aligning CPM cell clustertrajectories. Red-dotted line represents the theoretical pre-diction of our fast-aligning active Brownian particle theory,i.e. eq. (18). Simulation parameters used: N cells = 4, γ = 1, κ = 5, and τ large = 1300 MCS. Data for the MSDs has beenobtained by time-ensemble averaging over 30 trajectories eachconsisting of 50000 MCS. Active Brownian Motion
Identical Particles
To provide a more general framework for our results, wenow seek to rationalize the observed benefits of cell-cellalignment in our CPM simulations using so-called activematter theory. One of the most widely used models inthis field [23, 41, 42, 48] is the active Brownian particle(ABP) model. Such particles undergo Brownian motionwith a ’thermal’ diffusion coefficient D t , while they si-multaneously self-propel with an absolute speed v alongtheir orientation axis, called the director e i ( t ). The in-dex i here labels each of the N individual ABPs that,together, form a cluster in our theory. The evolution intime t of the position r i ( t ) = [ x i ( t ) , y i ( t )] of each parti-cle i is captured by the stochastic differential Langevinequation [23, 41–44, 48] d r i ( t ) dt = v e i ( t ) + (cid:112) D t ξ i ( t ) . (6) Here ξ i = [ ξ x i , ξ y i ] where ξ α ( α = x i , y i ) representsan independent white noise stochastic process with zeromean, (cid:104) ξ α ( t ) (cid:105) = 0, and delta-correlations (cid:104) ξ α ( t (cid:48) ) ξ β ( t ) (cid:105) = δ ( t (cid:48) − t ) δ α,β .Similar to the polarity vector of each CPM cell p σ , thedirector of each ABP is parametrised by the polar an-gle φ i ( t ) ∈ [0 , π ), i.e. e i ( t ) = [cos φ i ( t ) , sin φ i ( t )]. Inour aligning ABP model, we assume that φ i evolves intime not only according to a stochastic rotational diffu-sion process, but also due to a potential U that encodesvelocity alignment, dφ i ( t ) dt = − η ∂U∂φ i + (cid:114) τ ξ φ i ( t ) , (7)with η > τ the single particle per-sistence time, and ξ φ i another independent white noisestochastic process. For the aligning potential U we write U ( { r i } , { φ i } ) = − (cid:88) | r i − r j |
2. We haveused ( φ − φ ) = π/
9; it can be seen that the approximationis still reasonably accurate. a whole. Moreover, it allows us to simplify the sumsin eq. (9). As all involved angles are almost equal, thedirectors e i will point in roughly the same direction. Wecan therefore approximate the sum of the N directorsby N vectors which all point in the average direction ofthe particles. In other words, we can replace (cid:80) Ni =1 e i → N e cm in eq. (9) where e cm = [cos( φ cm ) , sin( φ cm )] and φ cm ≡ N (cid:80) Ni =1 φ i . A visualisation of this approximationfor N = 2 is shown in fig. 6.Additionally, since the zero mean stochastic processes ξ α are independent and delta-correlated, we can ef-fectively replace a sum of these variables by a sin-gle one via (cid:80) Ni =1 ξ α i → √ N ξ α cm . The factor √ N isadded to ensure that the correlation remains consistent,i.e. (cid:68)(cid:80) Ni =1 ξ α i ( t (cid:48) ) (cid:80) Ni =1 ξ α i ( t ) (cid:69) = N (cid:104) ξ α cm ( t (cid:48) ) ξ α cm ( t ) (cid:105) .Overall this allows us to simplify eq. (9) as d R ( t ) dt = v e cm ( t ) + (cid:114) D t N ξ cm ( t ) , (10)where ξ cm = [ ξ x cm , ξ y cm ] represents a new vector of inde-pendent stochastic processes with zero mean and delta-correlations.The time evolution of the average direction of the par-ticles (and thus of the cluster) can be formulated usingeq. (7), i.e. dφ cm ( t ) dt = 1 N (cid:114) τ N (cid:88) i =1 ξ φ i ( t ) , (11)where, due to symmetry, all alignment terms cancelagainst each other. As already mentioned, we can re-place a sum of the stochastic noise terms by a single one. Introducing a new stochastic process ξ φ cm with zero meanand delta-correlations we arrive at dφ cm ( t ) dt = (cid:114) N τ ξ φ cm ( t ) . (12)Interestingly, these resulting equations that govern themotion of the center of mass (and thus of the entirecluster of N particles) [eqs. (10) and (12)], are identi-cal in form to the equations that describe a single non-interacting ABP, i.e. eqs. (6) and (7) with η = 0. Theonly difference lies in the fact that the center-of-masspersistence time has increased in proportion to the num-ber of particles τ → N τ , while the ’thermal’ diffusioncoefficient has decreased as D t → D t /N . Summarising,we conclude that the center-of-mass motion of a clusterof N fast-aligning ABPs follows a PRW [eq. (3)] that ischaracterised by a cluster persistence time τ cm , thermaldiffusion coefficient D cmt , and active diffusion coefficient D cma given by τ cm = N τ, D cmt = D t /N, D cma = v N τ / . (13)Relating our ABPs to the CPM cells by interpreting N as N cells , we see that these theoretical results are in goodagreement with the ones from our CPM simulations (seefig. 4) and thus provide a theoretical underpinning for theincreased cluster persistence due to cell (particle) align-ment observed earlier. The prediction of the ’thermal’diffusion coefficient scaling as N − deviates slightly fromthe observed power of − . N ), we only require that all alignment contri-butions to the time evolution of φ cm will cancel out. Inother words, for all potentials that satisfy N (cid:88) i =1 ∂U∂φ i = 0 , (14)our argument and the results should be valid. Non-Identical Particles
To incorporate cluster heterogeneity, i.e. to account forthe fact that single-cell properties within a CTC clusterare generally not the same [22], we can extend our theoryanalysis to a set of N quickly aligning non-identical ABPsin several ways. One way would be to let all particlestravel at different speeds v → v ,i ; one could also makethe degree of alignment explicitly particle-dependent by0substituting η → η i and µ → µ ij (with µ ij = µ ji ) ineqs. (7) and (8) respectively. However, we refrain fromexploring these options in too much detail for severalreasons. Firstly, letting the particles travel at differentspeeds will eventually lead to particle separations thatare larger than the interaction range of the particles.This scenario is obviously incorrect for CTC clusters inwhich the cells stick together and travel at approximatelyequal velocities, and furthermore would invalidate the as-sumption of fast alignment between particles. Secondly,introducing µ ij will not change the equations governingthe center-of-mass motion: all alignment terms are stillcanceling out, i.e. eq. (14) still applies. Conversely, byintroducing a particle-dependent relaxation constant η i (e.g. to account for different cell sizes with different fric-tion constants), eq. (14) will not be valid anymore. Inparticular, the time evolution of φ cm will then also con-tain terms that are proportional to ( η i − η j ) sin( φ i − φ j ).Nonetheless, for strong enough alignment | φ i − φ j | willremain sufficiently small, allowing us to neglect theseterms and recover the same results as for identical par-ticles. In other words, a larger variety in relaxationconstants η i will result in a more narrow fast-alignmentregime for the cluster, but within this regime it does notqualitatively change its motion.The most relevant unexplored option for introducingheterogeneity in our aligning ABP model is therefore thescenario we have also studied numerically: to have eachparticle move with a different persistence time, i.e. toreplace τ → τ i in eq. (7). This implies that instead ofeq. (12) we have dφ cm ( t ) dt = (cid:16)(cid:80) Ni =1 1 τ i (cid:17) N / ξ φ cm ( t ) ≡ (cid:114) τ cm ξ φ cm ( t ) , (15)where the cluster persistence time is now given by τ cm = N (cid:32) N (cid:88) i =1 τ i (cid:33) − . (16)Note that again we have replaced the sum of stochasticterms by a single one, but due to the particle-dependent τ i this is less straightforward, i.e. (cid:80) Ni =1 (cid:113) τ i ξ α i → (cid:16)(cid:80) Ni =1 1 τ i (cid:17) / ξ α cm . When all persistence times are equal( τ i = τ ), we recover the linear increase of the persistencetime with the number of particles ( τ cm = N τ ).It is now interesting to see how the behavior of a clusterwith a distribution of persistence times compares to thecase in which all particles are identical. In fact, fromeq. (16) it can be (straightforwardly) derived (see [47]for details) that for each set of N persistence times { τ i } we have τ cm ≤ N (cid:104) τ (cid:105) (17)with (cid:104) τ (cid:105) = N (cid:80) Ni =1 τ i the average persistence time of theset and the equal sign corresponding to a constant τ i = τ . This shows that broadening the distribution of individualpersistence times of particles (keeping a constant aver-age) will always decrease the center-of-mass persistencetime τ cm . In other words, alignment is always less effec-tive in terms of increasing τ cm when individual particlestravel with a different persistence.Finally, let us verify one last numerical observation:that even a single less persistent particle among the align-ing particles can cancel out the benefits of alignment ofthe total cluster. Suppose we have N fast aligning parti-cles, N − τ large , and one particle has a smaller persistence time τ small < τ large . The persistence time of the center-of-mass motion is then given by τ cm = N τ small τ large ( N − τ small + τ large , (18)which agrees very well with our CPM simulation results(see fig. 5) and confirms that one particle can substan-tially decrease the mobility benefits of alignment. Partic-ularly, the effect of collective alignment will be canceledwhen τ cm = τ large . In that case, we find τ small = τ large N − N + 1 . (19)This value for τ small thus presents a critical value be-low which the cluster moves with less persistence thanthe individual particles. Note that for large N we have τ small ∼ τ large N → N is typically not large for CTC clusters, this effectis not negligible and the inclusion of a rapidly reorientingcell to a CTC cluster can, at least in principle, stronglysuppress its directional movement. Durotaxis
In the discussion above, we have demonstrated numer-ically, and explained theoretically, how mutual velocityalignment can increase the collective persistence of themotion of cell clusters. So far, however, we have assumedthe environment of the cell clusters (typically, the ECM)to be homogeneous, and we have included its interactionswith the cells only via constant values of the persistencetime ( τ σ ), the active speed ( κ σ ) and the adhesion coeffi-cient ( J σ, ). To establish proof-of-principle for the effectsof alignment in an inhomogeneous environment, we alsoextend the CPM simulation setup to include durotaxis(migration in a stiffness gradient); a behavior that maybe closely linked to cell persistence. Experimental resultssuggest that cells on substrates with higher stiffnessestend to exhibit greater persistence (longer persistencetimes) [19, 20] and simulations of persistently movingpoint particles have shown that a gradient in persistencetime, in itself, is sufficient to generate durotactic motion[18]. Moreover, experiments on a larger length scale have1 Time (MCS) › x fi ( a ) dτ/dx = 10 dτ/dx = 20 dτ/dx = 50 dτ/dx = 100 Time (MCS) D I x dτ/dx = 10 dτ/dx = 20 dτ/dx = 50 dτ/dx = 100 Time (MCS) D I x N cells = 1 N cells = 4 N cells = 9 − w w x τ min τ max τ A B C
FIG. 7. Plots of (A) the average displacement in the x -direction and (B,C) the x -component of the durotactic vector index as afunction of time for a single cell (A,B) and fast-aligning different-sized clusters (C) experiencing a linear gradient in persistencetime from τ min = 200 MCS ( ∼ . τ max = 2000 MCS ( ∼ . dτ /dx = 20 MCS /a . Gradients are controlled by the width w of the gradient region (see inset). Simulationparameters used: γ = 1 and κ = 5. Averages taken over 10 trajectories. indicated that in multicellular settings, entire cell mono-layers may exhibit much stronger durotaxis compared totheir isolated constituents under the same circumstances;an observation that has in part been attributed to thefact that a larger cellular collective experiences a largerstiffness differential between its leading and trailing edge[17]. This suggests that enhanced durotaxis might alsooccur in a smaller aggregate like the typical CTC cluster.Motivated by these experiments, and following the ap-proach in [18], we model durotaxis by implementing aposition-dependent single-cell persistence time τ σ = τ ( x )which, for convenience, is the same for all simulated cellsand only depends on their center-of-mass position alongthe x -axis. We then let τ ( x ) increase linearly from a min-imum value τ min to a maximum value τ max over a region x ∈ [ − w, w ] around the origin, with the cluster (or sin-gle cell) center of mass always starting in the middle ofthe gradient. Beyond the gradient region the parametersremain constant, such that τ ( x ) = τ min for x ≤ w and τ ( x ) = τ max for x ≥ w (see inset fig. 7A). This meansthat the width w effectively controls the steepness of thegradient dτ /dx . The stiffness gradient is always along thepositive x -direction, and as in most experimental setupsthe gradient only occupies part of the system [18], con-necting two regions of approximately constant stiffnessor persistence time.To first test whether our implementation of durotac-tic motion is consistent with earlier simulation work onpoint particles [18], we have studied single cell CPM sim-ulations for different gradients dτ /dx between τ min =200 MCS ( ∼ . τ max = 2000 MCS ( ∼ . x -components of theaverage cell (cluster) displacement (cid:104) x ( t ) (cid:105) and the duro-tactic vector index DI x ( t ) ≡ (cid:104) x ( t ) (cid:105) /v t respectively. Thelatter provides the fraction of the average drift velocityof the cell (cluster) along the gradient relative to its ab-solute speed v and allows us to quantify the drift upthe stiffness gradient [18, 51]. Note that y -componentsare not reported, since there is no gradient along this axis and thus no drift. The results are consistent withthe fact that a gradient in persistence time suffices toproduce a single cell flux towards the stiff region of thedomain (positive (cid:104) x ( t ) (cid:105) ), leading to a form of durotaxis.We also observe an increase (on the investigated timescale) of the drift velocity for increasing gradients (largervalues of DI x ( t )). Additionally, note that DI x ( t ) peaksand afterwards seems to decrease in the long time limit,which is a result of cells leaving the gradient region. Amapping of the retrieved results to the ones obtained forpoint particles presented in [18] shows that they are alsoquantitatively the same, making our work fully consistentwith literature and extending the point-particle results tocells with a finite area in the CPM.Having confirmed our implementation of single celldurotaxis, we now proceed by placing different-sizedaligned clusters in a fixed gradient dτ /dx = 20 MCS /a between the extremes τ min and τ max (that is, identicalenvironments but different cluster sizes). As one can seein fig. 7C, larger clusters indeed show stronger durotaxis,in the sense that the maximal durotactic index is largerfor clusters consisting of more cells. We also see thatDI x ( t ) takes longer to peak for larger clusters. Compar-ing cluster behavior to that of single cells, the enhanceddurotaxis can be attributed to the enhanced persistenceof clusters in combination with the fact that the cluster—simply because it is larger—spans a wider gradient regionand thus experiences a larger persistence differential be-tween its front and rear end. As a result, durotaxis of astrongly aligning cluster may also effectively be treatedas that of a single, fixed-size particle that moves in anincreasingly steep gradient as N increases.A prediction that follows from this observation is thatincreasing the distance between the leading and trailingedge of our model cell clusters, keeping the persistencegradient the same, will generally enhance collective duro-taxis. This is indeed what is seen in the experiments re-ported in [17], and suggests that even in the absence oflong-range force transmission any cluster (with sufficient2cell-cell adhesion to maintain cluster integrity) will showenhanced durotactic efficiency as they navigate stiffnessgradients. Conclusion
We have studied and compared the motility of single cellsand small cell clusters in the context of tumor cell clus-ters using a combination of cellular Potts modeling andanalytical active-matter theory. CTC clusters have beensuggested to pose a far more serious threat than singlecells in terms of metastatic potential, and while metas-tasis is far more involved than motility alone, differencesbetween the ways that small clusters and single cells nav-igate their environment might play an important part inthis striking difference. Our primary aim has been togain more insight into how the motile behavior of smallclusters is affected by their size, and in particular wehave focused on whether the effect of cell-cell alignmentprovides a possible mechanism for enhanced directionalmotion of clusters.We have first carried out CPM simulations to studycell motion in a homogeneous environment. Our sin-gle cell simulations show excellent agreement (evidencedby their mean squared displacement) with the persistentrandom walk, and with experimental results. Moreover,the emergent cell speed and manifested persistence timeare directly controlled by model parameters. Extendingthe simulations to small cell clusters, we have examinedthe effect of cell-cell alignment by adding a Vicsek-liketerm to the CPM. Our results demonstrate that align-ment enhances the persistence time of cell clusters; theenhancement scales linearly with the number of cells.This allows the cluster to cover more distance than asingle cell, which may play some part in its potential toinvade the extracellular matrix in the early stages of themetastatic cascade. Within our CPM description, how-ever, this advantage can be suppressed, partly or com-pletely, by adding only one rapidly reorienting cell to thecluster. In addition, we have found that reducing thestrength of alignment beyond a critical point results inrapid disintegration of the cluster.To explain the CPM results, we have proposed a theo-retical model in which CTCs are represented by a clusterof fast-aligning active Brownian particles (ABPs). Ouranalysis reveals that fast velocity alignment increases thepersistence time of the ABP cluster, yielding, consistentwith the CPM simulations, a linear scaling with the num-ber of particles. The added effect of alignment on theoverall cluster mobility is strongest for identical particles,and is also counteracted by adding a rapidly reorientingparticle to the cluster.As a first attempt to investigate the consequences ofcell-cell alignment in a more biologically relevant, inho-mogeneous environment, we have investigated durotaxis,i.e. the migration up a stiffness gradient. Within ourCPM simulations we have implemented such a stiffness gradient as a linear gradient in the persistence time. Inthis scenario, we have shown that in a fixed gradientthere indeed exists a durotactic benefit for larger clus-ters, which may be attributed to the overall persistencedifferential between the leading and the trailing edge ofthe cluster and the enhanced cluster persistence due tocell-cell alignment.Overall, we have shown that, in the presence of genericvelocity alignment, single-cell and cluster migration canbe significantly different. In particular, enhanced di-rectional migration is exhibited by larger clusters whenalignment is fast compared to a typical persistence timein the system. Since these persistence times for livingcells are generally on the order of several hours, the con-dition of rapid alignment may be quite broadly met.Our results offer specific predictions for the scaling ofboth the persistence time, as well as the random motion,of clusters of cells as a function of cluster size. These re-sults fill in a previously uncharted regime between single-cell behavior and large-scale collective motility in conflu-ent cell sheets, a physiologically very relevant regime forwhich our predictions should be directly observable inexperiments.
Acknowledgments
We acknowledge the Netherlands Organisation for Sci-entific Research (NWO) for financial support through aSTART-UP grant (V.E.D.).
A Aligned Cell Cluster Motion in 3D
To extend the CPM to a 3D system of fast-aligningcells we require new updating rules for the polarityvector p σ , which is now described by the sphericalangles φ σ ( t ) ∈ [0 , π ) and θ σ ( t ) ∈ [0 , π ): p σ =[sin( θ σ ) cos( φ σ ) , sin( θ σ ) cos( φ σ ) , cos( θ σ )]. Discretisingthe angular Langevin Dynamics for a 3D active Brow-nian particle we obtain [44] θ σ ( t + ∆ t ) = θ σ ( t ) + (cid:114) τ σ Γ(∆ t ) + ∆ t τ tan( θ σ ) ,φ σ ( t + ∆ t ) = φ σ ( t ) + (cid:114) τ σ θ σ ) Γ(∆ t ) , (A1)3 Time (MCS) -4 -3 -2 -1 D i n = M S D / t( a / M C S ) N cells = 1 N cells = 2 N cells = 4 N cells = 8 N cells -4 -2 D ( a / M C S ) , τ p ( M C S ) D a Simulation D t Simulation τ p SimulationTheory / Fit
A B
FIG. 8. (A) Plots of the instant diffusion coefficient D in = MSD / t for fast-aligning ( γ = 1) 3D cell clusters consisting of avariable amount of N cells identical cells (markers). The results have been fitted with a PRW (lines). (B) Cluster persistencetime τ p , active diffusion coefficient D a , and thermal diffusion coefficient D t obtained from the fits shown in (A). The obtainedvalues have been compared to or fitted with the derived results of the 2D fast-aligning active Brownian particle theory, i.e.eq. (13). Simulation parameters used: κ = 30 and τ = 500 MCS. Data has been obtained by time-ensemble averaging over 20trajectories each consisting of 50000 MCS. which in turn, as a result of Vicsek alignment, are ex-tended to θ σ ( t + ∆ t ) = arg θ (cid:32) γ p σ ( t ) + (cid:88) σ (cid:48) p σ (cid:48) ( t ) (cid:33) + (cid:114) τ σ Γ(∆ t ) + ∆ t τ σ tan( θ σ ) ,φ σ ( t + ∆ t ) = arg φ (cid:32) γ p σ ( t ) + (cid:88) σ (cid:48) p σ (cid:48) ( t ) (cid:33) + (cid:114) τ σ θ σ ) Γ(∆ t ) . (A2)Here arg θ ( a ) and arg φ ( a ) denote the spherical coordi-nates θ and φ of a vector a respectively.Using these updating rules we have calculated theMSDs for 3D fast-aligning ( γ = 1) clusters consistingof different numbers of N cells identical cells. Realisingthat the MSD of 2D and 3D ABPs are identical up to achange 4 D t → D t , we have again fitted the results to aPRW. Plots for D in = MSD / t including these fits andthe respective fit parameters ( τ p , D a , D t ) are shown infig. 8. Comparing with the 2D results (fig. 4) we see al-most the same behavior, i.e. D a increasing linearly with N cells , D t decreasing (almost) linearly with N cells , and τ p = N cells τ . This shows that the effect of our Vicsekalignment on the dynamics is qualitatively the same for2D and 3D CPM cell clusters. B Circular Shape Constraint
As shown in the Hamiltonian H , i.e. eq. (1), we modelcell-cell attachment via the adhesion coefficient J σ,σ (cid:48) . In particular, by setting the adhesion coefficient betweencells J cell − cell to a sufficiently small value relative to theone between cells and the medium J cell − medium , it be-comes energetically more favorable for cells to form asurface with other cells instead of with the medium. How-ever, when this difference becomes too large or the cell-cell adhesion becomes negative, the cells will be able toeasily create interfacial area with the other cells, whichcan lead to a disintegration of the cell shape. This isclearly unphysical behavior. To prevent it from happen-ing, we impose a shape constraint on the cells that forcesthe cells to have a circular or spherical shape. We caninterpret the constraint as a bending rigidity of the cellsand formulate it in the form of an energy bias given by[52]∆ H r = λ r (cid:0) r σ ( x i ) − (cid:12)(cid:12) x i − R σ ( x i ) (cid:12)(cid:12)(cid:1) (1 − δ σ ( x i ) , ) − λ r (cid:0) r σ ( x j ) − (cid:12)(cid:12) x i − R σ ( x j ) (cid:12)(cid:12)(cid:1) (1 − δ σ ( x j ) , ) . (B1)Here λ r denotes the relative strength of the constraintand r σ is the preferred radius of cell σ so that its areaor volume fits precisely in a circle or sphere respectively.The scalar | x i − R σ | denotes the length of the vector thatpoints from the center of mass of cell σ , i.e. R σ , to thelocation of the candidate site x i and can be seen as alocal cell radius. Note that the function only applies tocells ( σ > σ ( x i ) by the value of its randomly chosenneighboring site σ ( x j ). This means that the candidatecell locally retracts at its location x i , while the neighborcell locally extends towards x i . The bias checks whetheror not the extension or retraction moves the local cellradius ( | x i − R σ | ) towards or from the preferred radiusof the cell r σ . It then gives a negative energy bias formoves towards the preferred radius, thus making them4more favorable. The strength of the energy bias scaleswith the difference between the local and preferred cell radius; that is, when this difference is large, the cell ismore deformed and is therefore more likely to move to-wards the preferred radius. [1] F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A.Torre, and A. Jemal, CA: A Cancer Journal for Clini-cians , 394 (2018).[2] S. H. Au, B. D. Storey, J. C. Moore, Q. Tang, Y.-L. Chen,S. Javaid, A. F. Sarioglu, R. Sullivan, M. W. Madden,R. O’Keefe, D. A. Haber, S. Maheswaran, D. M. Lan-genau, S. L. Stott, and M. Toner, Proceedings of theNational Academy of Sciences , 4947 (2016).[3] Y. Hong, F. Fang, and Q. Zhang, International Journalof Oncology , 2206 (2016).[4] J. Massagu´e and A. C. Obenauf, Nature , 298 (2016).[5] R. L. Siegel, K. D. Miller, and A. Jemal, CA: A CancerJournal for Clinicians , 5 (2015).[6] K. J. Cheung and A. J. Ewald, Science , 167 (2016).[7] M. Giuliano, A. Shaikh, H. C. Lo, G. Arpino,S. De Placido, X. H. Zhang, M. Cristofanilli, R. Schiff,and M. V. Trivedi, Cancer Research , 845 (2018).[8] N. Aceto, A. Bardia, D. Miyamoto, M. Donaldson,B. Wittner, J. Spencer, M. Yu, A. Pely, A. Engstrom,H. Zhu, B. Brannigan, R. Kapur, S. Stott, T. Shioda,S. Ramaswamy, D. Ting, C. Lin, M. Toner, D. Haber,and S. Maheswaran, Cell , 1110 (2014).[9] R. Maddipati and B. Z. Stanger, Cancer Discovery ,1086 (2015).[10] K. J. Cheung, V. Padmanaban, V. Silvestri, K. Schipper,J. D. Cohen, A. N. Fairchild, M. A. Gorin, J. E. Verdone,K. J. Pienta, J. S. Bader, and A. J. Ewald, Proceedingsof the National Academy of Sciences , E854 (2016).[11] V. Padmanaban, I. Krol, Y. Suhail, B. M. Szczerba,N. Aceto, J. S. Bader, and A. J. Ewald, Nature ,439 (2019).[12] D. Cai, W. Dai, M. Prasad, J. Luo, N. S. Gov, andD. J. Montell, Proceedings of the National Academy ofSciences , E2134 (2016).[13] J. Kolega, Journal of Cell Science , 15 (1981).[14] B. A. Camley, Journal of Physics: Condensed Matter ,223001 (2018).[15] M. L. Lalli and A. R. Asthagiri, Cellular and MolecularBioengineering , 247 (2015).[16] M. R. King, K. G. Phillips, A. Mitrugno, T.-R. Lee,A. M. E. de Guillebon, S. Chandrasekaran, M. J.McGuire, R. T. Carr, S. M. Baker-Groberg, R. A. Rigg,A. Kolatkar, M. Luttgen, K. Bethel, P. Kuhn, P. Decuzzi,and O. J. T. McCarty, American Journal of Physiology-Cell Physiology , C792 (2015).[17] R. Sunyer, V. Conte, J. Escribano, A. Elosegui-Artola,A. Labernadie, L. Valon, D. Navajas, J. M. Garc´ıa-Aznar, J. J. Mu˜noz, P. Roca-Cusachs, and X. Trepat,Science , 1157 (2016).[18] E. A. Novikova, M. Raab, D. E. Discher, and C. Storm,Physical Review Letters , 078103 (2017).[19] M. Raab, J. Swift, P. D. P. Dingal, P. Shah, J.-W. Shin,and D. E. Discher, Journal of Cell Biology , 669(2012).[20] D. Missirlis and J. P. Spatz, Biomacromolecules , 195(2014). [21] L. Li, S. F. Nørrelykke, and E. C. Cox, PLOS ONE ,e2093 (2008).[22] P.-H. Wu, A. Giri, S. X. Sun, and D. Wirtz, Proceedingsof the National Academy of Sciences , 3949 (2014).[23] D. Campos, V. M´endez, and I. Llopis, Journal of Theo-retical Biology , 526 (2010).[24] B. A. Camley and W.-J. Rappel, Physical Review E ,062705 (2014).[25] A. Szab´o, R. Unnep, E. M´ehes, W. O. Twal, W. S.Argraves, Y. Cao, and A. Czir´ok, Physical Biology ,046007 (2010).[26] E. L˚ang, A. Po(cid:32)le´c, A. L˚ang, M. Valk, P. Blicher, A. D.Rowe, K. A. Tønseth, C. J. Jackson, T. P. Utheim,L. M. C. Janssen, J. Eriksson, and S. O. Bøe, NatureCommunications , 3665 (2018).[27] J. Zimmermann, B. A. Camley, W.-J. Rap-pel, and H. Levine, Proceedings of the Na-tional Academy of Sciences , 1226 (1995).[29] J. A. Glazier and F. Graner, Physical Review E , 2128(1993).[30] F. Graner and J. A. Glazier, Physical Review Letters ,2013 (1992).[31] C. Giverso, M. Scianna, L. Preziosi, N. Lo Buono, andA. Funaro, Mathematical Modelling of Natural Phenom-ena , 203 (2010).[32] R. M. H. Merks, S. V. Brodsky, M. S. Goligorksy, S. A.Newman, and J. A. Glazier, Developmental Biology ,44 (2006).[33] N. B. Ouchi, J. A. Glazier, J.-P. Rieu, A. Upadhyaya,and Y. Sawada, Physica A: Statistical Mechanics and itsApplications , 451 (2003).[34] B. M. Rubenstein and L. J. Kaufman, Biophysical Jour-nal , 5661 (2008).[35] R. Allena, M. Scianna, and L. Preziosi, MathematicalBiosciences , 57 (2016).[36] R. B. Potts, Mathematical Proceedings of the CambridgePhilosophical Society , 106109 (1952).[37] M. Durand and E. Guesnet, Computer Physics Commu-nications , 54 (2016).[38] A. Szab´o and R. M. H. Merks, Frontiers in Oncology ,87 (2013).[39] N. Guisoni, K. I. Mazzitello, and L. Diambra, Frontiersin Physics , 61 (2018).[40] A. Szab´o, K. Varga, T. Garay, B. Heged˝us, andA. Czir´ok, Physical Biology , 016010 (2012).[41] M. C. Marchetti, Y. Fily, S. Henkes, A. Patch, and D. Yl-lanes, Current Opinion in Colloid & Interface Science ,34 (2016).[42] C. Bechinger, R. Di Leonardo, H. L¨owen, C. Reichhardt,G. Volpe, and G. Volpe, Reviews of Modern Physics ,045006 (2016).[43] G. Szamel, Journal of Chemical Physics , 124901(2019). [44] R. Großmann, F. Peruani, and M. B¨ar, The EuropeanPhysical Journal Special Topics , 1377 (2015).[45] E. Mhes and T. Vicsek, Integrative Biology , 831 (2014).[46] K. W. Marco Scianna, Luigi Preziosi, Mathematical Bio-sciences and Engineering , 235 (2013).[47] V. E. Debets, Collective cell dynamics in cancer metasta-sis , Master’s thesis, Eindhoven University of Technology(2019).[48] A. Zttl and H. Stark, Journal of Physics: CondensedMatter , 253001 (2016). [49] F. Peruani, A. Deutsch, and M. B¨ar, The EuropeanPhysical Journal Special Topics , 111 (2008).[50] F. Peruani, F. Ginelli, M. Br, and H. Chat´e, Journal ofPhysics: Conference Series , 012014 (2011).[51] H. G. Othmer, S. R. Dunbar, and W. Alt, Journal ofMathematical Biology , 263 (1988).[52] Y. E. S´anchez-Corrales, M. Hartley, J. van Rooij,A. F. Mar´ee, and V. A. Grieneisen, Development145