Estimating parameters of a stochastic cell invasion model with fluorescent cell cycle labelling using approximate Bayesian computation

3.1. Simulation model

We adopt the discrete random walk model developed by Simpson et al. [7] on a 2D hexagonal lattice. Each lattice site has diameter Δ = 20 μm, which is the average cell diameter [21], and is associated with a set of unique Cartesian coordinates,

where i and j are the respective row and column indices. To mimic scratch assay experiments, cells in G1 phase are represented by red agents, cells in eS phase are represented by yellow agents, and cells in S/G2/M phase are represented by green agents. Agents are permitted to transition through phases of the cell cycle and undergo a nearest neighbour random walk by simulating from a Markov process using the Gillespie algorithm [22] where the time between events is exponentially distributed. The algorithm is presented in §1 of the electronic supplementary material.

To simulate cell migration, agents undergo a nearest neighbour random walk at rates M_r, M_y, M_g per hour for red, yellow and green agents, respectively (figure 2a–f). Potential movement events involve randomly selecting the target site from the set of six nearest-neighbouring lattice sites, with the movement event being successful only if the target site is vacant. In this way, crowding effects are simply accommodated. To simulate transitions through the cell cycle, red agents are allowed to transition into yellow agents at rate R_r per hour (figure 2h,i), yellow agents to green agents at rate R_y per hour (figure 2i,j) and green agents into two red daughter agents at rate R_g per hour (figure 2j,k). While we assume that the red-to-yellow and yellow-to-green transitions are unaffected by crowding, we model crowding effects for the green-to-red transition by aborting transitions where the additional red daughter agent would be placed onto an occupied lattice site. By prohibiting multiple agents from occupying the same lattice site, we are able to realistically incorporate crowding effects [23,24].

• Download figure
• Open in new tab
• Download PowerPoint

The Simpson et al. [7] model is dependent on the initial geometry, boundary conditions, the lattice spacing Δ, and the cell cycle transition and motility rates. Since we have reasonable estimates for Δ [21] and we calibrate the initial geometry and boundary conditions to the experimental data, our study is concerned with estimating the unknown cell cycle transition and motility parameters. In a Bayesian setting, the unknown model parameters, θ = (R_r, R_y, R_g, M_r, M_y, M_g), and the uncertainty around them can be quantified by the posterior distribution, which is dependent on the likelihood and the prior distribution. However, while the Markov process model can capture the stochastic nature of cell proliferation and migration, when the dimension of the generator matrix (a matrix of rate parameters which describe the rate of transitioning between states) is too high the likelihood function consequently becomes intractable due to the computational cost of computing the matrix exponential (see [25–28]). Since conventional Bayesian approaches to parameter estimation are no longer feasible, we are motivated to use likelihood-free methods.

3.2. Approximate Bayesian computation

Using a Bayesian framework, the uncertainty about the unknown parameter θ with respect to the data y can be quantified by sampling from the posterior distribution π(θ|y) ∝ π(y|θ)π(θ), where π(y|θ) is the likelihood function and π(θ) is the prior. However, the likelihood function for sufficiently complex models becomes intractable (see examples in biology [12,24,29], in ecology [30,31] and in cosmology [32]). Rather than reverting to simpler models with tractable likelihoods, these types of problems can be instead analysed using likelihood-free methods that avoid evaluating the likelihood function.

One popular likelihood-free approach is ABC [33]. ABC involves simulating data from the model x ∼ π(· |θ) instead of evaluating the intractable likelihood; accepting configurations of θ which produce simulated data x that is close to the observed data y. It can be impractical to compare the full data sets of x and y, so ABC often relies on reducing the full data sets to summary statistics by some summarizing function S(·), where the summary statistics for x and y are denoted S_x = S(x) and S_y = S(y), respectively. Provided the summary statistics are highly informative about the model parameters, then π(θ|y) ≈ π(θ|S_y) is a good approximation or exact if sufficient statistics are used [34]. However, the latter are usually difficult to attain in practice and so this study focuses on the use of summary statistics. In effect, ABC samples from the approximate posterior:

