Systems pharmacology models can be used to understand complex pharmacokinetic-pharmacodynamic behavior: an example using 5-lipoxygenase inhibitors.

Zileuton, a 5-lipoxygenase (5LO) inhibitor, displays complex pharmaokinetic (PK)-pharmacodynamic (PD) behavior. Available clinical data indicate a lack of dose-bronchodilatory response during initial treatment, with a dose response developing after ~1-2 weeks. We developed a quantitative systems pharmacology (QSP) model to understand the mechanism behind this phenomenon. The model described the release, maturation, and trafficking of eosinophils into the airways, leukotriene synthesis by the 5LO enzyme, leukotriene signaling and bronchodilation, and the PK of zileuton. The model provided a plausible explanation for the two-phase bronchodilatory effect of zileuton-the short-term bronchodilation was due to leukotriene inhibition and the long-term bronchodilation was due to inflammatory cell infiltration blockade. The model also indicated that the theoretical maximum bronchodilation of both 5LO inhibition and leukotriene receptor blockade is likely similar. QSP modeling provided interesting insights into the effects of leukotriene modulation.CPT: Pharmacometrics & Systems Pharmacology (2013) 2, e74; doi:10.1038/psp.2013.49; advance online publication 11 September 2013.

Arachidonic acid (AA) metabolites such as leukotrienes (LT) are key inflammatory mediators, whose overactivity causes bronchoconstriction in asthmatic patients. A number of therapies targeting this pathway are already marketed or in development–montelukast (Singulair; 10 mg once daily (q.d.)), a LT receptor (LTR) antagonist and zileuton (Zyflo IR/CR; 600 mg four times a day or 1,200 mg twice a day), a 5-lipoxygenase (5LO) redox enzyme inhibitor, are both marketed agents and 5LO activation protein inhibitors are in development.[1] In spite of extensive clinical experience with this pathway, some gaps in our understanding of this mechanism still exist. For example, the sparse clinical data[2,3,4] on the dose response of zileuton is interesting; on acute dosing, both 400 and 600 mg four times a day (q.i.d.) result in similar bronchodilatory response as measured by forced expiratory volume in 1 s (FEV1); however, when dosing is continued beyond 1–2 weeks, dose response emerges. Also, the relative bronchodilatory potential of different interventions along the AA pathway–e.g., LTR blockade vs. 5LO inhibition vs. 5LO activation protein inhibition is not known. Furthermore, it is not known whether zileuton doses higher than 600 mg q.i.d. would result in even higher efficacy. Understanding these properties would be critical to exploiting the full therapeutic potential of this pathway–for example, to develop a non-redox 5LO inhibitor with less frequent dosing regimen than zileuton (e.g., once daily). An empirical approach such as a pharmacokinetic-pharmacodynamic (PKPD) modeling approach is not suitable for answering these questions; a descriptive model would be unable to provide a mechanistic explanation of the observed PKPD data and a model developed based on 5LO inhibitor data would have limitations in predicting the effects of LTR blockade. Using a quantitative systems pharmacology (QSP) model, we have previously shown[5] that a non-redox 5LO inhibitor could have the same efficacy as a redox inhibitor. We further developed this model[6,7] using literature-published data to help answer these additional questions.

Several mathematical models describing various parts of the 5LO system have already been published.[8,9,10] These models describe basic aspects of the enzyme function–binditng AA in catalytic and regulatory sites of 5LO and its transformation to products–they do not take into account other important features such as LTA4 synthesis, reduction of 5-hydroperoxyeicosatetraenoic acid (HPETE) and other peroxides to activate 5LO (pseudoperoxidase reaction), reversible inactivation of 5LO by 5-hydroxyeicosatetraenoic acid,[11] and irreversible inactivation by LTA4.[12,13] Models of AA metabolism that include both lipoxygenase and cyclooxygenase pathways resulting in formation of LTs and prostanoids have also been reported.[12,13] The influence of various inhibitors of 5LO and cyclooxygenase on AA metabolism have been studied in these models. However, rate laws for enzymes involved in the AA metabolism have been described in semi-empirical manner using Michaelis–Menten type equations and, consequently, inhibition of 5LO by AA, and activation with product HPETE and pseudoperoxidase activity of 5LO have not been taken into account. Another effort to construct quantitative description of allergic airway inflammation has been performed by Walsh et al.[14] The computational model represents Boolean network including eosinophils, T cells, and other cell types as well as a variety of interleukins (IL), and is trained against literature data on C57BL/6 and BALB/c mice. However, this type of modeling does not allow us to understand mechanisms underlying complex PKPD behavior after administration of antiasthmatic drugs.

