Immunotherapy of cancer has had a long history of development, starting from pioneering efforts in using coley toxins to treat patients—a therapeutic approach named after Dr. William Coley (1). Even though these earlier efforts never turned into a standard treatment, further investigations on the relationships between tumor cells and the immune system led to discoveries which unveiled fundamental principles underlying cancer progression, such as immune surveillance (2, 3), cancer dormancy (4), cancer immuno-editing (5), and the cancer immunity cycle (6). These discoveries were foundational for clinical successes and corresponding regulatory approvals in recent years, of therapies targeting the CTLA-4, PD-1, and PD-L1 immune checkpoints. In the wake of these successes, there has been an explosion in the development of immuno-modulating, anti-cancer pharmacological modalities, leading to the initiation of, literally, thousands of clinical trials (7, 8). However, from the discovery phase to late-stage clinical testing and regulatory approval, challenges in the development of immuno-oncology (IO) drugs are multi-fold and complex (9), with related complexities in the design of clinical trials; if unaddressed, these may lead to a decreased probability of success (10). Some of these challenges can be mapped to an incomplete mechanistic understanding of immune response dynamics and the interplay of such immune responses with tumor infiltration processes and tumor cell growth (11). These quantitative knowledge gaps hinder: (i) effective translation of novel promising therapeutic approaches into the clinic, (ii) identification of predictive response biomarkers, and (iii) search of therapeutic drug combinations which may overcome intrinsic or acquired resistance to existing standards of care (12). This mini-review focuses on quantitative, mechanistically-oriented modeling approaches which have been sought in IO, to address, at least partially, the abovementioned challenges and knowledge gaps.

Application of mathematical modeling in support of preclinical and clinical research, as well as decision-making in Oncology, has a long-standing history covering multiple problems and addressing a variety of research questions—today often referred to as computational oncology (13–15). Historical milestones include adaptations of the Gompertz model for treatment outcomes in breast cancer (16). These earlier efforts started from models with a simplistic empirical structure, based on an ordinary differential equation (ODE) describing tumor size growth using an exponential or sigmoidal function (17). Such a model, however, would not adequately describe the interplay between tumor cells and tissue vs. the immune system, since it entirely ignores the immune component (18). It is nevertheless valuable to mathematically describe treatment response effects following various chemotherapies, which are adequately captured by generalized Gompertzian kinetics (19). In fact, such modeling results provided a basis for the use of specific “dose-dense” chemotherapeutic regimens, which subsequently showed favorable outcomes in the treatment of breast cancer (20). Additionally, such empirical considerations allowed for a gradual evolution of modeling concepts, which today can be grounded in mechanistically-oriented principles, including for tumor vs. immune system interactions (Figure 1).

Earlier efforts to describe tumor vs. immune system relationships via a general mathematical description appeared in the 1980's, following the pioneering IO work that introduced the concept of immune surveillance (2, 3). These mathematical models considered the addition of a second variable describing the dynamics of cytotoxic immune cells, which are able to attack tumor cells (22–24). The resultant “two-ODE” model actually follows a typical “predator-prey” model introduced by Alfred Lotka and Vito Volterra, in much earlier days, at the turn of the 20th century. In such a model, tumor cells may be interpreted as the “prey,” whereas cytotoxic immune cells may be viewed as the “predator”: their dynamic interplay may result in one possible system behavior reflective of cancer dormancy (4). Given the relative simplicity of such a “two-ODE” model and since the behavior of such a model could be assessed analytically, it gained immense popularity within the oncology modeling community and led to several theoretical hypotheses underlying fundamental principles of cancer progression. For example, it was shown, through modeling, that key parameters controlling tumor re-growth under steady-state conditions of cancer dormancy were those relating to activities of the immune system (25). A corollary result was that it is a reduction in the probability of achieving tumor cell kill, rather than a reduction in the probability of tumor cells being recognized by cytotoxic cells, which best explained immune evasion by tumor cells (26). Interestingly, this key result, derived theoretically at the time, has recently been supported by elegant modeling work linking high-level immunological and epidemiological data, which suggests that age-related decline in T cell output correlates better with risk of cancer diagnosis vs. age-related accumulation of somatic mutations in tumor cells (27).

