A mechanistic, multiscale mathematical model of immunogenicity for therapeutic proteins: part 2-model applications.

A mechanistic, multiscale mathematical model of immunogenicity for therapeutic proteins was built by recapitulating key underlying known biological processes for immunogenicity. The model is able to simulate immune responses based on protein-specific antigenic properties (e.g., number of T-epitopes and their major histocompatibility complex (MHC)-II binding affinities) and host-specific immunological/physiological characteristics (e.g., MHC-II allele genotype, drug clearance rate). Preliminary validation was performed using mouse studies with antigens such as ovalbumin (OVA) or OVA-derived peptide. Further, using adalimumab as an example therapeutic protein, the model is able to simulate immune responses against adalimumab in individual subjects and in a population, and also provides estimations of immunogenicity incidence and drug exposure reduction that can be validated experimentally. This is a first attempt at modeling immunogenicity of biologics, so the model simulations should be used to help understand the immunogenicity mechanisms and impacting factors, rather than making direct predictions. This prototype model needs to be subjected to extensive experimental validation and refinement before fulfilling its ultimate mission of predicting immunogenicity. Nevertheless, the current model could potentially set up the starting framework to integrate various in silico, in vitro, in vivo, and clinical immunogenicity assessment results to help meet the challenge of immunogenicity prediction.

The development of unwanted immune responses (immunogenicity) against a therapeutic protein, particularly the induction of antidrug antibodies (ADA), involves complex immunological mechanisms, such as antigen presentation, activation of immune cells, and cytokine production. Understanding immunogenicity of therapeutic proteins involves assessing their antigenic properties, while putting them into the context of the host's immunological environment. To predict immunogenicity, various techniques have been developed to assess their antigenic potentials. For example, in silico prediction tools are available for predicting the T-cell or B-cell epitopes based on protein sequences or structures.[1,2,3,4,5,6,7] Experimental approaches, such as in vitro major histocompatibility complex (MHC)-peptide binding assays,[8,9] T-cell proliferation assays[10,11] and humanized mice,[12,13] are being explored to assess the immunogenicity risk. Due to the complicated mechanisms for immunogenicity and the large number of impacting factors, it is often difficult to quantitatively integrate results for immunogenicity prediction.

Mathematical modeling may serve as a helpful tool for this purpose, since it can quantitatively recapitulate complicated mechanisms and incorporate the effect of multiple influencing factors. By mathematically describing the current knowledge of immunogenicity development, a multiscale, mechanistic model was developed. While many mathematical models were developed to describe immune system dynamics, none of them were applied to the development of immunogenicity in a therapeutic setting.[14,15,16] We developed a multiscale model of immunogenicity, described in detail in a companion report (Part 1). The current model is inherently compatible with parametric inputs informed by experimental results that correspond to various impacting factors for immunogenicity. For example, the model includes antigen presentation, during which the processing of antigenic protein into T-epitopes, and the binding between T-epitopes and MHC-II, take place. This model component allows for the integration of protein-specific information, particularly the number and MHC-II binding affinities of T-epitopes, which can be obtained through in silico or in vitro experiments. This component also permits the incorporation of patient-specific information, such as MHC allele genotype, which is known to be a crucial factor for the immune response. Many other potential impacting factors for immunogenicity, e.g., initial number of naive T and B cells and number and binding affinity of B-cell epitopes, are designed as integral parts of the model structure; these can also be conceivably informed by conducting appropriate experiments.

In this work, we applied the mathematical model to the simulation of in vivo immune response in mouse and human using selected case studies. The model is able to simulate immunological responses to therapeutic proteins based on protein-specific characteristics (e.g., T-cell epitope, B-cell epitope) and host-specific characteristics (e.g., MHC-II genotype). Model simulations include kinetics of immune cells, antigenic protein and ADA profiles, antibody affinity maturation profile, etc. Importantly, when certain population characteristics, e.g., MHC-II allele frequency, are known, the model can ultimately be used to simulate immunogenicity incidence within that population.

Results

Simulation of immune response against OVA in mouse

