SIAM News Blog
Research

Optimizing Immune Checkpoint Blockades with Complementary Modeling Strategies

By

Many types of cancer are notoriously difficult to treat due to the complex nature of interactions between cancer and immune cells. A healthy immune system utilizes checkpoint proteins on T-cells to prevent an extreme immune response, maintaining a system of checks and balances that allows the body to appropriately respond to threats such as cancer or infection without triggering autoimmunity. However, some cancers can exploit these proteins and “trick” the T-cells into thinking that cancerous cells are actually healthy. Affected T-cells then become trapped in a downregulated or “exhausted” state that hinders their ability to stop the cancer’s progression. Standard-of-care treatments for these types of cancers include immune checkpoint blockades: a form of immunotherapy that is meant to block immune cells’ inhibitory receptors (i.e., checkpoint proteins), preventing this downregulation and in some cases permitting exhausted T-cells to regain their ability to recognize and fight cancer cells [1, 6, 7].

However, there is certainly room for improvement in this strategy. Immune checkpoint blockades are commonly applied to treat melanoma—where response rates range from \(25\) to \(45\) percent—and non-small cell lung cancer, with response rates of \(8\) to \(40\) percent [4, 5, 9, 13]. Our recent work focuses on a complementary dynamical system and agent-based modeling approach to optimize the use of checkpoint blockades and improve response rates in cancer patients. 

We began by implementing a system of ordinary differential equations (ODEs) to evaluate the initial impact of a checkpoint blockade on a simulated tumor. Our ODE system describes interactions between T-cells \((T),\) cancer cells \((C),\) and non-cancerous antigen-presenting cells \((A)\): 

\[\begin{eqnarray}
  \frac{dT}{dt} &=& b_T+\frac{r_1TC}{k_1+C}-a(A+C)T-d_TT \tag{1}\\\\
  \frac{dC}{dt} &=& b_CC\bigg{(}1-\frac{C}{k_2}\bigg{)}-r_2CT-d_CC \tag{2}\\\\
  \frac{dA}{dT} &=& b_A-d_AA, \tag{3}\\\\
  \textrm{where}\; a &=& \frac{1}{n}IR_1 R_2.\tag{4}
\end{eqnarray}\]

We modified the above equations from the existing literature [3,10] by introducing a term for immune exhaustion control \((I)\) that explicitly considers blockade efficiency. We also accounted for T-cell robustness by sweeping through the parameter \(n,\) which represents the minimum number of interactions before a cell becomes exhausted. More information about parameter interpretations and sampled values is available in the literature [11].

<strong>Figure 1.</strong> Ordinary differential equation simulation for tumor-immune dynamics and select local sensitivity analysis. <strong>1a.</strong> Phase plane that represents cancer cell versus T-cell count in the simulated tumor over the course of \(100\) simulation runs. Each trajectory represents one tumor. The blue dashed lines with arrows trace one trajectory that illustrates each fate: progressive disease (high cancer, low T-cell in the left region of the plane) or immune control (high T-cell, low cancer in the bottom region of the plane and magnified by the inset). <strong>1b.</strong> Histogram of terminal cancer cell count from the simulation in 1a that highlights progressive disease or stable disease steady state as a binary outcome. <strong>1c.</strong> Surface that represents the effect of varying the nondimensionalized cell birth-immunotherapy parameter subset on cancer cells. The surface discontinuity again underscores a simulated binary outcome of progressive or stable disease. <strong>1d.</strong> Surface that represents the effect of varying the nondimensionalized cell death-immunotherapy parameter subset on cancer cells. The surface discontinuity again underscores a simulated binary outcome of progressive or stable disease. Figure adapted from [11] under the<a href="https://creativecommons.org/licenses/by/4.0" target="_blank"> Creative Commons Attribution 4.0 International Deed</a>.
Figure 1. Ordinary differential equation simulation for tumor-immune dynamics and select local sensitivity analysis. 1a. Phase plane that represents cancer cell versus T-cell count in the simulated tumor over the course of \(100\) simulation runs. Each trajectory represents one tumor. The blue dashed lines with arrows trace one trajectory that illustrates each fate: progressive disease (high cancer, low T-cell in the left region of the plane) or immune control (high T-cell, low cancer in the bottom region of the plane and magnified by the inset). 1b. Histogram of terminal cancer cell count from the simulation in 1a that highlights progressive disease or stable disease steady state as a binary outcome. 1c. Surface that represents the effect of varying the nondimensionalized cell birth-immunotherapy parameter subset on cancer cells. The surface discontinuity again underscores a simulated binary outcome of progressive or stable disease. 1d. Surface that represents the effect of varying the nondimensionalized cell death-immunotherapy parameter subset on cancer cells. The surface discontinuity again underscores a simulated binary outcome of progressive or stable disease. Figure adapted from [11] under the Creative Commons Attribution 4.0 International Deed.

