Institute of Fundamental Technological Research

An overwhelming majority of mathematical models of regulatory pathways, including the intensively studied NF-κB pathway, remains non-identifiable, meaning that their parameters may not be determined by existing data. The existing NF-κB models that are capable of reproducing experimental data contain non-identifiable parameters, whereas simplified models with a smaller number of parameters exhibit dynamics that differs from that observed in experiments. Here, we reduced an existing model of the canonical NF-κB pathway by decreasing the number of equations from 15 to 6. The reduced model retains two negative feedback loops mediated by IκBα and A20, and in response to both tonic and pulsatile TNF stimulation exhibits dynamics that closely follow that of the original model. We carried out the sensitivity-based linear analysis and Monte Carlo-based analysis to demonstrate that the resulting model is both structurally and practically identifiable given measurements of 5 model variables from a simple TNF stimulation protocol. The reduced model is capable of reproducing different types of responses that are characteristic to regulatory motifs controlled by negative feedback loops: nearly-perfect adaptation as well as damped and sustained oscillations. It can serve as a building block of more comprehensive models of the immune response and cancer, where NF-κB plays a decisive role. Our approach, although may not be automatically generalized, suggests that models of other regulatory pathways can be transformed to identifiable, while retaining their dynamical features.

Resolving practical non-identifiability of computational models typically requires either additional data or non-algorithmic model reduction, which frequently results in models containing parameters lacking direct interpretation. Here, instead of reducing models, we explore an alternative, Bayesian approach, and quantify the predictive power of non-identifiable models. We considered an example biochemical signalling cascade model as well as its mechanical analogue. For these models, we demonstrated that by measuring a single variable in response to a properly chosen stimulation protocol, the dimensionality of the parameter space is reduced, which allows for predicting the measured variable’s trajectory in response to different stimulation protocols even if all model parameters remain unidentified. Moreover, one can predict how such a trajectory will transform in the case of a multiplicative change of an arbitrary model parameter. Successive measurements of remaining variables further reduce the dimensionality of the parameter space and enable new predictions. We analysed potential pitfalls of the proposed approach that can arise when the investigated model is oversimplified, incorrect, or when the training protocol is inadequate. The main advantage of the suggested iterative approach is that the predictive power of the model can be assessed and practically utilised at each step.

Type I interferons (IFNs) are key coordinators of the innate immune response to viral infection, which, through activation of the transcriptional regulators STAT1 and STAT2 (STAT1/2) in bystander cells, induce the expression of IFN-stimulated genes (ISGs). Here, we showed that in cells transfected with poly(I:C), an analog of viral RNA, the transcriptional activity of STAT1/2 was terminated because of depletion of the interferon-β (IFN-β) receptor, IFNAR. Activation of RNase L and PKR, products of two ISGs, not only hindered the replenishment of IFNAR but also suppressed negative regulators of IRF3 and NF-κB, consequently promoting IFNB transcription. We incorporated these findings into a mathematical model of innate immunity. By coupling signaling through the IRF3–NF-κB and STAT1/2 pathways with the activities of RNase L and PKR, the model explains how poly(I:C) switches the transcriptional program from being STAT1/2 induced to being IRF3 and NF-κB induced, which converts IFN-β–responding cells to IFN-β–secreting cells.

When infected with a virus, cells may secrete interferons (IFNs) that prompt nearby cells to prepare for upcoming infection. Reciprocally, viral proteins often interfere with IFN synthesis and IFN-induced signaling. We modeled the crosstalk between the propagating virus and the innate immune response using an agent-based stochastic approach. By analyzing immunofluorescence microscopy images we observed that the mutual antagonism between the respiratory syncytial virus (RSV) and infected A549 cells leads to dichotomous responses at the single-cell level and complex spatial patterns of cell signaling states. Our analysis indicates that RSV blocks innate responses at three levels: by inhibition of IRF3 activation, inhibition of IFN synthesis, and inhibition of STAT1/2 activation. In turn, proteins coded by IFN-stimulated (STAT1/2-activated) genes inhibit the synthesis of viral RNA and viral proteins. The striking consequence of these inhibitions is a lack of coincidence of viral proteins and IFN expression within single cells. The model enables investigation of the impact of immunostimulatory defective viral particles and signaling network perturbations that could potentially facilitate containment or clearance of the viral infection.

Living cells utilize signaling pathways to sense, transduce, and process information. As the extracellular stimulation often has rich temporal characteristics which may govern dynamic cellular responses, it is important to quantify the rate of information flow through the signaling pathways. In this study, we used an epithelial cell line expressing a light-activatable FGF receptor and an ERK activity reporter to assess the ability of the MAPK/ERK pathway to transduce signal encoded in a sequence of pulses. By stimulating the cells with random light pulse trains, we demonstrated that the MAPK/ERK channel capacity is at least 6 bits per hour. The input reconstruction algorithm detects the light pulses with 1-min accuracy 5 min after their occurrence. The high information transmission rate may enable the pathway to coordinate multiple processes including cell movement and respond to rapidly varying stimuli such as chemoattracting gradients created by other cells.

A cell constantly receives signals and takes different fates accordingly. Given the uncertainty rendered by signal transduction noise, a cell may incorrectly perceive these signals. It may mistakenly behave as if there is a signal, although there is none, or may miss the presence of a signal that actually exists. In this paper, we consider a signaling system with two outputs, and introduce and develop methods to model and compute key cell decision-making parameters based on the two outputs and in response to the input signal. In the considered system, the tumor necrosis factor (TNF) regulates the two transcription factors, the nuclear factor κB (NFκB) and the activating transcription factor-2 (ATF-2). These two system outputs are involved in important physiological functions such as cell death and survival, viral replication, and pathological conditions, such as autoimmune diseases and different types of cancer. Using the introduced methods, we compute and show what the decision thresholds are, based on the single-cell measured concentration levels of NFκB and ATF-2. We also define and compute the decision error probabilities, i.e., false alarm and miss probabilities, based on the concentration levels of the two outputs. By considering the joint response of the two outputs of the signaling system, one can learn more about complex cellular decision-making processes, the corresponding decision error rates, and their possible involvement in the development of some pathological conditions.

cell decision making, decision theory, molecular signaling systems, signal transduction noise, cellular decision error probabilities

Nucleosomes are recognized as key regulators of transcription. However, the relationship between slow nucleosome unwrapping dynamics and bulk transcriptional properties has not been thoroughly explored. Here, an agent-based model that we call the dynamic defect Totally Asymmetric Simple Exclusion Process (ddTASEP) was constructed to investigate the effects of nucleosome-induced pausing on transcriptional dynamics. Pausing due to slow nucleosome dynamics induced RNAPII convoy formation, which would cooperatively prevent nucleosome rebinding leading to bursts of transcription. The mean first passage time (MFPT) and the variance of first passage time (VFPT) were analytically expressed in terms of the nucleosome rate constants, allowing for the direct quantification of the effects of nucleosome-induced pausing on pioneering polymerase dynamics. The mean first passage elongation rate γ(hc, ho) is inversely proportional to the MFPT and can be considered to be a new axis of the ddTASEP phase diagram, orthogonal to the classical αβ-plane (where α and β are the initiation and termination rates). Subsequently, we showed that, for β = 1, there is a novel jamming transition in the αγ-plane that separates the ddTASEP dynamics into initiation-limited and nucleosome pausing-limited regions. We propose analytical estimates for the RNAPII density ρ, average elongation rate v, and transcription flux J and verified them numerically. We demonstrate that the intra-burst RNAPII waiting times tin follow the time-headway distribution of a max flux TASEP and that the average inter-burst interval correlates with the index of dispersion De. In the limit γ→0, the average burst size reaches a maximum set by the closing rate hc. When α≪1, the burst sizes are geometrically distributed, allowing large bursts even while the average burst size is small. Last, preliminary results on the relative effects of static and dynamic defects are presented to show that dynamic defects can induce equal or greater pausing than static bottle necks.

We observed the interference between two prevalent respiratory viruses, respiratory syncytial virus (RSV) and influenza A virus (IAV) (H1N1), and characterized its molecular underpinnings in alveolar epithelial cells (A549). We found that RSV induces higher levels of interferon beta (IFN-β) production than IAV and that IFN-β priming confers higher-level protection against infection with IAV than with RSV. Consequently, we focused on the sequential infection scheme of RSV and then IAV. Using A549 wild-type (WT), IFNAR1 knockout (KO), IFNLR1 KO, and IFNAR1-IFNLR1 double-KO cell lines, we found that both IFN-β and IFN-λ are necessary for maximum protection against subsequent infection. Immunostaining revealed that preinfection with RSV partitions the cell population into a subpopulation susceptible to subsequent infection with IAV and an IAV-proof subpopulation. Strikingly, the susceptible cells turned out to be those already compromised and efficiently expressing RSV, whereas the bystander, interferon-primed cells are resistant to IAV infection. Thus, virus-virus exclusion at the cell population level is not realized through direct competition for a shared ecological niche (single cell) but rather is achieved with the involvement of specific cytokines induced by the host’s innate immune response.

influenza A virus,innate immunity,interferon beta,interferon lambda,respiratory syncytial virus,single cell,viral interference

Persistent epidemic can lead to the emergence of new virus strains due to virus mutations. This work is devoted to the SIR model with strain-dependence of infected individuals due to virus mutations and a continuous strain variable. Characterization of epidemic progression is obtained for a strain-dependent infectivity function in numerical simulations and with some analytical estimates. Taking into account limited infection-induced immunity, transition from recovered to susceptible compartment is considered. Different scenarios of epidemic progression are identified. In the case of a monotonically growing transmission rate as a function of strain, only one epidemic outbreak is observed. However, if the transmission rate is a non-monotonic function, then multiple outbreaks can occur.

epidemic model, strain-dependent transmission rate, reinfection

Omicron, the novel highly mutated SARS-CoV-2 Variant of Concern (VOC, Pango lineage B.1.1.529) was first collected in early November 2021 in South Africa. By the end of November 2021, it had spread and approached fixation in South Africa, and had been detected on all continents. We analyzed the exponential growth of Omicron over four-week periods in the two most populated of South Africa’s provinces, Gauteng and KwaZulu-Natal, arriving at the doubling time estimates of, respectively, 3.3 days (95% CI: 3.2–3.4 days) and 2.7 days (95% CI: 2.3–3.3 days). Similar or even shorter doubling times were observed in other locations: Australia (3.0 days), New York State (2.5 days), UK (2.4 days), and Denmark (2.0 days). Log–linear regression suggests that the spread began in Gauteng around 11 October 2021; however, due to presumable stochasticity in the initial spread, this estimate can be inaccurate. Phylogenetics-based analysis indicates that the Omicron strain started to diverge between 6 October and 29 October 2021. We estimated that the weekly growth of the ratio of Omicron to Delta is in the range of 7.2–10.2, considerably higher than the growth of the ratio of Delta to Alpha (estimated to be in in the range of 2.5–4.2), and Alpha to pre-existing strains (estimated to be in the range of 1.8–2.7). High relative growth does not necessarily imply higher Omicron infectivity. A two-strain SEIR model suggests that the growth advantage of Omicron may stem from immune evasion, which permits this VOC to infect both recovered and fully vaccinated individuals. As we demonstrated within the model, immune evasion is more concerning than increased transmissibility, because it can facilitate larger epidemic outbreaks.

