Breast cancer invasion into adipose tissue strongly influences disease progression and metastasis. The degree of cancer cell invasion into adipose tissue depends on both biochemical signaling and the mechanical properties of cancer cells, adipocytes, and other key components of adipose tissue. We model breast cancer invasion into adipose tissue using discrete element method simulations of active, cohesive spherical particles (cancer cells) invading into confluent packings of deformable polyhedra (adipocytes). We quantify the degree of invasion by calculating the interfacial area *A _{t}* between cancer cells and adipocytes. We determine the long-time value of

*A*vs the activity and strength of the cohesion between cancer cells, as well as the mechanical properties of the adipocytes and extracellular matrix in which adipocytes are embedded. We show that the degree of invasion collapses onto a master curve as a function of the dimensionless energy scale

_{t}*E*, which grows linearly with the cancer cell velocity persistence time and fluctuations, is inversely proportional to the system pressure, and is offset by the cancer cell cohesive energy. When $ E c > 1$, cancer cells will invade the adipose tissue, whereas for $ E c < 1$, cancer cells and adipocytes remain de-mixed. We also show that

_{c}*A*decreases when the adipocytes are constrained by the ECM by an amount that depends on the spatial heterogeneity of the adipose tissue.

_{t}## INTRODUCTION

Breast cancer cell invasion into mammary adipose tissue is a key step toward disease progression and metastasis.^{1,2} Gaining insight into the biochemical and biophysical processes by which cancer cells invade adipose tissue is crucial for advancing our knowledge of breast cancer and improving treatments. During invasion, cancer cells are exposed to environmental cues as they migrate through the narrow, collagen-rich interstitial space between individual adipocytes that comprise adipose tissue.^{3,4} Indeed, the dynamics of breast cancer cell invasion are known to depend on important biochemical factors in the tumor microenvironment (TME).^{5–7} For example, cancer cells secrete chemical signals that promote the growth of blood vessels, which in turn provide oxygen and other nutrients to the tumor.^{6} Oscillations in the oxygen concentration in the TME caused by collective diffusion affect tumor growth and invasion.^{5} Moreover, genetic and chemical cues can trigger the epithelial–mesenchymal transition (EMT), which decreases the cohesion between cancer cells and increases their motility.^{7} In the context of adipose tissue, interactions between cancer cells and adipocytes cause lipid loss and de-differentiation in adipocytes, which promote invasion by altering cancer cell metabolism and reconfiguring the stroma.^{4,8}

However, the *physical* properties of the cancer cells and adipose tissue also play important roles in mediating breast cancer invasion.^{4,9–14} In adipose tissue, ECM fiber alignment is influenced by adipocyte properties and regulates tumor cell invasion by controlling migration persistence.^{14,15} More persistently moving cancer cells, in turn, can invade adipose tissue more rapidly. Cancer cell cohesion also influences tumor growth and invasion.^{9,10} Highly cohesive cancer cells do not disconnect from each other and thus do not as easily leave the primary tumor. With decreasing cohesion, cancer cells can more readily detach from the tumor and invade distant tissues.^{9,12} During invasion, cancer cells navigate between adipocytes, which form dense polyhedral packings. Thus, the packing fraction and shape deformation of adipocytes, as well as the ECM network between adipocytes, also physically influence breast cancer invasion.^{14,16}

In this work, we investigate the key physical variables that influence cancer cell invasion into adipose tissue. We leverage the recently developed deformable particle model (DPM)^{17–23} and carry out DPM simulations of cancer cells invading adipose tissue, which we model as dense packings of deformable polyhedra that are tethered to each other to mimic adhesion to the ECM. Cancer cells are modeled as soft, self-propelled, cohesive spherical particles since the magnitude of cancer cell shape fluctuations is much smaller than that of adipocytes. In addition, the cancer cells in our model are allowed to slightly overlap with the adipocytes as they invade into the interstitial regions between adipocytes, mimicking the elongation of cancer cells moving in confined regions. The systems are de-mixed when initialized, with a smooth interface between the cancer cells and adipose tissue. During the simulations, we quantify the invasion process by calculating the interfacial area shared by the cancer cells and adipocytes as a function of time. In our computational studies, we do not consider pressure gradients caused by cell proliferation but instead focus on invasion due to cancer cell migration. These novel simulations allow us to study the degree of tumor invasion as a function of the cohesion and motility of the cancer cells, as well as the structural and mechanical properties of the adipocytes and ECM.

As shown in previous studies, we find that a strong cohesion between cancer cells hinders tumor invasion, while increased cell motility promotes invasion. We identify a characteristic dimensionless energy scale *E _{c}* that controls the degree of invasion.

