A Mathematical Model to Optimize the Neoadjuvant Chemotherapy Treatment Sequence for Triple-Negative Locally Advanced Breast Cancer

Abstract

Background: Triple-negative locally advanced breast cancer is an aggressive tumor type. Currently, the standard sequence treatment is applied, administering anthracyclines first and then a taxane plus platinum. Clinical studies for all possible treatment combinations are not practical or affordable, but mathematical modeling of the active mitotic cell population is possible. Our study aims to show the regions with the tumor’s most substantial cellular population variation by utilizing all possible values of the parameters $\left({\alpha }_s^i\right)$ that define the annihilatory drug capacity according to the proposed treatment. Method: A piecewise linear mathematical model was used to analyze the cell population growth by applying four treatments: standard sequences of 21 days (SS21) and 14 days (SS14), administering anthracyclines first, followed by a taxane plus platinum, and inverted sequences of 21 days (IS21) and 14 days (IS14), administering a taxane plus platinum first then anthracyclines. Results: The simulation showed a higher effect of IS14 over SS14 when the rate of drug resistance was larger in the cell population during DNA synthesis (G1 and S) compared to cells in mitosis (G2 and M). However, if the proportion of resistant cells in both populations was equivalent, then treatments did not differ. Conclusions: When resistance is considerable, IS14 is more efficient than SS14, reducing the tumor population to a minimum.

1. Introduction

During the last 50 years, different mathematical models have been developed to illustrate how a cancerous tumor’s growth begins, its dynamics, and treatments. Some models include the organism’s resistance to the different drugs that are usually administered as the first step to reduce the size of the tumor (neoadjuvant therapy) before the primary treatment, usually consisting of surgery for all types of cancer [1,2,3,4] but especially for breast cancer [5,6,7,8,9,10,11].

Because breast cancer is one of the most studied pathologies in the world, many efforts have been devoted to understanding the relationships between immune cells, tumor cells, and certain adjuvant treatments [6,12,13,14]. To cite an example, Jarett et al. [12] discussed the effects of trastuzumab on the overexpression of the HER2 gene, which is characteristic of breast cancer, using a mathematical model and integrating the experimental results [12,13,14]. Other examples include those where neoadjuvant treatment affects the immune system and its interaction with cancer cells [12,14,15].

When cancer is in its early stages, surgery is performed first and then radiotherapy. Therefore, implementing the reverse sequence (radiotherapy first and later surgery) has been widely ruled out for many cancers. For instance, Lopez-Alonso and Poleszczuk J. et al. [16,17] reversed the traditional treatment (radiotherapy first to induce antitumor immunity and then surgery) for different cancers and evaluated the overall survival (OS) and the disease-free survival (DFS), obtaining better results. We were interested in inquiring if a reverse or different sequence would produce better results compared with a traditional one as a research question.

Triple-negative breast cancer (no expression of estrogen receptor, progesterone receptor, and growth factor 2 receptor (Her2)) represents 20% of breast cancer cases in young women, with early recurrence and dissemination to the viscera and central nervous system. In patients with triple-negative, locally advanced breast cancer (TN-LABC; T3-T4, N1-N3), unlike other molecular subtypes, the complete pathological response (pCR) is a surrogate for more remarkable overall survival and DFS compared with that of patients with residual disease [18,19].

The order of chemotherapy administration before or after surgery may change the outcome. For example, the National Surgical Adjuvant Breast and Bowel Project (NSABP) Study B-18 demonstrated a 13% pCR rate when doxorubicin (DX) (60 mg/m²) plus cyclophosphamide (CPh) (600 mg/m²) was administered every 21 days for four cycles during the pre-surgery period. This chemotherapy slightly improved the survival and disease-free time compared with the same chemotherapy administered in the post-surgery period. However, the NSABP Study B-27 later showed that adding docetaxel every 21 days for four cycles after the standard neoadjuvant sequence with DX and CPh increased the pCR to 26% (p < 0.001) [20].

Adding carboplatin (CP) to the sequential neoadjuvant chemotherapy scheme increases the probability of the pCR in patients with TN-LABC by 13% [21]. However, other long-term results of randomized phase II studies, CALGB 40603 (Alliance) [22] and GeparSixto [23], showed that CP increases the pCR but does not improve the event-free interval. The low number of patients may explain these controversial findings.

The morphological phenotypes in TN-LABC might determine a low response to systemic therapy in TN-LABC; other poor-prognosis factors are inflammatory carcinoma in individuals older than 40 years and a high proliferation index [24]. Different susceptibilities to treatment depend on their genotype (tumor markers), e.g., basal-like (BL) tumors carrying pathogenic mutations in the BRCA1 gene (57%) or BRCA2 gene (23%) respond to treatment with platinum agents [25]. In a study with 290 patients with TN-LABC, 47% of BL-1 tumors responded to chemotherapy compared with 28% of BL-2 tumors [26]. Additionally, patients with tumors expressing p53 and Ki-67 protein markers treated with taxanes might respond better to neoadjuvant chemotherapy [27,28].

