A multiscale model of interleukin-6-mediated immune regulation in Crohn's disease and its application in drug discovery and development.

In this study, we have developed a multiscale systems model of interleukin (IL)-6-mediated immune regulation in Crohn's disease, by integrating intracellular signaling with organ-level dynamics of pharmacological markers underlying the disease. This model was linked to a general pharmacokinetic model for therapeutic monoclonal antibodies and used to comparatively study various biotherapeutic strategies targeting IL-6-mediated signaling in Crohn's disease. Our work illustrates techniques to develop mechanistic models of disease biology to study drug-system interaction. Despite a sparse training data set, predictions of the model were qualitatively validated by clinical biomarker data from a pilot trial with tocilizumab. Model-based analysis suggests that strategies targeting IL-6, IL-6Rα, or the IL-6/sIL-6Rα complex are less effective at suppressing pharmacological markers of Crohn's than dual targeting the IL-6/sIL-6Rα complex in addition to IL-6 or IL-6Rα. The potential value of multiscale system pharmacology modeling in drug discovery and development is also discussed.CPT: Pharmacometrics & Systems Pharmacology (2014) 3, e89; doi:10.1038/psp.2013.64; advance online publication 8 January 2014.

Introduction

Among the factors impeding discovery of new drugs is a lack of understanding of the complex biological systems underlying diseases and systemic repercussions of drug candidate–target interactions.[1] Model-based drug discovery and development has been increasingly utilized to address this challenge[2]; in recent years, the multidisciplinary approach of systems pharmacology that combines systems biology with pharmacokinetic/pharmacodynamic (PK/PD) modeling principles has progressively drawn attention in the pharmaceutical industry.[3] Systems pharmacology holds promise as a discipline that can help develop holistic understanding of interactions between drugs and complex biological systems underlying disease and has been termed the “next iteration” of model-based drug discovery and development.[4,5,6] In this study, we have developed a multiscale systems model of Crohn's disease to investigate candidate therapeutic strategies and to demonstrate the potential of quantitative systems pharmacology in drug discovery and development.

Crohn's disease is an inflammatory bowel disease that is thought to be caused by a combination of genetic and environmental factors,[7] resulting in a proinflammatory environment in the mucosal layer of the gastrointestinal (GI) tract.[8] Along with several other cytokines that are perturbed in Crohn's disease, interleukin (IL)-6 has been shown to be increased in intestinal mucosa and in peripheral blood.[9,10,11,12] IL-6 signaling is thought to be an important player in Crohn's disease contributing to enhanced T-cell survival and apoptosis resistance in the lamina propria along with elevated chemokine secretion.[13]

IL-6 signaling can occur via membrane-bound IL-6 receptor (IL-6Rα)–mediated classical pathway as well as soluble IL-6 receptor (sIL-6Rα)–mediated trans-signaling pathway (). The membrane-bound IL-6/IL-6Rα dimer or the circulating preformed complex of IL-6/sIL-6Rα can recruit a membrane-bound gp130 coreceptor, forming a heterotrimer which subsequently dimerizes, resulting in an active hexameric IL-6R receptor complex.[14,15,16] This complex initiates phosphorylation of gp130-bound Janus kinase (Jak) family proteins and subsequent signal transducer and activator of transcription 3 (STAT3) phosphorylation.[17] Trans-signaling may be important in disease because transmembrane IL-6Rα is expressed only in a small number of cell types, e.g., hepatocytes and a subset of leukocytes, whereas gp130 is almost ubiquitously expressed.[18] Alternative gene splicing and shedding of transmembrane IL-6Rα are thought to be responsible for production of sIL-6Rα.[19] IL-6Rα shedding was suggested to be accelerated by proteins including C-reactive protein (CRP).[20] The coreceptor gp130 can also exist in a soluble form (sgp130), acting as a natural inhibitor of trans-signaling by sequestering the IL-6/sIL-6Rα complex[16] (). The downstream effects of IL-6 signaling are multifarious, ranging from differentiation of lymphocytes to expression of acute-phase proteins such as CRP.[18] The reported increase in serum concentrations of CRP,[21] sIL-6Rα,[22] and pSTAT3 in colon biopsies[23] of Crohn's disease patients may be driven by increased IL-6, acting through the complex network of IL-6 signaling ().

