Description Usage Arguments Details Value See Also Examples
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.
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)
|
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. |
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.
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. |
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)
}
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.