## A parallelized cellular Potts model that enables simulations at tissue
scale

Shabaz Sultan1,2, Sapna Devi3, Scott N. Mueller3, Johannes Textor1,2

1 Data Science group, Institute for Computing and Information Sciences,
Radboud University, Nijmegen, The Netherlands

2 Department of Medical BioSciences, Radboud University Medical Center,
Nijmegen, The Netherlands

3 Department of Microbiology and Immunology, The University of Melbourne, at
the Peter Doherty Institute for Infection and Immunity, Melbourne, Victoria,
Australia

## Summary

The Cellular Potts Model (CPM) is a widely used simulation paradigm for
systems of interacting cells that has been used to study scenarios ranging
from plant development to morphogenesis, tumor growth and cell migration.
Despite their wide use, CPM simulations are considered too computationally
intensive for three-dimensional (3D) models at organ scale. CPMs have been
difficult to parallelise because of their inherently sequential update scheme.
Here, we present a Graphical Processing Unit (GPU)-based parallelisation
scheme that preserves local update statistics and is up to 3-4 orders of
magnitude faster than serial implementations. We show several examples where
our scheme preserves simulation behaviors that are drastically altered by
existing parallelisation methods. We use our framework to construct tissue-
scale models of liver and lymph node environments containing millions of cells
that are directly based on microscopy-imaged tissue structures. Thus, our GPU-
based CPM framework enables _in silico_ studies of multicellular systems of
unprecedented scale.

## Introduction

The CPM [1] is a powerful simulation method that is known for its ability to
simulate a wide range of biological scenarios. When modelling tissue
environments, many researchers wish to consider both the complex morphology of
single cells, and the behavior of these cells as part of larger collectives.
The appeal of the CPM lies partly in its ability to integrate both these
aspects in a single model. As such, it has been applied to diverse problems
including tumor growth [2], vascularization [3], wound healing [4], and
morphogenesis [5, 6]. Many extensions and hybrid model variants are available
in the literature to, for example, add support for cell motility [7, 8, 9],
cell signalling [10], and finite-element based mechanical interactions [11].
Thus, the CPM can be described as a generic platform upon which a wide range
of models can be built, possibly explaining its popularity in the tissue
modelling field.

For computational reasons, CPM models are often built using a small number of
cells (perhaps in the thousands) in simplified two-dimensional (2D)
geometries. Some research questions in tissue modelling however demand a
realistically large number of cells, and accurate 3D representations of
tissues. For instance, it is widely appreciated that the efficiency of space
exploration of motile cells depends critically on the dimensionality of the
environment, leading researchers to use 3D CPM simulations [9, 12, 13].
Assessment of tumor growth behavior and invasive patterns of tumors likewise
changes as a function of dimension, resulting in the use of 3D CPM simulations
in the tumor modelling field as well [2]. While there is thus a clear research
demand for realistically large 3D CPM simulations, scaling up such simulations
is generally considered computationally infeasible [14, 15, 16]. Removing this
computational barrier would allow a host of CPM based modelling research to be
scaled to what is biologically realistic, instead of scaling to what technical
constraints allow.

The CPM is an on-lattice method: simulations are executed on a grid. A common
strategy to speed up grid-based models is to cut the simulation into smaller
subsections that are run in parallel. However, for the CPM, model behavior is
intimately tied to locations on the grid being updated sequentially, in a
single-location update Markov chain Monte Carlo process. As such,
parallelisation has been cautioned against, as it risks drastic changes in CPM
“system kinetics” [5]. Nevertheless, several attempts have been made to employ
subsection-based parallelisation. Chen _et al._ proposed a checkerboard scheme
that divides the CPM lattice up over a computer cluster, with sections chosen
much larger than cells, preventing cell-level synchronization issues [17].
Tapia _et al._ propose a similar scheme for GPUs [18]: because the sections
are much smaller than a cell on the GPU, they use a per-cell locking mechanism
to prevent simultaneous updates of a cell. Neither approach has been
extensively checked for correct model behavior. To our knowledge, no publicly
available implementation of subsection-based parallelisation exists for the
CPM. The software “Morpheus” [19] includes a non-subsection-based Central
Processing Unit (CPU) parallelisation strategy; this method is less suited to
fully utilise parallel accelerator hardware such as GPUs.

In this paper, we aim to make several contributions: First, we introduce a
highly optimized parallelized implementation of the CPM using several
optimization strategies not previously explored in the literature. Second, we
conduct an empirical analysis of potential pitfalls in parallelising a model
dependent on serial updates of single locations and show that our approach
mitigates these issues while preserving performance benefits. Third, we
demonstrate that our approach retains the behavior of cell motility
simulations that are particularly prone to changed model kinetics upon
parallelisation. And finally, we build two large 3D image-based simulations of
T cell migration in liver and lymph nodes, demonstrating how complex tissue
environments can be analysed at a larger scale and realism than was possible
before.

## Results

### Serial Cellular Potts Models are limited to sub-millimetre scale

Figure 1: Serial CPM implementations do not scale well to large simulations
and big ensembles. (A) Computational cost of simulating a 2563$\mu$m3 volume
(1$\mu$m3 per voxel) for 10 simulated minutes (1s per simulation step). With a
high-performance CPU implementation, this takes 1.2 hours of wall-clock time.
Extrapolation to a full day of simulated time for a 1mm3 simulation estimates
1.4 years of wall-clock time. Time measurements (black) are averages of 5
simulations, projections are based on linear fits of simulations $\geq$
643$\mu$m3 in log-log space (grey). (B) Building CPMs typically requires
tuning parameters. Here we perform a multidimensional grid search for
parameters that yield motile cells using the Act extension of the CPM [7]; a
modest single parameter grid search where 10 values are tried takes 15.9
seconds; larger grid searches take exponentially more time as function of
parameters searched, with a five parameter search with 30 values is
extrapolated to 1.1 years of wall-clock time. Time measurements (black) are
averaged over 5 simulations, projections fitted on linear trends in linear-log
space. CPM parameters are shown in Table 1 and Table 2.

The CPM’s serial update schedule – only one location at a time is updated on a
grid – is central to its behavior, and thus not trivially changed. We
therefore first aimed to quantify the computational limitations for models
based on the serial CPM: how big can our computational experiments be, and
what timescale can they span, without making fundamental changes to the CPM?

We set up 3D CPM tissue experiments at different spatial and temporal scales
and measured the required wall-clock time. We choose spatial and temporal
simulation resolutions that allow reasonably accurate representations of
highly motile cells within tissues. Lymphocytes are among the smallest
nucleated cells in the mammalian body at 8-10μm diameter for B cells and T
cells. At speeds of up to a few body lengths per minute, they are also among
the fastest movers within tissues. We chose 1μm3 per voxel and 1 second per
simulation step, or Monte Carlo Step (MCS), as our simulation resolutions,
corresponding to Beltman _et al._ [9].

To create realistically large models of, for example, a complete murine Lymph
Node (LN), or scale a tumor model to the size of an avascular tumor spheroid
[20], experiments at mm3 scale would be needed. Certain organ scale dynamics
require timescales of at least a day, such as T cell transit through LNs [21].
More modest simulation experiments at 2563μm3 (Figure 1A) at 10 simulated
minutes require ~1h of wall-clock time on a highly optimized CPU
implementation of the CPM (Methods). Such experiments are sufficient to see
minimal movement of motile cells such as T cells and can capture organ
compartments such as a LN paracortical T cell zone, or a micrometastasis for
cancer modelling. For mm3 scale, one-day experiments however, researchers
would have to wait at least a year – making them impractical for all but the
most patient of researchers.

A single experiment is rarely sufficient, in both empirical and simulation
research. Beyond the need for experimental replicates, simulations can have
parameters for which associated model behavior is non-obvious, requiring
parameterized simulation ensembles. In Figure 1B, we explore a cell motility
extension to the CPM that coarsely models actin dynamics [7]. A single two-
dimension simulation is relatively cheap, but when exploring even a modest
number of parameters the execution time goes up exponentially.

Thus, serial 3D CPM simulations at a scale of 1 mm3 are currently not
feasible, and even 1003μm3 scale simulations were extremely challenging in our
experiments. While ensemble runs can be accelerated with parallel simulations
without changing the CPM, the high computational costs of running large
ensembles further motivates the need for a faster CPM.

### Designing a parallel execution scheme for the CPM

To make large scale simulations practical, we propose a parallelisation scheme
for the CPM that uses highly parallel GPU hardware. A natural way to
parallelise on-lattice models, such as the CPM, is to divide the simulation
grid in smaller sections; each section runs a local simulation in parallel. As
GPUs support up to thousands of calculations in parallel, these sections will
be quite small. The CPM also updates information outside of the grid, such as
sizes of simulated cells. We first consider three issues when designing our
parallelisation scheme: adjacent sections with active simulation interfering
with each other during updates of the simulation grid, non-grid information
such as cell sizes being updated simultaneously, and stale information being
used by one parallel simulation after said information got updated by another
simulation.

Figure 2: Designing a parallel CPM based on checkerboarding and atomic updates
for GPUs. (A) We divide simulation space into a subgrid, with alternating
regions being active in a “checkerboard” pattern, preventing update conflicts
on the grid. (B) We implement two strategies to prevent global variable update
conflicts (e.g. cell size tracking): a mutex locking strategy (Tapia _et al._
[18]), and a novel atomic operation strategy. (C) Combining interleaved
storage and atomic operations on adjacent data allows variables to be kept in
sync with each other. (D) A four-step execution scheme, where active
checkerboard areas are transferred to fast shared memory, accelerates copy
attempts. Thread barriers synchronize active threads after transfers. (E) A
multi-GPU execution scheme, where regions of the simulation are assigned to
different GPUs. Borders are synchronized between relevant adjacent regions,
and global variables are synchronized through delta tracking.

A CPM simulation proceeds by choosing a random position on the simulation
grid, and stochastically copying the position’s state to a neighbour. If
adjacent subsections of the grid each have an active parallel simulation
running, we run the risk of simulations “racing” to write to the same location
on the sections’ boundary. Such “race conditions” run the risk of information
getting lost, as one section overwrites an almost simultaneous update –
without incorporating information of the overwritten update. This is a
standard problem in parallelisation, with a straight-forward solution: never
allow two adjacent sections to be active at the same time. A spaced out
“checkerboard” pattern of sections is only ever active; active sections are
sufficiently spaced to prevent overlapping updates (Figure 2A). Regions cycle
periodically between active and dormant state.

The CPM also tracks information involving more than one grid location, such as
the current sizes of simulated cells. We face the same “race condition”
problem here, as it is possible that multiple parallel updates touch a single
cell (Figure 2B, left column). GPU parallelisation requires a high number of
small active subsections, making such concurrent updates more likely. In Tapia
_et al._[18] a per-cell locking scheme is proposed, that rejects all but the
first concurrent update attempt on a variable (Figure 2B, middle column). We
instead opt to use atomic operations (Figure 2B, right column), which equally
prevents update information from getting lost, as the hardware will eventually
integrate all updates for such operations.

In serial CPMs, global information travels instantaneously: as only one grid
position is updated at a time, every update attempt has access to the results
of all previous updates. This instantaneous information access is relaxed in
our proposed scheme: when information such as cell size is accessed, this
information may be out-of-date due to still in-flight updates to the accessed
information. The cell-locking scheme from Tapia _et al._ prevents this from
happening, by outright rejecting any updates that could have used out-of-date
information; we explored a less strict scheme, reasoning that the stochastic
CPM copy attempts may tolerate out-of-date information. For CPM extensions
that do not tolerate their global variables being slightly out sync with each
other (e.g. variables used in center-of-mass “centroid” tracking) we designed
an additional scheme that still will keep variables in sync with each other,
and sits on top of our atomics based scheme (Figure 2C, Methods). Otherwise,
due to atomic updating, our scheme still gives strong guarantees on tracked
data being correct – for example, cell sizes will be tracked correctly – and
the added flexibility allows us to explore the effects of allowing the use of
slightly outdated information.

