A Mathematical Model of Chemotaxis in Cancer Treatment

Belgacem Al-Azem, Ali Nadim

Institute of Mathematical Sciences, Claremont Graduate University, Claremont, CA, USA.

DOI: 10.4236/am.2025.1611041                

Abstract

Chemotaxis is the directed movement of cells or organisms in response to a chemical gradient. This process is critical in many biological systems. In the context of cancer, chemotaxis plays a significant role in how cancer cells migrate and invade surrounding tissues. Our goal is to design a mathematical model that utilizes chemo-attractants to draw cancer cells towards the core of the tumor domain. We include therapeutic agents to eradicate cancer cells concentrated at the core of the tumor domain. To that aim, we set up a coupled system of partial differential equations that describes the interactions and dynamics of the following components: cancer cells concentration, chemo-attractants concentration, and cancer killer agents concentration. The model components are designed to manipulate the spatial distribution and eradication of cancer cells within the tumor microenvironment over time. Each term in the equations of those components plays a specific role to that end. Our simulations confirm that the combined therapy is highly effective in reducing tumor burden. Counter-intuitively, we find that chemotactic drift alone can achieve up to a 78% reduction in tumor mass within one unit of dimensionless time by inducing overcrowding at the tumor core, while the addition of therapy further enhances eradication and accelerates tumor suppression.

Keywords

Chemotaxis, Cancer Cell Migration, Mathematical Modeling, Partial Differential Equations (PDEs), Tumor Microenvironment

1. Introduction

Chemotaxis is a critical process in many biological systems, including immune responses, wound healing, and cancer progression [1]-[5]. It directs cells to move in a desired direction. The movement can be of two types: Positive chemotaxis, where cells move towards higher concentrations of specific chemical attractants; these are substances that attract cancer cells. They can include growth factors, cytokines, or other molecules that promote cancer cell survival and proliferation [6] [7]. For instance, immune cells move towards sites of infection or inflammation where signaling molecules are released [8] [9].

The other type is negative chemotaxis which refers to the process by which cells detect and move away from high concentrations of repulsive chemical signals. In cancer biology, these are substances that repel cancer cells. They can include certain therapeutic agents or signals that induce cell death or inhibit cell migration [10]. Generally, negative chemotaxis can help cells avoid harmful substances or unfavorable environments. Although the detailed mechanisms and signaling pathways remain under active investigation, several studies over the past few years have explored how these chemo-repulsive agents modulate cell behavior [11]-[13].

Since our purpose here is the mathematical modeling of chemotaxis, it is useful to highlight some key concepts that can help our understanding of the model application: Chemotactic Gradient is the spatial variation in the concentration of a chemical substance that cells detect, leading them to move in a specific direction. Chemotactic Receptors are proteins on the cell surface that bind to chemoattractants or chemo-repellents, initiating intracellular signaling pathways that guide cell movement [11]. Cell Polarization refers to the establishment of spatial asymmetry in a cell. This asymmetry is essential for effective cell migration because it organizes cellular components into distinct front (leading edge) and rear (trailing edge) domains. The process involves the reorganization of the cyto-skeleton along with the redistribution of organelles, membrane receptors, and signaling molecules, thereby enabling the cell to interpret and respond to directional cues from chemical gradients [14]-[16]. This coordination is critical not only for physiological processes like wound healing and immune responses but also plays a significant role in pathological contexts such as cancer metastasis.

In the context of cancer, chemotaxis plays a significant role in how cancer cells migrate and invade surrounding tissues. Negative chemotaxis can be exploited to counteract metastasis by repelling tumor cells, outlining potential strategies for designing anti-metastatic agents [17]. Therefore, understanding and manipulating chemotaxis can constitute a crucial tool for developing effective cancer treatments. For example, by designing therapies that alter chemotactic gradients, researchers may be able to potentially control the movement of cancer cells, concentrating them in specific local areas of the tumor where they can be more effectively targeted by therapeutic agents. We are not biologists, so we do not make any feasibility or validity claim. Nevertheless, the following possible specific manipulations were the basis of our mathematical model and, of course, we would be happy if those ideas find some validity and resonance in the field of cancer treatment:

Chemo-attractants and Chemo-repellents: Recent developments in repelling cancer cells and impeding tumor spread through molecular inhibition [18] highlight potential therapeutic approaches that target negative chemotaxis. By administering or inhibiting certain chemicals, the movement of cancer cells can be directed. For example, increasing the concentration of a chemo-attractant in the tumor core can pull cancer cells away from critical areas.

Therapeutic Agents: The balance between chemo-attractive and chemo-repulsive signals in the tumor microenvironment has implications for new treatment modalities that aim to tip this balance against cancer progression [19]-[23]. Combining chemo-tactic manipulation with therapeutic agents can enhance treatment efficacy. Agents can be designed to be more effective in regions where cancer cells are concentrated due to chemo-tactic manipulation.

Micro-environment Modulation: Understanding the cellular signaling pathways that mediate chemorepulsion in cancer cells offers insights into how modulation of these pathways can influence cancer cell migration and dissemination [24] [25]. Therefore, altering the tumor microenvironment to change chemo-tactic gradients can help control cancer spread and improve the response to therapy.

Overall, we believe that chemotaxis is a fundamental process that could be strategically manipulated to improve cancer treatment outcomes without damaging unaffected or critically sensitive non-targeted parts of tissues.

2. The Mathematical Model

We introduce a comprehensive mathematical model for the evolution of cancerous tumors accounting for the effects of chemo-attractants and therapeutic agents. The model focuses on the growth and migration of cancer cells, as well as on their death due to the presence of therapeutic agents. In our model, we introduce sources of chemo-attractants and therapeutic agents near the center or core of the tumor to see whether one can cause the live cancer cells to migrate inwards toward the core, where the higher concentration of therapeutic agents can bring about their death.

