QSP‐IO: A Quantitative Systems Pharmacology Toolbox for Mechanistic Multiscale Modeling for Immuno‐Oncology Applications

Abstract

Immunotherapy has shown great potential in the treatment of cancer; however, only a fraction of patients respond to treatment, and many experience autoimmune‐related side effects. The pharmaceutical industry has relied on mathematical models to study the behavior of candidate drugs and more recently, complex, whole‐body, quantitative systems pharmacology (QSP) models have become increasingly popular for discovery and development. QSP modeling has the potential to discover novel predictive biomarkers as well as test the efficacy of treatment plans and combination therapies through virtual clinical trials. In this work, we present a QSP modeling platform for immuno‐oncology (IO) that incorporates detailed mechanisms for important immune interactions. This modular platform allows for the construction of QSP models of IO with varying degrees of complexity based on the research questions. Finally, we demonstrate the use of the platform through two example applications of immune checkpoint therapy.

During the past decade, the development of immunotherapy has become one of the most exciting breakthroughs in cancer research. Extensive preclinical and clinical research efforts that investigated the relationship between immune response and tumor growth have demonstrated that our body's immune system, upon proper education and activation, does have the full potential to completely eradicate malignant tumors.^1,2 Motivated by these discoveries, numerous therapeutics have been designed to modulate the activity of major immune cells toward enhanced recognition and killing of cancer cells through distinct mechanisms, the most successful being the immune checkpoint blockers, namely, antibodies against programmed cell death protein 1 (PD‐1), programmed death ligand 1(PD‐L1) and cytotoxic T‐lymphocyte‐associated protein 4 (CTLA‐4) pathways, which have now revolutionized the treatment standards for multiple cancer indications.^3–5

Besides these three commonly targeted checkpoints, numerous immune‐targeting drugs have been widely investigated in clinical trials against cancers. Examples include drugs that target agents such as indoleamine 2,3‐dioxygenase and interleukins⁶ or other immune cells (e.g., myeloid cells, dendritic cells).^7,8 In addition, various treatment combinations of immunotherapies with other anticancer therapeutics (e.g., chemotherapy, antiangiogenic therapies, targeted tyrosine kinase inhibitors) have also been investigated.^9–11 However, the large number of registered trials is by no means an indicator of clinical success (more than 3,000 immunotherapy trials in the United States alone, and more than 1,700 active trials testing PD‐1/PD‐L1 blockers with other anticancer therapies^12,13). In fact, the general response rates of immunotherapies (e.g., US Food and Drug Administration–approved PD‐1/PD‐L1 antibodies) in cancer patients are highly variable, and for many cancer types only a small subset of patients would respond favorably even in biomarker‐selected cohorts.^10,14 Therefore, it is of significant translational value to develop an integrative mechanistic understanding of cancer biology, interpatient variability, and drug–target interactions using state‐of‐the‐art, quantitative systems–level models to help answer critical questions throughout the immuno‐oncology (IO) drug development life cycle (e.g., target selection, biomarker identification, trial optimization), with the ultimate goal to reduce clinical trial failures and improve treatment efficacies in patients.¹⁵

In the pharmaceutical industry, the concept of modeling and simulating biological processes has long been implemented to investigate the pharmacokinetic (PK) and pharmacodynamic (PD) properties of candidate drugs. In recent years, the emergence of complex mechanism‐based, quantitative systems pharmacology (QSP) models has attracted increasing interest in translational research and has now become an essential component in the practice of model‐informed drug discovery and development.^16–21 Centered around the mechanisms of human physiology and pathogenesis, endogenous cellular dynamics and pharmacology, and disease‐associated gene/protein signatures, which can all be translated quantitatively into patient‐specific and disease‐specific parameters to generate virtual patient cohorts and simulate trial results, these fit‐for‐purpose or platform‐based QSP models have already proven to be useful throughout the drug development process in a wide range of therapeutic areas.^18,22–25

In the field of IO, a number of computational models that aligned closely with the concepts of QSP with varying degrees of physiological details about the immune system have been developed to characterize the complex cellular dynamics that regulate the antitumor immune response elicited by pharmacological agents (many of the fit‐for‐purpose IO models are reviewed in ref. 15). Regarding the recent developments in platform‐based IO models, Wang et al. constructed a detailed compartmentalized systems pharmacology model that considers both molecular‐level (e.g., multiple ligand/receptor binding pairs) and cellular‐level kinetics (e.g., the chain of events in immune‐cancer interactions) to simulate personalized treatment outcomes in response to checkpoint inhibitors for patients with metastatic breast cancer.²⁶ Milberg et al., based on the similar platform backbone, expanded the reaction dynamics of cell–cell interactions and reparameterized the new model using clinical trial data in melanoma; the authors then used their model to assess different treatment regimens of checkpoint inhibitors for simulated melanoma patient cohorts.²⁷ Jafarnejad et al. further enriched the platform with a mechanistic module of antigen processing and presentation; the extended model was then parameterized to investigate potential biomarkers of tumor regression and predict patient‐specific response to neo‐adjuvant anti‐PD‐1 therapy in non‐small cell lung cancer using clinical data.²⁸ This platform has also been used to model combination therapies in triple negative breast cancer²⁹ and colorectal cancer.³⁰ These QSP modeling studies from our group suggested the potential translational values and efficiency of platform‐based model formulation in the in silico investigation of complex human diseases, especially in IO where new mechanisms are being discovered and novel treatment modalities are being tested at a very rapid pace.