*E*increases with the speed and persistence of cancer cell migration and decreases with the pressure in the system. To model the adhesion of adipocytes to the ECM, we add springs connecting neighboring adipocytes to constrain their motion. We show that tethering the adipocytes decreases invasion compared to untethered adipocytes by an amount that is controlled by the packing fraction of the adipose tissue. Finally, we show that spatial heterogeneity in the mechanical properties of adipose tissue impedes invasion relative to adipose tissue with uniform mechanical properties.

_{c}We describe these findings in more detail in the remaining three sections of the article. In the Results section, we define the interfacial area *A _{t}* shared by the cancer cells and adipocytes and show how it depends on the dimensionless energy scale

*E*. We then show that

_{c}*A*decreases with tethering and heterogeneity in the mechanical properties of the adipocytes. In the Conclusions and Discussion section, we emphasize that the physical properties of adipose tissue, not only cancer cells, influence breast cancer invasion. We also propose future experimental studies of breast cancer cell invasion into soft spherical particles embedded in collagen fibers that will allow us to determine values of

_{t}*E*for systems that mimic the physical properties of adipose tissue. In the Methods section, we outline the DPM shape-energy function for adipocytes and the soft particle model for the cancer cells, the interactions between adipocytes, between cancer cells, and between adipocytes and cancer cells, as well as the integration of the equations of motion for the cancer cells and adipocytes. We also include three Appendices. In Appendix A, we show that our chosen reference shape parameter $ A 0$ for adipocytes is similar to those found in recent experimental studies of adipose tissue. In Appendix B, we further describe the interactions between cancer cells and adipocytes in simulations. Finally, in Appendix C, we describe how the interfacial area

_{c}*A*between cancer cells and adipocytes is normalized to remove the effects of dilation at large

_{t}*E*.

_{c}## RESULTS

In this section, we describe the results from discrete element method simulations of breast cancer invasion into adipose tissue. First, we describe how we quantify the degree of invasion using the interfacial area *A _{t}* between cancer cells and adipocytes. We then calculate

*A*as a function of time as the cancer cells and adipocytes evolve from a de-mixed to mixed state over a wide range of pressures and temperatures. If the system is above the glass transition temperature, the purely repulsive, passive cancer cells evolve toward a mixed (invaded) state. For attractive cancer cells, the long-time degree of invasion depends on the strength of the cancer cell cohesion. We identify a dimensionless energy scale, where $ E c > 1$ indicates the onset of cancer cell invasion. We find that

_{t}*E*is inversely proportional to the persistence of cancer cell motion and increases with the strength of attraction. Next, we study the effects of the ECM on cancer invasion by adding linear springs between neighboring adipocytes to constrain their motion. We find that constraining adipocytes does not alter the onset of invasion but decreases the value of

_{c}*A*relative to untethered adipocytes. Finally, we show that heterogeneity in the local structural and mechanical properties of tethered adipocytes restricts cancer cell invasion.

_{t}### Quantifying the degree of cancer cell invasion

As discussed in the Methods section, cancer cells and adipocytes are initialized in a de-mixed state. We then integrate their equations of motion [i.e., Eqs. (11) or (15) for cancer cells and Eq. (16) for adipocytes] to determine their spatiotemporal trajectories over a wide range of temperatures and pressures. To illustrate how we quantify the degree of cancer cell invasion, we first focus on the passive model with purely repulsive interactions for cancer cells [Eq. (15)]. We quantify the invasion process by calculating the interfacial area *A _{t}* between the cancer cells and adipocytes (as well as between cancer cells and the walls) as a function of time. As discussed in Appendix C,

*A*is determined by performing a Voronoi tessellation and identifying the faces of the Voronoi polyhedra that are shared by cancer cells and adipocytes (and cancer cells and the walls). In Fig. 1(a), we show that the time dependence of

_{t}*A*strongly varies with $ k b T / ( P \sigma 3 )$. The minimum possible value of

_{t}*A*is $ A t min = 470 \sigma 2$, which corresponds to the initial, de-mixed value [Fig. 1(b)], and the maximum possible value of

_{t}*A*is $ A t max = 2000 \sigma 2$, which corresponds to the surface area of the

_{t}*N*adipocytes [Fig. 1(d)]. Thus,

_{a}*A*must exist between these two values during the cancer cell invasion process. To enable comparisons of the invasion process for different system parameters, we normalize

_{t}*A*by its minimum and maximum values while accounting for the dilation of the system at high temperatures, as discussed in Appendix C. In Fig. 2(a), we show the normalized interfacial area

