Global dynamics of a colorectal cancer treatment model with cancer stem cells.

We present and analyze a mathematical model of the treatment of colorectal cancer using a system of nonlinear ordinary differential equations. The model describes the effectiveness of immunotherapy and chemotherapy for treatment of tumor cells and cancer stem cells (CSCs). The effects of CD8+T cells, natural killer cells, and interleukin proteins on tumor cells and CSCs under the influence of treatment are also illustrated. Using the method of localization of compact invariant sets, we present conditions on treatment parameters to guarantee a globally attracting tumor clearance state. Numerical simulations using estimated parameters from the literature are included to showcase various global dynamics of the model.

Introduction

Cancer is a leading cause of morbidity and mortality worldwide according to the World Health Organization (WHO) [1]. The WHO predicts that within the next two decades, the number of cancer cases will rise by approximately 70%. Of the 72 types of cancer documented by the American Cancer Society, colorectal cancer is reported to be the third most common cancer and third leading cause of cancer related death [2]. Though it can occur at any age, the risk for colorectal cancer increases after the age of fifty [3].

Numerous techniques have been implemented in modeling the growth and treatment of different types of cancer. Because the immune system has such a large impact on the growth of tumors, immune cell classes are usually considered separately in models to accurately capture the different interactions each immune cell type will have with the cancerous cells [4], [5], [6], [7]. A variety of treatment models, including both immunotherapy and chemotherapy, have also been implemented in a variety of publications [6], [8], [9]. In previous works, the key role of mutations and cell population dynamics for colorectal cancer have also been explored [10], [11]. The role of cancer stem cells (CSCs) in solid tumors continues to be a subject of interest in the medical literature, and several models have included them in an attempt to explain cancer recurrence [12], [13].

The cancer stem cell hypothesis, a recent theory for tumor growth, suggests a small population of cancerous cells known as cancer stem cells have stem cell-like qualities. These cells possess the ability to self-renew, resist apoptosis, differentiate, and resist certain treatments, such as chemotherapy. CSCs help explain cancer recurrence due to their ability to replenish the tumor cell population [14], [15], [16]. Because the CSC hypothesis is relatively new in application, incorporating CSCs into well-studied and tested models, and analyzing the global dynamics of such revised models, may provide new biological insights.

Modeling population dynamics at the cellular level is often accomplished using a system of nonlinear ordinary differential equations, and long-term behavior of such systems is analyzed using various local and global asymptotic stability methods. However, as the complexity of the model increases, standard approaches such as local linearization or Lyapunov functions become infeasible due to the large number of parameters involved, and only partial analytical results are possible. Recently, Starkov et al. [17], [18], [19], [20], [21], [22], [23], [24], [25] have successfully employed the method of localization of compact invariant sets to study a variety of cancer models in which standard approaches prove intractable.

The purpose of this paper is twofold: 1) to extend a colorectal cancer treatment model from de Pillis et al. [26] to include the dynamics of cancer stem cells, and 2) to study the global dynamics of the resulting system using the method of localization of compact invariant sets. We focus on the development of sufficient conditions on treatment terms to guarantee the existence of a globally attractive tumor clearance state. It is important to note that the global dynamics of the full system found in [26] have yet to be explored, but appear here as a subset of our extended model. The organization of the paper is as follows: Section 2 introduces the model and justifies new terms, Section 3 analyzes the global dynamics of the full model for a tumor clearance state, Section 4 explores numerical simulations of the model and discusses implications, and Section 5 closes the paper with concluding remarks.

Model

System of equations

In order to study the dynamics of the CSC hypothesis in colorectal cancer, we extend the model presented by de Pillis et al. [26] to include CSCs in the treatment of metastasized colorectal cancer. Since the focus of this paper is the global dynamics of the resulting nonlinear system of differential equations, we refer the reader to [26] for a full description of the original model. We briefly summarize below our additions to the de Pillis model and include references to justify the structure of each new term. In the model presented by de Pillis et al., the immune cell population is differentiated into various subgroups of: NK cells, ; CD+T cells, ; circulating lymphocytes, ; and interleukin proteins, . We shall invoke the same principle on the cancerous cell population by differentiating the tumor cell and the CSC populations, which we will denote by and respectively. Treatment of chemotherapy is represented by and immunotherapy by . All new terms are differentiated with bold text.

