Expokit provides a set of routines aimed at computing matrix exponentials. More precisely, it computes either a small matrix exponential in full, the action of a large sparse matrix exponential on an operand vector, or the solution of a system of linear OBEs with constant inhomogeneity. The backbone of the sparse routines consists of matrix-free Krylov subspace projection methods (Arnoldi and Lanczos processes), and that is why the toolkit is capable of coping with sparse matrices of large dimension. The software handles real and complex matrices and provides specific routines for symmetric and Hermitian matrices. The computation of matrix exponentials is a numerical issue of critical importance in the area of Markov chains and furthermore, the computed solution is subject to probabilistic constraints. In addition to addressing general matrix exponentials, a distinct attention is assigned to the computation of transient states of Markov chains.

Abstract. We present numerical methods for both the direct solution and simulation of the chemical master equation (CME), and, compared to popular methods in current use, such as the Gillespie stochastic simulation algorithm (SSA) and τ -Leap approximations, this new approach has the advantage of being able to detect when the system has settled down to equilibrium. This improved performance is due to the incorporation of information from the associated CME, a valuable complementary approach to the SSA that has often been felt to be too computationally inefficient. Hybrid methods, that combine these complementary approaches and so are able to detect equilibrium while maintaining the efficiency of the leap methods, are also presented. Amongst CME-solvers the recently suggested finite state projection algorithm is especially well suited to this purpose and has been adapted here for the task, leading to a type of "exact τ -Leap." It is also observed that a CMEsolver is often more efficient than an SSA or even a τ -Leap approach for computing moments of the solution such as the mean and variance. These techniques are demonstrated on a test suite of five biologically inspired models, namely, stochastic models of the genetic toggle, receptor oligomerization, the Schlögl reactions, Goutsias' model of regulated gene transcription, and a decaying-dimerizing reaction set. For the gene toggle it is observed that important experimentally measurable traits such as the percentage of cells that undergo so-called switching may also be more efficiently approximated via a CME-based approach. Key words. chemical master equation, stochastic simulation algorithm, systems biology AMS subject classifications. 60H35, 65C40, 65F10DOI. 10.1137/060678154 1. Introduction. Recent and rapid advances in molecular biology promise great benefits in areas such as medicine and agriculture, and computational biology is seen as having an important role to play in these advances by modeling cell biology at a systems level [10]. Modeling such vastly complicated systems as living cells is inherently a multiscale exercise due to the vast range of spatial and temporal scales on which these processes occur [8]. This paper focuses on the development of multiscale computational methods for coping with this complexity. In particular, this paper focuses on computational methods for models of gene regulatory networks (GRNs) as continuous-time, discrete-state, Markov processes. We note in passing that there are many other kinds of modeling frameworks for GRNs, such as partial differential equations (PDEs) or Kauffman's Boolean networks [8,9].GRNs have been successfully modeled by Markov processes, a notable example being that of the bacteriophage λ life cycle [2], and, in this setting, a GRN is modeled via the collection of biochemical reactions of which it is composed. Although under some circumstances, chemical kinetics have been well modeled by ordinary differential equations (ODEs), under many circumstances this is not appropriate. For example, there are on...

Abstract. The suprathermal particles, electrons and protons, coming from the magnetosphere and precipitating into the high-latitude atmosphere are an energy source of the Earth's ionosphere. They interact with ambient thermal gas through inelastic and elastic collisions. The physical quantities perturbed by these precipitations, such as the heating rate, the electron production rate, or the emission intensities, can be provided in solving the kinetic stationary Boltzmann equation. This equation yields particle fluxes as a function of altitude, energy, and pitch angle. While this equation has been solved through different ways for the electron transport and fully tested, the proton transport is more complicated. Because of charge-changing reactions, the latter is a set of two-coupled transport equations that must be solved: one for protons and the other for H atoms. We present here a new approach that solves the multistream proton/hydrogen transport equations encompassing the collision angular redistributions and the magnetic mirroring effect. In order to validate our model we discuss the energy conservation and we compare with another model under the same inputs and with rocket observations. The influence of the angular redistributions is discussed in a forthcoming paper.

Recently the application of the quasi-steady-state approximation (QSSA) to the stochastic simulation algorithm (SSA) was suggested for the purpose of speeding up stochastic simulations of chemical systems that involve both relatively fast and slow chemical reactions [Rao and Arkin, J. Chem. Phys. 118, 4999 (2003)] and further work has led to the nested and slow-scale SSA. Improved numerical efficiency is obtained by respecting the vastly different time scales characterizing the system and then by advancing only the slow reactions exactly, based on a suitable approximation to the fast reactions. We considerably extend these works by applying the QSSA to numerical methods for the direct solution of the chemical master equation (CME) and, in particular, to the finite state projection algorithm [Munsky and Khammash, J. Chem. Phys. 124, 044104 (2006)], in conjunction with Krylov methods. In addition, we point out some important connections to the literature on the (deterministic) total QSSA (tQSSA) and place the stochastic analogue of the QSSA within the more general framework of aggregation of Markov processes. We demonstrate the new methods on four examples: Michaelis-Menten enzyme kinetics, double phosphorylation, the Goldbeter-Koshland switch, and the mitogen activated protein kinase cascade. Overall, we report dramatic improvements by applying the tQSSA to the CME solver.

No abstract