A PDE Model of Breast Tumor Progression in MMTV-PyMT Mice.

The evolution of breast tumors greatly depends on the interaction network among different cell types, including immune cells and cancer cells in the tumor. This study takes advantage of newly collected rich spatio-temporal mouse data to develop a data-driven mathematical model of breast tumors that considers cells' location and key interactions in the tumor. The results show that cancer cells have a minor presence in the area with the most overall immune cells, and the number of activated immune cells in the tumor is depleted over time when there is no influx of immune cells. Interestingly, in the case of the influx of immune cells, the highest concentrations of both T cells and cancer cells are in the boundary of the tumor, as we use the Robin boundary condition to model the influx of immune cells. In other words, the influx of immune cells causes a dominant outward advection for cancer cells. We also investigate the effect of cells' diffusion and immune cells' influx rates in the dynamics of cells in the tumor micro-environment. Sensitivity analyses indicate that cancer cells and adipocytes' diffusion rates are the most sensitive parameters, followed by influx and diffusion rates of cytotoxic T cells, implying that targeting them is a possible treatment strategy for breast cancer.

1. Introduction

Breast cancer is the most common type of cancer in women and a major public health problem [1]. It accounts for 25% of the new female cancer cases around the world [2], and it was the cause of 43,600 deaths in the United States alone in 2021 [3]. Breast cancer has three major subtypes: human epidermal growth factor 2 positive (HER2+) (70% of patients), HER2- (15–20%), and triple-negative (tumors lacking all three standard molecular markers; 15%) [4]. Based on the stage of the disease, different treatment options are practiced, such as chemotherapy, radiation therapy, surgical removal, or a combination of these [4,5]. However, these treatments can be straining and do not always work. A comprehensive understanding of the biology of cancer as a complex system of interactions is essential for obtaining effective treatments.

The tumor micro-environment has been the focus of many studies, including therapeutic targeting [6,7] and spatial heterogeneity [8,9,10], as many scientists believe that the key to curing cancer lies within tumor micro-environment interactions [11,12]. However, to find effective therapies, many factors should be considered, and in-vivo examinations of the tumor micro-environment and its response to treatments can be expensive and straining for both patients and investigators. Therefore, mice models of breast cancer are prevalent for testing different hypotheses and treatment combinations.

To establish mouse models, many different approaches have been used, including transgenic [13,14,15,16], gene targeting [17,18,19], and RNA interference [20,21,22]. Transgenic mouse models (in which oncogenes can be expressed while other factors, such as tumor-suppressor genes, are muted), are used to study breast cancer and improve our understanding of cancer initiation and progression [23,24,25,26]. For example, using these mouse models, it has been shown that macrophages regulate the cancer progression and formation of a high-density vessel network [27], and the high mobility group box-1 (HMGB-1) proteins play roles in promoting angiogenesis and tumor migration [28]. One of the most commonly used mouse models in cancer research is the mammary-specific polyomavirus middle T antigen overexpression mouse model (MMTV-PyMT) [29,30,31]. Specifically for breast cancer, this mouse model has been utilized successfully and qualifies as a transgenic approach [32,33].

At the same time, mathematical modeling of cancer development and tumor micro-environment offers insights and can be used in discovering new treatments [34,35,36,37,38,39,40,41,42,43,44,45,46,47]. As the complex spatial cell-to-cell interactions in the tumor micro-environment has attracted many experimental studies, a more thorough mathematical model can help scientists gain a better insight into the mechanisms of cancer growth. Many mathematical models using partial differential equations (PDE), which allow the integration of spatial information into the governing differential equations, have been developed to study health problems, such as atherosclerosis, Alzheimer’s disease, and COVID-19 spread [48,49,50,51,52]. Since the pioneering work of Iwata et al. [53] that introduced a PDE dynamical model for the metastatic evolution of an untreated tumor, researchers have developed several PDE models for cancer [54,55,56,57]. For example, for breast cancer, a PDE model with a 1D domain was proposed in which a predetermined population of tumor cells interacts with a fixed population of macrophages in a tissue-culture experiment [58]. In another study, a 2D spatial hybrid model was suggested, where an ordinary differential equation (ODE) described the dynamics of each immune cell, and a PDE described the evolution of the average substances released by the tumor cells [59]. Lai et al. also developed a system with eight variables, including cell populations of macrophages and T cells, to model the combination therapy for breast cancer [60].

Biological tissues have specific mechanical properties which can affect the growth and deformation caused by cell proliferation and movement in addition to the important cellular and molecular interactions. Since the pioneering works of Y.C. Fung on modern biomechanics [61,62,63,64], many specific tissue constitutive models have been proposed [65,66,67,68]. As mechanical properties of tumors can differ depending on the hosting tissue, different mechanical models have been developed for different cancer types. In particular, breast tumors have been treated as porous media and fluid-like tissue. For example, Frieboes et al. modeled breast tumor tissues as porous media and implemented a model that incorporates the interplay between the local drug, oxygen, and nutrient concentrations [69]. On the other hand, Friedman and Hu argued that due to a high content of the extracellular fluid in mammary glands, breast tumors can be modeled by the Stokes equation [70], and their method has been successfully applied to many studies [71,72,73,74].

