Probing Cellular and Molecular Mechanisms of Cigarette Smoke-Induced Immune Response in the Progression of Chronic Obstructive Pulmonary Disease Using Multiscale Network Modeling

Chronic obstructive pulmonary disease (COPD) is a chronic inflammatory disorder characterized by progressive destruction of lung tissues and airway obstruction. COPD is currently the third leading cause of death worldwide and there is no curative treatment available so far. Cigarette smoke (CS) is the major risk factor for COPD. Yet, only a relatively small percentage of smokers develop the disease, showing that disease susceptibility varies significantly among smokers. As smoking cessation can prevent the disease in some smokers, quitting smoking cannot halt the progression of COPD in others. Despite extensive research efforts, cellular and molecular mechanisms of COPD remain elusive. In particular, the disease susceptibility and smoking cessation effects are poorly understood. To address these issues in this work, we develop a multiscale network model that consists of nodes, which represent molecular mediators, immune cells and lung tissues, and edges describing the interactions between the nodes. Our model study identifies several positive feedback loops and network elements playing a determinant role in the CS-induced immune response and COPD progression. The results are in agreement with clinic and laboratory measurements, offering novel insight into the cellular and molecular mechanisms of COPD. The study in this work also provides a rationale for targeted therapy and personalized medicine for the disease in future.


