In this paper, we aim to investigate the mechanisms involved in cancer spheroid growth and invasion into the extracellular matrix. Spheroid growth refers to the expansion of the central spheroid mass and its volume change over time as a result of cell proliferation. On the other hand, invasion describes the penetration of single cells or broad multicellular protrusions into the surrounding ECM. We build the model based on the in vitro experiments conducted by Jahin et al. (2023) studying the effect of ribose-induced ECM stiffening on cancer spheroid growth and invasion of non-invasive MCF7 and invasive HCC1954 breast cancer cell lines.

In their study, Jahin et al. (2023) formed tumour spheroids of 200 μ m in diameter using the hanging drop method and subsequently embedded them in a collagen matrix with varying ribose concentrations of 0 mM, 50 mM and 200 mM as in Phillips et al. (2023). They used the non-invasive parental MCF7 and invasive parental HCC1954 cells, both human breast carcinoma cell lines with epithelial-like morphology. To model the ECM, they chose collagen I, derived from rat tail tendons, as it is the most abundant protein component in the extracellular matrix surrounding solid tumours, with supplementary fibronectin also included to allow enhanced cell attachment. During collagen hydrogel formation, ribose, a cross-linker used for non-enzymatic glycation to induce gel stiffening in vitro models, was added at appropriate concentrations to increase hydrogel stiffness. More cross-linking between collagen fibres increases the ECM stiffness without altering the matrix organisation and the ligand binding sites for cell-ECM adhesion. This allows for the investigation of the effect of stiffness alone on cancer spheroid growth and invasion. The spheroids were imaged by combining Z-slices at 10 μ m intervals, covering the whole spheroid thickness. The images were captured daily over at least 96 h and were used to track the spheroid invasion.

To compare our results with the experimental data, we calculated spheroid area growth relative to the initial time, cell count and Delaunay mean distance between cells in Python version 3.10.12. Given the stochastic nature of our model, we ran 10 replicates for each simulation and computed the mean and 25th/75th percentile of spheroid area growth relative to the initial time, cell count and Delaunay mean distance every 60 min. Simulations were performed on the University of Birmingham’s high-performance computing (HPC) cluster, BlueBEAR. We utilised a single node ( 2 × 56 -core Intel^® Xeon^® Platinum 8570) and ran batches of 20 simulations simultaneously, each allocated 4 GB of RAM. Completion times varied with the number of agents, ranging from approximately 45 s to 5 min, with most simulations finishing within 2–3 min.

We computed the spheroid area growth relative to the initial time by calculating the area covered by the cells at each time point and dividing it by the area covered at the initial time t = 0 min. The spheroid area was approximated by dividing the entire domain into a 5000 × 5000 , initially setting all grid elements to a value of 0. This baseline value represents unoccupied space. We then drew disks of value 1 at the coordinates of each cell’s centre with their corresponding radius (Figure 2A). Overlaps were ignored, as grid elements covered by multiple cells are only counted once. To calculate the total spheroid area, we summed the grid elements with value 1 and rescaled to the original domain size to obtain the spheroid area in μ m 2 . This process is analogous to the method used for computing spheroid invasion relative to the initial time in Jahin et al. (2023).

For the cell count, we tallied the total number of cells in the simulations at each time point, corresponding to the number of nuclei in a slice of the experimental data.

Finally, the Delaunay mean distance measures the proximity of cells. It utilises the Delaunay triangulation, the dual graph of the Voronoi diagram, of a set of points. The edges of the graph form triangles whose circumscribed circles do not contain any other point. We used the cell centres as the input nodes of the network and generated the Delaunay network using the spacial algorithm Delaunay from the Python library SciPy version 1.11.1 (Figure 2B). We then computed the mean edge length between the nodes to find the Delaunay mean distance.

