Team:BCCS-Bristol/BSim/Features/ODE solvers

GRNs
BSim 2008: "Each of the modelling approaches [GRNs and agent-based modelling] have been considered in separate contexts, mainly due to the differing aspects of the system they are concerned with. Now, having working models for each, it would be possible to bring these together with the aim of improving simulation accuracy and allowing for the internal cellular dynamics to be studied in an ever changing physical environment. Such a hybrid model may also help shed light on the critical aspects of project as a whole."

BSim 2009 provides a robust implementation of the second and fourth order Runge-Kutta methods for systems of ordinary differential equations. It is possible to easily specify systems of ODEs as objects within the simulation. These ODE systems can be "attached" to objects in the simulation if necessary and can be used to simulate any aspect of the environment to which they are coupled, depending on the user's requirements. An example would be attaching an ODE system to each bacterium and coupling these systems via an external chemical field. See the overview of our ongoing quorum-coupled repressilators simulation for an example application of this.

As a result of the modular nature of the solver implementation it would also be possible to implement stochastic ODEs, and delay differential equations in a similar manner. These features are likely to be implemented soon to assist with the modelling of more complex GRN systems across a population.

Solvers
Currently we can specify custom ODEs by creating a class that implements BSimOdeSingle or BSimOdeSystem.

The custom classes can be passed directly to any solver, however the solvers are only single-step to make them more adaptable, therefore implementing actual iteration along a time series is left up to the programmer.

The interfaces are designed such that a user defined set of equations can be an external class or in fact an inner class. In general it is useful to define GRNs or other ODEs as inner classes as part of the simulation definition, mostly for reasons of convenience when coupling systems or accessing variables. However, in general the choice of implementation of a specific ODE or system of equations is left up to the user.

BSimOdeSingle
This interface allows you to specify a single ODE. The interface requires two methods, one to define the derivatives and one to return the initial conditions for the initial value problem.

An Example ODE
To define the problem

dy/dt = t

One would create a new class defining the ODE with the following code:

//Package, Imports etc here

// An example ODE public class Linear implements BSimOdeSingle { public double derivative(double x, double y){ dy = x;		return dy; }

public double getIC{ return 0.0; } }

The class Linear would then be passed to a solver along with current x (in this case x is in fact t) and y values. However in complex biological systems such as GRNs which oscillate, there are likely to be multiple ODEs, therefore it is reccommended that BSimOdeSystem be employed in general.

BSimOdeSystem
With this interface it is possible to define a system of multiple (potentially coupled) ODEs. In a similar manner to BSimOdeSingle, it is necessary to specify the derivatives and the initial conditions, however they must now be in the form of arrays as there are multiple values to consider. There is also a third method which must be specified which returns the number of equations present in the system. The implementation of these methods is left up to the user.

An Example ODE System
As an example, we show here how one would define a simple (damped) harmonic oscillator.

This is a second order ODE and so must be simplified to a system of first order equations via a substitution. One can then specify the initial conditions.

The Java implementation of this system is illustrated in the following segment of code:

//Package, Imports etc here

// An example ODE system - Damped Simple Harmonic Motion public class SHM implements BSimOdeSytem{ private int numEq = 2;

public int getNumEq{return numEq ;}

public double[] derivativeSystem(double x, double[] y){ double[] dy = new double[numEq]; dy[0] = y[1]; dy[1] = -2*y[0] - 0.1*y[1]; return dy; }	public double[] getICs{ double[] y0 = new double[numEq]; y0[0] = 0; y0[1] = 1; return y0; } }

The class SHM can then be passed to the solver of choice at each time step, and the values of y extracted and either logged or displayed to the user.