Mathematical Modeling and Parameter Estimation of Intracellular Signaling Pathway : Application to LPS-induced NF κ B Activation and TNF α Production in Macrophages

Due to the intrinsic stochasticity, the signaling dynamics in a clonal population of cells exhibit cell-to-cell variability at the single-cell level, which is distinct from the population-average dynamics. Frequently, flow cytometry is widely used to acquire the single-cell level measurements by blocking cytokine secretion with reagents such as GolgiplugTM. However, GolgiplugTM can alter the signaling dynamics, causing measurements to be misleading. Hence, we developed a mathematical model to infer the average single-cell dynamics based on the flow cytometry measurements in the presence of GolgiplugTM with lipopolysaccharide (LPS)-induced NFκB signaling as an example. First, a mathematical model was developed based on the prior knowledge. Then, average single-cell dynamics of two key molecules (TNFα and IκBα) in the NFκB signaling pathway were measured through flow cytometry in the presence of GolgiplugTM to validate the model and maximize its prediction accuracy. Specifically, a parameter selection and estimation scheme selected key model parameters and estimated their values. Unsatisfactory results from the parameter estimation guided subsequent experiments and appropriate model improvements, and the refined model was calibrated again through the parameter estimation. The inferred model was able to make predictions that were consistent with the experimental measurements, which will be used to construct a semi-stochastic model in the future.