A preliminary model validation/data fitting was performed using two mouse studies monitoring immune responses against an immunogenic protein, ovalbumin (OVA), or OVA-derived peptide. Simulations of mouse immune response overlaid with experimentally determined data are illustrated in ,. In the first study, by injecting OVA323–339, a well-known T-epitope peptide in OVA, significant T-cell response was elicited, with a dramatic increase of total T-cell number, and the generation of a large number of memory T cells.[17] Using parameters specific to the antigen (OVA323–339) and the host (C57BL/6 mice), e.g., dose and MHC-II binding affinity, the model simulation was reasonably consistent with the experimental results.

The second study measured the total plasma cell number after secondary immunization with OVA protein. To best describe the experimental data, we adjusted two parameter values, including g2 (percentage for activated B cells to differentiate to short-lived plasma cells) and CCN (carrying capacity for a functional T cell to stimulate the activation and proliferation of target naive B cells), which are not available from the literature. The parameter g2 is necessary to account for the fact that there are both short-lived and long-lived plasma cells for ADA secretion.[18] Another parameter, CCN is a parsimony parameter for modeling the activation of B cells via interacting with active T-helper cells.[19] By using the values of g2 = 0.4, and CCN = 10, the model simulations are in reasonable agreement with the literature results, thus increasing our confidence in model structure and parameter values.

Simulation of immune response against adalimumab in 1,000 North American subjects

One potential application of the mechanistic model is the simulation of immune responses against therapeutic proteins in a human subject, and eventually in a heterogeneous human population. Adalimumab (Humira), one of the anti-TNF-α treatments for various inflammatory and autoimmune diseases, is a fully human IgG1 monoclonal antibody that was reported to induce varying degree of ADA responses in a subset of patients.[20,21,22]

The immune response profiles for one simulated subject are illustrated in as an example. This subject carries an MHC-II allele (DRB1*04:11) with strong binding affinity (57.44 and 101.5 nmol/l) to the two predicted T-epitopes. Under the clinical regimen, this virtual patient gradually developed ADA against adalimumab, resulting in the reduction of drug exposure at later time points (). During antigen presentation (,), the T-epitopes are efficiently presented onto DRB1*04:11 (binder), resulting in significant T-cell activation and proliferation (). Subsequently, naive B cells are efficiently activated upon strong T-cell help (), leading to the generation of a high number of plasma cells, which are responsible for secreting ADA. The dynamics of the antibody affinity maturation is also captured by the model (,), showing that ADA with higher antigen-binding affinities (clone 7 to 10, dashed lines) are preferentially produced over time, leading to a gradually increased average binding affinity.

A more powerful potential application of the model is the simulation of expected immunogenicity incidence in a human population against a therapeutic protein. This is illustrated by simulating the immune response against chronically dosed adalimumab in 1,000 North American patients (). These virtual subjects carry different MHC-II alleles, whose frequency distribution follows the NCBI MHC database record. Another patient-specific variable illustrated here is the elimination rate of adalimumab, whose distribution was obtained from a population pharmacokinetics (PK) study. As discussed in Part 1 of this paper, the genetic background (MHC-II genotype) may have a strong impact on immunogenicity; therefore, the patients can be stratified according to the number of strong T-epitope-MHC-binding pairs. For example, for a patient with DRB1*04:01 and DRB1*04:03, the two T-epitopes can bind tightly to both 04:01 (123 and 85 nmol/l) and 04:03 (78.52 and 147.85 nmol/l). Therefore, the number of T-epitope-MHC pairs for this patient is 4. By stratifying patients using the number of binding T-epitope-MHC pairs, a significant impact of T-epitope on ADA response can be observed. , provide an overview of the projected ADA and adalimumab concentration for the 1,000 subjects at the end of simulation. In , there is a statistically significant difference in the ADA concentration between subjects who carry 0, 1, or 2 T-epitope-MHC binding pairs. With more than three pairs, the ADA response appears to reach a plateau. The higher ADA response in group 3 is not statistically significant, since the patient number is very small (n = 6, dictated by the MHC allele distribution frequency), and the variance is thus large. For Ag (adalimumab) concentration (), when the subject carries at least 1 binding pair, Ag concentration is significantly decreased. The time profiles of ADA and adalimumab are presented in ,.

