## 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 [

1

, 2

, 3

]. 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 [[4]

]. 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 [5

, 6

]. 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 [

2

, 3

, 7

, 8

]. 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 [

7

, 9

].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 [

[10]

].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 [

11

, 12

]. 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 [- 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

[13]

]. 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. [

7

, 14

].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 [

[15]

].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 [

[16]

]. 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 [[17]

]. 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 [18

, 19

, 20

]. However, it seems that adaptive interferon-inducible expression of PD-L1 is more common than the constitute expression in many cancers [21

, 22

]. 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 [

[2]

]. 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 [[23]

].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 [

24

, , 26

]. 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 [27

, 28

, 29

, 30

], 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. [

where $a$ is the linear growth rate (fitness) of tumor cells, and $b$ is the inverse carrying capacity. $\mu $ 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) [

[31]

]. 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 ($\lambda $). 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 [[31]

]:$\begin{array}{cc}\hfill \dot{x}& =\lambda -\mu \xb7x+k\xb7xy\hfill \\ \hfill \dot{y}& =ay\xb7(1-by)-m\xb7xy\hfill \end{array}$

(1)

where $a$ is the linear growth rate (fitness) of tumor cells, and $b$ is the inverse carrying capacity. $\mu $ 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) [

[31]

].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 ${y}_{0}$ (${y}_{1}$) and PD-1- (+) effector cells with ${x}_{0}$ (${x}_{1}$), respectively. ${y}_{i}$ or ${x}_{i}$ 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 $\lambda $.
- •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, ${q}_{0}$, PD-L1- tumor cells transform into PD-L1+ cells. Similarly, the effector cell population transforms from PD-1- into PD-1+ with a rate ${p}_{0}$.
- •Any T-cell subtype (PD-1+/-) can in interact with any tumor cell subtype (PD-L1+/-). Their interaction strengths are modeled into a matrix $\left\{{k}_{ij}\right\}$. Parameter ${k}_{ij}$ 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 (${k}_{00}>0$), while PD-L1+ tumor cells suppress anti-tumor response of PD-1+ T cells (${k}_{11}<0$).

Therefore, we set ${k}_{00}$, ${k}_{01}$ and ${k}_{10}>0$ while ${k}_{11}<0$ represents the immune suppression. The anti-tumor activity of the immune cell population is modeled with a matrix of killing rate, $\left\{{m}_{ij}\right\}$. ${m}_{ij}$ 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, ${x}_{0}$ and ${x}_{1}$, and PD-L1-/+ tumor cells, ${y}_{0}$ and ${y}_{1}$:

Here, $\lambda $ is the recruiting rate of the effector cells from the lymph nodes to the local tumor microenvironment. $\mu $ 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 ${p}_{0}$ while PD-L1- cell become PD-L+ with rate ${q}_{0}$. Values of ${p}_{0}$ and ${q}_{0}$ together quantify how fast the immune resistance is developed against an anti-tumor immune response. The above mechanisms are schematically presented in Fig. 2.

$\begin{array}{cc}{\dot{x}}_{0}\hfill & =\lambda -\mu {x}_{0}+({k}_{00}{y}_{0}+{k}_{01}{y}_{1}){x}_{0}-{p}_{0}{x}_{0},\hfill \\ {\dot{x}}_{1}\hfill & =-\mu {x}_{1}+({k}_{10}{y}_{0}+{k}_{11}{y}_{1}){x}_{1}+{p}_{0}{x}_{0},\hfill \\ {\dot{y}}_{0}\hfill & =a{y}_{0}(1-b({y}_{0}+{y}_{1}))-({m}_{00}{x}_{0}+{m}_{01}{x}_{1}){y}_{0}-{q}_{0}{y}_{0},\hfill \\ {\dot{y}}_{1}\hfill & =a{y}_{1}(1-b({y}_{0}+{y}_{1}))-({m}_{10}{x}_{0}+{m}_{11}{x}_{1}){y}_{1}+{q}_{0}{y}_{0}.\hfill \end{array}$

(2)

Here, $\lambda $ is the recruiting rate of the effector cells from the lymph nodes to the local tumor microenvironment. $\mu $ 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 ${p}_{0}$ while PD-L1- cell become PD-L+ with rate ${q}_{0}$. Values of ${p}_{0}$ and ${q}_{0}$ together quantify how fast the immune resistance is developed against an anti-tumor immune response. The above mechanisms are schematically presented in Fig. 2.

We can rewrite the Eq. (2) by rescaling tumor population with carrying capacity ${b}^{-1}$ and immune population with $\lambda $.

