text:
Mathematical modeling identifies Lck as a potential mediator for PD-1 induced inhibition of early TCR signaling

Programmed cell death-1 (PD-1) is an inhibitory immune checkpoint receptor that negatively regulates the functioning of T cell. Although the direct targets of PD-1 were not identified, its inhibitory action on the TCR signaling pathway was known much earlier. Recent experiments suggest that the PD-1 inhibits the TCR and CD28 signaling pathways at a very early stage ─ at the level of phosphorylation of the cytoplasmic domain of TCR and CD28 receptors. Here, we develop a mathematical model to investigate the influence of inhibitory effect of PD-1 on the activation of early TCR and CD28 signaling molecules. Proposed model recaptures several quantitative experimental observations of PD-1 mediated inhibition. Model simulations show that PD-1 imposes a net inhibitory effect on the Lck kinase. Further, the inhibitory effect of PD-1 on the activation of TCR signaling molecules such as Zap70 and SLP76 is significantly enhanced by the PD-1 mediated inhibition of Lck. These results suggest a critical role for Lck as a mediator for PD-1 induced inhibition of TCR signaling network. Multi parametric sensitivity analysis explores the effect of parameter uncertainty on model simulations.

Introduction
Activation and subsequent proliferation of T cell are crucial events preceding pathogen clearance. However, proper functioning of the immune system also relies on the ability of T cells to promote self-tolerance. Hence, these processes are tightly controlled at multiple levels by regulatory mechanisms [1]. T cells have co-stimulatory and co-inhibitory receptors that coordinate to modulate its response [2]. TCR (T cell receptor) activation is primarily responsible for the activation of effector functions of T cells and its full activation needs co-stimulation by CD28 (Cluster of Differentiation 28) receptor [3,4]. Induction of TCR and CD28 signaling pathways result in T cell proliferation, increased glucose uptake and production of cytokines [5]. On the other hand, inhibitory receptors CTLA-4 (Cytotoxic T-lymphocyte-associated antigen 4) and PD-1 (Programmed Cell Death-1) negatively regulate the T cell response. Activation of PD-1 receptor has been shown to negatively affect several processes upregulated by the TCR and its associated co-stimulatory signaling pathways [6,7].Knockouts of the genes encoding these a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 inhibitory receptors have produced autoimmune phenotypes in the animal models suggesting their role in preventing autoimmune diseases [8][9][10].
The finding that cancer cells can be recognized and destroyed by the immune system, has established the field of cancer immunology and the interaction between cancer cells and immune system is being studied extensively [11,12]. Cancer cells are found to evade the immune system by employing numerous mechanisms and one such mechanism is the activation of negative regulators, PD-1 and CTLA-4 [13]. High expressions of ligands that are specific to the negative regulatory receptors have been detected on the cancer and immune cells in the tumor microenvironment [14,15]. Further IFN-γ produced by the T cell induces the expression of these inhibitory ligands on the cells of the tumor microenvironment [16][17][18]. Consequently, T cells receiving high level of inhibitory signals become inactive and have suppressed effector functions. PD-1 and CTLA-4 are extensively being studied and are considered as potential targets for activating the tumor infiltrating T cells that remain inactive in the immunosuppressive tumor microenvironment [19,20]. Antibodies against these receptors have shown exceptional efficacy and are considered as promising drugs that could potentially revolutionize cancer treatment. A few of the antibodies for instance, Nivolumab and Pembrolizumab targeting PD-1 receptor have been approved by the FDA (Food and Drug administration) for the treatment of melanoma [2]. However, administration of these immune checkpoint inhibitor drugs has numerous adverse effects and the treatment remains ineffective for a significant proportion of patients [21]. Apart from its role in inducing tumor immune escape, its role in several viral infections such as HIV (Human immunodeficiency virus), HCV (Hepatitis C virus) and HBV (Hepatitis B virus) are also demonstrated [22]. Exhaustion of T cells due to persistent TCR stimulation is observed during chronic viral infections [23,24]. Hence, an understanding of how the PD-1 receptor influences the T cell response is crucial for the development of effective treatment against cancer, autoimmunity and several other diseases.
Mathematical models have been an integral part in understanding complex biological phenomena such as apoptosis [25], cell cycle [26,27], NF-κB oscillations [28], cellular differentiation [29], cell signaling [30]. Mathematical modelling tools have become popular in explaining various aspects of immune systems [31] such as discrimination of self and non-self antigen [32,33], T cell activation [34][35][36], cytokine signaling pathways [37][38][39], T cell differentiation [40]. With the accumulation of quantitative and semi quantitative experimental results, modeling the TCR signaling network is increasingly being explored [41]. Protein-protein docking, molecular dynamics and mathematical modeling studies on interaction of PD-1 with its ligands have provided insights into the atomistic details and factors affecting these interactions at cellular interfaces [42,43]. Mathematical models were developed to predict tumor response to immune checkpoint inhibitors, and immune checkpoint therapy in combination with radiotherapy [44,45]. However, until now there is no mathematical model available to understand the influence of PD-1 on TCR signaling molecules.
To better understand the molecular basis of PD-1 induced inactivation of T cell signaling molecules, we constructed a deterministic mathematical model of PD-1 regulatory pathway. The proposed model investigates the role of feedback regulatory mechanisms in inhibition of early TCR signaling by PD-1. The model simulations provides the mechanistic basis of several features of PD-1 mediated inhibition that were observed in the recent experiments done using the reconstitution system by Hui et al [46]. phosphorylation of cytoplasmic domain of the TCR complex at the ITAMs (Immunoreceptor tyrosine-based activation motifs) by the Src family tyrosine kinase, Lck (Lymphocyte specific protein tyrosine kinase). Phosphorylated CD3 ITAMs recruit Zap70 (Zeta-chain-associated protein of 70 KDa), a Syk family kinase. Studies on ZAP70 phosphorylation have shown that Y315 and Y319 are initially phosphorylated by Lck and this in turn allows the Y493 present in the activation loop of the Zap70 to get phosphorylated. Phosphorylation of Y493 has shown to be crucial for the complete activation of the kinase Zap70 [47,48]. Activated Zap70 phosphorylates TCR signaling molecules such as the adaptor molecule, LAT (Linker for activation of T cells) and Slp76 (SH2-domain-containing leukocyte protein of 76 KDa) [49]. Phosphorylated LAT recruits signaling proteins such as Gads (Grb2-related adapter protein 2), Grb2 which in turn interact with other signaling proteins. Gads associates with Slp76 and this brings Slp76 in close proximity to activated Zap70 bound to the cytoplasmic domain of TCR complex. Several other molecules are recruited to form a complex called as signalosome and the signal is transmitted to the downstream effector molecules. Signal is transduced to the nucleus by the activation of the transcription factors AP1, NFAT and NF-KB resulting in changes of the gene expression [50][51][52][53]. CD28 receptor is activated upon binding to its ligands B7-1 and B7-2 inducing a conformational change in its cytosolic domain. This recruits PI3K and it transmits signal by activating other signaling molecules [54].
PD-1 receptor is stimulated by binding to its ligand, PD-L1 or PD-L2 expressed on the surface of APCs (Antigen Presenting Cells) [55]. Studies have demonstrated the inhibitory effect of PD-1 stimulation on several of the TCR signaling components [56,57]. PD-1 receptor upon activation recruits the cytosolic tyrosine phosphatase Shp2, which in turn dephosphorylates the TCR and CD28 signaling components. Recently, a FRET (Fluorescence Resonance Energy Transfer) based assay was done by Hui et al. [46] using a biochemical reconstitution system, to understand the molecular basis of suppression of T cell functions by PD-1. Few of the components of the T cell signaling network were reconstituted to understand the effect of PD-1 activation on the sequence of events that occur after the phosphorylation of cytosolic domain of the receptors TCR and CD28. In the in vitro reconstitution system, the cytosolic domains of receptors and membrane bound proteins were tethered to the surface of LUVs (large unilamellar vesicles) and the cytosolic signaling components were added to the solution. FRET based techniques were used to quantify the phosphorylation state of the receptors involved. These experiments demonstrated that Shp2 recruited by the cytosolic domain of PD-1 receptor can directly dephosphorylate the cytosolic tails of CD3z and CD28. The key findings of the experiments are that the CD28 is preferentially dephosphorylated in response to increasing PD-1 concentration whereas the CD3z chain and the other components of the TCR signaling tested are less sensitive to the inhibitory effect of PD-1. Further these findings were also confirmed using a cell based experiment employing Jurkat T cells and Raji B cells.

