R code help: how to get p1 vs time plot rather than frequency time plot from this code. i am a beginner in coding

General
1 
library(GillespieSSA)
parms <- c(kmu=0.5, kmo = 0.0005, kp=0.167, gamma_m=0.005776, gamma_p=0.001155, kr=1.0, ku1=224.0, ku2=9.0)
tf <- 10000                                     # Final time
simName <- "REPRESSILATOR"  # Name
x0 <- c(m1=10, p1=10, m2=10, p2=10, m3=10, p3=10, n1=0, n2=0, n3=0)

nu <- matrix(c(1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
               0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0,
               0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
               0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1,
               0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
               0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 1, 0, 0,
               0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1, 0,
               0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1,
               0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, -1, 0, 0),
               nrow=9,byrow=TRUE)
          
a  <- c("kmu*(n3 == 0) + kmo*(n3 > 0)", 
        "kmu*(n1 == 0) + kmo*(n1 > 0)", 
        "kmu*(n2 == 0) + kmo*(n2 > 0)", 
        "kp * m1",
        "kp * m2",
        "kp * m3",
        "gamma_m * m1",
        "gamma_m * m2",
        "gamma_m * m3",
        "gamma_p * p1",
        "gamma_p * p2",
        "gamma_p * p3",
        "gamma_p * n1",
        "gamma_p * n2",
        "gamma_p * n3",
        "kr * p3 * (n3 < 2)",
        "kr * p1 * (n1 < 2)",
        "kr * p2 * (n2 < 2)",
        "ku1 * (n3 == 1) + 2 * ku2 * (n3 == 2)",
        "ku1 * (n1 == 1) + 2 * ku2 * (n1 == 2)",
        "ku1 * (n2 == 1) + 2 * ku2 * (n2 == 2)")
set.seed(1)
out <- ssa (
  x0 = x0,
  a = a,
  nu = nu,
  parms = parms,
  tf = tf,
  method = ssa.d(),
  simName = simName,
  verbose = FALSE,
  consoleInterval = 10
) 
ssa.plot(out,show.title = TRUE,show.legend = FALSE)
2

Your code doesn't work for me. It seems to be hanging in ssa(). But maybe it just has a long runtime?

3

I had the same experience as @woodward with tf = 10000, so I trimmed it to permit the example to run.

suppressPackageStartupMessages({
 library(dplyr)
 library(ggplot2)
 library(GillespieSSA)
})

parms <- c(kmu = 0.5, kmo = 0.0005, kp = 0.167, gamma_m = 0.005776, gamma_p = 0.001155, kr = 1.0, ku1 = 224.0, ku2 = 9.0)

# REDUCED FROM 10000 to permit reasonal runtime

tf <- 100 # Final time

simName <- "REPRESSILATOR" # Name

x0 <- c(m1 = 10, p1 = 10, m2 = 10, p2 = 10, m3 = 10, p3 = 10, n1 = 0, n2 = 0, n3 = 0)

nu <- matrix(c(
 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0,
 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1,
 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 1, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, -1, 0, 0
),
nrow = 9, byrow = TRUE
)

a <- c(
 "kmu*(n3 == 0) + kmo*(n3 > 0)",
 "kmu*(n1 == 0) + kmo*(n1 > 0)",
 "kmu*(n2 == 0) + kmo*(n2 > 0)",
 "kp * m1",
 "kp * m2",
 "kp * m3",
 "gamma_m * m1",
 "gamma_m * m2",
 "gamma_m * m3",
 "gamma_p * p1",
 "gamma_p * p2",
 "gamma_p * p3",
 "gamma_p * n1",
 "gamma_p * n2",
 "gamma_p * n3",
 "kr * p3 * (n3 < 2)",
 "kr * p1 * (n1 < 2)",
 "kr * p2 * (n2 < 2)",
 "ku1 * (n3 == 1) + 2 * ku2 * (n3 == 2)",
 "ku1 * (n1 == 1) + 2 * ku2 * (n1 == 2)",
 "ku1 * (n2 == 1) + 2 * ku2 * (n2 == 2)"
)

set.seed(1)

out <- ssa(
 x0 = x0,
 a = a,
 nu = nu,
 parms = parms,
 tf = tf,
 method = ssa.d(),
 simName = simName,
 verbose = FALSE,
 consoleInterval = 10
)