When the tumor’s interstitial fluid pressure (IFP) increases, the bioavailability of antineoplastic drugs decreases. Therefore, the drugs do not reach the tumor cells at a sufficient concentration. When paclitaxel (PX) is administered after anthracyclines, the IFP increases. However, PX is administered before anthracyclines. In that case, the IFP decreases [29,30], and the order of the treatment sequence in TN-LABC patients may increase reduction rates and favor increases in complete pathological responses.

Mathematical models help to explain the response probability to treatments and their combinations. The logarithmic death or log-kill is a cancer model with a constant exponential growth fraction per time step. However, the presence of effective cancer drugs also decreases the tumor size by a constant fraction. For example, if a drug eliminates 90% of a tumor cell population and a second drug destroys 90%, in the end, they may eliminate 99% of the tumor cells [31]. A modification to the log-kill is the Norton–Simon hypothesis, which relates the treatment effects to the growth rate and the tumor size [32]. For example, a small tumor with a high growth rate can be eliminated using a particular medication at a specific dose; however, the same drug has a weaker effect on treating tumors of a larger size with a low growth rate.

A piecewise linear mathematical model (PLMM) based on Roe-Dale R. et al. [13] was proposed to analyze four treatment schedules (standard sequence, 21 and 14 days, and inverted sequence, 21 and 14 days) [22]. The novelty of our PLMM is a wide spectrum of growth rates of the active mitotic cell population (AMCP) for treatment sequences and chemotherapy cell resistance.

2. Materials and Methods

Four chemotherapy sequences were described using standard doses of DX, CPh, CP, and PX (defined in Section 2.1); then, the evolution of the AMCP was described using linear ordinary differential equations (ODEs). Finally, we determined which treatment sequence had a higher probability of success for neoadjuvant chemotherapy patients with TN-LABC to reach a maximum tumor reduction that could direct a better pCR and more conservative breast surgeries.

2.1. Standard Doses

The standard dose scheme we used for our model is detailed below:

(a)

The DX dose was 60 mg/m² body surface area every 21 or 14 days for 4 cycles.

(b)

The CPh dose was 600 mg/m² body surface area every 21 or 14 days for 4 cycles. The treatment administered every 14 days is proposed by the NCCN guidelines in the United States of America; however, some countries continue administering treatment every 21 days as described in previous studies [20]. The treatment administration every 14 days is associated with unacceptable hematological toxicity. That must be balanced by administering granulocyte colony-stimulating factor (GCSF) [33].

(c)

The PX dose was 80 mg/m² of body surface area intravenously every week for 12 weeks.

(d)

The CP dose was obtained through the Calvert formula, $Dose=\left(IFG+25\right)\times AUC$, where IFG is the glomerular function index or creatinine clearance; this is the volume of fluid filtered per time unit from renal glomerular capillaries into Bowman’s capsule, usually measured in milliliters per minute. The IFG varies for each patient. For the CP every three weeks, which was our schedule, the AUC = 6 (6 units equals 6 mg. min/mL), where the AUC is the area under the curve of free plasma carboplatin concentration versus time. This is a method used to reduce toxicity based on renal clearance values calculated using the age and health condition of the patient.

2.2. Standard Sequence of 21-Day Cycle (SS21)

DX plus CPh was administered every 21 days for four cycles during the first phase. In the second phase, PX was issued weekly on days 1, 8, and 15. Moreover, CP was administered on day 1 of each of the four cycles, as shown in Figure 1a.

2.3. Standard Sequence of 14-Day Cycle (SS14)

The chemotherapy sequence was the same as the standard dose but with a DX and CPh administration every 14 days for four cycles, as shown in Figure 1b.

Neoadjuvant therapy with SS21 (Figure 1a) is considered a safe and efficient option for treating TN-LABC. However, because this therapy can induce drug resistance [34], treatment with SS14 (Figure 1b) prevents the repopulation of tumor cells but at the cost of higher toxicity. Therefore, GCSF was administered in patients with breast cancer receiving neoadjuvant chemotherapy per cycle of chemotherapy.

2.4. Inverted Sequence of 21-Day Cycle (IS21)

In the inverted sequence, phase 2 was used as the first treatment; that is, PX was administered on days 1, 8, and 15, and CP on day 1 of each cycle for four cycles. Next, DX and CPh were sequentially administered using standard doses every 21 days, as shown in Figure 2a.

2.5. Inverted Sequence of 14-Day Cycle (IS14)

This chemotherapy was identical in dosage to that described above, but DX and CPh were administered every 14 days for four cycles, as shown in Figure 2b.