Model construction
We have constructed an ordinary differential equation (ODE) based mathematical model for the early signaling events following the TCR, CD28 and PD-1 stimulation. The full network diagram of the model is given in Fig 1 and a simplified influence diagram of the model can be found in Figure A in S1 File.
Ligated PD-1 receptor is doubly phosphorylated in distributive manner by the active Lck kinase (Lck active ). Shp2 binds to both the singly and doubly phosphorylated PD-1 to form the complexes referred as CP 1 and CP 2 respectively. Shp2 bound to the phosphorylated PD-1 (CP 2 and CP 1 ) distributively self-dephosphorylate the PD-1 to return either CP 1 or mono/ unphosphorylated PD-1 and thereby releasing Shp2 from the complex. CD28 and CD3z receptors are phosphorylated by the active Lck kinase and dephosphorylated by the Shp2 bound to the phosphorylated PD-1 (CP 1 and CP 2 ). We must mention here that free Shp2 does not have any catalytic activity. We used Michaelis-Menten kinetics to model these phosphorylation and dephosphorylation reactions.
Phosphorylation of CD28 and CD3z lead to engagement of free PI3K and free Zap70 to form PI3K b andZap70 i complexes, respectively. Zap70 bound to CD3z (Zap70 i ) is phosphorylated by active Lck at tyrosine Y315 initially to form Zap70 a1 and then at tyrosine Y493 to form Zap70 a2 . Zap70 a2 is considered as the active form of Zap70 and it phosphorylates LAT, referred as LAT i , to form phosphorylated LAT, referred as LAT a . Although LAT i is phosphorylated at several sites by active Zap70, the individual phosphorylation events of LAT molecule are not considered separately in the model to reduce the complexity and therefore these phosphorylations are considered as processive in nature. LAT a is considered to be the active form of LAT and it binds to the adaptor protein Gads to form the complex Gads a . Gads a in turn interacts  Table 1 and Table 2. The description of model variables are listed in Table A  with free Slp76 to form LAT a -Gads-Slp76 complex referred as Slp76 i . Slp76 bound to the LAT a -Gads complex is phosphorylated by Zap70 a2 to form the complex Slp76 a .
Experimental and molecular dynamics studies have shown that phosphorylation of Lck at Y394 has an activating effect and whereas Y505 is inhibitory [66]. Further it is believed that when the Lck is phosphorylated at Y505 first, it acquires a closed conformation and subsequent phosphorylation at Y394 does not change its activity [67]. Therefore five different forms of Lck are considered in the model depending on the phosphorylation state of the two tyrosine amino acid residues (Y394 and Y505): a. unphosphorylated at both tyrosines (Lck i ), b. phosphorylated only at Y394 (Lck ya ), c. phosphorylated only at Y505 (Lck yi ), d.Y394 and Y505 phosphorylated sequentially (Lck yiya ) and e. Y505 and Y394 phosphorylated consecutively (Lck pi ). Consequently, among all forms of Lck in the model, only Lck ya and Lck yiya were considered to be active and are capable of phosphorylating their substrates. All these reactions are modeled using mass action rate law of chemical reactions.
Experimental results by Hui et al [46] on the recruitment of Shp2 at various doses of PD-1 demonstrated that Shp2 recruitment to the PD-1 saturates beyond a certain concentration of PD-1. As the mechanistic basis to this observation is not known, a correction factor was introduced in the equations for the phosphorylation of PD-1 to account for this observation. This correction factor ensures the saturation of Shp2 recruitment at increasing concentration of PD-1, depending on the concentration of the Lck in the system.

Model parameterization and simulation
The model consists of 20 ODEs and 36 parameters. The model equations, the parameters values and the biological description of variables of the model are listed in the Table 1, Table 2 and Table A in S1 File respectively. Among 36 parameters we managed to get values for 13 parameters from experimental literature. Rest of the rate constants were parameterized to reproduce the experimental data of Hui et al [46]. Shp2 may bind to doubly phosphorylated PD-1 with higher affinity as compared to the singly phosphorylated PD-1. However due to lack of experimental binding data and also to reduce number of parameters we preferred same binding parameters for both the phosphorylated forms. CP 1 and CP 2 are assumed to have same catalytic activity in dephosphorylation of all of its substrates. For doubly phosphorylated Lck, dephosphorylation of Y394 is faster than that of Y505. Such preferential dephosphorylation of activating tyrosine by SHP-1 has been observed earlier [68]. CP 1 and CP 2 are also considered to be equally stable. However, rate constants for the auto-phosphorylation of different forms of Lck were assumed to be different. Lck auto-phosphorylation rate constants have been estimated and employed in a published mathematical model by Rohrs et al [62].  [46]. Simulations for a few experiments for which the concentrations are not mentioned, are done by using guessed concentrations. However, changing the concentration in most cases did not alter the qualitative behavior of simulated responses. Initial concentration for the following Lck states-Lck yiya , Lck ya , Lck yi and Lck i are taken as 25% of the total Lck concentration. This is consistent with the relative proportion of different forms of Lck measured in Jurkat T cells [67]. The initial condition of Lck pi was chosen to be 0 as it is believed that majority of doubly phosphorylated Lck is derived by the phosphorylation of Lck ya at Y505 [67].
Units of concentration and time in the model are nM and seconds, respectively. Hence, initial concentrations of species were provided in nM. In few cases, nM concentrations in the simulation results were converted to molecules/μm 2 for easier quantitative comparison with published experimental results. Michaelis-Menten constant of phosphorylation of CD3z by Lck obtained from the literature had unit of μm -2 and was converted to nM before employing in the model. To interconvert the concentration units between nM and molecules/μm 2 we used 1 nM = 2.9 molecules/μm 2 as conversion factor. Concentration as surface density (molecules/μm 2 ) is estimated by dividing the number of molecules of a species (calculated from concentration in nM) by the total area of exposed vesicle membrane [63]. The volume, Avogadro's number, fraction of exposed lipid and the area of each lipid head. For a liposome of diameter 200 nm,~52.6% of total lipids were found to present on the outer membrane of liposome [63]. Thus with a lipid concentration of 1 mM (as used by Hui et al. [46]), the protein concentration of 1 nM becomes equivalent to surface density of~2.9 molecules/μm 2 .

