# Anatomically aware simulation of patient-specific
glioblastoma xenografts

Adam A. Malik  Mathematical Sciences, Chalmers University of Technology,
Gothenburg, Sweden.  Cecilia Krona  Department of Immunology, Genetics and
Pathology, Uppsala University, Uppsala, Sweden.  Soumi Kundu  Department of
Immunology, Genetics and Pathology, Uppsala University, Uppsala, Sweden.
Philip Gerlee  Mathematical Sciences, Chalmers University of Technology,
Gothenburg, Sweden.  Sven Nelander  Department of Immunology, Genetics and
Pathology, Uppsala University, Uppsala, Sweden.

###### Abstract

Patient-derived cells (PDC) mouse xenografts are increasingly important tools
in glioblastoma (GBM) research, essential to investigate case-specific growth
patterns and treatment responses. Despite the central role of xenograft models
in the field, few good simulation models are available to probe the dynamics
of tumor growth and to support therapy design. We therefore propose a new
framework for the patient-specific simulation of GBM in the mouse brain.
Unlike existing methods, our simulations leverage a high-resolution map of the
mouse brain anatomy to yield patient-specific results that are in good
agreement with experimental observations. To facilitate the fitting of our
model to histological data, we use Approximate Bayesian Computation. Because
our model uses few parameters, reflecting growth, invasion and niche
dependencies, it is well suited for case comparisons and for probing treatment
effects. We demonstrate how our model can be used to simulate different
treatment by perturbing the different model parameters. We expect in silico
replicates of mouse xenograft tumors can improve the assessment of therapeutic
outcomes and boost the statistical power of preclinical GBM studies.

##  1 Introduction

Glioblastoma (GBM), the most common primary malignant brain tumor in adults,
poses a significant clinical challenge. Despite various treatment modalities,
such as surgical intervention, radiation therapy, chemotherapy, or their
combinations, patients with GBM face a grim prognosis, with a median survival
of just over a year after diagnosis [2, 17]. The primary obstacle in the
surgical management of high-grade glioma lies in its infiltrative nature, as
tumor cells have the capacity to migrate extensively through healthy brain
tissue, infiltrating regions critical for patient survival [17].

In the late 1930s, Scherer’s pioneering studies of GBM pathology delineated
four patterns of invasion, encompassing subpial spread, accumulation near
neurons or blood vessels, and migration along white matter tracts [28]. Today,
we know that the dynamic interaction between GBM cells and different
anatomical niches is crucial in mediating, facilitating, and shaping spread of
the tumor. For instance, invasion along blood vessels presents notable
advantages for glioma cells [35]; The perivascular space, characterized by its
ample fluid content, offers minimal physical hindrance. Additionally, the
basement membrane showcases distinctive extracellular matrix (ECM) proteins
that facilitate the movement of glioma cells. Of note, invasion along
preexisting microvessels in the perivascular space is likely a VEGF-
independent mechanism of tumor vascularization , meaning that this process
likely cannot be suppressed by VEGF inhibitors [3, 34]. In addition to the
perivascular niche, the interaction with neurons and nerve fibers is also a
crucial mediator of tumor spread. Radiological data strongly support that GBMs
follow white matter fibers and data from animal models support distinct
invasion rates in gray matter. Jointly, the tumor cells’ interactions with
blood vessels and white matter cause a high degree of anisotropy of the tumor,
and contributes to the difficulty in determining a region for surgical
resection or radiotherapy [14, 39]. Amounting evidence in preclinical models
also suggests that distinct molecular mechanisms may be at play in mediating
invasion along different routes [24, 22, 13, 9, 10]. Jointly, these
observations motivate investigation into the quantitative details of GBM
invasion. Towards this goal, it is important to develop a quantitative
understanding of how diverse invasion outcomes arise in GBM.

This paper seeks to align two distinct approaches to brain tumor invasion:
mouse patient-derived cell (PDC) xenografts and patient-specifc, agent-based
simulation. Orthotopic xenografts from patient-derived primary cells replicate
many aspects of human GBM, including the formation of a core lesion with
increased mitotic index, pleiomorphic cells, and neoangiogenesis [38, 40, 26].
Also, in similarity with the original tumors, xenografts primarily spread
along the vasculature, through white matter structures, or diffusively through
the parenchyma [23, 40, 15]. The hypoxic niche and the perivascular niche are
also important in driving tumor growth and invasion [29, 27, 15]. Comparisons
of orthotopic xenografts derived from different source patients have shown
case-specific survival times and growth patterns [37, 40], providing a
possible way to elucidate mechanisms of GBM-brain interaction.

In recent years, mathematical and computational models have become
increasingly prevalent in cancer research, particularly in the context of
understanding tumor growth. One common modeling approach for describing the
anisotropic growth of gliomas involves using a continuum description of a
population of cells subject to diffusion and proliferation. For example, the
Proliferation-Invasion model [33], published in 2000, uses two different
diffusion constants in white and gray matter. More sophisticated continuum
models have since been developed, with some incorporating adhesion mechanisms
between glioma cells and ECM components [12, 1]. This approach to modeling has
also been extended to incorporate advanced imaging data such as Diffusion
Tensor Imaging to support the simulations and assess therapeutic effects [1].

Another group of mathematical models are those which consider individual cells
as the entities of interest. These models can be based on cellular automata
and its extensions or random walks. For example, a lattice-gas cellular
automaton was developed to model glioma cell invasion in the brain,
incorporating diffusion tensor data [16]. Another group of agent-based models
involves random walks, with some models using diffusion tensor data to guide
cell migration [20, 11]. A benefit of such agent-based models is that they are
based on explicit assumptions regarding cell behavior, rather than
phenomenological principles of macroscopic behavior. It is possible to derive
macroscopic descriptions from microscopic descriptions in certain cases [12,
32], but it is a non-trivial task and it has its own limitations.

While anisotropic spread induced by white matter invasion has been extensively
studied, the influence of perivascular invasion on glioma growth has not
received the same level of attention. In [3], an agent-based model was
developed and compared against experimental data obtained from mouse xenograft
experiments, where cells were assumed to be attracted to blood vessels.
However, no study has analyzed the influence of both white matter and the
brain vasculature simultaneously.

In this work we fill this gap by developing a model that explicitly
incorporates both white matter and blood vasculature. The model has four
parameters which we fit to data: migration rate, proliferation rate, white
matter preference and blood vessel preference. By fitting these model
parameters to data from 8 different xenografts, we have the goal to quantify
the different invasion phenotypes observed in experiments. In particular, we
will aim to answer the following questions. (i) Does adding explicit blood
vasculature into the model add anything in terms of predictive power? (ii) Can
our model capture the growth patterns seen in mice? We end by demonstrating
how the model can be used to simulate different types of treatment scenarios.

Figure 1: Figure illustrating our workflow. (A) PDCs have been grown in mice,
followed by histological sectioning and staining. Perivascular growth and
invasion along white matter tracts is present. (B) Our mathematical model
incorporates 3D maps of vasculature and white matter tracts. (C) The
mathematical model is an agent-based model where cell migration is influenced
by white matter and vasculature. (D) Data is segmented before being used for
model fitting (E) using Approximate Bayesian Computation.

##  2 Results

###  2.1 Anatomy-aware simulation of GBM xenografts

Histological sections from xenografted GBM typically present with a central
tumor lesion, together with evidence of spread in the perivascular and white
matter compartments (Figure 1A). Yet, relatively little is known about the
relationship between such phenotypes and the underlying dynamic processes,
such as the rate of cell proliferation and migration. To clarify this, we
propose a simulation strategy based on three main components:

First, the simulation is done on a 3D scaffold (map) that describes the normal
mouse brain anatomy (Figure 1B). This was obtained by magnetic resonance
Diffusion Tensor Imaging (DTI) data from normal adult C57BL/6J mouse brains
with a resolution of 43 $\mu$m [18]. Using this data, we calculated white
matter regions by thresholding the fractional anisotropy (FA) at 0.5.
Intuitively, FA is a DTI-derived score that captures to what degree the water
molecules in each small region of the brain are free to move, ranging from
free movement in all directions (FA=0) to being constrained to a single
direction (FA=1). We also incorporated brain vasculature data from a previous
study that combined tissue-clearing microscopy with image processing to
capture the vasculature down to the capillary level [36].

Second, we developed a mathematical model using an agent-based approach
(Figure 1C). In this model, cells migrate and proliferate within a three-
dimensional lattice, and changes occur at discrete time steps. The total
number of cells is denoted by $N_{t}$, and each lattice site has a carrying
capacity of $K$ cells. Cells have probabilities of proliferating ($p_{p}$) and
migrating ($p_{m}$) at each time step, with their migration influenced by
their microenvironment, particularly white matter and blood vessels. We
introduced parameters, $w_{wm}$ and $w_{bv}$, representing the strength of
attraction toward white matter and blood vessels, respectively. The
probability of migrating in different directions depends on these parameters.

