Differences in mechanical properties lead to anomalous phase separation in a model cell co-culture

TL;DR

This study uses simulations to explore how differences in mechanical properties of cells influence phase separation, revealing unique clustering behaviors relevant to tissue development and disease.

Contribution

It introduces a simulation model showing how mechanical heterogeneity drives phase separation and emergent structures in cell co-cultures.

Findings

Cells with different stiffness and adhesion segregate into distinct clusters.

Cluster growth and structure deviate from classical binary fluid phase separation.

Results align with experimental observations of cell co-cultures.

Abstract

During the morphogenesis of tissues and tumors, cells often interact with neighbors with different mechanical properties, but the understanding of its role is lacking. We use active Brownian Dynamics simulations to study a model co-culture consisting of two types of cells with the same size and self-propulsion speed, but different mechanical stiffness and cell-cell adhesion. As time evolves, the system phase separates out into clusters with distinct morphologies and transport properties for the two cell types. The density structure functions and the growth of cell clusters deviate from behavior characteristic of the phase separation in binary fluids. Our results capture emergent structure and motility previously observed in co-culture experiments and provide mechanistic insights into intercellular phase separation during development and disease.

Tables

Table 1. Table 1: Simulation parameters

$v_0$	1 $\mu$m min-1Smeets et al. (2016)	$D_r$	0.05 min-1Smeets et al. (2016)	${\mu }^{-1}$	8 nN$\mu$m min-1
$Y_h$	1.5 nN$\mu$m-2Lekka et al. (2012)	$Y_c$	0.4 $Y_h$Lee et al. (2012)	$\nu$	0.5
${\sigma }_{hh}$	0.5 nN$\mu$m-1Maître and Heisenberg (2013)	${\sigma }_{cc}$	0.3 ${\sigma }_{hh}$	${\sigma }_{hc}$	0.1 ${\sigma }_{hh}$
$d$	20 $\mu$m	$L$	256	$\varphi$	1.2

Equations

$$\frac{d{\vec{r}}_i}{dt}=v_0{\hat{n}}_{\theta ,i}+\mu \sum _{j\in n_i}\vec{F}(\overrightarrow{r_i},\overrightarrow{r_j}),\frac{d{\theta }_i}{dt}=\sqrt{2D_r}{\eta }_{{}_i},$$

$$\displaystyle\!\!\vec{F}^{\rm JKR}_{ij}=\frac{1}{2}\bigg{[}Y^{\rm eff}_{ij}{d}^{1/2}\xi_{ij}^{3/2}\!-\!\sqrt{3\pi\sigma_{ij}Y^{\rm eff}_{ij}}{d}^{3/4}\xi_{ij}^{3/4}\bigg{]}\frac{\vec{r}_{ij}}{|\vec{r}_{ij}|},$$

$$S\rho \rho (k,t)$$

$$P^{\alpha }(s,t)=\frac{1}{\mathcal{A}(t{)}^2}{f}^{\alpha }(s/\mathcal{A}(t)),$$

Full text

Self-organization and phase separation of cells are critical to tissue morphogenesis, whether in the context of formation of embryos or the formation and progression of tumors. During these processes, cells are often surrounded by neighbors with different mechanical properties. For example, recent studies have shown that for many types of tumors, cancer cells are mechanically softer and more deformable than the corresponding non-cancerous cells Suresh (2007); Lekka et al. (2012); Lee et al. (2012). Furthermore, while healthy, epithelial cells express specific proteins such as E-cadherin that facilitate cell-to-cell adhesion, in most cancer cells and other mesenchymal cells, the expression of this protein is significantly down regulated Oka et al. (1993); Wijnhoven et al. (2000); Song et al. (2016); Gamboa Castro et al. (2016). This leads to the following questions: how do the variations in the mechanical properties of cells influence their collective dynamics and self-organization? Can some cells harness these differences for enhanced organization and migration better than others? Addressing these questions will help to understand the role of mechanical heterogeneity in cellular organization, dynamics, and functions, and may provide insights into how changes in cells’ mechanical properties during tumorigenesis influence cancer invasion and metastasis.

