Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (3): 359-366

#### The article information

Hai-shan Yu, Zhi-chao Pan, Jie-lou Liao

Network Modeling of Inflammatory Dynamics Induced by Biomass Smoke Leading to Chronic Obstructive Pulmonary Disease

Chinese Journal of Chemical Physics, 2018, 31(3): 359-366

http://dx.doi.org/10.1063/1674-0068/31/cjcp1801001

### Article history

Accepted on: February 26, 2018
Network Modeling of Inflammatory Dynamics Induced by Biomass Smoke Leading to Chronic Obstructive Pulmonary Disease
Hai-shan Yu, Zhi-chao Pan, Jie-lou Liao
Dated: Received on January 9, 2018; Accepted on February 26, 2018
Department of Chemical Physics, University of Science and Technology of China, Hefei 230026, China
*Author to whom correspondence should be addressed. Jie-lou Liao, E-mail:liaojl@ustc.edu.cn
Abstract: Chronic obstructive pulmonary disease (COPD) is a chronic inflammatory disorder characterized by airflow obstruction and progressive damage of lung tissues. As currently more than 3 billion people use biomass fuel for cooking and heating worldwide, exposure to biomass smoke (BS) is recognized as a significant risk factor for COPD. Recent clinical data have shown that BS-COPD patients have a Th2-type inflammatory profile significantly different from that in COPD induced by cigarette smoke. As COPD is essentially proinflammatory, however, the mechanism underlying this Th2-type anti-inflammatory profile remains elusive. In this work, a network model is applied to study BS-induced inflammatory dynamics. The network model involves several positive feedback loops, activations of which are responsible for different mechanisms by which clinical phenotypes of COPD are produced. Our modeling study in this work has identified a subset of BS-COPD patients with a mixed M1-and Th2-type inflammatory profile. The model's prediction is in good agreement with clinical experiments and our in silico knockout simulations have demonstrated several important network components that play an important role in the disease. Our modeling study provides novel insight into BS-COPD progression, offering a rationale for targeted therapy and personalized medicine for treatment of the disease in future.
Key words: Network model    Inflammatory dynamics    Positive feedback loops    Biomass smoke    Chronic obstructive pulmonary disease
Ⅰ. INTRODUCTION

Chronic obstructive pulmonary disease (COPD) is a progressive inflammatory condition characterized by airflow limitation due to small airway obstruction (bronchitis) and destruction of the lung parenchyma (emphysema) [1, 2]. COPD has currently become the third leading cause of death globally [2]. There is no cure available for COPD and the present drug treatments are mainly effective in the improvement of symptoms and exacerbation but generally do not slow down the disease progression [3]. The disease poses a huge public health burden worldwide. Thus, it is of great importance to elucidate the detailed mechanism of the disease for the development of effective therapies.

Although inhalation of cigarette smoke (CS) is the primary cause of COPD (CS-COPD), a growing number of other risk factors such as exposure to air pollution particulate matter (PM) contribute to the pathology of the disease [4]. In particular, PM$_{2.5}$, defined as fine particulate matter with an aerodynamic diameter less than 2.5 ${\rm{\mu }}$m, plays a detrimental role in the pathogenesis of COPD, as it can readily penetrate into the small airway and alveoli of the lung [5]. There are several sources including biomass smoke (BS), a major contributor to PM$_{2.5}$ [6], as more than 3 billion people use biomass fuel for cooking and heating worldwide [7]. Exposure to BS has been widely recognized as a significant risk factor for COPD [1, 8]. BS-induced COPD (BS-COPD) has been considered mainly in developing regions of the world, but biomass burning is also recognized as a significant cause of the disease in industrialized countries [8].

BS contains many compounds similar to those in CS [9]. Both BS and CS activate inflammatory responses of the lung. The BS- and CS-induced inflammatory responses involve both innate and adaptive immunity, which are mediated through a complex network consisting of multiple immune cell types, molecular mediators, and lung tissues [1]. It would be expected that BS induces an inflammatory pattern similar to that of CS [1, 9]. Recent studies, however, have shown that immune response caused by exposure to BS [10] appears to be different from that induced by CS in COPD patients [11]. Clinical data have demonstrated that a set of patients with BS-COPD have higher levels of IL4-producing Th2 cells, which predominate over Th1 and Th17 [10], while subjects with CS-COPD develop a proinflammatory Th17-type immune profile, in which the Th17 population is significantly higher than that of Th2 [11]. However, a more recent study showed that another set of patients with BS-COPD exhibited higher levels of proinflammatory biomarkers than healthy subjects similar to those with CS-COPD, and there was no predominant Th2-type inflammation observed in BS-COPD subjects [8]. Although the discrepancies between these studies [8, 10] might be attributable to the high heterogeneity of COPD and a different distribution of males and females in each COPD group selected [8], the underlying mechanism is not completely clear.

