|Home | About | Journals | Submit | Contact Us | Français|
Pharmacokinetic (PK) and pharmacodynamic (PD) models seek to describe the temporal pattern of drug exposures and their associated pharmacological effects produced at micro- and macro-scales of organization. Antibody-based drugs have been developed for a large variety of diseases, with effects exhibited through a comprehensive range of mechanisms of action. Mechanism-based PK/PD and systems pharmacology models can play a major role in elucidating and integrating complex antibody pharmacological properties, such as nonlinear disposition and dynamical intracellular signaling pathways triggered by ligation to their cognate targets. Such complexities can be addressed through the use of robust computational modeling techniques that have proven powerful tools for pragmatic characterization of experimental data and for theoretical exploration of antibody efficacy and adverse effects. The primary objectives of such multi-scale mathematical models are to generate and test competing hypotheses and to predict clinical outcomes. In this review, relevant systems pharmacology and enhanced PD (ePD) models that are used as predictive tools for antibody-based drug action are reported. Their common conceptual features are highlighted, along with approaches used for modeling preclinical and clinically available data. Key examples illustrate how systems pharmacology and ePD models codify the interplay among complex biology, drug concentrations, and pharmacological effects. New hybrid modeling concepts that bridge cutting-edge systems pharmacology models with established PK/ePD models will be needed to anticipate antibody effects on disease in subpopulations and individual patients.
Monoclonal antibodies (mAbs) are utilized as primary therapeutics for a broad array of diseases, especially in the treatment of chronic conditions such as cancer, inflammation, and autoimmunity. Therapeutic applications of mAbs are diverse and include: 1) immunotoxicotherapy, in which mAbs alter the pharmacokinetic (PK) properties and activity of soluble ligands (e.g., drugs, cytokines, toxins, and others) to prevent their associated adverse effects, 2) cancer therapy, in which mAbs induce killing of target cells (e.g., antibody-dependent cell-mediated cytotoxicity), 3) impairment of a specific cellular function (e.g., platelet aggregation), and 4) use as vectors for targeted delivery of conjugated cytotoxic agents.1-3 This biological diversity in the pharmacological action of mAbs requires an extension of general PK and pharmacodynamic (PD) principles and the construction of more detailed mathematical models based on the major patho-physiological and drug-specific processes regulating drug action.4
A plethora of PK/PD models reported in the literature aim to describe the time-course of mAb exposures and effects, and quantify drug- and system-specific parameters and biological rate-limiting steps that control the activities of mAbs. A general scheme depicting the 4 major areas in which systems models for mAbs have the ability to substantially influence drug development is shown in Fig. 1. The first area depicts disease processes and PK behavior, such as target-mediated disposition, in which models are applied to describe how mAbs bind to their pharmacological targets and trigger or block various bio-signals mediated through networks of intracellular proteins. The second area involves cellular signal transduction, in which systems models of cellular responses can serve as a bridge linking drug exposure to macro-scale outcomes. The third area relates to translational physiology, in which models translate the pharmacological activities of mAbs in preclinical systems to simulate clinical outcomes. The fourth area is that of disease progression, with models describing the population variability in the temporal evolution of biomarkers, clinical symptoms, or disease scores. Altogether, integration across these 4 scales necessitates novel prospective strategies to develop accurate and unbiased mathematical models that can anticipate drug action in patient subpopulations and individuals.
Quantitative systems pharmacology (QSP) is a bourgeoning modeling field that seeks to better understand pathophysiology and pharmacological mechanisms from the levels of genetic molecular pathways (i.e., –omics), regulatory networks, cells, tissues, organs, and ultimately the whole organism.5 QSP is a multi-disciplinary endeavor, borrowing from: 1) genomics, which studies the DNA-sequence of genes and their interactions within the genome; 2) transciptomics, which investigates gene expression profiles; 3) proteomics, which analyzes the structure, function, and interaction of proteins; 4) metabolomics, which identifies metabolic changes in biochemical pathways, complex biomarker combinations, toxic effects, and disease progression; and 5) signaling network modeling, which describes cellular communication and functions through intracellular protein dynamics.6,7 As a result, computational methods for system network analyses have been developed and provide a powerful means for analyzing experimental findings, interpreting biological behaviors, testing hypotheses, and designing experiments by integrating large amounts of experimental data into mathematical models of regulatory networks and cell behavior.8,9
Systems network analyses are typically structural or dynamical in nature.10 Structural analysis provides insights into the general properties of the network through its dissection into functional units and the identification of key pathways that influence network behavior.11 Dynamical networks combine mathematical models with experimental data to characterize propagation of signals through the network.12 The type and extent of network models depend upon several factors, such as the desired outcome to be captured, the expected level of mechanistic detail of the system, and details of the biological system itself. Models extending from molecular to organismal systems and from minimal to high-level mechanistic details have been applied, including: 1) regulatory networks featuring inter-communication between system components; 2) Bayesian networks representing probabilistic relationships between system variables (e.g., genes and proteins) and their conditional dependencies (e.g., disease states); 3) Boolean networks that describe system behavior under different combinations of perturbations; 4) nonlinear ordinary differential equations (ODE), which also describe the dynamical features of the system; and 5) stochastic models that capture random behavior of potential outcomes caused by random variation in the inputs over time.13-16
The main objective of the current work is to review major systems PK/PD models in use for therapeutic mAbs. We highlight important model characteristics, the approaches used to integrate multiple system complexities into these models, the challenges of defining redundant and robust pathways in sufficient detail, and the reconciliation of biological scales into a multi-scale model architecture that vertically combines molecular, cellular, and macroscopic scales.
The absorption and disposition of antibodies are complex, challenging the characterization of their PK and the overall development of new antibody-based drug candidates. Common variables regulating antibody disposition include: 1) the location site of the pharmacological target (e.g., soluble or membrane bound), 2) immunogenicity, by generation of endogenous antibodies against the therapeutic mAb, 3) concomitant medications, and 4) patient covariates and disease status.17 Whereas many small molecule drugs are cleared through kidney filtration, this elimination pathway remains relatively insignificant for therapeutic mAbs due to their large size (molecular weight ~150 kDa), restricting their passage through the glomerulus. The three main processes by which mAbs are distributed and removed from the body are the salvage elimination pathway involving the neonatal Fc receptor (FcRn), tissue catabolism, and target-mediated drug disposition (TMDD).1,17
The FcRn is a ubiquitous membrane protein that is also found on the surface of certain specialized immune cells. The FcRn-mediated recycling pathway was first suggested by Brambell et al.18 It was hypothesized that the binding of immunoglobulin molecules (IgG) to a saturable target protects them from degradation in the circulation and that this physiological process results in a relatively extended half-life for IgG (average 25 days).19 Brambell's hypothesis was eventually supported by several independent studies, and the FcRn recycling pathway was identified as Brambell's saturable salvage pathway.20-22 FcRn binds exclusively to the constant domain Fc (fragment crystallizable) region of therapeutic mAbs and the mAb-FcRn complex may trigger intracellular catabolism. Generally, it is well accepted that elimination from the FcRn pathway is secondary for monomeric mAbs, considering the homologies between the Fc domain of endogenous IgG and mAbs for the binding to FcRn, their high affinities to FcRn, and the relatively high endogenous IgG plasma concentrations. In contrast, soluble immune complexes of mAbs with FcRn and mAbs attached to their cell-surface targets results in primary drug clearance via antibody opsonization and phagocytosis.
Once mAb-FcRn complexes are endocytosed, a significant fraction of mAbs is sorted toward the cell membrane surface via the FcRn-mediated recycling process, then released into the extracellular space. Antibodies bind to FcRn with a pH-dependent affinity.18,20 Thus, in the early endosomes where pH is low (i.e., acidic environment), antibodies are tightly attached to the receptor. As the pH increases to physiological values, the affinity of the antibodies toward FcRn gradually diminishes, resulting in their dissociation from the receptor and release into plasma or interstitial fluids. In knockout mice with impaired expression of the FcRn, the clearance of IgG was increased 10-fold, corresponding to a recycling efficiency of 90%,20 confirming the role of FcRn-mediated recycling in systemic exposures. However, since the expression of FcRn on cell surfaces is limited, the FcRn mediated-recycling process is capacity limited. Elevated antibody concentrations following very high doses (e.g., IVIG therapy) in plasma or extravascular fluids are capable of saturating the FcRn-mediated recycling process, resulting in a low recycling efficiency and an increase in the rate of their intracellular catabolism and clearance. Increased IgG elimination can lead to a decrease in endogenous immunoglobulins, which could be of benefit in mAb treatment of autoimmune diseases, in which effects may be achieved by decreasing plasma concentrations of endogenous auto-antibodies. Examples of systems PK models integrating FcRn-mediated recycling are provided in section (a.2) below.
The binding of mAbs to their pharmacological target, either soluble or on the cell surface, can trigger receptor-mediated endocytosis (RME) and intracellular catabolism. For TMDD systems,23,24 the pharmacological target influences the kinetics of mAb distribution and elimination. Due to the finite number of targets, the TMDD pathway is of limited capacity, which in part explains the nonlinear PK behavior of the majority of therapeutic mAbs. Several factors can influence the rate of uptake and elimination of mAbs by the TMDD pathway, including the level of the injected dose, the total capacity of the target, and the kinetics of receptor turnover (i.e., internalization and degradation).25 Mager and Jusko25 derived and emphasized the general characteristics of TMDD models (Fig 2). Since the original model was published, a myriad of examples have been reported, showcasing state-of-the-art mathematical models of mAbs exhibiting TMDD.26 Meijer et al.27 showed that T cell depletion by modulation of their CD3 surface antigens with an anti-CD3 antibody results in a decreased clearance of the mAb at subsequent doses. Ng et al.28 described the PK of TRX1, an anti-CD4 mAb, simultaneously modeling the time-course profiles of free drug exposure and free and total CD4 in healthy subjects.
One of the computational challenges encountered while using TMDD models is the reliable estimation of the drug-binding micro-constants parameters (i.e., kon and koff) from routine exposure time-course data. To overcome this challenge, several approximations to the general TMDD model have been developed. A quasi-equilibrium (QE) solution to the general TMDD model was derived and validated.29 It assumes rapid binding of the drug to its target, and the drug-binding micro-constants (kon and koff) are substituted with the equilibrium dissociation constant (KD = koff/kon). This model was successfully applied by Hayashi et al.30 to describe the nonlinear disposition of omalizumab, an anti-IgE antibody in the treatment of asthma. Other approximations of the original TMDD model were derived by Gibiansky et al.,31 in particular the quasi-steady-state (QSS) model, in which measurements of the total drug receptor and drug concentrations may allow for estimation of the unbound receptor concentrations. The model successfully described the PK time-course profiles of several mAbs, such as anti-IL632 and anti-neuropilin-1.33 Mathematical equations for the QE and QSS models are identical, except for their corresponding constants KD and Kss, with Kss = KD + Kint/Kon. Additionally, it has been demonstrated that when the amount of the drug significantly exceeds the amount of the target, models with Michaelis-Menten elimination may be sufficient to describe TMDD processes.34,35
TMDD models for systems with multiple targets have been derived,36 and the 2-target TMDD model can be approximated by equations resembling the inclusion of 2 Michaelis-Menten elimination terms in parallel (with different Vmax and KM parameters). Gibiansky et al.37 extended the application of TMDD models, including the QE, Qss, and the MM approximations, to describe the PK of antibody-drug conjugates (ADCs) with model-based simulations for ado-trastuzumab emtansine as an example ADC.
Physiologically-based pharmacokinetic (PBPK) models describe mAbs disposition using a more mechanistic physiological compartmental system. The major structural elements of PBPK models are the actual organs or tissues connected in an anatomical manner by blood (and sometimes lymphatic) circulation. The structure of a PBPK model is predetermined and is derived from anatomical and physiological structure of the animal species of interest. The parameter set of PBPK models include 2 distinct subsets: a drug-independent subset, which is derived from the underlying physiological processes, and a drug-specific subset that describes the PK properties of the mAbs.38 PBPK models for mAbs have the potential for interspecies scaling,39,40 predicting the PK of the mAbs prior to in vivo studies,41,42 and characterizing the physiological and biochemical basis for altered PK and associated toxicities.43,44
Growing interest in the disposition of antibodies is leading to the development of advanced mathematical models to account for drug and system complexities and to allow for improved predictions of mAbs PK.45 PBPK models are uniquely suited for this purpose as they facilitate the incorporation of processes such as convective transport, FcRn binding and recycling, TMDD, and tissue catabolism, resulting in a more accurate and informative description of mAbs disposition. A PBPK model could aid in understanding the factors that determine mAbs distribution in the body, and in identifying functions that are responsible for its efficient site-specific localization, leading to an overall improvement in the desired therapeutic outcome. PBPK models provide a systems approach to describing mAb PK, which may be used to guide the development of antibody-based compounds and explore biophysical properties that might translate into improved plasma and tissue specific PK. The physiological basis of PBPK models (e.g., organ volumes and blood flow rates) provides a natural framework for understanding inter-species differences in mAb PK and for the scale-up of preclinical biodistribution data to predict antibody disposition and to recommend a first dose in humans.38
PBPK models have been successfully developed to describe the disposition of mAbs in plasma and in various tissues.46 Detailed PBPK models of macromolecules typically feature the sub-compartmentalization of tissue compartments into vascular, intracellular, endosomal, and interstitial spaces to capture the permeability-limited distribution and non-uniform distribution within tissues. FcRn plays a major role in mAb disposition,47 and Garg et al.48 incorporated the influence of FcRn on IgG disposition through an FcRn-mediated transcytosis across vascular endosomal cells in a PBPK model for an IgG1 anti-platelet mAb. The model predicted well the IgG plasma and tissue kinetics in control and FcRn knock-out mice. The endosomal space has been further compartmentalized into a catenary series49 that integrates several features, including the transit time of IgG molecules through endosomes, the gradual acidification of endosomal contents, and the association and dissociation of the IgG-FcRn complex under changing pH conditions. This model also incorporates tissue-specific FcRn expression. The simplifying assumptions in this model include: pH conditions are fixed within the vascular endothelial sub-compartments, a fixed transit time for all sub-compartments, and the FcRn binding parameters are identical to those from in vitro experiments. Model-based simulations were compared to those of a non-catenary models,48,50 in which equilibrium binding of IgG and FcRn are assumed. Although the model predicted modest changes in the mAb half-life, it also suggested that increased rates of IgG–FcRn association provides a greater benefit to extending plasma half-life than decreases in the rate of IgG–FcRn dissociation. This catenary structure also accurately describes the PK of pH-sensitive mAbs or so-called ‘catch and release’ mAbs.51
The development of immunogenicity, or the production of anti-drug antibodies (ADA) following the administration of a therapeutic mAb, can alter the PK or PD characteristics of mAbs. It may originate from biological molecular structure recognition, suboptimal dosing, or temporary treatment discontinuation. Regulatory documents and literature reports have provided guidance to meet the regulatory requirements of characterizing and quantifying immune responses to mAb therapeutics.52,53
Circulating ADA can exhibit different specificities and binding affinities and have different effects on the PK/PD of mAbs: 1) a neutralizing effect on drug activity by interfering with mAb binding to its pharmacological target, 2) a non-neutralizing effect on drug activity and a sustaining effect on the PK, and 3) non-neutralizing effect on drug activity, but an enhanced elimination of the mAb. The mechanisms by which ADA alter the PK of mAbs are not completely known. Both binding and neutralizing ADAs have the potential to alter mAb clearance.54 It has been suggested that co-medication with immune-suppressants might suppress the production of ADAs and decrease the risk of immunogenicity.55 Adjusting the dosing of mAbs to target excess circulatory concentrations might also prevent immunogenicity. Although the use of ADA titers is challenging, quantifying the extent of the ADA response might improve the characterization of altered PK/PD properties in the presence of immunogenic reactions.56-58
For most cases, mAb plasma exposure and subsequent pharmacological effects are temporarily disconnected. Underlying physiological processes, such as interactions of the drug with the turnover of endogenous targets, make the indirect response models (IDR)59 essential in characterizing the PD of many mAbs. In general, IDR models describe the turnover process of pharmacological biomarkers by incorporating a zero-order production and first-order removal rate constants. Depending on their mechanism of action, the pharmacological activity of mAbs alter either of these 2 processes through a stimulatory or inhibitory capacity-limited (i.e., Hill-type) function. Several variants to the basic IDR models have been utilized to characterize the PD of mAbs in multiple disease conditions. In this section we describe some selected literature examples of mechanism-based PD models for several FDA-approved therapeutic mAbs. All mAbs referenced in this review article are summarized in Table 1.
Preclinical efficacy of TDM1: Jumbe et al.60 used a population-based PK/PD modeling approach to describe the anti-tumor activity of TDM1, a humanized monoclonal ADC specific to the human epidermal growth factor receptor 2 (HER2). A two-compartment model captured well TDM1 PK in mice bearing HER2 xenografts. A cell-cycle-phase non-specific tumor cell kill model61 including transit compartments62 described well the growth of xenograft tumors and the anti-tumor activity of TDM1. The mathematical model was useful at identifying key determinants of efficacy for TDM1, such as the growth rate of tumor cells, the anti-tumor activity of TDM1 was independent of dosing regimen, and the tumor response was linked to the ratio of exposure to a concentration required for tumor stasis.
Toxicity of TDM1: Thrombocytopenia (TCP) is a dose-limiting toxicity of TDM1. A semi-mechanistic population PK/PD model was developed by Bender et al.63 to characterize the effect of TDM1 on patient platelet counts. The PD model is consistent with a general model of myelosuppression developed by Friberg et al.64 The average plasma concentration of TDM1 was used as a driver for the blockade of megakaryocyte production in the bone marrow. The model captured well the platelet profiles and predicted the incidence of grade ≥3 TCP. The major findings from this model were the identification of the dose of 3.6 mg/kg given every 3 weeks as being less thrombocytopenic, and, in a subgroup of patients exhibiting variable degrees of downward drifting platelet profiles, the TCP was shown to stabilize above grade 3 after the eighth treatment cycle with TDM1. These detailed mathematical modeling examples of TDM1 preclinical and clinical drug effects identified the major determinants of efficacy and safety.
Shah et al.65 developed a multi-scale PK/PD model for brentuximab vedotin. This model serves as another excellent example of preclinical to clinical translation of ADC efficacy. Experimental data were extracted from various sources and used to build and qualify the final PK/PD model. The ADC cytotoxic effects were integrated into the model through a capacity limited function driven by the tumor concentration of the cytotoxic payload. The major finding of this model was the identification of tumor concentrations of monomethyl auristatin E, brentuximab vedotin's payload, as the driver for the preclinical efficacy of the ADC. In addition, the model was successfully translated to humans, and the final PK/PD model predicted observed clinical outcomes of brentuximab vedotin, such as progression free survival rates and complete response rates.
Marathe et al.66 modeled the complex dynamic system of bone remodeling in multiple myeloma (MM) patients under denosumab treatment. Denosumab is an IgG2 mAb directed against the receptor activator of nuclear factor-κB (RANK) ligand (RANKL), and the overexpression of RANKL is partly responsible for the increased bone degradation and decreased bone formation in MM. The PK of denosumab was well characterized with a simple TMDD model. The pharmacological effects of denosumab on a serum N-telopeptide (NTX) biomarker and the cellular interactions between osteoblast-osteoclasts were driven by free denosumab plasma concentrations. The complex interplay between active osteoblasts and active osteoclasts, as well as several critical regulatory proteins in the RANK-RANKL-osteoprotegrin (OPG) signaling pathway, was achieved through the adaptation of a prior mechanistic model by Lemaire et al.67 In this proof of concept study, denosumab PK was integrated, RANK occupancy was recalculated, and the NTX biomarker was successfully linked to a state variable in the model. The major findings of this mathematical model is that it not only assembles the fundamental mechanisms involved in the regulation of bone turnover, but also allows the assessment of the role of key cellular variables and the effect of denosumab exposure on the RANK-RANKL-OPG pathway. For example, increased RANKL expression with disease progression might explain diminished long-term drug responsiveness, and such hypotheses could be further tested and used to guide the design of effective therapeutic strategies.
The PD of rituximab (RTX) suggest that tumor burden can influence its concentration-effect relationships.68 The anti-cancer effect of RTX was described with a stimulatory capacity-limited function on a zero-order rate constant. The most important result from this model is that it was able to predict that high RTX plasma concentrations are associated with better anti-tumor responses and longer mouse survival. In addition, it showed that high tumor burden on Day 13 was associated with low RTX concentrations, an impairment of RTX concentration-effect relationship with low RTX efficacy, and both decreased tumor response and shorter mouse survival. Ternant et al.69 developed a PK/PD model for RTX and predicted the progression-free survival (PFS) of patients treated with RTX alone and in combination with chemotherapy (R-CHOP) in relapsed/resistant patients with FL. The impact of the polymorphism in the FCGR3A-F gene encoding the FcγRIIIa receptor was integrated in clinical efficacy of RTX. This model predicted a potential benefit of maintenance dose of RTX (1500 mg/m2) for both RTX monotherapy and R-CHOP, and it also showed that the PFS of FCGR3A-F carriers remains lower than that of homozygous FCGR3A-VV patients, even with markedly increased RTX dose levels.
Alemtuzumab targets CD52, a protein present on the surface of mature lymphocytes B, and drug-target binding initiates antibody-dependent cell-mediated cytotoxicity. Alemtuzumab is used as second-line therapy for B-CLL and was approved by the FDA for CLL patients who have been treated with alkylating agents and have failed fludarabine therapy. Mold et al.70 developed a population PK/PD model for alemtuzumab in this patient population. The PD of alemtuzumab was measured by whole blood cells, and the time-course of response was modeled with a stimulatory loss indirect response model. The final model integrated a direct relationship between maximal trough concentrations and clinical outcomes, and identified that the probability of achieving a complete or partial response was ≥ 50% when the maximal trough concentration exceeded a defined threshold.
Efalizumab is a humanized mAb directed against CD11a on the cell surface of T cells. CD11a binding to efalizumab inhibits CD11a-mediated cell signaling in T cells, and thus acts as an immunosuppressant by blocking lymphocyte activation and cell migration into tissues. Efalizumab was originally indicated for the treatment of chronic moderate to severe plaque psoriasis; however, due to drug-associated fatal brain infections in treated patients, efalizumab was withdrawn from the market in 2009. A PD model for efalizumab in subjects with psoriasis was first developed by Bauer et al.71 An indirect response relationship to describe CD11a turnover with stimulation of the removal process via a Hill function driven by the free plasma concentrations of the antibody successfully captured the time-course of CD11a cells. Later, Ng et al.28 developed a mechanism-based PK/PD model that characterized the dynamic interaction of efalizumab binding to CD11a followed by efalizumab effects on CD11a expression on T cells and on the psoriasis area and severity index, an efficacy end point for psoriasis. The model described the data in psoriasis patients reasonably well, and could be used to suggest an equivalent achievable efficacy with less frequent dosing of the mAb.
Mold et al. derived a semi-mechanistic PK/PD model of a cynomolgus macaque monkey-human chimeric anti-CD4 mAb, clenoliximab, in patients with moderate to severe rheumatoid arthritis (RA).72 The model was developed sequentially: first, the PK of mAb was defined with a 2-compartment model with Michaelis-Menten elimination from the central compartment, and the binding of clenoliximab to CD4 on T lymphocyte membranes was incorporated into the model. Second, the subsequent profiles of the drug-receptor complex were generated by simulation and used as a driver for the inhibitory effects of clenoliximab on the production of CD4 receptor density. The model was useful in that it provided a framework to examine relationships between clenoliximab dosing regimens and the initiation of durable immune tolerance and to assess the impact of tolerance on the treatment of RA. In a comparative study conducted by Sharma et al.,73 the PK and PD of clenoliximab and keliximab, another anti-CD4 mAb, were described with one unique structural model. A two-compartment model with a saturable clearance pathway captured the PK of both antibodies, and an IDR model with stimulation of CD4+ cell removal described well the PD. A limitation of this model is that, although plasma concentration-time data for clenoliximab and keliximab were captured with a unique set of parameter estimates, their pharmacological immunosuppressant activities differed, with system parameters were specific to each mAb. The major finding from this model is that it showed a greater efficacy for keliximab in depleting CD4+ cells than clenoliximab. Furthermore, the findings from this preclinical model were in agreement with results from clinical trials at comparable doses.
A population PK/PD model was developed for canakinumab,74 an anti-interleukin-1β (IL-1β) mAb, to examine the pharmacological responses of this drug in RA patients. The PK of total canakinumab was well captured with a standard 2-compartment model with a linear absorption process from a subcutaneous depot compartment, whereas drug binding to the IL-1β was well characterized with a quasi-equilibrium TMDD model. Normalized plasma concentrations of free IL-1β served as a stimulatory driver for the production of the C-reactive protein (CRP), a common biomarker for inflammation, which was described with a series of 3 transit compartments to account for the time delays in the CRP dynamical response. The final model for the graded American College of Rheumatology (ACR) scores used a modified model connecting a continuous latent variable, termed ACR Latent (ACRL), to multiple related binary outcome variables.75 The disease condition was modeled using a single transit compartment with identical first-order production and loss rate constants, and the difference of free IL-1β from its baseline was incorporated into a sigmoidal Emax function to indirectly represent the inhibitory effect of the drug on disease processes. The latent variable was related to the probability of achieving 20, 50, or 70% improvement from baseline conditions using a logit transform. The main conclusions from this canakinumab model were drawn from model-based simulations, which confirmed that 150mg every 4 weeks improved the ACR scores in patients with RA, but there is no additional benefit to the patients treated with higher doses or more frequent dose administration.
Furuya et al.76 developed a PK/PD model for infliximab, a chimeric murine-human monoclonal IgG1 antibody that targets tumor necrosis factor (TNF) and is used as a therapeutic agent for CD. The PK of infliximab was described with a TMDD model. The final model linked unbound TNF concentrations to the effect of interest, and it was assumed that the antibody-ligand complex was cleared with the same fractional catabolic rate as the unbound ligand (TNF). However, model fittings suggested a greater estimate for the plasma half-life of TNF (~40 days) as opposed to the physiological value of < 1h. Overall, the model captured the time-courses of infliximab serum concentrations and the changes of the CD activity index (CDAI) reasonably well. Furthermore, the model helped establish a rationale for individualized dosage regimens for infliximab in patients with CD. The value of the CDAI obtained after the first dosing event was used to generate individual patient PD parameters and to determine an appropriate dosage for subsequent drug administration.
Deng et al.77 developed a PK/PD model that links the exposure of an anti-platelet antibody, anti-CD41, to the time-course of platelet counts in a mouse model of sustained ITP. The model combines the effects of intravenous IgG (IVIG) on the clearance of anti-CD41 and platelet counts in anti-CD41-induced ITP mice. The stimulatory effects of anti-CD41 on the removal of circulating platelets were captured with an IDR model coupled with a series of transit compartments mimicking a cascade of events preceding platelet production (i.e., megakaryopoiesis) to account for the observable time-delay in the PD response (i.e., increase in platelet counts). A feedback-regulation loop from circulating platelets on their own production was integrated in the model to capture the apparent tolerance effect. The IVIG treatment enhanced anti-CD41 clearance through saturation of FcRn and decreased the elimination of opsonized platelets. The model suggested that the PK effect of IVIG (i.e., increasing anti-CD41 elimination) is responsible for 43% of the overall effect of IVIG on anti-CD41-induced ITP.
Abciximab (anti-glycoprotein IIb/IIIa) and daclizumab (anti-IL2) are examples of US Food and Drug Administration (FDA)-approved and marketed antibodies that achieve their effects by altering cell signaling pathways. Abciximab is an antigen-binding fragment that targets the GPIIb/IIIa receptor on platelets. Upon its binding to the receptor, abciximab inhibits platelet binding to fibrinogen and von Willebrand factor, thereby inhibiting platelet aggregation. An inhibitory Emax model78 was developed to describe the effects of abciximab on ex vivo platelet aggregation in patients undergoing coronary angioplasty. The PD model assumed that abciximab concentrations in plasma were directly related to its inhibitory effects. The model provided a good characterization of the data, and serves as a rare example of applying a simple direct effects PD model to antibody PD.
Diao L et al.79,80 developed a population PK/PD model for daclizumab in healthy volunteers and multiple sclerosis patients. The model integrated the PD of 3 biomarkers, including CD25 (or IL2 receptor) occupancy, CD56bright natural killer cell counts, and regulatory T cell (Treg) counts from 4 clinical trials. Data were analyzed using non-linear mixed effects modeling, and an Emax model captured CD25 saturation temporal profiles, which occurred very shortly after drug administration. An IDR model captured well the cellular PD responses. The main finding from this model was the that average maximum reduction in Treg cells of about 60% occurred ~4 days after a 150 mg subcutaneous dose. After the last dose, Tregs were projected to return to their initial baseline levels in ~20 weeks.
To illustrate the relevance of systems biology to understanding mechanisms of action of mAbs and their downstream pharmacological effects on diseases, we reviewed several examples, including recombinant immunotoxin (RIT) in solid tumors, etanercept in RA, PTEN involvement in resistance to trastuzumab treatment, HER2 signaling pathway in de novo resistance to trastuzumab treatment, cetuximab in metastatic colorectal cancer, and rituximab signaling pathway in non-Hodgkin's lymphoma.
The interest in RIT as an effective mAb-based cancer therapy stems from their ability to bind to membrane proteins that are unique to malignant cells, subsequently resulting in tumor regression. Chen at al.,81 developed a mathematical model that incorporates the so-called “binding site barrier concept” to explore RIT biological features on overall drug efficacy. The binding site barrier concept corresponds to the phenomenon of competition between binding and diffusion of mAbs. Thus, affinity that is too great between the RIT and its pharmacological target may lead to a reduced cytotoxicity. Very little is currently known about the key factors that maximize RIT antitumor activity, such as the optimal dose range for RIT to achieve the desired binding affinity, extracellular diffusivity, and target binding site density. The proposed mathematical model explores the dose-dependent cytotoxic activity of Pseudomonas exotoxin (PE)-based RIT against either A431 carcinoma or CD46 Burkitt's lymphoma cell xenografts, and employs the previously described QE TMDD model, with additional biological processes such as recycling of receptors to the cell surface membrane, random sorting and translocation of the toxin to cytosol, and the de novo synthesis of receptors. A naïve pooled data analysis was used and yielded the estimation of 38 parameters. A sensitivity analysis was conducted and revealed that tumor response decreases with the size ratio of the tumor to the blood vessel. In addition, the agreement between model simulations and experimental data of human carcinoma tumors in xenograft mice reinforced the binding site barrier concept, and showed that bidirectional transport across the capillaries reduces the amount of free RIT within the tumor extracellular space. The optimal antitumor activity is achieved through a logarithmic dependent relationship between the binding affinity and the target binding site density.
RA is a complex systemic, chronic, inflammatory and multifactorial autoimmune disease. A range of system-related genes have been identified in the pathogenesis of RA; however, a comprehensive understanding of how the disorder develops is still not well established. Thus, identifying the structure and dynamics of molecular networks is a key concept to unraveling such a complex disease. Computational network analyses may support the investigation of biological features of RA, and also guide optimization of its treatment. Etanercept is an Fc fusion protein for RA that neutralizes TNF and blocks its downstream molecular network. Etanercept therapy has been proven to slow the progression of RA, but its exact underlying molecular mechanisms remain unknown. Hecker et al.82 investigated the therapeutic effects of etanercept on the transcriptional regulation of many genes in RA using a gene regulatory network (GRN) analysis. Genome wide transcriptional profiles (GWTP) using DNA microarrays were obtained for 19 RA patients following their first week of treatment with etanercept. GWTP revealed a set of drug-responsive genes, among which transcriptional factors (TF) and TF binding sites (TFBS) were overexpressed. Mathematical modeling of the underlying GRN was conducted using a system of linear equations, along with least angle regression (LAR) and ordinary least squares regression to explain the controlling factors between TF and TFBS. This integrative modeling strategy is called TFBS-integrating LARS (TILAR). Overall, the TILAR method identified a set of 83 TF-genes involved in RA, among which 48 genes were downregulated and 35 genes were upregulated. Most of these genes are known to control immune responses. The constructed GRN model included 95 nodes and 106 edges, and included 22 parameters. The TILAR method was able to derive gene regulatory interactions from gene expression data by integrating TF and TFBS predictions. Furthermore, this technique distinguishes between genes, TF, TFBS, and identifies the connection between them, which provides testable hypotheses of etanercept mechanisms of action.
Trastuzumab targets the HER2 receptor tyrosine kinase (RTK), which is clinically effective in HER2-positive breast and ovarian cancers.83,84 However, resistance to trastuzumab is a common clinical problem, mainly because of the impaired expression of HER2 and the over-activation of cell signaling pathways.85 Thus, there is a considerable need to identify novel biomarkers and therapeutic targets to overcome trastuzumab resistance and predict therapeutic response. In this section, we review 2 mechanisms of resistance to trastuzumab treatment from a computational systems biology perspective: 1) the role of tumor suppressor protein phosphatase and tensin homolog (PTEN) in the mitogen-activated protein kinase (MAPK) and phosphatidylinositol 3-kinase (PI3K) signaling pathways,86 and 2) the overexpression of the transmembrane tyrosine kinase HER2.87
A dynamic mathematical model of a MAPK/PI3K signaling network developed by Faratian et al.86 integrates the binding of trastuzumab to the HER2, HER2/HER3 dimerization and inhibition, cross talk of AKT and MAPK pathways, and the dynamic and regulatory features of PTEN. Using experimental dynamic parameters extracted from literature, the model satisfactorily described the competitive processes between lipid and protein phosphatase activities of PTEN, its auto-phosphorylation, and the equilibrium between its active (PTEN) and inactive (pPTEN) forms. Model fittings to in vitro generated data from ovarian and cancer cell lines dosed with pertuzumab (an RTK inhibitor similar to trastuzumab) thoroughly characterized the dynamic responses of HER2/HER3 heterodimerization, phospho-ERK (pERK) and pAKT activation, and the inhibitory effects of pertuzumab. Furthermore, a resistance factor (γ) derived from the model indicates the balance between the activities of PTEN and PI3K and governs the sensitivity of AKT activation of pertuzumab, such that high γ values predict a sustained inhibition of AKT stimulation, and low γ values predict the opposite. Computational results from the proposed model predicted that signal transduction through the PI3K/PTEN cycle is sensitive to inhibition by pertuzumab and that pAKT signaling is decreased by 30–70%. Therefore, the model suggests that resistance to RTK inhibitors is dependent on the quantitative balance between PTEN and PI3K, which are both commonly deregulated in breast and ovarian cancers.88,89
HER2 is an overly activated protein in breast cancer patients. Trastuzumab disturbs downstream signaling events resulting in cell cycle arrest and cytotoxicity, but the response rate to trastuzumab is relatively low (<30%) and over half of patients treated with the drug develop resistance. Sahin et al.87 constructed a protein interaction network in which they connected the HER2 signaling pathway with the G1/S transition of the cell cycle through 2 critical transcription factors. Their modeling strategy consisted first of the development of a Boolean logical framework of the biological system by gathering information from literature. The resulting network consisted of 18 proteins, including the RTK insulin-like growth factor 1 receptor, 2 transcription factors (ER-α and c-Myc), signaling factors (AKT and MEK), and G1/S transition cyclins CDK and CDK inhibitors. Second, in silico simulations of protein function loss were performed, with each protein in the network considered either active (value 1) or inactive (value 0), and resulted in the prediction of significant decreases in phosphorylated RB (pRB), the marker for G1/S transition for Cyclin E1, ER-α, and c-Myc knockdowns. Third, model predictions were further confirmed by data from a series of in vitro experiments involving quantitative measurements of protein abundance and activation states using reverse phase protein arrays. Experiments indicated that c-Myc, ER-α, and Cyclin D1 were involved in de novo mechanisms of resistance to trastuzumab. Moreover, targeting c-Myc with a specific inhibitor alone or concomitant with trastuzumab also resulted in a strong reduction in pRB. Experimental results were further confirmed using model predictions that ultimately selected c-Myc as a candidate to be tested in in vitro and in vivo models for future treatments of breast cancer.
Cetuximab (CTX) targets epidermal growth factor receptor (EGFR) and is commonly used in the treatment of metastatic colorectal cancer (mCRC), in which mutations of the proto-oncogene KRAS gene are predictive of poor prognosis.90 Understanding the mechanisms underlying response failure and acquired resistance to CTX treatment is critical and may significantly improve the clinical management of mCRC. Oliveras-Ferraros et al.91 studied the transcriptomic signature triggered by CTX in a KRAS wild-type A341 tumor cell line. EGFR pathway activation and long-term inactivation of EGFR/RAS/MAPK pathways were investigated. A gene set enrichment analysis screening of the Kyoto encyclopedia of genes and genomes pathway database was used to assess genome-wide analyses of whole human genome arrays. Functioning of the CTX transduction pathway was shown to depend upon several processes, including activation of the EGFR pathway, the regulation of the MAPK pathway activation by dual-specificity phosphatases (DUSP), and the transcriptional status of gene pathways controlling the epithelial-to-mesenchymal transition (EMT) and its reversal (MET) program. Quantitative real-time PCR, flow cytometry analysis, and high-content immunostaining confirmed that CTX efficacy relies on 3 factors: 1) its ability to increase cell-to-cell contact by upregulating the expression of the epithelial markers E-cadherin and occluding; 2) downregulation of epithelial transcriptional repressors, such as Zeb, Snail, and Slug along with the restoration of cortical F-actin; and 3) complete suppression of the CD44pos/CD24neg/low mesenchymal immune-phenotype. As a result, the prediction of CTX efficacy in mCRC from gene transcripts of EGFR ligands/MAPK phosphatases relies on the ability of CTX to concomitantly reverse the EMT status known to be an important characteristic of migrating cancer stem cells.
The delta-like ligand 4 (DII4)-mediated Notch signaling pathway is essential for angiogenesis. Dll4 is a transmembrane protein acting downstream of vascular endothelial growth factor (VEGF) signaling and has shown substantial negative regulatory effects on VEGF-induced angiogenesis. Yan et al.92 reported that antagonizing DII4 via an anti-DII4 mAb results in tumor regression, especially in tumors that are resistant to VEGF inhibition. The proposed mechanism of action is that the newly developed tumor microvasculature after anti-DII4 treatment is highly dependent on VEGF, the removal of which results in apoptosis of endothelial cells and vascular regression preventing tumor progression. Thus, Dll4 blockade results in the formation of a nonfunctional vasculature that is unable to support tumor growth.
Using specific mAbs, the Bcr-Abl protein interaction network was characterized in CML K562 cells by Brehme et al.93 Among the 18 protein candidates that were identified to highly interact with Bcr-Abl, statistical analysis revealed that only a set of 7 protein “cores” were tightly connected around Bcr-Abl, including Grb2, Shc1, Crk-I, c-Cbl, p85, Sts-1, and SHIP-2, as well as their relationships to other signaling pathways. In addition, inhibition of Bcr-Abl kinase activity by tyrosine kinase inhibitors, such as nilotinib and desatinib, resulted in complete impairment of some protein interactions and suggested AP-2 as an adapter complex mediating intracellular receptor tyrosine kinase trafficking linked to Bcr-Abl.
Rituximab is the first mAb approved as a treatment for indolent B-cell non-Hodgkin's lymphoma (NHL). The use of rituximab, both as a single therapeutic antibody or together with chemotherapy, has advanced the treatment of patients with B-NHL. Rituximab inhibits survival pathways (MAPK, ERK-1/2, NF-kB, and Akt signaling pathways) in B-NHL cells after binding to CD20, which is an unglycosylated phosphoprotein specifically expressed within B-cells.94-96 Furthermore, rituximab treatment activates the caspase cascade through YY and FAS as a result of the inhibition of the survival pathways.97 Understanding the complexity of rituximab-mediated cell signaling and its effects are essential to guide the development of novel therapeutic interventions, as well as the identification of biomarkers that can help prognosis/diagnosis and improve treatment. In the following section, we review 2 modeling approaches that may be used to identify potential intervention targets and biomarkers: 1) Boolean networks and 2) quantitative system pharmacology models.
Boolean network analysis is a powerful tool to understand and predict system behavior based on qualitative data. They are qualitative representations of biological systems in which nodes are different components of the system and edges represent interactions and relationships among nodes (Fig. 3a). The nodes are characterized by either active (1 or ON) or inactive (0 or OFF) states. The relationships among nodes are described by either activation or inhibition. For example, A → B (activation) illustrates that if A is active (ON) and B is inactive (OFF) at timestep (t), A will remain in the active state and B will be active (ON) at the next timestep (t+1). Similarly, A –| B (inhibition) illustrates that if A is inactive (OFF) and B is inactive (OFF) at timestep (t), B will remain in an inactive state and A will be active (ON) at the next timestep (t+1). When there are multiple nodes affecting another node, the relationship can be described by either (AND) and (OR) statements. For example, both Caspase 8 and p53 should be (ON) in order to activate BID (Fig. 3a). In the absence of any of the 2 components, the activation of BID does not take place. This activation is defined by an (AND) statement. Odefy, a publicly available MATLAB→ toolbox that transforms Boolean relationships to ordinary differential equations, was used to simulate this simple Boolean network. The transformed model is described by 3 parameters (Fig. 3b), and details of Odefy and the transformed model have been described elsewhere.98 Literature data from Ramos cells94-97,99 were used to develop Boolean relationships that describe RTX-mediated cell signaling effects on survival pathways (Fig. 3a). In addition, the Boolean relationships of the caspase cascade were adapted from the literature.100 The steady-states of the nodes were considered in order to validate the Boolean network, which can be used to address potential intervention targets and identify biomarkers.
Algorithms for identifying both intervention points and biomarkers of the Boolean network are summarized in Fig. 3c. Intervention targets of the network were defined as the nodes in which inhibition or stimulation might result in the desired outcome. In the present case, the desired outcome was an active steady-state of DNA damage event in the absence of RTX. The intervention strategy may be constant inhibition or constant stimulation of a protein or node. Both intervention strategies can be modeled in the Boolean network by manipulating corresponding Boolean relationships and initial conditions. The initial condition of a given node was set to 0 or 1 for constant deactivation or constant deactivation, then the change of a node state over time was set to 0 so that a new network is created. Next, the new network is simulated and the steady-state of the output (DNA damage event) was compared to the steady-state of output of the original network. Fig. 3d illustrates the change of states for selected nodes in the original network. If the output of the new network reached active steady-state, then the node associated with the manipulated Boolean relationship and initial condition was considered significant. This node could be set as a potential target for combinatorial therapy. In addition to intervention points, biomarkers of a biological network are defined as the nodes that are crucial for network robustness. To identify biomarkers of a biological network, a given node is deleted (if there is no self-loop), and the Boolean relationships are recreated by following an algorithm described elsewhere.101 The main point of reducing a Boolean network is to preserve the wiring diagram and the dynamical properties of the system. When a node is deleted, the deleted node and its effect on other nodes of the network are replaced with the relationships that affect the deleted node. For example, when b is deleted in (a–| b → c), the recreated network is the following: (a–|c). After recreating and simulating the new network, the steady-states of all the nodes in the new network are compared to the original network. For a node to be significant, the steady-states of the original network could not be preserved in the new network.
The results of a simple intervention point and biomarker analysis for a simple network of RTX-mediated cell signaling102 are summarized in Fig. 3e. The inhibition of the NKFB-IKK-IKB protein set (modeled as constant deactivation in the Boolean network), YY, XIAP and the stimulation (modeled constant activation in the Boolean network) of FAS, Caspase8, Caspase3, Caspase 6, and APC were identified as potential intervention targets for combinatorial therapy. The identified intervention targets are biologically relevant. For example, rhApo2L/TRAIL was hypothesized to activate FAS and it has been shown that the combination of rhApo2L/TRAIL and RTX therapy attenuates xenograft tumor growth.99 Another identified potential target for combination therapy is the NKFB-IKK-IKB protein set, which is considered a significant target for anti-cancer agents.103 In addition, specific components of the caspase cascade are identified as intervention points. The biomarkers of RTX-mediated cell signaling were identified as RKIP, NIK, YY, p53 and Mitochondria. RKIP and NIK are upstream of RTX-mediated cell signaling, and the change in both RKIP and NIK is essential to RTX treatment. Although p53 is a widely studied cancer related protein, its involvement in Ramos cells and RTX-mediated cell signaling needs to be studied further and experimentally verified.
Harrold et al.102 derived an integrated QSP model for RTX and a small CD20 signaling pathway in NHL. The model accounted for critical factors determining anti-tumor efficacy of RTX as a single agent and in combination with fenretinide or recombinant human Apo2 ligand (rhApo2L). The model incorporated key regulatory mechanisms induced by CD20 occupancy, such as 1) the target-mediated disposition of RTX, 2) the regulation of apoptosis through the CD20 intracellular signaling, and 3) the relative efficacy of death receptor isoforms. The tumor responses from each agent were linked to their mechanism of action in vivo, and the PD temporal changes of some critical proteins, e.g., Bcl-xL, Fas, were linked to these in vivo responses to explain the apparent synergy of these drug combinations. Predictions of tumor progression profiles using the developed multiscale model agreed well with cell-based measurements and tumor growth kinetics in mice bearing xenograft tumors. The synergistic anti-tumor responses of these agents were captured with high fidelity, thus providing a mechanism-based platform for exploring new regimens with CD20 agonist combinations.
There are a large number of published mechanism-based PK/PD models for mAbs that present insightful descriptions of the drug and the biological systems of interest. However, their application may be limited by experimental design and to the specific drug. Such models may lack generalizability and are rarely sufficient for understanding complex clinical phenotypes. Systems pharmacology models may overcome these complexities by assembling together, in a quantitative framework, diverse sources of prior knowledge from reductionist and traditional biological experiments. The inclusion of more mechanistic components might impart improved projections of system behavior under different conditions and perturbations. The use of mathematical models is critical to the development of successful mAbs therapeutics, and they are needed to address scientific questions arising from the drug or the system. The combinatorial complexity of new antibody-based drugs in drug discovery is formidable, and it is not feasible to explore such complexity by experimentation alone. Two of the biggest challenges in developing translational mathematical models for mAbs, especially in complex diseases such as cancer and some autoimmune diseases, are the seamless integration of big data and making them readily accessible for exploratory analyses.
Advances in mathematical modeling of mAb drug action primarily emerge from the combination of 4 factors: 1) improved analytical methodologies, which allow for the measurement of nanomolar concentrations of biological entities at the cellular and molecular levels; 2) advanced computational approaches and tools, including computer hardware and software; 3) increased interest in delineating mechanisms underlying physiological and pathological processes of complex systems to provide insights into designing novel drug candidates and delivery systems to optimize exposure-response relationships and therapeutic efficacy; and 4) an increased scientific knowledge base. Successful translational models should be useful for predicting mAb drug action across phases of development, and to explain interspecies differences in drug effects. Furthermore, strategies are needed to involve sophisticated systems pharmacology models, often developed in early phases of drug development in a complementary manner with conventional PK/PD models of mAb mechanisms of action. Future efforts should focus on drug safety, identifying and validating novel biomarkers and surrogates of drug effects, resistance mechanisms, and novel combinatorial regimens.
No potential conflicts of interest were disclosed.