If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Corresponding author at: Department of Mathematics, Dartmouth College, 27 N. Main Street, 6188 Kemeny Hall Hanover, NH 03755, USA.
Affiliations
Department of Mathematics, Dartmouth College, 27 N. Main Street, 6188 Kemeny Hall Hanover, NH 03755, USADepartment of Biomedical Data Science, Geisel School of Medicine at Dartmouth, Lebanon, NH 03756, USA
Immune checkpoint therapy is one of the most promising immunotherapeutic methods that are likely able to give rise to durable treatment response for various cancer types. Despite much progress in the past decade, there are still critical open questions with particular regards to quantifying and predicting the efficacy of treatment and potential optimal regimens for combining different immune checkpoint blockades. To shed light on this issue, here we develop clinically-relevant, dynamical systems models of cancer immunotherapy with a focus on the immune checkpoint PD-1/PD-L1 blockades. Our model allows the acquisition of adaptive immune resistance in the absence of treatment, whereas immune checkpoint blockades can reverse such resistance and boost anti-tumor activities of effector cells. Our numerical analysis predicts that anti-PD-1 agents are commonly less effective than anti-PD-L1 agents for a wide range of model parameters. We also observe that combination treatment of anti-PD-1 and anti-PD-L1 blockades leads to a desirable synergistic effect. Our modeling framework lays the ground for future data-driven analysis on combination therapeutics of immune checkpoint treatment regimes and thorough investigation of optimized treatment on a patient-by-patient basis.
]. This has been identified long before a detailed understanding of the mechanisms and components of immune response were understood. Physicians had noticed that solid tumors regressed or even disappeared in patients with skin infections such as Erysipelas caused by Streptococci [
]. While immune system has a vast potential to mount a response against a malignant tumor, in most cases tumor evolves and gains the capability to escape the anti-tumor response. Cancer cells gain the capability to become ‘invisible’ (not recognized by) to the immune cells. Alternatively, they might be able to suppress or reverse the immune response [
]. This can be due to mutations or epigenetic/adaptive changes in cancer cells that increase the expression of immune suppressive pathways.
Immunotherapeutic techniques in cancer are, in essence, a category of methods to reverse the above mechanism of immune resistance (or escape) established by tumor cells. If, through therapeutic interventions, one can stop tumor cells from evading immune response, or inhibit tumor cells from suppressing it, then the immune system, in principle, is capable of eradicating the tumor population. Arguably, one of the most promising strategies in immunotherapy has been the development of immune checkpoint blockade antibodies [
]. Immune checkpoint blockade therapies are a new therapeutic method to make tumor cells ‘visible’ to the immune cells and re-activate and strengthen tumor-specific immune response. This is achieved by blocking the signals/proteins or corresponding ligands that contribute in the immune recognition/activation or immune suppression pathways.
There have been several classes of immunotherapy drugs that function as immune checkpoint blockades. Antibodies that block cytotoxic T lymphocyte-associates protein 4 (CTLA-4) or programmed cell death-1 (PD-1) or its ligand-1 (PD-L1) are the most well-known immune checkpoint treatments. There has been success in patient survival in different cancers mainly, melanoma, non-small cell lung cancer and non-hodgekins lymphoma, among others. Currently, there are several immune checkpoint blockades that are approved by the FDA, such as Ipilimumab, Tremelimumab (anti-CTAL-4), Nivolumab (anti-PD-1) and Avelumab (anti-PD-L1). There are more immune checkpoint antibodies under development for a broader range of tumor types [
The field of immune checkpoint therapy has joined the ranks of surgery, radiation, chemotherapy, and targeted therapy as a pillar of cancer therapy. These drugs represent a radical and disruptive change in cancer therapy. With the exception of anti-PD-L1 agents, these drugs do not directly attack tumor cells. But instead, they target immune cells and immune system at large. More importantly, instead of activating the immune system, they suppress inhibitory pathways that induce immune resistance and block effective anti-tumor response by the immune system. Immune checkpoint therapy, with anti-CTLA-4 having longer follow-up than other agents, leads to durable clinical responses that can last a decade and more, but only in a fraction of patients [
Unlike anti-PD-1 and anti-CTLA4 antibodies, anti-PD-L1 agents do target tumor cells. While anti-PD-1 antibodies block the pathways that suppress anti-tumor response in T-cells, anti-PD-L1 antibodies block similar pathways in tumor cells. They make tumor cells vulnerable to anti-tumor response. Successes have been seen in preclinical and clinical trials based on the combination of anti-CTLA4 and anti-PD-1 agents [
]. As such, there are high expectations that a highly effective immunotherapeutic strategy can be devised using potential combinations of these different agents. There are ongoing studies to identify immune biomarkers with which one can predict the efficacy of treatment for select patients [
]. But the complexity of the immune system has made this task more difficult so far.
For an adaptive anti-tumor response to initiate, cytotoxic T cells that are responsible for killing tumor cells need to be activated in the first place. T cell activation is commonly specific to tumors with particular genetic makeup. Cytotoxic T cells are denoted with their protein marker CD8 (vs CD4 for helper T-cells). Upon interaction with tumor cells – in the tumor microenvironment – T cells can be activated. This also can happen indirectly when T cells encounter other immune cells that carry the tumor antigen. These are so-called antigen presenting cells or APC. APC’s are commonly dendritic cells that have absorbed antigens from dying tumor cells, but they can be other immune cell type as well. [
Immune activation is done through a two-signal model (see Fig. 1). The first and main signal is the binding of T cell receptor (TCR) by the MHC protein on APC or tumor cells. Second signals, or co-stimulating signals, are commonly CD28 on T cells that interact with B7-1/2 on tumor cells. One major function of a second signal is to stabilize the first signal binding. Upon activation T cells undergo clonal expansion with help of the IL-2 cytokines, which is a growth and differentiation factor. Beside clonal expansion of T cells, tumor-specific APC cells are amplified in numbers as well. This is done through an intricate set of mechanisms that involve CD8 and CD4 T cells as well regulatory T cell (Tregs). This response can last by establishing memory T cells [
Fig. 1Schematic illustration of key tumor-immune molecular interactions. TCR/MHC is the first immune activation signal and CD28/B7-1/2 is the second signal. CTLA-4 competitive inhibition with CD28 and PD-1/PD-L1 engagement are inhibitory signals for T cells.
The mechanism of adaptive (or constitutive) immune resistance by tumor is understood by the over-expression of PD-L1 ligand on tumor cells. Upon engaging with its PD-L1 ligand, PD-1 acts as an immune suppressive signal [
]. When tumor antigen-specific T cells recognize their cognate antigen expressed by cancer cells, signaling through TCR leads to production of interferons and expression of regulatory receptors such as PD-1. Interferons tune the expression of immune-suppressive factors on tumor cells, including PD-L1 ligand. Upon engaging with PD-L1, PD-1 protein on a T cell works towards suppression of the adaptive immune response [
]. Cancer cells use these adaptive immune suppression programs that are set to limit the immune and inflammatory responses, to their benefit. The mechanism of immune resistance can be constitutive due to activation of oncogenic pathways that also lead to expression of PD-L1 ligand [
]. These key tumor-immune molecular interactions are depicted in Fig. 1.
Anti-PD-1/PD-L1 immune checkpoint blockade therapies are aimed to reverse the immune resistance mechanism by making tumor cells visible to immune cells. This effectively leads to re-launching an immune response that is suppressed by PD-L1+ tumor cells. As described above, interaction of programmed cell death-1 (PD-1) protein with its ligand on a tumor cell sets off inhibition programs that suppresses the adaptive immune response. Immune checkpoint PD-1/PD-L1 blockades effectively inhibit PD-1/PD-L1 engagement.
As aforementioned, another category of immune checkpoint blockades targets the CTLA-4 protein. The CTLA-4 is homologous to the T cell co-stimulatory protein, CD28, and both molecules bind B7-1/2 proteins. CTLA4 can outcompete CD28 in binding. It transmits an inhibitor signal for T cell activation, whereas CD28 transmits a stimulatory signal. CTLA-4 inhibition happens in T cell priming sites, i.e., lymph nodes [
]. Blocking CTLA-4 on T cells increases the immune response and amplifies the activation process. CTLA-4 blockade, furthermore, increases T cell motility and renders more T cells to move into the tumor microenvironment as well. Such an increase in motility can significantly help in T cell infiltration process and anti-tumor response [
In what follows, we aim to address fundamental open questions regarding the dynamics of the immune checkpoint blockades using a quantitative modeling approach. Specifically, potential biomarker mechanisms behind the choice of therapeutic strategies using either of anti-PD-1/PD-L1 agents or both (monotherapy versus combination therapy) are not well-understood [
]. Moreover, there are still critical open questions with particular regards to quantifying and predicting the efficacy of treatment and potential optimal regimens for combining different immune checkpoint inhibitors. Our approach follows the mathematical oncology paradigm [
], and our dynamical systems modeling can be used to inform the rational and personalized development of cancer immunotherapy using checkpoint blockades and their potential combinations in order to improve response rate and reduce resistance.
2. Dynamical systems model of adaptive immune resistance through PD-1/PD-L1 axis
A minimal dynamical model of immune-tumor interaction is recently suggested in Ref. [
Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: a mathematical and computational examination.
]. The model considers a population of tumor cells, and immune cells (effector cells ), in the tumor microenvironment. We assume a constant supply of (primed) effector cells from lymph nodes (). Upon interaction with tumor cells (either direct or indirect through APC cells) immune cells receive stimulatory signals to increase their proliferation. We assume this to happen with a rate per unit time and per capita. Tumor cells are killed with a rate upon interaction with effector cells. Assuming a logistic growth for tumor population in the absence of immune interaction, we have [
Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: a mathematical and computational examination.
where is the linear growth rate (fitness) of tumor cells, and is the inverse carrying capacity. is the average death rate of effector cells. The solution for the above equation shows bi-stability as the treatment parameters changes, and one can show the exact combinations of and values that are needed for effective immune control of tumor growth (characterized by low abundance of tumor cells or complete eradication) [
Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: a mathematical and computational examination.
In the following we generalize the above model for immune checkpoint combination therapies. In this paper, we focus on anti-PD-1 and anti-PD-L1 antibodies and their effects on suppressing the immune resistance.
We divide the population of tumor infiltrating T cells into low- and high- expression of PD-1 protein, namely, PD-1- and PD-1+ T cells. Similarly, we focus on two subpopulations of PD-L1+ and PD-L1- tumor cells. We denote tumor PD-L1- (+) with () and PD-1- (+) effector cells with (), respectively. or for are dynamical variables and their values can change in time depending on how immune-tumor interaction dynamics is going and how other tumor microenvironmental factors, including drug efficacies, affect them.
To be concrete, we use the following set of biologically plausible assumptions to construct our model:
•
The primed T-cells are supplemented at a constant rate from the lymph nodes into the tumor microenvironment. This rate is denoted with .
•
The majority tumor cell population is initially PD-L1- and sensitive to the anti-tumor immune response.
•
We assume a mechanism of immune resistance where with rate, , PD-L1- tumor cells transform into PD-L1+ cells. Similarly, the effector cell population transforms from PD-1- into PD-1+ with a rate .
•
Any T-cell subtype (PD-1+/-) can in interact with any tumor cell subtype (PD-L1+/-). Their interaction strengths are modeled into a matrix . Parameter represents the stimulation/inhibition level of type effector cells by type tumor cells. For PD-1- cells and for PD-1+ ones . Using similar notations, applies to PD-L1-/+ subtypes.
•
PD-L1+ subpopulation in the tumor drives the immune resistance. Upon interaction with PD-1- immune cells a PD-L1- tumor cell stimulates the response (), while PD-L1+ tumor cells suppress anti-tumor response of PD-1+ T cells ().
Therefore, we set , and while represents the immune suppression. The anti-tumor activity of the immune cell population is modeled with a matrix of killing rate, . is the death rate induced by a type immune cell on a type tumor cell (). The PD-L1+ subpopulation is assumed to have lesser death rate when encountering effector cells, especially PD-1+ T cells. The above dynamics is described by the following system of equations for four subpopulations of PD-1-/+ immune cells, and , and PD-L1-/+ tumor cells, and :
(2)
Here, is the recruiting rate of the effector cells from the lymph nodes to the local tumor microenvironment. are the death rates of PD-1+/- effector cells. is the linear growth rate of tumor and is the inverse carrying capacity for tumor cells. PD-1- T cells transform into PD-1+ with rate while PD-L1- cell become PD-L+ with rate . Values of and together quantify how fast the immune resistance is developed against an anti-tumor immune response. The above mechanisms are schematically presented in Fig. 2.
Fig. 2Model scheme of tumor-immune interactions in the microenvironment. Population of effector cells are divided into two subpopulation of PD-1+ and PD-1- types. Two tumor subpopulations of PD-L1+ and PD-L1- are distinguished in a similar fashion. Tumor and immune cell subpopulations interact based on matrices and . () determine how the population of tumor cells affect the immune response. () quantify the death rate induced by effector cells on tumor cells. The interaction-matrix elements, and are not the same. For example, PD-L1+ tumor cells are able to suppress the immune response by PD-1+ immune cells, and thus the strength of can be negligible. The mechanism of adaptive immune resistance is modeled as the increase in PD-1 and PD-L1 expression levels, which can be quantified by the transition rates and .
We can rewrite the Eq. (2) by rescaling tumor population with carrying capacity and immune population with .
(3)
This leads to rescaling of model parameters:
(4)
where the subscript index . Eq. (2) can be rewritten in terms of rescaled parameters as,
(5)
Such rescaling reduces the effective number of model parameters by two and thus facilitates our sensitivity analysis of model parameters. We could further rescale time, , in the derivative, and write it in units of average life time of immune cells, . This would reduce the model parameters to one less. However, we keep variable as we want to keep the time units in specific units such as days (or proper tumor cell generation times).
We summarize the key parameters and their biological interpretations in our model as follows:
Pre-existing immune resistance parameters, (, ). These basically quantify how impaired the immune response is due to -interferon pathway activation and up-regulation of PD-1 on T cells and PD-L1 on tumor cells.
Tumor growth parameters (, ). We already re-scaled the model system to absorb the tumor inverse carrying capacity into other parameters. Thus the main indicator of tumor growth and aggressiveness is the value of its linear growth rate .
Tumor-immune interaction matrices (). The value of identifies the set of death rates values for the tumor cells due to immune killing (). While determines how the strength of immune response (abundance of effector cells) is mediated by the interaction of PD-1+/- effector cells with PD-L1+/- tumor cells.
At this stage the above model can describe how quickly the tumor develops resistance to an anti-tumor immune response. Subpopulation of PD-L1- tumor cells is vulnerable to immune response and can be effectively suppressed by the presence of PD-1- effector cells. This means that PD-L1- cells can stimulate T cells response and also are lysed after an encounter with tumor specific T cells. PD-L1+ cells, however, have the capability to suppress the immune response. This effect in fact is not uniform and PD-L1+ cells effectively suppress immune response in PD-1+ cells and less so in PD-1- immune cells.
The main mechanism behind the above dynamics is that in the interactions between PD-1- T cells and PD-L1- tumor cells, the tumor cells are vulnerable to the tumor-specific response from PD-1- effector cells and thus immune response is effective. In contrast, owing to the checkpoint pathway engagement via PD-1+/PD-L1+, the immune response can be greatly impaired over time when PD-L1- tumor cells gain the capability to switch into PD-L1+ status and thus will be able to resist the tumor-specific response and even suppress it.
3. Results
3.1 Evolution of adaptive immune resistance
We first focus on the tumor-immune interaction dynamics in the absence of anti-PD-1 and anti-PD-L1 treatments. The dynamical model, as given in Eq. (5), describes the tumor growth over time, while effector cells respond to suppress the tumor growth. The anti-tumor activity of effector cells increases upon recognizing PD-L1- tumor cells while PD-L1+ tumor cells suppress PD-1+ effector cell activity upon interaction. The transformation from PD-L1- to PD-L1+ compartments (through up-regulation of PD-L1 ligand) results in tumor immune resistance. The parameter quantifies how tumor population transforms into immune-suppressive subpopulation of PD-L1+. In the meantime, at rate , PD-1- immune cells become PD-1+ cells with poor anti-tumor characteristic. The finite value of describes a gradual change of PD-L1- to PD-L1+, thereby indicating the level of PD-L1 expression. As a consequence of adaptive immune resistance, anti-tumor response is suppressed due to high proportion of PD-L1+ tumor cells and low abundance of PD-1- T cells (Fig. 3).
Fig. 3The impact of adaptive immune resistance on tumor burden and immune cell counts. The top panel (a) shows the total number of effector cells (tumor infiltrating cells) as a sum of PD-1+ and PD-1- cells. The bottom panel (b) shows the tumor population (total of PD-L1+ and PD-L1- cells). For different plots we have changed the value of while the value of is kept the same as , . The abundance of effector cells and tumor burden are shown as rescaled according to Eq. (3), and relevant model parameter values are rescaled according to Eq. (4). Rescaled model parameters are , , , , and , .
Without loss of generality, we assume the values of immune stimulation/inhibition, are the same order of magnitude. We model suppressive effect of PD-L1+ tumor cells on immune cells by setting to a negative value. Anti-tumor activity of effector cells is parametrized by prescribing values of . In the example shown in Fig. 3, we use and . The value of indicates that PD-1+ T-cells has no effector activity on PD-L1+ tumor cells. The tumor-immune interaction parameter, are set to , and . The positive (negative) values of mean stimulatory (or inhibitory) response due to immune-tumor interaction, respectively. For the rest of parameters we used and values of . In Fig. 3, we plot the total effector cell count, , and the total tumor size, as a function of time, using the initial condition . The higher the values of and are, the more pronounced adaptive immune resistance prior to treatment; tumor-immune interactions are largely driven by PD-1+/L1+ interactions with sizable tumor burden. Specifically, the equilibrium tumor size (obtained when the above model reaches the steady state) as a function of immune resistance parameters is shown in Fig. 4. Same model parameters as in Fig. 3 are used, and values of and are varied independently. As adaptive immune resistance is modulated by the -interferon pathway activation and up-regulation of PD-1 on T cells and PD-L1 on tumor cells, complete immune escape occurs only for high levels of PD-1 and PD-L1 expression, that is, large values of and (see Fig. 4).
Fig. 4Heatmap plot of tumor steady-state size as a function of immune resistance parameters, , and . We assume no treatment is applied, . Tumor burden is shown as rescaled according to Eq. (3), and relevant model parameter values are rescaled according to Eq. (4). The immune inhibition/suppression parameters are and , and the tumor killing rates by PD-1+/- cells are and . And we set and the same as in Fig. 3.
To obtain closed-form results beside numerical simulations, we can simplify the immune-tumor interaction matrix by writing four components in terms of a single parameter : . A simple analytical expression can be found for the steady-state tumor size along the diagonal :
(6)
For small-, the formula above can be expressed in a Taylor expansion:
(7)
Notably, this simplified formula provides us a clear and intuitive picture about the relationship between tumor burden and levels of adaptive immune resistance. The equilibrium tumor burden monotonically increases with the adaptive immune resistance parameters as the coefficient of the first order expansion in is positive. High expression levels of PD-1 and PD-L1 lead to severely impaired immune response and as a consequence, yield high tumor burden (as shown in Fig. 3).
3.2 Anti-PD-1 versus anti-PD-L1 monotherapy
We extend the above model of adaptive immune resistance by incorporating immune checkpoint therapy using PD-1 and PD-L1 blockades (see Fig. 5). Put simply, the effect of immune checkpoint blockade treatments is exerted on the PD-1+ (PD-L1+) subpopulation by neutralizing and blocking PD-1 (and PD-L1) receptors and further reversing them into PD-1- (PD-L1-) subpopulation, respectively. Such conversion happens after immune cells (tumor cells) are bound with anti-PD-1 (anti-PD-L1) antibodies and thus the immune checkpoint pathway is successfully blocked. In Eq. (5) developing adaptive immune resistance is modeled through transition rates and , which determine how quickly PD-1 (PD-L1) expression levels of cells transition from lower (-) to higher values (+). As such, we account for the effect of immune checkpoint therapy by corresponding reverse transition rates. That is, in the presence of anti-PD-1 antibodies, PD-1+ immune cells transition to PD-1- compartment with the rate . The bulk parameter reflects the overall neutralization (blockade) efficacy determined by the underlying concentration-dependent binding kinetics [
]. Similarly, we denote by the reverse transition rate for tumor cells transitioning from PD-L1+ to PD-L1- compartments in the presence of anti-PD-L1 treatment.
Fig. 5Model scheme of tumor-immune interactions in the presence of immune checkpoint blockades. Similar to Fig. 2, the effector cell population is divided into PD-1+ and PD-1- compartments and tumor cells are divided into two compartments of PD-L1+ and PD-L-. Immune-escape mechanisms cause transition from PD-L1- to PD-L1+ compartment, with rate . We assume similar transition from PD-1- to PD-1+ for effector cells with rate . Anti-PD-1 and anti-PD-L1 antibodies cause a reversal of this with corresponding rates and .
Incorporating the anti-PD-1/anti-PD-L1 treatment into the Eq. (5), we obtain the following system of differential equations:
(8)
Closed-form analytical solutions can be obtained for the steady-state tumor size, (that is, ), if we apply the same assumption as in Eq. (6), , as well as for anti-PD-1 treatment. We have,
(9)
Similar results can be obtained in the case of PD-L1 treatment.
It is worthy noting that while treatment efficacy is quantified mainly by model parameters and in this paper, the parameter characterizing the tumor-immune interaction via the checkpoint pathway PD-1+/PD-L1+ can be affected to some extent during treatment. PD-1+ T cells are the primary target of anti-PD-1 therapy, and after PD-1 blockade, they effectively become PD-1- cells that are responsive again to tumor-immune stimulatory signals. As PD-1 blockade antibody can also competitively weaken the binding activities between PD-1+/PD-L1+, we assume anti-PD-1 treatment is able to reduce the magnitude of immune suppression signaling . Similarly, we assume anti-PD-L1 therapy targets PD-L1+ tumor cells, and meanwhile can also reduce the value of tumor-immune interaction matrix . Hence, the parameter value of during treatment is uniformly altered as a secondary effect of monotherapy or combination therapy.
We perform bifurcation analysis of equilibrium tumor burden, , with respect to changes in the anti-tumor immune activity, . In Fig. 6, we plot the steady-state solutions (fixed points) for anti-PD-1 and anti-PD-L1 treatments, respectively. We also vary values of while other immune-tumor interaction parameters are kept constant . Values of tumor killing rates are parametrized as . In other words, the PD-1+/PD-L1+ interaction leads to much weaker cytotoxic activity from T-cells. For both anti-PD-1 and anti-PD-L1 monotherapies, we observe an interesting bistability behavior for and for intermediate levels of anti-tumor activity (Fig. 6). One potential implication of bistablity phenomena is that the existence of an interior unstable equilibrium allows the occurrence of hysteresis. Namely, the rebound of tumor mass after treatment stops has to follow the lower branch instead of the upper branch along which it plummeted during treatment, and thus tumor can be controlled at low abundance until the anti-tumor activity is decayed below a certain threshold.
Fig. 6Bifurcation analysis of equilibrium tumor burden with respect to immune killing rate under monotherapy with anti-PD-1 and anti-PD-L1. tumor fixed point, , for various values of (tumor killing rate) for PD-L1 mono-therapy as immune-tumor interaction parameter is varied due to treatment. Tumor burden is shown as rescaled according to Eq. (3), and relevant model parameter values are rescaled according to Eq. (4). Model parameters are set as , , , , , (a) anti-PD-1: , and (b) anti-PD-L1: . And we set and the same as in Fig. 3.
To further reveal the subtle difference between anti-PD-1 and anti-PD-L1 monotherapies, we plot the time evolution of tumor burden during treatment in Fig. 7. For anti-PD-1 monotherapy, we set while . Similarly, for anti-PD-L1 monotherapy, while . To compare the efficacy of the two treatments we keep the remaining model parameters exactly the same. In the examples shown in Fig. 7, we note that the treatment outcome of anti-PD-L1 therapy is better than anti-PD-1 therapy for the specific parameter choices ( and ) that indicate high PD-1 expressions prior to treatment (also see Fig. 8a). More generally, these exists an interesting crossover of the residual tumor curves (as shown in Fig. 8b and c): anti-PD-1 is slightly more effective for low neutralizing efficacies, but anti-PD-L1 works better for high neutralizing efficacies. These theoretical results suggest that PD-L1 expression level alone is insufficient to determine which anti-PD-1 or anti-PD-L1 blockade to be more effective [
]. However, for a wide range of model parameters, it appears that anti-PD-L1 therapy is more effective than anti-PD-1 therapy even for equally potent PD-1 and PD-L1 blockades [
Fig. 7Comparison of anti-PD-L1 (upper panel) and anti-PD-1 (lower panel) monotherapy (lower panel). The relative tumor size (to carrying capacity) is used as a quantitative measure of response dynamics and treatment outcomes. Treatment begins with tumor at its steady state in the absence of immune checkpoint therapy. In this example, anti-PD-L1 monotherapy seems to be more effective (i.e., yielding smaller residual tumor burden) than anti-PD-1 monotherapy, because of higher expression level of PD-1 than that of PD-L1 (). Tumor burden is shown as rescaled according to Eq. (3), and relevant model parameter values are rescaled according to Eq. (4). Other model parameters used are: , , , , (a) anti-PD-L1: , and (b) anti-PD-1: . And we set and the same as in Fig. 3.
Fig. 8Residual tumor burden as a function of treatment efficacy of monotherapy, anti-PD-1 versus anti-PD-L1 . In panels (a)-(c), anti-PD-1 and anti-PD-L1 monotherapies are compared by assuming the same degree of efficacy expressed in terms of neutralizing PD-1+ and PD-L1+ cells: for anti-PD-L1 we have for anti-PD-1 we have , where denotes the x-axis. Depending on the expression levels of PD-1 and PD-L1 status (immune resistance parameters and ), there exists an interesting crossover of the two tumor burden curves as shown in panels (b) and (c): anti-PD-1 is slightly more effective for low neutralizing efficacies, but anti-PD-L1 seems to work better for high neutralizing efficacies. Rescaled model parameters are, , , , , , , (a) , (b) , (c) .
3.3 Anti-PD-1 and anti-PD-L1 combination immunotherapy
It is straightforward to investigate combination immunotherapy based on Eq. (8). Immune checkpoint therapy combining anti-PD-1 and anti-PD-L1 blockades is charaterized by the parameter space , which quantifies the treatment efficacy of each checkpoint blockade alone. A natural question is whether there exists Loewe synergy between anti-PD-1 and anti-PD-L1 treatments [
In Fig. 9, we show a contour plot of residual tumor burden as a function of . Along each contour line, the residual tumor burden is constant for the combination of values of . Let us denote by () the critical efficacy of monotherapy which is needed to reach a fixed value of residual tumor burden. If anti-PD-1 and anti-PD-L1 treatments satisfy the Loewe additivity, that is, for any combination of and satisfying
(10)
we have the exactly same treatment effect where residual tumor burden always reaches the same value. However, the contour lines in Fig. 9 are not straight but instead convex downward, which suggest synergy between anti-PD-1 and anti-PD-L1 blockades. Figure 10 further demonstrates that neither of monotherapies is able to yield greater tumor shrinkage than combination therapy.
Fig. 9Synergistic effect of anti-PD-L1 and anti-PD-1 combination immunotherapy. Residual tumor burden is shown as a contour plot of the parameter space , which quantifies anti-PD-1 and anti-PD-L1 treatment efficacies. The Loewe additivity is indicated by the straight line , where and are the critical values for each corresponding monotherapy to have the same treatment effect (namely, residual tumor burden). The contour lines are convex downward, suggesting synergy between anti-PD-L1 and anti-PD-1 treatments. Tumor burden is shown as rescaled according to Eq. (3), and relevant model parameter values are rescaled according to Eq. (4). Model parameters are .
Fig. 10Immunotherapy combining PD-1 and PD-L1 blockades renders better tumor shrinkage than monotherapy. For the sake of proper comparison in accordance with the Loewe synergy
, we restrict parameter choices of and for combination therapy along the line (see Fig. 9), , where and are the critical values for each corresponding monotherapy to have the same treatment effect (namely, residual tumor burden). We assume tumor reaches steady state prior to treatment, and the relative tumor size is plotted as a function of time after the treatment starts at . Monotherapy: anti-PD-1 , anti-PD-L1 , combination therapy: , other rescaled model parameters: , , , , , , , .
]. While these potentially curative cancer therapies are rapidly being developed and tested, a major barrier is the lack of quantitative models of their efficacy. To address this issue, mathematical modeling of cancer-immune interactions and of immunotherapy has been a topic of interest and primary significance [
]. Building on and integrating these aforementioned models, here we focus on quantifying and predicting the efficacy of immune checkpoint PD-1/PL1 blockades treatment and their combinations.
An open question in the field is, other things being equal, which monotherapy, anti-PD-1 or anti-PD-L1 is more effective in tumor control and eradication provided they have the same checkpoint blockade efficacy [
Sequential blockade of PD-1 and PD-L1 causes fulminant cardiotoxicity: from case report to mice model validation.
Ann Oncol.2018; 29
]. Our results show that depending on PD-1 and PD-L1 expression levels, anti-PD-1 or anti-PD-L1 treatment can be more effective than the counterpart. However, in most cases anti-PD-L1 treatment is seen to be a more effective treatment in a wide range of model parameters. Altogether, our work demonstrates that PD-L1 expression level alone is insufficient to determine which anti-PD-1 or anti-PD-L1 blockade to be more effective [
]. Our theoretical modeling predicts synergy between anti-PD-1 and anti-PD-L1 treatments. This prediction is in line with a recent work, which shows that PD-1/PD-L1 blockade combination leads to better treatment outcomes among metastatic breast cancer patients [
], and develop theory-informed combination therapies to overcome immunotherapy resistance and improve response.
In this work, the switching rate of effector cells, (and similarly tumor cells, ) to increase the expression level of PD-1 (PD-L1 ligand, respectively) is assumed to be constant as a proof-of-principle model with a parsimonious number of parameters. It is promising for future work to incorporate detailed molecular reaction kinetics and dynamics associated with the PD-1/PD-L1 regulation pathway [
In silico simulation of a clinical trial with anti-CTLA-4 and anti-PD-L1 immunotherapies in metastatic breast cancer using a systems pharmacology model.
A QSP model for predicting clinical responses to monotherapy, combination and sequential therapy following CTLA-4, PD-1, and PD-L1 checkpoint blockade.
]. These meaningful extensions may further improve our understanding of the subtle differences in treatment outcomes between anti-PD-1 and anti-PD-L1 immune checkpoint blockades.
Recently, there has been growing interest and effort in advancing personalized medicine in cancer immunotherapy [
], it is particularly important for immunotherapy to account for individual heterogeneity for rational and informed treatment decisions, for example, regarding which checkpoint inhibitors to use, as well as when and how to use their potential combinations so as to improve response and reduce resistance. As such, treatment decisions and regimes ought to be highly tailored to individual patients based on their cancer specific and immune-related biomarkers. In this context, mathematical models will have an important role in furthering personalized immunotherapy [
]: parameters of mechanistic models have clear biological meanings and can be tuned and calibrated using real patient data in a truly personalized fashion; they can be used to quantify for the uncertainty in health outcomes due to stochasticity; they can be readily refined by taking advantage of available data and knowledge to ensure their clinic relevance [
]. Synergistically combining in-silico modeling with real patient data will greatly facilitate the way novel cancer immunotherapeutic strategies are conceived, tested, and understood.
In summary, we have proposed a proof-of-concept, dynamical systems model of tumor-immune interactions along the PD-1/PD-L1 axis. Based on this model, we have studied the evolution of adaptive immune resistance and subsequent response to anti-PD-1/PD-L1 treatments. Depending on model parameters describing immune resistance, tumor growth and treatment efficacies, our modeling framework provides mechanistic insights into quantifying and characterizing conditions for the success (complete or partial response) or failure of immunotherapy using checkpoint inhibitors and their potential combinations. Further parameterized with patient-specific tumor and immune biomarker data [
This work is supported by the NIH COBRE Program (grant no. 1P20GM130454). F.F. is grateful for the generous financial support by the Bill & Melinda Gates Foundation (award no. OPP1217336) and the Neukom CompX Faculty Grant.
Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: a mathematical and computational examination.
In silico simulation of a clinical trial with anti-CTLA-4 and anti-PD-L1 immunotherapies in metastatic breast cancer using a systems pharmacology model.
A QSP model for predicting clinical responses to monotherapy, combination and sequential therapy following CTLA-4, PD-1, and PD-L1 checkpoint blockade.