Comparing theory based and higher-order reduced models for fusion simulation data

1.   Introduction

Cancer is a group of diseases involving abnormal cell growth ^[1,2]. Currently, despite great progresses of many newly developed therapeutic methods ^[3,4,5], chemotherapy is still a common treatment method for many cancers. Patients are often administrated with high-dose of cytotoxic chemotherapeutics trying to eliminate tumor cells as much as possible ^{[6,7,8,9,10,11]}. Nevertheless, it is difficult to determine the proper dosage of chemotherapy, low dosage is ineffective in killing tumor cells, whereas excessive dosage may result in additional toxicity that is intolerable to patients ^[6]. Clinically, patients are often dosed at maximum or near maximum tolerated dose, which is carefully determined in phase Ⅰ studies ^[7,12]. However, high level doses often induce series side-effects, increasing the chemotherapy dose (also the treatment cost) would not yield the decreasing of the recurrence rate ^[13,14]. The recurrent tumors often show drug resistance, which is a major cause of treatment failure in chemotherapeutic drugs ^[15,16,17].

Drug resistance has been a major challenge in cancer therapy. The mechanisms of drug resistance are complex, and many reasons are involved, including cellular plasticity ^[18], heterogenous tumor cells ^[19], or therapy induced gene mutations ^[20,21]. In this study, we focus on a mechanism of drug resistance due to chemotherapy-induced genome instability. Chemotherapy agents are cytotoxic by means of interfering with cell division in a way of damaging or stressing cells, and leading to cell death through apoptosis. During the early stage of apoptosis, apoptotic chromosome fragmentation (C-Frag) are produced through the cleavage by caspase-3 activated DNase (CAD) ^[22]. Nevertheless, C-Frag does not always result in cell death, sometimes the chromosome fragments can randomly rejoin to form genome chaos so that the cells survive from crisis ^[23,24]. These survived cells carry non-clonal chromosome aberrations (NCCAs), the major form of genome variation and the key index of genome instability in cancer cells ^{[25,26,27,28]}. It was proposed that such genome instability induced by chemotherapy is a source of drug resistance in cancer therapy ^[26,27]; quantitative control of drug dosages is important for the long-term clinical effects.

Mathematical modelling approaches have been widely used in cancer research from different aspects ^[2,29]. A variety of models have been established to study the mechanisms of drug resistance after chemotherapy ^{[30,31,32,33,34]}. However, there are rare quantitative studies on how tumor cells population change in response to chemotherapy and drug resistance due to therapy-induced chromosome recombination. The roles of NCCAs in tumor recurrence is still controversial ^[35].

In this study, we intend to investigate how NCCAs may affect tumor growth, and present a mathematical model for chromosome recombination-induced drug resistance in cancer therapy. The model extends the previously well studied G0 cell cycle model ^{[36,37,38,39]}, and includes cell survival from C-Frag ^[23,24]. We mainly study cell population responses to various doses of chemotherapy, and show that there is an optimal dose (within the maximum tolerated dose) so that the steady state tumor cell number is relative low after chemotherapy. Moreover, the model implies that persistent extreme high dose therapy may induce oscillations in cell number, which is clinically inappropriate and should be avoided.

2.   Model formulation

In this study, we model the process of tumor growth through formulations of stem cell regeneration. Here, we mainly consider chemotherapy for leukemia, and mathematical models of hematopoiesis are referred in our modeling. We refer the classical G0 cell cycle model that has long been studied in literatures ^[40,41,42] (Figure. 1(a)). In the model, leukemia stem cells are classified as either resting or proliferative phase cells. Resting phase cells ($Q$, $\mathrm{cells/kg}$) either enter the proliferative phase at a rate $\beta$ ($\mathrm{day}^{-1}$), or be removed from the pool of resting phase due to differentiation, senescence, or death, at a rate $\kappa$ ($\mathrm{day}^{-1}$). The cells at the proliferative phase undergo apoptosis in a rate $\mu$ ($\mathrm{day}^{-1}$), and the duration of the proliferative phase is $\tau$ ($\mathrm{days}$), each survived cell divides into two daughter cells through mitosis at the end of the proliferative phase. These processes can be described by a delay differential equation model ^{[36,37,38,39]}:

Here, the subscript means the time delay, i.e., $Q_\tau = Q(t-\tau)$. The equation (2.1) has been widely applied in the study of hematopoietic stem cell regeneration dynamics and blood disease ^{[43,44,45,46]}, as well as the hematopoietic responses to chemotherapy ^[47,48].

