Immunogenicity, the Mounting of an Anti-tumor Response and its Resistance by the Tumor

The immunogenicity of a tumor, which is the ability of an immune response to be mounted against it, is a function of its tumor burden (i.e., the abundance of tumor antigens recognized by the immune cells) and its antigen mutational landscape (i.e., the variety and strength of those antigens)⁷. Tumors with antigens more abundant, stronger and varied in sequence, are more immunogenic and may activate more T cells against the tumor^8–10.

Tumor-associated or tumor-specific antigens (TAA or TSA, respectively) are usually proteins, peptides or other factors that are products of the cancer cells¹¹. In larger tumors (with higher tumor burden), more TAA and TSA are released through natural death and decay of cancer cells during growth, and during killing of cancer cells by lymphocytes. TAA (used from now on) can be picked up by antigen presenting cells (APCs), which mature (to mAPCs) and present the TAA to naïve T cells during priming to induce their differentiation into activated tumor-targeting effector T cells, which then can eliminate the tumor through cytotoxic activity.

During the process of antigen presentation, each naïve T cell can recognize a specific TAA by its unique T cell receptor (TCR). A T cell with a TCR specific to a TAA is called a clone. There are usually several copies of each clone specific to each TAA available for priming and activation. The more unique TAA a tumor releases, the more T cell clones can be primed and activated to target the tumor. Additionally, the strengths of the different antigens have been classified from weak to strong; the exact reason for why one antigen is stronger or weaker than another is not entirely clear and has been linked to several factors, including the affinity (K_d) and binding ability between the antigen presented on the major histocompatibility complex (MHC) and its cognate TCR during antigen presentation. Upon the successful elimination (death) of cancer cells by the effector T cells, TAA associated with those cancer cells are released into the environment and thus can be used to prime more effector cells to mount a subsequent attack against the tumor; some non-specific cytotoxic activity by the T cells is also thought to occur by infiltrating effector T cells¹².

As described, an effective immune response is dictated by the number of T cell clones, tumor burden (or antigen abundance, which is assumed proportional to tumor size) and antigen strength. Considering that an effective immune response can be generated, to eliminate the tumor it must overcome its growth rate and resistance mechanisms. The focus of this study is on the resistance to cancer elimination by immune checkpoints Cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), Programmed cell death protein 1 (PD-1), and Programmed death-ligand 1 (PD-L1) and how their blockade by immune checkpoint inhibitors can lead to a range of patient responses through variations in specific physiological parameters, such as those described above (i.e., numbers of T cell clones, tumor size and antigen strength), encoded in a mechanistic framework. Other key parameters mechanistically considered are the expression levels of the immune checkpoints (particularly PD-L1 on cancer cells, which is a known biomarker of anti-PD-1 efficacy), and inhibitory immune cells, such as T-regulatory cells present in the lymph nodes and the tumor microenvironment. While immune checkpoints in recent years have been found to be prime targets for cancer immunotherapy in melanoma, non-small cell lung cancer (NSCLC), bladder cancer and other cancer types, we particularly examine them here for melanoma^13,14.

Immune checkpoint blockade therapies

Both PD-1 and PD-L1 are expressed on effector T cells, cancer cells, as well as the mentioned regulatory immune cells in the tumor microenvironment (TME)^15–17. When PD-1 interacts with PD-L1, the effector T cells can be prevented from cytotoxic activity and may become anergic; in this way the tumor hijacks the signaling pathways for its own survival. Furthermore, CTLA-4 checkpoints have been found on T cells and play an important role during priming in the activation of effector T cells; their blockade can help boost effector T cell numbers^18,19. Other checkpoints have also been identified: those found that are most relevant to this study include the ones mentioned above, as well as the co-stimulatory B7 receptors (CD80 and CD86) and CD28, and PD-L2, which also interacts with PD-1^20–23.

Mathematical models for immunotherapy and their utility