Model validation
We benchmarked our model by performing several simulations corresponding to the experiments done by Hui et al on the biochemical reconstitution system. It includes the time course of recruitment of Shp2, Zap70 and PI3k to PD-1, CD3z and CD28 respectively; comparison of phosphorylation activity of Lck on CD3z and CD28; dissociation of Zap70 and PI3k from their corresponding receptors in the absence of Lck activity; the dephosphorylation of TCR signaling molecules and CD28 due to PD-1.
We investigated the recruitment of Shp2 to the PD-1 receptor upon its phosphorylation by Lck. Upon introduction of Lck at time zero, Shp2 rapidly binds to PD-1 and the extent of binding increases with PD-1 concentration (Fig 2A). Here we did not allow Shp2 to dephosphorylate the receptor after binding in order to recapture the experiments consisting of only binding domain of Shp-2 lacking catalytic activity. In order to determine the effect of self-dephosphorylation of PD-1 receptor bound to Shp2 on the kinetics of Shp2 recruitment, we allowed selfdephosphorylations of CP 1 and CP 2 complexes. Results from Fig 2B indicate that after a rapid initial engagement, Shp2 slowly dissociates from the PD-1 due to the dephosphorylations. Further there is a weak dependence of Lck on the dissociation kinetics of Shp2 as the slope of the line increases slightly with the increased Lck concentrations.
Next we investigated the kinetics of Lck phosphorylation and PD-1 mediated dephosphorylation of receptors in absence of each other. We determined the kinetics of engagement of PI3K and Zap70 to their respective receptors upon Lck mediated phosphorylation of CD28 and CD3z in absence of PD-1 (Fig 2C). Similar to the experimental profile, simulation results show that phosphorylation of CD3z and subsequent recruitment of Zap70 proceeds at a much higher rate as compared to the CD28 phosphorylation and subsequent PI3K recruitment. Further we investigated the PD-1 receptor mediated dephosphorylation dynamics of CD28 and CD3z by calculating the percentage of free PI3K and Zap70 in absence of Lck (Fig 2D). Here we started simulations with a certain initial concentrations of Zap70 i and PI3K b considering that all of the CD28 and CD3z are already phosphorylated leading to full engagement of PI3K and Zap70 respectively. The initial concentrations of Lck and PD-1 were also set to zero. Then to determine the dephosphorylation kinetics we used a non-zero concentration of CP 2 complex as a proxy for membrane bound Shp2 used in the experimental protocol. In order to be consistent with the experimental protocol we did not allow self-dephosphorylation, dissociation and the degradation of CP 2 complex. In experiments Hui et al used full length Shp2 protein and it was directly attached to the membrane of vesicles. As observed in the experiments, simulations results showed that the percentage of free PI3K and Zap70 increases monotonically with time with similar rates.
Till now we showed the kinetics of phosphorylations or dephosphorylations in absence of either PD-1 or Lck. Now we investigated the kinetics when both of them are present essentially simulating the whole network. Simulation results showed that in presence of PD-1 the steady state values of free PI3K and Zap70 are higher than in absence of PD-1 ( Figure B in S1 File). This is due to dephosphorylations of CD28 and CD3z resulting release of respecting signaling molecules. The kinetics of PI3K indicates that the effect of PD-1 on CD28 is drastic as compared to CD3z (Fig 2E and 2F, Figure B in S1 File). We must point that all the simulation results of Fig 2 are in excellent qualitative and quantitative agreement with the experimental plots reported in Hui et al [46].
We now turn our focus towards the steady state results of model simulations and their comparison with reported experimental literature. We have carried out PD-1 dose response simulations of the full model where at a given concentration of PD-1 we record response at 30 min (Fig 3) as done in the experiments. Concentration of phosphorylated forms of each species was calculated and normalized with respect to the phosphorylation in the absence of PD-1. Phosphorylated CD3z was calculated by adding the CD3 a , Zap70 i , Zap70 a1 and Zap70 a2 in the model. Similarly, phosphorylated CD28 includes CD28 a and PI3K b ; phosphorylated LAT  [46]. Concentration of phosphorylated species is calculated as: Phosphorylated CD3z = CD3 a +Zap70 i +Zap70 a1 +Zap70 a2 , Phosphorylated CD28 = CD28 a +PI3K b , ZAP70 phosphorylated at Y315 = Zap70 a1 +Zap70 a2 , Zap70 phosphorylated at Y493 = Zap70 a2 , Lck phosphorylated at Y505 = Lck yi +Lck yiya +Lck pi , Lck phosphorylated at Y394 = Lck ya +Lck yiya +Lck pi , Phosphorylated LAT = LAT a + Gads a +Slp76 i +Slp76 a and Phosphorylated SLP76 = Slp76 a . Consistent with the experimental results PD-1 leads to dramatic decrease in CD28 phosphorylation. Whereas in the same concentration range of PD-1, dephosphorylation of CD3z is not significant (Fig 3). We have calculated the IC 50 and Hill coefficient of the response by fitting the dose response curves with a standard Hill function (k � x n =ðIC n 50 þ x n Þ where k is the maximum response). See Figure C in S1 File for fitting of dose responses with Hill function. The IC 50 value of CD28 is smaller in several orders of magnitude than that of CD3z highlighting the sensitivity of CD28 towards PD-1 (

Insights and predictions
In Figs 2 and 3 we have shown results from the model that quantitatively recaptures the experimental observations reported in Hui et al and thus validating the model. Hui et al suggested that the net dissociation of Shp2 from PD-1 complexes seen in the experiments, is due to the dephosphorylation of PD-1 by the PD-1-Shp2 complex. To test this, we performed simulations in the presence and absence of self dephosphorylation activity of CP 1 and CP 2 complexes for a large range of PD-1 and Lck concentrations. We find that till 500 nM of Lck there is a net decrease in Shp2 binding with increase in PD-1 concentration (Fig 4) as compared to simulations lacking self-dephosphorylation of PD-1 receptor by Shp2.
We generated similar surface plots to determine the dependence of phosphorylation of CD28 and CD3z and the subsequent recruitment of PI3K and Zap70 by the respective phosphorylated receptors in absence of PD-1 ( Figure D in S1 File). We find CD28 and CD3z are completely phosphorylated irrespective of the substrate concentration for Lck concentrations above 10 nM due to the large catalytic activity of Lck. This also suggests that Lck can completely phosphorylate all the CD28 and CD3z at physiological concentration of these components. Consequently PI3K recruitment and Zap70 recruitment are almost independent of Lck concentration. In the case of Zap70, 100% recruitment is achieved at approximately 500nM of CD3z concentration. Percentage of active Zap70 increases with increasing CD3z and Lck concentrations thus showing stronger dependence on CD3z and Lck. As opposed to this, complete recruitment and activation of Slp76 is achieved at low Lck and CD3z concentrations ( Figure D in S1 File). Now to assess the effect of PD-1, on the same variables we ran simulations with 300 nM PD-1 and 300nM Shp2 (Fig 5). Introduction of PD-1 leads to a dramatic reduction in CD28 phosphorylation and subsequent reduction in PI3K recruitment as compared to their response without PD-1 ( Figure D in S1 File). Whereas CD3z phosphorylation is reduced only at low Lck concentration highlighting the differential sensitivity of CD28 and CD3z towards PD-1. Lck at higher concentrations could reduce the inhibitory effect of PD-1 on the phosphorylation of these receptors. Due to the less sensitive nature of CD3z to PD-1 the effect on association of downstream signaling molecule Zap70 is minimal (Fig 5D).
In order to find out the effect of PD-1 for a range of its concentration, we performed a two dimensional scan of PD-1 and either CD28 or CD3z as appropriate (Fig 6) with 300 nM Shp2 and 100 nM Lck. PD-1 fully reverses the effect of Lck on CD28 even at its moderate concentration (~300 nM) resulting a complete loss of engaged PI3K (Fig 6A & 6C). Inhibition of CD3z phosphorylation happens at PD-1 concentrations greater than 500nM only when the CD3z concentration is also higher (Fig 6B). The maximal inhibition happens only at very high concentrations of PD-1 and CD3z. Recruitment of Zap70, its subsequent activation (Zap70 a2 ) and activation of Slp76 are affected only at very high concentration of PD-1 (Fig 6D-6F). The effect of PD-1 on Slp76 is very weak as Slp76 is not a direct target of PD-1. Thus PD-1 causes a dual effect by a strong inhibition of CD28 pathway and a moderate inhibition of CD3z.
PD-1 dose response simulations (and experiments) showed that PD-1 dephosphorylates both the activating (Y394) and inhibitory (Y505) tyrosine amino acid residues of Lck (Fig 3). Although, both the tyrosine amino acids are dephosphorylated to similar extent by PD-1 (similar profiles in the PD-1 dose response curve in Fig 3), due to the presence of a doubly phosphorylated species of Lck which retains catalytic activity, it is unclear if the dephosphorylation of Lck has a net activating or inhibitory effect. In the model Lck yiya and Lck ya are considered as active Lck forms capable of phosphorylating its substrates PD-1, ZAP70, CD28 and CD3z, and all other forms of Lck are considered as inactive. In the absence of PD-1, both the active and inactive Lck forms are present in equal proportion but in the presence of 300nM PD-1, proportion of inactive Lck increases over time (Fig 7A). A two dimensional variation of PD-1 and Mathematical modeling of PD-1 inhibitory pathway on early TCR signaling Lck shows (Fig 7B) that above 100 nM Lck concentration the proportion of active Lck decreases systematically with increasing PD-1 indicating its net inhibitory effect on Lck. At Lck concentrations below 100nM, the decrease in active Lck proportion with increasing PD-1 concentration is less and this is due to the requirement of Lck for PD-1 activation. Thus the model predicts that PD-1 imposes a net inhibitory effect on Lck and thereby causing an indirect inhibition to CD28 and CD3z signaling.
These observations raise the question, how significant is the inactivation of Lck in inhibiting the TCR and CD28 early signaling molecules. To explore this, we performed simulations where CP 1 and CP 2 (PD-1-Shp2 complexes) do not dephosphorylate Lck but dephosphorylate its usual target CD28 and CD3z. In the absence of Lck dephosphorylation by PD-1/Shp2 complexes, the inhibitory effect on Zap70 and Slp76 molecules are shown in Fig 7C and 7D. Although activation of CD3z and CD28 receptors were not greatly affected (not shown) but the activation of Zap70 and Slp76 are significantly affected without Lck dephosphorylation. Hence, the inhibitory effect of PD-1 on the downstream components such as the Zap70 and Slp76 are achieved indirectly via Lck dephosphorylation.
We have shown that after rapid engagement Shp2 dissociates from the PD-1 at later time ( Fig 2B). To test if self-dephosphorylation is solely responsible for this net dissociation or loss of Lck activity due to Lck dephosphorylation affects the net dissociation, we again carried out perturbation simulations of the model. As shown before without self-dephosphorylation and with Lck dephosphorylation there is no recovery of Shp2 (Fig 8B). However in the reverse  Mathematical modeling of PD-1 inhibitory pathway on early TCR signaling condition also the model predicts no net dissociation of Shp2 from PD-1 (Fig 8A). These simulations show that both the self dephosphorylation and Lck dephosphorylation are collectively responsible for the observed net dissociation of Shp2 from PD-1. Previous simulations in Fig 7 show that the Lck is inactivated by this PD-1 induced dephosphorylation. Hence, the self dephosphorylation activity of CP 1 and CP 2 complexes, decrease their stability and the inactivation of Lck decreases the phosphorylation of PD-1 thus inhibiting the formation of these inhibitory complexes.

Sensitivity analysis
Predictions and output of a mathematical model rely largely on parameter values. Hence, it is crucial to test the robustness of the model output to changes in parameters. Parameter sensitivity analysis is a technique extensively used for mathematical models to determine the fluctuations in the model output due to uncertainty in the parameters used. We used global parameter sensitivity analysis where all the chosen parameters are varied simultaneously to explore their collective effect on the dynamics of the system (see Materials and methods). Global sensitivity analysis offers advantage over the local parameter sensitivity analysis, as it explores the impact of the interaction between uncertainties in different parameter values and is considered more reliable [69]. Owing to the large number of optimized parameters, global sensitivity analysis was performed to test the robustness of Shp2 recruitment, PI3K recruitment and Slp76 activation in a 'modular' manner. 1) Shp2 recruitment by PD-1: Dissimilarity measure, K-S statistic (see Materials and Methods) was calculated for parameters involved in the activation of PD-1 by phosphorylation, formation of CP 1 and CP 2 complexes and parameters for Lck auto-phosphorylation and dephosphorylation. This module includes 15 parameters (all other parameters were kept constant) and representative cumulative frequency distributions for most and least sensitive parameters were provided in Figure E in S1 File. The K-S statistic estimates for all of the 15 parameters are summarized in the bar plot given in Fig 9A. The results suggest that the percentage of Shp2 recruitment is influenced mainly by changes in the parameters k p,pd1 and k a,shp which are the phosphorylation rate constant of PD-1 and association rate constant of Shp2 to phosphorylated PD-1. Their K-S statistic values are approximately 0.13. Shp2 recruitment is relatively robust to changes in the phosphorylation and dephosphorylation rate constants of Lck.
2) PI3K recruitment by phosphorylated CD28: Sensitivities of PI3K recruitment to 21 different parameters were calculated keeping other parameters constant. Out of all the parameters tested, the parameter k a,pi3k , rate constant for association of PI3K to phosphorylated CD28 was the most influential, with a K-S statistic value Mathematical modeling of PD-1 inhibitory pathway on early TCR signaling of~0.3, in affecting PI3K recruitment. Parameter k d,pi3k , dissociation rate constant of PI3K from CD28, had a K-S statistic value of 0.1. However, the K-S statistic for other parameters were less than 0.05, suggesting that PI3K recruitment is insensitive to a wider range of these parameter values around their optimized values (Fig 9B).
3) Slp76 activation: Effect of uncertainties in 29 parameters on Slp76 activation was tested. The remaining 7 parameters involving CD28-PI3K module were kept constant. Among all the 29 parameters, K-S statistic for the parameters k p1,zap , k p2,zap , k p,lat and k p,slp are high (K-S statistic values are approximately 0.35, 0.2, 0.15 and 0.15 respectively), highlighting their relative importance in determining the Slp76 activation ( Fig 9C). These parameters are the rate constants of phosphorylation of Zap70 by Lck, LAT and Slp76 by activated Zap70. Rest of the parameters had a K-S statistic value less than 0.05.

