Optimal control model of tumor treatment in the context of cancer stem cell

1.   Introduction

Recent advances in oncology have demonstrated the existence of a subpopulation of tumor cells called cancer stem cells (CSCs) that are crucial in the initiation, metastasis, and resistance of tumors ^[1,2]. They have been identified in many types of cancers. For example, in immunodeficient mice, orthotopic xenografts of CSCs isolated from brain tumors were performed and demonstrated tumorigenesis properties ^[3,4,5].

CSCs constitute a small cell population within a tumor, divide extremely slowly, are capable of self-renewal, and can explain the genesis of a tumor. The heterogeneity observed in tumors is explained as the result of a partial differentiation gradient of different CSCs. These cells possess distinct genetic and epigenetic characteristics that can resist therapy and may later initiate metastatic processes ^[6,7]. Previous studies have proven that increased genetic instabilities in normal stem cells may lead to the formation of CSCs, suggesting that stemness is attained from additional genetic modifications. Currently, biological markers are available to identify these stem cells such as CD34+/ CD38 in leukemia, CD44+/CD24 in breast cancer, CD133+ in brain cancer, colorectal lung cancer, and endometrial cancer, CD44+/CD24+ in pancreatic cancer, CD44+/CD117+ in ovarian cancer, CD44+/CD271+ in throat cancer, CD90 in liver cancer, CD105 in kidney cancer, and CD271 in melanoma ^{[8,9,10,11,12]}.

Conventional therapies treat all tumor cells similarly. However, CSCs are not always affected. Indeed, they are generally, highly chemo-resistant, and radio-resistant; therefore, recurrence is almost inevitable in the near or long term. Targeting CSCs is a major challenge in cancer treatment and might limit tumor growth and prevent metastasis, offering new hope in anti-cancer treatment (trials have been conclusive in the fight against melanoma ^[13]). Several approaches are used in that measure, like targeted therapy that aims to sensitize CSCs to radiotherapy ^[1] or also using metronomic scheduling (MSAT) to affect both the cells and their micro-environment ^[14,15,16]. The problem can also be treated with conventional therapies only with new administration schedules ^[17]. This approach can be considered very close to the MSAT just in a more general way. Targeted therapy has also presented a way of enhancing the results of conventional treatments ^[18]

In this context, based on the CSC niche model of Jinzhi Lei ^[19], we propose a CSC mathematical model that describes the growth of a heterogeneous population of cancer cells and its treatment response as an optimal control problem. Our model contains two compartments of cancer cell populations. The first describes the dynamics of different types of CSCs. The second describes the dynamics of different types of non-stem cancer cells (CCs). Based on the hypothesis that stem cells are at the origin of cancer, we have assumed that each type of CC is an early progenitor of a certain type of CSC. These CCs can divide rapidly but in a limited manner in opposition to CSCs. Our goal is subsequently to identify an optimal multidrug scheduling that aims to stabilize the entire cancer cell population at a certain competitive equilibrium, by targeting CSCs. This approach allows transforming cancer into a chronic disease rather than a lethal one. The idea is based on using a tumor's own heterogeneity against it. In fact, stabilizing cancer around certain equilibrium that are not lethal but are sufficient to maintain competition presents an interesting way to treat cancer, whether from a health pointy of view or an economic one. A mathematical model was used to investigate this, as mathematical modeling has always been present in cancer studies ^[20,21].

The paper is organized as follows. In Section 2, we describe the model in detail and subsequently analyze it using optimal control theory, to identify optimal drug scheduling that may control the tumor by stabilizing and targeting CSCs. In Section 3, we present our numerical computations using an indirect method of optimal control theory. Finally, in Section 4, we discuss our results and present some perspectives and eventual extensions of the model.

2.   Model description and mathematical analysis