Invasion and migration of cancer cells from a cancer spheroid into the surrounding extracellular matrix depends strongly on how they interact with the ECM (Yamada et al., 2022). The cancer cells need to remodel the ECM to enable invasion, and can then use the adhesion sites on the collagen fibres to propel themselves and invade further. In our model, to study such cell-ECM interactions, we consider different biophysical parameters representing cell-ECM cross-talk and specific properties of the ECM. We control how quickly the agent cells reduce locally the ECM density by changing the degradation rate ( r deg,0 , Equation 1). In turn, the ECM density affects the cell’s speed (Equations 4, 6), which reaches its maximum cell-ECM interaction speed ( S rib ) when the ECM density is equal to 0.5. The ECM density also affects the cell’s proliferation rate ( r div ) through the inhibition of proliferation function (Equation 8). In this section, we present the results of our analysis on the degradation rate ( r deg,0 ) and maximum cell-ECM interaction speed ( S 0 ) without any ribose at different proliferation rates ( r div ) . Then we analyse the δ and σ parameters, which determine how strongly ribose affects the degradation rate ( r deg,rib , Equation 2) and the maximum cell-ECM interaction speed ( S rib , Equation 5) respectively at different ribose concentrations ( r i b ) .

We initiated all simulations with a spheroid of cancer cells of 200 μ m in diameter and homogeneous ECM density with value ρ = 1 throughout, except at the spheroid’s location, where the ECM density is zero. We set the overcrowding threshold N m a x used in the inhibition of proliferation function (Equations 5–8), equivalent to a cell fully surrounded by the other cells (hexagonal packing) (Metzcar et al., 2025). All of the parameters used in the simulations are also listed in the Supplementary Tables S1–S5 and other PhysiCell specific parameters are set to their default values as used in PhysiCell 1.12.0 (Ghaffarizadeh et al., 2018). We examined the impact of the degradation rate ( r deg,0 ) and the maximum cell-ECM interaction speed ( S 0 ) on spheroid area growth relative to the initial time and Delaunay mean distance (explained in Section 2.3), holding the ribose concentration at 0 mM. The analysis was conducted for three proliferation rates: r div = 0.0004 min − 1 , 0.0006 min − 1 and 0.0008 min − 1 (Figures 3A, B).

Our results demonstrate that an increase in proliferation rate ( r div ) enhances spheroid area growth relative to the initial time (Figure 3A) and reduces the Delaunay mean distance (Figure 3B). This is expected, as faster cell division leads to a denser cell population, which is correlated to a lower average cell-cell distance.

Further, increasing the maximum cell-ECM interaction speed ( S 0 ) leads to enhanced spheroid growth and a higher Delaunay mean distance when the degradation rate r deg,0 is above 0.0001 min − 1 , and to little change when equal to 0.0001 min − 1 (Figures 3A, B). A higher S 0 allows more cells to detach from the spheroid, reducing the number of neighbours, and consequently avoiding proliferation arrest. It also enables cells to access and degrade more areas of the ECM, further promoting proliferation. Additionally, with faster cell migration, cells at the spheroid’s edge become more dispersed, contributing to the increase in Delaunay mean distance.

Finally, we observed that increasing the degradation rate r deg,0 does not always induce a monotonic increase in spheroid area growth relative to the initial time and Delaunay mean distance. Generally, increasing the degradation rate enhances spheroid growth (Figure 3A) and increases cell-cell distance, leading to a higher Delaunay mean distance (Figure 3B). This occurs because a higher degradation rate reduces ECM density around the spheroid, which in turn positively affects proliferation and invasion. However, when the degradation rate is excessively high ( r deg,0 = 0.0128 min − 1 ), the ECM is degraded too quickly, which inhibits migration and invasion. Therefore, for fixed proliferation rate r div and maximum cell-ECM interaction speed S 0 , as the degradation rate r deg,0 increases, the spheroid growth slows down (Figure 3A) and the Delaunay mean distance decreases (Figure 3B).

Figure 3C shows simulation images at the final time point (96 h) for S 0 = 0.2, 0.5 and 0.8 μ m ⋅ min − 1 , r deg,0 = 0.0002, 0.0008 and 0.00064 min − 1 and r div = 0.0004, 0.0006 and 0.0008 min − 1 . These images illustrate that higher S 0 values lead to greater cell dispersion at the spheroid’s edge. When combined with higher ECM degradation r deg,0 more single cells are observed migrating away from the spheroid. Thus, increased degradation rate and maximum cell-ECM interaction speed contribute to the formation of protrusions in the spheroid. In contrast, lower degradation levels and migration speed limit spheroid growth, resulting in a more rounded spheroid shape.

