Seminars (SDB)

Co-authors: Kevin Leder (University of Minnesota), Marc Ryser (Duke University), Walter Lee (Duke University)

The accumulation and spatial spread of mutations during carcinogenesis leads to the production of premalignant 'fields' in many epithelial cancers. I will discuss a spatial stochastic model of this process, based on the biased voter model, and analysis of the extent and geometry of these fields. We apply this model to head and neck squamous cell carcinoma and use it to understand population-level incidence and recurrence patterns.

Co-authors: David F. Anderson (University of Wisconsin-Madison), Yu Sun (University of Wisconsin-Madison)

I will analyze and compare the computational complexity of different simulation strategies for continuous time Markov chains. I consider the task of approximating the expected value of some functional of the state of the system over a compact time interval. This task is a bottleneck in many large-scale computations arising in biochemical kinetics and cell biology. In this context, the terms 'Gillespie's method', 'The Stochastic Simulation Algorithm' and 'The Next Reaction Method' are widely used to describe exact simulation methods. For example, Google Scholar records more than 6,000 citations to Gillespie's seminal 1977 paper. I will look at the use of standard Monte Carlo when samples are produced by exact simulation and by approximation with tau-leaping or an Euler-Maruyama discretization of a diffusion approximation. In particular, I will point out some possible pitfalls when computational complexity is analysed. Appropriate modifications o f recently proposed multilevel Monte Carlo algorithms will then be studied for the tau-leaping and Euler-Maruyama approaches. I will pay particular attention to a parameterization of the problem that, in the mass action chemical kinetics setting, corresponds to the classical system size scaling.

Related Links

• http://personal.strath.ac.uk/d.j.higham/

Within the field of systems biology there is increasing interest in developing computational models which simulate the dynamics of intra-cellular biochemical reaction networks and incorporate the stochasticity inherent in such processes. These models can often be represented as nonlinear multivariate Markov processes. Analysing such models, comparing competing models and fitting model parameters to experimental data are all challenging problems. This talk will provide an overview of a Bayesian approach to the problem. Since the models are typically intractable, use is often made of algorithms exploiting forward simulation from the model in order to render the analysis "likelihood free". There have been a number of recent developments in the literature relevant to this problem, involving a mixture of sequential and Markov chain Monte Carlo methods. Particular emphasis will be placed on the problem of Bayesian parameter inference for the rate constants of stochastic b iochemical network models, using noisy, partial high-resolution time course data, such as that obtained from single-cell fluorescence microscopy studies.

Co-authors: Brian Drawert (UC Santa Barbara), Michael Lawson (Uppsala University), Tau-Mu Yi (UC Santa Barbara), Mustafa Khammash (ETH Zurich), Otger Campas (UC Santa Barbara), Michael Trogdon (UC Santa Barbara) AbstractPolarization is an essential behavior of living cells, yet the dynamics of this symmetry-breaking process are not fully understood. Previously, noise was thought to interfere with this process; however, we show that stochastic dynamics plan an essential role in robust cell polarization and the dynamic response to changing cues. 

To further our understanding of polarization, we have developed a spatial stochastic model of cellular polarization during mating of Saccharomyces cerevisiae. Specifically we investigate the ability of yeast cells to sense a spatial gradient of mating pheromone and respond by forming a projection in the direction of the mating partner. Our mechanistic model integrates three components of the polarization process: the G-protein cycle activated by pheromone bound receptors, the focusing of a Cdc42 polarization cap, and the formation of the tight localization of proteins on the membrane known as the polarisome. 

Our results demonstrate that higher levels of stochastic noise result in increased robustness, giving support to a cellular model where noise and spatial heterogeneity combine to achieve robust biological function. Additionally, our simulations predict that two positive feedback loops are required to generate the spatial amplification to produce focal polarization. We combined our modeling with experiments to explore the critical role of the polarisome scaffold protein Spa2 during yeast mating, and as a result, have characterized a novel positive feedback loop critical to focal polarization via the stabilization of actin cables. 

All intracellular processes involve components present in low numbers, creating spontaneous fluctuations that in turn can enslave the components present in high numbers. The mechanisms often appear complex, with reaction rates that depend nonlinearly on concentrations, indirect feedback loops, and distributed delays. Most systems are also sparsely characterized, with a few steps known in detail but many important interactions not even identified. In the first half of the talk, I will discuss mathematical approaches that exploit natural fluctuations to more reliably analyze data and to make predictions about what complex biological networks cannot do. In the second half I will discuss some of our recent experimental results on the role of fluctuations in cells, e.g. in the segregation of mitochondria, oscillations of synthetic genetic networks, bacterial cell fate decisions, and DNA repair.

Accurate modeling of reaction kinetics is important for understanding the functionality of biological cells. Depending on the particle concentrations and on the relation between particle mobility and reaction rate constants, different mathematical models are appropriate. In the limit of slow diffusion and small concentrations, both discrete particle numbers and spatial inhomogeneity must be taken into account. The most detailed model consists of particle-based reaction-diffusion dynamics, where all individual particles are explicitly resolved in time and space, and particle positions are propagated by diffusion equations, and reaction events may occur only when reactive species are adjacent. For rapid diffusion or large concentrations, the model may be coarse-grained in different ways. Rapid diffusion leads to mixing and implies that spatial resolution is not needed below a certain length scale. This permits the system to be modeled via a spatiotemporal chemical Master equation (STCME), i.e. a coupled set of chemical Master equations acting on spatial sub-volumes. The talk will discuss these different models; in particular, we will see how the STCME description can be derived from particle-based reaction-diffusion dynamics.

Joint work with Stefanie Winkelmann (FU Berlin) 

The last twenty years have seen a revolution in tracking data of biological agents across unprecedented spatial and temporal scales. An important observation from these studies is that path trajectories of living organisms can appear random, but are often poorly described by classical Brownian motion. The analysis of this data can be controversial because practitioners tend to rely on summary statistics that can be produced by multiple, distinct stochastic process models. Furthermore, these summary statistics inappropriately compress the data, destroying details of non-Brownian characteristics that contain vital clues to mechanisms of transport and interaction. In this talk, I will survey the mathematical and statistical challenges that have arisen from recent work on the movement of foreign agents, including viruses, antibodies, and synthetic microparticle probes, in human mucus. 