In recent years, there has been tremendous focus on phase separation within cells, the associated physical mechanisms, and the consequences for biological function Hyman et al. (2014); Brangwynne et al. (2009); Shin and Brangwynne (2017); Alberti et al. (2019). When it comes to phase separation of multicellular systems, while the fully segregated final states can be understood in terms of the Differential Adhesion Hypothesis (DAH) Steinberg (1962); Steinberg and Takeichi (1994); Foty and Steinberg (2005); Manning et al. (2010); Gonzalez-Rodriguez et al. (2012); Song et al. (2016) as corresponding to the minimum in the free energy of a mixture of immiscible liquids with different surface energies, the dynamics via which these states are reached is not well understood. Furthermore, the DAH focuses on the differences in cell adhesion and does not account for the differences in cell stiffness. During the past decade, there have also been many computational studies of phase separation in single species active colloidal systems Buttinoni et al. (2013); Palacci et al. (2013); Fily and Marchetti (2012); Redner et al. (2013a); Cates and Tailleur (2015) and in binary systems motivated by mixtures of cells with different motility or adhesion Belmonte et al. (2008); Méhes et al. (2012); Beatrici et al. (2017); Graner and Glazier (1992); Kabla (2012); Nakajima and Ishihara (2011); Osborne et al. (2017). In the latter, the interactions between cells were phenomenological, and an understanding of the role of differences in cell mechanical properties is still lacking. Here we ask how do differences in cell stiffness and adhesion impact their phase separation in a co-culture and whether the structure and dynamics of the phase separated system have characteristics similar to binary liquids.

To address these questions, we study a model co-culture consisting of two different types of cells that have the same size and self-propulsion speed but different Young’s moduli and surface energies. We model this system as a 2D suspension of spherical overdamped self-propelled particles with mechano-adhesive contact interactions. The equation of motion of a cell $i$ $\in\{1,N\}$ ($N$ being the total number of cell) with position vector $\vec{r}_{i}$ and the orientation angle $\theta_{i}$ of self-propulsion is given by Fily and Marchetti (2012); Redner et al. (2013a),

[TABLE]

where $\hat{n}_{\theta,i}$ and $v_{0}$ are the unit vector and the speed associated with the propulsion velocity, and $\mu$ is the mobility. The intercellular force between cells $i$ and $j$, $\vec{F}(\vec{r_{i}},\vec{r_{j}})$, is non-zero only when the two cells are in contact, and the sum is over all contacting neighbors of cell $i$. Here, we neglect the translation noise as in Szabó et al. (2006); Belmonte et al. (2008); Smeets et al. (2016). The random variable $\eta_{i}(t)$ has zero mean and is $\delta-$ correlated i.e., $\langle\eta_{i}(t)\rangle=0$ and $\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=2\,\delta_{i,j}\delta(t-t^{\prime})$ and $D_{r}$ is the angular diffusivity. Furthermore, there is no built-in alignment interaction between cells; these interactions may impact the collective dynamics of some cell types Méhes and Vicsek (2013), and will be incorporated in future work.

We model the contact force between cells using the well known Johnson Kendall Roberts (JKR) model Johnson et al. (1971). The JKR model was originally proposed for soft sticky spherical objects in contact Johnson et al. (1971), and has been experimentally shown to work well for several cell types Chu et al. (2005); Zhu et al. (2016). Given constitutive cell mechanical properties such as cell stiffness and interfacial energy of cell-cell contacts, the JKR model gives the intercellular contact force as a function of the overlap distance between cells, and is therefore better suited for studying differences in mechanical properties than commonly used phenomenological models which assume the contact force to be a piecewise linear function of the intercellular distance Szabó et al. (2006); Belmonte et al. (2008); Smeets et al. (2016).

For a pair of overlapping spherical particles $i$ and $j$, the JKR force can be formulated in terms of the virtual overlap $\xi_{ij}=d-|\vec{r}_{ij}|$, where $\vec{r}_{ij}=\vec{r}_{j}-\vec{r}_{i}$ and $d$ is the particle diameter (Fig. 1A). For computational simplicity, we use an approximate form of this force which assumes weak adhesion and is given by Brilliantov and Pöschel (2005); Beyer and Meyer-Hermann (2009),

[TABLE]