Treatment was inverted beginning with phase 2, starting with PX treatment, then CP was added on day 1 of each cycle, followed by phase 1, DX plus CPh. This was used to prevent cross-resistance (when acquired resistance induced by a drug treatment results in resistance to other drugs). This significantly increases the possibility of pCR and more conservative breast surgery in TN-LABC patients [21].

2.6. Glossary of Parameters Used in Models

N_i: N₁ and N₂ are the numbers of cells in compartments 1 and 2, respectively.

G0, G1, S, G2, and M are cellular phases: G₁ and S are in compartment N₁, and G₂ and M are in compartment N₂.

• ${\lambda }_i$: is the exchange between compartments N₁ and N₂.
• ${\alpha }_s^i$: is the tumor cell survival proportion after applying drugs s (DX plus CPh).
• $\gamma$: is the effect of PX on mitosis.
• ${\kappa }_i$: represents the portion of susceptible cells resistant because of drug s.
• ${\beta }_s^i$: parameter that indicates the decrease in the resistant population.

2.7. Quantitative Model

For the quantitative comparison of the chemotherapy sequencing strategies, we used a dynamic system describing the evolution of the different cell populations that contain the tumor. Based on the PLMM [13], we described the development of cell populations using a linear ODE, including the cancer treatment cycle. Although our model included the primary tumor and possible positive axillary nodes in TN-LABC, both were considered a single tumor volume, and cell proliferation was not limited to a specific geometry. Therefore, the model will not work with metastatic disease because cells spreading to other tissues have higher proliferation rates, and the response to the SS21 treatment was low.

Five phases of the cell cycle were described: quiescent cells (G0), cells that start synthesizing RNA and proteins (G1), DNA synthesis replication (S), proteins and RNA continue to be synthesized (G2), and mitosis (M). The G0 state is not properly part of the division process; thus, this phase was not considered in the PLMM. Instead, phases G1 and S were grouped into the $N_1$ compartment and phases G2 and M into $N_2$ (in this study, we defined $N_i, i=1, 2$ as the AMCP). The cell cycle followed an exponential growth pattern in each compartment, represented by a first-order ODE system:

$$\frac{d{N}_1}{dt}=-{\lambda }_1{N}_1+2{\lambda }_2{N}_2;\hspace{1em}\frac{d{N}_2}{dt}=-{\lambda }_1{N}_1-{\lambda }_2{N}_2$$

(1)

where ${\lambda }_1$ and ${\lambda }_2$ are the rate of exchange between compartments $N_1$ and $N_2$, and $2{\lambda }_2$ represents the M phase of the cell cycle. Under the Norton–Simon hypothesis [32], our model corresponded to the exponential growth phase (that is, we were far from the saturation point where medications have low effectiveness because of the tumor size). Therefore, the matrix representation of Equation (1) is defined as:

$$N=\left(\begin{array}{c}{N}_1\\ N_2\end{array}\right)\in {ℝ}^2$$

(2)

$$\frac{dN}{dt}=\frac{d}{dt}\left(\begin{array}{c}{N}_1\\ N_2\end{array}\right)=\left(\begin{array}{cc}-{\lambda }_1& 2{\lambda }_2\\ {\lambda }_1& -{\lambda }_2\end{array}\right)N=GN;\hspace{1em}G=\left(\begin{array}{cc}-{\lambda }_1& 2{\lambda }_2\\ {\lambda }_1& -{\lambda }_2\end{array}\right)$$

(3)

The solution for the whole linear system of the autonomous ODE is given by $exp\left(G\tau \right)$. The product of this matrix by the initial condition $N\left(0\right)$ provides the value of the population in period $\tau$:

$$N\left(\tau \right)=exp\left(G\tau \right)N\left(0\right)$$

(4)

where $exp\left(G\tau \right)={\displaystyle {\sum }_{i=0}^{\infty }}\frac{1}{i!}{\left(G\tau \right)}^i$. The growth in the $N_1$ and $N_2$ populations depends on parameters ${\lambda }_1$ and ${\lambda }_2$. Determining which value of these parameters corresponds to each subject is challenging. These can change because of numerous circumstances owing to individual differences in the cell proliferation speed (usually measured by the proliferation index Ki-67) [27,28]. Nevertheless, the growth of these cell populations can be estimated by sweeping a values interval of these parameters using clinical knowledge from published studies and comparing the development of the tumor cell population.

2.8. Doxorubicin and Cyclophosphamide Effect

A drug’s effect on treatment is described as an instant change in the AMCP when the dose is administered. The application of matrix $D_s$ to population N represents this behavior. The subscript s represents DX and CPh drugs: the type of chemotherapy that inhibits the replication of tumor cells. Following the PLMM, the matrix instantly acts on the N population; the matrix is diagonal because drugs independently modify each population.