In this tutorial, we present QSP‐IO, a modularized QSP modeling platform based in MATLAB (MathWorks, Natick, MA) that allows efficient formulation of physiology‐based IO models and quantitative simulation of patient outcomes in response to different immunotherapy combinations. The equations implemented in QSP‐IO are based on those implemented in Jafarnejad et al. with minor improvements outlined in the Modules section. This tutorial assumes the reader is able to run and edit scripts in MATLAB and has a basic understanding of the SimBiology toolbox. This tutorial does not require any specialized knowledge of QSP, although the interested reader is referred to the following articles: refs. 22,31,32

We first give an outline the platform, detailing the overall structure and workflow of modeling with QSP‐IO, followed by a description of the biological rationale and mathematical representation of the different mechanistic modules within the integrative QSP‐IO model platform. Then, through two examples that characterize the complex interplay between immune cells, cancer cells, checkpoints, antigen presentation, and targeted immunotherapy, we demonstrate how to construct mechanistic IO models from scratch by selecting and connecting modules in QSP‐IO, and how such models, based on physiological and clinical constraints, are mathematically initialized and simulated to predict patient response. Finally, we discuss how QSP‐IO as a modeling platform can be further enriched to provide insights into rising issues in IO from the computational perspective.

SOFTWARE PACKAGE OVERVIEW

This work presents a modular toolbox for assembling IO‐based QSP models in MATLAB using SimBiology (MATLAB 2018b, MathWorks). The code can be downloaded from our GitHub page (www.github.com/popellab/qspio). Overall, our QSP platform is a modular, multicompartmental dynamical systems model of IO that aims to predict various biological features important to antitumor immune physiology; similar models were implemented in our recent work.^26–28 QSP‐IO allows the user to create a SimBiology model object in a MATLAB script by calling functions from our platform. This section describes the directory structure (Directory Structure section), the parameter structures (Parameter Structure and Input section), our initial conditions procedure (Initial Conditions section), and the overall workflow (Workflow section).

Directory structure

The repository is organized into the following five main directories: @struct, model, parameters, scripts, and utils. The model directory contains all the functions necessary to create and manipulate the SimBiology model object. The parameters directory contains the parameter files (see the Parameter Structure and Input section for more details), and the scripts directory contains the user‐defined scripts for creating models. The utils directory contains useful functions for accessing and manipulating components of the model, and the @struct directory contains overloaded MATLAB functions (mpower, mrdivide, mtimes, plus, and simplify) to manipulate the parameter structure (see the Parameter Structure and Input section); note that the @struct directory must be added to the MATLAB path.

Parameter structure and input

Model parameters are stored as MATLAB structures that have the following three attributes: Value, Units, and Notes. Value and Units store the parameter's value and units, respectively, and the Notes attribute stores comments about that parameter that will get incorporated into the SimBiology object.

The model parameters are inputted into the model using a JavaScript Object Notation (JSON) file that can be edited with any plain text editor. All parameters have the fields of name, value, units, description, and source. The name field refers to the name of the parameter used by QSP‐IO; a complete list of the parameter names required by each module is included in the supplemental material. The value, units, and description fields correspond to the Value, Units, and Notes properties of the model parameter structures in MATLAB, respectively. The source field is used to indicate the reference where the parameter value was found. If the model parameter is determined by other parameters, the source field can be set to "derived"; these parameters must have the fields derived_from and expression. The derived_from field is an array of strings naming the parameters on which it is dependent. The expression field is a string giving the relationship between the named parameters in the derived_from property. An example for two different parameters is given next. Notice that in the expression field on line 15, p(1) and p(2) correspond to the first and second entries of the derived_from array, respectively; the parameters number of lymph nodes and diameter of lymph node must also be defined in the same file.

Listing 1: Sample JSON parameter file entries

{ “name”:“V_C”, “value”:5, “units”:“liter”, “description”:“Central compartment volume”, “source”:“doi:10.1002/psp4.12040” }, { “name”:“V_LN”, “value”:null, “units”:“liter”, “description”:“Lymph Node compartment volume”, “source”:“derived”, “derived_from”:[“number of lymph nodes”,“diameter of lymph node”], “expression”:“p (1) *4/3* pi *(p (2)/2)^3” },

Initial conditions

One challenge when facing large systems of ordinary differential equations is coming up with a set of realistic initial conditions to represent the state of the system at the beginning of the simulation. It is standard practice to treat the initial conditions as unknown parameters; however, in large systems, this introduces a large number of parameters. Milberg et al. overcame this issue by initializing the tumor to the desired size and set all other state variables to zero.²⁷ Although this solution overcomes the difficulty of having many unknown parameters, it results in a large tumor growth rate initially attributed to the delay in immune activation. Jafarnejad et al. overcame the latter difficulty by fixing the tumor to the desired size, then let the other state variables evolve until they reach steady state, then initialized the system with those values.²⁸ One concern with existing methods is that they could result in a potentially unrealistic behavior, for example, a realization where the tumor decreases in size artificially in the absence of therapy.

In our platform, the number of cancer cells is initialized to one cell, all other state variables are initialized to zero, and the simulation is run for a long period of time. The initial conditions are taken at the time where the volume reaches the target volume (specified by the user as a diameter; the target volume is calculated by assuming a spherical tumor). This procedure is meant to mimic cancer biology, where a normal cell mutates and becomes cancerous. Furthermore, any parameterization that would not reach the user‐specified initial size would not represent a realistic cancer patient, thus these subjects are eliminated from further simulations.

Workflow

Briefly, the platform incorporates the description of cancer cell growth with a death rate dependent on the number of T cells present in the tumor compartment. The platform has detailed mechanisms for immune activation by modeling the release of antigen from dying cancer cells; antigen uptake by antigen presenting cells; the presentation of antigen to naïve T cells; and the activation, proliferation, and transport of T cells between the four compartments. QSP‐IO also implements PK and PD models for immune checkpoint inhibitors of PD‐1 and PD‐L1. A summary of the model interactions is shown in Figure 1a. A detailed description of each module is presented in the following section.