Implementation of our scheme on real hardware requires consideration of
several technical details. On GPUs, memory access almost always bottlenecks
performance – but the CPM’s random access pattern prevents easy use of a GPU’s
caching hierarchy [18]. To maximize performance, we designed a data transfer
scheme that carefully operates in lock-step with checkerboard work splitting
(Figure 2D, Methods): we cache local data needed by the regions that are
becoming active, and overlap memory transfer and simulation calculations to
fully utilise GPU hardware. To increase simulation volume, and gain further
speed ups, we also prototyped a multi-GPU scheme (Figure 2E, Methods).

In short, our proposed parallel CPM scheme is tailored to make optimal use of
GPUs: hardware specialized in running parallel code. It still gives strong
guarantees on simulation state consistency – for example, a global variable
tracking a cell’s size will always equal the number of pixels on the grid
occupied by that cell. Yet, we do allow the use of slightly outdated values of
such global variables in copy attempts to gain flexibility and further
maximise performance gains.

### Aligning the lattice update statistics of parallel CPMs to serial CPMs

Figure 3: Checkerboarding can drastically alter local update statistics. (A)
To prevent threads from racing ahead of others too much, thread barriers are
imposed after each update. (B) Checkerboard size and switching frequencies
determine the distribution of waiting times between consecutive updates of the
same locations. (C) Small checkerboard areas or infrequent switches can alter
the waiting time distribution substantially compared to the single-update
scheme. Update frequency histograms were generated by simulating one MCS (256
updates) on a 16x16 grid. Waiting time histograms were simulated on the same
grid using 100 simulated MCSs.

Our novel parallel execution scheme guarantees simulation data consistency,
but drastically changes the way the simulation lattice is updated. We
investigated statistical properties of copy attempts on the simulation lattice
for different lattice positions; specifically, we analysed the relative
ordering of copy attempts, frequency of copy attempts, and waiting times
between copy attempts. This allowed us to quantify the divergence of parallel
lattice update kinetics from the serial CPM. Such divergence could lead to
changed simulation behavior in complex models of cells and tissues.

Without guardrails, performing independent CPM updates in parallel would allow
different parts of a simulated tissue to have been alive for different lengths
of simulated time (Figure 3A). Modern GPUs use thread barriers to synchronize
parallel threads without strong computational performance implications. We set
thread barriers after ever single copy attempt, so for every currently active
checkerboard region there is only a maximum time divergence of one copy
attempt between active regions.

Our parallel execution scheme allows us to freely choose the size of
subregions and frequency of active region switching (Figure 3B), without
impacting data consistency guarantees. We measure frequency and waiting times
between copy attempts per position on the lattice in the standard, serial
implementation of the CPM (Figure 3C, left column). For our initial parameter
choice, we found that underlying update statistics are drastically changed
(Figure 3C, middle column): waiting times between copy attempts became much
more predictable due to the much more even distribution of copy attempts
compared to the serial CPM. Using large parallel subsections with frequent
active region switching allows non-parallel update statistics to be closely
recapitulated in our parallel scheme (Figure 3C, right column).

We thus find that parallelisation parameters have major influence on copy
attempt statistics on the lattice. By switching active regions of the
checkerboard frequently, and avoiding too small areas per thread, we can
obtain statistics that closely resemble the serial CPM.

### Emergent cell and tissue behavior is maintained in a parallel CPM

Figure 4: Atomic updates lead to correct behavior in a simulation of cell
sorting. A cell sorting simulation based on Graner _et al._ [1] is implemented
with two strategies to prevent race conditions for cell size variables. The
mutex locking strategy proposed by Tapia _et al._ [18] leads to severe
starvation in our tests. Using an atomic update approach recovers the desired
sorting behavior, closely matching sorting speed of the original CPM. The
simulation was run with a 4x4 area per thread and 2x2 checkerboards, but other
sizes were tested and did not make a significant difference for this
simulation. Sorting trends of different CPM implementations are averaged over
5 simulations. Simulation parameters can be found in Table 1. Figure 5:
Motility behavior matches serial CPM behavior when parallel CPM local
statistics are similar. (A) Cells in a circular chemokine gradient simulation
move highly regular in a serial CPM (low per-pixel entropy, shaded: 1SD);
parallel simulations diverge slowly from reference behavior when copy attempt
statistics are preserved (low per-pixel cross-entropy). (B) In the Act model
of intra-cell polarization [7], stronger polarization strength leads to higher
directional autocorrelation (measured for 750 MCSs at 50 MCS intervals,
averaged over 50 simulations, shading: $\pm$1SD); directional persistence only
emerges when update waiting time distributions are preserved. (C) Phase
transition simulation in 3D where synchronous motility emerges at a density
around 0.5. Cell order indices after 2000 MCSs are shown (averaged over 5
simulation runs, shading: $\pm$1SD). The phase transition only occurs when
update waiting time statistics are preserved. Simulation parameters can be
found in Table 1 and Table 2.

Our goal is to build parallel CPMs that are much faster than previous
implementations but maintain cell and tissue kinetics consistent with existing
serial CPMs. We hypothesized that aligning the waiting times between update
attempts of our parallel CPM to the serial version should suffice to achieve
such consistency. To investigate this, we ran several simulation scenarios
with different CPM variants, measuring emergent behaviors on the cell and
tissue levels and comparing these to serial CPM versions.

We first investigated the classic cell sorting simulation [1]: two mixed
populations of cells with differential adhesion (i.e., cells of the same type
have stronger adhesion to each other than cells of different types) sort
themselves in two homogeneous regions over time (Figure 4, grey lines). The
mixing kinetics can be measured by tracking homogeneous and heterogeneous cell
neighbours over time [1]. Strikingly, the mutex based parallelisation approach
proposed by Tapia _et al._[18] – while still being able to run a simulation –
leads to a lack of sorting behavior; the simulation encounters so-called
“starvation”, because only one concurrent copy attempt is allowed per cell
under this scheme (Figure 4, red lines). Our less restrictive scheme allows
expected sorting behavior to re-emerge (Figure 4, black lines). The sorting
behavior remained consistent for all parallelisation parameters we examined
(Supplementary Figure S1).

We next examined three motility model variants of the CPM, where additional
biases are placed on copy attempts to obtain cells that move persistently; in
the standard CPM movement only occurs as a function of random membrane
fluctuations. Because motile cells have a leading and a retracting edge and
need to roughly maintain their target size, copy attempts need to sufficiently
alternate between these cell sections for movement to occur. We therefore
speculated that movement models might be more affected by changes to the
waiting times between copy attempts.

The first motility model we consider adds a chemotactic vector field to the
CPM; copy attempts are more likely to succeed when aligned with the underlying
vector field (see Methods). Because the directional bias is fixed in this
model – the field does not change over time – cells move in a highly regular
fashion, and simulations with the same initial conditions closely align even
at a pixel level (Figure 5A, left column). We quantified this alignment using
pixel-wise (cross-)entropy: we measure the probability that a pixel contains a
cell in multiple simulations, where low entropy of pixel probabilities
indicated high agreement between simulations, and low cross-entropy indicates
a high level of agreement of parallel simulations with reference serial runs.
In simulations with highly changed lattice update statistics, cross-entropy
increased quickly (Figure 5A, middle column); when the update statistics are
preserved cross-entropy rises much slower (Figure 5A, right column).

Next, we analysed a model that simulates localized cell polarization [7],
where recent pixel updates at the border of a cell makes further copy attempts
involving those pixels more likely to succeed. The waiting times between copy
attempts are very important in this model, as cells start to move when many
pixels in a small neighbourhood successfully expand at the same time; this is
less likely when copy attempts are spaced out evenly. A polarization strength
parameter controls how much polarization bias is added to a cell; varying the
polarization strength interpolates motility between Brownian and directionally
persistent motility modes (Figure 5B, left column). In parallel simulations
with changed update frequency and waiting time statistics this full range of
behaviors does not appear (Figure 5B, middle column), but they re-emerge when
update statistics are preserved (Figure 5B, right column).

Finally, we examined a simulation of a tightly packed 3D collective of immune
cells [9], where each cell has a preferred movement direction vector; copy
attempts for a cell are more likely to succeed when they align with this
vector, and the vector is periodically updated to align with recent movement
of the cell. This motility model involves information tracked on the lattice,
and additional information tracked external to the simulation grid, placing
additional demands on a parallel CPM. In a similar simulation model, low
density cell collective exhibit Brownian movement, transitioning sharply to
ordered collective movement as density is increased [22]; this characteristic
phase transition is reproduced in the CPM (Figure 5C, left column). In a
parallel CPM this phase transition disappears when update statistics are
changed (Figure 5C, middle column), but re-emerges when update statistics are
preserved (Figure 5C, right column).

Taken together, these results demonstrate that a strict update scheme
preventing parallel updates of the same cell is not necessary for maintaining
valid CPM behavior and can even lead to severe “starvation” where simulations
do not progress anymore. On the other hand, preserving waiting times between
copy attempts involving the same pixel can be very important to maintain model
behavior, and was critical in all three motility scenarios we examined.

###  GPU-based CPM parallelisation accelerates models by several orders of
magnitude

Having confirmed that our parallel CPM indeed produced correct model behavior
in a range of examples, we next investigated the performance improvements
resulting from our implementation choices when compared to a serial CPM.

We first examined the performance gains resulting from our use of a shared
memory caching scheme (Figure 2D, Methods) using a 2D sorting simulation
benchmark (Figure 4). On its own, shared memory caching achieved only minimal
speedup (Figure 6A). However, for further performance improvement, we pre-
generated different versions of our code based on which CPM extensions were in
use (Methods). Combining these specialized kernels with the caching scheme
achieved a significant 10-fold speedup, more than expected from either
optimization alone (Figure 6A). Testing our multi-GPU scheme (Figure 2E,
Methods) in the same scenario achieved a further 2.3x speedup by distributing
the simulation across four GPUs (Figure 6B). Hence, careful consideration of
GPU-specific implementation choices resulted in significant performance
benefits.

Figure 6: Parallel CPMs outperform serial implementations by orders of
magnitude. (A) Performance benefits are obtained by a shared memory scheme,
and a pre-compiled kernel specialization scheme, based on a sorting simulation
extended to densely fill simulation space (Figure 4, all sorting benchmarks
are average of 5 simulations of 100 MCSs, 100 MCS burn-in). (B) Distributing
the lattice over 4 GPUs speeds up the simulation further. (C) Performance of
sorting simulations are significantly higher relative to two serial
implementations; relative speedup shown for 2048x2048 simulations. (D)
Performance of 3D cell collective simulations is also improved (Figure 5C,
average of 2 simulations of 100 MCSs, 100 MCS burn-in); relative speedup shown
for 2563 simulations. (E) Trade-offs between speed and accuracy for chemokine
field simulations (Figure 5A, accuracy: 1 - cross-entropy of parallel ensemble
relative to serial ensemble, 250 simulations per ensemble). CPM parameter
settings can be found in Table 1 and Table 2.

We next compared our parallel CPM (on a single GPU) to our own optimized CPU
implementation (Methods), and to the state-of-the-art software “Morpheus”
[19]. Using the same 2D cell sorting benchmark, we observed a ¿3,500-fold
speedup compared to our own CPU implementation, and a ¿25,000-fold improvement
relative to Morpheus (Figure 6C). On the 3D cell collective simulation (Figure
5C), a ¿600-fold speedup was achieved relative to serial CPMs (Figure 6D). It
is important to note that our benchmark hardware did not permit use of shared
memory caching for the 3D model; future hardware with more shared memory may
improve 3D performance.