_{t}*A*as a function of $ k b T / ( P \sigma 3 )$, where $ 0 < A n < 1$, for purely repulsive, passive cancer cells for different times following initialization of the de-mixed state. For each time

_{n}*t*,

*A*is sigmoidal when plotted vs $ k b T / ( P \sigma 3 )$, but shifts to smaller $ k b T / ( P \sigma 3 )$ as

_{n}*t*increases. These results indicate that purely repulsive, passive cancer cells will mix with adipocytes, reaching

*A*= 1 for all $ k b T / ( P \sigma 3 )$ that are above the glass transition temperature, where structural relaxation times diverge.

_{n}### Effects of cancer cell cohesion and persistence on invasion

*A*is time dependent, evolving toward the mixed (invaded) state for purely repulsive, passive cancer cells. Next, we examined how

_{n}*A*is affected by cancer cell cohesion and persistence. In Fig. 2(b), we show that

_{n}*A*vs $ k b T / ( P \sigma 3 )$ for systems with attractive, active cancer cells does not possess strong time dependence. Moreover, the value of

_{n}*A*does not depend on the starting condition, reaching a steady-state value for times $ t / m \sigma 2 / \epsilon c < 5 \xd7 10 5$ as shown in the inset to Fig. 2(b). In Fig. 3(a), we plot the steady-state

_{n}*A*vs $ k b T / ( P \sigma 3 )$ for three values of the attractive strength

_{n}*β*and five values of the persistence time

*τ*is sigmoidal vs $ k b T / ( P \sigma 3 )$ for all

_{P}. A_{n}*β*and

*τ*combinations studied. Thus, we can describe

_{P}*A*using

_{n}*A*occurs and $ 0.5 < b < 2.5$ controls the rate of this increase. We demonstrate high-quality fits of

_{n}*A*to Eq. (1) in Fig. 3(a). Thus, when we plot

_{n}*A*as a function of $ E c b$, where $ E c \u2261 k b T / ( a P \sigma 3 )$, the data collapse as shown in Fig. 3(b). We also varied the damping coefficient

_{n}*γ*to verify that Eq. (1) collapses the data for systems with both underdamped and overdamped dynamics. We next determine the functional form for $ a ( \beta , \tau p )$. First,

*a*for passive cancer cells ( $ \tau p = 0$) at a given

*T*and

*P*is similar to that for active cancer cells at small

*τ*. [See black/purple circles and black/purple diamonds in Fig. 3(a)]. Next,

_{p}*a*decreases with increasing

*τ*for large

_{p}*τ*, emphasizing that cancer cell persistence enhances invasion. Thus, a reasonable ansatz for

_{p}*a*at fixed

*β*is $ a \u223c 1 / \tau d$, where

*τ*is the decorrelation time of the velocity autocorrelation function for the cancer cells,

_{d}*τ*on cancer cell packing fraction and

_{d}*τ*, we show

_{p}*τ*for simplified systems composed of only repulsive, active cancer cells in periodic boundary conditions in Fig. 4(b). $ \tau d \u223c \tau \gamma $ in the $ \tau p \u2192 0$ limit for all $ \varphi c$, where $ \tau \gamma = m \gamma $. For intermediate values of

_{d}*τ*, $ \tau d \u223c \tau p$, which agrees with the analytical result at $ \varphi c = 0$,

_{p}^{24}

*τ*,

_{p}*τ*deviates from the analytical result as $ \varphi c$ increases. In the inset to Fig. 3(b), we verify that $ a \u223c \tau \gamma / \tau d$.

_{d}At fixed *τ _{p}*,

*A*shifts to larger $ k b T / ( P \sigma 3 )$ with increasing attraction strength $ \beta \epsilon c$ since cohesion makes it more difficult for cancer cells to invade the adipocytes. In the $ \beta \u2192 0$ limit,

_{n}*A*converges to a time-dependent value for purely repulsive cancer cells. Thus, we expect $ a \u223c c 1 ( t ) + c \beta \epsilon c$, where $ c 1 ( t ) \u223c k b T g / ( P \sigma 3 )$ in the long time limit and

_{n}*a*is not time dependent when $ c \beta \u226b c 1 ( t )$. We show in Fig. 3(b) that $ a = c 1 ( 1 + c 2 \beta \epsilon c / ( P \sigma 3 ) ) \tau \gamma / \tau d$, where $ c 1 \u2248 0.06$ enforces $ A n = 1 / 2$ at

*E*= 1 and $ c 2 \u2248 0.02$.

_{c}### Tethering of neighboring adipocytes

*i*and

*j*. Adipocytes are considered neighbors when the distance between the centers of mass of adipocyte pairs is smaller than a threshold value