With the explosive growth of experimental data surrounding the complexity of tumor vs. immune system interplay, “two-ODE” models experienced a further evolution with additional biological entities and mechanisms being taken into mathematical consideration. At this point and looking forward, many biological candidates were tested as the “third modeling variable,” representing either specific immune cells or cytokines that modulate cytotoxic T lymphocyte (CTL) function (28). Such models were initially focused on including IL-2 function and effects, reflective of the potential importance of this cytokine and its associated dynamics in long-term tumor relapse (29). In further work, de Pillis et al. used a “three-ODE” model to reveal a difference between the dynamics of CD8⁺ CTLs vs. natural killer cells, which supported the importance of considering multiple cell types in the overall anti-tumor immune activity (30). More recently, CD4⁺ T helper cells were considered as the third component, in a quantitative, model-based investigation of adoptive cellular immunotherapy (31).

“Three-ODE” models, however, exhibit one significant structural limitation, namely they completely lack (an) immuno-suppressive component(s), which would be crucial when considering immune evasion mechanisms (32). Therefore, embedding a fourth variable into such models, to describe immuno-suppression, would seem rather natural; however, choices for the most appropriate candidate in this role are multi-fold. Several types of immuno-suppressive cells or molecules could be suitable candidates, including regulatory T cells (Tregs), myeloid-derived suppressor cells (MDSCs), or Type 2 tumor-associated macrophages, as well as cytokines such as TGFβ or IL-10. Thus, Arciero et al. chose TGFβ as the fourth model variable (33), while de Pillis et al. used Tregs as the principal immuno-suppressive component in their model (34). While these two modeling examples focused on immuno-suppressive effectors, other “four-ODE” models abound, declining a vast variety of immune players or tumor cell clones (35–41). Such a variety in potential key immuno-modulating factors made the generalization of any “three-ODE” or “four-ODE” model an overly difficult process, since any one of the models cited above can be challenged with newly generated experimental data featuring the importance of one vs. another immune factor. This may also explain, at least partially, the relatively minimal recognition, to date, of quantitative modeling approaches by immuno-oncologists (28, 42, 43).

On one hand, some of the biological complexities which compose the IO cycle, as summarized in recent reviews (6, 44, 45) clearly indicate the limitations of oversimplified models such as “prey-predator” models, which appear to be too remote from experimental reality and would not be applicable or of use for the majority of research relevant questions. On the other hand, increasing model complexity with additional mechanistic insights always comes with challenges of model calibration, as depicted in this famous quote by John von Neumann, “with four parameters, I can fit an elephant and with five, I can make him wiggle his trunk” [see in Dyson (46)]—pointing to the necessity of avoiding overparameterized “metastatic” models with unreliable extensions and loss of predictive power. Achieving such a balance in capturing necessary (not over-simplified) yet sufficient (not over-developed) features, and as constrained by the available data, is arguably one of the most difficult challenges in fit-for-purpose, parsimonious mechanistic model building and calibration. Overparameterization can easily negate all benefits brought forward by the incorporation of exquisite biological details of the system under consideration (47); models which attempt to explain everything may in fact not be useful, their predictive power remaining a question mark (48).

To address this challenge, part of the solution may reside in the combining of modeling methodologies developed previously and in other disciplines (49). This would result in a repository of prior information and knowledge validated elsewhere, to build mechanistic models in immuno-oncology which, on one hand, incorporate increasing system complexity and, on the other hand, avoid overparameterization based on newly generated data—thereby resulting, using terminology of a Bayesian mindset, in a posterior model based on existing, established quantitative priors (50).

If so, the question then becomes, “where to find such established prior models?” One obvious domain is quantitative immunology (51, 52), where the use of various modeling techniques by experimentalists has already gotten significantly more traction, arguably, than in other fundamental biological disciplines (53). For example, modeling has provided quantitative “inference frameworks” for immunology basics and fundamentals such as T cell activation, homeostasis or self / non-self recognition (54–57), immune receptor signaling (58), and understanding of T cell immunological memory (59). Prior models from quantitative immunology may then be combined with prior models from quantitative pharmacology (60–62), another field where modeling has provided quantitative “inference frameworks” (63) In the next section, we will discuss selected works which considered the combining of modeling methodologies, in attempts to develop pharmacologically-modulated posterior models, which were then used to prospectively address questions in the development of IO therapies (Table 1).

