Skip to contents

Formula-based simulation constraints

Source: vignettes/run_until.Rmd
run_until.Rmd

Note: This article presents advanced topics on tissue simulation. Refer to this article for an introduction to the subject.

Disclaimer: ProCESS/CLONES internally implements the probability distributions using the C++11 random number distribution classes. The standard does not specify their algorithms, and the class implementations are left free for the compiler. Thus, the simulation output depends on the compiler used to compile CLONES, and because of that, the results reported in this article may differ from those obtained by the reader.

ProCESS implements a first-order unquantified logic having variables representing the cardinality of the species, the number of events fired in a species (being duplications, deaths, or switches), and the simulation time. These variables and real values are summed by +, subtracted by -, and multiplied by * and form expressions. The expressions are then compared with the standard semantics by >, >=, ==, !=, <=, and < to form relations. A formula in this language is either a relation, the conjunction of two formulas (&), or the disjunction of two formulas (|).

Any formula in the above language expresses a condition on the simulation status. It can be used as a parameter of the method TissueSimulation$run_until() to let the simulation evolve until the condition no longer holds.

Variables

The variables represent one of the following quantities:

  • the cardinality of a species;
  • the number of fired deaths, duplications, or switches in a species;
  • the elapse simulation time.

All the above variables can be built by using the method TissueSimulation$var(). When the parameter is the string "Time", the elapsed simulation time variable is returned.

library(ProCESS)

# set the seed of the random number generator
set.seed(0)

# build a tissue simulation and add two species to it
sim <- TissueSimulation()
sim$add_mutant(name = "A",
               epigenetic_rates = c("+-" = 0.01, "-+" = 0.01),
               growth_rates = c("+" = 0.2, "-" = 0.08),
               death_rates = c("+" = 0.1, "-" = 0.01))

# get the variable representing the simulation time
v_time <- sim$var("Time")

v_time
#> Time

When the parameter is the name of a species, a variable representing the cardinality of the species is built.

# get the variable representing the cardinality of A+ in sim
va_p <- sim$var("A+")
va_p
#> |A+|

# get the variable representing the cardinality of A- in sim
va_m <- sim$var("A-")
va_m
#> |A-|

Finally, when the parameter is the name of a species followed by a . and the name of an event among deaths, duplications, or switches, TissueSimulation$var() returns the variable associated with the number of the corresponding event in the specified species.

# get the variable representing the number of epigenetic
# switches from A+
va_ps <- sim$var("A+.switches")
va_ps
#> |A+.switches|

# get the variable representing the number of duplications in A+
sim$var("A+.duplications")
#> |A+.duplications|

# get the variable representing the number of deaths in A+
sim$var("A+.deaths")
#> |A+.deaths|

Expressions and Formulas

An expression is one of the following objects:

  • a variable, e.g., sim$var("A+");
  • a numeric value, e.g., 3.4;
  • the sum of two expressions, e.g., sim$var("A+") + 3.4;
  • the subtraction of two expressions, e.g., sim$var("A+") - 3.4;
  • the multiplication of two expressions, e.g., sim$var("A+") * 3.4.

Two expression can be related by <=, <, ==, !=, > and >=.

A formula is:

  • a relation among two expressions, e.g., sim$var("A+")>=2;
  • the conjunction of two formulas, e.g., sim$var("A+")>=2 & sim$var("A+")<=500;
  • the disjunction of two formulas, e.g., sim$var("A+")>=2 | sim$var("A+")<=500.

The method TissueSimulation$run_until()

The method TissueSimulation$run_until() takes a formula as a parameter and runs the simulation until the formula no longer holds.

# build a condition stating that the cardinality of A+ doubles
# that of A-
c1 <- va_p >= 2 * va_m
c1
#> |A+|>=2*|A-|

# build a condition that holds when there are more than
# 100000 live cells of mutant A
c2 <- va_p + va_m > 1e5
c2
#> |A+|+|A-|>100000

# build a condition that holds when less than 4000
# epigenetic switches from the species A+ have occurred
c3 <- va_ps < 4000
c3
#> 4000>|A+.switches|

# build a condition that holds when 40 time unit have been
# simulated at least
c4 <- v_time >= 40
c4
#> Time>=40

# build a condition that holds when c4 and at least one
# among c1, c2, and c3 hold
c5 <- c4 & (c1 | c2 | c3)
c5
#> Time>=40 and (|A+|>=2*|A-| or |A+|+|A-|>100000 or 4000>|A+.switches|)

# place the initial cell
sim$place_cell("A+", 500, 500)

# run the simulation while c5 does not hold
sim$run_until(c5)
#>  [████████████████████████████████████████] 100% [00m:00s] Saving snapshot

sim
#> ──  ProCESS   D   S   M  /var/folders/tb/jqmdpgxs2t5129bny6pb96680000gn/T/ProCES
#> 
#> ── Species: 2, with epigenetics
#>    
#>    =======  ====  ====  ====  ======  ========
#>    species   λ      δ    ε    counts     %    
#>    =======  ====  ====  ====  ======  ========
#>         A-  0.08  0.01  0.01    32    76.19048
#>         A+  0.20  0.10  0.01    10    23.80952
#>    =======  ====  ====  ====  ======  ========
#> 
#> ── Firings: 74 total
#> 
#>  Species [A-]:  4 (deaths), 37 (duplications) and  5 (switches)
#>  Species [A+]:  8 (deaths), 16 (duplications) and  4 (switches)
#> ✖ The simulation has no samples yet!
sim$get_clock()
#> [1] 40.01942