A pilot clinical trial with anti–IL-6R antibody, tocilizumab, has shown positive clinical activity in Crohn's disease,[24] suggesting inhibition of IL-6 signaling as a therapeutic strategy. Several other strategies that modulate IL-6 signaling have been proposed, which include targeting IL-6[25] or IL-6R,[24] using sgp130 as a natural antagonist of IL-6 trans-signaling,[15,16] or simultaneously inhibiting classical (target IL-6 or IL-6R) and trans (target IL-6/sIL-6R complex) IL-6 signaling.[15] Antibodies against these targets in Crohn's disease are in early clinical development; however, clinical biomarker or efficacy data are currently unavailable.

Systems biology models of IL-6–mediated immune signaling have been published previously,[26,27,28,29] but these are single-scale models at the cellular level, and their translation to clinical end points has not been established. At a very different level of detail, PK/PD relationships of Ab–receptor interaction and biomarker or efficacy modulations are characterized in vitro and in vivo in preclinical and clinical studies. However, the PK/PD understanding is often limited by models based on limited data and lack of incorporation of biology. With various new IL-6–related therapeutic agents in the pipeline, we have developed a multiscale system model that integrates current knowledge about the biology of IL-6–mediated immune response in Crohn's and used it for mechanistic assessment and comparison of several proposed therapeutic strategies. We have further used the model to explore the value of quantitative system pharmacology modeling in discovery and development, particularly in the areas of target selection, candidate optimization, and dose selection.

Results

Construction of the multiscale model of IL-6 signaling in Crohn's disease

The multiscale model comprises three overlapping structural modules spanning spatial scales from cellular to organ levels. The first module describes events at the cellular level and consists of a reduced model of IL-6–mediated signal transduction. The second module is made up of target organs relevant in Crohn's disease, and the first module (the cell signaling model) is embedded within these organs. The third and final module is a general PK model for monoclonal antibodies, including specific target binding that is linked to the first two modules to study the action of drug on the system. Each of these modules is described below.

Module 1: IL-6–mediated cellular signal transduction. Detailed models of IL–mediated Janus kinase-Signal Transducer and Activator of Transcription (Jak-STAT) signaling, including IL-6/STAT3 signaling, have been described previously.[26,30,31,32] Starting from the canonical structure of these models (schematic diagram in ), we reduced the signaling network using a rationale-based approach to obtain a more parsimonious description while retaining the essential features of the IL-6 signaling pathway (). Simplifying assumptions used to reduce the model are listed in Supplementary Text S1.

Module 2: integrating cell signaling with organ-level kinetics. Instead of constructing a large-scale, physiologically based model accounting for all major organs in the body, we chose a parsimonious description by identifying the two most important target organs in Crohn's disease and modeling them as separate, communicating compartments (). The first of these compartments represents the GI tract, which is the organ affected directly by Crohn's disease.[8,13] The second compartment models the liver, which was included for its prominent role in the upregulation of acute-phase proteins associated with Crohn's disease, such as CRP.[27] We further assumed that each of these organs was separately made up of homogeneous cell populations with identical signaling dynamics with minor differences (see Supplementary Text S1 and S2). Assuming that the IL-6 signaling network topology and reaction rate parameters remain conserved across different organs, we embedded the IL-6 signaling module into the cells of these organs. The organs were allowed to exchange soluble components through blood serum, which was modeled as a separate compartment representing the circulation. Briefly, starting from a single cell model, we scaled up to organs, modeling them as compartmentalized collections of cells and connected the organs through blood circulation (). Assumptions accompanying this module are listed in Supplementary Text S1.

Module 3: PK of monoclonal antibodies. To study the effects of different therapeutic agents, we developed a module to incorporate drug PK. The antibody PK model comprises target-binding, dynamic exchange of free and bound drug between organ compartments, and clearance of free as well as bound antibody (Supplementary Figure S1). This model of drug PK describes the equilibration between multiple compartments incorporating both linear and target-mediated disposition of the antibody.