Several mathematical models of the immune system have been developed thus far; only a small number of them focus on immuno-oncology and immunotherapies, as described here^24–26. Previous models of immunotherapies have focused on single chemotherapies, combinations of chemotherapies, radiation therapy, a cancer vaccine and IL-2 immunotherapy, or focused on a specific molecule or pathway^27–31. Quantitative Systems Pharmacology (QSP) approaches to cancer immunotherapy are currently being pursued in academia and industry³². It is our belief that computational models, such as the one presented here, will play a substantial role in reducing costs and expediting the drug development and approval processes to bring therapies to market quicker, especially for combination therapies³³.

Much is still unknown about how and why only some patients benefit from particular therapeutic regimens^34,35. In order to make personalized cancer immunotherapy a reality, it is necessary to gain a better understanding of the physiological system dynamics at play on both the cellular and molecular levels for each cancer patient, which may ultimately help determine the optimal personalized dosing, regimens and combination strategies. Here, we aim to link how the mentioned physiological parameters that we described biologically (i.e., numbers of T cells clones, tumor size, average antigen strength, expression of PD-L1 on cancer cells, as well as the levels of inhibitory cells in the tumor and lymph node) when put into a well-mixed mechanistic framework and their values varied across defined ranges, can allow us to simulate a variety of anti-tumor responses when the dose level and regimen of each therapy is kept constant. We further show that the median values of our model parameters can be used to predict median responses to various combinations of CTLA-4, PD-1 and PD-L1 blockade therapies that have been reported in the literature specific to the doses and regimens. Lastly, we show the significances of the parameters we have varied in terms of how they may correlate to different types of patient responders, in addition to how variations in these parameters can lead to different types of tumor dynamics. We use this model to make predictions on outcomes of combination therapies and believe it can be used to predict various dosing regimens for optimization purposes within the given framework.

Simulations and predictions of anti-CTLA-4 monotherapy

Prior to utilizing the model to make predictions, we analyzed the simulated response for anti-CTLA-4 therapy (Fig. ^S1) and compared it to our understanding of response on a qualitative level. Response was similar to that seen in⁴⁴ with a common regimen represented by Ku et al., (ClinicalTrials.gov number, NCT00324155)⁴⁵; the boost in peripheral effector T cells was the cause of the response in the model, as also shown clinically⁴⁶. We show a clear and wide dose response between 0.3 and 10 mg/kg (red to green) in Fig. ^S1A, and for the 3 mg/kg case, simulations show that approximately 20% of patients would be responders (Fig. ^S1C, left). This would be slightly below or around 20–28.5% ORR seen with CTLA-4 therapy reported for advanced melanoma⁴⁷. From this point on, we utilized the model to make predictions via simulations and compared them to clinical outcomes of published data.

Simulations and predictions of anti-CTLA-4/anti-PD-1 combination therapy

In Fig. 3, we show our predictions for anti-CTLA-4/anti-PD-1 combination therapy, based on the trial regimen as carried out in Wolchok et al. (ClinicalTrials.gov number, NCT01024231)⁴⁸, and represented by the PK in Fig. 3A. The trial was simulated using 47 VPs, as done in the clinical trial with 1 mg/kg anti-PD-1 and 3 mg/kg anti-CTLA-4. Figure 3B shows that a higher proportion of VPs can be expected to be PR/CR responders, although some could have progressive disease or SD (left). In Fig. 3B, right, we show that the median of the clinical trial data and the median from our simulated virtual clinical trial overlap to a large extent within the narrow confines of the 90% confidence interval. We further see that the clinical trial also had patients who had stable and progressive disease, although to a much lower extent than those who showed partial or complete responses. Of our VP population, 85% were PR/CR responders (Fig. 3b, right), compared to clinical results showing 79% having robust responses, and with further analysis we indicate that the number of T cell clones and an increase in maximum effector T cell density in the tumor to be the main significant driving factors for response in our virtual population; similar to anti-PD-1 monotherapy case, except starting tumor diameter and %PD-L1 at therapy were less significant and the density of effector T cells in the tumor showed a wider range for responders.