2.1. Model Definitions

In order to include some effects of geometry of the tumor domain, we choose to work in a cylindrical polar coordinate system with axial symmetry. While tumors are likely more spherical than cylindrical, because in our simulation software (COMSOL) the axisymmetric polar system is already built-in, this choice makes the computations more straightforward.

We thus take the tumor domain Ω to be a disk of radius R (polar coordinates system with 1 − D axi-symmetry as defined by COMSOL Multi-Physics system). Within that disk, we somewhat arbitrarily divide the domain into three subregions: Core, Central, and Peripheral, with the core region nearest the symmetry axis (r = 0), and the peripheral region near to the outer boundary (r = R), and the central region in between those two. To summarize the numerical results later, we compute the mass of the tumor within those three regions to see how the combination of chemotaxis and therapeutics redistributes and potentially eliminates the cancer cells in the tumor.

The dependent variables in our model consist of: C(r,t), the cancer cell concentration; A(r,t), the chemo-attractant concentration; and T(r,t), the therapeutic agent concentration. Here, r is the radial variable in polar coordinates and t is time.

2.2. The Governing Equations

We start with the partial differential equation (PDE) that describes the evolution of the cancer cell concentration C(r,t). Essentially it is a conservation equation of the form:

$\frac{\partial C}{\partial t}+\nabla \cdot J=S$

in which $J$ is the flux of cancer cells and S is their net rate of production within the region. The flux $J=vC-D_c\nabla C$ includes an advective term characterized by velocity v and a diffusive term characterized by diffusivity $D_c$ . In the absence of any bulk flow, the advection velocity is purely due to chemotaxis and would be given by $v={\chi }_a\nabla A$ in the presence of chemo-attractant with concentration A. The source term S includes expressions for growth rate of cancer cells, as well as their removal or death rates (with negative signs). More specifically, the PDE for cancer cell concentration becomes:

$\frac{\partial C}{\partial t}=D_c{\nabla }^{2}C-{\chi }_a\nabla \cdot \left(C\nabla A\right)+\alpha C\left(1-\frac{C}{K}\right)-\beta C-\delta CT$ (1)

The individual terms in this equation can be interpreted as follows. The term $\partial C/\partial t$ is the time rate-of-change of cancer cell concentration resulting from all the transport and growth and removal processes on the right-hand side. The term $D_c{\nabla }^{2}C$ represents the transport of cancer cells by diffusion, with $D_c$ being the diffusion coefficient due to random motility of cancer cells within the tumor microenvironment. Cancer cells exhibit random movement due to interactions with the extracellular matrix and other cells. The term ${\chi }_a\nabla \cdot \left(C\nabla A\right)$ represents transport due to chemotaxis towards higher concentrations of chemo-attractants. Here, ${\chi }_a$ is the chemotactic sensitivity coefficient, A is the chemoattractant concentration, and $\nabla A$ is its gradient. Cancer cells are known to migrate towards higher concentrations of certain chemicals that promote their growth and survival such as growth factors or nutrients.

The next set of terms describe the growth and death rates of cancer cells, comprising the source term in the original conservation equation. The term $\alpha C\left(1-C/K\right)$ is an assumed logistic growth rate for cancer cells; $\alpha$ is the proliferation rate, and K is the carrying capacity of the cancer microenvironment. Cancer cell proliferation follows a logistic growth pattern where growth slows down as the population approaches the carrying capacity of the environment. Proliferation of cancer cells can be limited by the finite availability of space and nutrients. The term $\beta C$ is the intrinsic (natural) cancer cells death rate, where $\beta$ is the death rate coefficient. Cancer cells are subject to finite lifetime. Finally, the term $\delta CT$ is a removal rate that represents the death of cancer cells due to the presence of the therapeutic agent, $\delta$ being the proportionality coefficient of the death rate with therapeutic agent concentration T, such as chemotherapy drugs that kill cancer cells upon contact. Therapeutic agents are designed to specifically target and kill cancer cells, reducing the overall tumor burden.

Clearly Equation (1) shows the coupling between the evolution of cancer cell concentration C with the concentrations of chemo-attractants A and therapeutic agents T. We thus also write transport equations for these two concentrations. Both of these are assumed to be governed by diffusion equations in which we also allow for degradation of the molecules and potentially a distributed source within the tumor domain. The respective PDEs for the chemo-attractant concentration and the therapeutic agent concentration thus read:

$\frac{\partial A}{\partial t}=D_a{\nabla }^{2}A-{\lambda }_{a}A$ (2)

$\frac{\partial T}{\partial t}=D_t{\nabla }^{2}T-{\lambda }_{t}T$ (3)

In these equations the left-hand sides are the rates of change of the respective concentrations and the right-hand sides include diffusive transport terms with diffusivities $D_a$ and $D_t$ , and decay rates with rate constants ${\lambda }_a$ and ${\lambda }_t$ . Unlike general formulations that include bulk source terms $S_a$ and $S_t$ , our model supplies the chemo-attractant and therapeutic agents exclusively through fluxes imposed at the inner boundary of the domain. We therefore omit bulk sources from the governing equations to streamline the presentation.

