Three-dimensional chemotaxis-driven aggregation of tumor cells

One of the most important steps in tumor progression involves the transformation from a differentiated epithelial phenotype to an aggressive, highly motile phenotype, where tumor cells invade neighboring tissues. Invasion can occur either by isolated mesenchymal cells or by aggregates that migrate collectively and do not lose completely the epithelial phenotype. Here, we show that, in a three-dimensional cancer cell culture, collective migration of cells eventually leads to aggregation in large clusters. We present quantitative measurements of cluster velocity, coalescence rates, and proliferation rates. These results cannot be explained in terms of random aggregation. Instead, a model of chemotaxis-driven aggregation – mediated by a diffusible attractant – is able to capture several quantitative aspects of our results. Experimental assays of chemotaxis towards culture conditioned media confirm this hypothesis. Theoretical and numerical results further suggest an important role for chemotactic-driven aggregation in spreading and survival of tumor cells.


Synopsis
Supplementary material presents more detailed information on the results discussed in the Main Text, following the plan described hereafter. In Sec. 1, we derive the chemotactic response that we have considered, by studying the microscopic ligand-receptor dynamics in presence of a chemoattractant gradient. In 2, we define the chemottractant dynamics and the cluster chemotactic response. We focus on the chemoattractant statistics and on the distribution of cluster velocities. We show that the chemoattractant concentration results from the contribution of all the clusters within the interaction distance and is uniform in space. Thus, the velocity probability distribution depends only on the statistics of the chemoattractant gradient. In particular, we show that the distribution of velocities has heavy tails due to the interaction with the nearest cluster.
In 3, we propose an aggregation-proliferation analytic model based on macroscopic quantities such as the density of clusters of given number of cells. By using the chemotactic response previously described, we solve the analytic model for relevant observables such as cluster and cell number density. Sec. 4 contains a theoretical argument to illustrate why CDA might represent a selective advantage over non-aggregating cells. Furthermore, we compare the results of the analytic model with a computational simulation that is described in 5.
In Sec. 6 we present further experimental evidence supporting the results of the main text and describe the corresponding Materials and Methods.

Dynamics of the chemoattractant receptor
The interaction between ligand and receptor is the first step in the signaling cascade that leads to cell response. The response to chemical stimuli is assumed to be proportional to the number of actively signaling receptors on the cell (see [1] and [2]).
In our model a gradient of chemoattractant, i.e. the ligand, modifies the distribution of bound receptors on cell surface. As a result of this asymmetry, the cell moves in the direction of the gradient.We approximate the cell with a sphere of radius a characterized by a surface density of receptors R s which Figure 1: Representation of the receptor dynamics. k b is the receptor-ligand binding rate, k un is the unbinding rate. A bound receptor can be internalized by the cell with an endocytosis rate k en and k ex is the number of free receptors per unit time that exocytosis brings to the surface of the cell. Adapted from [1] is the sum of a surface density R un of unbound receptors and R b of bound receptors. The surrounding cells produce a ligand that diffuses with a diffusivity D c and its concentration profile rapidly reaches the equilibrium. We consider a gradient of ligand g = |∇c|. The concentration at large distances from the other clusters is: where we have used spherical coordinates with zenith direction z and inclination θ, and approximated the gradient by means of Taylor expansion. The system is invariant for rotations along the azimuth φ. We suppose that the gradient g is small, so that ga, i.e. the variation of concentration within the cell length, is small in comparison with the concentration c 0 . We consider a first-order dynamics where k b is the receptor-ligand binding rate and k un is the unbinding rate. A bound receptor can be internalized by the cell with an endocytosis rate k en . Inside the cell a complex system of vesicles and enzymes leads to unbinding of receptors, to their recycling or degradation, possibly with the synthesis of new receptors. We simplify the description assuming that exocytosis brings k ex free receptors per unit time to the surface of the cell. The whole receptor dynamics is graphically sketched in Fig.S1. The mass-action kinetics is given by the system of differential equations where we have introduced a diffusivity D s of the receptors on the surface of the cell. Since we suppose that the gradient g varies slowly with respect to the receptor dynamics timescales, we study the equilibrium solution. The solution of the diffusion equation with the boundary condition eq. (1) is where the coefficients A 1 and A 2 are set by the conservation of ligand at cell surface In order to simplify the equations, we rescale the variables with the total number of receptors in the cell R tot . Then, we consider the variables r b = R b /R tot and r s = R s /R tot . We introduce the characteristic concentrations We solve this system of ordinary differential equations perturbatively. We study the solutions at the lowest order in ga c0 1 that is The solution in absence of gradient is isotropic and the system shows perfect adaptation, i.e. the equilibrium surface density of bound receptors is independent of the chemoattractant concentration: and the total surface density of receptors is Since in general r b < r, one has that K 2 < c 0 in order to ensure a valid solution. Solving the equations for the first-order perturbation, one has For K 2 c 0 K 3 and K d K1 K 2 c 0 the response to the gradient is Finally, assuming that the chemotactic response is proportional to the perturbation of the density of bound receptors, one has