$$D_s=\left(\begin{array}{cc}{\alpha }_s^1& 0\\ 0& {\alpha }_s^2\end{array}\right)$$

(5)

Here, ${\alpha }_s^i$, i = 1, 2, is the tumor cell survival proportion after applying drug s to AMCPs $N_1$ and $N_2$. The cell population affected by drug s in each treatment within period $\tau$ is expressed as follows:

$$N\left(\tau \right)=D_{s}exp\left(G\tau \right)N\left(0\right)$$

(6)

The effect of matrix $D_s$ depends on the time it is applied. Figure 3 shows that when $D_s$ is used at different times (${\tau }_1$, ${\tau }_2$, and ${\tau }_3$), the population change $\Delta N_i$, i = 1, 2 after chemotherapy can be negative, zero, or positive.

The time interval was limited in size because of the toxicity effect of $D_s$. Therefore, $\tau$ was considered, given this restriction. According to our chemotherapy sequencing model with s drugs, the minimum time interval to prevent toxicity was 14 days, reinforced with the GCSF.

Note that when m doses of drug s are administered, the resulting number of AMCPs is as follows:

$$N\left(m\tau \right)=(D_{s}exp \left(G\tau \right))(D_{s}exp \left(G\tau \right))\dots (D_{s}exp \left(G\tau \right))N\left(0\right)$$

$$N\left(m\tau \right)={(D_{s}exp\left(G\tau \right))}^{m}N\left(0\right)$$

(7)

where $N\left(m\tau \right)$ denotes the administration of different chemotherapy sequences by the corresponding matrix multiplication for ${\tau }_1=21-\mathrm{day}$ and ${\tau }_2=14-\mathrm{day}$ cycles (as defined in the chemotherapy sequences in Figure 1 and Figure 2).

2.9. Qualitative Behavior of Paclitaxel

One of the mechanisms of the action of PX is the inhibition of mitotic spindle formation during cell division, blocking the mitosis process. The addition of PX as matrix P, a modification of matrix G, was represented. PX influences the cell cycle growth matrix; therefore, this drug has continuous effects over time during its application.

Therefore, the AMCP N now takes the following form:

$$\frac{dN}{dt}=PN$$

(8)

We represented matrix P as follows:

$$P=\left(\begin{array}{cc}-{\lambda }_1& \gamma {\lambda }_2\\ {\lambda }_1& -{\lambda }_2\end{array}\right)$$

(9)

where coefficient $\gamma$ is the effect of PX on mitosis; $\gamma \in \left[0, 2\right]$ is considered. For $\gamma =2$, we have cell mitosis; for $\gamma =1,$ a bifurcation point exists; and for $\gamma <1,$ the eigenvalue is negative (see Appendix A) as the tumor size decreases to zero. Therefore, in this case, the effectiveness of PX was sufficient for the tumor to disappear. Additionally, the addition of CP on day 1 of each cycle is shown in Figure 1 and Figure 2, and the action mechanism through which PX destroyed the tumor cells was enhanced, for which the $\gamma$ value was even more reduced. Therefore, we implemented the CP effect in the mathematical model to increase parameter $\gamma$, which was a robust approximation of the effect of CP and PX in each simulation cycle.

We contrasted cell growth for SS21, SS14, IS21, and IS14 of the s drugs, including PX. The results of the combination of chemotherapy with PX are represented in Figure 4. The objective was to minimize the active mitotic cell fraction $N_i\left(\tau \right), i=1, 2$.

Figure 4a shows that in period τ, exponential growth was such that after applying $D_s$, it could not return to the initial population state, in which $\Delta N_i$ was positive, referring to tumor growth at a slower rate. However, PX could reduce the maximum growth of the $exp\left(G\tau \right),$ resulting in a new matrix $exp\left(P\tau \right)$, where the effect of chemotherapy reached the objective of reducing the initial population, for which $\Delta N_i$ was negative, as shown in Figure 4b. The above demonstrated that the tumor would have been almost or entirely reduced.

When $\Delta N_i$ is negative, repeating this chemotherapy process allowed the tumor to be effectively reduced, as shown in Figure 4c, where the reduction is shown after $m{D}_s$.

2.10. Resistance Model

When the drug resistance effect was incorporated into the model, the cell population was divided into chemotherapy-resistant and susceptible. Furthermore, a pharmacological effect was defined as when a cell group transitions from susceptible to resistant, and such conversion occurred instantly in our simulations. Following the same PLMM scheme, we have $\overline{N}\left(t\right),$ the AMCP with a vector of four components, that is, in ${ℝ}^4$.

$$\frac{d\overline{N}}{dt}=\left(\begin{array}{cc}G& 0\\ 0& G\end{array}\right)\overline{N}\left(t\right)=\tilde{G}\overline{N}\left(t\right)$$