where ρ(S_y, S_x) is the discrepancy function which measures the difference between the two data sets and $K_{ϵ}(\cdot )$ is the kernel weighting function which weighs ρ(S_y, S_x) conditional on the tolerance ε. A common choice for the discrepancy function is the Euclidean distance, $\rho (S_y,S_x)={\parallel S_y-S_x\parallel }_2$, and for the kernel weighting function is the indicator function, $\mathbb{1}(\cdot )$, which is equal to one if ρ(S_y, S_x) ≤ ε and is zero otherwise. The approximate posterior in equation (3.2) converges to the posterior conditional on the observed summary (often referred to as the partial posterior) in the limit as ε → 0 [35].

To sample from the approximate posterior, commonly ABC-rejection [36,37], MCMC-ABC [38], or sequential Monte Carlo ABC (SMC-ABC) (e.g. [39–41]) algorithms are used. ABC-rejection samples particles from the prior distribution and accepts particles with a discrepancy measure ρ(S_y, S_x) less than the desired tolerance ε. In cases when the prior distribution is relatively diffuse compared to the posterior density (such as our application), lower acceptance rates are common because particles are predominantly sampled in regions of low posterior density [39]. To increase efficiency, one could instead use MCMC-ABC which constructs a Markov chain with a stationary distribution identical to the approximate posterior by proposing particles from a carefully-tuned proposal distribution, θⁱ ∼ q(· |θⁱ⁻¹), and accepting those with probability

which is based on the Metropolis–Hastings ratio [42,43]. While MCMC-ABC tends to be more computationally efficient compared to ABC-rejection [38], it is possible for the Markov chain to spend many iterations in areas of low posterior probability. In our application, we found MCMC-ABC to take a considerable effort to tune the proposal distribution while still being computationally cumbersome. However, SMC-ABC or more specifically the SMC-ABC replenishment algorithm [40] requires very little tuning comparatively and allows for simulations to be performed in parallel to increase computational efficiency.

The SMC-ABC replenishment algorithm traverses a set of distributions defined by T non-increasing tolerance levels ε₁ ≥ · · · ≥ ε_T to sample from the approximate posterior:

where the first target distribution is constructed by sampling from the prior distribution to attain a collection of parameter values (called particles) and their discrepancies, ${\{{\theta }^i,{\rho }^i\}}_{i=1}^N$. The first tolerance threshold, ε₁, is set as the maximum of the set of discrepancies. Thereafter, to propagate particles through the sequence of target distributions, particles are first sorted in ascending order by their discrepancy and the new tolerance is set as ${ϵ}_t={\rho }^{N-N_{\alpha }}$ where $N_{\alpha }=\lfloor \alpha N\rfloor$, α is the proportion of particles discarded and $\lfloor \cdot \rfloor$ is the floor function. Particles, ${\{{\theta }^i\}}_{i=N-N_{\alpha }+1}^N$, which do not satisfy the new tolerance are discarded and resampled, with replacement from the remaining particles to replenish the population. To prevent sample degeneracy (too many duplicated particles), resampled particles are then perturbed according to an MCMC kernel R_t times with an invariant distribution given by the current approximate posterior ${\pi }_{{ϵ}_t}(\theta |S_y)$. In each of the R_t iterations, the proposed particles are drawn from an automatically tuned proposal distribution θⁱ ∼ q( · |θⁱ⁻¹) and accepted with probability $p_{\mathrm{acc}}^{i,j}$ (equation (3.3)) where i denotes the ith particle and j the jth MCMC iteration. To ensure sample diversity, R_t can be dynamically set based on the overall MCMC acceptance rate, $\overline{\hspace{0.17em}{p}_{\mathrm{acc}}}=1/(N_{\alpha }\times R_t)\sum _{i=1}^{N_{\alpha }}$ $\sum _{\hspace{0.17em}j=1}^{R_t}{p}_{\mathrm{acc}}^{i,j}$, such that there is a 1 − c chance that all particles are moved at least once and is given by where the ceiling function $\lceil \cdot \rceil$ is used to be conservative and an estimate for $\overline{\hspace{0.17em}{p}_{\mathrm{acc}}}$ is calculated from R_t−1/2 pilot MCMC iterations. A popular choice for the proposal distribution is the multivariate normal distribution, $q({\theta }^i|{\theta }^{i-1})=\mathcal{N}({\theta }^i;{\theta }^{i-1},\Sigma )$, where the covariance matrix $\Sigma$ is the tuning parameter. To create a more efficient proposal distribution, we can adaptively tune $\Sigma$ by computing the empirical covariance matrix of the ${\{\theta \}}_{i=1}^{N-N_{\alpha }}$ particles which are already distributed according to the current target distribution. In this application, we do not worry about scaling the covariance matrix as there are only six parameters. The algorithm finally stops once the overall MCMC acceptance rate, $\overline{\hspace{0.17em}{p}_{\mathrm{acc}}}$, is unreasonably low (less than or equal to 1%) or the desired tolerance threshold is reached. For the two tuning parameters, Drovandi & Pettitt [40] suggest setting α = 0.5 and c = 0.01. The SMC-ABC replenishment algorithm is presented in algorithm 1 and is hereafter referred to as SMC-ABC.