The proliferation rate $\beta$ is a function of the resting phase cell number $Q$, indicating the regulation of cell proliferation through cytokines secreted from all stem cells. Normally, the proliferation rate is a decrease function of the cell number, and approaches $0$ when the cell number $Q$ is large enough ^[50]. Nevertheless, in the situation of cancer, the function can be non-monotonic because cancer cells can evade growth suppressors, and produce self-sustaining proliferative signaling ^[51].

Now, we consider the effects of chemotherapy, which often promote cell death during the proliferative phrase due to the toxicity ^[49]. Hence, we write the apoptosis rate as $\mu = \mu_0 + \mu_1$, where $\mu_0$ represents the baseline apoptosis rate in the absence of chemotherapy, and $\mu_1$ the extra apoptosis rate due to treatment stress (the rate $\mu_1$ is often increase with the chemotherapy dose, hence we also refer $\mu_1$ as the dose for short). At the early apoptosis stage, the cells undergo chromosome fragmentation (C-Frag) (cell number given by $(1-e^{-\mu\tau}) \beta(Q_\tau)Q_\tau$). However, in a small population of cells with C-Frag, the fragments can rejoin to yield chromosome recombination and the cells survive from crisis and re-enter the G0 phase; other cells continue the apoptosis process (Figure. 1(b)). We assume that the probability of chromosome recombination, $q(\mu_1)\ (0\leq q(\mu_1) < 1)$, is dependent on the chemotherapy dose, so that $q(\mu_1)$ is an increase function. Therefore, we modify the above G0 cell cycle model (2.1) to include chromosome recombination, the cell number $Q$ satisfies the following delay differential equation:

where

Hereafter, we always assume that $\beta(Q)$ is a decrease function, and $q(\mu_1)$ is an increase function.

For model simulation, we refer the classical models for hematopoietic stem cells (Table 1), and take the proliferation rate $\beta(Q)$ and the probability of chromosome recombination $q(\mu_{1})$ as Hill type functions ^[44,50]:

3.   Results

3.1.   Dependence of the steady state with drug dosage

From (2.2), the steady state $Q(t)\equiv Q^*$ is given by the equation

Obviously, there is a zero solution $Q^* = 0$. When the proliferation rate $\beta(Q)$ is a decrease function, there is a positive steady state $Q^* > 0$ if and only if the condition

is satisfied, and the steady state is given by the root of

From (3.1), there exists a positive steady state when the recombination rate $q(\mu_1)$ satisfies

Specifically, if

the system has only zero solution steady state, and the zero solution is global stable which means the situation of no cells; and if

the zero solution is unstable and all positive solutions approach to $+\infty$, which means the situation of uncontrolled cell growth.

Next, we consider the stability of the steady states. Let $x(t) = Q(t)-Q^*$ and linearize the equation (2.2) at $x = 0$, we obtain

where

The zero solution of equation (3.6) is stable if and only if the coefficients $a$ and $b$ take values from the region $S$ defined as

For the zero solution $Q(t)\equiv Q^{*} = 0$, we have

Hence, the zero solution is stable if and only if $b < a$, i.e., the inequality (3.4) is satisfied.

For the positive steady state $Q(t) = Q^{*} > 0$, we have

where $\bar{\beta} = \beta(Q^*) > 0$ and $\hat{\beta} = -\beta'(Q^*)Q^* > 0$. Hence, when the condition (3.1) is satisfied, the equation (3.7) gives the Hopf bifurcation curve for the positive steady state

where $\omega$ can be solved from the second equation of (3.9), and

When $\hat{\beta} < \hat{\beta}_{\mathrm{crit}}$, the positive steady state solution of equation (2.2) is stable, and when $\hat{\beta} > \hat{\beta}_{\mathrm{crit}}$, the steady state becomes unstable. Biologically, the Hopf bifurcation curve (3.9) gives the critical proliferation rate $\hat{\beta}_{\mathrm{crit}}$, so that the steady state is stable when the proliferation rate is less than the critical rate, but when the proliferation rate is larger than the critical rate, the steady state becomes unstable and the system exhibits oscillatory dynamics, which can be a source of dynamical blood diseases ^[42].

In summary, we have the following conclusion:

Theorem 1. Consider the model equation (2.2), where $\beta(Q)\ (\beta_0 > \beta(Q) > \beta_1)$ is a decrease function, and $0\leq q(\mu_1) < 1$. The equation always has zero steady state $Q(t) \equiv 0$; if and only if the condition

is satisfied, equation (2.2) has a unique positive steady state solution $Q(t) = Q^*$, which is given by

Moreover,

(1) the zero steady state is stable if and only if

(2) when (3.11) is satisfied and let $\bar{\beta} = \beta(Q^*)$, $\hat{\beta} = - \beta'(Q^*) Q^*$, for any $\mu_1 > 0$, there is a critical proliferation rate $\hat{\beta}_{\mathrm{crit}} > 0$, defined by (3.9), so that the positive steady state is stable if and only if $\hat{\beta} < \hat{\beta}_{\mathrm{crit}}$.

When the functions $\beta(Q)$ and $q(\mu_1)$ are defined by (2.4), equation (3.2) gives

and

Applying the default parameter values in Table 1, when $\mu_1$ varies from $0$ to $1.2$, the critical proliferation rate $\hat{\beta}_{\mathrm{crit}}$ (black line, given by (3.9)) and proliferation rate at the steady state $\hat{\beta}$ (red line, given by (3.15)) are shown in Figure 2(a). Here, $\hat{\beta} > \hat{\beta}_{\mathrm{crit}}$ implies the parameter region with unstable steady state, and there are oscillatory solutions due to Hopf bifurcation (Figure 2(b)).

From Figure 2(a), when the drug dose is low, we have $\hat{\beta} < \hat{\beta}_{\mathrm{crit}}$, and (2.2) has a stable positive steady state (Figure 2(b), blue line). The steady state cell number depends on the dose $\mu_1$, and reach a local minimum when $\mu_1$ takes intermediate values (green dot in Figure 2(a), also referred to Figure 3). When the drug dose further increase to a high level (black dot), the steady state of (2.2) becomes unstable, and the cell populations show oscillatory dynamics (Figure 2(b), black line).

3.2.   Optimal dose of chemotherapy

From the above analyses, when

the equation (2.2) has a unique positive steady state $Q(t) = Q^{*}$. Moreover, let $\bar{\beta} = \beta(Q^*)$, then $\bar{\beta} = \kappa/f$, and

Hence, the steady state can be expressed explicitly as

Equation (3.17) shows that the steady state $Q^*$ is dependent on the chemotherapy dose parameter $\mu_1$ through $\bar{\beta} = \kappa/f$. The following theorem shows that under certain conditions, there is an optimal dose (within a tolerated dose) so that the steady state cell number reach a local minimum.

Theorem 3.2. Consider the equation (2.2), if the following conditions are satisfied

(1) the functions $\beta(Q)$ and $q(\mu_1)$ are continuous, and satisfy

(2) the parameters ($\beta_0, \beta_1, \kappa, \mu_0$, and $\tau$) satisfy the condition (3.16) (here we note $\mu = \mu_0 + \mu_1$),

(3) both $q(0)$ and $q'(0)$ are sufficiently small,

(4) there exists $\mu_1^* > 0$ so that conditions (1)-(2) are satisfied when $\mu_1 \in (0, \mu_1^*)$, and

the steady state $Q^*$ reaches a local minimum value at an optimal dose $\mu_1 = \hat{\mu}_1 \in (0, \mu_1^*)$.

Proof. When the conditions (1) and (2) are satisfied, equation (2.2) has a unique positive steady state $Q(t)\equiv Q^*$, which is given by (3.17). Hence, we have

and

Hence, the sign of $\dfrac{d Q^*}{d \mu_1}$ is determined by the sign of $\dfrac{d f}{d \mu_1}$.

The condition (3) implies $\left.\dfrac{df}{d\mu_1}\right|_{\mu_1 = 0} < 0$, and condition (4) implies $\left.\dfrac{df}{d\mu_1}\right|_{\mu_1 = \mu_1^*} > 0$. Hence, there exists an optimal value $\hat{\mu}_1 \in (0, \mu_1^*)$, so that $\left.\dfrac{d Q^*}{d \mu_1}\right|_{\mu_1 = \hat{\mu}_1} = 0$, thus $Q^*$ reaches a local minimum value at $\mu_1 = \hat{\mu}_1$.