The estimation of immunogenicity incidence, particularly the percentage of ADA+ patients, can vary significantly depending on sensitivity and drug tolerance of the ADA assay. We arbitrarily set up a threshold for the absolute ADA concentration (250 ng/ml, the sensitivity recommended by the US Food and Drug Administration),[23] and a threshold for drug tolerance (molar ratio of ADA over Ag at 1). A patient would be reported as ADA+ only when the two thresholds are exceeded. The predicted time course of ADA+ incidence () indicates that a majority of this population would develop ADA after 3 months of treatment, with a suggested immunogenicity incidence of 75.3% after 574 days of treatment. To investigate the reduction of drug exposure due to ADA emergence, cutoff values for adalimumab trough concentration were chosen at 1/2, 1/5, 1/10, 1/20, and 1/50 of the average drug concentration in patients who carry no binding pairs (and thus do not develop ADA). The percentages of patients with drug concentration lower than the cutoff values were plotted in . The model predicts that up to 70% of the patients can have a twofold reduction, and 26.4% of the patients can experience severe (50-fold) reduction of the drug trough concentration after 574 days.

A preliminary comparison for adalimumab trough concentrations under chronic dosing was conducted by replicating a published clinical trial using model simulation[24] (). The overall trend for lower drug concentrations in ADA+ patients compared to ADA− patients is in agreement between model simulation and the trial observation, despite differences in absolute values. The trend for decreasing drug exposure over time in ADA+ patients is also captured by the model.

Discussion

Successful prediction of immunogenicity in human is still mostly in its infancy, due to the intrinsic complexity of the immune system, lack of translation from preclinical species to human, patient heterogeneity, and product-related factors such as aggregates and adjuvant-like contaminants. Although many platforms have been established for immunogenicity assessment, they usually provide partial information regarding one or two impacting factors.[25,26,27,28]

Mathematical models are naturally suitable for describing complex systems, integrating information from various sources, and eventually generating simulations or predictions that cannot be readily and intuitively processed. In this manuscript, we have used case studies to illustrate potential applications for a multiscale, mechanistic immunogenicity model, including simulations for immune responses and population immunogenicity incidence.

In the first case study, by using OVA-specific T-epitope information and mouse strain-specific MHC alleles, the simulated kinetics for T cell and plasma cell are in reasonable agreement with the experimental data, provided two parameters (whose values could not be obtained from the literature) are adjusted. Considering the model contains 88 parameters, the current framework is probably reasonable as a starting point. A future plan would be to perform controlled animal experiments with model antigenic proteins, to generate more immune response data, e.g., immune cell kinetics, antigenic protein PK, ADA kinetics, which can then be compared with model simulations. Upon the availability of more informative datasets, the model parameters can be calibrated by data fitting, thus increasing confidence in the model.

The second case study, a simulation of immune responses against adalimumab, aims to illustrate the potential of the current model to integrate antigen- and patient-specific information for projecting a population response. Antigen-specific information, particularly number of T-epitopes and their binding affinities to MHC, are incorporated. Patient-specific information includes the MHC allele population frequency and the population distribution of antigen elimination rate (drug PK). illustrates immune response simulations in a patient who carries an MHC that binds strongly to the T-epitopes. By simulating 1,000 patients, population-level immune responses against adalimumab were obtained (). The model suggests that T-epitopes have a strong impact on the ADA response, since more T-epitope-MHC binding pairs generally lead to higher ADA response, until a plateau is reached (). ADA development has a significant impact on drug level, as the patients who developed ADA have significantly lower drug concentration (), and this result agrees with clinical observations.[24]