(10)

Here, G is as defined in Equation (3), 0 is the null matrix, and $\overline{N}\left(t\right)$ is defined as

$$\overline{N}\left(t\right)=\left(\begin{array}{c}{N}_1^{se}\\ N_2^{se}\\ N_1^{re}\\ N_2^{re}\end{array}\right)$$

(11)

where $N_1^{se}$ and $N_2^{se}$ are the susceptible AMCPs and $N_1^{re}$ and $N_2^{re}$ are the resistant AMCPs. The fundamental matrix of this system is $exp\left(\tilde{G}\tau \right)$, and the AMCP with the effect of drugs ${\tilde{D}}_s$ after period τ is as follows:

$$\overline{N}\left(\tau \right)={\tilde{D}}_{s}exp\left(\tilde{G}\tau \right)\overline{N}\left(0\right)$$

(12)

In this case, matrix ${\tilde{D}}_s$ represents the population’s instantaneous change through the susceptible cell’s annihilation by drugs and the cell transformation from susceptible to resistant. Therefore, matrix ${\tilde{D}}_s$ takes the following form:

$${\tilde{D}}_s=\left(\begin{array}{cccc}\left(1-k_1\right){\alpha }_s^1& 0& 0& 0\\ 0& \left(1-k_2\right){\alpha }_s^2& 0& 0\\ k_1{\alpha }_s^1& 0& {\beta }_s^1& 0\\ 0& k_2{\alpha }_s^2& 0& {\beta }_s^2\end{array}\right)$$

(13)

where the ${\alpha }_s^i$ parameters are the same as those described in Section 2.8; parameters ${\beta }_s^i$, where i =1, 2, indicate a decrease in the resistant population (they represent the survival rate after applying drugs s to the N₁ and N₂ AMCPs); and the parameter $k_i$, where i = 1, 2, represents the rate of cell exchange speed for the susceptible to the resistant population. For example, the resistance model we described in Equation (12) depends on the $k_i$ parameters, where k₁ is the parameter of the $N_1$ population, and k₂ is the parameter of the $N_2$ population because of s drugs. As two inputs, ${\tilde{D}}_s\left(3, 3\right)={\tilde{D}}_s\left(4, 4\right),$ represent the effect of drugs on the resistant cells, and because such cells do not change owing to this effect, their numerical value is equal to 1. Therefore, the matrix of Equation (13) is as follows:

$${\tilde{D}}_S=\left(\begin{array}{cccc}\left(1-k_1\right){\alpha }_S^1& 0& 0& 0\\ 0& \left(1-k_2\right){\alpha }_S^2& 0& 0\\ k_1{\alpha }_S^1& 0& 1& 0\\ 0& k_2{\alpha }_S^2& 0& 1\end{array}\right)$$

(14)

In the results, Section 3.4 shows the simulations of numerical experiments with resistant and susceptible AMCPs under different parameters.

3. Results

3.1. Growth Rate Simulations of Tumor Cell Population

Figure 5 illustrates the rationale behind the lambda values that we used in our simulations and shows three graphs of the growth rate of the total tumor cell population for parameters ${\lambda }_i\in \left[0, 0.12\right],$ i = 1, 2. The unit of ${\lambda }_i$ is 1/day. We considered the interval reasonable compared with the possible growth. Using the linear ODE model, we presented the growth rate of the total AMCP ($N_f/N_i$, where $N_f=N_1\left(\tau \right)+N_2\left(\tau \right)$ was the final cell population, and $N_i=N_1\left(0\right)+N_2\left(0\right)$ was the initial AMCP) for end time of $\tau =25 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$. The purpose of timing was to observe any significant difference in AMCP growth within a time interval.

As shown in the graphs, the population growth rate $\left(N_1\left(\tau \right)+N_2\left(\tau \right)\right)/\left(N_1\left(0\right)+N_2\left(0\right)\right)$ was relatively homogeneous regarding parameters ${\lambda }_i$. ${\lambda }_2$ was more susceptible to population growth than ${\lambda }_1$ because ${\lambda }_2$ was the parameter that described final active mitosis. The yellow curve represents the ${\lambda }_i$ for which $N_f/N_i=2$ for $\tau =25 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$. Total tumor doubling time ranged from 60 to 175 days [35,36,37]; however, our model represented the growth rate of the tumor’s AMCP (10% to 20%) [38,39,40,41]. Appendix B shows that if the tumor doubling time was 100 days, then the AMCP rate of the tumor should double every 28.9 days. Because only the AMCP was considered in this study, we used a doubling time of 25 days. In these simulations, we used the same initial conditions: $N_1=500$ and $N_2=300$. Appendix C shows that the results obtained were independent of the initial conditions.