A crucial limitation of ABC methods is the curse of dimensionality, where despite the addition of more data, the approximation to the posterior can become distorted as a result of the discrepancy between observed data and simulated data ρ(S_y, S_x) naturally increasing with the dimension [35]. In applications where increasing the dimension of the summary statistics cannot be avoided, the discrepancy between observed and summary statistics can be accounted for, at least approximately, with regression adjustment [34,35]. Regression adjustment involves explicitly modelling the parameters against the discrepancy between observed and simulated data. Assume for the moment that θ is a scalar parameter. Consider the following regression model

where i = 1, …, N is the parameter sample index, β is the regression coefficients, β₀ is the intercept and ɛⁱ is the error term. Estimates for β₀ and β can be computed by minimizing the weighted least squares criterion $\sum _{i=1}^N{w}^i{({\theta }^i-{\beta }_0-{(S_x^i-S_y)}^{\mathrm{\top }}\beta )}^2$. Here, we choose to use the popular Epanechnikov weighting function [44], defined as $w^i=0.75(1-{({\rho }^i/max({\{{\rho }^i\}}_{i=1}^N))}^2)$, but other weighting functions could also be used. Using the estimated regression coefficients $\hat{\beta }$, we then make the adjustment The adjusted sample ${\{{\theta }^{i^{\ast }}\}}_{i=1}^N$ can often give a more accurate approximation of the posterior. To ensure that the adjusted parameters remain within the support of the prior distribution (if bounded), Hamilton et al. [45] suggest transforming parameter values before applying the regression adjustment. We use a logit transformation, $\tilde{\theta }=\log ((\theta -a)/(b-\theta ))$, where a and b are the respective lower and upper bounds of the prior. Given that we have a vector of parameters, we apply a regression adjustment to each component of the parameter vector separately.

3.3. Prior knowledge

The cell cycle time, which is related to the doubling time, can be thought of as the summation of the time spent in each phase of the cell cycle. Estimates of the cell doubling time for melanoma cells range from 16 to 47 h [4,9]. Furthermore, Simpson et al. [9] estimate the average time 1205Lu FUCCI-transduced melanoma cells spend in the G1 to be between 8 and 30 h. This means that the transition from G1 to S/G2/M phase is approximately 1/30−1/8 h⁻¹. Additionally, the duration spent in S/G2/M phases were reported as 8–17 h. This means that the transition from S/G2/M to G1 phase is approximately 1/17–1/8 h⁻¹. Therefore, we propose that our prior information of the transition rates is uniform over the range 0–1 h⁻¹ to be conservative.