Abstract Over the last decade or so, there has been an increased focus on the appearance of mixed-mode oscillations in a variety of applications, including cell dynamics and the environment. Taking a broad definition of mixed-mode oscillations, the talk will cover a number of biological phenonema that have mixed-mode oscillations as their basis.These include more familiar dynamics related to canards and coherence resonance, as well as new mechanisms appearing in systems with delays and nonsmooth dynamics. We concentrate on applications in which these mixed mode oscillations are in fact noise-induced, and discuss how the fundamental mechanisms in these examples appear broadly. These ideas point to a number of opportunities in transfering ideas between different areas of life sciences, as well as bringing new ideas to biological modeling from other areas of science and engineering. 

The development of technology to read and write DNA quickly and cheaply is enabling new opportunities for programming biological systems. One example of this is DNA computing, a field devoted to implementing computation in purely biological materials. The hope is that this would enable computation to be performed inside cells, which could pave the way for so-called “smart therapeutics”. Naturally, what we have learned in computer science can be applied to DNA computing systems, and has enabled the implementation of a wide variety of examples of performing computation. Examples include DNA circuits for computing a square root, implementing artificial neural networks, and a general scheme for describing arbitrary chemical reaction networks (CRNs), which itself can be thought of as a compiler.

We have used such a CRN compiler of DNA circuitry to implement the approximate majority (AM) algorithm, which seeks to determine the initial majority of a population of agents holding different beliefs. In its simplest form, the algorithm can be described by three chemical reactions. In this talk, I will describe how we implemented, characterized and modelled a purely DNA implementation of the AM reactions. Along the way, I will demonstrate our software platform for programming biological computation. The platform brings together a variety of stochastic methods that are relevant for both programming and understanding biochemical systems, including stochastic simulation, integration of the chemical master equation, a linear noise approximation, and Markov chain Monte Carlo methods for parameter inference. I will also show preliminary work on synthesizing CRNs with specified probabilistic behaviours.

Related Links

• http://research.microsoft.com/en-us/people/ndalchau/ - Personal website
• http://research.microsoft.com/en-us/projects/dna/ - 'Programming DNA circuits' project

This is a fast 45 minutes tutorial about modeling of stochatic chemical-reaction theory using Markov chains.

Introduction to various escape problems and presentation of the futures lectures.
-Ito calculus,
-stochatic integral and
-The Langevin equation. 
(D. Holcman)

It is well known that stochasticity can play a fundamental role in various biochemical processes, such as cell regulatory networks and enzyme cascades. Isothermal, well-mixed systems can be adequately modelled by Markov processes and, for such systems, methods such as Gillespie's algorithm are typically employed. While such schemes are easy to implement and are exact, the computational cost of simulating such systems can become prohibitive as the frequency of the reaction events increases. This has motivated numerous coarse grained schemes, where the ``fast'' reactions are approximated either using Langevin dynamics or deterministically. While such approaches provide a good approximation for systems where all reactants are present in large concentrations, the approximation breaks down when the fast chemical species exist in small concentrations, giving rise to significant errors in the simulation. This is particularly problematic when using such methods to compute statistics of extinction times for chemical species, as well as computing observables of cell cycle models. In this talk, we present a hybrid scheme for simulating well-mixed stochastic kinetics, using Gillepsie--type dynamics to simulate the network in regions of low reactant concentration, and chemical langevin dynamics when the concentrations of all species is large. These two regimes are coupled via an intermediate region in which a ``blended'' jump-diffusion model is introduced. Examples of gene regulatory networks involving reactions occurring at multiple scales, as well as a cell-cycle model are simulated, using the exact and hybrid scheme, and compared, both in terms weak error, as well as computational cost.  

Collective animal behaviour is often modelled by individual-based models which assume that each individual alters its behaviour according to signals in its neighborhoods. Basic self-propelled particle models can successfully explain some experimentally observed group-level properties, but additional conjectures have to be hypothesized at the individual-level to fully explain experimental data .  In this talk, we will discuss the influence of time delay and noise on the collective  behavior of self-propelled particles system. Firstly, we consider the directional switching of a self-driven particle model with constant, time-varying and random delay times, respectively. The presented analytical and numerical results have demonstrated that time delay can facilitate coherence in self-driven interacting particle systems. Then we discuss the effect of noise on the convergence speed of a stochastic  Cucker-Smale system. We show that noise can accelerate the emergence of flocking.  

Many signalling reactions depend on the signalling protein first localising near its substrate before catalysing a reaction. A common example is the binding of cytosolic enzymes to the unstructured cytoplasmic tails of immune receptors. The tethering (binding) of these enzymes to the tails of immune receptors influences the local concentration of substrate that they experience. In contrast to cytosolic reactions, we currently do not have experimental tools to study tethered signalling reactions. In this work-in-progress talk, I will present an experimental assay that we have been using to study tethered signalling. I will show that the PDE model fails to fit the data (which is averaged over moles of protein) whereas the equivalent stochastic simulation can fit the data. We derived a modified PDE model based on a pair density formalism that can fit the data and the stochastic simulation. Ultimately, we are able to recover 4 parameters that characterise a tethered signalling reaction from the data: binding rates, catalytic rate, and a parameter that determines the properties of the tether.

Asymptotic computation for the escape time through a barrier. 
Singular perturbation of the second order operator.
D. Holcman

Intracellular networks of interacting bio-molecules carry out many essential functions in living cells, but the molecular events underlying the functioning of such networks are ubiquitously random. Stochastic modelling provides an indispensable tool for understanding how cells control, exploit and tolerant the biological noise. A common challenge of stochastic modelling is to calibrate a large number of model parameters against the experimental data. Another difficulty is to study how the behaviours of a stochastic model depends on its parameters, i.e. whether a change in model parameters can lead to a significant qualitative change in model behaviours (bifurcation). One fundamental reason for these challenges is that the existing computational approaches are susceptible to the curse of dimensionality, i.e., the exponential growth in memory and computational requirements in the dimension (number of species and parameters). Herein, we have developed a tensor-based computational framew  ork to address this computational challenge. It is based on recently proposed low-parametric, separable tensor-structured representations of classical matrices and vectors. The framework covers the whole process from solving the underlying equations to automated parametric analysis of the stochastic models such that the high cost of working in high dimensions is avoided. One notable advantage of the proposed approach lies in its ability to capture all probabilistic information of stochastic models all over the parameter space into one single tensor-formatted solution, in a way that allows linear scaling of basic operations with respect to the number of dimensions. Within such framework, the existing algorithms commonly used in the deterministic framework can be directly used in stochastic models, including parameter inference, robustness analysis, sensitivity analysis, and stochastic bifurcation analysis.