Introduction
Chronic obstructive pulmonary disease (COPD) is characterized by airflow limitation caused by destruction of the lung parenchyma and/or airway obstruction [1][2][3]. COPD is currently the third leading cause of death worldwide and poses a major public health burden globally [4]. Moreover, COPD is associated with the development of lung cancer [5]. There is no cure available for COPD and current drugs are mainly effective in improving symptoms and exacerbations but generally do not slow down the progression of the disease [6]. Therefore, it is important to understand the cellular and molecular mechanisms of COPD for developing effective treatments of the disease.
COPD is a chronic inflammatory disease caused by inhalation of toxic particles and gases, mostly cigarette smoke (CS) [1][2][3]7]. Despite the fact that CS is the major risk factor for COPD, many chronic smokers maintain normal lung function (so-called resistant smokers) [2], so do some smokers even after more than 40 pack years of smoking [8], while only~20-30% of chronic smokers develop the disease [1,2,7,9]. This suggests that the susceptibility of smokers to COPD can vary significantly [1,2,8,9]. However, the cellular and molecular basis for the disease susceptibility remains to be elucidated albeit genetic or environmental factors may play a role [1,2]. As chronic cigarette smokers with normal lung function also have increased pulmonary inflammation, this inflammation seems to be magnified in COPD. Understanding of the amplification of inflammation is not yet complete [1]. Cigarette smoking cessation is considered currently as the most important intervention to reduce COPD progression [10]. While quitting smoking can prevent the COPD progression in some patients, who are referred as (reversibly) susceptible smokers, cigarette smoking cessation fails to slow or preclude the COPD progression in others (referred as severely susceptible smokers) [2,11]. The precise understanding of different effects of smoking cessation has not yet been fully achieved [1][2].
The CS-induced inflammatory response in COPD progression involving both innate and adaptive immunity [1,2] is mediated via a complex network that encompasses multiple immune cell types, molecular mediators and lung tissues. Several different types of immune cells and molecular mediators are found to accumulate in the lungs of patients with COPD [1-3, 5-7, 12]. Important immune cells include macrophages, neutrophils, dentritic cells, and T lymphocytes and molecular mediators include cytokines, chemokines, and protein proteases such as metalloproteases (MMPs). There exists an enormous amount of literature regarding these individual network elements. However, little is known about combined interactions between these elements or the associated pathways in the network. In particular, while COPD progression is a multistage and dynamic process, studies on the temporal sequence of inflammation in the disease are lacking [2]. It is not clear how immune cells and molecular mediators are dynamically linked and which of these elements are determinants in the disease progression. This is particularly important for identification of biomarkers in the disease [6,[13][14][15][16][17]. For example, the levels of proinflammatory cytokines, TNF-α and IL-1β are increased in the lungs of COPD patients and were suggested as potential targets [6]. However, inhibition of TNF-α or IL-1β has been proved to be unsuccessful in clinical trials of patients with COPD [6]. A pilot study on patients with COPD revealed no change in levels of inflammatory markers following inhibition of TNF-α [18], but the underlying reason remains to be elucidated [6,19]. Recent studies have shown that the levels of IL-6 as well as the associated proteins, C-reactive protein (CRP) and fibrinogen, are significantly increased in COPD patients compared to those in smokers with normal lung function and healthy non-smokers [16]. IL-6 is considered to be a potential biomarker for COPD, but the detailed mechanism of IL-6 action in the disease progression is not yet fully understood.
In this study, we develop a network model for CS-induced immune response in the lung to address the issues described above. Our model analysis identifies several positive feedback loops that play a determinant role in the CS-induced immune response and COPD progression, providing novel insight into the cellular and molecular mechanisms of the disease. IL-6, forming a positive feedback loop, IL-6!Th17!IL-17!TD!IL-6 [39]. Th17 cells also produce IL-21 for the differentiation of CD8 + T cells from naïve CD8 cytotoxic T lymphocytes (T0) [37,42]. While CD8 + T cells produce IFN-γ to enhance the M1 inflammatory activities, they also release granzyme B and perforins, causing apoptosis/necrosis of targeted cells and leading to TD further [37]. In addition, IL-6 can down-regulate the activation of Treg that secretes IL-10 to inhibit Th17 [43]. Consequently, a positive feedback loop, IL-6┫Treg!IL-10 ┫Th17!IL-17!TD ! IL-6, is formed.
As proposed earlier, the immune cells, cytokines and TD discussed above can be treated as network nodes. The aforementioned interactions between these nodes are then integrated into the network model shown in Fig 1 in the next section. The constructed network bears a multiple timescale feature. Cytokine regulation of cell function through signal transduction usually occurs on a sub-second timescale, for example, whereas cell production of cytokines takes minutes to hours [44]. The cytokine regulation activity can thus be considered to be at steady state in the equations that describes the slow timescale activities of the cells. In this way a positive or a negative input can be modeled using an increasing or decreasing Hill function, respectively [44,45]. The dynamics of the cytokines, the immune cells and TD can thus be described using a set of ordinary differential equations (ODEs). In this study, the system of ODEs is solved numerically using MATLAB (version R2013a Mathworks) with a variable order and multistep solver, ode15s. MATLAB is also used to plot the simulation data to generate the figures presented below.

Network Model
In the present network model shown in Fig 1, M1, DC, Th1, CD8 + T and Th17 cells along with their corresponding cytokines, TNF-α, IL-6, IL-12, IFN-γ, and IL-17 form multiple  Th2 and Treg cells with the associated cytokines,  IL-4, TGF-β, 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-β and IL-10 (Fig 1). These pathways eventually converge at the TD node that represents the tissue damage. Here, we focus on the immunologic aspects of COPD and the TD node is highly coarse-grained, involving neutrophil-induced tissue damage, epithelial and endothelial cell injury and extracellular matrix degradation etc. As discussed above, for example, M1 produces chemokines to recruit neutrophils (not shown in Fig  1) for the release of neutrophil elastase and ROS, contributing to TD. The representation of this process is included indirectly in the M1!TD motif shown in Fig 1. In addition, while M1 secretes MMPs to cause TD, TD also produces EFs to recruit monocytes for generating M1. This process is represented in the M1!TD!M1 positive feedback loop where MMPs and EFs are included indirectly (Fig 1). As TD is the major feature of CODP in the lung parenchyma, the dynamics of TD is used to characterize the progression of COPD in this work.
As mentioned above, the dynamics of the network elements can be described by a set of ODEs. Here, the ODEs involve the following 18 variables: M 1 , M 2 , D C , T 1 , T 2 , T 8 , T 17 , and T g represent the densities of M1, M2, DC, Th1, Th2, CD8 + T, 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 α , I γ , and I β denote concentrations of the cytokines, IL-4, IL-6, IL-10, IL-12, IL-17, IL-21,TNF-α, IFN-γ, and TGF-β. The variable, T D , is defined as a ratio (in terms of a percentage) of damaged tissue to the whole lung parenchymal tissue (normal tissue plus damaged tissue). This T D definition is similar to that of destructive index (DI) [46][47] and can be measured experimentally. While a T D value of 0-30% is considered normal, a T D value larger than 30% is associated with COPD [48].

Immune Cell Dynamics
The M1 population originates from a constant source of M0 [49] whose differentiation is stimulated by the external stimulus, CS [7]. The M0 differentiation process, which is inhibited by I 10 , is up-regulated by I γ and I α at maximal rate k 2 [50]. As shown in Fig 1, TD also increases the M1 population and this process is also down-regulated by I 10 . Equation for describing the population dynamics of M1 with a decay rate, d M1 , is then given by where and k 0 2 are rate constants), and the regulation of the M 0 !M 1 differentiation process is characterized by the Hill function whose coefficients are often chosen to be 2 to allow sufficient nonlinearity [45]. In Eq 1, S is cigarette smoking intensity, whose values are given below. M 2 is differentiated from M 0 and this process is up-regulated by I 4 or I 10 [26]. With a decay rate, d M2 , the M 2 dynamics can be described by The population of dendritic cells, D C , originates from a constant source of immature D 0 activated and up-regulated by the external stimulus, S and the damage signal from T D [7]. This process is inhibited by I 10 . With a decay rate coefficient, d D C , equation for the D C dynamics is given by where T cells derive from a constant source of immature T-cells (Th0 or T 0 ), which are recruited and activated by D C in combination with specific cytokines as mentioned previously. For T 1 , T 2 , T 8 , and T 17 , the associated T-cell differentiation processes are down-regulated by I 10 . With inclusion of a decay rate, equations describing the dynamics of these T-cell populations can be expressed by where T i represents T 1 , T 2 , or T 8 , k i = k 7 , k 8 , or k 9 ; I i = I 12 , I 4 , or I 21 , and k p i = k p 1 , k p 2 , or k p 8 , which are the maximum rates for the self-proliferation of T 1 , T 2 , and T 8 , respectively. The T 17 dynamics, which is up-regulated by both I 6 and I β and is down-regulated by I 10 , can be described by Equation describing the T g dynamics, which is up-regulated by I β and inhibited by I 6 , can be written as Cytokine Dynamics As shown in Fig 1, the cytokines are secreted from their associated immune cells. In particular, T D also contributes to the I 4 and I 6 productions, respectively, as discussed previously. The cytokine secretion processes are often down-regulated by I 10 . Accordingly, the population dynamics of the cytokines can be described by where C j is the density of the j-th cell, which produces I i with a maxima rate k i,j , k I i ;T D is the rate constant of the I i production by T D , f j (I 10 ) is the regulation function, i.e., the Hill function The Hill function, f D j (I 10 ), in Eq 7 describes the I 10 down-regulation of the process for I 4 or I 6 production by T D (Fig 1). Thus, for I 4 and I 6 , whereas f D (I 10 ) = 0 for other cytokines, whose population dynamics are governed by Equations A-G given in S1 File. For example, Eqs 10-11 describing the dynamics of I 6 and I 10 are given by In Eq 10, I 6 comes from M 1 [26] and T D [37], respectively. These production processes are down-regulated by IL-10 [50][51][52]. I 10 in Eq 11 is produced by M 2 [26] and T g [36], and can be inhibited by itself [49,51].

