1 Introduction

Longitudinal image data based on fluorescent proteins play a crucial role for both in vivo and in vitro analysis of various biological processes such as gene expression and cell lineage fate. Assessing the growth patterns of different cell types within a heterogeneous population and monitoring their interactions enables biomedical researchers to determine the role of different cell types in important biological processes such as organ development and regeneration, malignant growth or immune responses under various experimental conditions. For example, tumor progression has been shown to be affected by bidirectional interactions among cancer cells or between cancer cells and cells from the microenvironment, including tumor-infiltrating immune cells [1]. Being able to study these interactions in a laboratory setting is therefore highly relevant, but is complicated by the difficulty of dissecting the effect of the different cell types as soon as the number of cell types exceeds two. In the present study we used longitudinal image data collected from multicolor live-cell imaging growth experiments of co-cultures of cancer cells and fibroblasts (a key cell type in the tumor microenvironment) as well as behaviourally distinct (cloned) cancer cells. Using a high-content imaging system, we were able to acquire characteristics for each individual cell at subsequent times, including fluorescent properties, spatial coordinates, and morphological features. The motivation of this work was to design a model allowing the determination of spatio-temporal growth interactions between these multiple cell populations.

In longitudinal growth experiments, the two important goals are to determine growth rates for different cell populations and to assess how interactions between cell types may affect their growth. Whilst a wide range of descriptive data analysis approaches have been used in applications, inference based on a comprehensive model of multicolor cell data is an open research area. The main challenges are related to the presence of complicated spatio-temporal interactions amongst cells and difficulties related to tracking individual cells across time from image data. Typical longitudinal experiments consist of a relatively small number of measurements (e.g. 5 to 20 images taken every few hours), which is adequate for monitoring cell growth. Tracking individual cells would typically require more frequent measurements, complicating the practicality of the experiments in terms of the storage cost of very large image files and the cytotoxicity induced by the imaging process.

Although tracking individual cell trajectories is difficult due to cell migration, overlapping cells, changes in cell morphology, image artifacts, cell death and division, obtaining cell counts by cell type (represented by a certain color) is straightforward and can be easily automated. To describe the spatial distribution for different cell types, we propose to divide an image into a number of contiguous regions (tiles) to form a regular lattice structure as shown in Figure 1(a). We then record the frequency of cells of different colors in each tile at subsequent time points, and based on which we model the spatial and temporal dependencies of the cell growth.

Figure 1:

(a) Microscope images for the cancer cell growth data obtained from a high-content imager (Operetta, Perkin Elmer) at the initial and final time points of the experiment. In each image, colors for non-fluorescent fibroblasts, as well as red and green fluorescent cancer cells are merged. (b) Illustration of the local structure for the model in (1). The two planes correspond to 3×3 tiles at times t and t+1. The average number of cells of color c in a given tile at time t+1 is assumed to depend on the number of cells of other colors in contiguous neighboring tiles at time t.

To model spatio-temporal data, one could choose to approximate the spatio-temporal process by a spatial process of time series, that is, to view the process as a multivariate spatial process where the multivariate dependencies are inherited from temporal dependencies. In other words, it can be seen as a temporal extension of spatial processes.

The most popular way of developing a spatial process is through the conditionally auto-regressive (CAR) model proposed by Besag [2]. Waller, Carlin, Xia, and Gelfand [3] extend the CAR model into a spatio-temporal setting by allowing spatial effects to vary across time. However, the model lacks a specification of temporal dependency, as also noted by Knorr-Held [4]. More recently, Quick, Waller, and Casper [5] proposed a multivariate space-time CAR (MSTCAR) model, which is essentially a multivariate CAR model, where both temporal and between group dependencies are modelled as multivariate dependencies. Other works related to spatial process of time series include Sans, Schmidt, Nobre, et al. [6] and Quick, Waller, and Casper [7].

Alternatively, one also think of the process as a time series of spatial process, or a spatial extension of time series. This is the approach we take in our spatio-temporal modelling. The underlying notion is that “the temporal dependence is more natural to model than the spatial dependence” [8].

