Simulation of a biochemical system

Share:

Description

These functions simulate a biochemical reacton system parameterized as a Petri Net. GillespieOptimDirect, GillespieDirectGB, GibsonBruck, and GillespieDirectCR performs pure stochastic simulations, RungeKuttaDormandPrince45 a pure deterministic integration, HaseltineRawlings a hybrid of the above. PartitionedLeaping a dynamic-repartitioning simulation. Multiple runs can be performed at once.

See init for a way of defining the model that is close to the way reactions are written.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
## Exact stochastic simulation: 
GillespieOptimDirect(model, timep, delta=1, runs=1)
GillespieDirectGB(model, timep, delta=1, runs=1)
GibsonBruck(model, timep, delta=1, runs=1)
GillespieDirectCR(model, timep, delta=1, runs=1)

## Pure deterministic:
RungeKuttaDormandPrince45(model, timep, delta=1, ect = 1e-09)

## Hybrid stochastic/deterministic:
HaseltineRawlings(model, timep, delta=1, runs=1, ect = 1e-09)

## Dynamic re-partitioning:
PartitionedLeaping(model, timep, delta=1, runs=1, ect = 1e-09)

Arguments

model

list containing named elements:

timep

It can be either a numeric, indicating for how long (in the same time units as the propensity constants) the process will run, or a functions (R or C), in which case can be used to change the protocol at time intervals. See details.

delta

Interval time at which the state will be saved.

runs

How many runs will be performed.

ect

Precision for the fast reactions.

Details

model is a list containing the following elements:

  • model$pre: pre matrix, with as many rows as transitions (reactions), and columns as places (reactants). It has the stoichiometrics of the left sides of the reactions.

  • model$post: post matrix, with as many rows as transitions, and columns as places (products). It has the stoichiometrics of the right sides of the reactions.

  • model$h: list of propensity constants or functions returning the propensity (with as many elements as transitions).

  • model$slow: vector of zeros for slow transitions and ones for fast transitions. Only needed for HaseltineRawlings. Ignored otherwise.

  • model$M: initial marking (state) of the system.

  • model$place: vector with names of the places.

  • model$transition: vector with names of the transitions.

Value

The functions return a list with the following elements:

place

vector with the names of the places if supplied. If not, the function creates names as follows: P1, P2, ...

transition

vector with the names of the transitions if supplied. If not, the function creates names as follows: T1, T2, ...

dt

vector containing the discretized times at which the state is saved (according to delta)

run

list with as many elements as runs. We will describe the first element, run[[1]], as the rest have exactly the same structure. It is also a list, with the following elements:

run[[1]]$M

list with as many elements as places, each of them containing the state of the system sampled according to delta.

run[[1]]$transitions

vector with as many elements as transitions, with the total of time each slow reaction fired.

run[[1]]$tot.transitions

numeric with the summ of run[[1]]$transitions.

See Also