Cell diffusivity, D, the measurement of motility rate for particles undergoing random diffusive migration, can be used to quantify the cell motility rate $M\hspace{0.17em}\in \hspace{0.17em}\{M_r,M_y,M_g\}$, by D = MΔ²/4 [46], where Δ is the cell diameter. Empirical evidence finds estimates for cell diffusivity to range from 0 to 3304 μm² h⁻¹ [13,14,21]. Furthermore, Simpson et al. [7] suggest the cell diffusivity to be approximately 400 μm² h⁻¹, so that the rates are approximately 4 h⁻¹. We propose that the prior information of the motility rates to be uniform over the range 0−10 h⁻¹; attributing the larger interval to the greater variation of cell diffusivity estimates in existing literature.

4.1. Developing summary statistics and validation with synthetic data

The accuracy and precision of ABC methods in approximating the posterior distribution is sensitive to the quality of the summary statistics used [35]. We first trial and validate different summary statistics with multiple synthetic datasets such that the true parameter values are known. In this way, we are able to compare the performance of different summary statistics and determine which are the most effective. While trying to replicate the environment of the experimental data as close as possible, such as domain size, boundary conditions and initial number of cells, we do not calibrate the initial location of cells but rather randomly distribute the cells within a 200 μm by 1745.35 μm region on either side of the scratch. We attain the initial cell counts of red, yellow and green cells by using the procedure outlined in §4.2 (steps 1–4) and report them here to be 119, 35 and 121, respectively.

In our analysis of the simulation model, agents with relatively higher transition rates were found to correspond to lower population sizes, and vice versa. Therefore, we use the number of agents in each population (N_r, N_y, N_g) at the end of the experiment as summary statistics that may be informative about the transition rates R_r, R_y, R_g, respectively. We test the suitability of this summary statistic on four synthetic data sets produced by varying the transition rates amongst biologically plausible values and keeping the motility rates known and constant. Estimates in existing literature of cell transition rates are similar [4,9] and the efficiency of the simulation model is dependent on the number of agents in the system (higher transition rates increase the overall proliferation rate and the frequency of events). Therefore, we choose to keep the transition rates rather close to the estimates of Haass et al. [4], instead of varying them over the extents of the prior domain. Thus, the four parameter configurations we choose to generate the synthetic datasets are $\theta \hspace{0.17em}\in \{(0.04,0.17,0.08,4,4,4),(0.25,0.15,0.22,4,4,4),(0.12,0.07,0.03,4,4,4),(0.3,0.36,0.28,4,4,4)\}$. In §2 of the electronic supplementary material, we present the marginal posterior distributions produced and confirm the suitability of this summary statistic.

For the motility parameters, we explore and compare two sets of summary statistics, namely cell density and cell trajectory data. Of these two datasets, cell density data are desirable due to less manual effort needed to generate the data while cell trajectory data could offer more information but is more challenging to collect. For the cell density data, we first segment the imaged region at the end of the experiment (t = 48 h) directly down the centre of the image in the y direction and calculate the median position and interquartile ranges of the red, yellow and green agent populations in the x-direction for cells on the left and right sides. For the cell trajectory data, we average the distance of multiple cell trajectories through each cell phase until the cell returns to the initial phase or the simulation is terminated. We select cells to be tracked provided that the cell is initially in G1 (red) phase and the cell is located on the leading edge of the cell monolayer toward the gap in the scratch assay. The reasoning behind beginning tracking from the G1 phases is due to a short period of fluorescent negativity in between S/G2/M phases and G1 phase [4] which makes tracking between these phases difficult. Therefore, by starting from the G1 phase, we avoid having to track between these two phases. Furthermore, cells were identified as being on the leading edge if their path towards the middle of the scratch was unhindered. The process of manually tracking cell trajectories can be time consuming. Thus, we are interested in finding the minimum number of cells to track such that sufficient information is acquired. In §2 of the electronic supplementary material, we draw samples from the posterior distribution using 10, 20, 30, 40 and 50 cell trajectories with four synthetic data sets generated from $\theta \hspace{0.17em}\in \{(0.04,0.17,0.08,4,4,4),(0.04,0.17,0.08,2,5,8),\hspace{0.17em}(0.04,0.17,\hspace{0.17em}0.08,8,2,5),(0.04,0.17,0.08,8,2,5),(0.04,0.17,0.08,5,8,2)\}$. Our analysis finds diminishing returns of parameter precision as the number of cell trajectories increases. We find that using 20 cell trajectories achieves a good balance between precision and the number of cell trajectories used. Using the same four synthetic datasets, we also attempt to draw samples from the posterior distribution using cell density data with the cell transition rates held constant in §2 of the electronic supplementary material. However, under these settings, we found the motility parameters to be non-identifiable.