Figure 5 demonstrates the results from homogeneity; parameters ${\lambda }_1=0.1$ and ${\lambda }_2=0.05$ were considered throughout the study to ensure that the active cell fraction duplicated after 25 days. Appendix C shows that the qualitative behavior of the AMCP growth evolution was similar for several parameters because of the linearity of Equation (4).

3.2. Simulations of Effects of Doxorubicin and Cyclophosphamide

To perform a comparative analysis of cell growth, we defined the values of ${\alpha }_s^i$ for the s drugs. However, because we did not have a specific value for these parameters, we swept the possible oncological values of these parameters, so ${\alpha }_s^i\in \left[0, 1\right]$, i = 1, 2.

${\alpha }_s^i=1$ indicates that the drug does not affect the cells; ${\alpha }_s^i=0$ indicates that the drug annihilates the whole AMCP, as shown in Figure 6. The ${\lambda }_i$ of the G matrix parameters and the initial conditions were the same as those chosen in the previous section.

When performing the numeric simulations with the cell cycles, we first let the tumor cell growth evolve with no drug for five days (corresponding to the zero time). We exposed the AMCP to m = 4 cycles with s drugs (as shown in Figure 1 and Figure 2). The result of the growth rate with the effect of the s drugs, standard doses (defined in Section 2.1), and ${\tau }_2=14 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ is shown in Figure 6a, where the color bar represents the rate $N_{SS14}\left(t_f\right)/N_{SS14}\left(0\right)=\left(N_{\left(1\right)SS14}\left(t_f\right)+N_{\left(2\right)SS14}\left(t_f\right)\right)/\left(N_{\left(1\right)SS14}\left(0\right)+N_{\left(2\right)SS14}\left(0\right)\right)$ for $t_f=m{\tau }_2+5$.

When the sensitivity to medications was high, administering the same dose in a shorter time ${\tau }_2$ had a benefit up to an intermediate point of ${\alpha }_s^i\approx 0.5$ values because tumor growth was prevented. Drugs were highly effective in killing tumor cells at low values of ${\alpha }_s^i$. AMCP growth was small compared with that shown by the color bar (dark colors) in Figure 6a. Mathematically, reducing each drug period was convenient, not reaching an excessive increase in the dose to avoid intolerable toxicity.

The darkest part (black color) in Figure 6b is diagonal, indicating a more significant reduction in AMCP growth. As expected, SS14 treatment reduced the AMCP more than SS21 treatment, but the values of the medication parameters modulated this reduction ${\alpha }_S^i$, whereas those of SS14 produced a better response.

The two simulations showed the importance of drug administration in an optimal period to prevent tumor growth. In addition, the best SS14 treatment used GCSF to prevent high toxicity.

3.3. Simulations of the Paclitaxel Effect

Figure 7 illustrates the active mitotic cell population (the y-axis represents logarithmic values) for specific ${\alpha }_s^i$ values as a function of time. The active mitotic cell population (N1 + N2) evolved as a function of ${\lambda }_i$ to double this population in 25 days. The growth rate without drugs (red line) was exponential from an initial population of 800 million cells to reach an order of 129,825 million cells after 184 days, which was a doubling approximately every 25 days. Concerning the final cell population after the treatment with SS14, the population was 158 million cells (yellow line followed by dark blue line); after SI14, it was 125 million cells (green line followed by light blue line) after 184 days. SI14 killed the AMCP 21% more effectively than SS14.

Figure 8 shows a comparison of the final AMCP after treatments of SS14 and IS14, doing a complete sweep of the ${\alpha }_s^i$ interval, using ${\lambda }_i$ values to double the AMCP in 25 days.

The simulation results shown in Figure 8a show the difference in the population growth rate between the SS14 and IS14 treatments. In both cases, ${\tau }_2=14-\mathrm{d}\mathrm{a}\mathrm{y}$ cycles were considered for s drugs. We chose γ = 1.5 because it is a common value and is far from the bifurcation point for the cases in which the PX cycle is applied. The final total population result with the SS14 treatment was compared with that of IS14 $\left[N_{SS14}\left(t_f\right)/N_{IS14}\left(t_f\right)\right]$. The final versus initial AMCP growth rate of SS14 $\left[N_{SS14}\left(t_f\right)/N_{SS14}\left(0\right)\right]$ is shown in Figure 8b and of IS14 $\left[N_{IS14}\left(t_f\right)/N_{IS14}\left(0\right)\right]$ in Figure 7c to compare the growth order. As in previous simulations, the AMCP evolved without drugs for 5 days and m = 4 cycles so that $t_f=m{\tau }_2+5.$ Including PX therapy, in four cycles of 15 days (as shown in the second phase in Figure 1 and the first phase in Figure 2), $t_f$ remained equal to $t_f=4{\tau }_2+4\times 15+5$.