Applications of mechanistic modeling in support of preclinical and clinical research, commonly referred to as pharmacokinetic (PK)/pharmacodynamic (PD) modeling, are traditionally centered around the optimization of treatment dosing and scheduling—the “dose” representing a critical component of any drug development program (82). Such modeling approaches have thus been used in the development of IO agents such as PD-1 and PD-L1 inhibitors (83–85). In particular, mechanistic PKPD modeling has been applied in support of first-in-human dose selection of pembrolizumab, an anti PD-1 agent (65); this resulted in a seamless clinical trial design with a model-informed dose justification, which the US FDA accepted in the process of an accelerated regulatory review (86). Label updates with flat dosing schedules were subsequently granted, for both nivolumab and pembrolizumab, strongly supported by model-based simulations (87, 88). PKPD modeling has also been used in the translation of preclinical data for a conjugated IL-2 therapy, in particular to gain a better understanding of such a therapy's downstream effects (89). PKPD modeling has been further used in the development of bispecific biologics. Chen et al. used it for the estimation of the minimally anticipated biological effect level (MABEL) of a bispecific antibody targeting CD3 and p-cadherin (66), while Ribba et al. used it for guided dose escalation study design of cergutuzumab amunaleukin, a fusion protein consisting of IL-2 and a carcinoembryonic antigen (CEA) human monoclonal antibody (64). Such models are great examples of a “fit-for-purpose” quantitative approach, focused on addressing a specific pharmacological question. However, they do not take into account details of the tumor vs. immune system interactions, which would be critical to gain a better understanding of mechanisms of action (MoA) of immunotherapies.

Progressively adding components of tumor vs. immune system interactions into such PKPD models may well support the addressing of questions around pharmacologically-modulated IO biology, a topic of paramount importance in, for example, the search for therapeutic IO drug combinations (90). Such a systems approach may become an indispensable quantitative tool supporting “go/no-go” decisions in development programs, especially if sufficient biological knowledge for viable generalization is considered in the model (91). This prior knowledge is generally derived from two sources: (i) connectivity information to determine the system structure, e.g., molecular & cellular interactions, and their integration into patho-physiological processes; and (ii) quantitative data, for the calibration of model parameters. As discussed in the previous section, an imbalance in structural vs. quantitative information will in one way or another complicate integration into, and practical use of a mathematical model. For example, Lai and Friedman developed an elegant, yet complex model which includes a high number of biological elements, and considered their dynamics in space and time using partial differential equations (PDEs), to better understand the potential synergy between PD-(L)1 antagonists and either a GVAX vaccination or BRAFi/MEKi targeted therapies (72, 73). However, assessing the predictive power of such a model is impractical, given insufficient experimental data for model validation. Serre et al. provided another example of an elegant, yet insufficiently validated mathematical model describing the potential synergy between radiotherapy (RT) and immune checkpoint blockade (70).