Chemotactic model
We propose a model in which aggregation is driven by chemotaxis. We suppose that each cell in a cluster produces a diffusible chemoattractant that generates a concentration profile around the cluster itself. Cells measure the chemoattractant concentration through receptors on their surface and respond by moving toward the source, along the concentration gradient.
We suppose that each cell in a cluster responds in the same way to chemical stimuli, so that the cluster moves as a whole under the action of the chemoattractant gradient. In general the chemotactic response depends on ambient concentration. Here we assume the following form for the velocity where χ 0 is a reference chemotactic responsivity and K on , K off are reference concentrations. This expression have been derived in 1.
In this case the response is proportional to the relative variation with respect to the background value.

The chemoattractant concentration field
We suppose that each cell produces independently the chemoattractant at a rate β, so that the production of each cluster is proportional to its number of cells. Moreover, since cell motility is relatively slow with respect to molecular diffusion, the concentration profile rapidly reaches equilibrium. The stationary diffusion equation then gives where c is the chemoattractant concentration, r i are the positions of the sources, β is the cell chemoattractant production rate, D c is the diffusion constant and μ is the chemoattractant degradation rate. We remark that in this context we use the variable m i to indicate the number of cells in the i-th cluster. Since the clusters are approximately spheroidal and their size is negligible with respect to their distance, the chemoattractant production is considered as concentrated in their center of mass. By solving equation (21) with the Green's function method, we obtain the concentration profile outside the clusters where λ is the interaction length-scale of chemotaxis. For example, assuming a diffusivity of the chemoattractant D c ≈ 10 μm 2 /s and a degradation time μ −1 of 10 hours, the interaction length is λ ≈ 600μm.

Distribution of the chemoattractant concentration
The spatial distribution of chemoattractant is determined by the position of clusters. We consider an uniform distribution of clusters in a spherical volume V of radius R where τ (m) is the distribution of cluster size, i.e. number of cell. We define the random variable g of the concentration field generated by a single cluster in the center of the volume Then, by summing up the contributions of all the uniformly distributed clusters: The mean value of the single contribution is where m is the average cluster size. The second moment of the variable g is Therefore c is the sum of N independent identically distributed variables with finite mean value and variance and the central limit theorem can be applied. As a result c is distributed as a Gaussian with mean value c = N w and variance σ 2 c = N ( w 2 − w 2 ). Finally, one can take the infinite volume limit R → +∞, while keeping constant the density of clusters n = N V , to obtain with M = n m density of cells and M 2 = n m 2 . The coefficient of variation of the concentration is Over the experimental ranges of densities of clusters, considering m 2 / m 2 10, the standard deviation is always much smaller than the mean value. The average concentration is therefore approximately constant in space for all clusters configurations and equal to eq. (29). It is worth noting that here we assumed each cluster to be filled with cells, while in many contexts cancer cell spheroids are found to be hollow either because of some residual glandular morphology reminiscent of the epithelial origin or because of necrosis. However, the assumption that clusters are hollow would simply rescale all terms m into 3 m 2/3 , without affecting the dependence of the chemoattractant concentration on the cluster density. The velocity of a cluster depends on the chemoattractant concentration that it senses, including the contribution by the cluster itself. This induces an additional term that can be evaluated as Clearly, as the density of cells increases the self-interaction becomes negligible and the total concentration of chemoattractant is simply proportional to the density of clusters. We conclude that the concentration level is mostly determined by the effect of distant clusters, as they are the majority. In the next section we will see that the concentration gradient is instead largely governed by the interaction with neighboring clusters.

