Timescale Separation of Positive and Negative Signaling Creates History-Dependent Responses to IgE Receptor Stimulation

The high-affinity receptor for IgE expressed on the surface of mast cells and basophils interacts with antigens, via bound IgE antibody, and triggers secretion of inflammatory mediators that contribute to allergic reactions. To understand how past inputs (memory) influence future inflammatory responses in mast cells, a microfluidic device was used to precisely control exposure of cells to alternating stimulatory and non-stimulatory inputs. We determined that the response to subsequent stimulation depends on the interval of signaling quiescence. For shorter intervals of signaling quiescence, the second response is blunted relative to the first response, whereas longer intervals of quiescence induce an enhanced second response. Through an iterative process of computational modeling and experimental tests, we found that these memory-like phenomena arise from a confluence of rapid, short-lived positive signals driven by the protein tyrosine kinase Syk; slow, long-lived negative signals driven by the lipid phosphatase Ship1; and slower degradation of Ship1 co-factors. This work advances our understanding of mast cell signaling and represents a generalizable approach for investigating the dynamics of signaling systems.

: DNP-lysine does not induce degranulation. Cells were first exposed to a 5-min pulse of DNP-lysine and then a 5-min pulse of DF3, and degranulation was measured during each of these periods. DNP-lysine did not induce appreciable degranulation. Figure 2: 5 min of DNP-lysine is sufficient to break up receptor aggregates and reduce degranulation. Degranulation was monitored during an initial 5 min pulse of DF3, and then during a pulse of DNP-lysine that lasted either 5 or 10 min. Note that 5 min and 10 min of DNP-lysine exposure result in similar degranulation levels. The impact of S1 duration on responses during S2. The duration of S1 was varied from 30 s to 5 min. The durations of I and S2 were both 5 min in all cases. Desensitization occurs to a similar extent for all S1 durations: the response during S2 is alike for all of the S1 durations considered.

Supplementary
Supplementary Figure 4: Desensitization is nonspecific. Cells were sensitized with a half-and-half mixture of anti-DNP IgE, and anti-dansyl IgE. Cells were exposed to DF3 or a dansyl-BSA (DNS-BSA) in the indicated order. For all of these experiments, S1 = 5 min, I = 5 min, and S2 = 5 min. The blue bar represents the response to the first pulse and the red bar represents the response to the second pulse. Exposure to one antigen promotes desensitization to the other. were incubated with 10 nM DF3 for 10 min. B) Cells were incubated with 10 nM DF3 for the time shown in red (pulse), then DF3 was removed and cells were incubated for the time indicated in blue (chase) in media without DF3 for a total of 10 min. After each pulse with DF3 (10, 8, 5, or 2 min), supernatant was collected to measure degranulation prior to chase (red diamond, average degranulation values indicated above), and after the chase (blue diamond, average degranulation values indicated above). Degranulation increases for continuous stimulation up to 10 min. The increase in degranulation over this time period indicates that desensitization is not the result of a shortage of granules. , or 20 µM SHIP1 inhibitor, 3AC, and then exposed to a stimulation pattern of S1 = 5 min, I = 5 min, and S2 = 5 min. The degranulation level of S2 relative to that of S1 for cells treated with 1 µM NSC, 0.18 ± 0.05, was similar to that of untreated cells (S2 relative to S1, 0.11 ± 0.03). In contrast, for cells treated with 3AC, there was clear desensitization with a much higher S2 relative to S1 (0.78 ± 0.04).

