-
Notifications
You must be signed in to change notification settings - Fork 1
/
README.Rmd
232 lines (164 loc) · 10.3 KB
/
README.Rmd
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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
---
title: Package pdmpsim
output:
github_document:
toc: true
toc_depth: 2
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, echo = FALSE, message = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-"
)
library(pdmpsim)
```
# Introduction
The goal of `pdmpsim` is to simulate [piecewise deterministic
Markov processes](https://www.researchgate.net/publication/316281383_Piecewise-deterministic_Markov_processes_A_general_class_of_non-diffusion_stochastic_models_and_Discussion) (PDMPs) within R and to provide methods
for analysing the simulation results.
It is possible to
* simulate PDMPs
* store multiple simulations in a convenient way
* calculate some statistics on them
* plot the results (there are different plot methods available)
* compute the generator numerically
The PDMPs can have multiple discrete and continuous variables.
They are not allowed to have boundaries or a varying number of continuous
variables (the number should be independent of the state of the
discrete variable).
You can install `pdmpsim` from GitHub with:
```{r gh-installation, eval = FALSE}
# install.packages("devtools")
devtools::install_github("CharlotteJana/pdmpsim")
```
# A simple example
This is a simple example modelling gene expression with positive feedback:
```{r}
examplePDMP <- new("pdmpModel",
descr = "Gene regulation with positive feedback",
parms = list(b = 0.5, a0 = 1, a1 = 3, k10 = 1, k01 = 0.5),
init = c(f = 1, d = 1),
discStates = list(d = 0:1),
dynfunc = function(t, x, parms) {
df <- with(as.list(c(x, parms)), {
switch(d+1, a0 - b*f, a1 - b*f)
})
return(c(df, 0))
},
ratefunc = function(t, x, parms) {
return(with(as.list(c(x, parms)), switch(d + 1, k01*f, k10)))
},
jumpfunc = function(t, x, parms, jtype) {
c(x[1], 1 - x[2])
},
times = c(from = 0, to = 100, by = 0.1),
solver = "lsodar")
```
The variable names and initial values are specified in slot **init**. In this case, the model has two variables, `f` and `d`. Slot **discStates** specifies that `d` is the discrete variable and can take the possible values 0 and 1.
Constant parameters of the model are listed in slot **parms**, a short description is given in slot **descr** and is usually used as title for the plot methods. Slot **dynfunc** returns the ODEs of the variables, their value depends on the state of the discrete variable `d`. Slot **ratefunc** determines the rate and therefore probability that a jump occurs (there is only one type of jumps). Slot **jumpfunc** returns the new value after a jump.
# Simulation
A single simulation of a PDMP can be calculated with function `sim`. The function takes the model as argument and returns the same model, with the simulation results stored in a special slot named `out`.
```{r}
out(examplePDMP) # no simulation
```
```{r}
examplePDMP <- sim(examplePDMP)
head(out(examplePDMP)) # simulation stored in slot out
```
```{r}
plot(examplePDMP) # plot the simulation results
```
## Multiple Simulations
Package `pdmpsim` provides two similar methods to perform and store a large number of different simulations of one PDMP.
Function `multSim` returns an S3-object of class
`multSim` which contains a list of simulation results, a list of time values storing
the time needed for the corresponding simulation, the model that was used for the simulations and a vector of numeric numbers. This vector is named `seeds`, its elements are used as argument to function `sim` and control the stochastic part of the model, making the simulation results reproducible. The vector `seeds` and the PDMP model are the only arguments needed for function `multSim`.
```{r, eval = FALSE}
simulations <- multSim(examplePDMP, seeds = 1:300) # 300 simulations, they may take some time
```
The second function available to store multiple simulations is called `multSimCsv`. It is only useful if the memory used by all simulations exceeds the working memory. In this case, it is not possible to store all simulations in one object. They are stored in multiple csv files instead.
# Plot methods
There are several functions available to plot the simulations of a PDMP. Most of them are generated with package `ggplot2` and can be modified afterwards.
## Plot single simulations
The trajectory of a single simulation can be plotted with function `plot`, if the simulation is stored in slot `out` of the model. The generated plot is a base plot, but a plot generated by `ggplot2` is also possible, depending on the Boolean argument `ggplot`.
For objects of class `multSim`, which contain multiple simulations, another function is necessary. It is called `plotSeeds` and can plot up to four different trajectories at a time.
```{r, eval = FALSE}
plotSeeds(simulations, seeds = c(1, 5))
```
![](https://raw.githubusercontent.com/CharlotteJana/pdmpsim/charlotte/man/figures/README-plotSeeds.png)
## Boxplot and violin plot
Function `plotTimes` allows to plot a box plot or violin plot for different time values over all simulations. Argument `plottype` determines the type of the plot.
```{r, eval = FALSE}
plotTimes(simulations, plottype = "violin", times = c(3, 50, 100))
```
![](https://raw.githubusercontent.com/CharlotteJana/pdmpsim/charlotte/man/figures/README-violins.png)
A special argument to function `plotTimes` is called `nolo` which stands for *number of labelled outliers*. This allows to label the most extreme outliers with the seed number which was used to create the simulation. If for example `nolo = 3`, a maximum of three outliers will be labelled for every time value.
```{r, eval = FALSE}
plotTimes(simulations, vars = "f", times = c(20, 50, 80, 100), nolo = 3)
```
![](https://raw.githubusercontent.com/CharlotteJana/pdmpsim/charlotte/man/figures/README-boxplot.png)
## Heatmap over all simulations
If `plot` is called on an object of class `multSim`, it creates
a heat map for every continuous variable and a (possibly smoothed) line plot for the discrete values of every discrete variable. This is the only plot method which gives an overview over all simulated data at every time value.
```{r, eval = FALSE}
plot(simulations, discPlot = "line")
```
![](https://raw.githubusercontent.com/CharlotteJana/pdmpsim/charlotte/man/figures/README-heatmap.png)
## Density plot and histogram
Function `hist` plots a histogram over all simulations for a given time value. Function `density` creates a density plot, allowing multiple time values as input.
```{r, eval = FALSE}
density(simulations, t = c(3, 50, 80, 100))
```
![](https://raw.githubusercontent.com/CharlotteJana/pdmpsim/charlotte/man/figures/README-density.png)
## Calculate and plot statistics
Use function `summarise_at` to apply functions as `mean`, `sd`, etc. on the simulated values. Function `plotStats` calls `summarise_at` internally and plots the results.
```{r, eval = FALSE}
plotStats(simulations, vars = "f", funs = c("min", "mean", "median", "max"))
```
![](https://raw.githubusercontent.com/CharlotteJana/pdmpsim/charlotte/man/figures/README-statistics.png)
# The generator
Package `pmdpsim` provides the function `generator` to numerically calculate the generator of a PDMP without boundaries.
```{r}
Q <- generator(examplePDMP)
```
Because the generator operates on functions, it hast to be applied to a function and afterwards to specific values.
```{r}
f <- function(d, f) d*f^2
Q(f)(t = 10, x = c("d" = 1, "f" = 4))
```
See `?generator` for correct usage.
# A more sophisticated example
The PDMP presented in this paragraph models a gene regulation mechanism with toggle switch between two genes. It is included in package `pdmpsim` as an example, a detailed description of the different slots can be loaded with `?toggleSwitch`.
We consider two genes $A$ and $B$, and the concentration of their gene products, $f_A$ and $f_B$. A product of gene $A$ can block the transcription of gene $B$ and therefore affects the concentration $f_B$. Conversely, a product of gene $B$ can block gene $A$. See the left part of the picture below.
![](https://raw.githubusercontent.com/CharlotteJana/pdmpsim/charlotte/man/figures/README-toggleSwitch.png)
The PDMP contains six different constant parameters: $k₀₁, k₁₀, a_A, a_B, b_A$ and $b_B$. It has two discrete variables $d_A$ and $d_B$, where $d_A$ represents the status of gene $A$ (inactive/active) and $d_B$ represents the status of gene $B$. There are two different type of jumps, the first describes a change in gene $B$, the second a change in gene $A$. The ODEs for $f_A$ and $f_B$ depend on the values of the discrete variables, as can be seen in the right part of the figure above.
The formulation as `pdmpModel` object is as follows:
```{r}
toggleSwitch <- pdmpModel(
descr = "Toggle switch with two promotors",
parms = list(bA = 0.5, bB = 0.5, aA = 2, aB = 4, k01 = 0.5, k10 = 2),
init = c(fA = 0.5, fB = 0.5, dA = 1.0, dB = 1.0),
discStates = list(dA = c(0, 1), dB = c(0, 1)),
times = c(from = 0, to = 100, by = 0.01),
dynfunc = function(t, x, parms) {
df <- with(as.list(c(x, parms)), c(-bA*fA + aA*dA, -bB*fB + aB*dB))
return(c(df, 0, 0))
},
ratefunc = function(t, x, parms) {
return(with(as.list(c(x, parms)), c(switch(dB+1, k01, k10*fA),
switch(dA+1, k01, k10*fB))))
},
jumpfunc = function(t, x, parms, jtype){
return(with(as.list(c(x, parms)), c(fA, fB, switch(jtype,
c(dA, 1-dB),
c(1-dA, dB)))))
})
```
```{r, warning=FALSE, message=FALSE}
plot(sim(toggleSwitch), ggplot = TRUE) +
ggplot2::scale_color_manual(values = c("orange", "palegreen3")) # change color
```
# License
Some classes introduced in package `pdmpsim` are based on code from package [simecol](http://simecol.r-forge.r-project.org/). Both packages are free open source software licensed under the [GNU Public License](https://www.gnu.org/licenses/#GPL) (GPL 2.0 or above). The software is provided as is and comes WITHOUT WARRANTY.