Far from equilibrium dynamics of tracer particles embedded in a growing multicellular spheroid

The interplay between short-range forces and non-equilibrium processes arising from cell division and apoptosis gives rise to unexpected dynamics in the collective migration of cancer cells Friedl and Gilmour (2009); Shaebani et al. (2020); Kumar and Weaver (2009); Desoize and Jardillier (2000); Walenta et al. (2000); Laurent et al. (2013); Valencia et al. (2015). An example is the invasion of cancer cells (CCs) in a growing multicellular spheroid (MCS), which is relevant in cancer metastasis Laurent et al. (2013); Valencia et al. (2015). Imaging experiments show that collective migration of a group of cells that maintain contact for a long period of time exhibits far from equilibrium characteristics Valencia et al. (2015); Richards et al. (2019); Martino et al. (2019); Han et al. (2020); Palamidessi et al. (2019). Simulations and theory explain some of the experimental observations Malmi-Kakkada et al. (2018); Samanta and Thirumalai (2019); Sinha et al. (2020). Dynamics in a growing MCS is reminiscent of the influence of active forces in abiotic systems Marchetti et al. (2013); Bechinger et al. (2016); Nandi et al. (2018). In a growing MCS the analogue of active forces are self-generated Sinha and Thirumalai (2020), arising from biological events characterized here by cell division and apoptosis.

Experiments that probe the local stresses or pressure on the CCs Rauzi and Lenne (2011); Hutson et al. (2003); Boucher et al. (1990); Fadnes et al. (1977); Campàs et al. (2014); Dolega et al. (2017) have provided insights into the mechanism by which the CCs invade the extracellular matrix. Recently, the stresses within MCSs were measured Dolega et al. (2017) by embedding micron-sized inert deformable polyacrylamide beads as probes. For a TP to be a reliable sensor, it should not alter the dynamics of the CCs and the microenvironment, as assessed by local pressure. However, the influence of TPs on the CCs is unknown – a gap that we fill here.

We used theory and simulations to probe the dynamics of the TPs in the MCS, and the impact of the TPs on the CC mobility as well as radial pressure. The central results are : (1) The TPs exhibit sub-diffusive motion in the intermediate time scale ($t\lesssim\frac{1}{k_{b}}$), with the mean squared displacement, $\Delta_{TP}(t)\sim t^{\beta_{TP}}$, with $\beta_{TP}$ less than unity. In the long time limit, ($t\gtrsim\frac{1}{k_{b}}$), the TPs undergo hyper-diffusive motion, $\Delta_{TP}(t)\sim t^{\alpha_{TP}}$ with $\alpha_{TP}$ greater than 2. In contrast, the CCs, which are jammed at $t\lesssim\frac{1}{k_{b}}$, exhibit super-diffusive dynamics. (2) The dynamics of the CCs in the intermediate time regime changes slightly as the size of the TPs are varied. Remarkably, the long time exponent ($\alpha_{CC}$) as well as the decrease in the pressure as a function of the distance from the center of the MCS remains unaffected in the presence of TPs.

Theory:  We consider the dynamics of the TPs in a growing MCS in a dissipative environment by neglecting inertial effects. The deformable but inert TPs, do not grow, divide or undergo apoptosis. The CCs grow and divide at a given rate, and undergo apoptosis (see Figure 1). The CCs and TPs experience systematic short-ranged attractive and repulsive forces arising from other CCs and TPs Malmi-Kakkada et al. (2018). To obtain the dynamics of the TPs and CCs, we introduce the density fields $\phi({\bf r},t)=\sum_{i}\delta[{\bf r}-{\bf r}_{i}(t)]$ for the CCs, and $\psi({\bf r},t)=\sum_{i}\delta[{\bf r}-{\bf r}_{i}(t)]$ for the TPs. A formally exact Langevin equation for $\phi({\bf r},t)$, and $\psi({\bf r},t)$ may be derived using the Dean’s method Dean (1996) accounting for diffusion and non-linear interactions.

[figure]style=plain,subcapbesideposition=top

