A mechanistic, multiscale mathematical model of immunogenicity for therapeutic proteins: part 1-theoretical model.

A mechanistic, multiscale mathematical model of immunogenicity for therapeutic proteins was formulated by recapitulating key biological mechanisms, including antigen presentation, activation, proliferation, and differentiation of immune cells, secretion of antidrug antibodies (ADA), as well as in vivo disposition of ADA and therapeutic proteins. This system-level model contains three scales: a subcellular level representing antigen presentation processes by dendritic cells; a cellular level accounting for cell kinetics during humoral immune response; and a whole-body level accounting for therapeutic protein in vivo disposition. The model simulations for in vivo responses against antigenic protein challenge are consistent with many known immunological observations. By simulating immune responses under various initial parameter conditions, the model suggests hypotheses for future experimental investigation and contributes to the mechanistic understanding of immunogenicity. With future experimental validation, this model may potentially provide a platform to generate and test hypotheses about immunogenicity risk assessment and ultimately aid in immunogenicity prediction.

With the rapid expansion of therapeutic proteins into an important class of medicines, the issue of unwanted immunogenicity has stimulated much research effort and regulatory attention. The consequences of immunogenicity, in particular the induction of antidrug antibodies (ADA), have the potential to become a serious issue during drug development, due to their impact on drug pharmacokinetics (PK), efficacy, and/or safety.[1] Immunogenicity involves complex biological mechanisms, which could span multiple system scales, from subcellular processing and cellular interaction, to organ and whole-body functions. Although various techniques have been developed to assess the immunogenicity risk of therapeutic proteins,[2,3,4,5,6,7] success in predicting immunogenicity is still not prevalent, due to the involvement of complicated mechanisms and large numbers of impacting factors.

Mathematical modeling may serve as a complementary approach to help understand immunogenicity, since it can quantitatively recapitulate, and especially integrate, complicated mechanisms. Mathematical models for the immune system mainly involve two categories of modeling techniques, differential equations (DEs) and agent-based models. DEs have a long history in modeling the immune system. For example, Bell[8] developed a mathematical model for B-cell clonal selection and antibody production as early as 1970. Recently, the adaptive immune response to influenza A virus infection was modeled.[9] Conversely, agent-based models are a more recent approach and model each entity (e.g., an immune cell) as an “agent,” which adapts its behaviors over time (e.g., movement and differentiation) based on rules that have stochastic components. Some recent examples include ImmunoGrid, an integrated large-scale agent-based model environment to simulate the human immune system,[10,11] C-ImmSim, an agent-based simulator that combines computational immunology with bioinformatics,[12,13] and the Basic Immune Simulator.[14,15] One limitation for agent-based models is that they tend to require larger number of parameters than their DE counterparts, so it is often difficult for sufficient experimental data to be acquired to inform the model.[16] Given the comparatively long experience with DE models, we developed our model using DEs, to minimize the number of required parameters. An added benefit of a DE model is that it can be easily integrated with downstream applications more traditional in drug discovery and development, such as PK/PD modeling.

The objective of this work was to establish a multiscale, mechanistic model that can capture the key underlying mechanisms for immunogenicity against antigenic therapeutic proteins. To focus on the essential model components, while having the potential for modular expansion, this model considers the antigen-presenting cells, CD4+ T helper cells, and B cells as the major immune cells. Since dendritic cells (DCs) are the most efficient antigen-presenting cells,[17] they were chosen to represent all antigen-presenting cells in the model. DC activation could be driven by maturation/danger signals that are either signs of pathogen presence, e.g., endotoxin[18] or by tissue damage upon drug administration. Due to the complexity of this process and the unavailability of many parameters, DC activation was simplified and modeled as being directly driven by endotoxin, particularly, lipopolysaccharide, which is widely used in immunology studies to activate DCs[19] and is known to be present in many therapeutic protein dose forms.[20] Once the DCs become activated, they uptake and process the therapeutic protein, in this context the antigen (Ag), and present the T-epitope from the Ag for subsequent T-cell activation. These processes are collectively called “antigen presentation,” a critical step for efficient activation of the adaptive immune system, which ultimately evokes ADA production and immune memory. Efficient antigen presentation eventually leads to the activation, proliferation, and differentiation of T and B cells, as well as the secretion of ADA that modify the in vivo disposition of Ag. Although B-cell activation can be T cell dependent or independent,[21] the current model focuses on the first, because it leads to more robust antibody response with affinity maturation and isotype switching and is associated with more impactful clinical observations, such as high and persistent antibody titer.[21]