One obvious way to improve model validation and hence model predictive power is to use rich experimental data, to rigorously constrain model parameters. This, however, requires the use of adequate statistical methods to properly quantify uncertainty and variability, which are inherent to any experimental biomedical and life sciences dataset (49, 92). In oncology drug development, quantitative data supporting MoA elucidation are typically generated at the preclinical stage. Parra-Guillen et al. for example, used a nonlinear mixed-effects (NLME) model and experimental data from syngeneic tumor models, to reveal the most influential immuno-adjuvant capable of boosting anti-tumor vaccination effects (21, 67). Such a modeling approach, which combines mechanistic features and mixed effects, allows one to incorporate individual-level data into the model, which may then describe not only mean trends, but also the full range of individual biomarker dynamics (93). A similar, combined mechanistic and mixed-effects approach was used to develop a model describing synergistic effects between RT and PD-(L)1 blockade in mice (68). This model, in fact, synthesizes a fit-for-purpose, yet sufficiently detailed mathematical description of the IO cycle, together with adequate model validation based on data from multiple experiments. As a result, this model can be used as a simulation tool for experimental study design, and is also adequate for determining optimal schedule and sequencing of RT + IO, and IO + IO treatment combinations (68, 69). Interestingly, despite the well-known challenges in translating oncology preclinical results into the clinic, simulation results from this preclinical modeling exercise were recently supported, in a qualitative sense, with clinical data and a corresponding meta-analysis (94, 95). For a quantitative translation, the Kosinsky et al. model would require adjustments for multiple quantitative differences that exist between mouse vs. human immune systems, e.g., appropriate expressions of immune checkpoints and turnover of specific T cells (96). Another modeling approach aimed at supporting the development of such an RT + IO combination therapy was proposed by Poleszczuk et al. who developed a physiologically-based model which considered a detailed incorporation of T cell trafficking and was used for the identification of an optimal site for RT administration, to maximally increase the probability of incremental anti-tumor immune effects (71). Predictions from such a comprehensive modeling effort were also recently supported by clinical results, which showed that RT administered to liver metastases triggered a higher immunological response (97). A mechanistic model has also been proposed by Peng et al. in the search of an optimal combination strategy against castration-resistant prostate cancer (74).

The modeling applications discussed to this point emphasize the importance of addressing multi-pronged questions, e.g., not only around dose finding, but also on the identification of an adequate time window for maximizing therapeutic benefits (98). This problem is particularly challenging in the development of combination therapies, where multiple options around which cancer indication, which combination agents, which scheduling per agent, and which sequencing of the agents make trial design enormously complex (99, 100). In recent years, platform design of clinical studies, driven by one master protocol, has gained momentum (101, 102)—a format which, in fact, benefits even further from a supportive quantitative mechanistic modeling approach (103).

A third problem which is highly relevant in the development of IO therapies is the identification of predictive biomarkers. Indeed, there still is a lot of room for improving numbers of responder patients in pivotal IO trials, even in immunologically-active indications (104). Several computational models focusing on the identification of predictive biomarkers, with applications to personalized treatment against glioblastoma and prostate cancer have been developed (75, 77). These approaches have yet to find a general use in clinical practice. Part of the challenge arises from the biological complexity in the IO field, although there also are significant limitations from an experimental standpoint, such as differences in fresh vs. archived samples, difficulties in obtaining multiple biopsies per patient, with related risk and cost issues (105). One approach to alleviate some of these problems is the development of novel combinatorial biomarkers (“signatures”) which may relate multiple, routinely measured markers with clinically meaningful biological phenotypes (106). In fact, such a consensus approach, “Immunoscore,” has recently been validated in a large international study of colon cancer (107).

Another complicating factor in the development and interpretation of mechanistic modeling of IO data is the tremendous heterogeneity in tumor cell clones and elements of the surrounding immune microenvironment (108). A rapid development of novel experimental techniques may overcome this challenge, at least partially. Thus, the identification of specific gene expression signatures may help in further validating existing immunoscores and related biomarkers, even increasing their discriminatory ability (109), as recently shown with a PD-L1 expression signature which outperformed a standard PD-L1 immunohistochemistry (IHC) assay (110). Multiple immune signatures have now been identified, which allow for a better characterization of various aspects of anti-tumor immunity (111–116). Recent technological breakthroughs such as cytometry by time-of-flight (CyTOF) and single-cell mRNA sequencing (scRNA-seq) may further advance the utility and robustness of these immune signatures (117, 118); these techniques may allow for a deeper, more granular profiling of tumor and immune cell phenotypes involved in response or resistance to immunotherapies, in multiple indications (119–123). The importance in using quantitative models toward the selection and qualification (within the chain of events, from dosing to patient response) of IO biomarker signatures cannot be over-emphasized (108): immune biomarkers involve a high number of molecular and cellular species, and often exhibit complex temporal and spatial dynamics; these need to be properly framed in the context of a quantitative model, especially if the purpose is to relate multi-variate biomarker signatures to IO treatment effects and clinical endpoints (124). Quantitative modeling may also support the development of biomarkers in context, by integrating different data types, and following a model-based qualification of biomarkers as surrogate measures of efficacy and response. Such an approach has been proposed, recently, in the evaluation of neoantigen fitness as a surrogate measure of immunogenic quality of the existing neoantigen pool (125, 126). The progressive integration of such consensus, multi-variate combinatorial biomarkers into a unified, quantitative and mechanistic modeling framework will help overcome some of the limitations in the clinical use of IO biomarkers (127, 128).