Introduction
To integrate of multiple signaling pathways, their canonical transcription factors and downstream effector genes is required for cells to respond to various signals they encounter in their micro-environment.Therefore, understanding how information is sensed and processed by cells and the signaling pathways that are engaged by different stimuli can help elucidate cellular behaviors and responses.Typically, cellular signal dynamics and the response to stimuli have been studied using a combination of mathematical modeling and experimental analysis [1,2].A majority of these studies has modeled cell signaling at the population level and used population-averaged measurements such as Western blots to infer the dynamics of different proteins in the signaling pathway, as well as the possible network structure of signaling pathways [1].However, with recent advances in the ability to measure gene and protein expression at the single-cell level (reviewed in [2,3]), it has become possible to analyze signaling dynamics at the single-cell level.In contrast to the observations from population-average studies, the single-cell studies have demonstrated that individual cells in a clonal population may respond differently to the same stimulus, and the population level measurements could mask the temporal dynamics of individual cells [2].This variability in the responses of individual cells poses a challenge to their implementation in biology and medicine [4].Therefore, it is important to understand the stochasticity and heterogeneity in the single-cell responses that might be missed in population-averaged measurements.
Advances in experimental tools for single-cell analysis have led to a significant increase in single-cell studies [2,3].Despite these advancements, it is still difficult to study the single-cell signaling dynamics due to complex interactions at multiple levels between different proteins that are involved in signal transduction [1].Computational modeling has been proposed as a complementary approach to overcome some of these limitations and gain insights that cannot be obtained solely through experiments [1,2].A viable and computationally efficient approach to study the cell-to-cell variability is to use a deterministic model with parameters that have distributions [2,[5][6][7].In this approach, the computational cost is generally reduced by simulating the signaling dynamics through a deterministic modeling approach while the stochasticity is preserved by assigning a set of different parameter values for each simulation based on predetermined parameter distributions.
In order to construct such models, an experimentally validated deterministic model, which can capture average signaling dynamics at the single-cell level, is required.Although various deterministic models have been proposed for several well-studied signal transduction pathways [1,8], many demonstrate good qualitative, but not quantitative, agreement with the experimental data.This has been attributed to, among other factors, the limited breadth of data used for training the model (e.g., models trained using one dataset with a single stimulus concentration), which makes the models unable to make robust predictions under different conditions.Moreover, the identifiability issue of model parameters [9], which arises due to the model structure as well as the limited availability of experimental data of intracellular proteins [10,11], is not always addressed, which may lead to a suboptimal estimation of model parameters [10,12].Additionally, many models have been constructed and validated based on experimental data obtained from the population-averaged measurements, which mask the signaling dynamics at the single-cell level [2,13,14].Consequently, these models are inadequate to predict the average signaling dynamics of single cells.
Motivated by the above considerations, we developed a deterministic model that can accurately predict the average signaling dynamics of single cells.We chose lipopolysaccharide (LPS)-induced nuclear factor κB (NFκB) signaling in mouse macrophages for our model system as it is an extensively studied and characterized signaling pathway [8,15,16].In order to address the issues discussed above, both computational and experimental approaches have been implemented.First, a rigorous numerical scheme is used to identify the most important parameters that are to be estimated in the parameter estimation [17].Specifically, the sensitivity analysis and the parameter selection method quantitatively assess the significance of each model parameter with respect to experimental measurements under different LPS concentrations and select parameters whose values could be uniquely estimated [10,18].Second, flow cytometry with intracellular staining is used to measure the average single-cell dynamics of key molecules involved in the NFκB signaling pathway in response to a broad range of LPS concentrations [19,20].In this study, the intracellular concentrations of the inhibitor of κB-α (IκBα) and tumor necrosis factor α (TNFα) were measured.IκBα is an inhibitor of NFκB activity, and therefore, the IκBα dynamics are inversely correlated with the NFκB dynamics.At the same time, the activated NFκB induces the transcription and translation of TNFα upon the stimulation of LPS; hence, the TNFα can also be used to infer the dynamics of the NFκB signaling pathway [16].The obtained average single-cell kinetics is used to quantitatively calibrate and validate the model.Third, the discrepancy between the experimental measurements and the model predictions reveals important, yet unconsidered mechanisms, which is validated experimentally afterwards and leads to the model refinement.Through this integrated model development methodology, predictions from the resultant model quantitatively agree with the experimental measurements.Therefore, the proposed model represents a first step towards the construction of single-cell semi-stochastic models to investigate the stochasticity of intracellular NFκB signaling in macrophages.

Flow Cytometry Analysis
The expression of TNFα and IκBα under different experimental conditions was determined using flow cytometry.RAW264.7 cells were seeded into round-bottomed 96-well plate and stimulated with different concentrations of LPS for the indicated time.Golgiplug ™ (BD Biosciences, San Jose, CA, USA) was added along with LPS for TNFα detection experiments to block secretion of TNFα.Cells were then stained with Alexa Flour 700 fluorescence-tagged TNFα antibody (BD Biosciences) and PE-conjugated IκBα antibody (Cell Signaling Technology, Danvers, MA, USA) using the manufacturer's suggested protocol.Stained cells were analyzed using a BD Fortessa flow cytometer (BD Biosciences) at the Texas A&M Health Science Center College of Medicine Cell Analysis Facility.Ten thousands events per sample were acquired, and the data were analyzed using FlowJo software (Tree Star, OR, USA).Cells were gated based on side scattered light (SSC) and forward scattered light (FSC) values to eliminate cell debris, and TNFα-and IκBα-positive cells were gated based on the antibody isotype (see Supplementary Materials Figures S1-S3).All experiments were repeated using at least three different cultures.

Model Development
The schematic diagram of the NFκB signaling pathway is illustrated in Figure 1.The model used in this study was adopted from Caldwell et al. [21], which takes the extracellular LPS concentration as an input to predict the kinetics of key biomolecules in the NFκB signaling pathway.In this model, by forming a complex with Toll-like receptor 4 (TLR4), LPS activates IκB kinase (IKK) through myeloid differentiation primary response 88 (MyD88)-or TIR (Toll/Interleukin-1 receptor)-domain-containing adaptor-inducing interferon-β (TRIF)-dependent activation of TNF receptor-associated factor 6 (TRAF6).The activated IKK in turn promotes the translocation of NFκB to the nucleus, where the nuclear NFκB induces the transcription of NFκB inhibitors (IκB-α, -β, -and A20), as well as TNFα.Once translated, these inhibitors inhibit the NFκB signaling pathway.In contrast, the translated TNFα is secreted to the extracellular medium, and some of the secreted TNFα proteins will bind with TNFα receptor (TNFR) on the cellular membrane to initiate the TNFα-induced NFκB signaling pathway (see [21][22][23] for details of the model).
Additionally, nonlinear functions proposed by Junkin et al. [24] were added to describe how the rates of TNFα production and secretion increase as the amount of activated TRIF complex increases.This model incorporates the TLR4-mediated NFκB dynamics induced by LPS, as well as the production of TNFα in macrophages (see [21,23] for details).For the purpose of this study, two modifications were made to the model presented by Caldwell et al. [21].First, transcription delays were ignored to facilitate the simplicity of subsequent calculations for sensitivity analysis and parameter estimation.Second, a new role of A20 protein, which was introduced in the previous model [21,23] as an inhibitor of the TNFα-induced NFκB signaling [23,25,26], was included in the modified model to downregulate the LPS-induced signaling through deubiquitinating of TRAF6 [27].For this study, the TNFα production at the single-cell level was measured using flow cytometry by adding Golgiplug ™ since brefeldin A, the active agent of Golgiplug ™ , causes the Golgi apparatus to merge with endoplasmic reticulum (ER) and inhibits protein export from the Golgi complex [28,29].Hence, the addition of Golgiplug ™ enabled us to measure average single-cell production of TNFα.On the other hand, because Golgiplug ™ interferes with the normal cellular processes, it inevitably affects the NFκB signaling dynamics.Specifically, Golgiplug ™ suppresses the expression of receptors on the cellular membrane, which negatively regulates the LPS-mediated NFκB signaling pathway in different ways.First, the addition of Golgiplug ™ can block the translocation of TLR4 and its accessory molecules from the Golgi complex, which leads to the termination of signaling as these receptors are not replenished after turnover [28,[30][31][32].Similarly, TNFR is also depleted from the cellular membrane due to Golgiplug ™ [33,34], which may inhibit subsequent TNFα autocrine and paracrine signaling [35][36][37].Second, Golgiplug ™ can hinder the membrane expression of the cluster of differentiation 14 (CD14), which regulates the endocytosis of LPS or the TLR4-LPS complex [38][39][40].Therefore, the TRIF-dependent pathway, which is initiated only after LPS or LPS-TLR4 is endocytosed into cytoplasm [5,41], can also be partially impaired.Lastly, the secretion of TNFα proteins translated in response to the NFκB activation will also be inhibited, which helps measure the TNFα production at the single-cell level.
Consequently, the dynamic effects of Golgiplug ™ were parameterized and included in the model by the following equations: where G is the normalized activity of Golgiplug ™ , t is the elapsed time from the addition of Golgiplug ™ , τ is the characteristic time associated with Golgiplug ™ activity, k s TNFR and k s TLR4 are the constitutive synthesis rates of TNFR and TLR4, respectively, in the absence of Golgiplug ™ , k en LPS and k en cp are the endocytosis rates of LPS and the LPS-TLR4 complex, respectively, in the absence of Golgiplug ™ , k sec is the TNFα secretion rate in the absence of Golgiplug ™ and k s TNFR ,m , k s TLR4 ,m , k en LPS ,m , k en cp ,m and k sec,m are the corresponding rates in the presence of Golgiplug ™ .After Golgiplug ™ is added to the cells at t = 0, G slowly increases from zero to one, which corresponds to no inhibition of protein export to complete inhibition of protein export from the Golgi apparatus in the presence of Golgiplug ™ .
Since the signaling kinetics under the stimulation of LPS in the presence of Golgiplug ™ were measured experimentally, the dynamic model that consists of the model presented in [21] and Equation ( 1) was used to simulate the dynamics of LPS-induced NFκB signaling in the presence of Golgiplug ™ .In general, the dynamic model that simulates the signaling pathway can be represented by a set of nonlinear ordinary differential equations as follows: where x represents the concentration of the biomolecules involved in the signaling pathway (i.e., a vector of states), θ is a vector of model parameters that describe the biochemical reaction rates in the process, u is the concentration of LPS added to the cells (i.e., the process input), and y is the model output (i.e., the experimental measurements predicted by the model).When Golgiplug ™ is added, Equation ( 1) is included in Equation ( 2), and the overall model consists of 49 states and 146 parameters (see Supplementary Materials Tables S1-S2 and Equations (S1)-(S60)).

Parameter Estimation
Since we added the Golgiplug ™ module to the model developed by Caldwell et al. [21], the integrated dynamic model (the model presented in [21] and Equation ( 1)) was quantitatively calibrated by estimating its parameters using experimental measurements in response to different LPS concentrations in the presence of Golgiplug ™ .
The model parameter values were estimated by minimizing the difference between the experimental measurements and the model predictions of the protein concentration.In this work, we used flow cytometry to measure two key molecules in the LPS-induced NFκB signaling pathway: TNF α and IκBα.Since flow cytometry does not provide direct measurements of protein concentration, the mean fluorescence intensity (MFI), which is a measure of the number of copies of the target molecule per cell, was used to infer the protein concentration by assuming a linear relationship between MFI and protein concentration.The experimental data and model prediction were compared based on fold changes of MFI, which are defined as follows: where y IκBα (t) and y TNFα (t) are the fold changes of the IκBα and TNFα concentration at time t, x IκBα , x IκBα n , x NFκB-IκBα , x NFκB-IκBα n and x TNFα are the cytoplasmic IκBα, nuclear IκBα, cytoplasmic IκBα-NFκB complex, nuclear IκBα-NFκB complex and intracellular TNFα concentration, respectively, x i,0 is the initial concentration of the corresponding biomolecules, I IκBα and I TNFα are the MFI of IκBα and intracellular TNFα, respectively, and I j,0 and I j,c , ∀ j = {IκBα, TNFα}, are the corresponding MFI at t = 0 and MFI of negative control, respectively.In each cell, IκBα can be part of four biomolecules (x IκBα , x IκBα n , x NFκB-IκBα , x NFκB-IκBα n ); however, flow cytometry measurements can only provide the total IκBα concentration in each cell.Therefore, the simulated concentrations of four IκBα-containing biomolecules were initially summed, and the fold change of the sum (i.e., y IκBα ) was computed to compare with the measurements in the subsequent parameter estimation procedure.One of the biggest challenges in estimating parameters of signaling pathways with a large number of parameters is the parameter identifiability issue [10].That is, the exact values of some model parameters cannot be uniquely determined from experimental measurements even if a large amount of experimental measurements are available [10,11].As the proposed model has a large number of parameters, not all the model parameters can be estimated.To this end, a subset of the model parameters, which can be uniquely estimated from the available experimental measurements, was identified through a parameter selection method [10,18].Only these parameters were estimated against the experimental data.
First, local sensitivity analysis [10,42] was performed to compute two different sensitivity matrices S 1 and S 2 to quantify the effect of each model parameter on y IκBα and y TNFα (i.e., the process outputs).S 1 and S 2 represent the sensitivity matrices of the model parameters with respect to y IκBα and y TNFα , respectively, when the cells were stimulated with LPS in the presence of Golgiplug ™ .Specifically, a sensitivity matrix is defined as: where n p is the number of parameters in θ in Equation ( 2), and ∂y i (t l )/∂θ j quantifies the effect of a parameter θ j on an output y i at t = t l , ∀l = 1, • • • , N t , where N t is the number of measurement instants.∂y i (t l )/∂θ j can be computed by the following equation: Additionally, the term ∂x/∂θ j in Equation ( 5) can be computed by integrating the following equation along with Equation (2): Second, the Gram-Schmidt orthogonalization method [10,18] was used to identify the p i most important model parameters to be estimated for each S i , ∀ i = 1, 2. Here, p i is the number of singular values of S i whose magnitudes are at least 5% of the largest singular value [17,18].As a result, the parameter subset to be estimated, θ s ∈ R p×1 where p ≤ p 1 + p 2 , is chosen as the union of the selected parameters from S 1 and S 2 .Third, the least-squares problem was solved to estimate the values of θ s by minimizing the difference between the model predictions and the experimental data of y TNFα and y IκBα while the values of the unselected parameters were fixed at their nominal values selected from the literature [21,24,43,44] with some modifications.
In this study, three LPS concentrations (10, 50 and 250 ng/mL) were used to stimulate cells, and the MFI of IκBα and TNFα were measured at t = 0, 10, 20, 30, 60, 120, 240 and 360 min after the addition of LPS with Golgiplug ™ (i.e., t l , ∀ l = 1, • • • , 7).Specifically, the MFI data from 10 and 250 ng/mL of LPS (i.e., u k , ∀ k = 1, 2) were used to estimate the parameter values, while the dataset from 50 ng/mL LPS was used to validate the model with the updated parameters.Then, the least-squares problem is formulated as follows: x lb ≤ x k,i ≤ x ub (10) where y IκBα,k,1 (t l ) and y TNFα,k,1 (t l ) are the simulated fold changes of IκBα and TNFα, respectively, through Equation ( 9) at t = t l under the initial LPS concentration of u k in the presence of Golgiplug ™ , ŷIκBα,k,1 and ŷTNFα,k,1 are the corresponding experimentally measured fold changes and x 0 is the vector of the initial conditions of x (see Supplementary Materials Table S1).
In the least-squares problem of Equations ( 7)-( 11), the objective function of Equation ( 7) computes the difference between model predictions and the experimental measurements of the proteins in the presence of Golgiplug ™ .As a whole, the objective function minimizes the difference by varying the values of θ s .While Equation ( 8) is integrated to compute the predicted protein concentration x k,i , f 1 , which includes Equation (1), is used if Golgiplug ™ is present; otherwise, f 2 , which does not involve Equation ( 1), is integrated.The initial condition of the model, x0 , is assumed based on a previous study [21,23].Equations ( 10)- (11) impose lower and upper bounds on the states and parameters, respectively, based on previous studies and underlying biological knowledge [5,21,23].
It should be noted that we preserved one set of the experimental measurements (one obtained under 50 ng/mL of LPS) to validate the parameter estimation results [45].As Equations ( 7)-( 11) are likely to be non-convex, the choice of the initial guesses is important.In this study, the initial guesses for the above least-squares problem were obtained from Caldwell et al. [21], which were validated experimentally by comparing with the population-level measurements.Therefore, the parameter values estimated by Caldwell et al. [21] were suitable initial guesses.At the same time, Equations ( 7)- (11) were solved multiple times with different initial values to avoid any suboptimal optima.Model simulations and the parameter estimation were performed in MATLAB via its functions ode15s and fmincon.The absolute and relative tolerance criterion for ode15s were set as 10 −9 , and fmincon was implemented with multistart to obtain a better result by solving Equations ( 7)-( 11) multiple times with different initial conditions.

Results
Profiles of de novo synthesized intracellular TNFα under the stimulation of LPS in the presence of Golgiplug ™ demonstrated that the TNFα production increased around one hour after the stimulation (Figure 2).At around the same time, the IκBα concentration reached its minimum, which is consistent with experimental observations in the literature [46][47][48].Subsequently, the IκBα concentration increased due to the induction of IκB transcript (IκB t ) by nuclear translation of NFκB, while the TNFα production rate slowed down beyond 4 h of LPS stimulation (Figure 2).It should be noted that no experiments were conducted beyond 6 h after LPS was added to the cell culture based on the manufacturer's guideline on Golgiplug ™ use.This is most likely based on the fact that Golgiplug ™ might induce the apoptosis of cells exposed to it for a long time [49,50].As a result, the calibrated model is more suitable to describe the early NFκB signaling pathway (≤6 h) upon the LPS stimulation.

Model Validation
Based on the criteria outlined above in the previous section, six parameters (Table 1) were selected for parameter estimation.Figure 2 shows the simulated profiles of intracellular TNFα and IκBα in macrophages under the stimulation of LPS in the presence of Golgiplug ™ after the parameter estimation.While the model predictions agreed well with the experimental data obtained for 250 ng/mL of LPS, less concordance was observed between simulations and experimental data for 10 ng/mL of LPS.Specifically, the simulated concentration of intracellular TNFα was one order of magnitude lower than the MFI data, while the simulated IκBα dynamics were qualitatively similar to the measured MFI values.Since the discrepancy between the model prediction and the experimental measurements was pronounced with 10 ng/mL of LPS, we hypothesized that the lack of agreement between the simulations and experimental data was because the effects of Golgiplug ™ addition were not adequately represented in the model structure and were more pronounced at the lower LPS concentration.

Golgiplug ™ -Induced ER Stress
One possible explanation for this discrepancy could be that the addition of Golgiplug ™ induced other signaling pathways, which altered NFκB signaling dynamics [51].As Golgiplug ™ prevents protein secretion by causing collapse of the Golgi apparatus into the ER, synthesized proteins will be redistributed from the Golgi complex into the ER [29].A direct consequence of Golgiplug ™ addition could be accumulation of newly synthesized proteins in ER, which may induce ER stress [51].It is well established that the ER stress leads to phosphorylation of eukaryotic initiation factor 2 α-subunit (eIF2α), which partially inhibits the translation of IκB in the NFκB signaling pathway [51][52][53][54].This could lead to a decrease in the overall kinetics of the LPS-induced NFκB signaling as the concentration of IκB proteins would be kept lower, leading to the aforementioned mismatch between the model predictions and experimental data.Since the low LPS concentration induces less IκBα and its isomers (IκBβ and IκB ) (Figure 2), the entire LPS-induced NFκB pathway dynamics would be affected more significantly by Golgiplug ™ at a lower LPS concentration than at a high LPS concentration if the translation of IκBα and its isomers is partially inhibited.If this is true, it can lead to the pronounced disagreement between the model prediction and the experimental measurement under the stimulation of 10 ng/mL LPS as shown in Figure 2.
Therefore, we examined whether the Golgiplug ™ addition could modulate IκB levels in macrophages.First, the fold change in IκBα MFI with the stimulation of LPS alone, Golgiplug ™ alone, and LPS and Golgiplug ™ in macrophages were compared.Figure 3a shows that Golgiplug ™ alone lowers the concentration of IκBα, and Figure 3b-d show that the IκBα kinetics were altered when the cells were stimulated with LPS and Golgiplug ™ .While IκBα levels initially decreased when cells were exposed to LPS alone, they recovered to pre-stimulation levels after 3 h of exposure.However, when LPS was added along with Golgiplug ™ , IκBα levels continued to be lower than pre-stimulation levels (Figure 3b-d).These results suggested that Golgiplug ™ could affect the IκBα kinetics (presumably through the eIF2α phosphorylation) [52,54].This also explains the observations in Figure 2a-c, where the intracellular TNFα concentration continued to increase since the Golgiplug ™ -induced response prolonged the NFκB activation by inhibiting the IκBα synthesis.

Model Refinement
In order to account for the Golgiplug ™ -induced translation inhibition, the following equation was considered in addition to Equation (1): where k tl i ,m and k tl i are the IκBi, ∀ i = α, β, , translation rates in the presence and absence of Golgiplug ™ , respectively, ν is a coefficient for the maximum translation inhibition and K is the Michaelis constant of the eIF2α phosphorylation.Equation ( 12) was included in Equation ( 2) along with Equation (1) for an accurate simulation.The process affected by this translation inhibition is shown in Figure 1 via a red arrow.
The proposed dynamic model was calibrated again using the parameter estimation procedure as described above.Since the additional measurements of the IκBα dynamics in the absence of Golgiplug ™ were obtained, an extra sensitivity matrix was calculated, and the following was added to the objective function (7): where ŷIκBα,k,2 and y IκBα,k,2 are the simulated and measured IκBα fold change, respectively, in the absence of Golgiplug ™ .

Final Model Validation
Based on the updated model structure and the available experimental data, the aforementioned parameter selection approach determined eight parameters, which could be uniquely estimated (Table 2).Most of the parameters selected by the proposed parameter selection procedure were relevant to the core NFκB-IκB feedback system such as Hill coefficients for IκB-α and -transcription, IKK deactivation, and IκBα transcript degradation rate.The remaining identified parameters are the TLR4 constitutive generation rate, C1 (TNFR complex [23]) deactivation rate and eIF2α phosphorylation coefficient, which are most relevant to the LPS-and TNFα-induced NFκB activation, as well as the effects of Golgiplug ™ , respectively.Hence, all major processes considered in this system, which included the LPS-and TNFα-induced NFκB signaling pathway in the presence of Golgiplug ™ , were quantitatively validated against the single-cell experimental data.Figure 4 shows simulated fold changes in IκBα and intracellular TNFα after parameter estimation.The simulated profiles were again compared with the experimental data.The normalized root-mean-squares of the parameter estimation results before and after the incorporation of the Golgiplug ™ model (Equation ( 12)) were 3.8 and 2.5, respectively, which demonstrated the improvement of the model fidelity.Overall, the model predictions were in qualitative and quantitative agreement with both training datasets and validation datasets, as well as the literature data, which validated the prediction capability of the calibrated model, as well as our hypothesis on the effect of Golgiplug ™ on the inhibition of IκB translation.The results demonstrated that the calibrated model is capable of predicting input-output responses in the NFκB pathway.Additionally, the predictions from the current model were compared with the model proposed by Caldwell et al. [21] (Figure 4g-i).The proposed model was able to predict the observed IκBα dynamics under all LPS concentrations more accurately than the previous model by Caldwell et al. [21], which again demonstrated that the predictive capability of the model was improved in terms of simulating the IκBα dynamics.
In order to further assess the predictive capability of the newly calibrated model, the simulated dynamics of nuclear NFκB levels (i.e., activated NFκB) in the absence of Golgiplug ™ were computed and plotted in Figure 5a.The maximum NFκB translocation to the nucleus occurred within 2 h of LPS addition, which was consistent with previous experimental studies [15,55,56].Moreover, as the LPS concentration increased, the nuclear NFκB levels reached their maximum value earlier (i.e., at 50, 60, 75 and 105 min after adding LPS), and the areas under the curves in Figure 5a, which were computed as indicators of the signal strength, were around 20 µM•min for different concentrations of LPS.Interestingly, a 25-fold change in the LPS concentration only resulted in less than a 100% change in the signal strength.This observation was consistent with single-cell studies by Tay et al. [14,57], where they observed a relatively constant peak intensity and decreasing response time of the NFκB signal in mouse 3T3 cells upon TNFα or LPS stimulation.
Figure 5b shows the predicted amount of TNFα secreted under different LPS stimulation conditions in macrophages.As the LPS concentration increased, the concentration of secreted TNFα increased, which was expected since the signal (area under the peak) became stronger.Furthermore, similar to previous studies [15,21,58], the TNFα concentration peaked around 5 h after stimulation and gradually declined thereafter; however, the rate of decline was slower than that reported by Maiti et al. [15] (Figure 5c), where they measured the TNFα secretion dynamics from RAW264.7 macrophages in response to LPS stimulation at the population level.This observation was consistent with the observation reported by Xue et al. [13], who observed using human monocyte-derived macrophages the amount of TNFα secreted to the medium from a single cell in a cell population was less than that from an isolated single-cell at 20 h after the LPS stimulation.This suggested that the simulated dynamics by the proposed model is qualitatively similar to the signaling dynamics of an isolated single-cell instead of population-averaged dynamics, which was expected since the kinetic data obtained under Golgiplug ™ were used to train the model.[15] in response to 100 ng/mL of LPS.The TNFα concentration at each point was normalized to the maximum value obtained.

Discussion
In this study, we have developed a dynamic model that can accurately simulate the average single-cell dynamics of the NFκB signaling pathway by combining the single-cell measurements and a numerical scheme with sensitivity analysis, parameter selection and parameter estimation.The dynamic model was built based on a previously developed NFκB model [21,23,24] and calibrated using the experimental data and the aforementioned numerical scheme.Predictions from the developed dynamic model are in good agreement with the experimental measurements under all LPS concentrations, which demonstrates that the model is capable of simulating the average single-cell dynamics.
Previous studies have used stochastic simulation algorithms such as Gillespie's algorithm [59] and approximate methods of Gillespie's algorithm [60][61][62][63] to study single-cell dynamics and investigate heterogeneity in signaling pathways at the single-cell level [2,5,64].For example, Lipniacki et al. [14,64] proposed a hybrid stochastic-deterministic model of the TNFα-induced NFκB signaling pathway that was able to reproduce the heterogeneous responses observed in the single-cell measurements [14,65] and identify possible origins of the heterogeneity.However, stochastic simulation algorithms are computationally expensive, and they are difficult to fit to experimental measurements for model validation [7,66,67].A more viable method is a semi-stochastic model, which uses deterministic modeling with model parameters that have distributions [5][6][7], to reduce the computational cost while still studying the cell-to-cell variability.The dynamic model developed here can accurately simulate average single-cell dynamics and is a first step towards building a semi-stochastic model of the NFκB signaling.
The development of such a deterministic model for building a semi-stochastic model requires accurate parameter estimation, where values of model parameters are estimated by solving an optimization problem (Equations ( 7)-( 11)).However, parameter estimation is a nontrivial problem due to, but not limited to, ill conditioning, over-fitting and the non-identifiability of model parameters [9,68,69].The ill-conditioning and over-fitting problems during parameter estimation are attributed to the fact that available experimental measurements are usually very limited and noisy, while mathematical models of signaling pathways are often very comprehensive and include a large number of parameters [9,10].As a result, the solution to the parameter estimation problem is likely to be non-unique or very sensitive to noise present in the experimental measurements.Furthermore, even if a large number of noise-free experimental measurements are available, the value of a parameter cannot be uniquely determined if the parameter is not identifiable [10,11]; hence, it is necessary to check the parameter identifiability a priori.
The model developed in this work contains 148 parameters with limited experimental data, and hence, parameter estimation is very likely to suffer from the aforementioned issues in the parameter estimation procedure.Therefore, we implemented an integrated method combining sensitivity analysis and parameter selection before parameter estimation.Specifically, the sensitivity analysis quantified the effects of each parameter on the measurements, and the parameter selection method selected identifiable parameters via Gram-Schmidt orthogonalization.Then, the values of only the selected parameters were estimated in the parameter estimation, while the values of remaining parameters were fixed at their nominal values, which effectively alleviated the ill-conditioning problem by reducing the degrees of freedom in Equations ( 7)-( 11) [9,10,68].
After parameter estimation, the simulated profiles of intracellular TNFα and IκBα exhibited reasonable agreement between the model predictions and the experimental measurements at all LPS concentrations (Figure 4).Furthermore, as shown in Figure 5c, model predictions after parameter estimation were distinct from that of a cell population as the simulated profiles were closer to the signaling dynamics of isolated single-cells.This was likely because the use of Golgiplug ™ inhibited secretion of cytokines [70] and hence minimized potential autocrine and paracrine signaling from the secreted cytokines.This is important as the autocrine and paracrine signaling has been proposed as a key component in determining the overall signaling dynamics of cells in a population [13,35,36,71].Therefore, the proposed model, which was trained by the single-cell dynamics from flow cytometry in the presence of Golgiplug ™ , was able to describe the single-cell NFκB dynamics under minimal cytokine feedback.
It should be noted that the current model simulates the LPS-induced NFκB signaling dynamics in a cell, but it does not consider the initiation of the NFκB signaling pathway by TNFα secreted by neighboring cells.Hence, the flow cytometry measurements obtained in the presence of Golgiplug ™ are appropriate to identify realistic parameter values to reproduce average single-cell dynamics.At the same time, as flow cytometry measures cellular responses from thousands of cells simultaneously, flow cytometry can provide distributions of the measurements (see Supplementary Materials Figures S1-S3).Based on this statistical information, one can estimate the distributions of the parameters by different methods such as Bayesian approaches [6,7] or generalized polynomial chaos [72].The model with the estimated parameter distributions is then the semi-stochastic model that can be used to study the heterogeneity in cellular responses.
The present study also suggests that cytokine production data acquired using flow cytometry in the presence of Golgiplug ™ should be interpreted cautiously.As Golgiplug ™ can block cytokine secretion, it is often used to assess the cytokine production at the single-cell level using flow cytometry [73][74][75].The data shown in the work suggest that the dynamics of transcription factors and other signaling intermediates may be altered by the addition of Golgiplug ™ (Figure 3).Therefore, data from studies using Golgiplug ™ need to be interpreted cautiously, and a model-based approach like the one presented here can be useful in eliminating the effects of Golgiplug ™ and extract true signaling dynamics from flow cytometry data.

Conclusions
We systemically extracted the average single-cell dynamics of the LPS-induced NFκB signaling pathway through the integration of sensitivity analysis and a parameter selection scheme with flow cytometry data of key protein intermediates.Based on the measurements and the model structure, key model parameters were identified and estimated to maximize the prediction accuracy of the calibrated model while avoiding overfitting.The mismatch between the model predictions and experimental observations even after the parameter estimation revealed the existence of a previously unconsidered, yet important, mechanism related to Golgiplug ™ , which was subsequently validated by experiments and led to the update of the proposed model.Then, the resultant model was validated, and the simulated profiles from the updated model were in good agreement with experimental datasets under three different LPS concentrations.This model can be used as the nominal model to construct a deterministic model that has parameters with distributions and can be used to study the stochasticity in signaling.

Figure 1 .
Figure 1.Schematic diagram for the LPS-NFκB-TNFα signaling pathway.Due to space limitation, TRIF-dependent regulation of TNFα production, IκBβ and IκB -dependent NFκB deactivation and eIF2α-induced translation inhibition are not illustrated.Furthermore, some states related to TNFα-induced activation of IKK kinase (IKKK) are not shown.Colored arrows indicate the processes affected by the addition of Golgiplug ™ (see the text for details).

Figure 2 .
Figure 2. Parameter estimation before considering the Golgiplug ™ -induced ER stress.(a-c) Measured (empty circle) and simulated (solid line) fold changes of intracellular TNFα concentrations over time were plotted under different LPS concentrations in the presence of Golgiplug ™ .(d-f) Measured (empty circle) and simulated (solid line) fold changes of IκBα concentrations over time were plotted under different LPS concentrations in the presence of Golgiplug ™ .Indicated amounts of LPS were used for experiments and simulations.Experimental data are given as means ± SEM (standard error of means) with at least n = 6.

Table 1 .Figure 3 .
Figure 3. Kinetics of IκBα fold changes when the cells were stimulated by (a) 0, (b) 10, (c) 50, and (d) 250 ng/mL of LPS in the presence (empty circles) or absence (x marks) of Golgiplug ™ .Data are given as means ± SEM with at least n = 3.

Figure 4 .Figure 5 .
Figure 4. Parameter estimation considering Golgiplug ™ -induced ER stress.(a-c) Measured (empty circle) and simulated (solid line) fold changes of intracellular TNFα concentrations over time were plotted in the presence of Golgiplug ™ .(d-f)Measured (empty circle) and simulated (solid line) fold changes of IκBα concentrations over time were plotted in the presence of Golgiplug ™ .(g-i) Measured (empty circle) and simulated (solid line) fold changes of IκBα concentrations over time were plotted in the absence of Golgiplug ™ .The IκBα dynamics predicted by the model in[21,24] were also plotted in (g-i) for comparison.Indicated amounts of LPS were used for experiments and simulations.

Table 2 .
Selected parameters and their newly estimated values for the final model.