Optimizing cancer therapy for individuals based on tumor-immune-drug system interaction

1.   Introduction

Cancer is a severe public health problem. The International Agency for Research on Cancer estimated approximately 19.3 million new cancer cases and nearly 10.0 million cancer deaths occurred in 2020 ^[1]. Chemotherapy can slow down the reproduction of tumor cells, but can also cause side effects such as myelosuppression, immunosuppression, mucositis, and alopecia ^[2,3]. An increasingly trending topic is the use of mathematical models to design dosing treatments ^[4]. Most modeling guidance projects have focused on correlating drugs and tumor growth to discover the appropriate drug dosage ^[5,6]. However, responses to drugs vary substantially in patients with the same variety of tumors, as well as the same stage of tumors ^[7,8]. Chemotherapy causes a certain degree of either damage or inhibitory effects to the immune system; therefore, the efficacy of chemotherapy varies among individuals with different immune states ^[9]. Factors such as the microenvironment and immune status can significantly influence disease progression and the response in patients with cancer ^[8,10]. Therefore, customized chemotherapy that considers tumor-immune-drug (TID, Figure 1) interactions can minimize side effects and improve the quality of life, thereby significantly reducing the social and medical burden ^[8,10,11].

Figure 1.  Illustration of the tumor-immune-drug network. Immune cells grow faster in the presence of tumor cells. Immune cells inhibit tumor growth. Chemotherapy drugs have inhibitory effects on tumor cells and immune cells. Both tumor cells and immune cells have an apoptotic process, and drugs have a decay process.

DownLoad: Full-Size Img PowerPoint

There are many mathematical models that can be applied to evaluate TID interaction. Ordinary differential equations (ODEs) are a classic example of deterministic models. ODE-based tumor models are analyzed and optimized through computer simulations. It can be applied to tumor chemotherapy and immunotherapy evaluation ^[12]. S. Sameen et al. ^[13] further refined and utilized an ODE model for the analysis of the drug resistance characteristics of colorectal cancer based on Kirsten rat sarcoma virus (KRAS) mutations. Researchers similar to Sharma and Samanta applied the ODE model to explore the theoretical basis for the interaction between cancer and obesity and the immune response ^[13,14,15]. Other researchers constructed unique models that effectively incorporated the essential characteristics of the ODE model combined with the inherent elements of different tumors. However, there is a lack of an integrated model that merges existing mathematical models to provide more accurate suggestions for drug treatments for patients with cancer.

Some researchers have developed mathematical models focusing on the interplay between tumor and immune cells, as well as the impact of therapeutic drugs, to optimize treatment strategies ^[15,16,17]. Others have explored the influence of stochastic elements and noise on tumor-immune system dynamics, revealing how random disturbances can alter the outcomes and effectiveness of treatments ^[18,19,20]. Another group of researchers has concentrated on the role of the immune system in fighting tumors, from the signaling mechanisms of cytokines to reviews and evaluations of existing mathematical models on immune factors ^[21,22,23]. In addition, one study takes a unique approach by delving into the etiological and epigenetic aspects of tumors, emphasizing the need for a new model of tumor growth ^[24]. In this study, we analyzed the dynamic scenarios of five distinct models in recent publications and developed a master mathematical model that incorporates the vibrant characteristics of all five models reported in previous publications. Then, we developed a matrix for determining an optimal treatment for an individual patient's condition and created a platform that would provide treatment suggestions for all stages of tumor and self-immunity. This mathematical platform would effectively give clinicians a quantitative recommendation for a gentle treatment that is optimized based on the patient's tumor and immune status.

2.   Materials and methods

2.1.   Model establishment

We compared and evaluated five distinctive models from recent publications ^[15] and discovered distinctions among them, mainly in their cell growth function and the interactions between tumor cells and immune cells. In each model, $T\left(t\right)$ represents the tumor cell biomass, $I\left(t\right)$ represents the immune cell biomass, and $D\left(t\right)$ represents the amount (or concentration) of chemotherapeutic drugs in the bloodstream at time $t$.

2.2.   Description of five distinctive models reported in the literature

2.2.1.   Model 1 ^{[15,16,18,19,20,21]}

The equation and curve for Model 1 are shown in Figure 2a. The parameters used in the model are presented in Table 1. The assumptions of the model are as follows:

Figure 2.  Dynamic scenarios from five models. Equations and curves of Model 1 (a), Model 2 (b), Model 3 (c), Model 4 (d), and Model 5 (e) under different conditions.

DownLoad: Full-Size Img PowerPoint

Table 1.  Parameters used in five distinctive models. Columns include the meaning, the value (if one parameter is shared by multiple models and has different values, we utilized the average value), the citation, and the unit of each parameter.