where $\sigma_{ij}$ is the interfacial energy of the contact between $i$ and $j$, and $Y^{\rm eff}_{ij}$ is their effective Young modulus, given by $\frac{1}{Y^{\rm eff}_{ij}}=\frac{3}{4}\left(\frac{1-\nu^{2}_{i}}{Y_{i}}+\frac{1-\nu^{2}_{j}}{Y_{j}}\right)$. Here, $Y_{i(j)}$ and $\nu_{i(j)}$ are the Young’s modulus and Poisson’s ratio of particle $i(j)$.

The first term in Eq. 2 corresponds to a repulsive contribution due to cell stiffness, and the second term to cell-cell adhesion. The equilibrium overlap distance between the two particles is given by the value of $\xi_{ij}$ where $F^{\rm JKR}=0$, and the force required to pull the two particles off each other is the maximum negative value of the JKR force (See Fig. 1A). Both these quantities are critical to the collective structure and dynamics of the system and depend on the elasticity of the particles and the interfacial energy of the contact. While the JKR model has been used in several other theoretical and computational studies of cell mechanics Drasdo et al. (2007); Byrne and Drasdo (2008); Brilliantov and Pöschel (2005); Beyer and Meyer-Hermann (2009), to our knowledge it has not been used to study demixing in cell co-cultures.

With these model ingredients, we consider two cell types, which we refer to as $h$ and $c$ cells, with distinct mechanical properties. The $h$ and $c$ cells are motivated by the healthy breast epithelial cell line MCF10A and the breast cancer cell line MDA-MB231 respectively, and have approximately the same Young’s moduli and ratio of surface energies as these cells (see Table 1). Therefore, $Y_{h}>Y_{c}$, and $\sigma_{hh}>\sigma_{cc}>\sigma_{hc}$. The resulting JKR force exhibits stronger adhesion at small overlap distances and consequently a larger pull-off force, as well as larger elastic repulsion at large overlap distances for healthy breast cells ($h-h$) than breast cancer cells ($c-c$) (Fig. 1A).

To simulate the dynamics of this model system, we first non-dimensionalize the dynamical equations ( Eqs. 1) using cell diameter $d$ as the characteristic length scale, and time $\tau=d/v_{0}$ needed to cross $d$ when moving with a self-propulsion speed $v_{0}$ as the time scale. We integrate these equations using the Euler method with a step size $\delta t=0.001\,\tau$. The simulations are performed in a 2D box of size $L\,d\times L\,d$ with equal numbers of h and c cells for area fraction $\phi$, using periodic boundary conditions.

Starting with a spatially homogeneous mixture, we find that as time evolves, the $h$ and $c$ cells separate out, forming clusters of individual cell types which coalesce and grow. The time evolution of a typical phase-separated pattern is shown in Fig. 1B-D. The clusters of the two cell types have distinct morphologies: while $h$ cell clusters are compact with well-defined boundaries, those of $c$ cells are dispersed with ill-defined boundaries. Similar morphologies were observed in recent experimental studies of co-cultures of MCF10A and MDA MB231 Gamboa Castro et al. (2016); Pawlizak et al. (2015).

Here we note that unlike Beatrici and Brunnet (2011) where phase separation is driven by a difference in self-propelled velocity, in our model all the cells have the same self-propulsion speed and there is no alignment interaction between them. Coming to comparisons with binary liquids, interfacial tension plays a critical role in the phase separation of our model co-culture as in these systems. However, the compact and dispersed morphologies seen here are not usual in binary liquids Bray (1994). As far as we are aware, they have not been reported in theoretical studies of phase separation before. The origin of such distinct cluster structures lies in the interplay of differences in mechanical properties of the cells. The combination of greater stiffness and greater cell-cell adhesion of $h$ cells (Fig. 1A) causes them to form clusters that are simultaneously more closely packed and have a well-defined boundary. On the other hand, $c$ cells are less adhesive and less stiff, leading to both weak attraction when they come in contact and weak repulsion when they push each other, resulting in dispersed clusters with irregular boundaries. Together with the barely attractive $h-c$ interactions, this leads to rough interfaces between the $h$ and $c$ cell domains.

