Likelihood-based inference, identifiability and prediction using count data from lattice-based random walk models

In this study, we use two types of mathematical modelling frameworks. First, we use a stochastic lattice-based random walk model that will be described in Section 2.1.1. The motivation for using a computationally expensive stochastic model is that it provides a high-fidelity means of generating noisy data that is compatible with the kind of data generated experimentally. Second, we use a computationally efficient continuum limit description of the stochastic model. The continuum limit description takes the form of a system PDEs that will be described in Section 2.1.2.

For this simulation study we use synthetic data generated by stochastic model described in Section 2.1.1 to generate count data that has similar properties to count data obtained in a scratch assay. With this data we will explore different options for efficient likelihood-based parameter estimation, parameter identifiability analysis and model prediction. We will focus on two different cases: (i) Case 1 involves generating data at a single time point involving biased motility in the context of working with a single homogeneous population with $S=1$; and (ii) Case 2 involves generating data at a single time point involving biased motility with a population composed of two subpopulations with $S=2$, All data are denoted using the vector $y^{\text{o}}$. For Case 1, $y^{\text{o}}$ is a vector of length $I$; for Case 2, $y^{\text{o}}$ is a vector of length $2I$.

Case 1 involves a biased single population, $S=1$, on a $200\times 20$ lattice where the motion of agents is biased, $(P_{1},\rho_{1})^{\top}=(1,0.1)^{\top}$ which corresponds to $(D_{1},v_{1})^{\top}=(0.25,0.05)^{\top}$ in the continuum model under idealised conditions where $J$ is sufficiently large. Since the motion is biased in the positive $x$-direction, the initial placement of agents involves fully occupying sites with $10\leq i\leq 40$, as shown in Figure 3(a). The placement of agents towards the left boundary of the domain allows us to observe the biased motion in the positive $x$-direction, and data collected after $k=300$ time steps gives rise to the snapshot shown in Figure 3(b). Count data in Figure 3(c)–(d) are given at the beginning and end of the simulation, respectively.

Case 2 involves a biased population composed of two subpopulations, $S=2$, on a $200\times 20$ lattice with $(P_{1},P_{2},\rho_{1},\rho_{2})^{\top}=(0.8,1,0.2,0.0)^{\top}$, corresponding to $(D_{1},D_{2},v_{1},v_{2})^{\top}=(0.2,0.25,0.08,0.0)^{\top}$ in the continuum model under idealised conditions where $J$ is sufficiently large. The initial placement of agents involves placing a group of agent from subpopulation 1 in the central region of the lattice so that all sites with $80\leq i\leq 120$ are completely occupied by agents from subpopulation 1. All remaining lattice sites are randomly occupied by agents from subpopulation 2 with probability 0.5, as shown in Figure 3(e). This means that agents from subpopulation 1 will attempt to spread out from their original location with a small bias in their motion in the positive $x$-direction. Since the motion of agents from subpopulation 1 are hindered by the presence of surrounding agents from subpopulation 2 we collect data after a larger number of time steps, $k=1000$, which gives rise to the snapshot in Figure 3(f). Count data for subpopulation 1 are given in Figure 3(g)–(h) whereas count data for subpopulation 2 are given in Figure 3(i)–(j).

Our count data, summarised in Figure 3, are deliberately generated with $J=20$ so that the count data involve clear fluctuations and so the mathematical relationship between the discrete count data and the solution of the corresponding continuum limit model, Equation (2.1.2), is not guaranteed to hold. Therefore, we introduce two measurement error models [29] that allow us to relate the solution of the continuum model to the noisy count data in a probabilistic sense. In the first instance we take a standard approach and work with an additive Gaussian measurement error model because this is the most standard approach in the biological physics and mathematical biology literature. Secondly, we introduce a multinomial measurement error model which, as we will show, is a natural choice for working with count data.