This paper proposes a system of PDEs consisting of 15 variables coupled with the Stokes equation governing the tumor micro-environment. The coupling is performed through the velocity involved in the convection terms of the PDE system. This study improves our previous model developed via a system of ODEs [75]. We incorporate the parameter estimations for the ODE system in this paper and extract the PDE-specific parameters from the literature. Spatially dependent initial conditions for each cell type are obtained from an MMTV-PyMT mouse model by using multiplexed immunohistochemistry and investigating certain combinations of biomarkers such as Epcam, CC3, CD45, CD3, CD4, CD8, CD11c, F4/80, CSF1R, and MHC-II in the experimental domain. The simulations are done via the finite element method (FEM). We first investigate a naive case similar to our experimental data for which we do not consider any immune cell influx. Later, we include an influx of specific immune cells. We study tumor growth, cell and molecules dynamics, the spatial significance of immune cell distributions, and the effect of immune cell influx on cancer and domain growth. Finally, we perform an adjoint-based sensitivity analysis to find the sensitivity of total cancer cells to all the parameters of our model. We comment on the similarity of our results with the ODE study and investigate the biological importance of the PDE-specific parameter sensitivities.

2. Materials and Methods

The interaction network consists of helper T cells, cytotoxic T cells, regulatory T cells, dendritic cells, macrophages, adipocytes, cancer cells, and necrotic cells and molecules, including HMGB1, interleukin-6, interleukin-10, and interleukin-12, see Figure 1.

Figure 1

Interaction network. Interaction between key cells and molecules.

The effect of nutrients and metabolites on tumor growth has been studied extensively. They can effect cancer cell growth directly or through the immune cells such as T cells [76,77,78]. However, to avoid too much complexity, we neglect the effect of nutrients and metabolites in this study.

2.1. PDE System

While ODE models give valuable information about cell and molecule dynamics, they cannot capture the spatial effect of cell-cell and cell-molecule interactions. For a more realistic model of the cancer evolution, we considered a Reaction-Diffusion-Advection (RDA) PDE system, in which the reaction terms come from the ODE system mentioned above [75]. The advection was directed by a velocity field acquired from the mechanical properties of the tumorous tissue and triggered by the cells’ growth in the cancer micro-environment. We assumed that naive T cells and naive macrophages get activated outside the tumor micro-environment [79]. Therefore, the region’s deformation and growth depend on the model’s other cell types: helper, cytotoxic and regulatory T cells, naive and activated dendritic cells, activated macrophages, cancer cells, necrotic cells, and adipocytes.

We used a common approach of modeling cells’ and molecules’ movements and diffusion rates in the tumor [48,80,81,82,83,84]. Given the significant difference between the size of cells and molecules, in these studies, it was assumed that the domain volumetric change is only affected by the movement of cells inside of the tumor microenvironment and not the molecules. We, therefore, only considered advective terms in the PDEs for cells in the tumor and not the molecules or cells outside of the tumor. These assumptions lead to the following system of PDEs: where for corresponds to the 15 variables form the ODE system given in [75]. Table 1 indicates the list of the variables of the ODE model and their corresponding names in the PDE model. The operators and Δ are the divergence and Laplacian, respectively. The vector represents the advection velocity, and  is the diffusion coefficient for the cell or molecule . The function represents the right-hand side (biochemical reactions) of the i-th ODE equation, and

Table 1

PDE and ODE variables. This table shows the relationship between the variables from (1) and the system of ODEs in [75].

Variable in PDE	Variable in ODE	Name
X_1	T_N	Naive T cells
X_2	T_h	Helper T cells
X_3	T_C	Cytotoxic cells
X_4	T_r	Regulatory T cells
X_5	D_N	Naive dendritic cells
X_6	D	Activated dendritic cells
X_7	M_N	Naive macrophages
X_8	M	Activated macrophages
X_9	C	Cancer cells
X_10	N	Necrotic cells
X_11	A	Cancer associated Adipocytes
X_12	H	HMGB1
X_13	IL_12	IL-12
X_14	IL_10	IL-10
X_15	IL_6	IL-6

To exclude the effect of molecules and outsider cells on volumetric changes, indices correspond to the naive T cells and naive macrophages, and  correspond to the molecules. Similarly,

Because we did not model the diffusion outside the tumor micro-environment, the equations for naive T cells and naive macrophages remained as ODEs. Furthermore, we assumed that all cells in the model have approximately the same volume and surface area, and hence the same diffusion coefficient [50,81,83].

These small constant diffusion rates for cells help with smoothness while avoiding too much complexity in the sensitivity analysis. This is a common simplification [48,80,81,82,83] that leads to advection-dominated PDEs for cells.