Distribution of the chemoattractant gradient
We will now consider the distribution of the gradient From the isotropy of the distribution of clusters it follows that and we can consider, without loss of generality, the component u = ∇ x c of the gradient, that is the sum of the contributions w due to each cluster The random variable w has infinite variance because of the divergence in r = 0 with w 2 ∼ r −4 and the central limit theorem cannot be applied. In this case it is necessary to resort to the generalized central limit theorem [3]. The distribution of u, sum of a large number of independent w variables, tends to a stable distribution, determined by the tails of p(w). A symmetric stable distribution p(w) has the asymptotic behavior and it is defined by the index α and the scale factor w * . The tail of p(w) is determined by the behavior of w Gm r 2 for r → 0 The velocity u thus follows a stable distribution of index α = 3 2 , known as the Holtsmark distribution. The width parameter u * grows with the number of clusters according to Comparing (36) and (37), from (38) one obtains where n is the density of clusters.

Large gradients are determined by nearest-neighbor interaction
In the previous section we have calculated the exact distribution of the gradient for an uniform distribution of clusters. We show here that the tails of the distribution p(u) are determined by the interaction with the closest cluster. The distribution of the position of the nearest neighbor is [4] p nn (r 1 , cos θ 1 ) = 2πnr 2 1 e − 4 3 πr 3 1 n From this distribution we calculate that the average distance of the nearest neighbor is For a density of clusters of 100/mm 3 , their average distance r 1 ≈ 100 μm is smaller than the interaction length λ ≈ 600μm calculated above. Then we can approximate the velocity as: Under this approximation, we obtain for the distribution of gradient induced by the nearest neighbor p nn (u) = 1 2Gm where γ(s, x) is the lower incomplete Gamma The asymptotic behavior of the distribution is that coincides exactly with the asymptotic behavior of p(u). This agreement means that large values of u are entirely due to the interaction with the nearest neighbors. As shown in the Main Text, this approximation ceases to be valid for intermediate and low gradients.

The contribution of farther clusters
Since the nearest neighbor approximation captures the large-velocity events only, it is necessary to take into account the other contributions in order to improve the approximation. The gradient can be split in two terms, the first one due to the nearest neighbor, the second one to the other clusters Gm cos θ k r k The clusters contributing toũ are uniformly distributed outside a sphere of radius r 1 , the position of the nearest neighbor. For the single contributioñ the mean value w k = 0 and the variance w k 2 is finite. Thus, the central limit theorem can be applied. The total contribution of the outer clusters is distributed as a Gaussian with zero mean value and variance The distribution of u is the convolution of the distributions of the two contributions p(u) = where the integrations overũ and cos θ 1 have been performed. The fact that the two terms contribute to different regions of p(u) suggests to approximate them as if they were independent.
The velocity u thus can be approximated as the the sum of the nearest neighbor interaction plus a Gaussian contribution. The variance of this second contribution is given by (49), where r 1 is replaced by the average position of the nearest neighbor. As shown in the Main Text this approximation gives a very good agreement over the whole range of gradients.