To model a growing MCS, we modify the density equation for the CCs phenomenologically by adding a non-linear source term, $\propto\phi(\phi_{0}-\phi)$, accounting for cell birth and apoptosis, and a non-equilibrium noise term that breaks the CC number conservation. The noise, $f_{\phi}$, satisfies $<f_{\phi}({\bf r},t)f_{\phi}({\bf r^{\prime}},t^{\prime})>=\delta({\bf r}-{\bf r}^{\prime})\delta(t-t^{\prime})$. The source term, $\propto\phi(\phi_{0}-\phi)$, represents the birth and apoptosis, with $\phi_{0}=\frac{2k_{b}}{k_{a}}$Doering et al. (2003); Gelimson and Golestanian (2015). The coefficient , $\sqrt{k_{b}\phi+\frac{k_{a}}{2}\phi^{2}}$ of $f_{\phi}$, is the noise strength, corresponding to number fluctuations.

The $\phi({\bf r},t)$ and $\psi({\bf r},t)$ fields obey,

where ${\bf J}=\int_{\bf r^{\prime}}[\psi({\bf r^{\prime}},t)+\phi({\bf r^{\prime}},t)]{\bf\nabla}U({\bf r-\bf{r^{\prime}}})$, $\tilde{\eta}_{\psi}({\bf r},t)={\bf\nabla}\cdot\left(\eta_{\psi}({\bf r},t)\psi^{1/2}({\bf r},t)\right)$, $\tilde{\eta}_{\phi}({\bf r},t)={\bf\nabla}\cdot\left(\eta_{\phi}({\bf r},t)\phi^{1/2}({\bf r},t)\right)$, and $\eta_{\phi,\psi}$ satisfies $<\eta_{\phi,\psi}({\bf r},t)\eta_{\phi,\psi}({\bf r^{\prime}},t^{\prime})>=\delta({\bf r}-{\bf r}^{\prime})\delta(t-t^{\prime})$. The second term in Eq.(1) accounts for the TP-TP interactions (${\bf\nabla}\cdot\left(\psi({\bf r},t)\int_{\bf r^{\prime}}\psi({\bf r^{\prime}},t){\bf\nabla}U({\bf r-\bf{r^{\prime}}})\right)$) and TP-CC interactions, (${\bf\nabla}\cdot\left(\psi({\bf r},t)\int_{\bf r^{\prime}}\phi({\bf r^{\prime}},t){\bf\nabla}U({\bf r-\bf{r^{\prime}}})\right)$). The influence of the CCs on the TP dynamics is reflected in the TP-CC coupling. The third term in Eq.(2) results from cell birth and apoptosis, and the fourth term in Eq.(2) is the non-equilibrium noise. Eq. 1 does not satisfy the fluctuation-dissipation theorem in the tumor growth phase.

To obtain the exponents describing the dependence of the mean-square displacement (MSD) on time, we consider a change in scale, ${\bf r}\rightarrow s{\bf r}$, and $t\rightarrow s^{z}t$ where $z$ is the dynamic exponent. We use the Parisi-Wu Parisi and Wu (1981); Samanta et al. (2006); Samanta and Bhattacharjee (2006) technique to calculate $z$. The main features of the method, which are not needed to understand the central results here, are relegated to section I in the SI.

Simulations: We simulated a 3D MCS with embedded TPs using an agent-based model Drasdo and Hohme (2005); Malmi-Kakkada et al. (2018); Sinha et al. (2020). The cells, treated as interacting soft deformable spheres, grow with time, and divide into two identical daughter cells upon reaching a critical radius ($R_{m}$). The mean cell cycle time, $\tau=\frac{1}{k_{b}}$ is 15 hrs is used if not mentioned explicitly. The CCs also undergo apoptosis at the rate $k_{a}<<k_{b}$. The sizes and the number of TPs are held constant.

To determine the TP dynamics, we model the CC-CC, CC-TP, and TP-TP interactions using two potentials (Gaussian and Hertz) to ensure that the qualitative results are robust (simulation details are in section II in the SI).

The equation of motion for the dynamics of TP and CCs is taken to be Malmi-Kakkada et al. (2018), $\dot{\vec{r}}_{i}=\frac{\vec{F}_{i}}{\gamma_{i}}$, where $\dot{\vec{r}}_{i}$ is the velocity of $i^{th}$ CC or TP, $\vec{F}_{i}$ is the force on $i^{th}$ CC/TP, and $\gamma_{i}$ is the damping term. We used a pressure inhibition mechanism to model the observed growth dynamics in solid tumors Malmi-Kakkada et al. (2018). Dormancy or the growth phase of the CCs depends on the local microenvironment, determined by the pressure on the $i^{th}$ cell (Figure 1). Cell division and the placement of the daughter cells add stochasticity in the dynamics of the CCs. We initiated the simulations with 100 TPs and 100 CCs. The spatial coordinates of the CCs and TPs were sampled using a normal distribution with mean zero and standard deviation, $50~{}\mu m$.