Boundary Conditions

Immune cell infiltration in cancer studies has important prognostic implications [85,86,87,88,89]. In order to control and analyze the rate and intensity of infiltration of these cells, we used Robin boundary conditions through the tumor’s outer boundary, which has been used to model the influx rate of cells in tumors [90,91]. Namely, where is the boundary of the tumor at time t, the vector is the outward unit normal vector to the boundary, and  is the influx rate and is only nonzero for immune cells; we considered zero flux boundary conditions for the rest of the cells and molecules. The quantity pertains to the maximum levels of immune cells in lymph nodes and blood.

2.2. Mechanical Model

Many mathematical models of tumors assume that the tissue is a porous medium [92,93,94]. Then, due to the high permeability of macrophages and other cells [95,96], one can treat them as a low-speed flow through the porous tissue with the advection velocity [90]. On the other hand, some studies consider the tumor as fluid without a solid structure interaction [72,80,97]. Especially in the case of breast cancer, they argue this approach is reasonable since the tumor is mainly confined in the mammary gland, which has a high content of the extracellular fluid. For this study, we considered the breast tumor as a fluid, and we followed the method introduced in [70]. We also assumed that the changes in the volume and surface area of the cells in the cancer micro-environment were negligible and the domain was an incompressible, continuous fluid, with no voids inside. Hence, the sum of the densities of all cells remained constant [48,80]. If we take (with indices 1 and 7 missing on purpose) to be the set of indices corresponding to the cells present in the tumor micro-environment, then

Since the breast tissue is mainly decomposed into water, lipid, and protein with corresponding mass densities of 1, 0.924, and 1.35 g/cm3 [98], for simplicity, we assume that the constant in (8) is on average 1. Hence, summing both sides of (1) over and applying (8) implies:

For modeling the fluid-like behavior of the breast tumor, we use the Stokes equation: where is the tumor region at the time t and with and p being viscosity and the hydrostatic pressure. We assume that cell-cell adhesion force on the boundary of a tumor keeps the domain connected [99,100]. Taking to be this force, to be the point-wise mean curvature of the boundary, and  to be the unit outward normal vector for the boundary , we have the following boundary condition:

Additionally, we assumed the kinematic boundary condition: where is the velocity of the free boundary in the direction of . Finally, since the problem has a 6 dimensional kernel in 3D (and a 3 dimensional one in 2D) of the form , with and being arbitrary vectors and being the deformation vector, we consider the constraints to exclude rigid body movements such as translation and rotation. The system consisting of Equations (9), (10) and (12)–(14) has a unique solution [101].

There are no precise reports of the values and in the literature. Rianna and Radmacher [102], for thyroid cancer 350,000, and Sancho et al. report an approximate value of for general cell-cell adhesion force () in-vitro [103]. For this study, we scaled the problem so we can take , and this scaling would result in . Since we worked with mouse data, we tended to go with the lower bound. We, therefore, took the scaled to be exactly . However, we acknowledge this ad-hoc estimation as a limitation of our study.

2.3. Data of the Model

2.3.1. Mouse Model and Experiments

Tumors from naïve mouse mammary tumor virus-polyoma middle tumor-antigen (MMTV-PyMT) mice were harvested at 3–5 mm in diameter (early), and 12–15 mm in diameter (late) size and were prepared as formalin-fixed, paraffin-embedded (FFPE) samples. The 4–5 μm tumor sections were stained using multiplex immunohistochemistry (mIHC)—a process of serial immunostaining, imaging, and stripping—to assess a range of markers with specific staining patterns being cross-validated by using cyclic immunofluorescence (cycIF). Each mIHC image was analyzed by segmenting individual cells and calculating marker positivity for each segmented cell. For this study, we developed a comprehensive mouse-specific readout panel including proteins such as Epcam, CC3, CD45, CD3, CD4, CD8, CD11c, F4/80, CSF1R, and MHC-II to interrogate a broad range of tumor and tumor microenvironment states and functions. Table 2 lists cell classification based on biomarker combination. Notice that since we consider all subtypes of activated macrophages as one variable, we have considered two combination biomarkers for it. For regulatory T cells, we assume T cells that are neither cytotoxic nor helper are regulatory. Moreover, the combination Epcam(−) CD45(−) may include fibroblasts, endothelial cells and pericytes, and adipocytes. However, here, we assume it only refers to adipocytes.

Table 2

Biomarker combinations. (+) means high expression and (−) means lack of expression of a protein at a certain location.