The third aspect of our methodology involves data preprocessing (Figure 1D)
and parameter estimation using Approximate Bayesian Computation [7] (ABC)
(Figure 1E). The key idea of ABC is to run a simulation across multiple,
randomly picked parameter values. By subsequently comparing the simulations to
the experimental data (histological slides of mouse brains), a range of
credible parameters is obtained. We used summary statistics based on the
geometric properties of tumors to compare simulated tumor growth with real
tumor data. These summary statistics include the number of connected
components, area, perimeter, filled area, and eccentricity. We performed
regression adjustment with a logit transform to ensure the adjusted posterior
samples were within the same range as the respective prior distributions [8].
The parameters include $p_{m}$, $p_{p}$, $w_{wm}$, and $w_{bv}$, defined in
Methods.

Overall, our approach involves data preprocessing, mathematical modeling, and
parameter estimation to develop a comprehensive model for glioblastoma growth.
This model incorporates both anatomical influences and agent-based simulation
for a more accurate representation of tumor behavior.

(a) 3D map of white matter. Derived from previously published diffusion tensor
imaging dataset. abcabcabcabcabcabcabcabc

(b) 3D map of blood vasculature. Derived from previously published dataset
obtained by multi-dye staining and DISCO tissue clearing.

###  2.2 Improved estimation of parameters by geometric similarity measure

The task of comparing the location, size and shape of different tumors,
obtained from histological sections or from simulations, is a challenging
endeavour. In order to make a quantitative comparison between tumors, a metric
has to be chosen. A frequently used metric is the Jaccard Index, which
measures fraction of overlap between images. Different metrics have their own
strengths and weaknesses. For this reason we have developed our own metric
that is based on geometric similarity, rather than spatial overlap. The goal
is to investigate if histological sections together with our geometric
similarity metric is sufficient to estimate the four growth and invasion
parameters ($p_{m}$, $p_{p}$, $w_{wm}$, and $w_{bv}$) of our simulation model.

To investigate this idea, we first generated synthetic datasets by varying
each parameter within specified ranges and conducting simulations (see Methods
for details). By testing two (low vs high) values for each of the four
parameters, we obtained $2^{4}=16$ simulations from which we extracted virtual
histology slices, with clear differences in tumor size and invasion phenotype
(Figure (3)). Given such ”ground truth” data, we evaluated three versions of
the ABC algorithm. First, in the simplest version, we matched simulations to
experiments using the Jaccard index (the fraction of tumor-containing pixels
that overlapped between experiment and simulation). We further proposed a
geometric similarity measure with regression adjustment (as explained in
Section 3.3), and a “direct” method without regression adjustment.

For each of the 16 datasets, we calculated the relative error between the true
parameters $p_{m},p_{p},w_{wm},w_{bv}$ and the estimated parameters for the 50
accepted parameters using the formula given in (4). The overall error for a
specific method is then the average error over the 50 accepted parameters, as
defined in (5). The averaged error $E$ over the 16 datasets was 0.857 for the
Jaccard measure, 0.856 for the geometric measure without regression
adjustment, and 0.512 when using the geometric measure with regression
adjustment. The errors for the 16 cases are visualized in Figure (3). Based on
these results, we opted to use the geometric measure with regression
adjustment to fit the model to the histological sections from the PDC
xenografts.

Figure 3: Caption

###  2.3 Model simulations reproduces overall size and tumor morphology

We fitted the full model that contains both white matter and blood vasculature
to a total of eight PDC xenograft experiments (four cell lines, 2 replicates
each). The results are reported in Table 1. The numbers represent the point
estimate followed by the 10% and 90% percentiles. We also present the largest
value of $\delta$, defined in (3), that allowed for exactly 50 parameters to
be accepted. This value gives some insights into how well the simulated tumors
fits the real tumors for each case. A large value of $\delta$ means that in
order to accept 50 parameters, the accepted error threshold for the summary
statistics was large, whereas a smaller value of $\delta$ implies that the
overall difference in summary statistics was smaller. Hence $\delta$ can be
viewed as a measure of the inverse quality of the fit. We also present the
best fit for each tumor in Figure (4).

Tumor number | Cell line |  $P_{m}$ [$\mu$m/h] |  $P_{p}$ [1/day] | $w_{wm}$ | $w_{bv}$ | $\delta$
---|---|---|---|---|---|---
1 | U3013MG |  9.12 (1.96, 22.82) |  0.32 (0.28, 0.37) |  0.31 (0.15, 0.49) |  0.69 (0.5, 0.87) | 0.98
2 | U3013MG |  20 (5.6, 41.5) |  0.19 (0.14, 0.24) |  0.38 (0.26, 0.51) |  0.62 (0.34, 0.8) | 0.58
3 | U3230MG |  6.67 (5.2, 7.95) |  0.22 (0.22, 0.23) |  0.49 (0.21, 0.73) |  0.51 (0.25, 0.7) | 0.9
4 | U3230MG |  3.9 (1.52, 5.1) |  0.23 (0.18, 0.26) |  0.64 (0.19, 0.89) |  0.12 (0.04, 0.24) | 0.47
5 | U3033MG |  4.28 (1.8, 7.2) |  0.09 (0.067, 0.13) |  0.15 (0.02, 0.46) |  0.27 (0.04, 0.49) | 0.38
6 | U3033MG |  1.87 (0.9, 2.7) |  0.12 (0.89, 0.16) |  0.15 (0.02, 0.28) |  0.32 (0.14, 0.54) | 0.42
7 | U3062MG |  1.72 (0.22, 4.4) |  0.035 (0.028, 0.04) |  0.39 (0.08, 0.65) |  0.24 (0.1, 0.53) | 0.53
8 | U3062MG |  4.5 (1.7, 8.6) |  0.019 (0.014, 0.025) |  0.16 (0.03, 0.3) |  0.34 (0.14, 0.52) | 0.7
Table 1: Estimated parameter values of the eight experiments, along with the
10% and 90% percentiles in parentheses.

(a) Tumor 1 (U3013MG)

(b) Tumor 2 (U3013MG)

(c) Tumor 3 (U3230MG)

(d) Tumor 4 (U3230MG)

(e) Tumor 5 (U3033MG)

(f) Tumor 6 (U3033MG)

(g) Tumor 7 (U3062MG)

(h) Tumor 8 (U3062MG)

Figure 4: Best fit for each of the eight tumors

We can see that the overall size, and whether the tumor is diffuse or has a
distinct boundary is captured well by the model simulations. However, the
exact spatial overlap does not always agree, which is apparent in e.g. Tumors
1, 7 and 8. This is an inherent feature of the metric used, which does not
take spatial overlap into account. The best fit of tumor 2 shows a very small
collection of cells, in agreement with the real tumor. However, it did not
capture the few single cells that have migrated far away in the real tumor.
The fits shows good agreement in terms of overall size, and in most cases also
general shape. It is noticeable in tumors 3, 5 and 6 that the growth is
limited at the top of the tumor edge. It appears to be another anatomical
structure which is not white matter, which has influenced the growth. Tumors
3-7 exhibit growth along corpus callosum, which is replicated in the
simulations in Tumors 3,4 and 6.

By inspecting the fitted parameter values in Table 1, we can see that there is
some agreement between the replicates for each of the four cell lines. For
U3013MG (Tumors 1,2), there is good agreement in terms of $w_{wm}$ and
$w_{bv}$, but less agreement for $P_{m}$ and $P_{p}$. However, it is worth
pointing out that tumor 2 appears to be spread over a larger area, and hence
expected to have a larger migration rate $P_{m}$.

For U3230MG, (tumors 3 and 4) both replicates agree well for all parameters
except preference for blood vasculature, $w_{bv}$. The cause for this is not
obvious. However we observe that their direction of spread appears similar
along the corpus callosum, but their morpholoigies differ significantly apart
from that.

The fitted parameters of U3033MG agree well across all four parameters.

For U3062MG (tumors 7 and 8) there is some agreement. However, the spatial
overlap between the simulated tumors and the real tumors is poor for both
tumors. The real tumors also show somewhat distinct qualities, with tumor 7
appearing to be compact with a rather distinct boundary, whereas tumor 8 is
more diffuse.

###  2.4 Addition of blood vasculature improves model fit

Since the combination of white matter and blood vasculature is a novel feature
of our mathematical model, we tested if this new extension adds to the
predictive capacity of the model. To evaluate this, we considered two
alternative models, one full model (as above), and one white matter-only model
(without vasculature). We performed 5,000 simulations for each model (using
the same parameter ranges) and pooled all simulations into a set of 10,000.
For each of the eight tumors, we went on to identify the 50 best parameter
sets. To compare the merit of each model (full vs white-matter-only) we
calculated the fraction of the 50 best simulations that came from either
model. (The rationale for this comparison comes from the theory of ABC; more
probable models will be over-represented, while less probable models will be
under-represented or absent [8]). In six of the eight cases the full model had
stronger support, and in the two cases where the white-matter only model was
preferred, the fraction of accepted parameters was close to 50% from each of
the two models. (Table 2).