We can see in Figure 8a that by doing the complete sweep over the parameters ${\alpha }_S^1$ and ${\alpha }_S^2$, the resulting population in SS14 was more significant than that obtained in IS14. Furthermore, the area covered by values ${\alpha }_s^i$, corresponded to the greater efficiency of the IS14, whereas the SS14 occupied approximately 90%, which meant that IS14 was more efficient in containing tumor growth. We had a similar behavior if the value of ${\tau }_1=21 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$.

Figure 8b,c shows the tumor growth for SS14 and IS14, respectively. The color level represents how much the tumor grew compared with the initial population. The IS14 treatment, including the PX and CP, was found to be one of the best chemotherapies for patients with TN-LABC. Its population growth was minimal compared with that under the SS14 treatment [22,23,29]. Notably, the toxicity of the therapy had to remain within a specific limit.

3.4. Simulations of Resistance Model

The first simulation group came from an initial AMCP, with ${\lambda }_1=0.1$ and ${\lambda }_2=0.05$, used in Figure 5. As in the previous simulations, the matrix (14) ${\alpha }_S^i\in \left[0, 1\right], i=1, 2$ parameters represented the drug effects, and the unit intervals were partitioned into 100 equidistant portions. Figure 8 shows the results of five simulations of the tumor population rate with the effect of resistance to s drugs. The color bar corresponds to the final population rate versus the initial population $\left(N_{SS14}^{se}\left(t_f\right)+N_{SS14}^{re}\left(t_f\right)\right)/\left(N_{SS14}^{se}\left(0\right)+N_{SS14}^{re}\left(0\right)\right)$; the simulation time was $t_f=m\tau +5=4\times 14+5$.

When the plots of Figure 9 were analyzed, we observed homogeneity in the growing proportion of the resistant and susceptible populations as the $k_i$ value increased. The yellow line in the five charts shows that the final population was equal to the initial population, for which it took a value of one; this line separates the areas of the ${\alpha }_S^1$ and ${\alpha }_S^2$ parameters, where the population increased or decreased. The yellow line in Figure 9d,e represents no growth or a decrease in the active AMCP, closer to the axis origin, indicating that a more aggressive therapy should be chosen to prevent the tumor from becoming resistant to the drug.

In the following numerical experiment, we performed some simulations to compare the efficiency of the standard sequence versus the inverse sequence, considering the effect of s drug resistance. Figure 9 shows the results of the simulations of the growth rate of the final AMCP, SS14 versus IS14 $\left[\left(N_{SS14}^{se}\left(t_f\right)+N_{SS14}^{re}\left(t_f\right)\right)/\left(N_{IS14}^{se}\left(t_f\right)+N_{IS14}^{re}\left(t_f\right)\right)\right]$, for the same simulation time $t_f$ before and $k_i$ values as in Figure 9b–d. The extreme values of $k_i$ in Figure 9a,e were not considered because the smallest value produced almost no change from susceptible to resistant cells. The highest value corresponded to a quick change from susceptible to resistant cells. The values of the parameters (${\lambda }_1=0.1$ and ${\lambda }_2=0.05$) were the same as the previous simulations, and $\gamma =1.5$ was the same as that used in Figure 7 in Section 3.3, which was dedicated to PX.

When Figure 10a–c was analyzed for $k_1=k_2$, we observed that the growth rate of the tumor AMCP for SS14 versus IS14 was similar, with a slight variation, as shown in Figure 10c by the color bar. However, when comparing the $k_1\ne k_2$ values, a significant difference was observed in the efficiency of the two treatments. The conversion rate from susceptible to resistant cells differed between the $N_1$ and $N_2$ populations. Figure 10d–f shows that the performance of the IS14 treatment was more efficient than that of SS14, provided that $k_1>k_2$. Note that the variation in the population proportion in the color bar is higher than in Figure 10a–c, where rates above one predominated. At the right end of graph (f), we can observe a population below the yellow line.

A more detailed analysis of Figure 10 is shown in Figure 11; this figure compares the population growth rate simulations of the tumor cells, SS14 versus IS14, for three cases: $k_1>k_2$, $k_1<k_2$_, and $k_1=k_2$ where the exchanges of $k_1$ and $k_2$ values were analyzed when they were close.

The yellow line represents a value of 1; because the final populations of the SS14 and IS14 treatments were equal, the line shifts to the left when $k_1<k_2$ and $k_1=k_2$, respectively. When $k_1>k_2$, IS14 was more efficient than SS14. On the contrary, the yellow line divides the graph when $k_1<k_2$, and SS14 was more efficient than IS14. In Figure 11a, the area covered by the highest efficiency of the IS14 treatment is noticeably larger. In Figure 11b, the highest efficiency of SS14 barely exceeded that of IS14, and when $k_1=k_2$, we could not differentiate which of the two sequences was more efficient. Notably, the results were qualitatively similar for the 21-day cycles.