The standard approach for relating data to solutions of differential equations is to assume that the data $y^{\text{o}}$ is normally distributed about the solution of the differential equation of interest. Throughout this work we will refer to this approach as working with an additive Gaussian measurement error model. For example, if we consider Case 1 where we have $y^{\text{o}}=C_{1}^{\text{o}}(i,k)$ for $i=1,2,3,\ldots,I$ and the solution of the continuum limit mathematical model is $c_{1}(x_{i},t)$ for $i=1,2,3,\ldots,I$, we make the standard assumption that $C_{1}^{\text{o}}(i,k)=\mathcal{N}(c_{1}(x_{i},t),\sigma_{1}^{2})$ for $i=1,2,3,\ldots,I$. Since we deal with $S$ counts across $I$ columns of the lattice, invoking a standard independence assumption gives us a log-likelihood function that can be written as

where $\phi(x\mid\mu,\sigma^{2})$ is the Gaussian probability density function with mean $\mu$ and variance $\sigma^{2}$, and $\theta$ refers to the parameters in the continuum limit PDE model.

Working with an additive Gaussian measurement error model means that we have introduced additional parameters in the noise model, $\sigma_{s}^{2}$ for $s=1,2,3,\ldots,S$. Throughout this work we will present results generated using the Gaussian additive noise model in terms of the standard deviation $\sigma_{s}$ for $s=1,2,3,\ldots,S$. As we will demonstrate, data will be used to estimate the parameters in the continuum limit model description of the stochastic model as well as simultaneously estimating parameters in the measurement error model. For example, in Case 1 where we have $S=1$ we have two parameters in the continuum limit model $D_{1}$ and $v_{1}$, and one parameter in the measurement error model $\sigma_{1}$. When working with the additive Gaussian measurement error model we will include both sets of parameters in the vector $\theta$ which means that with $S=1$ we have $\theta=(D_{1},v_{1},\sigma_{1})^{\top}$ whereas with $S=2$ we have $\theta=(D_{1},D_{2},v_{1},v_{2},\sigma_{1},\sigma_{2})^{\top}$.

Instead of simply assuming that the data is normally distributed about the solution of the appropriate continuum limit model, we can take advantage of the structure of our count data and derive a more mechanistically-motivated measurement error model that we will refer to as a multinomial measurement error model. For example, when dealing with just a single population, $S=1$, our count data simply consists of counts of agents within each column of the lattice, which implicitly defines a count of the vacant lattice sites in each column, meaning that our data takes the form $C_{1}^{\text{o}}(i,k)$ and $E^{\text{o}}(i,k)=J-C_{1}^{\text{o}}(i,k)$ for $i=1,2,3,\dots,I$. If we interpret $C_{1}^{\text{o}}(i,k)$ at any given spatial location as corresponding to a finite number of independent samples from some underlying stochastic process, these samples can be thought of as an approximation of a continuous measurement of agent occupancy within that column. These data can be considered as samples from a distribution where the expected occupancy fraction is given by the solution of the continuum limit model $c_{s}(x_{i},t)$. Together, these measurements imply a binomial likelihood,

Assuming all column counts, $i=1,2,3,\ldots,I$, are conditionally independent, given the continuum-limit solution, and taking the logarithm of the likelihood function gives us a log-likelihood function based on the binomial measurement error model for a single population with $S=1$ that can be written as

Similar arguments lead to a multinomial-based log-likelihood function for $S>1$. Under these conditions our count data consists of counts of agents from each subpopulation within each column, $C_{s}^{\text{o}}(i,k)$ for $i=1,2,3,\ldots,I$ and $s=1,2,3,\ldots,S$. The geometry of the lattice means that counts of vacant sites is given by $\displaystyle{E^{\text{o}}(i,k)=J-\sum_{s=1}^{S}C_{s}^{\text{o}}(i,k)}$ for $i=1,2,3,\dots,I$. Again, we can interpret $C_{s}^{\text{o}}(i,k)$ at any given spatial location as being generated by a finite number of independent samples from some underlying stochastic process, these samples can be thought of as an approximation of expected, noise-free measure of agent occupancy from the $s$th subpopulation within that column. These data can be considered as samples from a distribution where the expected occupancy fraction is given by the solution of the continuum limit model, $c_{s}(x_{i},t)$. Following the same arguments as above for $S=1$, these measurements imply multinomial log-likelihood function which can be written as