Our model was developed for mouse and human using species-specific parameters. When applied to simulate primary and secondary immune responses, the model was able to reasonably reproduce many immunological phenomena, such as enhanced secondary immune response and antibody affinity maturation. It was also applied to generate hypotheses regarding potential important factors for immunogenicity development, including T-epitope affinity, Ag dose, and naïve T- and B-cell number, and suggested directions for future investigation. We believe that its current structure and features make this system model suitable to mechanistically study the immunogenicity of therapeutic proteins.

Results

Mathematical model

The multiscale model structure can be represented by the subcellular, cellular, and whole-body levels (,,): The subcellular model, capturing antigen presentation processes in mature DCs (), is built upon the work by Agrawal and Linderman.[22] It is assumed that, when the mature DCs are first generated, no antigen presentation has happened, and all free major histocompatibility complex (MHC) molecules are located in the endosome.[23] Antigenic proteins are then endocytosed into the endosomes of DCs and are digested into T-epitope peptides, which subsequently bind to MHC-II to appear on DC surfaces as T-epitope–MHC complexes for T-cell activation. The processing and presentation of endogenous competing proteins by DCs are also modeled, to reflect the fact that endogenous peptides can compete with T-epitopes for antigen presentation. To reproduce the physiology of MHC-II, our model includes six classical human MHC-II (2 HLA-DR, 2 HLA-DP, and 2 HLA-DQ) and two classical mouse MHC-II (H-2A(I-A) and H-2E(I-E) subclasses). The idea to model antigen presentation in detail is motivated by the availability of parameters for these processes, including, for example, in silico T-epitope prediction and experimental measurements of their MHC-binding affinities.[24] These Ag-specific parameters can be fed to the model to simulate Ag-specific immune response. In addition, by capturing the genotype of MHC-II in the model, it is possible to simulate host/subject-specific immune responses based on the MHC-II alleles each subject carries.

The cellular level accounts for immune cells kinetics. At the cellular level (), DCs serve as a connection between the innate and the adaptive immune system, by presenting T-epitopes for the activation of corresponding naïve helper T cells. Once the Ag-specific naïve T cells are activated by the DCs, they proliferate and differentiate into functional T cells to facilitate downstream B-cell activation or give rise to memory T cells that can be immediately activated via interacting with antigen-presenting DCs. Subsequently, antigen-specific naïve B cells, which recognize Ag through B-cell receptor, were fully activated upon stimulations from functional T cells and Ag-bound B-cell receptor. Activated B cells proliferate and differentiate into short-lived and long-lived plasma cells, as well as memory B cells. The plasma cells secrete ADA that modifies antigenic protein elimination, while the memory B cells behave similarly as naïve B cells and immediately react to an antigen challenge.

The whole-body level accounts for therapeutic protein in vivo disposition (). The PK profile of therapeutic proteins is usually modeled by empirical compartmental models, which lump rapidly equilibrated tissues for drug distribution into a single compartment. In the current model, plasma is modeled as the space for the immune cells to reside and for the antigenic proteins to interact with the immune system, by assuming that lymphocyte move fast between lymph organs and the blood. To accurately describe the PK of the antigen, while modeling interaction between the Ag and the immune system, a modified compartmental model was constructed. An extended two-compartment model is illustrated as an example. A peripheral compartment includes slowly perfused tissues, as is usually done. The central compartment generally includes plasma and other rapidly perfused tissues and therefore can be regarded as the sum of plasma volume (Vp) and any extra volume in the central compartment (Vec). It is assumed that plasma and the extracentral compartment equilibrate very fast. By fitting the PK profile of the antigenic protein using this modified compartment model, the rate constants (k13, k31, and kel) and Vec can be estimated (once plasma volume and plasma to extracentral rate constants are assumed). Then, the modified compartmental model can be integrated into the rest of the mechanistic model to specifically model the Ag in vivo disposition.