Note since we have differentiated the cancerous cell population into two cell groups, we include the term in multiple equations to denote the total cancerous cell population. The meanings and values of parameters are presented in Supplementary Materials A and B, respectively.

Justification of new terms

We model the structure of the CSC equation, and subsequent interactions between CSC and other cell populations, after the tumor cell dynamics from [26]. Let it be noted that CSCs develop at a faster rate than that of tumor cells but have a lower carrying capacity [15], [16]. Once the CSCs have reached their carrying capacity, they differentiate to maintain the tumor cell population [27]. NK cells preferentially attack CSCs [28]. The immunotherapy drug Cetuximab is more effective on CSCs than tumor cells, leading to a faster rate of NK induced death through ADCC and a faster rate of mAb-induced cell death. We also note that CD+T cells are preferential to binding to CSCs [29]. Due to the effectiveness of Cetuximab, the rate of mAb cell formation for CSCs is faster than tumor cells [29]. All these above CSC behaviors yield the following inequalities:

Definitive values for these CSC parameters are unknown at the moment. We base our values for numerical analysis on the parameters used by de Pillis, et al. [26]. We note that we are primarily interested in showcasing the possible dynamics of the model and will leave the fine-tuning of parameters as a subject for future research.

Also note that CSCs are more resistant to chemotherapy than normal tumor cells [16]. Published results have not yet provided a definitive value for the effectiveness of Irinotecan on CSCs. To account for this, we add the parameter τ to represent the percent effectiveness of the chemotherapy agent on CSCs in relation to tumor cells. This parameter will play a central role in the conditions for tumor clearance and, as such, is explored more fully in numerical simulations.

As the CSC hypothesis suggests, the growth of tumor cells is dependent on the CSC population. The term captures both that the tumor cells have limited growth potential and that CSCs drive their proliferation. Furthermore, the additional factor slows the CSC production of tumor cells and prioritizes their own regeneration. While there are likely many different potential structures for this growth term, we note that the one chosen here captures the qualitative behavior implicit in the CSC hypothesis while also admitting a complete global stability analysis as seen below. We include a natural death rate of ζ for tumors as well.

We also alter the immune system equations to account for the differentiation of the cancerous cell population. We note that both cancerous cell populations affect all the immune cells, therefore, we will consider the total population in terms of death due to exhaustion of cancer killing resources [30], [31].

It can clearly be seen that the vector field associated with (1) is continuously differentiable. Hence, given any set of initial conditions, there exists an unique solution. It is also clear that the non-negative orthant is invariant under (1).

Analysis

In this section, we determine the conditions under which we can guarantee a globally attractive tumor clearance state for our system of equations. To simplify our exploration of the dynamics of the model, we make the assumption that treatment is constant for both immunotherapy and chemotherapy. Although constant treatment is highly unrealistic, this assumption allows us to proceed with studying the global dynamics of a highly nonlinear system of differential equations. We consider more realistic treatments, such as periodic treatment, below in numerical simulations.

With the assumption of constant treatment terms, we denote the infusions of immunotherapy and chemotherapy as and respectively. To study the global dynamics of (1), we use the method of the localization of compact invariant sets as presented in Krishchenko [32] and Valle, et al. [17]. Before we proceed, we present some preliminary results that can be found in [17] and are included here for the benefit of the reader.

Consider a nonlinear system of the form where is a continuously differentiable vector field, and let Γ be a continuously differentiable scalar function that is not a first integral of (2). Further, let and denote and . We then define the localizing set of Γ to be , and we call Γ a localizing function of the system (2). The following theorem in [17] will be utilized in the argument below.

Given a localizing function Γ of the system (2), the localizing set contains all compact invariant sets of (2).

We note that examples of compact invariant sets include equilibrium points, limit cycles, chaotic attractors, homoclinic and heteroclinic orbits, and so on. Thus the construction of localizing sets gives useful information regarding the possible long-term dynamics of the model. We now present our main result:

Under the following conditions, (1) has a globally attractive tumor clearance state:

,

,

,

,

,

,

.

We shall construct these sufficient conditions for a globally attractive tumor clearance state of (1) using the method of the localization of compact invariant sets as follows. We first construct a localizing set containing all compact invariant sets of (1). We then establish the existence of a larger absorbing region that contains the localizing set in which all trajectories enter in finite time. Finally, we provide conditions on treatment parameters that guarantee is an attractive plane.