[figure]style=plain,subcapbesideposition=top

Results: Theory and simulations predict that, in the limit $t<\frac{1}{k_{b}}$, the TPs exhibit sub-diffusive behavior (Table I). The non-linear terms for the TP-TP and TP-CC interactions determine the scaling laws. A power counting analysis shows that the dynamical exponent $z=2+\frac{d}{2}$. As a consequence, the MSD of a TP behaves as, $\Delta(t)_{TP}=<[r(t)-r(0)]^{2}>\sim t^{2/z}=t^{\beta_{TP}}$. Because $\beta_{TP}$ is less than unity we surmise that the TPs are jammed.

We also calculated $\Delta_{TP}(t)=\langle[{\bf r}(t)-{\bf r}(0)]^{2}\rangle$, by averaging over $\approx 2,000$ trajectories, for the Gaussian (Figure 2) and Hertz (Figure S2a in the SI) potentials. Similar behavior is found in both cases. The duration of the plateau (Figure 2) increases as the cell cycle time increases. Although jamming behavior predicted theoretically is consistent with the simulations the $\beta_{TP}$ exponents differ (Table I).

[figure]style=plain,subcapbesideposition=top

In the $t>\frac{1}{k_{b}}$ limit, theory and simulations predict hyper-diffusive dynamics, $\Delta_{TP}\sim t^{\alpha_{TP}}$ with $\alpha_{TP}>2$. Non-linearities, due to the TP-CC interactions together with cell division and apoptosis, determine the long-time behavior of the TPs. Using the scaling analysis as before, the theory predicts hyper-diffusive dynamics for the TPs (Table I). Variations in $k_{b}$ do not change the value of $\alpha_{TP}$. It merely changes the coefficients of the linear term and the value of $\frac{2k_{b}}{k_{a}}$. Therefore, $\alpha_{TP}$ is independent of the cell cycle time in the long time limit. Simulation results are in excellent agreement with the theoretical predictions. For times exceeding the cell cycle time, ($t>\tau$), the TPs exhibit hyper-diffusive behavior (Figure 2). For the Gaussian potential, we obtain $\alpha_{TP}\approx$ 2.3 ( see the inset in Figure 2) for $0.5\tau$, $\tau$ and $2\tau$. Thus, $\alpha_{TP}$ is independent of cell cycle time. Note that TP-TP interactions play an insignificant role in the dynamics of TPs or the CCs (Figure S3 in the SI). They merely alter the amplitude of $\Delta_{TP}$ in the intermediate time. They do not affect the long time dynamics

[figure]style=plain,subcapbesideposition=top

Figure 3 shows that changing the tracer size affects only the amplitude of $\Delta_{CC}$ in the intermediate time without altering the $\alpha_{TP}$ values (see also figures S4 and S5 in the SI). This is because the TP size only changes the nature of the short-range interactions without introducing any new scale. Since, $\alpha_{CC}$ is a consequence of the long-range spatial and temporal correlations that emerge because of the inequality $k_{b}>>k_{a}$, it is independent of the tracer size (section IIIE of SI).

In the absence of the TPs, birth and apoptosis determine the CC dynamics in the long time regime, yielding $\alpha_{CC}=1.33$ Malmi-Kakkada et al. (2018); Samanta and Thirumalai (2019). When TPs are present, the CCs continues to exhibit super-diffusive motion with a modest increase in $\alpha_{CC}$ (Table I). The CC pair-correlation, $g(r)=\frac{V}{4\pi r^{2}N^{2}}\sum_{i=1}^{N}\sum_{j\neq i}^{N}\delta(r-|{\bf r}_{i}-{\bf r}_{j}|)\sim r^{-0.5}$, in the presence and absence of TPs at $t\approx 8\tau$ (Figure 4a). The dynamically-induced CC correlations is independent of the TPs, thus explaining the insignificant effect of the TPs on the CC dynamics.