Conclusion
We propose a mathematical model of PD-1 pathway that negatively regulates the TCR activation pathway. The model has been parametrized and benchmarked to quantitatively recapture kinetic and dose response experimental results from recently published paper by Hui et al [46]. Further it explored two dimensional dose responses of various signaling molecules for a wide concentration range of PD-1 and CD28 or CD3z. The model provides molecular insights into the inhibitory effect of PD-1 on several key regulators such as CD28, CD3z, PI3K, Zap70 and Slp76. A key finding of the model is that PD-1-Shp2 complex targets TCR pathway both directly and indirectly. On the direct path, it leads to dephosphorylation of CD28 and CD3z resulting a decrease in binding of PI3K and Zap70, respectively. On the indirect path, it dephosphorylates Lck leading to a net inhibitory effect on Lck and thereby it indirectly downregulates activation of Zap70, LAT, Slp76 in the TCR-CD3 pathway whose activation requires phosphorylated Lck. Therefore PD-1 causes a dual effect to the TCR activation pathway by downregulating CD28, CD3z directly and Zap70, LAT, Slp76 indirectly via Lck dephosphorylation. The model highlights the importance of Lck dephosphorylation by PD-1-Shp2 complex in downregulating the TCR pathway. Global parameter sensitivity analysis of the model finds crucial parameters in PD-1 mediated dephosphorylation of several key molecules. Higher parameter sensitivity of Lck mediated phosphorylation rate of Zap70 in the Slp76 activation pathway further suggests that the inhibition of Lck could have a relatively large impact on the early TCR signaling. Collectively, simulation results point out that the Lck could be a potential target employed by PD-1 pathway to inhibit the activation of several signaling components or alter several signaling pathways. Here we used manual approach in parameterizing our model to benchmark with the experimental data. We preferred the manual approach as our primary goal was to model the network of PD-1 pathway and not to estimate the parameter values of an established pathway. However one can use a parameter estimation tool box to optimize the model parameters. Our proposed model is deterministic in nature considering all the cells in a population as identical in every aspect. Therefore the physiological basis of dose response curves with varying quantity/dose of signaling molecules such as Lck, PD-1, CD28 and CD3z in our model needs to be highlighted. A clonal population of cell exhibits large cell-to-cell variation in cellular content, shape, size and cell cycle phases due to intrinsic and extrinsic source of heterogeneity. The variations in expression of key signaling components, both in upstream and downstream pathways, ultimately lead to heterogeneous response across the population [70][71][72][73][74]. Although our model does not consider cell-to-cell variation of proteins, however the deterministic dose response curves provides an avenue to estimate the dynamic range of response and the critical amount of signal needed for the population to respond (IC 50 ) for a heterogeneous population of cells. Therefore in a way our model sets the background for predicting the outcome of the system in presence of population heterogeneity.
A key feature in PD-1 pathway is the presence of a negative feedback loop involving PD-1 and Shp2 (Fig 1). In this motif activation of PD-1 leads to recruitment of Shp2 that ultimately deactivates PD-1 via dephosphorylation. A negative feedback loop is well known for its adaptation property that allows the system to respond to the external signal and reset back to its original state even in presence of persistent signal [75]. The resetting of the system is very crucial for T cell to avoid any unwanted over-or under-reaction of the immune system. Further for coherent and collective response of a population of T cell against pathogen the effect of inevitable noise must be minimized by the regulatory pathway. Previous experimental and computational studies have demonstrated that network architecture, in particular negative feedback loop and feed forward loop, play a crucial role in reducing molecular noise [76][77][78][79][80][81]. Therefore these network motifs in the PD-1 pathway may have been selected over evolution due to their role in reducing noise in addition to their usual deterministic properties.
Early TCR signaling has numerous potential implications. Recently, a study has shown a direct role of early TCR signaling in the activation of the enzyme pyruvate dehydrogenase kinase 1 (PDHK1). Consequently, mitochondrial import of pyruvate is inhibited and aerobic glycolysis is promoted [82]. Hence, early TCR signaling is important in regulating the aerobic glycolysis, which is characteristic of activated T cells. Model developed here could be employed to understand the inhibition of several early TCR signaling molecules due to PD-1 in the case of T cell activation or T cells exhausted due to chronic viral infections.
Experiments done on membrane based reconstitution system are considered physiologically more realistic when compared to traditional solution based experiments for signaling studies [62]. Although the model recaptures several of these observations, in vivo scenario could be considerably different from the network modeled here. Future improvements in the model could potentially result in a better prediction and reproduction of in vivo scenario. Although the impact of PD-1 on Lck activity was explored here, regulation of Lck activity is very complex due to the action of several other kinases and phosphatases such as Csk and CD45 [83]. Model simulation could be made more realistic by considering the timescale of separation of TCR and PD-1 activation. This is because PD-1 expression itself is under the control of TCR signaling. In fact, recent studies have shown that the strength of TCR signaling influences PD-1 expression [84,85]. Moreover, binding dynamics of receptors on the T cell surface to the ligands present on the surface of antigen presenting cells could significantly impact downstream activation. However, this model serves as a base to which such effects could be incorporated provided reliable experimental results or kinetic measurements are available.
Finally, the proposed model is based on a set of deterministic dynamical equations. Although our deterministic model recaptures many key experimental observations however being deterministic in nature it is not capable of explaining single cell data which requires stochastic modeling. It is established that in the early T cell receptor activation pathway population heterogeneity of key signaling proteins, such as Lck, Zap70, play a crucial role [70,86]. However in this paper our main objective was to establish a mechanistic mathematical model of PD-1 pathway and stochastic calculations of the current model is definitely a scope of the future.

Parameter sensitivity analysis
To implement global sensitivity analysis, we make use of multi parametric sensitivity analysis (MPSA) where multiple parameters are selected together for testing the sensitivity. Each parameter was picked from a uniform distribution with a range ±50% of the optimized parameter value used in the model ( Table 2). We used Latin hyper cube sampling technique to create a sample of 20000 values for each parameter. Parameter samples were permuted randomly and a combination of parameters was taken as a parameter set, thus generating 20000 parameter sets. We calculated the parameter sensitivity on the time course data of various quantities. 11 different PD-1 concentration values were generated that include 0 nM and 10 logarithmically spaced between 1 and 1000 nM. In each PD-1 concentration, we calculated the sum squared deviation over 10 different time points with respect to the time course of optimized parameter. Finally we estimated the overall error by summing the squared deviation over all PD-1 concentrations. In order to determine acceptable parameter combination, we set a threshold by taking the average of the overall sum of squared error determined over 20000 parameter sets. Any parameter combination that results overall error below the threshold was considered as acceptable and otherwise it was considered unacceptable [87]. We calculated cumulative frequency distribution for individual parameter values from the acceptable and unacceptable parameter sets. See Figure E in S1 File for representative cumulative frequency distributions. To determine the sensitivity of a parameter we calculated Kolmogorov-Smirnov (K-S) statistic which quantifies the dissimilarity of two probability distributions by measuring the maximum perpendicular distance between their respective cumulative distribution functions. The more dissimilar the two distributions are the distance between the two distributions would be higher. Hence, higher the K-S statistic, higher is the sensitivity of that particular parameter [88]. Concentration of components used in the sensitivity analysis is close to physiological concentrations. Sensitivity analysis was performed separately for three different measures in the model-percentages of Shp2 recruitment, PI3K recruitment and active form of Slp76 pertaining to sensitivities of PD-1-Shp2 complex formation, early CD28 signaling and early TCR signaling, respectively.


extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Medicine

text:
A new artificial diet for western corn rootworm larvae is compatible with and detects resistance to all current Bt toxins

Insect resistance to transgenic crops is a growing concern for farmers, regulatory agencies, the seed industry, and researchers. Since 2009, instances of field-evolved Bt resistance or cross resistance have been documented for each of the four Bt proteins available for western corn rootworm (WCR), a major insect pest. To characterize resistance, WCR populations causing unexpected damage to Bt maize are evaluated in plant and/or diet toxicity assays. Currently, it is not possible to make direct comparisons of data from different Bt proteins due to differing proprietary artificial diets. Our group has developed a new, publicly available diet (WCRMO-1) with improved nutrition for WCR larvae. For the current manuscript, we tested the compatibility of all Bt proteins currently marketed for WCR on the WCRMO-1 diet and specific proprietary diets corresponding to each toxin using a susceptible colony of WCR. We also tested WCR colonies selected for resistance to each protein to assess the ability of the diet toxicity assay to detect Bt resistance. The WCRMO-1 diet is compatible with each of the proteins and can differentiate resistant colonies from susceptible colonies for each protein. Our diet allows researchers to monitor resistance without the confounding nutritional differences present between diets.

Western corn rootworm (Diabrotica virgifera virgifera LeConte, WCR) has been a challenge for maize (Zea mays L.) farmers throughout much of the United States of America (USA) for decades. In 1909, this pest was discovered to attack maize roots in Colorado 1 . Since then, it has expanded its geographic distribution to most of the maize growing regions of North America. Multiple introductions of WCR to Europe 2 have increased its global importance. More than 30 years ago, it was estimated that the species caused over $1 billion (USD) in economic losses stemming from yield loss and control costs 3 .
Early attempts at managing WCR in the USA focused solely on crop rotation to a non-host such as soybean (Glycine max (L.) Merrill) or sorghum 1 . When they became available, management tactics included the application of insecticides for larval or adult management 4,5 . Recently, biotechnology has allowed farmers to grow maize which expresses one or more proteins from Bacillus thuringiensis Berliner (Bt). Except for banded applications of soil insecticides, each of these options has failed in one or more regions found within the distribution of this pest [6][7][8][9][10][11][12][13][14] .
When the first Bt events targeting lepidopterans were registered, insect resistance management (IRM) programs were required. Some believe these programs are especially important for products expressing toxins at less than high-dose or species with lesser susceptibility to the proteins 15 , since the likelihood of survival increases when products express toxins at less than high-dose, or when some targeted species are less susceptible to the toxins. Recently, biotechnology has allowed farmers to grow maize that expresses one or more proteins from Bacillus thuringiensis Berliner (Bt) for WCR control. However, modeling efforts and laboratory selection studies SCIEnTIfIC RePoRts | (2018) 8:5379 | DOI: 10.1038/s41598-018-23738-z for current proteins targeting WCR suggested products would lose efficacy within a decade of first use [16][17][18][19][20][21] , and field data support this 10,13,14 . Bt maize products targeting WCR have not yet met the high-dose criterion 22 , and this likely is a primary reason that the refuge strategy designed to delay resistance has been ineffective with this pest 23 .
Resistance monitoring efforts to determine whether a shift in the susceptibility of randomly sampled populations was required as part of the registration process and to date has used proprietary artificial diets developed by product registrants [24][25][26] . While this effort meets the requirement for IRM programs, there have been issues with the diet toxicity assays. In addition to nutritional differences between diets, high levels of contamination often occurred. Changes in assay methodology and diet formulations have also occurred over the years for some Bt proteins, likely creating additional variability in the data. Some academic researchers have tended to use on-plant assays as opposed to proprietary diets and purified Bt proteins, which require special agreements with the owners. This barrier can be overcome with relationships between industry and academia, but often takes considerable time and negotiation. As a result, several studies have evaluated Bt resistance with on-plant assays only [10][11][12]21,[27][28][29] , and some with both on-plant and diet toxicity assays 13,14,17,19,[30][31][32] . Following initial reports of field-evolved Bt resistance using on-plant assays, the Environmental Protection Agency (EPA) proposed on-plant assays be used as a replacement for diet toxicity assays 33 and the proposed changes for resistance monitoring have been implemented 34 . Instead of diet toxicity assays, registrants will be required to use on-plant assays to determine shifts in susceptibility, though they may continue to conduct diet toxicity assays in conjunction with on-plant assays.
Efforts have been underway to improve WCR diets and the first new publicly available diet was recently published 35 . Referred hereafter as "WCRMO-1", this diet has improved characteristics compared to the only other publicly-available diet 36 and all the proprietary diets used in the current study (Hibbard/Coudron labs unpublished data). Efforts are continuing to further improve the WCRMO-1 diet for weight gain and developmental rate. An optimum diet for WCR should match maize for survival, weight, and developmental rate parameters. If these parameters can be optimized, or can no longer be improved, then we believe a single diet should be used as a universal diet in all public diet toxicity assays. As a step toward this longer-term goal, we compared the WCRMO-1 diet to the appropriate proprietary diets for compatibility with Bt proteins. We also tested susceptible and resistant WCR colonies on the WCRMO-1 diet to evaluate the ability of the diet to detect differences in susceptibility to Bt proteins in these two phenotypes.