Localization of compact invariant sets

We begin developing our localizing set by applying the localizing function and note that implying

We then apply the same process, with , to find the bounds on the tumor cell population. We note that, using our upper bound on S from above, Solving this inequality for T yields . Hence,

We now consider our first treatment method, immunotherapy, with a localizing function of . It follows that and By setting the directional derivative to zero, we find

Defining , we see that due to the uncoupled nature of the M equation, it is immediate that where .

Using the bounds obtained on M above, we let and observe that which implies It follows that where .

Defining , we have Hence, where .

Applying yields Using bounds derived above, we arrive at the inequality Solving this inequality for L, yields This allows us to thus conclude that with denoted by where is given above.

Finally, we use the localization function and note that implies Thus,

We conclude that all compact invariant sets of (1) lie in the localization set

Absorbing region

We now provide conditions under which our localizing set is absorbing. We first note that the equations governing the dynamics of chemotherapy (M) and circulating lymphocytes (C) are uncoupled from the remainder of the system. Due to the simplicity of these equations, we can easily solve them for their globally attracting equilibrium solutions. From this, we can conclude that, in finite time, M will enter the set , where and . Likewise, C will enter the set in finite time.

Continuing with our argument that the localizing set is absorbing, we consider the candidate Lyapunov function and compute Utilizing the bounds for M and C obtained above, yields Completing the square for the first four terms and grouping strategically, we obtain

Thus, inequalities (H1)–(H6) provide sufficient conditions in which , for sufficiently large values of , and L.

Attractivity to plane

With the conditions given above, we have proven the existence of a compact global attractor. In order to conclude that we have tumor clearance, we must now argue that the planes and are attractive. Once again employing , we see that If (H7) is met, that is then and we can conclude that is attractive.

Note that this constraint is dependent on τ, suggesting that the effectiveness of the chemotherapy agent on CSCs is critical in the application of this combination treatment. Furthermore, with , using , it is easy to see that no further conditions are needed to ensure . Thus, we conclude that (1) contains a globally attractive tumor clearance state.

Results & discussion

In this section we examine numerical simulations of two treatment protocols. We first explore constant treatment protocols to demonstrate the existence of a globally attractive tumor clearance state as guaranteed by Theorem 2. We then further explore how the effectiveness of the chemotherapy agent on the CSC population affects how quickly this tumor clearance is achieved. Next, we turn our attention to more realistic periodic treatment methods. While rigorous analytical results are more difficult to establish for such protocols, we may still gain insights on possible stable tumor clearance states by exploring numerical simulations in these cases. It should be noted that all simulations in this section were performed using Mathematica 10.

Constant treatment

CSC chemotherapy effectiveness and tumor clearance

Recall that the parameter τ represents the effectiveness of the chemotherapy agent on the CSC population. While it is known CSCs are resistant to some chemotherapy drugs, such as Irinotecan, the degree of resistivity is currently unknown. From Theorem 2, our model supports that when chemotherapy is ineffective on CSCs, tumor clearance relies on the effectiveness of immunotherapy. As the effectiveness of chemotherapy on CSCs increases, the amount of immunotherapy required for tumor clearance decreases. As depicted in Figure 1, in order to reach a cure state in 50 days with no immunotherapy and a constant chemotherapy infusion rate [26], we must meet the condition , meaning chemotherapy must be at least 95% as effective on CSCs as on tumor cells.

Figure 1

Log plot of cell count vs. CSC chemotherapy effectiveness using constant chemotherapy at an infusion rate of v = 3.6217 mg/L.

However, with the inclusion of immunotherapy, the required level of effectiveness of chemotherapy on CSCs to guarantee tumor clearance decreases. Using a baseline treatment level [26], we see in Figure 2 that with combination therapy, the chemotherapy agent needs only be 80% effective on CSCs.

Figure 2

Log plot of cell count vs. CSC chemotherapy effectiveness using constant chemotherapy at an infusion rate of v = 3.6217 mg/L and immunotherapy at an infusion rate of v = 0.76515 mg/L.

We note that the reduced dependence on the effectiveness of chemotherapy on CSCs when immunotherapy is introduced directly relates to the attractivity argument in Section 3.3. In particular, condition (H7) reveals that when chemotherapy is not effective on CSCs (), guaranteeing tumor clearance relies solely on the dosage of immunotherapy. However, when chemotherapy is effective on CSCs (), a lower dosage of immunotherapy is required to satisfy the inequality.