In our previous study [12, 13], a network model was proposed to describe the dynamics of CS-induced immune responses in COPD progression. This network model is then extended to study the immune response dynamics in the progression of ulcerative colitis [14]. Our network model studies have identified several positive feedback loops (PFLs), activation of which plays a determinant role in several different mechanisms (endotypes) in COPD developments [12-14]. Inhibition of key elements in the activated PFLs could provide a possible therapeutic approach for COPD treatment. Moreover, lessons learned from our previous study are that a similar clinical phenotype of COPD patients may originate from different endotypes. Therefore, it is suggested that personalized medicine is required for an effective treatment of COPD [12, 13]. Our previous modeling predictions [13] were in good agreement with the clinical data for CS-induced COPD, in which Th17 predominates over Th2 [12]. As COPD is essentially proinflammatory [1, 12, 13], how an inflammatory profile of a Th2-type, which is anti-inflammatory, is achieved in a subset of patients with BS-COPD has not been addressed in our previous study [12] and remains largely unclear so far. In the present work, the network model developed in our previous study [12] is applied to elucidate BS-induced inflammatory dynamics leading to COPD to address the above issue.

Ⅱ. MATERIALS AND METHODS A. Network model

The immune system associated with COPD is highly complex, involving many molecular mediators, immune cells, and lung tissues. For the sake of simplicity, a network model is developed by treating important components as network nodes [12]. In this network model, there exist two types of inputs initiated from a node. A positive or an up-regulation input (denoted by $\rightarrow$) represents that an increasing of the concentration of the tail node will result in an increasing of the head node or an up-regulation of the process when the input arrow ends at an edge between two nodes, and vice versa for a negative or a down-regulation (inhibition) input (denoted by $\dashv$).

The network model (FIG. 1) used in this work is based on that presented in our previous study [12]. Upon exposure to external stimuli such as BS, immature macrophages (M0) are polarized to the inflammatory type, M1 macrophages, and initiate an inflammatory cascade [7]. M1 produces inflammatory cytokines such as tumor necrosis factor-$\alpha$ (TNF-$\alpha$), which activates M1 conversely, IL-6 and IL-12 [15]. M1 can cause tissue damage (TD) in the lung by releasing reactive oxygen species (ROS) leading to oxidative stress, proteases such as macrophage elastase and metalloproteases (MMPs) to ingest pathogens and apoptotic cells, and chemokines including IL-8 to recruit neutrophils into the lung [16, 17, 18]. Neutrophils, which are short-lived and subsequently cleared by macrophages [19], contribute to TD in a manner similar to M1 macrophages [17]. In addition, the inflammation in COPD is often described as neutrophilic [1]. Furthermore, the tissues damaged by M1 can produce elastin fragments (EFs) as strong attractors to recruit monocytes (precursors of macrophages, M0) into the lung from circulation [19]. Subsequently, these M0 cells are differentiated into M1, thus forming a positive feedback loop, M1$\rightarrow$TD$\rightarrow$M1.

 FIG. 1 Network model for BS-induced inflammatory response. The interactions between various nodes that represent immune cells, cytokines and tissues (TD) are depicted (see the main text for details). For clarity, only a few inhibitory interactions of IL-10 are shown (see details in Ref.[12]).

Lung tissue damage triggers an early wounding healing process by producing IL-4 that alternatively activates macrophages (M2) [18, 20, 21]. M2 secretes IL-10, which activates M2 inversely, and transforming growth factor, TGF-$\beta$ [22]. While IL-10 is a potent anti-inflammatory cytokine that down-regulates almost all important proinflammatory and TD-related processes, TGF-$\beta$ is a multi-functional growth factor. In the lung parenchyma, TGF-$\beta$ down-regulates tissue damage through inhibition of MMP-12 and MMP-9, whereas in small airways, TGF-$\beta$ is a potent inducer for extracellular matrix target genes such as collagens, and fibroblast proliferation and activation which both are key events in the fibrogenic process [23].

