A couple of years ago, I was interested in a pharmacokinetic/pharmacodynamics (PK/PD) modeling and simulation project developed at the institute where I worked. It seemed that Matlab and its SimBiology toolbox was the tool everybody around me went with. While I worked with Matlab before and consider it a great tool, it is not free. Trying to develop a model that may be of interest to many researchers, but limiting its use use to only those who own a Matlab licence, kind of disagreed with my research vision, so I started looking for the open source and free ways to do it.

First, I remembered the xCellerator package, an open source package for biological modeling, that we used in one of my grad school classes. While this package provides an automatic conversion of chemical reactions into ordinary differential equations (ODEs), ability to solve them via numerical integration, and many other great and useful options, it requires installation of Mathematica, which is not free. So I had to look for another solution.

I often use R for my other research projects, so the next logical step was to check R package repositories for any type of biological system modeling and simulation packages. I was surprised when I discovered that there are not many packages for systems biology or PK/PD modeling and none that were general and would allow creating a model from high level reactions equations. This was little bit disappointing. However, I found many examples of how to simulate a system modeled in a form of ODEs and I was able to successfully run a small simulation of my own model. This was nice, but I was not sure how convenient this approach is for a bigger systems – one would still need to convert all reactions into equations manually – it would take time, but also would be prone to errors. I kept wondering how hard would it be to create functions that would allow me to make a model from high level reactions equations and to solve them. So I decided to try to do that (and finally found some time to finish it).

In this post, I will give a high level description of the main functions I developed and show examples of how to create and simulate a model of a system. The functions are packed in a form of a package called sysBio (I eventually plan to release this as a package), but a number of things still need to be added or changed, so keep that in mind if you decide to use it as it is. I did a number of tests, but additional testing also needs to be done. I also plan to expand functionalities, but those will probably come after the package is cleared up and finalized (but probably before it is released).

So…

Molecular reactions can be described by differential equation that define the rate of change in molecular concentrations. There are multiple ways in which reaction can be transformed into a differential equation. I decided to use the law of mass action because of its simplicity and the fact that it has been validated experimentally many times. Given a reaction:

,

where A and B represent reactants and AB represents a product (in general, we can call A, B, and AB species), r represents the reaction rate, and k1, k2, and k3 represent stoichometric coefficients, we can use the law of mass action to determine the change in the amount of each species as:

.

Or in general, given a reaction , we can represent changes in the amount of each species as .

sysBio allows users to write reactions in two formats – as forward reactions, e.g., , and as reversible reactions, e.g.,. Such reversible reactions would be transformed and interpreted as two forward reactions and . sysBio uses the “->” symbol for forward reaction and the “=” symbol for the reversible reaction.

Reaction rate can be defined as a fixed numerical value (e.g., r=l, where l represents a constant parameter, l=0.9, or just as r=0.9) or initial assignments (e.g., r=l*A, where l represents a constant parameter and A represents a species). In the current version, we assume that parameters, if defined, will remain constant during simulation. Values that represent stoichometric coefficients should be explicitly included in the reaction definition (e.g., 2*A + B -> 1.5*AB).

The sysBio package also allows user to define rules. Rules can be used to specify values for model species based on the values of other components (mainly species) in the model. Current version allows only for the ODEs type of rules.

Model must contain at least one species. Initial values for all species, as well as all rates and parameters have to be defined. sysBio uses a term null to describe a null species (source or sink).

sysBio allows users to simulate the model (i.e., solve ODEs) using a differential equation solver based on a combination of backward differentiation formula and a direct linear system solution method (using the “daspk” function from the “deSolve” package) and using Gillespie stochastic simulation algorithm (using the “ssa” function from the “GillespieSSA” package).

If you want to try the sysBio package, you need to install it:

devtools::install_github("Vessy/sysBio") library("sysBio")

Now let’s try a few examples.

Let’s start with a simple reaction such as synthesis, which can be modeled as a zero order reaction.
# Synthesis - Zero order reaction exmp1 <- newModel("Synthesis") addMAreaction(exmp1, "null -> P", "k") addSpecies(exmp1, "P", 0) addMAreactRate(exmp1, "k", "fixed", 1) makeModel(exmp1) simResults <- simulateModel(exmp1) plotResults(simResults, title="Synthesis")