Therefore, a more detailed model of 5LO system which describes intracellular LT synthesis, extracellular LT transformation, and cell dynamics of eosinophils is required to fully understand the therapeutic potential of this pathway in asthma. This article describes the development and evaluation of such a model.

Results

In this manuscript, we define the term “5LO system” as a set of processes including intracellular biosynthesis of cysteinyl LT (CysLT) from AA by eosinophils, extracellular transformation of CysLTs, CysLT- and IL-5–mediated eosinophils maturation, migration and activation, production of IL-5 and histamine by eosinophils and mast cells, and effect of CysLTs and histamine on airway smooth muscles contraction. A schematic of the model is presented in . The model has been developed using data published in literature. These data and some of the key underlying assumptions are presented in . Briefly, the process described by the model is as follows: after maturation in the bone marrow in IL-5–dependent manner, eosinophils migrate to blood, where they are activated by CysLTs. Activated eosinophils produce CysLTs and IL-5, which further accelerates the process (positive feedback). Both activated and non-activated eosinophils produce histamine. Increase in IL-5 and CysLT blood concentration stimulates migration and accumulation of both activated and non-activated eosinophils to airways. Histamine and CysLT action in the airways results in bronchoconstriction. Zileuton inhibits intracellular production of CysLTs from AA; montelukast inhibits LT binding to its receptor to both activated and non-activated eosinophils and to cause bronchoconstriction.

The systems model describes the complex PKPD behavior of zileuton

To understand the complex PKPD relationship of zileuton observed in the clinic[2] and to understand the relative efficacy potential of 5LO inhibition vs. LTR inhibition with montelukast, we have simulated 400 and 600 mg q.i.d. zileuton doses and 10 (clinically marketed dose of montelukast as Singulair) and 50 mg montelukast doses q.d. – show the effect of these doses on FEV1, airway eosinophils and extracellular LT concentrations simulated with our model and clinically measured. These figures show that 5LO inhibition causes bronchodilation in two phases. In the first “acute” phase (0–1 days), bronchodilation is due to direct inhibition of airway LTD4 by ~60% (). Minimal changes in airway eosinophils are predicted in this phase (). Continued dosing at 400 mg results in equilibrium in the system between the airway eosinophils and the LTD4 and histamine they produce, which in turn attract eosinophils. At higher zileuton doses, the lower levels of LT and histamine achieved in the airways results in an unsustainable cell population in the airways, and hence a continued decrease in the cell numbers and subsequent bronchodilation. A final new steady state corresponding to higher FEV1 () is reached with low airway cells and LT/histamine concentrations at this higher dose after ~3 weeks of dosing.

Zileuton 600 mg q.i.d. has higher bronchodilation than montelukast 10 mg q.d.

Our simulations indicate that 5LO inhibition with zileuton at 600 mg q.i.d. has higher mean efficacy than LTR inhibition with montelukast at 10 and 50 mg q.d. The maximum FEV1 increase predicted by the model for zileuton 600 mg is 23% and that for montelukast 10 mg dosing is 5% (). Montelukast dosing 10 mg q.d. is predicted to result in essentially unchanged airway cell populations in peripheral tissues (see for comparison with clinically measured data) due to low blockade of LT action, similar to that reported by Reiss et al.[15] The marginal effect of montelukast on airway cell population (predicted by the model) was found to be due to high competition from airway LT for the CysLT receptors.

Doses of zileuton of 600 mg q.i.d. achieves maximum bronchodilation possible with this mechanism