Overview of the Model and Simulation Approach
We used the BioNetGen language (BNGL) to formulate a mathematical model for IgE receptor (Fc RI) signaling with the goal of representing the opposing effects of Syk and Ship1 on antigen-stimulated secretory responses of antigen-sensitized mast cells. Syk is a protein tyrosine kinase recruited to activated Fc RI that plays an important role in generating positive signals for mast cell degranulation, whereas Ship1 is a lipid phosphatase that plays an important role in generating negative signals. BNGL is a formal language for defining models (1; 2), which is akin to a domain-specific programming language. BNGL can be used to define not only models but also simulation protocols and simulation outputs. BNGL is compatible with the BioNetGen software package (1; 3), which includes various tools for deterministic and stochastic simulation of well-mixed chemical kinetics. In this study, we used the deterministic simulation capabilities of BioNetGen. The model is available as a plain-text BioNetGen input file named model.bngl. A listing of this file is provided in the Appendix. An electronic version is available as online supplementary material. The file model.bngl consists of declarations that introduce 1) molecule types (e.g., formal representations of Syk and Ship1), 2) reaction rules and associated rate laws for molecular interactions and biochemical processes (e.g., a rule for recruitment of Syk to phosphorylated Fc RI), 3) parameters (e.g., rate constants) and parameter settings, and 4) a list of seed species and their concentrations, which we took to define initial conditions for an initial value problem. The file model.bngl also includes a definition of simulation outputs, called observables, and simulation protocols/commands, called actions, such as generate_network and simulate.
Parameter fitting was performed using BioNetFit (4), a fitting program designed to work with BioNetGen. BioNetFit uses a genetic algorithm to stochastically find a best-fit parameter set for a model. BioNetFit was run on the plain text input file model_fit.conf, which configures BioNetFit with the appropriate model file (model_tofit.bngl) and experimental data files (p1_5.exp, p3_5.exp, p3_30.exp, p3_60.exp, p1_120.exp, p1_240.exp), as well as settings for fitting. As output from BioNetFit, we obtain a modified version of model_tofit.bngl, with the best-fit parameter values added to the beginning of the file.
The file model.bngl was processed by BioNetGen (1) to find the individual chemical species and reactions implied by the specified rules and seed species, as well as the corresponding set of coupled ordinary differential equations (ODEs) for well-mixed chemical kinetics. These ODEs were then integrated numerically. Numerical integration was performed by using CVODE in the SUNDIALS software package (5) with BioNetGen's default algorithmic settings, which are appropriate for stiff problems.
Below, we provide further information about the model and the computational methods used to analyze the model.

Molecule types
In the rule-based formulation of the model, we consider the following seven molecule types: The functional components of these molecules (i.e., their sites) are annotated in Table S1. Sites Yb and Yg in R, s in X, and loc in H are taken to have internal states. The internal states of Yb and Yg indicate if the sites are phosphorylated (P) or unmodified (0). The internal state of s indicates if X is competent (on) or incompetent for binding to Ship1 (off). The internal state of loc indicates the location of H: inside granules (in) or outside the cell as a consequence of secretion (out). Pairs of cognate binding sites are identified in Table S2. These pairs of sites identify the four binary interactions considered in the model.

Rules
In the model, molecular interactions and other processes (e.g., proteasomal degradation of X) are represented by the following 14 rules, which are associated with 17 rate constants. For reversible rules, we define the forward rate constant, followed by the reverse rate constant.
kpten k pten Basal clearance of PIP3 X(s∼on)->0 kdegX k degX Degradation of active free X by the proteasome X(s∼on!1).Ship1(x!1)-> kdegX k degX Degradation of active Ship1-bound X Ship1(x) by the proteasome PIP3()+H(loc∼in)-> kdegran k degran Degranulation caused by the PIP3()+H(loc∼out) presence of PIP3 At this point, we should make note of the following simplifying assumptions. The receptor is taken to be monoclonal anti-DNP IgE in complex with Fc RI, the high-affinity IgE receptor found on mast cells. The IgE-Fc RI complex is long lived, on the order of hours to days. In the model, we represent the abundance of trivalent DNP ligand DF3 with an effective concentration of isolated DNP sites. When a receptor is not bound to a DNP site, we assume that it is unaggregated, and therefore, not competent to undergo receptor phosphorylation. In contrast, when a receptor is bound to a DNP site, we assume that it is aggregated and, accordingly, competent for phosphorylation. Receptor phosphorylation depends on receptor aggregation and is mediated predominantly by the protein tyrosine kinase Lyn, which we consider implicitly. We further assume that phosphorylation simultaneously occurs at the Ship1 binding site on the β chain and the Syk binding sites on the γ chain. Finally, we assume that the concentration of PIP3 directly relates to the rate of degranulation. Downstream events, including the conversion of PIP3 to inositol trisphosphate (IP3), binding of IP3 to the IP3 receptor, and the subsequent mobilization of calcium, are considered implicitly.