Tumor number | Cell line |  Fraction of accepted parameters Full model |  Fraction of accepted parameters White matter-only model
---|---|---|---
1 | U3013MG | 0.66 | 0.34
2 | U3013MG | 0.48 | 0.52
3 | U3230MG | 0.58 | 0.42
4 | U3230MG | 0.84 | 0.16
5 | U3033MG | 0.72 | 0.28
6 | U3033MG | 0.6 | 0.4
7 | U3062MG | 0.82 | 0.18
8 | U3062MG | 0.46 | 0.54
Table 2: Model support calculated as the fraction of accepted parameters from
each of the two models.

###  2.5 Simulating treatment response

In this section we conduct simulation experiments which are designed to
resemble the outcome when treating cells with hypothetical drugs that can
reduce migration or proliferation. The aim is to determine if there is an
inter-tumor heterogeneity in the response to different combinations of drugs.

To this end we perform the following type of experiment. We assume that there
is a treatment that knocks out the mechanism (proliferation or migration) for
a proportion of all cells. For example, we assume that each cell either stops
migrating, or stops proliferating. The distribution of cells into these
subpopulations varies between a 100% - 0% distribution (only anti migration
drug), then a 90% - 10% distribution, until we reach the final distribution of
0% - 100% (only anti proliferation drug). We implement the drug action by
assuming that, in each time step, each cell either does not migrate or does
not proliferate, and the probabilities are given by the assumed drug
distribution. For example, when we simulate a distribution corresponding to
70% anti migration and 30% anti proliferation treatment, the cell fails to
move with probability 0.7, and if it did not fail to move, it fails to
proliferate.

Our experimental design serves the purpose of incorporating a type of
effect/side-effect trade-off. Provided that there is some limit to the allowed
or tolerated amount of side effects, there may be cases where a tumor can be
more efficiently treated by reducing the dose of one treatment (and hence the
effect), and instead allowing the introduction of a secondary one. We use as a
starting point the fitted parameters from tumors 1-8 and initiate a tumor at
day 0 and let it grow for 30 days unhindered. At day 30 the treatment is
introduced and remains for another 60 days. To compare the treatments we
compare the final tumors to the case of an untreated control tumor growing for
90 days.

To measure the outcome of the different drug combinations on the growth of
tumors we use two different measures of tumor size. The first is the total
cell count, and the second is the convex hull of the tumor. We chose to use
the convex hull as a proxy for the spatial distribution of cells, which the
cell count does not account for. We illustrate our results by comparing the
treated tumor’s cell count and convex hull relative to the untreated tumor’s
cell count and convex hull in Figures 5 and Figure 6.

We can see that different tumors respond differently to the different drug
combinations. For tumors 1-4 and 6, it appears that the cell count is lowest
when the tumor is treated with only anti migration or only anti proliferation
drugs, with a maximum somewhere between the two extreme cases. The convex hull
on the other hand seem to be lowest when treated only with an anti migration
drug. This is expected to some extent, as the convex hull is capturing the
spatial spread of the tumor. On the other hand, in the cases of tumors 1-4 and
6, we see that there exists a point where the convex hull starts to decrease
again, with increasing doses of the anti proliferation drug. Interestingly,
tumors 5,7 and 8 does not show the same qualitative behavior. It appears that
the convex hull more or less increases monotonically with increasing
proportions of anti proliferation drug. In addition, the count has a minimum
when the treatment consists of only anti proliferation drug.

(a) Tumor 1 (U3013MG)

(b) Tumor 2 (U3013MG)

(c) Tumor 3 (U3230MG)

(d) Tumor 4 (U3230MG)

Figure 5: Cell count and convex hull of the treated tumor relative to the
untreated case for different combinations of anti migration and anti
proliferation drugs for tumors 1-4. The combination 100%/0% corresponds only
anti migration drug, and 0%/100% corresponds to only anti proliferation drug.

(a) Tumor 5 (U3033MG)

(b) Tumor 6 (U3033MG)

(c) Tumor 7 (U3062MG)

(d) Tumor 8 (U3062MG)

Figure 6: Cell count and convex hull of the treated tumor relative to the
untreated case for different combinations of anti migration and anti
proliferation drugs for tumors 5-8. The combination 100%/0% corresponds only
anti migration drug, and 0%/100% corresponds to only anti proliferation drug.

##  3 Material and Methods

###  3.1 data