We use a command newModel to create a new model called exmp1. A type of model object is a list and it contains (or will contain after we populate and make model) information about model name, reactions, species, rates, parameters, rules, models, ODEs, stochastic matrix, and stochastic model (propensity function). After we created a model, we can start populating it. Using the addMAreaction function, we can add one or more reactions into the model - these reactions will be interpreted using the law of mass action. Next, we have to specify reaction rates and parameters (if needed).  The addMAreactRate and addParameters functions are used to specify information about the reaction rates and parameters involved in the model. Finally, we need to define species using the addSpecies function. Once we defined all components of a model, we can use function makeModel to create a mathematical representation of the model. This function will transform reactions into corresponding ODEs, and create stochastic matrix and propensity function. Next, we use the simulateModel function to run simulation (solve ODEs). This function calls the validateModel function that checks if all components of the models have been defined. Finally, we use plotResults( function to visualize simulation results. Figure one shows the simulation results.

Figure 1: Synthesis (zero order reaction) model

Now let's try a little bit more complex model, a gene regulation model. Gene regulation process can be described by six reactions: transcription, translation, binding, unbinding, mRNA degradation, and protein degradation (Figure 2a).
# Create a model exmp3 <- newModel("Gene Regulation") # Add reactions addMAreaction(exmp3, "DNA -> DNA + mRNA", r1="v1", name="Transcription") addMAreaction(exmp3, "mRNA -> mRNA + protein", r1="v2", name="Translation") addMAreaction(exmp3, "DNA + protein -> DNA_protein", r1="v3", name="Binding") addMAreaction(exmp3, "DNA_protein -> DNA + protein", r1="v4", name="Unbinding") addMAreaction(exmp3, "mRNA -> null", r1="v5", name="Degradation: mRNA") addMAreaction(exmp3, "protein -> null", r1="v6", name="Degradation:protein") # Add reaction rates addMAreactRate(exmp3, "v1", "assigned", "k1*DNA") addMAreactRate(exmp3, "v2", "assigned", "k2*mRNA") addMAreactRate(exmp3, "v3", "assigned", "k3*DNA*protein") addMAreactRate(exmp3, "v4", "assigned", "k3r*DNA_protein") addMAreactRate(exmp3, "v5", "assigned", "k4*mRNA") addMAreactRate(exmp3, "v6", "assigned", "k5*protein") # Add species addSpecies(exmp3, "DNA", 50) addSpecies(exmp3, "mRNA", 0) addSpecies(exmp3, "protein", 0) addSpecies(exmp3, "DNA_protein", 0) # Add parameters addParameters(exmp3, "k1", 0.2) addParameters(exmp3, "k2", 20) addParameters(exmp3, "k3", 0.2) addParameters(exmp3, "k3r", 1) addParameters(exmp3, "k4", 1.5) addParameters(exmp3, "k5", 1) # Make a model makeModel(exmp3) #printInfo(exmp3) # Run model simulation and plot the results simResults <-simulateModel(exmp3) plotResults(simResults, title="Gene regulation")

Figure 2b shows the simulation results. To see more details about the model, one can use the printInfo command:
> printInfo(exmp3) [1] "Model name: Gene Regulation" [1] "Number of reactions: 6" [1] "Reaction type(s): Mass Action" [1] "Number of species: 4" [1] "Model can be represented in form of 4 differential equation(s)"

To see the mathematical representation of the model, one can print it directly from the model object (after the makeModel command was executed):
> exmp3$model [1] "d_DNA <- -v3*DNA*protein + v4*DNA_protein" [2] "d_DNA_protein <- v3*DNA*protein - v4*DNA_protein" [3] "d_mRNA <- v1*DNA - v5*mRNA" [4] "d_protein <- v2*mRNA - v3*DNA*protein + v4*DNA_protein - v6*protein"

To run a stochastic simulation on the created model, we will use solveStoch function. Figure 2c shows the results of stochastic simulation.
simResults.stoch <- solveStoch(exmp3,10) plotResults(simResults.stoch, title="Gene regulation (stochastic simulation)")

Figure 2: Gene regulation

More model examples and details about package functions can be found at sysBio github wiki page, as well as in the package vignette and help files.

At the end, a few words about future plans. The priority will be given to cleaning the current code and finishing testing. After that, I'll work on adding new features. Features that I plan to add are additional rules to transform reaction into a differential equation (probably Michaelis-Menten rule and Hill equation), ability to process reactions in the form of cascade, events, additional types of reaction rates, additional types of rules, steady state calculation, sensitivity analysis, and parameter estimation. I also plan to include an option for different ODEs solver and I may need to look for or implement my own version of stochastic simulation that would support events and different types of rules. If I get really motivated and have extra time, I may also add an option for multiple compartments within a model and/or provide some type of visualization (probably using Shiny).