Constant	Description	Value	Units	Applied in
${r}_{1}$	Innate growth rate of tumor cells	5 [15,16,17,18,19,20,21,22,23,24]	hr-1	Models 1–5
${p}_{1}$	Mutual carrying capacity for tumor cells	0.01 [15,16,18,19,20,21,22,23,24]	ml-1	Models 1–4
${v}_{2}$	Constant influx rate of immune cells	1 [15,16,18,19,20,21,22,23,24]	ml*hr-1	Models 1–3, 5
${\mu }_{1}$	Regression term from immune to tumor	10 [15,16,17,18,19,20,21,24]	ml-1*hr-1	Models 1, 4, 5
${\mu }_{1}$	Regression term from immune to tumor	20 [22,23]	hr-1	Models 2, 3
${d}_{22}$	Loss of the immune cells due to tumor	0.1 [15,16,18,19,20,21,23]	ml-1*hr-1	Models 1, 3
$\text{­}{\mu }_{2}$	Promotion term from tumor to immune	0.5 [15,16,18,19,20,21,22,24]	hr-1	Models 1, 2, 5
$n$	Steepness coefficient of the immune cell recruitment curve by tumor cells	10 [15,16,18,19,20,21,22,23]	-	Models 1–3
${d}_{2}$	Per capita decay rate of immune cells without tumor cells	1 [15,16,18,19,20,21,22,23,24]	hr-1	Models 1–3, 5
${d}_{3}$	Decay rate of the drug.	0.1 [15,16,17,18,19,20,21,22,23,24]	hr-1	Models 1–5
$v\left(t\right)$	Dose of given drug.	0 [15,16,17,18,19,20,21,22,23,24]	ng*hr-1	Models 1–5
${d}_{11}$	Response to the drug for tumor cells	0.2 [15,16,17,18,19,20,21,22,23,24]	ng-1*hr-1	Models 1–5
${d}_{21}$	Response to the drug for immune cells	0.05 [15,16,17,18,19,20,21,22,23,24]	ng-1*hr-1	Models 1–5
${\rho }_{2}$	Antigenicity of the tumor	0.1 [22]	hr-1	Model 2
$g$	Constant	10 [22,23]	ml	Models 2, 3
${r}_{2}$	Immune cells growth rate	6 [17]	hr-1	Model 4
${p}_{2}$	Carrying capacity of B cells	0.2 [17]	ml-1	Model 4
$\beta$	Constant	0.00264 [24]	ml-1	Model 5

---

 | Show Table

DownLoad: CSV

1) In the absence of foreign drugs and immune cells, tumor cells grow logistically. The growth term of the tumor is considered to be tumor cell apoptosis, as represented by ${r}_{1}T(1-{p}_{1}T)$;

2) The repression of tumor cells by immune cells is proportional to the multiplication of immune cells and tumor cells, as represented by ${\mu }_{1}IT$;

3) Chemotherapy drugs have damaging effects on tumor cells and immune cells, as represented by ${d}_{11}TD$ and ${d}_{21}ID$, respectively;

4) The increase rate of immune cells is dictated by the influx rate of immune cells in the bloodstream, as represented by ${v}_{2}$;

5) As tumor cells stimulate the immune response, immune cells produce $\frac{{\mu }_{2}{T}^{n}I}{h+{T}^{n}}$, which leads to the promotion of immune cell activity;

6) Tumor cells can reduce immune cell cytotoxicity, as represented by ${d}_{22}TI$; and

7) The dose-response curve of a chemotherapeutic drug is dictated by the dose of the drug and its natural decay, as represented by $v\left(t\right)-{d}_{3}D$.

2.2.2.   Model 2 ^[22]

The equation and curve for Model 2 are shown in Figure 2b. The assumptions of model 2 are as follows:

1) Non-self-antigenicity of the tumor causes an increase of immune cells, as represented by ${\rho }_{2}T$;

2) Immune cells repress tumor growth, as represented by ${\mu }_{1}\frac{T{I}^{{n}_{1}}}{g+{I}^{{n}_{1}}}$; and

3) All other terms used for this model were described for Model 1.

2.2.3.   Model 3 ^[23]

The equation and curve used for Model 3 are shown in Figure 2c. The assumptions of the model are as follows:

1) Immune cells inhibit tumor growth, as represented by ${\mu }_{1}\frac{I{T}^{{n}_{1}}}{g+{T}^{{n}_{1}}}$. Notice this repression term is different from that used for the other models; and

2) All other terms used for this model were described for Model 1.

2.2.4.   Model 4 ^[17]

The equation and curve used for Model 4 are shown in Figure 2d. The assumptions of the model are as follows:

1) Immune cell growth restriction can be simulated by a logistic model, as represented by ${r}_{2}I(1-{p}_{2}I)$; and

2) All other terms used for this model were described for Model 1.

