Abstract
Therapeutic protocols in immunotherapy are usually proposed following the intuition and experience of the therapist. In order to deduce such protocols mathematical modeling, optimal control and simulations are used instead of the therapist’s experience. Clinical efficacy of dendritic cell (DC) vaccines to cancer treatment is still unclear, since dendritic cells face several obstacles in the host environment, such as immunosuppression and poor transference to the lymph nodes reducing the vaccine effect. In view of that, we have created a mathematical murine model to measure the effects of dendritic cell injections admitting such obstacles. In addition, the model considers a therapy given by bolus injections of small duration as opposed to a continual dose. Doses timing defines the therapeutic protocols, which in turn are improved to minimize the tumor mass by an optimal control algorithm. We intend to supplement therapist’s experience and intuition in the protocol’s implementation. Experimental results made on mice infected with melanoma with and without therapy agree with the model. It is shown that the dendritic cells’ percentage that manages to reach the lymph nodes has a crucial impact on the therapy outcome. This suggests that efforts in finding better methods to deliver DC vaccines should be pursued.
1. Introduction
Immunotherapy based on antigen-treated dendritic cells (DCs) is a promising treatment against certain types of cancer [1]. This kind of therapy is often regarded as a safe option alone or in concurrence with other therapies [2]. Nonetheless, clinical evidence of its efficacy is still unclear and researchers are still designing new experiments in mice and humans to decipher mechanisms that could lead to a successful cancer treatment. Murine models are a current way to test hypothetical therapy protocols by cultivating dendritic cells in vitro from mouse bone marrow to be inserted back into the mouse after being treated with antigens. For example, modest clinical success has been observed in C57BL/6S mice inoculated with B16 melanoma cells [3] after being treated with antigen loaded DCs.
Therapy schedules consist in previously planned protocols of injection times with their respective vaccine quantity. Such protocols in immunotherapy are usually proposed following therapist traits such as intuition and experience. We aim to provide a rational therapy planning not only relying on therapist’s experience and intuition but also by the help of mathematical modeling, optimal control, and simulations. Optimal control has a history as a tool to improve tumor therapy schedules. It is usually used to rationalize issues like infusion times and drug quantities. For instance, the optimal control for a model with mixed immunotherapy and chemotherapy is addressed in [4]. Numerical and analytic control techniques for continuous and bang-bang controls can be consulted in [5]. Therapeutic protocol improvement by simulations have been made in [6] for Cytotoxic T Cell (CTL) and in [7] for dendritic cell transfections.
Maturation of antigen-primed DCs reduces phagocytic capacities in an interchange for improving and presenting migration capabilities [2] to lymph tissues increasing expression of chemokines and adhesion molecules. DCs are characterized by improving the immune response activating natural killer cells and naive and memory B cells. Moreover, DCs create cytokines which help T cells differentiation. Several methods to generate and mature DCs are able to produce antigen-specific T cell response after DCs inoculation. Maturation is essential since mature DCs have a greater capacity to create cytokines and chemokines that activate T cells in contrast to immature DCs.
Tumor cell antigens need to be loaded into Major Histocompatibility Complex (MHC) molecules of the DCs in a process called antigen loading [2]. For such an endeavor DCs are curated with peptides and proteins from tumor cells. Another method to load antigens to DCs is by using virus vectors. Such viruses are used to load genes with encoded Tumor Associated Antigens (TAA) into the dendritic cells. Also, the natural immune response against virus vectors serves as an immunostimulant towards the TAA.
DCs are regarded as the best at activating cytotoxic T cells. They have a history of being used as the therapy to improve the host immune capabilities. Nonetheless, many clinical results of dendritic immunotherapy are not quite successful. The purpose of this work is to propose a model that accounts for the obstacles that dendritic cells could face when they are used as a medium to transmit antigens to the host immune system and create immune response. Such obstacles are as follows:(i)Only 5% of the DC vaccine reaches the lymph node [8].(ii)Injected DCs do not directly activate cells; they only transmit the antigens to the endogenous Antigen Presenting Cells (APCs) which are responsible for activating the T cells [9]. We suppose that DCs only interact with T cells and IL-2 cytokine.
We also include the immunosuppressive effects of the TGF-β in the cytotoxic T cells. Then, the model is used together with optimal control to overcome such obstacles. In this case, optimal control is used to predict satisfactory dendritic injection times with an optimal control technique developed by Castiglione and Piccoli [10]. Another consideration is regarding the kind of DCs vaccine. We do not consider a continuous infusion of DCs, but a short period or bolus injection that can be considered instantaneous. This resembles more accurately the procedure following the immunotherapy injection in [3] for which we adjusted our model.
The outline of this work is as follows. The explanation of the mathematical model and its terms are given in Section 2, followed by numerical simulations that agree with experimental data from mice bearing melanoma cells in Section 3. Section 4 contains the optimal control problem and the optimization algorithm, followed by Section 5 that explores the consequences of changing the percentage of dendritic cells arriving at the lymph node. Finally, Section 6 gives some biological implications of the results.
2. The Murine Mathematical Model
The clinical efficacy of DC vaccines by simulations of mathematical models is a topic of ongoing research in mathematical oncology [1]. One of the most used tools for modeling the intricate interactions between proteins, lymphocytes, and tumor cells is ordinary differential equations (ODEs). Other kinds of models include delay differential equations (DDEs), which are a system of ODEs where an unknown depends on past values. This is usually represented in the independent variable as (e.g., and is changed for ) where is called the delay. Such delay can elegantly represent a gestation time or a transport delay, hiding some of the intricacies of the phenomena. The negative side is that it is often more difficult to solve a DDE analytically and numerically than an ODE. Several DDE models have been used to hide the transport delay of DCs from the injection time to the interaction time with effector and helper T cells. For example, a compartment model between the spleen, blood, and tumor with a delay between the DCs injection and spleen arriving is given in [11]. A DDE with interactions between TGF-β and T cells is made in [7].
For the model considers an inherent delay h, which is equal to the time that DCs take to arrive at lymph nodes but only affects the times at which the simulations assume an injection and is not implicitly adding . Then, our model is an ODE. Also, we are considering that therapy consists of bolus injections with very small timespan in comparison with overall therapy. The details of implementation are explained in The Simulation.
2.1. The ODE System
The model of the present paper having been inspired by the works of [6, 7, 10] is given by the ODE system:which has state variables: (i), the tumor cells.(ii), the T helper cells.(iii), the T or cytotoxic cells.(iv), the antigen loaded dendritic cells.(v), the IL-2 Interleukin-2 cytokine.(vi), the TFG-β T cell inhibitor.(vii), the IFN-γ which upregulates MHC class 1.(viii), is the number of MHC class 1 receptors per melanoma cell.
The first term of (1) is a Gompertz growth which adjusted very well to experimental data; see Figures 1 and 2. The expression in the second term represents the tumor cell eliminations by the CD8 T cells, which in turn is suppressed by the cytokines TFG-β and the efficiency of MHC class 1 receptors. Both suppressing effects follow Michaelis-Menten saturation dynamics; see [6].
Equations (2) and (3) represent the dynamics of CD4 and CD8 T cells, respectively. The parameters and represent the production of CD4 and CD8 T cells. Both follow a logistic growth law. As mentioned before, it is supposed that the injections of DCs only interact with CD4 T cells since it was reported in [9] that CD8 activation is not made by DC vaccines. So, in the term we consider that activation of tumor-specific CD8 is only promoted by IL-2.
The term in (4) is a control variable representing the DCs injections. Since we are considering a bolus injection therapy, is the sum of impulses at injection times: . Where the therapy consists of doses of size given at times , but because there is a delay in the DCs arriving time at lymph nodes the impulse occurs at .
Equation (5) have interleukin IL-2 produced from the interactions of DCs and CD4 T cells by and is consumed by CD8 T cells at . Interleukin IL-2 has a major role in the system, prolonging persistence of CD8 T cells. It is necessary for the production of new tumor antigen-specific CD8 T.
Equation (6) describes the cytokine dynamics of TGF-β, an inhibitor of T cells activity, with secretion proportional to the number of tumor cells and degradation rate .
Equation (7) describes the cytokine IFN-γ, a weapon of the CD8 T cells to upregulate the MHC class 1 in the melanoma cells. The production is proportional to the CD8 T cells and degradation rate .
MHC 1 class dynamics is in (8). It has a basal production rate per melanoma cell . The first term shows a growing rate stimulated by IFN-γ which follows a Michaelis-Menten Kinetics with maximal effect and Michaelis parameter .
This model assumes the following:(1)The tumor mass is homogeneous. The present model assumes that each tumor cell has the same sensitivity to DCs vaccine.(2)It is monoclonal, that is, only T cells are able to recognize the TAA.
3. The Simulation
The model is calibrated to fit the data supplied in [7] made at the Laboratory of Immunotherapy and Tissue Engineering of UNAM (Universidad Nacional Autonoma de México), Mexico [3]. The immunotherapy consisted of three doses of DCs activated with MAGE-AX separated by hours, in Table 1. The tumor cells population is calculated from mice tumor diameter assuming a spherical form. The used parameters can be seen in Table 2. The initial condition is . A comparison between the experiment and simulations with and without therapy is shown in Figures 1 and 2.
As mentioned before, in the simulation there is an inherent delay of h, which was adjusted to fit the experimental data. Hence, we implement to affect the impulse time of the control variable . Therefore, for example, the therapy consists of a single injection of DCs at h; we havewhere is the Dirac delta function. Observe that this implies that the DC injection at h has an effect on the system until . So, we integrate our system from , and then restart the integration from with initial condition .
4. Optimal Control
Hypothetical protocols can be found simulating the model with several combinations of dose timing and size. This could lead to the implementation of new protocol therapies which could potentially replace the originally proposed by medics. For example, see [7]. Now in this section, we show that such hypothetical protocols can be derived or improved using an optimization algorithm based on optimal control.
It is considered that a therapeutic protocol consists of doses of size given at times . Now, let a schedule of injections be given by , where is the fixed time horizon. Let be the space of schedules, then for a particular schedule the control variable takes the form: where is the percentage of DC vaccine that actually reaches the lymph nodes. Notice, we are approximating the doses as an impulse given by the Dirac function . Now let , so Equations (1)–(8) take the form,
The optimal control problem consists in the following:(P)Determine the schedule of injections that minimize the final tumor mass of the trajectory given by (11) with initial condition .
To solve (P) we use the optimization algorithm described in [10], which is repeated here for clarity, with a slight modification.
4.1. Optimization Algorithm
Algorithm 1 (optimization). (S0)Fix the time horizon , the number of vaccine administrations, the value of vaccine quantity, an initial value of cells population, and an initial schedule .(S1)Integrate system (1)–(8) with initial value to obtain the trajectory . Solve for each of the schedule : in the interval with initial condition where is the fourth coordinate vector.(S2)Update the schedule with . The step size should be of the order to avoid an excessive big step. Go to step 1.
The correct choice of step size is important. If is too big, we could get a wrong optimization step. If is too small the optimization could improve very slowly. Final tumor mass was decreased in each optimization step choosing of the order . Also, to find the best an search method could be used such as the Golden Section Search (GSS) [12]. Using GSS we achieve a good optimization with steps in about seconds. In contrast to a fixed step it needs optimization steps in 60 seconds.
We tested Algorithm 1 with an initial random schedule and 100 optimization steps with . Also, in order to get at step 2 we used the same initial condition: as that used to test the model (Figure 2). This evolved into one big injection of at h and one of at h. Observe that a total dose of DCs is used, whereas the original protocol used DCs (Figure 2).
The difference of initial and optimized therapy is shown in Figure 5. Moreover, Figure 5(b) shows a maximum of roughly tumor cells and final tumor mass of around ; substantially less than in the experiment. The maximum and final tumor mass is and respectively, see Figure 2.
Figure 4 shows that a fast final tumor mass decay occurs in the first optimization steps with little improvement afterwards. Nonetheless, Figure 3 shows a substantial therapy evolution until about the 80th optimization step.
(a) The starting therapy
(b) Therapy after the optimization
4.2. Optimal Therapies in Large Time Horizons
A very slow convergence and unsatisfactory results of Algorithm 1 are observed searching for protocols far beyond h (e.g h). In view of such problems, we divide the time horizon in equal intervals of size and apply Algorithm 1 at every interval for in an iterative way. At the end of each iteration step the state is saved to be used as an initial condition in the following iteration. Also, in each iteration the optimized schedule is saved. Finally, the system is solved with the schedule .
Algorithm 2 (optimization of the intervals). (S0)Fix the time horizon , number of intervals , the number of vaccine administrations on each interval, the value of vaccine quantity, an initial value of cells population and an initial schedule . Let and .(S1)If equals , go to Step 3. Otherwise apply Algorithm 1 at interval with initial condition . Get the schedule .(S2)Let . Save the next initial condition , Let . Go to Step 1.(S3)Solve the system with the schedule .
5. Results and Discussion
The objective of this section is to show that the percentage of DC injection which arrives at lymph nodes has a great impact on the tumor mass behavior. Three therapy cases are shown, each starting with a random therapy which is improved by Algorithm 2.
For ethical reasons, melanoma inoculated mice are allowed to live maximum of 1000 h. Then, in order to find a therapy after such time we apply Algorithm 2 for timespans: 20,000, 10,000, and 5,000 hours; in addition, we change from to and to . In every case tumor oscillations with amplitude that stabilize after some time are obtained. This stabilization time seems to be smaller as grows. It is worth mentioning that , the maximum number of melanoma cells that mice can bare before dying cannot be greater than [7]. So, the immunotherapy protocol with shown in Figure 6 can be regarded as a successful therapy to control the tumor under .
When is increased, Figures 7 and 8 show a therapy with better tumor-killing performance than therapy with . Using (Figure 7) reduces the maximal oscillation amplitude to and reduces the amplitude of the stable oscillations to . Finally, in Figure 8 reduces the amplitude of the stable oscillations to the order which in some instances is regarded as clinically not detectable.
The increments of can be implemented in experiments using a secondary therapy that improves the migration capacity of the DCs. Also, alternative delivering methods for DCs injections could be used. Nevertheless, alternative delivering routes such as intradermally and intranodally showed a comparable response to the intravenously route [8].
5.1. Further Research
The overall schema followed in Algorithm 1 is called the steepest descent. This relies on following the correct gradient in each step until a minimum of the cost function is apparent. The gradient for our problem is given by (12) and only considers the cost function [10]. Despite this, the simulations in Figure 5 show improvements over the originally used therapy. However, there are many cost functions which can be minimized. For example, a cost function that admits the total dose given over the optimization can be proposed. This could help therapists to economize on materials which cost a lot of money and effort in cultivating and transferring antigens to the DCs. Therefore it is observed that optimization in Figure 5 results in a lesser total dose. For further cost functions associated with gradients on immunotherapy optimization we refer to [13].
Melanoma is well known for its immunogenic properties usually used for immunotherapy testing. For example, one of the most cited studies on cancer immunotherapy [14] consisted in a vaccine trial of 440 patients where of the patients had melanoma, although only of the patients responded to the vaccine and the majority was from melanoma. Corresponding patients can be found suffering from lung, childhood, and kidney cancer. This could motivate the creation of mathematical models for therapy optimization on such types of cancer.
It can be observed in Figures 6, 7, and 8, a periodically pulsed immunotherapy can be used to control the tumor. This results in what seems to be periodical tumor solutions. Are these periodical solutions stable? This can be investigated creating a map where is the solution of (1)–(8) at with initial condition . Then, the existence of periodical solutions is equivalent to solve the equation
The stability is given by the Jacobian of the map where is a periodic solution. These stable periodic solutions can be regarded as more robust than the unstable ones.
6. Conclusions
A model has been created that includes several of the obstacles that DCs face in the host environment such as immunosuppression and poor transference to the lymph nodes. The model shows an agreement with experimental data from mice inoculated with melanoma. This gives us more confidence about the optimizations outcome.
We used optimal control as a tool to rationalize the creation of immunotherapy protocols. These protocols can serve the therapist to complement their intuition and experience. Simulations show that optimized protocols outperform that proposed by the therapists in [7] regarding the total dose and final tumor size.
Although the real clinical efficacy of DC therapy is still under discussion, the results of the present work should be tested in experiments. The simulations show that an increment in the percentage of DCs that manage to arrive to the lymph nodes () has a huge impact on the amplitude of the oscillations made by the periodic therapy. Even an increment from to reduces the amplitude from the order of to the order of as is shown in Figure 8. This could be achieved using better methods to deliver DC vaccines or combining DC immunotherapy with treatments that enhance the DCs migration capacity.
Can immunotherapy be a cure for cancer? That is too early to be answered. What is sure is that most of the current therapies rely on prolonging the survival of the patient and not the cure. The simulations shown in Figures 6, 7, and 8 support the idea of DCs immunotherapy as a medium to tumor control.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors wish to acknowledge “Consejo Nacional de Ciencia y Tecnología” (CONACyT), “Comisión de Operación y Fomento de Actividades Académicas del Instituto Politécnico Nacional” (COFAA-IPN), and “Estímulos al Desempeño de los Investigadores del Instituto Politécnico Nacional” (EDI-IPN) for the support given for this work.