Patient-derived cells (PDCs) obtained from the HGCC biobank [40, 19] were
labeled with a GFP-luciferase construct (pBMN(CMV-copGFP-Luc2-Puro, Addgene
plasmid #80389, a kind gift from Prof. Magnus Essand, Uppsala University), to
enable monitoring of tumor growth by in vivo bioluminescence imaging on a
NightOWL imaging system (Berthold Technologies). Following puromycin
selection, 100.000 cells were injected into 6-11 weeks old female
NOD/MrkBomTac-Prkdcscid (Taconic), NOD.CB17-Prkdcscid/scid/Rj (Janvier Labs),
NOD.Cg-Prkdcscid
Il2rgtm1Sug/JicTac (Taconic), and Rj:NMRI-Foxn1 nu/nu (Janvier Labs) under
2-4% isoflurane anesthesia. The cells were injected into a precise location of
the striatum using a stereotactic frame (KOPF model 940) and bregma as
reference point 0: AP 0, ML 1.5 (right), DV -3.0. Mice were given 5 mg/kg
carprofen, (Orion Pharma Animal Health, Sweden) subcutaneously before
injection and again the following day. All mice were monitored closely until
the endpoint defined by the occurrence of the first of the following three
events; (1) luciferase levels increase by a magnitude of 100-1000, (2) mice
display neurological symptoms or weight loss exceeding 10% of their maximum
weight, or (3) 40 weeks has passed since injection of cells. All mouse
experiments have been pre-approved by the regional animal research ethical
committee (permits C41/14, 5.8.18-02571-2017).

At the experimental endpoint, mice were euthanized using a gradually
increasing concentration of CO2 in air, followed by cervical dislocation.
Intact brains were harvested and saved in 4% buffered formaldehyde (Histolab
products AB, Sweden) after brief washing in PBS. For paraffin embedding, each
brain was first sliced coronally into five equal-thickness pieces (2 mm) and
processed for dehydration and paraffin embedding using an automated tissue
processor system (TPC15 DUO, Medite Medizintecknik, Germany). The tissue
blocks were prepared using a consistent pattern alignment with the five slices
to generate comparable 3-um thick sections after cutting and mounting on glass
slides.

To visualize tumor cells, sections were stained with anti-human specific NuMA
(ab 97585, Abcam) and Stem121 antibodies (Y40410, Takara) following standard
procedures and visualized with DAB Quanto (Thermo Fisher Scientific) following
incubation with horseradish peroxidase-conjugated secondary antibodies; goat
anti-mouse IgG (AP308P, Millipore) and goat anti-rabbit IgG (AP307P,
Millipore). To generate high-resolution images, sections were scanned in a
pyramidal structure with the highest magnification of 20x and resolution of
0.5 $\mu$m using an Aperio ScanScope XT Slide Scanner (Aperio Technologies,
Vista, CA, USA) at the Swedish Science for Life Laboratory (SciLifeLab) Tissue
Profiling facility at Uppsala University.

The clinical information of the patient derived cell lines used in this study
are summarized in Table 3.

Cell line |  Tumor number | Subtype |  Patient age |  Patient sex |  Survival (days) |  Experiment duration (days)
---|---|---|---|---|---|---
U3013MG | 1, 2 | Proneural | 78 | Female | 122 | 60, 66
U3230MG | 3, 4 | Mesenchymal | 61 | Male | 712 | 119, 95
U3033MG | 5, 6 | Classical | 58 | Male | 178 | 165, 136
U3062MG | 7, 8 | Mesenchymal | 76 | Female | 166 | 282, 275
Table 3: Clinical data for the four patient-derived cell lines used in this
study.

###  3.2 Mathematical model

We consider an agent-based model where cells are migrating and proliferating
in a three-dimensional lattice with spacing $\Delta x$, $\Delta y$ and $\Delta
z$, and where changes to the population occur at discrete time steps. We
denote by $N_{t}$ the total number of cells at time step $t=0,1,\ldots,T$ of
duration $\Delta t$. Each lattice site has a carrying capacity and can
accommodate at most $K$ cells at the same time. In each time step, $N_{t}$
agents are chosen independently at random, one at a time. When an agent is
chosen it is given the opportunity to proliferate, with probability $p_{p}$,
provided the number of cells is less than $K$ in the current voxel. The
offspring is placed in the proliferating cell’s voxel, but the offspring
cannot proliferate or move in the current time step. After a possible
proliferation event, the cell attempts to move, with probability $p_{m}$. If
the cell performs a motility event, it changes position to one of the six
neighboring lattice sites in its von Neumann neighborhood. In the case of
unbiased migration all six sites are equally likely. However, we assume that
the cells can sense their microenvironment consisting of white matter and
blood vessels. The details of how the environmental cues influences cell
migration is discussed in detail in the next section. The cell moves into the
target site only if the site is below its carrying capacity. If the lattice
site is full, the cell aborts its attempt to move. Once $N_{t}$ agents have
been given the opportunity to proliferate and migrate, the simulation time is
updated, $N_{t+1}$ computed, and the procedure is repeated for $T$ time steps.

As initial conditions we seed cells in a spherical region with radius 8 and
center at the location corresponding to the injection site. In our 3D lattice
this is the site located at $(94,133,115)$. In each voxel we place 3 cells,
resulting in a total of 6327 cells. The boundary conditions are such that if a
cell attempts to move out of the brain, the migration attempt is aborted and
the cell remains in place.

####  3.2.1 Influence of white matter and blood vasculature on migration

In order to model the preferential motion of cells towards white matter and
blood vessels, we introduce a bias to the random migration. We do that by
assuming that each cell can perceive its surrounding region, and becomes more
likely to move towards voxels containing white matter or blood vasculature.
Let us consider the case of how cells are attracted towards white matter
first. The aim is to compute how much the cell is influenced by white matter
which is within a distance $d$ from the cell. This is done by first finding
all voxels within a distance $d$ to the cell, which contains white matter. Let
a particular cell be located at $\mathbf{x}_{cell}$. Assume that there are a
total of $i$ voxels containing white matter, located at
$\mathbf{x}^{(i)}_{wm}$, which our cell can sense. We then compute the vector
pointing from $\mathbf{x}_{cell}$ towards $\mathbf{x}^{(i)}_{wm}$ and divide
by the Euclidean distance:

$\mathbf{v}^{(i)}=\frac{\mathbf{x}^{(i)}_{wm}-\mathbf{x}_{cell}}{||\mathbf{x}^{(i)}_{wm}-\mathbf{x}_{cell}||_{2}}=(x_{i},y_{i},z_{i}).$
(1)

The vector $\mathbf{v}^{(i)}$ points in the direction of the white matter
voxel, but has unit length. Recall that a migrating cell can move in either of
the 6 directions in its von Neumann neighborhood. We therefore have to convert
the direction vector $\mathbf{v}^{(i)}$ into a vector $\mathbf{u}^{(i)}$ of
length 6 that represents how $\mathbf{v}^{(i)}$ influences the probability of
migrating in each of the 6 directions. For instance, if
$\mathbf{v}^{(i)}=(1,0,0)$ we want the migration to be towards the positive
$x-$direction, whereas if $\mathbf{v}^{(i)}=(-1,0,0)$ we want the migration to
be in the negative $x-$direction. We therefore define $\mathbf{u}^{(i)}$ have
6 elements corresponding to the positive and negative $x-$axis, the positive
and negative $y-$axis, as well as the positive and negative $z-$axis. We then
calculate

$\mathbf{u}^{i}=\begin{pmatrix}\max(x_{i},0)\\\ \min(x_{i},0)\\\
\max(y_{i},0)\\\ \min(y_{i},0)\\\ \max(z_{i},0)\\\
\min(z_{i},0)\end{pmatrix}.$

Finally, we add up all the contributions from the $\mathbf{u}^{(i)}$’s and
normalize to 1, so that we obtain vector containing probabilities to move in
each of the 6 directions:

$\mathbf{P}_{wm}=\frac{\sum_{i}\mathbf{u}^{(i)}}{||\sum_{i}\mathbf{u}^{(i)}||_{1}}$

where $i$ is the index of number of neighboring voxels containing white
matter. The exact same procedure is performed for the blood vasculature
voxels, resulting in a vector of length 6, $\mathbf{P}_{bv}\in\mathbb{R}^{6}$.
For both white matter and vasculature we use a sensing radius of $d=5$ voxels.
The interpretation should be that $\mathbf{P}_{wm}$ contains the probabilities
to migrate in each of the 6 directions if the migration is entirely governed
by white matter. However, different cell lines show different preferential
movement towards different anatomical structures, so we introduce two key
parameters representing the strength of attraction towards white matter,
$w_{wm}$ and blood vessels, $w_{bv}$, respectively. Both parameters takes
values between 0 and 1, and must satisfy $w_{wm}+w_{bv}\leq 1$. The
probability to migrate in each of the 6 directions is given by

$\mathds{1}\frac{1}{6}(1-w_{wm}-w_{bv})+w_{wm}\mathbf{P}_{wm}+w_{bv}\mathbf{P}_{bv},$

where $\mathds{1}$ is a vector of ones, of length 6. The first term represents
the isotropic part, and the second and third terms the contributions from
white matter and vasculature respectively.

###  3.3 Parameter estimation using Approximate Bayesian Computation

In order to estimate the four model parameters $p_{m}$ (probability to
migrate), $p_{p}$ (probability to proliferate), $w_{wm}$ (white matter
preference) and $w_{bv}$ (blood vessel preference), we use a version of the
Approximate Bayesian Computation method [31] (ABC). ABC methods are often used
when the computation of a likelihood function is expensive or intractable. In
short, we specify prior distributions for our four parameters which are
assumed to be uniformly distributed, and then sample parameter values, perform
a simulation and compare how well such a simulation matches the data at hand.
We perform 5000 simulations, and choose to accept the top 1% best performing
simulations.

One of the major challenges with fitting model parameters to data is how to
choose a summary statistic used to compare a simulated tumor with a real
tumor. One of the most common metrics used for comparing spatial data is the
Jaccard index, defined for sets $A$ and $B$ through

$J(A,B)=\frac{|A\cap B|}{|A\cup B|},$

which for binary images of tumors measures the amount of overlap relative to
the total area from both tumors. However, such a measure performs poorly for
certain morphologies, such as diffusely growing tumors with low cell density.
Consider for example a checkerboard type of pattern in two images, one
containing the white squares and one the black ones. Although by visual
inspection two such patterns look the same, their Jaccard index will be 0.
This is also true for metrics such as the $L^{2}$-norm applied to individual
pixels. For this reason we have decided to use summary statistics based on
geometric properties of the tumors, which are not based on spatial location of
the tumors. For binary images we first find all connected components, using
MATLABs built-in function bwconncomp, and then we compute the Number of
connected components, Area, Perimeter, FilledArea and Eccentricity. These
measures are obtained for each connected component in the binary image. Once
calculated we define the summary statistics as

$\mathbf{s}=\begin{pmatrix}\mbox{Num. connected components}\\\
\mbox{Mean(Area)}\\\ \mbox{std(Area)}\\\ \mbox{mean(Eccentricity)}\\\
\mbox{std(Eccentricity)}\\\ \mbox{std(Perimeter)}\\\ \mbox{max(Perimeter)}\\\
\mbox{max(FilledArea)}\end{pmatrix},$ (2)

where “std” is the standard deviation. We denote the summary statistics of
simulation $i$ by $\mathbf{s}_{i}$, and the one from a experimental tumor
image by $\mathbf{s}$. In ABC it is common practice to define a threshold
$\delta$ so that parameters satisfying

$||\mathbf{s}_{i}-\mathbf{s}||<\delta$ (3)

are accepted, or by accepting a fixed proportion (1% or 0.1%) of simulations
[30]. In this work we do the latter, and accept 1% of the 5000 simulations.

After having obtained a sample of 50 parameters of each type we perform local
linear weighted regression [7] with an Epanechnikov kernel (the bandwidth is
chosen to be the smallest value of $\delta$ that permits exactly 50
simulations) and a logit transform to ensure that the adjusted posterior
samples lie within the same range as the respective prior [30]. Once an
adjusted sample has been obtained we perform kernel smoothing using an
Epanechnikov kernel with bandwidth $1/10$, and use the mean value of the
resulting estimated posterior density function as an estimate of the parameter
value.

The parameters are randomly sampled from the uniform prior distributions

$\displaystyle p_{m}\sim U(0,1),$ $\displaystyle p_{p}\sim U(0,0.015),$
$\displaystyle w_{wm}\sim U(0,1),$ $\displaystyle w_{bv}\sim U(0,1).$

We choose a carrying capacity of $K=3$ and perform 1800 time steps. At the end
time we save the final simulation configuration, i.e. a 3D lattice where each
voxel is populated with 0,1,2 or 3 cells. Although a slightly larger value
would be more realistic, it would induce significant computational cost.

The mouse experiments run for different durations due to differences in tumor
formation rate and growth rate, whereas each simulation consist of 1800 time
steps. Because of this, a simulation time step corresponds to different
durations in real time in the eight experiments. Once the model has been
fitted to one image from each tumor, we rescale the fitted parameters to
measure migration rate $P_{m}$ ($\mu m/\mbox{hour})$ and proliferation rate
$P_{p}$ ($1/\mbox{day}$), i.e. the same units (per hour and per day) in all of
the eight experiments. Prior to applying the method to experimental data we
apply it to synthetic data to demonstrate that our summary statistics perform
on par with the Jaccard index on synthetic datasets when no regression
adjustment is applied, and considerably outperforms it after applying
regression adjustment.

### Simulation study

In our simulation study to assess different similarity metrics, we ran 16 test
cases with each of the four parameters at a high and a low value. Values used
were
$p_{m}\in\\{0.2,1\\},p_{p}\in\\{0.0005,0.005\\},w_{wm}\in\\{0.1,0.45\\},w_{bv}\in\\{0.1,0.45\\}$.
We thus obtained datasets, and the combination of parameter values are shown
in Figure (3).

For each dataset we test three different methods to estimate the parameters.
Each method result in a set of 50 accepted parameters. The three methods we
compare are the Jaccard index and our geometric similarity measure with and
without regression adjustment (see Supplementary Information Section 3.3 for
details). When no regression adjustment is used, meaning that parameter
estimates are chosen to be the mean of the accepted parameters, we refer to it
as being a “direct” method.

For each of the 16 datasets we compute the relative error between true
parameters $p_{m},p_{p},w_{wm},w_{bv}$ and the estimated parameters
$\hat{p}^{(i)}_{m},\hat{p}^{(i)}_{p},\hat{w}^{(i)}_{wm},\hat{w}^{(i)}_{bv}$,
for the 50 accepted parameters, $i=1,2,\ldots,50$, using

$E_{i}=\frac{1}{4}\left(\left|\frac{p_{m}-\hat{p}^{(i)}_{m}}{p_{m}}\right|+\left|\frac{p_{p}-\hat{p}^{(i)}_{p}}{p_{p}}\right|+\left|\frac{w_{wm}-\hat{w}^{(i)}_{wm}}{w_{wm}}\right|+\left|\frac{w_{bv}-\hat{w}^{(i)}_{bv}}{w_{bv}}\right|\right),$
(4)

which is the absolute relative error, averaged over the four model parameters.
The error of a given method is then given as the average error over the 50
accepted parameters:

$E=\frac{1}{50}\sum_{i=1}^{50}E_{i}.$ (5)

We find that the error $E$ defined in (5) averaged over the 16 datasets is
0.857 using the Jaccard measure, 0.856 when using the geometric measure, and
0.512 when using the geometric measure followed by regression adjustment. The
errors for the 16 cases is shown in Figure (3).

Based on these results we decided to use the geometric measure with regression
adjustment of similarity when fitting the model to the histological sections
obtained from the PDC xenografts.

##  4 Discussion

We have proposed a novel model of glioblastoma growth in mice that explicitly
incorporates both white matter and blood vasculature, in order to be able to
capture a wider range of tumor morphologies. By fitting the model to eight
images obtained from mouse experiments we find that the model replicates the
overall tumor size and shape (diffusive vs. bulky) very well in most cases.
Moreover, it is able to reproduce the growth along the corpus callosum that is
a common feature of glioblastoma growth in mouse models and patients. Because
we use a geometric measure of similarity rather than a similarity measure
based on spatial overlap, like the Jaccard index, the exact position of the
simulated tumors are not always accurate (see e.g. tumor 3 in Figure 4), but
instead we obtain a larger morphological similarity between real and simulated
tumours.

The choice of similarity measure used to compare tumors is of utmost
importance, and remains one of the greatest challenges when comparing tumor
images in general. To some extent, the choice of metric will be influenced by
the application at hand. In some cases the exact position of the tumor or its
boundary would have higher priority, and in other applications the overall
size and structure (bulky vs. diffuse, alignment with anatomical structures or
not) will be of interest. Yet another consideration is whether the initial
tumor location is known. If the model simulation is initialized at the
incorrect site in the brain or if image misalignment is present, a measure
based solely on spatial overlap will perform poorly. Based on tests on
synthetic data (see Section (2.2)), we found that the geometric measure with
regression adjustment performs better than the Jaccard index on average for
our model, even in the case where the previously mentioned sources of error
are not present.

Our investigation into the response to different hypothetical drug
combinations demonstrated a considerable inter-tumor heterogeneity, where some
tumors exhibited a clear trade-off while others were monotone in their
response. Furthermore, these results highlight the utility of patient-specific
models to elucidate the most fruitful treatment options. However, for such
models to be clinically relevant they would need to be validated on much
larger datasets.

Throughout the process of modelling the growth of tumor cells injected into
mice, a number of sources of error and uncertainty are introduced. First, the
DTI dataset and the blood vasculature are obtained from different mice, and
hence needs to be aligned. Second, once a tumor has been simulated, a slice
has to be chosen for comparison to the corresponding histology section.
Currently this was done manually by visual inspection of the section and the
white matter structures present. Third, when removing the needle during the
tumor transplantation procedure, a small number of cells can follow the needle
and initiate tumor growth at other locations than the injection site. Fourth,
mice of different strains may not be anatomically identical.

When setting up a mathematical model, a number of simplifying assumptions have
to be made. In our model the rules of motion of cells is such an example. We
assumed that cells can sense and respond to presence of white matter and blood
vasculature within a spherical region in its vicinity, and that the attraction
is stronger the closer the cell is to the region. The exact details of how
cells in actuality respond to such cues is not fully known. However, in an
agent-based model like ours it is a straightforward task to test other
assumptions of motion. It is also possible to incorporate other mechanisms
such as cell-cell adhesion, chemotaxis or subpopulations of cells with
different behavior.

In this work we chose to develop an agent based model, which is but one of
many types of models that could be used. We made this choice due to their
flexibility and ease with which one can implement different rules that dictate
the cells’ behavior and interactions with their environment. Another group of
models are the continuum models in the form of partial differential equations
describing the spatial and temporal evolution of the density of cells. These
types of models can be phenomenologically motivated or derived from agent
based models. Agent based models are best suited when the entity of interest
is the cell, and where stochastic effects may play a role, whereas continuum
models are more suitable for large populations of cells.

A number of important topics warrants further research. First, our model
should be systematically benchmarked against other models to determine its
performance, in particular after increasing the number of simulations used in
the ABC method. Second, it should be determined whether the estimated
parameters correlates with known markers of white matter invasion and
perivascular invasion. To do that the model would have to be fitted to a
larger number of tumors which are also characterised on the genomic and
proteomic scale.

##  5 Supplementary information

###  5.1 Diffusion tensor imaging

Diffusion tensor imaging was first introduced in clinical neuroradiology in
1994 [6, 5] and offers a noninvasive method for studying the tissue anatomy
and the neural connectivity of the brain. The effective diffusion tensor is
estimated by applying a sequence of magnetic field gradients and measuring the
echo signal [4], which gives information about the diffusion of water
molecules in the brain. Since anatomical structures impede the motion of
water, the movement of water molecules can be used to infer the location and
direction of white matter tracts [25]. The resulting data from the DTI
measurements is a tensor field, where a rank 2 tensor is associated with each
voxel in the spatial domain. Each tensor is symmetric and has six independent
elements. To obtain a better intuitive understanding of the diffusion tensor
one can consider the following thought experiment as described in [4]. A
collection of water molecules diffuse away from the center of a voxel in an
isotropic medium. After some time of diffusion, the distribution of water
molecules is spherically symmetric and the surfaces of constant probability
are concentric spheres. In an anisotropic medium however, the rate of
diffusion is not equal in all directions, and the resulting surfaces of
constant probability are instead ellipsoids. Such three-dimensional diffusion
ellipsoids are fully described by six parameters, which determine its size,
shape and orientation. The diffusion ellipsoid can also be characterized by
three eigenvalues $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ and their
corresponding eigenvectors $\mathbf{v}_{1}$, $\mathbf{v}_{2}$ and
$\mathbf{v}_{3}$, where $\lambda_{1}\mathbf{v}_{1}$,
$\lambda_{2}\mathbf{v}_{2}$ and $\lambda_{3}\mathbf{v}_{3}$ correspond to the
perpendicular axes of symmetry of the ellipsoid. The eigenvector
$\mathbf{v}_{1}$ corresponding to the largest eigenvalue $\lambda_{1}$ is
called the principal eigenvector.

To describe the characteristic size of the diffusion ellipsoid in a given
voxel, one can compute the trace of the tensor $D$, which is independent of
its orientation and shape, and given by the sum of its eigenvalues:

$Tr(D)=D_{xx}+D_{yy}+D_{zz}=\lambda_{1}+\lambda_{2}+\lambda_{3},$ (6)

where $D_{xx}$, $D_{yy}$ and $D_{zz}$ are the diagonal elements of the tensor.
The average of the three eigenvalues is a quantity referred to as the
directionally averaged diffusivity [25],
$D_{avg}=(\lambda_{1}+\lambda_{2}+\lambda_{3})/3$. A commonly used measure of
the amount of anisotropy within a voxel is given by the fractional anisotropy
[25]:

$FA=\frac{\sqrt{(\lambda_{1}-\lambda_{2})^{2}+(\lambda_{2}-\lambda_{3})^{2}+(\lambda_{3}-\lambda_{1})^{2}}}{\sqrt{2}\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2}}}$
(7)