*d*, such that the number of neighbors per adipocyte is 6. The potential energy associated with the spring connecting adipocytes

_{c}*i*and

*j*is

*i*share the interaction with adipocyte

*j*equally, and thus the force on the vertex

*μ*on adipocyte

*i*from adipocyte

*j*is given by

*A*vs $ k b T / ( P \sigma 3 )$ for several values of

_{n}*K*and two values of

_{ecm}*τ*and

_{p}*β*. As discussed above,

*A*vs $ k b T / ( P \sigma 3 )$ is normalized between 0 and 1 using the interfacial area obtained for untethered (

_{n}*K*= 0) systems. Tethering of adipocytes decreases the degree of invasion but not the onset value of $ k b T / ( P \sigma 3 )$ for invasion. We find that the plateau values $ A n max < 1$ at large $ k b T / ( P \sigma 3 )$ for all

_{ecm}*K*> 0, but the values of $ k b T / ( P \sigma 3 )$ at which $ A n > 0$ does not vary strongly with

_{ecm}*K*. In the insets to Fig. 5, we show that $ A n max$ decreases monotonically from 1 to lower values that depend on

_{ecm}*τ*and

_{p}*β*. Thus, the tethering of adipocytes by the ECM reduces the degree of cancer cell invasion into adipose tissue.

### Packing fraction heterogeneity

The mechanical properties of adipose tissue are affected by many factors including obesity, inflammation, and lipolysis.^{4,14,25} Previous experimental studies of cancer cell migration through collagen networks demonstrated that geometrical confinement strongly affects cancer cell motion.^{16,26} Thus, we expect the packing fraction of adipocytes to influence the degree of invasion. Adipocytes are initialized using athermal, quasistatic compression to reach a total packing fraction $ \varphi = 0.72$ as described in the Methods section. We now vary the adipocyte packing fraction over the range $ 0.4 < \varphi a < 0.9$ by changing the length of the simulation box and shifting the positions of the cancer cells and adipocytes affinely, followed by energy minimization. We then add linear springs between neighboring adipocytes, setting $ l i j 0 = R i j$ and $ K ecm = 0.04 \epsilon c$, and calculate *A _{n}* vs $ \varphi a$ for $ \tau p / m \sigma 2 \epsilon c = 25$ and $ \beta = 10 \u2212 2$. In Fig. 6(a), we confirm that

*A*decreases with increasing adipocyte packing fraction and vanishes for $ \varphi a \u2273 0.85$.

_{n}In addition to packing fraction, the local structural and mechanical properties of adipose tissue are often heterogeneous.^{27} We model structural and mechanical heterogeneity in adipose tissue by resetting $ l i j 0$ of the tethering springs for half of the adipocytes in the *z*-direction to $ l i j 0 = ( 1 \u2212 \lambda ) R i j$ and to $ l i j 0 = ( 1 + \lambda ) R i j$ for the other half of the adipocytes. In this idealized system, half of the adipocytes are tightly coupled with smaller equilibrium lengths, and half are loosely coupled with larger equilibrium lengths. This heterogeneity in equilibrium lengths gives rise to root mean square fluctuations in the local packing fraction that increase with *λ*. As shown in Fig. 6(b), the maximum values of $ \Delta \varphi $ are $ \u2248 0.1$, 0.2, and 0.25 for *λ* = 0, 0.3, and 1. In Fig. 6(a), we show that invasion is less pronounced for heterogeneous systems compared to more uniform systems at the same adipocyte packing fraction. This result can be explained given that for sufficiently large *λ*, tightly coupled adipocytes exhibit smaller gaps between adipocytes that are less than the diameter of the cancer cells, reducing the degree of invasion in half of the system. Since the adipocytes are tethered, cancer cells invading the loosely coupled adipocytes do not affect invasion in the tightly coupled half of the system.

## CONCLUSIONS AND FUTURE DIRECTIONS

In this work, we aimed to understand how the physical properties of cancer cells and adipose tissue influence breast cancer invasion. To this end, we developed discrete element method simulations of tumor invasion, modeling the cancer cells as cohesive, active, soft particles and the adipocytes as deformable, polyhedral particles. We initialized the system in a de-mixed state and quantified the degree of invasion as the interfacial area *A _{t}* shared between the cancer cells and adipocytes.

*A*is bounded by the area of the interface in the initial de-mixed state $ A t min$ and the total surface area of the adipocytes $ A t max$, and we normalized the interfacial area by these two values so that $ 0 < A n < 1$. We then calculated

_{t}*A*as a function of temperature, pressure, and the cohesive strength and persistence of cancer cell motion. For both underdamped and overdamped dynamics of the cancer cells, we showed that