We herein propose a mathematical model structured into two compartments: one for the CSC population and the other one for the non-stem CCs. We assume that the CC compartment growth is stimulated by the differentiation of CSCs. Indeed, during the cell cycle, CSCs may undergo apoptosis or divide to yield two daughter cells indefinitely (no natural mortality) or differentiate to a CC. However, CCs do not proliferate indefinitely owing to the aging process (a growth rate less than a mortality rate), do not differentiate, and are not subjected to apoptosis.

Suppose that the tumor heterogeneity is due to CSC population heterogeneity. Thus, the CSC population is structured by $n$ types. Through the differentiation process, each CSC of type $i$ yield a CC of the same type. It should also be emphasized that, in some cases, the CC can acquire the ability to de-differentiate and regain the CSC population with all stem cell properties ^[22]. However, we ignore this de-differentiation process and assume that the CSC compartment is independent of the CC compartment (see Figure 1 for the conceptual model overview). In the following, we present the description of the two compartments separately.

2.1.   CSC compartment model

Let $\forall\ i \in \{1, .., n \}$, and $M_i(t)$ denotes the number of resting phase cancer stem cells of type $i$ at time $t$. During the cell cycle, each cell of type $i$ can either proliferate with $\beta_i$ rate or differentiate with $\delta_i$ rate.

Subsequently, the CSC of type $i$ that has entered the proliferative state and escaped apoptosis with a rate ($1- \mu_i$), divide and yield two daughter CSCs. Thus, the total number of daughter CSCs from type $i$ is

Considering that the daughter cells of type $i$ may come from a mother cell of type $j$ with a probability $p(i, j)$, the quantity of cells added to type $i$ is then

The evolution of the different classes of CSCs is therefore given, $\forall \ i \in \{1, .., n \}$, by

with $M = (M_1, \ldots, M_n)$ representing the state vector of all CSC populations.

The CSC compartment model as an optimal control problem. It is typical in cancer treatment to deliver not just one drug, but a combination of therapeutic agents that act on different pathways. We herein propose a combination of drugs, each targeting one specific CSC type. Basically, the sensitivity of cells to treatment is not absolute, but it depends significantly on the drugs administered. Moreover, cell culture experiments suggest that the treatment-induced cell response rate (dose-response relationship) is proportional to the drug concentration ^[23,24]. Subsequently, let $u_i(t)$ be the cell-type specific response rate owing to drug concentration at time $t$.

We therefore formulate an optimal control problem for the system (2.1), where we consider $n$ controls $u = (u_1, \ldots, u_n)^T$ corresponding to $n$ dose rate. Drugs administered are performed under a homologous stem cell treatment (cytapheresis). Thus, $u_i(t)$ might be either positive (i.e., $u_i(t) > 0$, corresponding to the death rate) or negative (i.e., $u_i(t) < 0$, corresponding to the recruitment term).

The evolution of the different classes of CSC population under treatment is therefore given, $\forall \ i \in \{1, \ldots, n \}$, by

where $M = (M_1, \ldots, M_n)$ the state vector of all CSC populations.

The treatment principle consists in maintaining a certain competition between the resistant and sensitive cancer cells in a way that keeps the tumor small enough for the patient to live with it. This principle aims at transforming cancer from a lethal disease to a chronic one. Subsequently, our goal is to obtain the control $\hat{u} = (\hat{u}_1, \ldots, \hat{u}_n)^T$, for each CSC type, to minimize the following functional of the form:

with $T$ the final time of the therapeutic protocol, $\lambda_j$ are positive weights relative to each CSC type, and $M_i^*$ a positive constant. This functional measures the difference between the CSC populations and the constant $M_i^*$ during the therapy interval. Therefore, minimizing $J$ results in the minimization of this distance $\left| M_i(t)-M_i^*\right|$. Thus, the primary objective is to stabilize the volume of the CSC population $\sum_{i = 1}^{n} M_i(t)$ at a certain competitive equilibrium $\sum_{i = 1}^{n} M_i^*$. We assume that the total amount of each drug to be given is specified a priori based on a medical assessment represented by constrained values. In fact, any kind of drug can become toxic for the patient if the dose crosses a certain maximal threshold and is of no benefit to the body if given below a certain minimal dose. This results in the following mathematical condition, $\forall \ i \in \{1, \ldots, n \}$:

Thus, the optimal control problem (OCP) can be stated as follows: minimize the functional $J$ over all Lebesgue measurable functions $u = (u_1, \ldots, u_n)^T$, such that $u_i:[0, T]\rightarrow [u_i^{min}, u_i^{max}]$, $i \in \{1, \ldots, n \}$, subject to the dynamic (2.2).

Existence of an optimal control. We examine sufficient conditions for the existence of a solution to the quadratic optimal control problem. The necessary conditions for existence are stated and verified as follows ^[25].

Theorem 2.1. There exists control functions $\hat{u} = (\hat{u}_i), \; i \in \{1, ..., n\}$ so that

where $U$ is the set of admissible controls.

Proof. To prove the existence of an optimal control pair, it is easy to verify that:

● The set of controls and corresponding state variables, $V = \{ (u, M_u(t)) \}$, where $u = (u_i)$ is the control vector and $M_u = (M_{i \, u})$ is the corresponding state variables vector, is nonempty, convex and closed.

● The right-hand side of the state system is a linear function of the state variables as well as the control.

● The integral $L = \sum\limits_{j = 1}^{n} \lambda_j \left(M_j(t) - M_j^* \right)^2$ in the objective functional is positive quantity in Eq (2.3) is convex. This condition holds, because the integral $L$ is quadratic.

● There exist constants $\omega_1 > 0$, $\omega_2 > 0$ and $\rho > 1$ such that the integral $L$ of the objective functional satisfies:

In fact, since $L > 0$, there exists $l > 0$ such that $L > l$, then, for $\omega_1 = \frac{l}{2\, n\, (u^{max})^2}$ where $u^{max} = \underset{1 \leq i \leq n}{max} \, u^{max}_i$ and $\omega_2 = \frac{l}{2}$ one has the previous equation.

These conditions provide the existence of an optimally controlled trajectory $(\widehat{M_i}, \widehat{u_i})$ for the problem (OCP) defined over the interval $[0, T]$, that minimizes Eq (2.3).

Characterization of the optimal control. In order to characterize the optimal control, we use the Pontryagin maximum principle ^[26,27]. This principle allows transforming the optimality problem into a mere study of the Hamiltonian variations. So, we start by giving the Hamiltonian defined as:

where $\psi = (\psi_1, ..., \psi_n)$ is the vector of co-state (adjoint) variables and $F(t, M(t), u(t)) = (F_1(t, M(t), u(t)), ..., F_n(t, M(t), u(t)))$ is the vector of state variables.

Theorem 2.2. Given optimal controls $\hat{u}_i$ and solutions $M_i$ of the corresponding state system, there exist co-state variables $\psi_i$ that satisfy

with terminal condition:

Furthermore, the optimal controls $\hat{u} _i$ are given by

where

Proof. Using the Pontryagin maximum principle conditions, we obtain the co-state variables equations given by

with terminal condition:

Furthermore, the optimality conditions state that the optimal control $\widehat{u}(.)$ maximizes the Hamiltonian function H. By using them, one has the switching function:

As $M_i$ is assumed positive, the sign of the switching function $\gamma_i$ is determined by the sign of the co-state vector $\psi_i$. According to the previous condition, the Hamiltonian function $H$ is linear with regard to $u_i$, thus, the maximization property implies that $\forall \ i \in \{1, \ldots, n \}$:

Therefore, the maximum rate $u_i^{max}$ of the drug should be administered whenever $\psi_i$ is negative, and the minimum rate $u_i^{min}$ of the drug should be injected whenever $\psi_i$ is positive. If we do not consider the singular solution, these types of controls are called bang-bang controls, while they assume values only in the boundary points of the control interval.