Using Green’s function, we  study the additive properties of the MFPT. We use path-integral approach to study Fokker-Planck with a killing term. We introduce a non Poissonian stochastic escape from a limit cycle (part I).D. Holcman

Co-authors: Jose Juan Tapia (University of Pittsburgh), John Sekar (University of Pittsburgh)

In this talk I will address some of the challenges faced in developing detailed models of biochemical networks, which encompass large numbers of interacting components. Although simpler coarse-grained models are often useful for gaining insight into biological mechanisms, such detailed models are necessary to understand how molecular components work in the network context and essential to developing the ability to manipulate such networks for practical benefits. The rule-based modeling (RBM) approach, in which biological molecules can be represented as structured objects whose interactions are governed by rules that describe their biochemical interactions, is the basis for addressing multiple scaling issues that arise in the development of large scale models. Currently available software tools for RBM, such as BioNetGen, Kappa, and Simmune, enable the specification and simulation of large scale models, and these tools are in widespread use by the modeling community. I will re view some of the developments that gave rise to those capabilities, and then I will describe our current efforts broaden the appeal of these tools as well as to better enable collaborative development of models through re-use of existing models and improving visual representations of models.

Related Links

• http://bionetgen.org - For more information about rule-based modeling tools described in this talk.

Co-authors: CHRISTIAN BAYER (WIAS, BERLIN), CHIHEB BEN HAMMOUDA (KAUST, THUWAL), ALVARO MORAES (ARAMCO, DAMMAM), FABRIZIO RUGGERI (IMATI, MILAN), PEDRO VILANOVA (KAUST, THUWAL) 

Stochastic Reaction Networks (SRNs), that are intended to describe the time evolution of interacting particle systems where one particle interacts with the others through a finite set of reaction channels. SRNs have been mainly developed to model biochemical reactions but they also have applications in neural networks, virus kinetics, and dynamics of social networks, among others. 

This talk is focused on novel fast simulation algorithms and statistical inference methods for SRNs. 

Regarding simulation, our novel Multi-level Monte Carlo (MLMC) hybrid methods provide accurate estimates of expected values of a given observable at a prescribed final time. They control the global approximation error up to a user-selected accuracy and up to a certain confidence level, with near optimal computational work. 

With respect to statistical inference, we first present a multi-scale approach, where we introduce a deterministic systematic way of using up-scaled likelihoods for parameter estimation. In a second approach, we derive a new forward-reverse representation for simulating stochastic bridges between consecutive observations. This allows us to use the well-known EM Algorithm to infer the reaction rates. 

Co-authors: Juraj Szavits-Nossan, Kym Eden, Ryan Morris, Martin Evans and Cait MacPhee

In small volumes, the kinetics of filamentous protein self-assembly is expected to show significant variability, arising from intrinsic molecular noise. We introduce a simple stochastic model including nucleation and autocatalytic growth via elongation and fragmentation, which allows us to predict the effects of molecular noise on the kinetics of autocatalytic self-assembly. We derive an analytic expression for the lag-time distribution, which agrees well with experimental results for the fibrillation of bovine insulin. Our analysis shows that significant lag-time variability can arise from both primary nucleation and from autocatalytic growth and should provide a way to extract mechanistic information on early-stage aggregation from small-volume experiments.

Co-author: David Soloveichik (University of Texas, Austin)

Suppose a chemical system requires a single molecule of a certain species $L$. Preparing a solution with just a single copy of $L$ is a difficult task to achieve with imprecise pipettors. Could we engineer artificial reactions (a chemical election algorithm, so to speak) that whittle down an initially large count of $L$ to 1? Yes, with the reaction $L+L \to L+F$: whenever two candidate leaders encounter each other, one drops out of the race. In volume $v$ convergence to a single $L$ requires expected time proportional to $v$; the final reaction --- two lone $L$'s seeking each other in the vast expanse of volume $v$ --- dominates the whole expected time.

One might hope that more cleverly designed reactions could elect a leader more quickly. We dash this hope: $L+L \to L+F$, despite its sloth, is the fastest chemical algorithm for leader election there is (subject to some reasonable constraints on the reactions). The techniques generalize to establish lower bounds on the time required to do other computational tasks, such as computing which of two species $X$ or $Y$ holds an initial majority.

Democracy works... but it's painstakingly slow.

Related Links

• http://web.cs.ucdavis.edu/~doty/papers/slepprlt.pdf - Paper

In intracellular signaling pathways, biochemical activation events are transmitted from one node within the signaling network to another.  Recent work examining the information capacity of signaling pathways has concluded that most signaling pathways have limited abilities to resolve different strengths of inputs.  However, these studies are based on data in which only a single signal is measured in each cell, in response to a given cell, with the limitation that transmission of a signal from one signaling node to another cannot be directly observed.  Other published data suggest that single cells may have a much higher capacity to transmit quantitative information, which is obscured by population heterogeneity.  To better understand the properties of information transmission through biochemical cascades in individual cells, we have developed a panel of live-cell reporters to monitor multiple signaling events in the cell proliferation and growth network (CPGN).  These reporters include activity biosensors for the kinases ERK, Akt, mTOR, and AMPK, and CRISPR-based reporters for ERK target gene expression.  Experimental analysis with these tools reveals the temporal and quantitative linkage properties between nodes of the CPGN.  I will discuss two studies currently underway in our lab.  The first examines the how the CPGN manages the interplay between ATP-producing and ATP-consuming processes during cell proliferation; we find that loss of Akt signaling results in unstable levels of ATP and NADH in proliferating cells.  The second project focuses on how variations in amplitude and duration of ERK activity control the expression of the target gene Fra-1, which is involved in metastasis; here, we show that cancer therapeutics directed at inhibiting this pathway create strikingly different kinetics of ERK activity at the single-cell level, with distinct effects on Fra-1 expression.

Co-authors: Anthony Bowman (Imperial College London), Xi-Ming Sun (MRC CSC), Samuel Marguerat (MRC CSC)