The TID dynamic was evaluated when the initial condition was set to the parameters ${(T}_{0} = 1, {I}_{0} = 1, {D}_{0} = 0$). Tumors grew fast during the plateau period. Next, the TID dynamic was evaluated when the initial condition was set to the parameters $({T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 0)$. The initial density of the tumor cells was much higher, as tumor growth decreased to the thickness observed during the plateau period. Then, the TID dynamic was evaluated when the initial condition was set to the parameters $({T}_{0} = 10, {I}_{0} = 1, {D}_{0} = 10)$. Initially, the drug inhibited the tumor's density; however, the tumor grew gradually to the thickness that was observed during the plateau period. Finally, the TID dynamic was evaluated when the initial condition was set to the parameters ${(T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$. The initial density of the tumor cells was much higher. Initially, the drug inhibited the density of the tumor (approximate density of 1); however, the tumor grew gradually to the thickness that was observed during the plateau period.

Model 4 indicated that, without a drug, the tumor multiplied out of control and failed to stimulate an immune response. However, the initial density of tumor cells was much higher, as tumor growth decreased to the thickness it maintained during the plateau period. The drug initially inhibited tumor growth; however, the tumor subsequently grew to the density it maintained during plateau period and became out of control, regardless of the initial thickness of tumor cells, without stimulating the immune response.

2.2.5.   Model 5 ^[24]

The equation and curve for Model 5 are shown in Figure 2e. The assumptions of the model are as follows:

1) By assuming a fixed limited bearing capacity, the immune response pair reduces the tumor volume; and

2) The immune system's response to tumors is state-dependent (i.e., small tumors stimulate the proliferation of immune cells, whereas large tumors inhibit the activity of the immune system), as represented by ${\mu }_{2}T(1-\beta T)I$.

3.   Results

3.1.   Assessment of five distinctive models reviewed from the literature

3.1.1.   Model 1 ^{[15,16,17,18,19,20]}

When assessing the tumor cells, immune cells, and drug density dynamics for different initial conditions in Model 1, we observed the following.