Now, we investigate the singular solution. This solution is considered if, over some interval $I_i \subset [0, \, T ]$, one has $\Gamma_i(t) = 0$ for all $t \in I_i$. Note that over this interval, it is assumed that $M_i > 0$ and $\psi_i = 0$ since $M_i = 0$ means that the patient is perfectly healthy and does not require any treatment. To compute the singular control, one needs to differentiate the switching function $\Gamma_i(t)$, it is well-known that the control can appear for the first time only in an even numbered derivative, and thus we need to compute the first two derivatives.

This expression obtained by replacing $\dot{M}, \ddot{M}, \dot{\psi}$ and $\ddot{\psi}$ by their dynamics expressions can be written as follows:

where

and

Since $\ddot{\Gamma}_i = 0$, the singular control is then given by:

under the following admissibility hypothesis:

In the following, we elaborate the CC compartment model and analyze the asymptotic behavior of the entire tumor population when the CSC equilibrium $M_i^*$ is reached for each $i \in \{1, \ldots, n \}$.

2.2.   CC compartment model

Let $\forall \ i \in \{1, \ldots, n \}$, $C_i$ denote the density of the CC of type $i$. Each type $i$ of CC is subject to natural mortality with a rate $\mu_{Li}$, and to proliferation with a rate $\beta_{Li}$. Moreover, the growth of each type $i$ of CC is stimulated by the differentiation of the type $i$ of CSC with rate $\delta_i$. Considering the aging process, we assume that for each type of CC, we have $2\beta_{Li}-\mu_{Li} < 0$. In fact, unlike CSCs which have very low apoptosis rate, fully differentiated CCs are more likely subject to natural mortality than to proliferation and the renewing of the population depends on CSCs. Therefore, the CC equations are as follows:

The solution of the Eq (2.8) is given by:

If $M_i$ is equal to $M_i^*$, subsequently the Eq (2.9) resumes to

Therefore, considering the aging process, i.e., $(2\beta_{Li}-\mu_{Li}) < 0$, we have

This proves that the optimal control of the CSC population results in the stabilization of the entire tumor.

In the following section, we describe the numerical computations with two type of CSCs. Our aim is to identify the switching times for the bang bang solutions with different scenarios of cancer aggressiveness.

3.   Numerical computations

For the numerical simulations, we consider our model as qualitative, rather than quantitative, such that the parameter values cannot be claimed as being entirely realistic. They are not based on medical data, but merely reflect some qualitative biological tumor properties.

We assume that the cells are structured by two types of CSCs: Type $1$ (type $2$) exhibits a high (low) proliferation rate with a lower (higher) differentiation and apoptosis rate. The treatment administered was performed under a homologous stem cell treatment (i.e cytapheresis). Time was normalized to set the boundary points of the control interval to $u_i^{max} = 1$ and $u_i^{min} = -1$. As explained above, we seek to determine a control that can stabilize CSC populations at $M_i^*$ equilibrium; therefore, we chose $M_1^* = 2500$ and $M_2^* = 1500$. Moreover, we considered the aging process for different types of CCs; therefore, $\forall \ i \in \{1, \ldots, n \}$, $2\beta_{Li}-\mu_{Li} < 0$.

We numerically solved the following Hamiltonian system using the single shooting method ^[27], as the problem is a two-point boundary value problem for which both the initial states of the state variables and the final states of the co-state variables are known. The control was chosen in accordance with the sign of the corresponding switching function, $\Gamma = - \psi_i M_i$ (i.e., we did not consider the MRAP solution):

Figures 2(a), (b) illustrate the therapeutic protocol adapted to both type $1$ and type $2$ CSCs. The protocol adapted starts with a constant maximal dose over a $1.2$ unit of time period for type 1 and $2.8$ for type 2; which results in both cases in a decreasing of the population density. This shot is followed by a recess in both cases that causes the populations to regrow. The cycle is then maintained until reaching the equilibrium values. The generated protocols succeeded in rebalancing these cell populations in the neighborhood of $M_1^*$ and $M_2^*$ (i.e., prescribing a lifetime treatment for the patient that will control the cancer without completely eliminating it). The principle is to maintain the competitive environment between resistant and sensitive cancer cells. The optimal treatments obtained by the simulations postulate that courses should be increasingly more dense (i.e., chemotherapy treatments are closer together) and progressively less intense (i.e., intake duration).