Gene expression is affected by both random timing of reactions (intrinsic noise) and interaction with global stochastic systems in the cells (extrinsic noise). A challenge in inferring parameters of gene expression using models of stochastic gene expression is that these models usually only inlcude intrinsic noise. However, experimental distributions of transcripts are strongly influenced by extrinsic effects including cell cycle and cell division. Here, we present a multi-scale approach in stochastic gene expression to deal with this problem. We apply our methodology to data obtained using single molecule Fish technique in fission yeast. The data suggests cell size influences transcription parameters. We use Approximate Bayesian Computation (ABC) along with sequential Monte Carlo to infer the dependence of gene expression parameters on cell size. Our analysis reveals a linear increase of transcription burst size during the cell cycle.

In order to analyse stochastic chemical systems, we solve the corresponding Fokker-Planck equation numerically. The dimension of this problem corresponds to the number of chemical species and the standard numerical methods fail for systems with already four or more chemical species due to the so called curse of dimensionality. Using tensor methods we succeeded to solve realistic problems in up to seven dimensions and an academic example of a reaction chain of 20 chemical species.

In the talk we will present the Fokker-Planck equation and discuss its well-posedness. We will describe its discretization based on the finite difference method and we will explain the curse of dimensionality. Then we provide the main idea of tensor methods. We will identify several types of errors of the presented numerical scheme, namely the modelling error, the domain truncation error, discretization error, tensor truncation error, and the algebraic error. We will present an idea that equilibration of these errors based on a posteriori error estimates yields considerable savings of the computational time.

Co-authors: Radek Erban (University of Oxford), Kostantinos Zygalakis (University of Edinburgh)

It is well known that stochasticity can play a fundamental role in various biochemical processes, such as cell regulatory networks and enzyme cascades. Isothermal, well-mixed systems can be adequately modelled by Markov processes and, for such systems, methods such as Gillespie's algorithm are typically employed. While such schemes are easy to implement and are exact, the computational cost of simulating such systems can become prohibitive as the frequency of the reaction events increases. This has motivated numerous coarse grained schemes, where the "fast" reactions are approximated either using Langevin dynamics or deterministically. While such approaches provide a good approximation for systems where all reactants are present in large concentrations, the approximation breaks down when the fast chemical species exist in small concentrations, giving rise to significant errors in the simulation. This is particularly problematic when using such methods to compute statistics of extinction times for chemical species, as well as computing observables of cell cycle models. In this talk, we present a hybrid scheme for simulating well-mixed stochastic kinetics, using Gillepsie-type dynamics to simulate the network in regions of low reactant concentration, and chemical Langevin dynamics when the concentrations of all species is large. These two regimes are coupled via an intermediate region in which a "blended"' jump-diffusion model is introduced. Examples of gene regulatory networks involving reactions occurring at multiple scales, as well as a cell-cycle model are simulated, using the exact and hybrid scheme, and compared, both in terms weak error, as well as computational cost.

Biological systems can emerge complexity from simple yet multitude of interactions. Capturing such biological phenomenon mathematically for predictions and inference is being actively researched. Computing systems where the interacting components are inherently stochastic demands large amounts of computational power. Recently, splitting the dynamics of the system into deterministic and stochastic components has been a new strategy for computing biological networks. This hybrid strategy drastically reduces the number of equations to solve, however, the new equations are naturally stiff and nonlinear. Hybrid models are a strong candidate as a numerical method for probing large biological networks with intrinsic stochasticity. In this talk we will take on a new mathematical and numerical perspective of hybrid models. Through many biological examples, we will aim to gain insight into the benefits and stumbling blocks of the hybrid framework.

Related Links

• https://github.com/vikramsunkara/PyME - PyME: Python base CME solver.

Co-author: Mustafa Khammash (ETH Zurich)

The internal dynamics of a living cell is generally very noisy. An important source of this noise is the intermittency of reactions among various molecular species in the cell. The role of this noise is commonly studied using stochastic models for reaction networks, where the dynamics is described by a continuous-time Markov chain whose states represent the molecular counts of various species. In this talk we will discuss how the long-term behavior of such Markov chains can be assessed using a blend of ideas from probability theory, linear algebra and optimisation theory. In particular we will describe how many biomolecular networks can be viewed as generalised birth-death networks, which leads to a simple computational framework for determining their stability properties such as ergodicity and convergence of moments. We demonstrate the wide-applicability of our framework using many examples from Systems and Synthetic Biology. We also discuss how our results can hel p in analysing regulatory circuits within cells and in understanding the entrainment properties of noisy biomolecular oscillators.

Stochasticity is a defining feature of biochemical reaction networks, with molecular fluctuations influencing cell physiology. In principle, master probability equations completely govern the dynamic and steady state behavior of stochastic reaction networks. In practice, a solution had been elusive for decades, when there are second or higher order reactions. A large community of scientists has then reverted to merely sampling the probability distribution of biological networks with stochastic simulation algorithms. Consequently, master equations, for all their promise, have not inspired biological discovery.

We recently presented a closure scheme that solves chemical master equations of nonlinear reaction networks [1]. The zero-information closure (ZI-closure) scheme is founded on the observation that although higher order probability moments are not numerically negligible, they contain little information to reconstruct the master probability [2]. Higher order moments are then related to lower order ones by maximizing the entropy of the network. Using several examples, we show that moment-closure techniques may afford the quick and accurate calculation of steady-state distributions of arbitrary reaction networks.

With the ZI-closure scheme, the stability of the systems around steady states can be quantitatively assessed computing eigenvalues of the moment Jacobian [3]. This is analogous to Lyapunov’s stability analysis of deterministic dynamics and it paves the way for a stability theory and the design of controllers of stochastic reacting systems [4, 5].

In this seminar, we will present the ZI-closure scheme, the calculation of steady state probability distributions, and discuss the stability of stochastic systems.

1. Smadbeck P, Kaznessis YN. A closure scheme for chemical master equations. Proc Natl Acad Sci U S A. 2013 Aug 27;110(35):14261-5.

2. Smadbeck P, Kaznessis YN. Efficient moment matrix generation for arbitrary chemical networks, Chem Eng Sci, 20

Co-author: Arnab Bandyopadhyay (University of Kansas)

Cytoplasmic regulatory networks often approximate well-mixed reaction kinetics in single cells, but with variability from cell to cell. As a result, inheritance dynamics and kin correlations have been implicated in effects on cell cycle, regulatory networks, and modulation of population growth rate. Based on an experimental result in our lab suggesting lineage correlations in bacterial growth arrest, we developed a cellular stochastic simulation framework to analyse the role of lineage in bacterial cells regulating growth rate by means of an intracellular molecular network. The simulation framework thus models both intrinsic and inherited noise sources while maintaining lineage data between cell agents assigned individual unique identifiers.

