Immune checkpoint therapy modeling of PD-1/PD-L1 blockades reveals subtle difference in their response dynamics and potential synergy in combination

  • Kamran Kaveh
    Affiliations
    Department of Mathematics, Dartmouth College, 27 N. Main Street, 6188 Kemeny Hall Hanover, NH 03755, USA
    Search for articles by this author
  • Feng Fu
    Correspondence
    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, USA

    Department of Biomedical Data Science, Geisel School of Medicine at Dartmouth, Lebanon, NH 03756, USA
    Search for articles by this author
Open AccessPublished:July 29, 2021DOI:https://doi.org/10.1016/j.immuno.2021.100004

      Abstract

      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.

      Graphical abstract

      Keywords

      1. Introduction

      Immune system is shown to have the potential to activate a response that can eradicate a tumor [
      • Ribas A.
      • et al.
      Releasing the brakes on cancer immunotherapy.
      ,
      • Ribas A.
      • Wolchok J.D.
      Cancer immunotherapy using checkpoint blockade.
      ,
      • Sharma P.
      • Allison J.P.
      Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential.
      ]. 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 [
      • Kucerova P.
      • Cervinkova M.
      Spontaneous regression of tumour and the role of microbial infection–possibilities for cancer treatment.
      ]. 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 [
      • Ribas A.
      Adaptive immune resistance: how cancer protects from immune attack.
      ,
      • Kalbasi A.
      • Ribas A.
      Tumour-intrinsic resistance to immune checkpoint blockade.
      ]. 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 [
      • Ribas A.
      • Wolchok J.D.
      Cancer immunotherapy using checkpoint blockade.
      ,
      • Sharma P.
      • Allison J.P.
      Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential.
      ,
      • Wei S.C.
      • Duffy C.R.
      • Allison J.P.
      Fundamental mechanisms of immune checkpoint blockade therapy.
      ,
      • Littman D.R.
      Releasing the brakes on cancer immunotherapy.
      ]. 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 [
      • Wei S.C.
      • Duffy C.R.
      • Allison J.P.
      Fundamental mechanisms of immune checkpoint blockade therapy.
      ,
      • Tang J.
      • Shalabi A.
      • Hubbard-Lucey V.
      Comprehensive analysis of the clinical immuno-oncology landscape.
      ].
      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 [
      • Sharma P.
      • Allison J.P.
      The future of immune checkpoint therapy.
      ].
      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 [
      • Wei S.C.
      • Levine J.H.
      • Cogdill A.P.
      • Zhao Y.
      • Anang N.-A.A.
      • Andrews M.C.
      • Sharma P.
      • Wang J.
      • Wargo J.A.
      • Pe?er D.
      Distinct cellular mechanisms underlie anti-CTLA-4 and anti-PD-1 checkpoint blockade.
      ,
      • Wei S.C.
      • Anang N.-A.A.
      • Sharma R.
      • Andrews M.C.
      • Reuben A.
      • Levine J.H.
      • Cogdill A.P.
      • Mancuso J.J.
      • Wargo J.A.
      • Pe?er D.
      • et al.
      Combination anti–CTLA-4 plus anti–PD-1 checkpoint blockade utilizes cellular mechanisms partially distinct from monotherapies.
      ]. 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 [

      Willis J.C., Lord G.M.. Immune biomarkers: the promises and pitfalls of personalized medicine. Nat Rev Immunol2015. 15, 5, 323–329,

      ]. 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. [
      • Wei S.C.
      • Duffy C.R.
      • Allison J.P.
      Fundamental mechanisms of immune checkpoint blockade therapy.
      ,
      • Chaplin D.D.
      Overview of the immune response.
      ].
      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 [
      • Nicholson L.B.
      The immune system.
      ].
      Fig. 1
      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 [
      • Johnson R.M.G.
      • Dong H.
      Functional expression of programmed death-ligand 1 (B7-H1) by immune cells and tumor cells.
      ]. 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 [
      • Wu Y.
      • Chen W.
      • Xu Z.P.G.
      • Gu W.
      PD-L1 distribution and perspective for cancer immunotherapy–blockade, knockdown, or inhibition.
      ]. 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 [
      • Parsa A.T.
      • Waldron J.S.
      • Panner A.
      • Crane C.A.
      • Parney I.F.
      • Barry J.J.
      • et al.
      Loss of tumor suppressor PTEN function increases B7-H1 expression and immunoresistance in glioma.
      ,
      • Akbay E.A.
      • Koyama S.
      • Carretero J.
      • Altabef A.
      • Tchaicha J.H.
      • Christensen C.L.
      • et al.
      Activation of the PD-1 pathway contributes to immune escape in EGFR-driven lung tumors.
      ,
      • Atefi M.
      • Avramis E.
      • Lassen A.
      • Wong D.J.
      • Robert L.
      • Foulad D.
      • et al.
      Effects of MAPK and PI3K pathways on PD-L1 expression in melanoma.
      ]. However, it seems that adaptive interferon-inducible expression of PD-L1 is more common than the constitute expression in many cancers [
      • Taube J.M.
      • Anders R.A.
      • Young G.D.
      • Xu H.
      • Sharma R.
      • McMiller T.L.
      • et al.
      Colocalization of inflammatory response with B7-H1 expression in human melanocytic lesions supports an adaptive resistance mechanism of immune escape.
      ,
      • Tumeh P.C.
      • Harview C.L.
      • Yearley J.H.
      • Shintaku I.P.
      • Taylor E.J.
      • Robert L.
      • Chmielowski B.
      • Spasic M.
      • Henry G.
      • Ciobanu V.
      • et al.
      PD-1 blockade induces responses by inhibiting adaptive immune resistance.
      ]. 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 [
      • Ribas A.
      • Wolchok J.D.
      Cancer immunotherapy using checkpoint blockade.
      ]. 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 [
      • Pentcheva-Hoang T.
      • Simpson T.R.
      • Montalvo-Ortiz W.
      • Allison J.P.
      Cytotoxic T lymphocyte antigen-4 blockade enhances antitumor immunity by stimulating melanoma-specific T-cell motility.
      ].
      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 [
      • Linhares A.D.S.
      • Battin C.
      • Jutz S.
      • Leitner J.
      • Hafner C.
      • Tobias J.
      • et al.
      Therapeutic PD-L1 antibodies are more effective than PD-1 antibodies in blocking PD-1/PD-L1 signaling.
      ,
      • Chen Y.-J.
      • Huang W.-C.
      • Liu S.-Y.
      • Ko C.C.
      Sequential blockade of PD-1 and PD-L1 causes fulminant cardiotoxicity: from case report to mice model validation.
      ,
      • Shen X.
      • Zhao B.
      Efficacy of PD-1 or PD-L1 inhibitors and PD-L1 expression status in cancer: meta-analysis.
      ]. 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 [
      • Basanta D.
      • Scott J.G.
      • Fishman M.N.
      • Ayala G.
      • Hayward S.W.
      • Anderson A.R.
      Investigating prostate cancer tumour–stroma interactions: clinical and biological insights from an evolutionary game.
      ,
      • Szeto G.L.
      • Finley S.D.
      Integrative approaches to cancer immunotherapy.
      ,
      • Altrock P.M.
      • Liu L.L.
      • Michor F.
      The mathematics of cancer: integrating quantitative models.
      ,
      • Brady R.
      • Enderling H.
      Mathematical models of cancer: when to predict novel therapies, and when not to.
      ], 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. [
      • Garcia V.
      • Bonhoeffer S.
      • Fu F.
      Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: a mathematical and computational examination.
      ]. The model considers a population of tumor cells, y and immune cells (effector cells x), 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 k per unit time and per capita. Tumor cells are killed with a rate m upon interaction with effector cells. Assuming a logistic growth for tumor population in the absence of immune interaction, we have [
      • Garcia V.
      • Bonhoeffer S.
      • Fu F.
      Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: a mathematical and computational examination.
      ]:
      x˙=λμ·x+k·xyy˙=ay·(1by)m·xy
      (1)


      where a is the linear growth rate (fitness) of tumor cells, and b 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 m and k values that are needed for effective immune control of tumor growth (characterized by low abundance of tumor cells or complete eradication)  [
      • Garcia V.
      • Bonhoeffer S.
      • Fu F.
      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 y0 (y1) and PD-1- (+) effector cells with x0 (x1), respectively. yi or xi for i=0,1 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, q0, PD-L1- tumor cells transform into PD-L1+ cells. Similarly, the effector cell population transforms from PD-1- into PD-1+ with a rate p0.
      • Any T-cell subtype (PD-1+/-) can in interact with any tumor cell subtype (PD-L1+/-). Their interaction strengths are modeled into a matrix {kij}. Parameter kij represents the stimulation/inhibition level of type i effector cells by type j tumor cells. For PD-1- cells i=0 and for PD-1+ ones i=1. Using similar notations, j=0,1 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 (k00>0), while PD-L1+ tumor cells suppress anti-tumor response of PD-1+ T cells (k11<0).
      Therefore, we set k00, k01 and k10>0 while k11<0 represents the immune suppression. The anti-tumor activity of the immune cell population is modeled with a matrix of killing rate, {mij}. mij is the death rate induced by a type j immune cell on a type i tumor cell (i,j={0,1}). 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, x0 and x1, and PD-L1-/+ tumor cells, y0 and y1:
      x˙0=λμx0+(k00y0+k01y1)x0p0x0,x˙1=μx1+(k10y0+k11y1)x1+p0x0,y˙0=ay0(1b(y0+y1))(m00x0+m01x1)y0q0y0,y˙1=ay1(1b(y0+y1))(m10x0+m11x1)y1+q0y0.
      (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. a is the linear growth rate of tumor and b is the inverse carrying capacity for tumor cells. PD-1- T cells transform into PD-1+ with rate p0 while PD-L1- cell become PD-L+ with rate q0. Values of p0 and q0 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. 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 kij and mij. kij (i,j=0,1) determine how the population of tumor cells affect the immune response. mij (i,j=0,1) quantify the death rate induced by effector cells on tumor cells. The interaction-matrix elements, mij and kij 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 m11 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 p0 and q0.
      We can rewrite the Eq. (2) by rescaling tumor population with carrying capacity b1 and immune population with λ.
      b·y0,1y0,1λ1·x0,1x0,1
      (3)


      This leads to rescaling of model parameters:
      b1kijkijλmijmij
      (4)


      where the subscript index i,j{0,1}. Eq. (2) can be rewritten in terms of rescaled parameters as,
      x˙0=1μx0+(k00y0+k01y1)x0p0x0x˙1=μx1+(k10y0+k11y1)x1+p0x0y˙0=a·y0(1(y0+y1))(m00x0+m01x1)y0q0y0y˙1=a·y1(1(y0+y1))(m10x0+m11x1)y1+q0y0
      (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, t, in the derivative, and write it in units of average life time of immune cells, μ1. 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, (p0, q0). 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 (a, b). We already re-scaled the model system to absorb the tumor inverse carrying capacity b into other parameters. Thus the main indicator of tumor growth and aggressiveness is the value of its linear growth rate a.
      • Tumor-immune interaction matrices ({kij},{mij}). The value of mij identifies the set of death rates values for the tumor cells due to immune killing (i,j={0,1}). While kij 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 q0 quantifies how tumor population transforms into immune-suppressive subpopulation of PD-L1+. In the meantime, at rate p0, PD-1- immune cells become PD-1+ cells with poor anti-tumor characteristic. The finite value of q0 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. 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 p0=1,2,3 while the value of q0 is kept the same as p0, q0=p0. The abundance of effector cells and tumor burden are shown as rescaled according to , and relevant model parameter values are rescaled according to . Rescaled model parameters are k00=k01=k10=1, k11=1, m00=m10=m01=0.1, m11=0, and a=0.5, μ=0.5.
      Without loss of generality, we assume the values of immune stimulation/inhibition, k00,k01,k10 are the same order of magnitude. We model suppressive effect of PD-L1+ tumor cells on immune cells by setting k11 to a negative value. Anti-tumor activity of effector cells is parametrized by prescribing values of mij. In the example shown in Fig. 3, we use m00=m01=m10=0.1 and m11=0. The value of m11 indicates that PD-1+ T-cells has no effector activity on PD-L1+ tumor cells. The tumor-immune interaction parameter, kij are set to k00=k10=k01=1, and k11=1. The positive (negative) values of kij mean stimulatory (or inhibitory) response due to immune-tumor interaction, respectively. For the rest of parameters we used a=0.5,μ=0.5 and values of p0=q0=1,2,3. In Fig. 3, we plot the total effector cell count, x0+x1, and the total tumor size, y0+y1 as a function of time, using the initial condition x1(0)=y1(0)=0. The higher the values of p0 and q0 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 (p0,q0) is shown in Fig. 4. Same model parameters as in Fig. 3 are used, and values of p0 and q0 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 p0 and q0 (see Fig. 4).
      Fig. 4
      Fig. 4Heatmap plot of tumor steady-state size as a function of immune resistance parameters, p0, and q0. We assume no treatment is applied, q1=p1=0. Tumor burden is shown as rescaled according to , and relevant model parameter values are rescaled according to . The immune inhibition/suppression parameters are k00=k01=k10=1 and k11=1, and the tumor killing rates by PD-1+/- cells are m00=m10=m01=0.1 and m11=0. And we set a=0.5 and μ=0.5 the same as in .
      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 β: k00=k01=k10=k11=β. A simple analytical expression can be found for the steady-state tumor size along the diagonal q0=p0=η:
      T0+T1=y0*+y1*=aηη2λμ+aη(aη)β
      (6)


      For small-η, the formula above can be expressed in a Taylor expansion:
      Taβ+a+a2(β1)2+4βλμ2aβ+12(aβ+a2(β1)2+4aβλμ+a)βa((β1)2a+4βλμ)η+O(η2)
      (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 T* monotonically increases with the adaptive immune resistance parameters p0=q0=η 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 p0 and q0, 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 p1. The bulk parameter p1 reflects the overall neutralization (blockade) efficacy determined by the underlying concentration-dependent binding kinetics [
      • Brown M.E.
      • Bedinger D.
      • Lilov A.
      • Rathanaswami P.
      • Vásquez M.
      • Durand S.
      • Wallace-Moyer I.
      • Zhong L.
      • Nett J.H.
      • Burnina I.
      • et al.
      Assessing the binding properties of the anti-PD-1 antibody landscape using label-free biosensors.
      ]. Similarly, we denote by q1 the reverse transition rate for tumor cells transitioning from PD-L1+ to PD-L1- compartments in the presence of anti-PD-L1 treatment.
      Fig. 5
      Fig. 5Model scheme of tumor-immune interactions in the presence of immune checkpoint blockades. Similar to , 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 q0. We assume similar transition from PD-1- to PD-1+ for effector cells with rate p0. Anti-PD-1 and anti-PD-L1 antibodies cause a reversal of this with corresponding rates p1 and q1.
      Incorporating the anti-PD-1/anti-PD-L1 treatment into the Eq. (5), we obtain the following system of differential equations:
      x˙0=1μx0+(k00y0+k01y1)x0p0x0+p1x1,x˙1=μx1+(k10y0+k11y1)x1+p0x0p1x1,y˙0=a·y0(1(y0+y1))(m00x0+m01x1)y0q0y0+q1y1,y˙1=a·y1(1(y0+y1))(m10x0+m11x1)y1+q0y0q1y1.
      (8)


      Closed-form analytical solutions can be obtained for the steady-state tumor size, T=T0+T1 (that is, y0*+y1*), if we apply the same assumption as in Eq. (6), k00=k01=k10=k11=β,q0=p0=η, as well as p1=ζ,q1=0 for anti-PD-1 treatment. We have,
      T=(ηζ)aλμ+a2ζ22a(aη+λμ2a)ζ+(aηλμ+2a)22aβ.
      (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 p1 and q1 in this paper, the parameter k11 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 k11. Similarly, we assume anti-PD-L1 therapy targets PD-L1+ tumor cells, and meanwhile can also reduce the value of tumor-immune interaction matrix k11. Hence, the parameter value of k11 during treatment is uniformly altered as a secondary effect of monotherapy or combination therapy.
      We perform bifurcation analysis of equilibrium tumor burden, T*, with respect to changes in the anti-tumor immune activity, mij. 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 k11=+1,5,10 while other immune-tumor interaction parameters are kept constant k00=k01=k10=5. Values of tumor killing rates mij are parametrized as m00=m01=m10=m,m11=0.01m. 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 k11<0 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. 6
      Fig. 6Bifurcation analysis of equilibrium tumor burden with respect to immune killing rate m under monotherapy with anti-PD-1 and anti-PD-L1. tumor fixed point, T=T0+T1, for various values of m (tumor killing rate) for PD-L1 mono-therapy as immune-tumor interaction parameter k11 is varied due to treatment. Tumor burden is shown as rescaled according to , and relevant model parameter values are rescaled according to . Model parameters are set as m00=m01=m10=m, m11=0.01m, q0=p0=5, k00=k10=k01=5, k11=10,5,1, (a) anti-PD-1: p1=1,q1=0, and (b) anti-PD-L1: p1=0,q1=1. And we set a=0.5 and μ=0.5 the same as in .
      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 q1=0 while p1>0. Similarly, for anti-PD-L1 monotherapy, q1>0 while p1=0. 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 (p0=5 and q0=1) 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 [
      • Shen X.
      • Zhao B.
      Efficacy of PD-1 or PD-L1 inhibitors and PD-L1 expression status in cancer: meta-analysis.
      ]. 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 [
      • Linhares A.D.S.
      • Battin C.
      • Jutz S.
      • Leitner J.
      • Hafner C.
      • Tobias J.
      • et al.
      Therapeutic PD-L1 antibodies are more effective than PD-1 antibodies in blocking PD-1/PD-L1 signaling.
      ].
      Fig. 7
      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 (p0=5,q0=1). Tumor burden is shown as rescaled according to , and relevant model parameter values are rescaled according to . Other model parameters used are: k00=k10=k10=10, k11=10,0,10, m00=m01=m10=0.1, m11=0, (a) anti-PD-L1: p1=0,q1=5, and (b) anti-PD-1: p1=5,q1=0. And we set a=0.5 and μ=0.5 the same as in .
      Fig. 8
      Fig. 8Residual tumor burden as a function of treatment efficacy of monotherapy, anti-PD-1 p1 versus anti-PD-L1 q1. 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 q1=p,p1=0 for anti-PD-1 we have q1=0,p1=p, where p denotes the x-axis. Depending on the expression levels of PD-1 and PD-L1 status (immune resistance parameters p0 and q0), 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, m00=m01=m10=0.1, m00=0.01, k00=k01=k10=7, k11=1, a=0.5, μ=0.5, (a) p0=5,q0=1, (b) p0=3,q0=1, (c) p0=1,q0=1.

      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 (p1,q1), 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 [
      • Loewe S.
      • Muischnek H.
      Über kombinationswirkungen.
      ].
      In Fig. 9, we show a contour plot of residual tumor burden as a function of (p1,q1). Along each contour line, the residual tumor burden is constant for the combination of values of (p1,q1). Let us denote by p1* (q1*) 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 p1 and q1 satisfying
      q1*p1+p1*q1=p1*q1*.
      (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. 9
      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 (p1,q1), which quantifies anti-PD-1 and anti-PD-L1 treatment efficacies. The Loewe additivity is indicated by the straight line p1q1+p1q1=p1q1, where q1 and p1 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 , and relevant model parameter values are rescaled according to . Model parameters are μ=0.5,m00=m01=m10=0.1,m11=0.01,k00=k01=k10=5,k11=5,a=0.5,p0=2,q0=1..
      Fig. 10
      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 
      [
      • Loewe S.
      • Muischnek H.
      Über kombinationswirkungen.
      ]
      , we restrict parameter choices of p1 and q1 for combination therapy along the line (see ), p1q1+p1q1=p1q1, where q1 and p1 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 t=50. Monotherapy: anti-PD-1 p1=2.5,q1=0, anti-PD-L1 p1=0,q1=2, combination therapy: p1=1.25,q1=1, other rescaled model parameters: m00=m01=m10=0.1, m00=0.01, k00=k01=k10=5, k11=5, p0=2, q0=1, a=0.5, μ=0.5.

      4. Discussion & conclusion

      In recent years, there has been a surge of interest in developing immune checkpoint inhibitors to treat cancer [
      • Tang J.
      • Shalabi A.
      • Hubbard-Lucey V.
      Comprehensive analysis of the clinical immuno-oncology landscape.
      ]. 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 [
      • DePillis L.
      • Eladdadi A.
      • Radunskaya A.
      Modeling cancer-immune responses to therapy.
      ,
      • de Pillis L.G.
      • Radunskaya A.E.
      • Wiseman C.L.
      A validated mathematical model of cell-mediated immune response to tumor growth.
      ,
      • Wilson S.
      • Levy D.
      A mathematical model of the enhancement of tumor vaccine efficacy by immunotherapy.
      ,
      • Yamamoto Y.
      • Offord C.P.
      • Kimura G.
      • Kuribayashi S.
      • Takeda H.
      • Tsuchiya S.
      • et al.
      Tumour and immune cell dynamics explain the PSA bounce after prostate cancer brachytherapy.
      ,
      • Serre R.
      • Benzekry S.
      • Padovani L.
      • Meille C.
      • André N.
      • Ciccolini J.
      • et al.
      Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy.
      ,
      • Castiglione F.
      • Piccoli B.
      Cancer immunotherapy, mathematical modeling and optimal control.
      ,
      • Owens K.L.
      • Bozic I.
      Modelling CAR T-cell therapy with patient preconditioning.
      ,
      • Kimmel G.J.
      • Locke F.L.
      • Altrock P.M.
      Response to car T cell therapy can be explained by ecological cell dynamics and stochastic extinction events.
      ]. Prior mathematical models of cancer-immune interactions [
      • Kirschner D.
      • Panetta J.C.
      Modeling immunotherapy of the tumor–immune interaction.
      ,
      • Eftimie R.
      • Bramson J.L.
      • Earn D.J.
      Interactions between the immune system and cancer: a brief review of non-spatial mathematical models.
      ] are mostly based on ordinary differential equations in combination with stochastic in-silico simulations [
      • Lakatos E.
      • Williams M.J.
      • Schenck R.O.
      • Cross W.C.
      • Househam J.
      • Zapata L.
      • et al.
      Evolutionary dynamics of neoantigens in growing tumors.
      ,
      • West J.
      • Robertson-Tessi M.
      • Luddy K.
      • Park D.S.
      • Williamson D.F.
      • Harmon C.
      • et al.
      The immune checkpoint kick start: optimization of neoadjuvant combination therapy using game theory.
      ]. 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 [
      • Linhares A.D.S.
      • Battin C.
      • Jutz S.
      • Leitner J.
      • Hafner C.
      • Tobias J.
      • et al.
      Therapeutic PD-L1 antibodies are more effective than PD-1 antibodies in blocking PD-1/PD-L1 signaling.
      ,
      • Chen Y.-J.
      • Huang W.-C.
      • Liu S.-Y.
      • Ko C.C.
      Sequential blockade of PD-1 and PD-L1 causes fulminant cardiotoxicity: from case report to mice model validation.
      ]. 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 [
      • Shen X.
      • Zhao B.
      Efficacy of PD-1 or PD-L1 inhibitors and PD-L1 expression status in cancer: meta-analysis.
      ].
      Randomized trials of monotherapy and combination therapy for treating advanced cancers are underway and some have been completed [
      • Sui H.
      • Ma N.
      • Wang Y.
      • Li H.
      • Liu X.
      • Su Y.
      • et al.
      Anti-PD-1/PD-L1 therapy for non-small-cell lung cancer: toward personalized medicine and combination strategies.
      ,
      • Sato H.
      • Okonogi N.
      • Nakano T.
      Rationale of combination of anti-PD-1/PD-L1 antibody therapy and radiotherapy for cancer treatment.
      ]. 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 [
      • Page D.B.
      • Bear H.
      • Prabhakaran S.
      • Gatti-Mays M.E.
      • Thomas A.
      • Cobain E.
      • McArthur H.
      • Balko J.M.
      • Gameiro S.R.
      • Nanda R.
      • et al.
      Two may be better than one: PD-1/PD-L1 blockade combination approaches in metastatic breast cancer.
      ]. Future work can incorporate into the current models with cancer and immune biomarkers measured before and during immunotherapy [

      Willis J.C., Lord G.M.. Immune biomarkers: the promises and pitfalls of personalized medicine. Nat Rev Immunol2015. 15, 5, 323–329,

      ], and develop theory-informed combination therapies to overcome immunotherapy resistance and improve response.
      In this work, the switching rate of effector cells, p0 (and similarly tumor cells, q0) 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 [
      • Wang H.
      • Milberg O.
      • Bartelink I.H.
      • Vicini P.
      • Wang B.
      • Narwal R.
      • et al.
      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.
      ,
      • Milberg O.
      • Gong C.
      • Jafarnejad M.
      • Bartelink I.H.
      • Wang B.
      • Vicini P.
      • et al.
      A QSP model for predicting clinical responses to monotherapy, combination and sequential therapy following CTLA-4, PD-1, and PD-L1 checkpoint blockade.
      ] and also spatial infiltration of T cells into tumor mass [
      • Gong C.
      • Milberg O.
      • Wang B.
      • Vicini P.
      • Narwal R.
      • Roskos L.
      • et al.
      A computational multiscale agent-based model for simulating spatio-temporal tumour immune response to PD1 and PDL1 inhibition.
      ]. 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 [
      • Sharma P.
      • Allison J.P.
      The future of immune checkpoint therapy.
      ]. Due to the intricate nature of cancer-immune interactions [
      • Szeto G.L.
      • Finley S.D.
      Integrative approaches to 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 [
      • Altrock P.M.
      • Liu L.L.
      • Michor F.
      The mathematics of cancer: integrating quantitative models.
      ]: 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 [
      • Brady R.
      • Enderling H.
      Mathematical models of cancer: when to predict novel therapies, and when not to.
      ]. 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 [
      • Brady R.
      • Enderling H.
      Mathematical models of cancer: when to predict novel therapies, and when not to.
      ], our models can be used to test hypothetical drug administration schedules in-silico and to optimize cancer treatment in a personalized fashion.

      CRediT authorship contribution statement

      Kamran Kaveh: Conceptualization, Formal analysis, Funding acquisition, Supervision, Writing – original draft, Writing – review & editing. Feng Fu: Conceptualization, Formal analysis, Funding acquisition, Writing – original draft, Writing – review & editing.

      Acknowledgments

      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.

      References

        • Ribas A.
        • et al.
        Releasing the brakes on cancer immunotherapy.
        N Engl J Med. 2015; 373: 1490-1492
        • Ribas A.
        • Wolchok J.D.
        Cancer immunotherapy using checkpoint blockade.
        Science. 2018; 359: 1350-1355
        • Sharma P.
        • Allison J.P.
        Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential.
        Cell. 2015; 161: 205-214
        • Kucerova P.
        • Cervinkova M.
        Spontaneous regression of tumour and the role of microbial infection–possibilities for cancer treatment.
        Anti-Cancer Drugs. 2016; 27: 269
        • Ribas A.
        Adaptive immune resistance: how cancer protects from immune attack.
        Cancer Discov. 2015; 5: 915-919
        • Kalbasi A.
        • Ribas A.
        Tumour-intrinsic resistance to immune checkpoint blockade.
        Nat Rev Immunol. 2019; : 1-15
        • Wei S.C.
        • Duffy C.R.
        • Allison J.P.
        Fundamental mechanisms of immune checkpoint blockade therapy.
        Cancer Discov. 2018; 8: 1069-1086
        • Littman D.R.
        Releasing the brakes on cancer immunotherapy.
        Cell. 2015; 162: 1186-1190
        • Tang J.
        • Shalabi A.
        • Hubbard-Lucey V.
        Comprehensive analysis of the clinical immuno-oncology landscape.
        Ann Oncol. 2018; 29: 84-91
        • Sharma P.
        • Allison J.P.
        The future of immune checkpoint therapy.
        Science. 2015; 348: 56-61
        • Wei S.C.
        • Levine J.H.
        • Cogdill A.P.
        • Zhao Y.
        • Anang N.-A.A.
        • Andrews M.C.
        • Sharma P.
        • Wang J.
        • Wargo J.A.
        • Pe?er D.
        Distinct cellular mechanisms underlie anti-CTLA-4 and anti-PD-1 checkpoint blockade.
        Cell. 2017; 170: 1120-1133
        • Wei S.C.
        • Anang N.-A.A.
        • Sharma R.
        • Andrews M.C.
        • Reuben A.
        • Levine J.H.
        • Cogdill A.P.
        • Mancuso J.J.
        • Wargo J.A.
        • Pe?er D.
        • et al.
        Combination anti–CTLA-4 plus anti–PD-1 checkpoint blockade utilizes cellular mechanisms partially distinct from monotherapies.
        Proceedings of the national academy of sciences. vol. 116. 2019: 22699-22709
      1. Willis J.C., Lord G.M.. Immune biomarkers: the promises and pitfalls of personalized medicine. Nat Rev Immunol2015. 15, 5, 323–329,

        • Chaplin D.D.
        Overview of the immune response.
        J Allergy Clin Immunol. 2010; 125: S3-S23
        • Nicholson L.B.
        The immune system.
        Essays Biochem. 2016; 60: 275-301
        • Johnson R.M.G.
        • Dong H.
        Functional expression of programmed death-ligand 1 (B7-H1) by immune cells and tumor cells.
        Front Immunol. 2017; 8: 961
        • Wu Y.
        • Chen W.
        • Xu Z.P.G.
        • Gu W.
        PD-L1 distribution and perspective for cancer immunotherapy–blockade, knockdown, or inhibition.
        Front Immunol. 2019; 10: 2022
        • Parsa A.T.
        • Waldron J.S.
        • Panner A.
        • Crane C.A.
        • Parney I.F.
        • Barry J.J.
        • et al.
        Loss of tumor suppressor PTEN function increases B7-H1 expression and immunoresistance in glioma.
        Nat Med. 2007; 13: 84
        • Akbay E.A.
        • Koyama S.
        • Carretero J.
        • Altabef A.
        • Tchaicha J.H.
        • Christensen C.L.
        • et al.
        Activation of the PD-1 pathway contributes to immune escape in EGFR-driven lung tumors.
        Cancer Discov. 2013; 3: 1355-1363
        • Atefi M.
        • Avramis E.
        • Lassen A.
        • Wong D.J.
        • Robert L.
        • Foulad D.
        • et al.
        Effects of MAPK and PI3K pathways on PD-L1 expression in melanoma.
        Clin Cancer Res. 2014; 20: 3446-3457
        • Taube J.M.
        • Anders R.A.
        • Young G.D.
        • Xu H.
        • Sharma R.
        • McMiller T.L.
        • et al.
        Colocalization of inflammatory response with B7-H1 expression in human melanocytic lesions supports an adaptive resistance mechanism of immune escape.
        Sci Transl Med. 2012; 4: 127ra37
        • Tumeh P.C.
        • Harview C.L.
        • Yearley J.H.
        • Shintaku I.P.
        • Taylor E.J.
        • Robert L.
        • Chmielowski B.
        • Spasic M.
        • Henry G.
        • Ciobanu V.
        • et al.
        PD-1 blockade induces responses by inhibiting adaptive immune resistance.
        Nature. 2014; 515: 568
        • Pentcheva-Hoang T.
        • Simpson T.R.
        • Montalvo-Ortiz W.
        • Allison J.P.
        Cytotoxic T lymphocyte antigen-4 blockade enhances antitumor immunity by stimulating melanoma-specific T-cell motility.
        Cancer Immunol Res. 2014; 2: 970-980
        • Linhares A.D.S.
        • Battin C.
        • Jutz S.
        • Leitner J.
        • Hafner C.
        • Tobias J.
        • et al.
        Therapeutic PD-L1 antibodies are more effective than PD-1 antibodies in blocking PD-1/PD-L1 signaling.
        Sci Rep. 2019; 9: 1-9
        • Chen Y.-J.
        • Huang W.-C.
        • Liu S.-Y.
        • Ko C.C.
        Sequential blockade of PD-1 and PD-L1 causes fulminant cardiotoxicity: from case report to mice model validation.
        Ann Oncol. 2018; 29
      2. Viii431
        • Shen X.
        • Zhao B.
        Efficacy of PD-1 or PD-L1 inhibitors and PD-L1 expression status in cancer: meta-analysis.
        Bmj. 2018; 362
        • Basanta D.
        • Scott J.G.
        • Fishman M.N.
        • Ayala G.
        • Hayward S.W.
        • Anderson A.R.
        Investigating prostate cancer tumour–stroma interactions: clinical and biological insights from an evolutionary game.
        Br J Cancer. 2012; 106: 174-181
        • Szeto G.L.
        • Finley S.D.
        Integrative approaches to cancer immunotherapy.
        Trends Cancer. 2019; 5: 400-410
        • Altrock P.M.
        • Liu L.L.
        • Michor F.
        The mathematics of cancer: integrating quantitative models.
        Nat Rev Cancer. 2015; 15: 730-745
        • Brady R.
        • Enderling H.
        Mathematical models of cancer: when to predict novel therapies, and when not to.
        Bull Math Biol. 2019; 81: 3722-3731
        • Garcia V.
        • Bonhoeffer S.
        • Fu F.
        Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: a mathematical and computational examination.
        J Theor Biol. 2020; 492: 110185
        • Brown M.E.
        • Bedinger D.
        • Lilov A.
        • Rathanaswami P.
        • Vásquez M.
        • Durand S.
        • Wallace-Moyer I.
        • Zhong L.
        • Nett J.H.
        • Burnina I.
        • et al.
        Assessing the binding properties of the anti-PD-1 antibody landscape using label-free biosensors.
        PLoS ONE. 2020; 15: e0229206
        • Loewe S.
        • Muischnek H.
        Über kombinationswirkungen.
        Naunyn-Schmiedebergs Archiv für Experimentelle Pathologie und Pharmakologie. 1926; 114: 313-326
        • DePillis L.
        • Eladdadi A.
        • Radunskaya A.
        Modeling cancer-immune responses to therapy.
        J Pharmacokinet Pharmacodyn. 2014; 41: 461-478
        • de Pillis L.G.
        • Radunskaya A.E.
        • Wiseman C.L.
        A validated mathematical model of cell-mediated immune response to tumor growth.
        Cancer Res. 2005; 65: 7950-7958
        • Wilson S.
        • Levy D.
        A mathematical model of the enhancement of tumor vaccine efficacy by immunotherapy.
        Bull Math Biol. 2012; 74: 1485-1500
        • Yamamoto Y.
        • Offord C.P.
        • Kimura G.
        • Kuribayashi S.
        • Takeda H.
        • Tsuchiya S.
        • et al.
        Tumour and immune cell dynamics explain the PSA bounce after prostate cancer brachytherapy.
        Br J Cancer. 2016; 115: 195-202
        • Serre R.
        • Benzekry S.
        • Padovani L.
        • Meille C.
        • André N.
        • Ciccolini J.
        • et al.
        Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy.
        Cancer Res. 2016; 76: 4931-4940
        • Castiglione F.
        • Piccoli B.
        Cancer immunotherapy, mathematical modeling and optimal control.
        J Theor Biol. 2007; 247: 723-732
        • Owens K.L.
        • Bozic I.
        Modelling CAR T-cell therapy with patient preconditioning.
        bioRxiv. 2020;
        • Kimmel G.J.
        • Locke F.L.
        • Altrock P.M.
        Response to car T cell therapy can be explained by ecological cell dynamics and stochastic extinction events.
        bioRxiv. 2020; : 717074
        • Kirschner D.
        • Panetta J.C.
        Modeling immunotherapy of the tumor–immune interaction.
        J Math Biol. 1998; 37: 235-252
        • Eftimie R.
        • Bramson J.L.
        • Earn D.J.
        Interactions between the immune system and cancer: a brief review of non-spatial mathematical models.
        Bull Math Biol. 2011; 73: 2-32
        • Lakatos E.
        • Williams M.J.
        • Schenck R.O.
        • Cross W.C.
        • Househam J.
        • Zapata L.
        • et al.
        Evolutionary dynamics of neoantigens in growing tumors.
        Nat Genet. 2020; 52: 1057-1066
        • West J.
        • Robertson-Tessi M.
        • Luddy K.
        • Park D.S.
        • Williamson D.F.
        • Harmon C.
        • et al.
        The immune checkpoint kick start: optimization of neoadjuvant combination therapy using game theory.
        JCO Clin Cancer Inf. 2019; 3: 1-12
        • Sui H.
        • Ma N.
        • Wang Y.
        • Li H.
        • Liu X.
        • Su Y.
        • et al.
        Anti-PD-1/PD-L1 therapy for non-small-cell lung cancer: toward personalized medicine and combination strategies.
        J Immunol Res. 2018; 2018
        • Sato H.
        • Okonogi N.
        • Nakano T.
        Rationale of combination of anti-PD-1/PD-L1 antibody therapy and radiotherapy for cancer treatment.
        Int J Clin Oncol. 2020; 25: 801-809
        • Page D.B.
        • Bear H.
        • Prabhakaran S.
        • Gatti-Mays M.E.
        • Thomas A.
        • Cobain E.
        • McArthur H.
        • Balko J.M.
        • Gameiro S.R.
        • Nanda R.
        • et al.
        Two may be better than one: PD-1/PD-L1 blockade combination approaches in metastatic breast cancer.
        NPJ Breast Cancer. 2019; 5: 1-9
        • Wang H.
        • Milberg O.
        • Bartelink I.H.
        • Vicini P.
        • Wang B.
        • Narwal R.
        • et al.
        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.
        R Soc Open Sci. 2019; 6: 190366
        • Milberg O.
        • Gong C.
        • Jafarnejad M.
        • Bartelink I.H.
        • Wang B.
        • Vicini P.
        • et al.
        A QSP model for predicting clinical responses to monotherapy, combination and sequential therapy following CTLA-4, PD-1, and PD-L1 checkpoint blockade.
        Sci Rep. 2019; 9: 1-17
        • Gong C.
        • Milberg O.
        • Wang B.
        • Vicini P.
        • Narwal R.
        • Roskos L.
        • et al.
        A computational multiscale agent-based model for simulating spatio-temporal tumour immune response to PD1 and PDL1 inhibition.
        J R Soc Interface. 2017; 14: 20170320