Abstract
We present a mathematical model developed to reproduce the immune response entitled with the combined administration of activated OT1 cytotoxic T lymphocytes (CTLs) and Anti-CD137 monoclonal antibodies. The treatment is directed against melanoma in B16 OVA mouse models exposed to a specific immunotherapy strategy. We model two compartments: the injection point compartment where the treatment is administered and the skin compartment where melanoma tumor cells proliferate. To model the migration of OT1 CTLs and antibodies from the injection to the skin compartment, we use delay differential equations (DDEs). The outcomes of the mathematical model are in good agreement with the in vivo results. Moreover, sensitivity analysis of the mathematical model underlines the key role of OT1 CTLs and suggests that a possible reduction of the number of injected antibodies should not affect substantially the treatment efficacy.
1. Introduction
Melanoma is a malignant tumor caused by the mutation of melanocytes, that is, the cells that produce the melanin and are responsible of the color of the skin. Despite intensive research, melanoma still represents one of the most aggressive malignant cancers [1]. Many experimental approaches are now focused on targeting cytotoxic T lymphocytes (CTLs) against cancer. A common strategy to enable CTL efficacy against tumor is to activate naïve CTLs in vitro through the use of cells engineered to present the tumor antigen, and to reinject them in the host. However, even if activated CTLs are able to infiltrate into tumor masses, in most cases they remain unable to contrast cancer growth [2]. As experimental evidence suggests, tumor-infiltrating lymphocytes are rendered ineffective by coinhibitory molecules expressed by tumor and stroma cells surfaces [3].
In order to gain complete rejection of tumors, injection and stimulation of CTLs is not sufficient and should be, therefore, coupled with complementary measures voted at boosting CTLs migration inside tumor masses, and conjugation and killing of target cells [4–6]. One way of boosting CTLs actions is represented by stimulation through the binding of costimulatory proteins expressed on CTLs surface. Among possible surface proteins, Anti-CD137, also known as 4-1BB, represents a valuable target. This protein is expressed by multiple IS cells such as activated T, NK, B-lymphocytes, dendritic cells and also by tumor endothelium cells [7]. Its natural ligand (CD137L) can be found on activated antigen-presenting cells surface [8].
The combined administration of monoclonal antibodies specifically targeted to bind Anti-CD137 proteins and in vitro activated-OT1 CTLs was demonstrated to be able to prevent the melanoma formation in B16-OVA mouse models [7]. Moreover, the combined treatment avoided appearing of undesired side effects like the hepatotoxicity, observed only under anti-CD137 only high-dosage treatment [9]. The IS stimulation mechanisms of Anti-CD137 immunostimulatory monoclonal antibodies are multilayered and include the improving of cytotoxicity, duplication rates, and chemotaxis sensitivity of activated-OT1 CTLs [6, 10–12].
To reproduce the dynamics of this biological process, a delay differential-equation-(DDE-) based model has been developed. The model reproduces two different compartments: the injection point compartment, where both antibodies and OT1 cells are injected and the skin compartment where melanoma develops.
2. Biological Background
The in vivo experiment is carried on B16-OVA mice, mice transduced with the chicken ovalbumin gene. The ovalbumin is used as a model tumor antigen. B16 melanoma cell line was derived from an aggressive spontaneous melanoma in pure C57BL6, and B16F10 was derived as a clonal variant from a lung metastasis of this cell line. In tumor immunology, these variants of melanoma are considered poorly immunogenic in the sense that immune-mediated rejections or growth retardations are difficult to achieve.
The experimental setup is oriented to model therapeutic synergy between anti-CD137 monoclonal antibodies and adoptive T cell therapy in melanoma. B16-OVA is a poorly immunogenic murine tumor. The treatment protocol includes a single injection of anti-CD137 mAb and adoptive T cell transfer of OVA-specific TCR-transgenic CD8 CTLs.
In vivo experiments have been executed by Professor Melero and coworkers at the University of Navarra [13]. Mice are divided in five different groups; all groups are composed by five individuals. Each group is treated with a different treatment: Untreated (control) mice, mice treated with naïve OT1 CTLs, mice treated with naïve OT1 CTLs and Anti-CD137 monoclonal antibodies, and mice treated with in vitro activated OT1 CTLs, mice treated with Anti-CD137 monoclonal antibodies, mice treated with in vitro activated OT1 CTLs and Anti-CD137 monoclonal antibodies. The experiment runs for 30 days. At day 0, all B16-OVA mice receive one injection of melanoma malignant cells. The therapeutic treatment used during in vivo experiments is composed by one single boost, and it is administered at day 3. Melanoma surface measurements () are taken at given times for each treatment and are used to estimate the efficacy of each vaccination strategy. We note here that in order to compare in vivo and in silico results we computed the estimated mean surfaces entitled with the use of each treatment. Among the tested treatments, only the combined administration of activated OT1 CTLs and antibodies was able to show complete depletion of the tumor burden, whereas the other treatments remained almost ineffective [13].
3. The Model
We realized a model with two compartments in order to reproduce the dynamics of the process. The first compartment is represented by the injection point compartment where the treatment is administered, whereas the second one is represented by the skin compartment where melanoma tumor cells proliferate and where the cancer-IS competition occurs. To this end, a system of seven delay differential equations has been set up. The model takes into account the following entities: injected activated OT1 CTLs () and injected antibodies () for the injection point compartment; melanoma cells (), tumor antigens (), activated OT1 CTLs and antibodies that have reached the skin ( and ) and naïve CTLs () for the skin compartment. It has also been assumed that Injected OT1 CTLs, and antibodies move from the injection point to the skin compartment only. Figure 1 shows the conceptual model for the biological problem; model entities are listed in Table 1.
Equation (1) refers to first compartment and represents the time evolution of injected activated OT1 CTLs. In (1), represents a known function that models the number of inoculated entities at the scheduled injection time . cells migrate from the injection point compartment to the skin compartment with given rates () and are subject to natural death ().Antibodies (injection point compartment):
Similary to (1), (2) refers to the first compartment and represents the time evolution of antibodies (). are injected at given times and at given quantities , according to the function , and can migrate and disappear from the system by natural degradation at given rates ( and , resp.).Activated OT1 CTLs (skin compartment):
Activated OT1 CTLs that reach the skin compartment () are modeled by (3). Antibodies in the skin compartment () have multiple positive effects on activated OT1 CTLs dynamics (). One of these effects is represented by the ability of promoting duplication. This is modeled through the Holling type II function , where is the maximum biological duplication rate of and is a tuned threshold. When the number of is high enough, the term tends towards 1, thus entitling maximum duplication rates for . The term is used to model migration of OT1 CTLs from the injection compartment to the skin compartment. We suppose here that migration from the first compartment () to the second one () entitles a time delay of and occurs with a given rate .
The term models the activation of naïve CTLs thanks to the presence of antigens () released by killed cancer cells. The biological process that explains the presence of this term is summarized as follows. Antigenic sequences released by killed melanoma cells may be captured by antigen presenting cells (APC) such as macrophages and dendritic cells. These cells process the antigens and present them on their cellular surface to naïve CTLs (). After this presentation process, naïve CTLs cells may be activated, and if some complementary biological steps are accomplished (i.e., stimulation by cytokines released by T helper cells), they can become able to kill tumor cells. This process is not modeled in depth since it involves the modeling of other entities that are not considered fundamental for the problem. The number of newly activated OT1 CTLs is instead directly estimated on the basis of the quantity of released antigens. The last term () is used to take into account natural death of OT1 CTLs.Melanoma cells:
Equation (4) describes the melanoma cells () behavior in the skin compartment. The first term (() represents a Gompertzian growth [14], whereas the second term denotes killing of by activated OT1 CTLs that are already in the skin compartment (). One of the most important actions accomplished by antibodies () is to boost chemotaxis sensitivity of CTLs () thus enabling better infiltration into the tumor mass. This translates into higher killing rates of melanoma cancer cells by activated CTLs (). This is modeled through the term , where is the maximum killing rate of by , and and () are tuned constants. When the number of is high enough, the term tends towards 1, thus entitling maximum killing rates for . In absence of antibodies, having , the term translates into , which involves lower killing rates.Antigens:
With (5), we describe the tumor-associated antigen () dynamics. Antigens are released in the skin compartment by killed melanoma cells () and are subject to natural degradation (). They can also be captured by APC, which will present the antigen to naïve CTLs. As already stated in (3), capturing of the antigen by APC is not modeled and the number of captured antigens is estimated on the basis of naïve CTLs () that are activated by APC ().Naïve CTLs:
Equation (6) models the behavior of naïve OT1 CTLs (). It is supposed here that these cells are already present in the skin compartment. The term is used to model homeostasis. is the number of circulating naïve CTLs under safe conditions given by the leukocyte formula. If switching of naïve CTLs to activated CTLs occurs, the number of naïve CTLs gets lower. As a consequence of that, the naïve population is repopulated with newborn cells and tends towards at a rate . The second term () models the CTLs state changing from naïve to activated (), thanks to presentation of the antigen by APC.Antibodies (skin compartment):
Antibodies that have reached the skin compartment () are modeled and described by (7). Antibodies in the skin compartment are supposed to be proportional to the number of antibodies in the injection point compartment () with a proportionality constant and a time delay of . They also disappear by stimulating OT1 cells activities and are subject to a natural degradation ( and ).
According to the considered cell populations and interactions, the mathematical model can be then represented by the following system of seven nonlinear delay differential equations:
Since we consider mainly populations that appear in the system as a consequence of treatment administrations (except for melanoma and naïve CTLs), the following Cauchy initial conditions have been set for the equations:
The physical time-step has been chosen equal to 8 hours, and the integration time has been then computed up to days. The reason of this choice is biological, and it is given by the fact that in the in vivo experiment it is not possible to observe relevant biological phenomena in smaller time intervals. In particular, the minimum time required for cell division, which represents one of the most important biological phenomena, is usually not lower than 6–8 hours [15]. This may be not true in other in vivo setups, such as in case of allergies, where the time scale varies from seconds to minutes.
Some parameters appearing in the equations have been estimated from the literature (see Table 2), and from measurements made during the in vivo experiment.
In particular, melanoma growth rates have been estimated from diameter measurements made at different times in melanomas in five untreated mice. Melanomas have two growth phases, radial and vertical [16]. The first phase is represented by a radial growth, in which malignant cells grow in a radial fashion in the epidermis. At later stages, most melanomas progress to the vertical growth phase, in which the malignant cells invade the dermis and develop the ability to metastasize. In this case, we supposed that only the radial growth phase is involved. By supposing a disk-shaped layout for the melanoma and having knowledge of the mean diameter of melanoma cancer cells, the number of cancer cells (in the observed melanomas) has been estimated for all the measurements made. This data has been used with a curve fitting procedure to estimate the unknown parameters needed to model the cancer growth kinetics under the hypothesis of a Gompertzian growth for the tumor, used successfully in our previous experience [13, 17, 18]. The Gompertz law, which is commonly considered suitable for describing populations growths, uses two factors: a constant growth term and a shrink term that increases in time and is related to antiangiogenic factors, giving as a result a sigmoid shape to the curve. The initial number of injected cancer cells , the number of injected OT1 CTLs and antibodies ( and ), and the number of antigen-specific naïve OT1 CTLs under safe conditions () have been chosen into reasonable ranges given by in vivo measurements.
Remaining parameters have been chosen in plausible biological ranges in such a way to reproduce the set of experimental data for the activated OT1 CTLs + Anti-CD37 Ab combined treatment, and counterchecked against the other vaccine scenarios taken into account (Table 2).
4. Results
To reproduce the dynamics of the model, we took into account four treatments that have been tested in vivo: untreated (control), in vitro activated OT1 CTLs, Anti-CD137 monoclonal antibodies, and in vitro activated OT1 CTLs + Anti-CD137 monoclonal antibodies.
Among such treatments, only the last one (treated with activated CTLs + Anti-CD137 monoclonal antibodies) showed complete eradication of the tumor burden, whereas the other treatments remained almost ineffective.
Entities behaviors for all the analyzed treatments are shown in Figure 2. In absence of therapy (blue dashed line), there is no induced immune response and thus the number of melanoma cells grows without any intervention from IS cells, whose plots remain flat.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
The same scenario arises when the Anti-CD137 monoclonal antibodies treatment is administered (see Figure 2, green dotted lines). The mechanism that triggers the IS response is driven by the presence of activated OT1 CTLs that can kill melanoma cells, which release antigens able to further stimulate the IS response.
This fact can be partially seen if the treatment based upon the administration of in vitro activated OT1 CTLs is used (Figure 2, red dot-dashed lines). The plots show some evidence of activated OT1 CTLs in the skin. These cells are able to kill melanoma cells (), which release antigens () that are captured by APCs and then used to promote the differentiation of newborn naïve CTLs () to activated OT1 CTLs (). However, in absence of Anti-CD137 antibodies, which promote both duplication and infiltration into tumor masses of activated OT1 CTLs, recruitment of newborn naïve CTLs is too bland and tardy to stimulate an IS response able to stop the melanoma.
Only the combined action of activated CTLs + Anti-CD137 monoclonal antibodies is able to to contrast the melanoma growth (see Figure 2, black solid line). In this case, both Anti-CD137 antibodies () and activated OT1 CTLs () are injected at day three and then migrate to the skin compartment. Activated OT1 CTLs () are stimulated to duplicate and to infiltrate into tumor masses thanks to the presence of antibodies in the skin compartment (). As a consequence of this, many more melanoma cells () are killed, and copious number of antigens () is released, thus promoting the activation of naïve CTLs () to activated OT1 CTLs (), which further act against melanoma.
To qualitatively compare the model results with the in vivo data, we estimated (using the number of cancer cells) the melanoma surface given by the mathematical model under the same assumptions made for estimating the Gompertz growth parameters. We then compared the time evolution of the melanoma surface given by the mathematical model with the mean surface observed in the in vivo experiment for the following cases: the untreated (control) case and the activated CTLs + Anti-CD137 monoclonal antibodies combined treatment case. The other cases are not presented since they give back the same scenario observed in the untreated case for both the experiments.
Comparison for the untreated case is shown in Figure 3(a). As expected, since we only observe the melanoma growth, the mathematical model perfectly reproduces the in vivo setup. Comparison for the combined treatment case is shown in Figure 3(b). The mean in vivo melanoma behavior is well reproduced by the mathematical model that is able to quantitatively and qualitatively represent the shape of the curve observed in the in vivo experiments. In Figures 3(c) and 3(d), we show absolute and relative differences between in silico and in vivo measurements for the untreated and activated OT1 CTLs + Anti-CD137 combined treatments, respectively. Relative differences are computed using the following metric: .
(a) Untreated
| (b) Act. OT1 + Anti-CD137 |
(c) Untreated
| (d) Act. OT1 + Anti-CD137 |
It must be noted here that the initial gap between the mathematical model and the in vivo experiment measurements visible in Figures 3(b) and 3(d) is mainly due to the different nature of measurements. In the simulations, the total number of cells is always known and it is used to estimate the melanoma surface. In the in vivo experiment, melanoma surface is instead measured on mice skin. Even if the injection of tumor cells is done at time 0, melanoma needs some time to arise and become visible even if the melanoma cells are all already present in mice.
5. Sensitivity Analysis
To understand how the system varies under different parameter values, it is important to analyze the sensitivity of the model to variation of parameters. Classical sensitivity analysis is done usually by varying a given parameter in reasonable ranges and keeping the other ones fixed. The results obtained are obviously dependent on the values of the fixed parameters, and different sets of values may entitle completely different results.
In order to overcome the limits of classical sensitivity analysis, many techniques have been developed. One of these is represented by Partial rank correlated coefficients (PRCC) [23], a statistical sensitivity analysis technique, which computes a partial correlation on rank-transformed data between input and output values, represented in this case by the model parameters (input) and the model entities behaviors (output).
It important to note here that the obtained correlation indexes do not depend on a given set of parameters, and it is therefore possible to estimate how the variations of a parameter may influence the results of the model, no matter what the value of the other parameters is. PRCC returned values varying in and estimated the correlation between input and output parameters. A value near 1 suggests a high (linear) positive correlation, whereas a value near −1 indicates a negative correlation. Values around 0 usually indicate little or no correlation.
Using this methodology in conjunction with the Latin Hypercube Sampling (LHS), which is used to sample the parameters’ space, we analyze the effects of the input parameters most influencing the growth of cancer cells (). We plot PRCCs versus the experiment time for the most important parameters to see how the sensitivity of parameters changes as the system dynamics progresses.
5.1. Impact Variation of the Treatment Quantities and on the Number of Melanoma Cells
The two parameters and refer to the number of activated OT1 effector cells () and Anti-CD137 antibodies () injected as a single boost treatment against the melanoma, according to the in vivo experiment. In Figure 4, we show the PRCC time plots for the two parameters. As expected, just after the injection of the treatment, both the two parameters show a negative correlation with the number of cancer cells. It is worth to note that shows good correlation in particular in the first ten days, whereas a weaker and almost constant correlation is related to . This may suggest that the effects of the treatment are mainly driven by the OT1 effector cells and, even if antibodies are needed to obtain protection (as discussed earlier and shown in Figure 2), a reduction in their quantity may entitle similar treatment effectiveness.
| (a) PRCC time-plot for |
| (b) PRCC time-plot for |
5.2. Impact of Naïve OT1 CTLs Initial Number () and Recovery Rate () on the Number of Cancer Cells
The parameter indicates the initial number of antigen-specific naïve CTLs in the host, and the parameter represents the rate of newborn naïve CTLs generated by thymus selection. The parameter negatively correlates over time with the number of melanoma cells (Figure 5(a)), thus confirming that the action of naïve OT1 T cells induced by the treatment administration is fundamental for treatment effectiveness. Surprisingly it seems that the variation of the rate of introduced newborn naïve CTLs does not influence the number of melanoma cells, since no correlation is shown over time (see Figure 5(b)).
| (a) PRCC time-plot for |
| (b) PRCC time-plot for |
5.3. Impact of Activated OT1-T-Related Parameters on the Number of Cancer Cells
The activated OT1 CTLs in the skin () represent the number of entities that are directly able to kill melanoma tumor cells. It is, therefore, trivial to see that and , which represent the rate of killed melanoma cells by activated OT1 effector cells and the half-life of OT1 effector cells, have a strong negative and positive correlation with the number of cancer cells, respectively, (see Figures 6(a) and 6(b)). A negative correlation is also present for , which represents the maximal duplication rate of () (see Figure 6(c)). However, in this case the correlation remains weaker, thus suggesting how small variations on the cells duplication rates, which can be associated with individual diversity, do not influence considerably the total behavior of cancer cells.
| (a) PRCC time-plot for |
| (b) PRCC time-plot for |
| (c) PRCC time-plot for |
5.4. Impact of Antibody-Related Parameters on the Number of Cancer Cells
The parameter is used to represent the rate of antigens released by killed , whereas the represents the specific antigen half-life. If the former has a strong negative PRCC correlation with the number of (Figure 7(a)), particularly in the first period just after treatment administration, the latter interestingly seems to be not correlated with the number of (Figure 7(b)).
| (a) PRCC time-plot for |
| (b) PRCC time-plot for |
5.5. Impact of the Tumor Growth and Shrink Parameters ( and ) on the Number of Cancer Cells
The and parameters represent the Gompertz growth parameters used to model the growth of melanoma. The range of variation for these two parameters has been set equal to the confidence range given by the curve fitting procedure used to estimate the parameters’ values from in vivo measurements. In this way, we reasonably take into account the possible melanoma growths that may also be observed in the in vivo experiments. PRCC plots over time are represented in Figures 8(a) and 8(b). Even if positive and negative PRCC correlations are somewhat expected, there is a small time window just after the treatment injection time where no correlation occurs. This may be explained by the fact that variations of the melanoma growth rate do not influence substantially (at least initially) the effectiveness of the treatment.
| (a) PRCC time-plot for |
| (b) PRCC time-plot for |
6. Conclusions
We presented a mathematical model, which reproduces the immune response against B16-melanoma induced by the combined administration of activated OT1 CTLs and Anti-CD 137 immunostimulatory monoclonal antibodies. The model uses delay differential equations to reproduce the presence of two different compartments: the injection point compartment where the treatment is administered, and the skin compartment, where the melanoma-Immune system competition occurs.
The model proved to be able to coherently reproduce the in vivo experiment results obtained with four vaccination strategies (untreated, only activated OT1 CTLs, only monoclonal , OT1 CTLs + monoclonal ). Moreover, the model is able to qualitatively and quantitatively reproduce the time dynamics of melanoma under the administration of the combined treatment. Results show that activated CTLs + Anti-CD137 monoclonal combined treatment acts in two ways: directly by activated OT1 CTLs that are able to kill melanoma and antibodies that boost CTLs activities and indirectly by promoting activation of naïve OT1 CTLs thanks to the releasing of melanoma cells antigens.
Among some useful findings, sensitivity analysis underlined the important role of activated OT1 cytotoxic treatment, suggesting that it would be in principle possible to obtain similar effectiveness lowering the number of administered antibodies, which, however, remain fundamental to gain effectiveness. Such kind of suggestions may be useful for optimizing treatment effectiveness and minimizing the risk of side effects.
Future work will be focused on studying analytically a simplified model without delay and in comparing the obtained results with Sim-B16, an agent-based model developed to reproduce the same in vivo experiment [13].
Acknowledgments
The paper was partially supported by the INDAM-GNFM-Young Researchers Project-prot.n.43-Mathematical Modelling for the Cancer-Immune System Competition Elicited by a Vaccine and by University of Catania under PRA grant. The author would like to thank Professor I. Melero and coworkers for supplying data.