TD Dynamics
As discussed above, T D can be generated by S [7], M 1 [24][25][26] and T 8 [37], respectively. I γ and I 17 also cause TD [37]. These TD generation processes are all down-regulated by I 10 [52]. As mentioned earlier, in the lung parenchyma the process for MMP-induced TD generation is also inhibited by I β [29]. Equation describing the T D dynamics thus is where T D is defined by damaged tissue/(normal tissue+damaged tissue)×100% and tissue repair is involved indirectly in the last term in Eq 12.

Sensitivity Analysis
The parameter values in the above equations were taken or estimated from literature in the following discussion. For those whose experimental data are not available, we performed a global sensitivity analysis to obtain order-of-magnitude estimates. Sensitivity analysis is useful to determine parameters that play a critical role in affecting the model outcome [49]. In this study, a global sensitivity analysis outlined by Marino et al. [53] is applied to assess the sensitivity of the model outcome, T D , at t = 4000 day in the stable stage of COPD to variations of all parameters in the model. Baseline values of these parameters are listed in Table A in S1 File. A range of 10% to 200% of the baseline values was specified for each parameter in a way similar to that in ref. 49. All parameters are assumed to be uniformly distributed in their corresponding interval and 2000 samples are generated for each parameter using the Latin Hypercube Sampling method [53]. The partial rank correlation coefficient (PRCC) as well as p-value for each parameter is then calculated. The calculated PPRCs with the p-values smaller than 0.01 are presented in Figure A in S1 File. The values of PRCC range between -1 and +1 with the sign determining whether an increase in the parameter will increase (+) or decrease (-) the T D output. Our sensitivity analysis shows that a set of parameters including k 13 , k 14 , k 15 , and d TD have relatively large PRCCs (>0.1, Figure A in S1 File), suggesting that TD outcomes are sensitive to these parameters. This set of parameters, whose values can be changed from individual to individual, are associated with different network elements, positive feedback loops and pathways, which are critical in the CS-induced immune response and COPD progression. Therefore, investigation of model outcome changes with variations in these parameters can be of great use to disclose the cellular and molecular mechanisms of COPD as discussed below.