Parameter estimation

We exploited the modular structure of our model to break down the parameter estimation problem into a two-step process. In the first step, the model was divided into structural modules, and parameters were either assigned values available from the literature or estimated for each module in isolation using any published data we could find (example in ). The model was reassembled by putting the modules together, and the parameters were assigned values obtained from the first step as initial estimates. It was expected that the estimates obtained from isolated modules would no longer be accurate with complex interactions between the modules. Accordingly, in the second step, the initial parameter estimates were further fine-tuned to fit the system-level behavior of the whole model to available data.

When optimizing the parameters in step 2, we hypothesized that observed differences in concentrations of molecular markers between healthy subjects and Crohn's disease patients are a systemic outcome of the differences in IL-6 concentrations. Therefore, although the model does not account for the mechanism of increased IL-6 in Crohn's patients, it explains differences in other molecular markers as a consequence of IL-6 signaling. Given the large variability in measurements of biomarkers in human subjects and the sparsity of data, we chose parameter values that produced reasonable fits across all data points rather than parameter sets that resulted in over fitting of some data points at the cost of others (Supplementary Text S3). Using baseline concentrations of a set of biomarkers available from healthy subjects and Crohn's disease patients, the model was optimized to fit its steady-state values to measured typical values, and the ratio of fitted value to measured value was close to 1 for all molecular markers ( and ).

Since our study was largely focused on assessing the effects of treatment with monoclonal antibodies, we calibrated the drug PK parameters to fit the serum concentration of drug to median PK of monoclonal antibodies as reported by Dirks and Meibohm[28] (,).

Model validation

Data from a previously published pilot trial of a humanized anti–IL-6Rα monoclonal antibody, tocilizumab, in Crohn's subjects[24] was used to validate the model. Using the optimized model, we simulated the pharmacological effects of treatment with placebo or 8 mg/kg anti–IL-6Rα antibody at intervals of 4 or 2 weeks (). The simulated anti–IL-6Rα antibody was assigned a Kd of 1 nmol/l as an order of magnitude estimate based on target-binding studies on tocilizumab.[29] The PK profile of anti–IL-6Rα antibody in our simulation was in agreement with the reported PK profile of tocilizumab[24] (Supplementary Figure S2). Comparing the simulation results with the mean CRP levels measured by Ito et al. under identical dose, we found that the model was successful at capturing the trend shown by the clinical study even though it could not reproduce the data quantitatively (). Similar to the clinical data, the model showed that an every 4-week intravenous (i.v.) administration of tocilizumab transiently reduced serum CRP, only to be restored to baseline level by the end of each 4-week cycle. For an every 2-week i.v. dosing schedule, both the model and clinical data showed a sustained reduction in serum CRP levels. However, under both dose schedules, the model underestimated the reduction in serum CRP when compared with the clinical data (compare and ).

Dose–response characteristics of anti–IL-6 and anti–IL-6R antibodies are distinct

Using the optimized model, we simulated the effects of anti–IL-6 and anti–IL-6Rα on PD biomarkers including serum CRP, pSTAT3 activity in the GI tract, and IL-6 levels in the GI tract and serum after every 4-week i.v. dosing for a 12-week period. To examine the role of drug dose and antibody target–binding affinity in determining drug efficacy, we ran multiple simulations while systematically varying the dose and binding affinities of the antibodies. When administered at sufficiently high doses, both anti–IL-6 and anti–IL-6Rα (assumed to target both transmembrane and soluble IL-6Rα) treatments lead to the suppression of serum CRP and gut pSTAT3 in a time-dependent manner (,, respectively). Measured in terms of CRP suppression at the end of the 12-week treatment, anti–IL-6 antibody had a graded dose response tending toward saturation at high dose (). Anti–IL-6Rα antibody, on the other hand, remained ineffective over a large dose range but suddenly transitioned to high efficacy at a certain dose threshold (). We found that this qualitative difference of graded vs. switch-like response was not CRP specific and was maintained across multiple biomarkers (Supplementary Figure S3). Furthermore, both antibodies showed increased effectiveness with increasing target-binding affinity.