1 Figure. (a) Diagram of immuno‐oncology quantitative systems pharmacology (QSP) model interactions. Naïve and mature antigen presenting cells are represented by APC and mAPC, respectively, naïve, proliferating and mature T cells are represented by nT, aT, Tcyt, respectively and regulatory T cells are represented by Treg. (b) Workflow for creating a model in QSP‐IO. The user starts by creating a SimBiology object using the QSP‐IO's initialization function. Module can be added to include the necessary detail in the model based on the research questions. The initial conditions are then generated based on the model parameters using the novel initial conditions procedure. Finally, the simulations can be run using standard SimBiology functions and visualization can be performed using MATLAB plot functions or using QSP‐IO's plotting function.

A summary of the workflow for using the toolbox is shown in Figure 1b. First, the workflow involves calling the initialization function that creates an empty SimBiology object with four compartments (called C, P, T, and LN, which stand for central, peripheral, tumor, and lymph node, respectively). The volumes of each compartment are constant, except for the tumor compartment whose volume is determined based on the number of cancer and T cells in the tumor compartment (see the Cancer Module section for more details). Second, the modules are called, which adds the relevant equations to the SimBiology object. After the model is created, the initial conditions routine is run, which creates the initial conditions for all the species in the model. Finally, the simulation can be run with the defined treatment scheme. The toolbox also includes routines for visualizing the results of the simulation. Overall, our toolbox consists of seven modules that can be called independently to control the complexity of the model to suit the specific research questions.

MODULES

The QSP‐IO platform currently consists of the following seven modules that may be added independently to build a QSP model: A cancer module, an antigen presentation module, an antigen module, a T cell module, a regulatory T cell (Treg) module, a checkpoint module, and a pharmacokinetics module. Only the cancer module is required, and all other modules are optional and can be omitted to decrease the complexity of the model. Because of the modular nature of QSP‐IO, additional modules can be created by following the structure of the existing modules; more details on extending the model are provided in the Conclusions.

In the current version, the cancer, T cell, and antigen modules can be called multiple times to simulate different clones of cancer, T cells, and antigen, respectively. Having multiple cancer cell clones may be used, for example, to model a heterogeneous tumor having populations of cells with different growth rates, different sensitivities to T cell killing, or having different neo‐antigens. Defining multiple T cell clones may be used to model subpopulations of T cells that recognize different neo‐antigens. Similarly, the antigen module can be called multiple times to simulate different antigens (including the effect of having competing self‐peptides). In addition, the pharmacokinetics module can be called to simulate administration of multiple drugs. Table 1 gives a summary of the modules included in QSP‐IO at the time of this writing; the modules are listed in the suggested order in which they should be called.

1 TableSummary of modules

Module	Required	Number of calls	Dependencies
Cancer	Yes	Unlimited	None
T cell	No	Unlimited	Cancer
Treg	No	1	Cancer, T cell
Antigen presenting cell	No	1	Cancer
Antigen	No	Unlimited	Cancer, T cell
Immune checkpoint	No	1	Cancer, T cell

Cancer module

The number of cancer cells in the tumor, C, is given by the following:[Image Omitted. See PDF]

where the first term on the right represents the rate of cancer cell proliferation and the second term represents the rate of cancer cell death. Cancer cell proliferation is modeled by logistic growth as done in refs. 33 and 34 with maximal growth rate, $k_{\text{growth}}$, and carrying capacity, $C_{\max }$. Cancer cell death is assumed to be attributed to two mechanisms: death attributed to the innate immune system using first‐order kinetics with rate constant, $k_{\mathrm{innate}}$, and death attributed to cytotoxic T cells simulated by Michaelis‐Menten kinetics with maximal death rate, $k_{T_{\text{cell}}}$, and the number of T cells in the tumor compartment, $T$.$H_{\text{PD}1}$ is the fraction of T cell inhibition due to the PD‐1 checkpoint, calculated according to:[Image Omitted. See PDF]

where $X=\text{Asyn}\left[Y{Y}_1\right]$ (see the Immune Checkpoint Module section) is the number of PD1/PDL1 complexes in the immune synapse between T cells and cancer cells, $PD{1}_{50}$ is the number of complexes for half‐maximal T cell inhibition from the PD‐1 checkpoint, and n is the Hill coefficient.

When modeling multiple cancer clones, Eq. 1 is used to calculate the number of cancer cells for each clone and the logistic growth term changes as follows:[Image Omitted. See PDF]

where C_i is the i^th cancer clone and C_total is the total number of cancer cells. The maximal growth rate constant can be different for each cancer clone. Note that this assumes that all cancer clones are in competition for the same limiting resources.

Similarly, when modeling multiple T cell clones, the rate of cancer cell death from T cells becomes[Image Omitted. See PDF]

where $T_i$ is the i^th T cell clone and $T_{\mathrm{total}}$ is the total number of T cells.

The tumor compartment volume, $V_T$, is then calculated based on the number of cancer cells and T cells in the tumor compartment according to:[Image Omitted. See PDF]

where $V_{\mathrm{cancer}}$ and $V_{T_{\text{cell}}}$ are the volumes of single cancer cells and T cells, respectively. $V_{{\mathbf{T}}_{min}}$ is the minimum volume of the tumor compartment, included to avoid a zero tumor compartment volume. Although this term has very little effect on the solution, it becomes more important at low tumor volumes; the tumor volume will tend toward $V_{{\mathbf{T}}_{min}}$ as the number of cells in the tumor go to zero. Biologically, $V_{{\mathbf{T}}_{min}}$ may be interpreted as the volume of the extracellular matrix or interstitial space around the tumor when the volume is small. The volume of other immune cells such as natural killer (NK) cells, tumor associated macrophages, and cancer‐associated fibroflasts are not considered, but may be addressed in future work.