We then analysed the impact of the parameters δ and σ on spheroid area growth relative to the initial time (Figure 4). For δ , σ = 0 the functions are constant, so r deg,rib = r deg,0 and S rib = S 0 , while for δ , σ > 0 the functions are monotonically decreasing with respect to the ribose concentration r i b , so r deg,rib < r deg,0 and S rib < S 0 for ribose r i b greater than zero (Equations 2, 5). For this analysis, we fixed the degradation rate at ribose concentration 0 mM ( r deg,0 = 0.0032 min − 1 ), the maximum cell-ECM interaction speed at ribose concentration 0 mM ( S 0 = 0.7 μ m ⋅ min − 1 ) and the proliferation rate ( r div = 0.00072 min − 1 ). The analysis was conducted for two ribose concentrations: r i b = 50 mM and 200 mM.

Figure 4 shows that higher ribose concentrations correspond to a decrease in spheroid growth, and increasing either δ or σ also results in reduced spheroid growth. Notably, at 200 mM ribose, when σ is greater than or equal to 0.015 mM − 1 , spheroid growth remains unchanged for fixed values of δ . This occurs because for σ = 0.015 mM − 1 the maximum cell-ECM interaction speed is significantly reduced, S 200 = S 0 e − σ 200 ≈ 0.035 μ m ⋅ min − 1 compared to S 50 = S 0 e − σ 50 ≈ 0.33 μ m ⋅ min − 1 . As a result, as the cell population becomes denser, proliferation is inhibited throughout the spheroid except at its boundary, where proliferation depends on the surrounding ECM density. Consequently, spheroid growth is determined by the rate at which cells degrade the ECM at the boundary, facilitating increased proliferation in this region.

To begin with, we replicated experiments from Jahin et al. (2023) that study the effect of ribose concentration on two different cell lines of parental breast cancer cells: MCF7 and HCC 1954. MCF7 cells are a non-invasive cell line, which correlates with weaker cell-ECM interactions (Comşa et al., 2015). On the other hand, HCC1954 cells are a more aggressive and invasive cell line, characterised by enhanced contractility, and therefore stronger interactions with the ECM fibres enabling migration, and further ECM remodelling (de Abreu Pereira et al., 2022; Jahin et al., 2023). The in vitro experiments showed that with increasing ribose concentration, and therefore collagen fibre stiffness, the spheroid invasion was inhibited for the invasive HCC1954 cells, while the non-invasive MCF7 cells did not invade for any of the ribose concentrations.

The cell-ECM interactions in our model depend on two parameters: ECM degradation rate, which controls the ECM remodelling by the cells (Equation 1), and maximum cell-ECM interaction speed, which affects the cell’s movement (Equation 4). Given the different invasiveness of the two cell lines, we assume that the invasive cells have a higher ECM degradation rate and maximum cell-ECM interaction speed than the non-invasive cells. From the wide range of values studied in Section 3.1 (Supplementary Figure S3), we choose the values that lead to simulations matching the experimental observations (Figures 5E, F). Hence, we use the following set of values: ECM degradation rate r deg,0 = 0.0001 min − 1 and maximum cell-ECM interaction speed S 0 = 0.1 μ m ⋅ min − 1 for non-invasive cells, and ECM degradation rate r deg,0 = 0.0032 min − 1 and maximum cell-ECM interaction speed S 0 = 0.7 μ m ⋅ min − 1 for invasive cells. Further, the strength of the effect of ribose concentration on the cell behaviour depends on the parameters δ for degradation rate r deg,rib (Equation 2) and σ for maximum cell-ECM interaction speed S rib (Equation 5). Following Section 3.1 (Figure 4), we use δ = 0.02 mM − 1 and σ = 0.035 mM − 1 to match the spheroid growth after 96 h of the invasive cells for ribose concentrations of 50 mM and 200 mM (Figure 5Di). Finally, we set the proliferation rate r div = 0.00072 min − 1 , which is the default parameter for the live cell cycle model of MCF10A breast cancer epithelial cells in PhysiCell (Ghaffarizadeh et al., 2018). Following the in vitro experiments performed in Jahin et al. (2023), we set the initial spheroid diameter to be 200 μ m (corresponding to 139 cells), with homogeneous ECM density with value ρ = 1 throughout, except at the spheroid’s location, where the ECM density is zero. We study the effect of ribose concentrations, 0 mM, 50 mM and 200 mM, on spheroid area growth relative to the initial time, cell count and Delaunay mean distance, as shown in Figures 5C, D. Simulation images at 24 h and 96 h for non-invasive and invasive cells are shown in Figures 5A, B respectively.