_{n}*A*can be collapsed when plotted as a function of a dimensionless energy scale

_{n}*E*, where $ E c \u226a 1$ indicates a de-mixed (non-invaded) state and $ E c \u226b 1$ indicates a mixed (invaded) state. We found that

_{c}*E*increases with the mean square fluctuations and persistence of cancer cell velocities and is inversely related to cancer cell cohesion and pressure in the system. We also investigated the physical effects of the extracellular matrix on cancer cell invasion. We demonstrated that tethering of neighboring adipocytes and denser packing of adipocytes inhibit invasion. In addition, we showed that heterogeneity in the local packing fraction for tethered adipocytes decreases

_{c}*A*relative to more uniform local packing.

_{n}The current computational studies represent an important first step toward modeling breast cancer invasion into adipose tissue as a physical process. However, there are numerous future research directions to pursue. First, we accounted for some of the effects of the ECM by tethering neighboring adipocytes to each other. However, alignment of ECM fibers also strongly influences the speed and persistence of cancer cell velocities during invasion.^{16,28} In preliminary experimental studies, we have embedded plastic spherical beads within collagen matrices and quantified the local density and alignment of the collagen fibers, as well as the positions of the spherical beads. These studies have varied the bead packing fraction and collagen density and correlated enhanced cancer cell motion to collagen fiber alignment.^{29} In the current work, we assumed that the cancer cells share the same migration speed determined by *k _{bT}* and persistence time

*τ*, and that these properties are static. In future computational studies, we can incorporate separate fields into the simulations to account for local variations in collagen density and alignment caused by cancer cell motion, couple these fields to cancer cell velocities, and quantitatively compare the simulated trajectories of the cancer cells with experimental results from

_{p}*in vitro*systems where collagen fiber alignment can be controlled. Second, in the current simulations, we modeled cancer cells as active, soft spherical particles that do not change shape. However, we know that cancer cells can significantly alter their shapes to squeeze through narrow gaps and their motion has been correlated with shape elongation.

^{16,30,31}Thus, in future studies, we aim to model cancer cells using active,

*deformable*particles. Studies have also suggested that breast cancer invasion into adipose tissue can induce lipolysis in adipocytes.

^{4}To model this effect, we can allow the equilibrium volume $ v i 0$ of adipocyte

*i*to decrease over time when in contact with cancer cells.

In our computational studies, we assumed that the cancer cell invasion time scale $ \tau inv$ is shorter than the time scale of the cell cycle $ \tau cell$, and thus we did not model cell growth, division, and apoptosis. As a result, de-mixed (non-invaded) systems in our current studies could become mixed (invaded) due to cancer cell proliferation on time scales longer than $ \tau cell$. Moreover, recent studies demonstrated that coupling cell growth rate and local pressure generates an interface instability with a fingering invasion pattern even when cancer cells do not migrate.^{32} This prior work emphasizes that systems with $ \tau cell \u226a \tau inv$ are also relevant to cancer invasion. In future computational studies, we will model cancer cell growth, division, and apoptosis and determine the mechanism of cancer cell invasion and shape of the invasion front as a function of $ \tau cell / \tau inv$.

## METHODS

In this section, we describe the methods for simulating breast cancer cell invasion into adipose tissue. First, we introduce the shape-energy function for the 3D deformable particle model (DPM), which is used to model adipocytes. We then describe the self-propelled, attractive soft particle model for the cancer cells and the repulsive contact interactions between cancer cells and adipocytes. We outline the equations of motion for the cancer cells and adipocytes and describe the methods we use to integrate them to determine the cell trajectories. Finally, we specify the initial configuration as a de-mixed system with a smooth, planar interface between cancer cells and adipocytes, which also includes physical walls in one direction to control pressure and periodic boundary conditions in the other two directions.

### Deformable particle model for adipocytes

*N*= 42 vertices belonging to a mesh with

_{v}*N*= 80 triangular faces, which is sufficient to describe the shape fluctuations observed in adipocytes. [See Fig. 7(f)]. The shape energy for the

_{f}*i*th adipocyte as a function of the positions of its vertices $ r i 1 , \u2026 , r i N v$ is

The first term imposes a harmonic energy penalty for changes in cell volume *v _{i}* from its preferred value $ v i 0$ and

*ε*controls the fluctuations in adipocyte volume. We set

_{v}*ε*to be much larger than the other energy scales in the system so that the adipocytes cannot change their volume. The second term imposes a harmonic energy penalty for deviations in the area

_{v}*a*of each triangular face from its preferred value $ a i k 0$, and