To explain this difference in response to the two antibodies, we examined their PK and found that despite identical parameters of distribution and linear clearance, the two types of antibodies could have very different PK time profiles depending on dose. For example, at 200 mg dose administered every 4 weeks, the anti–IL-6 antibody concentration was sustained at high enough levels to suppress IL-6 signaling (), whereas the anti–IL-6Rα antibody was severely depleted in 2 weeks (). Increasing the dose to 500 mg sustained anti–IL-6Rα over the dose cycle in a manner similar to anti–IL-6, resulting in higher efficacy (,). Increasing the dosing frequency also resulted in greater efficacy at lower doses of anti–IL-6Rα (Supplementary Figure S4).

Targeting the IL-6/sIL-6R complex in addition to free IL-6 or IL-6R as a therapeutic approach

In all the simulations presented so far, it was assumed that the anti–IL-6 and anti–IL-6Rα antibodies bind only to free IL-6 and IL-6Rα, respectively. Given the hypothesized importance of the role of IL-6/sIL-6Rα complex in Crohn's disease,[15] we ran simulations to assess the effects of targeting this complex in addition to free IL-6 or IL-6Rα.

When anti–IL-6 was allowed to bind to IL-6 and IL-6/sIL-6Rα complex with the same affinity, the dose–response characteristics of the antibody showed improved pharmacological efficacy over a range of doses and binding affinities in a dose-dependent manner (). However, similar improvement was not observed in the case of anti–IL-6Rα treatment, and the dose response remained largely unchanged ().

Evaluating alternative therapeutic approaches using the model

We performed global sensitivity analysis on the optimized model to assess the sensitivity of steady-state values of a panel of biomarkers to parameters of the model (see Methods section). The parameters to which the output was found to be least and most sensitive are shown in (see Supplementary Figure S5 for a complete list). The parameters that most strongly influenced the output biomarkers were rate of synthesis of IL-6 in the GI tract, rate of degradation of soluble receptor, parameters associated with IL-6–mediated production of CRP and CRP degradation. The production and degradation rates of sgp130 were among the least influential parameters. The sensitivity analysis suggests that while IL-6 and sIL-6Rα are among the best targets, sgp130-based therapies may not be as effective.

A fusion protein combining the extracellular portion of gp130 with the Fc region of human IgG1 (sgp130Fc) has been shown to inhibit IL-6 trans-signaling in cultured cells.[16] To study the potential of sgp130 as a therapeutic approach in Crohn's disease, we simulated i.v. dosing of sgp130Fc at 4-week intervals. We assumed that sgp130Fc binds the IL-6/sIL-6Rα complex with the same affinity as natural sgp130 and that the PK of sgp130Fc is identical to that of a typical monoclonal antibody with linear clearance. Dose–response curves in terms of CRP suppression for several dose schedules are shown in . Similar dose–response curves were obtained for pSTAT3 suppression in gut (Supplementary Figure S6). In agreement with the sensitivity analysis which suggested relative insensitivity of the output to sgp130, simulations showed sgp130Fc to be effective only at very high and frequent doses.

Discussion

We have constructed a system-level, multiscale, deterministic model of IL-6 signaling in Crohn's disease that elucidates how ideas from traditional pharmacometrics can be combined with strategies from systems biology to design disease models in systems pharmacology. Our work advances existing systems pharmacology modeling methodologies and illustrates model development strategies with particular focus on integrating the role of cell signaling networks in disease. Although highly detailed quantitative models of system dynamics are desirable for several reasons,[33] the construction of such models is often hindered by the sparsity of pertinent molecular-level kinetic measurements. Therefore, a balance needs to be struck between the level of detail in the model against the quality and quantity of data available to reliably incorporate such detail. To attain this balance, whenever possible, we have chosen to favor parsimony in describing the biological processes in our model. Notably, despite the simplifying assumptions and limited amount of data available to train the model, we were able to qualitatively validate the model against an independent clinical data set (). Particularly, even though the CRP time course was not exactly replicated quantitatively by the model in the validation step, the model was successful at qualitatively predicting the long-term behavior of CRP following drug treatment in the typical Crohn's disease patient showing the model's potential as a tool to study the effects of perturbation through various treatments.