By summarizing the population ADA responses, an immunogenicity incidence of adalimumab in the North American population was calculated to be 75.3% after 19 months of treatment (). Most of the patients are predicted to become ADA+ after about 3 months. Compared to the immunogenicity incidence observed in clinical trials (5–89%),[20,21,22] our simulation result lies at the higher end. Several potential confounding factors may explain the large variation in immunogenicity incidence between clinical studies, and between the clinical observations and model simulations. From a technical perspective, the sizes, duration, and dosing regimens of the clinical studies, together with the ADA assay format, its sensitivity and drug tolerance, may potentially affect the immunogenicity assessment. The inaccuracy in Immune Epitope Database (IEDB) prediction for T-epitopes and their binding affinities may also account for errors in the model prediction, since the in silico prediction of class II T-epitope is not absolutely correct.[29] On the other hand, patient-related variability, such as population genetic background (e.g., MHC allele frequencies), disease indication, immune status, and comedications, can further complicate the immune response outcome. For instance, certain MHC alleles have shown genetic association to diseases,[30] so that patients with autoimmune diseases (some indications for adalimumab) may have a different distribution for MHC alleles compared to normal patients, which are simulated currently. In addition, immunosuppressive compounds such as methotrexate are known to reduce the immune responses to many therapeutic proteins,[31,32,33] and may also partially contribute to the discrepancy between the clinical trials and our simulation.

The reduction in drug exposure due to ADA emergence was also simulated by the model (). Depending on the disease indication, the reduction in drug exposure may ultimately result in variable degrees of loss of efficacy in the patients. Upon the availability of suitable PK/pharmacodynamic (PD) models for adalimumab's efficacy end point in specific disease populations, the current mechanistic model may be flexibly extended for predicting loss of efficacy due to ADA development.

By replicating a clinical trial with adalimumab using simulation, we obtained qualitative comparison for drug trough concentrations (). For ADA+ patients, drug concentrations in the trial are lower than the model prediction. We speculate the lower drug concentrations could be confounded by the ADA assay, which has low drug tolerance.[24] Only the ADA+ patients who show considerable reduction in drug exposure may be captured by the ADA assay. This is probably why the percentage of ADA+ patients in the trial is much smaller compared to the model prediction (see data labels). For ADA− patients, the model prediction is close to the observations at week 4, but underestimates the measurements at later time points; however, the published data are limited to 67 patients where adalimumab trough exposure was measured after 4 weeks. Overall, despite these differences, the high-level trends are reasonably consistent between model simulation and the measurements. ADA generation greatly reduced drug exposure in ADA+ patients. Considering the presence of many potential confounding factors, e.g., small numbers of patients with ADA measurements and the lack of patient-specific information such as MHC genotype, a qualitative agreement constitutes a reasonable start.

In summary, a mechanistic model for immunogenicity was constructed, and its potential applications were illustrated using case studies. Preliminary model validation was performed using mouse immunization studies. The capability of the model to simulate population responses was demonstrated by simulating immune responses against adalimumab in a North American population as well as in a European clinical trial. By integrating protein-specific antigenic properties and a subset of subject-specific characteristics, the model provides predictions of the relevant immune response features that are subjected to validation, including drug and ADA exposure, as well as immunogenicity incidence. It is worthwhile to point out that the current model is a prototype model, which needs to be further improved and vigorously validated before attempting to predict immunogenicity. The examples are provided to illustrate its potential applications, and to highlight the impacting factors for immunogenicity, such as the MHC genotype, ADA assay sensitivity and drug tolerance, etc.

The main limitation of the model is that it is inherently difficult to validate some of the model predictions, both for mouse and human, and thus this initial development is hypothesis generating. A direct comparison of study end points, such as ADA titer, immunogenicity incidence, or loss of clinical efficacy, against model simulations is not straightforward. For instance, the concentration and affinity of ADA from model simulation cannot directly translate into titer values, due to the heterogeneous nature of ADA and the lack of a reference standard in the titer assay. Our hope is that experimental efforts can be directed to independent model validation by us as well as the broad scientific community.

Model validation using clinical data will be critical. Since many confounding factors can greatly impact immunogenicity incidence, more granular information, such as patient genetic background, immune status, use of comedication, needs to be collected in clinical studies. A suitable clinical trial can then be applied toward model validation by comparing ADA status, drug concentration, etc. In addition, model improvements are expected by incorporating other impacting factors, for example comedication, aggregation, excipient in the formulation, or patient immune status. Upon vigorous validation and improvement, the current model can potentially be applied to integrate preclinical and clinical data, and ultimately aid the risk assessment of immunogenicity for therapeutic proteins.