In summary, Equations (1) - (3) constitute a system of three coupled PDEs that govern the transport and evolution of cancer cell, chemo-attractant and therapeutic agent concentrations within the tumor. The initial and boundary conditions imposed on this system are specified in Section 3 below in the context of the numerical simulations. This system and its associated initial and boundary conditions constitute our complete mathematical model. Each term in the governing equations corresponds to specific physical or biological processes that are crucial for modeling cancer cell behavior in the tumor microenvironment. By incorporating these terms, the model can simulate how cancer cells evolve, proliferate, and respond to treatments and environmental cues, providing a comprehensive framework for understanding and optimizing cancer therapy. In the simulations, we set up the boundary conditions so as to ensure that chemo-attractants draw cancer cells towards the core of the tumor where we expect them to be eradicated by therapeutic agents that are also introduced mainly in that location.

3. Numerical Simulations

To test our main hypothesis, we work on three variants of the model describing tumor evolution. Our working hypothesis is that if chemo-attractants and therapeutics are introduced into the tumor core, they can induce cancer cells to migrate toward the core where targeted therapy can lead to a significant reduction or complete eradication of cancer. The three variants of the model loosely correspond to three of the main stages of the hypothetical tumor. One application aims at capturing the dynamics of a typical stage I tumor, a second attempts to describe the evolution of a typical stage II tumor, and a third describes the dynamics of a typical stage III tumor. Note that we are excluding stage zero and metastatic tumor (stage IV). Before presenting and analyzing the simulations we have performed, we briefly describe the different stages of cancer tumor.

Cancer staging provides a standardized method to describe the extent of disease progression. It is crucial for determining prognosis and guiding treatment decisions. One of the most widely used systems is the TNM system, which classifies cancer based on three key components: T (Tumor) refers to the size and/or local extent of the primary tumor, N (Node) refers to the involvement of regional lymph nodes, M (Metastasis), refers to the presence of distant spread [5] [4] [11] [23]. This staging system is continuously refined based on accumulating clinical data and evolving understanding of tumor biology, making it a pivotal tool in both clinical practice and research.

A stage I tumor is small and limited to its site of origin, without spreading to lymph nodes or distant sites. It is easier to treat and it is often curable. Left unchecked, a stage I tumor can evolve to stage II, a larger tumor that may spread to nearby tissues or lymph nodes but not to distant organs. This stage requires a more extensive treatment compared to stage I. A stage II tumor, in turn, can evolve to stage III and become more aggressive and extensively invasive of a larger number of regional lymph nodes. At this stage, there is still no distant metastasis, but the disease is more aggressive and can evolve to the next stage. At stage IV, the cancer can spread to distant parts of the body (e.g., lungs, liver, bones), and becomes much harder to treat and often considered incurable but manageable with palliative care or systemic therapies. Understanding these stages is crucial for tailoring treatment and predicting outcomes.

In our simulations, while we do not account for spread of cancer to nearby regions, we distinguish between stages I through III using different initial conditions corresponding to the initial extent of cancer cells within the domain of interest, with stage I starting with a small cancerous region, while for stage III, the cancer cells initially fill the entire computational domain.

We emphasize that this use of “Stage I,” “Stage II” and “Stage III” is a modeling simplification based solely on initial spatial distributions. It does not correspond directly to the clinical TNM staging system, which incorporates tumor size, lymph node involvement, and metastasis.

3.1. Dimensionless Equations

For computational and visualization purposes, we now proceed to non-dimensionalize the mathematical model described in Equations (1) - (3). To that end, we choose as our length scale the radius R of the disk-like domain and for our time scale, the diffusion time of cancer cells over that distance: ${R^2/D}_c$ . We thus define a dimensionless radial coordinate $\widehat{r}=r/R$ and a dimensionless time $\widehat{t}=D_{c}t/R^2$ , and we scale all concentrations with the carrying capacity K, taking $\widehat{C}=C/K$ , $\widehat{A}=A/K$ , and $\widehat{T}=T/K$ . Upon dropping the hats for notational clarity, we obtain:

$\frac{\partial C}{\partial t}={\nabla }^{2}C-{\beta }_1\nabla \cdot \left(C\nabla A\right)+{\beta }_{2}C\left(1-C\right)-{\beta }_{3}C-{\beta }_{4}CT$ (4)

$\frac{\partial C}{\partial t}={\beta }_5{\nabla }^{2}A-{\beta }_{6}A$ (5)

$\frac{\partial C}{\partial t}={\beta }_7{\nabla }^{2}T-{\beta }_{8}T$ (6)

Where the eight dimensionless parameters appearing therein are given by

${\beta }_1=\frac{{\chi }_{a}K}{D_c},{\beta }_2=\frac{\alpha R^2}{D_c},{\beta }_3=\frac{\beta R^2}{D_c},{\beta }_4=\frac{\delta R^{2}K}{D_c}$ (7)

${\beta }_5=\frac{D_a}{D_c},{\beta }_6=\frac{\lambda \alpha R^2}{D_c},{\beta }_7=\frac{D_t}{D_c},{\beta }_8=\frac{{\lambda }_T{R}^2}{D_c}$ (8)

Note that we model the introduction of attractants and therapeutic agents as a boundary flux, corresponding to localized delivery methods such as intratumoral injection or infusion devices. Instead of having distributed sources of attractants and therapeutic agents, we introduce them into the tumor through Neumann boundary conditions.

3.2. Notes on COMSOL Multiphysics Software

The tumor domain Ω that we have in mind is a sphere with radial symmetry. However, for convenience, we model a disk-shaped tumor in polar coordinates with axisymmetry. This means that the properties of the tumor are invariant under rotation about a central axis, with no dependence on the angular coordinate. Thus all quantities depend only on the radial coordinate r. In COMSOL this geometry is built-in and referred to as “1-D axisymmetric”.