We compared four different strategies that have been proposed to treat Crohn's disease by inhibition of IL-6 signaling through targeting IL-6, IL-6Rα, IL-6/sIL-6Rα in addition to IL-6 or IL-6Rα, or using sgp130Fc as an antagonist of IL-6 trans-signaling. The model predicts that two of these strategies, targeting IL-6 or IL-6Rα, can effectively suppress markers of inflammation (CRP) and IL-6 signaling in target organs. Using sgp130Fc to exploit the natural ability of sgp130 to inhibit IL-6 signaling was predicted to be effective only at very high and frequent doses. Simulated dose–response curves as well as sensitivity analysis done on the model suggest that targeting IL-6 requires lower and less frequent dosing as compared with anti–IL-6Rα to achieve similar reduction of the systemic effects of IL-6 signaling. This statement must be qualified, however, by our observation that the model underestimated the effect of anti–IL-6Rα treatment when compared with clinical data (). Nonetheless, the dose–response characteristics of anti–IL-6Rα are clearly distinct from those of anti–IL-6 (,). This difference arises due to the immense disparity between the levels of free IL-6 (order of 10 pg/ml ≈ 10−4 nmol/l)[12,13] and free sIL-6Rα (order of 100 ng/ml ≈ 1 nmol/l) in circulation.[22] At doses below a certain threshold, the antibody becomes the limiting factor in the anti–IL-6Rα/sIL-6Rα interaction as most of the antibody is sequestered by sIL-6Rα (, 200 mg curve). Once this threshold is crossed, the antibody overwhelms sIL-6Rα production and accumulates in the system in sufficient amount (, 500 mg curve) to suppress IL-6 signaling, resulting in the switch-like dose–response curve seen in .

The PK profiles of anti–IL-6 and anti–IL-6Rα antibodies are different in our model (compare and ) due to differences in the kinetics of their interactions with their targets. Nonlinear clearance very similar to that predicted by the model for the anti–IL-6Rα antibody has been observed for tocilizumab at therapeutic doses in clinical studies in rheumatoid arthritis.[34,35] This provides further validation to our modeling approach and shows that a well-constructed model of the underlying biological system can automatically accommodate target-mediated disposition of the antibody, reducing the need for case-by-case optimization of drug PK.

The IL-6/sIL-6Rα complex is an important intermediate in IL-6 signaling and has been proposed as a therapeutic target.[15] We tested the relative gain in the efficacy of biomarker modulation from combined targeting of the complex in addition to IL-6 or IL-6Rα. Greater dose-dependent modulation of the immunological biomarkers was suggested in the case of anti–IL-6 but not for anti–IL-6Rα. As argued above, at small drug doses, substantial concentration of sIL-6Rα remains in circulation and any IL-6/sIL-6Rα complex blocked by the antibody is quickly replenished, negating the effect of targeting the complex. At larger drug doses when sIL-6Rα is effectively suppressed by the antibody, IL-6/sIL-6Rα complex is naturally depleted because reduced sIL-6Rα concentration drives the IL-6/sIL-6Rα binding equilibrium toward the unbound state. Taken together, at both high and low concentrations, the effect of targeting the complex does not add much to the efficacy of IL-6Rα targeting.

Targeting IL-6/sIL-6Rα alone to inhibit trans-signaling by sgp130Fc was found to be effective only at very high and frequent doses (). This result was supported by sensitivity analysis using the model, which indicated that altering the levels of sgp130 has very little effect on steady-state values of a set of output parameters. This is explained by the high baseline level of sgp130 (≈ 300 ng/ml or 3 nmol/l) as compared with free IL-6 (≈ 10 pg/ml or 4 × 10−4 nmol/l).[10,11,12] This means that sgp130 is basally present in large excess when compared with the IL-6/sIL-6Rα complex, which is limited by the low IL-6 concentration. The implication of this concentration difference is that small perturbations in sgp130 level have minimal effect on system dynamics, and any observable effects require large changes in sgp130.