In Eq. 2, $C_{\mathrm{total}}$ and $T_{\mathrm{total}}$ include dead cancer and T cells, because the volume of the tumor would not change instantly following cell death. The dead cells are assumed to clear with first‐order kinetics with rate constant $k_{\mathrm{clear}}$. Including the T cells and dead cells in the volume calculation allows for the model to capture the pseudo‐progression phenomenon observed in clinical settings³⁵ where the tumor appears to be growing while the number of viable cancer cells are decreasing (see Figure 2). In this simulation, the number of viable cancer cells reach a maximum in 50 days, whereas the tumor reaches its maximum volume at 294 days; this apparent discrepancy is attributed to the presence of T cells and dead cancer/T cells in the tumor. Once the number of viable cancer cells is low, the T cells deplete from the tumor. The tumor volume begins to decrease when number of cancer cells are low enough that the rate of death is less than the rate of dead cell clearance.

2 Figure. (a) Tumor volume as a function of time. (b) Number of viable (blue) and dead (red) cancer cells as a function of time. (c) Cell count for viable (blue) and exhausted (red) T cells as a function of time. The parameters from example 1 were used with cell clearance rate (kclear) = 0.001 day−1 and rate of cancer death from T cells (kTcell) = 8.7 day−1 to demonstrate the pseudo‐progression phenomenon.

Antigen presenting cell module

The number of antigen presenting cells (APCs) are simulated in the tumor and the tumor‐draining lymph node compartments. In our platform, we consider both naïve APCs and mature APCs. The APCs are assumed to mature in the tumor compartment where they are exposed to maturation cytokines, c. Mature APCs can migrate from the tumor to the lymph node,^36–38 where they will present the antigen to activate naïve T cells.

The number of naïve APCs, $D_i$, are calculated in both compartments (i = T, LN) by assuming first‐order kinetics about a target density, ${\rho }_i^D$, with rate constant, $k_D$, according to:[Image Omitted. See PDF][Image Omitted. See PDF]

where the last term in Eq. 3a represents the rate of maturation; with maximal rate, $k_{\text{mat}}$, and Michaelis‐Menten constant, ${[c]}_{50}$.^39,40 The mature APCs, ${\widehat{D}}_i$, are described by the following,[Image Omitted. See PDF][Image Omitted. See PDF]

where $k_{\text{mig}}$ is the first‐order rate constant for the rate of migration to the lymph nodes from the tumor.

The maturation cytokines, $c$, are calculated according to:[Image Omitted. See PDF]

where the first term on the right represents the first‐order kinetics about ${[c]}_0$ with rate constant, $k_c$.³⁹ The second term represents the increased rate of cytokine release attributed to inflammation, measured by the rate of T cell killing given in Eq. 1 with the constant, $x_c$, representing the amount of cytokines released per cancer cell death because of T cells.

Antigen module

The antigen module models the release of antigen from lysed cancer cells, the uptake by mature APCs, degradation of antigen in the APC endosomes, and the presentation of antigen‐major histocompatibility complex (MHC) on the surface of APCs.^39–41 The equation governing the antigen concentration, ${[P]}_{\mathbf{T}}$, in the tumor is given by[Image Omitted. See PDF]

where the first term on the right is the number of antigen clones represented, $n_{\mathrm{clones}}$, times the concentration of a single antigen clone in a cancer cell, ${[P]}_0$, times the rate of cancer cell death. The second term on the right is the first‐order uptake of antigen by APCs, with rate constant $k_{\mathrm{up}}$ and the first‐order degradation in the extracellular space with rate constant, $k_{\mathrm{deg}}$.

The equations governing the concentration of antigen, ${\left[P\right]}_e$, and epitope, ${\left[p\right]}_e$, in the endosomes of APCs are given by[Image Omitted. See PDF][Image Omitted. See PDF]where $V_e$ and $A_e$ are the volume and surface area of the endosomal compartment, respectively; $k_{deg}^P$ and $k_{deg}^p$ are the rate constants for antigen and epitope degradation, respectively; and $k_{\text{on}}$ and $k_{\text{off}}$ are the rate constants for antigen‐MHC binding kinetics.

The following equations govern the binding of antigen epitope, $p$, to MHC, $M$, to form the antigen‐MHC complex, $Mp$, in the endosomes (subscript e) and on the surface of APCs (subscript s).[Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF]

where $A_e$ and $A_s$ are the surface areas of the endosomes and APCs, respectively, and $k_{\mathrm{in}}$ and $k_{\mathrm{out}}$ are the rate constants for MHC internalization and externalization, respectively. Note that in Eqs. 7 and 8, the concentrations of $M$ and $Mp$ are per unit area.

T cell module

In the T cell module, we model the activation of naïve CD8+ T cells, the proliferation of activated T cells and transport of mature, activated, cytotoxic T cells. This module is also used by the Treg module to model the activation, proliferation, and transport of regulatory T cells (see the Regulatory T Cells section).

Naïve T cells

The number of naïve T cells are modeled in the central, peripheral, and tumor‐draining lymph node compartments.^42,43 In the central compartment, naïve T cells, ${\mathcal{N}}_{\mathbf{C}}$, are modeled according to:[Image Omitted. See PDF]

where the first term on the right represents the rate of release of naïve T cells from the thymus described by zero‐order kinetics with rate constant, $\sigma$, the second term represents the rate of proliferation modeled by Michaelis‐Menten kinetics with maximum proliferation rate, $k_{\text{prolif}}$, and Michaelis constant, $K_m$. The third term on the right represents the transport between compartments using first‐order kinetics with rate constants $q_{\mathbf{P}}^{\text{in}}$, $q_{\mathbf{LN}}^{\text{in}}$, $q_{\mathbf{P}}^{\text{out}}$, and $q_{\mathbf{LN}}^{\text{out}}$ representing the rates for entry and exit of the peripheral and lymph node compartments. The last term on the right is the first‐order death term.