When the initial condition was set to $({T}_{0} = 1, {I}_{0} = 1, {D}_{0} = 0)$, the TID dynamic involved an initial increase and then decrease of the tumor cell density, maintaining a viscosity similar to that of immune cells without drugs. When the initial condition was set to ${(T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 0)$, the TID dynamic involved a rapid increase of tumor cell density until it reached the plateau stage, while the immune cell density remained low without drugs. When the initial condition was set to ${(T}_{0} = 10, {I}_{0} = 1, {D}_{0} = 10)$ and drug was added, the TID dynamic involved a decrease of tumor cell density from a high position until it maintained a viscosity similar to that of the immune cells. When the initial condition was set to ${(T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$, the TID dynamic involved a gradual increase in tumor cell density until it reached a high plateau, when the drug lost its effect, while immune cell density remained low.

Model 1 revealed that only a small number of tumor cells stimulated the immune response to achieve a balance between tumor and immune cells (Figure 2a). When there was a large number of tumor cells, the immune response was suppressed, and the tumor cells grew out of control (Figure 2a). If the number of tumor cells increased to a certain extent, chemotherapy drugs inhibited the tumor cells, stimulated the immune response, and balanced the proportion of tumor and immune cells (Figure 2a). If the number of tumor cells was too large, even chemotherapy drugs could not inhibit tumor growth from proliferating out of control (Figure 2a).

3.1.2.   Model 2 ^[21]

When assessing the tumor cells, immune cells, and drug density dynamics for different initial conditions in Model 2, we observed the following.

When the initial condition was set to (${T}_{0} = 1, {I}_{0} = 1, {D}_{0} = 0$) and (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 0$), the TID dynamic indicated that, without drugs, tumor cell density fluctuated in correlation with immune cell density. When immune cell viscosity was high, the tumor cell density was low. When the immune cells were depleted to a low viscosity, the tumor cell density increased again. When the initial condition was set to (${T}_{0} = 10, {I}_{0} = 1, {D}_{0} = 10$) and (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10$), the TID dynamic indicated that with the addition of the drug, the immune response was initially suppressed to a certain extent. Later, the density of subsequent tumor cells fluctuated in correlation with the density of the immune cells. Finally, all four groups reached a balance between tumor and immune cells.

The addition of drugs inhibited the immune response to a certain extent, which indicated that regardless of the initial density, tumor cells stimulated the immune response and, in turn, immune cells inhibited tumor growth in a wave-like dynamic process. The final immune cell density was higher than the tumor cell density, forming a steady-state (Figure 2b). However, in the dynamics of the system, the tumor cell density fluctuated with immune cell density. Eventually, the tumor cell density became proportionately lower than the immune cell density (Figure 2b).

3.1.3.   Model 3 ^[22]

When assessing the tumor cells, immune cells, and drug density dynamics for different initial conditions in Model 3, we observed the following.

When the initial condition was set to (${T}_{0} = 1, {I}_{0} = 1, {D}_{0} = 0$) and (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 0$), the TID dynamic indicated that, without drugs, tumor cells proliferate out of control, and the density of immune cells is always low. When the initial condition was set to (${T}_{0} = 10, {I}_{0} = 1, {D}_{0} = 10$) and (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10$), the TID dynamic indicated that with the addition of the drug, tumor cells grew slightly during the initial stage and then gradually increased to a high plateau, where they proliferated out of control.

Model 3 indicated that, regardless of the initial tumor cell density, without drugs, tumor cells proliferated out of control while failing to stimulate the immune response. Initially, the drug inhibited tumor growth to a certain degree, but it gradually failed to be as effective as when the tumors proliferated out of control.

3.1.4.   Model 4 ^[24]

When assessing the tumor cells, immune cells, and drug density dynamics for different initial conditions in Model 4, we observed the following.

When the initial condition was set to (${T}_{0} = 1, {I}_{0} = 1, {D}_{0} = 0)$, the TID dynamic indicated that the tumor cells proliferated rapidly during the plateau period. When the initial condition was set to ${(T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 0)$, the TID dynamic indicated that the initial density of tumor cells was much higher, but decreased during the plateau period. When the initial condition was set to (${T}_{0} = 10, {I}_{0} = 1, {D}_{0} = 10)$, the TID dynamic indicated that the drug initially inhibited the tumor's density; however, the tumor cell density gradually grew in thickness during the plateau period. When the initial condition was set to (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$, the TID dynamic indicated that the initial density of tumor cells was much higher. Although the drug initially inhibited the density of the tumor (approximate density 1), the tumor cell density gradually grew in thickness during the plateau period.

Model 4 indicates that, without a drug, the tumor multiplied out of control and failed to stimulate the immune response. However, the initial density of tumor cells was much higher, and tumor growth decreased to the density observed during the plateau period. The drug initially inhibited tumor growth; however, subsequently, the tumor density increased to a level that was observed during the plateau period and proliferated out of control, regardless of the initial density of the tumor cells. Meanwhile, the immune response was not stimulated.

3.1.5.   Model 5 ^[23]

When assessing the tumor cells, immune cells, and drug density dynamics for different initial conditions in Model 5, we observed the following.

When the initial condition was set to (${T}_{0} = 1, {I}_{0} = 1, {D}_{0} = 0)$, the TID dynamic indicated that the tumors stimulated an immune response, and the immune cells inhibited tumor growth, finally reaching a balance ($T = 1, I = 10$). When the initial condition was set to (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 0)$, the TID dynamic indicated that an initial tumor cell density of 50 proliferated rapidly until reaching the plateau period, when it proliferated out of control. When the initial condition was set to (${T}_{0} = 10, {I}_{0} = 1, {D}_{0} = 10)$, the TID dynamic indicated that drugs inhibited the initial immune response, finally reaching a balance ($T = 1, I = 10$). When the initial condition was set to (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$, the TID dynamic indicated that drugs inhibited the initial immune response and the tumor at the initial density of 50, finally reaching a balance ($T = 1, I = 10$).

Model 5 shows that a small number of tumor cells stimulated a robust immune response and reached a balance between tumor and immune cells. When there was a large number of tumor cells, the immune response was suppressed and tumor cells proliferated out of control. If the number of tumor cells increased to a certain extent, chemotherapy drugs promoted a stronger immune response, finally achieving a balance between tumor and immune cells. If there was a large number of tumor cells, chemotherapy drugs could inhibit tumor growth and the immune response, eventually reaching a balance between tumor and immune cells.

3.2.   Establishment of the master model

Based on the above five model assumptions (Table 1), applicable scenarios, and dynamics, we concatenated all five distinct equations into a master equation to capture all possible tumor-immune-drug relationships.

The master equation is shown in Figure 3a. The assumption of the master equation combines the assumptions of all five models. For simplification, the TID dynamic in all five models are all the same due to the small molecule dynamic as described by pharmaco-dynamic and pharmaco-kinetic measurements. Therefore, the master equation incorporates a drug input profile and drug degradation rate to describe its dynamics. We present six conditional curves in Figure 3b:

Figure 3.  Dynamic scenarios of the master equation and numerical simulation of the TID model. a. Master equation; b. Capturing dynamics in five distinctive published TID models; c. Taking drug effect function; d. The effect from drug influx rate to TID dynamic. TID, tumor-immune-drug.

DownLoad: Full-Size Img PowerPoint

a) The TID dynamic when the initial condition was set to (${T}_{0} = 1, {I}_{0} = 1, {D}_{0} = 0)$. The TID dynamic is similar to that of Model 1 under the same initial condition.

b) The TID dynamic when the initial condition was set to (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10$). The system dynamic is similar to that of Model 1 under the same initial condition.

c) The TID dynamic when the initial condition was set to (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$. The TID dynamic is similar to that of Model 2 under the same initial condition.

d) The TID dynamic when the initial condition was set to (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$. The TID dynamic is similar to that of Model 3 under the same initial condition.

e) The TID dynamic when the initial condition was set to (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$. The TID dynamic is similar to that of Model 4 under the same initial condition.

f) The TID dynamic when the initial condition was set to (${T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$. The TID dynamic is similar to that of Model 5 under the same initial condition.

As Figure 3b shows that the master equation can capture the dynamics of five distinctive published TID systems. This means that the master equation is not heavily assumption-based. When creating a customized chemotherapy treatment plan without this master model, one would need to determine which TID model describes the patient's situation the best, and then design the treatment according to that particular model. The master model eliminates this step of choosing which model to use. Regardless of the initial model assumption, the master model will describe various TID dynamics, thereby allowing the master equation to develop an optimized drug treatment.

3.3.   Numerical simulation of the TID model

After establishing the master equation, we investigated how we could use the model to design the best drug treatment to eliminate the tumor in the shortest timeframe. It is assumed that the chemotherapeutic drug is given to the patient as a continuous Dirac delta Function (see Figure 3c).

The initial condition of TID was fixed at ${(T}_{0} = 50, {I}_{0} = 1, {D}_{0} = 10)$. while ${t}_{1}, {t}_{2},$ and D_const were variables. Interestingly, when the drug influx rate was too low (${D}_{const} = 1$), the tumor growth rate slowed significantly, but still increased within an oscillation (Figure 3d). However, when ${t}_{1} = 2$, ${D}_{const} = 10$, tumor cells were eliminated from the system within three drug cycles (Figure 3d).

To conduct a more precise search, we searched thousands of different combinations of ${t}_{1}$, ${t}_{2}$, and ${D}_{const}$. We set an initial range for the values of ${t}_{1}$, ${t}_{2}$, and ${D}_{const}$ (e.g., $0\le {t}_{1}\le 2$; $\text{10≤}{t}_{2}\le 20$; and $0\le {D}_{const}\le 20$). We looked for the following:

● The minimum timeframe for tumor density to decrease under 0.01.

● The area of ${t}_{1}*{D}_{const}*$ number of drug cycles to eliminate tumor density to 0.01 as a minimum.

● ${t}_{1}$ and ${D}_{const}$ should be as short and small as possible, whereas ${t}_{2}$ should be as long as possible.

3.4.   Optimization of drug treatment

To test which drug treatment yields the best effect, we investigated four different drug treatment profiles (delta direct function, normal distribution function, sin function, and trapezoid function) as a chemotherapeutic drug treatment (see Figure 4a). To make a fair comparison, all parameters of one single dosage of each medicine, which was the quantity of the drug, were kept the same. The amount of the drug was 16. ${t}_{1}$ represents the time interval when the drug was added. ${t}_{2}$ represents the combination of ${t}_{1}$ and the time interval between each drug treatment. In the delta direct function graph, the concentration of the drug was represented ${D}_{const}$ when the time is in ${t}_{1}$; for all other time intervals, the concentration of the drug was 0.

Figure 4.  Optimization of drug treatment. a. The effect plots of different treatments include direct function, normal distribution function, sin function, and trapezoid function. Plots are generated through the matplotlib Python package. b. Effects on tumor density by different treatment profiles during the early stages of cancer, given T₀ = 1, 10, 50. T₀ is small. c. Effects on tumor density by different treatment profiles at the advanced stages of cancer, given T₀ = 100,300, and 500. T₀ is large.

DownLoad: Full-Size Img PowerPoint

In the typical distribution function graph, since the range of normal distribution is infinite, the range was considered to be three standard deviations from the median of normal distribution. In other words, ${t}_{1}$ had the length of six standard deviations. Only some parts of normal distribution within this range were maintained. Since the total area of normal distribution was 1,

We multiplied the value from a normal distribution with ${D}_{const\_n}$ to obtain the regular distribution function graph.

In the sin function graph, we transformed the sin function from 0 to pi, which was the positive part of the sin function, into the range of ${t}_{1}$, proportionally. We derived ${D}_{const\_s}$ from integration, which allowed the area of the graph to equal the total drug quantity, 16. We multiplied the value from the sin function with ${D}_{const\_s}$ to obtain the sin function graph.

In a trapezoidal graph, the corresponding lengths of two hypotenuses on the x-axis are two variables, ${t}_{1\_1}$ and ${t}_{1\_3}$. In this graph, they are 0.2 and 0.5, respectively. The bottom edge has a length of ${t}_{1}$. We derived the length of the top edge, ${t}_{1\_2}$, from the difference of ${t}_{1}$ and the sum of ${t}_{1\_1}$ and ${t}_{1\_3}$. Based on the trapezoid area formula,

we derived ${D}_{const\_t}$ and obtained the trapezoidal graph.

To examine the effect on cancer treatment imposed by different drug treatment profiles, we arbitrarily chose two sets of initial values–one representing an early stage of cancer (${T}_{0}$ is small tumor) and the other representing an advanced stage of cancer (${T}_{0}$ is large tumor). We drew the tumor density after using different drug treatment functions. In each graph, the initial values were the same. Different colors represented other drug treatment functions.

Figure 4b simulates the early stage of cancer. ${T}_{0}$ is the initial value of a tumor cell, and ${I}_{0}$ is the initial value of an immune cell. It can be seen that during the early stage of cancer, the trapezoidal drug treatment function (black curve) has the best inhibitory effect on tumor cell density. The curve corresponding to the trapezoidal drug treatment function has the lowest tumor cell density at almost all times.

Figure 4c simulates the advanced stage of cancer, which is stimulated by a very high ${T}_{0}$ value. It can be concluded that during the advanced stage of cancer, the effect of each drug treatment function is approximately the same. Each curve has almost the same tumor cell density. This conclusion aligns with the current medical observation that when a tumor is found at a later stage, it is too late for any drug treatment.

To measure the different drug treatment regimens more precisely, we established a matrix to quantitatively measure how well the drug can decrease the tumor density (treatment effect). The score considers the following:

● The duration of one treatment (the shorter the treatment time, the higher the score) accounts for 40% of the total score.

● The time between each treatment (the longer the time between each treatment, the higher the score) accounts for 40% of the total score.

● The number of drugs used to control the tumor under a threshold of $T < 0.01$ (the fewer drugs, the better) accounts for 10% of the total score.

● The time it takes to control the tumor under a threshold of $T < 0.01$ (the shorter the time, the higher the score) accounts for 5% of the total score.

● The total number of treatments patients received to control the tumor under a threshold of $T < 0.01$ (the fewer treatments, the better) accounts for 5% of the total score.

Since each drug treatment function performs approximately the same during the advanced phase of cancer, and the trapezoidal drug treatment function performs best for early-stage tumors, the trapezoidal drug treatment function was chosen to optimize future treatment designs.

To further analyze the contribution of each parameter to the score regimen, we generated 100,000 random combinations of parameters. We graphed them into polygons by using the seismic colormap and connecting points corresponding to parameters in each variety.

We used the initial condition of ${T}_{0} = 50, {I}_{0} = 1$ and generated 100,000 random combinations of ${t}_{1}, {t}_{1\_1, }{t}_{1\_2}, {t}_{1\_3}, {t}_{2, }{D}_{const\_t}$:

-${t}_{1}$ varies between 0 and 0.5,

-${t}_{1\_1}$ varies between 0 and ${t}_{1}$,

-${t}_{1\_2}$ varies between 0 and ${(t}_{1}-{t}_{{1}_{1}})$,

-${t}_{1\_3} = {t}_{1}-{t}_{1\_1}-{t}_{1\_2}$,

-${t}_{2}$ varies between ${(t}_{1}+1)$ and ${(t}_{1}+2)$, and

-${D}_{const\_t}$ is a random number between 0 and 50.

Black lines represent a score of 0 and 1. White bars represent a score of 0.5. The higher the score of the combination, the redder a line is. The lower the score, the bluer a line is. To leverage the seismic colormap to demonstrate the score distribution, each line's score was squared before graphing the line. If the square of the score of a line was higher than 0.5, then the line was graphed in this polygon. The color of the line was adjusted based on the seismic colormap and square of the score. A high score reflects a score whose square was greater than and equal to 0.5, whereas a low score reflects a square that was less than 0.5. Each corner represents the following parameters, respectively: ${t}_{1\_1}$, ${t}_{1\_2}$, ${t}_{1\_3}$, ${t}_{2}$, and ${D}_{const\_t}$. The center represents the minimum values of parameters, and the corners represent the maximum values of parameters. The value increases linearly as points move from the center to the vertex. The max and min values of ${t}_{1\_1}$, ${t}_{1\_2}$, and ${t}_{1\_3}$ are 0 and 0.5, respectively. The max and min values of ${t}_{2}$ are 1 and 2.5, respectively. The max and min values of ${D}_{const\_t}$ are 0 and 50, respectively.

3.5.   Application of the master model

We plotted two polygons. Parameter combinations that successfully eliminated tumors are graphed separately into the following polygons (Figure 5a). There were 3348 combinations graphed into a blue polygon and 31,387 combinations graphed into a red polygon. The red lines in the right graph are based on the seismic colormap.

Figure 5.  Application of the master model. a. Polygons with all successful treatments were separated by high score and low score. The polygon in blue represents parameter combinations with low scores, whereas the polygon in red represents parameter combinations with high scores. Polygons demonstrate successful treatments and their parameters related to immune and tumor. b. The initial condition of low immune cells and early-stage tumors. c. The initial condition of low immune cells and medium-stage tumors. d. The initial condition of low immune cells and advanced-stage tumors. e. The initial condition of high immune cells and early-stage tumors. f. The initial condition of high immune cells and medium-stage tumors. g. The initial condition of the high immune and advanced-stage tumors. High immunity is defined as the top 50% of ${T}_{0}, {\mu }_{2}, {r}_{2}$, and vice versa. The early-stage tumor is defined as the bottom 33% of ${T}_{0}, {\mu }_{2}, {r}_{2}$. The medium stage of tumor is defined as ${T}_{0}, {\mu }_{2}, {r}_{1}$, whose values are between the parameters of early-stage and advanced-stage tumor progression. The advanced stage of the tumor is defined as the top 33% of ${T}_{0}, {\mu }_{2}, {r}_{2}$. Polygons demonstrate successful high-score treatments and their parameters related to immune and tumor cells. h. An initial condition of low immune cells and early-stage tumor. i. An initial condition of low immune cells and medium-stage tumor. j. An initial condition of low immune cells and advanced-stage tumor. k. An initial condition of high immune cells and early-stage tumor. l. An initial condition of high immune cells and medium-stage tumor. m. An initial condition of high immune cells and advanced-stage tumor.

DownLoad: Full-Size Img PowerPoint

Based on two polygons, it can be concluded that combinations with high scores and low scores have similar ${t}_{1\_1}$, ${t}_{1\_2}$, ${t}_{1\_3}$ values. ${t}_{1\_1}$, ${t}_{1\_2}$, ${t}_{1\_3}$ values cannot be increased simultaneously. This is reasonable, because the sum of ${t}_{1\_1}$, ${t}_{1\_2}$, ${t}_{1\_3}$ is equal to t₁. The value of ${t}_{1\_2}$ is inversely proportional to the value of ${t}_{1\_1}$ and ${t}_{1\_3}$. Both ${t}_{1\_1}$ and ${t}_{1\_3}$ tend to be low. This indicates that the drug concentration should change in a short period, regardless of whether the change is either an increase or decrease. However, this is trivial because this pattern is shown in both graphs.

Combinations with high scores generally have a slightly higher ${t}_{2}$, and they have a tremendously higher ${D}_{const\_t}$ than combinations with lower scores. The minimum drug concentration was 12.01 in the high score graph, as shown by the green dot. The chart was empty for a ${D}_{const\_t}$ below 12.01, indicating that no combinations with a ${D}_{const\_t}$ lower than 12.02 can have a good score. When designing optimized treatments, the drug concentration must be greater than 12.01, regardless of other parameters.

The green dot had the most significant ${t}_{2}$ value in the blue polygon when the ${D}_{const\_t}$ was at its max. An orange line is plotted to show that this dot has the max ${D}_{const\_t}$ value. Other points in the red polygon have a more considerable ${t}_{2}$ value when the ${D}_{const\_t}$ is the max, compared to the green dot in the blue polygon. This reveals that the ${t}_{2}$ value of the green dot in the blue polygon was the min ${t}_{2}$ value for a high-score treatment. When ${t}_{2}$ was more extensive than this value, there was a greater possibility that the treatment would be optimized with a high score.

It is recognized that each patient has a different tumor and immune status. Therefore, parameters indicating tumor and immune statuses were included in the graphs (Figure 5b–g). To simplify the diagram, ${t}_{1\_1}$, ${t}_{1\_2}$, ${t}_{1\_3}$ were not graphed, whereas ${t}_{1}$ was graphed, as it is the sum of ${t}_{1\_1}$, ${t}_{1\_2}$, ${t}_{1\_3}$.

Blue lines represent parameter combinations with low scores, and red lines represent high scores. A high score is a score whose square is more significant than or equal to 0.5, whereas a low score is a score whose square is less significant than 0.5. Each corner represents the following parameters, respectively: ${t}_{1}$, ${t}_{2}$, and ${D}_{const\_t}$, ${T}_{0}$, ${r}_{1}$, ${\mu }_{2}$, ${I}_{0}$, ${r}_{2}$, ${\mu }_{1}$. The center represents the minimum values of parameters, whereas the corners represent the maximum values of parameters. The max and min value of ${t}_{1}$ are 0 and 0.5, respectively. The max and min values of ${t}_{2}$ are 1 and 2.5, respectively. The max and min values of ${D}_{const\_t}$ are 0 and 50, respectively. The max and min values of ${T}_{0}$ are 1 and 500, respectively. The max and min values of ${I}_{0}$ and ${\mu }_{1}$ are 0.1 and 5, respectively. The max and min values of ${r}_{1}$, ${\mu }_{2}$, ${r}_{2}$ are 0.1 and 10, respectively. Figure 5d, g have more blue lines than Figure 5b, e. This suggests that it would be more challenging to find effective parameter combinations for treating a patient when the patient has an advanced stage tumor. This is intuitive because more tumor cells exist in patients' bodies in a developed location of cancer. Figure 5e–g have more red lines than Figure 5b–d. This indicates that when the immunity is strong, there are more efficient parameter combinations than those when the immunity is weak. Therefore, it is essential to increase patients' self-immunity, because it would make treatment more effective.

Figure 5h–m show combinations with high scores. The two lines in yellow and cyan in Figure 5m each represent parameter combinations with high scores. However, the yellow line has a high ${t}_{2}$ and low ${D}_{const\_t}$, whereas the cyan line has a high ${D}_{const\_t}$ and low ${t}_{2}$. The drug concentration and time between each treatment are much longer for the yellow regimen than the cyan regimen, which is a mild, favorable treatment. This shows that instead of qualitative suggestions, this platform can effectively give quantitative recommendations for discovering a gentle treatment that has been optimized and individualized for a patient's personal tumor and immune systems.

4.   Discussion

After investigating five distinct models in the current publications in this project, we incorporated their essential features into a master mathematical model for dynamic programming of TID system interaction to optimize cancer therapy. We thoroughly studied the model assumptions and system dynamics for all five published models. First, we analyzed publications about the TID model that were published after 2010 and summarized the five most representative models that can cover all models in recent publications. Then, we built a master model that can cover the dynamic characters of all five summarized models. Now, we have one single set of equations that can describe all TID model dynamic publications published after 2010, which is essential for designing an optimized, customized drug treatment plan.

We investigated four different profiles of chemotherapeutic drug treatments (direct function, normal distribution function, sin function, and trapezoid function) to test which treatment yielded the best effect. The trapezoidal drug treatment function had the best inhibitory effect on tumor cell density during the early stage of cancer. We modeled the tumor density after using different drug treatment functions with the same initial values. The curve corresponding to the trapezoidal drug treatment function had the lowest tumor cell density at almost all times. However, in the advanced stage of cancer, the effect of each drug treatment function was approximately the same.

Knowing the dynamics of the tumor-immune system is essential and instructive for designing future drug treatments. We investigated how the master model can be applied to create an optimal drug treatment. In this project, we developed a scoring matrix for an optimal treatment that is effective and painless for patients. The trapezoidal drug treatment function was used in this step to establish an optimized treatment.

We began assessing a customized cancer therapy based on an individual case first. Given a patient with early-stage cancer with weak self-immunity, we first generated 10,000 random combinations of parameters. We plotted two polygon graphs using a seismic colormap mapping non-effective vs. effective treatment. The results showed that when designing optimized treatments, the concentration of drugs should change over a short period of time, the drug concentration must be above 12.01 regardless of other parameters and t₂ has a min value for a high-score treatment. When t₂ is more extensive than this value, the treatment will likely be optimized with a high score.

Then, we generated successful treatments for broader situations. Namely, we established approximate effective treatment plans at different stages of cancer and levels of self-immunity. The results showed that there were more efficient parameter combinations when the immunity was strong compared to when the immunity was weak. Therefore, it is imperative to increase patients' self-immunity, because it makes treatment more effective. This platform can provide quantitative recommendations for establishing a suitable, personalized treatment based on patients' tumors and immune systems instead of qualitative suggestions.

Other than designing a new drug, our method focused on improving treatment efficiency by improving the treatment regime using existing drugs. This methodology has been proven ^[25]. Meredith et al. ^[25] developed a metric to guide the design of a dosing protocol to conquer antibiotic resistance. The protocol intends to optimize treatment efficacy for any antibiotic-pathogen combination. The wet lab experimental data proved the validity of designing a more appropriate the treatment regime to improve treatment outcomes. We believe by combining web lab data with our model, anti-tumor drug treatment regimens can be better designed and customized.

Our research presents certain limitations. First, the spatial evolution cannot be accomplished solely through the application of ODEs. Second, the proposed computational model has not yet been substantiated with clinical validation.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Author contributions

Xin Chen downloaded and processed the data, and wrote the paper; Yangxiaolu Cao provided some advice for this research; Tengda Li edited the paper; Yangxiaolu Cao discussed with Xin Chen, and provided suggestions for the designation of this study.

Availability of data and materials

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

Conflicts of interest

The authors declare they have no competing conflicts of interest.