Initial Conditions
Simulations began with initial quantities of each inactive molecule type (called seed species in BNGL). R, Syk, and Ship1 were each initiated at 3e5 copies per cell. H(loc∼in) had an initial count of 1e6. PIP3 had an initial count of 0. The initial copy number of X was left as a free parameter, and was estimated through fitting to have a value of 1.84e6 copies per cell.

Simulation Protocol
Simulated time courses were run with conditions corresponding to experimental data. In the first phase of the simulation, the count of Ag was set to 3e6 copies, corresponding to a concentration of 10 nM in an extracellular volume of 0.5 nL per cell. This count of 3e6 was held constant even as Ag bound and unbound R (because of the approximately constant bulk concentration of antigen in the flow channel of the microfluidic device). The first phase of simulation was run for 5 min. The degranulation was reported as the final abundance of H(loc∼out), as a percentage of the total 1e6 copies of H per cell. In the second phase of simulation, Ag count was set to 0, to represent exchange of DF3 for excess DNP-lysine in the flow channel. This was run for a variable time interval of 5 to 240 min. Before the final phase of simulation, H(loc∼out) was reset to 0, to measure only degranulation occurring in the final phase. Then Ag was reset to 3e6 copies, the simulation was run for another 5 min, and the final H(loc∼out) count was recorded.
To simulate proteasome inhibition, we disabled the two rules corresponding to degradation of active free and Ship1-bound X. To simulate inhibition of Ship1 by 3AC, we multiplied the rate k deg by a factor ranging from 0.1 to 0.9 (corresponding to varying possible levels of Ship1 inhibitor efficacy). To simulate a Shc1 knockdown, assuming X represents Shc1, we multiplied the initial concentration of X by 0.32, corresponding to the knockdown efficiency shown in Fig. 6C.

Parameters
We used BioNetFit to fit the 17 rate constants in Table S3, as well as the initial copy number of X. Parameters were fit to the experimental, on-device degranulation data shown in Fig. 4, which was tabulated in BioNetFit-readable format in the plain text files p1_5.exp, p3_5.exp, p3_30.exp, p3_60.exp, p1_120.exp, and p1_240.exp. Fitting was performed using a genetic algorithm with 300 generations, with 80 parameter sets per generation. The initial population of parameter sets was chosen according to a log uniform random distribution spanning 4 orders of magnitude around an initial guess for each parameter (see model_fit.conf for specific values). Simulation results for each parameter set were compared to the data with a chi-squared objective function.
In an initial fitting run, rate constants k on and k deg were set to values larger than a reasonable upper limit of 1e7 /M/s. These parameters were then set equal to 1e7 /M/s, and the fitting was repeated with the remaining 16 parameters free to vary.
The final fit parameters are given in the file model.bngl.