As shown in Figure 3 and Figure 5, checkerboard size and switching frequency
influences model behavior: larger divisions and more frequent switching leads
to behavior concordant with serial CPMs. These update patterns also interact
with caching, further impacting performance. As the sorting simulation has low
sensitivity to parallelisation parameters (Supplementary Figure S1), we
instead benchmarked our 2D chemotaxis simulation (Figure 5A). We found that
larger subdivisions indeed slowed simulations down considerably, whereas
frequent switching led to a less marked slowdown (Figure 6E). These results
suggested a trade-off between performance and model accuracy.

Hence, our parallel CPM represents a significant speed improvement over even
well optimized serial CPMs such as Morpheus. This improvement was to a
substantial extent a result of GPU-specific optimizations and led to an
overall speedup of several orders of magnitude.

### Building an image-based model of T cell surveillance of the liver

Our parallel CPM was motivated by the need for tissue simulations at realistic
scale with realistic cell morphology. While our approach increases engineering
complexity of the model, it can also simplify model building, as it enables
direct use of imaging data to initialize the simulation environment. To
determine the feasibility of directly using real-world imaging data in our
system, we constructed a model of T cells moving in the sinusoidal regions of
the liver.

Figure 7: T cells in a liver sinusoid simulation exhibit a transition from a
non-motile to a motile state. (A) A 1024$\times$1024$\times$182 μm3 scan of a
liver sinusoid is directly imported into our simulation environment (64μm
slice shown). (B) Simulated T cells explore liver sinusoids (5000 cell
simulation shown). (C) Simulations at different cell densities explore the
environment at different speeds (single simulation of 24hr, with 1 min
snapshots). (D) An increase of explored volume per cell occurs when cell
density rises from 1,000 to 10,000 cells; 100,000 cells fill $\sim$1/3 of
sinusoid volume, causing exploration efficiency to saturate. (E) At low cell
densities most simulated cells are in an un-polarized, non-motile state; at
higher densities collisions cause cells to transition in to a polarized,
motile state (1 simulation of 30 min per density, arrest coefficient: fraction
cells ¡ 2μm/min[23], shaded: Wilson score interval based on the number of
2-minute subtracks). CPM parameter settings can be found in Table 3. (F)
Average speed of the simulated T cells at different densities.

Tissue resident memory T cells (Trms) are confined to specific organs, where
they sample their environment for signs of disease [24] – in particular, liver
Trms are key in providing sterilising immunity for malaria [25]. Trms display
a distinct motility pattern in the liver, and efficacy of malaria clearing
depends on the number of Trms present [26]. We constructed a 3D simulation
environment of liver sinusoids – capillaries that Trms reside in and migrate
through – and combined it with a polarized cell motility model [7] to
represent the Trms (Figure 7B). Trm motility parameters were derived from a LN
model constructed in parallel, and our parameter choice will be explained in
more detail with the LN data. These parameter settings lead to motility modes
also observed in intravital liver data: intermittent movement and pausing of
Trms. The goal of our simulation was to predict how much of the liver volume
would be explored by a given number of T cells within a set amount of time;
this same question was addressed by Fernandez-Ruiz _et al._ using a more
abstract mathematical model [26].

As can be expected, we found that liver exploration by the Trm population
increased with the number of exploring cells, and with the amount of
exploration time (Figure 7C). Remarkably, the sinusoid volume explored by an
individual cell increased dramatically as the total Trm population was raized
from 1,000 to 10,000 cells (Figure 7D). Further investigation revealed the
source of this phase transition: cells in our simulations can be in an
unpolarized state, in which they do not move persistently; collisions can
trigger polarization of stationary cells. Since such collisions occur more
frequently at higher densities, this causes a larger fraction of cells to be
motile (Figure 7E). Thus, if Trms can help each other to move, a larger Trm
population would screen the environment more efficiently than expected from
the population size alone.

Notably, the mean speed of Trms peak at $\sim$10.5 μm/min (Figure 7F),
matching average speeds in our lymph node simulations – but only at densities
where the arrest coefficient is close to 0. In two-photon experiments, liver
Trm speed is reported to be significantly lower than T cells in lymph nodes
[27, 28]; we can reproduce this in our simulations with moderate cell density,
which leads to both a realistic slowdown, and higher arrest coefficient. The
simulated liver Trm speed and arrest coefficient closely matches previously
reported values in literature [29], where they were used to quantify observed
characteristic patrolling behavior of Trms.

These liver experiments demonstrated that our system allowed us to combine 3D
imaging data with a relatively realistic model of polarizing cell motility.
The resulting simulation produced an interesting phase transition of Trm
screening efficiency. While we do not claim that such an effect is present in
the real liver – more detailed tuning of the motility model to match the
characteristics of real sinusoidal Trm cells would be out of scope for this
simple proof of concept – this does illustrate how detailed models such as the
CPM generate emergent behavior that is not an immediately obvious consequence
of the model ingredients and can generate interesting biological hypotheses
worthy of further exploration.

### A T cell collective moves robustly in a realistic, organ-scale model of
the lymph node

In LNs, millions of T cells move in search of antigen-presenting cells [21,
30], and finding the cognate antigen is a key early step in adaptive immune
responses [31]. T cell motility in the LN is shaped by the high deformability
of the cells [32] and their interactions with a complex stromal environment,
the fibroblastic reticular cell network [33, 34]. Given the LN’s unique
structure as an organ consisting mainly of moving cells, it presents a
challenging test case for any tissue simulation system. Indeed, many models of
LN T cell dynamics already exist, but they were often limited in either cell
count and simulation volume, or in realism of cell morphology and LN
environment. For example, T cell interactions with a synthetic stromal
structure have been investigated using 3D CPM simulations of thousands of
cells [9], while collective motility at a more realistic scale of hundreds of
thousands of cells has been investigated with non-morphological agent-based
models [35, 36]. We were therefore interested in our system’s ability to
simulate a realistically large collective of millions of realistically shaped
T cells. As in the previous use case, we used microscopy data to initialize –
and parameterize – our model.

Figure 8: Reproducing _in vivo_ T cell motility characteristics in a CPM. We
compare intravital data from Mandl _et al._ [21] to the motion of 50,135 CPM
cells moving in a spherical structure (256μm diameter). (A) A 130x130x110 μm3
scan of the stromal structure (fibroblastic reticular cells) in the
paracortical T cell zone of a murine LN. The images show average intensity
projections of different depth ranges. After segmentation and mirroring, these
data are used to fill a spherical simulation volume. (B) Speeds observed in
simulated data are used to construct a probability distribution (upper panel,
solid), which we compare to the empirical distribution (upper panel, dashed)
using the KL divergence metric. To this end, we bin the speed distributions,
compute the divergence in each bin, and sum the values for each bin (lower
panel). (C) Increased adhesion of T cells to stromal cells causes the speed
distributions to be more right-skewed (top) and match real data more closely
(bottom). (D) Increased adhesion can slow T cell paths non-uniformly by
creating longer and/or more frequent pauses (red) but does not prevent rapid
motion altogether. Two representative tracks are shown. (E) A grid search of
motility parameters for cell polarization reveals which parameter ranges best
fit the intravital data. Both low and high adhesion simulations have similar
fit landscapes; the best fit is found in a high adhesion simulation, although
the difference to the low-adhesion best fit is small (1 simulation per
parameter combination, best fit: red box/red line). CPM parameter settings can
be found in Table 3.

We constructed a LN simulation using the work of Beltman _et al._ [9] as a
starting point. To improve realism, we used an imaged LN section to initialize
the stromal structure (Figure 8A), we used the Act-CPM extension for cell
motility [7], and we fit the two Act motility parameters directly to
intravital imaging data (Figure 8B-E). When fitting these parameters, we aimed
to reproduce the intermittent stop-and-go motion [37] observed in real data
whilst also quantitatively approximate the real speed distribution. Our system
allowed fitting to be performed on realistically large simulations (256μm
diameter simulation volume, containing 50,135 T cells, 20 minutes simulated
time). Crucially, we did not use periodic boundaries which, in dense
simulations, can lead to artificial global synchronisation of motion in the
cell collective.

We used the KL divergence to compare simulated to real speed distributions
(Figure 8B). Increased adhesion of T cells to the stromal structure made the
speed distribution more right-skewed (Figure 8C, top), which was more
consistent with real data (Figure 8C bottom). For most examined parameter
settings, we found the pronounced variations in cell speed along cell
trajectories that are also observed in real data (Figure 8D). When we
performed grid searches of motility parameters at two different stromal
adhesion strengths (Figure 8E), both grid searches had broadly similar fit
landscapes, and the overall lowest KL divergence values were very similar.
Qualitatively, we found more right-skewed and bimodal speed distributions in
high-adhesion than low-adhesion simulations for a subset of the motility
parameters (Supplementary Figure S2). Based on these fitting results, we chose
the best fit motility parameters of the high adhesion setting for our further
simulations.

Figure 9: Mixing of motile T cell collectives in a realistic-scale LN
simulation. (A) Spherical LN simulation at 1mm diameter containing 3,359,752 T
cells and stromal cells (grey; Figure 8A). For clearer visualization, only 10%
of the simulated cells are shown, and red and green cells are equally
parameterized although coloured differently. (B) The simulation is initialized
with two equally parameterized populations each starting at one half of the
sphere. Mixing is quantified by tracking the fraction of cells having at least
one neighbour of the other population. (C) Global mixing is slower in the
high-adhesion simulation (red) than at low adhesion (black), although the
difference is minor. Both simulations reach 50% mixing in ~1h. (D) Cells do
not synchronize their directionality much at longer spatial scales, and this
absence of long-range synchronization is the same for low and high adhesion.
CPM parameter settings can be found in Table 3.

We next investigated organ-scale collective T cell motility by scaling up our
simulated volume to 1mm diameter, now containing 3,359,752 T cells (Figure
9A). While this is not entirely realistic – a real LN has a more complex
compartmental structure, while we simulate only a T cell zone – our main
intention was to determine the spatial range of collective behavior that can
now be investigated with our system. We simulated one virtual day in 1.9 days
of wall clock time – where serial implementations were projected to have taken
a full 1.4 years (Figure 1A).

We used a mixing simulation of two differentially labelled populations as an
easily observable proxy for collective search (Figure 9B). In the simulation,
it takes about 1 hour to reach a state where 50% of the cells are next to at
least one cell of the other population (Figure 9C). While stromal adhesion
influences intra-cellular stop-and-go motility, it has minimal influence on
collective motion: high and low stromal adhesion simulations have broadly
similar mixing speeds. While collective motility is slightly ordered at cell-
scale distances, no ordering persists at organ-scale (Figure 9D).

In summary, we used an organ-scale LN simulation as a challenging test case:
an environment known for complex, large scale collective motility dynamics.
Our system allowed us to integrate of microscopy data directly into our
simulation, and to fit motility parameters directly on two-photon data. We
investigated the behavior of the fitted cell population at organ-scale
(millions of cells). We hope that this will inspire further research on large
cell collectives – integrating various sources of empirical data, at a high
level of realism, and at realistic scale.

## Discussion

We have implemented a parallel version of the CPM that, in some experiments,
was 3 to 4 orders of magnitude faster on GPUs than a state-of-the-art serial
CPM. To achieve this, we designed novel schemes for lattice work splitting and
synchronization, maintaining global state, efficient GPU memory transfer, and
multi-GPU execution. We found that preservation of local lattice update
statistics was sufficient to maintain high-level system dynamics, despite
removing the sequential nature of the original CPM algorithm.