However, the equilibrium $M_2^*$ for the type-$2$ CSC is reached faster compared to the other type (see Figure 2(b)). Indeed, we observed that the dynamic of type-$2$ CSC indicates less important oscillations than that of type $1$. This can be explained by the fact that type-$1$ CSCs are more aggressive than type-$2$ CSCs, considering our parameter configuration. Moreover, the ratio $\frac{\rm{Cytotoxic\; drug\; amount}}{\rm{cytapheresis\; drug\; amount}}$ for type-$2$ CSC is larger than that for type $1$. This is consistent with the hypothesis, in that with not extremely aggressive cancers, relativity large doses of cytotoxic treatment can be administered without an extended rest period.

Monotherapy against multidrug therapy. Monotherapy refers to medical treatment that relies on only one drug administration (i.e., one agent) at a time, unlike a multidrug therapy that is a treatment comprising several agents simultaneously. Although the response rates and survivals observed in multidrug therapy programs appear to be higher than those observed with monotherapy, few randomized trials have been conducted for this purpose ^[28,29].

Figure 3 shows that the two optimally generated scheduling of the two drugs overlap several times (i.e., the two drugs should be injected simultaneously). Indeed, at 45% of the treatment period, we have to administer both types of drugs simultaneously. However, for 52% of this treatment period, we alternate between these two drugs. This suggests that combining the two strategies: multidrug therapies and monotherapies, is essential for heterogeneous cancers. Moreover, the efficiency of these two strategies depends on their scheduling.

Influence of CSC population treatment on the entire tumor. We herein propose to study the effect of targeting CSCs and CCs on the whole tumor. We study three scenarios depending on the sign of $\mu_{Li} -(\delta_i+2\beta_{Li})$, for $i = 1, 2$, thus ensuring that at the beginning (i.e., initial conditions), the density of the CC is much higher than that of the CSC. Figure 4 represents the densities of CSCs and their corresponding CCs in three different cases. In all cases (see Figure 4), we observed that targeting the CSC using the optimal drug scheduling can control and stabilize all the tumors.

However, we observed that we can even reduce the size of the CC populations to less than those of the CSC in both cases where $\mu_{Li} > \delta_i+2\beta_{Li}$, for all $i = 1, 2$, and $\mu_{Li} > \delta_i+2\beta_{Li}$, for $i = 1$ or $i = 2$ (see figure 4c). In the latter case, we can reduce the size of only the non-aggressive CC population, allowing the size of the aggressive CC population to be relatively high and subsequently rendering the tumor homogenous. Nonetheless, if $\mu_{Li} < \delta_i+2\beta_{Li}$, for all $i = 1, 2$, we cannot reduce the size of CC populations to less than those of the CSC. In this case, it would be advisable to add a cytotoxic treatment to kill the remaining CCs (see figure 4b). Moreover, we observed in Figure 4 a significant delay between the CSC and CC. This is because the stabilization of CSC induces CC stabilization.

4.   Discussion

In this study, we created a controlled mathematical model for the treatment of CSCs and CCs where drug doses are the control. Optimal control allows determining the optimal dose and schedule for the delivery of treatment. The mathematical model includes an $n$ ordinary differential equations of first order for both the CSCs and the CCs.The model describes the dynamics of CSCs and CCs when applying the optimal combination. Other than describing the dynamics, the model yields the best schedule according to the parameters used. The mathematical analysis of the problem was followed by numerical simulations using the single shooting method. The optimal drug scheduling generated by our optimal control problem suggests that drugs should be administered continuously over relatively long periods without extended rest periods and in a chronic manner. This optimal scheduling joins strongly the new modality of drug administration called "metronomic chemotherapy" that creates reference to the schedule ^[14,16]. It is a new therapeutic approach that began to evolve in the early 2000s, when Fidler and Ellis reported, "Cancer is a chronic disease and should be treated like other chronic disease ^[15]". This approach appeared after a reflection regarding the side effects caused by traditional therapies. Because targeting a CSC is not a typical practice, we subsequently compared the optimal drug scheduling generated by our model with traditional ones.