To explain the finding, $\alpha_{TP}>\alpha_{CC}$, we calculated the force-force autocorrelation function, FFA=$\langle{\bf F}(t+t_{d})\cdot{\bf F}(t)\rangle_{t}$. Here, $\langle...\rangle_{t}$ is the time average and $t_{d}$ is a delay time. Since, the TPs (CCs) exhibit hyper-diffusion (super-diffusion), we expect that the FFA of TPs should decay slower relative to the CCs, which is confirmed in Figure 4b, which shows that the FFA for the TPs (CCs) decays as $t_{d}^{-0.7}$ ($t_{d}^{-1.2}$). Thus, the TP motion is significantly more persistent than the CCs, which explains the hyper-diffusive nature of the TPs.

[figure]style=plain,subcapbesideposition=top

A mechanistic explanation for $\alpha_{TP}>\alpha_{CC}$ can be gleaned from the simulations. To illustrate the difference between the dynamics of the TPs and CCs, we calculated the angle ($\theta$) between two consecutive time steps in a trajectory. We define $\cos(\theta(t,\delta t))=\frac{\delta{\bf r}(t+\delta t)\cdot\delta{\bf r}(t)}{|\delta{\bf r}(t+\delta t)||\delta{\bf r}(t))|}$, where $\delta{\bf r}(t)={\bf r}(t+\delta t)-{\bf r}(t)$. Figure 5a shows the ensemble and time averaged $\cos\theta$ distribution for $\frac{\delta t}{\tau}=1$. If the motion of the CCs and TPs were diffusive, the distribution of $\cos\theta$ would be uniformly distributed from -1 to 1. However, Figure 5a shows that $P(\cos\theta)$ is skewed towards unity implying that the motion of both the TPs and CCs are persistent. Interestingly, the skewness is more pronounced for the TPs compared to the CCs (Figure 5a). During every cell division, the motion of CCs is randomized and hence the persistence is small compared to TPs. This explains the hyper-diffusive (super-diffusive) for the TPs (CCs) and thus $\alpha_{TP}>\alpha_{CC}$ in the long time limit (see Figure S6 which shows TPs move much more straighter than CCs).

To assess the effectiveness of the TPs as sensors of the local microenvironment, we calculated the radial pressure (Eq S28 in the SI) profile at $t=7.5\tau$, which has been measured Dolega et al. (2017). We measured pressure profiles using Eq. S28 (in the SI) for the system with both CCs and TPs. We also calculated pressure experienced just by the TPs. Finally, we computed pressure in a system consisting of solely the CCs. Two pertinent remarks about Figure (5b are worth making. (i) All three curves qualitatively capture the pressure profile found in the experiments. Pressure decreases roughly by factor of four, as the distance $r$ from the center of the tumor increases. The pressure is almost constant in the core, with a decrease that can be fit using the logistic function, as the boundary of the tumor is reached (Figure (5b). The high core pressure is due to small number of cell divisions. As a consequence, the CCs are jammed, leading to high internal pressure. As $r$ increases, the CCs proliferate, resulting in a decrease in self generated stress, and consequently a decreases in the pressure (Figure (5b).(ii) Most importantly, the CC pressure profiles are unaltered even in the presence of the TPs (compare blue circles and greenish yellow squares in Figure (5b), which shows that the latter does faithfully report the microenvironment of the CCs.

We conclude with a few comments. (i) The excellent agreement for $\alpha_{TP}$ between simulations and theory suggests that the CC dynamics is determined by the overall tumor growth, and not the details of the TP-CC interactions. Although $\alpha_{TP}>\alpha_{CC}$, the CCs move a larger distance than the TPs (compare Figure 2 and Figure 3), which implies that the effective diffusion constant of the CCs is greater than the TPs. (ii) The radial dependence of the pressure (Figure 5b), with the core experiencing higher values than the periphery cells, is another manifestation of non-equilibrium that is induced by an imbalance in division and apoptosis rates of the CCs. (iii) Most importantly, we find that the long time dynamics of CCs and the pressure profiles are unchanged when the TPs are embedded in the MCS. Therefore, it is justified to use the TPs as tumor stress sensors.

Acknowledgements We thank Mauro Mugnai, Hung Nguyen, Rytota Takaki and Davin Jeong for useful discussions and comments on the manuscript. This work is supported by the National Science Foundation (PHY 17-08128), and the Collie-Welch Chair through the Welch Foundation (F-0019).