which is 0 in the case of isotropic diffusion
($\lambda_{1}=\lambda_{2}=\lambda_{3}$), and 1 in the case where one
eigenvalue is nonzero, and the other two zero.

###  5.2 Data processing

####  5.2.1 Processing the DTI data

The DTI data was obtained upon request by the authors of [18], the original
publication of the dataset. The DTI datasets comes from normal adult C57BL/6J
mouse brains, and acquired with an isotropic Nyquist limited resolution of 43
$\mu m$. We used the diffusion tensor labelled by N56014 for this study. In
order to calculate white matter regions we first calculated the fractional
anisotropy (FA) of each voxel, resulting in a $256\times 420\times 256$ 3D
array. We finally calculate the white matter regions by thresholding at $0.5$,
meaning that all voxels where the FA is above 0.5 are considered white matter,
and voxels where it is below 0.5 non-white matter. In addition to this we also
find the brain volume and brain shell by binarizing the FA at 0, meaning that
every voxel which has an FA greater than 0 is considered to belong to the
brain. This is used to align the brain volume obtained from the DTI dataset to
the brain volume obtained from the blood vasculature. The alignment process is
described in Section 5.2.3.

####  5.2.2 Processing the blood vasculature dataset

The dataset containing the blood vasculature map was published in [36] and is
publicly available for download. We used the dataset named _BALBc-no2 iso3um
stitched segmentation subsamples_. The dataset was reduced in size by a factor
2 (every other voxel removed) to be able to fit into 16 GB of ram. The
resulting array is of size $1555\times 2392\times 991$.