Following Cox et al. [9], it is useful to distinguish two modelling approaches for the analysis of time series data commonly seen in spatial-temporal modelling literature: the parameter-driven and observation-driven model. In a parameter-driven model, the dependence between subsequent observations is modelled by a latent stochastic process, which evolves independently of the past history of the observation process. In contrast, in an observation-driven model, time dependence arises because the conditional expectation of the outcome given the past depends explicitly on the past values.

For multivariate count data, the advantage of parameter-driven models is that one can easily assume that the conditional expectation of the observed process (on log-scale), as a latent process, is (multivariate) normal. There are extensive works related to latent spatio-temporal models under the Bayesian framework, including models with Gaussian data modelled by (multivariate) Gaussian process with an additive error [10, 11, 12, 13], Poisson data with conditional expectation modelled by Gaussian latent process ([14, 15] and Chapter 7 of [8]) and Poisson data with multivariate log-gamma latent process [16]. However, estimation of parameters in parameter-driven models requires considerable computational effort, as does prediction of the latent process.

On the other hand, in observation-driven models, inference is possible in a (penalized) maximum likelihood framework and therefore can be easily fitted even for quite complex regression models [17]. Schrödle, Held, and Rue [18] proposed a parameter-driven spatio-temporal model and compared it with a similar observation-driven model proposed by Paul, Held, and Toschke [19]. They conclude that the parameter-driven models perform slightly better in terms of prediction in some cases, however, while the computation time for the observation-driven model is mostly less than a second, fitting a parameter-driven model takes several hours if it ever converges, because of the complexity with the latent autoregressive process. Besides, their model contains only five parameters, while in our application, the number of parameters of interest grows quadratically with the number of cell populations, which will make the parameter-driven models intractable even with a moderate number of cell populations.

Therefore, we choose to work with a spatial extension of observation-driven time series. Zeger and Qaqish [20] review various observation-driven time series models with a quasi-likelihood estimation. Fokianos and Tjøstheim [21] develop and study the probabilistic properties of a log-linear autoregressive time series model for Poisson data, as an extension of the model considered by Fokianos, Rahbek, and Tjøstheim [22]. See Scott, et al. and Kedem and Fokianos [23, 24] for a complete review.

Literature about observation-driven spatio-temporal models, however, is relatively sparse. Held, Höhle, and Hofmann [25] propose a multivariate time series model where parameters are allowed to vary across space. Paul et al. [19] extended the model such that spatial dependences are captured by additional parameters that quantify the “directed influence” of neighboring areas at previous time points on the observation of interest. Paul and Held [26] further extend the model by introducing random effects. Note that these approaches model directly the conditional expectation of the count data, meaning they are using an identity link function, instead of the canonical log-link. Thus, it is required that the parameters are positive to ensure that the resulting conditional expectation is positive. Knorr-Held and Richardson [27] propose a space-time model for surveillance data, apart from separate seasonal and spatial components, they include an autoregressive term with a latent indicator.

In this paper, we develop a conditional spatial-temporal model for multivariate count data on tiled images, and provide its application on tiled images in the context of longitudinal cancer cell monitoring experiments. Our model enables us to measure the effect on the growth rate of each cell population and changes due to local cross-population interactions. Specifically, we consider a multivariate Poisson model with intensity modeled as a log-linear form similar to those in [27] and [21], and we quantify spatio-temporal impacts of different cell populations in neighboring tiles through model parameters, as illustrated in Figure 1(b). Impacts are allowed to be positive or negative, and unlike those models that describe between group dependence through a covariance matrix, influences do not have to be symmetrical in our model. Another main advantage of the proposed framework is that it enables one to accommodate spatio-temporal cell interactions for heterogeneous cell populations within a relatively parsimonious statistical model.