Dendritic cells (DCs) are antigen-presenting, playing a critical role in linking the innate to the adaptive immune response [24]. Immature dendritic cells (DC0) near the epithelial surface are activated directly by BS or dangerous signals generated from TD [1, 11]. DC undergoes a maturation process and migrates towards the local lymph nodes. Naïve, quiescent T cells cannot enter the lung parenchyma. But once activated by matured DC [25, 26], they can move into the lung and differentiate into Th1, Th2, Th17, T-regulatory (Treg) and CD8$^{+}$ T cells in their corresponding cytokine environments, e.g., in the presence of IL-12 secreted by M1 (as well as DC), naïve CD4$^{+}$ T cells (Th0) differentiate into T helper 1 (Th1) cells [27, 28, 29]. Th1 secretes interferon-$\gamma$ (IFN-$\gamma$) to up-regulate the polarization process from M0 to M1 [30]. A multi-node positive feedback loop, M1$\rightarrow$IL12$\rightarrow$Th1$\rightarrow$IFN$\rightarrow$$\gamma$$\rightarrow$M1, is thus created. In contrast to Th1, Th2 is polarized from Th0 in the presence of IL-4. Release of IL-4 further enhances the production of IL-10 and TGF-$\beta$ by M2 (FIG. 1).

In the presence of TGF-$\beta$, Th0 cells differentiate into Treg, which secretes IL-10 [31]. TGF-$\beta$ and IL-6 together induce Th17 differentiation, leading to the production of IL-17 [32, 33, 34, 35]. While IL-17 acts on epithelial cells to recruit neutrophils to cause TD further, the activated epithelial cells in TD secrete IL-6, forming another positive feedback loop, IL-6$\rightarrow$Th17$\rightarrow$IL-17$\rightarrow$TD$\rightarrow$IL-6 [27, 35]. Th17 cells also produce IL-21 for the differentiation of CD8$^{+}$ T cells from naïve CD8$^{+}$ cytotoxic T lymphocytes (T0) [27, 36]. While CD8$^{+}$ T cells produce IFN-$\gamma$ to enhance the M1 inflammatory activities, they also release granzyme B and perforins, causing apoptosis/necrosis of targeted cells and leading to TD further [27]. In addition, IL-6 can down-regulate the activation of Treg that secretes IL-10 inhibiting Th17 [37]. Consequently, a positive feedback loop, IL-6$\dashv$ Treg$\rightarrow$IL-10$\dashv$ Th17$\rightarrow$IL-17$\rightarrow$TD$\rightarrow$IL-6, is generated. The aforementioned molecular mediators, immune cells, and TD are then treated as nodes, and are integrated into the network model presented in FIG. 1.

In this network model M1, DC, Th1, Th17, and CD8$^{+}$ T with their associated cytokines, TNF-$\alpha$, IL-6, IL-12, IFN-$\gamma$, and IL-17 form multiple proinflammatory pathways, whereas M2, Th2, and Treg with the related cytokines, IL-4, TGF-$\beta$, and IL-10, form anti-inflammatory/regulatory pathways. The inflammatory and anti-inflammatory/regulatory pathways are interlinked with each other through several nodes representing molecular mediators such as IL-6, TGF-$\beta$, IL-10, IL-4, and IFN-$\gamma$ (FIG. 1). These pathways eventually converge at the TD node. As we focus on the immunologic aspects of BS-COPD, the TD node of the network is highly coarse-grained, involving neutrophil-induced tissue damage, epithelial cell injury and extracellular matrix degradation etc. [12]. For example, epithelial cell injury in lung tissue (TD) can release molecular mediators including IL-4, IL-33 and thymic stromal lymphopoitin (TSLP) to up-regulate Th2 and type Ⅱ innate lymphoid cells (ILCs), both of which secrete IL-4, IL-5, and IL-13 [1, 38]. Dysregulated expressions of these cytokines in the airway smooth muscle are associated with asthma [1]. For simplicity, ILCs and TSLP, IL-5, and IL-13 are not included in the present network model.

B. Network dynamics