Cell Type	Biomarker Combination
Helper T cells (T_h)	Epcam(−) CD45(+) CD3(+) CD4(+) CD8(−)
Cytotoxic T cells (T_C)	Epcam(−) CD45(+) CD3(+) CD4(−) CD8(+)
Naive dendritic cells (D_N)	Epcam(−) CD45(+) F4/80(−) CD11C(+)
Dendritic cells (D)	Epcam(−) CD45(+) F4/80(−) CD11C(+) MHC-II(+)
Activated macrophages (M)	Epcam(−) CD45(+) F4/80(+) CD11C(−) CSF1R(+) or Epcam(−) CD45(+) F4/80(+) CD11C(−) CSF1R(−) MHC-II(+)
Cancer cells (C)	Epcam(+) CD45(−)
Necrotic cells (N)	CC3(+)

The tumors were derived from immunocompetent MMTV-PyMT mice with spontaneously growing tumors that mirrored the morphology and aspects of progression of human breast cancer [30]. For the details of the experiments, please see Appendix B. Table 2 explains the biomarker combinations.

2.3.2. Preparation of Initial Conditions

To avoid the instability caused by an irregular boundary, we used an elliptical domain containing all the model’s cell types. Figure 2 shows the mathematical region superimposed on the experimental domain. Single cell-type locations were determined by the bio-marker combinations given in Table 2.

Figure 2

Nine figures showing the position of the chosen elliptical domain compared to each cell type. Blue dots represent a single cell of the corresponding cell type, and gray dots are the rest.

We used a triangular mesh on the mentioned elliptical domain. We assigned a discontinuous Galerkin function space of degree zero to this mesh. Then, for each cell type, we defined piece-wise functions as follows: where is the number of the cell type in the triangle . The function for a domain point x is a characteristic function defined by:

We project the function defined in (15) onto a function space with linear Lagrangian elements to get a continuous representation for the initial state of each cell type. However, this projection might result in non-smoothness. Figure 3 shows a one-dimensional case of this issue.

Figure 3

Solid red: A discontinuous function. Dashed blue: The projection onto a finite element space with linear Lagrangian bases.

To prevent the propagation of such anomalies in our simulation, we flatten the negative values and then introduce a primary diffusion step to smooth things out. After this step, we non-dimensionalize each field by dividing it by its maximum value across the domain. This concludes the initial condition preparation. See Appendix C and Figure A2 for a visualization of these steps and the final products, which we use as the initial conditions for solving the system (1).

Figure A2

Column 1: Discontinuous fields. Column 2: Projection onto a function space with linear Lagrangian elements. Column 3: Smoothened fields via diffusion. Column 4: Non-dimensionalized fields.

The parameters involved in the reaction term of Equation (1) are evaluated based on an estimation method proposed in a previous study [75]. In that study, we had time-course data for three PyMT mice, and we performed a least-squares optimization to obtain the reaction parameters. Here, we consider constant initial conditions for HMGB1, IL-12, IL-10, and IL-6; these constants are taken to be the initial values of mouse 1 from the ODE model presented in [75]. Finally, as mentioned before, Equations (1)–(3) for naive T cells and naive macrophages are actually ODEs. Therefore, we take their initial conditions to be a constant of 1.

3. Results

3.1. No Influx

We used a finite element method to simulate our results (see the Appendix D). We started by investigating a case with no influx source for the immune cells. In other words, we tooko in (7). Using the initial conditions from the fourth column of Figure A2, constant initial conditions for cytokines (taken from [75]) and naive T cells and naive macrophages (taken to be 1), we solved the discrete mechanical and biological problem discussed in Appendix D. Figure 4 shows the change in the length and width of the bounding box containing the domain. For this case, we carried out our simulations for 600 h (25 days). It is worth pointing out that the extractions of early and late tumors mentioned in Section 2.3.1 were about three weeks apart. The almost circular form of the domain is due to the effect of the boundary condition (12), and the deformation generally happens in the direction of the mean curvature. Even though for this problem, we have conveniently picked an elliptic reference region; cancer scientists commonly observe a blob-shaped final results [104,105,106]. The small deviations from a fully circular region in Figure 4 and the rest of this paper is the result of the competition between deformation by the reaction term (9) and the boundary term (12) throughout the region.

Figure 4

Comparison between the dimensions of the tumor at h versus h. The graphs show the time evolution of the bounding box dimensions.

Figure 5A describes the spatial and evolutionary behavior of cytokines, and Figure 5B shows the evolution of two naive cell types excluded from the tumor micro-environment. Figure 5C shows the spatial distribution of each cell type in the tumor micro-environment next to their maximum, average, and minimum over the whole domain with respect to time.

Figure 5

Results with no flux of immune cells. (A) Column 1: Spatial distribution of cytokines. Column 2: Maximum, average, and minimum concentration (ng/mL) of each cytokine over the whole domain with respect to time. (B) Evolution of naive T cells and naive macrophages. (C) Column 1: Spatial distribution of cell types. Column 2: Maximum, average, and minimum number of each cell type over the whole domain with respect to time.

