Markov chain aggregation and its application to rule-based modelling

Rule-based modelling allows to represent molecular interactions in a compact and natural way. The underlying molecular dynamics, by the laws of stochastic chemical kinetics, behaves as a continuous-time Markov chain. However, this Markov chain enumerates all possible reaction mixtures, rendering the analysis of the chain computationally demanding and often prohibitive in practice. We here describe how it is possible to efficiently find a smaller, aggregate chain, which preserves certain properties of the original one. Formal methods and lumpability notions are used to define algorithms for automated and efficient construction of such smaller chains (without ever constructing the original ones). We here illustrate the method on an example and we discuss the applicability of the method in the context of modelling large signalling pathways

Syntactic Markovian Bisimulation for Chemical Reaction Networks

In chemical reaction networks (CRNs) with stochastic semantics based on continuous-time Markov chains (CTMCs), the typically large populations of species cause combinatorially large state spaces. This makes the analysis very difficult in practice and represents the major bottleneck for the applicability of minimization techniques based, for instance, on lumpability. In this paper we present syntactic Markovian bisimulation (SMB), a notion of bisimulation developed in the Larsen-Skou style of probabilistic bisimulation, defined over the structure of a CRN rather than over its underlying CTMC. SMB identifies a lumpable partition of the CTMC state space a priori, in the sense that it is an equivalence relation over species implying that two CTMC states are lumpable when they are invariant with respect to the total population of species within the same equivalence class. We develop an efficient partition-refinement algorithm which computes the largest SMB of a CRN in polynomial time in the number of species and reactions. We also provide an algorithm for obtaining a quotient network from an SMB that induces the lumped CTMC directly, thus avoiding the generation of the state space of the original CRN altogether. In practice, we show that SMB allows significant reductions in a number of models from the literature. Finally, we study SMB with respect to the deterministic semantics of CRNs based on ordinary differential equations (ODEs), where each equation gives the time-course evolution of the concentration of a species. SMB implies forward CRN bisimulation, a recently developed behavioral notion of equivalence for the ODE semantics, in an analogous sense: it yields a smaller ODE system that keeps track of the sums of the solutions for equivalent species.Comment: Extended version (with proofs), of the corresponding paper published at KimFest 2017 (http://kimfest.cs.aau.dk/

Exact Hybrid Particle/Population Simulation of Rule-Based Models of Biochemical Systems

Detailed modeling and simulation of biochemical systems is complicated by the problem of combinatorial complexity, an explosion in the number of species and reactions due to myriad protein-protein interactions and post-translational modifications. Rule-based modeling overcomes this problem by representing molecules as structured objects and encoding their interactions as pattern-based rules. This greatly simplifies the process of model specification, avoiding the tedious and error prone task of manually enumerating all species and reactions that can potentially exist in a system. From a simulation perspective, rule-based models can be expanded algorithmically into fully-enumerated reaction networks and simulated using a variety of network-based simulation methods, such as ordinary differential equations or Gillespie's algorithm, provided that the network is not exceedingly large. Alternatively, rule-based models can be simulated directly using particle-based kinetic Monte Carlo methods. This "network-free" approach produces exact stochastic trajectories with a computational cost that is independent of network size. However, memory and run time costs increase with the number of particles, limiting the size of system that can be feasibly simulated. Here, we present a hybrid particle/population simulation method that combines the best attributes of both the network-based and network-free approaches. The method takes as input a rule-based model and a user-specified subset of species to treat as population variables rather than as particles. The model is then transformed by a process of "partial network expansion" into a dynamically equivalent form that can be simulated using a population-adapted network-free simulator. The transformation method has been implemented within the open-source rule-based modeling platform BioNetGen, and resulting hybrid models can be simulated using the particle-based simulator NFsim. Performance tests show that significant memory savings can be achieved using the new approach and a monetary cost analysis provides a practical measure of its utility. © 2014 Hogg et al

Computational and Mathematical Modelling of the EGF Receptor System

This chapter gives an overview of computational and mathematical modelling of the EGF receptor system. It begins with a survey of motivations for producing such models, then describes the main approaches that are taken to carrying out such modelling, viz. differential equations and individual-based modelling. Finally, a number of projects that applying modelling and simulation techniques to various aspects of the EGF receptor system are described

Efficient Syntax-Driven Lumping of Differential Equations

We present an algorithm to compute exact aggregations of a class of systems of ordinary differential equations (ODEs). Our approach consists in an extension of Paige and Tarjan’s seminal solution to the coarsest refinement problem by encoding an ODE system into a suitable discrete-state representation. In particular, we consider a simple extension of the syntax of elementary chemical reaction networks because (i) it can express ODEs with derivatives given by polynomials of degree at most two, which are relevant in many applications in natural sciences and engineering; and (ii) we can build on two recently introduced bisimulations, which yield two complementary notions of ODE lumping. Our algorithm computes the largest bisimulations in O(r⋅s⋅logs)O(r⋅s⋅log⁡s) time, where r is the number of monomials and s is the number of variables in the ODEs. Numerical experiments on real-world models from biochemistry, electrical engineering, and structural mechanics show that our prototype is able to handle ODEs with millions of variables and monomials, providing significant model reductions

Search for the Decays B^0 -> D^{(*)+} D^{(*)-}

Using the CLEO-II data set we have searched for the Cabibbo-suppressed decays B^0 -> D^{(*)+} D^{(*)-}. For the decay B^0 -> D^{*+} D^{*-}, we observe one candidate signal event, with an expected background of 0.022 +/- 0.011 events. This yield corresponds to a branching fraction of Br(B^0 -> D^{*+} D^{*-}) = (5.3^{+7.1}_{-3.7}(stat) +/- 1.0(syst)) x 10^{-4} and an upper limit of Br(B^0 -> D^{*+} D^{*-}) D^{*\pm} D^\mp and B^0 -> D^+ D^-, no significant excess of signal above the expected background level is seen, and we calculate the 90% CL upper limits on the branching fractions to be Br(B^0 -> D^{*\pm} D^\mp) D^+ D^-) < 1.2 x 10^{-3}.Comment: 12 page postscript file also available through http://w4.lns.cornell.edu/public/CLNS, submitted to Physical Review Letter