Using this framework, we observed that simulated tumors are more likely to enter a controlled steady state, or remission, when \(80\) to \(90\) percent of inhibitory receptors are blocked on the T-cell surface. At this range of blockade efficiencies, more than \(50\) percent of simulated tumors enter the controlled (i.e., high T-cell, low cancer) steady state (see Figures 1a and 1b). In contrast, less efficient blockades resulted in a greater likelihood that simulated tumors would enter the progressive disease (i.e., high cancer, low T-cell) steady state. Local sensitivity analysis also revealed a relationship between immune checkpoint blockade efficiency and T-cell birth/death rates, which indicates a need to maintain a critical threshold of immune activity to achieve tumor control (see Figures 1c and 1d). Less robust immune responses require greater blockade efficiency to achieve the same level of tumor control.

Next, we implemented an agent-based model (ABM) to investigate the individual interactions and behaviors that are associated with immune checkpoint blockades at the cellular level. Specifically, we defined and initialized tumor and T-cells on a grid. The simulated cells were advanced off lattice through time and space according to the Lennard-Jones potential and Verlet algorithm (see Figures 2a and 2b). We tallied interactions with cancer cells for each T-cell and assigned an exhaustion probability according to the simulated blockade, which yielded an exhaustion score for every cell (see Figures 2c and 2d). This score indicates the potential for local patterning, i.e., clusters of cells with higher or lower degrees of exhaustion potential. The initial T-cell exhaustion process typically occurs rapidly—as early as six hours after infiltration [8]—so this stage of our model focuses primarily on early time points.

<strong>Figure 2.</strong> Agent-based model that represents solid tumor and infiltrating T-cell exhaustion trajectories in a subset of the solid tumor microenvironment. <strong> 2a.</strong> Initialization of the tumor-cell grid (in magenta) and T-cell layer (in blue). <strong> 2b.</strong> Zoomed-in screen capture of cell positions that depicts T-cell infiltration in the cancer cell layers and tumor expansion due to mechanical forces \((\textrm{time step} = 300).\) <strong> 2c. </strong> Number of exhausting interactions in the absence of an immune checkpoint blockade, indicated by the color bar. <strong> 2d.</strong> Predicted number of exhausting interactions in the presence of a blockade that is \(85\) percent effective, indicated by the color bar. This simulation highlights spatial regions of projected T-cell exhaustion over time in both the presence and absence of a checkpoint blockade. Figure adapted from [11] under the<a href="https://creativecommons.org/licenses/by/4.0" target="_blank"> Creative Commons Attribution 4.0 International Deed</a>.
Figure 2. Agent-based model that represents solid tumor and infiltrating T-cell exhaustion trajectories in a subset of the solid tumor microenvironment. 2a. Initialization of the tumor-cell grid (in magenta) and T-cell layer (in blue). 2b. Zoomed-in screen capture of cell positions that depicts T-cell infiltration in the cancer cell layers and tumor expansion due to mechanical forces \((\textrm{time step} = 300).\) 2c. Number of exhausting interactions in the absence of an immune checkpoint blockade, indicated by the color bar. 2d. Predicted number of exhausting interactions in the presence of a blockade that is \(85\) percent effective, indicated by the color bar. This simulation highlights spatial regions of projected T-cell exhaustion over time in both the presence and absence of a checkpoint blockade. Figure adapted from [11] under the Creative Commons Attribution 4.0 International Deed.

Finally, we implemented immersive visualization in a cave automatic virtual environment system at the National Institute of Standards and Technology  to better understand cell interaction behaviors at the tumor’s center [12] (see Figure 3). This technique allowed us to zoom into local regions of the solid tumor to determine interactions that cause exhaustion scores that are higher or lower than expected. Through immersive visualization, we were able to identify unexpected behaviors—such as unique infiltration trajectories—that shaped both the simulation’s development and the potential for biological insights.

<strong>Figure 3.</strong> Immersive rendering of the virtual tumor microenvironment through the cave automatic virtual environment system at the National Institute of Standards and Technology. Figure adapted from [12].
Figure 3. Immersive rendering of the virtual tumor microenvironment through the cave automatic virtual environment system at the National Institute of Standards and Technology. Figure adapted from [12].