Methods

The model structure has been discussed in detail in the companion manuscript (Part 1). The current manuscript intends to demonstrate practical applications for the mechanistic model, using some case studies as examples. The software Matlab (The MathWorks, Natick, MA) was used for model implementation and for the simulations described below.

Simulation of immune responses against OVA in mouse

For model validation/data fitting purposes, two sets of experimental data were obtained from the literature on mouse immunization studies, where OVA or OVA-derived peptide was used as an antigen. We chose OVA because it is a well-studied antigenic protein for various animal species, and the antigen-specific parameters (e.g., T-epitope and its MHC-II affinity) can be readily obtained from the literature.

The first study determined the number of T-helper cells specific for the OVA peptide (OVA323–339) before and after immunization with OVA323–339 in mice.[17] As described in the original reference, C57BL/6 mice were given i.v. injection of 50 μg of OVA323–339 plus 5 μg LPS. Mean total number (±SD, n = 2–8) of OVA323–339-specific helper T cells were determined using MHC-peptide tetramers, which are composed of four identical biotinylated OVA peptide:I-Ab MHC molecules complexed to a fluorochrome-labeled streptavidin core. Kinetics of total T-helper cell number after OVA323–339 challenge in of the ref. 17 was extracted using DigitizeIt (ShareIt, Eden Prairie, MN). For model simulation, antigen-specific parameters, including the MHC-II binding affinity and the elimination rate of OVA323–339 in mouse, were obtained from the literature. Since MHC-II genotype is an important determinant for binding affinity between T-epitope and MHC, the MHC-II genotype of C57BL/6 mice, I-Ab, is used. The binding affinity for OVA323–339 peptide against I-Ab was reported to be 400 nmol/l.[34] The elimination rate of OVA323–339 peptide was assumed to be the same as OVA due to lack of available data in the literature, and was calculated to be 2.088 day−1 by fitting the PK data with a one-compartment model[35] using the exponential curve fitting function in Microsoft Excel. The antigen-specific parameters, together with the dosing regimen, were fed to the mouse model for simulation. The numbers of total T-helper cells, including naive, activated, functional, and memory T cells were calculated from the simulation results and compared with the experimental result.

The second study determined the absolute number of OVA-specific plasma cells (PCs) in the spleen and bone marrow following secondary immunization with OVA.[18] Briefly, for primary immunization, BALB/c mice received 100 µg alum-precipitated OVA through i.p. injection. After 3–5 weeks, primed mice were boosted by i.v. injection of 100 µg OVA. The OVA-specific PCs were defined as surface Ig low/intracellular OVA-binding blasts. Data were extracted using DigitizeIt (ShareIt) from in ref. 18. The whole body total number of PCs was calculated based on the equation: total PCs = spleen PCs × 4 + bone marrow PCs. The scalar “4” reflects the fact that the spleen contains about ¼ of the PCs in the body.[36] For model simulation, antigen-specific parameters were obtained from the literature. Since OVA323–339 is the only widely and consistently reported strong T-epitope on OVA, it is assumed that OVA323-339 is the only T-epitope on OVA. The binding affinity between OVA323–339 and I-Ad (MHC-II genotype for BALB/c mice) is 150 nmol/l,[37] and the elimination rate for OVA was calculated to be 2.088 day−1 (ref. 35) The simulated total PC number was calculated by adding up short-lived and long-lived PCs, and was compared against the experimental data. Some adjustment to the model parameters was necessary to account for the data (see Results section).

Simulation of immune responses against adalimumab in a human population

To demonstrate the potential of the model to simulate immune responses in a human population, adalimumab (Humira) was used as a representative therapeutic protein. Although it is a fully human monoclonal antibody, adalimumab elicits immune response in a subset of patients. The model was used to simulate human immune responses against adalimumab under clinical dosing regimen in 1,000 virtual patients in North America. The simulation intends to predict the immune response under chronic dosing; therefore, adalimumab was dosed subcutaneously at 40 mg every 2 weeks (41 doses) in the simulation for 574 days. The simulated immune responses depend on (i) protein-specific antigenic characteristics, e.g., the number of T-cell epitope, their MHC-II binding affinities, and therapeutic protein dosing regimen; and (ii) patient-specific characteristics, e.g., patient-specific MHC-II genotype and individual clearance rate for adalimumab.