####  5.2.3 Alignment of DTI and blood vasculature datasets

The two binary 3D arrays containing the white matter and the blood vasculature
needs to be aligned into a single brain region, containing both white matter
and blood vasculature. The first step we take is to manually crop both
volumes, so that as little empty space surrounding the regions remains. After
that we manually rotate both volumes so that their orientation in 3D space is
the same. Next we resize the blood vasculature brain to have the same size as
the DTI brain, which is the smaller of the two. Then we apply a finite
iterative closest point algorithm [21] on the point cloud making up the brain
shells, which provides a transformation which maps one brain onto the other.
Once the two brains are aligned, some misalignment will remain. In order to
avoid white matter or blood vasculature to stretch outside the brain boundary,
we define the final brain volume to be the union of the two aligned brain
volumes. The final brain domain, which constitutes the domain of the
stochastic simulations, is a 3D array of size $250\times 346\times 182$, and
contains both white matter and blood vasculature.

####  5.2.4 Comparing simulations and experimental data

All simulations result in a 3D array of size $250\times 346\times 182$. In
order to compare a microscopy image of a tumor to a simulated tumor we do the
following steps.

  1. 1.

Decide which slice in the simulated brain that corresponds to the microscopy
slice.

  2. 2.

Binarization step

  3. 3.

Resize the microscopy image to have the same size as the simulated tumor. In
this case $250\times 182$, corresponding to a coronal section.

  4. 4.

Align the brain regions using the finite iterative closest point method [21].

  5. 5.

Apply the transformation obtained from the brain regions to align the tumors.

  6. 6.

Binarize both tumors using Matlabs built in function im2bw. A threshold of 0.1
was used on the real tumor, and 0.05 on the simulated tumor. The thresholds
were manually selected.

(a) Segmented microscopy image

(b) Binary microscopy image

(c) Original simulation image

(d) Binary simulation image

Figure 7: Images illustrating the original tumor/simulated tumor, and their
binary counterparts.

