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 At between cancer cells and adipocytes. We determine the long-time value of At 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 Ec, 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 At decreases when the adipocytes are constrained by the ECM by an amount that depends on the spatial heterogeneity of the adipose tissue.

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 Ec that controls the degree of invasion. Ec 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.

We describe these findings in more detail in the remaining three sections of the article. In the Results section, we define the interfacial area At shared by the cancer cells and adipocytes and show how it depends on the dimensionless energy scale Ec. We then show that At 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 Ec 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 At between cancer cells and adipocytes is normalized to remove the effects of dilation at large Ec.

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 At between cancer cells and adipocytes. We then calculate At 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 Ec 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 At relative to untethered adipocytes. Finally, we show that heterogeneity in the local structural and mechanical properties of tethered adipocytes restricts 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 At 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, At 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 At strongly varies with k b T / ( P σ 3 ). The minimum possible value of At is A t min = 470 σ 2, which corresponds to the initial, de-mixed value [Fig. 1(b)], and the maximum possible value of At is A t max = 2000 σ 2, which corresponds to the surface area of the Na adipocytes [Fig. 1(d)]. Thus, At 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 At 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 An as a function of k b T / ( P σ 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 t, An is sigmoidal when plotted vs k b T / ( P σ 3 ), but shifts to smaller k b T / ( P σ 3 ) as t increases. These results indicate that purely repulsive, passive cancer cells will mix with adipocytes, reaching An = 1 for all k b T / ( P σ 3 ) that are above the glass transition temperature, where structural relaxation times diverge.

FIG. 1.

(a) The extent of tumor invasion is quantified by plotting the shared interfacial surface area At (in units of σ 2) between the cancer cells and adipocytes as a function of time (in units of m σ 2 / ε c). We show At for five values of k b T / ( P σ 3 ), representing tumor cells with varying motility increasing from blue to red. Simulation snapshots for different At values and times from (a) are shown in (b)–(d): (b) A t / σ 2 = 610 (c) 1280, and (d) 1720. The time of each snapshot is indicated by the corresponding shapes shown in (a): (b) triangle, (c) square, and (d) circle.

FIG. 1.

(a) The extent of tumor invasion is quantified by plotting the shared interfacial surface area At (in units of σ 2) between the cancer cells and adipocytes as a function of time (in units of m σ 2 / ε c). We show At for five values of k b T / ( P σ 3 ), representing tumor cells with varying motility increasing from blue to red. Simulation snapshots for different At values and times from (a) are shown in (b)–(d): (b) A t / σ 2 = 610 (c) 1280, and (d) 1720. The time of each snapshot is indicated by the corresponding shapes shown in (a): (b) triangle, (c) square, and (d) circle.

Close modal
FIG. 2.

Normalized interfacial area An between the cancer cells and adipocytes (defined in  Appendix B) plotted as a function of k b T / ( P σ 3 ) at different times t, over a wide range of temperatures T, and where pressure P = 1.3 × 10 3 ε c σ 3 for (a) purely repulsive, passive cancer cells with attraction strength β = 0 and persistence time τ p = 0 and (b) attractive, active cancer cells with β = 10 2 and τ p / m σ 2 ε c = 25. In (a), the vertical line indicates the glass transition temperature k b T g / ( P σ 3 ) 0.03 below which the structural relaxation τ r . In (b), the inset shows that the equilibrium value of A n 0.55 (black dashed line) for k b T / ( P σ 3 ) 0.3 is independent of the initial state since the same value of An is obtained when the system is initialized in a de-mixed state (blue line) and a mixed state (red line). In both (a) and (b), time increases from green to blue.

FIG. 2.

Normalized interfacial area An between the cancer cells and adipocytes (defined in  Appendix B) plotted as a function of k b T / ( P σ 3 ) at different times t, over a wide range of temperatures T, and where pressure P = 1.3 × 10 3 ε c σ 3 for (a) purely repulsive, passive cancer cells with attraction strength β = 0 and persistence time τ p = 0 and (b) attractive, active cancer cells with β = 10 2 and τ p / m σ 2 ε c = 25. In (a), the vertical line indicates the glass transition temperature k b T g / ( P σ 3 ) 0.03 below which the structural relaxation τ r . In (b), the inset shows that the equilibrium value of A n 0.55 (black dashed line) for k b T / ( P σ 3 ) 0.3 is independent of the initial state since the same value of An is obtained when the system is initialized in a de-mixed state (blue line) and a mixed state (red line). In both (a) and (b), time increases from green to blue.

Close modal
In the previous section, “Quantifying the degree of cancer cell invasion,” we found that the degree of invasion An is time dependent, evolving toward the mixed (invaded) state for purely repulsive, passive cancer cells. Next, we examined how An is affected by cancer cell cohesion and persistence. In Fig. 2(b), we show that An vs k b T / ( P σ 3 ) for systems with attractive, active cancer cells does not possess strong time dependence. Moreover, the value of An does not depend on the starting condition, reaching a steady-state value for times t / m σ 2 / ε c < 5 × 10 5 as shown in the inset to Fig. 2(b). In Fig. 3(a), we plot the steady-state An vs k b T / ( P σ 3 ) for three values of the attractive strength β and five values of the persistence time τP. An is sigmoidal vs k b T / ( P σ 3 ) for all β and τP combinations studied. Thus, we can describe An using
A n = 1 2 ( tanh [ log 10 ( k b T a P σ 3 ) b ] + 1 ) ,
(1)
where a ( β , τ p ) controls the value of k b T / ( P σ 3 ) at which the rapid increase in An occurs and 0.5 < b < 2.5 controls the rate of this increase. We demonstrate high-quality fits of An to Eq. (1) in Fig. 3(a). Thus, when we plot An as a function of E c b, where E c k b T / ( a P σ 3 ), the data collapse as shown in Fig. 3(b). We also varied the damping coefficient γ to verify that Eq. (1) collapses the data for systems with both underdamped and overdamped dynamics. We next determine the functional form for a ( β , τ p ). First, a for passive cancer cells ( τ p = 0) at a given T and P is similar to that for active cancer cells at small τp. [See black/purple circles and black/purple diamonds in Fig. 3(a)]. Next, a decreases with increasing τp for large τp, emphasizing that cancer cell persistence enhances invasion. Thus, a reasonable ansatz for a at fixed β is a 1 / τ d, where τd is the decorrelation time of the velocity autocorrelation function for the cancer cells,
C v v ( t ) = v q ( t ) · v q ( t + t ) v q 2 ( t ) q , t ,
(2)
and
C v v ( τ d ) = e 1 .
(3)
[See Fig. 4(a)]. To understand the dependence of τd on cancer cell packing fraction and τp, we show τd for simplified systems composed of only repulsive, active cancer cells in periodic boundary conditions in Fig. 4(b). τ d τ γ in the τ p 0 limit for all ϕ c, where τ γ = m γ. For intermediate values of τp, τ d τ p, which agrees with the analytical result at ϕ c = 0,24 
C v v ( t ) = τ γ e t τ γ τ p e t τ p τ γ τ p .
(4)
At large τp, τd deviates from the analytical result as ϕ c increases. In the inset to Fig. 3(b), we verify that a τ γ / τ d.
FIG. 3.

(a) Normalized interfacial area An is plotted vs k b T P σ 3 for several values of attraction strength β and persistence time τp of the cancer cells. The symbols indicate the attraction strength: β = 10 4 (triangles), 10 3 (squares), and 10 2 (diamonds). The colors indicate the persistence time: τ p / m σ 2 ε c = 0.02 (purple), 25 (blue), 250 (magenta), and 2500 (orange). We also show data for passive cancer cells with τ p = 0 and β = 10 2 as filled black diamonds. The solid lines are fits to Eq. (1). We include data for both underdamped ( γ = 0.2 m ε c / σ; filled symbols and solid lines) and overdamped ( γ = 20 m ε c / σ; open symbols and dashed lines) dynamics. (b) An is plotted vs the dimensionless energy scale E c b for the same data in (a). The inset shows the relationship between the parameter a in Eq. (1) and c 1 ( 1 + c 2 β ε c P σ 3 ) τ γ τ d. The solid black line indicates a = c 1 ( 1 + c 2 β c P σ 3 ) τ γ τ d.

FIG. 3.

(a) Normalized interfacial area An is plotted vs k b T P σ 3 for several values of attraction strength β and persistence time τp of the cancer cells. The symbols indicate the attraction strength: β = 10 4 (triangles), 10 3 (squares), and 10 2 (diamonds). The colors indicate the persistence time: τ p / m σ 2 ε c = 0.02 (purple), 25 (blue), 250 (magenta), and 2500 (orange). We also show data for passive cancer cells with τ p = 0 and β = 10 2 as filled black diamonds. The solid lines are fits to Eq. (1). We include data for both underdamped ( γ = 0.2 m ε c / σ; filled symbols and solid lines) and overdamped ( γ = 20 m ε c / σ; open symbols and dashed lines) dynamics. (b) An is plotted vs the dimensionless energy scale E c b for the same data in (a). The inset shows the relationship between the parameter a in Eq. (1) and c 1 ( 1 + c 2 β ε c P σ 3 ) τ γ τ d. The solid black line indicates a = c 1 ( 1 + c 2 β c P σ 3 ) τ γ τ d.

Close modal
FIG. 4.

(a) The velocity autocorrelation function C v v ( t ) at different persistence times τp and cancer cell packing fraction ϕ c 0. The filled circles indicate the velocity decorrelation times τd at which C v v ( τ d ) = e 1. The colors from red to blue indicate short to long τp. Open stars represent Eq. (4). (b) The decorrelation time τ d / m σ 2 ε c plotted vs τ p / m σ 2 ε c. The filled symbols indicate τd from simulations of repulsive, active cancer cells without adipocytes at ϕ c = 0 (circles), 0.3 (triangles), 0.4 (squares), and 0.5 (diamonds). The analytical solution for C v v ( τ d ) = e 1 when ϕ c = 0 [Eq. (4)] is shown as a blue solid line.24 The horizontal dashed line indicates τ d = τ γ.

FIG. 4.

(a) The velocity autocorrelation function C v v ( t ) at different persistence times τp and cancer cell packing fraction ϕ c 0. The filled circles indicate the velocity decorrelation times τd at which C v v ( τ d ) = e 1. The colors from red to blue indicate short to long τp. Open stars represent Eq. (4). (b) The decorrelation time τ d / m σ 2 ε c plotted vs τ p / m σ 2 ε c. The filled symbols indicate τd from simulations of repulsive, active cancer cells without adipocytes at ϕ c = 0 (circles), 0.3 (triangles), 0.4 (squares), and 0.5 (diamonds). The analytical solution for C v v ( τ d ) = e 1 when ϕ c = 0 [Eq. (4)] is shown as a blue solid line.24 The horizontal dashed line indicates τ d = τ γ.

Close modal

At fixed τp, An shifts to larger k b T / ( P σ 3 ) with increasing attraction strength β ε c since cohesion makes it more difficult for cancer cells to invade the adipocytes. In the β 0 limit, An converges to a time-dependent value for purely repulsive cancer cells. Thus, we expect a c 1 ( t ) + c β ε c, where c 1 ( t ) k b T g / ( P σ 3 ) in the long time limit and a is not time dependent when c β c 1 ( t ). We show in Fig. 3(b) that a = c 1 ( 1 + c 2 β ε c / ( P σ 3 ) ) τ γ / τ d, where c 1 0.06 enforces A n = 1 / 2 at Ec = 1 and c 2 0.02.

The ECM surrounding the adipocytes provides structural support for adipose tissue. In this section, we investigate the effects of the ECM on breast cancer cell invasion. To model the structural support provided by the ECM, we add linear springs between the centers of mass, R i = 1 N v μ = 1 N v r i μ, 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 dc, such that the number of neighbors per adipocyte is 6. The potential energy associated with the spring connecting adipocytes i and j is
U i j ECM = K ecm 2 ( R i j l i j 0 σ ) 2 ,
(5)
where K ecm / σ 2 is the stiffness of the spring, R i j = | R i R j |, and l i j 0 is the equilibrium length, which is set equal to the adipocyte separations in the initial configuration. All vertices on adipocyte i share the interaction with adipocyte j equally, and thus the force on the vertex μ on adipocyte i from adipocyte j is given by
F i μ j ECM = 1 N v R i U i j ECM .
(6)
In Figs. 5(a) and 5(b), we calculate the interfacial area An vs k b T / ( P σ 3 ) for several values of Kecm and two values of τp and β. As discussed above, An vs k b T / ( P σ 3 ) is normalized between 0 and 1 using the interfacial area obtained for untethered (Kecm = 0) systems. Tethering of adipocytes decreases the degree of invasion but not the onset value of k b T / ( P σ 3 ) for invasion. We find that the plateau values A n max < 1 at large k b T / ( P σ 3 ) for all Kecm > 0, but the values of k b T / ( P σ 3 ) at which A n > 0 does not vary strongly with Kecm. In the insets to Fig. 5, we show that A n max decreases monotonically from 1 to lower values that depend on τp and β. Thus, the tethering of adipocytes by the ECM reduces the degree of cancer cell invasion into adipose tissue.
FIG. 5.

An is plotted vs the cancer cell motility k b T P σ 3 for persistence times and attractive strengths: (a) τ p = 25 m σ 2 ε c and β = 10 2, and (b) τ p = 2500 m σ 2 ε c and β = 10 3. Stiffness of the ECM increases from dark to light green. Systems with Kecm = 0 are indicated by black filled circles. The dashed horizontal lines indicate A n max for different values of Kecm. Insets: A n max is plotted as a function of K ecm / ε c for the data in the main panels. The black dashed horizontal lines indicate A n max in the K ecm / ε c limit.

FIG. 5.

An is plotted vs the cancer cell motility k b T P σ 3 for persistence times and attractive strengths: (a) τ p = 25 m σ 2 ε c and β = 10 2, and (b) τ p = 2500 m σ 2 ε c and β = 10 3. Stiffness of the ECM increases from dark to light green. Systems with Kecm = 0 are indicated by black filled circles. The dashed horizontal lines indicate A n max for different values of Kecm. Insets: A n max is plotted as a function of K ecm / ε c for the data in the main panels. The black dashed horizontal lines indicate A n max in the K ecm / ε c limit.

Close modal

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 ϕ = 0.72 as described in the Methods section. We now vary the adipocyte packing fraction over the range 0.4 < ϕ 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 ε c, and calculate An vs ϕ a for τ p / m σ 2 ε c = 25 and β = 10 2. In Fig. 6(a), we confirm that An decreases with increasing adipocyte packing fraction and vanishes for ϕ a 0.85.

FIG. 6.

(a) Standard deviation in the local packing fraction of adipocytes Δϕa plotted vs ϕa. Open dark green circles indicate results for uniformly distributed adipocytes at a total packing fraction ϕ = 0.72. (b) The degree of invasion An plotted vs the packing fraction of adipocytes ϕ a at τ P = 25 m σ 2 ε c , β = 10 2, and k b T / ( P σ 3 ) = 0.37. λ, which determines the packing fraction heterogeneity of adipose tissue, increases from dark to light green.

FIG. 6.

(a) Standard deviation in the local packing fraction of adipocytes Δϕa plotted vs ϕa. Open dark green circles indicate results for uniformly distributed adipocytes at a total packing fraction ϕ = 0.72. (b) The degree of invasion An plotted vs the packing fraction of adipocytes ϕ a at τ P = 25 m σ 2 ε c , β = 10 2, and k b T / ( P σ 3 ) = 0.37. λ, which determines the packing fraction heterogeneity of adipose tissue, increases from dark to light green.

Close modal

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 λ ) R i j and to l i j 0 = ( 1 + λ ) 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 Δ ϕ are 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.

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 At shared between the cancer cells and adipocytes. At 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 An 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 An can be collapsed when plotted as a function of a dimensionless energy scale Ec, where E c 1 indicates a de-mixed (non-invaded) state and E c 1 indicates a mixed (invaded) state. We found that Ec 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 An relative to more uniform local packing.

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 kbT and persistence time τp, 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 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 τ inv is shorter than the time scale of the cell cycle τ 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 τ 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 τ cell τ 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 τ cell / τ inv.

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.

