Stochastic dynamics of Type-I interferon responses

Interferon (IFN) activates the transcription of several hundred of IFN stimulated genes (ISGs) that constitute a highly effective antiviral defense program. Cell-to-cell variability in the induction of ISGs is well documented, but its source and effects are not completely understood. The molecular mechanisms behind this heterogeneity have been related to randomness in molecular events taking place during the JAK-STAT signaling pathway. Here, we study the sources of variability in the induction of the IFN-alpha response by using MxA and IFIT1 activation as read-out. To this end, we integrate time-resolved flow cytometry data and stochastic modeling of the JAK-STAT signaling pathway. The complexity of the IFN response was matched by fitting probability distributions to time-course flow cytometry snapshots. Both, experimental data and simulations confirmed that the MxA and IFIT1 induction circuits generate graded responses rather than all-or-none responses. Subsequently, we quantify the size of the intrinsic variability at different steps in the pathway. We found that stochastic effects are transiently strong during the ligand-receptor activation steps and the formation of the ISGF3 complex, but negligible for the final induction of the studied ISGs. We conclude that the JAK-STAT signaling pathway is a robust biological circuit that efficiently transmits information under stochastic environments.

Transcription.  The reaction rates are given in Table A.  Table A. Reaction rates considered in the model k 11.2 · ST AT 2 − IRF 9 n f 12. 1 k 12 · I irf 9 · ISGF 3 n f 12.2 k 12 · I irf 9 · ST AT 2 − IRF 9 n f 13. 1 k 13.1 · irf 9 f 13.2 k 13.2 · irf 9 * f 14.1 k 14 · I socs · ISGF 3 n f 14  represented by the following general scheme: for j = 1, ..., M , where v ij =v ij −v ij are the reaction stoichiometries. The stoichiome-34 try matrix v is composed of all the reaction stoichiometries.

35
Given the microscopic random processes that govern chemical reactions, it is possible 37 to describe the evolution of X(t) as a homogeneous Markov process in continuous time.

38
Considering this framework, the Chemical Master Equation (CME) describes the prob-39 ability that the system has a specific copy number of each S i at a given point in the 40 future: where the two terms within brackets give the rate at which the probability of being in 42 state X increases or decreases over time because of reactions into or out of state X, 43 respectively. a j indicates the probability that the reaction R j will occur in the next 44 infinitesimal interval [t, t + dt). weighted by the its own probability, that is: The second statistical moment is also known as the variance, and higher moments can 70 be computed, but this will not be discussed in the manuscript.

SI-10
C Parameter estimation strategy 72 The aim of the parameter estimation strategy is to find a unique parameter set that  Table 1, initial conditions are given in Table 3 of the main section.
where φ 1 is an scaling factor.

105
Given that different complexes are detected by the used antibody against phosphory-106 lated STAT1, the experimental measurements of cytoplasmic phosphorylated STAT1 107 (pST AT 1 † ) were mapped with the different cytoplasmatic complexes of pSTAT1 that 108 are described in the model as follows: where φ 2 is an scaling factor.

111
The experimental measurements of nuclear IRF9 (IRF 9 † ) were mapped with the nuclear 112 complexes involving IRF9 in the model as follows: where φ 3 is an scaling factor.   The initial particle numbers of MXA and IFIT1 of the model given in Table 3 of the 138 main section were calculated as follows: where φ 4 and φ 5 are scaling factors. Note that for the stochastic simulations the re-141 porter's initial conditions were sampled from a log-normal distribution which adequately 142 reproduces experimental data (Fig F)    The process to fit the stochastic system to flow cytometry data was developed based on Time series data of all flow cytometry measurements of unstimulated cells (0 UI/mL IFN) were combined and its mean value and standard deviation calculated to determine the formula of the logarithmic normal distribution for basal expression. The lognormal distribution for MxA is characterized by a log mean value of 8.85 and a a standard deviation on the log scale of 0.45, while IFIT1's initial particle number is sampled from a lognormal distribution with a log mean of 9.48 and a log sd of 0.45. The IFIT initial particle number was adjusted with 65 molecules times scaling factor 5 to adress for IFIT1 expression through cross-talk (see Section S3.1.2). The Kolmogorov-Smirnov distance D K S between experiment and log-normal distribution is below 0.05. In the plot, the x-axis represents the fluorescence level in arbitrary units of fluorescence. The integral over the density (area under the curve) is normalized so that it equals one. Prior to that, it was tested whether the conditions for assuming a normal distribution were given.  Genetic algorithms are stochastic search algorithms that resemble natural selection and sexual reproduction by mimicking the biological mechanisms of selection, recombination and mutation. This algorithm is made of a population of individuals (parameter sets), and each contains a genome that is defined by the number of parameters to optimize. The individuals are ranked after solving the objective function, and a population of parental individuals is selected according to an elitism rate (υ). New individuals (offspring) are generated by pairing and recombining the parental genomes (cross-over). Variability is introduced in the population by adding mutations in the new individuals according to a given mutation rate (µ). By the continuous process of selecting the best parameters after each generation, the algorithm evolves towards a minimum in parameter space. Our optimization strategy is based on Aguilera et al. [12]. The proposed method improves its performance by selecting parameters values after comparing the similitude between the first statistical moment of the system and the first statistical moment in the experimental data distribution. By this pre-selection of parameter values most of the original parameters are rejected and the algorithms focus on the finding of parameters that reproduce the observed distribution dynamics. This pre-step significantly reduces the computational cost.    MxA promoter contains two functional IRES sites (marked in red in the following se-220 quence) that have been experimentally proved to be activated by ISGF3. [16,17].  The range of the distributions is indicated by the light gray ribbons, while dark gray ribbons represent 50% KI intervals. B) Variability in the JAK-STAT signaling pathway was measured during different time points and for all the elements that form the pathway. The effects of extrinsic noise in the system were calculated by the coefficient of variation (cv = σ S /(µ S + 0.1) , where the subindex s represents the species in the pathway). In the plot the colorbar varies between 0 (white color) and larger than 4 (blue color), dark colors represent high variability in the dynamics of the studied species. The plots are consistent with those where a N:C ratio of 13.5% was assumed (Fig 5B). C) Stochastic simulations of different time points after IFN stimulation displaying a correct agreement in shape and location for multiple nuclear-cytoplasmic ratios. Simulated time-dependent distributions were computed by solving our model under stochastic dynamics using a distribution of values as initial conditions and repeating the simulations 1,000 times assuming a N:C ratio of 13.5% (black histogram) and 27% (red histogram), respectively. The red histogram is scaled by 24% to compensate for shifted steady state conditions.