In their experiments with non-invasive MCF7 cells, Jahin et al. (2023) observed that increasing ribose concentration did not affect spheroid growth (Figure 5Ei), consistent with previous reports (Ziegler et al., 2014). Our findings also indicate that spheroid growth was inhibited for the non-invasive cells, with the spheroid area growth at the final time point remaining below 2 for all ribose concentrations (Figure 5Ci). However, while the invasion of MCF7 cells was low, Jahin et al. (2023) observed an increase in the number of nuclei, particularly at 0 mM ribose (Figure 5Eii). Similarly, our simulations showed a larger increase in cell count at 0 mM ribose, with the cell count reaching ∼ 280 after 96 h (Figure 5Cii), which matches the corresponding mean number of nuclei observed in vitro (Figure 5Eii). Interestingly, they found that the Delaunay mean distance increased for ribose concentration of 0 mM and it decreased at 50 mM and 200 mM, indicating a denser spheroid for higher ribose concentrations (Figure 5Eiii). In our simulations, we observed that the Delaunay mean distance between 24 h and 96 h slightly increased at ribose concentration 0 mM and decreased for ribose 200 mM, but remained constant at 50 mM (Figure 5Ciii). This means that, in our simulations, the spheroid at 50 mM is less dense than in the experimental data. This can be due to more rapid degradation of the ECM or slower cell proliferation than in the experiments. We also observed a larger decrease in Delaunay mean distance for ribose 200 mM in our simulations compared to the experimental data. This difference could be attributed to the absence of cell death in our in silico model, whereas the experimental data show a reduction in the number of nuclei over time likely due to cell death (Figure 5Eii). The higher cell count maintained in the simulations likely results in a denser spheroid and therefore a larger decrease in Delaunay mean distance.

In contrast, the invasive HCC1954 cells exhibited a reduction in spheroid growth and the number of nuclei with a ribose concentration increase in vitro (Figures 5Fi, ii). Our simulations closely matched the experimental data, resulting in a spheroid growth of ∼ 6 after 96 h and a cell count of ∼ 570 after 72 h for ribose 0 mM, spheroid growth of ∼ 3.4 after 96 h with a cell count of ∼ 370 after 72 h for ribose 50 mM, and spheroid growth of ∼ 1.8 after 96 h and a cell count of ∼ 250 after 72 h for ribose 200 mM (Figures 5Di, ii). Further, Jahin et al. (2023) found that the Delaunay mean distance at 72 h was higher than that at 24 h, for all ribose concentrations considered. However, it decreased with an increase in the ribose concentration (Figure 5Fiii). Our simulations showed that the Delaunay mean distance at 24 h and 72 h is almost the same for each ribose concentration, but lowers as the ribose increases (Figure 5Diii). However, we see that the Delaunay mean distance does not remain constant over time (Supplementary Figure S4Biii). We noticed that the Delaunay mean distance at ribose 0 mM initially increases and peaks between 24 h and 48 h, before decreasing. While we see an overall increase in Delaunay mean distance for ribose concentration 0 mM, as seen experimentally, our model predicts a rapid increase of Delaunay mean distance and has a similar value at 24 h and 72 h, in contrast to the experiments. With only two experimental time points, we cannot capture the dynamics in the first 24 h and between 24 h and 72 h. Including more time points would allow us to better understand variations over time and a more accurate representation of how the spheroid evolves.