Simulations of Dynamics of Immune Responses to CS at Different Dose Levels
Eqs 1-12 are applied in the following discussion with the parameters given in Table A in S1 File. The time course of the simulations is marked in days and units of cell populations and cytokine concentrations are in terms of cells per milliliter (cells/ml) and pmol per liter (pmol/ L). The initial conditions of the variables for all activated cells and secreted cytokines are set to zero. CS-induced immune response and COPD progression are dependent upon the dose of CS inhalation [54]. Figs 2-4 respectively show the time courses of the changes in the immune cells, cytokines and tissue damage (T D ) in response to CS at a relatively high level of the cigarette smoking intensity, S (S = 1.67). Here, S is defined as the ratio of a CS dose to the minimal CS dose required to cause COPD with the parameters given in Table A in S1 File.
As shown in Fig 2(B), M 1 along with I α and I 6 [Fig 3(B)] rises to a peak at 12 days of CS exposure and then decays until day 60 because of the suppression of IL-10 secreted mainly by M2, showing an acute inflammatory response to the CS exposure. These results are qualitatively consistent with experiments in mice [47]) (see Figure B in S1 File for comparison with mice experiments). This period of time can be referred as step 1 in the progression of COPD as proposed by Agusti et al. [2][3]. After step 1, however, the inhibitory effect of IL-10 on M1 (as well as proinflammatory cytokines) is counteracted by the production of M1 up-regulated by TNF-α, TD, and INF-γ (Fig 1). It turns out that M 1 is increased (but slowly) up to day 180 [ Fig  2(B)]. This time period is referred as step 2 in the COPD progression [2,3]. During this period, D C along with I 12 (data not shown), I 6 and I β [Fig 3(B)], and I 21 (data not shown) increases slowly and gradually, resulting in the slow productions of T 1 , T 17 , and T 8 , respectively. Consequently, T D rises gradually during this period. After step 2, innate proinflammatory network elements start to go up quickly as M1 predominates over M2 again [Fig 2(A)]. In particular, I 6 is increased more rapidly than I α [Fig 3(A)], enhancing the production of T 17 and T 8 and allowing the adaptive immunity to play an increasingly important role in response to the CS exposure and in tissue damage (Fig 4). Notably, T 8 goes up more quickly than T 1 and T 17 [ Fig  4(A)], dominating in the late phase in the COPD progression. Eventually, T D together with the immune cells and cytokines goes to a steady state (stable COPD) as shown in Fig 4(B). The results for the above immune cells and cytokines at the steady state are listed in Table B in S1 File. Our results are in agreement with experiments [55][56][57][58][59][60][61][62][63].   The results for the cases of S = 1.67 and S = 0.7 indicate that during the transition from the innate to the adaptive immunity, when M 1 predominates over M 2 , the system would proceed to high-grade chronic inflammation and eventually toward COPD (Figs 2-4); while M 2 (Fig 5) [Treg (Fig 7A)] is predominant, the acute inflammation turns into the low-grade chronic inflammation (Fig 6), and COPD does not occur [Fig 7(B)].