ssa.plot(out, show.title = TRUE, show.legend = FALSE)


# inspect the out object

# str(out)

# List of 3
# $ data : num [1:20037, 1:10] 0 0.0269 0.0273 0.0302 0.0318 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : NULL
# .. ..$ : chr [1:10] "t" "m1" "p1" "m2"

# out is a matrix (which is either all character or all numeric)
# t is presumably time, in column 1, and p1 is in column 3

head(out$data)
#>               t m1 p1 m2 p2 m3 p3 n1 n2 n3
#> [1,] 0.00000000 10 10 10 10 10 10  0  0  0
#> [2,] 0.02692223 10  9 10 10 10 10  1  0  0
#> [3,] 0.02729366 10 10 10 10 10 10  0  0  0
#> [4,] 0.03021189 10  9 10 10 10 10  1  0  0
#> [5,] 0.03181018 10  8 10 10 10 10  2  0  0
#> [6,] 0.09477203 10  8 10  9 10 10  2  1  0

# subset t and p1

sub_out <- data.frame(out$data[, c(1, 3)]) # all rows and columns 1 and 3

head(sub_out)
#>            t p1
#> 1 0.00000000 10
#> 2 0.02692223  9
#> 3 0.02729366 10
#> 4 0.03021189  9
#> 5 0.03181018  8
#> 6 0.09477203  8

p <- ggplot(sub_out, aes(t, p1))
p + geom_point() + theme_minimal()

For a newcomer, there's a lot happening, so come back with questions.

2 Likes
4

sir thanks for the reply...
how to get the p1 p2 p3 plots in the same figure....

5

Here's a start. Come back with a reprex if you need to tweak the color.

suppressPackageStartupMessages({
  library(dplyr)
  library(ggplot2)
  library(GillespieSSA)
})

parms <- c(kmu = 0.5, kmo = 0.0005, kp = 0.167, gamma_m = 0.005776, gamma_p = 0.001155, kr = 1.0, ku1 = 224.0, ku2 = 9.0)

# REDUCED FROM 10000 to permit reasonable runtime

tf <- 100 # Final time

simName <- "REPRESSILATOR" # Name

x0 <- c(m1 = 10, p1 = 10, m2 = 10, p2 = 10, m3 = 10, p3 = 10, n1 = 0, n2 = 0, n3 = 0)

nu <- matrix(c(
  1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0,
  0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1,
  0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 1, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, -1,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1, 0, 0, -1, 0, 0
),
nrow = 9, byrow = TRUE
)

a <- c(
  "kmu*(n3 == 0) + kmo*(n3 > 0)",
  "kmu*(n1 == 0) + kmo*(n1 > 0)",
  "kmu*(n2 == 0) + kmo*(n2 > 0)",
  "kp * m1",
  "kp * m2",
  "kp * m3",
  "gamma_m * m1",
  "gamma_m * m2",
  "gamma_m * m3",
  "gamma_p * p1",
  "gamma_p * p2",
  "gamma_p * p3",
  "gamma_p * n1",
  "gamma_p * n2",
  "gamma_p * n3",
  "kr * p3 * (n3 < 2)",
  "kr * p1 * (n1 < 2)",
  "kr * p2 * (n2 < 2)",
  "ku1 * (n3 == 1) + 2 * ku2 * (n3 == 2)",
  "ku1 * (n1 == 1) + 2 * ku2 * (n1 == 2)",
  "ku1 * (n2 == 1) + 2 * ku2 * (n2 == 2)"
)

set.seed(1)

out <- ssa(
  x0 = x0,
  a = a,
  nu = nu,
  parms = parms,
  tf = tf,
  method = ssa.d(),
  simName = simName,
  verbose = FALSE,
  consoleInterval = 10
)

sub_out <- data.frame(out$data[, c(1,3,5,7)]) # all rows and columns 1 and 3,5,7
p <- ggplot(sub_out)
p + 
  geom_point(aes(t,p1, color = p1)) +
  geom_point(aes(t,p2, color = p2)) +
  geom_point(aes(t,p3, color = p3)) +
  theme_minimal()

1 Like
6

Sir,
I need different colours for the plot.. and to check if am getting a reasonable output..i need to run it for a minimum of 9000secs...what to do for that. TIA.

7

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.