The above constructed network has a multiple timescale nature. For example, while cytokine regulation of cellular function via signal transduction usually takes place on a sub-second timescale, cell production of cytokines takes minutes to hours [39]. Therefore, the cytokine regulation activity can be considered to be at steady state in the equation that describes the slow timescale activities of the cells. Thereby, a positive or a negative input can be modeled using an increasing or decreasing Hill function [12]. In this work, we use a set of ordinary differential equations (ODEs), which are similar to those in our previous study [12], to describe the dynamics of the above network components (see text in Supplementary materials). These ODEs involve 18 variables, i.e., M$_{1}$, M$_{2}$, D$_{\textrm{C}}$, T$_{1}$, T$_{2}$, T$_{8}$, T$_{17}$, and T$_{\textrm{g}}$ represent the densities of M1, M2, DC, Th1, Th2, CD8$^{+}$, Th17, and Treg cells (in units of cell numbers in a cubic millimeter of tissue), respectively, whereas I$_{4}$, I$_{6}$, I$_{10}$, I$_{12}$, I$_{17}$, I$_{21}$, I$_{\alpha}$, I$_{\gamma}$, and I$_{\beta}$ denote concentrations of the cytokines, IL-4, IL-6, IL-10, IL-12, IL-17, IL-21, TNF-$\alpha$, IFN-$\gamma$, and TGF-$\beta$. The variable, T$_{\textrm{D}}$, represents the tissue damage (in terms of a percentage) [12]. Given the network model (FIG. 1) as well as the ODEs (see text in the Supplementary materials), the values of the parameters in the ODEs determine the immune response dynamics and inflammatory profile of a subject exposed to BS. In the following discussion, the values of the parameters (Table S1 in the Supplementary materials) in the ODEs (see text in the Supplementary materials) were adopted or estimated from experimental data (Table S2 in the Supplementary materials). A subset of the parameters including k$_{9}$ (3.36$\times$10$^{4}$/day), k$_{26}$ (8.00 pmol/(cell day)) and k$_{27}$ (2.00 pmol/(cell day)), which govern the dynamics of Th2 and IL-4, respectively, have different values from their counterparts, k$_{8}$ (0.41$\times$10$^{4}$/day), k$_{\textrm{I4, TD}}$ (1.56 pmol/(cell day)) and k$_{\textrm{I4, T2}}$ (0.83 pmol/(cell day)) in Table S1 in Ref.[12] for patients with CS-COPD studied previously.

As mentioned above, 18 ODEs, which involve 18 variables, M$_{1}$, M$_{2}$, D$_{\textrm{C}}$, T$_{1}$, T$_{2}$, T$_{8}$, T$_{17}$, T$_{\textrm{g}}$, I$_{4}$, I$_{6}$, I$_{10}$, I$_{12}$, I$_{17}$, I$_{21}$, I$_{\alpha}$, I$_{\gamma}$, I$_{\beta}$, and T$_{\textrm{D}}$, are used to describe the network dynamics. These ODEs are listed in Supplementary materials. The system of ODEs is solved numerically using MATLAB (version R2013a Mathworks) with a variable order and multistep solver, ode15s, and the parameters used in the simulations are listed in Table S1. MATLAB is also used to plot the simulation data to generate the figures presented below.

Ⅲ. RESULTS A. Dynamics of BS-induced immune response and inflammatory profile

FIGs. 2-4 present the population dynamics of the immune cells and cytokines in response to BS, respectively. The M1 population (FIG. 2(b)) along with TNF-$\alpha$, IL-6, and IL-12 (FIG. 3(b)) ascends to a peak after BS exposure for $\sim$13 days, and then goes downuntil day 20 due to the down-regulation of IL-10 mainly produced by M2 (FIG. 1(a)), exhibiting an acute inflammatory response to BS exposure. As discussed in the previous study, this time period is referred to as phase Ⅰ in COPD progression [14]. Thereafter, the IL-10 inhibitory effect on M1 is countervailed by the M1 production up-regulated by TNF-$\alpha$, TD and IFN-$\gamma$ (see FIG. 1). M1 is then raised again up to day of 180. This time of period is referred to as phase Ⅱ, which bridges the innate and adaptive immunity in the progression of COPD [14]. During phase Ⅱ, DC along with IL-12, IL-6, and TGF-$\beta$ (FIG. 3(b)), and IL-21 (data not shown) increases slowly and gradually, leading to the slow productions of Th1, Th17, and CD8$^{+}$, respectively. After phase Ⅱ, TD (FIG. 5(a)) along with the immune cells and molecular mediators eventually reaches a steady state (stable COPD, see results in Table S2). Our modeling results are consistent with laboratory and clinical experiments (see Table S2 in the Supplementary materials).

 FIG. 2 Population dynamics of M$_{1}$, M$_{2}$ and D$_{\textrm{C}}$ over a time period of (a) 6000 days and (b) 180 days (the dashed square region in (a)) in BS-COPD progression.
 FIG. 3 Population dynamics of I$_{\alpha}$, I$_{6}$, I$_{12}$, I$_{17}$, I$_{\gamma}$, I$_{4}$, I$_{10}$ and I$_{\beta}$ over a time period of (a) 6000 days and (b) 180 days (the dashed square region in (a)) in BS-COPD progression.
 FIG. 4 Population dynamics of T$_{1}$, T$_{2}$, T$_{8}$, T$_{17}$, and T$_{\textrm{g}}$ over a time period of (a) 6000 days and (b) 180 days (the dashed square region in (a)) in BS-COPD progression.