We now combine the summary statistics formulated to estimate the cell cycle transition and motility rates together with four synthetic data sets. Due to the similarity in estimates for cell cycle transition rates in existing literature (see [4,9]), we adopt estimates for the cell cycle transition rates from Haass et al. [4] for all four parameter configurations. Since estimates for motility rates have been reported to vary by two orders of magnitude (see [13,14,21]), we choose to vary the motility rates over the range of the prior for the four parameter configurations. That is, we generate four synthetic data sets with θ ∈ {(0.04, 0.17, 0.08, 4, 4, 4), (0.04, 0.17, 0.08, 2, 5, 8), (0.04, 0.17, 0.08, 8, 2, 5), (0.04, 0.17, 0.08, 5, 8, 2)}. In figure 3, we present the marginal posterior distributions produced when using the number of cells in each subpopulation and cell density data as summary statistics. In figure 4, we present the marginal posterior distributions when using cell trajectory data in place of cell density data. Again, we see that the motility parameters are practically non-identifiable when cell density data are used while both cell cycle transition and motility parameters are practically identifiable when cell trajectory data are included. Furthermore, it is clear from the concentration of the marginal posterior distributions around the true parameter values (dashed line) in figure 4 that the cell count and cell trajectory data are highly informative about the transition and motility parameters, respectively. We note that the precision of these distributions is greater for the cell cycle transition parameters than the motility parameters. Importantly, these results show for the first time that practical parameter inference on both transition and motility parameters of a FUCCI scratch assay experiment using Bayesian inference techniques is possible. These results justify the choice of the Markov process model compared to simpler continuum models which do not give insight into cell trajectory data (see [9]).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

4.2. Image analysis of experimental data

We analyse the experimental images using ImageJ [47] to record the Cartesian coordinates of cells. Of primary interest is processing the initial frame such that we can replicate the experimental settings as accurately as possible in the simulation but we also repeat this procedure for the final frame to retrieve the final cell counts and cell density data, which we use as summary statistics. The process is as follows:

• Step 1: Read in image: File > Open > select image (figure 5a).

• Step 2: Convert image to 8-bit: Image > Type > 8 − bit (figure 5b).

• Step 3: Identify cell edges: Convert the image to black and white ($\mathrm{Process}>\mathrm{Binary}>\mathrm{convert}\hspace{0.17em}\mathrm{to}\hspace{0.17em}\mathrm{mask}$) and then distinguish conjoined cells (Process > Binary > Watershed) (figure 5c).

• Step 4: Compute Cartesian coordinates: Analyse > Analyse particles…> OK.

• Download figure
• Open in new tab
• Download PowerPoint

