Visual Overview
Abstract
Crohn's disease (CD) is a complex inflammatory bowel disease whose pathogenesis appears to involve several immunologic defects causing functional impairment of the gut. Its complexity and the reported loss of effectiveness over time of standard of care together with the increase in its worldwide incidence require the application of techniques aiming to find new therapeutic strategies. Currently, systems pharmacology modeling has been gaining importance as it integrates the available knowledge of the system into a single computational model. In this work, the following workflow for robust application of systems pharmacology modeling was followed: 1) scope definition; 2) species selection and circulating plasma levels based on a search in the literature; 3) representation of model topology and parametrization of the interactions, after literature data extraction and curation, and the implementation of ordinary differential equations in SimBiology (MATLAB version R2018b); and 4) model curation and evaluation by visual comparison of simulated interleukin (IL) concentrations with the reported levels in plasma, and sensitivity analysis performed to confirm model robustness and identify the most influential parameters. Finally, 5) exposure to two dose levels of recombinant human IL10 was evaluated by simulation and comparison with reported clinical study results. In summary, we present a quantitative systems pharmacology model for the main ILs involved in CD developed using a standardized methodology and supported by a comprehensive repository summarizing the most relevant literature in the field. However, it has to be taken into account that external validation is still pending as available clinical data were primarily used for model training.
SIGNIFICANCE STATEMENT Crohn's disease (CD) is a complex heterogeneous inflammatory bowel disorder. Systems pharmacology modeling offers a great opportunity for integration of the available knowledge on the disease using a computational framework. As a result of this work, a comprehensive repository along with a quantitative systems pharmacology model for the main interleukins involved in CD is provided. This model is useful for the in silico evaluation of biomarkers and potential therapeutic targets and can be adapted to address research gaps regarding CD.
Introduction
Crohn's disease (CD) is a complex inflammatory bowel disease (IBD) with a controversial autoimmune nature which causes a functional impairment of the gut wall leading to abdominal pain, severe diarrhea, fatigue, weight loss, and malnutrition (Kaser et al., 2010; Nóbrega et al., 2016). The disease is characterized by discontinuous ulceration and inflammation throughout the entire gastrointestinal tract, with the ileum and colon being the most affected areas (Doherty et al., 2018).
CD is a heterogeneous disorder with a multifactorial etiology in which both genetics and environment play a key role (Gajendran et al., 2018). The pathogenesis of CD appears to involve a primary defect in innate immune mechanisms leading to an aggressive acquired (T-cell) immune response to a subset of commensal enteric bacteria, resulting in dysregulated expressions of molecules involved in proinflammatory and anti-inflammatory processes (Sartor, 2006; Park et al., 2017). Some of these upregulated molecules, named interleukins (ILs), are tumor necrosis factor alpha (TNFα), IL6, IL1β, IL17, IL21, IL12, IL23, interferon gamma (IFNγ), transforming growth factor beta (TGFβ), and IL10 (Sartor, 2006). These altered expressions are associated with impaired mucosal barrier functionality and/or bacterial clearance at the epithelial interface causing an alteration in microbiota homeostasis, increasing luminal gut antigens, and leading to a life-long risk of inadequate and recurrent innate and adaptive immune activation, which implies the development of chronic tissue damage (Buttó and Haller, 2016).
Despite the remarkable accumulation of knowledge regarding the disease mechanisms associated with CD development in recent decades and the wide range of therapeutic options available, there is still a significant number of patients that do not respond to primary therapy or present lower response throughout treatment, e.g., anti-TNFα, which presents in up to 30%–40% of primary nonresponders and similar percentages of patients undergoing loss of response (Roda et al., 2016; Qiu et al., 2017; Adegbola et al., 2018). It is worth noting that currently available biologic therapies against CD target individual signaling pathways. Taking into account the large number of immune components and signaling pathways involved in CD pathogenesis, a single treatment may not be feasible for all of the patients with CD. In line with this, there is evidence suggesting that the optimal treatment of CD should involve a combination of different drugs targeting several pathways (Colombel et al., 2010; Grevenitis et al., 2015; Coskun et al., 2017).
The increasing global incidence of CD (Kaplan and Ng, 2017; Ng et al., 2018) together with the reported loss of effectiveness over time of standard of care requires a change in the paradigm of drug discovery and development by applying techniques that can help to find new targets and promising drugs to treat patients (Iyengar et al., 2012). Taking into account that many molecular pathways in both the immune system and CD’s etiopathology have already been described individually, modeling becomes an attractive option as it enables the integration of available data into a computational framework to describe the dynamics of the disease, support target identification, and allow the in silico evaluation of different therapeutic strategies.
Recently, a network model integrating key elements of IBD and the immune system has been developed in our group (Balbas-Martinez et al., 2018). Despite the capability of the model to predict a variety of real scenarios, it does not account for the magnitude or the temporality of the response to the simulated therapies given its Boolean nature.
Consequently, herein we aim to build a quantitative systems pharmacology (QSP) model able to characterize the dynamic behavior of plasma interleukins in CD using the previous network as a topological representation and starting scaffold (Balbas-Martinez et al., 2018) and incorporating those molecular pathways known to be relevant in CD for which human quantitative information is available. This modeling effort contributes to and expands the limited current arsenal of systems pharmacology models developed for IBD covering the as-yet small number of interleukins (Yates et al., 2004; Wendelsdorf et al., 2010; Lo et al., 2013; Dwivedi et al., 2014; Mei et al., 2015; Hontecillas et al., 2016; Peters et al., 2017).
In this work, the workflow for robust application of systems pharmacology modeling proposed by Gadkar et al. (2016) was adapted to the CD scenario (Fig. 1). Once the objective (mentioned above) was formulated, the steps considering scope definition and mathematical representation of biology were established as detailed in the Materials and Methods section. The Results section presents the ability of the model to capture the biologic behavior as well as the outcome from the sensitivity analysis (SA). Finally, the last step of the workflow is addressed in the Discussion section.
The iterative workflow followed for quantitative systems pharmacology modeling of CD. The text included in the blue section describes the conceptual objective of each stage, and the corresponding technical objective is summarized in the adjacent gray squares. The workflow is iterative, and model-based insights of different nature and degrees of robustness can be obtained at each stage. Figure adapted from Gadkar et al. (2016).
Materials and Methods
Scope Definition
Plasma steady-state levels of each of the ILs and cells included in the previously developed network (Balbas-Martinez et al., 2018) were obtained for healthy subjects and patients with CD based on a literature search. Medical Subject Headings terms were used during the review in different search engines such as PubMed, clinicaltrials.gov, and Google Scholar. Interestingly, a recent publication (Rogers et al., 2018) provided information for several of the species under study, which was complemented with published clinical trials (Supplemental Table 1). The average plasma concentration (
) was calculated by weighting the different reported average concentrations by the number of individuals per study (eq. 1).
(1)where 
and 
refer to the average plasma concentration of the species of interest and the number of individuals reported in the ith study, respectively.
Parameters characterizing model interactions among components were obtained either directly from the literature, when available, or estimated using human in vitro or in vivo data presented in tables or figures from the literature, always when a clear description of experimental conditions was provided. A detailed description of the data used for the parameter optimization of each interleukin is summarized in the Supplemental Methodology and Results. For modeling purposes, data shown in figures were extracted using the WebPlotDigitalizer version 3.8 software.
Representing the Biology of CD
Visual Note-Taking Diagrams.
Visual note-taking diagrams (i.e., topological representations) were generated for each of the interleukins making up the identified interactions (Balbas-Martinez et al., 2018). Then, these diagrams were used to develop the corresponding mathematical representation.
Model Parametrization and Ordinary Differential Equation Definition.
The general ordinary differential equation (ODE) describing the interleukins dynamics is represented in eq. 2.
(2)where dILi/dt is the rate of change of the ith interleukin, its degradation is governed by the first-order rate constant 
, and the synthesis is controlled by the zero- or first-order rate constants, 
and 
, respectively. The choice of zero- or first-order rate constant was dependent upon whether the production was considered to be constant or proportional to cell levels, respectively.
The synthesis processes can also be modulated by other interleukins in the environment (jth interleukin) as represented by the vectors 
and 
with the next general expression 
, where the structure of 
refers to either a stimulation or an inhibition of synthesis following linear, sigmoidal EMAX, or exponential dynamics by the jth interleukin.
The values used for 
were derived (
) from the human half-life values reported in the literature for the different ILs. In the case of IL22, IL17 and IL18 no human data were found, thus their half-lives were allometrically scaled from animals according to the following equation: 
(Boxenbaum, 1982).
Parameter estimation and selection of the structure for 
were accomplished using the nonlinear regression function available in R version 3.5.0 and the Akaike information criterion (Hirotogu, 1998). Each tested model was selected according to the richness of the available data, accounting for parameter identifiability. For instance, in the case where only two data points were available, only a linear model was tested. On the other hand, when richer data were available, different submodels were tested, and the best one was selected. A detailed description of the data used for the parameter estimation of each interleukin is summarized in the Supplemental Methodology and Results.
The rate constants of synthesis were derived at the steady-state condition (dILi/dt = 0) using the steady-state levels extracted from the literature (Table 1), the 
values, and the parameter estimates for the selected 
functions, except for the synthesis rate of IFNγ elicited by CD4+IFNγ+ T (Th1) and CD4+IL17+ T (Th17) cells, the values of which were found in the literature (Kamada et al., 2008; Axtell et al., 2010; Liu et al., 2011).
List of the mean value, with S.D. or range, of the model species’ circulating levels for healthy subjects and patients with CD
Healthy condition.
The dynamic of the interleukins was characterized as explained, by processes of synthesis and degradation and the initial conditions reported in Table 1. In the healthy state, levels of activated antigen-presenting cells (aAPCs) were considered negligible, and therefore the aAPC-mediated processes were not considered at this step.
Crohn's disease condition.
For most of the ILs in the model, the rate constants of synthesis were derived using the initial condition for aAPCs and ILs reported for patients with CD with the rest of the parameters obtained in healthy subjects. On the other hand, for the case of IL4, IL17, and IL22, whose synthesis processes are not modulated by aAPCs, the initial conditions were calculated using the CD reported levels of the corresponding producing cells and the synthesis and degradation parameters used for the healthy condition.
To summarize, the types of data considered were 1) baseline levels of model components reported in literature (mean and range), 2) in vitro dynamics extracted from literature used to obtain model parameter estimates, and 3) overall noncontinuous measures of clinical response expressed as CD activity index score and safety, used to be compared on qualitative means with the outcome obtained from model simulations.
Capturing Behavior
Model Curation.
The systems pharmacology model was implemented in the MathWorks SimBiology toolbox (MATLAB version R2018b) (MATLAB, 2018). The model, which is not fit-for-purpose, was challenged by running unperturbed deterministic simulations and comparing the model predicted steady-state levels of the different interleukins with those obtained from the literature (Table 1) for healthy and CD conditions.
Additionally, 1000 simulations were performed in which the initial conditions of the different components (ILs and cells) were randomly generated using log-normal distributions with the mean and S.D. provided in the literature (Table 1). Agreement between simulated and reported values (i.e., means and range) was visually evaluated.
Sensitivity Analysis.
A local SA evaluating the change in steady-state levels of each of the ILs as a function of variations in the model parameters was performed in MATLAB version R2018b, using the complex-step approximation method (Zhou et al., 2015; MATLAB, 2018). The integrals of the fully normalized time-dependent sensitivities from day 0 to day 28 were computed and reported as a heatmap.
Therapy Simulation.
Recombinant human IL10 (rhuIL10) was chosen to be tested as therapy due to the high relevance of IL10 in many of the model ILs. The levels of the different ILs in the model were simulated and evaluated during a virtual treatment with IL10 administered as an intravenous bolus once daily over 7 consecutive days. The simulated doses allowed maximum concentrations of IL10 of 1 and 6 ng/ml to be achieved, thus mimicking the exposure to rhuIL10 reached in a clinical study in patients with CD (Fedorak et al., 2000).
Results
Species Selection and Blood Levels
A total of 21 species representative of the innate and adaptive immune response in CD were included in the model considering both cells and interleukins. The cells included were aAPCs [activated macrophages (MACRact) and activated dendritic cells (DCact)], activated CD4+ T (Th0act) cells, CD4+IL4+ T (Th2) cells, Th1 cells, Th17 cells, and CD4+CD25+FOXP3+ T (Treg) cells. On the other hand, the ILs included were IL4, IL10, TGFβ1, IL12, IFNγ, IL18, IL2, IL17, IL23, IL6, IL22, IL15, TNFα, and IL1β. The average plasma concentration, range, and demographics for each species for healthy subjects and patients with CD along with the assumptions made and the corresponding set of references are reported in Table 1 and Supplemental Table 1.
Representing the Biology of CD
Visual Note-Taking Diagrams.
Figure 2 provides an overall overview of the model components and their relationships consisting of 21 species (seven cells and 14 interleukins). Furthermore, Fig. 3 shows in detail the dynamics of IL12 as an example, where its synthesis (both baseline and mediated by aAPCs) is promoted by IFNγ and inhibited by TNFα, IL4, IL10, or TGFβ1. The topological representations for all the IL submodels considered, including model assumptions, are presented in the Supplemental Material and have served as a basis for the development and implementation of the corresponding ODEs and integration into the entire QSP model.
Graphical representation of the CD model. Cells, represented by colored figures adapted from Servier Medical Art (smart.servier.com), are in charge of the synthesis of the interleukins included in the model, represented as small spheres. The interleukins promoting stimulation or inhibition of the synthesis process are denoted in green or red, respectively. DC, dendritic cell; MACR, macrophage.
Graphical representation of IL12 kinetics. DC and MACR represent activated dendritic cells and macrophages, respectively. Synthesis stimulation is represented by green arrows, and inhibitions are represented by red lines with a round edge.
Model Parametrization and ODE Definition.
The QSP model includes 14 ODEs, 1 for each IL, and 84 parameters accounting for the different processes. From the 84 total parameters, 45 were estimated from in vitro data (Supplemental Methodology and Results), 14 were directly obtained from the literature (degradation constants), and the remaining 24 were derived as previously explained. Supplemental Table 2 lists the values for all the model parameters. A full description of the development and parameterization of the ODEs for all the ILs of the model is given in the Supplemental Material. Below, the process followed for one of the ILs of the model (IL12) is detailed as an example. Figure 4 shows the submodels fitted for the IL12 influencing components, where the raw data and the model predictions are presented, along with the constructed ODE that characterizes IL12 dynamics. The IL12 model was developed based on the following information and assumptions: 1) synthesis is inhibited by TNFα, IL4, IL10, and TGFβ1 (Bonder et al., 1999; Ma et al., 2000; Xia and Kao, 2003; Fogel-Petrovic et al., 2007) and stimulated by IFNγ (Xia and Kao, 2003); 2) aAPCs (DCact and MACRact) are responsible for the increased levels of IL12 in patients with CD with respect to the healthy condition in which aAPCs are considered negligible (Bonder et al., 1999; Ma et al., 2000; Xia and Kao, 2003); and 3) the magnitude of the effect of different ILs on the synthesis mediated by the two types of aAPCs is considered the same.
Graphical and mathematical representation of the models characterizing the influence of the different ILs on IL12 synthesis. The dashed blue lines refer to the model predictions, and black points to the in vitro data obtained from the literature after being normalized to represent the effect of each specific interleukin in picograms of IL12 synthesis per milliliter day and cell. (A) IL4 inhibition (Bonder et al., 1999), characterized by a linear model. (B) IL10 inhibition (Xia and Kao, 2003), characterized by a sigmoidal IC50 model. (C) TGFβ1 inhibition (Fogel-Petrovic et al., 2007), characterized by an exponential model. (D) TNFα inhibition (Ma et al., 2000), characterized by a sigmoidal Imax model. (E) IFNγ stimulation (Xia and Kao, 2003), characterized by a linear model.
Additionally, the full model structure is provided in a sbml and sbproj format (CD_Interleukins.xml and CD_Interleukins.sbproj) ensuring reusability and reproducibility.
Capturing Behavior
Model Curation.
Results from the deterministic simulation (in both healthy and CD conditions) showed that the simulated steady-state levels matched well the levels of ILs reported in the literature and used as initial conditions to derive ILs dynamics (data not shown). Additionally, the simulation results showed that randomly varying the different components of the model within their reported physiologic ranges did not produce nonphysiologic steady-state levels for any of the ILs in both CD (Fig. 5A) and the healthy condition (data not shown). The agreements between the mean predicted and mean observed ILs values are shown in Fig. 5B, where it can be observed that for all the entities comprising the model the corresponding data pairs fell within the 1.33-fold limits of acceptance.
Results of 1000 simulations assigning random values for each model component. In (A) the simulated interleukin distribution is represented in light green color with the mean (solid line) and the 2.5–97.5th percentile of the simulations (dashed lines). Red color represents the average (solid line) and the distribution (dashed lines) of the reported concentrations in the literature for patients with CD (Table 1). Left panel of (B) indicates the predicted versus observed average values. Solid line represents the unity line and dashed lines represent the 1.33-fold prediction line. The table included lists the actual mean values used to construct the graphic, together with the absolute prediction error (PE%) calculated as 
Sensitivity Analysis.
Figure 6 summarizes the results from the local SA. None of the ILs showed a high degree of sensitivity to changes in any of the parameters, supporting biologic stability. The most relevant parameters affecting interleukins dynamics were those directly related to degradation (intense dark red colored cells in Fig. 6). Consequently, IFNγ, IL10, and TNFα degradation-related parameters seemed to be the parameters exerting the greatest impact in the model, as they interfere in the synthesis of many interleukins. As expected, interleukins whose synthesis is highly regulated by other interleukins (TNFα, IFNγ, and IL12 rows) showed a higher influence on the parameters governing their turnover.
Local sensitivity analysis. The heatmap indicates the effect of the alteration of a single parameter (columns) in every interleukin of the model (rows). The color in each cell corresponds to the result of the sensitivity of the parameter alteration for each interleukin, ranging from null (purple color, close to 0) to high influence (red color, close to 28).
Therapy Simulation.
The simulated profiles of the 14 ILs under the healthy condition and nontreated and treated CD condition after administration of two virtual doses of rhuIL10 as described previously are presented in Fig. 7. With the lowest and highest doses, equivalent to 5 and 20 μg/kg, respectively, a minor impact on IL12 and IL18 was observed as a reduction in their levels toward the healthy condition. Regarding IFNγ and TNFα, the average concentration oscillated around the healthy level for the lower dose, which could be correlated to the optimal response. Moreover, further decreases in the average ILs steady-state levels (below the healthy condition levels) are shown for the higher simulated doses (20 μg/kg) (Fedorak et al., 2000).
Results from the rhuIL10 simulated therapy. The blue and orange lines represent the IL simulated profiles after administration of 5 and 20 µg/kg rhuIL10 doses, respectively, given as a bolus for 7 consecutive days, starting at day 7. The red and green lines represent the Crohn's disease and healthy population average IL concentration obtained from the literature (Table 1).
Discussion
In the current study, we present a QSP model characterizing the dynamics of the relevant ILs in CD in plasma. Our analysis appears timely, as CD incidence has dramatically increased worldwide (Ng et al., 2018) and the disorder has recently attracted more attention, as demonstrated by the increasing number of CD research studies and drug development programs (Peters et al., 2017; Gregoire et al., 2018; Li et al., 2018; Klaassen et al., 2019; Ye et al., 2019).
Different efforts have been undertaken in the past to describe quantitatively the different parts of the IBD system (Wendelsdorf et al., 2010; Lo et al., 2013; Dwivedi et al., 2014; Mei et al., 2015; Hontecillas et al., 2016). The considerable number of ILs considered and the interaction among them, along with the focus just on the CD condition of IBD and our effort to obtain parameter estimates exclusively of human origin, represent the main contributions of our work with respect to others. Furthermore, the amount of updated literature compiled for the development of our model not only represents a detailed open and useful repository, but also highlights current knowledge gaps on the interaction between cytokines, suggesting future experiments that could be performed to feed and refine the existing model.
The fact that the information gathered from the literature was obtained under very different experimental designs/conditions represents a challenge with respect to integration for model building. Therefore, to mitigate the possible sources of variability, the methodology employed followed the workflow proposed by Gadkar (Gadkar et al., 2016), which aimed to standardize the building of QSP models.
A number of assumptions were made throughout the process of model development, with the following being the most critical. 1) CD is a disease that occurs in the gastrointestinal tract. We assume that the plasma circulating interleukins are correlated to those in the affected areas, as previously demonstrated for some of the ILs considered in the current model (Braegger et al., 1992; Arnott et al., 2001; León et al., 2009). 2) We used in vitro data as a surrogate of in vivo system dynamics. 3) No interaction (i.e., antagonism or synergism) was assumed between the different ILs modulating a particular synthesis process, and solely additivity scenarios were considered. 4) Linear, exponential, or EMAX models were used to describe IL effects on the different synthesis processes. We should acknowledge that although some of these relationships might not be physiologically plausible, they can certainly appear within the biologic ranges of the cytokines. Moreover, no indication of model misspecification or implausible simulated values was observed when challenging the model. Given the abovementioned assumptions, deterministic and stochastic simulations were performed to challenge the QSP model and evaluate the consistency of the model structure and the parameter values. The simulated steady-state values obtained after randomly and simultaneously varying the initial conditions of the different ILs were located within reported physiologic plasma ranges for CD condition. Nonetheless, the predicted ranges were in general narrower than those reported in the literature. This could be due to the fact that observations from several studies using different analytical methods and degrees of disease were pooled, but also to the fact that the proposed QSP model does not account for the dynamics of the different cells yet, which could also influence (increasing) predicted variability. Additionally, neither interindividual variability nor residual error was added to the different parameters, which would otherwise have resulted in an increase in the simulated range. Finally, the local SA (Fig. 6) demonstrated that the impact of changing the parameters of the interleukin levels was not large, thus confirming the stability of the modeled biologic system. Overall, these results support model robustness and indicate that the assumptions made at this stage were reasonable. Unfortunately, taking into account that the sparsely available clinical data were primarily used for model training, no external validation was performed, and predictions based on the model should be taken cautiously.
Regarding cell components, the model has been established in steady-state conditions where the levels of the IL-expressing cells (Th0act, Th1, Th2, Th17, Treg, and aAPCs) were assumed to be constant in the two conditions studied (the CD and healthy state). Therefore, the transition from the healthy to CD status and vice versa is not accounted for dynamically by the model in its present form but is rather a consequence of modifying the cells’ values based on the condition to emulate. However, this structural characteristic does not prevent the use of the proposed model to mimic IL-based therapies already tested in CD.
In line with this, the Crohn's disease activity index score combined with biomarkers, such as C-reactive protein or fecal calprotectin (Fengming and Jianbing, 2014), are currently used to monitor disease status and treatment efficacy. Interestingly, the current model provides the opportunity to use the dynamics of additional components associated with CD that could contribute to evaluate therapeutic effects and disease status. In fact, the results from the review performed by Rogers (Rogers et al., 2018) showed differences in circulating levels of the different ILs considered in our model between healthy patients and patients with CD, opening up the possibility to at least consider a wider set of potential biomarkers.
Indeed, the proposed model was used to test a virtual therapy with rhuIL10 as shown in the Results section (Fedorak et al., 2000), where the results obtained are in concordance with the general outcome of the clinical study. The simulation of a virtual dose equivalent to 5 and 20 µg/kg with our model provided return of main IL values to healthy levels. On the other hand, the levels for IFNγ and TNFα resulted in further decreases in their steady-state levels (below the healthy condition levels) for both doses, resulting in lower levels for the higher simulated dose (20 μg/kg). Accordingly, the authors concluded that rhuIL10 showed clinical remission and endoscopic improvement at the dose of 5 µg/kg, whereas higher doses showed a reduced response and limited safety (Fedorak et al., 2000). It might well be that the lack of efficacy reported by the authors could be the result of a higher incidence of adverse effects, as some of these are part of the primary and secondary measurement of efficacy, thus jeopardizing efficacy. Other therapies such as TNFα, IL12–23, or the IL17 or IFNγ monoclonal antibodies are approved or have already been tested in a clinical setting against CD (Hommes et al., 2006; Hueber et al., 2012; European Medicines Agency and European Medicines Evaluation Agency, 2016, 2018). Even though the targeted interleukins of these therapies are implemented in the model, their pathways are more related to cell dynamics (T-cell differentiation processes or natural killer activation/response (Balbas-Martinez et al., 2018)), which are mechanisms beyond the scope of the current IL QSP model. Therefore, the present model brings certain confidence to simulations involving those pathways highly involved in IL synthesis, such as IL10, and needs further development to improve the confidence in therapy simulations related to T-cell differentiation or natural killer activation/response pathways.
On the other hand, one of the characteristics that modeling offers is reusability, which is even more relevant in the QSP arena as it is linked with the concept of translational research. In the current CD model, the parameter values were obtained under healthy and CD conditions, and therefore, it is likely that those values will change when evaluating other diseases where the immune system also plays a relevant role.
In summary, a QSP model describing the dynamics of circulating ILs in humans with CD has been successfully developed and implemented computationally in SimBiology. This manuscript contributes with detailed submodels for each interleukin and related parameter values to be applied to CD, thus providing a repository of great value by integrating the most relevant up-to-date literature. As the model is not fit-for-purpose, and taking into account that the sparsely available clinical data were primarily used for model training, no external validation was performed, and predictions based on the model should be taken cautiously. Finally, as any QSP model, the current framework is a living project in which the dynamics of the cells are being considered, thus expanding the applicability of the network to evaluate therapeutics acting on those specific pathways.
Authorship Contributions
Participated in research design: Balbas-Martinez, Asin-Prieto, Troconiz.
Conducted experiments: Balbas-Martinez.
Contributed new reagents or analytic tools: Balbas-Martinez, Asin-Prieto.
Performed data analysis: Balbas-Martinez, Asin-Prieto, Parra-Guillen, Troconiz.
Wrote or contributed to the writing of the manuscript: Balbas-Martinez, Asin-Prieto, Parra-Guillen, Troconiz.
Footnotes
- Received May 30, 2019.
- Accepted December 2, 2019.
This work was supported by a fellowship grant to Violeta Balbas-Martinez from the Navarra Government, Spain [944/2016: BON N° 18, de 26 de enero de 2017].
↵
This article has supplemental material available at jpet.aspetjournals.org.
Abbreviations
- aAPC
- activated antigen-presenting cell
- CD
- Crohn's disease
- DCact
- activated dendritic cells (in the presence of antigen)
- IBD
- inflammatory bowel disease
- IFNγ
- interferon gamma
- IL
- interleukin
- MACRact
- activated macrophages (in the presence of antigen)
- ODE
- ordinary differential equation
- QSP
- quantitative systems pharmacology
- rhuIL10
- recombinant human IL10
- SA
- sensitivity analysis
- TGFβ1
- transforming growth factor beta 1
- Th0act
- activated CD4+ T cells
- Th1
- CD4+IFNγ+ T cells
- Th17
- CD4+IL17+ T cells
- Th2
- CD4+IL4+ T cells
- TNFα
- tumor necrosis factor alpha
- Treg
- CD4+CD25+FOXP3+ T cells
- Copyright © 2020 by The American Society for Pharmacology and Experimental Therapeutics