## References

  * [1] Abramo Agosti, Clara Cattaneo, Chiara Giverso, Davide Ambrosi, and Pasquale Ciarletta.  A computational framework for the personalized clinical treatment of glioblastoma multiforme.  ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 98(12):2307–2327, 2018.
  * [2] Komal Anjum, Bibi Ibtesam Shagufta, Syed Qamar Abbas, Seema Patel, Ishrat Khan, Sayed Asmat Ali Shah, Najeeb Akhter, and Syed Shams ul Hassan.  Current status and future therapeutic perspectives of glioblastoma multiforme (gbm) therapy: a review.  Biomedicine & Pharmacotherapy, 92:681–689, 2017.
  * [3] Gregory J Baker, Viveka Nand Yadav, Sebastien Motsch, Carl Koschmann, Anda-Alexandra Calinescu, Yohei Mineharu, Sandra Ines Camelo-Piragua, Daniel Orringer, Serguei Bannykh, Wesley S Nichols, et al.  Mechanisms of glioma formation: iterative perivascular glioma growth and invasion leads to tumor progression, vegf-independent vascularization, and resistance to antiangiogenic therapy.  Neoplasia, 16(7):543–561, 2014.
  * [4] Peter J Basser.  Diffusion and diffusion tensor mr imaging: fundamentals.  Magnetic resonance imaging of the brain and spine, pages 1752–1767, 2008.
  * [5] Peter J Basser, James Mattiello, and Denis LeBihan.  Estimation of the effective self-diffusion tensor from the nmr spin echo.  Journal of Magnetic Resonance, Series B, 103(3):247–254, 1994.
  * [6] Peter J Basser, James Mattiello, and Denis LeBihan.  Mr diffusion tensor spectroscopy and imaging.  Biophysical journal, 66(1):259–267, 1994.
  * [7] Mark A Beaumont, Wenyang Zhang, and David J Balding.  Approximate bayesian computation in population genetics.  Genetics, 162(4):2025–2035, 2002.
  * [8] Giorgio Bertorelle, Andrea Benazzo, and S Mona.  Abc as a flexible framework to estimate demography over space and time: some cons, many pros.  Molecular ecology, 19(13):2609–2625, 2010.
  * [9] L. J. Brooks, M. P. Clements, J. J. Burden, D. Kocher, L. Richards, S. C. Devesa, L. Zakka, M. Woodberry, M. Ellis, Z. Jaunmuktane, S. Brandner, G. Morrison, S. M. Pollard, P. B. Dirks, S. Marguerat, and S. Parrinello.  The white matter is a pro-differentiative niche for glioblastoma.  Nat Commun, 12(1):2184, 04 2021.
  * [10] M. S. Carro, W. K. Lim, M. J. Alvarez, R. J. Bollo, X. Zhao, E. Y. Snyder, E. P. Sulman, S. L. Anne, F. Doetsch, H. Colman, A. Lasorella, K. Aldape, A. Califano, and A. Iavarone.  The transcriptional network for mesenchymal transformation of brain tumours.  Nature, 463(7279):318–325, Jan 2010.
  * [11] L Leon Chen, Stephan Ulmer, and Thomas S Deisboeck.  An agent-based model identifies mri regions of probable tumor invasion in a patient with glioblastoma.  Physics in Medicine & Biology, 55(2):329, 2009.
  * [12] Christian Engwer, Thomas Hillen, Markus Knappitsch, and Christina Surulescu.  Glioma follow white matter tracts: a multiscale dti-based model.  Journal of mathematical biology, 71(3):551–582, 2015.
  * [13] E. Eskilsson, G. V. Røsland, G. Solecki, Q. Wang, P. N. Harter, G. Graziani, R. G. W. Verhaak, F. Winkler, R. Bjerkvig, and H. Miletic.  EGFR heterogeneity and implications for therapeutic intervention in glioblastoma.  Neuro Oncol, 20(6):743–752, 05 2018.
  * [14] Morteza Esmaeili, Anne Line Stensjøen, Erik Magnus Berntsen, Ole Solheim, and Ingerid Reinertsen.  The direction of tumour growth in glioblastoma patients.  Scientific reports, 8(1):1–6, 2018.
  * [15] Rajesh Kumar Gupta, Mia Niklasson, Tobias Bergström, Anna Segerman, Christer Betsholtz, and Bengt Westermark.  Tumor-specific migration routes of xenotransplanted human glioblastoma cells in mouse brain.  Scientific Reports, 14(1):864, 2024.
  * [16] Haralambos Hatzikirou and Andreas Deutsch.  Cellular automata as microscopic models of cell migration in heterogeneous environments.  Current topics in developmental biology, 81:401–434, 2008.
  * [17] Eric C Holland.  Glioblastoma multiforme: the terminator.  Proceedings of the National Academy of Sciences, 97(12):6242–6244, 2000.
  * [18] Yi Jiang and G Allan Johnson.  Microscopic diffusion tensor atlas of the mouse brain.  Neuroimage, 56(3):1235–1243, 2011.
  * [19] Patrik Johansson, Cecilia Krona, Soumi Kundu, Milena Doroszko, Sathishkumar Baskaran, Linnéa Schmidt, Claire Vinel, Elin Almstedt, Ramy Elgendy, Ludmila Elfineh, et al.  A patient-derived cell atlas informs precision targeting of glioblastoma.  Cell reports, 32(2):107897, 2020.
  * [20] Anitha Priya Krishnan, Delphine Davis, Paul Okunieff, and Walter O’Dell.  Random walk model based on dti for predicting the microscopic spread of gliomas.  In 2008 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pages 891–894. IEEE, 2008.
  * [21] Dirk-Jan Kroon.  Finite iterative closest point.
  * [22] B. Krusche, C. Ottone, M. P. Clements, E. R. Johnstone, K. Goetsch, H. Lieven, S. G. Mota, P. Singh, S. Khadayate, A. Ashraf, T. Davies, S. M. Pollard, V. De Paola, F. Roncaroli, J. Martinez-Torrecuadrada, P. Bertone, and S. Parrinello.  EphrinB2 drives perivascular invasion and proliferation of glioblastoma stem-like cells.  Elife, 5, 06 2016.
  * [23] Jeongwu Lee, Svetlana Kotliarova, Yuri Kotliarov, Aiguo Li, Qin Su, Nicholas M Donin, Sandra Pastorino, Benjamin W Purow, Neil Christopher, Wei Zhang, et al.  Tumor stem cells derived from glioblastomas cultured in bfgf and egf more closely mirror the phenotype and genotype of primary tumors than do serum-cultured cell lines.  Cancer cell, 9(5):391–403, 2006.
  * [24] H. Miao, N. W. Gale, H. Guo, J. Qian, A. Petty, J. Kaspar, A. J. Murphy, D. M. Valenzuela, G. Yancopoulos, D. Hambardzumyan, J. D. Lathia, J. N. Rich, J. Lee, and B. Wang.  EphA2 promotes infiltrative invasion of glioma stem cells in vivo through cross-talk with Akt and regulates stem cell properties.  Oncogene, 34(5):558–567, Jan 2015.
  * [25] P Mukherjee, JI Berman, SW Chung, CP Hess, and RG Henry.  Diffusion tensor mr imaging and fiber tractography: theoretic underpinnings.  American journal of neuroradiology, 29(4):632–641, 2008.
  * [26] Michele Patrizii, Monica Bartucci, Sharon R Pine, and Hatem E Sabaawy.  Utility of glioblastoma patient-derived orthotopic xenografts in drug discovery and personalized therapy.  Frontiers in oncology, 8:23, 2018.
  * [27] Gong Peng, Yin Wang, Pengfei Ge, Christopher Bailey, Peng Zhang, Di Zhang, Zhaoli Meng, Chong Qi, Qian Chen, Jingtao Chen, et al.  The hif1$\alpha$-pdgfd-pdgfr$\alpha$ axis controls glioblastoma growth at normoxia/mild-hypoxia and confers sensitivity to targeted therapy by echinomycin.  Journal of Experimental & Clinical Cancer Research, 40:1–16, 2021\.
  * [28] HJ Scherer.  Structural development in gliomas.  The American Journal of Cancer, 34(3):333–351, 1938.
  * [29] Sheng Si et al.  Notch signaling regulates metabolic heterogeneity in glioblastoma stem cells.  2017\.
  * [30] Scott A Sisson, Yanan Fan, and Mark Beaumont.  Handbook of approximate Bayesian computation.  CRC Press, 2018.
  * [31] Mikael Sunnåker, Alberto Giovanni Busetto, Elina Numminen, Jukka Corander, Matthieu Foll, and Christophe Dessimoz.  Approximate bayesian computation.  PLoS computational biology, 9(1):e1002803, 2013.
  * [32] Amanda Swan, Thomas Hillen, John C Bowman, and Albert D Murtha.  A patient-specific anisotropic diffusion model for brain tumour spread.  Bulletin of mathematical biology, 80(5):1259–1291, 2018.
  * [33] Kristin R Swanson, Ellsworth C Alvord Jr, and JD Murray.  A quantitative model for differential motility of gliomas in grey and white matter.  Cell proliferation, 33(5):317–329, 2000.
  * [34] Aaron C Tan, David M Ashley, Giselle Y López, Michael Malinzak, Henry S Friedman, and Mustafa Khasraw.  Management of glioblastoma: State of the art and future directions.  CA: a cancer journal for clinicians, 70(4):299–312, 2020.
  * [35] Emily G Thompson and Harald Sontheimer.  A role for ion channels in perivascular glioma invasion.  European Biophysics Journal, 45(7):635–648, 2016.
  * [36] Mihail Ivilinov Todorov, Johannes Christian Paetzold, Oliver Schoppe, Giles Tetteh, Suprosanna Shit, Velizar Efremov, Katalin Todorov-Völgyi, Marco Düring, Martin Dichgans, Marie Piraud, et al.  Machine learning analysis of whole mouse brain vasculature.  Nature methods, 17(4):442–449, 2020.
  * [37] Einar Osland Vik-Mo, Cecilie Sandberg, Havard Olstorn, Mercy Varghese, Petter Brandal, Jon Ramm-Pettersen, Wayne Murrell, and Iver Arne Langmoen.  Brain tumor stem cells maintain overall phenotype and tumorigenicity after in vitro culturing in serum-free conditions.  Neuro-oncology, 12(12):1220–1230, 2010.
  * [38] Hiroaki Wakimoto, Santosh Kesari, Christopher J Farrell, William T Curry Jr, Cecile Zaupa, Manish Aghi, Toshihiko Kuroda, Anat Stemmer-Rachamimov, Khalid Shah, Ta-Chiang Liu, et al.  Human glioblastoma–derived cancer stem cells: Establishment of invasive glioma models and treatment with oncolytic herpes simplex virus vectors.  Cancer research, 69(8):3472–3481, 2009.
  * [39] Jun Wang, Liang Yi, Qing-mei Kang, Ji Zhou, Tian-qing Chen, Jean-philippe Hugnot, and Shi-cang Yu.  Glioma invasion along white matter tracts: A dilemma for neurosurgeons.  Cancer letters, 526:103–111, 2022.
  * [40] Yuan Xie, Tobias Bergström, Yiwen Jiang, Patrik Johansson, Voichita Dana Marinescu, Nanna Lindberg, Anna Segerman, Grzegorz Wicher, Mia Niklasson, Sathishkumar Baskaran, et al.  The human glioblastoma cell culture resource: validated cell models representing all molecular subtypes.  EBioMedicine, 2(10):1351–1363, 2015.


SPL

    Software Product Line

/tikz/mandatory/.style= circle,fill=drawColor, draw=drawColor, inner
sep=0.25em, /tikz/optional/.style= circle, fill=white, draw=drawColor, inner
sep=0.25em, featureDiagram/.style= for tree= parent anchor = south, child
anchor = north, draw = drawColor, edge = draw=drawColor, ,
/tikz/abstract/.style= fill = blue!85!cyan!5, draw = drawColor ,
/tikz/concrete/.style= fill = blue!85!cyan!20, draw = drawColor ,
mandatory/.style= edge label=node [mandatory] , optional/.style= edge
label=node [optional] , or/.style= tikz+= .parent)coordinateA) – (!u.children)
coordinate (B) – (!ul.parent) coordinate (C) pic[fill=drawColor, angle
radius=0.8em]angle; , /tikz/or/.style= , alternative/.style= tikz+=
.parent)coordinateA) – (!u.children) coordinate (B) – (!ul.parent) coordinate
(C) pic[draw=drawColor, angle radius=0.8em]angle; , /tikz/alternative/.style=
, /tikz/placeholder/.style= , collapsed/.style= rounded corners, no edge, for
tree= fill opacity=0, draw opacity=0, l = 0em, , /tikz/hiddenNodes/.style=
midway, rounded corners, draw=drawColor, fill=white, minimum size = 1.2em,
minimum width = 0.8em, scale=0.9 ,

11affiliationtext: University of Ulm, Ulm, Germany22affiliationtext:
University of Magdeburg, Magdeburg, Germany33affiliationtext: TU Darmstadt,
Darmstadt, Germany44affiliationtext: University of Siegen, Siegen,
Germany55affiliationtext: Harz University of Applied Sciences, Wernigerode,
Germany

# T-Wise Presence Condition Coverage and Sampling for Configurable Systems