The model details are provided in the Supplementary Materials, with the following organization:

DC activation and differentiation;

Antigen presentation by DCs;

T-helper cell activation and differentiation;

B-cell activation and differentiation;

ADA in vivo disposition;

Antigenic protein in vivo disposition.

The definitions of all model variables and parameters are provided in Supplementary Tables S1 and S2. By using the same model structure, but species-specific parameters, the model can be used to simulate immune responses in mouse and human. Parameters were collected from the literature as much as possible. More parameters are available for mouse compared with human: when there is no reference for human-specific parameters, the current model assumes that they are the same as the mouse parameters. The software Matlab (The MathWorks, Natick, MA) was used for model implementation and for the simulations described below.

Simulation of immune responses in human against a hypothetical antigen

By inputting antigen-specific parameters for a hypothetical antigenic therapeutic protein, the basic model could simulate the human kinetic profiles for DCs, T cells, B cells, ADA, and average ADA affinity upon the injection of antigenic protein (). In addition, the antigen presentation processes, including the internalization, processing, and presentation of antigenic protein and competing protein, are also captured by the model (Supplementary Figure S1). For this particular antigenic protein, the model predicts the generation of strong immune responses, reflected by the rapid activation and proliferation of immune cells, and the production of high level of ADA. The model also predicts greatly enhanced immune response upon secondary immunization, consistent with the immune memory phenomenon. The average antigen-binding affinity of ADA was also calculated to increase over time, which agrees with the antibody affinity maturation phenomenon observed after immunization.

Simulation of immune responses in human under various immune system conditions

Several potential impacting factors for immunogenicity development, including naïve T-cell number, naïve B-cell number, binding affinity of T-epitopes to MHC-II, antigen dose, were tested for their impact on ADA generation by changing the parameter values over physiologically plausible range. The ADA kinetic profiles and the cumulative amount of ADA under various parameter conditions were summarized in . The naïve T-cell number has a significant impact on the magnitude and the timing of ADA generation. The cumulative ADA response as a function of naïve T-cell number exhibited a bell shape (). In addition, a greater naïve T-cell number generally leads to a faster rise of the ADA (), with the earliest rising time around 40 days. In contrast, naïve B-cell number has less impact on ADA generation (). Within the normal parameter range (from 10–7 % to 101 % of total B cells), the cumulative ADA amounts varied only within twofold. Only when excessive number of naïve B cells are present (10–100% of total B cells), the amount of ADA drops dramatically. The binding affinities for T-epitopes (dissociation constant, Kd) demonstrate a strong impact on ADA magnitude, with the most sensitive range of 1 nmol/l to 1 µmol/l (). Ag dose also exhibited a bell-shaped relationship with the amount of ADA generated (,).

Sensitivity analysis

A sensitivity analysis for the Ag amount and the total ADA amount was performed against model parameters. The top 10 sensitive parameters are listed in . Six parameters are shared between Ag and ADA, including rate constants for the proliferation and death of activated T cells (ρAT, βAT), differentiation fraction of T cells (f1) and B cells (g1, g2), and elimination rate of the Ag (βAg). The other four parameters that have high sensitivity on ADA response are involved in T-cell activation (δNT), B-cell proliferation and death (ρAB_N, βAB), and T-cell–B-cell interaction (CCN). The other four sensitive parameters for Ag response include Ag PK parameters (Ag0, VP), number of T-epitope–MHC on DCs for half-maximum T-cell activation (KpM,M), and death rate for functional T cell (βFT).

Discussion

By integrating key biological mechanisms for immunogenicity development, a mathematical model was developed for simulating humoral immune responses against antigenic therapeutic proteins. The model can reasonably reproduce many immunological phenomena, such as the enhanced secondary immune response, antibody affinity maturation, and antigen presentation. illustrates the simulation results for immune responses in human against two doses of a therapeutic protein. Notably, upon the secondary immune challenge, the magnitude of the immune response (number of immune cells, amount of ADA) is much higher than the primary response, consistent with the notion of immune memory. The model behavior for enhanced secondary immune response was achieved by the following model features: (i) incorporating memory T and memory B cells, which can get activated, proliferated, and differentiated into functional cells upon antigen stimulation; (ii) assigning lower activation threshold for memory cells compared to naïve cells[25,26]; (iii) assigning higher proliferation rate for activated memory B cells compared to activated naïve B cells.[27] Another interesting result is the time profile of the average ADA affinity (), which gradually increases over time. This model behavior is achieved by treating the B-cell lineages, and the corresponding ADA, as polyclonal populations with different antigen-binding affinities. The model assumes that the activation, proliferation, and differentiation of B cells are dependent on the B-cell receptor occupancy, which is affected by their antigen-binding affinity. Therefore, when the amount of antigen is limited, B cells that have higher receptor-binding affinity to the antigen will be preferentially activated, producing ADA with higher average affinities over time.