_{ik}*ε*controls the fluctuations in adipocyte surface area. We quantify the adipocyte shape using the dimensionless shape parameter $ A = ( \u2211 k = 1 N f a i k ) 3 2 6 \pi v i$. The shape parameter for a sphere is 1 and $ A > 1$ for all nonspherical shapes. We set the reference shape parameter $ A 0 = 1.1$ for adipocytes based on experimental observations shown in Appendix A.

_{a}Modeling adipocytes using the deformable particle model has several advantages. In general, 3D deformable particles without bending energy are floppy and can change shape without increasing energy. In particular, deformable particles with *N _{v}* = 42 vertices and

*N*= 80 faces have $ \Delta N = ( 3 N v \u2212 6 ) \u2212 ( N f + 1 ) = 39$ zero-energy modes of the dynamical matrix for the shape-energy function in Eq. (7), allowing the adipocytes to change shape without energy cost, while maintaining their volume and surface area.

_{f}^{18}In addition, using the DPM allows packings of adipocytes to be nearly confluent without overlaps between adipocytes.

### Interactions between adipocytes, between cancer cells, and between adipocytes and cancer cells

*μ*on adipocyte

*i*and vertex

*ν*on adipocyte

*j*is purely repulsive,

*ε*controls the strength of the repulsive forces, $ r \u0302 ij \mu \nu = r ij \mu \nu / r ij \mu \nu , \u2009 r ij \mu \nu = r i \mu \u2212 r j \nu $ is the separation vector, $ r ij \mu \nu = | r ij \mu \nu |$ is the distance between vertices, and $ \sigma ij \mu \nu = 1 2 ( \sigma i \mu + \sigma j \nu )$ is their average diameter. (See Fig. 8 in Appendix B). $ \sigma i \mu $ is approximately two times the average cancer cell diameter, which allows us to cover the surface of the adipocytes with vertices to prevent interpenetration by cancer cells and other adipocyte vertices. Since vertex neighbor lists can be employed to calculate forces between adipocytes, the computational complexity of the DPM is only $ \u223c N v N a$.

_{c}^{33}while cancer cells are much smaller, ranging between 10 and $ 20 \u2009 \mu m$.

^{34}Since the magnitude of cancer cell shape fluctuations is much smaller than that for adipocytes, we do not consider shape fluctuations of cancer cells and instead model them as soft spheres. We assume that the force on vertex

*μ*on adipocyte

*i*from cancer cell

*q*is purely repulsive,

*μ*on adipocyte

*i*and cancer cell

*q*, $ \sigma i \mu q = 1 2 ( \sigma i \mu + \sigma q )$, and

*σ*is the diameter of the

_{q}*q*th cancer cell.

*q*from cancer cell

*s*obeys

*q*and

*s*, and $ \sigma q s = 1 2 ( \sigma q + \sigma s )$ is the average diameter of the cancer cells. As shown in Fig. 9, the pair force between cancer cells is repulsive for $ r q s < \sigma q s$, and attractive for $ 1 < r q s / \sigma q s < r \alpha / \sigma q s$ with maximum attractive force $ F q s = \epsilon c \beta / \sigma q s$ at $ r q s = r \beta $, where $ r \alpha / \sigma q s = 1 + \alpha $ and $ r \beta / \sigma q s = 1 + \beta $. We vary the depth of the attractive potential

*β*from $ 10 \u2212 4$ to $ 10 \u2212 2$ and set the range of attractive interactions at $ \alpha = 0.2$. To prevent partial crystallization of the cancer cells during invasion, we assume that the cancer cells are bidisperse in size where half of the cells have a smaller diameter $ 0.9 \sigma $ and the other half have a larger diameter $ 1.1 \sigma $, where

*σ*is the average diameter of the cancer cells. The ratio of the adipocyte diameter (including the vertex radii) and the average cancer cell diameter

*σ*is approximately 6, which matches recent experimental studies.

^{33,34}

### Equations of motion for cancer cells and adipocytes

**r**

_{q}of cancer cell

*q*is

*m*is the mass of the cancer cells, $ F i \mu q$ is the pair force between cancer cells and the vertices of the adipocytes [Eq. (9)], and $ F q s$ is the pair force between cancer cells [Eq. (10)]. $ F w q$ indicates the repulsive forces between the physical walls and the

*q*th cancer cell and is defined in Eq. (18). The equation of motion for active cancer cells also includes a damping force from the surrounding environment that is proportional to the cancer cell velocity

**v**

_{q}with damping coefficient $ \gamma \sigma / m \epsilon c = 0.2$ and a self-propulsion term with active force magnitude

*f*

_{0}. The direction of the active force $ n \u0302 q$ obeys rotational diffusion,

^{35–37}

^{,}

*x*-,