To obtain the protein-specific parameters (particularly, T-epitope information) for adalimumab, the IEDB prediction tool for MHC-II (http://tools.immuneepitope.org/analyze/html/mhc_II_binding.html) was applied to predict strong T-epitopes.[1,29] It was assumed that T-epitopes on adalimumab only reside within the complementarity determining region (CDR): adalimumab is a fully human antibody, thus other regions would be conceivably well tolerated by the immune system. We defined the predicted T-epitopes as the top 2% of hits based on IEDB ranking. Two strong promiscuous T-epitopes are predicted, including “AKVSYLSTASSLDYW” on the heavy chain CDR-3 region, and “KLLIYAASTLQSGVP” on the light chain CDR-2 region. The binding affinities of the epitopes against individual MHC-II alleles were then predicted by the IEDB class II prediction tool, by using the IC50 values from SMM_align or NetMHCIIpan method. The two T-epitopes were predicted to bind to some MHC-II DR alleles, but not to any DP and DQ alleles. For MHC-II alleles that are predicted to have no binding, a binding affinity of 4 µM is assumed, same as the binding affinity between endogenous competing peptides and MHC-II.[38] The binding affinity predictions are summarized in .

Patient-specific characteristics currently include population MHC-II allele frequency and population elimination rate for adalimumab. The MHC-II allele frequency for HLA-DRB1 in North America was obtained from the database “dbMHC” (http://www.ncbi.nlm.nih.gov/projects/gv/mhc/). The population distribution of the elimination rate of adalimumab was based on a population PK study,[39] reporting a log-normally distributed half-life with variability estimated at 54.6% coefficient of variation. The PK profile of adalimumab after subcutaneous injection was modeled using a modified two-compartment model (manuscript Part 1), based on a single i.v. dose PK study[40] and the Tmax (time to reach maximum serum concentration, reported to be 131 h) after subcutaneous dosing.[41]

Upon the availability of the protein- and patient-specific information, simulations were conducted for immune responses against adalimumab in 1,000 virtual subjects (). Since the two T-epitopes were predicted to bind to some MHC-II DR alleles, but not to any DP and DQ allele, only specific DR alleles were simulated for each virtual patient, and the DP and DQ alleles were simplified as generic nonbinding alleles with binding affinity of 4 µmol/l. To generate the genotype of the DR alleles for a virtual patient, two DR alleles were randomly chosen from the candidate alleles listed in . The probability of choosing each allele is based on the allele frequency. Once the two DR alleles were determined, the binding affinities between T-epitopes and DR alleles were obtained from the table. Another patient-specific variability is added by randomly sampling elimination rates of adalimumab from the predefined distribution determined from the population PK study.[39] The subject-specific MHC-II binding affinities and elimination rate of adalimumab were then applied for simulating this subject. By repeating the above processes, simulations were performed for 1,000 subjects. All the parameters used for this simulation are in the model code and provided as Supplementary Material online.

Simulation of immune responses against adalimumab in a clinical trial

A clinical trial of adalimumab therapy[24] was used for comparison with the model simulation. The patient number (N = 168) and dosing regimen were replicated in the simulation as in the reference. In the trial, dose escalation was given to overcome the loss of clinical efficacy. We tried to replicate the dose escalation in of the ref. 24 by raising the dose in patients who had the lowest drug concentration at designated time points. The European population MHC-II allele frequency was used since the trial was conducted in Europe. Other variables in the simulation, such as T-epitopes and the population eliminavation rate of adalimumab, were the same as in the previous simulation for a North American population.

Author Contributions

X.C., T.P.H., and P.V. wrote the manuscript. X.C., T.P.H., and P.V. designed the research. X.C., T.P.H., and P.V. performed the research.

Conflict of Interest

X.C., T.P.H., and P.V. are employed by Pfizer. X.C., T.P.H., and P.V. hold stock in Pfizer. Associate Editor P.V. was not involved in the review or decision process for this paper.

Study Highlights