Bayesian Parameter Analysis
To assess our degree of uncertainty in the parameter values found with BioNetFit, we performed Bayesian Markov Chain Monte Carlo (MCMC) simulations to sample the credible region of parameter space. We followed the method described in (6). Briefly, for a parameter set Θ used to predict experimental data y, we seek to sample the posterior probability distribution P (Θ|y). By Bayes' theorem, For the log likelihood ln P (y|Θ), we use the negative of a χ 2 objective function, and for the prior P (Θ), we use a log (base 10) normal distribution, with mean equal to our value from BioNetFit and standard deviation of 1. P (y) can be disregarded, as it remains constant for all Θ. Thus we can compute P (Θ|y) by evaluating the prior and objective function for any parameter set Θ.
To perform a step of the MCMC simulation, we propose a new parameter set Θ * (by perturbing Θ by a distance of 0.1 in a random direction in log (base 10) space), and compute P (Θ * |y). The proposed move is accepted with probability min(1, P (Θ * |y)/P (Θ|y)). Our MCMC simulations consisted of 30 independent runs of 200000 steps each. Θ was recorded every 20 steps after an initial 20000-step burn-in period.
From this analysis, we obtain the marginal probability distribution of each parameter value (Fig. S10). This distribution indicates the range of reasonable values for each parameter; a narrower distribution indicates that the model is more sensitive to the value of that parameter.
Additionally, we assessed the robustness of the predictions of our model. For this analysis, we arbitrarily chose 540 parameter sets that were sampled during the MCMC simulation. We consider the predicted degranulation data, as shown in Fig. 4, for each of the sampled parameter sets. At each time point, the middle 68% of the predictions gives our 68% credible interval. In Fig. S11, we present this credible interval for each time interval tested. The credible interval forms a narrow band around each predicted value, which indicates that the degranulation predictions are quite consistent throughout the sampled parameter space.
Similarly, we evaluated the Syk and Ship1 activation time courses, as shown in Fig. 5a-b. At each time point, we considered the Syk or Ship1 activation level with each of our 540 sampled parameter sets, and took the middle 68% of those values to define the 68% credible interval at that time point. The credible intervals for the entire time courses are shown in Fig.  S12. We find that these credible intervals preserve the qualitative features of the time course, indicating that this qualitative behavior is robust. Figure S10: Marginal probability distributions for each parameter of the model. X-axes indicate the log base 10 of the parameter values. Red lines indicate the parameter value found by BioNetFit, which was used in the model for the rest of the study, and served as a prior for the Bayesian algorithm. Units are s -1 for first order rate constants, and M -1 s -1 for second order rate constants. X tot is expressed as the ratio of X concentration to Ship1 concentration.

of 25
Harmon B, Chylek LA, Liu Y et. al The chemical reaction network implied by rules and seed species Table S4 lists the 23 individual chemical species that can be populated according to the rules of the model. Table S4. Abbreviated names of chemical species. Symbol Representation in BNGL Network generation by BioNetGen yields the following list of 23 differential equations that define the model. These equations were numerically integrated to obtain simulation results.
+k mShip x 19 x 11 = +k of f x 10 − k on x 1 x 11 − k pase x 11 − k pSyk x 11 x 3 −k pShip x 11 x 4 + k mSyk x 15 − k pShip x 11 x 18 + k mShip x 16 + k mShip x 21 x 12 = +k pSyk x 10 x 3 − k of f x 12 − k mSyk x 12 − k pShip x 12 x 4 +k on x 1 x 15 − k pShip x 12 x 18 + k mShip x 17 + k mShip x 22 x 13 = +k pShip x 10 x 4 − k of f x 13 − k pSyk x 13 x 3 − k mShip x 13 −k pX x 14 x 13 + k on x 1 x 16 + k mSyk x 17 + k mX x 19 + k degX x 19 x 14 = +k Xon x 10 x 7 + k Xon x 11 x 7 + k Xon x 12 x 7 + k Xon Molecules Xall X() # total abundance of hypothetical Ship1 cofactor X Molecules X_on_free X(s∼on) # abundance of free X in activated state Molecules X_on_free_or_bound X(s∼on!?) # abundance of X (bound or free) in # activated state Molecules XShip1 X(s∼on!1).Ship1(x!1) # abundance of Ship1 bound to (activated) # cofactor X end observables begin reaction rules # ligand-receptor binding # As a simplification, we consider a one-step binding mechanism. # Thus, the rate constants are effective parameters that reflect # capture and release of (multivalent) antigen as well as the effects # of antigen-mediated receptor aggregation on the residence time of antigen