Impact of τ on long-term cell dynamics with insufficient immunotherapy

As shown in Figure 2, the value of τ plays a significant role in how effective the combination treatment strategy is at eradicating cancer. We further explore the role that τ plays by simulating constant combination treatment and noting the effects changes in τ have on long-term cancer dynamics. Continuing with the infusion rates from above, Figure 3 demonstrates that a clearance state is never reached when chemotherapy is ineffective on CSCs ().

Figure 3

Log plot with τ = 0 using constant treatment with infusions of v = 3.6217 mg/L for chemotherapy and v = 0.76515 mg/L for immunotherapy.

We observe that since chemotherapy is only effective on the tumor cells, the tumor cell population is lower than the CSC population. Note that CSCs reach their carrying capacity at approximately day 10. We also note that the tumor cell population repopulates following day 6 and stabilizes at the same time that the CSCs reach their carrying capacity. The tumor cells never reach their carrying capacity because the constant chemotherapy treatment continuously limits their population.

To contrast this behavior with a chemotherapy drug that is effective on CSCs, we choose a representative value of τ above the 80% threshold seen in Figure 2 that will lead to a cure state. For instance, when chemotherapy is 85% as effective on CSCs, a cure state will be reached, as depicted in Figure 4. We observe that the NK cell, interleukin, and circulating lymphocyte populations are unchanged by the effectiveness of the chemotherapy. With , a cancer clearance state is reached at approximately day 33. Since chemotherapy is still not as effective on CSCs as on tumor cells, it takes additional 12 days for the elimination of CSCs following the eradication of tumor cells. Thus, if the immunotherapy treatment is insufficient, the effectiveness of chemotherapy on CSCs is vital in reaching a cure state.

Figure 4

Log plot with τ = 0.85 using constant treatment with infusions of v = 3.6217 mg/L for chemotherapy and v = 0.76515 mg/L for immunotherapy.

Impact of τ on long-term cell dynamics with sufficient immunotherapy

Let us now consider the case in which chemotherapy is ineffective on CSCs, but the patient is receiving a sufficient amount of immunotherapy to guarantee a cure state, as depicted in Figure 5. We now reach a tumor clearance state at approximately day 49. The NK cell, interleukin, and circulating lymphocyte populations are unaffected by immunotherapy treatment. The two different cancerous cell populations are eradicated approximately three weeks apart.

Figure 5

Log plot with τ = 0 using constant treatment at an infusion rate of v = 3.6217 mg/L for chemotherapy and v = 3.4 mg/L for immunotherapy.

Next, we consider increasing the value of τ so that chemotherapy is 30% as effective on CSCs as on tumor cells, depicted in Figure 6.

Figure 6

Log plot with τ = 0.3 using constant treatment at an infusion rate of v = 3.6217 mg/L for chemotherapy and v = 3.4 mg/L for immunotherapy.

Increasing the value of τ only affects the two cancerous cell populations. The time to reach a cancer clearance state is drastically decreased to approximately 24 days, with only a six day difference between the eradication of tumor and CSC populations. Hence, we may conclude that the more effective the chemotherapy is on CSCs, the faster a cure state is achieved. The time between the eradication of the two cancerous cell populations also decreases with greater values of τ, reducing the likelihood of cancer recurrence.

Periodic treatment

Our globally attractive tumor clearance state above has been found analytically with constant treatment. However, in practice, most treatments are given out in shorter doses on a regular basis. In this subsection, we explore cancer dynamics under periodic treatment. For our model, we will suppose Irinotecan treatment, , will be given weekly for two hours in a dose of 57.947 mg/L/day. We will also suppose that Cetuximab treatment, , will be given for two hours with a loading dose of 139.072 mg/L/day followed by 173.840 mg/L/day given weekly. The treatment values chosen are the highest dosages that have been seen clinically to ensure a cure state [26]. Both treatment regimens are given for four weeks, followed by a two week rest period before introducing the next round of treatment. For practical purposes we assume that the cancer cell population is fully eradicated when the number of cells drops below 1.

Periodic combination treatment with low CSC chemotherapy effectiveness

