Modeling and simulation have developed as important tools for rational decision making in drug development and use. Appropriate models can predict the time course of drug exposure, response, and adverse effects for different dose regimens. Of modeling techniques, the widespread use of population modeling methods has provided a quantitative framework for understanding individual variability in drug exposure and response.

Population modeling is a tool for identifying relationships between observed drug exposure or response and a subject’s physiologic characteristics. It originated from the population pharmacokinetics (PK) modeling introduced in 1972 by Sheiner, et al.1 Initially developed to analyze sparse PK data collected from routine clinics,2 this approach was expanded to modeling drug responses [e.g., pharmacodynamics (PD)] in order to investigate the relationship between PK and PD quantitatively.3 Modeling, especially population modeling, has become an important tool in drug development. The method developed by Sheiner, et al.1 is based on a mixed effect modeling framework that estimates population mean parameters and between-individual variability by pooling sparse data from many individuals. The method also estimates the covariate effects associated with variability in drug exposure and response. With such growing importance of modeling approaches in drug development, in recent years, pharmacokinetic-pharmacodynamic (PKPD) modeling has also become a key tool for modernizing oncologic drug development through early identification and quantification of dose–response relationships.

One of the difficulties with oncologic drug development is that the full dose–response relationship is difficult to be characterized as typically only one or two doses are given to patient population and a placebo group is rarely used. Also, the therapeutic index is typically narrow because drug concentrations causing tumor shrinkage can also cause adverse effects (AEs).

Given such practical difficulties, PKPD modeling could facilitate more efficient oncology drug therapy by systematically assessing various tumor metrics, which can be potential predictors for oncology drug trials.

This review summarizes works published on population PK PD modeling to describe the time course of tumor size, tumor marker and biomarker responses, and adverse effects, as well as model-based predicted tumor metrics which were then subsequently used in survival analyses as predictors of survival along with baseline patient factors.

Fig. 1 depicts a model-based framework for oncology drug development. First, a PK model is developed to describe the relationship between dose and PK metrics, such as trough concentration, area under concentration (AUC), and time course of concentration [C(t)]. Then, a PKPD model is developed to describe the relationship between PK and PD metrics, where PK metrics obtained from PK model are used as inputs in the PKPD model for various PD metrics, such as tumor size metrics {e.g., tumor size ratio [TSR], tumor size [SIZE(t)], time to tumor growth [TTG], tumor growth rate constant [K_grow]}, tumor marker or biomarker metrics, survival time [e.g., progression-free survival (PFS), overall survival (OS)], and AE metrics [e.g., decrease of white blood cell (WBC) count].

Endpoints of clinical efficacy in anticancer drug treatment include objective response rate (ORR), PFS, and OS. Of these, early evaluation of antitumor activity is typically based on ORR, which is defined by Response Evaluation Criteria in Solid Tumors (RECIST) categories:4 complete response, partial response, stable disease, or progressive disease.

The above categorization is based on the following criteria: reduction in tumor size, measured as the sum of longest diameters (SLD) of target lesions and assessment of non-target legions and the appearance of new lesions. RECIST-based categorization of antitumor response may not be appropriate for drugs whose mode of action is to delay disease progression without noticeable tumor shrinkage, where disease stabilization may be more related to survival benefit. Thus, ORR does not always reflect improvement in survival, and many clinical studies, which were effective as assessed by ORR in early-stage of study, have failed to show efficacy in late stages.

PFS is defined as the time from disease occurrence until diagnosis of progressive disease and is used in later stages of studies. This endpoint, which requires a shorter follow-up period than OS, is dependent on the assessment time of tumor response, not always used as a valid surrogate endpoint for OS. According to meta-analyses, while PFS improvement has translated into OS improvement in some cancers (e.g., colorectal),56 no such correlation has been found in others (e.g., breast cancer).78 Both ORR and PFS reduce a set of time-course data of tumor response to a single summary measure of categorized tumor progression status and survival time, thereby ignoring the time-varying longitudinal nature of tumor response. OS represents the time from disease occurrence until death and remains the universally accepted gold standard.

In solid-tumor clinical trials, imaging techniques, such as computed tomography scan and X-ray, are typically used to measure tumor size, which is recorded according to the RECIST criteria4 as SLD measured on a limited number of organs across targeted lesions. Then, according to RECIST, tumor size, originally measured on a continuous scale, is transformed into the 4 categories defined above. Accordingly, this step of categorization facilitates clinical interpretation of measured tumor sizes. However, this method of assessing drug effectiveness based on a point estimate of categorized response has its limitations:910 first, transforming a continuous variable into a categorical variable results in the loss of information. Second, RECIST criteria are evaluated at a set of discrete time points selected a priori, thereby all the dynamic characteristics related to tumor progression, including natural tumor growth, treatment-related tumor shrinkage, and treatment-related resistance development, are ignored.