We also ran simulations to evaluate whether higher doses of zileuton could result in higher efficacy. As shown in , 600 mg zileuton q.i.d. results in maximum efficacy. Higher doses result in faster approach to maximum efficacy, but no higher bronchodilation is predicted. Unlike zileuton, higher doses of montelukast are predicted to result in greater efficacy (), with doses >250 mg resulting in FEV1 increases similar to that of zileuton at the higher FEV1 steady state.

Discussion

We have illustrated how a systems pharmacology modeling approach can be used to understand the complex pharmacological behavior of the 5LO inhibitor zileuton and to compare potential efficacy of intervention at different points along the AA pathway–i.e., 5LO inhibition vs. LTR blockade.

The final model of the 5LO system is an ensemble model of six submodels each describing different aspects of the PKPD system. Each component submodel is a system of ordinary differential equations (ODE), which represent a simplified quantitative representation of the interactions within the subsystem. For example, the airway trafficking of eosinophils in our model is assumed to be entirely due to IL-5 and LTE4. However, other cytokines such as granulocyte-macrophage colony-stimulating factor and stem cell factor are also known to be involved in eosinophil maturation, release, and airway trafficking.[16] This is illustrated by the fact that complete inhibition of IL-5 activity reduced lung eosinophilia incompletely by mepolizumab.[17] Similarly, in our model LTs are produced only by eosinophils whereas there is evidence that basophils and mast cells may also be involved in LT production.[16] These important assumptions in our model are made to ensure that only the minimum number of variables is included in our system of equations in order to answer the main questions of interest. Therefore, while each component of the model is quantitative, the overall result should be considered semi-quantitative only. Within these constraints, the model predictions are roughly in line with observations on eosinophil changes observed in the clinic. For example, our model predicts complete reduction in airway eosinophil count at 600 mg q.i.d. dosing of zileuton, whereas Hasday et al.[18] reported ~70% reduction in bronchial lavage eosinophils count and preclinical data indicate a maximum reduction in airway eosinophils and CysLT of 80%.[19] For montelukast also, the predicted blood eosinophil reduction is in line with observations; 10% predicted vs. ~10% observed.[15,20]

Also, it should be noted that in our biophysical model, placebo effects on FEV1 changes are not accounted for; however, for the sake of qualitative description of bronchodilation by zileuton and to compare the relative effects of 5LO inhibition and LTR blockade, this level of abstraction is sufficient.

All submodel parameters were estimated with adequate precision (Supplementary Data J online). Each submodel provided adequate description of the underlying data. The model predicted the two-phase FEV1 change and the transition time from the first to the second phase after zileuton administration.

Our model provides a plausible explanation for the unusual PKPD behavior of zileuton. In the acute-dosing phase, bronchodilation is achieved by inhibiting LT action on the airways. However, airway eosinophils continue to produce LTs, which continue to attract and activate eosinophils in the airways. Due to the positive feedback in this system, higher level of LT blockade is required to disrupt the trafficking process and therefore, higher doses of zileuton–i.e., doses ≥600 mg–are required to achieve additional bronchodilation. The time required for this additional bronchodilation to be achieved corresponds to the lifespan of activated eosinophils in the airways (3–4 days). Intriguingly, the model suggests that this phenomenon is also true for montelukast, albeit at doses >250 mg q.d., which have not been reported as tested in clinical trials. For montelukast, direct blockade of LTR provides acute bronchodilation; however, competition from high levels of airway LT limits its ability to block eosinophil trafficking, except at doses >250 mg. To the best of the authors' knowledge, this is the first instance of this mechanism being shown to be a plausible explanation of the observed data. Apart from 5LO inhibitors and LTR antagonists, this mechanism is also likely applicable for other anti-inflammatory bronchodilation agents including steroids, which have delayed action.