In our work we aimed to ensure that important properties of CPM simulation
trajectories were maintained. Some previous CPM parallelisation efforts
instead focused on equilibrium state distributions as important quantities.
For example, Chen _et al._[17] argue that their parallel CPM maintains certain
equilibrium properties. While such arguments are frequently applied in the
analysis of energy-based models in physics, Voss-Böhme _et al._ showed that
these are less relevant for the CPM because the Metropolis algorithm used in
the CPM drastically alters longtime behavior of the model [38].

While we have shown that the same range of behaviors is present in our
parallel CPM as in a serial CPM in our benchmark scenarios, this would still
need to be carefully checked when re-using simulation parameters that were
originally fitted on a serial implementation. Specifically, researchers might
want to re-fit their simulation parameters when using our system instead of a
serial version – advice that is sensible in general when moving between
different implementations of the same model formalism, especially with complex
models such as the CPM where implementation details might anyway differ
subtly. Furthermore, the observed relation between model kinetics and
performance (Figure 6E) also suggests it might be fruitful to explore more
aggressive parallelisation settings that might yield further speed-up whilst
maintaining desired simulated behavior.

Through accelerating the CPM on accessible GPU hardware, we have increased the
range of feasible simulation sizes – because our scheme fully utilises all GPU
cores, simulation size is now mostly limited by the available GPU memory. As
memory per GPU increases, our scheme will easily scale to larger simulations;
our multi-GPU prototype scheme will allow for further increase of simulation
volume. Our multi-GPU scheme could be further extended to include inter-GPU
communication over MPI [39], allowing our system to be used as a fundamental
building block for a GPU cluster-based version of the CPM – potentially
scaling our 1mm3 simulations further to 1cm3 scale.

Large-scale multicellular simulations are often performed using more coarse-
grained agent-based model environments [35, 36, 14, 15]. Our high-performance
CPM now allows researchers to perform simulations that are morphologically
much more detailed and allow the investigation of new questions that require
taking cell shape into account, such as migration in complex environments or
morphogenesis. We hope that our example simulations of liver and LN
environments will serve as inspiration for many such spatially resolved tissue
and organ simulations in the future. Similar approaches to the one shown here
could be useful for parallelisation of CPM variants or other sequential
lattice-based models [8].

## Tables

Table 1:  Model settings for 2D simulations. For all 2D simulations we use $V_{\mathrm{target}}=2500$, $A_{\mathrm{target}}=600$, $\lambda_{\mathrm{perimeter}}=2$, $J_{\mathrm{cell},\mathrm{ECM}}=10$ and $T=10$. The vector field in the 2D circular chemotaxis simulations (Figure 5A) consists of vectors at a 90 degree angle to radial vectors from points in the field to the simulation center, with all vectors in the field normalized to length 1.  | $\lambda_{\mathrm{volume}}$ | $\lambda_{\mathrm{act}}$ | $\mathrm{Max}_{\mathrm{act}}$ | $\lambda_{\mathrm{chemotaxis}}$
---|---|---|---|---
Fig 1B | 20, 30, 40 | 200, 250, 300, 350, 400 | 50 | —
Fig 5A, Fig 6E | 20 | — | — | 10
Fig 5B | 20 | 50-800, in increments of 50 | 50 | —
Table 2:  Model settings for 3D simulations. For all simulations $V_{\mathrm{target}}=150$, $\lambda_{\mathrm{volume}}=25$, $A_{\mathrm{target}}=1500$, $\lambda_{\mathrm{perimeter}}=0.2$, $J_{\mathrm{cell},\mathrm{cell}}=10$ and $T=20$. The time resolution is 1 second per Monte Carlo step.  | $\mu$ | $\Delta t$ | $\lambda_{\mathrm{persistence}}$ | $J_{\mathrm{cell},\mathrm{ECM}}$
---|---|---|---|---
Fig 1A | — | — | — | 10
Fig 5C, Fig 6D | 0.95 | 15 | 40 | 0
Table 3:  Model settings for liver and lymph node simulations. For all simulations $V_{\mathrm{target}}=150$, $\lambda_{\mathrm{volume}}=25$, $A_{\mathrm{target}}=1500$, $\lambda_{\mathrm{perimeter}}=0.2$, $J_{\mathrm{cell},\mathrm{cell}}=10$, $J_{\mathrm{cell},\mathrm{ECM}}=0$, and $T=20$. The time resolution is 1 second per Monte Carlo step.  | $\lambda_{\mathrm{act}}$ | $\mathrm{Max}_{\mathrm{act}}$ | $J_{\mathrm{cell},\mathrm{FRC}}$
---|---|---|---
Fig 7 | 2600 | 10 | 0
Fig 8B | 2600 | 10 | 0
Fig 8C | 2600 | 10 | -400, -300, -200, -100, 0
Fig 8D | 2600 | 10 | -400, 10
Fig 8E | 200-3000, in increments of 200 | 10, 15, 20, 25, 30, 40, 50, 60, 70 | -3000, 10
Fig 9C | 2600 | 10 | -800, 0
Fig 9D | 2600 | 10 | -800

## Acknowledgments

This research was supported by Young Investigator Grant 10620 from KWF
Kankerbestrijding and HFSP program grant RGP0053/2020 (to JT). JT was also
supported by NWO Vidi grant VI.Vidi.192.084. We thank Franka Buytenhuijs for
giving valuable feedback on this manuscript, and we thank L.W.A. Stuurman and
Pedro Aroca Lara for testing early versions of this simulation system. We also
gratefully acknowledge the support of NVIDIA Corporation with the donation of
the Titan V GPU used for this research. The funders had no role in design,
execution, and interpretation of the research reported in this manuscript.

## Author contributions

SS and JT conceived the study. SS designed the parallel CPM algorithms. SS
performed the research and analysed the data. SS and JT wrote the manuscript.
JT supervized the project. SM provided critical feedback and provided the
liver and lymph node imaging data. SD performed image acquisition of the liver
and LN data. All authors critically revized the manuscript for intellectual
content.

## Declaration of interests

The authors declare no competing interests.

## Materials and Methods

### Cellular Potts model and model variants

The CPM is a lattice based formalism for simulating biological cells.
Simulations proceeds with a Markov chain Monte Carlo process, where a single
position at a time can transition to one of the states of its neighbours. The
transition probability is based on the change of a Hamiltonian energy
function. The original CPM has

$\displaystyle\mathcal{H}=\sum_{i,j}J_{\tau(\sigma_{i}),\tau(\sigma_{j})}(1-\delta_{\sigma_{i},\sigma_{j}})+\sum_{\sigma}\lambda_{\mathrm{volume}}(\tau_{\sigma})(V(\sigma)-V_{\mathrm{target}}(\tau_{\sigma}))^{2}$
(1)

as its Hamiltonian, with $i,j$ being possible neighbour locations in the
lattice, $\sigma_{i}$ the cell present at location $i$, $\tau$ is the type of
a particular cell, $\delta$ is the Kronecker delta, $V(\sigma)$ is the number
of lattice sites a cell $\sigma$ occupies, and
$\lambda_{\mathrm{volume}}(\tau)$ how strongly a volume constraint contributes
to the Hamiltonian for a particular cell type. In all our simulations a Moore
neighbourhood was used. Transitions stochastically succeed based on the change
in energy $\Delta\mathcal{H}$, using a Boltzmann probability distribution

$\displaystyle P(\Delta\mathcal{H})=\begin{cases}1&\text{if
$\Delta\mathcal{H}<0$,}\\\ e^{-\Delta\mathcal{H}/T}&\text{if
$\Delta\mathcal{H}\geq 0.$}\end{cases}$ (2)

By adding extra terms to the Hamiltonian, constraints can be added to cells
and new types of behavior introduced to the model. When cell shape needs to be
parameterized beyond just cell size, a perimeter constraint [40] is often
introduced,

$\displaystyle\mathcal{H}_{\mathrm{perimeter}}=\sum_{\sigma}\lambda_{\mathrm{perimeter}}(\tau_{\sigma})(A(\sigma)-A_{\mathrm{target}}(\tau_{\sigma}))^{2}.$
(3)

Here $A(\sigma)$ represents the number of interfaces between lattice sites of
a cell with sites not part of the same cell.

We also implemented several motility models. The first allows a vector field
to bias the energy change $\Delta\mathcal{H}$ during a copy attempt

$\displaystyle\Delta\mathcal{H}_{\mathrm{chemotaxis}}=\lambda_{\mathrm{chemotaxis}}(\tau_{\sigma})\hat{x}_{i\rightarrow
j}\cdot\vec{C_{i}}.$ (4)

Here $\hat{x}_{i\rightarrow j}$ is the unit vector pointing in the direction
of a copy attempt from position $i$ to $j$, and $C_{i}$ the value of a vector
field at position $i$. In this paper, $\vec{C_{i}}$ is always a field of unit
vectors.

For the model where cells exhibit persistent movement based on one vector per
cell, we add

$\displaystyle\Delta\mathcal{H}_{\mathrm{persistence}}=\lambda_{\mathrm{persistence}}(\tau_{\sigma})\hat{x}_{i\rightarrow
j}\cdot\hat{p}_{\sigma,t}$ (5)

to the Hamiltonian delta. The vector $\hat{p}_{\sigma,t}$ is the persistence
direction of cell $\sigma$ at time $t$, which is updated iteratively at every
MCS with

$\displaystyle\vec{p}_{\sigma,t}=(1-\mu_{\tau})\frac{\vec{X}_{\sigma,t}-\vec{X}_{\sigma,t-\Delta
t}}{\lVert\vec{X}_{\sigma,t}-\vec{X}_{\sigma,t-\Delta
t}\rVert}+\mu_{\tau}\hat{p}_{\sigma,t-1},$ (6)
$\displaystyle\hat{p}_{\sigma,t}=\frac{\vec{p}_{\sigma,t}}{\lVert\vec{p}_{\sigma,t}\rVert}.$
(7)

Here $\mu_{\tau}$ is the strength of directional persistence for cell type
$\tau$, $\vec{X}_{\sigma,t}$ is the centroid location of cell $\sigma$ at time
$t$ and $\Delta t$ is the number of MCSs we look back for to determine the
current direction of a cell.

Finally, we implement the Act model [7] that allows for localized
polarization. For each position $i$ we track an activity value each time said
location is updated, here labelled as $\mathrm{Act}(i)$. This localized
‘memory’ is initialized at $\mathrm{Max}_{\mathrm{act}}$, and decays by 1 each
timestep until it reaches 0. The Hamiltonian delta is adapted by

$\displaystyle\Delta\mathcal{H}_{\mathrm{act}}=\frac{\lambda_{\mathrm{act}}(\tau_{\sigma})}{\mathrm{Max}_{\mathrm{act}}(\tau_{\sigma})}(\mathrm{GM}_{\mathrm{act}}(i)-\mathrm{GM}_{\mathrm{act}}(j)),$
(8)

making copies that originate from areas of a cell with a high activity value
more likely. The activity is calculated for a full Moore neighbourhood using a
geometric mean

$\displaystyle\mathrm{GM}_{\mathrm{act}}(i)=(\prod_{n\in
N(i)}\mathrm{Act}(n))^{1/\lvert N(i)\rvert}.$ (9)

$N(i)$ are all neighbours of pixel or voxel $i$.

### Optimized CPU implementation of CPM