COVID-19 pandemic, SARS-CoV-2, Omicron variant, genome sequencing, mutation

Countries worldwide have adopted various strategies to minimize the socio-economic impact of the ongoing COVID-19 pandemic. Stringency of imposed measures universally reflects the standpoint from which protecting public health and avoiding damage to economy are seen as contradictory objectives. Based on epidemic trajectories of 25 highly developed countries and 10 US states in the (mobility reduction)–(reproduction number) plane we showed that delay in imposition of nation-wide quarantine elevates the number of infections and deaths, surge of which inevitably has to be suppressed by stringent and sustained lockdown. As a consequence, cumulative mobility reduction and population-normalized cumulative number of COVID-19-associated deaths are significantly correlated and this correlation increases with time. Overall, we demonstrated that, as long as epidemic suppression is the aim, the trade-off between the death toll and economic loss is illusory: high death toll correlates with deep and long-lasting lockdown causing a severe economic downturn.

The novel SARS-CoV-2 Variant of Concern (VOC)-202012/01 (also known as B.1.1.7), first collected in United Kingdom on 20 September 2020, is a rapidly growing lineage that in January 2021 constituted 86% of all SARS-CoV-2 genomes sequenced in England. The VOC has been detected in 40 out of 46 countries that reported at least 50 genomes in January 2021. We have estimated that the replicative advantage of the VOC is in the range 1.83–2.18 [95% CI: 1.71–2.40] with respect to the 20A.EU1 variant that dominated in England in November 2020, and in range 1.65–1.72 [95% CI: 1.46–2.04] in Wales, Scotland, Denmark, and USA. As the VOC strain will likely spread globally towards fixation, it is important to monitor its molecular evolution. We have estimated growth rates of expanding mutations acquired by the VOC lineage to find that the L18F substitution in spike has initiated a fast growing VOC substrain. The L18F substitution is of significance because it has been found to compromise binding of neutralizing antibodies. Of concern are immune escape mutations acquired by the VOC: E484K, F490S, S494P (in the receptor binding motif of spike) and Q677H, Q675H (in the proximity of the polybasic cleavage site at the S1/S2 boundary). These mutants may hinder efficiency of existing vaccines and expand in response to the increasing after-infection or vaccine-induced seroprevalence.

COVID-19 pandemic, SARS-CoV-2, spike protein, VOC-202012/01, spike L18F, genome sequencing, mutation

Alternations in the p53 regulatory network may render cancer cells resistant to the radiation-induced apoptosis. In this theoretical study we search for the best protocols combining targeted therapy with radiation to treat cancers with wild-type p53, but having downregulated expression of PTEN or overexpression of Wip1 resulting in resistance to radiation monotherapy. Instead of using the maximum tolerated dose paradigm, we exploit stochastic computational model of the p53 regulatory network to calculate apoptotic fractions for both normal and cancer cells. We consider combination protocols, with irradiations repeated every 12, 18, 24, or 36 h to find that timing between Mdm2 inhibitor delivery and irradiation significantly influences the apoptotic cell fractions. We assume that uptake of the inhibitor is higher by cancer than by normal cells and that cancer cells receive higher irradiation doses from intersecting beams. These two assumptions were found necessary for the existence of protocols inducing massive apoptosis in cancer cells without killing large fraction of normal cells neighboring tumor. The best found protocols have irradiations repeated every 24 or 36 h with two inhibitor doses per irradiation cycle, and allow to induce apoptosis in more than 95% of cancer cells, killing less than 10% of normal cells.

In this study, we established a dynamic micromodel of urinary tract infection to analyze the impact of UT-segment-specific urinary outflow on the persistence of E. coli colonization. We found that the adherence of Dr+ E. coli to bladder T24 transitional cells and type IV collagen is maximal at lowest shear stress and is reduced by any increase in flow velocity. The analyzed adherence was effective in the whole spectrum of physiological shear stress and was almost irreversible over the entire range of generated shear force. Once Dr+ E. coli bound to host cells or collagen, they did not detach even in the presence of elevated shear stress or of chloramphenicol, a competitive inhibitor of binding. Investigating the role of epithelial surface architecture, we showed that the presence of budding cells–a model microarchitectural obstacle–promotes colonization of the urinary tract by E. coli. We report a previously undescribed phenomenon of epithelial cell "rolling-shedding" colonization, in which the detached epithelial cells reattach to the underlying cell line through a layer of adherent Dr+ E. coli. This rolling-shedding colonization progressed continuously due to "refilling" induced by the flow-perturbing obstacle. The shear stress of fluid containing free-floating bacteria fueled the rolling, while providing an uninterrupted supply of new bacteria to be trapped by the rolling cell. The progressive rolling allows for transfer of briefly attached bacteria onto the underlying monolayer in a repeating cascading event.

The basic reproduction number R0 of the coronavirus disease 2019 has been estimated to range between 2 and 4. Here, we used an SEIR model that properly accounts for the distribution of the latent period and, based on empirical estimates of the doubling time in the near-exponential phases of epidemic progression in China, Italy, Spain, France, UK, Germany, Switzerland and New York State, we estimated that R0 lies in the range 4.7-11.4. We explained this discrepancy by performing stochastic simulations of model dynamics in a population with a small proportion of super-spreaders. The simulations revealed two-phase dynamics, in which an initial phase of relatively slow epidemic progression diverts to a faster phase upon appearance of infectious super-spreaders. Early estimates obtained for this initial phase may suggest lower R0.

COVID-19, reproduction number

Characterization of decision-making in cells in response to received signals is of importance for understanding how cell fate is determined. The problem becomes multi-faceted and complex when we consider cellular heterogeneity and dynamics of biochemical processes. In this paper, we present a unified set of decision-theoretic, machine learning and statistical signal processing methods and metrics to model the precision of signaling decisions, in the presence of uncertainty, using single cell data. First, we introduce erroneous decisions that may result from signaling processes and identify false alarms and miss events associated with such decisions. Then, we present an optimal decision strategy which minimizes the total decision error probability. Additionally, we demonstrate how graphing receiver operating characteristic curves conveniently reveals the trade-off between false alarm and miss probabilities associated with different cell responses. Furthermore, we extend the introduced framework to incorporate the dynamics of biochemical processes and reactions in a cell, using multi-time point measurements and multi-dimensional outcome analysis and decision-making algorithms. The introduced multivariate signaling outcome modeling framework can be used to analyze several molecular species measured at the same or different time instants. We also show how the developed binary outcome analysis and decision-making approach can be extended to more than two possible outcomes. As an example and to show how the introduced methods can be used in practice, we apply them to single cell data of PTEN, an important intracellular regulatory molecule in a p53 system, in wild-type and abnormal cells. The unified signaling outcome modeling framework presented here can be applied to various organisms ranging from viruses, bacteria, yeast and lower metazoans to more complex organisms such as mammalian cells. Ultimately, this signaling outcome modeling approach can be utilized to better understand the transition from physiological to pathological conditions such as inflammation, various cancers and autoimmune diseases.

Cell decision making, noise, decision theory, machine learning, signaling errors, p53 system

We analyze heteroclinic traveling waves propagating on two dimensional manifolds to show that the geometric modification of the front velocity is proportional to the geodesic curvature of the frontline. As a result, on surfaces of concave domains, stable standing fronts can be formed on lines of constant geodesic curvature. These lines minimize the geometric functional describing the system's energy, consisting of terms proportional to the front line-length and to the inclosed surface area. Front stabilization at portions of surface with negative Gaussian curvature, provides a mechanismof pattern formation. In contrast to the mechanism associated with the Turing instability, the proposed mechanism requires only a single scalar bistable reaction–diffusion equation and connects the intrinsic surface geometry with the arising pattern. By considering a system of equations modeling boundary-volume interactions, we show that polarization of the boundary may induce a corresponding polarization in the volume.

heteroclinic traveling waves, standing fronts, geodesic curvature, negative Gaussian curvature, domain polarization, pattern formation

This paper presents the development and comparison of label-free electrochemical immunosensors based on screen-printed gold and glassy carbon (GC) disc electrodes for efficient and rapid detection of respiratory syncytial virus (RSV). Briefly, the antibody specific to the F protein of RSV was successfully immobilized on modified electrodes. Antibody coupling on the Au surface was conducted via 4-aminothiophenol (4-ATP) and glutaraldehyde (GA). The GC surface was modified with poly-L-lysine (PLL) for direct anti-RSV conjugation after EDC/NHS (1-Ethyl-3-(3-dimethylaminopropyl)carbodiimide/N-Hydroxysuccinimide) activation. Electrochemical characterizations of the immunosensors were carried out by cyclic voltammetry (CV) and electrochemical impedance spectroscopy (EIS). GC-based immunosensors show a dynamic range of antigen detection from 1.0 × 10^5 PFU/mL to 1.5×10^7 PFU/mL, more than 1.0 × 10^5 PFU/mL to 1.0 × 10^7 PFU/mL for the Au-based sensor. However, the GC platform is less sensitive and shows a higher detection limit (LOD) for RSV. The limit of detection of the Au immunosensor is 1.1 × 10^3 PFU/mL, three orders of magnitude lower than 2.85 × 10^6 PFU/mL for GC. Thus, the Au-based immunosensor has better analytical performance for virus detection than a carbon-based platform due to high sensitivity and very low RSV detection, obtained with good reproducibility.

respiratory syncytial virus, cyclic voltammetry, electrochemical impedance spectroscopy, sensor, gold electrode, glassy carbon

We constructed a simple Susceptible−Exposed–Infectious–Removed model of the spread of COVID-19. The model is parametrised only by the average incubation period, τ, and two rate parameters: contact rate, β, and exclusion rate, γ. The rates depend on nontherapeutic interventions and determine the basic reproduction number, R0 = β/γ, and, together with τ, the daily multiplication coefficient in the early exponential phase, θ. Initial R0 determines the reduction of β required to contain the spread of the epidemic. We demonstrate that introduction of a cascade of multiple exposed states enables the model to reproduce the distributions of the incubation period and the serial interval reported by epidemiologists. Using the model, we consider a hypothetical scenario in which β is modulated solely by anticipated changes of social behaviours: first, β decreases in response to a surge of daily new cases, pressuring people to self-isolate, and then, over longer time scale, β increases as people gradually accept the risk. In this scenario, initial abrupt epidemic spread is followed by a plateau and slow regression, which, although economically and socially devastating, grants time to develop and deploy vaccine or at least limit daily cases to a manageable number.