Cessation of Cigarette Smoking
Cigarette smoking cessation is considered a most important intervention to reduce COPD progression [10]. The dynamics of T D is shown in Fig 8 as cigarette smoking cessation occurs after different days of CS exposure (S = 1.67) with the parameters in Table A in S1 File. Our simulations demonstrate that when early smoking cessation happens before 920 days of CS exposure, T D falls to the baseline and COPD is prevented (Fig 8). Smoking cessation starting after day 920 leads to the reduction of T D to some extent, but COPD and inflammatory responses persist (Figure C in S1 File). Our results are qualitatively consistent with experimental and clinical observations [64][65][66]. Analysis of positive feedback loops will elucidate further these effects of smoking cessation as discussed later.

Susceptibility to COPD
As mentioned above, although CS is the major risk factor for COPD, only about 20-30% of chronic smokers are susceptible to the disease, suggesting that susceptibility of smokers to COPD varies [1,7,9]. This implies that parameters governing the dynamics of T D as well as other important network elements in the model can change among smokers with different levels of COPD susceptibility. Importantly, effects of cigarette smoking cessation can be different among smokers with variable susceptibility [65]. While quitting smoking can prevent COPD Cellular and Molecular Mechanisms of CS-Induced Immune Response in COPD Using Network Modeling in some patients (i.e., reversible susceptible smokers), smoking cessation fails to slow or stop COPD progression in others (severely susceptible smokers) [2,66].
Our global sensitivity analysis shows that the T D outcome is sensitive to a subset of parameters including k 13 , k 14 , and k 15 and d TD in Eq 12.  Table A in S1 File are used. The results show that when k 13 is less than 2.6×10 −2 ml/(cell day),T D remains at a low-level (<30%), exhibiting a COPD resistant feature seen in Fig 9(A). While k 13 lies in a value range of 2.6×10 −2 and 0.31 ml/(cell day), COPD occurs but smoking cessation leads to T D decreasing to the baseline shown in Fig 9(B). In this case, COPD is reversible [2]. When k 13 is larger than 0.31 ml/(cell day), T D at the steady state is reduced to some extent but still remains at a high level (>30%) after smoking cessation at day 2500. A COPD patient in this case is severely susceptible [11]. Interestingly, for a severely susceptible smoker whose k 13 is relatively large, the M1-induced destruction of lung tissue predominates. In this case, M1 is sufficient for the progression of COPD, consistent with experiments in mice [24,25].
A well-recognized example for severely susceptibility is the severe deficiency of α-1 antitrypsin, which is the major inhibitor of neutrophil elastase. The severe α-1 antitrypsin deficiency is present in only 1-2% of individuals with COPD [67]. As neutrophil elastase functions to degrade macrophage elastase inhibitor of metalloproteinase-1, the severe α-1 antitrypsin deficiency significantly increases the tissue destruction activity of macrophage elastase [68]. Therefore, the overall effect of the severe α-1 antitrypsin deficiency is to significantly enhance the M1-induced T D generation [68], corresponding to the case where k 13 adopts a relatively large value for a severely susceptible smoker (Fig 9). Similar results were obtained by varying k 14 and k 15 shown in Figures D-E in S1 File, respectively. Intriguingly, the changes in some parameters such as k I 6 ;M 1 in Eq 10, which do not affect directly on the T D outcome, also lead to similar results for T D dynamics, shown in Figure F in S1 File. Here, our modeling results demonstrate that similar T D outcomes can result from variations of different parameters. As will be discussed later, these parameters may be associated with different mechanisms of CS-induced immune responses in COPD progression, implying that similar COPD phenotypes can be caused by different mechanisms (endotypes). This finding could be important for developing therapeutic methods to treat COPD, as also discussed below.