This multinomial log-likelihood function relaxes to the binomial log-likelihood function, derived previously, when $S=1$, and here we have written the data using the compact vector notation $y^{\text{o}}$. In practice, when we evaluate the multinomial log-likelihood function we simply set the proportionality constant to unity. Unlike the additive Gaussian measurement error model, working with the multinomial measurement error model does not involve introducing any new parameters. When we work with $S=1$ we have $\theta=(D_{1},v_{1})^{\top}$, whereas $S=2$ involves $\theta=(D_{1},D_{2},v_{1},v_{2})^{\top}$.

Given a set of count data $y^{\text{o}}$ together with an appropriate process model and measurement error model with all unknown parameters summarised in the vector $\theta$, we have access to a log-likelihood function, $\ell(\theta\,|y^{\text{o}})$. Consequently, the choice of parameters that gives the best match to the data is given by $\hat{\theta}$ which maximises the log-likelihood, giving rise to the maximum likelihood estimate (MLE),

Throughout this work we always use numerical optimization calculate $\hat{\theta}$. In particular all numerical optimization calculations are performed using the Nelder-Mead algorithm with simple bound constraints [32] implemented within the NLopt routine [33] in Julia. In general we find that our numerical optimization calculations are insensitive to the initial estimate and that the algorithm converges using default stopping criteria.

Given our estimate of $\hat{\theta}$, we use then the profile likelihood to quantify the precision of our parameter estimates by examining the curvature of the log-likelihood function [24, 25, 26, 27, 28]. To compare our results against asymptotic thresholds we work with a normalised log-likelihood function

so that $\bar{\ell}(\hat{\theta}\mid y^{\text{o}})=0$. To proceed we partition the full parameter $\theta$ into interest parameters $\psi$, and nuisance parameters $\omega$, so that $\theta=(\psi,\omega)$. In this study we restrict our attention to univariate profile likelihood functions which means that our interest parameter is always a single parameter. For a set of data $y^{\text{o}}$, the profile log-likelihood for the interest parameter $\psi$ given the partition $(\psi,\omega)$ is

which implicitly defines a function $\omega^{*}(\psi)$ of optimal values of $\omega$ for each value of $\psi$. As for calculating the MLE, all profile likelihood functions in this work are calculated using the Nelder-Mead numerical optimization algorithm with the same bound constraints used to calculate $\hat{\theta}$. As a concrete example, if we work with the continuum limit model for a single population using the additive Gaussian measurement error model we have $\theta=(D_{1},v_{1},\sigma_{1})^{\top}$. In this scenario we can compute three univariate profile likelihood functions by choosing the interest parameter to be: (i) the diffusivity, $\psi=D_{1}$ and $\omega=(v_{1},\sigma_{1})^{\top}$; (i) the drift velocity, $\psi=v_{1}$ and $\omega=(D_{1},\sigma_{1})^{\top}$; and (iii) the standard deviation in the measurement error model, $\psi=\sigma_{1}$ and $\omega=(D_{1},v_{1})^{\top}$. For all univariate profile likelihood calculations we work with a uniformly-discretized interval that contains the MLE. For example, if our interest parameter is $D_{1}$ we identify an interval, $D_{1}^{-}<\hat{D}_{1}<D_{1}^{+}$ and evaluate $\bar{\ell}_{p}$ across a uniform discretization of the interval to give a simple univariate function that we call the profile likelihood. The degree of curvature of the profile likelihood function provides information about the practical identifiability of the interest parameter. For example, if the profile log-likelihood function is flat then the interest parameter is non-identifiable. In contrast, when the profile log-likelihood function is curved the degree of curvature indicates inferential precision and we may determine likelihood-based confidence intervals where $\bar{\ell}_{p}<-\Delta_{q,n}/2$ where $\Delta_{q,n}$ denotes the $q$th quantile of the $\chi^{2}$ distribution with $n$ degrees of freedom, which we take to be the relevant number of unknown parameters [34]. For example, identifying the interval where $\bar{\ell}_{p}<-\Delta_{0.95,1}/2\approx-1.9207$ allows us to identify the asymptotic 95% confidence interval [24, 29] for a univariate profile likelihood function with one free parameter. This procedure gives us a simple way of identifying the width of the interval $[D_{1}^{-},D_{1}^{+}]$. A simple approach to determine suitable choices of $D_{1}^{-}$ and $D_{1}^{+}$ is to compute $\bar{\ell}_{p}$ across some initial estimate of the interval and if the values of the profile log-likelihood do not intersect the relevant asymptotic threshold then we simply continue to compute the profile log-likelihood function on a wider interval. In contrast, if we compute the profile log-likelihood on some interval and find that this interval is very wide compared to the relevant asymptotic threshold we can simply re-compute the profile log-likelihood function on a narrower interval. Since these computations involve numerical optimisation calculations where we always have a reasonably good estimate to start the iterative calculations (e.g. $\hat{\theta}$) we find that these computations are efficient.