Similarly, naïve T cells in the peripheral and tumor‐draining lymph node compartments are given by the following:[Image Omitted. See PDF][Image Omitted. See PDF]

where $n_{\mathrm{sites}}$ is the maximum number of T cell binding sites on the APCs. The rate of activation is assumed to be dependent on the number of APCs in the tumor and the number of antigen‐MHC molecules presented on the surface of APCs. $H_p$ is the fraction of activation due to T cell receptor (TCR) binding; see the T Cell Activation section for details.

T cell activation

The activation of naïve T cells by mature APCs is modeled by considering the interaction between TCRs on naïve T cells and antigen‐MHC complexes on mature APCs. This engagement is modeled based on kinetic proofreading with limited signaling that was shown by Lever et al.⁴⁴ to best represent our understanding of T cell activation. Kinetic proofreading with limited signaling captures the presence of an optimum in T cell activation with respect to binding affinity of the TCR with antigen‐MHC complex, which is observed at all levels of antigen‐MHC.⁴⁴ The activation of the TCR is modeled by:[Image Omitted. See PDF][Image Omitted. See PDF]

where $M{p}_s$ is the total amount of antigen‐MHC complex on the surface of mature APCs, ${\mathrm{TCR}}_{\mathrm{tot}}$ is the total amount of TCR on naïve T cells, K_D is the binding affinity ($k_{\text{on}}^{\text{TCR}}/k_{\text{off}}^{\text{TCR}}$) of TCR binding to antigen‐MHC complex, $C_{\mathrm{tot}}$ is the total TCR‐antigen‐MHC complex, $k_p$ is the modification rate, $\varphi$ is the modification rate to the nonsignaling state, and $m$ is the number of intermediate states. ${\text{TCR}}_{\text{active}}$ was then related to the activation rate of naïve T cells by mature APCs using the following Hill function,[Image Omitted. See PDF]where $K_{p,50}$ is the half‐maximal level for T cell activation.

T cell proliferation

The number of first generation of proliferating T cells, $\mathcal{P}$, are calculated according to:[Image Omitted. See PDF]

where the first term represents the rate of activation and the second term represents the rate of proliferation.

Assuming the T cells have a division destiny of N and that proliferation is in quasi‐steady state, the number of proliferating T cells in the N^th generation is given by $2^N\mathcal{P}$. Recent work suggests that the division destiny is dependent on TCR binding, costimulatory signals, and the presence of interleukin‐2 (IL2) and that each of these contributions are additive⁴⁵ such that,[Image Omitted. See PDF][Image Omitted. See PDF]

In the present work, $N_{\text{TCR}}$ and $N_{\text{costim}}$ are assumed to be constant; however, future work will include variable $N_{\text{TCR}}$ dependent on $H_p$ from the kinetic proofreading model described in the T Cell Activation section.

IL2 is assumed to be secreted from proliferating T cells and is consumed by both proliferating regulatory T cells and proliferating cytotoxic T cells. The equation governing IL2 concentration is given by,[Image Omitted. See PDF]

where T_LN and $T_{\mathbf{LN}}^{\text{Treg}}$ are the number of cytotoxic and regulatory T cells in the lymph node compartment, respectively.

Activated T cell distribution

T cell transport between compartments is modeled by first‐order dynamics in the same way as with naïve T cells.⁴⁶ The rate of transport between the central and tumor compartments was assumed to scale with tumor volume. The equations governing the number of mature, activated, cytotoxic T cells in the four compartments are given by,[Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF]

where $k_{T\text{reg}}$ is the rate constant for T cell exhaustion by Tregs, $k_C$ is the rate constant for T cell exhaustion by cancer cells. As with the cancer module, the user may define multiple clones of cytotoxic T cells, in which case the equations will be introduced for each clone, and the new clones will be added to T_total.

Regulatory T cells

The Treg module models regulatory T cells in all four compartments using the same equations as the T cell module with minor differences. First, the induction of Tregs in the tumor‐draining lymph nodes are modeled using Eqs. 10b and 13 with the difference that Tregs are induced by presentation of self‐reactive peptides on immature APCs,⁴⁷ thus the activation term becomes:[Image Omitted. See PDF]

Second, the exhaustion terms from Tregs and cancer cells in Eq. 16c become zero. In addition, both natural and induced Tregs are accounted for in this module, but are both pooled together as Tregs. Natural Tregs are directly derived in the thymus and enter the blood.⁴⁸ The dynamics of the natural Tregs are modeled using an additional zero‐order term in the central compartment. Finally, the number of Tregs is added to the T cell killing of cancer cells term used in Eqs. 1, 5 and 6 as follows:[Image Omitted. See PDF]

Immune checkpoint module

Interactions between PD‐1 on T cell with PD‐L1 and PD‐L2 on cancer cell are modeled within an additional compartment that represents the immunological synapse between the two cells.⁴⁹ Different antibodies against these immune checkpoint molecules can be modeled. In this example, we have represented the currently approved bivalent antibodies against PD‐1 and PD‐L1. The equations governing the PD‐1/PD‐L1/PD‐L2 checkpoints and the antibodies targeting them are as follows:[Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF]