Open in a separate window

Figure 3

Dose response and clinical validation of anti-PD-1/anti-CTLA-4 combination-therapy. (A) From top to bottom: Tumor response to combination therapy at doses of 0.3, 1, 3 and 10 mg/kg of anti-CTLA-4, as represented by the colors in the bottom figure in ascending order and 3 mg/kg for anti-PD-1 was used for all simulations; the blue line indicates no therapy (top figure), and orange indicates only anti-PD-1. Then, effector T cell density in the tumor (second from the top), mAPC density in the lymph nodes (third from the top) and finally, the PK of anti-PD-1 and lastly, anti-CTLA-4 at the given doses. For all following figures, 1 mg/kg anti-PD-1 and 3 mg/kg anti-CTLA-4 were used, following the same regimen. (B) Diversity of tumor response (left), prediction of median clinical response data (right). (C) Waterfall plot of VPs (left) and pie chart (right) with percent of virtual non-responders (NR), stable disease (SD) and partial or complete responders (PR/CR). (D) Bar graph comparison parameters varied in model for each responder type (left) and box plots of significant differentiators (right). (E) Max effector T cell density in the tumor (left) and average mAPC density in the lymph nodes (right) for each responder category.

Simulations and predictions of anti-PD-L1 monotherapy

We further employed the same kind of analysis to anti-PD-L1 monotherapy in comparison to a trial run by Brahmer et al., (ClinicalTrials.gov number, NCT00729664)⁴⁹, where 20 mg/kg anti-PD-L1 was administered Q2W. Figure 4A shows that much like for anti-PD-1 therapy, the dose response is largely flat. As shown previously, we demonstrate that a diversity of responses across a VP population of 200 is possible (Fig. 4B, left) and that the simulated median within a 30% confidence interval largely overlaps with the clinically determined median (Fig. 4B, right). We further predict that overall clinical response would be lower for anti-PD-L1 than that shown for anti-PD-1, where 35% (Fig. 4C) showed overall response compared with 43% in our VP population, respectively, which is in line with what was seen qualitatively in the mentioned trials. Interestingly, we determine PD-L1 therapy as lower in efficacy compared to PD-1 therapy for melanoma, which may be a consequence of the associated checkpoints (i.e., CD80, PD-1 and PD-L2, expressed at literature reported averages for melanoma). In fact, PD-L2 expression can help overcome PD-L1 blockade, since the latter cannot prevent the former’s interaction with PD-1; while anti-PD-1 blocks both interactions of PD-1 with PD-L1 and PD-L2 and CD80 is minimally expressed in comparison. The effects of varying CD80, PD-1, PD-L1 and PD-L2 expression levels on response during each therapy are examined in a sensitivity analysis, described later. In a similar manner, as for anti-PD-1 therapy, we find that number of T cell clones, %PD-L1 expression and tumor diameter at therapy are significantly different between non-responders and those showing partial or complete responses (Fig. 4D). We further show the same type of trend for maximum effector T cell density in the tumor and average mAPC density in the lymph nodes as for anti-PD-1 therapy.

Open in a separate window

Figure 4

Dose response and clinical validation of anti-PD-L1 monotherapy. (A) From top to bottom: Tumor response to anti-PD-L1 therapy at doses of 0.3, 1, 3, 10, 15 and 20 mg/kg as represented by the colors in the bottom figure in ascending order; the blue line indicates no therapy (top figure). Then, effector T cell density in the tumor (second from the top), mAPC density in the lymph nodes (third from the top) and finally, the PK of anti-PD-L1 at the given doses. For all following figures, 20 mg/kg anti-PD-L1 was used, following the same regimen. (B) Diversity of tumor response (left), prediction of median clinical response data (right). (C) Waterfall plot of VPs (left) and pie chart (right) with percent of virtual non-responders (NR), stable disease (SD) and partial or complete responders (PR/CR). (D) Bar graph comparison parameters varied in model for each responder type (left) and box plots of significant differentiators (right). (E) Max effector T cell density in the tumor (left) and average mAPC density in the lymph nodes (right) for each responder category.