Results
Lethal effects. The WCRMO-1 diet was tested with each Bt protein alongside the proprietary diet of the corresponding registrant using a Bt-susceptible colony (Brookings-ND, Table 1). Average mortality for the buffer dose (control) on WCRMO-1 was less than 10 percent for all but one instance (Fig. 1d). The average mortality for the buffer dose on proprietary diet C was between 12 and 13 percent ( Fig. 1b and c), while the buffer dose for proprietary diets corresponding to Cry34/35Ab1 (proprietary diet A) and Cry3Bb1 (proprietary diet B) toxins was less than five percent ( Fig. 1a and d). Based on overlapping confidence intervals between concentrations required to kill 50 percent of insects (LC 50 ), WCRMO-1 provided similar data as the proprietary diet for each protein except for one ( Table 2). WCRMO-1 produced a significantly lower LC 50 value for eCry3.1Ab compared to the corresponding proprietary diet, labelled as "proprietary diet C" ( Table 2). This means less eCry3.1Ab was required on WCRMO-1 to kill half of the susceptible colony tested when compared to proprietary diet C.
When Bt-selected colonies were evaluated on WCRMO-1, significant differences in LC 50 values were found for each Bt protein between the Brookings-ND colony and the resistant colonies (Table 2). This was the case for each selected colony on each of the four Bt toxins evaluated. Overall, resistant colonies showed minimal mortality in response to toxins. For Cry3Bb1, mCry3A, and eCry3.1, all resistant colonies had maximum mortality under 14 percent regardless of the dose (Fig. 1f-h). Only one of the DAS-59122-selected colonies reached maximum mortality 50% while the other was 41% following exposure to the highest dose Cry34/35Ab1 toxins (Fig. 1e). Sublethal effects. Dry weight was the only parameter where a significant difference was observed between susceptible larvae fed Bt toxin on WCRMO-1 versus a proprietary diet. For eCry3.1Ab, the concentration required to cause a 50 percent reduction in dry weight (EC 50 ) for the Brookings-ND colony fed WCRMO-1 diet was significantly greater than when fed proprietary diet C ( Table 2). No other significant differences were found in EC 50 values between the WCRMO-1 diet and the respective proprietary diet for any of the other toxins  Table 2). Both DAS-59122-7-selected colonies had similar EC 50 values compared to the Brookings-ND colony when exposed to the Cry34/35Ab1 proteins on the WCRMO-1 diet ( Table 2). Both MON88017-selected colonies had significantly greater EC 50 values for Cry3Bb1 protein on the WCRMO-1 when compared to the susceptible colony ( Table 2). The Brookings-ND colony had a significantly lower EC 50 value for eCry3.1Ab protein on the WCRMO-1 diet compared to the 5307-S colony (Table 2). Lastly, the MIR604-S colony had a significantly greater EC 50 value for mCry3A protein when compared to the susceptible colony on the WCRMO-1 diet ( Table 2). There were large differences in susceptible WCR larval dry weight on the buffer dose between the WCRMO-1 and proprietary diet C ( Fig. 2b and c). For both toxins, the insects reared on proprietary diet C weighed one-fifth as much as insects reared on WCRMO-1 after 10 d (Fig. 2b and c). Insects from Brookings-ND exposed to the buffer dose of Cry34/35Ab1 weighed more on WCRMO-1 than those reared on proprietary diet A. All but one colony (DAS-59122-7-S1) had greater dry weight than the Brookings-ND insects on the buffer dose ( Fig. 2e-h). Additionally, some colonies had greater dry weight at the lowest Bt dose than for the buffer (Fig. 2a,d-f).
Molting was significantly inhibited for all susceptible colonies on all Bt proteins ( Fig. 3a-d). The dose that causes 50 percent molting inhibition (MIC 50 ) was significantly greater for all the Cry3-selected colonies than the Brookings-ND colony ( Table 2). The MIC 50 value was not significantly different for the two Cry34/35Ab1-selected colonies and the Brookings-ND colony ( Table 2).