This result shows that most of the cell types deplete in time except cancer, necrotic and naive dendritic cells (Figure 5). The qualitative behavior of most of the results is comparable with the ODE case, especially mouse 1; see Figure A1 in Appendix A. The most significant difference is the behavior of naive macrophages. The naive macrophages have a strictly decreasing dynamic in the ODE paper, unlike here. The reason might be connected to the slower depletion of cytokines IL10, IL12, and helper T cells, which are the main contributors to the inhibition of naive macrophages. Additionally, HMGB1 has a sharp increase, followed by a decrease in our ODE results. However, the decrease is not monotonic, and in mouse 2, it stops for a short while. Additionally, the timescale of this simulation is much smaller than the ODE paper. To summarize, naive macrophages and HMGB1 are the only fundamentally different results from the behaviors observed in our ODE paper. For the rest of the variables, the changes in the total populations/densities are similar to their dynamics observed in the ODE model. Biologically, the observed decreasing behavior of the immune cells is attributed to cancer cells’ ability to evade identification and invasion by host immune responses in later stages of cancer [107,108,109].

Figure A1

Cells and molecule dynamics from the ODE model [75].

An ODE model fails to capture the significance of the spatial distribution of cells and cytokines. To see how that affects cells and molecules, we studied their level-curve plots. Figure 6 show the contours for cells (left) and molecules (right) separating the regions with values above and under the average of their corresponding field at . Colored areas have the highest number of intersections within the plot. Comparing Figure 5C and Figure 6, we can see that the cancer cells have a minor presence in the area with the most overall immune cells (region ). There is also an intersection between regulatory and cytotoxic T cells with macrophages in the region, which might be the reason for the slightly decreased cancer population in the corresponding location. The intersection between molecules (region ) happens too close to the shaded region to leave a discernible spatial footprint.

Figure 6

(Left): Level-curves indicating the mean value of each cell type at . (Right): Level-curves indicating the mean value of each molecule at . Areas and correspond to the regions with the most and second-most immune cell intersections, respectively. Area corresponds to the region with the highest cytokine intersections.

We compared the simulation results with the experiments. The late foldout of the mouse model used in this study is given in Figure 7. The late foldout is extracted under a naive regime as well, i.e., no treatment has been applied to the mouse model (see Section 2.3.1). Therefore, immune recruitment and infiltration are negligible. Even though it is impossible to show exactly where our simulation domain is located due to the lack of a common reference frame, we can see a good agreement in the overall qualitative behavior of the cell populations. We can see a significant increase in the number of cancer cells compared to Figure 2, while most of the other cell types are depleted. Adipocytes are settled at lower levels, just like the mathematical model’s results. The major differences between the mathematical model’s results and the late foldout are the naive dendritic and necrotic cells. Because of the simplicity and lack of data, the mathematical model does not include metabolites and nutrients, which play a crucial role in necrosis. Here, the only source of necrosis is the death of cancer cells, so they inevitably follow the same trend as cancer cells. Additionally, since HMGB1 acts as a major inhibitor of naive dendritic cells, the mouse model may have a higher level of this molecule.

Figure 7

Nine figures showing each cell type in the mouse model. Blue dots represent a single cell of the corresponding cell type, and gray dots represent the rest.

3.2. Immune Cell Influx

In this section, we consider the influx rates and non-dimensional influx sources for the helper, cytotoxic and regulatory T cells, naive dendritic cells, and macrophages. We assume no-flux boundary conditions for other cells and molecules. Since activated dendritic cells are differentiated from the naive ones, and this activation happens mostly inside of the tumor micro-environment [110], we assume a no-flux boundary condition for the activated dendritic cells. Figure 8 shows the domain after h. Compared to Figure 4, the domain is smaller.

Figure 8

Dimensions of the tumor subject to immune cells influx at h. The curves show the time evolution of the bounding box dimensions.

Figure 9 shows the spatial and evolutionary behavior of cytokines and cells. Naive T cells and macrophages do not change much since they are not spatially dependent. However, due to the Robin type boundary condition, we see either a stationary maximum value (, , and M) or a stationary minimum value () close to the boundary. The former is because the field values are depleting across the domain and the influx tends to bring it up, and the latter is precisely the opposite. The distribution of dendritic cells follows the same trend as the naive dendritic cells. However, it decreases quickly, just like the no-influx case. Due to the near boundary focus of the cells, the molecules are more intense closer to the boundary. As for cancer cells, it seems that their distribution inside of the region and away from the boundary is similar to what Figure 5 shows but at much lower values. On the other hand, there is a higher intensity close to the boundary. Note that this high-intensity region covers a very tiny area, and in total, there are fewer cancer cells present in the region due to the immune cell influx. This is more evident when we compare the integral of cancer cells over the domain for the two cases (Figure 10).

Figure 9

Results with flux of immune cells. (A) Column 1: Spatial distribution of cytokines. Column 2: Maximum, average, and minimum concentration (ng/mL) of each cytokine over the whole domain with respect to time. (B) Evolution of naive T cells and naive macrophages. (C) Column 1: Spatial distribution of cell types. Column 2: Maximum, average, and minimum number of each cell type over the whole domain with respect to time.