With our choice of parameters, we observe that low degradation rate and maximum cell-ECM interaction speed inhibit both invasion and proliferation of the spheroid. This is mainly due to the cell-ECM repulsive velocity v cmr (Equation 6) and the inhibition of proliferation function (Equation 8). As observed in the parameter analysis in Section 3.1, low degradation rate and maximum speed make the spheroid denser thanks to the ECM acting as a wall because of the repulsive velocity. Proliferation is inhibited at the centre of the spheroid due to the high number of neighbours and at the edge due to the high ECM density. The chosen values for δ and σ in our simulations give low degradation rates and maximum cell-ECM interaction speeds for the invasive cells at ribose concentrations of 50 mM and 200 mM. The degradation rates at ribose concentrations 50 mM and 200 mM are r deg,50 ≈ 0.001 min − 1 and r deg,200 ≈ 0.00006 min − 1 , while the maximum cell-ECM interaction speeds are S 50 ≈ 0.1 μ m ⋅ min − 1 and S 200 ≈ 0.0006 μ m ⋅ min − 1 . This indicates that both ECM degradation and cell speed are substantially reduced for the invasive cells as the ribose concentration increases, which lowers tumour invasion and proliferation. Thus, our simulations are in line with the observation by Jahin et al. (2023), that ribose-induced cross-linking of collagen possibly reduces ECM remodelling and migration, slowing spheroid growth and invasion.

Cancer cells remodel the extracellular matrix mechanically and proteolytically when invading, creating paths that facilitate the migration of nearby attached cancer cells (Walker et al., 2018). The cells mechanically apply forces to the ECM fibres by pushing or pulling the fibres when adhering to ligand binding sites, resulting in fibre displacement and orientation changes. Fibre orientation can direct migration, and the realignment of the collagen fibres has been associated with higher invasion (Provenzano et al., 2006). On the other hand, proteolytic remodelling involves enzymatic degradation of ECM fibres through the activity of matrix metalloproteinases (MMPs). It has been shown that ECM degradation by MMPs significantly contributes to cell invasion as it facilitates migration and realignment of the fibres (Itoh, 2015). Jahin et al. (2023) investigated the role of MMPs in fibre alignment and invasion by treating the invasive HCC1954 cells with the pan-MMP inhibitor GM6001. They observed that MMPs inhibition significantly reduces invasion for ribose 0 mM, but not for ribose 50 mM and 200 mM (Figure 6C).

In our model, ECM remodelling is represented by the degradation rate parameter r deg,0 , which is responsible for changes in ECM density. We compared the experimental results in control conditions (DMSO) with the invasive cells simulations discussed in Section 3.2. Following the parameter analysis from Section 3.1 (Supplementary Figure S3), we choose the degradation rate matching the experimental results of growth relative to the initial time after 72 h of the invasive HCC1954 cell line with the addition of the pan-MMP inhibitor GM6001 (Figure 6C). We selected the degradation rate r deg,0 = 0.0004 min − 1 . This non-zero value arises because the pan-MMP inhibitor only blocks ECM degradation, but not mechanical remodelling of the fibres, which makes it possible for the cells to locally change the density of the ECM even without degrading the fibres. Simulation images at 72 h for invasive cells with high (Control) and low (GM6001) degradation rates are shown in Figure 6A.

Similarly to the results from the experimental data by Jahin et al. (2023) shown in Figure 6C, spheroid area growth relative to the initial time after 72 h was significantly reduced for ribose 0 mM going from ∼ 4 in the control conditions simulation (Control) to ∼ 2 with the reduced degradation rate (GM6001) (Figure 6B). At ribose concentration of 50 mM spheroid growth was less affected by the MMP inhibition, both in vitro and in the simulations. Finally, at 200 mM, both the Control and the reduced degradation rate (GM6001) conditions have similar spheroid sizes in the simulations after 72 h, which is in line with the experimental data in the DMSO and GM6001 conditions (Figure 6C).