Discussion
Here, we demonstrate for the first time that a single artificial diet can be used to detect resistance with all Bt proteins currently targeting WCR. Previously, different proprietary proteins were evaluated on the proprietary diet from the company that produced the protein. Use of a single diet allows for direct comparisons of toxicity between proteins without the confounding effects of nutritional differences. In general, performance of susceptible WCR was similar between the proprietary diets and the WCRMO-1 diet ( Table 2, Figs 1-3). We also documented that differences between WCR colonies selected for resistance to each Bt protein and a control WCR population can be detected in diet toxicity assays with the WCRMO-1 diet ( This study adds to the limited number of studies which have documented a significant difference in LC 50 values between colonies or populations in Cry34/35Ab1 diet toxicity assays 13,14 . Resistance ratios in the previous studies were relatively minor compared to the more than 40 fold difference in LC 50 values of the two DAS-59122-7-selected populations and the control population evaluated here (Table 2). Interestingly, a colony originally described by Lefko et al. 30 and Nowatzki et al. 37 showed no difference in LC 50 values when compared to a control colony after being selected for more than 30 generations on DAS-59122-7 maize (which produces the Cry34/35Ab1 toxin) using proprietary diet A (pers. communication from Stephen Thompson, DuPont Pioneer). Using the WCRMO-1 diet, significant differences in LC 50 values between both DAS-59122-7-selected colonies and control colonies were documented with less than 20 generations of selection ( Table 2, Figs 1-3). The difference in LC 50 values did not translate to significant differences in EC 50 values; however, the DAS-59122-7-S2 colony did have a significantly greater MIC 50 value than the control WCR colony ( Table 2). These results suggest a possibility for more sensitive screening of Cry34/35Ab1 proteins with WCRMO-1. However, there is also a chance that the difference in LC 50 values detected on the WCRMO-1 diet is unique to the DAS-59122-7-selected colonies in this study. Additional comparisons of proprietary diet A and WCRMO-1 diets are needed to determine whether the basis of this detection is the result of nutritional or genetic factors, or a combination of the two factors.
Although we do not have the formulations for the proprietary diets tested here, life history parameters indicate that nutritive qualities varied among the diets. We previously found significant differences in molting rate, weight gain, and survival among the proprietary and WCRMO-1 diets (unpublished data). Research with lepidopteran insects has shown that the toxicity of Bt proteins can vary when protein to carbohydrate ratios are manipulated 38,39 . While lepidopteran and coleopteran insects are quite different, we may be able to draw parallels to this study. One proprietary artificial diet appears to lack significant nutritive qualities and/or attractiveness as a food source. Insects on proprietary diet C with buffer alone for eCry3.1Ab and mCry3A proteins weighed less than 18 percent of insects reared on the WCRMO-1 diet with buffer alone (Fig. 2b and c). Such a difference in nutritive qualities could explain the significant difference in dry weight observed between larvae fed WCRMO-1 or the proprietary diet C. However, this difference in weight gain between diets did not cause a significant difference in LC 50 for the mCry3A protein. Under conditions where an insect ingests food at a normal rate (i.e. rate on maize), then the LC 50 is representative of the amount needed to kill targeted insects in the field. However, when a diet is less attractive, then the LC 50 will likely shift towards a greater value as less diet is ingested. Larvae were capable of molting in the 10-day time period used for these experiments, unlike previously reported assays with WCR which were terminated at five or six days depending on the protein. Previously reported assays had shorter time frames primarily due to contamination (EPA 2016). Longer assays for slower acting products (e.g. dsRNA) are possible using techniques described here and in Huynh et al. 35 . Consequently, we were able use molting as a measure of resistance, similar to studies with Ostrinia nubilalis (Hübner) 40 . Here, exuviae were clearly visible through the sealing film and provided a non-invasive determination of sublethal effects that could be recorded over time, whereas dry weight could only be collected at the end of each experiment. Recording molting data as a measure of resistance may allow for additional studies where researchers may look at genetic factors, microbes, or other variables. This approach may provide a more sensitive measure for detecting resistance in populations where no differences in LC 50 or EC 50 exist. For example, the DAS-59122-7-S2 colony had a significantly greater MIC 50 value when compared to the Brookings-ND colony while the DAS-59122-7-S1 showed no significant difference (Table 2). This was the only difference noted between the two DAS-59122-7 selected colonies. Future studies should include these data, when possible. Before molting data can be collected, a diet must have adequate nutrition to allow first instar larvae to molt into the second instar. Proprietary diet C did not produce second instar larvae, even on diet with buffer alone (Fig. 3b and c). All other diets produced second instar larvae on diet with buffer alone (Fig. 3a and d).
Nutritional improvements of an artificial diet 35 and validation of compatibility with all current Bt proteins documented here, support the utility of the diet toxicity assay in resistance monitoring efforts. We recommend a standardized approach to how these diet toxicity assays are conducted to reduce data variability in IRM programs. Some of this variability is likely due to the alterations of some diets over the years for various reasons. Additionally, the number of days for observation varies between proteins complicating comparisons of Bt proteins. These factors likely explain some of the variability that has been reported within diet toxicity assays over the years. Shifts in susceptibility may be masked by these factors. Comparisons between current and future products could help to evaluate novel modes of action and physiological effects of protein exposure.
WCR has evaded management tactics for more than a century. Some populations of this species have developed resistance to crop rotation 9 , broadcast insecticides 6-8 , and transgenic maize expressing Bt proteins 10,11,13,14,27,29 . Even with the advent of biotechnology, WCR continues to be a pest of maize. Current and future products targeting the species will require IRM programs to ensure product viability and longevity. Both plant and diet assays have value in detecting resistance 13 .

Materials and Methods
Artificial Diet and Bt Protein. Lyophilized Cry34/35Ab1 proteins and the corresponding proprietary diet were provided by Dow AgroSciences. Cry3Bb1 protein was provided in solution by Monsanto Company in addition to their proprietary diet. The lyophilized mCry3A and eCry3.1Ab proteins and their proprietary diet were provided by Syngenta Biotechnology. The WCRMO-1 diet was produced by the USDA-ARS Biological Control of Insects Research Laboratory (BCIRL) in Columbia, MO 35 . Buffers to suspend the protein in solution were prepared at BCIRL (Table 1). The sodium citrate buffer was stored at 4 °C immediately after preparation. The sodium carbonate bicarbonate buffer was stored at room temperature in a container which blocked all light.

Insects.
Eggs from the non-diapausing, susceptible colony, hereafter known as Brookings-ND, were provided by the USDA-ARS laboratory in Brookings, SD. This colony is derived from a field population collected prior to the release of transgenic crops 41 , so this colony should be susceptible to all Bt proteins active against WCR.
Several resistant colonies were used in this study, all of which have been described in other studies. Two MON88017-selected colonies, MON88017-S1 and MON88017-S2, exhibited resistance to Cry3Bb1-expressing plants (event MON88017) and protein in Zukoff et al. 13 . MON88017-S1 was derived from the Canby population while MON88017-S2 was derived from Hills population in Zukoff et al. 13 . Both MON88017-selected colonies were reciprocally crossed with adults from the Brookings-ND colony and then reared on maize containing event MON88017 for eleven generations post Zukoff et al. 13 at the time Cry3Bb1 diet toxicity assays were conducted.
Two DAS-59122-7-selected colonies, hereafter known as DAS-59122-7-S1 and DAS-59122-7-S2, were also used in this study. The DAS-59122-7-S1 colony was derived from a field population described in Ludwick et al. 14 . The field population had been exposed to SmartStax, a maize product which expresses both Cry34/35Ab1 (event DAS-59122-7) and Cry3Bb1 proteins. DAS-59122-7-S1 had incomplete resistance to both Cry3Bb1 and Cry34/35Ab1 proteins 14 . Eggs from this population were reared on non-Bt corn and resulting adults were then reciprocally crossed with the Brookings-ND colony to produce a non-diapausing DAS-59122-7-S1 colony. The DAS-59122-7-S2 colony was started from pooled eggs of three field populations reared on non-Bt maize. The resulting adults were reciprocally crossed with Brookings-ND at differing points earlier and then selected on maize containing event DAS-59122-7 similarly as in Deitloff et al. 21 . The DAS-59122-7-S1 and DAS-59122-7-S2 colonies were selected on DAS-59122-7 plants for six and 18 generations under laboratory/greenhouse conditions, respectively, at the time Cry34/35Ab1 diet toxicity assays occurred.
Two additional colonies with resistance to mCry3A 31 and eCry3.1Ab 19 were also tested. The colony selected on eCry3.1Ab-expressing corn will hereafter be known as 5307-S. We tested this strain after 35 generations of selection on maize containing event 5307 or an earlier event also expressing eCry3.1Ab. The MIR604-selected colony will hereafter be known as MIR604-S. The MIR604-S colony was evaluated after more than 50 generations of selection on MIR604 maize which expresses mCry3A protein.
Equipment Sterilization. The WCRMO-1 diet was prepared in new 96-well immunoassay plates (product #3370, Corning Inc., Corning, NY) using sterile techniques described in Huynh et al. 35 All proprietary diets were prepared in 96-well immunoassay plates by the registrants. For shipment, proprietary diets were placed in sterile SCIENtIfIC REPORTS | (2018) 8:5379 | DOI:10.1038/s41598-018-23738-z packaging and shipped to the BCIRL in styrofoam containers with ice packs. Diets were stored in a 4 °C refrigerator and infested within two weeks of arrival. Lids for diet were only removed inside a UV-sterilized biosafety cabinet. Paintbrushes, insect pins, deli containers, coffee filters, and beakers were all sterilized by UV lights in a biosafety cabinet prior to use. A spray bottle with autoclaved DI water was used to keep the coffee filter and neonate larvae moist.
Egg Sterilization. The methods used to surface sterilize eggs were modified from descriptions by Pleau et al. 36 . Briefly, eggs from each colony were incubated in Petri dishes with 70 mesh sieved soil within an incubator at 25 °C in complete darkness until hatch started. Once approximately 10 percent of the eggs hatched, each dish was washed through a 60 mesh sieve. Remaining eggs were submerged in a beaker with water where debris (soil, fungus, hatched egg shells) floated to the water's surface and were then decanted. As much excess water as possible was poured off before the eggs were submerged in 20 mL of disinfectant (undiluted Lysol ® ) for three minutes. The disinfectant was then decanted and the eggs were triple rinsed with distilled water. Approximately 20 mL of 10 percent buffered zinc formalin (10 percent formaldehyde) covered the eggs for three minutes. Again, a triple rinse with distilled water was used to remove any chemical residues. Eggs were then dispensed onto a UV-sterilized coffee filter with a 1.5 ml transfer pipette. The coffee filter was then placed in a UV-sterilized 16 oz. Solo ® deli container with holes (#0 insect pin) punched in the lid to allow for air exchange. Incubation, typically less than three days, occurred until enough eggs hatched for infestation. If sufficient numbers of larvae did not hatch, then the coffee filter with eggs were transferred to a recently UV-sterilized deli container each additional day and placed back in the incubator. Contamination tends to increase when eggs are used more than three days after sterilization (personal observation). The deli container with larvae on its wall was placed in a sterilized biosafety cabinet. Neonate larvae were then transferred to diet using sterile equipment to prevent any possible contamination. Sealing film (Excel scientific, Inc., Thermalseal RTS RM , TSS-RTQ-100) was placed over the entire plate to prevent escape. A single hole (#0 insect pin) was poked into the sealing film of each infested well for diffusion of oxygen.
Diet Toxicity Assays. All diet toxicity assays were conducted for a length of 10 days. Lyophilized proteins were dissolved in solution and then a serial dilution was made (Table 1). These solutions were then overlaid on artificial diet and allowed to dry before larvae were placed on the diet. Each row (12 wells) per plate received one dose and one colony, thereby constituting a single replication. Plates with more than 25 percent mortality on the control (buffer dose, 0 µg/cm 2 ) were excluded from the study. The number of replications included in this study varied between diets and proteins ( Table 2). Survival and molting were recorded prior to the collection of larvae. Surviving larvae were placed in ethanol for each dose and plate and then dried in an oven at 50 °C for one week (Blue M Therm Dry Bacteriological Incubator, Model #602752). After drying, larvae were weighed on a digital microbalance (Sartorius TM Cubis TM , 6.6S), which weighed samples to one thousandth of a milligram. Statistical Analyses. Mortality was calculated by dividing the number of dead larvae by the initial number of larvae infested per dose. An average mortality value of the buffer dose was calculated. Then, mortality at each subsequent dose was divided by the average mortality of the buffer dose and multiplied by 100 42 . These mortality percentages were analyzed with a probit analysis to generate LC 50 and 95 percent confidence intervals using SAS software version 9.2 (SAS Institute). Values were considered significantly different when 95 percent confidence intervals did not overlap.
Dry weight was calculated by dividing the total weight per dose by the initial number of larvae infested, effectively giving all dead larvae a weight of 0. The average dry weight per larva recovered from the buffer dose was then averaged across replications to generate a mean value. The dry weight of larvae recovered from each dose was then divided by the buffer's mean dry weight and multiplied by 100 to generate a percentage relative to the buffer dose. These data were then analyzed with a nonlinear probit analysis described by Marçon et al. 43 in SAS. This analysis calculated the EC 50 and 95 percent confidence intervals.
Molting was calculated by dividing the number of molted individuals per dose by the initial number of larvae infested and multiplying by 100 to obtain a percentage. The average molting rate for each colony on the buffer dose was calculated. Each subsequent dose was divided by the buffer dose's average molting rate for the corresponding colony 42 . These molting rates were then analyzed in SAS with a probit analysis to generate MIC 50 values and 95 percent confidence intervals. MIC 50 values between colonies on the same protein were considered significantly different when confidence intervals did not overlap. Data Availability. All pertinent data are found in the figures and tables. Requests for data and additional information should be submitted to the corresponding author.


s2fieldsofstudy: Biology Chemistry Medicine

text:
PM2.5 downregulates miR-194-3p and accelerates apoptosis in cigarette-inflamed bronchial epithelium by targeting death-associated protein kinase 1

Background Persistent exposure to cigarette smoke or biomass fuels induces oxidative stress and apoptosis in bronchial epithelium, which is one of the most important pathogenic mechanisms of chronic obstructive pulmonary disease (COPD). Fine particulate matter (PM2.5) is an aggravating risk factor of COPD exacerbation. Animal evidence showed PM2.5accelerated lung inflammation and oxidative stress in COPD mice, but the mechanism is still not clear. Recently, we found that miR-194-3p is a novel biomarker of both COPD and PM2.5 exposure, and miR-194 family has been reported to be involved in cell proliferation and apoptosis. Thus, we propose a hypothesis: PM2.5 can accelerate apoptotic response of airway epithelial cells in COPD and miR-194 is a potential involved regulator. Materials and methods Human bronchial epithelial cells (HBEpiCs) were treated with normal media, cigarette smoke solution (CSS) and PM2.5-CSS for 24 h. miR-194-3p mimics, inhibitors and scrambled controls were non-transfected or pre-transfected into HBEpiCs for 48 h. MircroRNAs and mRNA expression were quantified by qRT-PCR. Protein expression was analyzed by western blotting. Caspase activities, mitochondrial membrane potential and TUNEL-positive cells were detected to analyze apoptosis. Bioinformatics and luciferase analysis were used to identify the predicted binding site of miR-194-3p and potential targets. Results In our study, we found that PM2.5 significantly aggravated apoptosis in cigarette-inflamed HBEpiCs. miR-194-3p was dramatically downregulated in PM2.5-CSS-treated HBEpiCs. Bioinformatics and luciferase experiments reported that death-associated protein kinase 1 (DAPK1), regulating caspase 3 activities in apoptosis, was directly targeted by miR-194-3p. Inhibition of miR-194-3p increased DAPK1 expression and apoptosis in normal HBEpiCs. Importantly, overexpression of miR-194-3p suppressed apoptosis in PM2.5-CSS HBEpiCs. Conclusion These results suggested that miR-194-3p was a protective regulator involved in apoptosis pathway and a potential therapeutic target for treatment of bronchial epithelial injury aggravation induced by PM2.5.

Introduction
Chronic obstructive pulmonary disease (COPD) is characterized by persistent and aggressive airflow limitation and lung inflammation, 1 and COPD has become the third highest cause of death in the world. 2 One of the most important pathogenic mechanisms of COPD is inflammatory response and abnormal apoptosis in airway epithelium, due to long-term exposure to cigarette smoke or biomass fuels. [3][4][5][6] Not only that, the inflammation and cell apoptosis induced by cigarette smoke remain in the lungs after quitting smoking. 7,8 Thus, COPD is not easily managed and causes severe economic burden. Air pollution is one of the most challenging problems in COPD management. According to the World Health Organization, approximately 14% of 3.7 million deaths in 2012 were associated with COPD and respiratory infection. 9 Particulate matter (PM) is a complex mixture of small particles and liquid droplets in the air, and fine particulate matter (PM 2.5 ) is the PM with a diameter of 2.5 μm or less. 10 Once inhaled, PM 2.5 deposits in lung tissues and diffuses in blood inducing airway and systematic inflammation. 11,12 In recent decades, epidemiological studies demonstrated that PM 2.5 is an aggravating risk factor of COPD exacerbation. The latest systemic review showed that every 10 μg/ m 3 increase of PM 2.5 was associated with a 2.5% increased risk of COPD-related emergency visits and hospital admissions. 13 Animal studies have shown that PM 2.5 accelerated the original lung inflammation in COPD mice, 14 but the mechanism is still not clear.
MicroRNAs (miRNAs) are a class of endogenous small noncoding RNAs regulating cell development, proliferation, differentiation and death. 15,16 miRNAs are confirmed to play an important role in the pathogenesis of COPD. For example, the downregulation of miR-146a 17 and the upregulation of miR-223 18 were associated with the inflammatory responses of COPD. Meanwhile, miRNAs are sensitive biomarkers to PM 2.5 . Expression of miRNAs such as miR-9, miR-10b and miR-21 were significantly altered after exposure to PM 2.5 , 19 suggesting that miRNAs might be ideal tools for understanding the mechanism of PM 2.5 in COPD.
In our unpublished data, we found that miR-194-3p was a novel biomarker of COPD, and miR-194-3p was also correlated with PM 2.5 -induced lung function decline (Zhou et al, unpublished material, 2018). miR-194 family participates in cell proliferation and senescence. [20][21][22][23] Recently, miR-194 was confirmed to express in non-small lung cancer cells. 24 However, miR-194 has not been discussed either in COPD or in PM 2.5 exposure. Thus, we suggested a hypothesis that PM 2.5 regulates the expression of miR-194-3p and accelerates apoptosis of cigarette-inflamed bronchial epithelial cells. In this study, we aimed to find out whether PM 2.5 accelerated apoptosis and altered the expression of miR-194-3p in cigarette-inflamed human bronchial epithelial cells (HBEpiCs). If so, we aimed to further explore the underlying mechanism.

Cigarette smoke solution (CSS)
According to the method of Blue and Janoff, 25 smoke from one cigarette was pumped slowly into a negative pressure tube that contained 10 mL of sterile medium via a 50 mL syringe. After being filtered through 22 μm filters, the concentration of the original CSS was considered to be 100%. The concentration of final CSS was 5%.

PM 2.5
The preparation of PM 2.5 has been already described. 26 Briefly, PM 2.5 was gathered by a PM 2.5 high volume sampler system (Staplex ® No PM-2.5 SSI; USA) set on the roof of Peking University First Hospital, which was located very close to one of the busiest roads in central Beijing. Glass fiber filter paper (20.3×25.4 cm) was used to collect PM 2.5 with an air flow rate of 1.13 m 3 /min. Then, the paper was sonicated three times with deionized water. Finally, the filtered products were lyophilized and refrigerated at −80°C. Final extracted PM 2.5 was suspended and sonicated in culture media to ensure that the stimulus dose for cells was 500 ng/mL.

Cell culture
HBEpiCs (ScienCell Research Laboratories, Inc., San Diego, CA, USA) were cultured in bronchial epithelial cell medium (ScienCell Research Laboratories, Inc.) at 37°C in a 5% CO 2 atmosphere. In total, 30 nM miR-194-3p mimics, inhibitors or scrambled controls were non-transfected or pre-transfected into HBEpiCs with DMEM with 3% fetal calf serum (Thermo Fisher Scientific, Waltham, MA, USA) for 48 h. Transfection of each reagent was performed using Lipofectamine ® 3000 (Thermo Fisher Scientific). The mimics and negative control of miR-194-3p were chemically synthesized as double-stranded RNAs (GenePharma Co., Ltd., Shanghai, China) using the following sequence (5′→3′): mimics sense CCAGUGGGGCUGCUGUUAUCUG and anti-sense GAUAACAGCAGCCCCACUGGUU; negative control sense UUCUCCGAACGUGUCACGUTT and antisense ACGUGACACGUUCGGAGAATT. The inhibitors and inhibitor-negative control were synthesized as singlestranded RNAs using the following sequence: inhibitor CAGAUAACAGCAGCCCCACUGG and inhibitor-negative control CAGUACUUUUGUGUAGUACAA. Before analysis, HBEpiCs were treated with normal media, CSS or PM 2.5 -CSS for 24 h.

Quantitative assessment of mirna and mrna
Total RNA was isolated from HBEpiCs via TRIzol ® Reagent (Thermo Fisher Scientific). Kit (Thermo Fisher Scientific), while for analysis of miRNA, 600 ng of total RNA was reverse transcribed by Mir-X miRNA First-Strand Synthesis Kit (TAKARA, Tokyo, Japan). qRT-PCR was performed on StepOnePlus Real-Time PCR System (Thermo Fisher Scientific) in a 20 μL reaction that contained 1 μL of cDNA and Power SYBR ® Green PCR Master Mix (Thermo Fisher Scientific). Relative expressions of death-associated protein kinase 1 (DAPK1) and miR-194-3p were determined. β-Actin and U6 were used as the reference genes, respectively. The sequences of primers for miRNA analysis were as follows (5′→3′): U6 from Mir-X miRNA First-Strand Synthesis Kit; hsa-miR-194-3p ATTATTCCAGTGGGGCTGCT. Gene primers for mRNA analysis: hsa-β-actin forward CTGTGGCATC CACGAAACTA, reverse GTGTTGGCGTACAGGTCTT; hsa-DAPK1 forward GTGGATGGTCATTGCAGTTTAAG, reverse TACTGGAGGATGAGAGATGGAG.

TUNEL analysis
TUNEL (Roche Molecular Systems, Inc., Basel, Switzerland) staining was used for measuring DNA fragmentation. HBEpiCs were fixed before detection. The procedures were carried out according to the manufacturer's instructions.

Caspase-3/7-positive cell detection and TMrM detection
The CellEvent™ Caspase-3/7 Green detection reagent (Thermo Fisher Scientific) labels nuclei of caspase 3/7-positive cells to report apoptosis. Tetramethylrhodamine, methyl ester reagent (TMRM; Thermo Fisher Scientific) was used to observe the mitochondrial membrane potential. HBEpiCs were loaded with 50 nM TMRM followed by 20 μM CellEvent™ detection reagent. Each reagent was incubated for 30 min, and fluorescence signals were observed in living cells.

Statistical analysis
All experiments were conducted independently at least 3 replicates. Continuous variables are presented as mean ± SD. Data were analyzed using SPSS 13.0 (SPSS Inc., Chicago, IL, USA), and bar graphs were protracted using Prism (version 5.0, GraphPad Software Ltd, San Diego, CA, USA). If the data distribution was normal, comparisons among three or more groups were estimated with an ANOVA test and comparisons between two groups were estimated by twotailed Student's t-test. If the data distribution was not normal, comparisons were estimated by Wilcoxon signed-rank test for nonparametric analysis. p,0.05 was considered statistically significant.

PM 2.5 increased the apoptotic response in cigarette-inflamed HBEpiCs
To verify the effect of PM 2.5 in cigarette-inflamed bronchial epithelium, HBEpiCs were cultured with normal media, CSS or PM 2.5 -CSS for 24 h. Typical features of apoptosis were evaluated including mitochondrial membrane potential loss, DNA fragmentation and caspase activation. Caspase activities and mitochondrial membrane potential were detected in living cells ( Figure 1A). Compared with normal HBEpiCs, CSS-and PM 2.5 -CSS-treated HBEpiCs showed a significant loss of mitochondrial membrane potential with AKT plays a critical role in cell survival and resistance, 27 while caspases play an essential role in apoptosis activation. 28,29 Expressions of these proteins were detected using western blot analysis ( Figure 1B). The results showed that AKT phosphorylation decreased, while cleavage of caspase 3 increased in CSS-treated HBEpiCs. In PM 2.5 -CSS-treated cells, AKT phosphorylation further decreased and cleavage of caspase 3 further increased. Cleavage of caspase 9 increased significantly in PM 2.5 -CSS cells but showed no significant change in CSS cells. The abovementioned results demonstrated that PM 2.5 aggravated apoptosis in CSS-treated HBEpiCs.  00±0.13, p=0.019). The difference between CSS-and PM 2.5 -CSS treated HBEpiCs was significant ( p=0.003).
DaPK1 was a potential target of mir-194-3p and upregulated in PM 2.

-CSS HBEpiCs
To explore the mechanism of miR-194-3p in PM 2.5 -CSS HBEpiCs, three bioinformatics algorithms, miRanda, TargetScan and miRWalk, were searched to find potential regulatory targets of miR-194-3p. All algorithms predicted DAPK1 as a target of miR-194-3p and the putative target sequence is located in 495-501nt of the 3′-UTR of human DAPK1 mRNA ( Figure 2B). Then, expressions of DAPK1 mRNA and protein were detected in normal, CSS-and PM 2.5 -CSS-treated HBEpiCs. There was no significant difference in the expression of DAPK1 mRNA (normal 1.00±0.06 vs CSS 0.95±0.04 vs PM 2.5 -CSS 0.96±0.06, p=0.433) ( Figure 2C). On the other hand, compared with control, the expression of DAPK1 protein showed a significant upregulation in PM 2.5 -CSS-treated HBEpiCs (normalized gray value: PM 2.5 -CSS 2.88±0.44 vs control 1.00±0.00, p,0.001), but it was not significantly different in CSS-treated HBEpiCs (normalized gray value: CSS 1.12±0.31 vs control 1.00±0.00, p=0.662). The difference in the expression of DAPK1 protein between CSS and PM 2.5 -CSS HBEpiCs was statistically significant ( p,0.001) ( Figure 2D).

DAPK1 mRNA was directly targeted by mir-194-3p
Next, luciferase analysis was used to identify the predicted binding site of miR-194-3p and DAPK1. The DAPK1 mRNA Decreased miR-194-3p directly upregulated the expression of DAPK1 protein and induced activation of caspase 3 To explore the potential regulation relationship between DAPK1 and miR-194-3p, normal HBEpiCs were transfected with miR-194-3p mimics, inhibitors and scrambled controls. The efficiency of transfection and the expression of DAPK1 mRNA were detected by qRT-PCR, and the expression of DAPK1 mRNA showed no significant difference among different HBEpiCs ( Figure 3C). The expression of DAPK1 protein was detected using western blot analysis. The expression of DAPK1 protein increased significantly after inhibition of miR-194-3p (normalized value: inhibitors 3.12±0.92 vs null 1.00±0.00, p,0.001), whereas the expression decreased significantly after overexpression of miR-194-3p (normalized value: mimics 0.40±0.16 vs null 1.00±0.00, p=0.006) ( Figure 3D). Previous research has reported that DAPK1 is an upstream regulator of caspase 3, which activates apoptosis. 30 p=0.005). These results demonstrated that overexpression of miR-194-3p inhibited the production of DAPK1 and suppressed apoptosis in PM 2.5 -CSS HBEpiCs.

Discussion
miRNAs are acknowledged as key regulators in posttranscriptional modification. After constituting RNA-induced silencing complex, miRNAs bind directly to the 3′-UTR of target mRNA leading to mRNA degradation (complete complementary pairing) or translation suppression (incomplete complementary pairing). The stable expression of miRNAs in specific tissues and the conservative role of their mechanism in different species indicate their importance, with an enormous research and therapeutic potential. 15,16 Several miRNA therapeutics such as mimics and inhibitors have been used in clinical trials or basic animal research. The most famous clinical appliance is miravirsen, an miR-122 inhibitor for the treatment of hepatitis C virus. 33 So far, several miRNAs have been identified as potential therapeutic targets of lung diseases. For example, miR-16 mimics induce apoptosis and suppress proliferation of non-small cell lung cancer cells. 34 miR-155 35 and miR-19a 36 are considered to be therapeutic targets of asthma. Nowadays, many experts have claimed that miRNAs are ideal therapeutic targets for COPD in the present and future. 37,38 Hundreds of miRNAs in serum, plasma, sputum and lung tissues have been profiled to be associated with COPD. Several miRNAs have been identified to have therapeutic effects on COPD in vitro. Osei et al 17 reported that miR-146-5p mimics reduced the release of IL-8 in COPD pulmonary fibroblasts. Ezzie et al 18 found that miR-15b was significantly upregulated in stage IV COPD, and overexpression of miR-15b directly targets SMAD7 and then downregulates TGF-β in human bronchial epithelial cell line BEAS-2B. Hassan et al 39 found that miR-199a-5p mimics effectively downregulated the unfolded protein response and decreased the release of cytokines in a1-antitrypsin-deficient monocytes.
In the meantime, epidemiological evidence has clarified that COPD patients are more vulnerable to PM 2.5 . 40 Once PM 2.5 is inhaled, COPD patients suffer from significant dyspnea and lung function decrease. 41 Meta-analysis demonstrated that a 10 μg/m 3 increase in PM 2.5 was associated with a 2.36% (95% CI 1.00%-3.73%) increased risk of COPD-associated hospital admissions. 42 Short-term exposure to outdoor air pollution increased the mortality rate of COPD by 1%, 1% and 6% in China, the European Union and the USA, respectively. Long-term exposure increased the mortality rate by 10%. 43 Oxidative stress is the main effect of PM 2.5 . Researchers suggest that the inflammatory response involved in PM 2.5 -induced lung injury is very similar to the inflammation involved in COPD. 9,44 For example, PM 2.5 activated NF-κB-dependent pathway inducing release of IL-6 and monocyte adhesion, 45 just like the reaction in COPD. Ye et al 46 found that PM 2.5 promoted bronchial smooth muscle cell migration and airway remodeling, just as in COPD. So, this suggests that the mechanism of PM 2.5 in COPD is that PM 2.5 aggravates the original lung inflammation of COPD. In an animal model, Gu et al 14 found that PM 2.5 accelerated the original Th1 and Th17 inflammation in lung tissues of COPD mice via Notch signal pathways.
miRNAs are also sensitive to air pollution, 19 and so miRNAs might be regulators in the mechanism of PM 2.5 in COPD. Previously, we found that 10 miRNAs (miR-691, miR-181a, miR-146a, miR-146b, miR-21a-5p, miR-129, miR-155, miR-139-5p, miR-21a-3p and miR-340) were upregulated in lung tissues of mice exposed to PM 2.5 . 47 However, because the design of these air pollution-related miRNA studies were different, there are still no consensus from different researches. Our unpublished preliminary study firstly identified that blood miR-194-3p was decreased significantly in untreated COPD patients. Then, we verified that miR-194-3p was decreased in healthy adults after PM 2.5 exposure, and miR-194-3p showed statistically positive correlation with lung function at lag1 and lag0. This inspired us to conclude that miR-194-3p might be an important regulator participating in the mechanism of PM 2.5 on COPD.
Members of the miR-194 family are novel markers that were first reported in the intestine and kidney. Previous evidence from tumor research showed that miR-194 was able to inhibit cell proliferation [20][21][22][23] and suppress inflammation 48,49 and epithelial mesenchymal transition. 21,50,51 In the miR-194 family, miR-194-3p is reported less than miR-194-5p. Alteration of miR-194-3p has been reported in chronic hepatitis B infection 52 and necrotizing enterocolitis. 53  In the present study, we found that PM 2.5 aggravated apoptosis of cigarette-inflamed HBEpiCs, and we verified that miR-194-3p declined significantly in PM 2.5 -CSS-treated HBEpiCs. Furthermore, we found that miR-194-3p directly targeted DAPK1. DAPK1, belonging to a family of Ser/Thr Kinase, is a regulator of cell apoptosis, and DAPK1 acts as a mediator in apoptosis. 55 Previous research found that DAPK1 mediates the pro-apoptotic activity via TGF-β death signals, 56 p53 pathway 57 and NF-κB signaling pathways. 58 There is sufficient evidence showing that DAPK1 is the upstream regulator of caspase 3 in apoptosis. Apoptosis in cells transfected with DAPK1 increased with activation of the caspase 3-dependent pathway, 30 while caspase 3 activation decreased in DAPK1 knockdown cells compared to control siRNA-transfected cells. 59 Our study found that inhibition of miR-194-3p upregulated the expression of DAPK1 and cleaved caspase 3 in normal HBEpiCs. Not only that, overexpression of miR-194-3p decreased DAPK1 and caspase 3 cleavage, resulting in the suppression of apoptosis in PM 2.5 -cigarette smoke exposed HBEpiCs. This suggested that miR-194-3p might be a protective miRNA and a potential therapeutic target.
There are some problems that need to be discussed in the future. First, as miR-194-3p decreased significantly only in PM 2.5 -CSS HBEpiCs not in CSS HBEpiCs, we did not discuss the role of miR-194-3p in CSS HBEpiCs. Whether miR-194-3p is a regulator in COPD should be verified in future using more appropriate cells from patients or animal models. Secondly, we found that miR-194-3p targeted DAPK1 directly and activated caspase 3 downstream; whether there were different apoptotic pathways involved was not analyzed in our study. Lastly, our study illustrated that miR-194-3p is a very promising therapeutic target of bronchial epithelial cells. Next, transfected animal studies should be prepared to find a reasonable dose and measurement of treatment. It is necessary to find the most appropriate way to deliver miR-194-3p mimics to bronchial epithelial cells and to reduce the nonspecific toxicity at the same time.

Conclusion
Our experiments proved that PM 2.5 downregulates miR-194-3p and accelerates apoptosis in cigarette-inflamed bronchial epithelium by directly targeting DAPK1. Reduction of miR-194-3p upregulated expression of DAPK1 and cleavage of caspase 3. Overexpression of miR-194-3p could suppress apoptosis in PM 2.5 -cigarette-exposed bronchial epithelial cells. These findings suggest that miR-194-3p could be a potential therapeutic target for the treatment of PM 2.5 -induced bronchial epithelial injury aggravation.

Publish your work in this journal
Submit your manuscript here: http://www.dovepress.com/international-journal-of-chronic-obstructive-pulmonary-disease-journal The International Journal of COPD is an international, peer-reviewed journal of therapeutics and pharmacology focusing on concise rapid reporting of clinical studies and reviews in COPD. Special focus is given to the pathophysiological processes underlying the disease, intervention programs, patient focused education, and self management protocols.
This journal is indexed on PubMed Central, MedLine and CAS. The manuscript management system is completely online and includes a very quick and fair peer-review system, which is all easy to use. Visit http://www.dovepress.com/testimonials.php to read real quotes from published authors.