$\begin{array}{cc}\hfill b\xb7{y}_{0,1}& \to {y}_{0,1}\hfill \\ \hfill {\lambda}^{-1}\xb7{x}_{0,1}& \to {x}_{0,1}\hfill \end{array}$

(3)

This leads to rescaling of model parameters:

where the subscript index $i,j\in \{0,1\}$. Eq. (2) can be rewritten in terms of rescaled parameters as,

$\begin{array}{c}{b}^{-1}{k}_{\mathit{ij}}\to {k}_{\mathit{ij}}\\ \lambda {m}_{\mathit{ij}}\to {m}_{\mathit{ij}}\end{array}$

(4)

where the subscript index $i,j\in \{0,1\}$. Eq. (2) can be rewritten in terms of rescaled parameters as,

$\begin{array}{cc}{\dot{x}}_{0}\hfill & =1-\mu {x}_{0}+({k}_{00}{y}_{0}+{k}_{01}{y}_{1}){x}_{0}-{p}_{0}{x}_{0}\hfill \\ {\dot{x}}_{1}\hfill & =\phantom{\rule{0.33em}{0ex}}-\mu {x}_{1}+({k}_{10}{y}_{0}+{k}_{11}{y}_{1}){x}_{1}+{p}_{0}{x}_{0}\hfill \\ {\dot{y}}_{0}\hfill & =a\xb7{y}_{0}(1-({y}_{0}+{y}_{1}\left)\right)-({m}_{00}{x}_{0}+{m}_{01}{x}_{1}){y}_{0}-{q}_{0}{y}_{0}\hfill \\ {\dot{y}}_{1}\hfill & =a\xb7{y}_{1}(1-({y}_{0}+{y}_{1}\left)\right)-({m}_{10}{x}_{0}+{m}_{11}{x}_{1}){y}_{1}+{q}_{0}{y}_{0}\hfill \end{array}$

(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, ${\mu}^{-1}$. This would reduce the model parameters to one less. However, we keep variable $\mu $ 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, (**${p}_{0}$, ${q}_{0}$**)**. These basically quantify how impaired the immune response is due to $\gamma $-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 (**$\left\{{k}_{ij}\right\},\left\{{m}_{ij}\right\}$**)**. The value of ${m}_{ij}$ identifies the set of death rates values for the tumor cells due to immune killing ($i,j=\{0,1\}$). While ${k}_{ij}$ 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 ${q}_{0}$ quantifies how tumor population transforms into immune-suppressive subpopulation of PD-L1+. In the meantime, at rate ${p}_{0}$, PD-1- immune cells become PD-1+ cells with poor anti-tumor characteristic. The finite value of ${q}_{0}$ 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).

Without loss of generality, we assume the values of immune stimulation/inhibition, ${k}_{00},{k}_{01},{k}_{10}$ are the same order of magnitude. We model suppressive effect of PD-L1+ tumor cells on immune cells by setting ${k}_{11}$ to a negative value. Anti-tumor activity of effector cells is parametrized by prescribing values of ${m}_{ij}$. In the example shown in Fig. 3, we use ${m}_{00}={m}_{01}={m}_{10}=0.1$ and ${m}_{11}=0$. The value of ${m}_{11}$ indicates that PD-1+ T-cells has no effector activity on PD-L1+ tumor cells. The tumor-immune interaction parameter, ${k}_{ij}$ are set to ${k}_{00}={k}_{10}={k}_{01}=1$, and ${k}_{11}=-1$. The positive (negative) values of ${k}_{ij}$ mean stimulatory (or inhibitory) response due to immune-tumor interaction, respectively. For the rest of parameters we used $a=0.5,\mu =0.5$ and values of ${p}_{0}={q}_{0}=1,2,3$. In Fig. 3, we plot the total effector cell count, ${x}_{0}+{x}_{1}$, and the total tumor size, ${y}_{0}+{y}_{1}$ as a function of time, using the initial condition ${x}_{1}\left(0\right)={y}_{1}\left(0\right)=0$. The higher the values of ${p}_{0}$ and ${q}_{0}$ 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 $({p}_{0},{q}_{0})$ is shown in Fig. 4. Same model parameters as in Fig. 3 are used, and values of ${p}_{0}$ and ${q}_{0}$ are varied independently. As adaptive immune resistance is modulated by the $\gamma $-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 ${p}_{0}$ and ${q}_{0}$ (see Fig. 4).

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 $\beta $: ${k}_{00}={k}_{01}={k}_{10}=-{k}_{11}=\beta $. A simple analytical expression can be found for the steady-state tumor size along the diagonal ${q}_{0}={p}_{0}=\eta $:

${T}_{0}^{\u2605}+{T}_{1}^{\u2605}={y}_{0}^{*}+{y}_{1}^{*}=\frac{a\eta -{\eta}^{2}-\lambda \mu +a-\eta}{(a-\eta )\beta}$

(6)

For small-$\eta $, the formula above can be expressed in a Taylor expansion:

$\begin{array}{ccc}\hfill {T}^{\u2605}& \approx & \frac{a\beta +a+\sqrt{{a}^{2}{(\beta -1)}^{2}+4\beta \lambda \mu}}{2a\beta}\hfill \\ & & +\phantom{\rule{0.16em}{0ex}}\frac{1}{2}\frac{(-a\beta +\sqrt{{a}^{2}{(\beta -1)}^{2}+4a\beta \lambda \mu}+a)}{\beta \sqrt{a({(\beta -1)}^{2}a+4\beta \lambda \mu )}}\eta +\mathcal{O}\left({\eta}^{2}\right)\hfill \end{array}$

(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 ${p}_{0}={q}_{0}=\eta $ as the coefficient of the first order expansion in $\eta $ 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 ${p}_{0}$ and ${q}_{0}$, 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 ${p}_{1}$. The bulk parameter ${p}_{1}$ reflects the overall neutralization (blockade) efficacy determined by the underlying concentration-dependent binding kinetics [

[32]

]. Similarly, we denote by ${q}_{1}$ the reverse transition rate for tumor cells transitioning from PD-L1+ to PD-L1- compartments in the presence of anti-PD-L1 treatment.Incorporating the anti-PD-1/anti-PD-L1 treatment into the Eq. (5), we obtain the following system of differential equations:

$\begin{array}{cc}{\dot{x}}_{0}\hfill & =1-\mu {x}_{0}+({k}_{00}{y}_{0}+{k}_{01}{y}_{1}){x}_{0}-{p}_{0}{x}_{0}+{p}_{1}{x}_{1},\hfill \\ {\dot{x}}_{1}\hfill & =\phantom{\rule{0.33em}{0ex}}-\mu {x}_{1}+({k}_{10}{y}_{0}+{k}_{11}{y}_{1}){x}_{1}+{p}_{0}{x}_{0}-{p}_{1}{x}_{1},\hfill \\ {\dot{y}}_{0}\hfill & =a\xb7{y}_{0}(1-({y}_{0}+{y}_{1}\left)\right)-({m}_{00}{x}_{0}+{m}_{01}{x}_{1}){y}_{0}-{q}_{0}{y}_{0}+{q}_{1}{y}_{1},\hfill \\ {\dot{y}}_{1}\hfill & =a\xb7{y}_{1}(1-({y}_{0}+{y}_{1}\left)\right)-({m}_{10}{x}_{0}+{m}_{11}{x}_{1}){y}_{1}+{q}_{0}{y}_{0}-{q}_{1}{y}_{1}.\hfill \end{array}$

(8)

Closed-form analytical solutions can be obtained for the steady-state tumor size, ${T}^{\u2605}={T}_{0}^{\u2605}+{T}_{1}^{\u2605}$ (that is, ${y}_{0}^{*}+{y}_{1}^{*}$), if we apply the same assumption as in Eq. (6), ${k}_{00}={k}_{01}={k}_{10}=-{k}_{11}=\beta ,{q}_{0}={p}_{0}=\eta $, as well as ${p}_{1}=\zeta ,{q}_{1}=0$ for anti-PD-1 treatment. We have,

${T}^{\u2605}=\frac{(\eta -\zeta )a-\lambda \mu +\sqrt{{a}^{2}{\zeta}^{2}-2a(a\eta +\lambda \mu -2a)\zeta +{(a\eta -\lambda \mu +2a)}^{2}}}{2a\beta}.$

(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 ${p}_{1}$ and ${q}_{1}$ in this paper, the parameter ${k}_{11}$ 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 ${k}_{11}$. Similarly, we assume anti-PD-L1 therapy targets PD-L1+ tumor cells, and meanwhile can also reduce the value of tumor-immune interaction matrix ${k}_{11}$. Hence, the parameter value of ${k}_{11}$ 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, ${m}_{ij}$. 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 ${k}_{11}=+1,-5,-10$ while other immune-tumor interaction parameters are kept constant ${k}_{00}={k}_{01}={k}_{10}=5$. Values of tumor killing rates ${m}_{ij}$ are parametrized as ${m}_{00}={m}_{01}={m}_{10}=m,{m}_{11}=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 ${k}_{11}<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.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 ${q}_{1}=0$ while ${p}_{1}>0$. Similarly, for anti-PD-L1 monotherapy, ${q}_{1}>0$ while ${p}_{1}=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 (${p}_{0}=5$ and ${q}_{0}=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 [

[26]

]. 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 [[24]

].### 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 $({p}_{1},{q}_{1})$, 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 [

[33]

].In Fig. 9, we show a contour plot of residual tumor burden as a function of $({p}_{1},{q}_{1})$. Along each contour line, the residual tumor burden is constant for the combination of values of $({p}_{1},{q}_{1})$. Let us denote by ${p}_{1}^{*}$ (${q}_{1}^{*}$) 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 ${p}_{1}$ and ${q}_{1}$ satisfying

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.

${q}_{1}^{*}{p}_{1}+{p}_{1}^{*}{q}_{1}={p}_{1}^{*}{q}_{1}^{*}.$

(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.

## 4. Discussion & conclusion

In recent years, there has been a surge of interest in developing immune checkpoint inhibitors to treat cancer [

[9]

]. 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 [34

, 35

, 36

, 37

, 38

, 39

, 40

, 41

]. Prior mathematical models of cancer-immune interactions [42

, 43

] are mostly based on ordinary differential equations in combination with stochastic *in-silico*simulations [44

, 45

]. 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 [

24

, ]. 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 [[26]

].Randomized trials of monotherapy and combination therapy for treating advanced cancers are underway and some have been completed [

46

, 47

]. 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 [[48]

]. Future work can incorporate into the current models with cancer and immune biomarkers measured before and during immunotherapy [[13]

], and develop theory-informed combination therapies to overcome immunotherapy resistance and improve response.In this work, the switching rate of effector cells, ${p}_{0}$ (and similarly tumor cells, ${q}_{0}$) 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 [

49

, 50

] and also spatial infiltration of T cells into tumor mass [[51]

]. 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 [

[10]

]. Due to the intricate nature of cancer-immune interactions [[28]

], 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 [[29]

]: 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 [[30]

]. 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 [

[30]

], 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

- Releasing the brakes on cancer immunotherapy.
*N Engl J Med.*2015; 373: 1490-1492 - Cancer immunotherapy using checkpoint blockade.
*Science.*2018; 359: 1350-1355 - Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential.
*Cell.*2015; 161: 205-214 - Spontaneous regression of tumour and the role of microbial infection–possibilities for cancer treatment.
*Anti-Cancer Drugs.*2016; 27: 269 - Adaptive immune resistance: how cancer protects from immune attack.
*Cancer Discov.*2015; 5: 915-919 - Tumour-intrinsic resistance to immune checkpoint blockade.
*Nat Rev Immunol.*2019; : 1-15 - Fundamental mechanisms of immune checkpoint blockade therapy.
*Cancer Discov.*2018; 8: 1069-1086 - Releasing the brakes on cancer immunotherapy.
*Cell.*2015; 162: 1186-1190 - Comprehensive analysis of the clinical immuno-oncology landscape.
*Ann Oncol.*2018; 29: 84-91 - The future of immune checkpoint therapy.
*Science.*2015; 348: 56-61 - Distinct cellular mechanisms underlie anti-CTLA-4 and anti-PD-1 checkpoint blockade.
*Cell.*2017; 170: 1120-1133 - 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
Willis J.C., Lord G.M.. Immune biomarkers: the promises and pitfalls of personalized medicine. Nat Rev Immunol2015. 15, 5, 323–329,

- Overview of the immune response.
*J Allergy Clin Immunol.*2010; 125: S3-S23 - The immune system.
*Essays Biochem.*2016; 60: 275-301 - Functional expression of programmed death-ligand 1 (B7-H1) by immune cells and tumor cells.
*Front Immunol.*2017; 8: 961 - PD-L1 distribution and perspective for cancer immunotherapy–blockade, knockdown, or inhibition.
*Front Immunol.*2019; 10: 2022 - Loss of tumor suppressor PTEN function increases B7-H1 expression and immunoresistance in glioma.
*Nat Med.*2007; 13: 84 - Activation of the PD-1 pathway contributes to immune escape in EGFR-driven lung tumors.
*Cancer Discov.*2013; 3: 1355-1363 - Effects of MAPK and PI3K pathways on PD-L1 expression in melanoma.
*Clin Cancer Res.*2014; 20: 3446-3457 - 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 - PD-1 blockade induces responses by inhibiting adaptive immune resistance.
*Nature.*2014; 515: 568 - Cytotoxic T lymphocyte antigen-4 blockade enhances antitumor immunity by stimulating melanoma-specific T-cell motility.
*Cancer Immunol Res.*2014; 2: 970-980 - Therapeutic PD-L1 antibodies are more effective than PD-1 antibodies in blocking PD-1/PD-L1 signaling.
*Sci Rep.*2019; 9: 1-9 - Sequential blockade of PD-1 and PD-L1 causes fulminant cardiotoxicity: from case report to mice model validation.
*Ann Oncol.*2018; 29 - Viii431
- Efficacy of PD-1 or PD-L1 inhibitors and PD-L1 expression status in cancer: meta-analysis.
*Bmj.*2018; 362 - Investigating prostate cancer tumour–stroma interactions: clinical and biological insights from an evolutionary game.
*Br J Cancer.*2012; 106: 174-181 - Integrative approaches to cancer immunotherapy.
*Trends Cancer.*2019; 5: 400-410 - The mathematics of cancer: integrating quantitative models.
*Nat Rev Cancer.*2015; 15: 730-745 - Mathematical models of cancer: when to predict novel therapies, and when not to.
*Bull Math Biol.*2019; 81: 3722-3731 - Cancer-induced immunosuppression can enable effectiveness of immunotherapy through bistability generation: a mathematical and computational examination.
*J Theor Biol.*2020; 492: 110185 - Assessing the binding properties of the anti-PD-1 antibody landscape using label-free biosensors.
*PLoS ONE.*2020; 15: e0229206 - Über kombinationswirkungen.
*Naunyn-Schmiedebergs Archiv für Experimentelle Pathologie und Pharmakologie.*1926; 114: 313-326 - Modeling cancer-immune responses to therapy.
*J Pharmacokinet Pharmacodyn.*2014; 41: 461-478 - A validated mathematical model of cell-mediated immune response to tumor growth.
*Cancer Res.*2005; 65: 7950-7958 - A mathematical model of the enhancement of tumor vaccine efficacy by immunotherapy.
*Bull Math Biol.*2012; 74: 1485-1500 - Tumour and immune cell dynamics explain the PSA bounce after prostate cancer brachytherapy.
*Br J Cancer.*2016; 115: 195-202 - Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy.
*Cancer Res.*2016; 76: 4931-4940 - Cancer immunotherapy, mathematical modeling and optimal control.
*J Theor Biol.*2007; 247: 723-732 - Modelling CAR T-cell therapy with patient preconditioning.
*bioRxiv.*2020; - Response to car T cell therapy can be explained by ecological cell dynamics and stochastic extinction events.
*bioRxiv.*2020; : 717074 - Modeling immunotherapy of the tumor–immune interaction.
*J Math Biol.*1998; 37: 235-252 - Interactions between the immune system and cancer: a brief review of non-spatial mathematical models.
*Bull Math Biol.*2011; 73: 2-32 - Evolutionary dynamics of neoantigens in growing tumors.
*Nat Genet.*2020; 52: 1057-1066 - The immune checkpoint kick start: optimization of neoadjuvant combination therapy using game theory.
*JCO Clin Cancer Inf.*2019; 3: 1-12 - Anti-PD-1/PD-L1 therapy for non-small-cell lung cancer: toward personalized medicine and combination strategies.
*J Immunol Res.*2018; 2018 - Rationale of combination of anti-PD-1/PD-L1 antibody therapy and radiotherapy for cancer treatment.
*Int J Clin Oncol.*2020; 25: 801-809 - Two may be better than one: PD-1/PD-L1 blockade combination approaches in metastatic breast cancer.
*NPJ Breast Cancer.*2019; 5: 1-9 - 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 - 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 - 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

## Article info

### Publication history

Published online: July 29, 2021

Accepted:
July 24,
2021

Received in revised form:
June 24,
2021

Received:
March 23,
2021

### Identification

### Copyright

© 2021 The Author(s). Published by Elsevier B.V.

### User license

Creative Commons Attribution – NonCommercial – NoDerivs (CC BY-NC-ND 4.0) | How you can reuse

Elsevier's open access license policy

Creative Commons Attribution – NonCommercial – NoDerivs (CC BY-NC-ND 4.0)

## Permitted

### For non-commercial purposes:

- Read, print & download
- Redistribute or republish the final article
- Text & data mine
- Translate the article (private use only, not for distribution)
- Reuse portions or extracts from the article in other works

## Not Permitted

- Sell or re-use for commercial purposes
- Distribute translations or adaptations of the article

Elsevier's open access license policy