Our simulation results show that overall, BS-induced dynamic behaviors of the innate immune cells, M1 and DC (FIG. 2) are similar to those induced by CS in COPD progression [12]. However, there exist significant differences in the dynamics of some important network components including IL-4 and Th2 between these two cases. As shown in FIG. 3 IL-4 (I$_{4}$) is increased dramatically compared to that in CS-COPD [14]. Th2 is also significantly enhanced, predominating over Th1 and Th17, respectively, during the COPD progression (FIG. 4). The BS-COPD patients, which were all female [12], studied in this work have an inflammatory profile in which the Th2 level is higher than that of Th17 (FIG. 5(b)), different from that of CS-COPD in which Th17 predominates over Th2 [13, 14] (see FIG. S1 in this work). Here, our simulations results are in good agreement with the clinical data [12], as shown in FIG. 5(b). Nevertheless, M1 still remains at a high level (FIG. 2) similar to that of CS-COPD [14] as the host is persistently exposed to BS and the M1-involved PFLs are constantly activated. As a result, our modeling study has identified a subset of COPD patients whose immune dynamics is dominated by a mixed M1- and Th2-type response.

 FIG. 5 (a) TD dynamics. (b) Calculated results (red) for Th1, Th2, Th17, and Treg compared with clinical data (black) for BS-COPD. The clinical data are taken from Ref.[12].
B. Knockout simulations

In the following discussion, in silico knockout simulations are performed to identify important network components for BS-COPD through deletion of a node by setting all parameters of the component and the rate to zero. The results for knockouts of the immune cells, M1, Th1, Th2, Th17, CD8$^{+}$, and the cytokines, TNF-$\alpha$, IL-6, IL-17, IFN-$\gamma$, and IL-4 are presented in FIG. 6 (a) and (b), respectively.

 FIG. 6 Simulations of TD dynamics for in silico knockout of (a) an immune cell, M$_{1}$ (red), T$_{1}$ (green), T$_{2}$ (blue), T$_{17}$ (cyan), or T$_{8}$ (pink) and (b) a cytokine, I$_{\alpha}$ (red slashed line), I$_{6}$ (cyan), I$_{17}$ (pink), I$_{\gamma}$ (green), or I$_{4}$ (blue). WT (wild-type) is represented by black solid line.

FIG. 6(a) shows that M1 knockout leads to a significant reduction in TD in spite of persistent exposure to BS. As discussed earlier, M1 not only produces proinflammatory cytokines such as IL-6 and IL-12 to activate the adaptive immune responses, but also induces TD by producing chemokines such as IL-8 to recruit neutrophils into the lung, ROS leading to oxidative stress, and elastolytic enzymes including MMPs. The M1 knockout result along with that shown in FIG. 3 demonstrates that M1 predominates in the immune response, consistent with experiments in which M1 has a determinant role for BS-COPD [1] similar to the case of CS-COPD [40]. While deletion of Th1 leads to a small change in BS-induced TD, knockouts of Th17 and CD8$^{+}$ result in a large amount of reduction in TD, respectively. Interestingly, deleting Th2 leads to an increase in TD (FIG. 6(a)), which is consistent with knockout of Th2-produicng IL-4 (FIG. 6(b)). These results are not surprising as Th2 and IL-4 execute their anti-inflammatory and wound healing effects on TD. Our modeling study indicates that anti-IL-4/Th2 strategy may not be effective in the treatment of BS-COPD, although both Th2 and IL-4 are significantly enhanced that might be associated with coexisting asthma [41, 42], in line with clinical experiments [38].