To characterizes the morphologies of the clusters and growth, we measure spatio-temporal density structure factors. First, we define a local coarse-grained density field for a cell-type $\alpha$, $\rho^{\alpha}(\vec{r},t)$ as the total number of $\alpha$ ($=h$ or $c$) cells inside a square box of unit area. Then, the density structure factor is given by $S^{\alpha}_{\rho\rho}(\vec{k},t)=\langle\tilde{\rho}^{\alpha}(\vec{k},t)\tilde{\rho}^{\alpha}(-\vec{k},t)\rangle$, where $\tilde{\rho}^{\alpha}(\vec{k},t)$ is the Fourier transform of $\rho^{\alpha}(\vec{r},t)$ at a wave vector $\vec{k}$ and the angular bracket represents average over many ensembles.

For a system undergoing phase ordering kinetics, the structure factor follows the scaling hypothesis Bray (1994):

[TABLE]

where $\mathcal{L}(t)$ is the growing length scale for the density field and $d$ is the dimensionality. The scaling function $g(y)$ shows a power law decay as $g(y)\sim y^{-\theta}$ for $y\gg 1$. For a system with a sharp, regular interface, the structure factor of a scaler order parameter, such as density, has an exponent $\theta=d+1$. This is the well known Porod law Porod (2 NY). Phase separation in binary fluids ordering towards equilibrium often show morphologies obeying this law Bray (1994). Deviations from the Porod law may indicate a rough interface or structural hierarchy and have been observed in non-equilibrium systems Das and Barma (2000); Shinde et al. (2007); Mishra and Ramaswamy (2006); Dey et al. (2012). The length scale of typical $\mathcal{L}(t)$ grows as $t^{\beta}$, where $\beta=1/3$ for a system with a conserved order parameter and is attributed to growth dominated by diffusion such as in diffusive coalescence of droplets in binary liquids Bray (1994); Siggia (1979); Furukawa (1985).

We measure the structure factor for $h$ and $c$ cell clusters, separately. In Fig. 2A and B, we show the structure factor data $S^{(h)}_{\rho\rho}$ for $h$ cells and $S^{(c)}_{\rho\rho}$ for $c$ cells respectively at four different times. We found that for $h$ cells, the data at different times show an excellent collapse according the scaling hypothesis (Eq. 3) with $\mathcal{L}(t)\sim t^{\beta_{h}}$, where $\beta_{h}\simeq 0.23\pm 0.02$ (Fig. 2C). The large $k\mathcal{L}(t)$ data of scaled structure factor, show a power law decay with exponent $\theta^{h}=2.6$. This is a small but clear deviation from the Porod law which predicts an exponent $3$ in $2D$. We believe it arises because the interfaces between the $h$ and $c$ domains are rough even though the $h$ domains themselves are compact. The structure factor for $c$ cell clusters are distinct from $h$ cell clusters. For $c$ cells, we are able to obtain a good scaling collapse of the structure factor data only at large wavevectors, and only when we replace $d$ in equation (Eq. 3) by $d/2$, or in other words when the structure factor is additionally scaled with $t^{-0.23}$, an unusual scaling which was previously observed in a freely cooling granular gas Shinde et al. (2007). At large wavevectors, this scaled structure factor shows a power law decay with an exponent $\theta^{c}=1.1$, indicating considerably stronger violation of the Porod law due to the dispersed nature of the $c$ cell clusters in addition to the rough $h-c$ phase boundaries.

Also, note the deviation of the cluster growth exponent $\beta$ from the scaling $1/3$ commonly observed in droplet coalescence in liquid-liquid phase separation. This may be qualitatively understood as follows. In binary liquids, droplet merging due to Brownian diffusion gives a timescale of growth $t\sim{\mathcal{L}(t)}^{2}/D$, where the diffusivity $D\sim K_{B}T/\eta\mathcal{L}(t)$, which in turn gives $\mathcal{L}(t)\propto t^{1/3}$ Bray (1990); Siggia (1979). In the above, $K_{B}$ is the Boltzmann constant, $T$ the temperature, and $\eta$ the viscosity. In our model co-culture, the effective diffusivity underlying the coarsening of clusters scales as $1/\mathcal{L}(t)^{d}$ in the absence of alignment interactions between cells Beatrici et al. (2017); Kolb (1984), which leads to a growth law $\mathcal{L}(t)\propto t^{\frac{1}{d+2}}$ ($t^{1/4}$ in 2D) Furukawa (1985). A growth exponent $\sim 1/4$ was previously reported in some passive Furukawa (1985); Corberi et al. (1995); Wolterink et al. (2006) and active phase separating systems Dey et al. (2012); Redner et al. (2013b); Beatrici et al. (2017); Stenhammar et al. (2014).