Since the model complexity can be potentially very large in the presence of many cell types, it is also important to address the question of how to select an appropriate model by retaining only the meaningful spatio-temporal interactions between cell populations We cary out model selection using the common model selection criteria for parametric models, the Akaike and the Bayesian information criteria (AIC and BIC).

The remainder of the paper is organized as follows. In Section 2, we introduce the conditional spatio-temporal lattice model for multivariate count data and develop maximum likelihood inference tools. In the same section, we discuss the asymptotic properties of our estimator and standard errors. In Section 3, we study the performance of our methodology using simulated data, and compare it to that of the multivariate conditional autoregressive (MCAR) model. In Section 4, we apply our method, as well as the MCAR model to analyze datasets from an in-vitro experiment, where cancer cells are co-cultured with fibroblasts. In Section 5, we conclude and give final remarks.

2 Methods

3 Monte Carlo simulations

In our Monte Carlo experiments, we generate data from a Poisson model as follows. At time t=0, we populate nL tiles using equal counts for cells of different colors. For t=1,…,T, observations are drawn from the multivariate Poisson model Yi,t(c)|Yt−1∼Poisson(λi,t(c)),c∈C. Recall that the rate λi,t(c) defined in Section 2.1 contains autoregressive coefficients β(c|c′), which are collected in the nC×nC matrix B.

We assess the performance of MLE under different settings concerning the size and sparsity of B. Consider the three models with the following choices of B:

Denote Model i as the model corresponding to Bi,i=1,2,3. In Model 1, all the effects in B have the same size; in Model 2, the effects have decreasing sizes; Model 3 is the same as Model 1, but with some interactions exactly equal to zero.

We set α(1)=⋯=α(nC)=−0.1 for all three models. The above parameter choices reflect the situation where the generated process Y has a moderate growth.

In Table 1 and Table 2, we show results based on 1000 Monte Carlo runs generated from Models 1-3, for n=25,nC=3 and T=10 and 25. In Table 1, we show Monte Carlo estimates of squared bias and variance of θˆ. Both squared bias and variance of our estimator are quite small in all three models, and decrease as T gets larger. The variances of Model 2 are slightly larger than those in the other two models due to the increasing difficulty in estimating parameters close to zero.

In Table 2, we report the coverage probability for symmetric confidence intervals of the form θˆ±z1−α/2sdˆ(θˆ), where zq is the q−quantile for a standard normal distribution, with α=0.01,0.05,0.10. The standard error, sdˆ(θˆ), is obtained by the squared root of diagonal elements of Vn(θˆ) and the parametric bootstrap estimate, Vˆest and Vˆboot, described in Section 2.3. The coverage probability of the confidence intervals are very close to the nominal level for both methods.

In Table 3, we show results for the model selection based on 1000 Monte Carlo samples from Model 3 using the AIC and the BIC given in Section 2 for n=25 and T=10,25. We report Type A error (a term is not selected when it actually belongs to the true model ) and Type B error (a term is selected when it is not in the true model ). For both AIC and BIC model selection is more accurate for large T. As expected AIC tends to over select, and BIC outperforms AIC, with zero Type A error, and very low Type B error.

Finally, we compare the performance of our model with the following Multivariate conditional autoregressive (MCAR) model proposed by Leroux, Lei, and Breslow [32]:

where Zi,i∈Ln are random effects with conditional distribution

where ρ is a spatial autocorrelation parameter, with ρ=0 corresponding to independence, while ρ=1 corresponds to the intrinsic model, and ΣZ is a nCT×nCT between variable covariance matrix, which is assumed to have no fixed structure, and ni is the number of tiles in a neighborhood of tile i as defined in Section 2.1. Let β=(αT,vec(B)T)T be a vector of regression parameters, where B is defined in Section 2.1 and α is the intercept. Let the covariate xi,t be a nC2-dimensional vector consists of nC vectors: (Si,t−1(1),…,Si,t−1(nC)), where Si,t−1(c) carries the information from the neighbouring tiles on the previous time point, defined in eq. (2).