Our initial application of the framework demonstrates the role of lineage in the probability of bacterial growth arrest controlled by an endogenous toxin from a toxin-antitoxin system. These systems have tight binding between toxin and antitoxin, so that there is a discrete critical threshold in the toxin:antitoxin ratio below which a cell is essentially toxin-free and growth is unrestricted, and above which toxin rapidly slows the growth rate. The subset of high-toxin cells crossing into the growth arrested state are associated with antibiotic persistence. Our implementation of a simple toxin-antitoxin system in the simulation framework revealed the statistical dependence of growth arrest on cellular lineage: after several generations of growth, the probability of cellular growth arrest began to depend on lineage distance. Clusters of closely related cell agents had a high probability of transitioning into growth arrest, while the rest of the lineage continued to grow withou t restriction.

We consider various quantities of interest in multiscale lineage simulations, and conclude that growth transitions in a cellular colony cannot be fully understood without quantitative knowledge of its lineage.

Antibiotic-resistant bacterial infections are significant causes of morbidity and mortality worldwide. These infections kill over 700,000 people globally every year. Because of the paucity of new antibiotics, this number may increase to over 10 million by 2050, surpassing deaths from cancer and diabetes combined (http://amr-review.org/). Using synthetic biology techniques, we engineer probiotic bacteria to produce antimicrobial peptides (AMPs) and deliver the potent antibiotic proteins in the gastrointestinal tract of hosts. We test this new technology and demonstrate that it safely reduces multi-drug resistant pathogenic bacteria in the gut of animals [1,2]. At the heart of biological engineering efforts are multiscale models that guide explanations and predictions of the antagonistic activity of recombinant LAB against pathogenic strains [3]. Models are developed to quantify how AMPs kill bacteria at distinct but tied scales. Using atomistic simulations the various interaction steps between peptides and cell membranes are explored. Mesoscopic models are developed to study ion transport and depolarization of membranes treated with AMPs [4]. Stochastic kinetic models are developed to quantify the strength of synthetic promoters and AMPs expression [5]. In this presentation, we will discuss how modeling facilitates biological engineering and stress important theoretical and numerical challenges.
References

1. Volzing K, Borrero J, Sadowsky MJ, Kaznessis YN. “Antimicrobial Peptides Targeting Gram-negative Pathogens, Produced and Delivered by Lactic Acid Bacteria.” ACS Synth Biol. 2013. 15;2(11):643-50. doi: 10.1021/sb4000367 .
2. Geldart K., Borrero J., Kaznessis Y.N., “A Chloride-Inducible Expression Vector for Delivery of Antimicrobial Peptides Against Antibiotic-Resistant Enterococcus faecium”, Applied and Environmental Microbiology, 2015 81:11 3889-3897.3.
3. Kaznessis Y, “Multiscale Models of Antibiotic Probiotics”, Curr Opin Chem Eng, 2014, 1;6:18-24
4. Bolintineanu D, Kaznessis YN. "Computational studies of protegrin antimicrobial peptides: a review." Peptides. 2011 Jan;32(1):188-201
5. Kaznessis YN. Computational methods in synthetic biology. Biotechnol J. 2009 Oct;4(10):1392-405. doi: 10.1002/biot.200900163. 

This lecture introduces novel concepts in asymptotic of second order PDE. It is a continuation of lecture 7. The motivation is coming from noisy dynamical systems, modelling neuronal networks with synaptic properties. The lecture presents a novel matched asymptotic, based on Mobius conformal mapping, Hopf normal form, to estimate the distribution of exit times and exit points (concentrated at one boundary point). The spectrum of the Fokker-Planck non-selfadjoint operator is computed. The role of the second eigenvalue is explained and generated the oscillation escape.

This lecture continues lecture 8 (asymptotics of second order PDEs) and escape through a limit cycle.  The second part introduces the Stochastic Narrow Escape theory. D. Holcman.

I will discuss an emerging framework where the extensive and powerful toolbox of deterministic dynamical systems can be used to study an important class of noise induced stochastic phenomena, called metastable dynamics. Consider two trajectories: one deterministic and one perturbed by weak noise, each having the same initial conditions. A single fluctuation is very unlikely to perturb the stochastic trajectory very far from the deterministic trajectory, and on short time scales the two remain close. On long time scales, both trajectories approach a stable steady state.  Given enough time, it is possible for a rare sequence of fluctuations to perturb the stochastic trajectory far enough that it moves toward a different steady state. Hence, a stable steady state for a deterministic system becomes a metastable state under the influence of weak noise. I will discuss two applications in biology.   1. The first application is epigenetic switching in gene circuits, the biological problem that first motivated my research. It has been known for decades that gene expression is strongly influenced by random forces within the cell. Rare events driven by noise can cause a dramatic shift in the way a gene is expressed, which can radically alter the state of a cell. One example is an altered state that imparts antibiotic resistance to e-coli bacteria. I will show that diffusion approximations, widely used to study gene expression systems, are inaccurate and unreliable for metastable phenomena. As an alternative, I will propose a far more accurate direct method that eliminates the need for a diffusion approximation.   2. The second application is spontaneous neural activity. Intrinsic noise from molecular fluctuations of voltage-gated ion channels cause spontaneous activity that propagates into and affects local network function.  A spontaneous action potential is a physical example of a new type of first-exit-time problem: the random time to initiate an excitable event in an excitable system with a single fixed point. Using a metastable phase plane analysis, I will show how noise induced excitable events in the stochastic Morris-Leccar model are initiated through a predictable sequence of events. In other words, a single mechanism explains how spontaneous activity is generated. Moreover, the generating mechanism contradicts the current understanding of this phenomena.  It is widely believed that spontaneous activity in most neurons is driven primarily by fast sodium channels, because these channels govern the fast initiation stage of an action potential. Potassium channels respond much more slowly and are responsible for reseting the membrane voltage at the final stage of the action potential. Contrary to the standard paradigm, metastable dynamics predicts that the primary driving force behind spontaneous initiation is the slow potassium channel noise.

This lecture introduces asymptotic computations to estimate the escape time through a cusp. The method uses conformal mapping in novel banana shaped domain.  The second part of the lecture concerns the computation of  the effective diffusion coefficient of a crowded membrane. D. Holcman.

I will present results from Molecular Dynamics simulations of pure andcrowded lipid bilayer systems, giving evidence of anomalous diffusion.While in the pure lipid bilayer this anomaly is very short ranged, theaddition of colesterols or the passage to the gel phase leads to exten-ded anomalous diffusion.  The character of the dynamics corresponds tothat of the fractional Brownian motion in the dilute bilayer, changingto non-Gaussian motion when the bilayer is crowded with proteins. Realbiological membranes show macroscopic anomalous diffusion, which is ev-idenced from superresolution microscopy experiments. In particular themotion becomes non-ergodic and ageing.

This lecture presents models and coarse-graining analysis of stochastic chemical reactions. The first time that k Brownian particles are bound is computed using Markov chains (Mean time to Threshold) . The conditional mean time and splitting probability are computed using certain paths in two dimensional Markov chains. D. Holcman.

The field of anomalous diffusion can be very polarising, with some saying that there is little basis for such studies in biology to those saying most of spatial crowding effects lead to these ideas.  This presentation  will be mainly extempore and off the top of my head but I will try to put some of these ideas in context, warts and all.  I will then branch off into such esoteric fields as fractional differential equations, fractional stochastic differential equations, Poisson tau leap with memory, fractional Poisson processes, Mittag Leffler waiting times and finish with a new Poisson tau leap method for anomalous diffusive effects. This will be very much seat of the pants and I hope there is lots of discussion.

This lecture presents the final computation of the mean time to threshold. The second part concerns modeling telomere length during cell division (computation of the shortest telomere and comparison with the second shortest). The last introduce Rouse and the beta polymer model.

This lecture presents recent methods to study the anomalous exponent of a monomer from a polymer (Rouse and beta). The second part is about asymptotic estimations of the mean first looping time using high dimensional perturbation method of the Laplace operator when removing a tubular neighbourhood of a low dimensional submanifold.

Gene transcription is a highly stochastic, dynamic process. As a result, the mRNA copy number of a given gene is heterogeneous both between cells and across time. I will present a framework to model gene transcription in populations of cells with time-varying (stochastic or deterministic) transcription and degradation rates. Such rates can be understood as upstream cellular drives representing the effect of different aspects of the cellular environment, e.g. external signalling, circadian rhythms, or chromatin remodelling. I will show that the full solution of the generic master equation for gene transcription contains two components: a model-specific, upstream effective drive, which encapsulates the effect of the cellular drives (e.g., entrainment, periodicity or promoter randomness), and a downstream transcriptional Poissonian part, which is common to all models. This analytical framework allows us to treat cell-to-cell and dynamic variability consistently, unifying several approaches in the literature.   Our general solution confers to us two broad advantages. The first is pragmatic: the theory provides us with new approaches for solving non-stationary gene transcription models, analysing single-cell snapshot and time-course data, and reducing the computational cost of sampling solutions via stochastic simulation. The second advantage is conceptual: studying the solution of a broad class of models in generality provides us with physical intuition for the sources of noise and their characteristics, and the ability to deduce which models are analytically solvable (along with the form and structure of their solutions). I will demonstrate some such benefits by applying our solution to several biologically-relevant examples.

We present path-space, information theory-based, sensitivity analysis, uncertainty quantification and variational inference methods for complex high-dimensional stochastic dynamics, including chemical reaction networks with hundreds of parameters, Langevin-type equations and lattice kinetic Monte Carlo. We establish their connections with goal-oriented methods in terms of new, sharp, uncertainty quantification inequalities that scale appropriately at both long times and for high dimensional state and parameter space.   The combination of proposed methodologies is capable to (a) tackle non-equilibrium processes, typically associated with coupled physicochemical mechanisms or boundary conditions, such as reaction-diffusion problems, and where even steady states are unknown altogether, e.g. do not have a Gibbs structure. The path-wise information theory tools,  (b) yield a surprisingly simple, tractable and easy-to-implement approach to quantify and rank parameter sensitivities, as well as  (c) provide reliable parameterizations for coarse-grained molecular systems based on fine-scale data, and rational model selection through path-space (dynamics-based) variational inference methods.

In the bond percolation model on a lattice, we colour vertices with n_c colours independently at random according to Bernoulli distributions. A vertex can receive multiple colours and each of these colours is individually observable. The colours colour the entire component into which they fall. Our goal is to estimate the n_c +1
parameters of the model: the probabilities of colouring of single vertices and the probability with which an edge is open. The input data is the configuration of colours once the complete components have been coloured, without the information which vertices were originally coloured or which edges are open.

We use a Monte Carlo method, the method of simulated moments to achieve this goal. We prove that this method is a strongly consistent estimator by proving a strong law of large numbers for the vertices' weakly dependent colour values. We evaluate the method in computer tests. The motivating application is cross-contamination rate estimation for digital PCR in lab-on-a-chip microfluidic devices.

Co-authors: Libin Abraham (University of British Columbia), Alejandra Herrera (University of British Columbia), Michael R Gold (University of British Columbia), Keng Chou (University of British Columbia), Reza Tafteh (University of British Columbia), Joshua Scurll (University of British Columbia), Henry Lu (University of British Columbia), Ki Woong Sung (University of British Columbia)

I will describe recent progress in interpretation and modelling of B cell signalling based on super-resolution light microscopy (primarily, STORM imaging and single particle tracking).

Co-author: Jonathan Chapman (University of Oxford)

We discuss nonlinear Fokker-Planck models describing diffusion processes with particle interactions. These models are motivated by the study of many particle systems in biology, and arise as the population-level description of a stochastic particle-based model. In particular, we consider a system of impenetrable diffusing spheres and use the method of matched asymptotic expansions to obtain a systematic model reduction. In the second part of the talk, we discuss how this method can be used to derive an effective transport equation in heterogeneous domains, such as porous media or crowded environments. A nice feature of this approach is that it can easily account for macroscopic gradients in porosity or crowding.

I will discuss methods for spatio-temporal modelling in molecular and cell biology, including all-atom and coarse-grained molecular dynamics (MD) and stochastic reaction-diffusion models, with the aim of developing and analysing multiscale methods which use MD simulations in parts of the computational domain and (less-detailed) stochastic reaction-diffusion approaches in the remainder of the domain. The main goal of this multiscale methodology is to use a detailed modelling approach in localized regions of particular interest (in which accuracy and microscopic details are important) and a less detailed model in other regions in which accuracy may be traded for simulation efficiency. Applications using all-atom MD include intracellular dynamics of ions and ion channels. Applications using coarse-grained MD include protein binding to receptors on the cellular membrane, where modern stochastic reaction-diffusion simulators of intracellular processes can be used in the bulk and a ccurately coupled with a (more detailed) MD model of protein binding which is used close to the membrane.

References:
[1] Radek Erban, "Coupling all-atom molecular dynamics simulations of ions in water with Brownian dynamics", Proceedings of the Royal Society A, Volume 472, Number 2186, 20150556 (2016)
[2] Radek Erban, "From molecular dynamics to Brownian dynamics", Proceedings of the Royal Society A, Volume 470, Number 2167, 20140036 (2014)

Related Links

• http://rspa.royalsocietypublishing.org/content/472/2186/20150556 - Reference [1]
• http://rspa.royalsocietypublishing.org/content/470/2167/20140036 - Reference [2]
• http://people.maths.ox.ac.uk/erban/ - Preprints 

I will describe a particle-based algorithm for reaction-diffusion problems that combines Brownian dynamics with a Markov reaction process. The microscopic model simulated by our Split Reactive Brownian Dynamics (SRBD) algorithm is based on the Doi or volume-reactivity model. This model applies only to reactions with at most two reactants, which is physically realistic. Let us consider the simple reaction A+B->product. In the Doi model, particles are independent spherical Brownian walkers (this can be relaxed to account for hydrodynamic interactions), and while an A and a B particle overlap, there is a Poisson process with a given microscopic reaction rate for the two particles to react and give a product. Our goal is to simulate this complex Markov process in dense systems of many particles, in the presence of multiple reaction channels.

Our algorithm is inspired by the Isotropic Direct Simulation Monte Carlo (I-DSMC) method and the next subvolume method. Strang splitting is used to separate diffusion and reaction; this is the only approximation made in our method. In order to process reactions without approximations, with the particles frozen in place, we use an event-driven algorithm. We divide the system into a grid of cells such that only particles in neighboring cells can react. Each cell schedules the next potential reaction to happen involving a particle in that cell and a particle in one of the neighboring cells, and an event queue is used to select the next cell in which a reaction may happen.

Note that, while a grid of cells is used to make the algorithm efficient, the results obtained by the SRBD method are grid-independent and thus free of grid artifacts, such as loss of Galilean invariance and sensitivity of the results to the grid spacing. I will compare our SRBD method with grid-based methods, such as (C)RDME and a variant of RDME that we call Split Brownian Dynamics with Reaction Master Equation (S-BD-RME), on a problem involving the spontaneous f

We have introduced interacting-particle reaction-diffusion (iPRD) simulations as a hybrid between molecular dynamics simulations and particle-based reaction-diffusion simulations. iPRD is a suitable modeling framework for many cellular signaling processes, especially such processes involving dense protein mixtures, supramolecular architectures or protein scaffolds at membranes. I will introduce to the theory behind iPRD and present computational algorithms.

One of our key application areas is clathrin-mediated endocytosis in neurons, and I will briefly elude to simulations of protein recruitment to clathrin-coated pits and the interplay between membrane-associated proteins and membrane deformation in endocytosis.

Related Links

• http://readdy-project.org

Microtubule (MT) acetylation is done post-translationally in the MT lumen by the acetyltransferase alpha-TAT1. A simple picture is that alpha-TAT1 enters the lumen at the open ends of MT, and diffuses inside while acetylating. Complicating this picture, alpha-TAT1 is approximately half the luminal diameter --- indicating that it may undergo single-file diffusion (SFD). Computationally, we explore the consequences of SFD in this system. Experimentally, in collaboration with the group of Guillaume Montagnac (Institut Gustave Roussy), immunofluorescence techniques allows us to measure the acetylation pattern of individual MT. We try to reconcile the computational and experimental approaches, and highlight open questions. 

In mathematical modeling one typically progresses from observations of the world (and some serious thinking!) to equations for a model, and then to the analysis of the model to make predictions. Good mathematical models give good predictions (and inaccurate ones do not) > - but the computational tools for analyzing them are the same: algorithms that are typically based on closed form equations. While the skeleton of the process remains the same, today we witness the development of mathematical techniques that operate directly on observations -data-, and "circumvent" the serious thinking that goes into selecting variables and parameters and writing equations. The process then may appear to the user a little like making predictions by  "looking into a crystal ball". Yet the "serious thinking" is still there and uses the same -and some new- mathematics: it goes into building algorithms that "jump directly" from data to the analysis of the model (which is never available in closed form) so as to make predictions. I will present a couple of efforts that illustrate this new path from data to predictions. It really is the same old path, but it is travelled by new means.

Co-authors: Michael J. Hallock (University of Illinois at Urbana-Champaign), Joseph R. Peterson (University of Illinois at Urbana-Champaign), John A. Cole (University of Illinois at Urbana-Champaign), Tyler M. Earnest (University of Illinois at Urbana-Champaign), John E. Stone (University of Illinois at Urbana-Champaign)

High-performance computing now allows integration of data from cryoelectron tomography, super resolution imaging, various –omics, and systems biology reaction studies into coherent computational models of cells and cellular processes functioning under in vivo conditions. Here we analyze the stochastic reaction-diffusion dynamics of ribosome biogenesis in slow growing bacterial cells undergoing DNA replication and probe the metabolic reprogramming that occurs within dense colonies of Escherichia coli cells over periods of hours. Using our GPU-based Lattice Microbe software, the some 1300 reactions and 250 species involved in transcription, translation and ribosome assembly are described in terms of reaction-diffusion master equations and simulated over a cell cycle of two hours. The ribosome biogenesis simulations account for DNA replication that takes place within the cell cycle, and the results are compared to super resolution imaging results. In the case of the c ell colony simulations, reaction-diffusion kinetics of the surrounding medium are coupled with the cellular metabolic networks to demonstrate how small colonies of interacting bacterial cells differentially respond to the competition for resources according to their position in the colony. The predicted metabolic reprogramming has been observed experimentally. Finally we will report on the progress we have achieved to date and how supercomputers will provide us a window into cellular dynamics within bacterial and eukaryotic cells.

Co-authors: Hepburn, Iain (OIST), Chen, Weiliang (OIST)

Spatial stochastic molecular simulations in biology are limited by the intense computation required to track molecules in space either by particle tracking or voxel-based methods, meaning that the serial limit has already been reached in sub-cellular models. This calls for parallel simulations that can take advantage of the power of modern supercomputers. GPU parallel implementations have been described for particle tracking methods [1,2] and for voxel-based methods [3], where good parallel performance gain up to 2 order of magnitude have been demonstrated but this depends strongly on model specificity. MPI parallel implementations have gained less attention than GPU implementations to date but offer several advantages including a greater range of platform support from personal computers to advanced supercomputer clusters. An initial MPI implementation for irregular grids has been described and almost ideal speedup demonstrated but only up to 4 cores [4], which indicates the potential for good scalability of such implementations.

We describe an operator splitting implementation for irregular grids with a novel method to improve accuracy over Lie-Trotter splitting that is somewhat comparable to tau-reduction but without the performance cost. We systematically investigate parallel performance for a range of models and mesh partitionings using the STEPS simulation platform [5]. Finally we introduce a whole cell parallel simulation of a published reaction-diffusion model [6] within a detailed, complete neuron morphology and demonstrate a speedup of 3 orders of magnitude over serial computations. 

[1] L Dematte 2012. IEEE/ACM Trans. Comput. Biol. Bioinf. 9: 655-667 [2] DV Gladkov et al. 2011. Proc. 19th High Perf. Comp. Symp. 151-158 [3] E Roberts, JE Stone, Z Luthey-Schulten 2013. J. Comp. Chem. 34: 245–255 [4] A Hellander et al. 2014. J. Comput. Phys. 266: 89-100 [5] I Hepburn et al. 2012. BMC Syst. Biol. 6:36 [6] H Anwar et al. 2013. J. Neurosci. 33: 15848-15867

Related Links

• http://steps.sourceforge.net/ - Software site 

Co-authors: Jean-Francois Rupprecht (National University of Singapore, Singapore), Olivier Bénichou (CNRS - UPMC, France), Raphael Voituriez (CNRS - UPMC, France)

We present an exact calculation of the mean first-passage time to a target on the surface of a 2D or 3D spherical domain, for a molecule alternating phases of surface diffusion on the domain boundary and phases of bulk diffusion. The presented approach is based on an integral equation which can be solved analytically. Explicit solutions are provided for normal and biased diffusion in a general annulus with an arbitrary number of regularly spaced targets on a partially reflecting surface. In the framework of this minimal model of surface-mediated reactions, we show analytically that the mean reaction time can be minimized as a function of the desorption rate from the surface. As a consequence, an intermittent exploration process may enhance search and reaction, as compared to pure surface diffusion or pure bulk diffusion. Our method is applicable to extended targets of arbitrary size (i.e., beyond the narrow escape limit). Higher-order moments and the probability distribution of the first-passage time can also be derived.

Related Links

• http://pmc.polytechnique.fr/pagesperso/dg/publi/publi_e.htm - Webpage with publications 

Cell biology processes such as motility, proliferation and death are essential to a host of phenomena such as development, wound healing and tumour invasion, and a huge number of different modelling approaches have been applied to study them. In this talk I will explore a suite of related models for the growth and invasion of cell populations. These models take into account different levels of detail on the spatial locations of cells and, as a result, their predictions can differ depending on the relative magnitudes of the various model parameters. To this end, I will discuss how one might determine the applicability of each of these models, and the extent to which inference techniques can be used to estimate their parameters, using both cell- and population-level quantitative data.

Co-author: Qing Nie (UC Irvine)

Morphogens provide positional information for spatial patterns of gene expression during development. However, stochastic effects such as local fluctuations in morphogen concentration and noise in signal transduction make it difficult for cells to respond to their positions accurately enough to generate sharp boundaries between gene expression domains. In this talk, I will present a novel noise attenuation mechanism during the development of spatial pattern. First, we investigate the boundary sharpening in the zebrafish hindbrain. Computational analyses of spatial stochastic models show, surprisingly, that a combination of noise in RA concentration and noise in hoxb1a/krox20 expression promotes sharpening of boundaries between adjacent segments. Second, we integrate spatial and temporal noise attenuation in the BMP-FGF signaling network in the dorsal telencephalon. We demonstrate that perturbing FGF signaling transiently will lead to a noisier boundary with loss of boundary s harpness and a global delay in development.

In this talk we review recent work with Sean Lawley on diffusion in randomly switching environments. One of the fundamental transport processes in biological cells is the exchange of ions, proteins and other macromolecules between subcellular domains, or between the interior and exterior of the cell, via stochastically gated membrane pores and channels. For example, the nucleus of eukaryotes is surrounded by a protective nuclear envelope within which are embedded nuclear pore complexes (NPCs). The NPCs are the sole mediators of exchange between the nucleus and cytoplasm, which requires the formation of complexes with chaperone molecules known as karyopherins. Other examples include the membrane transport of particles via voltage-gated and ligand-gated ion channels, and intercellular gap-junction coupling. One example at the more macroscopic level is the passive diffusion of oxygen during insect respiration. We show how each of these systems can be modeled in terms of diffusio n in a bounded domain with (partially) switching boundaries, and use a combination of PDE theory and probabilistic methods to determine statistical properties of the system. We highlight important differences between cases where the diffusing particles switch conformational state and cases where the boundary physically switches.

Co-authors: Andrew Seeber (Friedrich Miescher Institute for Biomedical Research), Susan M. Gasser (Friedrich Miescher Institute for Biomedical Research), David Holcman (Ecole Normale Supérieure)

Double-strand break (DSB) repair by homologous recombination (HR) requires an efficient and timely search for a homologous template. Here we developed a statistical method of analysis based on single-particle trajectory data which allows us to extract forces acting on chromatin at DSBs. We can differentiate between extrinsic forces from the actin cytoskeleton and intrinsic alterations on the nucleosomal level at the cleaved MAT locus in budding yeast. Using polymer models we show that reduced tethering forces lead to local decondensation near DSBs, which reduces the mean first encounter time by two orders of magnitude. Local decondensation, likely stems from loss of internal mechanical constraints and a local redistribution of nucleosomes that depends on chromatin remodelers. Simulations verify that local changes in inter-nucleosomal contacts near DSBs would shorten drastically the time required for a long-range homology search.