In summary, we determined a predicted threshold for blockade efficiency that can sufficiently reduce exhaustion and falls within the range of blockade efficiency values that are currently under investigation for the standard-of-care strategy that inhibits programmed cell death protein 1 (anti-PD-1):  an immune checkpoint receptor on T-cells that serves as a checkpoint for immune system regulation [2]. We also identified a relationship between immune checkpoint blockade efficiency and the T-cell activity that is responsible for determining the fate of a simulated tumor (either progressive disease or remission). This work lays an important foundation for future model iterations and subsequent insights in immuno-oncology that can account for longer time scales and associated dosing schedules, alternative checkpoints and associated pathways, parameterization that is specific to standard-of-care therapeutics, scaled-up simulations in time and increased spatial resolution, and integrated complementary ODE and ABM model components. These and other modeling efforts will be invaluable during the optimization of cancer drug design and related development processes.

Anne Talkington delivered a minisymposium presentation on this research at the 2025 SIAM Conference on Computational Science and Engineering, which took place last year in Fort Worth, Texas. She received funding to attend CSE25 through a SIAM Early Career Travel Award. To learn more about Early Career Travel Awards and submit an application, visit the online page.

References 
[1] He, X., & Xu, C. (2020). Immune checkpoint signaling and cancer immunotherapy. Cell Res., 30(8), 660-669. 
[2] Kuah, C.Y., Monfries, R., Quartagno, M., Seckl, M.J., & Ghorani, E. (2023). What is the optimal duration, dose and frequency for anti-PD-1 therapy of non-small cell lung cancer? Ther. Adv. Med. Oncol., 15, 17588359231210271. 
[3] Kuznetsov, V.A., Makalkin, I.A., Taylor, M.A., & Perelson, A.S. (1994). Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis. Bull. Math. Biol., 56(2), 295-321. 
[4] Man, J., Millican, J., Mulvey, A., Gebski, V., & Hui, R. (2021). Response rate and survival at key timepoints with PD-1 blockade vs chemotherapy in PD-L1 subgroups: Meta-analysis of metastatic NSCLC trials. JNCI Cancer Spectr., 5(3), pkab012. 
[5] Merck & Co., Inc. (2014). Highlights of prescribing information (Keytruda). Retrieved from https://www.merck.com/product/usa/pi_circulars/k/keytruda/keytruda_pi.pdf.  
[6] Pardoll, D.M. (2012). The blockade of immune checkpoints in cancer immunotherapy. Nat. Rev. Cancer, 12(4), 252-264.
[7] Postow, M.A., Callahan, M.K., & Wolchok, J.D. (2015). Immune checkpoint blockade in cancer therapy. J Clin. Oncol., 33(17), 1974-1982. 
[8] Rudloff, M.W., Zumbo, P., Favret, N.R., Roetman, J.J., Detrés Román, C.R., Erwin, M.M., … Philip, M. (2023). Hallmarks of CD8+ T cell dysfunction are established within hours of tumor antigen encounter before cell division. Nat. Immunol., 24(9), 1527-1539. 
[9] Sun, J.-Y., Zhang, D., Wu, S., Xu, M., Zhou, X., Lu, X.-J., & Ji, J. (2020). Resistance to PD-1/PD-L1 blockade cancer immunotherapy: Mechanisms, predictive factors, and future perspectives. Biomark. Res., 8, 35. 
[10] Talkington, A., Dantoin, C., & Durrett, R. (2018). Ordinary differential equation models for adoptive immunotherapy. Bull. Math. Biol., 80(5), 1059-1083. 
[11] Talkington, A.M., & Kearsley, A.J. (2025). Optimization of immune checkpoint blockade via a multiscale model system. Comput. Syst. Oncol., 5(2), e70007. 
[12] Talkington, A.M., Satterfield, S.G., & Kearsley, A.J. (2026). Immersive visualization of a 3D model of tumor-immune dynamics. In M. Golubitsky, S. Boccaletti, & C.M.A. Pinto (Eds.), Mathematical approaches to challenges in biology and biomedicine (ICMASC 2024), Springer proceedings in mathematics & statistics (Vol. 507). Porto, Portugal: Springer.
[13] Zhang, T., Xie, J., Arai, S., Wang, L., Shi, X., Shi, N., … Zhang, Y. (2016). The efficacy and safety of anti-PD-1/PD-L1 antibodies for treatment of advanced or refractory cancers: A meta-analysis. Oncotarget., 7(45), 73068-73079. 

About the Authors

Anne M. Talkington

Assistant professor, University at Buffalo

Anne M. Talkington is an assistant professor at the University at Buffalo Department of Pharmaceutical Sciences, where she works at the intersection of mathematical modeling, immunology, oncology, and pharmacology. Her research group integrates modeling and simulation with experimental data to predict optimal drug design and treatment strategies.

Anthony J. Kearsley

Staff research mathematician, National Institute of Standards and Technology

Anthony J. Kearsley is a staff research mathematician in the Applied and Computational Mathematics Division at the National Institute of Standards and Technology.