Models of tumor dynamics are expressed in two different approaches. One approach uses an analytic model described as algebraic equations, and the other uses a set of ordinary differential equations. For an analytic model, Stein, et al.11 proposed the following algebraic equation for tumor size y over time t, where the term "-1" is used due to the initial condition of y(0)=y₀:12

In the above equation, tumor size is assumed to increase exponentially with a net growth rate constant g and to decrease exponentially with the drug-induced decay rate constant d. Other models using algebraic equations are found elsewhere.

Analytical models present some advantages from an implementation point of view. Due to mathematical simplicity, they are easy to implement in classical software programs, and allow for very quick computations. However, they have several disadvantages. First, varying dosing information, such as dose de-escalation and modification, cannot be taken into account. Second, the models are purely empirical in nature and thus cannot be extended to extrapolating the developed model to account for dosing regimen changes within the same study, as well as different dosing regimens in a different design and/or study. Lastly, they have limitations with formulating new hypotheses for mechanisms pertaining to tumor growth and response to drug treatments.

Models of tumor dynamics using ordinary differential equations are generally expressed as12 where "dy/dt" denotes the derivative of y with respect to t, that is, the change of tumor size over time; "growth_{_net}" denotes a function of net growth, that is the difference between natural tumor growth and natural tumor death; and "decay_{_drug}" denotes drug-induced decay processes.

The "growth_{_net}" term can take different forms, including linear growth, exponential growth, and more complex forms, such as logistic and Gompertz growths that separate growth and natural (i.e., non-drug induced) decay as shown below:

The "decay_{_drug}" term in Eq. (2) represents the effect of an anticancer drug and is assumed to follow a first-order decay process as below, which ensures non-negativity of the solution for y.

The "effect" term represents a constant or a function of a PK metric (e.g., plasma concentration), e.g., where C(t) denotes the drug plasma concentration. The "effect" term can also be nonlinearly related to C(t), say, through an Emax model. The "β" term can be a constant or a function of time to reflect the decay with time as below

This expression in Eq. (6) represents the loss of drug effect over time due to the emergence of "resistance."

From 2008 up to 2014, 13 papers have been published in eight different therapeutic areas to propose models for the time course of tumor size in patients: colorectal cancer,1314 nonsmall cell lung cancer (SCLC),1516 renal cell carcinoma,17181920 thyroid cancer,21 metastatic breast cancer,22 prostate cancer,11 gastrointestinal stromal tumor1723 and low-grade glioma.24 See Ribba, et al.12 for a summary of these models.

Claret, et al.25 developed the tumor growth inhibition (TGI) model using data collected from colorectal cancer patients who received capecitabine with a schedule of 2 weeks on followed by 1 week off or 5 fluorouracil daily for 5 consecutive days every 4 weeks. This model has been applied to several cancer types and drugs by other investigators. The TGI model is described in Eq. (7) and (8).

In Eq. (7) and (8), y(t) is the tumor size (=SLD) at time t, y₀ is the baseline tumor size, K_L is the tumor growth rate, K_D(t) is the drug-constant cell kill rate that decreases exponentially with time (according to λ) from an initial value of K_D,0 to account for the progressive development of resistance, and Exposure(t) is the drug exposure at time t. Due to no concentration data available in their work, daily dose was used as the driving force of drug effect.

The key to this model is that it accounts for three important clinical features of tumor progression in anticancer drug treatment (the dynamics of tumor growth, antitumor drug effect, and resistance to drug effect) in one model, based on a previously published simulation model.26 Describing tumor size as a function of time and drug exposure, it accounts for the natural tumor growth and the drug action on the tumor (i.e., tumor cell kill driven by drug exposure). With a first-order tumor growth rate, the model incorporates a resistance process to describe the regrowth of tumor.

However, it should be noted that, because dose is not included in the model, it cannot be used to predict tumor response under different dosing regimens.

Tumor markers are produced by cancer or other cells in the body in response to cancer. Examples are prostate specific antigen (PSA) in prostate cancer, M-protein in myeloma, cancer antigen 125 (CA125) in ovarian cancer, and carcinoembryonic antigen (CEA) in colorectal cancer.27 Since these tumor markers can be readily measured in blood, they may represent the total body burden of cancer better than tumor SLD. According to RECIST 1.1, SLD measurements of tumor only assess a maximum of five target lesions,4 and conventional two dimensional scans, on which RECIST is based, do not capture changes in tumor density.28 Also, tumor SLD measurements are not objective, costly, and the assessment is performed every 6 to 8 weeks.