where $Y$, $Y_1$ and $Y_2$ represent PD‐1, PD‐L1, and PD‐L2, respectively, and $A$ and $A_1$ represent the antibodies to PD‐1 and PD‐L1, respectively. $k_{\mathrm{on}}$ and $k_{\mathrm{off}}$ are the binding constants, ${\gamma }_{\mathbf{T}}$ is volume fraction of interstitial space in tumor, $\chi$ is the intrinsic antibody cross‐arm binding efficiency, $A_{\mathrm{syn}}$ is the surface area of the synapse, $d_{\mathrm{syn}}$ is the thickness of the confinement space between the two cells, and $N_A$ is Avogadro's number. In this model, the effect of diffusion of the checkpoint molecules in and out of the synapse has been incorporated by overestimation of the size of the synapse.^28,49 In addition, target‐mediated drug disposition is neglected as it is thought that the antibodies of interest are in general in excess. For specific antibodies that are dosed in low concentrations, it might be necessary to implement the local degradation through target‐mediated drug disposition.

PK module

The PK module aims to capture the population‐averaged plasma concentrations of therapeutics, which are predicted in published model‐based population PK analyses, and predict their concentrations in the tumor.²⁸ This module is called directly by the corresponding PD module and does not need to be called by the user; at the time of this writing, the immune‐checkpoint module is the only PD module implemented. The permeabilities between the central and the other compartments are estimated based on molecular weights of the therapeutics. The other PK parameters, including clearance rate and volume fractions of interstitial space available to the therapeutics, were determined for nivolumab using optimization using pattern search in the Global Optimization Toolbox in MATLAB to fit the plasma concentration^28,50; the PK parameters for other drugs will be added in future releases of QSP‐IO.

The equations governing the PK of antibodies used in this tutorial are described as follows,[Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF][Image Omitted. See PDF]where $Q_{LD}$ refers to the lymph flow rate from the tumor to the blood through lymph nodes^38,51 $k_{cl}$ is the clearance rate from central compartment; and ${\gamma }_{\mathbf{C}}$, ${\gamma }_{\mathbf{P}}$, ${\gamma }_{\mathbf{T}}$, and ${\gamma }_{\mathbf{LN}}$ are the volume fractions of interstitial space available to the antibody in each compartment, respectively. The volumetric flow rates between central and other compartments in volume per time, $Q_{\mathbf{P}}$, $Q_{\mathbf{T}}$, and $Q_{\mathbf{LN}}$, are calculated to be the product of the estimated permeability of the antibody and the total surface area of capillary walls.

EXAMPLE APPLICATIONS

In this section, we outline the procedure of developing a QSP model through two examples; the scripts for generating the examples can be found in the scripts directory included with this package. In both examples, we consider immune checkpoint therapy against PD‐1 using the nivolumab antibody. The parameter values used in these examples were taken from the study by Jafarnejad et al.²⁸ with minor modifications.

SimBiology models are created using MATLAB scripts, which call various functions from QSP‐IO. In general, scripts should perform the five following steps:

1. Load the parameters into the MATLAB workspace from the JSON file
2. Initialize the SimBiology model object
3. Call the individual modules to customize the model
4. Create the SimBiology dose object
5. Generate the initial conditions

In step 1, the parameters are loaded into the workspace by calling the load_parameters function that takes the file name as input. Next, in step 2, the SimBiology object is created by calling the simbio_init function. This function requires the model name, an array of time points at which to evaluate the simulation and the parameter structure that was created in step 1. Step 3 is where the structure of the model is defined by calling various modules; a summary of the modules currently included in QSP‐IO and their dependencies is summarized in Table 1. In step 4, the SimBiology dose object is created; this object dictates which drugs are delivered, the time at which they are administered, and the dose of each drug. The dose object is created using the schedule_dosing function, which takes the names of drugs as input using a cell array of strings; at the time of this writing, nivolumab, durvalumab, and ipilimumab are currently supported; however, support for other common immunotherapy drugs will be added in the future. Optionally, the user may specify the patient weight, dose, and dose schedule through name–value pairs; the default information is hard‐coded into that function as a lookup table and new drugs can be added by users by following the same format. Listing 2 shows example calls to the schedule_dosing function. Finally, the model is completed in step 5 by calling the initial_conditions function, which takes the model object as input and calculates the initial values for all the species in the model. Once the model object is created, it can be manipulated with any of the SimBiology functions. For example, the simulation can be run by calling the sbiosimulate function, which takes the model and dose objects as inputs. The results can be visualized using MATLAB's plot function or QSP‐IO's simbio_plot function.

Listing 2: Examples of calling the schedule dose function

ds = schedule_dosing ({‘nivolumab’}, ‘patientWeight’, 90); % specify weight in kg ds = schedule_dosing ({‘nivolumab’, ‘ipilimumab’}); % specify two drugs % change dosing from default (3 mg/kg every two weeks ) to 1 mg/kg every week for 10 weeks sched = [0 ,7 ,10]; % format [ start_time , interval , number_of_repeats] ds = schedule_dosing ({‘nivolumab’}, ‘nivolumab_dose’, 1, ‘nivolumab_schedule’, sched);

Example 1