Sequential therapy of anti-CTLA-4 first then anti-PD-1 vs. the reverse for induction

Our final comparison to a human clinical trial was performed on a more qualitative basis. We assessed the model’s ability to emulate sequential therapy responses based on the induction therapy trial published by Weber et al., (ClinicalTrials.gov number, NCT01783938)⁵⁰. Here, anti-CTLA-4 administered first followed by anti-PD-1 therapy was compared to the reverse sequence (Fig. 5). Overall, we found our results to indicate that on average (and for individual cases) anti-PD-1 administered first followed by anti-CTLA-4 (Fig. 5A, blue) produced a greater effect on tumor response than anti-CTLA-4 first followed by anti-PD-1 (Fig. 5A, red), consistent with clinical results. It was interesting to see these results, even though in both scenarios the ranges and medians of maximum effector T cells and average mAPC densities (Fig. 5B,C, bottom) were about equal. Two significantly differentiating parameters, however, were the number of T cell clones and %PD-L1 expression on the cancer cells (Fig. 5B,C, second from bottom).

Open in a separate window

Figure 5

Sequential therapy of anti-PD-1 and anti-CTLA-4. (A) From top to bottom: Tumor response to sequential therapy with dosing as listed in phase II CheckMate 064 trial (NCT01783938) with regimens shown in bottom two figures. Following response (top), effector T cell density in the tumor (second from the top), mAPC density in the lymph nodes (third from the top) and finally, the PK of anti-PD-1 first regimen and lastly, anti-CTLA-4 first regimen. (B) Anti-PD-1 therapy first: Comparison of variability in responses (top), waterfall plots of responders in virtual trial (second from top), Bar graph comparison parameters varied in model for each responder type (third from top), and boxplots of Max effector T cell density in the tumor and average mAPC density in the lymph nodes (bottom) for each responder category: virtual non-responders (NR), stable disease (SD) and partial or complete responders (CR/PR). (C) Same as (B), except for anti-CTLA-4 therapy administered first.

Priming and activation of effector T cells

The model consists of several compartments for effector T cell trafficking once they exit the TDLN following priming and activation, as described previously by Zhu et al. Here, we built the model encompassing the following trafficking compartments: lymph node, blood, tumor, lung, GI tract, spleen, and liver, while the remainder of the organs were lumped into the peripheral compartment³⁸. For each organ, the vascular and interstitial volumes are represented, as well as the blood and lymph flow into and out of the organ. As in the physiologically-based model, effector cells can bind reversibly to the vasculature in the space of the lung, liver, spleen and tumor, attach irreversibly, becoming arrested and extravasate into the interstitial space from where they can circulate through the lymph and back into peripheral circulation. In the tumor, however, once the effector cells enter the interstitial space, they stay there to interact with the tumor cells. Elimination is set to occur in the lungs, as described in the mentioned study, and in the tumor, based on interactions with cancer and supporting cells. The antibodies are added to the central (blood) compartment and follow a minimal-PBPK layout. The cellular dynamics in the model are the prime elements that define the outcomes of the immunotherapy regimens, while the molecular dynamics are the defining factors for the cellular dynamics. Overall, the model consists of 282 ordinary differential equations (ODEs) and 218 algebraic equations and is implemented using the SimBiology plugin in MATLAB (MathWorks, Natick, MA). To ensure reproducibility of the model, the complete set of governing equations and parameters are presented in the supplement and the SBML code is provided.