Sebastian Krieter  Thomas Thüm  Sandro Schulze  Sebastian Ruland  Malte Lochau
Gunter Saake  Thomas Leich

###### Abstract

Sampling techniques, such as $t$-wise interaction sampling are used to enable
efficient testing for configurable systems. This is achieved by generating a
small yet representative sample of configurations for a system, which
circumvents testing the entire solution space. However, by design, most recent
approaches for $t$-wise interaction sampling only consider combinations of
configuration options from a configurable system’s variability model and do
not take into account their mapping onto the solution space, thus potentially
leaving critical implementation artifacts untested. Tartler et al. address
this problem by considering presence conditions of implementation artifacts
rather than pure configuration options, but do not consider the possible
interactions between these artifacts. In this paper, we introduce $t$-wise
presence condition coverage, which extends the approach of Tartler et al. by
using presence conditions extracted from the code as basis to cover $t$-wise
interactions. This ensures that all $t$-wise interactions of implementation
artifacts are included in the sample and that the chance of detecting
combinations of faulty configuration options is increased. We evaluate our
approach in terms of testing efficiency and testing effectiveness by comparing
the approach to existing $t$-wise interaction sampling techniques. We show
that $t$-wise presence condition sampling is able to produce mostly smaller
samples compared to $t$-wise interaction sampling, while guaranteeing a
$t$-wise presence condition coverage of 100%.

##  1 Introduction

Testing is an important task in software engineering to detect faults and to
check intended behavior [ammann2016introduction, M10]. However, exhaustive
testing may be impossible and binds resources that could be used in other
phases of development. This is especially an issue when testing highly-
configurable systems, such as Software Product Lines. A product of an SPL can
be derived from a configuration, which consist of a list of configuration
options (i.e.,features) that can be set either to _true_ or _false_. The
entirety of all possible configurations for an SPL is called the
_configuration space_ , which typically grows exponentially with the number of
features [ER:IST11, LKL:SPLC12].

A straight-forward testing strategy for SPLs is product-based testing. For
this purpose, a set of products is derived from a set of different
configurations and then test cases are run on each selected product
respectively [TAK+:CSUR14]. As testing every possible configuration in a
product-based manner is usually not feasible due to an enormous configuration
space, sampling strategies have been defined to generate a small yet
representative set of configurations to test (i.e., a configuration sample).
One such sampling strategy is _$t$ -wise interaction sampling_, which aims to
generate a preferably small sample that covers all possible interactions of
features of degree $t$ [CDS:TSE08, MGSH:SPLC13]. Using $t$-wise interaction
sampling, developers can ensure that all interactions of at most $t$ features
(e.g., all selected, none selected, only one selected, etc.) are indeed
contained in at least one configuration in the generated sample. $T$-wise
interaction sampling has shown to be a feasible trade-off between testing
effectiveness and testing efficiency, as for small values of t (i.e.,
$t\in\\{2,3\\}$) it usually returns a relatively small sample, while also
achieving reliable test results [AMS+:TOSEM18, KAR+:SPLC09, MGSH:SPLC13].

A property of $t$-wise interaction sampling is that it works purely on the
problem space of an SPL. Thus, it is a black-box approach that does not take
into account the mapping between features and actual implementation artifacts,
such as source code, models, and test cases. This can lead to a number of
issues. First, a traditional t-wise sampling algorithm may create samples that
contain configurations that generate similar products, which can decrease
overall testing efficiency. Second, when choosing the parameter for $t$ the
degree of faulty interactions is not known to the developers. Thus, developers
are incentivized to choose a low value for $t$ in order to keep the testing
effort feasible. However, some faults may require a feature interaction of a
high degree to be included in a configuration, and thus may not be included at
all in a sample, which potentially decreases testing effectiveness.

Due to the black-box nature of $t$-wise interaction sampling, it may not
reveal certain faults resulting from feature interaction beyond $t$ or
requires too many configurations to even reach a certain code coverage or
fault-detection rate. To overcome this limitation, Tartler et al. [TDS+:ATC14]
propose _statement coverage_ , a white-box approach that considers _presence
conditions_ (i.e., the selection of features for which an artifact is included
in a product) to ensure that every artifact is present in at least one
configuration in the sample and gets a chance of being tested. However, simply
including an implementation artifact in a product does not grantee that is
indeed being tested. Therefore, Ruland et al. [RLB+:GPCE18] argue that, in
addition to including each artifact, at least one test case must also cover
each artifact in order to properly test it. While the approach of Tartler et
al. [TDS+:ATC14] can generate a relatively small sample, the approach of
Ruland et al. [RLB+:GPCE18] produces a relatively large sample, but guarantees
that each artifact will indeed be tested. Therefore, with this work, we aim to
find a reasonable trade-off between both approaches by building on the
approach of Tartler et al. [TDS+:ATC14] and try to increase testing
effectiveness by combining it with t-wise interaction sampling.

We propose an extension of Tartler et al.’s coverage criterion $t$-wise
presence condition coverage, which combines $t$-wise interaction coverage with
presence condition coverage of implementation artifacts. Rather than counting
only interactions between features, our new criterion considers interactions
of presence conditions of implementation artifacts. A sample with a $t$-wise
presence condition coverage of 100% ensures that every _$t$ -wise interactions
of all implementation artifacts_ is contained in at least one product and can
thus be tested. In addition, we present an algorithm for $t$-wise presence
condition sampling, which generates samples with a 100% $t$-wise presence
condition coverage for a given configurable system and a given $t$. This
sampling algorithm works independently from the employed variability
mechanisms such as preprocessors, build systems, feature-oriented programming,
or aspect-oriented programming. We implemented our algorithm within a
prototype called _PRESICE (PRESence condItion CoveragE)_ to investigate its
testing effectiveness, testing efficiency, and sampling efficiency compared to
traditional $t$-wise interaction sampling algorithms. Our current prototype
supports extraction of presence conditions from SPLs that use the C
preprocessor and the Kbuild build tool (e.g., BusyBox and Linux).
Additionally, we extensively evaluate both, the coverage criterion and the
algorithm using 27 real-world systems from different domains. In summary, we
contribute the following:

  * •

We propose a novel coverage criterion for product-based testing: $t$-wise
presence condition coverage.

  * •

We present a sampling algorithm for $t$-wise presence condition coverage.

  * •

We provide an open-source implementation named named PRESICE as part of
FeatureIDE111https://featureide.github.io/.

  * •

We evaluate our coverage criterion and algorithm using 27 real-world systems.

  * •

We publish all data from our
experiments222https://github.com/skrieter/evaluation-pc-sampling.

##  2 Running Example and Problem Statement

In the following, introduce a running example, which we use throughout the
paper. We use this example to describe the potential problems with the current
approaches for $t$-wise interaction sampling and motivate our solution.

###  2.1 Running Example

To illustrate the challenge of finding effective and efficient samples, we
show a slightly edited code snippet from the system BusyBox in Figure 1, which
uses the C preprocessor [stallman1987c] to implement its variability
[LAL+:ICSE10]. The example is taken from the file tftp.c, which handles
client-server communication via tftp. The code contains five features, TFTP
$(T)$, TFTP_GET $(G)$, TFTP_PUT $(P)$, TFTP_BLOCKSIZE $(B)$, and TFTP_DEBUG
$(D)$, which each can be set to either true or false. We display their
dependencies in LABEL:fig_busbox_fm using an excerpt of the BusyBox
variability model, represented as a feature diagram. The additional features
BusyBox_TFTP $(BB)$ and TFTPD $(TD)$ do not appear in the code snippet, but
are necessary to visualize the feature diagram hierarchy. In the remainder of
the paper, we use the provided abbreviations of the feature names to ease the
readability of all propositional formulas using these features.

In the example, we changed the statement in Line 671 such that a compilation
error occurs for certain products. The variable blksize is declared in Line
626, which is dependent on feature $B$. Then, blksize is used in Line 671,
which is dependent on feature $D$. Thus, if $D$ is selected in a
configuration, but $B$ is not, the generated product will be syntactically
incorrect.

⬇

1#if TFTP_GET || TFTP_PUT // $G\vee P$

// ...

621# if TFTP // $T$

622int tftp_main(int argc, char **argv) {

// ...

625 # if TFTP_BLOCKSIZE // $B$

626 const char *blksize_str=TFTP_BLKSIZE_DEFAULT;

// ...

649 int blksize = tftp_blksize_check(blksize_str, 65564);

650 if (blksize < 0) return EXIT_FAILURE;

651 # endif

// ...

670 # if TFTP_DEBUG // $D$

671 printf("blksize = %d\n", blksize); // changed

672 # endif

// ...

690}

691# endif

// ...

827#endif

Figure 1: Example adopted from BusyBox.

  *[SPL]: Software Product Line
  *[SPLs]: Software Product Line