cell {BSL}  R Documentation 
Cell biology example
Description
This example estimates the probabilities of cell motility and cell proliferation for a discretetime stochastic model of cell spreading. We provide the data and tuning parameters required to reproduce the results in An et al. (2019).
Usage
data(ma2)
cell_sim(theta, Yinit, rows, cols, sim_iters, num_obs)
cell_sum(Y, Yinit)
cell_prior(theta)
Arguments
theta 
A vector of proposed model parameters,

Yinit 
The initial matrix of cell presences of size 
rows 
The number of rows in the lattice (rows in the cell location matrix). 
cols 
The number of columns in the lattice (columns in the cell location matrix). 
sim_iters 
The number of discretisation steps to get to when an
observation is actually taken. For example, if observations are taken every
5 minutes but the discretisation level is 2.5 minutes, then

num_obs 
The total number of images taken after initialisation. 
Y 
A 
Details
Cell motility (movement) and proliferation (reproduction) cause tumors to spread and wounds to heal. If we can measure cell proliferation and cell motility under different situations, then we may be able to use this information to determine the efficacy of different medical treatments.
A common method for measuring in vitro cell movement and proliferation is
the scratch assay. Cells form a layer on an assay and, once they are
completely covering the assay, a scratch is made to separate the cells.
Images of the cells are taken until the scratch has closed up and the cells
are in contact again. Each image can be converted to a binary matrix by
forming a lattice and recording the binary matrix (of size rows
\times
cols
) of cell presences.
The model that we consider is a random walk model with parameters for the
probability of cell movement
(P_m
) and the probability
of cell proliferation
(P_p
) and it has no
tractable likelihood function. We use the vague priors
P_m \sim U(0,1)
and P_p \sim U(0,1)
.
We have a total of 145 summary statistics, which are made up of the Hamming distances between the binary matrices for each time point and the total number of cells at the final time.
Details about the types of cells that this model is suitable for and other information can be found in Price et al. (2018) and An et al. (2019). Johnston et al. (2014) use a different ABC method and different summary statistics for a similar example.
Functions

cell_sim
: The functioncell_sim(theta, Yinit, rows, cols, sim_iters, num_obs)
simulates data from the model, using C++ in the backend. 
cell_sum
: The functioncell_sum(Y,sum_options)
calculates the summary statistics for this example. 
cell_prior
: The functioncell_prior(theta)
evaluates the log prior density at the parameter value\theta
.
A simulated dataset
An example “observed” dataset and the tuning parameters relevant to that
example can be obtained using data(cell)
. This “observed” data is
a simulated dataset with P_m = 0.35
and
P_p = 0.001
. The lattice has 27 rows
and 36
cols
and there are num_obs = 144
observations after time 0
(to mimic images being taken every 5 minutes for 12 hours). The simulation
is based on there initially being 110 cells in the assay.
Further information about the specific choices of tuning parameters used in BSL and BSLasso can be found in An et al. (2019).