To model dense packings of adipocytes, which form faceted cell–cell contacts under compression, we employ the recently developed 3D deformable particle model. Specifically, we consider adipocytes as deformable polyhedra composed of Nv = 42 vertices belonging to a mesh with Nf = 80 triangular faces, which is sufficient to describe the shape fluctuations observed in adipocytes. [See Fig. 7(f)]. The shape energy for the ith adipocyte as a function of the positions of its vertices r i 1 , , r i N v is
U i = ε v 2 ( 1 v i v i 0 ) 2 + ε a 2 k = 1 N f ( 1 a i k a i k 0 ) 2 .
(7)
FIG. 7.

(a) and (b) Histological slices of the interface between a mammary tumor and adipose tissue in mice stained with Masson's trichrome. Cancer cell nuclei are stained dark purple, adipocytes are the large unstained white regions, and collagen fibers are stained blue. (a) Non-invaded adipose tissue, characterized by a smooth interface between the tumor and adipose tissue. (b) Invaded adipose tissue with a rough interface. (c) and (d) 2D slices in the xy-plane from 3D DPM simulations in (e), both (c) before and (d) after invasion. Red circles represent cancer cells and adipocytes are shown as smooth, cyan-colored deformable polygons. Dashed lines indicate periodic boundaries in the y-direction. (e) Initial configuration of de-mixed cancer cells (red spheres) and adipocytes (smooth, multi-colored deformable polyhedra) in simulations with Na = 128 and Nc = 7000. (f) Adipocytes are modeled as deformable polyhedra with vertices of diameter σ i μ and a triangular mesh with individual triangle area aik. Cancer cells are modeled as soft spheres with diameter σq.