Besides a high-performance GPU implementation, we also developed our own CPU
implementation. This served both as a reference implementation for model
behavior and internal implementation details for us, and as a reference
performance benchmark. We strived to reduce the major performance bottlenecks
we could see in profiling. The implementation tracks lattice sites on cell
boundaries, and only tries copy attempts involving such positions; a common
optimisation, which is also implemented in Morpheus. The tracking
implementation uses a high-performance hashmap implementation [41]. We use a
custom high-performance pseudorandom number generator, xoshiro256 [42]. The
lattice data is stored with a Z-order space filling curve, allowing data that
is close in the simulation space also to be close in the linear memory space,
maximizing cache efficiency. Memory addresses are calculated using intrinsic
bit interleaving instructions from the x86 bit manipulation instruction set
extension.

###  GPU implementation details

Our CPM was implemented using CUDA; concepts equivalent to ones described here
exist in non-CUDA GPU hardware. A GPU consists of multiple streaming
multiprocessors (SMs). Each SM runs multiple threads in parallel; contrasting
with CPU threads, threads on the same SM execute in lockstep, in a Single
Instruction Multiple Threads approach. A single GPU process – called a kernel
– can consist of up to millions of threads.

GPUs have relatively slow global memory available to all threads, and much
faster shared memory local to SMs. Global memory transfer speed is often the
main performance bottleneck for GPU programs. Threads are bundled in thread
blocks (TBs); each TB is bound to a specific SM. By caching data on their SM
in shared memory, TBs can greatly reduce global memory transfer bottlenecks.

Memory transfer bottlenecks can be further reduced by two strategies: read
coalescing and high SM occupancy. Read coalescing has sequential threads in a
TB do contiguous reads of global memory; contiguous reads are needed to fully
utilise the memory bus bandwidth. High occupancy is achieved when a high
number of threads are scheduled to execute per SM; the SM scheduler will
switch out threads when they are blocked on memory transfers, ensuring a high
number of in-flight memory instructions at a time.

We propose a scheme of memory transfer, local SM execution and synchronization
(Figure 2D). We focus on caching just the grid of cell ids, as this data is
always needed, and can be used to quickly reject copy attempts internal to a
cell. The simulation grid is divided in subsections, with regions larger than
earlier checkerboard subdivisions; each subsection is assigned to a TB,
subsection size is chosen so active regions fit in shared memory. Each TB
determines which checkerboard sections are active in its region; active
sections are cached in shared memory using coalesced reads.

After caching – intended to speed up reads – the simulation is run, and
successful copy attempts write new values to both shared and global memory.
Alternatively, changes could be cached in shared memory and written back
coalesced to global memory; in our testing, an immediate shared/global write
yielded better performance. If a TB is done with simulation execution, it
immediately starts a memory request for the next active region, ensuring
maximum occupancy. All threads in the kernel are synchronized after new memory
transfers are finished, ensuring maximum overlap of code execution and memory
transfer.

CPM literature offers a wide range of model variants and extensions, of which
we implement a modest subset. This comes with potential performance issues:
even when extensions are disabled through branching statements, on GPUs such
branching code still incurs a significant performance loss compared to CPUs.
We therefore pre-generated specialized versions of the kernel for a range of
different extensions.

We use a custom high-performance random number generator on the GPU [43].

### Composite variable tracking

Several CPM extensions introduce additional global variables. Some of these
variable need to be kept in-sync with each other – one variable being slightly
out of date with others is not tolerated. One example is tracking of a center-
of-mass centroid; for 2D this can be tracked as a sum of x and y values, and
corresponding count of lattice locations occupied (Figure 2C). These are
integer values: tracking centroid locations as floating-point would cause
values to drift due to rounding errors. The centroid value is calculated by

$\displaystyle\\{x,y\\}=\frac{1}{\lvert\bullet\rvert}\times\\{x_{\mathrm{sum}},y_{\mathrm{sum}}\\}.$
(10)

Online updates change both a cell’s centroid position and size. When a new
lattice position $\\{x_{\mathrm{position}},y_{\mathrm{position}}\\}$ is added
to a cell, its centroid’s $x$ value can be updated with

$\displaystyle
x_{\mathrm{new}}=\frac{x_{\mathrm{old}}\times\lvert\bullet\rvert_{\mathrm{old}}+x_{\mathrm{position}}}{\lvert\bullet\rvert_{\mathrm{new}}}.$
(11)

Cell sizes are stored interleaved with sum values; sizes and sums are both
stored as 32-bit unsigned integers. The interleaved storage enables use of
64-bit operations to update both: a copy-and-swap (CAS) operation is used,
which atomically updates sum and cell size values simultaneously. This scheme
also handles periodic boundaries, by adding or subtracting
$\mathrm{dimension}\times\lvert\bullet\rvert$ to the sum values before
executing the CAS as needed, to keep centroids within simulation bounds.

### Multi-GPU execution

We devized a scheme to execute our CPM on multiple GPUs (Figure 2E); this
scheme is intended as a prototype, and as such not as optimized as other parts
of our CPM. Simulation space is cut into equal regions, each region is
assigned to a GPU. A region’s lattice data, such as the cell ids, are sent to
each GPU, with an additional border layer. Global variables such as cell sizes
are tracked on each GPU as deltas of the value since last inter-GPU data
synchronisation. Every four MCSs, these deltas are gathered from all GPUs,
deltas are integrated, and updated values distributed to each GPU. Multi-GPU
speedup is mostly hampered by frequent or slow inter-GPU communication; inter-
GPU synchronization after every MCS negates most speedup, and a faster GPU-to-
GPU interlink would likely increase performance.

### Hardware

The CPU-based simulations, using our custom optimized CPU implementation and
Morpheus 2.2.6 (the last non-parallelized version of Morpheus), and most GPU
simulations were run on a system with a Nvidia Titan V, an AMD 2950x CPU and
128 GB of RAM at 2133 MHz. The only exceptions to this are the multi-GPU
simulations (Figure 2E), which were run on a system with four Nvidia Geforce
RTX 2080 Ti GPUs, an Intel i9-9820x CPU at 3.3 GHz, and 128 GB of RAM at 2133
MHz.

### Imaging of thick lymph node sections

Naïve inguinal LNs were harvested from CCL19-EYFP mice [44] and fixed with 4%
paraformaldehyde (PFA, Electron Microscopy Sciences) for 1 hour at 4°C.
Following a brief wash with PBS, LNs were embedded in 2% agarose and 250-300
µm tissue sections were generated using a vibratome (Leica). LN sections were
subsequently immersed in a Special Wash Buffer (Baxter water, 0.2% BSA, 10%
10X PBS, 0.2% Tween-20, 0.2% Triton X-100, 0.2% 10% SDS) at 4°C overnight,
followed by staining with a rabbit anti-GFP Ab (Thermo Scientific, 1:100,
diluted in blocking buffer made with Baxter water, 0.2% BSA, 10% 10X PBS, 0.1%
Tween-20, 0.2% 10% SDS) at 4°C for three days. After incubation, LN sections
were thoroughly washed in PBS and mounted in 85% glycerol in PBS with glass
coverslips. Thick LN sections were imaged using a LSM880 AiryscanFast confocal
microscope (Carl Zeiss) with a x63/1.4NA oil immersion lens and images
processed with ZEN 2010 software (Carl Zeiss).

### Imaging of liver sinusoids

For whole-liver imaging, AlexaFluor 647 conjugated anti-CD31 rat monoclonal
antibody (Biolegend) was intravenously injected into wild-type C57BL/6 mice at
5ug for 15 minutes prior to harvesting the liver. This was followed by
overnight 4% PFA fixation at 4°C. Livers were then thoroughly washed with PBS
and small cubes of liver were prepared before immersion in FUnGI [45] for two
days before mounting in 85% glycerol in PBS. Thick liver sections were imaged
using FVMPE-RS multiphoton microscope (Olympus) with a x20/1.05 NA water
immersion lens.

## Data availability

The liver imaging data and the segmented version are available on Zenodo [46].
The lymph node imaging data and the segmented version are available on Zenodo
[47]. Source data and code for the figures will be published together with the
final, archival version of this manuscript. In the meantime, these data and
code are available upon request to the corresponding author and are made
available to reviewers.

## Code availability

GPU CPM simulation code is available on GitHub at
https://github.com/shabaz/gpu-cpm and is made available under the GNU General
Public License, version 3. CPU CPM simulation code is available on GitHub at
https://github.com/shabaz/cpm and is made available under the GNU General
Public License, version 3.