data
: Therows
\times
cols
\times
num_obs
array of the cell presences at times 1:144. 
sim_args
: Values ofsim_args
relevant to this example. 
sum_args
: Values ofsum_args
relevant to this example, i.e. just the value ofYinit
. 
start
: A vector of suitable initial values of the parameters for MCMC. 
cov
: The covariance matrix of a multivariate normal random walk proposal distribution used in the MCMC, in the form of a 2\times
2 matrix.
Author(s)
Ziwen An, Leah F. South and Christopher Drovandi
References
An Z, South LF, Nott DJ, Drovandi CC (2019).
“Accelerating Bayesian Synthetic Likelihood With the Graphical Lasso.”
Journal of Computational and Graphical Statistics, 28(2), 471–475.
doi: 10.1080/10618600.2018.1537928.
Johnston ST, Simpson MJ, McElwain DLS, Binder BJ, Ross JV (2014).
“Interpreting scratch assays using pair density dynamics and approximate Bayesian computation.”
Open Biology, 4(9), 140097.
doi: 10.1098/rsob.140097.
Price LF, Drovandi CC, Lee A, Nott DJ (2018).
“Bayesian Synthetic Likelihood.”
Journal of Computational and Graphical Statistics, 27, 1–11.
doi: 10.1080/10618600.2017.1302882.
Examples
## Not run:
require(doParallel) # You can use a different package to set up the parallel backend
# Loading the data for this example
data(cell)
model < newModel(fnSim = cell_sim, fnSum = cell_sum, simArgs = cell$sim_args,
sumArgs = cell$sum_args, theta0 = cell$start, fnLogPrior = cell_prior,
thetaNames = expression(P[m], P[p]))
thetaExact < c(0.35, 0.001)
# Performing BSL (reduce the number of iterations M if desired)
# Opening up the parallel pools using doParallel
cl < makeCluster(min(detectCores()  1,2))
registerDoParallel(cl)
resultCellBSL < bsl(cell$data, n = 5000, M = 10000, model = model, covRandWalk = cell$cov,
parallel = TRUE, verbose = 1L)
stopCluster(cl)
registerDoSEQ()
show(resultCellBSL)
summary(resultCellBSL)
plot(resultCellBSL, thetaTrue = thetaExact, thin = 20)
# Performing uBSL (reduce the number of iterations M if desired)
# Opening up the parallel pools using doParallel
cl < makeCluster(min(detectCores()  1,2))
registerDoParallel(cl)
resultCelluBSL < bsl(cell$data, n = 5000, M = 10000, model = model, covRandWalk = cell$cov,
method = "uBSL", parallel = TRUE, verbose = 1L)
stopCluster(cl)
registerDoSEQ()
show(resultCelluBSL)
summary(resultCelluBSL)
plot(resultCelluBSL, thetaTrue = thetaExact, thin = 20)
# Performing tuning for BSLasso
ssy < cell_sum(cell$data, cell$sum_args$Yinit)
lambda_all < list(exp(seq(0.5,2.5,length.out=20)), exp(seq(0,2,length.out=20)),
exp(seq(1,1,length.out=20)), exp(seq(1,1,length.out=20)))
# Opening up the parallel pools using doParallel
cl < makeCluster(min(detectCores()  1,2))
registerDoParallel(cl)
set.seed(100)
sp_cell < selectPenalty(ssy, n = c(500, 1000, 1500, 2000), lambda_all, theta = thetaExact,
M = 100, sigma = 1.5, model = model, method = "BSL", shrinkage = "glasso",
parallelSim = TRUE, parallelMain = FALSE)
stopCluster(cl)
registerDoSEQ()
sp_cell
plot(sp_cell)
# Performing BSLasso with a fixed penalty (reduce the number of iterations M if desired)
# Opening up the parallel pools using doParallel
cl < makeCluster(min(detectCores()  1,2))
registerDoParallel(cl)
resultCellBSLasso < bsl(cell$data, n = 1500, M = 10000, model = model, covRandWalk = cell$cov,
shrinkage = "glasso", penalty = 1.3, parallel = TRUE, verbose = 1L)
stopCluster(cl)
registerDoSEQ()
show(resultCellBSLasso)
summary(resultCellBSLasso)
plot(resultCellBSLasso, thetaTrue = thetaExact, thin = 20)
# Performing semiBSL (reduce the number of iterations M if desired)
# Opening up the parallel pools using doParallel
cl < makeCluster(min(detectCores()  1,2))
registerDoParallel(cl)
resultCellSemiBSL < bsl(cell$data, n = 5000, M = 10000, model = model, covRandWalk = cell$cov,
method = "semiBSL", parallel = TRUE, verbose = 1L)
stopCluster(cl)
registerDoSEQ()
show(resultCellSemiBSL)
summary(resultCellSemiBSL)
plot(resultCellSemiBSL, thetaTrue = thetaExact, thin = 20)
# Plotting the results together for comparison
# plot using the R default plot function
oldpar < par()
par(mar = c(5, 4, 1, 2), oma = c(0, 1, 2, 0))
combinePlotsBSL(list(resultCellBSL, resultCelluBSL, resultCellBSLasso, resultCellSemiBSL),
which = 1, thetaTrue = thetaExact, thin = 20, label = c("bsl", "ubsl", "bslasso", "semiBSL"),
col = 1:4, lty = 1:4, lwd = 1)
mtext("Approximate Univariate Posteriors", outer = TRUE, cex = 1.5)
par(mar = oldpar$mar, oma = oldpar$oma)
## End(Not run)