In Silico Knockout Simulations for Identification of Important Network Elements in COPD Progression
To identify important network elements for CS-induced COPD, in silico knockout simulations for the proinflammatory elements, e.g., M1, DC, Th1, Th17, CD8 + T cells and TNF-α, IL-6, IFN-γ, and IL-17 were conducted in the following discussion. Deletion of a network element in the model was performed by setting all parameters of the element and the rate to zero [69]. Here, the parameters used in the simulations are listed in Table A  As seen in Fig 10, M1 knockout results in significant reductions of T D and proinflammatory signals to rather low levels despite continuous CS exposure. M1 not only secretes proinflammatory cytokines such as TNF-α, IL-6 and IL-12, but also causes tissue damage as discussed above (Fig 1). Our results are in line with experiments demonstrating that M1 is a determinant for CS-induced COPD [24].
DC knockout leads to the low-level outcomes of T D and the adaptive immune signals such as I 17 shown in Fig 10. The innate immune signals, e.g. I α and I 6 , are reduced partially but still remain at a certain level. Upon deletion of Th1 there are no significant changes in the CS-induced TD and proinflammatory outcomes. This result is in agreement with clinical data for low Th1 population density [63], and mice experiments showing that CS induced significantly the production of Th17 rather than Th1 [70]. Deletion of Th17 or CD8 + T cells leads to a significant reduction of T D so that the progression of COPD is halted seen in Fig  10(A). However, the CS-induced proinflammatory signals such as I α and I 6 are reduced but still maintain at a relatively high level compared to those in silico experiment of M1 knockout. Different from the CD8 + T deletion, the knockout of Th-17 leads to a reduction of I 17 to the baseline, as shown in Fig 10(D). These results are in line with experiments in mice [70,71].
TNF-α knockout simulations show that there is no significant reduction in T D and other proinflammatory signals (Fig 10). This result is in line with clinical data and pharmaceutical studies discussed previously [6,18,19]. The IL-6 deletion allows the adaptive immune signals, T 17 , T 8 , I 21 (data not shown) and I 17 to maintain at the baseline level, and results in a significant reduction of T D (<10%) at the steady state (Fig 10), whereas innate proinflammatory signals such as I α maintains at a relatively high level (blue triangles in Fig 10). While the IL-17 knockout leads to the results similar to those in the deletions of Th17 and CD8 + T, the IFN-γ deletion results in T D (at the steady state) at a relatively high level (~32%, black asterisks in Fig 10).

Loop Breaking Simulations to Identify Important Positive Feedback Loops in COPD Progression
To investigate the importance of the positive feedback loops shown in Fig 1 in COPD progression, a feedback loop breaking approach is applied in the present study. Here, four aforementioned positive feedback loops, M1!TD!M1 (Loop 1), IL-6!Th17! IL-17!TD!IL-6 (Loop 2), M1!IL-12!Th1!IFN-γ!M1 (Loop 3), and IL-6 ┫Treg!IL-10 ┫Th17!IL-17!TD!IL-6 (Loop 4), are explored by setting, e.g., k 3 in Eq (A), k I6, TD in Eq (J), and k I12, M1 in S1 File to zero, leading to the breaking of Loops 1-3 on TD!M1, TD!IL-6, and M1!IL-12, respectively. We set K Tg, I6 in Eq 6 to a very large value, e.g., 10 4 to break Loop 4 on IL-6 ┫Treg in the calculation. Simulations were performed using our model with the parameters in Table A in S1 File (d TD = 2.9×10 −3 /day). The results reveal that breaking Loop 1, 2, or 4 allows T D to remain less than 10% (Fig 11), indicating that COPD does not occur. However, breaking Loop 3 has no profound effects on COPD progression (cyan line in Fig 11). This finding is consistent with virtual knockout of Th1 discussed above.
Loops 1, 2 and 4, which all include the TD node, can amplify the immune response to CS and enhance the TD outcome when they are activated. As Loop 1, 2 or 4 is necessary for COPD progression under the condition given in Table A in S1 File discussed above, we also carried out simulations to investigate whether one of these loops is sufficient for COPD when other positive feedback loops are broken. Here, k 3 = 1.9×10 6 cell/(ml day), k I6,TD = 22.0pmol/(cell day), and K Tg,I6 = 2.3 pmol/L are used for Loop 1, 2, and 4, respectively. The other parameter values are given in Table A in S1 File (d TD = 2.9×10 −3 /day). The results show that activation of Loop 1, 2, or 4 can lead to COPD [Fig 12(A)]. Again, Loop 3 alone is not able to drive COPD progression by modification of k I12, M1 .
It is of interest to note that Loops 1, 2 and 4 play important but different roles in the CSinduced immune response in COPD progression as shown in Fig 12. The activation of Loop 1 results in a marked predominance of M1 over M2, thus responsible for the innate immune response [Fig 12(B)]. In this case, T g is predominant over T 8 or T 17 [Fig 12 (C)]. The activation of Loop 4, as well as Loop 2 (data not shown), drives Th17 and CD8 + to predominate over Treg, responsible for the adaptive immune response [Fig 12(D)]. In this case, both M 1 and M 2 are relatively low [Fig 12(B)]. These results demonstrate that there exist different cellular and molecular mechanisms (endotypes) in COPD progression. In particular, these different endotypes (mechanisms) can lead to similar T D (CODP) outcomes, highlighting the heterogenous nature of COPD. This finding would be of particular importance in treatment of COPD using personalized medicine and target therapy [72] as discussed below.

