Department of Information Technology

# Past research endeavours of the Computational Systems Biology group

## Internal state model

When molecules move by diffusion in a cytosol in a cell, they are sometimes found to behave anomalously. The mean square displacement of the molecule does not grow linearly with time as expected from Brownian motion. The molecule is slower and we observe subdiffusion or anomalous diffusion. A stochastic, mesoscopic model for
subdiffusion is found in [1] with internal, hidden states. At a macroscopic level, the concentrations of the molecules satisfy a fractional partial differential equation. The internal state model is extended to chemical reactions in [2]. Different models for the reactions are possible. In some of these models, the reaction kinetics is also anomalous in a certain time interval.

1. E W Mommer, D Lebiedz, SIAM J Appl Math, 70
(2009), 112-132.
2. E Blanc, S Engblom, A Hellander, P Lötstedt,
Mesoscopic modeling of stochastic reaction-diffusion
kinetics in the subdiffusive regime, manuscript, 2015.

## Macroscopic, mesoscopic, and microscopic models in molecular biology

The most accurate models for the chemical reactions in biological cells are based on differential equations. The reaction rate equations for the concentrations of the molecular species are a system of nonlinear ordinary differential equations and the macroscopic model equations. They can be solved using standard mathematical software but they are too inaccurate for certain problems e.g. involving a small number of molecules of each kind. This is often the situation in molecular biology. For such systems a more precise mesoscopic, stochastic modelling is necessary. The master equation of chemical kinetics is an equation for the probability density of the distribution of molecules. A computational difficulty with the master equation is the growth in the number of dimensions of the equations with the number of molecular species involved in the reactions. Standard numerical methods for solving these equations suffer from an exponential growth in computational work and memory.

A harder problem is the solution of the reaction-diffusion master equation where the species are not well-stirred but have a space dependence. There the master equation may depend on a state vector with thousands of dimensions making computational solution of the equation impossible. A code for simulation of such systems on unstructured meshes with a Monte Carlo method is found here. The method has been generalized to active transport problems and an adaptive switch between macroscopic and mesoscopic diffusion. A microscopic model where single molecules are tracked is necessary for some systems. Then the molecules follow Brownian dynamics and react with a certain probability with each other when they are close. An efficient method avoiding the small time steps in straightforward Brownian motion is GFRD [1]. The method is simplified and made more flexible by introducing operator splitting.

The Stochastic Simulation Algorithm (SSA) due to Gillespie [2] is widely used to simulate biochemical reaction networks modeled as a continuous-time discrete-space Markov process. CellMC is an XSLT-based, automated SBML (Systems Biology Markup Language) model compiler capable of producing very efficient SSA executables for the Cell/BE or multicore x86 PCs. CellMC was developed by Emmet Caulfield as a part of his master's thesis: "CellMC: An XSLT-based SBML model compiler for Cell/BE and IA32'.

1. J. S. van Zon and P. R. ten Wolde, Green's-function reaction dynamics: A particle-based approach for simulating biochemical networks in time and space, J. Chem. Phys., 123 (2005), art. 234910.
2. D. T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, J. Comput. Phys., 22 (1976), pp. 403-434.

Living organisms have to adapt their behavior to different periodic changes in their environment e.g. the daily variation of light and the annual variation of temperature. Many of them have developed molecular clocks as their internal time-keeping mechanisms. An example of such an oscillatory process is the circadian rhythm with a period of about 24 h [3]. A model for this rhythm has been developed in [4] and is simplified in [5]. This model has two molecular species and the corresponding partial differential equation is two dimensional.

PDE approximation Monte Carlo
The Fokker-Planck equation of the simplified model is solved with a numerical method for partial differential equations with space adaptivity. A description of this method is found here. The chemical master equation is solved by the stochastic simulation algorithm (SSA) [1]
These two movies show the behavior of the probability density function in one period of the oscillation computed by two different numerical methods. A comparison between the methods is found here.

The original model has nine species and a reduction of the complexity is necessary with standard discretization methods. By combining the reaction rate equations on the macroscale with the Fokker-Planck equation on the mesoscale a two dimensional problem is obtained.

Dimension reduction
This movie shows one period of the projection of the probability function of the reduced problem on the same subspace as in the example above. This method is described here.

3. A. Golbeter, Computational approaches to cellular rhythms, Nature, 420 (2002), pp. 238-245.
4. N. Barkai, S. Leibler, Circadian clocks limited by noise, Nature, 403 (2000), pp. 267-268.
5. J. M. G. Vilar, H. Y. Kueh, N. Barkai, S. Leibler, Mechanisms of noise-resistance in genetic oscillators, Proc. Natl. Acad. Sci., 99 (2002), pp. 5988-5992.

Reaction-diffusion process
The cell with a nucleus is discretized by a tetrahedral mesh. A snapshot of the spatial distribution of one species in the cytosol (left) and another species in the nucleus (right). Red represents high copy number and blue low. The method for simulation of stochastic reaction-diffusion processes is described here.
Updated  2017-02-04 14:30:40 by Kurt Otto.