An independent Gaussian prior, N(0,100000), is specified for each regression parameter in β. A uniform prior on the unit interval, U(0,1), is specified for ρ. For covariance matrix ΣZ, assume an inverse Wishart distribution with identity scale matrix and nCT degree of freedom.

To evaluate the performance of MLE under our model and estimators obtained by the MCAR model, we generate 1000 set of data from Model 1. Estimation of the MCAR model is done by MCMC sampling, using R package CARBayes by Lee [33]. Table 4 show Monte Carlo estimates of squared bias, variance, the coverage probability of 95% confidence intervals and computation time for n,T∈{10,25} and nC=1,2,3. Two of the settings are the same as those shown for Model 1 in Table 1: n=25,nC=3,T=10 and n=25,nC=3,T=25 . In estimation of MCAR, we also show results of two MCMC settings: 1. MCAR1: 1000 MCMC samples generated and 200 discarded as the burn-in period; 2. MCAR2: 5000 samples with 100 discarded. Coverage probabilities of our model is computed as θˆ±z0.975sdˆ(θˆ), where zq is the q−quantile for a standard normal distribution. The standard error, sdˆ(θˆ), is obtained by taking the squared root of diagonal elements of Vn(θˆ) described in Section 2.3.

In overall, our method performs better than MCAR at analysing the kind of data that we generate, especially when n and/or T is small, with much smaller bias and variance, as well as computation time. The performance of MCAR improves significantly as the model gets more complicated (i.e. larger nC), and when n and T increases. In the case where n=25,T=25 and nC=3, it almost performs equally well with our model, however, it takes almost an hour to obtain the estimates, while our method requires less than a minute. Besides, for the coverage probabilities to reach the nominal level, it seems that MCAR requires larger MCMC sample size as the model gets more complicated, while those of our model has been stable and close to the nominal level in all cases.

4 Analysis of the cancer cell growth data

Cancer cell behaviour is believed to be determined by several factors including genetic profile and differentiation state. However, the presence of other cancer cells and non-cancer cells has also been shown to have a great impact on overall tumor behaviour [34, 35]. It is therefore important to be able to dissect and quantify these interactions in complex culture systems. The data sets in this section represent a cancer cell-fibroblast co-culture experiment. The data sets analyzed consist of counts of cell types (different cancer cell populations expressing different fluorescent proteins, and non-fluorescent fibroblasts) from 9 subsequent images taken at an 8-hour frequency over a period of 3 days using the Operetta high-content imager (Perkin Elmer). Information regarding cell type (fluorescent profile) and spatial coordinates for each individual cell were extracted using the associated software (Harmony, Perkin Elmer).

Each image was subsequently tiled using a 25×25 regular grid.

We choose the number of tiles for a balance between the fit of the model and capturing the local impact between cell populations. More specifically, decreasing tile sizes enables one to detect local impacts between cell populations, which is one of the objectives of our analysis. However, if the tiles are too small, we will end up with mostly no cells in most tiles. In this situation the conditional Poisson model would not fit well the data. On the other hand, when the tiles are too large the model would fit the data well (the conditional Poisson would be approximately a conditional normal model), but we lose information on local impacts. We recommend 0 to 20 average cells per tile, since for such choice our diagnostic and goodness-of-fit analyses suggest that the conditional Poisson model fits well the data whilst enabling us to measure local correlation effects between populations.

4.1 Cancer cell-fibroblast co-culture experiment

In this experiment, cancer cells are co-cultured with fibroblasts, a predominant cell type in the tumor microenvironment, believed to affect tumor progression, partly due to interactions with and activation by cancer cells [34]. In this experiment, fibroblasts (F) are non-fluorescent whereas cancer cells fluoresce either in the red (R) or green (G) channels due to the experimental expression of mCherry or GFP proteins, respectively. Cells were initially seeded at a ratio of 1:1:2 (R:G:F).