Figure 10

Evolution of the integral of cancer over the domain with and without immune cell influx.

In the model, is considered the primary inhibitor of the cancer cells, and  and A are the main contributors to their production. Since is more intense close to the boundary and and A have the same behavior as in the previous case, one might expect to see fewer cancer cells at the boundary. We hypothesize that, even though reactions are decisive in the model, this phenomenon is more because of the cells’ advection direction in the presence of immune cell influx. In other words, cancer cells tend to leave the region. This can be a dangerous trait, given that usually mitotic regions are close to the boundary of the tumor [106,111]. This might explain why some types of breast cancers still metastasize despite therapy and significant immune recruitment.

Similar to the no-influx case, the necrotic cells follow the same pattern as cancer cells. Finally, adipocytes show precisely the same behavior as before. Even though they influence other cells and molecules, their behavior is independent of the other variables. Thus, we observe the same behavior.

3.3. Sensitivity Analysis

To investigate the impact of parameters on the cancer population, we performed a sensitivity analysis on all the model parameters for the case with immune cell influx. We applied an adjoint sensitivity analysis for the functional

We used this method because the model contains a large set of parameters, and the computational cost of adjoint-based sensitivity analysis is almost independent of the number of input parameters. Hence, we can compute sensitivities with respect to numerous parameters or even entire functions. For more information on the mathematical detail of the method, please see [112]. To carry out the sensitivity analysis using adjoint methods, we took advantage of the dolfin-adjoint package [113]. The benefit of this package is that it can be mounted on FEniCS and can record each step of the simulation. Once the adjoint-based sensitivity starts, it automatically steps backward and calculates the sensitivities. Figure 11 shows the top five sensitivity values of (16) to four category of parameters: diffusion rates, influx rates, influx sources, and the rest of the reaction parameters. For a full report of the sensitivity analysis, parameter notations, and definitions, see Appendix E and Table A1.

Figure 11

The sensitivity of at to four categories of parameters: diffusion rates, influx rates, influx sources, and the reaction parameters.

Table A1

Full sensitivity report. This table contains the sensitivity value of to all the parameters used in the model. The rows are ordered in a decreasing fashion based on the absolute value of their sensitivity values.