In this first example, we simulate the response of a non‐small cell lung cancer tumor to anti‐PD‐1 treatment. For this, we use the model in the study by Jafarnejad et al.²⁸ Their model considers cancer cells, T cells, Tregs, neo‐antigens, self‐antigens, antigen presentation, and the immune checkpoints (PD‐1, PD‐L1 and PD‐L2). The script to create the model in this example is included with QSP‐IO in the example1.m script file and is included in Listing 3 for convenience. In step 1, the parameters for this model are loaded; the parameter values were taken from the study by Jafarnejad et al.²⁸ with minor modifications to account for small differences in the model structure. Steps 2, 4, and 5 follow exactly from the outline in Example Applications. For this example, in step 3, the cancer, T cell, Treg, APC, and antigen and immune checkpoint modules are called. The antigen module is called twice, once to create the neo‐antigens and once to create a competing self‐peptide. All of the modules take the SimBiology model object and a structure containing all of the model parameters as inputs. The cancer module requires an additional input that specifies a unique species name; this is required to allow for multiple subclones of cancer cells to be defined. The T cell module also requires two additional inputs: one that specifies a unique identifier for that specific T cell clone and one that specifies the name of the cancer cell subclones that it can recognize. The antigen module requires two additional inputs: a unique identifier and an antigen object that is created using the create_antigen function; the identifier is required to match that of a T cell clone. To define an antigen as a self‐peptide, the identifier should be set to 0. The create_antigen function requires cancer subclone names that contain the antigen, the intracellular antigen concentration in the cancer cell and optionally, the user may specify an identifier that uses a lookup table to assign the binding rate constants for the antigen to MHC molecules. The immune checkpoint module requires two additional inputs: one that specifies the name of the T cell clone and one that specifies the name of the cancer subclone. Optionally, the user may specify the names of antibodies to PD‐1 or PD‐L1 (see line 30 of Listing 3); if the user does not specify any antibodies, no antibodies are simulated.

Listing 3: Example 1 script

% Model Settings model_name = ‘Immune Oncology Model’; start_time = 0.0; % solution start time [days] % Load Model Parameters ( Step 1) parameter_filename = ‘parameters/example1_parameters . json’; params = load_parameters (parameter_filename); Model Initialization (Step 2) time = start_time : time_step : end_time; model = simbio_init (model_name, time, params); % Add Modules (Step 3) % Cancer Module model = cancer_module (model , ‘C1’, params); % T Cell Module model = Tcell_module (model, ‘1’, params, {‘C1’}); % T Reg Module model = Treg_module (model, params); % APC Module model = APC_module (model, params); % Antigen Module self_antigen = create_antigen ({‘C1’}, 1e -8, ‘antigenID’,0); neo_antigen = create_antigen ({‘C1’}, 1e -8, ‘antigenID’,1); model = antigen_module (model , ‘0’, params, self_antigen); model = antigen_module (model, ‘1’, params, neo_antigen); % PD1 Immune Checkpoint Module model = PD1_module (model, params, ‘T1’, ‘C1’, ‘drugName’, ‘nivolumab’); % Dose Schedule (Step 4) dose_schedule = schedule_dosing ({‘nivolumab’}); % Generate Initial Conditions (Step 5) model = initial_conditions (model);

Figure 3 shows the results of anti‐PD‐1 treatment and no treatment for this example model. Both simulations use the same parameterization, the only difference being that no dose object is specified for the no‐treatment simulation. The parameters for the example were hand tuned to give a parameterization that simulates a virtual patient that responds to treatment. The tumor in the no‐treatment case grows exponentially, whereas in the treatment case, the tumor decreases in size. The number of cancer cells in the treatment case decreases faster than the tumor volume because of the presence of dead cancer cells and T cells in the tumor. The amount of antigen present in the tumor compartment is much larger in the no‐treatment case because the rate of cancer cell death is higher. The density of T cells and Tregs in the tumor reach maximum and decrease because the naïve T cells/Tregs in the blood are limited by their thymic output.

3 Figure. Simulation results as a function of time for control (blue) and anti‐PD‐1 treatment (red). (a) Tumor volume. (b) Concentration of free antigen in the tumor. (c) Naïve T cell density in the blood. (d) Number of cancer cells. (e) Number of antigen‐MHC complex molecules per APC binding site. (f) Activated Treg density in the blood. (g) Number of APCs. (h) Concentration of nivolumab in the blood. (i) Activated T cell density in the tumor. APC, antigen presenting cell; LN, lymph node compartment; mAPC, mature antigen presenting cell; PD‐1, programmed cell death protein 1; Tregs, regulatory T cells.

To show the variability in the dynamics of the model, we simulated 100 parameterizations of the model using a Latin hypercube sampling of 10 parameters (Figure 4). The range for each parameter was chosen to be physiologically realistic; see Table 2 for more details. Figure 4 also shows the simulation results in terms of the percent change in tumor diameter both as a function time for 400 days (Figure 4b) and after 60 days of treatment (Figure 4d) because these are frequently reported in clinical studies and are used for assessing patient response to therapy. The Response Evaluation Criteria In Solid Tumors (RECIST), often used as a metric for evaluating tumor response to therapy, is shown as a function of time (Figure 4c); this demonstrates the sensitivity of this metric on the time of measurement. In this study, RECIST criteria are calculated assuming a single lesion for which the diameter is assumed to be the diameter of sphere with the same volume as the tumor; this approach was taken in previous work.^26–28 Note that because we sampled the parameter space uniformly, the distribution of virtual patient responses shown in Figure 4 does not reflect clinical outcomes. To simulate a virtual clinical trial, the distribution from which the parameters must be sampled needs to be determined. Potential approaches for determining parameter distributions for simulating clinical trials are implemented in refs. 52–54

4 Figure. (a) Tumor volume as a function of time for a Latin hypercube sampling (LHS) (n = 100) of a subset of the parameter space; 5 of the 100 parameterizations did not reach the initial tumor diameter. Ten parameters were varied over a physiological range; see Table 2. (b) Percent change in tumor diameter as a function of time for the LHS. Dashed lines indicate the threshold for the Response Evaluation Criteria In Solid Tumors (RECIST); regions are labeled with PR/CR for partial/complete response, SD for stable disease, and PD for progressive disease. (c) RECIST values as a function of time for the LHS. (d) Waterfall plot showing the percent change in diameter at 60 days following the start of treatment. Each bar represents a simulation, with height representing the percent change in tumor diameter.

2 TableParameter ranges used to generate the Latin hypercube sampling of parameter space for Figure 4