To understand the above morphologies and growth laws better, we measure the distribution of cluster sizes $P^{(\alpha)}(s,t)$ at different times, where $s$ is the total number cells in a cluster of cell type $\alpha$. In Fig. 3, we plot the probability distribution for cluster sizes for different times for $h$ and $c$ cells. We observe a data collapse according to,

[TABLE]

where $\mathcal{A}(t)=\mathcal{L}(t)^{d}=t^{2\beta}=t^{0.46}$ which implies that cluster sizes for both $h$ and $c$ cells grow as $t^{0.23}$, consistent with the structure factor measurement. For $c$ cells, we observe that the cluster sizes are distributed as a power law over many decades (Fig. 3D), suggesting hierarchical clusters as the origin of the non-Porod behavior. We find the decay exponent for cluster size distribution is $-1.8$. The implication of this particular exponent will be explored in future work. For $h$ cells, we do not observe a clear power law decay of clusters sizes.

Finally, we turn our attention to the emergent motility of $h$ and $c$ cells. We measure the mean squared displacement (MSD) for the $h$ and $c$ cells once they reach a “steady state” at time $t_{0}$. As the phase separation dynamics is very slow, for computational efficiency, we calculate MSDs for a relatively small system size $L=64\,d$, and $t_{0}=8192$ (Fig. 4). The MSD for a single isolated cell with self-propelled speed $v_{0}$ and angular noise strength $D_{r}$ is given by, $\langle\delta r^{2}(t)\rangle=4D_{0}\left[t+1/D_{r}({\rm e}^{D_{r}t}-1\right)]$, where $D_{0}=v^{2}_{0}/(2D_{r})$ Szabó et al. (2010). The motion is ballistic at small times scales and diffusive at large timescales. We observe a similar ballistic to diffusive crossover in the model co-culture, but with a smaller effective diffusion coefficient and effective velocity than for a single isolated cell. More interestingly, the $h$ and $c$ cells follow two different MSD curves although their self-propulsion speeds are the same, with the MSD of the $c$ cells being greater than that of the $h$ cells at any given time. This suggests that the effective motility of $c$ cells in a co-culture is greater than that of the $h$ cells, consistent with experimental observation Lee et al. (2012), and can arise solely from the differences in mechanical properties of the two cell types.

In summary, we have studied the phase separation of a model co-culture of two cell types $h$ and $c$, where $c$ cells are mechanically softer and show weaker cell-cell adhesion than $h$ cells. We used previously reported experimental values of breast epithelial cells and breast cancer cells for the mechanical properties of the $h$ and $c$ cells respectively. In numerical simulations of this system where we modeled the cells as active particles interacting via mechanical contact forces, we observed distinct morphologies of the phase separated domains qualitatively similar to those observed in experiments Gamboa Castro et al. (2016); Pawlizak et al. (2015); while $h$ cells formed compact clusters, $c$ cell domains were much more dispersed. We found that the structure functions of both $h$ and $c$ cell clusters deviate from the Porod behavior. The deviation was particularly striking for the $c$ cell clusters and arise due to a hierarchy of cluster sizes. We note that the non-Porod behavior has been observed in other biological contexts such as the microstructure of goat femur Mayya et al. (2013). The $c$ cells also exhibited enhanced emergent motility aided by bursty movement consistent with previous experimental observations Lee et al. (2012); Gamboa Castro et al. (2016). The domain coarsening follows a $t^{1/4}$ growth law due to slower coalescence of clusters. Our results demonstrate how mechanical heterogeneity can lead to phase separation and differential cell motility in co-cultures of epithelial and mesenchymal cells, and may provide insights into the roughness of tumor spheroids and enhanced migration of cancer cells during tumor progression.