## References

  * [1] Graner Fm and Glazier JA.  Simulation of biological cell sorting using a two-dimensional extended potts model.  Phys. Rev. Lett., 69:2013–2016, 1992.  doi:10.1103/PhysRevLett.69.2013.
  * [2] Szabó A and Merks RM.  Cellular potts modeling of tumor growth, tumor invasion, and tumor evolution.  Frontiers in Oncology, 3:87, 2013.  doi:10.3389/fonc.2013.00087.
  * [3] Merks RMH, Perryn ED, Shirinifard A, and Glazier JA.  Contact-inhibited chemotaxis in de novo and sprouting blood-vessel growth.  PLOS Computational Biology, 4(9):1–16, 2008.  doi:10.1371/journal.pcbi.1000163.
  * [4] Scianna M.  An extended cellular potts model analyzing a wound healing assay.  Computers in Biology and Medicine, 62:33–54, 2015.  doi:https://doi.org/10.1016/j.compbiomed.2015.04.009.
  * [5] Hirashima T, Rens EG, and Merks RMH.  Cellular potts modeling of complex multicellular behaviors in tissue morphogenesis.  Development, Growth & Differentiation, 59(5):329–339, 2017.  doi:https://doi.org/10.1111/dgd.12358.
  * [6] Marée AFM and Hogeweg P.  How amoeboids self-organize into a fruiting body: Multicellular coordination in dictyostelium discoideum.  Proceedings of the National Academy of Sciences, 98(7):3879–3883, 2001.  doi:10.1073/pnas.061535198.
  * [7] Niculescu I, Textor J, and de Boer RJ.  Crawling and gliding: A computational model for shape-driven cell migration.  PLOS Computational Biology, 11(10):1–22, 2015.  doi:10.1371/journal.pcbi.1004280.
  * [8] Thüroff F, Goychuk A, Reiter M, and Frey E.  Bridging the gap between single-cell migration and collective dynamics.  eLife, 8:e46842, 2019.  doi:10.7554/eLife.46842.
  * [9] Beltman JB, Marée AF, Lynch JN, Miller MJ, and de Boer RJ.  Lymph node topology dictates T cell migration behavior .  Journal of Experimental Medicine, 204(4):771–780, 2007.  doi:10.1084/jem.20061278.
  * [10] Vega R, Carretero M, Travasso RDM, and Bonilla LL.  Notch signaling and taxis mechanisms regulate early stage angiogenesis: A mathematical and computational model.  PLOS Computational Biology, 16(1):1–31, 2020.  doi:10.1371/journal.pcbi.1006919.
  * [11] Rens EG and Merks RM.  Cell shape and durotaxis explained from cell-extracellular matrix forces and focal adhesion dynamics.  iScience, 23(9):101488, 2020.  doi:https://doi.org/10.1016/j.isci.2020.101488.
  * [12] Marco Scianna KW Luigi Preziosi.  A cellular potts model simulating cell migration on and in matrix environments.  Mathematical Biosciences & Engineering, 10(1):235–261, 2013.
  * [13] Vroomans RMA, Marée AFM, de Boer RJ, and Beltman JB.  Chemotactic migration of t cells towards dendritic cells promotes the detection of rare antigens.  PLOS Computational Biology, 8(11):1–13, 2012.  doi:10.1371/journal.pcbi.1002763.
  * [14] Van Liedekerke P, Palm MM, Jagiella N, and Drasdo D.  Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results, volume 2.  Springer International Publishing, 2015.
  * [15] Metzcar J, Wang Y, Heiland R, and Macklin P.  A review of cell-based computational modeling in cancer biology.  JCO Clinical Cancer Informatics, (3):1–13, 2019.  doi:10.1200/CCI.18.00069.  PMID: 30715927.
  * [16] Osborne JM, Fletcher AG, Pitt-Francis JM, Maini PK, and Gavaghan DJ.  Comparing individual-based approaches to modelling the self-organization of multicellular tissues.  PLOS Computational Biology, 13(2):1–34, 2017.  doi:10.1371/journal.pcbi.1005387.
  * [17] Chen N, Glazier JA, Izaguirre JA, and Alber MS.  A parallel implementation of the cellular potts model for simulation of cell-based morphogenesis.  Computer Physics Communications, 176(11):670 – 681, 2007.  doi:https://doi.org/10.1016/j.cpc.2007.03.007.
  * [18] Tapia JJ and D’Souza RM.  Parallelizing the cellular potts model on graphics processing units.  Computer Physics Communications, 182(4):857 – 865, 2011.  doi:https://doi.org/10.1016/j.cpc.2010.12.011.
  * [19] Starruß J, de Back W, Brusch L, and Deutsch A.  Morpheus: a user-friendly modeling environment for multiscale and multicellular systems biology.  Bioinformatics, 30(9):1331–1332, 2014.  doi:10.1093/bioinformatics/btt772.
  * [20] Shirinifard A, Gens JS, Zaitlen BL, Popławski NJ, Swat M, and Glazier JA.  3d multi-cell simulation of tumor growth and angiogenesis.  PLOS ONE, 4(10):1–11, 2009.  doi:10.1371/journal.pone.0007190.
  * [21] Mandl JN, Liou R, Klauschen F, Vrisekoop N, Monteiro JP, Yates AJ, Huang AY, and Germain RN.  Quantification of lymph node transit times reveals differences in antigen surveillance strategies of naïve cd4+ and cd8+ t cells.  Proceedings of the National Academy of Sciences, 109(44):18036–18041, 2012.  doi:10.1073/pnas.1211717109.
  * [22] Szabó B, Szöllösi GJ, Gönci B, Jurányi Z, Selmeczi D, and Vicsek T.  Phase transition in the collective migration of tissue cells: Experiment and model.  Phys. Rev. E, 74:061908, 2006.  doi:10.1103/PhysRevE.74.061908.
  * [23] Hugues S, Fetler L, Bonifaz L, Helft J, Amblard F, and Amigorena S.  Distinct t cell dynamics in lymph nodes during the induction of tolerance and immunity.  Nature Immunology, 5(12):1235–1242, 2004.  doi:10.1038/ni1134.
  * [24] Masopust D and Soerens AG.  Tissue-resident t cells and other resident leukocytes.  Annual Review of Immunology, 37(1):521–546, 2019.  doi:10.1146/annurev-immunol-042617-053214.  PMID: 30726153.
  * [25] Lefebvre MN and Harty JT.  You shall not pass: Memory cd8 t cells in liver-stage malaria.  Trends in Parasitology, 36(2):147–157, 2020.  doi:10.1016/j.pt.2019.11.004.
  * [26] Fernandez-Ruiz D, Ng WY, Holz LE, Ma JZ, Zaid A, Wong YC, Lau LS, Mollard V, Cozijnsen A, Collins N, Li J, Davey GM, Kato Y, Devi S, Skandari R, Pauley M, Manton JH, Godfrey DI, Braun A, Tay SS, Tan PS, Bowen DG, Koch-Nolte F, Rissiek B, Carbone FR, Crabb BS, Lahoud M, Cockburn IA, Mueller SN, Bertolino P, McFadden GI, Caminschi I, and Heath WR.  Liver-resident memory cd8+ t cells form a front-line defense against malaria liver-stage infection.  Immunity, 45(4):889–902, 2016.  doi:10.1016/j.immuni.2016.08.011.
  * [27] Holz LE, Prier JE, Freestone D, Steiner TM, English K, Johnson DN, Mollard V, Cozijnsen A, Davey GM, Godfrey DI, et al.  Cd8+ t cell activation leads to constitutive formation of liver tissue-resident memory t cells that seed a large and flexible niche in the liver.  Cell reports, 25(1):68–79, 2018.
  * [28] Hor JL, Whitney PG, Zaid A, Brooks AG, Heath WR, and Mueller SN.  Spatiotemporally distinct interactions with dendritic cell subsets facilitates cd4+ and cd8+ t cell activation to localized viral infection.  Immunity, 43(3):554–565, 2015.
  * [29] McNamara H, Cai Y, Wagle M, Sontani Y, Roots C, Miosge L, O’Connor J, Sutton H, Ganusov V, Heath WR, et al.  Up-regulation of lfa-1 allows liver-resident memory t cells to patrol and remain in the hepatic sinusoids.  Science immunology, 2(9):eaaj1996, 2017.
  * [30] Textor J, Henrickson SE, Mandl JN, von Andrian UH, Westermann J, de Boer RJ, and Beltman JB.  Random Migration and Signal Integration Promote Rapid and Robust T Cell Recruitment.  PLoS Computational Biology, 10(8), 2014.  doi:10.1371/journal.pcbi.1003752.
  * [31] Gasteiger G, Ataide M, and Kastenmüller W.  Lymph node - An organ for T-cell activation and pathogen defense.  Immunological Reviews, 271(1):200–220, 2016.  doi:10.1111/imr.12399.
  * [32] Dupré L, Houmadi R, Tang C, and Rey-Barroso J.  T lymphocyte migration: An action movie starring the actin and associated actors.  Frontiers in Immunology, 6:586, 2015.  doi:10.3389/fimmu.2015.00586.
  * [33] Bajénoff M, Egen JG, Koo LY, Laugier JP, Brau F, Glaichenhaus N, and Germain RN.  Stromal cell networks regulate lymphocyte entry, migration, and territoriality in lymph nodes.  Immunity, 25(6):989–1001, 2006.
  * [34] Textor J, Mandl JN, and de Boer RJ.  The reticular cell network: A robust backbone for immune responses.  PLOS Biology, 14(10):1–6, 2016.  doi:10.1371/journal.pbio.2000827.
  * [35] Bogle G and Dunbar PR.  Agent-based simulation of T-cell activation and proliferation within a lymph node.  Immunology and Cell Biology, 88(2):172–179, 2010.  doi:10.1038/icb.2009.78.
  * [36] Gong C, Mattila JT, Miller M, Flynn JAL, Linderman JJ, and Kirschner D.  Predicting lymph node output efficiency using systems biology.  Journal of Theoretical Biology, 335:169–184, 2013.  doi:10.1016/j.jtbi.2013.06.016.
  * [37] Miller MJ, Wei SH, Cahalan MD, and Parker I.  Autonomous t cell trafficking examined in vivo with intravital two-photon microscopy.  Proceedings of the National Academy of Sciences, 100(5):2604–2609, 2003.
  * [38] Voss-Böhme A.  Multi-scale modeling in morphogenesis: A critical analysis of the cellular potts model.  PLOS ONE, 7(9):1–14, 2012.  doi:10.1371/journal.pone.0042852.
  * [39] Gropp W, Lusk E, and Skjellum A.  Using MPI: portable parallel programming with the message-passing interface, volume 1.  MIT press, 1999.
  * [40] Ouchi NB, Glazier JA, Rieu JP, Upadhyaya A, and Sawada Y.  Improving the realism of the cellular potts model in simulations of biological cells.  Physica A: Statistical Mechanics and its Applications, 329(3):451 – 458, 2003.  doi:https://doi.org/10.1016/S0378-4371(03)00574-0.
  * [41] Skarupke M.  A new fast hash table in response to google’s new fast hash table.  Available from https://probablydance.com/2018/05/28/a-new-fast-hash-table-in-response-to-googles-new-fast-hash-table/, 2018\.
  * [42] Blackman D and Vigna S.  Scrambled linear pseudorandom number generators.  ACM Transactions on Mathematical Software, 47(4):1–32, 2021\.  doi:10.1145/3460772.
  * [43] Langdon WB.  A fast high quality pseudo random number generator for nvidia CUDA.  In Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference: Late Breaking Papers, GECCO ’09, page 2511–2514, New York, NY, USA, 2009. Association for Computing Machinery.
  * [44] Chai Q, Onder L, Scandella E, Gil-Cruz C, Perez-Shibayama C, Cupovic J, Danuser R, Sparwasser T, Luther SA, Thiel V, Rülicke T, Stein JV, Hehlgans T, and Ludewig B.  Maturation of lymph node fibroblastic reticular cells from myofibroblastic precursors is critical for antiviral immunity.  Immunity, 38(5):1013–1024, 2013.  doi:10.1016/j.immuni.2013.03.012.
  * [45] Rios AC, Capaldo BD, Vaillant F, Pal B, van Ineveld R, Dawson CA, Chen Y, Nolan E, Fu NY, Jackling FC, Devi S, Clouston D, Whitehead L, Smyth GK, Mueller SN, Lindeman GJ, and Visvader JE.  Intraclonal plasticity in mammary tumors revealed through large-scale single-cell resolution 3d imaging.  Cancer Cell, 35(4):618–632.e6, 2019.  doi:10.1016/j.ccell.2019.02.010.
  * [46] Devi S, Mueller S, Sultan S, and Textor J.  Liver sinusoids.  Zenodo, 2023.  doi:10.5281/zenodo.8273333.
  * [47] Devi S, Mueller S, Sultan S, and Textor J.  FRC network in cleared lymph node.  Zenodo, 2023.  doi:10.5281/zenodo.8273496.

## List of abbreviations

2D

    two-dimensional
3D

    three-dimensional
CPM

    Cellular Potts Model
GPU

    Graphical Processing Unit
CPU

    Central Processing Unit
MCS

    Monte Carlo Step
LN

    Lymph Node
KL

    Kullback-Leibler

## Supplementary Figures

Supplementary Figure S1:  Checkerboard parallelisation scheme does not have
major influence on sorting model accuracy. We tested sorting simulations with
atomic updates, using different options of parallelisation schemes, varying
space assigned per thread and switching frequency of active checkerboard
areas. Accuracy is measured as the root mean square error of the GPU
simulation runs with 5 CPU reference runs and 5 GPU runs per parallelisation
setting (median measurement marked with line), measuring
homogeneous/heterogeneous cell border interfaces over time at t=0, 1, 10, 100,
1000, 10000, and 100000 seen in Figure 4. Absolute error values remain low,
and there are no pronounced differences between errors for the different GPU
settings (two-way ANOVA, p=0.43).

Supplementary Figure S2:  Cell speed distributions obtained when fitting LN
motility data. Each plot corresponds to one cell in the grid search shown in
Figure 8E. Polarization decay refers to the parameter
$\mathrm{Max}_{\mathrm{act}}$; polarization strength refers to the parameter
$\lambda_{\mathrm{act}}$.

  *[CPM]: Cellular Potts Model
  *[3D]: three-dimensional
  *[CPMs]: Cellular Potts Model
  *[GPU]: Graphical Processing Unit
  *[2D]: two-dimensional
  *[GPUs]: Graphical Processing Unit
  *[CPU]: Central Processing Unit
  *[MCS]: Monte Carlo Step
  *[LN]: Lymph Node
  *[LNs]: Lymph Node
  *[MCSs]: Monte Carlo Step
  *[KL]: Kullback-Leibler


# Latent Aspect Detection from Online Unsolicited Customer Reviews