Using the initial conditions found in de Pillis et al., , , , , , , , we model a healthy patient with a large tumor [26]. We assume CSCs are 2% of the total cancerous population for colorectal cancer, as seen in Todaro et al. [16]. When the chemotherapy is 30% as effective on CSCs as on tumor cells (), it takes around 30 days for the CSC population to be eliminated, as can be seen in Figure 7.

Figure 7

Log plot with τ = 0.3 and periodic treatment of v(t)=57.947 mg/L/day for chemotherapy and v(t)=139.072 mg/L/day as a starter dose for immunotherapy and v(t)=173.840 mg/L/day for proceeding dosages.

We have shown analytically that the tumor clearance state for this model using constant treatment is dependent on τ. This relationship is present in periodic treatment simulations as well. Tumor cells appear to plateau between treatments, while CSCs are able to slightly grow back before declining rapidly once tumor cells are eradicated. The different cell populations within the immune system seem to be approaching periodic solutions.

Periodic combination treatment with high CSC chemotherapy effectiveness

As we increase the efficiency of chemotherapy, as seen in Figure 8, the cancerous population reaches a cure state more rapidly.

Figure 8

Log plot with τ = 1 and periodic treatment of v(t)=57.947 mg/L/day for chemotherapy and v(t)=139.072 mg/L/day as a starter dose for immunotherapy and v(t)=173.840 mg/L/day for proceeding dosages.

When chemotherapy is as effective on CSCs as on tumor cells, the period between the eradication of the two cancerous cell populations is slightly decreased, but the total time to clear cancer from the system is drastically decreased. It should be noted that a decrease in the amount of CSCs, perhaps caused by a more effective chemotherapy agent, would also aid in depleting the tumor cell population, since the growth of tumor cells relies on CSCs. We have seen, analytically and numerically, that the speed at which a combination treatment eradicates the cancerous population relies on the effectiveness of the chemotherapy agent on CSCs. This shows the importance of discovering chemotherapy agents that are effective at attacking CSCs as well as tumor cells. Chemotherapy agents that are effective at fighting CSCs could lead to lower necessary doses of treatments, as well as quicker cancer eradication.

Conclusion

The model presented here aims to provide a framework for understanding the role that CSCs play in the dynamics and treatment of colorectal cancer. While analyzing the global dynamics of such a highly nonlinear model is usually quite difficult, we note that the ability to construct a globally attractive tumor clearance state of (1) showcases the usefulness of the method of localization of compact invariant sets developed in Starkov et al. [17], [18], [19], [20], [21], [22], [23], [24], [25]. To the authors' knowledge, there have been no attempts to consider analytically the global dynamics of the full model presented in [26]. Here, the complexity of system (1) prohibits the identification of specific globally asymptotically stable equilibrium solutions, yet the method we employ still allows for the creation of sufficient conditions on model parameters to guarantee the existence of a general tumor clearance state. We leave open the resolution of all possible dynamics on the remaining subsystem of (1) with , but we are hopeful that the method presented here can serve as a starting point for exploring the global dynamics of other biologically complex models which have not previously admitted sufficient analysis. Furthermore, our numerical simulations reveal the potential significance of designing CSC-targeted chemotherapy agents to guarantee tumor clearance and avoid cancer recurrence.

This work can be furthered by seeking conditions under which periodic solutions exist and are globally stable. Additionally, since we have shown the significance that τ has on balancing levels of chemotherapy and immunotherapy to achieve tumor clearance, incorporating a general healthy cell population to monitor the detrimental effect chemotherapy has on the overall health of a patient is an important extension of the work presented here. Finally, consideration of different submodels for possible cancer persistent states or other specific chemotherapy and immunotherapy agents could also be numerically assessed with this model.

Declarations

Author contribution statement

Kristen Abernathy, Zachary Abernathy, Kelsey Brown, Claire Burgess, Rebecca Hoehne: Conceived and designed the analysis; Analyzed and interpreted the data; Contributed analysis tools or data; Wrote the paper.

Funding statement

This work was supported by the National Science Foundation (Grant No. DMS-1358534). This work was supported by an Institutional Development Award (IDeA) from the National Institute of General Medical Sciences from the National Institutes of Health (2 P20 GM103499 15).

Competing interest statement

The authors declare no conflict of interest.

Additional information

Supplementary content related to this article has been published online at 10.1016/j.heliyon.2017.e00247.

No additional information is available for this paper.

a < b	K1 < K2
c < ν	ξT < ξS
ψT < ψS	sS < sT
λT < λS.