Figure 3 shows the dependence of steady state with the dose $\mu_1$ when we take parameters from Table 1, which exhibits a local minimum at about $\hat{\mu}_1 = 0.3$. Here we note that the steady state cell number $Q^*$ decrease with further increasing of $\mu_1 > 0.5$ and approaches zero when $\mu_1$ is large enough. However, clinically, higher level $\mu_1$ may exceed the maximum tolerated dose, and have a risk to induce oscillatory dynamics (Figure 2, and the discussion below). Hence, there is an optimal dose within a certain region (for example, $\mu_1 < 0.5$) so that the cell number reaches a relative low level.

3.3.   Chemotherapy-induced oscillations

Above analyses show that chemotherapy may induce oscillatory dynamics through Hopf bifurcation when $\mu_1$ takes proper values. Here, we further analyze the condition for Hopf bifurcation, and identify how clinical conditions may affect the chemotherapy-induced oscillations ^[46]. Hematopoiesis can exhibit oscillations in one or several circulating cell types and show symptoms of periodic hematological diseases ^[43,53]. Here, we show the conditions to induce oscillations, and try to identify the criteria to avoid therapy-induced oscillations in cancer treatments.

Using the previously introduced notations $\bar{\beta}$ and $f$ in (3.14), the bifurcation curve (3.9) can be rewritten as

and hence

Moreover, from equation (3.15), we have

Hence, from (3.19)–(3.21), and the bifurcation condition $\hat{\beta}_{\mathrm{crit}} = \hat{\beta}$, we obtain the bifurcation curve in terms of $f$ and $\kappa$ (defined by a parameter $\omega > 0$) as

Thus, given the parameter $\beta_0, \beta_1, n, \tau$, equation (3.22) defines a bifurcation curve in the $\kappa$-$f$ plane when $\omega$ varies over the region in (3.10).

Figure 4 shows the bifurcation curves obtained from (3.22). For any $\kappa > 0$, there exist an upper bound $\overline{f}$ and a lower bound $\underline{f}$, so that the positive steady state (if exists) is stable when $\underline{f} < f < \overline{f}$. Moreover, when the self-sustained proliferation rate $\beta_1$ increases, the lower bound $\underline{f}$ is independent of $\beta_1$, while the upper bound $\overline{f}$ decreases with $\beta_1$. Here, we note that $f$ depends on the dose parameter $\mu_1$ through (3.14), these results provide a strategy to quantitatively control the chemotherapy dose to avoid therapy-induced oscillation, i.e., try to take values from the gray region for different differentiation rate $\kappa$.

4.   Discussion

The era of chemotherapy has been ruled by the routine use of dose-intense protocols based on the "maximum-tolerate dose" concept. This protocol plays a prominent role in veterinary oncology, however there are many debates on using the high dose chemotherapy because of the side effects and recurrence rates ^[54]. Chromosome recombination induced by chemotherapy is a source to produce chaotic genome in cancer cells and may lead to drug resistance. Here, we present a mathematical model to consider cell population dynamics in response to chemotherapy with occasional occurrence of chromosome recombination. Model analyses show that maximum-tolerate dose may not always result in the best outcome when the probability of chromosome recombination is dependent on the dose; there is an optimal dose within a tolerated dose so that the steady state tumor cell number reaches a relative low level. Moreover, sustained administration of high dose chemotherapy may induced oscillatory cell population dynamics through Hopf bifurcation. Clinically, the long period oscillation can be considered as the frequent recurrence of tumor cells, which is a result of drug-resistant or chemotherapy-induced dynamical disease. For instance, the periodic chronic myelogenous leukemia show obvious oscillations in circulating blood cells ^[55]. We identify the parameter regions for the occurrence of oscillation dynamics, which is valuable when we try to avoid the frequent recurrence. To our knowledge, this is the first try to quantitatively study the tumor cell population dynamics in response to chemotherapy when chromosome recombination is induced. Genome chaotic has been widely reported in studies of genomic structure of cancer cells, however how genome chaotic play roles in cancer development and drug resistance is not well documented. The current study is the first try to consider this issue. Here, the proposed model only consider the main effect of cell survival from C-Frag through chromosome recombination. Nevertheless, many other effects should be included for a complete understanding of drug resistance due to the chromosome changes induced by cancer therapy. Moreover, the cell heterogeneity can also play important roles in drug resistance. These effects should be extended to the current model in order to obtain an optimal dose of chemotherapy.

Acknowledgments

This work was supported by grant from National Natural Science Foundation of China (91730301, 11762011).

Conflict of interest

No potential conflicts of interest were disclosed.