One of the limitations of the current model is that it underestimates the CRP response under anti–IL-6R treatment. The reasons for this are twofold. First, there was limited data available to both train and validate the model, meaning that estimates of some parameters could be inaccurate. Although it may ultimately be possible to simulate first use of a drug in humans in silico, achieving this routinely and with confidence requires more two-way flow of information from clinical trial back to the model for model training, and from model-generated hypothesis to the clinic for testing. Currently, we do not have sufficient dose–response clinical data available in the literature to warrant model refinement based on the data. The model will be refined to incorporate emerging data from ongoing anti-IL-6–related clinical studies to improve its predictive power. Second, the structure of our model is limited by our biological knowledge about the system, which leaves room for improvement. For example, while we have singled out IL-6 signaling to study the Crohn's disease system, other cytokines such as IL-10 and tumor necrosis factor-α are known to be important in this disease. Future models could also include distinct subpopulations of cells, such as T-lymphocyte subsets in the GI tract which differentially influence the disease state.

While admitting scope for improvement, the present work demonstrates the potential of the systems pharmacology approach in drug discovery and development. As demonstrated by model-based comparison of the competing IL-6 signaling–related therapeutic strategies, the system model can be tremendously valuable at the drug discovery stage. Whereas intuition, simplified experimental systems, and single-scale computational models are limited in their scope, a multiscale systems pharmacology model is better able to predict outcomes arising from system-level interactions and combination treatment, and has enormous potential in facilitating target selection and candidate optimization in drug discovery. The applications of systems pharmacology are not limited to the discovery stage. Phase I studies in patients are particularly relevant to the development of a quantitative systems pharmacology model, because it is at this stage that the human PK and PD information of the new therapeutic agent becomes available. Once a system model structure is in place, integration of phase I biomarker data can facilitate model refinement, resulting in more quantitative prediction of clinical efficacy and/or safety responses based on correlation with biomarker response. The refined model can aid optimal proof-of-concept phase II design and support an optimal dosing regimen. It is at this stage that one can truly test the correct selection of targets, drugs, dosing, and therapeutic strategy. Model development during phase II may be used to address different questions related to the understanding of drug effect for both safety and efficacy. Overall continuous integration of system-level models with data from multiple sources and various development stages would permit more complete exploration and interpretation of clinical response data for new and existing therapeutic agents and could enhance decision making throughout the development process.

Methods

Model implementation and parameter optimization. All modeling, simulation and analysis was done in Matlab (R2012b). The SimBiology toolbox of Matlab was used to implement the model (see Supplementary file 1 CD model and example Matlab script to run the model in Supplementary file 2 Sample Matlab code). To facilitate platform-independent exchange, Systems Biology Markup Language (SBML) versions of the models are also provided (Supplementary Files 3–6). Enzymatic reactions in the signaling pathway and pSTAT3-induced protein expression were modeled using Michaelis–Menten kinetics. All other events were modeled using mass action kinetics. All first-order rates have units of 1/hr, second-order rates are in l/(nmol·hr), and all concentrations are defined in nmol/l. All simulations were performed by first running the model to steady state and then simulating the effects of exogenous perturbations. Parameter estimation was performed using a combination of fitting by hand and functions from the optimization toolbox of Matlab. Further details of the fitting process are described in Supplementary Text S3.

Global sensitivity analysis. Partial rank correlation coefficients between the steady-state output and input parameters were calculated (Supplementary Text S4).[36] To get a single sensitivity measure per parameter, the partial rank correlation coefficientss for each parameter were summed over the output species and divided by a constant to scale the highest sensitivity to 1 (see Supplementary Figure S5).

Author Contributions

C.L. and G.D. wrote the manuscript. C.L. and G.D. designed the research. C.L., G.D., L.F., M.H., J.H., and A.H. performed the research. G.D. analyzed the data.

Conflict of Interest

The authors declared no conflict of interest.

Study Highlights

Table 1

Concentrations of biomarkers used to optimize model parameters