Computational modeling of the physical features that influence breast cancer invasion into adipose tissue

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 Ec>1, cancer cells will invade the adipose tissue, whereas for Ec<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.


INTRODUCTION
Breast cancer cell invasion into mammary adipose tissue is a key step toward disease progression and metastasis. 1,2aining 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.4][5] 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. 4Oscillations in the oxygen concentration in the TME caused by collective diffusion affect tumor growth and invasion. 3Moreover, genetic and chemical cues can trigger the epithelial-mesenchymal transition (EMT), which decreases the cohesion between cancer cells and increases their motility. 5In the context of adipose tissue, interactions between cancer cells and adipocytes cause lipid loss and dedifferentiation in adipocytes, which promote invasion by altering cancer cell metabolism and reconfiguring the stroma. 6,7][10][11][12][13] In adipose tissue, ECM fiber alignment is influenced by adipocyte properties and regulates tumor cell invasion by controlling migration persistence. 13,14re persistently moving cancer cells, in turn, can invade adipose tissue more rapidly.Cancer cell cohesion also influences tumor growth and invasion. 8,9Highly cohesive cancer cells do not disconnect from each other and thus do not as easily invade out of the primary tumor.With decreasing cohesion, cancer cells can more readily detach from the tumor and invade distant tissues. 8,11During invasion, cancer cells navigate between adipocytes which form dense polyhedral packings.Thus, the packing fraction and shape deformation of adipocytes, as well as the porous ECM network between adipocytes also physically influence breast cancer invasion. 13,15n 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) [16][17][18][19][20][21][22] and carry out DPM simulations of self-propelled, cohesive spherical particles (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.The systems are de-mixed when initialized, with a smooth interface between the cancer cells and adipose tissue.During 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 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 c 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 Methods section, we outline the DPM shape-energy function for adipocytes and the soft particle model for the cancer cells, the interac-tions 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.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 c .We then show that A t 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 E c for systems that mimic the physical properties of adipose tissue.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 A t between cancer cells and adipocytes is normalized to remove the effects of dilation at large E c .FIG. 2. The pairwise interaction force F plotted as a function of the separation r between different cell types.For adipocyte-adipocyte interactions, F = F i jµν , r = r i jµν , and σ = σ i jµν , where i, j = 1, . . ., N a are the adipocyte indices and µ, ν = 1, . . ., N v are the vertex indices.N a and N v are the number of adipocytes in the simulation box and the number of vertices on an adipocyte, respectively.For adipocytecancer cell interactions, F = F iµq , r = r iµq , and σ = σ iµq , where q = 1, . . ., N c .N c 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 qs , r = r qs , and σ = σ qs .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 dashdotted line).

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, which is used to model dense packings of 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 calculate 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
To model dense packings of adipocytes, which form faceted cell-cell contacts under compression, we leverage the recently developed 3D deformable particle model.Specifically, we consider adipocytes as deformable polyhedra composed of N v = 42 vertices belonging to a mesh with N f = 80 triangular faces, which is sufficient to describe the shape fluctuations observed in adipocytes.(See Fig. 1 (f).)The shape energy for the ith adipocyte as a function of the positions of its vertices r i1 , . . ., r iN v is: The first term imposes a harmonic energy penalty for changes in cell volume v i from its preferred value v i0 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 a ik of each triangular face from its preferred value a ik0 and ε a controls the fluctuations in the adipocyte surface area.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 f = 80 faces have 39 zero-energy modes of the dynamical matrix for the shape-energy function in Eq. 1. 17 ) We quantify adipocyte shape using the dimensionless shape parameter . 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.