In the TDLN compartment, the immune response is generated through the activation of tumor-targeting effector T cells resulting from the presentation of antigens by mature antigen presenting cells (mAPCs). Activated effector T cells migrate out of the TDLN and are trafficked as described above. In total, the priming process has been shown by intravital microscopy to occur in three distinct stages in the LN, which are incorporated in the model⁵⁴. The possibility of a successful priming interaction between mAPCs and T cells is represented by an algebraic equation that was extrapolated from an agent-based model, which was coupled with experiments to calculate the probability of at least one successful priming interaction to occur over the designated time⁵⁵. The final equation and an example of how it is used is shown in Eqs ^1.1 and ^1.2

Equation ^1.1 shows that the association of Naïve T cells (NT) with non-interacting mAPC (mAPC_nInt) during the first priming stage to become a single state variable (NT1) representing the associated state, is determined by an association rate (k_Assoc,P1 defined by measured literature values), the number of NT undergoing the interaction, the strength of the priming co-receptor (where SigCD28 represents the CD80/86 interaction with CD28 on mAPCs based on receptor occupancy of 0–1.0), and the probability of a successful priming interaction (P1). In Eq. ^1.2, P1 considers the number of mAPCs available to associate with the NT cells (mAPC_nInt), as well as other factors that define a value function representing the possibility of a successful priming interaction to occur. These include: The number of NT cells available for priming (calculated from the total number of T cells in the TDLN, T, and the fraction which are available for priming (ϕ_NT)), their migration or diffusion rate in the TDLN (Diff_T), the T cell clonality (TClonality), the strength of the antigen landscape during presentation (Antigen), a scaling factor for the possibility of a successful priming interaction between an mAPC and an NT cell (S_{NT_mAPC}), the geometry of the TDLN, the diameter of each T cell (D_T) and the total time over which priming occurs (t_max). The number of mAPCs available for priming at any given time is algebraically calculated. Other equations during priming govern the dissociation of mAPCs and NT during the first priming stage and the association of primed NT with free mAPCs for the second priming stage; the equations are listed in the supplement. Thus, the more mAPCs associate with each state (up to the maximum possible) and the stronger the antigen is, the more cells mature to the next stages and ultimately the higher the number of activated effector T cells are generated; this is the mAPC dose response of the T cells during priming⁵⁶.

A successful priming interaction is one that leads to maturation of the NT cells, whereas an unsuccessful one has been shown to lead to a relatively quick dissociation between the NT and mAPCs without the maturation of the NT cells⁵⁷. Once both stages of priming have occurred in a similar manner, the proliferation parameters were optimized to represent division of the cells three times a day over about five days at the maximum levels to produce effector T cells. When disregarding inhibitory checkpoint signals in total, the priming stage takes about a week as Naïve T cells mature to become active effector T cells^54,58.

Once primed and activated, effector T cells leave the lymph nodes via the efferent lymphatics and make their way into the blood, from which they access the tumor and other tissues⁵⁹. In a human body, there are about 500–600 lymph nodes, but since we consider a tumor localized to a particular region of the body, we assume that only a fraction of total lymph nodes are involved in producing an anti-tumor immune response^60,61. In total, we approximated 20 distinct regions that have clusters of lymph nodes in the human body, and therefore we estimated that about 30–35 TDLN are involved in the anti-tumor response in a particular given region⁶². Since one TDLN is considered in the compartment in the model, in order to simulate multiple TDLN, the number of effector T cells generated in the TDLN compartment that enter the blood (or the central compartment) are multiplied by the number of TDLN considered in the model. In this manner, we assume that each TDLN functions in the same manner. Furthermore, we assume that each TDLN receives the same number of mAPCs and TAA that are transported from the tumor. Once in the blood, a chemokine factor further promotes effector T cell infiltration into the tumor⁶³. In fact, immunotherapies to enhance lymphocyte-recruiting chemokine production by the tumor can also be explored in this model⁶⁴. Furthermore, there is evidence that tumor lymphangiogenesis promotes T cell infiltration, which can be factored into the simulations through the chemokine factor⁶⁵.

Interactions between effector and target cells, and elimination of cancer cells