*y*-, and

*z*-components that are Gaussian random numbers with zero mean and satisfy

**I**is the identity matrix, $ \delta q s = 1$ when

*q*=

*s*and 0 otherwise, and $ \delta ( . )$ is the Dirac

*δ*-function. $ \tau p$ is the persistence time of the active force as defined through the autocorrelation function of $ n \u0302 q$,

^{28}Experimental studies demonstrate that tumor cells locally align ECM fibers during migration by generating traction forces that cause plastic deformation of ECM fibers.

^{38,39}More aligned ECM increases the migration speed and persistence of subsequent tumor cells. In this model, we include a background temperature

*k*that captures the random motion of the cancer cells and persistence time

_{bT}*τ*that captures the correlation between cancer cell velocity and local ECM fiber alignment. Both

_{p}*k*and

_{bT}*τ*describe the migration properties of cancer cells, and implicitly model the degree of alignment of the ECM.

_{p}*q*th cancer cell,

**v**

_{q}of cancer cell

*q*with damping coefficient

*γ*. The last term is the noise term with strength $ 2 k b T 0 \gamma m$ that sets the target temperature

*T*

_{0}. $ \eta ( t )$ is a vector with

*x*-,

*y*-, and

*z*-components that are Gaussian random numbers with zero mean and unit variance.

*μ*on the

*i*th adipocyte obeys

*m*is the mass and

*γ*is the damping constant associated with each adipocyte vertex. The first and second terms give the repulsive pair forces between vertices

*μ*and

*ν*on different adipocytes

*i*and

*j*and between vertex

*μ*on adipocyte

*i*and cancer cell

*q*. $ F w i \mu $ is the repulsive force between the physical walls and vertex

*μ*on the

*i*th adipocyte, which is defined in Eq. (19). The fifth term gives the mechanical forces generated from the shape-energy function for the

*i*th adipocyte in Eq. (7).

### Initialization and boundary conditions for invasion simulations

We initialize the cancer cells and adipocytes in a dilute, de-mixed state. The simulation box is cubic with initial side length $ L 0 = 70 \sigma $, *N _{a}* = 28 adipocytes,

*N*= 1500 cancer cells, and initial total packing fraction $ \varphi = 0.01$. (Note the total packing fraction can be written in terms of the packing fractions for cancer cells and adipocytes, separately, $ \varphi = \varphi c + \varphi a$.) The cancer cells and adipocytes are placed randomly in the simulation box such that the

_{c}*x*-components of the cancer cell centers of mass satisfy $ 2 \sigma < x q < 12 \sigma $ and the

*x*-components of the adipocyte centers of mass satisfy $ 16 \sigma < X i < 68 \sigma $. We then perform athermal, isotropic compression, applying successive compression steps $ \Delta \varphi / \varphi = 0.03$ to generate a jammed packing of adipocytes. Compression is followed by energy minimization until the final packing fraction of the system reaches $ \varphi = ( \pi 6 \u2211 q = 1 N c \sigma q 3 + \u2211 i = 1 N a V i ) / L f 3 = 0.72$, where $ L f = 34 \sigma $ is the final side length of the box, and

*V*is the total volume of the

_{i}*i*th adipocyte, i.e.,

*v*plus the additional volumes from the vertices that do not overlap the polyhedron. [See Fig. 7(e)].

_{i}*y*- and

*z*-directions. In the

*x*-direction, we include two physical walls, a “right” wall at

*x*=

_{r}*L*and a “left” wall at

_{f}*x*= 0. The right wall is stationary, while the left wall moves in the

_{l}*x*-direction to maintain fixed pressure

*P*. The equation of the motion for the left wall position is

*m*is the mass and

*γ*is the damping coefficient for the left wall,

*v*is the speed of the left wall. The forces that act on the wall from the

_{l}*q*th cancer cell and

*μ*th vertex on the

*i*th adipocyte are

*w*=

*l*,

*r*indicates the left and right walls, respectively.) We quantify the temperature of the system using the kinetic energy of the cancer cells,

For the passive model of cancer cells, the temperature of the cancer cells *T _{c}* and adipocytes

*T*equilibrate. For the active model of cancer cells, the temperature of the cancer cells is higher than that of the adipocytes. In the systems considered here, $ 1 < T c / T a < 1.2$.

_{a}## ACKNOWLEDGMENTS

We acknowledge support from the National Institutes of Health under Grant No. R01CA276392 (Y.Z., D.W., and C.S.O.), the National Science Foundation under Grant No. DGE1650441 (G.B.), the National Cancer Institute under Grant No. F31CA278410 (G.B.), the National Cancer Institute under Grant No. R01CA259195 (C.F.), and the Center on the Physics of Cancer Metabolism under Grant No. 1U54CA210184 (C.F.). This work was also supported by the High Performance Computing facilities operated by Yale's Center for Research Computing.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Ethics Approval