Although the TNF-$\alpha$ level (FIG. 3) is significantly increased in BS-COPD, no significant reduction in TD is found in the TNF-$\alpha$ knockout simulation (FIG. 6(a)). This result is consistent with clinical data showing that TNF-targeted therapy is ineffective in COPD treatment [3, 43]. As shown in FIG. 6(b), the deletions of IL-17 and IFN-$\gamma$ lead to a relatively small decrease (${\rm{<}}$$\sim$${\rm{10\%}}$) in TD, respectively. These low-level reductions are not contradictory to the effects of CD8$^{+}$ and Th17 knockouts, as the latter effects come from both cytokine deletion and elimination of the cytotoxicity of CD8$^{+}$ that produce granzyme B and perforins to cause TD further as discussed earlier. Despite a relatively low level of IL-6 in BS-COPD, the IL-6 knockout still results in a large reduction of TD at the steady state, demonstrating an important role of IL-6 in bridging the innate and adaptive immunity. Our knockout simulations indicate that IL-6 is a promising anti-inflammatory target for an effective treatment of COPD [3].

Ⅳ. DISCUSSION

COPD is associated with chronic inflammation that affects predominantly the lung parenchyma and airways leading to airway limitation. This inflammation is amplified in patients with COPD and persists even after exposure to smoke is stopped [1]. However, the precise mechanisms for the inflammatory amplification and persistence are not clear [1]. Recently, we proposed a network model to probe the cellular and molecular mechanisms of CS-induced COPD [12]. Our modeling study has identified several positive feedback loops, activations of which are responsible for such inflammatory amplification and persistence in COPD patients and have an important role in several distinct mechanisms (endotypes) by which clinical COPD phenotypes are developed [12]. Our previous modeling results indicate that similar clinical phenotypes of COPD patients can come from different endotypes, suggesting that personalized medicine is required for COPD treatment.

Recent clinical data have shown that a set of patients with BS-COPD have a Th2-type inflammatory profile, in which the levels of Th2 are higher than those of Th1 and Th17, respectively [10], significantly different from that with CS-COPD where Th17 predominates over Th2 [11]. As the immune response dynamics of CS-COPD was extensively studied in our previous work [12], the underlying mechanism by which BS-COPD patients develop a Th2-type inflammatory profile has remained elusive [1, 9-11]. In the present work, we employed the network model developed in our previous study [12] to investigate BS-induced inflammatory response in COPD progression. Our modeling study has identified a subset of patients with BS-COPD, whose immune response is of mixed M1- and Th2- type, in which M1 dominates M2, whereas Th2 predominates over Th1 and Th17, respectively, in good agreement with clinical data (FIG. 5(b)) [10]. In silico knockout simulations in this work have demonstrated several important network components that have an important role in BS-COPD.

It is of interest to note that BS-COPD occurs mostly in women [10] and while female BS-COPD patients have a Th2-type inflammatory profile [10], male patients with BS-COPD have an immune response similar to that of CS-COPD [8]. This modeling study has identified a subset of COPD patients (i.e., females) who are likely more sensitive to tissue damage in the lung and have developed a protective mechanism by which IL-4 is significantly enhanced in response to the signal from TD (see FIG. 1) in the early phase as shown in FIG. 3(b) [12]. The positive feedback loop, IL-4$\rightarrow$Th2$\rightarrow$IL-4 (FIG. 1), is then activated so that the Th2 level is significantly enhanced, predominating over Th1 and Th17, respectively (see FIG. 4). IL-4 also alternatively activates M2 that secretes TGF-$\beta$ and IL-10. While TGF-$\beta$ promotes tissue repair in the lung parenchyma as mentioned above, IL-10 attenuates the proinflammatory responses and reduces tissue damage. Intriguingly, in silico knockout of IL-4 or Th2 in this study leads to an increase of TD although IL-4 and Th2 are both significantly enhanced in patients with BS-COPD, implicating that anti-IL-4/Th2 therapy may not be effective in the treatment of BS-COPD. Our modeling study provides novel insight into the cellular and molecular mechanism of BS-COPD with a Th2 profile, providing a rationale for targeted therapy and the personalized medicine treatment of COPD in future.

Supplementary materials: Equations (S1)-(S7) are given for the population dynamics of cytokines, I$_{\alpha}$ (TNF-$\alpha$), I$_{12}$ (IL-12), I$_{\gamma}$ (IFN-$\gamma$), I$_{17}$ (IL-17), I$_{21}$ (IL-21), I$_{4}$ (IL-4), and I$_{\beta}$ (TGF-$\beta$). FIG. S1 shows calculated results (red) for Th1, Th2, Th17, and Treg compared with clinical data (black) for CS-induced COPD. TABLE S1 lists parameters for the equations describing network dynamics of BS-induced immune response. TABLE S2 lists cell density and cytokine concentrations at the steady state from the simulations compared to experiments.