4. Discussion

In this study, we used a PLMM to describe the parameters that defined the order of four proposed oncological treatments. In addition, a uniform numerical sweep was performed for the parameters ${\lambda }_i\in \left[0, 0.12\right]$ and ${\alpha }_s^i\in \left[0, 1\right]$ to obtain a realistic tumor growth function. This approach produced an extensive qualitative and quantitative vision of the problem to explain the effects of the drugs on the AMCP.

Cancerous tumors actively divide cells during the G2 and M phases (active reproduction phase). In TN-LABC, the AMCP is between 20% and 60% of the tumor mass (43,44), which is used when calculating the doubling time of the tumor. In this study, the simulations of the tumor duplication were obtained for reasonable clinical times, $\tau =25 \mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ using the initial conditions of $N_1=500$ and $N_2=300$, as shown in Figure 5. Appendix C shows that the results obtained were independent of the initial conditions.

The matrix $D_s$, from Equation (5), represents the annihilating effect of s drugs on the AMCP. The drug effect was considered instantaneous concerning the timescales involved in the AMCP growth evolution. These timescales involved in sudden cell death were compared with the timeframe of AMCP growth without drugs. Proposing the killing effect of the medications as instantaneous was a suitable approximation, as illustrated in Figure 4.

We used a simplified mathematical model based on Roe et al.’s model (16), as mentioned throughout the text. One strength of this analysis were the exhaustive scenarios considered regarding ratios between cell populations and cellular dynamics by rates after treatments. However, one of this model’s limitations was using two populations instead of more populations that corresponded to the cell cycle, as Roe et al. did in their model. The AMCPs grouped the G1 (protein and RNA synthesis) and S (replication synthesis) cycles in population $N_1$ and the G2 (protein and RNA synthesis continues) and M (mitosis) cycles in population $N_2$; in this process, the state G0 (quiescent cells) was not considered. The reason for grouping the populations this way was that the drugs acted approximately homogeneously on the G1 and S cycles and the G2 and M cycles, which allowed us to handle a limited number of parameters. Although this was a limitation of the model, it was an advantage because we could study the global behavior of the model for a wide range of these parameters. Therefore, we obtained a broad qualitative and quantitative overview of the model.

Regarding the complexity of the tumor, owing to the model representing generic populations at any point in the tumor, by varying the model’s parameters, we could explore a wide range of cell types that considered this complexity. In the next phase of the study, we will add the effects of the immune system to the mathematical model.

A randomized phase II trial demonstrated that adding CP to weekly PX followed by DX and CPh significantly increased the pCR rate in stages II-III of TN-LABC [22]. Strategies to reduce the toxicity of sequential neoadjuvant chemotherapy in TN-LABC include adjusting the weekly CP dose, compensated with GCEF [29], adjusting by the age and comorbidity of patients to define the best treatment, determine adequate premedication, and prevent severe adverse reactions to PX.

The findings of Shepherd et al. [22] supported our mathematical model because the authors obtained higher PCR rates (administering first CP plus PX and then DX plus CPh) and provided helpful information for training our mathematical model with the specific parameters λ and α. As a result, our mathematical model was flexible enough to increase prediction accuracy with future adjustments using the poor clinical prognosis features of TN-LABC. These characteristics may include age less than 40 years, tumor size greater than 5 cm, presence of palpable axillary lymph nodes, inflammatory changes in the breast, histological characteristics (high histological grade tumors), and the immunohistochemical characteristics (triple-negative) of elevated Ki67, p53-positivity, and tumor-infiltrating lymphocytes ≥30. Comparing the findings of this study with our mathematical model, and we verified the consistency in the results obtained in the clinic, as observed in Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11, which showed that the IS14 treatment was more efficient than the SS14 treatment.

This type of PLMM swept all values of the computed parameters, helping to examine a wide range of drugs that can be used in clinical trials.

5. Conclusions

The sweep of all parameters ${\lambda }_i$ and ${\alpha }_s^i$ allowed us to visualize the variability in cell growth concerning variations in the amount of drug applied. Here, the inverted treatment sequence outperformed the standard treatment, allowing maximum AMCP and tumor size reductions. In addition, IS14 was more effective than IS21 but with a higher degree of toxicity; however, this can be compensated for with GCSF because some TN-LABC patients (with the basal-like subtype, BRCA1 gene mutations, p53 and ki-67, inflammatory carcinoma in individuals under 40 years of age, and other tumors with a high proliferation index) will need it, especially those who do not respond to IS21 treatment. The numerical results obtained in this study are consistent with the clinical results referenced in the manuscript. This study provides a complementary approach to that used in the clinic for optimal control.