Order	Notation	Sensitivity	Definition	Order	Notation	Sensitivity	Definition
1	D_C	2242.835	Diffusion coefficient of C	44	δ_T_h	−2.87×10^−5	Death rate of T_h
2	D_A	56.46447	Diffusion coefficient of A	45	λ_MT_h	2.10×10^−5	Activation rate of M by T_h
3	δ_C	−3.16066	Death rate of C	46	A_D_N	−1.48×10^−5	Independent production rate of D_N
4	λ_C	2.989047	Proliferation rate of C	47	δ_T_N	1.41×10^−5	Death rate of T_N
5	α_T_C	−0.31011	Influx rate of T_C	48	λ_T_hH	1.03×10^−5	Activation rate of T_h by H
6	D_T_C	−0.2628	Diffusion coefficient of T_C	49	λ_DC	−6.42×10^−6	Activation rate of D by C
7	λ_CIL_6	0.138431	Proliferation rate of C by IL_6	50	λ_T_hIL_12	6.27×10^−6	Activation rate of T_h by IL_12
8	λ_CA	0.107539	Proliferation rate of C by A	51	A_M_N	3.50×10^−6	Independent production rate of M_N
9	δ_CT_C	−0.06665	Inhibition rate of C by T_C	52	λ_T_CD	−3.46×10^−6	Activation rate of T_C by D
10	δ_A	−0.03611	Death rate of A	53	δ_D_N	3.26×10^−6	Death rate of D_N
11	λ_A	0.03531	Proliferation rate of A	54	λ_IL_10D	3.14×10^−6	Production rate of IL_10 by D
12	α_M	0.025664	Influx rate of M	55	δ_T_hIL_10	−1.71×10^−6	Inhibition rate of T_h by IL_10
13	δ_IL_6	−0.01152	Decay rate of IL_6	56	δ_T_hT_r	−1.46×10^−6	Inhibition rate of T_h by T_r
14	λ_IL_6A	0.009239	Production rate of IL_6 by A	57	λ_DH	−1.37×10^−6	Activation rate of D by H
15	λ_IL_6M	0.007813	Production rate of IL_6 by M	58	A_0	9.75×10^−7	Carrying capacity of A
16	α_T_h	0.002154	Influx rate of T_h	59	δ_M_N	−3.58×10^−7	Death rate of M_N
17	δ_T_C	0.001508	Death rate of T_C	60	λ_IL_12M	−3.55×10^−7	Production rate of IL_12 by M
18	α_T_r	0.001441	Influx rate of T_r	61	δ_IL_12	2.87×10^−7	Decay rate of IL_12
19	D_D_N	−0.00094	Diffusion coefficient of D_N	62	λ_IL_12T_h	−2.31×10^−7	Production rate of IL_12 by T_h
20	D_N	−0.00094	Diffusion coefficient of N	63	λ_IL_12T_C	−2.18×10^−7	Production rate of IL_12 by T_C
21	C_0	0.000541	Carrying capacity of C	64	D_IL_12	9.50×106−8	Diffusion coefficient of IL_12
22	D_T_r	−0.00044	Diffusion coefficient of T_r	65	λ_T_hD	7.78×10^−8	Activation rate of T_h by D
23	δ_IL_10	−0.00039	Decay rate of IL_10	66	λ_T_rD	6.93×10^−8	Activation rate of T_r by D
24	λ_IL_10C	0.000386	Production rate of IL_10 by C	67	δ_D	−3.70×10^−8	Death rate of D
25	λ_T_CIL_12	−0.00035	Activation rate of T_C by IL_12	68	δ_H	−1.49×10^−8	Decay rate of H
26	D_M	−0.00034	Diffusion coefficient of M	69	λ_HM	8.05×10^−9	Production rate of H by M
27	λ_IL_10T_r	0.000319	Production rate of IL_10 by T_r	70	λ_HT_C	5.35×10^−9	Production rate of H by T_C
28	δ_M	−0.00029	Death rate of M	71	λ_HC	5.00×10^−9	Production rate of H by C
29	λ_IL_10M	0.00024	Production rate of IL_10 by M	72	δ_DC	3.46×10^−9	Activation rate of D by C
30	λ_IL_10T_h	0.00018	Production rate of IL_10 by T_h	73	λ_IL_12D	−2.50×10^−9	Production rate of IL_12 by D
31	λ_IL_10T_C	0.000153	Production rate of IL10 by T_C	74	λ_HN	4.23×10^−10	Production rate of H by N
32	D_T_h	−0.00015	Diffusion coefficient of T_h	75	λ_HD	1.96×10^−10	Production rate of H by D
33	D_IL_6	0.000126	Diffusion coefficient of IL_6	76	δ_N	1.14×10^−11	Death rate of N
34	A_T_N	−9.73×10^−5	Independent production rate of T_N	77	D_H	1.90×10^−13	Diffusion coefficient of H
35	α_D_N	9.34×10^−5	Influx rate of D_N	78	α_NC	8.36×10^−15	C to N conversion fraction
36	δ_T_CT_r	8.92×10^−5	Inhibition rate of T_C by T_r	79	TC*	−3.56×10^−16	T_C influx source
37	D_D	8.69×10^−5	Diffusion coefficient of D	80	M^*	3.44×10^−17	M influx source
38	δ_T_CIL_10	7.88×10^−5	Inhibition rate of T_C by IL_10	81	Th*	2.87×10^−18	T_h influx source
39	λ_IL_6D	7.41×10^−5	Production rate of IL_6 by D	82	Tr*	2.47×10^−18	T_r influx source
40	D_IL_10	−5.72×10^−5	Diffusion coefficient of IL_10	83	DN*	−8.55×10^−20	D_N influx source
41	λ_MIL_10	5.35×10^−5	Activation rate of M by IL_10
42	δ_T_r	−4.66×10^−5	Death rate of T_r
43	λ_MIL_12	3.41×10^−5	Activation rate of M by IL_12

Interestingly, the top four most sensitive reaction parameters were captured in the same order in the ODE sensitivity analysis [75]. However, despite acknowledging as a sensitive parameter, it was not this high up in the ODE study. This is interesting because -related PDE parameters also show significant sensitivities, which are negative. This is due to the model’s higher levels of after imposing an influx for this cell type.

Moreover, the diffusion rate of cancer cells is one of the most sensitive parameters. In other words, the more motile the cancer cells get, the larger their population becomes. This is because cancer cells will interact more with cells and molecules that promote their proliferation. The diffusion rate of adipocytes is the following most sensitive parameter in this model. Adipocytes directly activate cancer cells, and increasing their motility means more cancer/adipocyte handshakes.

The macrophage-related parameter, , has a relatively high positive sensitivity value. This means that a greater influx of macrophages leads to more cancer cells. This is due to the fact that macrophages produce and , which, respectively, cause cancer cell proliferation and inhibition of cytotoxic T cells. Traditionally, macrophages are the immune system soldiers in charge of clearing target cells. However, many studies have investigated the tumor-promoting capabilities of macrophages [114,115,116,117].

Figure 12A shows the result of 10% perturbation of the top 20 most sensitive parameters given in Table A1. We can see a significantly lower number of cancer cells at the lower 5%. Moreover, Figure 12B shows a drastic change in the size of the tumor when the most sensitive parameters are varied.

Figure 12

(A) The variation of the total number of cancer cells as a result of 10% perturbation of the most sensitive parameters. (B) The leftmost circle corresponds to the lower bound, the middle circle corresponds to the thick solid red curve and the right circle corresponds to the upper bound of the graph in (A).