Distribution of velocities
On the basis of the previous calculations we can now consider the velocity of the clusters As discussed above, the distribution of c is strongly peaked and the concentration can be taken uniform and not fluctuating. Since ∇c results the only fluctuating variable, we can consider v = χ( c )∇c (55) where χ( c ) is the chemotactic coefficient where c and u * are respectively given by (29) and (39). We remark that since c ∝ n and u * ∝ n 2/3 , in the range K on c K off the characteristic velocity v * is proportional to n −1/3 .

Cluster proliferation-aggregation model
A complete description of cluster aggregation should include the spatial and geometrical properties of the clusters, such as their position, shape, size and velocity. Here we focus initially on a mean field model whose relevant variables are the number n i of clusters per unit volume of a given number of cells i.
The master equation describing the evolution of the system is a generalization of the Smoluchowski aggregation equation [5] The aggregation process is determined by the kernel K i,j , i.e. the average rate at which two clusters of i and j cells coalesce. We have considered a constant cell replication rate α.
The most important quantities for our purposes are the total density of clusters N = i n i and the total density of cells M = i n i i evolving as

Chemotactic kernel
for short distances with respect to the interaction length λ. This assumption is justified by the fact that the interactions leading to collision take place between neighboring clusters. We choose a reference system in the center of the cluster of size i and we impose at large distance a density n j (∞) of clusters of size j. The equation for the the spatial distribution of the clusters of size j is where v = v i + v j is the relative velocity of the clusters. Since we want to calculate the average collision rate, we study the stationary solution. By integrating the stationary equation on a spherical shell around the cluster i, one obtains that the flux of clusters across each spherical surface is conserved and given by Therefore, we conclude that the aggregation rate of two clusters of number of cells i and j is

Chemotactic aggregation
The chemotactic kernel (64) is now inserted into the master equation (58) to study the aggregative dynamics. In the range of concentration K on c K off it takes the form where eq. (29) has been used. The equation for the density of clusters N (59) becomes with solution: where N 0 is the initial density of clusters. Then, for this kernel, the aggregation time is independent of the initial density of clusters. The equation for M (t) has already been obtained in eq. (60). In the case of a constant replication rate α it has the solution It is important to stress that, thanks to the independence of the density of clusters from cell density, the aggregative dynamics does not depend on the replication rate.

Advantage of aggregating vs non-aggregating cells.
An interesting question is whether aggregation might represent an advantage over cells simply producing chemoattractant but without coalescing. This question can be answered by the following calculations. The flux of molecules secreted by a spherical cluster of radius R is where D c is the diffusion coefficient. This expression follows from the solution of the diffusion equation D c ∇ 2 c = 0 outside the sphere, with vanishing concentration at infinity. Now assume that at the surface of the cluster there will be production of chemoattractant at a rate Π + , and consumption at a rate Π − . These two contributions read: where n p and n a are the surface densities of cells that are producing or absorbing the ligand, respectively. The other parameters are k b , the absorption rate, and β, the production rate of ligand, per cell.
At the stationary state, the flux equates the balance between production and consumption As a result, the consumption rate of ligand at the surface per cell is given by which grows with the cluster size R and saturates for R R * = D c /(k b n a ) to the value βn p /n a . Therefore, as long as the cluster's size is smaller or comparable to R * the absorption of ligand per cell increases with the cluster's size.
The physical reason for such an enhancement resides in the properties of diffusion. Secreted molecules tend to bounce back to the surface with great probability: this results in the well known linear dependence of φ on the size R as opposed to the quadratic dependence for ballistic motion. Conversely, production and consumption do grow proportionally to the surface, i.e. as R 2 , in this case. Increasing R mitigates the effect of ligand dispersion relative to production and thereby increases the efficiency of recapture. Indeed, as R R * the ratio of consumption to production approaches unity.
This argument can be generalized to the case when production and consumption takes place also inside the volume. The calculations are more cumbersome and will not be reproduced here, but, as the intuition suggest, the argument above holds as well and the result is qualitatively similar.
Finally, it is only left to estimate R * in order to check that is sufficiently large with respect to the cell size. For this purpose we use the following values: where for n −1 a we have used an estimation of about half the surface of a cell with a 10μm diameter, while D c has the same estimation as in the paper and for k b we refer to the values cited in ref. [1]. With these values, we obtain R * ∼ 10 6 which is four orders of magnitudes larger than the radius of a single cell.