The methods outlined so far focus on taking noisy count data and exploring how to find the MLE parameter estimates, and to take an asymptotic log-likelihood threshold to define a parameter confidence set so that we can understand how variability in count data corresponds to variability in parameter estimates. In this section we will now explore how variability in parameter estimates maps to variability in model predictions. In particular, our focus will be to compare likelihood-based predictions using the traditional additive noise model with the multinomial noise model and explore differences in these approaches.

Given a set of count data, $y^{\textrm{o}}$, and an associated normalised log-likelihood function, $\bar{\ell}(\theta\mid y^{\textrm{o}})$, we proceed by identifying the asymptotic 95% log-likelihood threshold $-\Delta_{0.95,n}/2$, where $n$ is the number of free parameters. Using rejection sampling we obtain $M$ parameter sets that lie within the 95% log-likelihood threshold so that we have $\theta_{m}$ that satisfy $\bar{\ell}(\theta_{m}\,|\,y^{\text{o}})\geq-\Delta_{0.95,n}/2$ for $m=1,2,3,\ldots,M$. For each $\theta_{m}$ within the confidence set we solve the relevant continuum model to give $M$ solutions that correspond to the mean trajectory prediction as defined by the relevant measurement error model [29]. Since all measurement error models considered in this work involve a parametric distribution, we can also construct various data realization predictions [29] by considering various measures of the width of that parameteric distribution about the mean. For simplicity we will consider the 5% and 95% quantiles of the relevant distributions as a simple measure of distribution width. To make a prediction using the single population model we denote the mean trajectory as $c_{1}^{(m)}(x,t)$ to denote the solution of Equation (5) using the $m$th parameter sample from within the parameter confidence set. To plot the solutions we discretize the spatial location so that we have $c_{1}^{(m)}(x_{r},t)$, where $x_{r}=(r-1)/10$ for $r=1,2,3,\ldots,2001$. For each of the mean trajectories we now have an interval $c^{(m)}_{1,0.05}(x_{r},t)<c_{1}^{(m)}(x_{r},t)<c^{(m)}_{1,0.95}(x_{r},t)$ for $r=1,2,3,\ldots,2001$, where the lower bound corresponds to the 5% quantile of the measurement error noise model and the upper bound corresponds to the 95% quantile of the measurement error noise model. For a fixed parameter $\theta_{m}$, the 5% and 95% quantiles of the additive Gaussian measurement error model are constants determined by $\sigma_{1}$. For the multinomial measurement error model the 5% and 95% quantiles are no longer constants but vary with the mean of the distribution. For $S=1$ the multinomial distribution relaxes to the binomial distribution and hence the bounds at location $x_{r}$ correspond to the 5% and 95% quantiles of the binomial distribution with $J$ trials with probability of success $c_{1}(x_{r},t)$. For the multinomial measurement error model with $S>1$ the bounds for subpopulation $s$ at location $x_{r}$ correspond to the 5% and 95% quantiles of the $s$th marginal distribution, again with $J$ trials with probability of success $c_{s}(x_{r},t)$. After calculating these $M$ intervals around each mean trajectory we then take the union to define $\left[\text{min}\left(c^{(m)}_{1,0.05}(x_{r},t)\right),\text{max}\left(c^{(m)}_{1,0.95}(x_{r},t)\right)\right]$. Here the minimum value of $c^{(m)}_{1,0.05}(x_{r},t)$ and the maximum value of $c^{(m)}_{1,0.95}(x_{r},t)$ are computed for each fixed value of $x_{r}$ across the set of $M$ different parameter values, $m=1,2,3,\ldots,M$. This simple computational procedure gives us a complete picture of the prediction interval that can be interpreted as a form of tolerance interval [35] accounting for both parameter uncertainty and data realization uncertainty [21].