Since the diffusion of cancer cells and adipocytes are the most sensitive parameters, the results indicate targeting cancer cells motility might be a treatment strategy, as suggested in some other studies [118,119]. Parameters number 3 and 4 from Table A1 are not informative, because they engulf production and death caused by reasons other than the ones we have included in the model. Parameters number 5, 6, and 9 emphasize the importance of cytotoxic T cells in cancer inhibition. Targeting CD8+ for immunotherapy has been discussed extensively in the literature [120,121,122]. As for parameters 7 and 8, we refer to the ODE paper for a detailed discussion about the importance of adipocytes and IL6 in cancer therapy [75].

4. Discussion

In this study, we investigated the spatial interaction network of key cells and molecules in the breast cancer tumors of a PyMT mouse model by developing a bio-mechanical system of PDEs. We adopted the critical reactions among cells and molecules from an ODE model of mice breast tumors [75], and the initial conditions were extracted from an MMTV-PyMT mouse model.

Since there was no treatment applied in the experimental mouse model, the recruitment of the immune cells was negligible. We, therefore, first investigated a case with no immune cell influx. We noticed that the domain grows and deforms into a larger blob-shaped region, a shape which is commonly observed in experiments [104,105,106].

The spatial distribution of the model state variables, which are in good agreement with the late foldout experimental results, shows that the regions with the higher number of immune cells have much fewer cancer and necrotic cells. The only differences between the late foldout cells distributions in the mathematical model and experimental results are the number of necrotic cells and naive dendritic cells. We hypothesize that the former is due to the lack of nutrients and metabolites in the mathematical model, which affects necrosis [123,124]. The latter is possibly due to the level of HMGB1 (a major activator of naive DCs) in the mouse model compared to the mathematical model.

Moreover, we modeled an influx for immune cells through Robin boundary conditions. The results indicate a lower growth rate of the domain with influx compared to the no-influx case, and there was a significant overall decrease in the number of cancer cells in the domain. As a result of these influxes, most cells and cytokines were focused near the boundary of the tumor. Interestingly, despite the significant presence of the immune cells near the boundary, cancer cells have a higher concentration there. We hypothesize that the interactions resulting from the immune cells’ influxes create an outward divergence for the velocity field. This drives many of the cancer cells to the boundary of the domain. Since mitotic regions are usually close to the tumor boundary [106,111], this might explain why tumor cells can still escape after immunotherapy.

We calculated the sensitivity of total cancer cells to all the model parameters. The most sensitive reaction parameters are in agreement with the ODE study. Importantly, we observe a significant sensitivity to the diffusion coefficient of cancer cells and adipocytes. The interaction between cancer cells and adipocytes promotes cancer proliferation, and increasing these values elevates their interaction, leading to more cancer cells.

The 10% perturbation of the top 20 most sensitive parameters shows a small cancer growth at the lower bound, primarily due to the diffusion rates of cancer cells and adipocytes. Therefore, controlling the motility of these cells can lead to better prognoses, and there are already studies targeting cancer cells’ motility [118,119]. The results also indicate that the influx rate and source of cytotoxic T cells have a high impact on the total cancer cells. In other words, increased influx of these cells leads to fewer cancer cells, which is consistent with clinical and experimental observations [120,121,122]. Another important observation is the positive sensitivity values for the influx rate and source of macrophages. This is because macrophages produce IL6 and IL10, promoting cancer proliferation and inhibiting cytotoxic T cell proliferation. This result is also in line with biological and biomedical findings on the tumor-promoting properties of macrophages [114,115,116].

The findings of this study have to be seen in light of some limitations. One of the main challenges of mathematical modeling of cancer is the lack of data for more reliable parameter estimation and validation. Although we used the reaction parameters from the earlier ODE study, we did not have access to time-course data for this model to obtain the PDE-specific parameters. Therefore, we had to refer to studies that had estimated these parameters based on certain assumptions. Moreover, the cell-cell adhesion force value for the mechanical problem was a crude estimate. Currently, there is no direct report of such values in the literature to the best of the authors’ knowledge. In addition, nutrients and metabolites can significantly affect the tumor shape, size and number of cells. In avascular tumor models, the type of nutrients and the competition for their consumption between normal cells and cancer cells will lead to different growth pattern tendencies: circular or papillary-like [125]. On the other hand, depletion of nutrients and metabolites in the tumor microenvironment can cause the immune cells to lose their functionality, which can lead to cancer cells’ growth and immune invasion behavior [126]. However, to reduce the complexity of the model, we did not consider metabolites such as oxygen and nutrients.

Nevertheless, the model unveils some spatial features of breast tumor growth and identifies the most sensitive parameters, including diffusion rates of cancer cells and adipocytes to control the tumor growth. With access to more initial spatial data, this model can make predictions about the effect of early immune cells’ infiltration patterns on cancer progression. Additionally, future research can build upon this model to overcome its limitations, such as by including more biological processes or integrating different treatment options.