Cluster proliferation-aggregation computational model
We propose a computational model in which clusters are treated as spheres of volume proportional to the number of cells. Initially clusters are uniformly distributed in a cubic box of side L. Clusters move as rigid bodies with a velocity v = χ 0 K on ∇c c where the concentration is measured in the center of mass of the clusters and it comprises the self-interaction term. At any time, the concentration measured by a cluster is as given by eq. (22). The r j are the positions of the clusters and a selfinteraction term is given by eq. (32). Coalescence occurs when two clusters overlap. In this case the two clusters coagulate in a single one, located in the center of mass of the pair. We have used a constant replication rate. The parameters of the models are defined in Table 1 in the main text.

Periodic boundary conditions
Chemotactic aggregation is a typical many-body problem with long-range interactions. The computational cost for a single time-step scales as O(N 2 ) where N = λ 3 ∼ 10 4 . These figures are clearly prohibitive, so that we opted for the alternative strategy of partitioning the interaction volume λ 3 in smaller boxes of volume L 3 and use periodic boundary conditions. The concentration field produced by the clusters is the solution of eq. (21). The periodic Green function is a sum over all the lattice sites r s where R = r − r . The concentration is obtained integrating over the source distribution ρ(r ) and in the case of point sources it is with r i position of the sources. As the series in eq. (79) is slowly convergent, we resorted to the Ewald method, originally employed in solid state physics [6]. This method relies on the transformations of the series in eq. (79) in an alternative form of extremely fast convergence, by dividing appropriately the summation in coordinate and reciprocal spaces: where is a summation over reciprocal lattice vectors k n and is a summation in real space. The parameter η defines the way the sum is divided in the two terms and it has to be optimized to obtain the fastest convergent series. The calculation of the integral of equation (83) leads to

Supplementary experimental procedures
Migration assays presented in the supplementary material have been performed by means of the xCELLigence RTCA DP instrument (ACEA, San Diego, CA, USA, details can be found in ref. [7] and refs therein) which was placed in a humidified incubator at 37 • C and 5% CO2. Briefly, this instrument exploits a well-established technique based on the measurement of the frequency dependent electrical impedance of cell-covered electrodes subject to a small alternate electric current. Cells adhering on the electrodes vary the impedance in a frequency dependent manner. Migration is measured with a porous membrane analogous to Transwells which bottom is covered with electrodes measuring the presence of cells that migrated through the membrane (CIM-plate, ACEA, San Diego, CA, USA). 40000 cells were seeded for this experiment. Bars in the plot represent final points at 3h after seeding.

Supplementary Movies legends SI Movie 1
PC3 cells were grown overnight in U-bottom multiwell plates and the culture media was supplemented with 1% methylcellulose, washed twice with PBS, resuspended in Matrigel and used in aggregation assay. The result of this procedure is to generate a broad distribution of cluster size.

SI Movie 2
Aggregation assay of PC3. Cells were detached from standard culture plates and resuspended as a single cell suspension into Matrigel at the desired concentration and put under the microscope. The result of this procedure is to generate a narrow size distribution, i.e. mostly single cells.

SI Movie 3
DU145 cells were processed as in SIMovie 2 and observed under the microscope.

SI Movie 4
LNCAP cells were processed as in SI Movie 2 and observed under the microscope.

SI Movie 5
MDA-MB-231 cells were processed as in SI Movie 2 and observed under the microscope.