The immune response is initiated from the tumor. While in the tumor compartment, the effector T cells interact with and kill the cancer cells through cytotoxic activity, and the dead cancer cells release TAA. In general, all effector-target interactions in the model are represented by a first interaction step (where two cell types become represented by one state variable), and subsequently the effect the interaction has on the target and effector cells over the course of the dissociation step. The interaction steps in each case are defined by a validated model for effector cell killing of cancer cells⁶⁶, and the rate of the subsequent step is defined by the measured rates of interaction between the different cell types from published literature referenced in the supplement. In Eq. ^1.3, E refers to the effector cells producing the response, while T refers to the target cells. Algebraic equations govern the accounting of the cells during this process in each case, and the effect on each cell type following interaction is dependent upon the immune checkpoint signaling. d is the rate of association for killing per day, λ is the exponent of fractional cell association for killing and s is the steepness coefficient

For cancer cells specifically, their death produces the release of TAA to a defined level per cell, as described previously. In parallel, cancer cells that naturally decay during tumor growth also release TAA. Monocytes migrate into the tumor, differentiate into APCs, which pick up the TAA and become mAPCs. The mAPCs then migrate back to the TDLN via the lymphatic vessels for presentation of the TAA to Naïve and Primed Naïve T cells, during priming for effector T cell activation, and the process begins again. Additionally, the TAA that are not picked up in the tumor by APCs can flow back to the TDLN through the lymphatic vessels, where they are picked up by resident APCs, which would subsequently be involved in the priming process, leading to effector T cell activation. Altogether, this process is a feedforward loop where continuous amplification of effector T cells occurs, and the cycle continues until there are no more cancer cells left to release enough TAA that is required for priming.

Immune checkpoint signaling and inhibitory effects

The presented scenario considers the rate of natural turnover of the different immune and cancer cell types (i.e. effector T cells, mAPCs and cancer cells) in the human body. In the model, we specifically focus on the molecular species that are either essential for our analysis of the immune response (e.g. CD28:CD80/CD86 binding, or the overall strength of the antigen landscape being presented during priming), or are the direct or indirect targets of the therapies we are testing in the model (e.g. the cognate receptors of CTLA-4, PD-1, and PD-L1 immune checkpoints).

The regulatory immune cells include the Tregs and MDSCs, which are present in the TME along with the cancer cells; Tregs are also present in the TDLN. These cells reside in their respective compartments and do not migrate from one compartment to another like the effector T cells do. Furthermore, their numbers are a function of the other cell types, which change with time within their respective compartments: In the TDLN, Tregs are defined as a certain percentage of the total number of T cell as determined experimentally by flow cytometry (5–35%), while in the tumor Treg are defined as a percentage (0–0.25%) of the number of cancer cells, and the number of MDSCs are two to four times the number of Tregs, based on literature evidence^67–69. Note that all numerical values of the parameters here and elsewhere in the paper are baseline parameters and can be varied; these values are not the principal part of the model. Overall, these regulatory cells are known to play a role in inhibiting the anti-tumor immune response, and they carry out their activities through various molecular mechanisms, including inhibitory checkpoint receptors and other factors that target effector T cells. Cancer cells also employ inhibitory checkpoint receptors and other mechanisms as a means of protecting themselves, which slow down or inhibit effector T cell activities. The model incorporates checkpoint receptors expressed on the surface of cells and interactions between cells within the volume of the immunological synapses that are formed during cell-cell interactions⁷⁰.