Ⅴ. ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China (No.21273209).

Reference
 [1] P. J. Barnes, J. Allergy Clin. Immuno. 138 , 16 (2016). DOI:10.1016/j.jaci.2016.05.011 [2] S. I. Rennard, and M. B. Drummond, Lancet 385 , 1778 (2015). DOI:10.1016/S0140-6736(15)60647-X [3] P. J. Barnes, Nat. Rev. Drug Discov. 12 , 543 (2013). DOI:10.1038/nrd4025 [4] F. He, B. Liao, J. Pu, C. Li, M. Zheng, L. Huang, Y. Zhou, D. Zhao, B. Li, and P. Ran, Sci. Rep. 7 , 45666 (2016). [5] S. L. Hwang, S. E. Guo, M. C. Chi, C. T. Chou, Y. C. Lin, C. M. Lin, and Y. L. Chou, Int. J. Environ. Res. 13 , 366 (2016). [6] J. L. Lopez-Campos, E. Marquez-Martin, and J. B. Soriano, Expert. Rev. Respir. Med. 9 , 1 (2016). [7] J. Olloquequi, and R. O. Silva, Innate Immun. 22 , 373 (2016). DOI:10.1177/1753425916650272 [8] R. Golpe, I. Martín-Robles, P. Sanjuán-López, L. Pérezde-Llano, C. González-Juanatey, J. L. López-Campos, and E. Arellano-Orden, Int. J. COPD 12 , 2639 (2017). DOI:10.2147/COPD [9] R. Silva, M. Oyarzún, and J. Olloquequi, Arch. Bronconeumol. 51 , 285 (2015). [10] H. Solleiro-Villavicencio, R. Quintana-Carrillo, R. Falfán-Valencia, and M. I. Vargas-Rojas, Clin. Immunol. 161 , 150 (2015). DOI:10.1016/j.clim.2015.07.009 [11] M. I. Vargas-Rojas, A. Ramirez-Venegas, L. Lim, ó nCamacho, L. LimónCamacho, R. Hernández-Zenteno, and R. H. Sansores, Respir. Med. 105 , 1648 (2011). DOI:10.1016/j.rmed.2011.05.017 [12] Z. Pan, H. Yu, and J. L. Liao, PLoS One 11 , e0163192 (2016). DOI:10.1371/journal.pone.0163192 [13] (a) J. L. Liao, J. Immuno. Biol. 2, 119(2017). (b) J. L. Liao, Immunotherapy 3, 107(2017). [14] D. Wu, H. Shan, and J. L. Liao, Chin. J. Chem. Phys. , in press (2018). [15] C. E. Boorsma, C. Draijer, and B. N. Melgert, Mediators Inflamm. 2013 , 1 (2013). [16] R. D. Hautamaki, D. K. Kobayashi, R. M. Senior, and S. D. Shapiro, Science 277 , 2002 (1997). DOI:10.1126/science.277.5334.2002 [17] A. F. Ofulue, and M. Ko, Am. J. Physiol. 277 , L97 (1999). [18] D. M. Mosser, and J. P. Edwards, Nat. Rev. Immunol. 8 , 958 (2008). DOI:10.1038/nri2448 [19] R. A. Holloway, and L. E. Donnelly, Curr. Opin. Pulm. Med. 19 , 95 (2013). DOI:10.1097/MCP.0b013e32835cfff5 [20] P. Loke, I. Gallagher, M. G. Nair, X. Zang, F. Brombacher, M. Mohrs, J. P. Allison, and J. E. Allen, J. Immunol. 179 , 3926 (2007). DOI:10.4049/jimmunol.179.6.3926 [21] E. Brandt, G. Woerly, A. B. Younes, S. Loiseau, and M. Capron, J. Leukoc. Biol. 68 , 125 (2000). [22] R. Faner, T. Cruz, and A. Agusti, Expert Rev. Clin. Immunol. 9 , 821 (2013). DOI:10.1586/1744666X.2013.828875 [23] M. Königshoff, N. Kneidinger, and O. Eickelberg, Swiss Med. WKLY 139 , 554 (2009). [24] P. Stoll, M. Ulrich, K. Bratke, K. Garbe, J. C. Virchow, and M. Lommatzsch, Respir. Res. 16 , 19 (2015). DOI:10.1186/s12931-015-0174-x [25] G. R. Van Pottelberge, K. R. Bracke, I. K. Demedts, K. D. Rijck, S. M. Reinartz, C. M. van Drunen, G. M. Verleden, F. E. Vermassen, G. F. Joos, and G. G. Brusselle, Respir. Res. 11 , 35 (2010). DOI:10.1186/1465-9921-11-35 [26] H. Torres-Aguilar, M. Blank, L. J. Jara, and Y. Shoenfeld, Autoimmun. Rev. 10 , 8 (2010). DOI:10.1016/j.autrev.2010.07.015 [27] P. J. Barnes, Nat. Rev. Immunol. 8 , 183 (2008). DOI:10.1038/nri2254 [28] G. Trinchieri, S. Pflanz, and R. A. Kastelein, Immunity 19 , 641 (2003). DOI:10.1016/S1074-7613(03)00296-6 [29] D. R. Milich, S. F. Wolf, J. L. Hughes, and J. E. Jones, Natl. Acad. Sci. USA 92 , 6847 (1995). DOI:10.1073/pnas.92.15.6847 [30] P. J. Barnes, Am. J. Respir. Cell Mol. Biol. 41 , 631 (2009). DOI:10.1165/rcmb.2009-0220TR [31] Y. Y. Wan, and R. A. Flavell, Proc. Am. Thorac. Soc. 4 , 271 (2007). DOI:10.1513/pats.200701-020AW [32] A. Kimura, and T. Kishimoto, Eur. J. Immunol. 40 , 2830 (2010). DOI:10.1002/eji.200940115 [33] E. Bettelli, Y. Carrier, W. Gao, T. Korn, T. B. Strom, M. Oukka, H. L. Weiner, and V. K. Kuchroo, Nature 441 , 235 (2006). DOI:10.1038/nature04753 [34] P. R. Mangan, L. E. Harrington, D. B. O' Quinn, W. S. Helms, D. C. Bullard, C. O. Elson, R. D. Hatton, S. M. Wahl, T. R. Schoeb, and C. T. Weaver, Nature 441 , 231 (2006). DOI:10.1038/nature04754 [35] H. Ogura, M. Murakami, Y. Okuyama, M. Tsuruoka, C. Kitabayashi, M. Kanamoto, M. Nishihara, Y. Iwakura, and T. Hirano, Immunity 29 , 628 (2008). DOI:10.1016/j.immuni.2008.07.018 [36] M. C. Duan, Y. Huang, X. N. Zhong, and H. J. Tang, Mediators Inflamm. 2012 , 1 (2012). [37] M. Fujimoto, M. Nakano, F. Terabe, H. Kawahata, T. Ohkawara, Y. Han, B. Ripley, S. Serada, T. Nishikawa, A. Kimura, S. Nomura, T. Kishimoto, and T. Naka, J. Immunol. 186 , 32 (2011). DOI:10.4049/jimmunol.0903314 [38] K. F. Chung, Lancet 386 , 1086 (2015). DOI:10.1016/S0140-6736(15)00157-9 [39] U. Alon, Math. Biosci. 215 , 193 (2008). DOI:10.1016/j.mbs.2008.07.002 [40] R. D. Hautamaki, D. K. Kobayashi, R. M. Senior, and S. D. Shapiro, Science 277 , 2002 (1997). DOI:10.1126/science.277.5334.2002 [41] R. Golpe, P. Sanjuán-López, E. Cano-Jiménez, O. Castro-Añón, and L. A. Pérez-de-Llano, Arch. Bronconeumol. 50 , 318 (2014). [42] B.G. Cosio, J. B. Soriano, J. L. Lopez-Campos, M. Calle-Rubio, J. J. Soler-Cataluna, J. P. de-Torres, J. M. Marín, C. Martínez-Gonzalez, P. de Lucas, I. Mir, G. Peces-Barba, N. Feu-Collado, I. Solanes, I. Alfageme, and C. Casanova, Chest 149 , 45 (2016). DOI:10.1378/chest.15-1055 [43] M. A. Dentener, E. C. Creutzberg, H. J. Pennings, G. T. Rijkers, E. Mercken, and E. F. Wouters, Respiration 76 , 275 (2008). DOI:10.1159/000117386