The model predicts that presently marketed doses of zileuton provide maximum possible bronchodilation available for this mechanism, but those of montelukast only provide bronchodilation due to direct inhibition of LT action. LTR antagonists with greatly improved receptor affinities may provide improved asthma treatment over montelukast. Our simulations also suggest that at marketed doses, zileuton has higher FEV1 increase (23%) than montelukast 10 mg (~5%). The predicted effect for montelukast is broadly in line with a meta-analysis,[21] which indicated bronchodilation in the 5–10% range for montelukast doses of 10–50 mg q.d. Reported bronchodilation for zileuton is lower than predicted[2,3,4] (15-20%), presumably due to the assumptions regarding the role of eosinophils as the sole inflammatory cells in the airways. Furthermore, recently Kubavat et al.[22] have shown in a head-to-head trial that zileuton 1,200 mg twice daily is a more effective bronchodilator (measured by peak expiratory flow rate), than montelukast 10 mg, supporting the model prediction. Therefore, even though the model predictions are semi-quantitative, in the absence of a large database of clinical data, this analysis was able to provide an initial comparison of these two targets' potential efficacy.

The AA system has been extensively studied; hence there are sufficient data in literature to help develop this model. However, for many novel targets, there is often limited information on the link among target modulation, biomarker changes, and ultimate clinical effect. Indeed, this has often been cited as a key reason for not investing in QSP model-based approaches early during drug research and development. Our model suggests that in such situations where mechanistic information is sparse, investment in in vitro, in vivo, or human tissue assays to generate mechanistic data can provide benefits by enabling such a model-based analysis. This approach is preferable to the often-existing philosophy of running experiments to “confirm” the efficacy of the target in poorly understood animal “disease models”. Development of such mechanistic information to enrich our knowledge of the pathways of interest should be an important goal of translational research and systems models such as these can provide a quantitative framework for making sense of such complex data.

In conclusion, we have developed a QSP model to help understand the complex behavior of zileuton a 5LO inhibitor. The model provides a number of insights into the AA–LT mechanism causing bronchodilation–it is able to explain the complex PKPD behavior of zileuton; it suggests that zileuton is more efficacious than montelukast at currently marketed doses; it also indicates that improvements in target affinity are more likely to be useful for LTR antagonists than for 5LO inhibition. Our results illustrate how QSP techniques can form a key component of quantitative drug discovery and development strategy.

Methods

The QSP model of the 5LO system is an ensemble model consisting of six components, describing (i) intracellular LT biosynthesis, (ii) extracellular conversion of LT, (iii) distribution model describing transport of molecules between compartments, (iv) cell dynamics model, (v) biophysical model describing the link between LT levels and bronchial tone, and (vi) PK model. The main stages of model construction are outlined in . The final system of ODEs, experimental facts and assumptions underlying the model structure, rate laws of each process, values of model parameters, and experimental/clinical data used to identify them are presented in Supplementary Materials online. The final system consisted of 33 ODEs, 64 rate laws, and 113 parameters, of which 26 were estimated using literature data and the rest directly taken from literature-reported estimates. Below, we will briefly describe each of the component submodel (i–vi). Submodels (i–iv) are shown in Supplementary Materials (section A, Figure A1) online; submodels (v) and (vi) are shown in the Supplementary Materials (sections F and G) online.

Submodel of intracellular LT biosynthesis. This submodel represents ODE system describing intracellular LTC4 biosynthesis from AA. The model includes reactions catalyzed by 5LO, cytosolic phospholipases A2, glutathione peroxidase, 5-hydroxyeicosanoid dehydrogenase, and processes of LTA4 degradation and LTC4 excretion from eosinophils located in blood plasma (processes 1–10) and airways interstitium (processes 31–40). Rate laws of the processes and parameter values of the submodel have been taken from Karelina et al.[5] (and references contained therein). Intracellular concentration of LTC4 in eosinophils of healthy subjects was obtained from Bandeira-Melo et al.[23] The ODE system describing intracellular LT biosynthesis, rate equations, and values of parameters are presented in Supplementary Materials (sections A and B) online.

Submodel of extracellular conversion of LT. This submodel represents ODE system describing conversion of LTC4 to LTD4 and further to LTE4 catalyzed by γ-glutamyl transpeptidase and dipeptidase located in blood plasma (see processes 11 and 12 in Supplementary Materials (section A, Figure A1) online) and airway interstitium (processes 41 and 42), correspondently. All three CysLTs are exposed to degradation in blood plasma (processes 13–15) and airway interstitium (processes 55–57). Parameters of the model have been identified on the basis of in vitro data,[24] ex vivo experimental data measured for blood of healthy subjects,[25] and in vivo experimental data (time series of radiolabeled LTs after intravenous infusion) measured for monkey.[26] The ODE system describing extracellular conversion of LT, rate equations, and values of parameters are presented in Supplementary Materials (sections A and C) online.