After confirming that the model behaviors are reasonable, we also investigated the impact of several parameters of interest on the ADA responses, by changing the parameter values over a wide, physiologically plausible range ().

Four parameters, including the number of naïve T cells, number of naïve B cells, MHC-binding affinity of T-epitope, and Ag dose, were studied in the simulation. Although not routinely assessed for immunogenicity risk, the naïve T- and B-cell numbers were probably important; this would be expected, since more naïve cells likely better prepare the immune system against antigen challenge. Interestingly, the model suggests that a bigger impact on ADA response may come from the number of naïve T cells, rather than naïve B cells (,). The results indicate that a sufficient number of naïve T cells needs to be present to efficiently initiate an immune response, while the initial naïve B-cell number is not as critical. Even a small number of naïve B cells can be activated and proliferate to a sufficient size for mounting ADA response, as long as enough stimulation (T cells and Ag) is present. Studies have suggested that the size of naïve T-cell repertoire can vary dramatically depending on the antigenic T-epitope.[28,29,30] Our result emphasizes the potential importance of measuring naïve T-cell number when assessing the immunogenicity risk of therapeutic proteins. The model simulations also indicate that excessive amounts of naïve T and B cells may cause the expected immune response to decline, rather than to increase (,). This model behavior is caused by competition for resources, such as competing for DC help among T cells or competing for T-cell help and for antigen stimulation among B cells. Under the current simulation settings, frequencies to achieve the peak response are around 102 per million for naïve T and around 104 per million for naïve B cells.

MHC-binding affinity of T-epitope also demonstrated strong impact on the ADA response (,). The simulation shows that higher MHC-binding affinity generally leads to greater ADA response, because of the higher number of peptide-MHC presented on DCs. Under the current settings, the sensitive affinity range for T-epitope is between 1 nmol/l and 1 µmol/l (dissociation constant, Kd), suggesting that T-epitopes with Kd weaker than 1 µmol/l are unlikely to induce strong immune response, while epitopes with very strong binding (Kd stronger than 1 nmol/l) probably have reached a plateau for inducing response. Interestingly, under the current simulation, the suggestion of 1 µmol/l as a weak binding affinity agrees with other publications' cutoff for binders vs. no-binders for class II T-epitope prediction and shows biological relevance.[31,32]

Lastly, the dose of the antigenic protein is a recognized impacting factor for immunogenicity development.[33] In our simulation, a bell-shaped antigen dose–ADA response curve was obtained (,), suggesting that the model can reasonably capture the immune tolerance phenomenon. As demonstrated above, the model can potentially be applied to formulate or test hypotheses related to immunogenicity development, by modifying the model components or changing the parameter values.

In order to identify sensitive model parameters for ADA and Ag kinetics, sensitivity analysis was performed for the model with human parameters (). Overall, ADA kinetics is much more sensitive to model parameters than Ag kinetics (except that, as expected, Ag kinetics is very sensitive to antigen elimination rate, βAg). This is reasonable when we consider that the model is constructed for simulating ADA response, so its parameter values have more direct impact on ADA kinetics, while the impact on Ag kinetics is mostly secondary, through the formation of immune complexes. The sensitivity analysis indicates that the activation, proliferation, and differentiation of helper T cells and B cells are critical processes in the development of immunogenicity, and these parameter values have a big impact on the model outcome. It would be helpful to further confirm these parameter values by doing controlled in vitro or in vivo experiments in the future.