Model selection and inference.   We applied our methodology to quantify the magnitude and direction of the impacts have on growth for the considered cell types. To select the relevant terms in the intensity expression (1), we carry out model selection using the BIC model selection criterion. In Table 5, we show estimated parameters for the full and the BIC models, with bootstrap 95% confidence intervals in parenthesis. Figure 2 illustrates estimated spatio-temporal impacts between cell types using a directed graph. The solid and dashed arrows represent respectively significant and not significant impacts between cell types at the 95% confidence level. Significant impacts coincide with parameters selected by BIC.

The interactions within each cell type (βˆ(c|c),c=R,G,F) are significant, which is consistent with healthy growing cells. As anticipated, the effects βˆ(c|c) for the cancer cells are larger than those for the slower growing fibroblasts. The validity of the estimated parameters is also supported by the similar sizes of the parameters for the green and red cancer cells. This is expected, since the red and green cancer cells are biologically identical except for the fluorescent protein they express. Interestingly, the size of the estimated effects within both types of cancer cells (βˆ(c|c),c=R,G) are larger than the impact they have on one another (βˆ(G|R) and βˆ(R|G)). This is not surprising, since βˆ(c|c)(c=R,G) reflects not only impacts between cells from the same cell population, but also cell proliferation. The fact that we are able to detect the impacts between the red and green cancer cells confirms that our methodology is sensitive enough to detect biologically relevant impacts even though no interactions were found between the cancer cells and the fibroblasts. This might be due to the fact that we used normal fibroblasts that had not previously been in contact with cancer cells and thus had not been activated to support tumor progression as is the case with cancer-activated fibroblasts.

Figure 2:

Directed graph showing fitted spatio-temporal interactions between GFP cancer cells (G), mCherry cancer cells (R) and fibroblasts (F). The solid and dashed arrows represent respectively the significant and not significant interactions between cell types at the 95% confidence level.

Goodness-of-fit and one-step ahead prediction   To illustrate the goodness-of-fit of the estimated model, we generate cell counts for each type in each tile, yˆi,t(c), from the Pois(λˆi,t(c)) distribution for t≥1, where λˆi,t(c) is computed using observations at time t−1, with parameters estimated from the entire dataset. In Figure 4, we compare the actually observed and generated cell counts for GFP cancer cells (G) and mCherry cancer cells (R) and fibroblasts (F) across the entire image. The solid and dashed curves for all cell types are close, suggesting that the model fits the data reasonably well. As anticipated, the overall growth rate for the red and green cancer cells are similar, and sensibly larger than the growth rate for fibroblasts.

To assess the prediction performance of our method, we consider one-step-ahead forecasting using parameters estimated from a moving window of five time points. In Figure 3, we show quantiles of observed cell counts against predicted counts for each tile. The upper and lower 95% confidence bounds are computed non-parametrically by taking Fˆ1−1(Fˆ0(yt(c))−0.95) and Fˆ1−1(Fˆ0(yt(c))+0.95), where Fˆ0 and Fˆ1 are the empirical distributions of the observations and predictions at time t respectively [36]. The identity line falls within the confidence bands in each plot, indicating a satisfactory prediction performance.

Figure 3:

QQ-plots for cell growth, comparing observed (horizontal axis) and one-time ahead predicted (vertical axis) cell counts per tile on the entire image at times t=6,7,8 for GFP cancer cells (G), mCherry cancer cells (R) and fibroblasts (F). One-time ahead predictions are based on the model fitted using a moving window of five time points.

Comparison with MCAR model   Next, we compare the estimates as well as the goodness-of-fit on the real data with the MCAR model. Parameter estimates are shown in Table 5, with 95% confidence intervals given in parenthesis. Results from both models are mostly consistent with each other, specifically, both models show that impacts within each cell type (βˆ(c|c),c=R,G,F) are significant, the effects βˆ(c|c) for cancer cells are larger than those for the slower growing fibroblasts, the green and red cancer cells have positive impact on each other, and cancer cells have no impact on fibroblasts. The only difference is, the MCAR model shows a negative impact of fibroblasts on the green cancer cells only, while our model detect no significant impact on either cancer cells. Since the red and green cancer cells are biologically identical except for the fluorescent protein they express, we expect a symmetrical result with both cancer cells.