In addition, mAPCs also exhibit a variety of mechanisms that regulate the activation of effector T cells, which can stymie or boost an effective anti-tumor immune response. These mechanisms include the types and origins of the antigens that they present, the ratios of mAPC-to-Naïve and Primed Naive T cells during priming (where, generally, the greater the antigen abundance, the higher the number of mAPCs), the presence of positive immune checkpoints (including co-receptors) and (like the regulatory cells) inhibitory immune checkpoints. During antigen presentation, the B7 co-receptors (CD80 and CD86) on the mAPCs must engage CD28 on the Naïve and Primed T cells, in addition to antigen presentation, for priming to occur. Additionally, a stronger recognized antigen and a higher ratio of mAPC-to-T cells results in a greater extent of priming; i.e., a larger number of activated effector T cells are produced during priming as a result of greater activation and proliferation. Lastly, the inhibitory immune checkpoints expressed on the mAPCs, which are recognized by the cognate checkpoints on the T cells, act to limit the extent of priming and ultimately limit the number of activated effector T cells generated in the model. The specific cell types and the corresponding immune checkpoints in the model are depicted in Fig. 1B. All of the interactions between the ligands and their cognate receptors that are represented in the figure take into account experimentally determined K_d values for their binding, receptor expression values measured by flow cytometry (or their approximate expression levels) and the volumes in which they interact (being either within the immunological synapses, or within the entire volume of the designated compartment where the interaction is occurring). Note that in the model, we do not account for the mechanisms of antigen processing and we do not know the variety of antigens present and their binding affinity values; therefore, we use a general antigen strength value from 0.1–1.0, which for melanoma is known to be on the higher end and is so set between 0.8 and 1.0 for all simulations and fitting, as mentioned.

Representations of tumor heterogeneity

While tumor heterogeneity is not represented on a spatial scale in this version of the model, it is represented in terms of checkpoint heterogeneity. The expression of the different checkpoints on different cancer cells within a tumor can be heterogeneously distributed, but more importantly in this analysis is the percent of cancer cells that express particular checkpoints. Considering that a cancer cell expresses a certain immune checkpoint, it would affect an effector T cell differently than if the cancer cell were to express two different immune checkpoints to the same extents. The problem translates to the fact that if we state that 60% of the cancer cells express PD-L1 for one tumor, for example, while 5% express PD-1 from an independent tumor study, we do not know if 5% of the cancer cells should express both immune checkpoints in our model (which represents a single unified tumor), or if less than 5% express both checkpoints; in other words, PD-L1 and PD-1 can be expressed on the same or distinct cancer cells and their expressions can overlap to varying extents across the tumor, which would affect the lymphocytes differently depending on which cells they interact with. Although the model is considered “well-mixed” in general, unlike spatially-resolved agent-based models⁷¹, local interactions are accounted for distinctly by separate compartments and submodules.

Therefore, in this model, the percent expression of each immune checkpoint is an input, and for the presented simulations we assumed that the immune checkpoints are equally and independently distributed across the tumor; this assumption can be easily relaxed. For example, if we consider checkpoints A and B, at 50% expression each across the tumor, then 25% would express A and B, 25% would express only A, 25% would express only B, and 25% would express neither A nor B. This distribution can be modified to fit whatever the experimental data determine is the correct distribution for the tumor; however, at this stage we have not seen any consensus on how the percentages of immune checkpoints across the tumor are distributed across the different cancer cells, which resulted in the stated assumption for modeling purposes. In the model, while the number of immune checkpoints considered can be expanded, we consider four different immune checkpoints that are expressed on the cancer cells that have overlapping interactions with the targets of our therapies, PD-1, PD-L1, PD-L2 and CD80, which translates to 2ⁿ different checkpoint expression combinations across the tumor (where n=4, the number of different checkpoints considered). This in turn translates to 16 unique cancer cell-effector T cell signaling interactions, assuming all effector T cells have the same expression of immune checkpoints.

The advantage of structuring the model in such a manner is that we can translate percentage of different checkpoints within the tumor that are found individually on cells into a unified single tumor model, and in the future, as more checkpoint expression data are captured from a single tumor regarding the percentage of checkpoint expressions relative to one another, they can be entered in the model. Thus, the distribution of immune checkpoints on the cancer cells in the model can be modified if need be. Furthermore, it is useful to note that in the current model and for the following analyses we assume that no particular cancer cell type which has a certain sequence of immune checkpoints expressed on is ever fully removed from the system in terms of the percent it is present in the tumor as a whole; instead, as cancer cells get killed in the model, the tumor as a whole retains the same constant heterogeneous checkpoint distribution, fraction-wise. In the future, it would be interesting to compare the results obtained here with those when certain cancer cells with specific checkpoint expressions are able to be eradicated first or with some stochastic variability in checkpoint expression based on certain criteria.