init, atr

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
## bioPN has been tested only on 64 bits machines.
## It may fail in 32 bits architecture.
if (.Machine$sizeof.pointer == 8) {

####### Reaction constants
H <- 10
K <- 6
r <- 0.25
c <- 3
b <- 2
#######

Gi <- 1
Ga <- 2
mRNA <- 3
Protein <- 4

model <- list(
    pre=matrix(c(1,0,0,0, 0,1,0,0, 0,1,0,0,
        0,0,1,0, 0,0,1,0, 0,0,0,1),
        ncol=4, byrow=TRUE),
    post=matrix(c(0,1,0,0, 1,0,0,0, 0,1,1,0,
        0,0,0,0, 0,0,1,1, 0,0,0,0),
        ncol=4, byrow=TRUE),
    h=list(c, b, H, 1, K, r),
    M=c(1,0,0,0))

timep <- 200
delta <- 1

##############################
## Completely Deterministic ##
##############################
Sim <- RungeKuttaDormandPrince45(model, timep, delta)

## Note, it could also be done as follows
## slow <- rep(0, transitions)
## Sim <- HaseltineRawlings(model, timep, delta, runs = 1)

mRNA.run <- Sim$run[[1]]$M[[mRNA]]
protein.run <- Sim$run[[1]]$M[[Protein]]

## Theoretical results (red lines in following plots)
Mean.mRNA <- c/(c+b)*H
Mean.protein <- Mean.mRNA * K/r

par(mfrow=c(1,2))
par(mar=c(2, 4, 2, 1) + 0.1)
plot(Sim$dt, mRNA.run,type="l", ylab="Mean",main="mRNA")
legend(x="bottom", paste("Deterministic run"))
abline(h=Mean.mRNA,col="red", lwd=1)
plot(Sim$dt, protein.run,type="l", ylab="Mean",main="Protein")
legend(x="bottom", paste("Deterministic run"))
abline(h=Mean.protein,col="red", lwd=1)


runs <- 100	## Increase to 10000 for better fit
###########################
## Completely Stochastic ##
###########################
set.seed(19761111)  ## Set a seed (for reproducible results)
Sim <- GillespieOptimDirect(model, timep, delta, runs)

## Note, it could also be done as follows
## slow <- rep(1, transitions)
## Sim <- HaseltineRawlings(model, timep, delta, runs)

mRNA.run <- sapply(Sim$run, function(run) {run$M[[mRNA]]})
protein.run <- sapply(Sim$run, function(run) {run$M[[Protein]]})

## Histograms of protein at different time points.
par(mfrow=c(2,2))
par(mar=c(2, 4, 2.5, 1) + 0.1)
hist(protein.run[Sim$dt == 1,], main="Protein Distribution at t=1sec")
hist(protein.run[Sim$dt == 2,], main="Protein Distribution at t=2sec")
hist(protein.run[Sim$dt == 10,], main="Protein Distribution at t=10sec")
hist(protein.run[Sim$dt == 200,], main="Protein Distribution at t=200sec")

## Theoretical results (red lines in following plots)
Mean.mRNA <- c/(c+b)*H
Var.mRNA <- b/(c*(1+c+b))*Mean.mRNA^2 + Mean.mRNA
Mean.protein <- Mean.mRNA * K/r
Var.protein <- r*b*(1+c+b+r)/(c*(1+r)*(1+c+b)*(r+c+b))*Mean.protein^2 +
  r/(1+r)*Mean.protein^2/Mean.mRNA + Mean.protein

if (runs > 1 ) {
  par(mfrow=c(2,2))
} else {
  par(mfrow=c(1,2))
}
par(mar=c(2, 4, 2, 1) + 0.1)
plot(Sim$dt, apply(mRNA.run,1,function(tpt) {mean(tpt)}),type="l", ylab="Mean",main="mRNA")
legend(x="bottom", paste("Gene, mRNA and Protein Stochastic\nRuns :", runs))
abline(h=Mean.mRNA,col="red", lwd=1)
plot(Sim$dt, apply(protein.run,1,function(tpt) {mean(tpt)}),type="l", ylab="Mean",main="Protein")
legend(x="bottom", paste("Gene, mRNA and Protein Stochastic\nRuns :", runs))
abline(h=Mean.protein,col="red", lwd=1)
if (runs > 1 ) {
  par(mar=c(2, 4, 0, 1) + 0.1)
  plot(Sim$dt, apply(mRNA.run,1,function(tpt) {var(tpt)}),type="l", ylab="Var")
  abline(h=Var.mRNA,col="red", lwd=1)
  plot(Sim$dt, apply(protein.run,1,function(tpt) {var(tpt)}),type="l", ylab="Var")
  abline(h=Var.protein,col="red", lwd=1)
}


######################################################################
## Hybrid: mRNA and protein fast, gene activation/inactivation slow ##
######################################################################
model$slow <- c(1,1,0,0,0,0)

Sim <- HaseltineRawlings(model, timep, delta, runs)

mRNA.run <- sapply(Sim$run, function(run) {run$M[[mRNA]]})
protein.run <- sapply(Sim$run, function(run) {run$M[[Protein]]})
    
Mean.mRNA <- c/(c+b)*H
Var.mRNA <- b/(c*(1+c+b))*Mean.mRNA^2
Mean.protein <- Mean.mRNA * K/r
Var.protein <- r*b*(1+c+b+r)/(c*(1+r)*(1+c+b)*(r+c+b))*Mean.protein^2

if (runs > 1 ) {
  par(mfrow=c(2,2))
} else {
  par(mfrow=c(1,2))
}
par(mar=c(2, 4, 2, 1) + 0.1)
plot(Sim$dt, apply(mRNA.run,1,function(tpt) {mean(tpt)}),type="l", ylab="Mean",main="mRNA")
legend(x="bottom", paste("Only Gene Stochastic\nRuns :", runs))
abline(h=Mean.mRNA,col="red", lwd=1)
plot(Sim$dt, apply(protein.run,1,function(tpt) {mean(tpt)}),type="l", ylab="Mean",main="Protein")
legend(x="bottom", paste("Only Gene Stochastic\nRuns :", runs))
abline(h=Mean.protein,col="red", lwd=1)
if (runs > 1 ) {
  par(mar=c(2, 4, 0, 1) + 0.1)
  plot(Sim$dt, apply(mRNA.run,1,function(tpt) {var(tpt)}),type="l", ylab="Var")
  abline(h=Var.mRNA,col="red", lwd=1)
  plot(Sim$dt, apply(protein.run,1,function(tpt) {var(tpt)}),type="l", ylab="Var")
  abline(h=Var.protein,col="red", lwd=1)
}
}