Table 5:

Estimated parameters for the full, the BIC models and the MCAR model based on the cancer cell growth data described in Section 4. Bootstrap 95% confidence intervals based on 50 bootstrap samples are given in parenthesis.

	Full model
c=	G	R	F
αˆ(c)	−0.99 (−1.19, −0.79)	−0.50 (−0.70, −0.30)	−0.26 (−0.45, −0.06)
βˆ(G|c)	1.23 (1.10, 1.35)	0.34 (0.21, 0.48)	0.12 (−0.03, 0.27)
βˆ(R|c)	0.28 (0.17, 0.38)	1.09 (0.96, 1.21)	0.02 (−0.09, 0.13)
βˆ(F|c)	0.10 (−0.01, 0.21)	0.02 (−0.07, 0.12)	0.92 (0.81, 1.03)
	BIC model
c=	G	R	F
αˆ(c)	−0.88 (−1.04, −0.72)	−0.49 (−0.66, −0.31)	−0.19 (−0.36, −0.02)
βˆ(G|c)	1.24 (1.11, 1.37)	0.35 (0.21, 0.48)	/
βˆ(R|c)	0.28 (0.17, 0.39)	1.09 (0.96, 1.21)	/
βˆ(F|c)	/	/	0.93 (0.82, 1.04)
	MCAR
c=	G	R	F
αˆ(c)	−0.45 (−0.54, −0.38)	−0.45 (−0.54, −0.38)	−0.45 (−0.54, −0.38)
βˆ(G|c)	1.06 (0.93, 1.16)	0.16 (0.09, 1.16)	−0.15 (−0.22, −0.09)
βˆ(R|c)	0.25 (0.15,0.31)	1.01 (0.92,1.08)	0.05 (−0.02,0.10)
βˆ(F|c)	0.03 (−0.06,0.20)	0.03 (−0.07,0.19)	0.96 (0.83,1.08)

In Figure 4, apart from the observed (solid curve) and generated (dashed curve) cell counts from our model, we also show the generated cell counts from the MCAR model (dotted curve) for the green cancer cells (G), red cancer cells (R) and fibroblasts (F) across the entire image. Compared to the dotted curves, the dashed curves are slightly closer to the solid ones, which means our model seems more appropriate for analysing this type of data than the MCAR model.

Figure 4:

Goodness-of-fit of the estimated models. Observed (solid) and predicted (dashed for our model and dotted for the MCAR model) number of GFP cancer cells (G), mCherry cancer cells (R) cancer cells and fibroblasts (F) for the entire image. Predicted cell counts for each cell type in each tile yˆi,t(c) is generated from the conditional Poisson model with intensity λˆi,t(c) defined in eqs. (1) and (2), where the coefficients βˆ(c|c′) are estimated from the entire dataset.

5 Conclusion and final remarks

In this paper, we introduced a conditional spatial autoregressive model and accompanying inference tools for multivariate spatio-temporal cell count data. The new methodology enables one to measure the overall cell growth rate in longitudinal experiments and spatio-temporal interactions with either homogeneous or heterogeneous cell populations. The proposed inference approach is computationally tractable and strikes a good balance between computational feasibility and statistical accuracy. Numerical findings from simulated and real data in Sections 3 and 4 confirm the validity of the proposed approach in terms of prediction, goodness-of-fit and estimation accuracy.