Evidence for an excess of B -> D(*) Tau Nu decays

Based on the full BaBar data sample, we report improved measurements of the ratios R(D(*)) = B(B -> D(*) Tau Nu)/B(B -> D(*) l Nu), where l is either e or mu. These ratios are sensitive to new physics contributions in the form of a charged Higgs boson. We measure R(D) = 0.440 +- 0.058 +- 0.042 and R(D*) = 0.332 +- 0.024 +- 0.018, which exceed the Standard Model expectations by 2.0 sigma and 2.7 sigma, respectively. Taken together, our results disagree with these expectations at the 3.4 sigma level. This excess cannot be explained by a charged Higgs boson in the type II two-Higgs-doublet model. We also report the observation of the decay B -> D Tau Nu, with a significance of 6.8 sigma.Comment: Expanded section on systematics, text corrections, improved the format of Figure 2 and included the effect of the change of the Tau polarization due to the charged Higg

A search for the decay modes B+/- to h+/- tau l

We present a search for the lepton flavor violating decay modes B+/- to h+/- tau l (h= K,pi; l= e,mu) using the BaBar data sample, which corresponds to 472 million BBbar pairs. The search uses events where one B meson is fully reconstructed in one of several hadronic final states. Using the momenta of the reconstructed B, h, and l candidates, we are able to fully determine the tau four-momentum. The resulting tau candidate mass is our main discriminant against combinatorial background. We see no evidence for B+/- to h+/- tau l decays and set a 90% confidence level upper limit on each branching fraction at the level of a few times 10^-5.Comment: 15 pages, 7 figures, submitted to Phys. Rev.

Study of the reaction e^{+}e^{-} -->J/psi\pi^{+}\pi^{-} via initial-state radiation at BaBar

We study the process $e^+e^-\to J/\psi\pi^{+}\pi^{-}$ with initial-state-radiation events produced at the PEP-II asymmetric-energy collider. The data were recorded with the BaBar detector at center-of-mass energies 10.58 and 10.54 GeV, and correspond to an integrated luminosity of 454 $\mathrm{fb^{-1}}$. We investigate the $J/\psi \pi^{+}\pi^{-}$ mass distribution in the region from 3.5 to 5.5 $\mathrm{GeV/c^{2}}$. Below 3.7 $\mathrm{GeV/c^{2}}$ the $\psi(2S)$ signal dominates, and above 4 $\mathrm{GeV/c^{2}}$ there is a significant peak due to the Y(4260). A fit to the data in the range 3.74 -- 5.50 $\mathrm{GeV/c^{2}}$ yields a mass value $4244 \pm 5$ (stat) $\pm 4$ (syst)$\mathrm{MeV/c^{2}}$ and a width value $114 ^{+16}_{-15}$ (stat)$\pm 7$(syst)$\mathrm{MeV}$ for this state. We do not confirm the report from the Belle collaboration of a broad structure at 4.01 $\mathrm{GeV/c^{2}}$. In addition, we investigate the $\pi^{+}\pi^{-}$ system which results from Y(4260) decay

First Observation of τ→3πηντ\tau\to 3\pi\eta\nu_{\tau} and τ→f1πντ\tau\to f_{1}\pi\nu_{\tau} Decays

We have observed new channels for $\tau$ decays with an $\eta$ in the final state. We study 3-prong tau decays, using the $\eta\to\gamma\gamma$ and \eta\to 3\piz decay modes and 1-prong decays with two \piz's using the $\eta\to\gamma\gamma$ channel. The measured branching fractions are \B(\tau^{-}\to \pi^{-}\pi^{-}\pi^{+}\eta\nu_{\tau}) =(3.4^{+0.6}_{-0.5}\pm0.6)\times10^{-4} and \B(\tau^{-}\to \pi^{-}2\piz\eta\nu_{\tau} =(1.4\pm0.6\pm0.3)\times10^{-4}. We observe clear evidence for $f_1\to\eta\pi\pi$ substructure and measure \B(\tau^{-}\to f_1\pi^{-}\nu_{\tau})=(5.8^{+1.4}_{-1.3}\pm1.8)\times10^{-4}. We have also searched for $\eta'(958)$ production and obtain 90% CL upper limits \B(\tau^{-}\to \pi^{-}\eta'\nu_\tau)<7.4\times10^{-5} and \B(\tau^{-}\to \pi^{-}\piz\eta'\nu_\tau)<8.0\times10^{-5}.Comment: 11 page postscript file, postscript file also available through http://w4.lns.cornell.edu/public/CLN