In dbarrows/bondr: Chemical reaction networks

require(svglite)
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = '#>',
    dev = 'svglite',
    fig.path = 'man/figures/README-',
    fig.showtext = TRUE,
    fig.width = 7,
    fig.height = 4,
    out.width = '100%',
    message = FALSE
)

bondr

Provides utilities and classes for working with reaction networks in R.

Installation

You can install the development version from GitHub with:

devtools::install_github('dbarrows/bondr')

Creating networks

You write systems of reactions using a natural syntax, and bondr will parse it and turn it into an S3 object.

devtools::load_all()

library(bondr)

(synthesis <- network('A + B -> C, 2.4e-5'))

Sources / sinks

You can specify sources / sinks using 0 as the species name.

network('0 -> A, 4')

Multiple reactions

Additional reactions can be entered on new lines.

network('
    S + E -> SE, 1.66e-3
    SE -> E + P, 1e-1
')

Bidirectional reactions

You can use <-> to indicate bidirectional reactions, with an additional rate specified at the end of the line.

network('A <-> B, 1e-1, 2.2')

Species orders

Prefixing a species name with a number will be interpreted as a reaction coefficient.

net <- network('2A -> B, 1e-1')
str(net$reactions[[1]])

Using networks

Propensity functions

You can generate the propensity functions for a reaction network.

net <- network('
         A -> B, 2.5
    2B + C -> A, 4e-2
')
(props <- propensities(net))

Note that dimerisations and multiple reactants are handled properly.

Propensity functions take a state vector of species quantities, ordered according to the output of the species function.

species(net)
state <- c(2, 5, 4) # Corresponds to A, B, C
(props[[1]](state))
(props[[2]](state))

Stoichiometric matrix

A matrix that conveys how the system updates when reactions fire. The columns correspond to reactions, and the rows to species.

mm_string <- network_string_examples('mm')
cat(mm_string)
(net <- network(mm_string))
stmat(net)

Solving

bondr also provides a function deriv to get a derivative function compatible with the deSolve R package, which contains a number of numerical integrators.

Using deSolve

Obtaining a deterministic solution to a system as in the Reaction Rate Equation can be done as follows.

library(deSolve)

y <- c(S = 300, E = 120, SE = 0, P = 0)
times <- seq(0, 30, length.out = 100)
func <- deriv(net)

sol <- ode(y, times, func)
head(sol)

Plotting

You can then plot the solution using a few tidyverse functions fairly easily.

library(tidyverse)
library(wplot)

sol %>%
    data.frame() %>%
    rename(Time = time) %>%
    pivot_longer(species(net), names_to = 'Species', values_to = 'Quantity') %>%
    ggplot(aes(x = Time, y = Quantity, colour = Species)) +
        geom_line() +
        theme_wc()

dbarrows/bondr documentation built on Feb. 7, 2022, 11:01 a.m.