A limitation of using the watershed tool is that we must convert the image to black and white. In doing so, we lose the cell phase identity associated with the cell coordinates recovered from ImageJ. To overcome this, we use R [48] to retrieve the RGB decimal colour code and Cartesian coordinates of pixels. Matching pixel coordinates recovered from R and the coordinates of the centroid of the cells recovered in ImageJ, we create a data set of cell coordinates and their associated RGB decimal codes. To classify the RGB coordinates into one of the three cell cycle phases we use the conditions outlined in table 1.

To extract summary statistics from the experimental data, we repeat the image processing procedure previously outlined above with the final frame (t = 48 h) and extract the final cell counts and cell density data. Additionally, we extract cell trajectory data by processing the entire sequence of still images in ImageJ with the ‘Multi-point’ tool to manually track cell coordinates between frames. We use a similar process as before to identify cell phases in these summary statistics using R and present them in table 2 and the cell trajectory data in figure 6.

• Download figure
• Open in new tab
• Download PowerPoint

Finally, we calibrate the hexagonal lattice used in the simulation model with the dataset of Cartesian coordinates recovered previously by rearranging equation (3.1) to find their associated lattice row and column indices denoted

where $\lfloor \cdot \rceil$ rounds to the nearest integer. We treat the rare instances ($<1\text{\%}$) where multiple coordinates are mapped to the same lattice space as duplicated values and omit them rather than place them on the next closest lattice site. The result from this translation of data is presented in figure 5d. We repeat this process for the initial frame of the cell trajectory data to identify the starting position. However, due to manually tracking cell trajectories, often the coordinate retrieved was not centred on the cell which in some cases caused the starting position to be mapped to an unoccupied lattice site. We intervene prior to transforming the starting position and adjust the coordinate values to the closest occupied lattice site which is chosen such that the radial distance between the coordinate and the lattice site is minimized.

4.3. Estimating model parameters with experimental data

After calibrating the simulation to the experimental data of WM983C FUCCI-transduced melanoma cells, we first attempt to sample from the posterior distribution using the number of cells in each subpopulation and cell density data (summary statistics S₁ to S₃ and S₇ to S₁₂ in table 2, respectively) and present the samples from the posterior distribution in figure 7. Consistent with results found in §4.1 and Simpson et al. [9], estimates for the motility rates are practically non-identifiable when cell density data are used.

• Download figure
• Open in new tab
• Download PowerPoint

Next, we attempt to sample from the posterior distribution using the number of cells in each subpopulation and cell trajectory data (summary statistics S₁ to S₃ and S₄ to S₆ in table 2, respectively). We present the marginal posterior distributions produced in figures 8a,b along with the mean, standard deviation, (2.5%, 50%, 97.5%) quantiles, and the coefficient of variation (CV) in table 3. The practical identifiability in the transition and motility parameters clearly shows the benefits of using cell tracking data as the distributions are unimodal and concentrated. We estimate the cell cycle transition rates to be between 0.0411−0.193 h⁻¹ which is consistent with estimates in existing literature [4,9]. Our estimates for cell motility were found to range between 0.316−1.12 h⁻¹ which corresponds to estimates of cell diffusivity between 31.6 and 112 μm² h⁻¹ which is reasonable considering the degree of uncertainty in existing estimates which can vary between 0 and 3304 μm² h⁻¹ [13,14,21]. The precision in parameter estimates can be quantified by the CV which is a standard measure for the dispersion of data around the mean. Using the CV, the dispersion in the transition rates range from 2.65 to 5.31% and the motility rates range from 10.9 to 18.4%. To validate the parameter estimates recovered, we also present the posterior predictive distributions for the summary statistics retained from each parameter value in the posterior in figures 8c,d. These distributions are formed by plotting the distribution of simulated summary statistics produced from the posterior samples and is compared to the observed summary statistics (dashed line). These results suggest that the Markov process model developed by Simpson et al. [7] is promising as it is able to recover the observed summary statistics of the experimental data. However, further model validation should be considered in future research to determine if the simulated cell trajectories produce similar paths to those which were observed.

• Download figure
• Open in new tab
• Download PowerPoint