Mohammad Forouhesh
Amirkabir University of Technology, Iran
<EMAIL_ADDRESS>
&Arash Mansouri
Press’nXPress, Canada
<EMAIL_ADDRESS>
Hossein Fani
University of Windsor, Canada
<EMAIL_ADDRESS>

###### Abstract

Within the context of review analytics, aspects are the features of products
and services at which customers target their opinions and sentiments. Aspect
detection helps product owners and service providers to identify shortcomings
and prioritize customers’ needs, and hence, maintain revenues and mitigate
customer churn. Existing methods focus on detecting the surface form of an
aspect by training supervised learning methods that fall short when aspects
are latent in reviews. In this paper, we propose an unsupervised method to
extract latent occurrences of aspects. Specifically, we assume that a customer
undergoes a two-stage hypothetical generative process when writing a review:
(1) deciding on an aspect amongst the set of aspects available for the product
or service, and (2) writing the opinion words that are more interrelated to
the chosen aspect from the set of all words available in a language. We employ
latent Dirichlet allocation to learn the latent aspects distributions for
generating the reviews. Experimental results on benchmark datasets show that
our proposed method is able to improve the state of the art when the aspects
are latent with no surface form in reviews.

##  1 Introduction

The key characteristic of customers’ opinions is what they target in order to
convey a sentiment, referred to as aspect (variously known as features,
properties, targets, attributes, or components). For example, in the review
‘advertise free wifi but charge you for that,’ the aspect is ‘advertisement’
toward which an opinion ‘misinformation’ is expressed with ‘negative’
sentiment. Detection of aspects is a crucial step in the analysis of
customers’ reviews for products or services in e-commerce and social platforms
as it helps product owners and service providers to find shortcomings and
prioritize the improvements according to customers needs, and hence,
maintaining revenues and decrease customer churn Napitu et al. (2017). As
such, aspect detection along with customers’ opinion mining and sentiment
analysis has long gained much attention Cortis and Davis (2021).

Reviews: [latent] vs. explicit aspects
---
-were given table far from river [management]
-the girl [staff] at the front desk was really nice
-they advertise free wifi but charge you for that
-menu is extensive and a bar with live music
Table 1: Reviews for restaurants from Google reviews.

Existing methods to detect aspects in reviews and their respective sentiment
with the opinions fall into three categories based on the level of human
supervision: i) rule-based, ii) supervised, and iii) unsupervised methods.
Rule-based methods use the co-occurrence association rule mining approach to
match aspects with explicit surface forms (words) Hai et al. (2011); Zeng and
Li (2013). However, such methods are not scalable to an overwhelming number of
combinations in reviews. Supervised methods employ supervised machine learning
techniques on labeled datasets where aspects’ surface forms are manually
annotated in reviews Peng et al. (2020); Zhang et al. (2020). To eschew
reliance on labeled data whose curation is time-consuming, prone to human
biases, and has high manpower cost, unsupervised methods have been adopted He
et al. (2017); Zong et al. (2021). Proposed unsupervised methods, however,
assume that aspects are expressed explicitly in customers’ reviews and forego
the latent occurrences of aspects in reviews. Indeed, when writing an online
unsolicited review, a customer may write about her opinion and overall rating
while overlooking mentioning the aspects’ surface form for a product or
service due to their trivialities or being common knowledge. In Table 1, some
samples from Google reviews have been shown. For instance, the restaurant
review ‘we were given table far from river view’ entails one or several
aspects ‘management’ or ‘staff’ that are latent and not explicitly mentioned
(do not have surface form) due to the background knowledge on restaurants
common services. Indeed, 35% of reviews on restaurants and electronics include
latent aspects Cai et al. (2021); Xu et al. (2015). The need for latent
aspects detection is even more pressing in online social review platforms
where reviews are short, noisy, informal, and mostly rely on social background
knowledge and context.

To this end, in this work, we consider the task of latent aspect detection.
Specifically, given a collection of reviews with no human supervision (no
labels for aspects, opinion, nor sentiment and whether of their explicit or
latent appearances), we aim at finding latent (hidden) aspects in customer
reviews. Although Marrese-Taylor et al. (2013); Poria et al. (2014); Chen and
Chen (2016); Wan et al. (2020) have attempted to solve the latent aspect
detection problem, they all focus on supervised methods on labeled datasets
provided, e.g., SemEval-2016 Pontiki et al. (2016). In this work, unlike the
existing studies, we have no direct access to the aspects’ annotations.

Inspired by LDA Blei et al. (2003), we propose to model aspects as latent
variables. Previous studies have extensively shown the effectiveness of the
LDA for unsupervised learning of the latent representation of topics Cheng et
al. (2014), social users Sun et al. (2017), and communities Xianghua et al.
(2013); Li et al. (2016), to name but a few. We are not, however, aware of
work on latent aspect modeling except for that of a variational autoencoder by
Fei et al. (2021) where a review $r$ and its sentences $s_{r}$ were encoded to
two sets of latent variables that follow Gaussian distribution: i) coarse-
grain review-level ${z_{r}}$, and ii) fine-grain sentence-level $z_{s_{r}}$,
followed by Transformer-based hierarchical classification for the task of
review’s sentiment classification. In their work, Fei et al. (2021) assume
that a review is formed by a bag of aspect and opinion words that are
explicitly mentioned while filtering out other words using opinion lexicons.
Successful for the task of sentiment classification, Fei et al. (2021)’s
unsupervised method detects explicit aspects and overlooks their latent
occurrences. However, our main goal is to propose an unsupervised method to
detect latent aspects when they are hidden with no surface form.

##  2 Proposed Method

Let $\mathcal{D}$ be the set of all unique words in a language. Given a review
$r_{p}=\\{w_{i}\\}\subseteq\mathcal{D}$ as a bag of its constituent words
about a product or a service $p$ concerning an aspect $a\in
A_{p}\subseteq\mathcal{D}$ which is not necessarily in $r_{p}$, that is,
$a\notin r_{p}$, the latent aspect detection aims at identifying aspect $a$ in
$r_{p}$. $A_{p}$ is the set of all aspects of the product or service $p$ and
$R_{p}$ is the set of all $p$’s reviews.

We propose to model aspects as latent variables. We employ latent Dirichlet
allocation to learn the latent aspects distributions for generating the
reviews. Specifically, we assume that a customer undergoes a two-stage
hypothetical generative process to write a review about an aspect of a product
or service: 1) deciding on an aspect $a$ of the product or service amongst the
set of aspects available for the product or service ${A_{p}}$ but with a
different probability, and 2) selecting the words from the set $\mathcal{D}$
of words in a language but with a different probability that is more
interrelated to $a$ to express her opinion about $a$. In this generative
process, a review is about all aspects of the product or service according to
a probability distribution, i.e., $r_{p}\sim\theta_{r_{p}}$, and an aspect $a$
is not an explicit word but a latent probability distribution of all words
with respect to $a$, i.e., $a\sim\varphi_{a}$, where ${\theta}$ and
${\varphi}$ follow Dirichlet probability distribution. Formally, the
generative process for a set of reviews $R_{p}=\\{r_{p}\\}$ is:

  1. 1.

For each $r_{p}\in R_{p}$, draw $\theta_{r_{p}}\sim\text{Dir}(\bm{\alpha})$,
where $\bm{\alpha}$ is a ${|A_{p}|}$-dimensional vector of positive reals
which typically is sparse since a review is assumed to be about one or a few
aspects of the product not all.

  2. 2.

For each $a\in A_{p}$, draw ${\varphi_{a}\sim\text{Dir}(\bm{\beta})}$, where
${\bm{\beta}}$ is a ${|\mathcal{D}|}$-dimensional vector of positive reals.
${\bm{\beta}}$ is also assumed to be sparse since an aspect is interrelated to
a few words and not to all words.

  3. 3.

Given ${r_{p}\in R_{p}}$ along with its $\theta_{r_{p}}$:
3.1. Draw an aspect ${a\sim\text{Multinomial}(\theta_{r_{p}})}$
3.2. For each of the words ${w_{i}}$ in ${r_{p}}$, draw ${w_{i}\sim}$
$\text{Multinomial}({{\varphi}_{a}})$

Given the observed review set of a product ${R_{p}}$, our task is to infer
$\theta_{r_{p}}$ for all reviews ${r_{p}\in R_{p}}$ and ${\varphi_{a}}$ for
all aspects ${a\in A_{p}}$ of the product and estimation of the
${\bm{\alpha}}$ and ${\bm{\beta}}$ parameters using Gibbs sampling Xiao and
Stibor (2010)111Our codebase and PxP dataset have been made publicly available
at https://github.com/MohammadForouhesh/latent-aspect-detection.

##  3 Experiment

We show the empirical performance of our proposed method for $p=restaurant$.
We trained our proposed model on datasets scraped from Google reviews of
restaurants1 (PxP) across North America, and evaluate it on SemEval datasets
Pontiki et al. (2014, 2015, 2016) as the golden standards. Table LABEL:stats
shows the detailed statistics of training and testing datasets.

|  | SemEval
---|---|---
Dataset | PxP1 | 2014 | 2015 | 2016
#reviews | 71,842 | 2,068 | 1,075 | 1,393
#unique words | 70,578 | 4,451 | 2,777 | 3,285
avg #sentences | 1.48 | 0.93 | 0.91 | 0.89
avg #words | 20.36 | 16.85 | 14.99 | 14.9

We removed numerical and non-English words as well as stop-words, emojis, and
punctuations from reviews. We used flair’s part-of-speech tagger Akbik et al.
(2019), and lemmatized extracted phrases according to their part-of-speech.

###  3.1 Evaluation Strategy

We estimate ${\theta}$ and ${\varphi}$ of our proposed model on Google reviews
of restaurants with no supervision about whether the reviews have surface form
for aspects. We trained the model for different numbers of aspects from 5 to
100 for restaurants and used the coherence score Mimno et al. (2011) to find
the optimum one. We set ${\bm{\alpha}}$ and ${\bm{\beta}}$ priors to 5.0
divided by the number of aspects and 0.01, respectively. As seen in Fig. 1,
the coherence reaches a plateau at 30 aspects for restaurants.

To evaluate our trained model, we used SemEval datasets whose reviews have
aspects with surface form (explicit) and are annotated. Given a review from
SemEval, we removed the aspect from the review (make it latent) and used our
model to infer the review’s latent aspect as the one with the highest
probability among all aspects of restaurants, i.e.,
${a=max(\theta\sim\text{Dir}(\bm{\alpha}))}$. The predicted latent aspect $a$
is a probability distribution of words that should be compared to the removed
aspect. A true prediction of latent aspect $a$ should have the surface form of
the removed aspect in the top-${k}$ most probable words of
${a\in\varphi_{a}}$.

For our training and test datasets are disjoint, we may face out of vocabulary
(OOV) at evaluation. For instance, the word ‘sushi’ has been labeled as an
aspect in SemEval reviews while our training dataset misses it. In this case,
a true prediction of latent aspect $a$ should have the most semantically
similar words to the removed explicit aspect (e.g., ‘sushi’) in the top-${k}$
most probable words of ${a\in\varphi_{a}}$ (e.g., ‘food’, ‘dish’). We used the
Resnik similarity score Resnik (1999) to calculate the inter-word semantic
similarities.

We report mean reciprocal rank (MRR), recall, nDCG, and success (hit) at 5 as
evaluation metrics for the performance of our trained model.

Figure 1: Coherence score for PxP restaurant dataset. Figure 2: Comparison
results for the baselines and our model across different SemEval test sets.

###  3.2 Baselines

We benchmarked our proposed model against the following baselines:

Random: is a naive method that chooses a random noun word as the latent
aspect.