Ethics approval is not required.

### Author Contributions

**Yitong Zheng:** Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Dong Wang:** Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Software (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). **Garrett Beeghly:** Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Writing – review & editing (supporting). **Claudia Fischbach:** Conceptualization (equal); Funding acquisition (equal); Resources (supporting); Supervision (supporting); Writing – review & editing (supporting). **Mark D. Shattuck:** Conceptualization (supporting); Formal analysis (supporting); Methodology (supporting); Supervision (supporting). **Corey S. O'Hern:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: SHAPE PARAMETER OF ADIPOCYTES

Due to the highly scattering nature of lipids, *in situ* 3D visualization of adipose tissue remains challenging. Instead, the extent of breast cancer invasion into adipose tissue is routinely assessed using 2D histology slices of tumors as shown in Figs. 7(a) and 7(b). We analyzed $ \u223c 1000$ adipocytes from five different mice with mammary tumors 11 days after implantation. For each slice, we identify the adipocytes and calculate the shape parameter of the cross sections, $ A 2 D = p 2 / ( 4 \pi s )$, where *p* is the perimeter and *s* is the area of the cross section. In Fig. 10, we compare the distribution $ P ( A 2 D )$ from the 2D histology slices to the shape parameter distributions from random 2D slices through adipocytes from the DEM invasion simulations at several reference shape parameters $ A 0$ in 3D. We find that using $ A 0 \u2248 1.1$ gives reasonable agreement for $ P ( A 2 D )$ between the experiments and simulations.

### APPENDIX B: INTERACTIONS BETWEEN CANCER CELLS AND ADIPOCYTES

In Eqs. (8)–(10) in the Methods section, we define the pair forces between vertices on contacting adipocytes, between a vertex on an adipocyte and a contacting cancer cell, and between contacting cancer cells. Key parameters in these expressions include the separations $ r ij \mu \nu , \u2009 r i \mu q$, and *r _{qs}* and the diameters $ \sigma i \mu $ and

*σ*, which are pictured in Fig. 8.

_{q}### APPENDIX C: NORMALIZATION OF INTERFACIAL AREA *A*_{n}

_{n}

In this Appendix, we describe how we calculate the interfacial area *A _{t}* shared by cancer cells and adipocytes and how we define the normalized interfacial area $ 0 < A n < 1$. For a given configuration of cancer cells, adipocytes, and locations of the confining walls in the

*x*-direction, we perform Voronoi tessellation of the positions of the cancer cells and vertices of the adipocytes, as shown in Figs. 11(a)–11(c). Next, we find all of the Voronoi polyhedral faces that are shared by a Voronoi polyhedron associated with a cancer cell and by a Voronoi polyhedron associated with a vertex of an adipocyte, as well as all of the faces of cancer cell Voronoi polyhedra that form the physical walls as shown in Figs. 11(d)–11(f). The interfacial area (between cancer cells and adipocytes and between cancer cells and the physical boundaries)

*A*is defined as the sum of the areas of these faces. The maximum value of

_{t}*A*is $ A t max = 2 L f 2 + \Sigma i = 1 N a S i$, where

_{t}*S*is the surface area of the Voronoi polyhedron for adipocyte

_{i}*i*. The relative interfacial area $ A t / A t max$ serves as an indicator of the degree of cancer cell invasion.

Note that the volume of the simulation box, as well as *A _{t}* and $ A t max$, can increase with $ k b T / ( P \sigma 3 )$, especially at large values when $ k b T$ overcomes the cohesive energy $ \beta \epsilon c$. In this regime, the ratio $ A t / A t max$ can vary with $ k b T / ( P \sigma 3 )$. Since we are not interested in changes in the relative interfacial area $ A t / A t max$ caused by the dilation of the simulation box, we define $ r A t = ( A t \u2212 \Delta A ) / ( A t max \u2212 \Delta A )$, where $ \Delta A = A t max \u2212 min ( A t max )$ and $ min ( A t max )$ is the value of $ A t max$ at which dilation occurs. We then define the normalized value $ A n = ( r A t \u2212 min ( r A t ) ) / ( max ( r A t ) \u2212 min ( r A t ) )$, where $ min ( r A t )$ and $ max ( r A t )$ are the minimum and maximum values of

*rA*, such that $ 0 < A n < 1$.

_{t}## REFERENCES

*Angiogenesis—A Self-Adapting Principle In Hypoxia*

*Methods in Enzymology*