The above sections focused on traditional deterministic models, which make use of ODEs and PDEs for the description of IO systems dynamics. Other modeling techniques can be used to describe tumor vs. immune interactions. For example, cellular automata and agent-based models (ABMs) (129), as well as various hybrid models which link continuous and discrete modeling elements have been developed (130). Such models may be useful in raising new hypotheses, which may arise from emergent properties of the system based on existing data, rather than generating bona fide forward predictions. For example, a lattice gas automata technique has been used to gain a better understanding of a vaccination treatment mechanism and its corresponding anti-tumor immune response dynamics (131, 132). ABMs also represent a popular modeling technique, since they are well-suited to describe stochastic processes which do occur at various stages of the IO cycle. For example, Gong et al. developed an ABM to reveal spatio-temporal characteristics of PD-L1 blockade (79). In another publication, Kather et al. presented an elegant 2D ABM framework for an improved understanding of the role of stromal cells in colorectal cancer (CRC) (80). These authors determined that malignant cells hiding in the stroma cannot be eradicated completely, while stromal cells, at the same time, would not allow for rapid tumor progression. Consequently, simulations of an immuno-therapy illustrated how stroma permeabilization, concomitantly with immune activation, were able to markedly increase response to therapy in silico. Additionally, it was shown that a stroma-targeted therapy with insufficient activation of tumor-specific CTLs can lead to rapid tumor escape and hyper-progression (80). More recently, this model has been extended and generalized to a 3D spatial description, incorporating macrophage effects; it accurately reproduced the tissue architecture typically observed in CRC and can be used, similarly to ODE systems models, for the identification of effective IO therapeutic combinations (81).

Following the approval, in recent years, of the first immune checkpoint inhibitors, the landscape of cancer treatment has changed dramatically and has shifted to a deep reconsideration of the role of the immune system in cancer progression and treatment. This led to an unprecedented number of clinical trials and generation of clinical data in the IO field. Clinical success rates, however, while improving significantly, are still relatively low. The observed imbalance, between the amount of biological and clinical data being generated vs. probability of trial success is not uncommon in biomedical disciplines, and calls for the development and updating of a companion, integrative, quantitative modeling framework with predictive value for MoAs and simulation value for study design purposes. As described by Sidney Brenner in his “Sequences and Consequences” landmark paper: “We should welcome with open arms everything that modern technology has to offer us but we must learn to use it in new ways. Biology urgently needs a theoretical basis to unify it and it is only theory that will allow us to convert data to knowledge” (133). We propose that quantitative, mechanistically-oriented modeling represents a means toward the establishment of such a “theoretical basis,” pending proper integration of prior knowledge gained from biology and clinical research. One of the main factors limiting a wider application of quantitative systems modeling is its demand for rich experimental data necessary for precise parameter estimation. Historically, generation of such datasets in oncology research has been challenging, due to translational limitations of experimental preclinical models and sparse collection of tissue samples in clinical settings. Also, in the IO field, another challenge is the lack of predictive power for univariate biomarkers (e.g., PD-L1 IHC status or tumor mutational burden taken in isolation), which may unequivocally link immunologically-driven therapeutic effects to clinical response; a multi-variate approach is clearly needed (128). Recent developments in multi-modality biomarkers and associated molecular signatures, together with innovative pharmacologies and clinical design under platform trials (134) will help in the progressive build-out and qualification of such a unified quantitative modeling framework, which in turn may help in predicting patient responses based on a given pharmacological intervention choice and multi-variate biomarker signatures.

KP, VV, and YK generated Figure 1. The Table was generated by IA, VV, and LC. All authors contributed to the writing of the manuscript.