basic reproduction number, novel coronavirus

In December 2019, the first case of infection with a new virus COVID-19 (SARS-CoV-2), named coronavirus, was reported in the city of Wuhan, China. At that time, almost nobody paid any attention to it. The new pathogen, however, fast proved to be extremely infectious and dangerous, resulting in about 3–5% mortality. Over the few months that followed, coronavirus has spread over entire world. At the end of March, the total number of infections is fast approaching the psychological threshold of one million, resulting so far in tens of thousands of deaths. Due to the high number of lives already lost and the virus high potential for further spread, and due to its huge overall impact on the economies and societies, it is widely admitted that coronavirus poses the biggest challenge to the humanity after the second World war. The COVID-19 epidemic is provoking numerous questions at all levels. It also shows that modern society is extremely vulnerable and unprepared to such events. A wide scientific and public discussion becomes urgent. Some possible directions of this discussion are suggested in this article.

COVID-19, epidemic progression, mathematical models, crisis management, open questions

We consider a singular limit of diffusion equations in 3D domains of thickness converging to zero. In the 2D limit the resulting reaction-diffusion equation has a source term resulting from the Robin-type boundary conditions imposed on boundaries of the original 3D domain. The proposed approach can be applied to constructing approximate solutions of diffusion problems in thin planar, cylindrical, or spherical layers between two membranes. As an example we refer to the problem of activation of B lymphocytes, which typically have large nuclei and a thin cytoplasmic layer which can be considered as a spherical shell. For this example, assuming additionally axial symmetry we provide a rigorous convergence theorem in the language of semigroups of operators.

semigroups of operators, degenerate convergence, singular perturbation, boundary conditions, thin layers, signaling pathways, phosphorylation–dephosphorylation cycle

Two important signalling pathways of NF-κB and ERK transmit merely 1 bit of information about the level of extracellular stimulation. It is thus unclear how such systems can coordinate complex cell responses to external cues. We analyse information transmission in the MAPK/ERK pathway that converts both constant and pulsatile EGF stimulation into pulses of ERK activity. Based on an experimentally verified computational model, we demonstrate that, when input consists of sequences of EGF pulses, transmitted information increases nearly linearly with time. Thus, pulse-interval transcoding allows more information to be relayed than the amplitude–amplitude transcoding considered previously for the ERK and NF-κB pathways. Moreover, the information channel capacity C, or simply bitrate, is not limited by the bandwidth B = 1/τ, where τ ≈ 1 h is the relaxation time. Specifically, when the input is provided in the form of sequences of short binary EGF pulses separated by intervals that are multiples of τ/n (but not shorter than τ), then for n = 2, C ≈ 1.39 bit/h^-1; and for n = 4, C ≈ 1.86 bit/h^-1. The capability to respond to random sequences of EGF pulses enables cells to propagate spontaneous ERK activity waves across tissue.

cellular signal transduction, pulsatile stimulation, pulse-interval transcoding, bandwidth, representation problem

The innate immune system processes pathogen-induced signals into cell fate decisions. How information is turned to decision remains unknown. By combining stochastic mathematical modelling and experimentation, we demonstrate that feedback interactions between the IRF3, NF-κB and STAT pathways lead to switch-like responses to a viral analogue, poly(I:C), in contrast to pulse-like responses to bacterial LPS. Poly(I:C) activates both IRF3 and NF-κB, a requirement for induction of IFNβ expression. Autocrine IFNβ initiates a JAK/STAT-mediated positive-feedback stabilising nuclear IRF3 and NF-κB in first responder cells. Paracrine IFNβ, in turn, sensitises second responder cells through a JAK/STAT-mediated positive feedforward pathway that upregulates the positive-feedback components: RIG-I, PKR and OAS1A. In these sensitised cells, the 'live-or-die' decision phase following poly(I:C) exposure is shorter—they rapidly produce antiviral responses and commit to apoptosis. The interlinked positive feedback and feedforward signalling is key for coordinating cell fate decisions in cellular populations restricting pathogen spread.

cellular signalling networks, innate immunity, regulatory networks, stochastic modelling

In bistable reaction–diffusion systems, transitions between stable states typically occur on timescales orders of magnitude longer than the chemical equilibration time. Estimation of transition rates within explicit Brownian dynamics simulations is computationally prohibitively costly. We present a method that exploits a single trajectory, generated by a prior simulation of diffusive motions of molecules, to sample chemical kinetic processes on timescales several orders of magnitude longer than the duration of the diffusive trajectory. In this approach, we „loop” the diffusive trajectory by transferring chemical states of the molecules from the last to the first time step of the trajectory. Trajectory looping can be applied to enhance sampling of rare events in biochemical systems in which the number of reacting molecules is constant, as in cellular signal transduction pathways. As an example, we consider a bistable system of autophosphorylating kinases, for which we calculate state-to-state transition rates and traveling wave velocities. We provide an open-source implementation of the method.