locLDA Brody and Elhadad (2010) is an unsupervised method to identify explicit
aspects. This method uses the standard implementation of LDA, but they failed
to consider that a review may have no surface form for aspects.

CMLA Wang et al. (2017) is a supervised aspect-opinion co-extraction model
that leverages attention mechanisms.

HAST Li et al. (2018) is a supervised model on the tasks of explicit aspect
detection. It takes into account the history of aspects with an attention
block using bi-directional LSTMs.

OTE-MTL Zhang et al. (2020) is a supervised multitask learning framework to
extract aspect and opinion words jointly and parse sentiment dependencies
between them.

PxP is our proposed model. In contrast to supervised baselines, our model is
trained without any direct human supervision. Also, unlike locLDA, it takes
into consideration that a review might have a latent aspect with no surface
form.

###  3.3 Results

In this section, we seek to answer two research questions when dealing with
noisy, short, unsolicited reviews in which aspects are latent:

RQ1. Does unsupervised modeling of aspects as latent variables lead to more
accurate aspect detection compared to an unsupervised baseline that detect
explicit aspects?

RQ2. Do unsupervised methods fare better in the task of latent aspect
extraction compared to supervised methods?

Fig. 2 summarizes the comparative results for customer reviews for restaurants
in terms of evaluation metrics. X-axis shows to what extent aspects are
removed, as explained in the evaluation strategy. For example, at 0.3, we hide
labeled explicit aspects in 30% of SemEval’s reviews and use the baselines to
infer them. As seen, all of the baselines suffer from masking of explicit
aspects in the restaurant domain but at different rates. In response to RQ1,
it is evident that our unsupervised model consistently outperforms the state-
of-the-art unsupervised baseline (loc-LDA) in terms of all evaluation metrics
and across three different test sets. As shown, modeling aspects as latent
variables helps with better detection of latent aspects in the context of
noisy and short restaurant reviews.

Regarding RQ2, our proposed unsupervised method fares better in the task of
latent aspect extraction compared to the state-of-the-art supervised method
(CMLA) as the number of latent aspects increases. For instance, on SemEval-14
dataset, our model (PxP) could consistently maintain higher MRR score when
aspects are latent in more than 20% of reviews. Our model predicts a ranked
list of possible aspects for a given review, even for a review whose words
were totally unseen in the training phase. The supervised baselines however
fall short and their performance decreases linearly as the number of latent
aspects increases.

##  4 Conclusion

In this paper, we proposed an unsupervised LDA-based method to detect latent
aspects of products and services in unsolicited reviews which are noisy, short
and aspects are not explicitly mentioned due to social background knowledge.
We assume that writing a review is a generative process of deciding on an
aspect amongst the set of all available aspects for a product or service, and
then writing the most interrelated words for the chosen aspect. In contrast to
supervised methods, our model is trained without direct human supervision, and
unlike existing unsupervised methods, it takes into consideration that a
review might have a latent aspect. Benchmark results show that our model could
outperform supervised and unsupervised baselines when reviews happen to have
latent aspects.

## References

  * Akbik et al. (2019) Alan Akbik, Tanja Bergmann, Duncan Blythe, Kashif Rasul, Stefan Schweter, and Roland Vollgraf. 2019.  Flair: An easy-to-use framework for state-of-the-art nlp.  In _NAACL 2019_ , pages 54–59.
  * Blei et al. (2003) David M. Blei, Andrew Y. Ng, and Michael I. Jordan. 2003.  Latent dirichlet allocation.  _J. Mach. Learn. Res._ , 3:993–1022.
  * Brody and Elhadad (2010) Samuel Brody and Noemie Elhadad. 2010.  An unsupervised aspect-sentiment model for online reviews.  In _NAACL 2010_ , pages 804–812.
  * Cai et al. (2021) Hongjie Cai, Rui Xia, and Jianfei Yu. 2021.  Aspect-category-opinion-sentiment quadruple extraction with implicit aspects and opinions.  In _ACL 2010_.
  * Chen and Chen (2016) Huan-Yuan Chen and Hsin-Hsi Chen. 2016.  Implicit polarity and implicit aspect recognition in opinion mining.  In _ACL 2016_.
  * Cheng et al. (2014) Xueqi Cheng, Xiaohui Yan, Yanyan Lan, and Jiafeng Guo. 2014.  BTM: Topic modeling over short texts.  _IEEE Trans. Knowl. Data Eng. Transactions on Knowledge and Data Engineering_ , 26(12):2928–2941.
  * Cortis and Davis (2021) Keith Cortis and Brian Davis. 2021.  Over a decade of social opinion mining: a systematic review.  _Artif Intell Rev_ , 54(7):4873–4965.
  * Fei et al. (2021) Hao Fei, Yafeng Ren, Shengqiong Wu, Bobo Li, and Donghong Ji. 2021.  Latent target-opinion as prior for document-level sentiment classification: A variational approach from fine-grained perspective.  In _Proceedings of the Web Conference 2021_. ACM.
  * Hai et al. (2011) Zhen Hai, Kuiyu Chang, and Jung jae Kim. 2011.  Implicit feature identification via co-occurrence association rule mining.  In _Computational Linguistics and Intelligent Text Processing_ , pages 393–404. Springer Berlin Heidelberg.
  * He et al. (2017) Ruidan He, Wee Sun Lee, Hwee Tou Ng, and Daniel Dahlmeier. 2017.  An unsupervised neural attention model for aspect extraction.  In _ACL 2017_.
  * Li et al. (2018) Xin Li, Lidong Bing, Piji Li, Wai Lam, and Zhimou Yang. 2018.  Aspect term extraction with history attention and selective transformation.  In _IJCAI 2018_ , pages 4194–4200.
  * Li et al. (2016) Yonggan Li, Xueguang Zhou, Yan Sun, and Huanguo Zhang. 2016.  Design and implementation of weibo sentiment analysis based on LDA and dependency parsing.  _China Commun._ , 13(11):91–105.
  * Marrese-Taylor et al. (2013) Edison Marrese-Taylor, Juan D. Velasquez, and Felipe Bravo-Marquez. 2013.  Opinion zoom: A modular tool to explore tourism opinions on the web.  In _2013 IEEE/WIC/ACM International Joint Conferences on Web Intelligence (WI) and Intelligent Agent Technologies (IAT)_. IEEE.
  * Mimno et al. (2011) David M. Mimno, Hanna M. Wallach, Edmund M. Talley, Miriam Leenders, and Andrew McCallum. 2011.  Optimizing semantic coherence in topic models.  In _EMNLP 2011_ , pages 262–272.
  * Napitu et al. (2017) Fiernad Napitu, Moch. Arif Bijaksana, Agung Trisetyarso, and Yaya Heryadi. 2017\.  Twitter opinion mining predicts broadband internet's customer churn rate.  In _2017 IEEE International Conference on Cybernetics and Computational Intelligence_. IEEE.
  * Peng et al. (2020) Haiyun Peng, Lu Xu, Lidong Bing, Fei Huang, Wei Lu, and Luo Si. 2020.  Knowing what, how and why: A near complete solution for aspect-based sentiment analysis.  _AAAI Conference on Artificial Intelligence_ , 34(05):8600–8607.
  * Pontiki et al. (2016) Maria Pontiki, Dimitris Galanis, Haris Papageorgiou, Ion Androutsopoulos, Suresh Manandhar, Mohammad AL-Smadi, Mahmoud Al-Ayyoub, Yanyan Zhao, Bing Qin, Orphée De Clercq, Véronique Hoste, Marianna Apidianaki, Xavier Tannier, Natalia Loukachevitch, Evgeniy Kotelnikov, Nuria Bel, Salud María Jiménez-Zafra, and Gülşen Eryiğit. 2016.  SemEval-2016 task 5: Aspect based sentiment analysis.  In _Proceedings of the 10th International Workshop on Semantic Evaluation (SemEval-2016)_ , pages 19–30, San Diego, California. Association for Computational Linguistics.
  * Pontiki et al. (2015) Maria Pontiki, Dimitris Galanis, Haris Papageorgiou, Suresh Manandhar, and Ion Androutsopoulos. 2015.  SemEval-2015 task 12: Aspect based sentiment analysis.  In _Proceedings of the 9th International Workshop on Semantic Evaluation (SemEval 2015)_ , pages 486–495, Denver, Colorado. Association for Computational Linguistics.
  * Pontiki et al. (2014) Maria Pontiki, Dimitris Galanis, John Pavlopoulos, Harris Papageorgiou, Ion Androutsopoulos, and Suresh Manandhar. 2014.  SemEval-2014 task 4: Aspect based sentiment analysis.  In _Proceedings of the 8th International Workshop on Semantic Evaluation (SemEval 2014)_. Association for Computational Linguistics.
  * Poria et al. (2014) Soujanya Poria, Erik Cambria, Lun-Wei Ku, Chen Gui, and Alexander Gelbukh. 2014\.  A rule-based approach to aspect extraction from product reviews.  In _Proceedings of the Second Workshop on Natural Language Processing for Social Media (SocialNLP)_. Association for Computational Linguistics and Dublin City University.
  * Resnik (1999) P. Resnik. 1999.  Semantic similarity in a taxonomy: An information-based measure and its application to problems of ambiguity in natural language.  _jair_ , 11:95–130.
  * Sun et al. (2017) Shiliang Sun, Chen Luo, and Junyu Chen. 2017.  A review of natural language processing techniques for opinion mining systems.  _Information Fusion_ , 36:10–25.
  * Wan et al. (2020) Hai Wan, Yufei Yang, Jianfeng Du, Yanan Liu, Kunxun Qi, and Jeff Z. Pan. 2020.  Target-aspect-sentiment joint detection for aspect-based sentiment analysis.  _AAAI Conference on Artificial Intelligence_ , 34(05):9122–9129.
  * Wang et al. (2017) Wenya Wang, Sinno Jialin Pan, Daniel Dahlmeier, and Xiaokui Xiao. 2017.  Coupled multi-layer attentions for co-extraction of aspect and opinion terms.  In _Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, February 4-9, 2017, San Francisco, California, USA_ , pages 3316–3322. AAAI Press.
  * Xianghua et al. (2013) Fu Xianghua, Liu Guo, Guo Yanyan, and Wang Zhiqiang. 2013.  Multi-aspect sentiment analysis for chinese online social reviews based on topic modeling and HowNet lexicon.  _Knowledge-Based Systems_ , 37:186–195.
  * Xiao and Stibor (2010) Han Xiao and Thomas Stibor. 2010.  Efficient collapsed gibbs sampling for latent dirichlet allocation.  In _Proceedings of the 2nd Asian Conference on Machine Learning, ACML 2010, Tokyo, Japan, November 8-10, 2010_ , volume 13 of _JMLR Proceedings_ , pages 63–78. JMLR.org.
  * Xu et al. (2015) Hua Xu, Fan Zhang, and Wei Wang. 2015.  Implicit feature identification in chinese reviews using explicit topic mining model.  _Knowledge-Based Systems_ , 76:166–175.
  * Zeng and Li (2013) Lingwei Zeng and Fang Li. 2013.  A classification-based approach for implicit feature identification.  In _Lecture Notes in Computer Science_ , pages 190–202. Springer Berlin Heidelberg.
  * Zhang et al. (2020) Chen Zhang, Qiuchi Li, Dawei Song, and Benyou Wang. 2020.  A multi-task learning framework for opinion triplet extraction.  In _Findings of the Association for Computational Linguistics: EMNLP 2020_. Association for Computational Linguistics.
  * Zong et al. (2021) Chengqing Zong, Rui Xia, and Jiajun Zhang. 2021.  Sentiment analysis and opinion mining.  In _Text Data Mining_ , pages 163–199. Springer Singapore.