At this point, an explanation of a small subtlety in COMSOL formulation of the so called 1-D axisymmetric geometry is in order: In many finite element (and other numerical) simulations involving cylindrical coordinates, the origin (r = 0) is a singular point. In COMSOL Multiphysics that point is not treated as a standard boundary but rather as an interior point. This modeling strategy avoids the issues associated with applying a boundary condition at a singular point. Instead, the underlying formulation takes advantage of the invariance of the physical laws under rotation about the axis, ensuring that the appropriate conditions are implicitly satisfied at the origin without requiring an explicit boundary condition [26]-[28]. In our case, however, we would like to impose flux-type boundary conditions at the core of the tumor domain. To get around this limitation, we modify our domain from the dimensionless range 0 ≤ r < 1 to a new one: 0.02/R < r < 1, starting away from the axis of symmetry. We can now impose Neumann boundary conditions at the point r = 0.02/R, to specify the fluxes of all three components in our model. In our simulations the numerical value of R was taken to be 2.56 (in centimeters), and the value 0.02 cm corresponds to 200 microns.

Using COMSOL Multiphysics software, we simulated several instances of the model by varying two key parameters: The Drift Coefficient β₁ and the Therapeutic Agent Coefficient β₄ of the cancer cells PDE (4). We expected that these particular parameters would play a pivotal role in determining the qualitative and numerical evolution of the body mass of cancer cells within the tumor microenvironment. Indeed, our overarching hypothesis is that the combined effect of convecting cancer cells towards the core of the tumor domain, coupled with targeted therapy in that same region, would dramatically reduce, if not eliminate, the tumor from the entire microenvironment.

After many trials, we settled on {0, 30, 60, 90} as the set of possible values for the drift coefficients β₁ and {0, 10, 20, 30, 40} as that of the therapeutic coefficient β₄. Against each of the 5 values in the therapeutic coefficient set, we ran four instances that correspond to the four values in the drift coefficients set, resulting in 20 total completed numerical simulations. Other special cases or variants of simulations were also performed for the purpose of further illustrations and validations of the results. For each instance of the above simulations, we used the built-in integration capability of COMSOL to compute the non-dimensional body mass of cancer cells in the three regions of the tumor: The core region (0.02/R < r < 0.52/R) , the central region (0.52/R < r < 0.54/R), and the peripheral region (0.52/R < r < 2.56/R). In polar coordinates, these integrals have the forms:

$\text{mass}=2\pi {\displaystyle {\int }_{ri}^{rj}rC\left(r,t\right)}dr.$

We also computed the total mass in the entire domain separately. We performed this process for each of the 3 stages. The results are tabulated and visualized below. Oncological staging is implemented by assigning differing initial values to the cancer concentration C for each stage.

3.3. Initial and Boundary Conditions

We performed numerical simulations on the evolution of the three concentration fields $C\left(r,t\right)$ , $A\left(r,t\right)$ , and $T\left(r,t\right)$ starting with three different initial conditions that were meant to represent three of the stages of cancer at diagnosis. For all three cases, we took the initial conditions on the attractant and therapeutic to be

$A\left(r,0\right)=0,T\left(r,0\right)=0$

which means that they were not initially present in the system. For the cancer cell concentration, however, we took the initial conditions to be one of:

$\text{Stage I}:C\left(r,0\right)=\{\begin{array}{l}0.1\left(0.02/R<r<0.52/R\right)\\ 0\text{ elsewhere}\end{array}$

$\text{Stage II}:C\left(r,0\right)=\{\begin{array}{l}1\left(0.02/R<r<0.52/R\right)\\ 0.1\left(0.52/R\le r<2.54/R\right)\\ 0\text{ elsewhere}\end{array}$

$\text{Stage III}:C\left(r,0\right)=1\text{ everywhere}$

Stage III thus represents a case where the cancer cells have reached the carrying capacity of unity throughout the domain at the point of diagnosis. Stage I starts with the cancer cells at one-tenth of the carrying capacity in a small region near the center. In stage II, we have cancer at the carrying capacity near the center and one-tenth of it in the central region, and zero in the peripheral region. Recall that the numerical value of R is 2.56 in our computations.

As far as the boundary conditions are concerned, in all cases, on the cancer cell concentration C, we impose no-flux boundary conditions at both ends of the interval. For the other two variables A and T, we impose no-flux conditions at the periphery (r = 1), but at the core (r = 0.02/R) we impose an incoming flux into the tumor domain, with a dimensionless value of unity. The boundary conditions are thus:

$\frac{\partial C}{\partial r}=0,\frac{\partial A}{\partial r}=0,\frac{\partial T}{\partial r}=0,\text{ at }r=1$

and

$\frac{\partial C}{\partial r}=0,\frac{\partial A}{\partial r}=-1,\frac{\partial T}{\partial r}=-1,\text{ at }r=0.02/R$

for all time t > 0.

This modeling choice reflects a localized delivery mechanism, such as direct intratumoral injection into the tumor core or infusion through an implanted device. Intratumoral injection has been investigated in experimental oncology as a means to achieve high local concentrations of therapeutic agents while minimizing systemic exposure [29] [30]. Representing the introduction of chemo-attractants and therapeutics as boundary fluxes thus captures the effect of sustained, spatially focused administration without assuming a uniform distribution throughout the tumor.

4. Results

We have simulated stage I, II, and III tumor evolutions according to the model described above using COMSOL. All three models reach the same final steady state for the sets of parameters we examined. Their transient dynamics are different because they start with different initial conditions. To avoid redundancy, we simply highlight the major differences among the 3 stages during their transient states, and focus the steady-state analysis on stage III only. If the model shows significant eradication of stage III tumor, our worst-case scenario, we surmise that it will succeed in the less acute cases of stage I and II (Figures 1-6).

Figure 1. The figure captures the tumor evolution during the transit state of stage 1, where the initial values of 10⁻¹ at the core and zero elsewhere still have effect. Initially the tumor grows because the carrying capacity of cancer cells is small, so the logistic term of (4) leads to initial growth. The graph shows that after about 5 dimensionless units of time the tumor mass declines significantly everywhere.

Figure 2. The figure captures the tumor evolution during the transit state where the initial values of 1 at the core, 10⁻¹ at the central region, and zero near the boundary still play a noticeable role. At the core region where the carrying capacity of the logistic term of equation (4) is maximum, the cancer mass decreases at all times. In the central and the peripheral region where the carrying capacity is small, the logistic term leads to initial growth. The graph shows that after about 5 dimensionless units of time, the tumor mass declines significantly everywhere.

Figure 3. The figure captures the tumor evolution during the transit state where the initial value is 1 everywhere. Since the carrying capacity of the logistic term of equation (4) is maximum, the cancer mass decreases everywhere and at all times. The graph shows that after about 5 dimensionless units of time, the tumor mass declines significantly everywhere.

Figure 4. The reduction in total tumor body mass in time for the Therapeutics coefficient T = 0.

Figure 5. The reduction in total tumor body mass in time for the Therapeutics coefficient T = 20.

Figure 6. The reduction in total tumor body mass in time for the Therapeutics coefficient T = 40.

Summary: The numerical simulations show that the only difference between Stage I, Stage II, and Stage III implementations is during the transit states of their respective evolutions.

During stage 1 transit state, the tumor evolves under the effect of the initial values of 10⁻¹ at the core and zero elsewhere. Initially, the tumor grows because the carrying capacity of cancer cells is small, so the logistic term of Equation (4) leads to initial growth.

During stage 2 transit state, the tumor evolves under the effect of the initial values of 1 at the core, 10⁻¹ at the central region, and zero elsewhere. At the core region where the carrying capacity of the logistic term of Equation (5) is maximum, the cancer mass decreases at all times. In the central and the peripheral region where the carrying capacity is small, the logistic term leads to initial growth.

During stage 3 transit state, the tumor evolves under the influence of the initial values of 1 in the entire tumor domain. Since the carrying capacity of the logistic term of Equation (6) is maximum, the cancer mass decreases everywhere and at all times. All the 3 graphs above show that after about 5 dimensionless units of time, the tumor mass declines significantly and everywhere.

Below we present tables and graphs that convey the results of our simulations. We only included a sample from the set of visualizations we performed to avoid redundancy and clutter. In order to further help the reader navigate the numerical results in the tables, we state the following guidelines:

1) To see the effect of the Drift, compare the values within each table as you scroll down any column in that table (Tables 1-6).

For example, if we want to see the effect of Drift on the Central region of the tumor domain, when the Therapeutic coefficient is 10, we scroll down the Central Mass column of Table 2.

Table 1. Stage III tumor with a Therapeutic Coefficient of 0.

Stage 3 Cancer Tumor Evolution with Therapeutics Coef. = 0
Drift Coef.	Core Mass	Central Mass	Peripheral Mass	Total Mass
0	0.1118	12.9160	0.3102	13.3380
30	0.0875	09.4640	0.2296	09.7811
60	0.0510	06.1093	0.1492	06.3095
90	0.0226	02.8087	0.0689	02.9002
Rate*	20.2%	21.7%	22.2%	21.7%

Note: *The “Rate” row shows the % of tumor mass remaining after applying a ’maximum’ drift coefficient (β₁ = 90) for each region at the therapeutics coefficient level (β₄ = 0) compared to the baseline drift coefficient (β₁ = 0), computed as $\left[1-\frac{\text{Coef}.0\text{ Mass}-\text{Coef}\text{. 90 Mass}}{\text{Coef}\text{. 0 Mass}}\right]\times 100\%$ . Similar interpretation applies to the tables below.

Table 2. Stage III tumor with a Therapeutic Coefficient of 10.

Stage 3 Cancer Tumor Evolution with Therapeutics Coef. = 10
Drift Coef.	Core Mass	Central Mass	Peripheral Mass	Total Mass
0	0.0986	11.7380	0.2832	12.1200
30	0.0759	08.3177	0.2026	08.5962
60	0.0410	04.9880	0.1222	05.1513
90	0.0137	01.7065	0.0420	01.7622
Ratio*	13.9%	14.5%	14.8%	14.5%

Table 3. Stage III tumor with a Therapeutic Coefficient of 20.

Stage 3 Cancer Tumor Evolution with Therapeutics Coef. = 20
Drift Coef.	Core Mass	Central Mass	Peripheral Mass	Total Mass
0	0.0855	10.5570	0.2562	10.8980
30	0.0644	07.1747	0.1756	07.4148
60	0.0313	03.8722	0.0953	03.9990
90	0.0051	00.6105	0.0151	00.6308
Ratio*	6.0%	5.8%	5.9%	5.8%

Table 4. Stage III tumor with a Therapeutic Coefficient of 30.

Stage 3 Cancer Tumor Evolution with Therapeutics Coef. = 30
Drift Coef.	Core Mass	Central Mass	Peripheral Mass	Total Mass
0	0.0725	9.3726	0.2291	9.6742
30	0.0533	6.0359	0.1486	6.2378
60	0.0221	2.7629	0.0683	2.8533
90	0.0003	0.0000	0.0000	0.0003
Ratio*	0.4%	0.0%	0.0%	0.0%

Table 5. Stage III tumor with a Therapeutic Coefficient of 40.

Stage 3 Cancer Tumor Evolution with Therapeutics Coef. = 40
Drift Coef.	Core Mass	Central Mass	Peripheral Mass	Total Mass
0	0.0597	8.1840	0.2019	8.4456
30	0.0425	4.9024	0.1215	5.0664
60	0.0133	1.6611	0.0413	1.7158
90	0.0002	0.0000	0.0000	0.0002
Ratio*	0.3%	0.0%	0.0%	0.0%

From the above table, one can glean an important piece of information-the significant relative rate of reduction of tumor mass occurring at the core region of the tumor domain. This fact is maintained by all other similar cases. Accordingly, if one succeeds in convecting most or all cancer cells to that region of interest, a significant reduction of mass must result. Bear in mind that the core region can be replaced by any other local region and we expect identical results.

In what follows we present the numerical results we obtained through COMSOL simulation when we fix the Drift Coefficient to a moderate level of 30 and vary the Therapeutics Coefficients.

Table 6. Stage III tumor with a Drift Coefficient of 30.

Stage 3 Cancer Tumor Evolution with Drift Coef. = 30
Therapeutics Coef.	Core Mass	Central Mass	Peripheral Mass	Total Mass
0	0.0875	9.4640	0.2296	9.7811
10	0.0759	8.3177	0.2026	8.5962
20	0.0644	7.1747	0.1756	7.4148
30	0.0533	6.0359	0.1486	6.2378
40	0.0425	4.9024	0.1215	5.0664
Rate*	51.4%	48.2%	47.1%	48.2%

A remarkably consistent trend immediately emerges. For each fixed value of the Therapeutic Agent and in all regions of the tumor, the cancer body mass decreases as the Drift Coefficient value increases. Even in the absence of therapy, the Drift of cancer cells to the core region alone leads to a significant decrease of the cancer body mass in the entire tumor domain. This unintuitive, curious phenomenon can be explained as follows:

By inspecting the cancer cells PDE Equation (4) we see that the net rate of growth of the cancer cells is the result of two competing terms: the logistics term and the intrinsic rate of death of cancer cells. The logistics term has the effect of reducing the rate of growth in overcrowded regions where cancer cells compete for finite resources. On the other hands, the rate of death remains constant since it is intrinsic. The net rate of growth at the core of the tumor where the cancer populations is overcrowded can, therefore, be negative. Since most of the cancer populations is drifting towards the core over time, where the growth rate can be negative, it is no surprise then that the drift alone can lead to reduction of the overall cancer population in the entire tumor domain and in the absence of therapy. This qualitative mathematical result has the following corresponding physical justification: If it were possible to drift most or all cancer cells to the core region of the tumor domain, or to any local region farther away from the vascular region, cancer cells would be kept away from vital resources of oxygen and various nutrients, and thus deprived from essential ingredients of survival. Cells must compete for finite resources—a phenomenon captured by the reaction part of the logistics term of the cancer density equation (4). As a consequence, the net growth rate of cancer cells declines and can even become negative. Indeed, Table 1 provides a numerical illustration of this phenomenon.

This finding points to a potentially original therapeutic avenue, which might be described as therapy by induced overcrowding (or equivalently, induced hypoxia therapy). The mechanism can be understood directly from the structure of the logistic growth law. When the local cancer cell density surpasses the carrying capacity (C > K), the logistic term turns negative, driving net decline rather than growth. In our model, chemotactic drift concentrates cells into the tumor core, where this overcrowding condition is inevitably satisfied. The result is a mathematically guaranteed suppression of growth in the central region, which in turn reduces the total tumor mass. Remarkably, this decline emerges without recourse to traditional cytotoxic interventions such as chemotherapy or radiation. Instead, it exploits the tumor’s own self-limiting growth dynamics under conditions of induced density stress. Such a strategy suggests a conceptual departure from conventional approaches, highlighting the possibility of non-toxic modalities that suppress tumor burden by engineering microenvironmental crowding or hypoxic stress, thereby harnessing intrinsic biological constraints rather than external destructive agents.

Our analysis, therefore, reveals a novel therapeutic paradigm: therapy by induced overcrowding (or induced hypoxia therapy).

As a side note, we mention an important fact related to the above discussion: Cancer cells, when forced towards the core region or strategically migrate there to evade the preying of the immune system, face intense competition due to overcrowding. This hypoxic environment favors the survival of the fittest cells, leading to the emergence of aggressive variants that are hypothesized to be primarily responsible for metastatic progression [31]. Research indicates that hypoxic conditions within the tumor core create a selection pressure favoring resilient, fitness-enhancing phenotypes and stemlike tumor cells [32]. These cells can migrate to the tumor boundary and contribute to metastatic dissemination [6] [33].

The tumor core’s inhospitable environment can indeed select for the most resilient cancer cells, which may then play a pivotal role in the metastatic progression of the tumor. This concept aligns with current understandings in oncology regarding tumor microenvironments and metastasis [34]. However, currently, our model will not incorporate the metastatic stage 4 of the tumor.

Summary: All simulations point to a significant decrease of cancer cells body mass and possible cancer eradication from the tumor microenvironment, leaving a controlled and treated tumor region with minimal to no remaining cancer cells. The above tables and their corresponding graphs illustrate this result.

5. Conclusions and Future Directions

Our overall strategy is to effectively manipulate and eradicate cancer cells within the tumor microenvironment by using the combined effect of chemo-attractants and therapeutic agents. The high concentration of chemo-attractants in the core region creates a gradient that attracts cancer cells towards the core. This results in the migration towards the core region of the tumor, keeping cancer cells at a minimum distance of 200 microns from the vascular region, thus achieving two goals at once: reducing the likelihood of metastasis and preventing the cells from accessing nutrients and oxygen from the blood vessels. Therapeutic agents at the tumor core target and kill cancer cells as cancer cells accumulate. The combined effect of the accumulation of cancer cells in a specific local region of the tumor due to chemo-attractants enhances the efficacy of the therapeutic agents.

There are several potential complicating factors that have not been accounted for in the current version of the model: Heterogeneity: Variations in cell responsiveness to chemo-attractants and chemorepellents could affect the uniformity of outcomes; Therapeutic Resistance: Some cancer cells might develop resistance to the therapeutic agent, necessitating adjustments in treatment strategy. Model Parameters: The effectiveness of the model depends on accurate parameter values for diffusion coefficients, sensitivity coefficients, decay rates, and production rates. For the time being, we have chosen the values of these parameters, through their dimensionless counterparts, rather arbitrarily. In the future, one can try adapting the model to actual clinical cases once the true values of the relevant physical parameters are known.

Notwithstanding the above limitations, this strategy of combining chemotaxis and targeted therapy, demonstrably offers an effective treatment paradigm of cancer tumor. Whether this paradigm finds clinical merit in the field of oncology remains to be seen.

Our model is fundamentally mathematical. We lay no claims of biological soundness. Our knowledge of the cancer tumor is essentially common knowledge. The tumor we are modeling is a solid tumor that may or may not correspond to any realistic one. Nevertheless, we believe that our model could be the beginning of a long-term endeavor to better understand this illusive disease. In deepening our understanding of its evolutions, we strive to formulate models that incorporate efficient treatments that are ontologically valid and mathematically optimal. Some possible improvements on the current model are to incorporate chemo-repellents to strengthen our ability to manipulate cancer cells into migrating to certain desired local regions in the tumor domain. Negative chemotaxis can be exploited by designing potential anti-metastatic agents [35]. Cellular signaling pathways that mediate chemorepulsion in cancer cells give insights into how modulation of these pathways can influence cancer cell migration and dissemination [36], and hence offer a potential therapeutic opportunity. Recent developments in repelling cancer cells and impeding tumor spread through molecular inhibition [12] are yet another opportunity for a novel therapeutic approach. As highlighted above, we believe that the combined effect of chemo-attractive and chemo-repulsive signals in the tumor microenvironment is a potentially promising modality that aims to tip the balance against cancer progression. Here we have attempted to give a mathematical visualization of the effect of chemo-attraction. As a possible enhancement, one can attempt to illustrate the dynamics of chemo-repulsion and the combined effect on the efficacy of therapy.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Griffith, J.W., Sokol, C.L. and Luster, A.D. (2014) Chemotaxis and Migration of Immune Cells in Health and Disease. Nature Reviews Immunology, 14, 514-528.
[2]	Heng, A.H.S., Han, C.W., Abbott, C., McColl, S.R. and Comerford, I. (2022) Chemokine-Driven Migration of Pro-Inflammatory CD4+ T Cells in CNS Autoimmune Disease. Frontiers in Immunology, 13, Article ID: 817473.[CrossRef] [PubMed]
[3]	Roussos, E.T., Condeelis, J.S. and Patsialou, A. (2013) Cell Migration and Chemotaxis in Tissue Repair. Current Opinion in Cell Biology, 25, 557-565.
[4]	Shellard, A. and Mayor, R. (2020) All Roads Lead to Directional Cell Migration. Trends in Cell Biology, 30, 852-868.[CrossRef] [PubMed]
[5]	Balkwill, F. (2004) Cancer and the Chemokine Network. Nature Reviews Cancer, 4, 540-550.[CrossRef] [PubMed]
[6]	Godet, I., Oza, H.H., Shi, Y., Joe, N.S., Weinstein, A.G., Johnson, J., et al. (2024) Hypoxia Induces ROS-Resistant Memory upon Reoxygenation in Vivo Promoting Metastasis in Part via MUC1-C. Nature Communications, 15, Article No. 8416.[CrossRef] [PubMed]
[7]	Spencer, J.V., Lockridge, K.M. and Barry, P.A. (2014) cmvIL-10 Stimulates the Invasive Potential of MDA-MB-231 Breast Cancer Cells. PLOS ONE, 9, e88789.
[8]	Van Haastert, P.J.M., et al. (2023) Competition between Chemoattractants Causes Unexpected Patterns of Chemotaxis and Chemorepulsion. Current Biology, 33, R374-R376.
[9]	Rempel, S.A., Dudas, S., Ge, S. and Gutierrez, J.A. (2000) Identification and Localization of the Cytokine SDF1 and Its Receptor, CXC Chemokine Receptor 4, to Regions of Necrosis and Angiogenesis in Human Glioblastoma. Clinical Cancer Research, 6, 102-111.
[10]	Smith, J.A. and Doe, R.C. (2020) Chemorepulsion: Negative Chemotaxis as a Mechanism to Inhibit Cancer Metastasis. Trends in Cancer Research, 6, 123-134.
[11]	Chan, E., Li, W. and Zhou, M. (2021) The Role of Chemorepellent Signals in Cancer Cell Migration and Metastasis. Nature Reviews Cancer, 21, 701-718.
[12]	Kumar, A., Patel, R. and Ng, S. (2023) Targeting Negative Chemotaxis: Therapeutic Opportunities in Inhibiting Cancer Cell Dissemination. Journal of Experimental Oncology, 15, 204-215.
[13]	Garcia, M.L. and Robinson, A. (2022) Balancing Chemoattraction and Chemorepulsion in Cellular Migration: Implications for Cancer Therapy. Cellular Mechanisms, 10, 301-317.
[14]	Lee, J. and Kim, S. (2020) Cell Polarization and Cytoskeletal Reorganization: Mechanisms of Directional Migration. Cellular Signalling, 68, 109-118.
[15]	Johnson, M.T. and Rivera, E. (2021) Asymmetric Organization of Cellular Components in Directed Migration: The Role of Cell Polarization. Journal of Cell Biology, 220, e202101045.
[16]	Rappel, W. and Edelstein-Keshet, L. (2017) Mechanisms of Cell Polarization. Current Opinion in Systems Biology, 3, 43-53.[CrossRef] [PubMed]
[17]	Dowdell, A., Paschke, P.I., Thomason, P.A., Tweedy, L. and Insall, R.H. (2023) Competition between Chemoattractants Causes Unexpected Complexity and Can Explain Negative Chemotaxis. Current Biology, 33, 1704-1715.e3.[CrossRef] [PubMed]
[18]	Do, H.T.T., Lee, C.H. and Cho, J. (2020) Chemokines and Their Receptors: Multifaceted Roles in Cancer Progression and Potential Value as Cancer Prognostic Markers. Cancers, 12, Article No. 287.[CrossRef] [PubMed]
[19]	Liu, M., Yang, J., Xu, B. and Zhang, X. (2021) Tumor Metastasis: Mechanistic Insights and Therapeutic Interventions. MedComm, 2, 587-617.[CrossRef] [PubMed]
[20]	Tang, L., Peng, S., Zhuang, X., He, Y., Song, Y., Nie, H., et al. (2025) Tumor Metastasis: Mechanistic Insights and Therapeutic Intervention. MedComm-Oncology, 4, e70012.[CrossRef]
[21]	Shi, X., Wang, X., Yao, W., Shi, D., Shao, X., Lu, Z., et al. (2024) Mechanism Insights and Therapeutic Intervention of Tumor Metastasis: Latest Developments and Perspectives. Signal Transduction and Targeted Therapy, 9, Article No. 192.[CrossRef] [PubMed]
[22]	Wang, F. (2009) The Signaling Mechanisms Underlying Cell Polarity and Chemotaxis. Cold Spring Harbor Perspectives in Biology, 1, a002980.[CrossRef] [PubMed]
[23]	Babar, Q., Saeed, A., Tabish, T.A., Sarwar, M. and Thorat, N.D. (2023) Targeting the Tumor Microenvironment: Potential Strategy for Cancer Therapeutics. Biochimica et Biophysica Acta (BBA)-Molecular Basis of Disease, 1869, Article ID: 166746.[CrossRef] [PubMed]
[24]	Liu, Y., Guo, J. and Huang, L. (2020) Modulation of Tumor Microenvironment for Immunotherapy: Focus on Nanomaterial-Based Strategies. Theranostics, 10, 3099-3117.[CrossRef] [PubMed]
[25]	Glaviano, A., Lau, H.S., Carter, L.M., Lee, E.H.C., Lam, H.Y., Okina, E., et al. (2025) Harnessing the Tumor Microenvironment: Targeted Cancer Therapies through Modulation of Epithelial-Mesenchymal Transition. Journal of Hematology & Oncology, 18, Article No. 6.[CrossRef] [PubMed]
[26]	COMSOL Inc. (2021) COMSOL Multiphysics® v. 5.6 User Guide. https://www.comsol.com
[27]	Farahani, B.V., Berardo, J.M., Drgas, R., César de Sá, J.M.A., Ferreira, A.J.M. and Belinha, J. (2015) The Axisymmetric Analysis of Circular Plates Using the Radial Point Interpolation Method. International Journal for Computational Methods in Engineering Science and Mechanics, 16, 336-353.[CrossRef]
[28]	Mohseni, K. and Colonius, T. (2000) Numerical Treatment of Polar Coordinate Singularities. Journal of Computational Physics, 157, 787-795.[CrossRef]
[29]	Zhang, B., Wang, H., Shen, S. and Wang, J. (2015) Intratumoral Delivery of Therapeutic Agents for Cancer Therapy. Journal of Controlled Release, 220, 585-594.
[30]	Goldberg, M.S. (2017) Intratumoral Drug Delivery: Future Directions and Challenges. Cancer Letters, 419, 77-84.
[31]	Schito, L. and Rey-Keim, S. (2023) Hypoxia Signaling and Metastatic Progression. Seminars in Cancer Biology, 97, 42-49.[CrossRef] [PubMed]
[32]	Emami Nejad, A., Najafgholian, S., Rostami, A., Sistani, A., Shojaeifar, S., Esparvarinha, M., et al. (2021) The Role of Hypoxia in the Tumor Microenvironment and Development of Cancer Stem Cell: A Novel Approach to Developing Treatment. Cancer Cell International, 21, Article No. 62.[CrossRef] [PubMed]
[33]	Belisario, D.C., Kopecka, J., Pasino, M., Akman, M., De Smaele, E., et al. (2020) Hypoxia Dictates Metabolic Rewiring of Tumors: Implications for Chemoresistance. Cells, 9, Article No. 2598.[CrossRef] [PubMed]
[34]	Chen, Z., Han, F., Du, Y., Shi, H. and Zhou, W. (2023) Hypoxic Microenvironment in Cancer: Molecular Mechanisms and Therapeutic Interventions. Signal Transduction and Targeted Therapy, 8, Article No. 70.[CrossRef] [PubMed]
[35]	Lai, W.Y. (2023) Investigations into Intracellular Mechanisms of Chemokine Signaling in Cancer. PhD Thesis, University of East Anglia.
[36]	Moore, M.A.S. (2001) The Role of Chemoattraction in Cancer Metastases. BioEssays, 23, 674-676.[CrossRef] [PubMed]