FIG. 7.

(a) and (b) Histological slices of the interface between a mammary tumor and adipose tissue in mice stained with Masson's trichrome. Cancer cell nuclei are stained dark purple, adipocytes are the large unstained white regions, and collagen fibers are stained blue. (a) Non-invaded adipose tissue, characterized by a smooth interface between the tumor and adipose tissue. (b) Invaded adipose tissue with a rough interface. (c) and (d) 2D slices in the xy-plane from 3D DPM simulations in (e), both (c) before and (d) after invasion. Red circles represent cancer cells and adipocytes are shown as smooth, cyan-colored deformable polygons. Dashed lines indicate periodic boundaries in the y-direction. (e) Initial configuration of de-mixed cancer cells (red spheres) and adipocytes (smooth, multi-colored deformable polyhedra) in simulations with Na = 128 and Nc = 7000. (f) Adipocytes are modeled as deformable polyhedra with vertices of diameter σ i μ and a triangular mesh with individual triangle area aik. Cancer cells are modeled as soft spheres with diameter σq.

Close modal

The first term imposes a harmonic energy penalty for changes in cell volume vi from its preferred value v i 0 and εv 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 aik of each triangular face from its preferred value a i k 0, and εa controls the fluctuations in adipocyte surface area. We quantify the adipocyte shape using the dimensionless shape parameter A = ( k = 1 N f a i k ) 3 2 6 π 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.

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 Nv = 42 vertices and Nf = 80 faces have Δ N = ( 3 N v 6 ) ( 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.18 In addition, using the DPM allows packings of adipocytes to be nearly confluent without overlaps between adipocytes.

The interactions between adipocytes are controlled by the pair forces between their vertices. In particular, the repulsive pair force between vertex μ on adipocyte i and vertex ν on adipocyte j is purely repulsive,
F ij μ ν = { ε c σ ij μ ν ( 1 r ij μ ν σ ij μ ν ) r ̂ ij μ ν , r ij μ ν σ ij μ ν , 0 , r ij μ ν > σ ij μ ν ,
(8)
where εc controls the strength of the repulsive forces, r ̂ ij μ ν = r ij μ ν / r ij μ ν , r ij μ ν = r i μ r j ν is the separation vector, r ij μ ν = | r ij μ ν | is the distance between vertices, and σ ij μ ν = 1 2 ( σ i μ + σ j ν ) is their average diameter. (See Fig. 8 in  Appendix B). σ i μ 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 N v N a.
FIG. 8.

(a) Distance r ij μ v between vertex μ on adipocyte i and vertex ν on adipocyte j. The average diameter of the vertices is σ ij μ ν = ( σ i μ + σ j ν ) / 2. (b) Distance r i μ q between vertex μ on adipocyte i (with diameter σ i μ) and cancer cell q (with diameter σq). (c) Distance rqs between cancer cells q and s.

FIG. 8.

(a) Distance r ij μ v between vertex μ on adipocyte i and vertex ν on adipocyte j. The average diameter of the vertices is σ ij μ ν = ( σ i μ + σ j ν ) / 2. (b) Distance r i μ q between vertex μ on adipocyte i (with diameter σ i μ) and cancer cell q (with diameter σq). (c) Distance rqs between cancer cells q and s.

Close modal
The diameter of adipocytes in human breast adipose tissue ranges from 50 to 150 μ m,33 while cancer cells are much smaller, ranging between 10 and 20 μ 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,
F i μ q = { ε c σ i μ q ( 1 r i μ q σ i μ q ) r ̂ iuq , r i μ q σ i μ q , 0 , r i μ q > σ i μ q ,
(9)
where r ̂ i μ q = r i μ q / r i μ q , r i μ q = r i μ r q is the separation vector, r i μ q = | r i μ q | is the distance between vertex μ on adipocyte i and cancer cell q, σ i μ q = 1 2 ( σ i μ + σ q ), and σq is the diameter of the qth cancer cell.
For the interactions between cancer cells, we include contact repulsions as well as short-range attractions. We assume that the pair force on cancer cell q from cancer cell s obeys
F q s = { ε c σ q s ( 1 r q s σ q s ) r ̂ q s , r q s r β , ε c σ q s r β σ q s r α r β r q s r α σ q s r ̂ q s , r β < r q s r α , 0 , r q s > r α ,
(10)
where r ̂ q s = r q s / r q s , r q s = r q r s is the separation vector, r q s = | r q s | is the distance between cancer cells q and s, and σ q s = 1 2 ( σ q + σ 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 < σ q s, and attractive for 1 < r q s / σ q s < r α / σ q s with maximum attractive force F q s = ε c β / σ q s at r q s = r β, where r α / σ q s = 1 + α and r β / σ q s = 1 + β. We vary the depth of the attractive potential β from 10 4 to 10 2 and set the range of attractive interactions at α = 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 σ and the other half have a larger diameter 1.1 σ, 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
FIG. 9.

The pairwise interaction force F ̃ between two cells determines cell–cell adhesion (attractive) and volume exclusion (repulsive) and varies with the separation r ̃ between the cells. For adipocyte–adipocyte interactions, F ̃ = F ij μ ν , r ̃ = r ij μ ν, and σ ̃ = σ ij μ ν, where i , j = 1 , , N a are the adipocyte indices and μ , ν = 1 , , N v are the vertex indices. Na and Nv are the number of adipocytes in the simulation box and the number of vertices on an adipocyte, respectively. For adipocyte–cancer cell interactions, F ̃ = F i μ q , r ̃ = r i μ q, and σ ̃ = σ i μ q, where q = 1 , , N c. Nc is the number of cancer cells in the simulation box. Both adipocytes and cancer cells interact via repulsive linear spring forces when they overlap (black solid line) and do not interact when they are out of contact (blue dashed and black dash-dotted lines). For the cancer cell–cancer cell interactions, F ̃ = F q s , r ̃ = r q s, and σ ̃ = σ q s. These interactions include linear repulsion when two cancer cells overlap (black solid line), short-range attraction for 1 < r ̃ / σ ̃ < r α / σ ̃ (red dotted line), and no interactions for r ̃ / σ ̃ > r α / σ ̃ (black dash-dotted line).

FIG. 9.

The pairwise interaction force F ̃ between two cells determines cell–cell adhesion (attractive) and volume exclusion (repulsive) and varies with the separation r ̃ between the cells. For adipocyte–adipocyte interactions, F ̃ = F ij μ ν , r ̃ = r ij μ ν, and σ ̃ = σ ij μ ν, where i , j = 1 , , N a are the adipocyte indices and μ , ν = 1 , , N v are the vertex indices. Na and Nv are the number of adipocytes in the simulation box and the number of vertices on an adipocyte, respectively. For adipocyte–cancer cell interactions, F ̃ = F i μ q , r ̃ = r i μ q, and σ ̃ = σ i μ q, where q = 1 , , N c. Nc is the number of cancer cells in the simulation box. Both adipocytes and cancer cells interact via repulsive linear spring forces when they overlap (black solid line) and do not interact when they are out of contact (blue dashed and black dash-dotted lines). For the cancer cell–cancer cell interactions, F ̃ = F q s , r ̃ = r q s, and σ ̃ = σ q s. These interactions include linear repulsion when two cancer cells overlap (black solid line), short-range attraction for 1 < r ̃ / σ ̃ < r α / σ ̃ (red dotted line), and no interactions for r ̃ / σ ̃ > r α / σ ̃ (black dash-dotted line).

Close modal
To determine the spatiotemporal trajectories of adipocytes and cancer cells, we integrate Newton's equations of motion. For the cancer cells, we consider both active and passive models. For the active model, the equation of motion for the position rq of cancer cell q is
m d 2 r q d t 2 = s q = 1 N c F q s i = 1 N a μ = 1 N v F i μ q w = l r F w q γ v q + f 0 n ̂ q ,
(11)
where m is the mass of the cancer cells, F i μ 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 qth 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 vq with damping coefficient γ σ / m ε c = 0.2 and a self-propulsion term with active force magnitude f0. The direction of the active force n ̂ q obeys rotational diffusion,35–37,
d n ̂ q d t = χ q × n ̂ q ,
(12)
where χ q is a vector with x-, y-, and z-components that are Gaussian random numbers with zero mean and satisfy
χ q ( t ) χ s ( t + t ) t = τ p 1 I δ q s δ ( t ) ,
(13)
where . t is an average over time origins, I is the identity matrix, δ q s = 1 when q = s and 0 otherwise, and δ ( . ) is the Dirac δ-function. τ p is the persistence time of the active force as defined through the autocorrelation function of n ̂ q,
n ̂ q ( t ) · n ̂ q ( t + t ) q , t = exp ( t τ p ) ,
(14)
where . q , t is an average over cancer cells and time origins. Persistence in cancer cell migration originates from the fact that the migration direction and speed are strongly correlated with the local alignment of ECM fibers.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 kbT that captures the random motion of the cancer cells and persistence time τp that captures the correlation between cancer cell velocity and local ECM fiber alignment. Both kbT and τp describe the migration properties of cancer cells, and implicitly model the degree of alignment of the ECM.
We also consider a passive model for cancer cell migration to better understand the role of persistent cancer cell motion on breast cancer invasion. For the passive model, we employ the Langevin equation of motion for the qth cancer cell,
m d 2 r q d t 2 = s q = 1 N c F q s i = 1 N a μ = 1 N v F i μ q w = l r F w q γ v q + 2 k b T 0 γ m η q ( t ) ,
(15)
where the first and second terms include interactions among cancer cells and between cancer cells and adipocytes. As in Eq. (11), the equation of motion for passive cancer cells includes a damping force proportional to the velocity vq of cancer cell q with damping coefficient γ. The last term is the noise term with strength 2 k b T 0 γ m that sets the target temperature T0. η ( t ) is a vector with x-, y-, and z-components that are Gaussian random numbers with zero mean and unit variance.
We employ damped equations of motion for the adipocytes, i.e., the position of vertex μ on the ith adipocyte obeys
m d 2 r i μ d t 2 = j i = 1 N a ν = 1 N v F ij μ ν + q = 1 N c F i μ q w = l r F w i μ γ v i μ r i μ U i ,
(16)
where 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 μ is the repulsive force between the physical walls and vertex μ on the ith adipocyte, which is defined in Eq. (19). The fifth term gives the mechanical forces generated from the shape-energy function for the ith adipocyte in Eq. (7).

We integrate the equations of motion for the adipocyte vertices and cancer cells [i.e., Eqs. (11), (15), and (16)] using the modified velocity Verlet integration scheme with time step Δ t = 5 × 10 2 m σ 2 ε c.

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 σ, Na = 28 adipocytes, Nc = 1500 cancer cells, and initial total packing fraction ϕ = 0.01. (Note the total packing fraction can be written in terms of the packing fractions for cancer cells and adipocytes, separately, ϕ = ϕ c + ϕ a.) The cancer cells and adipocytes are placed randomly in the simulation box such that the x-components of the cancer cell centers of mass satisfy 2 σ < x q < 12 σ and the x-components of the adipocyte centers of mass satisfy 16 σ < X i < 68 σ. We then perform athermal, isotropic compression, applying successive compression steps Δ ϕ / ϕ = 0.03 to generate a jammed packing of adipocytes. Compression is followed by energy minimization until the final packing fraction of the system reaches ϕ = ( π 6 q = 1 N c σ q 3 + i = 1 N a V i ) / L f 3 = 0.72, where L f = 34 σ is the final side length of the box, and Vi is the total volume of the ith adipocyte, i.e., vi plus the additional volumes from the vertices that do not overlap the polyhedron. [See Fig. 7(e)].

During compression, we apply periodic boundary conditions in the y- and z-directions. In the x-direction, we include two physical walls, a “right” wall at xr = Lf and a “left” wall at xl = 0. The right wall is stationary, while the left wall moves in the x-direction to maintain fixed pressure P. The equation of the motion for the left wall position is
m d 2 x l d t 2 = q = 1 N c F l q · x ̂ + i = 1 N a μ = 1 N v F l i μ · x ̂ + P L f 2 γ v l ,
(17)
where m is the mass and γ is the damping coefficient for the left wall, vl is the speed of the left wall. The forces that act on the wall from the qth cancer cell and μth vertex on the ith adipocyte are
F w q = { ε c σ q ( 1 x w q σ q ) x ̂ , x w q σ q , 0 , x w q > σ q ,
(18)
and
F w i μ = { ε c σ i μ ( 1 x w i μ σ i μ ) x ̂ , x w i μ σ i μ , 0 , x w i μ > σ i μ ,
(19)
where x w q = x w x q and x w i μ = x w x i μ. (Note that the subscript 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,
k b T = 2 m 3 N c q = 1 N c v q 2 .
(20)

For the passive model of cancer cells, the temperature of the cancer cells Tc and adipocytes Ta 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.

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.

The authors have no conflicts to disclose.

Ethics approval is not required.

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).

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

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 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 π 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 1.1 gives reasonable agreement for P ( A 2 D ) between the experiments and simulations.

FIG. 10.

(a) A representative histology slice of mammary adipose tissue from a C57BL/6 mouse. (b) Cell boundaries are identified through standard image processing techniques. Cells are colored by A 2 D = p 2 / ( 4 π s ), where p is the perimeter and s is the area. (c) Probability distribution of A 2 D , P ( A 2 D ), for random 2D cross sections of adipocytes in 3D DEM invasion simulations with reference shape parameters A 0 = 1.0 (dashed line), 1.1 (dotted line), and 1.3 (dot-dashed line) and from histological slices (solid line). The experimental data include 1000 adipocytes collected from histology samples from 5 mice.

FIG. 10.

(a) A representative histology slice of mammary adipose tissue from a C57BL/6 mouse. (b) Cell boundaries are identified through standard image processing techniques. Cells are colored by A 2 D = p 2 / ( 4 π s ), where p is the perimeter and s is the area. (c) Probability distribution of A 2 D , P ( A 2 D ), for random 2D cross sections of adipocytes in 3D DEM invasion simulations with reference shape parameters A 0 = 1.0 (dashed line), 1.1 (dotted line), and 1.3 (dot-dashed line) and from histological slices (solid line). The experimental data include 1000 adipocytes collected from histology samples from 5 mice.

Close modal

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 μ ν , r i μ q, and rqs and the diameters σ i μ and σq, which are pictured in Fig. 8.

In this  Appendix, we describe how we calculate the interfacial area At 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) At is defined as the sum of the areas of these faces. The maximum value of At is A t max = 2 L f 2 + Σ i = 1 N a S i, where Si is the surface area of the Voronoi polyhedron for adipocyte i. The relative interfacial area A t / A t max serves as an indicator of the degree of cancer cell invasion.

FIG. 11.

(a)–(c) Voronoi tessellations of the centers of cancer cells (red shaded Voronoi polyhedra) and adipocytes (blue shaded Voronoi polyhedra) during cancer cell invasion into adipose tissue at several values of shared interfacial area A t / σ 2: (a) 610, (b) 1280, and (c) 1720. (d)–(f) The Voronoi polyhedral faces that are shared between the cancer cells and adipocytes for panels (a)–(c). To improve visualization, we do not include the faces of cancer cell Voronoi polyhedra that form the physical boundaries in the x-direction.

FIG. 11.

(a)–(c) Voronoi tessellations of the centers of cancer cells (red shaded Voronoi polyhedra) and adipocytes (blue shaded Voronoi polyhedra) during cancer cell invasion into adipose tissue at several values of shared interfacial area A t / σ 2: (a) 610, (b) 1280, and (c) 1720. (d)–(f) The Voronoi polyhedral faces that are shared between the cancer cells and adipocytes for panels (a)–(c). To improve visualization, we do not include the faces of cancer cell Voronoi polyhedra that form the physical boundaries in the x-direction.

Close modal

Note that the volume of the simulation box, as well as At and A t max, can increase with k b T / ( P σ 3 ), especially at large values when k b T overcomes the cohesive energy β ε c. In this regime, the ratio A t / A t max can vary with k b T / ( P σ 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 Δ A ) / ( A t max Δ A ), where Δ A = A t max 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 min ( r A t ) ) / ( max ( r A t ) min ( r A t ) ), where min ( r A t ) and max ( r A t ) are the minimum and maximum values of rAt, such that 0 < A n < 1.

1.
J.
Majidpoor
and
K.
Mortezaee
, “
Steps in metastasis: An updated review
,”
Med. Oncol.
38
,
3
(
2021
).
2.
N. M.
Novikov
,
S. Y.
Zolotaryova
,
A. M.
Gautreau
, and
E. V.
Denisov
, “
Mutational drivers of cancer cell migration and invasion
,”
Br. J. Cancer
124
,
102
114
(
2021
).
3.
O.
Ilina
,
L.
Campanello
,
P. G.
Gritsenko
,
M.
Vullings
,
C.
Wang
,
P.
Bult
,
W.
Losert
, and
P.
Friedl
, “
Intravital microscopy of collective invasion plasticity in breast cancer
,”
Dis. Models Mech.
11
,
dmm034330
(
2018
).
4.
Q.
Zhu
,
Y.
Zhu
,
C.
Hepler
,
Q.
Zhang
,
J.
Park
,
C.
Gliniak
,
G. H.
Henry
,
C.
Crewe
,
D.
Bu
,
Z.
Zhang
,
S.
Zhao
,
T.
Morley
,
N.
Li
,
D. S.
Kim
,
D.
Strand
,
Y.
Deng
,
J. J.
Robino
,
O.
Varlamov
,
R.
Gordillo
,
M. G.
Kolonin
,
C. M.
Kusminski
,
R. K.
Gupta
, and
P. E.
Scherer
, “
Adipocyte mesenchymal transition contributes to mammary tumor progression
,”
Cell Rep.
40
,
111362
(
2022
).
5.
E.
Milotti
,
S.
Stella
, and
R.
Chignola
, “
Pulsation-limited oxygen diffusion in the tumour microenvironment
,”
Sci. Rep.
7
,
39762
(
2017
).
6.
H. H.
Marti
,
Angiogenesis—A Self-Adapting Principle In Hypoxia
(
Springer
,
Basel
,
2005
) pp.
163
180
.
7.
F.
Sauer
,
S.
Grosser
,
M.
Shahryari
,
A.
Hayn
,
J.
Guo
,
J.
Braun
,
S.
Briest
,
B.
Wolf
,
B.
Aktas
,
L.
Horn
,
I.
Sack
, and
J. A.
Käs
, “
Changes in tissue fluidity predict tumor aggressiveness in vivo
,”
Adv. Sci.
10
,
e2303523
(
2023
).
8.
M.
Zhang
,
J. S.
Di Martino
,
R. L.
Bowman
,
N. R.
Campbell
,
S. C.
Baksh
,
T.
Simon-Vermot
,
I. S.
Kim
,
P.
Haldeman
,
C.
Mondal
,
V.
Yong-Gonzales
,
M.
Abu-Akeel
,
R. T.
Merghoub
,
D.
Jones
,
X. G.
Zhu
,
A.
Arora
,
C. E.
Ariyan
,
K.
Birsoy
,
J. D.
Wolchok
,
K. S.
Panageas
,
T.
Hollmann
,
J. J.
Bravo-Cordero
, and
R. M.
White
, “
Adipocyte-derived lipids mediate melanoma progression via FATP proteins
,”
Cancer Discovery
8
,
1006
1025
(
2018
).
9.
P.
Friedl
,
J.
Locker
,
E.
Sahai
, and
J. E.
Segall
, “
Classifying collective cancer cell invasion
,”
Nat. Cell Biol.
14
,
777
783
(
2012
).
10.
M.
Janiszewska
,
M. C.
Primi
, and
T.
Izard
, “
Cell adhesion in cancer: Beyond the migration of single cells
,”
J. Biol. Chem.
295
,
2495
2505
(
2020
).
11.
R.
Li
,
J. D.
Hebert
,
T. A.
Lee
,
H.
Xing
,
A.
Boussommier-Calleja
,
R. O.
Hynes
,
D. A.
Lauffenburger
, and
R. D.
Kamm
, “
Macrophage-secreted TNFα and TGFβ1 influence migration speed and persistence of cancer cells in 3D tissue culture via independent pathways
,”
Cancer Res.
77
,
279
290
(
2017
).
12.
N. V.
Krakhmal
,
M. V.
Zavyalova
,
E. V.
Denisov
,
S. V.
Vtorushin
, and
V. M.
Perelmuter
, “
Cancer invasion: Patterns and mechanisms
,”
Acta Nat.
7
,
17
28
(
2015
).
13.
I. G.
Gonçalves
and
J. M.
Garcia-Aznar
, “
Extracellular matrix density regulates the formation of tumour spheroids through cell migration
,”
PLOS Comput. Biol.
17
,
e1008764
(
2021
).
14.
B. R.
Seo
,
P.
Bhardwaj
,
S.
Choi
,
J.
Gonzalez
,
R. C. A.
Eguiluz
,
K.
Wang
,
S.
Mohanan
,
P. G.
Morris
,
B.
Du
,
X. K.
Zhou
,
L. T.
Vahdat
,
A.
Verma
,
O.
Elemento
,
C. A.
Hudis
,
R. M.
Williams
,
D.
Gourdon
,
A. J.
Dannenberg
, and
C.
Fischbach
, “
Obesity-dependent changes in interstitial ECM mechanics promote breast tumorigenesis
,”
Sci. Transl. Med.
7
,
301ra130
(
2015
).
15.
A. A.
Shimpi
,
E. D.
Williams
,
L.
Ling
,
T.
Tamir
,
F. M.
White
, and
C.
Fischbach
, “
Phosphoproteomic changes induced by cell-derived matrix and their effect on tumor cell migration and cytoskeleton remodeling
,”
ACS Biomater. Sci. Eng.
9
,
6835
6848
(
2023
).
16.
K.
Wolf
,
M.
Te Lindert
,
M.
Krause
,
S.
Alexander
,
J.
Te Riet
,
A. L.
Willis
,
R. M.
Hoffman
,
C. G.
Figdor
,
S. J.
Weiss
, and
P.
Friedl
, “
Physical limits of cell migration: Control by ECM space and nuclear deformation and tuning by proteolysis and traction force
,”
J. Cell Biol.
201
,
1069
1084
(
2013
).
17.
A. T.
Ton
,
A. K.
MacKeith
,
M. D.
Shattuck
, and
C. S.
O'Hern
, “
Mechanical plasticity of cell membranes enhances epithelial wound closure
,”
Phys. Rev. Res.
6
,
L012036
(
2024
).
18.
D.
Wang
,
J. D.
Treado
,
A.
Boromand
,
B.
Norwick
,
M. P.
Murrell
,
M. D.
Shattuck
, and
C. S.
O'Hern
, “
The structural, vibrational, and mechanical properties of jammed packings of deformable particles in three dimensions
,”
Soft Matter
17
,
9901
9915
(
2021
).
19.
Y.
Cheng
,
J. D.
Treado
,
B. F.
Lonial
,
P.
Habdas
,
E. R.
Weeks
,
M. D.
Shattuck
, and
C. S.
O'Hern
, “
Hopper flows of deformable particles
,”
Soft Matter
18
,
8071
8086
(
2022
).
20.
A.
Boromand
,
A.
Signoriello
,
F.
Ye
,
C. S.
O'Hern
, and
M. D.
Shattuck
, “
Jamming of deformable polygons
,”
Phys. Rev. Lett.
121
,
248003
(
2018
).
21.
J. D.
Treado
,
A. B.
Roddy
,
G.
Théroux-Rancourt
,
L.
Zhang
,
C.
Ambrose
,
C. R.
Brodersen
,
M. D.
Shattuck
, and
C. S.
O'Hern
, “
Localized growth and remodelling drives spongy mesophyll morphogenesis
,”
J. R. Soc. Interface
19
,
20220602
(
2022
).
22.
J. D.
Treado
,
D.
Wang
,
A.
Boromand
,
M. P.
Murrell
,
M. D.
Shattuck
, and
C. S.
O'Hern
, “
Bridging particle deformability and collective response in soft solids
,”
Phys. Rev. Mater.
5
,
055605
(
2021
).
23.
M. L.
Manning
, “
Essay: Collections of deformable particles present exciting challenges for soft matter and biological physics
,”
Phys. Rev. Lett.
130
,
130002
(
2023
).
24.
L.
Caprini
and
U. M. B.
Marconi
, “
Inertial self-propelled particles
,”
J. Chem. Phys.
154
,
024902
(
2021
).
25.
M.
Alsaggar
,
S.
Bdour
,
Q.
Ababneh
,
T.
El-Elimat
,
N.
Qinna
, and
K. H.
Alzoubi
, “
Silibinin attenuates adipose tissue inflammation and reverses obesity and its complications in diet-induced obesity model in mice
,”
BMC Pharmacol. Toxicol.
21
,
1
8
(
2020
).
26.
O.
Ilina
,
P. G.
Gritsenko
,
S.
Syga
,
J.
Lippoldt
,
C. A. L.
Porta
,
O.
Chepizhko
,
S.
Grosser
,
M.
Vullings
,
G. J.
Bakker
,
J.
Starruß
,
P.
Bult
,
S.
Zapperi
,
J. A.
Käs
,
A.
Deutsch
, and
P.
Friedl
, “
Cell–cell adhesion and 3D matrix confinement determine jamming transitions in breast cancer invasion
,”
Nat. Cell Biol.
22
,
1103
1115
(
2020
).
27.
N.
Alkhouli
,
J.
Mansfield
,
E.
Green
,
J.
Bell
,
B.
Knight
,
N.
Liversedge
,
J. C.
Tham
,
R.
Welbourn
,
A. C.
Shore
,
K.
Kos
, and
C. P.
Winlove
, “
The mechanical properties of human adipose tissues and their relationships to the structure and composition of the extracellular matrix
,”
Am. J. Physiol.
305
,
E1427
E1435
(
2013
).
28.
W. Y.
Wang
,
A. T.
Pearson
,
M. L.
Kutys
,
C. K.
Choi
,
M. A.
Wozniak
,
B. M.
Baker
, and
C. S.
Chen
, “
Extracellular matrix alignment dictates the organization of focal adhesions and directs uniaxial cell migration
,”
APL Bioeng.
2
,
046107
(
2018
).
29.
D.
Sun
,
Y.
Liu
,
H.
Wang
,
F.
Deng
,
Y.
Zhang
,
S.
Zhao
,
X.
Ma
,
H.
Wu
, and
G.
Sun
, “
Novel decellularized liver matrix-alginate hybrid gel beads for the 3D culture of hepatocellular carcinoma cells
,”
Int. J. Biol. Macromol.
109
,
1154
1163
(
2018
).
30.
M. S. U.
Rahman
,
J.
Wu
,
H.
Chen
,
C.
Sun
,
Y.
Liu
, and
S.
Xu
, “
Matrix mechanophysical factor: Pore size governs the cell behavior in cancer
,”
Adv. Phys.
8
,
2153624
(
2023
).
31.
S.
Grosser
,
J.
Lippoldt
,
L.
Oswald
,
M.
Merkel
,
D. M.
Sussman
,
F.
Renner
,
P.
Gottheil
,
E. W.
Morawetz
,
T.
Fuhs
,
X.
Xie
,
S.
Pawlizak
,
A. W.
Fritsch
,
B.
Wolf
,
L.-C.
Horn
,
S.
Briest
,
B.
Aktas
,
M. L.
Manning
, and
J. A.
Käs
, “
Cell and nucleus shape as an indicator of tissue fluidity in carcinoma
,”
Phys. Rev. X
11
,
011033
(
2021
).
32.
Y.
Ye
and
J.
Lin
, “
Fingering instability accelerates population growth of a proliferating cell collective
,”
Phys. Rev. Lett.
132
,
018402
(
2024
).
33.
S. D.
Parlee
,
S. I.
Lentz
,
H.
Mori
, and
O. A.
MacDougald
, “
Quantifying size and number of adipocytes in adipose tissue
,” in
Methods in Enzymology
(
Elsevier
,
Cambridge
,
2014
), Vol.
537
, pp.
93
122
.
34.
T. N.
Truongvo
,
R. M.
Kennedy
,
H.
Chen
,
A.
Chen
,
A.
Berndt
,
M.
Agarwal
,
L.
Zhu
,
H.
Nakshatri
,
J.
Wallace
,
S.
Na
,
H.
Yokota
, and
J. E.
Ryu
, “
Microfluidic channel for characterizing normal and breast cancer cells
,”
J. Micromech. Microeng.
27
,
035017
(
2017
).
35.
V. E.
Debets
,
X. M. D.
Wit
, and
L. M.
Janssen
, “
Cage length controls the nonmonotonic dynamics of active glassy matter
,”
Phys. Rev. Lett.
127
,
278002
(
2021
).
36.
M. R.
Shaebani
,
A.
Wysocki
,
R. G.
Winkler
,
G.
Gompper
, and
H.
Rieger
, “
Computational models for active matter
,”
Nat. Rev. Phys.
2
,
181
199
(
2020
).
37.
S. C.
Takatori
and
J. F.
Brady
, “
Towards a thermodynamics of active matter
,”
Phys. Rev. E
91
,
032117
(
2015
).
38.
N.
Gjorevski
,
A. S.
Piotrowski
,
V. D.
Varner
, and
C. M.
Nelson
, “
Dynamic tensile forces drive collective cell migration through three-dimensional extracellular matrices
,”
Sci. Rep.
5
,
11458
(
2015
).
39.
K. M.
Wisdom
,
K.
Adebowale
,
J.
Chang
,
J. Y.
Lee
,
S.
Nam
,
R.
Desai
,
N. S.
Rossen
,
M.
Rafat
,
R. B.
West
,
L.
Hodgson
et al, “
Matrix mechanical plasticity regulates cancer cell migration through confining microenvironments
,”
Nat. Commun.
9
,
4144
(
2018
).