The current model represents a starting point and has potential for future development. For example, while plasma is assumed as the space for the immune response, ideally the site of action should also include lymphoid organs such as lymph nodes, spleen, and bone marrow, where various immune reactions such as antigen presentation, clonal selection, and long-lived plasma cell survival take place. DC activation is modeled rather empirically in the current settings. Upon the availability of more information to expand this module, the LPS activation signal could be replaced or complemented by other danger signals, e.g., tissue damage on injection. Some human-specific parameter values are directly taken from mouse values, and they may be properly scaled from mouse to human or be determined by performing appropriate experiments.

Several features of this model should be highlighted. First, the model recapitulates and integrates essential biological components and mechanisms for the development of immunogenicity and provides a wide spectrum of predictions for the kinetics of the immune system. Second, the model simulations agree with many known immunological phenomena, such as immune memory, antibody affinity maturation, and immune tolerance. Last but not least, the model was developed by incorporating both antigen- and host-specific characteristics, since immunogenicity ultimately depends on both the antigenic properties of the protein and the immunological environment of the host. Antigen-specific characteristics, such as T- and B-epitope content and affinity, and host-specific characteristics, such as MHC-II genotype and naïve T- and B-cell numbers, are designed as potential inputs for simulation. Since many antigen or host characteristics can be obtained from experimental approaches, the model can potentially be applied to quantitatively integrate results from immunogenicity assessment studies and provide simulations/predictions about the putative immune responses in human.

In summary, a novel multiscale, mechanistic model was constructed for the development of immunogenicity against antigenic/therapeutic proteins, by mathematically describing key underlying mechanisms and collecting physiologically plausible parameters from scientific publications. As a result, the current model simulations reasonably agree with many accepted immunological phenomena. Here, we used the model to formulate and test some hypotheses on influential factors for immunogenicity. We believe that the model can be explored to potentially serve as a supplemental tool to aid immunogenicity risk assessment, by quantitatively integrating protein- or patient-specific information from in vitro and in vivo experiments or from clinical studies.

Methods

The in vivo immune responses in human against a hypothetical antigenic protein were simulated. The protein, with molecular weight of 150 kD, and a half-life of 2 days, was injected i.v. at 50 mg/kg on day 0 and day 14. It is assumed that one promiscuous T-epitope is present on this protein, which binds to three of the six MHC-II molecules for one human subject, with a strong binding affinity of 150 nmol/l. Assumptions are also made for antigen-specific naïve T and B cells, with naïve T-cell frequency at 0.33 per million and B-cell frequency at 10 per million. Other parameters are kept the same as the generic human parameter set. The simulation was run for 100 days from day 0, and the kinetics for cell populations, antigen presentation, antigenic protein, ADA, and immune complex were recorded.

To explore the model behavior while changing the values of several factors of interest, the corresponding parameters were varied within a large range. The factors tested included naïve T-cell number (NT0), naïve B-cell number (NB0), antigenic protein dose (Ag0), and T-epitope–binding affinity against MHC-II (Kd). The naïve T-cell frequency ranges from 3.3 10–4 to 3.3 105 per million, and the naïve B-cell frequency ranges from 10–3 to 106 per million. The antigenic protein doses are from 3.3 × 101 to 3.3 × 1010 pmol, while the MHC-II–binding affinities (Kd = koff/kon) of T-epitope are from 10–3 to 106 nmol/l. The other model parameters were adapted from the model for human response against a hypothetical protein. The ADA profiles and cumulative ADA over 154 days (area under the curve for ADA-time profile) were simulated.

We performed a sensitivity analysis to evaluate the impact of individual parameter values on the interested state variables. The sensitivity calculation was performed for two key state variables (amount of therapeutic protein and amount of ADA) against the system parameters. The parameter set used for sensitivity analysis is the human set.

Control coefficients are used as indicators for sensitivity. The control coefficients of variable x to parameter p () were calculated using the following equation[34,35]:

The finite difference method was used to approximate the derivative in the above equation, where h is the step size:

The step size h was chosen as h = 1%·p, because 1% increase in the value of p has been shown to yield numerically stable results.[34] We reported CCmax, the value with the maximum absolute value over the time course of the system. A positive CCmax suggests that the parameter increase results in a increase in x value, while a negative CCmax means that the parameter increase causes a decrease in x value.

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

Table 1

Sensitivity analysis: top 10 sensitive parameters