Interactions between adipocytes, between cancer cells, and between adipocytes and cancer cells
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: where ε c controls the strength of the repulsive forces, ri jµν = r i jµν /r i jµν , r i jµν = r iµ − r jν is the separation vector, r i jµν = |r i jµν | is the distance between vertices, and σ i jµν = 1 2 (σ iµ + σ jν ) is their average diameter.(See Fig. 10 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.
The diameter of adipocytes in human breast adipose tissue ranges from 50 to 150µm, 23 while cancer cells are much smaller, ranging between 10 and 20µm. 24.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: where riµ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: where rqs = r qs /r qs , r qs = r q − r s is the separation vector, r qs = |r qs | is the distance between cancer cells q and s, and σ qs = 1 2 (σ q + σ s ) is the average diameter of the cancer cells.As shown in Fig. 2, the pair force between cancer cells is repulsive for r qs < σ qs , and attractive for 1 < r qs /σ qs < r α /σ qs with maximum attractive force F qs = ε c β /σ qs at r qs = r β , where r α /σ qs = 1 + α and r β /σ qs = 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. 23,24

Equations of motion for cancer cells and adipocytes
To determine the spatiotemporal trajectories of adipocytes and cancer cells, we integrate Newtwon'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 r q of cancer cell q is 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.3), and F qs is the pair force between cancer cells (Eq.4).F wq indicates the repulsive forces between physical walls and the qth cancer cell and is defined in Eq. 12.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 γσ / √ mε c = 0.2 and a self-propulsion term with active force magnitude f 0 .The direction of the active force nq obeys rotational diffusion 25,26 , where χ q is a vector with x-, y-, and z-components that are Gaussian random numbers with zero mean and satisfy where ⟨.⟩ t ′ is an average over time origins, I is the identity matrix, δ qs = 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 nq : 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 27 .In this model, we include a background temperature 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.
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: where the first and second terms include interactions among cancer cells and between cancer cells and adipocytes.As in Eq. 5, the equation of motion for passive cancer cells includes a damping force proportional to the velocity v q of cancer cell q with damping coefficient γ.The last term is the noise term with strength 2k b T 0 γm that sets the target temperature T 0 .η(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: 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 wiµ is the repulsive force between the physical walls and vertex µ on the ith adipocyte, which is defined in Eq. 13.The fifth term gives the mechanical forces generated from the shape-energy function for the ith adipocyte in Eq.

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σ , N a = 28 adipocytes, N c = 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.Compression is followed by energy minimization until the final packing fraction of the system reaches φ = ( π 6 ∑ 72, where L f = 34σ is the final side length of the box, and V i is the total volume of the ith adipocyte, i.e. v i plus the additional volumes from the vertices that do not overlap the polyhedron.(See Fig. 1

(e).)
During compression, we apply periodic boundary conditions in the yand z-directions.In the x-direction, we include two physical walls, a "right" wall at x r = L f and a "left" wall at x l = 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: where m is the mass and γ is the damping coefficient for the left wall, v l 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: and where x wq = x w − x q and x wiµ = 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: For the passive model of cancer cells, the temperature of the cancer cells T c and adipocytes T a equilibrate.For the active model of cancer cells, the temperature of the cancer cells is larger than that for the adipocytes.In the systems considered here, 1 < T c /T a < 1.2.

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 t 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 in- dicates the onset of cancer cell invasion.We find that E c 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 A t relative to untethered adipocytes.Finally, we show that heterogeneity in the local structural and mechanical properties of tethered adipocytes restricts cancer cell invasion.

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 inte-grate their equations of motion (i.e.Eqs. 5 or 9 for cancer cells and Eq. 10 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.9).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 t 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. 3 (a), we show that the time dependence of A t strongly varies with k b T /(Pσ 3 ).The minimum possible value of A t is A min t = 470σ 2 , which corresponds to the initial, de-mixed value (Fig. 3 (b)), and the maximum possible value of A t is A max t = 2000σ 2 , which corresponds to the surface area of the N a adipocytes (Fig. 3 (d)).Thus, A t 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 A t by its minimum and maximum values, while accounting for dilation of the system at high temperatures, as discussed in Appendix C. In Fig. 4 (a), we show the normalized interfacial area A n 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, A n is sigmoidal when plotted versus 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 A n = 1 for all k b T /(Pσ 3 ) that are above the glass transition temperature, where structural relaxation times diverge.

Effects of cancer cell cohesion and persistence on invasion
In the previous section, we found that the degree of invasion A n is time-dependent, evolving toward the mixed (invaded) state for purely repulsive, passive cancer cells.Next, we examined how A n is affected by cancer cell cohesion and persistence.In Fig. 4 (b), we show that A n versus k b T /(Pσ 3 ) for systems with attractive, active cancer cells does not possess strong time dependence.Moreover, the value of A n 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. 4 (b).In Fig. 5 (a), we plot the steady-state A n versus k b T /(Pσ 3 ) for three values of the attractive strength β and five values of the persistence time τ P .A n is sigmoidal versus k b T /(Pσ 3 ) for all β and τ P combinations studied.Thus, we can describe A n using where a(β , τ p ) controls the value of k b T /(Pσ 3 ) at which the rapid increase in A n occurs and 0.5 < b < 2.5 controls the rate of this increase.We demonstrate high-quality fits of A n to Eq. 15 in Fig. 5 (a).Thus, when we plot A n as a function of

E b c
, where E c ≡ k b T /(aPσ 3 ), the data collapses as shown in Fig. 5 (b).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. 5 (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, and At large τ p , τ d deviates from the analytical result as φ c increases.In the inset to Fig. 5 (b), we verify that a ∼ τ γ /τ d .At fixed τ p , A n 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, A n 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 timedependent when cβ ≫ c 1 (t).We show in Fig. 5 (b) that a = c 1 1 + c 2 β ε c /(Pσ 3 ) τ γ /τ d , where c 1 ≈ 0.06 enforces A n = 1/2 at E c = 1 and c 2 ≈ 0.02.

Tethering of neighboring adipocytes
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, µ=1 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 d c , such that the number of neighbors per adipocyte is 6.The potential energy associated with the spring connecting adipocytes i and j is where K ecm /σ 2 is the stiffness of the spring, R i j = |R i − R j |, and l i j0 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 In Fig. 7 (a) and (b), we calculate the interfacial area A n versus k b T /(Pσ 3 ) for several values of K ecm and two values of τ p and β .As discussed above, A n versus k b T /(Pσ 3 ) is normalized between 0 and 1 using the interfacial area obtained for untethered (K ecm = 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 max n < 1 at large k b T /(Pσ 3 ) for all K ecm > 0, but the values of k b T /(Pσ 3 ) at which A n > 0 does not vary strongly with K ecm .In the insets to Fig. 7, we show that A max n 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.

Packing fraction heterogeneity
The mechanical properties of adipose tissue are affected by many factors including obesity, inflammation, and lipolysis. 6,13,29Previous experimental studies of cancer cell migration through collagen networks demonstrated that geometrical confinement strongly affects cancer cell motion 15,30 .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 j0 = R i j and K ecm = 0.04ε c , and calculate A n versus φ a for τ p / mσ 2 ε c = 25 and β = 10 −2 .In Fig. 8 (a), we confirm that A n decreases with increasing adipocyte packing fraction and vanishes for φ a ≳ 0.85.
In addition to packing fraction, the local structural and mechanical properties of adipose tissue are often heterogeneous. 31We model structural and mechanical heterogeneity in adipose tissue by resetting l i j0 of the tethering springs for half of the adipocytes in the z-direction to l i j0 = (1 − λ )R i j and to l i j0 = (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-meansquare fluctuations in the local packing fraction that increase with λ .As shown in Fig. 8 (b), the maximum values of ∆φ are ≈ 0.1, 0.2, and 0.25 for λ = 0, 0.3, and 1.In Fig. 8 (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 t is bounded by the area of the interface in the initial de-mixed state A min t and the total surface area of the adipocytes A max t , and we normalized the interfacial area by these two values so that 0 < A n < 1.We then calculated A n as a function of temperature, pressure, and the cohesive strength and persistence of cancer cell motion.We showed that A n can be collapsed when plotted as a function of a dimensionless energy scale E c , where E c ≪ 1 indicates a de-mixed (noninvaded) state and E c ≫ 1 indicates a mixed (invaded) state.We found that E c 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 A n relative to more uniform local packing.
The current computational studies represent a first step toward modeling breast cancer invasion into adipose tissue as a physical process.However, there are numerous important 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. 15,27In 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. 32In future computational studies, we can incorporate separate fields into the simulations to account for local variations in collagen density and alignment, couple these fields to cancer cell velocities, and quantitatively compare 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. 15,33,34Thus, 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. 6To model this effect, we can allow the equilibrium volume v i0 of adipocyte i to decrease over time when in contact with cancer cells.
Finally, in our computational studies, we assumed that the cancer cell invasion time scale τ inv is shorter than the timescale 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 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. 35This 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 .(c) Distance r qs between cancer cells q and s.

FIG. 1 .
FIG. 1.(a) and (b) Histological slices of the interface between a mammary tumor and adipose tissue in mice.Cancer cell nuclei are colored reddish-purple and the large white regions are adipocytes.(a) Noninvaded 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 N a = 128 and N c = 7000.(f) Adipocytes are modeled as deformable polyhedra with vertices of diameter σ iµ and a triangular mesh with individual triangle area a ik .Cancer cells are modeled as soft spheres with diameter σ q .

2 FIG. 3 .
FIG. 3. (a) Interfacial surface area A t (in units of σ 2 ) shared between the cancer cells and adipocytes plotted as a function of time (in units of mσ 2 /ε c ).We show A t for five values of k b T Pσ 3 = 1.25 × 10 −3 , 7.89 × 10 −2 , 1.25 × 10 −1 , 1.98 × 10 −1 , and 1.25 increasing from blue to red.Simulation snapshots for different A t 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.

1 FIG. 4 .
FIG. 4. Normalized interfacial area A n 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 β = 0 and τ 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 A n 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. 5 .
FIG. 5. (a) Normalized interfacial area A n is plotted versus 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. 15.(b) A n is plotted versus the dimensionless energy scale E b c for the same data in (a).The inset shows the relationship between the parameter a in Eq. 15 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 .

1 FIG. 7 .
FIG. 7. A n is plotted versus k b T Pσ 3 for cancer cells with (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 K ecm = 0 are indicated by black filled circles.The dashed horizontal lines indicate A max n for different values of K ecm .Insets: A max n is plotted as a function of K ecm /ε c for the data in the main panels.The black dashed horizontal lines indicate A max n in the K ecm /ε c → ∞ limit.

(
See Fig. 6 (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. 6 (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 28 :

FIG. 8 .
FIG. 8. (a) The degree of invasion A n plotted versus the packing fraction of adipocytes φ a at τ P = 25 mσ 2 ε c , β = 10 −2 and k b T /(Pσ 3 ) = 0.37.λ increases from dark to light green.(b) Standard deviation in the local packing fraction of adipocytes ∆φ a plotted versus φ a .Open dark green circles indicate results for uniformly distributed adipocytes at a total packing fraction φ = 0.72.

FIG. 9 .FIG. 10 .
FIG. 9. (a) A representative histology slice of adipose tissue from a C57BL6 mouse.(b) Cell boundaries are identified through standard image processing techniques.Cells are colored by A 2D = p 2 /(4πs), where p is the perimeter and s is the area.(c) Probability distribution of A 2D P(A 2D ) for random 2D cross-sections of adipocytes in 3D from 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.The experimental data includes ∼ 1000 adipocytes collected from histology samples from 5 mice.