Discussion
The molecular and cellular mechanisms of the CS-induced immune response in COPD progression are not completely understood. In particular, the issues regarding the dynamics of CSinduced immune response in COPD, the effects of cigarette smoking cessation, and the disease susceptibility remain to be elucidated [1][2]. As COPD is a chronic and progressive inflammatory disease whose dynamical time scale is usually very long (over 20 years of CS exposure [11]), it would be extremely difficult for real-time measurements in the clinic or the laboratory. The main objective of this paper is to use computer modeling and simulation to address these issues. As discussed above, the nodes in our network model, which bears a multiscale nature, represent the cytokines, immune cells and damaged tissues (TD) whose dynamics characterizes the progression of COPD. The interactions between these nodes are often highly nonlinear and can be described using the Hill functions. The population dynamics of the network elements are described by a set of DOEs with parameters, whose values are known or estimated from established literature. For those no experimental data are available, we performed the global sensitivity analysis to obtain order-of-magnitude estimates.
The results in this work demonstrate that CS-induced COPD development is a multi-step process involving both innate and adaptive immune responses. In the early acute phase of CS exposure, innate immune response predominates. During the transition from the innate to the adaptive immunity, if M 1 predominates over M 2 , the system proceeds to high-grade chronic inflammation and eventually toward COPD where the adaptive immunity play a dominant role (Figs 2-4). However, when M 2 (Treg) is predominant over M1 (Th17 and CD8 + T), the acute inflammation turns into the low-grade chronic inflammation, and COPD does not occur (Fig 5-7).
CS inhalation has been considered the major risk factor for COPD, but only about 20-30% chronic smokers develop the disease, suggesting that cigarette smokers have different levels of COPD susceptibility [1,9]. Our simulations disclose that there exist three types of smokers according to their COPD susceptibilities, i.e., resistant, reversely susceptible and severely susceptible smokers. While long-term CS inhalation can cause just low levels of chronic inflammation in resistant smokers but without COPD, susceptible smokers can develop COPD eventually under the same CS exposure conditions as shown in Fig 9. After cigarette smoking cessation, COPD can be prevented in reversible susceptible smokers, but the disease cannot be fully reversed in severely susceptible smokers (Fig 9). The sensitivity analysis in this work has identified a subset of parameters in the model that govern the dynamical behaviors of COPD and describe the different disease susceptibilities of different smokers as shown in Fig 9 and Figures D-F in S1 File.
The network model in this paper involves multiple immune cells, cytokines and lung tissues, forming multiple proinflammatory and anti-inflammatory/regulatory pathways with several positive feedback loops. The in silico knockout experiments performed in this study have identified several important proinflammatory elements including IL-6, IL-17 cytokines, and M1, DC, Th17, CD8 + T cells. It is intriguing to note that despite high concentrations of TNF-α often found in COPD smokers, our present study shows that this cytokine may not play a significant role in the progression of COPD. This finding is consistent with biopharmaceutic and clinical studies [6,18,19]. The feedback loop breaking simulations demonstrate that Loops 1, 2 and 4, which all involve the TD node, play important but different roles in the COPD progression. Activation of Loop 1, which enhances the M1 and TD productions, can promote the M1/ M2-type COPD while Loops 2 and 4 contribute to the (Th17+CD8 + T)/Treg-type disease where IL-6 and IL-17 are key molecules for the disease progression. These results indicate that COPD can be heterogeneous and can result from different molecular and cellular mechanisms.
The reason why inflammation persists in COPD patients even after long-term smoking cessation is currently unknown [2]. The above loop breaking results may provide novel insight into the aforementioned smoking cessation effects. After cigarette smoking cessation, the disease remains if these positive feedback loops continuously work to cause the tissue damage. But when TD is reduced due to smoking cessation to such an extent that these positive feedback loops are unable to drive tissue damage further, COPD will be suppressed.
The feedback loop breaking simulations are consistent with the in silico knockout experiments in this study, together identifying the network inflammatory determinants for CODP progression. For example, IL-6 is not only a product from Loop 1, but also a key component of both Loops 2 and 4. IL-6 is an important proinflammatory factor for synthesis of acute phase proteins such as C-reactive protein, which is associated with several acute and chronic inflammatory diseases including COPD [15,73,74]. This cytokine is also identified as a major regulator of the balance between Treg and Th17 (CD8 + T) cells as seen in our above discussions. As IL-6 plays an important role in COPD progression, it has recently been recognized as a potential target for COPD [73,74]. Another example is that IL-17 as well as Th-17 is identified as a key component for COPD. Targeting IL-17 and Th-17 has become a promising strategy for treatment of the disease [75]. However, COPD is a complex and heterogeneous disorder, e.g. similar clinical phenotypes can come from different endotypes, which are associated with different molecular and cellular mechanisms, as shown in our above modeling analysis. It is critical to identify the molecular and cellular disease mechanisms, by which subtypes of COPD patients are defined. As such, future treatment options would target the identified endotypes using available anti-inflammatory drugs [72]. For this purpose, specific biomarkers of these endotypes would be particularly useful. Our modeling study offers a possible approach to probe endotype biomarkers for COPD and provides novel insight into this personalized medicine strategy as discussed above.
As we focus on the CS-induced immune response in the progression towards stable COPD in this work, other pathogenic mechanisms including exaggerated proteolytic activity, resulting from an imbalance between protease and antiprotease, and excessive oxidative stress from an oxidant-antioxidant imbalance [7] are implicitly involved, for example, in the representation of the M1!TD process as discussed above. Disruption of the balance between cell death and repair is also included indirectly in the TD component. These highly coarse-grained representations can be extended in more detail. Indeed, such studies have been undertaken recently in references [76][77]. However, the adaptive immunity has not been investigated and the insight into the cellular and molecular mechanisms of COPD is incomplete in these studies [76,77].
Although cigarette smoking is the main risk factor for COPD, other factors can also play an important role in disease development and progression [78]. For example, as almost every component of the immune system undergoes age-associated changes, aging is a risk factor for developing COPD [78]. An exacerbation of COPD (ECOPD) is defined as an acute event characterized by a worsening of the patient's respiratory symptoms that is beyond normal day-today variations and leads to a change in medication, responsible for substantial COPD mortality [79]. Viral or bacterial infections are the main causes of ECOPD, leading to an acute flare-up of inflammation in the lung with stable COPD. Interestingly, mechanical forces of lung tissues have been shown to contribute to COPD progression by allowing rupture of tissue elements which directly leads to increased airspaces and this progression can go on even after smoking cessation [80]. It is important to point out that our current network model focusing on CSinduced COPD does not take these factors into account. Elucidating the roles of these factors in COPD progression is beyond the scope of this article. Possible extensions of the present work to incorporate these factors into the network model will be investigated in future studies.
It is of great interest to note that COPD and lung cancer are frequently induced by cigarette smoking, but these two disorders show opposite phenotypes [81,82]. While COPD is featured by excessive lung injury and airway epithelial cell death, lung cancer is caused by unregulated proliferation of epithelial cells [82]. Numerous epidemiological studies have linked the presence of COPD with increased lung cancer incidence, however, the molecular and cellular links between these two diseases remain obscure [5,[82][83][84]. Our network modeling research along this line remains for future.

Conclusion
In this work, we develop a network model for the dynamics of CS-induced immune response in COPD progression. Using this model, computer simulations are performed for the investigations of smoking cessation effects and susceptibility of smokers to COPD, and the identification of important network elements in COPD progression. Our modeling study identifies several positive feedback loops that play important but different roles in COPD progression. The computational results in this study are consistent with laboratory and clinical observations, providing novel insight into the cellular and molecular mechanisms of CS-induced COPD.