Rule-based modeling is a powerful approach for studying biomolecular site dynamics. Here, we present SPATKIN, a general-purpose simulator for rule-based modeling in two spatial dimensions. The simulation algorithm is a lattice-based method that tracks Brownian motion of individual molecules and the stochastic firing of rule-defined reaction events. Because rules are used as event generators, the algorithm is network-free, meaning that it does not require to generate the complete reaction network implied by rules prior to simulation. In a simulation, each molecule (or complex of molecules) is taken to occupy a single lattice site that cannot be shared with another molecule (or complex). SPATKIN is capable of simulating a wide array of membrane-associated processes, including adsorption, desorption and crowding. Models are specified using an extension of the BioNetGen language, which allows to account for spatial features of the simulated process. AVAILABILITY AND IMPLEMENTATION: The C ++ source code for SPATKIN is distributed freely under the terms of the GNU GPLv3 license. The source code can be compiled for execution on popular platforms (Windows, Mac and Linux). An installer for 64-bit Windows and a macOS app are available. The source code and precompiled binaries are available at the SPATKIN Web site (http://pmbm.ippt.pan.pl/software/spatkin).

In this study a new computational method is developed to quantify decision making errors in cells, caused by noise and signaling failures. Analysis of tumor necrosis factor (TNF) signaling pathway which regulates the transcription factor Nuclear Factor κB (NF-κB) using this method identifies two types of incorrect cell decisions called false alarm and miss. These two events represent, respectively, declaring a signal which is not present and missing a signal that does exist. Using single cell experimental data and the developed method, we compute false alarm and miss error probabilities in wild-type cells and provide a formulation which shows how these metrics depend on the signal transduction noise level. We also show that in the presence of abnormalities in a cell, decision making processes can be significantly affected, compared to a wild-type cell, and the method is able to model and measure such effects. In the TNF—NF-κB pathway, the method computes and reveals changes in false alarm and miss probabilities in A20-deficient cells, caused by cell’s inability to inhibit TNF-induced NF-κB response. In biological terms, a higher false alarm metric in this abnormal TNF signaling system indicates perceiving more cytokine signals which in fact do not exist at the system input, whereas a higher miss metric indicates that it is highly likely to miss signals that actually exist. Overall, this study demonstrates the ability of the developed method for modeling cell decision making errors under normal and abnormal conditions, and in the presence of transduction noise uncertainty. Compared to the previously reported pathway capacity metric, our results suggest that the introduced decision error metrics characterize signaling failures more accurately. This is mainly because while capacity is a useful metric to study information transmission in signaling pathways, it does not capture the overlap between TNF-induced noisy response curves.

Decision making, Radar, Probability distribution, Transcription factors, Signal processing, Signal transduction, Signaling networks, Statistical signal processing

Pattern formation is one of the most fundamental yet puzzling phenomena in physics and biology. We propose that traveling front pinning into concave portions of the boundary of 3-dimensional domains can serve as a generic gradient-maintaining mechanism. Such a mechanism of domain polarization arises even for scalar bistable reaction-diffusion equations, and, depending on geometry, a number of stationary fronts may be formed leading to complex spatial patterns. The main advantage of the pinning mechanism, with respect to the Turing bifurcation, is that it allows for maintaining gradients in the specific regions of the domain. By linking the instant domain shape with the spatial pattern, the mechanism can be responsible for cellular polarization and differentiation.

We formulated a computational model for a MAPK signaling cascade downstream of the EGF receptor to investigate how interlinked positive and negative feedback loops process EGF signals into ERK pulses of constant amplitude but dose-dependent duration and frequency. A positive feedback loop involving RAS and SOS, which leads to bistability and allows for switch-like responses to inputs, is nested within a negative feedback loop that encompasses RAS and RAF, MEK, and ERK that inhibits SOS via phosphorylation. This negative feedback, operating on a longer time scale, changes switch-like behavior into oscillations having a period of 1 hour or longer. Two auxiliary negative feedback loops, from ERK to MEK and RAF, placed downstream of the positive feedback, shape the temporal ERK activity profile but are dispensable for oscillations. Thus, the positive feedback introduces a hierarchy among negative feedback loops, such that the effect of a negative feedback depends on its position with respect to the positive feedback loop. Furthermore, a combination of the fast positive feedback involving slow-diffusing membrane components with slower negative feedbacks involving faster diffusing cytoplasmic components leads to local excitation/global inhibition dynamics, which allows the MAPK cascade to transmit paracrine EGF signals into spatially non-uniform ERK activity pulses.

MAPK signaling, oscillations, mathematical modelling

The NF-κB pathway is known to transmit merely 1 bit of information about stimulus level. We combined experimentation with mathematical modeling to elucidate how information about TNF concentration is turned into a binary decision. Using Kolmogorov-Smirnov distance, we quantified the cell's ability to discern 8 TNF concentrations at each step of the NF-κB pathway, to find that input discernibility decreases as signal propagates along the pathway. Discernibility of low TNF concentrations is restricted by noise at the TNF receptor level, whereas discernibility of high TNF concentrations it is restricted by saturation/depletion of downstream signaling components. Consequently, signal discernibility is highest between 0.03 and 1 ng/ml TNF. Simultaneous exposure to TNF or LPS and a translation inhibitor, cycloheximide, leads to prolonged NF-κB activation and a marked increase of transcript levels of NF-κB inhibitors, IκBα and A20. The impact of cycloheximide becomes apparent after the first peak of nuclear NF-κB translocation, meaning that the NF-κB network not only relays 1 bit of information to coordinate the all-or-nothing expression of early genes, but also over a longer time course integrates information about other stimuli. The NF-κB system should be thus perceived as a feedback-controlled decision-making module rather than a simple information transmission channel.

cellular signaling networks, innate immunity, stress signaling

Downstream of growth factor receptors and of the guanine triphosphatase (GTPase) RAS, heterodimers of the serine/threonine kinases BRAF and RAF1 are critical upstream kinases and activators of the mitogen-activated protein kinase (MAPK) module containing the mitogen-activated and extracellular signal–regulated kinase kinase (MEK) and their targets, the extracellular signal–regulated kinase (ERK) family. Either direct or scaffold protein–mediated interactions among the components of the ERK module (the MAPKKKs BRAF and RAF1, MEK, and ERK) facilitate signal transmission. RAF1 also has essential functions in the control of tumorigenesis and migration that are mediated through its interaction with the kinase ROKα, an effector of the GTPase RHO and regulator of cytoskeletal rearrangements. We combined mutational and kinetic analysis with mathematical modeling to show that the interaction of RAF1 with ROKα is coordinated with the role of RAF1 in the ERK pathway. We found that the phosphorylated form of RAF1 that interacted with and inhibited ROKα was generated during the interaction of RAF1 with the ERK module. This mechanism adds plasticity to the ERK pathway, enabling signal diversification at the level of both ERK and RAF. Furthermore, by connecting ERK activation with the regulation of ROKα and cytoskeletal rearrangements by RAF1, this mechanism has the potential to precisely coordinate the proper timing of proliferation with changes in cell shape, adhesion, or motility.

MAPK pathway, kinase RAF, protein isoform, phosphorylation, mathematical modeling

Background: In favorable conditions bacterial doubling time is less than 20 min, shorter than DNA replication time. In E. coli a single round of genome replication lasts about 40 min and it must be accomplished about 20 min before cell division. To achieve such fast growth rates bacteria perform multiple replication rounds simultaneously. As a result, when the division time is as short as 20 min E. coli has about 8 copies of origin of replication (ori) and the average copy number of the genes situated close to ori can be 4 times larger than those near the terminus of replication (ter). It implies that shortening of cell cycle may influence dynamics of regulatory pathways involving genes placed at distant loci. Results: We analyze this effect in a model of a genetic toggle switch, i.e. a system of two mutually repressing genes, one localized in the vicinity of ori and the other localized in the vicinity of ter. Using a stochastic model that accounts for cell growth and divisions we demonstrate that shortening of the cell cycle can induce switching of the toggle to the state in which expression of the gene placed near ter is suppressed. The toggle bistability causes that the ratio of expression of the competing genes changes more than two orders of magnitude for a two-fold change of the doubling time. The increasing stability of the two toggle states enhances system sensitivity but also its reaction time. Conclusions: By fusing the competing genes with fluorescent tags this mechanism could be tested and employed to create an indicator of the doubling time. By manipulating copy numbers of the competing genes and locus of the gene situated near ter, one can obtain equal average expression of both genes for any doubling time T between 20 and 120 min. Such a toggle would accurately report departures of the doubling time from T.

Mathematical modeling, Stochastic simulations, Regulatory pathways, Bistability, DNA replication, Gene copy number

The p53 transcription factor is a regulator of key cellular processes including DNA repair, cell cycle arrest, and apoptosis. In this theoretical study, we investigate how the complex circuitry of the p53 network allows for stochastic yet unambiguous cell fate decision-making. The proposed Markov chain model consists of the regulatory core and two subordinated bistable modules responsible for cell cycle arrest and apoptosis. The regulatory core is controlled by two negative feedback loops (regulated by Mdm2 and Wip1) responsible for oscillations, and two antagonistic positive feedback loops (regulated by phosphatases Wip1 and PTEN) responsible for bistability. By means of bifurcation analysis of the deterministic approximation we capture the recurrent solutions (i.e., steady states and limit cycles) that delineate temporal responses of the stochastic system. Direct switching from the limit-cycle oscillations to the “apoptotic” steady state is enabled by the existence of a subcritical Neimark—Sacker bifurcation in which the limit cycle loses its stability by merging with an unstable invariant torus. Our analysis provides an explanation why cancer cell lines known to have vastly diverse expression levels of Wip1 and PTEN exhibit a broad spectrum of responses to DNA damage: from a fast transition to a high level of p53 killer (a p53 phosphoform which promotes commitment to apoptosis) in cells characterized by high PTEN and low Wip1 levels to long-lasting p53 level oscillations in cells having PTEN promoter methylated (as in, e.g., MCF-7 cell line).

Apoptosis, Cell cycle and cell division, DNA damage, DNA repair, Phosphorylation, Biochemical simulations, Cell cycle inhibitors, Transcription factors

BACKGROUND:
Importins and exportins influence gene expression by enabling nucleocytoplasmic shuttling of transcription factors. A key transcription factor of innate immunity, NF-κB, is sequestered in the cytoplasm by its inhibitor, IκBα, which masks nuclear localization sequence of NF-κB. In response to TNFα or LPS, IκBα is degraded, which allows importins to bind NF-κB and shepherd it across nuclear pores. NF-κB nuclear activity is terminated when newly synthesized IκBα enters the nucleus, binds NF-κB and exportin which directs the complex to the cytoplasm. Although importins/exportins are known to regulate spatiotemporal kinetics of NF-κB and other transcription factors governing innate immunity, the mechanistic details of these interactions have not been elucidated and mathematically modelled.
RESULTS:
Based on our quantitative experimental data, we pursue NF-κB system modelling by explicitly including NF-κB-importin and IκBα-exportin binding to show that the competition between importins and IκBα enables NF-κB nuclear translocation despite high levels of IκBα. These interactions reduce the effective relaxation time and allow the NF-κB regulatory pathway to respond to recurrent TNFα pulses of 45-min period, which is about twice shorter than the characteristic period of NF-κB oscillations. By stochastic simulations of model dynamics we demonstrate that randomly appearing, short TNFα pulses can be converted to essentially digital pulses of NF-κB activity, provided that intervals between input pulses are not shorter than 1 h.
CONCLUSIONS:
By including interactions involving importin-α and exportin we bring the modelling of spatiotemporal kinetics of transcription factors to a more mechanistic level. Basing on the analysis of the pursued model we estimated the information transmission rate of the NF-κB pathway as 1 bit per hour.

Karyopherins, Nucleocytoplasmic transport, Negative feedback, Channel information capacity, Mathematical modelling

Digital signaling enhances robustness of cellular decisions in noisy environments, but it is unclear how digital systems transmit temporal information about a stimulus. To understand how temporal input information is encoded and decoded by the NF-κB system, we studied transcription factor dynamics and gene regulation under dose- and duration-modulated inflammatory inputs. Mathematical modeling predicted and microfluidic single-cell experiments confirmed that integral of the stimulus (or area, concentration × duration) controls the fraction of cells that activate NF-κB in the population. However, stimulus temporal profile determined NF-κB dynamics, cell-to-cell variability, and gene expression phenotype. A sustained, weak stimulation lead to heterogeneous activation and delayed timing that is transmitted to gene expression. In contrast, a transient, strong stimulus with the same area caused rapid and uniform dynamics. These results show that digital NF-κB signaling enables multidimensional control of cellular phenotype via input profile, allowing parallel and independent control of single-cell activation probability and population heterogeneity.

We constructed a mechanistic, computational model for regulation of (macro)autophagy and protein synthesis (at the level of translation). The model was formulated to study the system-level consequences of interactions among the following proteins: two key components of MTOR complex 1 (MTORC1), namely the protein kinase MTOR (mechanistic target of rapamycin) and the scaffold protein RPTOR; the autophagy-initiating protein kinase ULK1; and the multimeric energy-sensing AMP-activated protein kinase (AMPK). Inputs of the model include intrinsic AMPK kinase activity, which is taken as an adjustable surrogate parameter for cellular energy level or AMP:ATP ratio, and rapamycin dose, which controls MTORC1 activity. Outputs of the model include the phosphorylation level of the translational repressor EIF4EBP1, a substrate of MTORC1, and the phosphorylation level of AMBRA1 (activating molecule in BECN1-regulated autophagy), a substrate of ULK1 critical for autophagosome formation. The model incorporates reciprocal regulation of mTORC1 and ULK1 by AMPK, mutual inhibition of MTORC1 and ULK1, and ULK1-mediated negative feedback regulation of AMPK. Through analysis of the model, we find that these processes may be responsible, depending on conditions, for graded responses to stress inputs, for bistable switching between autophagy and protein synthesis, or relaxation oscillations, comprising alternating periods of autophagy and protein synthesis. A sensitivity analysis indicates that the prediction of oscillatory behavior is robust to changes of the parameter values of the model. The model provides testable predictions about the behavior of the AMPK-MTORC1-ULK1 network, which plays a central role in maintaining cellular energy and nutrient homeostasis.

Biological signals in cells are transmitted with the use of reaction cycles, such as the phosphorylation-dephosphorylation cycle, in which substrate is modified by antagonistic enzymes. An appreciable share of such reactions takes place in crowded environments of two-dimensional structures, such as plasma membrane or intracellular membranes, and is expected to be diffusion-controlled. In this work, starting from the microscopic bimolecular reaction rate constants and using estimates of the mean first-passage time for an enzyme–substrate encounter, we derive diffusion-dependent effective macroscopic reaction rate coefficients (EMRRC) for a generic reaction cycle. Each EMRRC was found to be half of the harmonic average of the microscopic rate constant (phosphorylation c or dephosphorylation d), and the effective (crowding-dependent) motility divided by a slowly decreasing logarithmic function of the sum of the enzyme concentrations. This implies that when c and d differ, the two EMRRCs scale differently with the motility, rendering the steady-state fraction of phosphorylated substrate molecules diffusion-dependent. Analytical predictions are verified using kinetic Monte Carlo simulations on the two-dimensional triangular lattice at the single-molecule resolution. It is demonstrated that the proposed formulas estimate the steady-state concentrations and effective reaction rates for different sets of microscopic reaction rates and concentrations of reactants, including a non-trivial example where with increasing diffusivity the fraction of phosphorylated substrate molecules changes from 10% to 90%.

Enzymes, Enzyme kinetics, Diffusion, Reaction rate constants, Membrane biochemistry

We investigate the kinetics of the ubiquitous phosphorylation-dephosphorylation cycle on biological membranes by means of kinetic Monte Carlo simulations on the triangular lattice. We establish the dependence of effective macroscopic reaction rate coefficients as well as the steady-state phosphorylated substrate fraction on the diffusion coefficient and concentrations of opposing enzymes: kinases and phosphatases. In the limits of zero and infinite diffusion, the numerical results agree with analytical predictions; these two limits give the lower and the upper bound for the macroscopic rate coefficients, respectively. In the zero-diffusion limit, which is important in the analysis of dense systems, phosphorylation and dephosphorylation reactions can convert only these substrates which remain in contact with opposing enzymes. In the most studied regime of nonzero but small diffusion, a contribution linearly proportional to the diffusion coefficient appears in the reaction rate. In this regime, the presence of opposing enzymes creates inhomogeneities in the (de)phosphorylated substrate distributions: The spatial correlation function shows that enzymes are surrounded by clouds of converted substrates. This effect becomes important at low enzyme concentrations, substantially lowering effective reaction rates. Effective reaction rates decrease with decreasing diffusion and this dependence is more pronounced for the less-abundant enzyme. Consequently, the steady-state fraction of phosphorylated substrates can increase or decrease with diffusion, depending on relative concentrations of both enzymes. Additionally, steady states are controlled by molecular crowders which, mostly by lowering the effective diffusion of reactants, favor the more abundant enzyme.

Background
Compared to engineering or physics problems, dynamical models in quantitative biology typically depend on a relatively large number of parameters. Progress in developing mathematics to manipulate such multi-parameter models and so enable their efficient interplay with experiments has been slow. Existing solutions are significantly limited by model size.

Results
In order to simplify analysis of multi-parameter models a method for clustering of model parameters is proposed. It is based on a derived statistically meaningful measure of similarity between groups of parameters. The measure quantifies to what extend changes in values of some parameters can be compensated by changes in values of other parameters. The proposed methodology provides a natural mathematical language to precisely communicate and visualise effects resulting from compensatory changes in values of parameters. As a results, a relevant insight into identifiability analysis and experimental planning can be obtained. Analysis of NF- κB and MAPK pathway models shows that highly compensative parameters constitute clusters consistent with the network topology. The method applied to examine an exceptionally rich set of published experiments on the NF- κB dynamics reveals that the experiments jointly ensure identifiability of only 60 % of model parameters. The method indicates which further experiments should be performed in order to increase the number of identifiable parameters.

Conclusions
We currently lack methods that simplify broadly understood analysis of multi-parameter models. The introduced tools depict mutually compensative effects between parameters to provide insight regarding role of individual parameters, identifiability and experimental design. The method can also find applications in related methodological areas of model simplification and parameters estimation.

We present an integrated dynamical cross-talk model of the epithelial innate immune reponse (IIR) incorporating RIG-I and TLR3 as the two major pattern recognition receptors (PRR) converging on the RelA and IRF3 transcriptional effectors. bioPN simulations reproduce biologically relevant gene-and protein abundance measurements in response to time course, gene silencing and dose-response perturbations both at the population and single cell level. Our computational predictions suggest that RelA and IRF3 are under auto- and cross-regulation. We predict, and confirm experimentally, that RIG-I mRNA expression is controlled by IRF7. We also predict the existence of a TLR3-dependent, IRF3-independent transcription factor (or factors) that control(s) expression of MAVS, IRF3 and members of the IKK family. Our model confirms the observed dsRNA dose-dependence of oscillatory patterns in single cells, with periods of 1–3 hr. Model fitting to time series, matched by knockdown data suggests that the NF-κB module operates in a different regime (with different coefficient values) than in the TNFα-stimulation experiments. In future studies, this model will serve as a foundation for identification of virus-encoded IIR antagonists and examination of stochastic effects of viral replication.

Our model generates simulated time series, which reproduce the noisy oscillatory patterns of activity (with 1–3 hour period) observed in individual cells. Our work supports the hypothesis that the IIR is a phenomenon that emerged by evolution despite highly variable responses at an individual cell level.

Bistable regulatory elements are important for nongenetic inheritance, increase of cell-to-cell heterogeneity allowing adaptation, and robust responses at the population level. Here, we study computationally the bistable genetic toggle switch—a small regulatory network consisting of a pair of mutual repressors—in growing and dividing bacteria. We show that as cells with an inhibited growth exhibit high stability of toggle states, cell growth and divisions lead to a dramatic increase of toggling rates. The toggling rates were found to increase with rate of cell growth, and can be up to six orders of magnitude larger for fast growing cells than for cells with the inhibited growth. The effect is caused mainly by the increase of protein and mRNA burst sizes associated with the faster growth. The observation that fast growth dramatically destabilizes toggle states implies that rapidly growing cells may vigorously explore the epigenetic landscape enabling nongenetic evolution, while cells with inhibited growth adhere to the local optima. This can be a clever population strategy that allows the slow growing (but stress resistant) cells to survive long periods of unfavorable conditions. Simultaneously, at favorable conditions, this stress resistant (but slowly growing—or not growing) subpopulation may be replenished due to a high switching rate from the fast growing population.

Gene expression, Bistability, Stochastic processes, Genetic toggle switch, Cell growth and division

NF-κB is a key transcription factor that regulates innate immune response. Its activity is tightly controlled by numerous feedback loops, including two negative loops mediated by NF-κB inducible inhibitors, IκBα and A20, which assure oscillatory responses, and by positive feedback loops arising due to the paracrine and autocrine regulation via TNFα, IL-1 and other cytokines. We study the NF-κB system of interlinked negative and positive feedback loops, combining bifurcation analysis of the deterministic approximation with stochastic numerical modeling. Positive feedback assures the existence of limit cycle oscillations in unstimulated wild-type cells and introduces bistability in A20-deficient cells. We demonstrated that cells of significant autocrine potential, i.e., cells characterized by high secretion of TNFα and its receptor TNFR1, may exhibit sustained cytoplasmic–nuclear NF-κB oscillations which start spontaneously due to stochastic fluctuations. In A20-deficient cells even a small TNFα expression rate qualitatively influences system kinetics, leading to long-lasting NF-κB activation in response to a short-pulsed TNFα stimulation. As a consequence, cells with impaired A20 expression or increased TNFα secretion rate are expected to have elevated NF-κB activity even in the absence of stimulation. This may lead to chronic inflammation and promote cancer due to the persistent activation of antiapoptotic genes induced by NF-κB. There is growing evidence that A20 mutations correlate with several types of lymphomas and elevated TNFα secretion is characteristic of many cancers. Interestingly, A20 loss or dysfunction also leaves the organism vulnerable to septic shock and massive apoptosis triggered by the uncontrolled TNFα secretion, which at high levels overcomes the antiapoptotic action of NF-κB. It is thus tempting to speculate that some cancers of deregulated NF-κB signaling may be prone to the pathogen-induced apoptosis.

Transitions between steady states of a multi-stable stochastic system in the perfectly mixed chemical reactor are possible only because of stochastic switching. In realistic cellular conditions, where diffusion is limited, transitions between steady states can also follow from the propagation of travelling waves. Here, we study the interplay between the two modes of transition for a prototype bistable system of kinase–phosphatase interactions on the plasma membrane. Within microscopic kinetic Monte Carlo simulations on the hexagonal lattice, we observed that for finite diffusion the behaviour of the spatially extended system differs qualitatively from the behaviour of the same system in the well-mixed regime. Even when a small isolated subcompartment remains mostly inactive, the chemical travelling wave may propagate, leading to the activation of a larger compartment. The activating wave can be induced after a small subdomain is activated as a result of a stochastic fluctuation. Such a spontaneous onset of activity is radically more probable in subdomains characterized by slower diffusion. Our results show that a local immobilization of substrates can lead to the global activation of membrane proteins by the mechanism that involves stochastic fluctuations followed by the propagation of a semi-deterministic travelling wave.

multi-stability, Markov process, spatially extended system, kinetic Monte Carlo on the lattice, kinase autophosphorylation, cell signalling

We perform a mathematical analysis of a spatially extended model describing mutual phosphorylation of cytosolic kinases and membrane receptors. The analyzed regulatory system is a part of signal transduction mechanisms, which enables communication of the cell with its extracellular environment or other cells. The mutual receptor-kinase interaction is characteristic for immune receptors and Src family kinases. From the mathematical viewpoint, the considered system is interesting because it couples differential equations defined in a domain $\Omega$ and on its boundary $\partial\Omega$ via nonlinear Robin boundary conditions. Assuming a spherically symmetric framework, our approach is to consider an auxiliary problem in which the Robin boundary condition on the external boundary of the spherical shell $\Omega$ is replaced by a uniform Dirichlet boundary condition. This method allows us to find the stationary spherically symmetric solutions, both stable and unstable. Interestingly, numerical computations suggest also the existence of nonspherically symmetric unstable stationary solutions to the spherically symmetric problem. These conjectured solutions appear to lie between super- and subsolutions that converge in time to two different stable spherically symmetric steady states. The study is completed by the proof of an existence theorem for a general version of the original system encompassing the earlier models. The theorem holds in domains of arbitrary shape, with a number of holes that may represent various organelles impenetrable to the considered kinases. The last result ensures that the problem in which the flux of active kinases (which replaces the source term) is determined by the Robin-type boundary condition at the cell membrane is well posed. The considered generalized model may thus serve as a template for intracellular signal transduction analysis.

reaction-diffusion systems, Robin-type boundary conditions, bistability, molecular signal transduction

The aim of this study is to demonstrate that in molecular dynamical systems with the underlying bi- or multistability, the type of noise determines the most strongly attracting steady state or stochastic attractor. As an example we consider a simple stochastic model of autoregulatory gene with a nonlinear positive feedback, which in the deterministic approximation has two stable steady state solutions. Three types of noise are considered: transcriptional and translational – due to the small number of gene product molecules and the gene switching noise – due to gene activation and inactivation transitions. We demonstrate that the type of noise in addition to the noise magnitude dictates the allocation of probability mass between the two stable steady states. In particular, we found that when the gene switching noise dominates over the transcriptional and translational noise (which is characteristic of eukaryotes), the gene preferentially activates, while in the opposite case, when the transcriptional noise dominates (which is characteristic of prokaryotes) the gene preferentially remains inactive. Moreover, even in the zero-noise limit, when the probability mass generically concentrates in the vicinity of one of two steady states, the choice of the most strongly attracting steady state is noise type-dependent. Although the epigenetic attractors are defined with the aid of the deterministic approximation of the stochastic regulatory process, their relative attractivity is controlled by the type of noise, in addition to noise magnitude. Since noise characteristics vary during the cell cycle and development, such mode of regulation can be potentially employed by cells to switch between alternative epigenetic attractors.

Gene expression, Bistability, Stochastic processes, Epigenetic attractors

Background
Apoptosis is a tightly regulated process: cellular survive-or-die decisions cannot be accidental and must be unambiguous. Since the suicide program may be initiated in response to numerous stress stimuli, signals transmitted through a number of checkpoints have to be eventually integrated.

Results
In order to analyze possible mechanisms of the integration of multiple pro-apoptotic signals, we constructed a simple model of the Bcl-2 family regulatory module. The module collects upstream signals and processes them into life-or-death decisions by employing interactions between proteins from three subgroups of the Bcl-2 family: pro-apoptotic multidomain effectors, pro-survival multidomain restrainers, and pro-apoptotic single domain BH3-only proteins. Although the model is based on ordinary differential equations (ODEs), it demonstrates that the Bcl-2 family module behaves akin to a Boolean logic gate of the type dependent on levels of BH3-only proteins (represented by Bad) and restrainers (represented by Bcl-xL). A low level of pro-apoptotic Bad or a high level of pro-survival Bcl-xL implies gate AND, which allows for the initiation of apoptosis only when two stress stimuli are simultaneously present: the rise of the p53 killer level and dephosphorylation of kinase Akt. In turn, a high level of Bad or a low level of Bcl-xL implies gate OR, for which any of these stimuli suffices for apoptosis.

Conclusions
Our study sheds light on possible signal integration mechanisms in cells, and spans a bridge between modeling approaches based on ODEs and on Boolean logic. In the proposed scheme, logic gates switching results from the change of relative abundances of interacting proteins in response to signals and involves system bistability. Consequently, the regulatory system may process two analogous inputs into a digital survive-or-die decision.

Apoptosis, Cell survival, Signaling pathway, Bcl-2 family, Bistability, Boolean logic, Ordinary differential equations

The Raf/MEK/ERK cascade is one of the most studied and important signal transduction pathways. However, existing models largely ignore the existence of isoforms of the constituent kinases and their interactions. Here, we propose a model of the ERK cascade that includes heretofore neglected differences between isoforms of MEK. In particular, MEK1 is subject to a negative feedback from activated ERK, which is further conferred to MEK2 via hetero-dimerization. Specifically, ERK phosphorylates MEK1 at the residue Thr292, hypothetically creating an additional phosphatase binding site, accelerating MEK1 and MEK2 dephosphorylation. We incorporated these recently discovered interactions into a mathematical model of the ERK cascade that reproduces the experimental results of Catalanotti et al (2009 Nature Struct. Mol. Biol. 16 294–303) and Kamioka et al (2010 J. Biol. Chem. 285 33540–8). Furthermore, the model allows for predictions regarding the differences in the catalytic activity and function of the MEK isoforms. We propose that the MEK1/MEK2 ratio regulates the duration of the response, which increases with the level of MEK2 and decreases with the level of MEK1. In turn, the amplitude of the response is controlled by the total amount of the two isoforms. We confirm the proposed model structure performing a random parameter sampling, which led us to the conclusion that the sampled parameters, selected to properly reproduce wild-type (WT) cell behavior, to allow for qualitative reproduction of differences in behavior WT cells and cell mutants studied experimentally.

Bistable regulatory elements enhance heterogeneity in cell populations and, in multicellular organisms, allow cells to specialize and specify their fate. Our study demonstrates that in a system of bistable genetic switch, the noise characteristics control in which of the two epigenetic attractors the cell population will settle. We focus on two types of noise: the gene switching noise and protein dimerization noise. We found that the change of magnitudes of these noise components for one of the two competing genes introduces a large asymmetry of the protein stationary probability distribution and changes the relative probability of individual gene activation. Interestingly, an increase of noise associated with a given gene can either promote or suppress the activation of the gene, depending on the type of noise. Namely, each gene is repressed by an increase of its gene switching noise and activated by an increase of its protein-product dimerization noise. The observed effect was found robust to the large, up to fivefold deviations of the model parameters. In summary, we demonstrated that noise itself may determine the relative strength of the epigenetic attractors, which may provide a unique mode of control of cell fate decisions.

Gene expression, Bistability, Stochastic processes, Genetic toggle switch

BCR signaling regulates the activities and fates of B cells. BCR signaling encompasses two feedback loops emanating from Lyn and Fyn, which are Src family protein tyrosine kinases (SFKs). Positive feedback arises from SFK-mediated trans phosphorylation of BCR and receptor-bound Lyn and Fyn, which increases the kinase activities of Lyn and Fyn. Negative feedback arises from SFK-mediated cis phosphorylation of the transmembrane adapter protein PAG1, which recruits the cytosolic protein tyrosine kinase Csk to the plasma membrane, where it acts to decrease the kinase activities of Lyn and Fyn. To study the effects of the positive and negative feedback loops on the dynamical stability of BCR signaling and the relative contributions of Lyn and Fyn to BCR signaling, we consider in this study a rule-based model for early events in BCR signaling that encompasses membrane-proximal interactions of six proteins, as follows: BCR, Lyn, Fyn, Csk, PAG1, and Syk, a cytosolic protein tyrosine kinase that is activated as a result of SFK-mediated phosphorylation of BCR. The model is consistent with known effects of Lyn and Fyn deletions. We find that BCR signaling can generate a single pulse or oscillations of Syk activation depending on the strength of Ag signal and the relative levels of Lyn and Fyn. We also show that bistability can arise in Lyn- or Csk-deficient cells.

The MAPK cascades are principal kinase transduction pathways in eukaryotic cells. This family includes RAF/ERK, JNK, and p38 pathways. In the MAPK cascade, the signal is transmitted through three layers of sequentially activated kinases, MAP3K, MAP2K, and MAPK. The latter two kinases require dual phosphorylation for activation. The dual phosphorylation requirement has been implicated in bringing about bistability and switch-like responses in the cascade. MAPK signaling has been known to involve scaffolds—multidomain proteins that can assemble protein complexes; in this case the three MAPK components. Scaffolds are thought to increase the specificity of signaling by pairing enzymes and substrates. Scaffolds have been shown to biphasically control the response (the level of activated MAPK) and amplify it at a certain scaffold concentration range. In order to understand the interplay of scaffolding and multisite phosphorylation, in this study we analyze simplified MAPK signaling models in which we assume that either mono- or double phosphorylation of MAP2K and MAPK is required for activation. We demonstrate that the requirement for double phosphorylation directs signaling through scaffolds. In the hypothetical scenario in which mono-phosphorylation suffices for kinase activity, the presence of scaffolds has little effect on the response. This suggests that double phosphorylation in MAPK pathways, although not as efficient as mono-phosphorylation, evolved together with scaffolds to assure the tighter control and higher specificity in signaling, by enabling scaffolds to function as response amplifiers.

Kinase cascades, Rule based modeling, Molecular pathways evolution

Living cells may be considered as biochemical reactors of multiple steady states. Transitions between these states are enabled by noise, or, in spatially extended systems, may occur due to the traveling wave propagation. We analyze a one-dimensional bistable stochastic birth–death process by means of potential and temperature fields. The potential is defined by the deterministic limit of the process, while the temperature field is governed by noise. The stable steady state in which the potential has its global minimum defines the global deterministic attractor. For the stochastic system, in the low noise limit, the stationary probability distribution becomes unimodal, concentrated in one of two stable steady states, defined in this study as the global stochastic attractor. Interestingly, these two attractors may be located in different steady states. This observation suggests that the asymptotic behavior of spatially extended stochastic systems depends on the substrate diffusivity and size of the reactor. We confirmed this hypothesis within kinetic Monte Carlo simulations of a bistable reaction–diffusion model on the hexagonal lattice. In particular, we found that although the kinase–phosphatase system remains inactive in a small domain, the activatory traveling wave may propagate when a larger domain is considered.

multi-stability, Markov process, spatially extended system, kinetic Monte Carlo on the lattice, cell signalling

We proposed a spatially extended model of early events of B cell receptors (BCR) activation, which is based on mutual kinase-receptor interactions that are characteristic for the immune receptors and the Src family kinases. These interactions lead to the positive feedback which, together with two nonlinearities resulting from the double phosphorylation of receptors and Michaelis-Menten dephosphorylation kinetics, are responsible for the system bistability. We demonstrated that B cell can be activated by a formation of a tiny cluster of receptors or displacement of the nucleus. The receptors and Src kinases are activated, first locally, in the locus of the receptor cluster or the region where the cytoplasm is the thinnest. Then the traveling wave of activation propagates until activity spreads over the whole cell membrane. In the models in which we assume that the kinases are free to diffuse in the cytoplasm, we found that the fraction of aggregated receptors, capable to initiate B cell activation decreases with the decreasing thickness of cytoplasm and decreasing kinase diffusion. When kinases are restricted to the cell membrane - which is the case for most of the Src family kinases - even a cluster consisting of a tiny fraction of total receptors becomes activatory. Interestingly, the system remains insensitive to the modest changes of total receptor level. The model provides a plausible mechanism of B cells activation due to the formation of small receptors clusters collocalized by binding of polyvalent antigens or arising during the immune synapse formation.

We demonstrate that a single reconnection of two quantum vortices can lead to the creation of a cascade of vortex rings. Our analysis involves localized induction approximation, high-resolution Biot-Savart and Gross-Pitaevskii simulations. The latter showed that the rings cascade starts on the atomic scale, with rings diameters orders of magnitude smaller than the characteristic line spacing in the tangle. Vortex rings created in the cascades may penetrate the tangle and annihilate on the boundaries. This provides an efficient decay mechanism for sparse or moderately dense vortex tangle at very low temperatures.

Rule-based modeling provides a means to represent cell signaling systems in a way that captures site-specific details of molecular interactions. For rule-based models to be more widely understood and (re)used, conventions for model visualization and annotation are needed. We have developed the concepts of an extended contact map and a model guide for illustrating and annotating rule-based models. An extended contact map represents the scope of a model by providing an illustration of each molecule, molecular component, direct physical interaction, post-translational modification, and enzyme–substrate relationship considered in a model. A map can also illustrate allosteric effects, structural relationships among molecular components, and compartmental locations of molecules. A model guide associates elements of a contact map with annotation and elements of an underlying model, which may be fully or partially specified. A guide can also serve to document the biological knowledge upon which a model is based. We provide examples of a map and guide for a published rule-based model that characterizes early events in IgE receptor (FceRI) signaling. We also provide examples of how to visualize a variety of processes that are common in cell signaling systems but not considered in the example model, such as ubiquitination. An extended contact map and an associated guide can document knowledge of a cell signaling system in a form that is visual as well as executable. As a tool for model annotation, a map and guide can communicate the content of a model clearly and with precision, even for large models.

B and Mast cells are activated by the aggregation of the immune receptors. Motivated by this phenomena we consider a simple spatially extended model of mutual interaction of kinases and membrane receptors. It is assumed that kinase activates membrane receptors and in turn the kinase molecules bound to the active receptors are activated by transphosphorylation. Such a type of interaction implies positive feedback and may lead to bistability. In this study we apply the Steklov eigenproblem theory to analyze the linearized model and find exact solutions in the case of non-uniformly distributed membrane receptors. This approach allows us to determine the critical value of receptor dephosphorylation rate at which cell activation (by arbitrary small perturbation of the inactive state) is possible. We found that cell sensitivity grows with decreasing kinase diffusion and increasing anisotropy of the receptor distribution. Moreover, these two effects are cooperating. We showed that the cell activity can be abruptly triggered by the formation of the receptor aggregate. Since the considered activation mechanism is not based on receptor crosslinking by polyvalent antigens, the proposed model can also explain B cell activation due to receptor aggregation following binding of monovalent antigens presented on the antigen presenting cell.

Cells operate in dynamic environments using extraordinary communication capabilities that emerge from the interactions of genetic circuitry. The mammalian immune response is a striking example of the coordination of different cell types1. Cell-to-cell communication is primarily mediated by signalling molecules that form spatiotemporal concentration gradients, requiring cells to respond to a wide range of signal intensities2. Here we use high-throughput microfluidic cell culture3 and fluorescence microscopy, quantitative gene expression analysis and mathematical modelling to investigate how single mammalian cells respond to different concentrations of the signalling molecule tumour-necrosis factor (TNF)-α, and relay information to the gene expression programs by means of the transcription factor nuclear factor (NF)-κB. We measured NF-κB activity in thousands of live cells under TNF-α doses covering four orders of magnitude. We find, in contrast to population-level studies with bulk assays2, that the activation is heterogeneous and is a digital process at the single-cell level with fewer cells responding at lower doses. Cells also encode a subtle set of analogue parameters to modulate the outcome; these parameters include NF-κB peak intensity, response time and number of oscillations. We developed a stochastic mathematical model that reproduces both the digital and analogue dynamics as well as most gene expression profiles at all measured conditions, constituting a broadly applicable model for TNF-α-induced NF-κB signalling in various types of cells. These results highlight the value of high-throughput quantitative measurements with single-cell resolution in understanding how biological systems operate.

Cell biology, Biophysics, Immunology, Genetics, Genomics

The spatiotemporal kinetics of proteins and other substrates regulate cell fate and signaling. In this study, we consider a reaction–diffusion model of interaction of membrane receptors with a two-step kinase cascade. The receptors activate the ‘up-stream’ kinase, which may diffuse over cell volume and activate the ‘down-stream’ kinase, which is also diffusing. Both kinase species and receptors are inactivated by uniformly distributed phosphatases. The positive feedback, key to the considered dynamics, arises since the up-stream kinase activates the receptors. Such a mutual interaction is characteristic for immune cell receptors. Based on the proposed model, we demonstrated that cell sensitivity (measured as a critical value of phosphatase activity at which cell maybe activated) increases with decreasing motility of receptor-interacting kinases and with increasing polarity of receptors distribution. These two effects are cooperating, the effect of receptors localisation close to one pole of the cell grows with the decreasing kinase diffusion and vanishes in the infinite diffusion limit. As the cell sensitivity increases with decreasing diffusion of receptor-interacting kinase, the overall activity of the down-stream kinase increases with its diffusion. In conclusion, the analysis of the proposed model shows that, for the fixed substrate interaction rates, spatial distribution of the surface receptors together with the motility of intracellular kinases control cell signalling and sensitivity to extracellular signals. The increase of the cell sensitivity can be achieved by (i) localisation of receptors in a small subdomain of the cell membrane, (ii) lowering the motility of receptor-interacting kinase, (iii) increasing the motility of down-stream kinases which distribute the signal over the whole cell.

In living cells proteins motilities regulate the spatiotemporal dynamics of molecular pathways. We consider here a reaction–diffusion model of mutual kinase–receptor activation showing that the strength of positive feedback is controlled by the kinase diffusion coefficient. For high diffusion, the activated kinase molecules quickly leave the vicinity of the cell membrane and cannot efficiently activate the receptors. As a result, in a broad range of parameters, the cell can be activated only if the kinase diffusion coefficient is sufficiently small. Our simple model shows that change in the motility of substrates may dramatically influence the cell responses.

Reaction–diffusion system, Signal transduction, Positive feedback, Kinase activation, Membrane receptors, Protein motility

A number of regulatory networks have the potential to generate sustained oscillations of irregular amplitude, but well conserved period. Single-cell experiments revealed that in p53 and nuclear factor (NF)-kB systems the oscillation period is homogenous in cell populations, insensitive to the strength of the stimulation, and is not influenced by the overexpression of p53 or NF-kB transcription factors. We propose a novel computational method of validation of molecular pathways models, based on the analysis of sensitivity of the oscillation period to the particular gene(s) copy number and the level of stimulation. Using this method, the authors demonstrate that existing p53 models, in which oscillations are borne at a saddle-node-on-invariantcircle or subcritical Hopf bifurcations (characteristic for systems with interlinked positive and negative feedbacks), are highly sensitive to gene copy number. Hence, these models cannot explain existing experiments. Analysing NF-kB system, the authors show the importance of saturation in transcription of the NF-kB inhibitor IkBa. Models without saturation predict that the oscillation period is a rapidly growing function of total NF-kB level, which is in disagreement with experiments. The study supports the hypothesis that oscillations of robust period are based on supercritical Hopf bifurcation, characteristic for dynamical systems involving negative feedback and time delay. We hypothesise that in the p53 system, the role of positive feedback is not sustaining oscillations, but terminating them in severely damaged cells in which the apoptotic programme should be initiated.

Nuclear factors p53 and NF-kB control many physiological processes including cell cycle arrest, DNA repair, apoptosis, death, innate and adaptive immune responses, and inflammation. There are numerous pathways linking these systems and there is a bulk of evidence for cooperation as well as for antagonisms between p53 and NF-kB. In this theoretical study, the authors use earlier models of p53 and NF-kB systems and construct a crosstalk model of p53–NF-kB network in order to explore the consequences of the two-way coupling, in which NF-kB upregulates the transcription of p53, whereas in turn p53 attenuates transcription of NF-kB inhibitors IkBa and A20. We consider a number of protocols in which cells are stimulated by tumour necrosis factor-a (TNFa) (that activates NF-kB pathway) and/or gamma irradiation (that activates p53 pathway). The authors demonstrate that NF-kB may have both anti- and pro-apoptotic roles. TNFa stimulation, preceding DNA damaging irradiation, makes cells more resistant to irradiation-induced apoptosis, whereas the same TNFa stimulation, when preceded by irradiation, increases the apoptotic cell fraction. The finding suggests that diverse roles of NF-kB in apoptosis and cancer could be related to the dynamical context of activation of p53 and NF-kB pathways.

The stochastic dynamics of T cell receptor (TCR) signaling are studied using a mathematical model intended to capture kinetic proofreading (sensitivity to ligand–receptor binding kinetics) and negative and positive feedback regulation mediated, respectively, by the phosphatase SHP1 and the MAP kinase ERK. The model incorporates protein–protein interactions involved in initiating TCR-mediated cellular responses and reproduces several experimental observations about the behavior of TCR signaling, including robust responses to as few as a handful of ligands (agonist peptide–MHC complexes on an antigen-presenting cell), distinct responses to ligands that bind TCR with different lifetimes, and antagonism. Analysis of the model indicates that TCR signaling dynamics are marked by significant stochastic fluctuations and bistability, which is caused by the competition between the positive and negative feedbacks. Stochastic fluctuations are such that single-cell trajectories differ qualitatively from the trajectory predicted in the deterministic approximation of the dynamics. Because of bistability, the average of single-cell trajectories differs markedly from the deterministic trajectory. Bistability combined with stochastic fluctuations allows for switch-like responses to signals, which may aid T cells in making committed cell-fate decisions.

T cell activation, Mathematical model, Kinetic proofreading, Hysteresis, Ordinary differential equations, Stochastic simulations

The p53 regulatory pathway controls cell responses, which include cell cycle arrest, DNA repair, apoptosis and cellular senescence. We propose a stochastic model of p53 regulation, which is based on two feedback loops: the negative, coupling p53 with its immediate downregulator Mdm2, and the positive, which involves PTEN, PIP3 and Akt. Existence of the negative feedback assures homeostasis of healthy cells and oscillatory responses of DNA-damaged cells, which are persistent when DNA repair is inefficient and the positive feedback loop is broken. The positive feedback destroys the negative coupling between Mdm2 and p53 by sequestering most of Mdm2 in cytoplasm, so it may no longer prime the nuclear p53 for degradation. It works as a clock, giving the cell some time for DNA repair. However, when DNA repair is inefficient, the active p53 rises to a high level and triggers transcription of proapoptotic genes. As a result, small DNA damage may be repaired and the cell may return to its initial ‘‘healthy’’ state, while the extended damage results in apoptosis. The stochasticity of p53 regulation, introduced at the levels of gene expression, DNA damage and repair, leads to high heterogeneity of cell responses and causes cell population split after irradiation into subpopulations of apoptotic and surviving cells, with fraction of apoptotic cells growing with the irradiation dose.

Signaling pathways, Positive and negative feedbacks, Mathematical modeling, Stochastic simulations, Apoptosis, p53, Mdm2, PTEN

Background
The NF-κ B regulatory network controls innate immune response by transducing variety of pathogen-derived and cytokine stimuli into well defined single-cell gene regulatory events.

Results
We analyze the network by means of the model combining a deterministic description for molecular species with large cellular concentrations with two classes of stochastic switches: cell-surface receptor activation by TNFα ligand, and Iκ Bα and A20 genes activation by NF-κ B molecules. Both stochastic switches are associated with amplification pathways capable of translating single molecular events into tens of thousands of synthesized or degraded proteins. Here, we show that at a low TNFα dose only a fraction of cells are activated, but in these activated cells the amplification mechanisms assure that the amplitude of NF-κ B nuclear translocation remains above a threshold. Similarly, the lower nuclear NF-κ B concentration only reduces the probability of gene activation, but does not reduce gene expression of those responding.

Conclusion
These two effects provide a particular stochastic robustness in cell regulation, allowing cells to respond differently to the same stimuli, but causing their individual responses to be unequivocal. Both effects are likely to be crucial in the early immune response: Diversity in cell responses causes that the tissue defense is harder to overcome by relatively simple programs coded in viruses and other pathogens. The more focused single-cell responses help cells to choose their individual fates such as apoptosis or proliferation. The model supports the hypothesis that binding of single TNFα ligands is sufficient to induce massive NF-κ B translocation and activation of NF-κ B dependent genes.

The paper concerns the problem of fitting mathematical models of cell signaling pathways. Such models frequently take the form of sets of nonlinear ordinary differential equations. While the model is continuous in time, the performance index used in the fitting procedure involves measurements taken at discrete time moments. Adjoint sensitivity analysis is a tool which can be used for finding the gradient of a performance index in the space of parameters of the model. In the paper, a structural formulation of adjoint sensitivity analysis called the Generalized Backpropagation Through Time (GBPTT) is used. The method is especially suited for hybrid, continuous-discrete time systems. As an example, we use the mathematical model of the NF-kB regulatory module, which plays a major role in the innate immune response in animals.

Biology and genetics, modeling, ordinary differential equations, parameter learni

The paper is devoted to a stochastic process introduced in the recent paper by Lipniacki et al. [T. Lipniacki, P. Paszek, A. Marciniak-Czochra, A.R. Brasier, M. Kimmel, Transcriptional stochasticity in gene expression, J. Theor. Biol. 238 (2006) 348–367] in modelling gene expression in eukaryotes. Starting from the full generator of the process we show that its distributions satisfy a (Fokker–Planck-type) system of partial differential equations. Then, we construct a c0 Markov semigroup in L1 space corresponding to this system. The main result of the paper is asymptotic stability of the involved semigroup in the set of densities.

Piece-wise deterministic process, Stochastic gene expression, Semigroups of operators, Feller semigroups, Dual semigroups, Markov semigroups, Asymptotic stability

In the article, we discuss the state of art and perspectives in deterministic and stochastic models of NFκB regulatory module. The NFκB is a transcription factor controlling various immune responses including inflammation and apoptosis. It is tightly regulated by at least two negative feedback loops involving IκBα and A20. This mode of regulation results in nucleus-to-cytoplasm oscillations in NFκB localization, which induce subsequent waves of NFκB responsive genes. Single cell experiments carried by several groups provided comprehensive evidence that stochastic effects play an important role in NFκB regulation. From modeling point of view, living cells might be considered noisy or stochastic biochemical reactors. In eukaryotic cells, in which the number of protein or mRNA molecules is relatively large, stochastic effects primarily originate in regulation of gene activity. Transcriptional activity of a gene can be initiated by trans-activator molecules binding to the specific regulatory site(s) in the target gene. The stochastic event of gene activation is amplified by transcription and translation, since it results in a burst of mRNA molecules, and each copy of mRNA then serves as a template for numerous protein molecules. Another potential source of variability can be receptors activation. At low-dose stimulation, important in cell-to-cell signaling, the number of active receptors can be low enough to introduce substantial noise to downstream signaling. Stochastic modeling confirms the large variability in cell responses and shows that no cell behaves like an “average” cell. This high cell-to-cell variability can be one of the weapons of the immune defense. Such non-deterministic defense may be harder to overcome by relatively simple programs coded in viruses and other pathogens.

Mathematical modeling, Systems biology, Immune response, NFκB, Stochastic regulation, Stochastic robustness, Single cell

The higher organisms, eukaryotes, are diploid and most of their genes have two homological copies (alleles). However, the number of alleles in a cell is not constant. In the S phase of the cell cycle all the genome is duplicated and then in the G2 phase and mitosis, which together last for several hours, most of the genes have four copies instead of two. Cancer development is, in many cases, associated with a change in allele number. Several genetic diseases are caused by haploinsufficiency: Lack of one of the alleles or its improper functioning. In the paper we consider the stochastic expression of a gene having a variable number of copies. We applied our previously developed method in which the reaction channels are split into slow (connected with change of gene state) and fast (connected with mRNA/protein synthesis/decay), the later being approximated by deterministic reaction rate equations. As a result we represent gene expression as a piecewise deterministic time-continuous Markov process, which is further related with a system of partial differential hyperbolic equations for probability density functions (pdfs) of protein distribution. The stationary pdfs are calculated analytically for haploidal gene or numerically for diploidal and tetraploidal ones. We distinguished nine classes of simultaneous activation of haploid, diploid and tetraploid genes. This allows for analysis of potential consequences of gene duplication or allele loss. We show that when gene activity is autoregulated by a positive feedback, the change in number of gene alleles may have dramatic consequences for its regulation and may not be compensated by the change of efficiency of mRNA synthesis per allele.

stochastic gene expression, feedback regulation, diploid genes, haploinsufficiency, piecewise deterministic time-continuous Markov process

Living cells may be considered noisy or stochastic biochemical reactors. In eukaryotic cells, in which the number of protein or mRNA molecules is relatively large, the stochastic effects originate primarily in regulation of gene activity. Transcriptional activity of a gene can be initiated by transactivator molecules binding to the specific regulatory site(s) in the target gene. The stochasticity of activator binding and dissociation is amplified by transcription and translation, since target gene activation results in a burst of mRNAs molecules, and each copy of mRNA then serves as a template for numerous protein molecules. In this article, we reformulate our model of the NF-κB regulatory module to analyze a single cell regulation. Ordinary differential equations, used for description of fast reaction channels of processes involving a large number of molecules, are combined with a stochastic switch to account for the activity of the genes involved. The stochasticity in gene transcription causes simulated cells to exhibit large variability. Moreover, none of them behaves like an average cell. Although the average mRNA and protein levels remain constant before tumor necrosis factor (TNF) stimulation, and stabilize after a prolonged TNF stimulation, in any single cell these levels oscillate stochastically in the absence of TNF and keep oscillating under the prolonged TNF stimulation. However, in a short period of ∼90 min, most cells are synchronized by the TNF signal, and exhibit similar kinetics. We hypothesize that this synchronization is crucial for proper activation of early genes controlling inflammation. Our theoretical predictions of single cell kinetics are supported by recent experimental studies of oscillations in NF-κB signaling made on single cells.

Due to the small number of copies of molecular species involved, such as DNA, mRNA and regulatory proteins, gene expression is a stochastic phenomenon. In eukaryotic cells, the stochastic effects primarily originate in regulation of gene activity. Transcription can be initiated by a single transcription factor binding to a specific regulatory site in the target gene. Stochasticity of transcription factor binding and dissociation is then amplified by transcription and translation, since target gene activation results in a burst of mRNA molecules, and each mRNA copy serves as a template for translating numerous protein molecules. In the present paper, we explore a mathematical approach to stochastic modeling. In this approach, the ordinary differential equations with a stochastic component for mRNA and protein levels in a single cells yield a system of first-order partial differential equations (PDEs) for two-dimensional probability density functions (pdf). We consider the following examples: Regulation of a single auto-repressing gene, and regulation of a system of two mutual repressors and of an activator–repressor system. The resulting PDEs are approximated by a system of many ordinary equations, which are then numerically solved.

Gene regulation, Transcription, Stochasticity, Probability density function, Transport-type equations

The paper explores the two scale approach to the incompressible dynamics of superfluid 4He. The fluid is described by a system of equations: Navier–Stokes and Euler equations for the macroscopic normal and superfluid velocity fields respectively. The two equations couple via a mutual friction force exerted on superfluid (quantum) vortices by the normal component. The magnitude of this force, calculated via the analysis of dynamics of quantum vortices in the microscopic scale, is proportional to the value of the counterflow (relative velocity of two helium components) and to the density of quantum vortices. The latter is determined by the generalized Vinen equation, adequate for flows having a net macroscopic vorticity. The generalized equation includes the drift of the anisotropic vortex tangle caused by Magnus force. The derived system of equations is applied first to the analysis of steady state solution of rotating turbulence, and then to the numerical analysis of formation of plane Couette flow between two infinite parallel material surfaces z=0 and z=D. For t<0 both surfaces and the fluid are at rest, then at t=0 one material surface starts moving along the x axis with velocity V0. The viscosity forces cause the motion of normal component and the counterflow which make the line-length density grow, and the two components become coupled by the mutual friction. The fact that superfluid velocity tends to match with the normal velocity makes that ωs≠0, which implies polarization and drift of the tangle. At a given temperature the dynamics depends solely on DV0.

Dynamics of 4He, Quantum turbulence, Anisotropy of quantum tangle, Rotating turbulence

The stochastic nature of gene regulation still remains not fully understood. In eukaryotes, the stochastic effects are primarily attributable to the binary nature of genes, which are considered either switched “on” or “off” due to the action of the transcription factors binding to the promoter. In the time period when the gene is activated, bursts of mRNA transcript are produced. In the present paper, we investigate regulation of gene expression at the single cell level. We propose a mechanism of gene regulation, which is able to explain the observed distinct transcription profiles assuming the number of co-regulatory activities, without attempting to identify the specific proteins involved. The model is motivated by our experiments on NF-κBκB-dependent genes in HeLa cells. Our experimental data shows that NF-κBκB-dependent genes can be stratified into three characteristic groups according to their expression profiles: early, intermediate and late having maximum of expression at about 1, 3 and 6 h, respectively, from the beging of TNF stimulation. We provide a tractable analytical approach, not only in the terms of expected expression profiles and their moments, which corresponds to the measurements on the cell population, but also in the terms of single cell behavior. Comparison between these two modes of description reveals that single cells behave qualitatively different from the cell population. This analysis provides insights useful for understanding of microarray experiments.

Stochastic gene regulation, Expression profiles, Single cell simulations, NF-κκB

p53 is the key transcription factor controlling cellular responses to oncogenic stimulation and DNA da mage. Its activity is tightly controlled by numerous feedback loops. In response to DNA damage, p53 promotes expression of proteins, which suppress cell cycle and activate DNA repair. If the damage is irreparable or the repair takes too long, the programmed cell death (apoptosis) is initiated. In the current study we analyze the apoptotic module, a part of our larger p53 pathway model. In the model, the apoptosis is triggered due to the suppression of Akt activity and/or elevated level of p53 killer.p53 killer, i.e. p53 form phosphorylated at Ser-46 (in addition to Ser-15 and Ser-20), promotes synthesis of pro- apoptotic protein Bax. In healthy cells Bax is inactive due to binding to Bcl-2, another member of Bcl-2 family proteins. Suppression of Akt activity leads to the dissociation of Bax:Bcl-2 complexes and release of Bax. Therefore, two signals may lead to the accumulation of free Bax: one coming from elevated level of p53 killer, the other resulting from decreased level of active Akt. We demonstrated that, depending on parameters, the apoptosis can be controlled by the logic gate ‘AND’ as well as gate ‘OR’. In the first case both signals are required simultaneously, while in the latter case any of the two signals suffices for the initiation of apoptosis.

The aim of this study is to demonstrate that in dynamical systems with underlying bistability the type of noise qualitatively influences the stationary probability distribution (SPD). Specifically, we consider a simplified model of gene expression with the nonlinear positive feedback, which in the deterministic approximation has two stable steady state solutions. Two types of noise are considered; transcriptional - due to the limited number of protein molecules, and gene switching noise - due to gene activation and inactivation. In the limit of zero noise, the SPD generically concentrates in the decreasing vicinity of one of the two stable steady states. We demonstrated that for a range of parameters the SPD corresponding to the system with transcriptional noise only concentrates around a different steady state than SPD corresponding to the system with gene switching noise only.

Gene expression, Bistability, Stochastic processes, Epigenetic attractors

The MAPK cascade is a principal kinase transduction pathway in eukaryotic cells. It transmits signals through three layers of sequentially activated kinases, RAF, MEK and ERK. The latter two kinases require dual phosphorylation for activation. Another property of MAPK signalling is its involvement of scaffolds - multidomain proteins that can assemble protein complexes; in this case the three MAPK components. In this study we analyze analytically and numerically four heuristic models of MAPK signalling in the presence or absence of a scaffold considering both real MEK and ERK kinases and their hypothetical isoforms that require only monophosphorylation. Based on this analysis we will demonstrate that double phosphorylation enforces signaling through scaffolds, which increases both versatility and specificity of the regulation.

We consider a kinase auto-activation model in which the number of activated kinases follows the timecontinuous Markov process. In the deterministic approximation the process is described by the single nonlinear ordinary differential equation, which may have two stable steady states. We found that for sufficiently large number of kinases, the stationary probability distribution given by the Markov process concentrates in the vicinity of the two stable steady states of the deterministic approximation. However, if the number of kinases diverges to the infinity (zero noise limit), the stationary probability distribution concentrates (generically) in only one of the two steady states.