https://orcid.org/0000-0001-8442-0536
Figures
Abstract
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.
Author summary
We investigate the impact of intrinsic and extrinsic noise on the reliability of interferon signaling. Information must be transduced robustly despite existing biochemical variability and at the same time the system has to allow for cellular variability to tune it against changing environments. Getting insights into stochasticity in signaling networks is crucial to understand cellular dynamics and decision-making processes. To this end, we developed a detailed stochastic computational model based on single cell data. We are able to show that reliability is achieved despite high noise at the receptor level.
Citation: Maier BD, Aguilera LU, Sahle S, Mutz P, Kalra P, Dächert C, et al. (2022) Stochastic dynamics of Type-I interferon responses. PLoS Comput Biol 18(10): e1010623. https://doi.org/10.1371/journal.pcbi.1010623
Editor: Christopher Rao, University of Illinois at Urbana-Champaign, UNITED STATES
Received: April 13, 2022; Accepted: September 30, 2022; Published: October 21, 2022
Copyright: © 2022 Maier et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All data and code is available via https://heibox.uni-heidelberg.de/d/3e7c5832e71e4f39a76f/.
Funding: The author(s) received no specific funding for this work.
Competing interests: The authors have declared that no competing interests exist.
Introduction
The interferon (IFN) system is the first line of innate immune defense. IFNs are polypeptides secreted by infected cells, inducing cell-intrinsic antimicrobial states that limit the spread of infectious agents. IFNs particularly act against viral pathogens in infected and neighboring cells. There are three distinct IFN families. The type-I IFN family comprises IFN-β and various subtypes of IFN-α. Most cell types produce IFN-β upon viral infection, whereas plasmacytoid dendritic cells are the predominant producers of IFN-α [1]. IFNγ is the only type-II IFN, and the type-III IFN family, the latest to be discovered, comprises IFNλ1–4. This study focuses on the widely studied type-I IFNs.
The IFN system is tightly regulated, as aberrant or excessive IFN responses can be detrimental and contribute to the pathogenesis of autoimmune diseases [2]. To prevent over-activation, multiple mechanisms determine the efficiency and flexibility in IFN-mediated responses, by amplifying or suppressing IFN-dependent signaling responses. Positive regulation of signaling includes the triggering of the JAK-STAT signaling pathway by binding of extracellular type-I IFN to the IFN-α receptor (IFNAR). Upon ligation, the two receptor chains, IFNAR1 and IFNAR2, heterodimerize and activate the receptor-associated protein tyrosine kinases Janus kinase 1 (JAK1) and tyrosine kinase 2 (TYK2), which in turn phosphorylate the signal transducers and activators of transcription 1 (STAT1) and STAT2. Phosphorylated STAT1 and STAT2 heterodimerize and together with IFN-regulatory factor 9 (IRF9) assemble the ternary complex called IFN-stimulated gene factor 3 (ISGF3) that translocates into the nucleus. ISGF3 acts as a transcription factor that directly activates the transcription of a set of several hundred IFN-stimulated genes (ISGs) and thus induces a cellular antiviral state [3]. IRF9 also constitutes a crucial positive feedback of JAK-STAT signaling, as its initial concentration is limiting and subsequently increased by ISGF3 activity [4]. Reciprocally, negative feedback regulation of the IFN response has also been well established and includes the induction of negative regulators in the JAK-STAT pathway, such as suppressor of cytokine signaling 1 (SOCS1) [3]. In addition to the receptor-dependent activation of ISGs, there is also an IFN-independent basal expression caused by constitutive formation of STAT2–IRF9 that translocates to the nucleus and initiates the transcription of many ISGs [5]. After IFN treatment, STAT2-IRF9 complexes are largely replaced by tyrosine-phosphorylated STAT1–STAT2 heterodimers, which have a lower dissociation rate and a higher quantity [5]. A graphical description of the IFN-induced response is given in Fig 1.
Free IFN binds to the IFNAR subunits 1 and 2 to form an active complex. After IFNAR activation the signal is transduced inside the cytoplasm; here STAT1 and STAT2 are phosphorylated. Phosphorylated forms of STAT1 and STAT2 form a heterodimer (dimerSTAT). dimerSTAT interacts with IRF9 to form ISGF3 (a trimolecular complex). ISGF3 translocates to the nucleus. In the nucleus, ISGF3 binds to free transcription factor binding sites (ISRE), inducing transcriptional activity leading to production of more IRF9 as a positive feedback and SOCS as a negative feedback given by the SOCS1 degradation of active receptors. ISGF3 induces the expression of around 350 different ISGs, including MxA and IFIT1. Additionally, there is a constitutive formation of STAT2-IRF9 heterodimers, that stimulate the expression of interferon-induced genes (ISGs) without a signaling requirement (basal expression). The model consists of 42 species in two compartments and 62 kinetic reactions and is fully described in the methods section. In the figure, boxes represent the chemical species, empty symbols represent degradation processes, arrows represent the reactions described in Section B.1 in S1 Text, initial conditions and parameters are given in Tables 1, 2 and 3, respectively. The pathway diagram was created with the Newt Editor [26, 27] following the Systems Biology Graphical Notation (SBGN) [28].
The expression of MxA and IFIT1 (previously referred to as ISG56) are gold standards for surrogate markers representing IFN-α’s biological activity in various experimental and clinical settings [6, 7]. MxA is one of the most commonly studied ISGs. Its expression is highly induced by IFN and it exerts strong antiviral activity against Orthomyxoviruses, most prominently Influenza A virus, but also other families of negative-strand RNA viruses such as Rhabdoviridae or Paramyxoviridae [8–10]. IFIT1 on the other hand was one of the first ISGs to be discovered, is the most prominent member of the viral stress-inducible gene family and is strongly induced in response to type-I IFN as well as various viruses [11–13].
As in many signalling pathways with low participating species concentrations, stochasticity certainly plays a role in the information processing. Even though our knowledge about the role of stochasticity in biological communication processes increased, the interpretation of stochastic signal transducing processes remains a challenge. The expression of many ISGs is believed to follow a digitalized (bi-modal all-or-none) expression pattern [14–18], however several studies are indicating that there are also ISGs that follow a graded (unimodal) pattern [19, 20]. Moreover, the number and percentage of ISGs displaying a bimodal expression varies between species and cell types [21–25], which is believed to reflect functionally important differences [15]. Despite evidence of all-or-nothing responses for MxA and IFIT1 gene induction [14, 15], the stochastic behavior of these genes cannot be fully explained yet. For instance, IFIT1 mRNA was shown to be bimodal in early time points before transitioning to a uniform response in Shalek et. al (2013) [15]. In order to answer this question and to better understand pathway information processing, the behaviour of these IFN activity markers has to be further verified. In the case of the IFN stimulated genes, neither the mechanisms nor the function of intrinsic variability in discrete reactions along the IFN activated pathway are well understood.
Previously, the dynamics of the JAK-STAT signaling pathway have been deterministically modeled and validated with time course data describing average dynamics in cell populations [4, 29–31]. Many questions could be answered with these studies. However, looking at the experimental data from single-cell measurements, it is obvious that stochastic influences play a major role. These data show distributions of behaviour rather than deterministic phenotypes. Obviously, the above efforts lack an exhaustive and systematic analysis of the influence of such stochastic fluctuations in the molecular reactions taking place along the pathway. The term “biochemical noise” is used to describe the inherently discrete and random interactions of molecules in the cellular environment [32]. Total biochemical noise is conventionally divided into two components: intrinsic and extrinsic noise. Intrinsic noise results from the probabilistic nature of molecular processes. Extrinsic noise results from (cell-to-cell) variations in the concentrations of cellular biomolecules such as ribosomes, enzymes, metabolites, overall proteins and nucleic acids [33]. Extrinsic and intrinsic noise have been proven to critically affect cellular decision-making processes [34]. Plenty of research has revealed genetic circuits that are subject to noise at the level of their components. In contrast, little is known about the effect of noise in cellular signaling pathways. However, previously, it has been shown that the information transduced in cellular signaling pathways is affected by noise [35]. Noise in cellular signaling pathways is conceptually different from the noise resulting from stochastic fluctuations of low-abundance genetic regulators. In cellular signaling pathways, noise can consist of low levels of enzyme activity through random transient interactions of receptors, binding of ligands to receptors in nonfunctional complexes, and background interactions along the pathway [36]. Noise plays an important role in IFN responses. It has been observed that within populations of nearly identical cells, not all cells that are infected by a virus or stimulated by IFN express ISGs [14, 37–41]. In all these cases the suggested mechanisms behind this heterogeneity has been related to stochastic events along the JAK-STAT signaling pathway.
Live-cell imaging data is indicating a multi-layered stochasticity (viral replication, IFN induction and IFN response), whereby heterogeneous IFN induction is largely determined by the translocation time of the key transcription factors and cell-intrinsic noise [14]. This result is challenged however by a recent integrated ChIP-seq and transcriptome analysis, which found that the molecular switch from an unstimulated homeostasis to a strong IFN-induced receptor-dependent response is caused by a rapid exchange of unphosphorylated transcription factors for the phosphorylated key transcription factors [5]. Furthermore, their results indicate that unstimulated and activated state are not only differing by intensity as previously assumed but also mechanistically as different signaling cascades are used, which has not been taken into account in earlier mathematical models and might be crucial to understand the above described (cell-intrinsic) noise. Hence, this study attempts to achieve greater detail than previous stochastic models, using more readily available flow cytometry data and novel literature findings.
Little is known about the sources and principles behind basal activity and expression in the JAK-STAT signaling pathway. Studies suggest that constitutive expression varies greatly between different cell types and occurs at different levels (receptor, signal transduction, gene expression) [42]. Basal STAT phosphorylation in diseased cells and tissues has been demonstrated by several groups [43–46]. Recently, it was shown that the basal expression of many interferon-induced genes (ISGs) is stimulated by STAT2–IRF9 complexes which form constitutively in the cytoplasm, shuttle between nucleus and cytoplasm and bind to the promoter region without a signaling requirement [5]. To which extent basal expression fluctuates within a population of genetically identical cells is, however, unknown so far.
Here, we aim to determine whether and how biochemical noise affects the information transduced in the JAK-STAT signaling pathway and investigate the transition between basal and activated state. To this end, we studied the stochastic responses of MxA and IFIT1 expression in Huh7.5 cells stimulated with IFN-α. Using fluorescent reporters under the control of the authentic promoter/enhancer region of IFIT1 and MxA we collected data displaying the differences between expressing and non-expressing cells for the marker genes in a time-course experiment. We hypothesize that the JAK-STAT signaling pathway efficiently transmits information under stochastic environments. To test our working hypothesis, we developed a detailed mathematical model using the obtained time-resolved flow cytometry data to describe the elements in the JAK-STAT signaling pathway at single-cell resolution. This model allowed us to systematically test the influence of intrinsic and extrinsic noise in the IFN response. We determined that the effects of intrinsic noise are particularly strong at the receptor level, during the formation of the transcription factor ISGF3 and during the transcription of the ISGs. We concluded that the JAK-STAT signaling pathway is a robust system that can filter extrinsic fluctuations.
Materials and methods
Experimental set up
The data described the expression of MxA and IFIT1 after IFN-α stimulation in a population of Huh7.5 cells. The experiments were done in the following way: Cells harboring a BAC (Bacterial Artificial Chromosome) containing the studied ISGs (mxa and ifit1) fused with the destabilized enhanced green fluorescent protein (deGFP) were seeded into a 6-well plate with 1.5x105 cells per well. Next day, cells were treated in a time-dependent manner with IFN-α (8, 12, 16, 20 and 32 h) and harvested at the same time point. Cells were fixed with 2% paraformaldehyde to stop further protein expression for 10 min, washed three times with PBS and applied for flow cytometry (BD Accuri C6). Measuring GFP fluorescence intensity served as marker for ISG induction. Data was evaluated with the FlowJo v10 software package. Human IFN-α (Alpha 2a) was obtained from Pbl Assay Science with a reported specific activity of 3.91x108 units/mg.
Mathematical model of the JAK-STAT signaling pathway upon IFN-α stimulation
The developed model consists of 42 species and 62 reactions (reactions m1 to m62). To name the variables in the model we used the following conventions: 1) variables referring to mRNA are denoted by m prefix. 2) Variables in phosphorylated use p as prefix. 3) Gene promoters are represented by the gene’s name in lowercase. 4) R1, R2, IR, AR and RC, represent the IFN receptor subunits, inactive, active and complex forms, respectively. 5) The compartment is superscripted to the species if the species exist in multiple compartments. A graphical representation of the interaction between variables in the model is given in Fig 1 and all reactions are listed in the Section B.1 in S1 Text. The model can be accessed via Biomodels: https://www.ebi.ac.uk/biomodels-main/ under reference: MODEL2210070001.
Mapping the system dynamics to experimental deterministic data
The measurements of time course data represent the average dynamics of proteins in a population of Huh7.5 cells [4]. Experimental data measurements were made comparable with the corresponding observable chemical species SO in the model by a function h that maps state SO(ti) to observation S†(ti) as follows:
(1)
at each measurement time point ti, i = 1, …, n; where n is the total number of observable time points.
Subsequently, average data measurements were made comparable with the corresponding observable chemical species in the model, considering a specific set of parameter values θ = {θ1, …, θd}, by the use of the following objective function:
(2)
Eq (2) compares the time courses for and
, and the difference between simulated and experimental data. A full description of the functions and scaling parameters used to map systems dynamics to experimental data is given in the Section C.1 in S1 Text.
Fitting the stochastic system to flow cytometry data
The process to fit the stochastic system to flow cytometry data was developed based on [47] and [48]. In short, empirical cumulative density functions (ECDF) were built for experimental data and simulations results. First, having nm repetitions of single-cell experimental data from flow cytometry measurements at I time points ti, i = 1, …, I, that is m(ti) = {m1(ti), …, mnm(ti)} ECDFs for the experimental data were built. In a similar way, considering a specific set of parameter values θ = {θ1, …, θd}, we performed ns repetitions of the stochastic simulations s(ti) = {s1(ti), …, sns(ti)}. The total of those stochastic simulations were used to build the ECDF for each ti that is
. To calculate the distance between
and
we used the Kolmogorov distance (DKS), that is the absolute value of the difference between the two distributions. For
and
their Kolmogorov distance is:
(3)
Parameter estimation strategy
Parameter searches consisted on optimization routines based on genetic algorithms implemented in R [49–51]. 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 [48]. The full strategy for parameter estimation is given in the Section C in S1 Text.
Numerical methods
Raw experimental data was analyzed using the function FCS data reader coded in Matlab [52]. Statistical moments were computed using R 4.1. The parameter estimation was performed using the genetic algorithm from the package “GA” coded in R [50, 51] and the COPASI R Connector [49]. Intrinsic noise was introduced by solving the model under stochastic dynamics using the adaptive explicit-implicit tau-leaping method with automatic tau selection [53] coded in COPASI 4.30 [54]. Extrinsic noise was simulated following the method suggested by Shahrezaei et al., [55] that proposes a modification of the standard SSA introducing extrinsic variability in the model parameters by random sampling a distribution. For our model, we focus our attention on understanding the effect of extrinsic variations in the initial conditions. For this, we randomly sampled a Normal distribution with values for μ given in Table 1 and one of three values of σ: 0, 0.3, 0.6. In a final step, a parameter-robustness analysis was carried out to determine to what extent the evaluated system functionality is preserved under considered perturbations. Parameters were altered between half and double their original values individually and after subsequent stochastic simulations, the effect on the ISG induction was quantified using the KS-distance.
Results
Heterogeneity in IFN responses
To analyze the dynamics of the IFN signaling pathway, we employed Huh7.5 cells stably transfected with bacterial artificial chromosomes (BACs) containing MxA or IFIT1, respectively, tagged with a destabilized enhanced green fluorescent protein (deGFP) in their genuine genomic context, including up- and downstream regulatory sequences as we described previously [38, 40]. This reporter system allowed for fluorescence activated flow cytometry (FACS) based quantification of the authentic induction of the respective ISG expression upon IFN stimulation. We implemented a work-flow to measure fluorescence from the reporter genes at a single-cell level by FACS and their distribution in the cell population, thereby capturing cell-to-cell variability of the IFN response (Fig 2A). Our results indicated that, both, MxA and IFIT1 sustained graded responses, and heterogeneity in the IFN response was observed by the deformation and shifting of an unimodal distribution. In this unimodal distribution, responder cells are in the upper percentile and the center of the curve, and the non-responder cells are in the lowest percentile. We observed that IFN stimulation shifted the whole distribution to higher fluorescence intensities, i.e. ISG concentrations (Fig 2B and 2C). This contrasts somewhat with previous publications indicating an all-or-none switch-like response with two discernible subpopulations of responder and non-responder cells, based on IRF7 that was used as the reporting ISG [14]. However, apart from the ISG, also the cell system and species (mouse) was different than in our study, suggesting a possible cell-type or species specific response behavior.
A: Experimental set-up showing the graded expression of MxA and IFIT1 after IFN-α stimulation in a population of Huh7.5 cells. A threshold is defined to differentiate responder vs. non-responder cells in the population. B: Distributions represent the flow cytometry measurements of MxA expression at different time points after IFN stimulation (blue distribution). C: Distributions represent the flow cytometry measurements of IFIT1 expression at different time points after IFN stimulation (red distribution). In figures B and C, the vertical dashed lines represents a threshold value (defined as the mean plus two standard deviations) calculated form the unstimulated populations. This threshold value is used to differentiate responser vs non-responder cells in a cell population with a graded response. In the observed treatments, the mean fluorescence level shifts from 8x103 a.u. (arbitrary units of fluorescence) to 4.5x104 a.u. for MxA, and from 2x104 to 9x104 a.u. for IFIT1.
We further tested in what way different concentrations of IFN affect the obtained graded response by comparing different doses of IFN (0, 10, 50, 250, 1250 UI/mL) and different time points (0, 8, 12, 16, 20 and 32 hours). We observed a dose-dependent shift of the population to higher signal intensities starting from 16 h post stimulation for MxA (Fig A in S1 Text) and from 12 h for IFIT1 (Fig B in S1 Text), but the graded responses were always maintained at different doses of IFN stimulation and for different time points. This may exclude the possibility of temporal monostable dynamics and then an all-or-none switch response. By plotting the mean fluorescence intensities (MFI) we analyzed time dependence of MxA and IFIT1 responses (Fig C in S1 Text). Our results indicated that IFIT1 reached its maximum response after 16 hours of IFN stimulation, whereas the maximum response for MxA was obtained at 32 hours of treatment. In summary, our results confirmed that MxA and IFIT1 only support graded responses, that those responses- as expected- depended on the dose of IFN used to stimulate the cells, and that MxA and IFIT1 exhibited distinct temporal induction dynamics.
Stochastic modeling of the JAK-STAT signaling pathway and parameter estimation
In order to develop a deeper understanding of the JAK-STAT signaling pathway we devised a new stochastic model. The developed mathematical model has single-cell resolution comprising two compartments, known key components, feedback responses, constitutive regulatory mechanisms and a basal expression. The model architecture is based on previous models that consider the well accepted topology of the pathway [4, 29–31]. Contrary to these previous models and a more recently published model [65], we constructed a stochastic model and implemented a detailed model of the basal expression. Additionally, the receptor-independent activation of interferon stimulated genes via STAT2-IRF9 heterodimers was included as recent studies identified them as key players for basal expression and rapid activation of the signaling pathway [5]. The developed model consists of 42 species and 62 kinetic reactions, the complete model is given in Fig 1. All enzymatic reactions are assumed to be mass action kinetics (Table A in S1 Text). Processes of gene induction were modeled following a two-state promoter scheme [66]. To solve the model, we used the Adaptive Stochastic Simulations Algorithm (SSA) with τ-leap approximation [53, 67], which was shown to lead to the same results as the standard SSA algorithm. This was also confirmed by us, as can be seen in Fig D in S1 Text. Extrinsic variability in protein concentration was introduced in the initial conditions by random sampling from a normal distribution following the scheme introduced by Shahrezaei et al. [55].
To identify the model parameters, we made use of our recently developed algorithm that fits single-cell data to stochastic models [48]. The method uses a comparison based on Kolmogorov-Smirnov (KS)-distances of the cumulative density functions that are computed based on Monte Carlo simulations and the experimental data. Multiple parameter values are iteratively evaluated using genetic algorithms (Fig 3A). The complete parameter estimation method is described in the Section C in S1 Text. To increase certainty in our analysis we performed a bibliographic search to narrow down parameter space with already determined parameter values (Table 2). The parametrized model was capable to reproduce the observed IFN response for the full probability distributions for MxA and IFIT1 at five different time points and different doses (Fig 3B and Fig H in S1 Text). Please note that at higher doses the variance is slightly underestimated which leads to an overestimation of the amplitude of the distribution, since total numbers have to be the same. The fit could possibly be further improved by including important elements whose implications in the JAK-STAT signaling pathway have only recently been described and which were not considered in this model (for example basal expression of phosphorylated forms of STAT proteins [64], receptor recycling and endocytic regulation [68], and more elements involved in the negative feedback [3]). Including those additional elements in our model is possible in principle. Nevertheless, this would increase the number of free parameters in the system, probably causing overfitting problems as described in the Section D in S1 Text. The overall good agreement indicated that the model can capture the heterogeneity in the IFN response. We also tested our model for resolving the dynamics of the pathway at different levels of the signaling cascade. For that purpose, we exploited previously published population-level data for the activation of the key signaling components JAK1, STAT1 and IRF9, all measured in Huh7.5 cells stimulated with IFN-α [4]. Notably, the average results of our stochastic model largely matched the temporal response of those signaling components (Fig 4A). Hence we conclude that our model (using a single parameter set) based on two different experimental data sources (flow cytometry data and immunoblotting measurements) can faithfully reproduce the system’s dynamics.
A: The parameter estimation strategy consist of optimization routines based on genetic algorithms. The proposed methods measure the similitude between experimental and simulated distributions using KS-distance. The best parameters are obtained by minimizing the KS-distance. The full strategy for parameter estimation is given in the Section C in S1 Text. B: Experimental time-dependent distributions were computed from the flow cytometry datasets (filled histograms). Simulated time-dependent distributions were computed by solving our model under stochastic dynamics and repeating the simulations 1,000 times (red density plots). In the plots, the y-axis represents the normalized cell count and the x-axis represents the fluorescence quantity (arbitrary units, a.u.) associated with the expression of the MxA and IFIT1 proteins at various time points after stimulating Huh7.5 cells with 250 UI/mL IFN (5,000 IFN molecules). For each distribution, the median (M) and variance (s) is given. The initial conditions are given in Table 1, fitted parameter values are given in Table 2 and compartment sizes are given in Table 3. See Figs H and I in S1 Text for fits using different IFN doses.
A: Cell population data describing the temporal dynamics of phosphorylated JAK1, pSTAT1 and nuclear IRF9. Experimental data describe quantitative immunoblotting measurements in Huh7.5 cells after stimulation with 500 UI/mL of IFN-α at different time points for a total time of 180 min [4]. Measurement error of 18% is represented in the figures as error bars. Model dynamics were obtained by repeating the stochastic simulation 1,000 times, each trajectory representing a single cell. The procedure for mapping the model variables and experimental data is given in Section C.1 in S1 Text. B: Time courses data describing the temporal dynamics of all species involved in the JAK-STAT signaling pathway. The plots represent the repetitions of the stochastic model. The y-axis has units of Molecules per Cell (M/C). Orange lines represent the median, while light gray display the range of values and dark gray ribbons contain 50% of the values. The ribbons may not be visible if there is almost no variation (e.g. IFN) or if there is no particle present for the majority of trajectories (e.g. mIRF9n).
Effect of noise in the ISG induction after IFN stimulation
Intrinsic noise is caused by the stochasticity of the involved processes meaning that there is variability even if the same initial conditions are used whereas extrinsic noise originates in the cell-to-cell variability of protein abundancies etc. In order to compare and analyse these two effects we computed repeated time-series for both scenarios.
Intrinsic noise in the JAK-STAT signaling pathway.
Here, we simply repeated time course simulations many times. Simulations were run for 32 hours and 1,000 repetitions were used to obtain a landscape of the heterogeneity in the population of trajectories. These different trajectories represent the effects of the intrinsic noise and are solely caused by the stochasticity of the individual processes. Our results suggest differences in the temporal responses: species involved in the ligand-receptor activation steps show fast but transient dynamics, whereas elements involved in positive (IRF9) and negative feedback (SOCS) loops display slower and more sustained dynamics (Fig 4B and Fig J in S1 Text). Fig 5B, first panel, depicts (for σ equal zero) the intrinsic variation of species over time in a different way. The figure indicates that the highest intrinsic variation lies in the short-lived RC. Otherwise, there is an expected trend that species with low particle numbers show higher variability (for instance SOCS) than species with higher particle numbers (for example STAT1c). However, there are some notable exceptions. Thus, mIRF9 and IRF9 itself which plays an important role in the system show low variability despite the fact that there are only small numbers of the respective molecules around. Of course, such exceptions can be explained by the fact that the stochasticity is not only influenced by particle numbers, but also by reaction rates and their ratios. Therefore, it is noteworthy that the used visualization offers a quick and not completely forseeable overview of the variability of the components in this system.
A: Extrinsic noise in the system was introduced by considering the effect of variability in the initial copy number of the proteins in the pathway. Initial conditions were generated by random sampling using a normal distribution with values for μ given in Table 1 and one of three values of σ: 0, 0.3, 0.6. B: Variability in the elements of the JAK-STAT signaling pathway over time. 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. Species were ordered based on the average coefficient of variation over all time points for systems with an extrinsic noise of σ = 0.3. The overall system dynamics under the influence of extrinsic noise can also be observed in Figs J to L in S1 Text. C: Stochastic simulations of different time points after IFN stimulation using a distribution of values as initial conditions. The results show a transient perturbation in the MxA and IFIT1 expression when extrinsic fluctuations are considered. As reference, experimental data for a system without extrinsic noise (σ = 0.0) are given (gray histograms). D: Distributions at multiple time points considering different strengths of extrinsic noise vs a system with only intrinsic noise are compared using the Kolmogorov-Smirnov distance.
Extrinsic noise in the JAK-STAT signaling pathway.
Extrinsic fluctuations can originate from cells undergoing different stages of their cell cycle, and cell-to-cell variability in the copy number of proteins inherited from parent cells during cell division [34]. We focused our attention on the identification of the effect of variability in the copy number of the proteins in the pathway. To this end, we simulated the effect of having a distribution of values as initial conditions instead of fixed values. Initial conditions were generated by random sampling using a normal distribution with values for μ given in Table 1 and one of three values of σ: 0 (again, only intrinsic noise), 0.3, 0.6. By introducing wider distributions as initial condition we tested noisier extrinsic effects in our system (Fig 5A). Subsequently, we quantified the size of the variability in the JAK-STAT signaling pathway by computing the coefficient of variations, which is defined as the ratio of the standard deviation, σs, to the mean μs, where the subindex s represents the species in the pathway. It is important to remark, that the coefficient of variation is very sensitive to small changes in the mean if the mean value is near zero. Given that there are three species present in the model that are expressed with a copy number of less than one on average (RC, mMxA and mIFIT1), uncertainty in low-probability events plays a major role. Hence, the formula for the coefficient of variations was modified as follows to compensate for high noise caused by low discrete particle numbers:
(4)
Note, that the effect of the modification is negligible for species with a mean value above one, while statistic uncertainty is reduced for species with a very low average expression.
In contrast to the findings for variability due to intrinsic noise as described above, effects of extrinsic noise are not governed as strongly by the abundance of species. Thus, species that are close to the beginning of the signalling cascade (Fig 5B right column) often show a high level of variability even if they are abundant like pSTAT1 and pSTAT2. The species that constitute the response of the pathway, especially MxA show a lower variability. This is also preseved when increasing the level of extrinsic noise. However, differences in local time-dependent variance are not huge and may not constitute the best way to look at the robustness of the response.
Therefore, we quantified the similarity between responses by measuring the KS-distance between distributions for MxA and IFIT1 under multiple values of extrinsic variability at different time points. Even though, there is clearly an effect of extrinsic noise, the response is conserved with a KS-distance of less than 0.1 for σ equal 0.3 and around 0.2 for σ equal 0.6 as can be observed in Fig 5C and 5D. In summary, our results indicate that despite variability of individual trajectories, there is an efficient robust transmission of the signal even in noisy environments.
It is noteworthy that fitting of the distribution of the experimental data (Fig 5C) seems to be better if moderate levels of extrinsic noise are added. This observation is in line with previous studies which found that models that incorporate comparable levels of extrinsic noise best resembles their experimental data [83].
Parameter robustness in the JAK-STAT signaling pathway.
To calculate the quality of the obtained parameters it is common practice to repeat the optimization method hundreds or thousands of times until a posteriori distribution for each estimated parameter are obtained [84]. This strategy is commonly used for deterministic systems or stochastic systems with a small number of variables. Calculating posteriori distributions for the estimated parameter values is used to drive conclusions referent to the sensitivity or robustness of each parameter. In our specific model, given the large number of variables, this strategy is computationally prohibitive. For this reason, we performed a parameter robustness assay to determine to what extent the evaluated system functionality is preserved under considered perturbations (Fig 6). Parameters were altered between half and double their original values individually and after subsequent stochastic simulations, the effect on the ISG induction was quantified using the KS-distance. Given the low KS-distance for all parameter alterations involved in receptor dynamics and signal transduction, it can be assumed that the JAK-STAT signaling pathway is robust to parameter variations. However, higher KS-values have been observed for parameter alterations of parameters involved in the gene expression and degradation of ISGs (transcription, translation, mRNA and protein degradation). Based on the same data, we investigated which parameters prevent reliable transduction of the signal from receptor to gene expression (Fig 7). Again, it was found that perturbing processes in the gene induction (transcription & translation) and decay of the ISGs MxA and IFIT1 altered the average abundance of ISGs while perturbations of other parts of the signal pathway did not cause considerable changes.
Parameters were altered between half (0.5 * k) and double (2.0 * k) their original values individually to determine to what extent the evaluated system functionality is preserved under considered perturbations. The effect on the ISG induction was quantified by repeating the stochastic simulation 600 times and computing the KS-distance to the unpertubated system.
Parameters were altered between half (0.5 * k) and double (2.0 * k) their original values individually to determine which parameters prevent reliable transduction of the signal from receptor to gene expression under considered perturbations. The effect on the ISG induction was quantified by repeating the stochastic simulation 600 times and computing the relative change in IFIT and MxA molecules to the unpertubated system.
Discussion
During the last two decades, important efforts have been made to increase the understanding of the IFN system by integrating experimental data sets and mathematical models. Existing models of the JAK-STAT signaling pathway describe the deterministic dynamics in cell populations [4, 29–31]. All those efforts neglect the influence of molecular noise in the biochemical reactions taking place along the pathway. Noise has been proven to be a fundamental force affecting multiple cellular pathways, and IFN is no exception. It has been well documented that populations of identical cells produce heterogeneous responses to IFN stimulation [17, 38–41]. The lack of stochastic models to describe the IFN system can be attributed to the scarcity of single-cell data and efficient methods to analyze and make use of this data. Live cell imaging and flow cytometry are among the methods that can be used to measure protein levels in high numbers of single cells. Given that flow cytometry is a relatively easy technique obtaining single-cell measurements of thousands of cells, our approach could be widely translated to other signaling pathways and contexts. Additionally, new methods have recently been introduced to efficiently fit stochastic models to such single-cell data [47, 48]. To the best of our knowledge, we here present the first mechanistic model of the JAK-STAT signaling pathway that is validated with single-cell data and can explain both basal (Fig N in S1 Text) and stimulated states (Fig J in S1 Text).
In summary, we studied the stochasticity of MxA and IFIT1 induction triggered by stimulation of Huh7.5 cells with IFN-α for a timescale of 32 hours. To test our working hypothesis, we developed a detailed mathematical model and calibrated it using time-course, single-cell flow cytometry data. This model allowed us to systematically test the influence of intrinsic and extrinsic noise in the IFN response. The major novelty of our model is the fact that it is mechanistic, incorporating the well known architecture of the pathway including major feedback systems and basal states. This contrasts with previous stochastic models that use minimalistic approaches and neglect unstimulated homeostatic states [14]. Unlike previous studies, we present a model that fits two experimental data sets that represent full probability distributions at different time points and immunoblotting measurements. It is important to remark that our model reproduces dynamics within the correct order of magnitude reported for mammalian mRNAs (reported average of 17 Molecules/Cell with range between 1 to 200 Molecules/Cell) and proteins (reported average of 50,000 Molecules/Cell with range between 100 to 108 Molecules/Cell) [71]. Comparing our molecule numbers for STAT1, STAT2 and IRF9-complexes to a more recent publication [65], we can also state that the numbers for STAT1 and STAT2 coincide well while our numbers for IRF9 are lower than measured in this study. However, they correspond to earlier measurements [4] by the same group.
Our results support graded (unimodal) dynamics for MxA and IFIT1. This clearly contrasts with previous reports based on single-cell data for IRF7 (IRF7 is an ISG with a central role in the immune response) where an all-or-none switch response (two subpopulations, responder and non-responder) was observed [14, 15]. We argue that differences in the expression of distinct ISGs may be explained by each ISG promoter architecture. For example, the IRF7 promoter contains two different transcription factor binding sites (ISRE and IRF-E) that are activated by ISGF3 and an IRF7 dimer, respectively. Bimodality in IRF7 expression can be justified by a circuit with a positive feedback loop and the non-linearity caused by the complex activation of its promoter as we showed earlier [48]. On the contrary, our experiments for MxA and IFIT1 show graded responses that were consistently obtained for all IFN doses and multiple time points. MxA and IFIT1 promoters only contain two binding sites for ISGF3 (Section F in S1 Text), and cooperativity has not been proven to take place during type-I IFN responses [85]. The lack of cooperative behavior in the MxA and IFIT1 promoters can explain the observed graded response [85]. It is also tempting to speculate that an all-or-none switch response vs graded responses for the ISG induction might be explained by the specific gene regulatory sequences encoded in the gene promoter (Fig M in S1 Text). Moreover, given the simple architecture of their transcriptional sites (only ISREs), very similar MxA and IFIT1 dose-response expression have been observed in different individuals [86, 87]. This suggests that MxA and IFIT1 response to IFN-α is robust and explains the graded response that was consistently obtained for all IFN-α doses at different time points. Maiwald et al. [4] demonstrated that the dynamic of ISG response and thus the expression of ISGs is largely dependent on the initial concentration of the signal-transduction molecules, which differs both between experimental set-ups and neighboring cells. Hence, contradicting results may be caused by the choice of the cell line, the initial states of the signal-transduction molecules in cells. Finally, we cannot rule out the fact that MxA and IFIT1 graded responses are cell-type specific. Still, we tested a different cell line A549. Our results confirm that MxA and IFIT1 show a graded response upon IFN stimulation and is also consistent with previous reports for MxA [39].
Most initial conditions in the model are assumed to be non-zero before IFN stimulation. This assumption can be justified by previous studies showing a basal expression (around a 10% of the total protein concentration) of phosphorylated forms of STAT proteins [64] and the formation of unphosphorylated STAT2-IRF9 transcription factors [5] even in the absence of IFN stimulation. Additionally, the immunoblotting measurements for IRF9, STAT1 and the active form of the IFN receptor show basal levels of expression even before IFN stimulation [4]. These basal expression levels can be caused by the low and transient enzymatic activities between elements in the pathway [36], or by paracrine responses in the cell culture even before the application of the external IFN stimulus [14].
Conclusion
By challenging our model we could determine a low intrinsic noise for all species except the receptor complex (RC) and mIFIT1, where the variability is caused by uncertainty in low-probability events. In contrast, high extrinsic noise was observed at the receptor level and during the formation of the transcription factor ISGF3. Nonetheless, extrinsic noise seems to play only a minor role for the JAK-STAT signaling pathway as strong and stable ISG expression was observed for all simulation conditions despite extrinsic noise affecting the dynamics of cellular constituents locally. Additionally, kinetic parameters in the JAK-STAT pathway were found to be resistent to parameter alterations, which supports our hypothesis that the JAK-STAT pathway is a robust system filtering extrinsic noise even further. From a biological point of view this makes perfect sense as multi-cellular organisms depend on a robustly functioning first line of defense against invading pathogens.
Acknowledgments
The authors thank Dr. Frank T. Bergmann (Heidelberg Germany) for software support. The authors are also very thankful to Dr. Jesus Rodriguez Gonzalez (Cinvestav Mexico) for valuable feedback on the manuscript.
References
- 1. Horner SM, Gale M Jr. Regulation of hepatic innate immunity by hepatitis C virus. Nature medicine. 2013;19(7):879–888. pmid:23836238
- 2. Hertzog PJ, Williams BR. Fine tuning type I interferon responses. Cytokine & growth factor reviews. 2013;24(3):217–225. pmid:23711406
- 3. Ivashkiv LB, Donlin LT. Regulation of type I interferon responses. Nature reviews Immunology. 2014;14(1):36–49. pmid:24362405
- 4. Maiwald T, Schneider A, Busch H, Sahle S, Gretz N, Weiss TS, et al. Combining theoretical analysis and experimental data generation reveals IRF9 as a crucial factor for accelerating interferon α-induced early antiviral signalling. FEBS journal. 2010;277(22):4741–4754. pmid:20964804
- 5. Platanitis E, Demiroz D, Schneller A, Fischer K, Capelle C, Hartl M, et al. A molecular switch from STAT2-IRF9 to ISGF3 underlies interferon-induced gene transcription. Nature Communications. 2019. pmid:31266943
- 6. Haller O, Kochs G. Human MxA protein: an interferon-induced dynamin-like GTPase with broad antiviral activity. Journal of Interferon & Cytokine Research. 2011;31(1):79–87. pmid:21166595
- 7. Levy D, Larner A, Chaudhuri A, Babiss LE, Darnell JE Jr. Interferon-stimulated transcription: isolation of an inducible gene and identification of its regulatory region. Proc Natl Acad Sci U S A. 1986;83(23):8929–8933. pmid:3466167
- 8. Haller O, Gao S, Von Der Malsburg A, Daumke O, Kochs G. Dynamin-like MxA GTPase: Structural insights into oligomerization and implications for antiviral activity; 2010.
- 9. Roers A, Hochkeppel HK, Horisberger MA, Hovanessian A, Haller O. MxA gene expression after live virus vaccination: A sensitive marker for endogenous type i interferon. Journal of Infectious Diseases. 1994;169(4):807–813. pmid:7510764
- 10. Gilli F, Marnetto F, Caldano M, Sala A, Malucchi S, Capobianco M, et al. Biological markers of interferon-beta therapy: comparison among interferon-stimulated genes MxA, TRAIL and XAF-1. Mult Scler. 2006;12(1):47–57. pmid:16459719
- 11. Guo J, Peters KL, Sen GC. Induction of the human protein P56 by interferon, double-stranded RNA, or virus infection. Virology. 2000;267(2):209–19. pmid:10662616
- 12. Fensterl V, Sen GC. The ISG56/IFIT1 gene family. Journal of interferon & cytokine research: the official journal of the International Society for Interferon and Cytokine Research. 2011;31(1):71–78. pmid:20950130
- 13. Wathelet MG, Clauss IM, Content J, Huez GA. Regulation of two interferon-inducible human genes by interferon, poly (rI)· poly (rC) and viruses. The FEBS Journal. 1988;174(2):323–329. pmid:2454816
- 14. Rand U, Rinas M, Schwerk J, Nöhren G, Linnes M, Kröger A, et al. Multi-layered stochasticity and paracrine signal propagation shape the type-I interferon response. Molecular systems biology. 2012;8(1):584. pmid:22617958
- 15. Shalek AK, Satija R, Adiconis X, Gertner RS, Gaublomme JT, Raychowdhury R, et al. Single-cell transcriptomics reveals bimodality in expression and splicing in immune cells. Nature. 2013;498(7453):236–240. pmid:23685454
- 16. Bhushal S, Wolfsmüller M, Selvakumar TA, Kemper L, Wirth D, Hornef MW, et al. Cell polarization and epigenetic status shape the heterogeneous response to type III interferons in intestinal epithelial cells. Front Immunol. 2017;8:671. pmid:28659914
- 17. El-Sherbiny YM, Psarras A, Md Yusof MY, Hensor EMA, Tooze R, Doody G, et al. A novel two-score system for interferon status segregates autoimmune diseases and correlates with clinical features. Scientific Reports. 2018;8(1):5793. pmid:29643425
- 18. Czerkies M, Korwek Z, Prus W, Kochańczyk M, Jaruszewicz-Błońska J, Tudelska K, et al. Cell fate in antiviral response arises in the crosstalk of IRF, NF-κB and JAK/STAT pathways. Nature Communications. 2018;9(1):493. pmid:29402958
- 19. Mostafavi S, Yoshida H, Moodley D, LeBoité H, Rothamel K, Raj T, et al. Parsing the Interferon Transcriptional Network and Its Disease Associations. Cell. 2016;164(3):564–578. pmid:26824662
- 20. O’Neal JT, Upadhyay AA, Wolabaugh A, Patel NB, Bosinger SE, Suthar MS, et al. West Nile Virus-Inclusive Single-Cell RNA Sequencing Reveals Heterogeneity in the Type I Interferon Response within Single Cells. Journal of Virology. 2019;93(6):e01778–18. pmid:30626670
- 21. Bar-Even A, Paulsson J, Maheshri N, Carmi M, O’Shea E, Pilpel Y, et al. Noise in protein expression scales with natural protein abundance. Nature Genetics. 2006;38(6):636–643. pmid:16715097
- 22. Sigal A, Milo R, Cohen A, Geva-Zatorsky N, Klein Y, Liron Y, et al. Variability and memory of protein levels in human cells. Nature. 2006;444(7119):643–646. pmid:17122776
- 23. Tang F, Barbacioru C, Wang Y, Nordman E, Lee C, Xu N, et al. mRNA-Seq whole-transcriptome analysis of a single cell. Nature Methods. 2009;6(5):377–382. pmid:19349980
- 24. Islam S, Kjällquist U, Moliner A, Zajac P, Fan JB, Lönnerberg P, et al. Characterization of the single-cell transcriptional landscape by highly multiplex RNA-seq. Genome Res. 2011;21(7):1160–1167. pmid:21543516
- 25. Hashimshony T, Wagner F, Sher N, Yanai I. CEL-Seq: Single-Cell RNA-Seq by Multiplexed Linear Amplification. Cell Reports. 2012;2(3):666–673. pmid:22939981
- 26. Sari M, Bahceci I, Dogrusoz U, Sumer SO, Aksoy BA, Babur Ã, et al. SBGNViz: A Tool for Visualization and Complexity Management of SBGN Process Description Maps. PLOS ONE. 2015;10(6):1–14. pmid:26030594
- 27. Balci H, Siper MC, Saleh N, Safarli I, Roy L, Kilicarslan M, et al. Newt: a comprehensive web-based tool for viewing, constructing and analyzing biological maps. Bioinformatics. 2020;37(10):1475–1477.
- 28. Le Novère N, Hucka M, Mi H, Moodie S, Schreiber F, Sorokin A, et al. The Systems Biology Graphical Notation. Nature Biotechnology. 2009;27(8):735–741. pmid:19668183
- 29. Rybiński M, Gambin A. Model-based selection of the robust JAK-STAT activation mechanism. Journal of theoretical biology. 2012;309:34–46. pmid:22677400
- 30. Smieja J, Jamaluddin M, Brasier AR, Kimmel M. Model-based analysis of interferon-β induced signaling pathway. Bioinformatics. 2008;24(20):2363–2369. pmid:18713791
- 31. Yamada S, Shiono S, Joo A, Yoshimura A. Control mechanism of JAK/STAT signal transduction pathway. FEBS letters. 2003;534(1):190–196. pmid:12527385
- 32.
Andrews SS, Dinh T, Arkin AP. Stochastic models of biological processes. Springer; 2009.
- 33. Elowitz MB, Levine AJ, Siggia ED, Swain PS. Stochastic gene expression in a single cell. Science. 2002;297(5584):1183–1186. pmid:12183631
- 34. Joo J, Plimpton SJ, Faulon JL. Statistical ensemble analysis for simulating extrinsic noise-driven response in NF-κB signaling networks. BMC systems biology. 2013;7(1):1. pmid:23742268
- 35. Cheong R, Rhee A, Wang CJ, Nemenman I, Levchenko A. Information transduction capacity of noisy biochemical signaling networks. science. 2011;334(6054):354–358. pmid:21921160
- 36. Ladbury JE, Arold ST. Noise in cellular signaling pathways: causes and effects. Trends in biochemical sciences. 2012;37(5):173–178. pmid:22341496
- 37. Hu J, Sealfon SC, Hayot F, Jayaprakash C, Kumar M, Pendleton AC, et al. Chromosome-specific and noisy IFNB1 transcription in individual virus-infected human primary dendritic cells. Nucleic Acids Research. 2007;35(15):5232–5241. pmid:17675303
- 38. Bauhofer O, Ruggieri A, Schmid B, Schirmacher P, Bartenschlager R. Persistence of HCV in quiescent hepatic cells under conditions of an interferon-induced antiviral response. Gastroenterology. 2012;143(2):429–438. pmid:22522091
- 39. von Recum-Knepper J, Sadewasser A, Weinheimer VK, Wolff T. Fluorescence-Activated Cell Sorting-Based Analysis Reveals an Asymmetric Induction of Interferon-Stimulated Genes in Response to Seasonal Influenza A Virus. Journal of virology. 2015;89(14):6982–6993. pmid:25903337
- 40. Schmid B, Rinas M, Ruggieri A, Acosta EG, Bartenschlager M, Reuter A, et al. Live Cell Analysis and Mathematical Modeling Identify Determinants of Attenuation of Dengue Virus 2’-O-Methylation Mutant. PLoS Pathog. 2015;11(12):e1005345. pmid:26720415
- 41. Zawatzky R, De Maeyer E, De Maeyer-Guignard J. Identification of individual interferon-producing cells by in situ hybridization. Proceedings of the National Academy of Sciences. 1985;82(4):1136–1140. pmid:3856251
- 42. Zurney J, Howard KE, Sherry B. Basal expression levels of IFNAR and Jak-STAT components are determinants of cell-type-specific differences in cardiac antiviral responses. Journal of Virology. 2007;81(24). pmid:17942530
- 43. Lin TS, Mahajan S, Frank DA. STAT signaling in the pathogenesis and treatment of leukemias. Oncogene. 2000; p. 2496–2504. pmid:10851048
- 44. Bowman T, Garcia R, Turkson J, Jove R. STATs in oncogenesis. Oncogene. 2000;19(21):2474–2488. pmid:10851046
- 45. Garcia R, Bowman T, Niu G, Yu H, Minton S, Muro-Cacho C, et al. Constitutive activation of Stat3 by the Src and JAK tyrosine kinases participates in growth regulation of human breast carcinoma cells. Oncogene. 2001;20(20):2499–2513. pmid:11420660
- 46. Mora LB, Buettner R, Seigne J, Diaz J, Ahmad N, Garcia Rea. Constitutive activation of Stat3 in human prostate tumors and cell lines: direct inhibition of Stat3 signaling induces apoptosis of prostate cancer cells. Cancer research. 2002;62(22):6659–6666. pmid:12438264
- 47. Lillacci G, Khammash M. The signal within the noise: efficient inference of stochastic gene regulation models using fluorescence histograms and stochastic simulations. Bioinformatics. 2013;29(18):2311–2319. pmid:23821649
- 48. Aguilera LU, Zimmer C, Ursula K. A New Efficient Approach to Fit Stochastic Models on the Basis of High-throughput Experimental Data Using a Model of IRF7 Gene Expression as Case Study. BMC Systems Biology. 2017;11(26). pmid:28219373
- 49. Förster J, Bergmann FT, Pahle J. Bioinformatics. 2020
- 50. Scrucca L. GA: A Package for Genetic Algorithms in R. Journal of Statistical Software. 2013;53(4):1–37.
- 51. Scrucca L. On Some Extensions to GA Package: Hybrid Optimisation, Parallelisation and Islands Evolution. The R Journal. 2017; p. 187–206.
- 52.
Balkay L. FCS data reader. MATLAB Central File Exchange. 2006.
- 53. Cao Y, Gillespie DT, Petzold LR. Adaptive explicit-implicit tau-leaping method with automatic tau selection. Journal of Chemical Physics. 2007;126(22). pmid:17581038
- 54. Hoops S, Sahle S, Gauges R, Lee C, Pahle J, Simus N, et al. COPASI–a COmplex PAthway SImulator. Bioinformatics (Oxford, England). 2006;22(24):3067–74. pmid:17032683
- 55. Shahrezaei V, Ollivier JF, Swain PS. Colored extrinsic fluctuations and stochastic gene expression. Molecular systems biology. 2008;4(1):196. pmid:18463620
- 56. Arnheiter H, Ohno M, Smith M, Gutte B, Zoon KC. Orientation of a human leukocyte interferon molecule on its cell surface receptor: carboxyl terminus remains accessible to a monoclonal antibody made against a synthetic interferon fragment. Proceedings of the National Academy of Sciences. 1983;80(9):2539–2543. pmid:6302694
- 57. Theofilopoulos AN, Baccala R, Beutler B, Kono DH. Type I interferons (α/β) in immunity and autoimmunity. Annu Rev Immunol. 2005;23:307–335. pmid:15771573
- 58.
DeMaeyer EM, Maeyer-Guignard D, et al. Interferons and other regulatory cytokines. Wiley; 1988.
- 59. Wilmes S, Beutel O, Li Z, Francois-Newton V, Richter CP, Janning D, et al. Receptor dimerization dynamics as a regulatory valve for plasticity of type I interferon signaling. The Journal of cell biology. 2015;209(4):579–593. pmid:26008745
- 60. Precious B, Carlos T, Goodbourn S, Randall R. Catalytic turnover of STAT1 allows PIV5 to dismantle the interferon-induced anti-viral state of cells. Virology. 2007;368(1):114–121. pmid:17640695
- 61. Wenta N, Strauss H, Meyer S, Vinkemeier U. Tyrosine phosphorylation regulates the partitioning of STAT1 between different dimer conformations. Proceedings of the National Academy of Sciences. 2008;105(27):9238–9243. pmid:18591661
- 62. Hottiger MO, Felzien LK, Nabel GJ. Modulation of cytokine-induced HIV gene expression by competitive binding of transcription factors to the coactivator p300. The EMBO Journal. 1998;17(11):3124–3134. pmid:9606194
- 63. Bachmann J, Raue A, Schilling M, Böhm ME, Kreutz C, Kaschek D, et al. Division of labor by dual feedback regulators controls JAK2/STAT5 signaling over broad ligand range. Molecular Systems Biology. 2011;7(1):516. pmid:21772264
- 64. Chatterjee-Kishore M, Wright KL, Ting JPY, Stark GR. How Stat1 mediates constitutive gene expression: a complex of unphosphorylated Stat1 and IRF1 supports transcription of the LMP2 gene. The EMBO journal. 2000;19(15):4111–4122. pmid:10921891
- 65. Kok F, Rosenblatt M, Teusel M, Nizharadze T, Gonçalves Magalhães V, Dächert C, et al. Disentangling molecular mechanisms regulating sensitization of interferon alpha signal transduction. Molecular Systems Biology. 2020;16(7):e8955. pmid:32696599
- 66. Kærn M, Elston TC, Blake WJ, Collins JJ. Stochasticity in gene expression: from theories to phenotypes. Nature Reviews Genetics. 2005;6(6):451–464. pmid:15883588
- 67. Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics. 1976;22(4):403–434.
- 68. Cendrowski J, Mamińska A, Miaczynska M. Endocytic regulation of cytokine receptor signaling. Cytokine & growth factor reviews. 2016;32:63–73. pmid:27461871
- 69. Darzacq X, Shav-Tal Y, De Turris V, Brody Y, Shenoy SM, Phair RD, et al. In vivo dynamics of RNA polymerase II transcription. Nature structural & molecular biology. 2007;14(9):796–806. pmid:17676063
- 70. Yang E, van Nimwegen E, Zavolan M, Rajewsky N, Schroeder M, Magnasco M, et al. Decay rates of human mRNAs: correlation with functional characteristics and sequence attributes. Genome research. 2003;13(8):1863–1872. pmid:12902380
- 71. Schwanhäusser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, et al. Global quantification of mammalian gene expression control. Nature. 2011;473(7347):337–42. pmid:21593866
- 72. Lee CK, Bluyssen HA, Levy DE. Regulation of interferon-α responsiveness by the duration of Janus kinase activity. Journal of Biological Chemistry. 1997;272(35):21872–21877. pmid:9268319
- 73. Siewert E, Müller-Esterl W, Starr R, Heinrich PC, Schaper F. Different protein turnover of interleukin-6-type cytokine signalling components. European Journal of Biochemistry. 1999;265(1):251–257. pmid:10491180
- 74. Feifei C, Ogawa K, Nagarajan RP, Zhang M, Kuang C, Yan c. Regulation of TG-interacting factor by transforming growth factor-beta. Biochemical Journal. 2003;371(2):257–263.
- 75. Branca A, Faltynek CR, D’Alessandro SB, Baglioni C. Interaction of interferon with cellular receptors. Internalization and degradation of cell-bound interferon. Journal of Biological Chemistry. 1982;257(22):13291–13296. pmid:6292184
- 76. Balakrishnan S, Mathad SS, Sharma G, Raju SR, Reddy UB, Das S, et al. A Nondimensional Model Reveals Alterations in Nuclear Mechanics upon Hepatitis C Virus Replication. Biophysical Journal. 2019;116(7):1328–1339. pmid:30879645
- 77. Kondo F, Wada K, Kondo Y. Morphometric analysis of hepatocellular carcinoma. Virchows Arch A Pathol Anat Histopathol. 1988;413(5):425–430. pmid:2845643
- 78. Gayotto LC, Scheuer PJ. Letter: Liver-cell mass and nuclear-cytoplasmic ratio in human liver. Journal of clinical pathology. 1975;28(7):599. pmid:1150901
- 79. Sainz B, Barretto N, Martin DN, Hiraga N, Imamura M, Hussain S, et al. Hepatitis C virus infection in phenotypically distinct Huh7 cell lines. PloS one. 2009;4(8). pmid:19668344
- 80. Kessler S, Hosseini K, Hussein U, Kim K, List M, Schultheiß C, et al. Hepatocellular Carcinoma and Nuclear Paraspeckles: Induction in Chemoresistance and Prediction for Poor Survival. Cell Physiol Biochem. 2019;52(4):787–801. pmid:30946555
- 81. Riccardi K, Ryu S, Lin J, Yates P, Tess D, Li R, et al. Comparison of Species and Cell-Type Differences in Fraction Unbound of Liver Tissues, Hepatocytes, and Cell Lines. Drug Metab Dispos. 2018;46:415–421. pmid:29437874
- 82.
Lodish H. Molecular cell biology. Macmillan; 2007.
- 83. Mitchell S, Roy K, Zangle TA, Hoffmann A. Nongenetic origins of cell-to-cell variability in B lymphocyte proliferation. Proceedings of the National Academy of Sciences. 2018;115(12):E2888–E2897. pmid:29514960
- 84. Munsky B, Trinh B, Khammash M. Listening to the noise: random fluctuations reveal gene network parameters. Molecular systems biology. 2009;5(1):318. pmid:19888213
- 85. Begitt A, Droescher M, Meyer T, Schmid CD, Baker M, Antunes F, et al. STAT1-cooperative DNA binding distinguishes type 1 from type 2 interferon signaling. Nature immunology. 2014;15(2):168–76. pmid:24413774
- 86. Ronni T, Melen K, Malygin A, Julkunen I. Control of IFN-inducible MxA gene expression in human cells. The Journal of Immunology. 1993;150(5):1715–1726. pmid:7679692
- 87. Pinto AK, Williams GD, Szretter KJ, White JP, Proença-Módena JL, Liu G, et al. Human and Murine IFIT1 Proteins Do Not Restrict Infection of Negative-Sense RNA Viruses of the Orthomyxoviridae, Bunyaviridae, and Filoviridae Families. Journal of Virology. 2015;89(18):9465–9476. pmid:26157117