Traditional therapies (i.e., chemotherapy) have been designed to kill as many tumor cells as possible by treating them with the "maximum tolerated doses" (MTDs) of cytotoxic agents. Many of these therapies are effective and result in regression, or at least stabilization, of the tumor. However, these responses are typically short-lived and relapses are often marked by more aggressive and resistant cancers. In addition, these traditional therapies seriously affect the quality of life of the patient. Moreover, the MTDs has been proven not to be optimal ^[30]

It is noteworthy that metronomic dosing has been performed in a clinical setting. Indeed, dozens of chemotherapeutic agents have been administered by metronomic infusion ^[17]. Nevertheless, it remains a rare practice owing to the lack of a standard planning of metronomic infusions. Hence, our study attempts to partially solve this problem, by focusing on the switching times between cytotoxic treatment and recovery time.

Another result of our model suggests that combining multidrug therapies and monotherapies is more efficient for heterogeneous cancers than using only either one. It contributes to the therapeutic controversy that aims to determine the best manner to administer drugs. In this context, the Goldie–Coldman hypothesis predicted that alternating two drugs (i.e., monotherapy) would be the most effective strategy. This hypothesis was tested negatively in Italy in women with mammary carcinoma ^[31]. This negative result is also supported by the experimental study performed by Alberto et al. ^[32]. Alberto et al. demonstrated an unquestionably higher objective response rate when four drugs were injected simultaneously, followed by the same four agents injected separately.

Moreover, our simulations confirm the fact that by controlling and targeting only the CSC, we manage to control the entire tumor, thus preventing possible tumor recurrences. We also noticed a certain delay between the CSC and CC; the CSC reached equilibrium faster than the CC, thus supporting the fact that CCs are dependent on CSCs. Therefore, an optimal protocol should target cells that are critical in the tumor ecosystem instead of the most abundant ones.

Few studies have been conducted on the development of molecules specifically targeting CSCs. However, it has been shown that the signaling pathways: Wnt, Notch, and Shh, are involved in the regulation of CSC differentiation and proliferation, and are necessary for the survival of these cells ^[33]. The Wnt signaling pathway has been shown to be selectively expressed in colon CSCs and promotes their maintenance and proliferation ^[34]. In addition, the over-expression of the Wnt pathway would promote the tolerance of DNA damages in mammary CSCs, thus resulting in their survival ^[1] from a transformation of stem cells into neoplastic cells ^[35]. In addition, the inhibition of Shh and Notch pathways in glioblastoma cell lines results in the increased sensitivity of CSCs to temozolimide ^[36].

Finally, the exact origin of these CSCs has not been clearly established. Genetic alterations could occur during different stages of cell maturation, resulting in the de-differentiation of different precursors to CSCs ^[2]. For example, it has been shown in mice that neural stem cells or mutated neural progenitors may be responsible for gliomas formation ^[37]. This will render the task of balancing the stem cell population even more difficult owing to the new term of recruitment to this population. Several hypotheses suggest that this term of recruitment could be induced by the treatment that explains the effect of resistance and tumor recurrence. Thus, it would be wise to extend our model, considering the assumption that in some cases the normal cancer cells (i.e., the non-stem cancer cells or poorly differentiated cells) acquire the ability to de-differentiate and regain its stemness. We will thus be able to study the resistance to treatment and the phenomenon of recurrence.

Acknowledgments

GH would like to thank PDI-MSC of IRD and SU for the PhD fellowship.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

Parameters