Men with prostate cancer or with other prostate disorders often show elevated PSA levels. You, et al.29 used a bi-exponential model to describe in prostate cancer patients the decline in PSA levels with time after prostatectomy. They also described using mono-exponential models alpha fetoprotein and human chorionic gonadotropin levels in germ cell tumor patients receiving conventional chemotherapy.30

In multiple myeloma, monoclonal immunoglobulin proteins (called M-protein) are produced in excessive amounts by malignant plasma cells. Using a drug exposure-driven model with the same model structure as in the TGI model (equation 7), Jonsson, et al.31 described M-protein levels over time in multiple myeloma patients receiving dexamethasone treatment.

Recently, Desmée, et al.32 using a simulation approach showed that joint modeling of PSA kinetics and survival time produces a precise estimation of PSA time-course and survival parameters, compared with two simplified alternatives, two-stage and joint sequential models. They suggested the developed method as a way to improve treatment prediction and evaluation in oncology.

Fig. 2 depicts a schematic diagram of the model used in the work by Desmée, et al.32 In the absence of treatment, it is assumed that prostatic cancer cells, C, proliferate with rate Kprol and eliminate with rate Kd. PSA is produced and secreted with rate Kp and eliminated from the blood with rate Ke. Then, it is supposed that a chemotherapy with time-varying effectiveness, e(t), acts by blocking cell proliferation, and thus, the proliferation rate under treatment is suppressed, becoming Kprol·[1-e(t)] with 0≤e(t)≤1. Mathematically,

Kanefendt, et al.33 fitted indirect response (IDR) models to the time courses of soluble vascular endothelial growth factor receptor 2 (sVEGFR-2) and soluble vascular endothelial growth factor receptor 3 (sVEGFR-3) data in metastatic colorectal cancer patients treated with sunitinib.

Hansson, et al.2334 used the IDR model to describe vascular endothelial growth factor (VEGF), sVEGFR-2, sVEGFR-3, and soluble stem cell factor receptor (sKIT) data in gastrointestinal stromal tumour (GIST) patients receiving sunitinib. Based on correlations between the biomarker responses analyzed, they proposed an analytical framework for investigating the relationships among model-predicted time courses of five metrics: drug exposure, biomarkers, tumor SLD, adverse effects, and survival.

These approaches used in sunitinib analyses support using population PKPD modeling technologies for biomarkers to 1) examine which biomarkers are worthwhile to be investigated, 2) better understand the action mechanism of the drug, 3) early evaluate treatment efficacy, and 4) predict long-term treatment outcomes.

Buil-Bruna, et al.35 modeled lactate dehydrogenase (LDH) and neuron specific enolase (NSE) concentrations in SCLC patients. In their work, a modeling framework was proposed that relates LDH and NSE, which are circulating biomarkers in the plasma and are easily obtained from patients, to tumor progression levels assessed by RECIST categories. LDH and NSE are known to be independent prognostic factors for SCLC. In their model, an underlying latent variable representing unobserved "disease level," corresponding to tumor size dynamics, was incorporated as a driving source of biomarker production influenced by exposure to treatment. They showed that model-based unobserved disease levels are strongly correlated with RECIST criteria-based disease progression measurements, suggesting the feasibility of circulating biomarkers in predicting treatment outcomes as powerful tools to monitor disease.

A schematic diagram of the final model is described in Fig. 3 and differential equations related to the model are as follows:

In the model diagram and equations, notations are as follows: C(t) is plasma drug concentration associated with chemotherapy, D(t) is a latent variable that represents disease progression and drives LDH and NSE production, and R(t) is drug resistance, modeled by linking cumulative drug exposure with a decrease in the drug effect, i.e., R(t)=exp(-γ·∫C(t)dt) and K_IN,LDH and K_IN,NSE are zero-order rate constants and K_D,LDH, K_OUT,LDH, K_D,NSE and K_OUT,NSE are first-order rate constants. Granulocyte colony-stimulating factor (GCSF) increases the physiological LDH synthesis.

Neutrophils and platelets, which are essential for fighting infections and blood clotting, respectively, have the characteristics of rapid proliferation and therefore are easily affected by chemotherapy. As a result, myelosuppression has been one of the most frequent drug adverse reactions encountered during chemotherapy. Using leukocyte and neutrophil data obtained from several chemotherapy drugs, Friberg, et al.36 developed a myelosuppression model, which has been the most widely used myelosuppression model so far.