Submodel of distribution. This submodel describes transport of LTC4, LTD4, LTE4, and histamine between blood and airway (processes 43, 44, 45, 52), and transport of IL-5 between airway and blood (process 54) and blood and bone marrow (process 28). The submodel describing distribution of the molecules is based on physiologically based PK approach,[27] on the one hand, and is simplified as much as possible to meet the compartmental structure of QSP model of 5LO system. All parameters of the submodel can be divided into three groups. First group consists of parameters and those values are calculated on the basis of physicochemical properties of the molecule (for example, tissue partition coefficients, logP, pKa, etc). Values of the parameters have been taken from literature.[28,29] Second group includes parameters characterizing physiological properties of the living system (for example, hematocrit, lymph flow rates, blood flow rates, etc). Values of the parameters have been taken from literature.[30] Third group includes only two parameters: effective permeability constants of blood to airway and blood to bone marrow transport of IL-5. The values of the parameters have been identified via fitting of the model against in vivo data measured in asthmatic patients.[31] The ODE system describing distribution of LTC4, LTD4, LTE4, histamine, and IL-5, rate equations, and values of parameters are presented in Supplementary Materials (sections A and D) online.

Submodel of cell dynamics. This submodel represents ODE system describing maturation, migration, death, and activation of the eosinophils with LTD4 and LTC4 as well as processes associated with the cells such as production/degradation of IL-5 and histamine (see processes 16–27, 29 and 30, 46–51, 53, 58 and 59 in Supplementary Materials (section A, Figure A1) online). Parameters of the model have been identified on the basis of the ex vivo and in vivo experimental data measured for asthmatic patients and healthy subjects.[18,31,32,33,34,35] ODE system describing cell dynamics of eosinophils, rate equations, and values of parameters are presented in Supplementary Materials (sections A and E) online.

Submodel of the link among LT, histamine levels, and bronchial tone. This submodel allows us to couple LTD4, LTC4, and histamine concentrations produced by eosinophils with clinically measured endpoint such as FEV1. The submodel represents several algebraic expressions (Supplementary Materials (section F) online) empirically describing level of bronchial smooth muscle contraction as a function of concentration of bronchoconstrictors such as LTD4, LTC4, and histamine. Parameters of the model have been identified on the basis of the ex vivo experimental data measured for guinea pigs and healthy subjects.[36,37,38]

The PK submodel. PK submodel parameters for zileuton and montelukast were estimated using data from refs. [39,40,41]. The plasma PK of zileuton was described by a simple two-compartment model. The plasma PK of montelukast was described by an one-compartment model. It has been assumed that intracellular concentration of zileuton is equal to that of free unbound concentration in plasma or airway interstitium. It has been also assumed that free unbound concentration of montelukast in airway interstitium is equal to that in plasma. The ODE system describing PK of zileuton and montelukast, rate equations, and values of parameters are presented in Supplementary Materials (sections A and G) online.

All parameters were estimated using Hooke–Jeeves method implemented in the DBSolve Optimum package (Institute for Systems Biology SPb, Moscow, Russia).[42] Parameter identification was performed individually for each submodel by fitting to literature data sets pertinent to the part of the system described by the submodel. The 95% confidential intervals were calculated for fitted parameters using method described by Motulsky[43] (see Supplementary Materials online).

Conflict of interest

The authors declared no conflict of interest.

Author contributions

B.M.A., O.D., and D.F. wrote the manuscript. B.M.A., D.F., and P.H. van der G. designed the research. O.D., T.K., D.S., E.M., G.S., and O.D., Jr. performed the research. O.D., D.S., E.M., G.S., and O.D., Jr. analyzed data. O.D. contributed new reagents/analytical tools.

Study Highlights

Table 1

Description of available experimental data and facts

Table 2

Description of stages of model construction