The data sets described in this paper serve as a proof-of-concept that the proposed methodology works. However, the potential applications and the relevant questions that the methodology can help to answer in cancer cell biology are plentiful. To build on from the examples given in this paper, the methodology can be used to study interactions between cancer cells and a wide range of cancer-relevant cell types such as cancer-activated fibroblasts, macrophages, and other immune cells when co-cultured. Since a substantial proportion of cancer cells in tumors are in close proximity to other cell types that have been shown to affect tumor progression, using these co-cultures is more representative of the situation in a patient compared to studying cancer cells on their own. In addition to just giving the final cell number, the presented approach can dissect which cell types affect the growth of others and to what extent in complex heterogeneous populations. This could be relevant in a drug discovery setting to determine if a drug affects cancer cell growth due to internal effects (on other cancer cells) or by interfering with the interaction between the cancer cells and other cell types. Finding drugs with different targets and mechanisms of action are particularly sought after as they provide a wider target profile, increasing the chance of patients responding as well as reducing the risk of tumors becoming resistant. The impact of different genes and associated pathways in different cell types in relation to inter-cellular interactions can also be studied by genetically modifying the cell type(s) in question before mixing the cells together. This could be beneficial to identify new potential drug targets. Our approach is also applicable in other kinds of studies where local spatial cell-cell interactions are believed to affect cell growth such as studies of neurodegenerative diseases [37] and wound healing/tissue re-generation [38]. In addition to evaluating cell growth, our approach can also be used to study transitions between cellular phenotypes upon interaction with other cell types, provided that the different phenotypes studied can be distinguished from one another based on the image data. Finally, it is worth noting that issues may arise when cells become too confluent/dense, this may lead to segmentation problems of the imaging system. If they become completely confluent, they are likely to progressively stop growing. If one wants to measure for longer period of time, experiments can be performed in larger wells/plates or with smaller starting cell numbers.

Our method offers several practical advantages to researchers interested in analysing multivariate count data on heterogeneous cell populations. First, the conditional Poisson model does not require tracking individual cells across time, a process that is often difficult to automate due to cell movement, morphology changes at subsequent time points, and additional complications related to storage of large data files. Second, we are able to quantify local spatio-temporal interactions between different cell populations from a very simple experimental set-up where the different cell populations are grown together in a single experimental condition (co-culture). An alternative, solely experimentally-based strategy would require monitoring the different cell types alone and together at different cell densities (number of cells per condition) in order to make inferences in terms of potential interactions. However, such an approach would give no possibility of evaluating the spatial relations in the co-culture conditions and would still restrict the number of simultaneously tested cell types to two.

In the future, we foresee several useful extensions of the current methodology, possibly enabling the treatment of more complex experimental settings. First, complex experiments involving a large number of cell populations, nC, would imply an over-parametrized model. Clearly, this large number of parameters would be detrimental to both statistical accuracy and reliable optimization of the likelihood objective function ℓn(θ) in (4). To address these issues, we plan to explore a penalized likelihood of form ℓn(θ)−penλ(θ), where pen(θ) is a nonnegative sparsity-inducing penalty function. For example, in a different likelihood setting, Bardic et al. [39] consider the L1-type penalty pen(θ)=λ∑|θ|,λ>0.

Second, for certain experiments, it would be desirable to modify the statistics in eq. (2) to include additional information on cell growth such as the distance between heterogeneous cells, and covariates describing cell morphology.

Thirdly, it would be useful to develop a more principled way to select the tile sizes/number, and consider tiling the microscope image into a hexagonal lattice, which is a more natural choice in real application, since the distance between neighboring tiles would be more even than that of a regular lattice.

Finally, although numerical results (results not reported here) show that our method are quite robust in the presence of mild outliers (with around 5% of contaminated data), for more severe situations, we expect that severe or numerous outliers will have some influence on the estimates since the Poisson score function is unbounded. To address this problem, the log-likelihood scores in eq. magenta (5) should be replaced by some other robust alternative. Following Ferrari and Vecchia [29] and La Vecchia et al. [30], robustness can be obtained by the so-called q-entropy estimation method simply obtained by replacing the usual logarithm in the log-likelihood estimating equation by the q-logarithm logarithm function logq(x)=(x1−q−1)/(1−q)if q≠1, and logq(u)=log(x) if q=1, for all x>0. This ensures a bounded influence function for the implied estimator and therefore guarantees control of the bias under contamination.