Other types of adverse drug effects induced by chemotherapy that have been investigated using model-based approaches include the following:27 1) the work by Agoram, et al.37 in chemotherapy-induced anemia on the hemoglobin time course whose production is stimulated by darbepoetin alfa; 2) the work by Fetterly, et al.38 in trabectedin-induced liver toxicity on the time course of alanine aminotransferase (ALT); 3) the works by Keizer, et al.39 in E7080 therapy and Hansson, et al.40 and Houk, et al.17 in sunitinib therapy on the time course of elevated diastolic blood pressure; 4) the work by Xie, et al.41 in irinotecan therapy on the occurrence of diarrhea; and 5) the works by Hénin, et al.42 in capecitabine therapy and Hansson, et al.40 in sunitinib therapy on the occurrence of hand-and-foot syndrome.

In the work by Friberg, et al.,36 chemotherapy-induced myelosuppression was described using a semi-mechanistic PKPD model, developed based on leukocyte and neutrophil data obtained in patients after administration of docetaxel, paclitaxel, and etoposide, which was then applied to myelosuppression data obtained from 2′-deoxy-2′-methylidenecytidine, irinotecan, and vinflunine administrations. The schematic diagram of the model is given in Fig. 4. In the model, Prol denotes a proliferating compartment that represents stem cells and progenitor cells sensitive to drugs; Kprol denotes a proliferation rate constant of cells in Prol, that is, rate constant of self-renewal or mitosis for generation of new cells in Prol, which was dependent on the number of cells in the compartment; Ktr denotes a rate constant between transit compartments, which was incorporated to allow a time delay between drug administration and the observed effect; n is the number of transit compartments, which was chosen to be 3; Circ denotes a compartment of observed circulating blood cells with Circ₀ being the baseline value; Kcirc denotes an elimination rate constant for blood cells from the circulating compartment; and E_Drug denotes drug effect, which was assumed to be either a linear (E_Drug=Slope·Conc) or an Emax model [E_Drug=Emax·Conc/(EC50+Conc)], where the drug concentration (Conc) was assumed to reduce the proliferation rate or induce cell loss by the function E_Drug. The assumption Kprol=Ktr was used to guarantee dProl/dt=0 at steady state, and Kcirc=Ktr was used to minimize the number of parameters to be estimated.

The differential equations were written as

As illustrated in Fig. 1, once models for tumor metrics, such as tumor size (or tumor size ratio or tumor growth rate, etc.), tumor markers [PSA, alpha fetoprotein (AFP), human chorionic gonadotropin (hCG), CA125, etc.], and biomarkers (LDH, VEGF-A, sVEGFR, etc.), are established, these tumor metrics can be implemented into parametric time-to-event (TTE) models to provide a model-based prediction for PFS and OS.

Central to implementing TTE models is the proper choice of a hazard model, a general form of which is described as below:

h(t) in equation (20) is the hazard function for the Cox proportional hazard model, representing the instantaneous rate at which an event (e.g., death) occurs. It consists of the baseline hazard function h₀(t) and the explanatory variables x_i, i=1,2,..., n. h₀(t) is defined by a set of estimated parameters and describes how the risk of event changes over time at baseline levels of covariates. x₁, x₂,..., x_n represent predictors (e.g., tumor metrics mentioned above or patient covariates that were related to PFS or OS). These predictors can be a baseline metric for each patient (e.g., baseline tumor size), an individual parameter estimate (e.g., tumor growth rate), or a time varying metric [e.g., tumor(t), tumor marker(t), or biomarker(t)]. Also, patient baseline characteristics (e.g., tumor stage, Eastern Cooperative Oncology Group performance (ECOG) status, number of lesions, etc.) are used as predictors. The size of the coefficients β₁, β₂, ..., β_n denotes the relevant importance of the corresponding predictor in the model.43

In the work by Claret, et al.,25 an OS model was described by baseline tumor size, denoting patient characteristics, and change in tumor size at week 7 relative to baseline, denoting drug effect. Poorer prognosis was associated with larger baseline tumor size and smaller tumor shrinkage at week 7.

In the work by Buil-Bruna, et al.,35 based on the Cox proportional hazards model in equation (20), PFS was found to be significantly influenced by the predicted change in individual disease level (=latent tumor size), which was defined as where $D_{t_{ji}}$ is the predicted disease level for j-th patient at i-th observation time point (i.e., the current CT scan), and $D_{t_{ji+1}}$ is the predicted disease level for j-th patient at (i+1)-th or subsequent observation time point (i.e., the following CT scan, approximately 8 weeks later).

The conventional approach of evaluating oncology drug response based on RECIST criteria is limited in that it results in the loss of information by transforming a continuous variable into a categorical variable and allowing tumor size evaluation only at a set of discrete time points selected a priori. Through a review of published works, this report intended to suggest a model-based approach as a promising tool to be used in oncology drug therapy to overcome limitations with the traditional RECIST-based approach.