Pharmacokinetics of immune checkpoint antibodies

Antibodies are delivered to the central compartment (blood) as a zero-order intravenous infusion over the course of one or one and a half hours, as prescribed for each antibody. The antibodies are also cleared from the central compartment through a first order elimination. The distribution of the antibody was designed using a minimal PBPK model³⁹ between the central, tight and leaky tissues and lymph node compartments. Distribution between the central and tumor, and between the tumor and lymph node compartments are calculated by considering the surface area of blood or lymph vessels in the tumor volume and permeability of antibodies across that surface area, as represented in Finley et al.⁷²; and the PK tumor compartment for the PK volume was varied with the true tumor size in the model, given this approach.

For the estimated PK parameters of all the antibodies used in this model, the concentrations of the antibodies in the central compartment were fit to available mean patient population serum concentrations taken over time, displayed in Fig. ^S4. For anti-PD-1, the data came from Nivolumab, for anti-PD-L1 it came from Durvalumab, and for anti-CTLA-4, it came from Ipilimumab published sources listed in the supplement. Fitting of each therapy to 3 mg/kg single dose administration is shown in Fig. ^S4, along with the PBPK outputs for the other tissues, including the tumor, which was found to have a concentration about 15–20% of the serum concentration within a few days, and having an increased ratio with time due to an enhanced permeability and retention (EPR) effect⁷³. The concentration in the lymph was about half of the tumor concentration, or less. While target-mediated drug disposition (TMDD) was only represented for CTLA-4, there is literature evidence that TMDD for anti-PD-1 does not occur⁷⁴ and we assumed that to be the case for anti-PD-L1 as well; although it could be easily incorporated should the need arises. The free fraction is reported as the concentration in each compartment. Once reaching their sites of activity in the TDLN or tumor, each antibody would bind to its respective target, which would influence the pharmacodynamics (PD) of the immune response and the tumor.

Immune checkpoint signal modeling

Multiple immune checkpoints were considered to have effects on cells at specific points. The combined effect of the signaling was calculated as follows in Eq. ^1.4

This equation was based on the generalized equations derived by Chou and Talalay to determine summation of effects (f_a) of two or more mutually exclusive drugs⁷⁵. In order to avoid having to calculate the potency of each checkpoint interaction, D in the original equation was replaced by receptor occupancy (RO) and instead of D50, the maximum RO possible during signaling for each respective receptor was used (RO_max). Lastly the weights of receptor signaling (w_i) were assumed to take on the fraction of total receptor occupancy for all the receptors, as shown in Eq. ^1.5; this turned out to be a consequence of the receptor expression levels, their binding affinities to the respective antibodies, any competing binding interactions and the concentration (or exposure) of the antibodies in the environment. For all signaling, the Hill-coefficient (m) was greater than 1 to represent the higher-order system, which denotes the sigmoidicity of the dose-effect curve. The same values for all the variables were used consistently throughout the model, with exception of the RO related values, which were calculated dynamically based upon receptor/ligand interactions.

Parameter optimization

The Multistart function along with lsqnonlin were run using the parallel computing toolbox in MATLAB to perform global optimization on single or multiple parameters in the model in a constrained manner. This allowed for multiple local minima to be found by lsqnonlin for each parameter by starting at a defined number of randomly dispersed initial values within the constrained search range. Then, the local minima that were found were compared to one another and the approximate global minimum was found.

Supported by grants from MedImmune and NIH R01 CA138264 to ASP. The authors thank members of the Popel laboratory Richard Sove and Hanwen Wang, for reading the manuscript and useful comments.