Description	Parameter	Variable name	Range
Cancer cell growth rate	$k_{\mathrm{growth}}$	k_C1_growth	0.001–0.05 day−1
Cancer cell capacity	$C_{\max }$	C_max	1011–1013 cells
Initial tumor diameter	$D_{\mathbf{T}}\left(0\right)$	initial_tumor_diameter	0.5–5.0 cm
Number of T cell clones	$n_{\mathrm{clones}}$	n_T1_clones	1–1000
Rate of thymic output	$\sigma$	Q_nT1_thym	108–1010
Cancer cell death rate from T cells	$k_{T_{\text{cell}}}$	k_C_T1	0.5–50 day−1
T cell cell death rate	$k_{\text{death}}^T$	k_T1_death	0.01–10 day−1
T cell transport rate out of the lymph nodes	$q_{\text{LN}}^{\text{out}}$	q_T1_LN_out	0.1–10 day−1
Cell clearance rate	$k_{\mathrm{clear}}$	k_cell_clear	0.001–0.04 day−1
Number of PD1‐PDL1 complex molecules for half‐maximal inhibition of T cell cytotoxicity	$PD{1}_{50}$	PD1_50	1–50 molecules

Example 2

In this second example, we consider a patient with two different cancer cell clones: one with a larger growth rate than the other. We also consider two different cytotoxic T cell clones: the first clone will only recognize one of the cancer clones, and the second clone recognizes both cancer cell clones. Finally, we simulate 1 year of therapy and follow the tumor progression after the end of treatment. To construct this model, we need to call the cancer module twice, modifying the growth rate for the second call, and the T cell module twice, modifying the rate at which the T cells kill the cancer cells for the second call. We also need to call the immune checkpoint module and create the dose object. For simplicity, this example will not consider antigen presentation or the presence of Tregs.

Figure 5 shows the simulation results for example 2. In the no‐treatment case, the tumor increases rapidly until it asymptotically approaches the cancer cell capacity dictated by $C_{\max }$. In the case of anti‐PD‐1 treatment, the tumor increases in size at the beginning of treatment, followed by a rapid decrease in size. Shortly after the end of the treatment, the C₁ cancer cells are eliminated, leaving only a small (on the order of 10⁵) number of C₂ cells, which go on to grow exponentially.

5 Figure. Simulation of a patient with two different cancer cell clones, C1 and C2, where C2 has a slower growth rate and lower rate of death from T cells. The simulation is run over approximately 8 years (3,000 days), with treatment given for 1 year. (a) Tumor volume as a function of time. (b) Number of cancer cells as a function of time. The vertical dotted line indicates the end of the treatment.

CONCLUSIONS

Following our experience in developing and publishing QSP models for immunotherapy applications, we identified sources that limit the implementation of QSP models for decision‐making in drug development and attempted to address this problem through the introduction of a new open‐source platform we called QSP‐IO. In this tutorial, we introduced the rationale behind the structure of the modeling framework focused on explaining the individual modules and finished by showing two examples to demonstrate the capabilities of QSP‐IO. The framework was designed to allow easy expansion of the model and facilitate adoption and modification of the modules by the industry and academic users. The assumptions behind mathematical representation of each biological process were explained in detail for each module.

Development of this open‐source, modular, and expandable toolbox that is being actively verified and validated through regular use and rigorous peer review could facilitate the resolution of a number of limitations facing wider spread use of QSP models. The reusability of the verified modules reduces the time from conception of the problem to the production of practical results for decision making. The modular platform allows the user to select and refine the minimal model required to address the questions at hand, which reduces the uncertainty from too many degrees of freedom and makes the model fit for purpose. The modular platform also allows the developers and users to expand the features of the model over time to make the toolbox suitable for a variety of targets relevant in the IO area.

Despite all the improvements provided by the QSP‐IO platform, there are a number of areas for further improvement in the toolbox. For example, the current QSP‐IO models can investigate the efficacy of the immune checkpoint‐blocking therapeutics, but toxicity modules are necessary to enable the optimization of the dosing schedule as implemented in refs. 55–57. Furthermore, to facilitate the translation of preclinical animal experimental findings to humans, a parallel parameterization of the model for mice is required to enable the users to exploit the relatively rich data sets from mice for prediction of the response in humans.

In addition, current modules could be improved to include more biological mechanisms. For instance, the immune checkpoint module could be augmented by allowing for checkpoint expression to be dependent on the tumor immune microenvironment, for example, PD‐L1 expression on cancer cells has been shown to be induced by the presence of interferon gamma.^58–61 In addition, other important checkpoint pathways could be included such as CTLA‐4 and tumor necrosis factor recepor superfamily member 4 (TNFRSF4, also known as OX‐40). Models of CTLA‐4 pathways have been implemented in refs. 62 and 63. A pathway model for OX‐40 can be developed using similar methodology based on our current understanding of the biology.^64–66

Thus far, our efforts have focused on immune‐checkpoint therapy, and therefore we have developed detailed mechanistic models of T cells; however, currently, QSP‐IO does not explicitly have modules for innate immune cells such as NK cells, macrophages, or myeloid‐derived suppressor cells. Models for NK cells^55,67–69 and macrophages^69,70 have been implemented in previous studies, and we plan to add modules for macrophages and NK cells based on these studies that will be available in future releases of QSP‐IO.

In conclusion, QSP‐IO serves as the first step toward a more comprehensive framework to run virtual clinical trials. A second follow‐up tutorial is planned for release of an update for the toolbox to enable the user to utilize the model developed by QSP‐IO to generate plausible patient populations and select a virtual patient population that captures the known variability in actual patient populations in the clinical trials.

Funding

This work was supported by National Institutes of Health Grants U01CA212007 and R01CA138264.

Conflict of Interest

The authors declared no competing interests for this work.