A stochastic individual-based model to explore the role of spatial interactions and antigen recognition in the immune response against solid tumours

TL;DR

This paper introduces a spatial stochastic individual-based model that captures antigen recognition and spatial interactions between immune and cancer cells, providing insights into immune response dynamics against solid tumours.

Contribution

It presents a novel model integrating antigen expression, recognition, and spatial movement of immune cells, advancing understanding of tumour-immune interactions.

Findings

Conditions for tumour clearance, dormancy, or escape identified

Antigen heterogeneity impacts immune efficacy

Spatial interactions influence immune response outcomes

Abstract

Spatial interactions between cancer and immune cells, as well as the recognition of tumour antigens by cells of the immune system, play a key role in the immune response against solid tumours. The existing mathematical models generally focus only on one of these key aspects. We present here a spatial stochastic individual-based model that explicitly captures antigen expression and recognition. In our model, each cancer cell is characterised by an antigen profile which can change over time due to either epimutations or mutations. The immune response against the cancer cells is initiated by the dendritic cells that recognise the tumour antigens and present them to the cytotoxic T cells. Consequently, T cells become activated against the tumour cells expressing such antigens. Moreover, the differences in movement between inactive and active immune cells are explicitly taken into account by the model. Computational simulations of our model clarify the conditions for the emergence of tumour clearance, dormancy or escape, and allow us to assess the impact of antigenic heterogeneity of cancer cells on the efficacy of immune action. Ultimately, our results highlight the complex interplay between spatial interactions and adaptive mechanisms that underpins the immune response against solid tumours, and suggest how this may be exploited to further develop cancer immunotherapies.

Tables

Table 1. Table 1: Model parameters and related values used in computational simulations. Note, standard deviation has been abbreviated to StD.

Symbol	Description	Value(s)	Reference
${\Delta }_t$	time-step	$1$ min	(Boissonnas et al., 2007)
${\Delta }_{x,y}$	grid spacing in the $x$ or $y$ direction	$10\mu$m	(Macfarlane et al., 2018)
$N_T\left(0\right)$	initial number of tumour cells	$400$ cells	(Christophe et al., 2015)
$N_C$	total number of CTLs	$400$ cells	(Christophe et al., 2015)
$N_D$	total number of DCs	$400$ cells	(Macfarlane et al., 2018)
$n$	Index identifier of each tumour cell	$n=1,\dots ,N_T$	-
$k$	Index identifier of each DC	$k=1,\dots ,N_D$	-
$m$	Index identifier of each CTL	$m=1,\dots ,N_C$	-
${𝐀}_{Tn}\left(t\right)$	Antigen profile of tumour cell $n$ at time $t$	values $\ge 0$	-
${𝐀}_{Dk}\left(t\right)$	Recognised antigen profile of DC $k$ at time $t$	values $\ge 0$	-
${𝐀}_{Cm}\left(t\right)$	Antigen receptor profile of CTL $m$ at time $t$	values of 0 or 1	-
$\lambda$	tumour cell division rate	$0.001{\text{min}}^{-1}$	(Christophe et al., 2015)
${\theta }_E^{\ast }$	Average probability of epimutations	$0.23$	(De Smet et al., 1996)
$\mu$	removal rate of tumour cells by CTLs	$0.03\text{cells}{\text{min}}^{-1}$	(Christophe et al., 2015)
$D_{Act}$	DC activation rate	$0.07\text{cells}{\text{min}}^{-1}$	(Bianca et al., 2012)
$C_{Act}$	CTL activation rate	$\approx 0.12\text{cells}{\text{min}}^{-1}$	(Engelhardt et al., 2012)
$\alpha$	Lévy walk exponent	1.15	(Harris et al., 2012)
$\beta$	T cell binding affinity	$[0,0.2]$	(Schmid et al., 2010)
$V_T$	StD of the initial tumour antigen profiles from the reference experimental profile $M$	$[0,1]$	-
$V_D$	StD of the antigen profiles recognised by DCs from the tumour antigen profiles	$[0,1]$	-
$V_{E/M}$	StD of the tumour antigen profiles from the previous profiles after epimutations/mutations	same as $V_T$	-

Table 2. Table 2: Average levels of expression of the MAGE-A antigens in oral squamous cell cancer cell lines. The experimental data were taken from  (Hartmann et al., 2016 ) and then normalised.

Antigen $\left(i\right)$	Mean Expression, $M_i$
MAGE-A1 (1)	0.10
MAGE-A2 (2)	0.25
MAGE-A3 (3)	0.41
MAGE-A4 (4)	0.24
MAGE-A5 (5)	0.36
MAGE-A6 (6)	0.35
MAGE-A8 (7)	0.17
MAGE-A9 (8)	0.16
MAGE-A10 (9)	0.38
MAGE-A11 (10)	0.06
MAGE-A12 (11)	0.32

Equations

$${\mathbf{A}}_{Tn}(t)=(A_{Tn}^{(1)}(t),\dots ,A_{Tn}^{(11)}(t)),$$

$$A_{Tn}^{(i)}(0)=(M_i+V_T{R}_i{)}_{+},i=1,\dots ,11.$$

$$A_{Tn}^{(i)}(t+{\Delta }_t)=(A_{Tn}^{(i)}(t)+V_E{R}_i{)}_{+},i=1,\dots ,11.$$

$$A_{Tn}^{(i)}(t+{\Delta }_t)=A_{Tn}^{(i)}(t),i=1,\dots ,11$$

$$A_{T{N}_T(t)+1}^{(i)}(t+{\Delta }_t)=(A_{Tn}^{(i)}(t)+V_M{R}_i{)}_{+},i=1,\dots ,11.$$

$${\mathbf{A}}_{Dk}(t)=(A_{Dk}^{(1)}(t),\dots ,A_{Dk}^{(11)}(t)).$$

$$A_{Dk}^{(i)}(0)=0,i=1,\dots ,11$$

$$A_{Dk}^{(i)}(t+{\Delta }_t)=(A_{Tn}^{(i)}(t)+V_D{R}_i{)}_{+},i=1,\dots ,11.$$

$${\mathbf{A}}_{Cm}(t)=(A_{Cm}^{(1)}(t),\dots ,A_{Cm}^{(11)}(t)).$$

$$A_{Cm}^{(i)}(0)=0,i=1,\dots ,11$$

$$A_{Cm}^{(i)}(t)=\{\begin{array}{l}1\text{mbox}fori=\hat{i},\\ 0\text{mbox}fori\ne \hat{i},\end{array}\text{with }\hat{i}=\underset{j}{\mathrm{argmax}}{A}_{Dk}^{(j)}(t),$$

$$A_{Tn}^{(i)}(t)\ge (M_i-\beta )\text{mbox}fori\text{ such that }{A}_{Cm}^{(i)}(t)=1.$$

$${\theta }_E^{\ast }=\mathrm{0.23.}$$

Full text

††journal: Journal of Theoretical Biology

1 Introduction

The immune system is a collection of cells, structures and processes that work together to remove harmful foreign material from the body. Understanding the mechanisms which underpin the cellular immune response to solid tumours is crucial for the development of new immunotherapy techniques, which can strengthen the natural human immune response to support the successful treatment of cancer (Couzin-Frankel, 2013; June et al., 2018; Mellman et al., 2011; Ribas and Wolchok, 2018).

Immune cells, specifically dendritic cells (DCs) and cytotoxic T lymphocytes (CTLs), detect and target solid tumours through recognising and processing the small peptides (i.e. the antigens) that are expressed by tumour cells (Messerschmidt et al., 2016). The antigenic composition of solid tumours can be heterogeneous, whereby each cell within the tumour mass may have an antigen profile characterised by different expression levels of tumour antigens. On the surface of a tumour cell, antigens are integrated with major histocompatibility complex (MHC) molecules allowing the DCs to recognise the antigens and present them to the CTLs. The CTLs then become activated and subsequently express the corresponding antigen receptor, which enables them to interact with the tumour cells expressing that particular antigen at a sufficiently high level (Coulie et al., 2014). Further interactions between CTLs and tumour cells can then trigger tumour cell death.

There are three distinct classes of tumour antigens: tumour associated antigens (TAAs), tumour specific antigens (TSAs) and cancer testis antigens (CTAs). These antigens are expressed, respectively, by: both normal and cancer cells, cancer cells only, and both cancer cells and human germ-line cells. Germ-line cells do not contain MHC molecules and, therefore, these cells cannot present antigens that can be recognised by T cells. This makes CTAs a viable target for immunotherapy as targeting these antigens implies a lower risk of autoimmune reactions (Coulie et al., 2014). One specific group of CTAs is represented by the melanoma associated antigen (MAGE) genes (Boon et al., 2006; Connerotte et al., 2008; Müller-Richter et al., 2009; Urosevic et al., 2005). Within this group, the MAGE-A family consists of eleven genes that are linked to poor prognosis in various types of cancer (Coulie et al., 2014; Hartmann et al., 2016; Zajac et al., 2017). In general, the expression of MAGE-A genes promotes tumour progression by enhancing tumour cell division and reducing apoptosis (Van Tongelen et al., 2017). Some of the specific functions of MAGE-A genes have been discovered and include: inducing chemoresistance through a decrease in p53 (Marcar et al., 2010; Monte et al., 2006; Yang et al., 2007), enhancing tumour proliferation (Costa et al., 2007) and decreasing the efficacy of anti-tumour drugs (Hartmann et al., 2013, 2014). The MAGE-A antigens are commonly expressed in melanomas (Boon et al., 2006; Connerotte et al., 2008; Coulie et al., 2014; Urosevic et al., 2005), oesophageal cancers (Zajac et al., 2017), lung, breast, prostate and colorectal carcinomas (Coulie et al., 2014), and head and neck cancers (Hartmann et al., 2016; Müller-Richter et al., 2009). Therefore, finding an effective way of targeting these antigens via immunotherapy may be beneficial to the treatment of multiple types of cancer (Zajac et al., 2017). For instance, clinical trials have shown that targeting MAGE-A3 can be a successful treatment option (Chinnasamy et al., 2011; Connerotte et al., 2008).

The antigen expression profiles of tumour cells can evolve over time due to epigenetic and genetic mechanisms. Amongst epigenetic mechanisms, spontaneous epimutations are ‘stochastic and heritable changes in gene expression that leave the sequence of bases in the DNA unaltered’ (Oey and Whitelaw, 2014). Through spontaneous epimutations, the MAGE-A antigen expression levels of tumour cells can change over time, which may result in variability between the antigen profiles of histological samples from the same patient or between patients (Urosevic et al., 2005).

One potential cause of spontaneous epimutations is DNA methylation, whereby methylation of specific promoter regions of the gene represses transcription. Consequently, if the gene is methylated the corresponding protein will not be expressed. MAGE genes are methylated in normal cells, but they can be de-methylated, and thus expressed, in cancer cells  (Boon et al., 2006; Chalitchagorn et al., 2004; Chinnasamy et al., 2011; Müller-Richter et al., 2009; Zajac et al., 2017). Demethylation of these genes can become more prominent during cancer progression, which suggests that demethylation of the MAGE genes may support tumour development (Coulie et al., 2014).

Cell-cell interactions involved in the immune response to cancer depend upon the spatial position of both immune cells and tumour cells within the tumour micro-environment (Chaplin, 2010; Messerschmidt et al., 2016). In particular, as a result of tumour growth, cells within the tumour can become more exposed to immune action depending on their location in relation to other cells of the tumour micro-environment (Hanahan and Weinberg, 2011). Moreover, stochastic antigenic variations can induce further spatial heterogeneity within the tumour (Boon et al., 2006; Yarchoan et al., 2017). Additionally, the movement of immune cells is dictated by the spatial distribution of tumour antigens within the tumour micro-environment (Boissonnas et al., 2007).

Mathematical models are a useful tool for simulating and investigating biological systems, and have been increasingly used to describe tumour antigen expression and tumour-immune interactions. Tumour antigen expression and the effects of epigenetic and genetic events have been modelled through differential equation models (Asatryan and Komarova, 2016; Lorenzi et al., 2016; Johnston et al., 2007; Tomasetti and Levy, 2010) and cellular-automaton (CA) models (Bouchnita et al., 2017; Manem et al., 2014). Traditionally, tumour antigen expression and recognition by the immune system have been implicitly modelled by tuning the rates of T cell recruitment, T cell proliferation or tumour cell removal (Arciero et al., 2004; Balea et al., 2014; Besse et al., 2018; De Boer et al., 1985; de Pillis et al., 2009; Köse et al., 2017; Mallet and de Pillis, 2006). More recently, these processes have been explicitly captured by mathematical models formulated in terms of either ordinary differential equations (Balachandran et al., 2017; d’Onofrio and Ciancio, 2011; Łuksza et al., 2017) or integro-differential equations (Delitala et al., 2013; Delitala and Lorenzi, 2013; Kolev et al., 2013; Lorenzi et al., 2015). However, these models rely on the assumption that cells are well-mixed and, as such, they do not reflect spatial aspects of tumour-immune competition. On the contrary, the spatial and temporal dynamics of tumour-immune interactions have been described through partial differential equation (PDE) models (Al-Tameemi et al., 2012; Matzavinos et al., 2004; Matzavinos and Chaplain, 2004) and hybrid PDE-CA models (de Pillis et al., 2006; Mallet and de Pillis, 2006), but these do not take into account the antigen expression and recognition processes.

In light of these considerations, we present here a spatial individual-based model of tumour-immune competition that explicitly captures antigen expression and recognition. In our model, each cancer cell is characterised by an antigen profile which can change over time due to either epimutations or mutations. The immune response against the cancer cells is initiated by the DCs that recognise the tumour antigens and present them to the CTLs. Consequently, T cells become activated against the tumour cells which express such antigens. Moreover, exploiting the modelling strategies that we have previously developed (Macfarlane et al., 2018), the differences in movement between inactive and active immune cells are explicitly taken into account. Computational simulations of this model clarify the conditions for the emergence of tumour clearance, dormancy or escape, and allow us to assess the impact of antigenic heterogeneity of cancer cells on the efficacy of immune action. Ultimately, our results highlight the complex interplay between spatial interactions and adaptive mechanisms that underpins the immune response against solid tumours, and suggest how this may be exploited to further develop cancer immunotherapies.

The remainder of this work is organised as follows. In Section 2, we present the individual-based model and detail how each biological mechanism is mathematically described. In Section 3, we parametrise the model and present the results of computational simulations. In Section 4, we discuss the results obtained and highlight their biological implications along with potential further applications of this work.

2 The mathematical model

Building upon our previous work (Macfarlane et al., 2018), we consider three cell types: tumour cells, dendritic cells and cytotoxic T lymphocytes. We use an on-lattice individual-based approach to describe the interactions between these three cell types. Our model is posed on a 2D spatial grid of spacing $\Delta_{x}$ in the $x$ direction and $\Delta_{y}$ in the $y$ direction, with the constraint that only one cell of any type is allowed at each grid-site, at any time-step of duration $\Delta_{t}$.

The system is initially composed of cancer cells and inactive immune cells only. The immune cells are randomly distributed on the spatial grid, while the cancer cells are tightly packed in a circular configuration positioned at the centre of the grid, to reproduce the geometry of a solid tumour. We let the tumour grow through cell division. A cancer cell divides at rate $\lambda$ into two progeny cells of which one occupies the position of the parent cell while the other is positioned at an unoccupied neighbouring grid-site. This ensures that only cancer cells with free grid-sites in their neighbourhood can divide. Inactive immune cells update their position according to a Lévy-like walk (Harris et al., 2012). This process allows them to move in a randomly chosen direction for a number of time-steps $s$ sampled from a Lévy distribution $L(s)$ with exponent $0\leq\alpha<2$, i.e. $L(s)\sim s^{-(\alpha+1)}$. DCs are activated at rate $D_{Act}$ upon contact with tumour cells, and CTLs become activated at rate $C_{Act}$ upon contact with active DCs. Upon activation, DCs and CTLs switch to Brownian motion, i.e. at each time-step they can move to any of the neighbouring grid sites with the same probability. Moreover, we let active CTLs remove tumour cells, upon contact, at rate $\mu$. Once activated, both DCs and CTLs remain active throughout the simulations. For simplicity, we omit natural death of tumour cells and proliferation of both DCs and CTLs (i.e. the total numbers of DCs and CTLs are constant over time). We refer the reader to our previous paper (Macfarlane et al., 2018) for a detailed description of these modelling strategies.

In this work, we develop each of these modelling strategies further to include the antigen profiles of cancer cells and their possible variation, the immune recognition of tumour antigens by DCs, and the targeted activation of CTLs against specific tumour antigens. The modelling strategies used to take into account such additional layers of biological complexity are described in detail in the following subsections, and are also schematically illustrated in Figure 1 and Figure 2.

2.1 Mathematical modelling of antigen expression

We denote by $N_{T}(t)$ the number of tumour cells in the system at time $t=h\,\Delta_{t}$, with $h\in\mathbb{N}_{0}$, and we label each cell by an index $n=1,\ldots,N_{T}(t)$. We incorporate antigen expression into our model by letting each tumour cell express eleven different antigens, to represent the eleven MAGE-A antigens that, as mentioned in Section 1, have a key role in tumour development (Coulie et al., 2014). These antigens are reported in Table 2 and we label them by an index $i=1,\ldots,11$. There can be high variability in each antigen’s expression between patients with the same type of cancer (Hartmann et al., 2016; Müller-Richter et al., 2009; Urosevic et al., 2005) and even within cancer cell samples from the same patient (Hartmann et al., 2016; Müller-Richter et al., 2009; Urosevic et al., 2005). Therefore, at each time instant $t$, we characterise the antigen profile of the $n^{th}$ tumour cell by means of a vector

[TABLE]

with $A_{Tn}^{(i)}(t)$ representing the expression level of antigen $i$. Biologically, there is evidence that there can be correlation between antigen expression in some cancers, e.g. epithelial ovarian cancer (Daudi et al., 2014), but not all cancers, e.g. hepatocelllular carcinoma (Roch et al., 2010). To consider a more generalised situation, we assume that the expression levels of each antigen $i$ can evolve independently from the others. As schematically illustrated in Figure 1, for each tumour cell $n$ we define the initial expression of the $i^{th}$ antigen as

[TABLE]

In equation (2.1), the parameter $M_{i}$ denotes a mean expression level of antigen $i$ taken from published experimental data, the values of which are reported in Table 2. The value of $R_{i}$ is sampled from a standard normal distribution centred at zero. In equation (2.1) we take the positive part of the right-hand side to ensure non-negativity of the antigen expression level. As $R_{i}$ is taken from a standard normal distribution, the parameter $V_{T}$ represents the standard deviation of the initial antigen profile from the experimental value $M_{i}$ (Kenney and Keeping, 1962). Therefore, the parameter $V_{T}$ determines how close the value of $A_{Tn}^{(i)}(0)$ will be to the value of $M_{i}$.

2.2 Modelling variations in antigen expression

At each time-step, we let the tumour cells divide at rate $\lambda$ [refer to the schemes in Figure 1b] and change their antigen profile either through epimutations or through mutations. We assume that epimutations can occur at any time during the life of a cell [refer to the schemes in Figure 1c and Figure 2a], whereas mutations take place during cell division and may cause the antigen profile of one progeny cell to be different from that of the parent cell [refer to the scheme in Figure 2b]. We allow epimutations and mutations to occur with probabilities $\theta_{E}$ and $\theta_{M}$, respectively. In the absence of changes in antigen expression (i.e. if $\theta_{E}=0$ and $\theta_{M}=0$), the antigen profiles of the tumour cells will remain constant over time, that is, ${\bf{A}}_{Tn}(t)={\bf{A}}_{Tn}(0)$ for all $t>0$. If antigenic changes do occur through epimutations or mutations, the antigen profiles of tumour cells are updated using the methods described in the following paragraphs.

Epimutations

A variation in the level of expression of the $i^{th}$ antigen of the $n^{th}$ tumour cell at the time instant $t$ due to an epimutation is modelled according to the following equation

[TABLE]

In equation (2.2), the value of $R_{i}$ is sampled from a normal distribution centred at zero and we take the positive part of the right-hand side to ensure non-negativity of the antigen expression level. Following on from the definition of the parameter $V_{T}$, since $R_{i}$ is taken from a standard normal distribution the parameter $V_{E}$ is the standard deviation of the updated antigen profile from the previous expression profile (Kenney and Keeping, 1962). Therefore, $V_{E}$ determines how close the value of $A_{Tn}^{(i)}(t+\Delta_{t})$ will be to the value of $A_{Tn}^{(i)}(t)$.

Mutations

Upon division at the time instant $t$, the $n^{th}$ tumour cell is replaced by two cells, one labelled by the index $n$ and the other one labelled by the index $N_{T}(t)+1$. If mutations do not occur, the progeny cells inherit the antigen profile of the parent cell, i.e. $A_{Tn}^{(i)}(t+\Delta_{t})=A_{Tn}^{(i)}(t)$ and $A_{T\,N_{T}(t)+1}^{(i)}(t+\Delta_{t})=A_{Tn}^{(i)}(t)$. Conversely, if a mutation occurs, the antigen profile of the progeny cells will be given by the following equations

[TABLE]

and

[TABLE]

Equations (2.3) and (2.4) rely on notation analogous to that of equation (2.2) and, therefore, $V_{M}$ is the standard deviation of the progeny antigen profiles from the parent antigen profiles (Kenney and Keeping, 1962). The value of $V_{M}$ determines how close the value of $A_{T\ N_{T}(t)+1}^{(i)}(t+\Delta_{t})$ will be to the value of $A_{Tn}^{(i)}(t)$.

2.3 Activation of immune cells

Activation of DCs

We denote by $N_{D}$ the number of DCs, which we assume to be constant, and we label each DC by an index $k=1,\ldots,N_{D}$. Activation of DCs occurs, at rate $D_{Act}$, through contact with tumour cells. At any time instant $t$, the $k^{th}$ DC is characterised by a recognised antigen profile

[TABLE]

We let all DCs be initially inactive and thus assume

[TABLE]

for all $k=1,\ldots,N_{D}$. There is biological evidence supporting heterogeneity within the antigen presentation process where a less prevalent antigen may be recognised and presented by the DCs (Boes et al., 2002; Fehres et al., 2014; Ljunggren et al., 1990). Additionally, it is known that the MAGE-A genes have similar homology (Roch et al., 2010; Zajac et al., 2017) and, therefore, there is a potential that they could be mis-recognised as each other (Graff-Dubois et al., 2002; Linette et al., 2013; Raman et al., 2016; Schueler-Furman et al., 1998; Tong et al., 2004). As schematically described by Figure 1d, we consider the case where there may be potential variation in the antigen recognition process. To capture this idea, upon activation through contact with the $n^{th}$ tumour cell at the time instant $t$, we assign the recognised antigen profile of the $k^{th}$ DC using the following equation

[TABLE]

In equation (2.5), the value of $R_{i}$ is sampled from a normal distribution centred at zero and we take the positive part of the right-hand side to ensure non-negativity of the antigen expression level. We can describe $V_{D}$ as the standard deviation of the antigen profile recognised by each DC from the actual tumour antigen profile (Kenney and Keeping, 1962). Therefore, $V_{D}$ determines how close the value of $A_{Dk}^{(i)}(t+\Delta_{t})$ will be to the value of $A_{Tn}^{(i)}(t)$. Following the method of our previous work (Macfarlane et al., 2018), once activated, we let the DCs remain activated against their recognised tumour antigen profile.

Activation of CTLs

We denote by $N_{C}$ the number of CTLs, which we assume to be constant, and we label each CTL by an index $m=1,\ldots,N_{C}$. As schematically described by Figure 1e, activation of CTLs occurs, at rate $C_{Act}$, through contact with activated DCs. At any time instant $t$, each CTL $m$ has a receptor profile

[TABLE]

We let all CTLs be initially inactive and thus assume

[TABLE]

for all $m=1,\ldots,N_{C}$. While DCs can recognise multiple types of antigens, CTLs can produce copies of one antigen receptor only (Brenner et al., 2008; Coico and Sunshine, 2015). This means that each CTL can only be activated against one of the eleven MAGE-A antigens. To capture this fact, upon activation through contact with the $k^{th}$ DC at the time instant $t$, we let the $m^{th}$ CTL become activated against the highest expressed antigen within the tumour antigen profile recognised by the $k^{th}$ DC, i.e. we assign the receptor profile of the $m^{th}$ CTL using the following equation

[TABLE]

where the index $\hat{i}$ specifies the target antigen of the activated CTL. Note that, if $|\operatorname*{arg\,max}_{j}A_{Dk}^{(j)}(t)|>1$, then we arbitrarily choose $\hat{i}=\min\operatorname*{arg\,max}_{j}A_{Dk}^{(j)}(t)$. Using the same assumptions as in our previous work, once activated, a CTL remains activated against the same tumour antigen.

2.4 Removal of tumour cells by activated CTLs

Upon contact, each activated CTL can induce death of the tumour cells which express a sufficiently high level of the CTL’s target antigen (Coulie et al., 2014; Stone et al., 2009), which we assume to be given by the mean antigen expression levels reported in Table 2. In particular, as schematically described by Figure 1f, when the $m^{th}$ CTL interacts with the $n^{th}$ tumour cell at time $t$, we compare the receptor profile ${\bf{A}}_{Cm}(t)$ with the antigen profile ${\bf{A}}_{Tn}(t)$ and let the tumour cell be removed from the system at rate $\mu$ provided that

[TABLE]

In equation (2.7), the parameter $\beta$ describes the binding affinity of the CTLs, which determines the range of tumour cells that each CTL can interact with. If $\beta$ is larger, then the CTL can recognise tumour cells with a lower level of expression of the antigen that they target. Independently of the outcome of the interaction, the CTL can subsequently interact with further tumour cells following the same process.

3 Computational simulations

3.1 Model parametrisation and simulation set up

We use a 2D spatial domain with $100$ grid sites, of length $\Delta_{x}=\Delta_{y}=10\ \mu$m, both in the $x$ and in the $y$ direction, which correspond to a domain of size $1\ \text{mm}^{2}$. Simulations were developed and run in Matlab, for an appropriate number of time-steps, with one time-step chosen to be $\Delta_{t}=1$ min, to allow for the resulting dynamics of the system to be investigated. All quantities in the results we report on in this section were obtained through averaging the results of $5$ simulations. We refer the interested reader to Macfarlane et al. (2018) for a detailed description of the parameterisation of the original model, and we describe here the way in which the additional components of the model were calibrated using the parameter values reported in Table 1.

Initial tumour antigen expression levels

Hartmann et al. (2016) investigated the levels of expression of the eleven MAGE-A antigens in oral squamous cell cancers of $38$ patients. The mean antigen expression levels were taken to be values between [math]-$12$ arbitrary units as an immune-reactivity score, which are normalised and reported in Table 2.

In our model, the initial expression level of each antigen for each tumour cell is defined by using equation (2.1) along with the values of $M_{i}$ from Table 2. To determine the value of the product $V_{T}\ R_{i}$ in equation (2.1) we consider the properties of $R_{i}$, which is a random value taken from a standard normal distribution. A standard normal distribution $\mathcal{N}(0,1)$, with mean 0 and standard deviation 1, has a 95$\%$ confidence interval of $\pm$1.96 (Kenney and Keeping, 1962). Therefore, for 95$\%$ of values we expect -1.96 $\leq R_{i}\leq$1.96 with the majority of the values being close to the mean value. We then use the parameter $V_{T}$ to control the minimum and maximum value of the product $V_{T}\ R_{i}$, i.e. when $V_{T}$=1 then $V_{T}\ R_{i}\in[-1.96,1.96]$, for most values. However, if $V_{T}$ is lower, e.g. $V_{T}$=0.1, then the values of the product $V_{T}\ R_{i}$ will also be lower, e.g. $V_{T}\ R_{i}\in[-0.196,0.196]$. To consider a wide range of biological situations corresponding to different initial levels of heterogeneity in tumour antigen profiles, we use a range of values, between 0 and 1, for the parameter $V_{T}$. As the experimental data in Table 2 are dimensionless, we consider also $V_{T}$ to be dimensionless.

Probabilities of epimutations and mutations

De Smet et al. (1996) found that cancer cell lines expressing the MAGE-A1 antigen were $23\%$ more likely to undergo demethylation events than tumour cell lines that did not express this antigen. Such a value is supported by other studies that consider the likelihood of DNA demethylation in various cancers (Chalitchagorn et al., 2004; Ehrlich, 2002). We make the assumption that this holds true for the other ten MAGE-A antigens, and let the average probability of epimutations in our model be

[TABLE]

In the model we consider the effect of increase or decreasing the probability of epimutations or mutations by setting $\theta_{E}$ and $\theta_{M}$, respectively, as multiples of $\theta_{E}^{*}$. The parameters $V_{E}$ and $V_{M}$ control how much the antigen expression of a tumour cell can change through epimutations or how much the antigen expression of a progeny tumour cell can change through mutations [refer to equations (2.2) and (2.3)]. The values of $V_{M}$ and $V_{E}$ are chosen to match the values of $V_{T}$. Since the experimental data in Table 2 are dimensionless, we also consider $V_{E}$ and $V_{M}$ to be dimensionless..

Antigen recognition process:

To consider a wide range of biological situations corresponding to different scenarios in terms of the number of CTLs activated against each antigen, we use a range of values, between 0 and 1, for the parameter $V_{D}$ [refer to equation (2.5)]. As the experimental data in Table 2 are dimensionless, we also consider $V_{D}$ to be dimensionless.

T cell binding affinity:

The binding affinity of a T cell is related to the association rate, that is the inverse of the dissociation rate $K_{D}$. In general this value is between $0.005\ \mu M^{-1}$ and $1\ \mu M^{-1}$ for all natural T cells (Davis et al., 1998; Slansky and Jordan, 2010). Furthermore, the MAGE-A T cell receptors association rates have been found to be even larger than this range, e.g. for MAGE-A3 the values are between $0.018\ \mu M^{-1}$ and $5.917\ \mu M^{-1}$ (Tan et al., 2015). However, Schmid et al. (2010) have shown that an association rate of $0.2\ \mu M^{-1}$ or higher did not improve the binding affinity and, therefore, higher binding affinities may have a limited effect. We take $M_{i}$ to be the minimal binding and allow the likelihood of binding to increase depending on $\beta$. In line with experimental evidence, we investigate a range of dimensionless values between [math] and $0.2$ for the parameter $\beta$ that models the T cell receptor binding affinity [refer to equation (2.7)].

3.2 Main Results

Variability in the initial tumour antigen profiles determines the effectiveness of the immune response

To investigate how the immune response is affected by variability in the initial tumour antigen profiles, we test for three increasing values of the parameter $V_{T}$ (i.e. $V_{T}=0.001$, $V_{T}=0.01$ and $V_{T}=0.1$). We choose $\beta=0.01$ and we let the antigen profiles of the tumour cells remain constant over time (i.e. we choose $\theta_{E}=\theta_{M}=0$). For each value of $V_{T}$ considered, we also explore the effect of increasing the value of the parameter $V_{D}$. In all cases, we carry out numerical simulations for $1000$ time-steps. As shown by Figure 3a, for a low value of $V_{T}$, very few tumour cells remain in the system after $1000$ time-steps for all considered values of $V_{D}$. Conversely, the results presented in Figure 3c show that if we set a relatively large value for $V_{T}$, there is a significant number of remaining tumour cells after $1000$ time-steps for all values of $V_{D}$. Moreover, as shown by Figure 3b, for an intermediate value of $V_{T}$, there appears to be a correlation between the number of tumour cells remaining after $1000$ time-steps and the parameter $V_{D}$. In particular, larger values of the parameter $V_{D}$ correspond to smaller numbers of the remaining tumour cells after $1000$ time-steps. These results suggest that for tumours characterised by intermediate levels of initial antigenic heterogeneity between tumour cells, higher deviation between the antigen profile recognised by DCs and the actual antigen profile of tumour cells may result in a more effective immune response. This is further illustrated by the computational results presented in the next paragraph.

Increasing variations between the antigen profile recognised by DCs and the actual antigen profile of tumour cells can result in immune escape, chronic dormancy or immune clearance of the tumour

The results discussed in the previous paragraph illustrate how different cell dynamics can be observed for increasing values of the parameter $V_{D}$. We test this further by using the parameter setting of Figure 3b (i.e. $V_{T}=0.01$, $\theta_{E}=0$ and $\beta=0.01$) and comparing the dynamics obtained for three different values of $V_{D}$ (i.e. $V_{D}=0.001$, $V_{D}=0.05$ and $V_{D}=0.1$). In Figure 4, we compare the average antigen profile of the tumour cells at the end of simulations with the average antigen profile recognised by the DCs, and we show the corresponding time evolution of the number of tumour cells. The insets also display the spatial cell distributions observed at the end of simulations to allow for a clearer understanding of the resulting dynamics. We observe that $V_{D}$ is a bifurcation parameter whereby three distinct situations result from choosing increasing values of $V_{D}$. In particular, Figures 4a,d refer to the case where the value of $V_{D}$ is relatively low (i.e. $V_{D}=0.001$), and show that there is very little difference between the average antigen profile of the tumour cells and the average recognised antigen profile at the end of simulations. Moreover, after an initial decrease, the tumour cell number increases steadily over time resulting in a relatively large final number of tumour cells. Furthermore, Figures 4b,e refer to the case where an intermediate value of $V_{D}$ is considered (i.e. $V_{D}=0.05$), and show that there a is larger variation between the average antigen profile of the tumour cells and the average recognised antigen profile at the end of simulations. Additionally, after a steep decrease, the tumour cell number remains at a low, almost constant, level for the remainder of the simulation time interval. Finally, Figures 4c,f refer to the case where the value of $V_{D}$ is relatively large (i.e. $V_{D}=0.1$), and show that the difference between the average antigen profile of the tumour cells. Moreover, the average recognised antigen profile at the end of simulations is even more varied than in the previous cases and the number of tumour cells decreases steadily over time until eventually the tumour is completely removed.

Increasing the T cell receptor binding affinity can benefit the immune system response to cancer

To explore the effect of altering the T cell receptor binding affinity, we test for three increasing values of the parameter $\beta$ (i.e. $\beta=0$, $\beta=0.0001$ and $\beta=0.001$). We choose $V_{T}=0.0001$ and we let the tumour cell antigen profiles remain constant over time (i.e. we choose $\theta_{E}=\theta_{M}=0$). For each value of $\beta$ considered, we also explore the effect of increasing the value of the parameter $V_{D}$. In all cases, we carry out numerical simulations for $1000$ time-steps. As shown by Figure 5a, for $\beta=0$, a considerable number of tumour cells remain inside the system at the end of simulations for all values of $V_{D}$. Conversely, the results presented in Figure 5c show that when $\beta$ is sufficiently high, very few tumour cells remain in the system after $1000$ time-steps for all values of $V_{D}$. Moreover, as shown by Figure 5b, for intermediate values of $\beta$, there appears to be a correlation between the number of tumour cells remaining after $1000$ time-steps and the parameter $V_{D}$. In particular, larger values of the parameter $V_{D}$ correspond to smaller numbers of the remaining tumour cells after $1000$ time-steps. These results suggest that the T cell receptor binding affinity plays a key role in the immune response to tumour cells. In the same way as Figure 4, Figure 6 shows that, under the parameter choice of the computational simulations related to Figure 5b, increasing the value of the parameter $V_{D}$ leads to immune escape or chronic dormancy or immune clearance of the tumour.

Increasing the probability of epimutations can lead to variations in the immune response to tumour cells

So far, we have considered only the situation where the antigen profile of each tumour cell remains constant over time (i.e. the probability with which epimutations and mutations leading to antigenic variations occur are $\theta_{E}=0$ and $\theta_{M}=0$). To investigate the effect of epimutations on the success of the immune response against tumour cells, we consider the parameter setting that we have used in the computational simulations shown either in Figures 4a,d or in Figures 6c,f but now we allow the antigen profiles of tumour cells to change through epimutations (i.e. we choose $\theta_{E}>0$). We consider eight distinct values of $\theta_{E}$ defined as fractions or multiples of the average probability of epimutations $\theta_{E}^{*}$, given in Table 1. The range we consider, $0\leq\theta_{E}\leq 0.92$, is chosen to include the range of De Smet et al. (1996), $012\leq\theta_{E}\leq 0.45$. In all cases, we carry out numerical simulations for $8640$ time-steps and we report on tumour cell numbers obtained as the average over $5$ simulations. Figure 7a displays the time evolution of the tumour cell number for the parameter setting of Figures 4a,d. These results show that increasing values of $\theta_{E}$ correspond to decreasing numbers of tumour cells inside the system at the end of simulations. In summary, by increasing the probability of epimutations the dynamics of tumour cells change from immune escape, through to chronic dormancy to immune clearance. On the other hand, Figure 7b displays the time evolution of the tumour cell number for the parameter setting of Figures 6c,f. These results show that for sufficiently small values of $\theta_{E}$ there are no tumour cells left inside the system at the end of simulations, whereas for larger values of $\theta_{E}$ a small number of tumour cells persist at the final time-step. Generally, by increasing the probability of epimutations the dynamics of tumour cells change from immune clearance to chronic dormancy.

Mutations have a weaker impact on the immune response to tumour cells compared to epimutations

We now compare the impact of mutations and epimutations on the immune response to tumour cells. Following what we have done in the previous paragraph, we consider the parameter setting used in the computational simulations shown either in Figures 4a,d or in Figures 6c,f but now we allow the antigen profiles of tumour cells to change through mutations (i.e. we choose $\theta_{M}>0$). We consider eight distinct values of $\theta_{M}$ defined as fractions or multiples of the average probability of epimutations $\theta_{E}^{*}$, given in Table 1. In all cases, we carry out numerical simulations for $8640$ time-steps and we report on tumour cell numbers obtained as the average over $5$ simulations. Figure 8a displays the time evolution of the tumour cell number for the parameter setting of Figures 4a,d and shows that immune escape occurs for all values of $\theta_{M}$ considered. On the other hand, Figure 8b displays the time evolution of the tumour cell number for the parameter setting of Figures 6c,f and shows that immune clearance occurs for all values of $\theta_{M}$ considered. Comparing these results with those displayed in Figure 7a and Figure 7b, respectively, we reach the conclusion that mutations have a weaker impact on the immune response to tumour cells compared to epimutations.

4 Discussion and conclusions

Spatial interactions between cancer and immune cells, as well as the recognition of tumour antigens by cells of the immune system, play a key role in the immune response against solid tumours. The existing mathematical models generally focus only on one of these key aspects. We have presented here a spatially explicit stochastic individual-based model that incorporates the adaptive processes driving tumour antigen recognition. Our model takes explicitly into account the dynamical heterogeneity of tumour antigen expression, and effectively captures the way in which this affects the immune response against the tumour.

Our computational simulation results show that the initial antigen expression profiles of cancer cells within the tumour have a crucial impact upon the outcome of the immune response [refer to Figure 3]. In the situation of an almost homogeneous tumour (i.e. where all tumour cells have a similar antigen profile), immune clearance occurs. Conversely, when the antigenic composition between cancer cells is highly heterogeneous the tumour may be able to escape the immune system response and continue growing. Interestingly, for moderate levels of initial antigenic heterogeneity our results demonstrate that the fate of the tumour is determined by the specificity of the cellular immune response.

The computational outcomes of our model indicate that the parameter controlling the specificity of the antigen recognition process of the dendritic cells (i.e. the parameter $V_{D}$) ultimately dictates which receptors are produced by the cytotoxic T lymphocytes [refer to Figures 3 and 4]. A larger value of this parameter brings about a more diverse receptor repertoire of the CTLs, and in turn results in a better immune response. This suggests that it is advantageous for the T cell pool to be multi-specific, whereby several different antigen receptors are simultaneously expressed by the CTL pool. In this respect, the outcomes of our model recapitulate the conclusions of experimental papers showing the success of a more diverse T cell repertoire in response to cancer (Carreno et al., 2015; Gerdemann et al., 2011; Ott et al., 2017; Sahin et al., 2017; Schumacher and Hacohen, 2016; Sharma and Allison, 2015). We remark that new experimental techniques have recently been developed to alter the specificity of T cell receptors (Smith et al., 2014). One particular approach is to use gene editing, in vitro, to modify which antigen the T cell receptors will respond to (Albers et al., 2019).

Moreover, our numerical results support the idea that varying the specificity of the immune response can result in three distinct scenarios, from immune escape, through to chronic dormancy to immune clearance of the tumour [refer to Figures 4 and 6]. The importance of tumour dormancy controlled by the immune system (i.e. immunological dormancy) has been highlighted by previous experimental and theoretical work (Kuznetsov et al., 1994; Lorenzi et al., 2015; Matzavinos et al., 2004; Wu et al., 2018). In particular, immunological dormancy can explain situations where there is an extended period of time before the occurrence of tumour relapse (Aguirre-Ghiso, 2007; Gomis and Gawrzak, 2017; Manjili, 2018; Teng et al., 2008; Wang and Lin, 2013; Yeh and Ramaswamy, 2015). In this regard, our model suggests the existence of a possible relationship between the specificity of the immune response and the emergence of prolonged immunological dormancy.

We have also explored the way in which altering the binding affinity of the CTLs to their corresponding tumour antigen may change the immune response to the tumour. Our results indicate that a stronger binding affinity leads to a more effective immune response, as the CTLs have a wider range of tumour cells that they can interact with [refer to Figure 5]. Previously, Gerdemann et al. (2011) found experimentally that a strong T cell binding affinity to tumour antigens played a key role in the overall immune response to the disease. Integrating the outcomes of our model and such experimental findings suggests that enhancing the binding affinity of T cells – e.g. through the modification of the receptors that the T cells of a patient can produce – could be a potential target of adoptive T cell therapy.

The results from our computational simulations suggest that changes in the antigenic expression of tumour cells due to epimutations can be either beneficial or detrimental to the immune response to a solid tumour [refer to Figure 7]. In more detail, we have found that in some cases increasing the probability of epimutations could transform situations of immune escape into tumour dormancy and eventually tumour removal. These findings are interesting in light of cancer therapy as they suggest that the efficacy of the immune response against solid tumours could be enhanced by increasing the probability of epimutations. In this respect, the loss of DNA methylation was the first epimutation to be identified in cancer cells (Feinberg and Tycko, 2004) and several experimental and clinical works found that the expression of MAGE antigens could be increased through demethylation (Chinnasamy et al., 2011; Gerdemann et al., 2011; Graff-Dubois et al., 2002; Wischnewski et al., 2006). Taken together, the outcomes of our model suggest that by combining a T cell therapy targeting multiple MAGE genes – e.g. using approaches similar to those of Gerdemann et al. (2011) – and increasing the probability of epimutations through demethylating agents – e.g. using methods similar to those of Chinnasamy et al. (2011) and Wischnewski et al. (2006) – a stronger immune response could be induced. However, in other cases, we have observed that increasing the probability of epimutations can turn instances of tumour removal into scenarios whereby a small number of tumour cells persisted over time. These contrasting results were also suggested previously through experimental research, where epimutations could be either beneficial or detrimental to tumour development (Chen and Mellman, 2017; Yarchoan et al., 2017).

We have additionally studied the effect of variations in the antigenic expression of tumour cells caused by mutations. In all parameter settings we have considered, increasing the probability of mutations did not change the resulting dynamics of the tumour-immune response [refer to Figure 8]. This suggests that mutations have a weaker effect on tumour-immune competition than antigenic variations caused by epimutations. This finding is coherent with experimental observations indicating that epimutations generally occur more frequently than mutations in tumour development (Feinberg, 2004; Peltomäki, 2012). Hence, our results demonstrate the importance of understanding the underlying causes of antigenic variations in tumour cells when considering tumour-immune competition.

Looking to the future, our individual-based model could be developed further in several ways. We could incorporate extended aspects of the tumour micro-environment, such as, proliferation of the immune cells and the interaction of tumour cells with abiotic factors (e.g. oxygen and glucose) that can affect the phenotypic composition of the tumour. Moreover, by posing the model on a 3D domain, further understanding could be obtained regarding the spatial dynamics of the tumour-immune response. The flexibility of our model would also allow for the inclusion of a wider range of antigens. Finally, although useful for investigating the qualitative and quantitative dynamics of a biological system, individual-based models like ours are not amenable to mathematical analysis. In a variety of contexts the benefit of relating stochastic individual-based models to deterministic continuum models has been highlighted (Champagnat et al., 2006; Chisholm et al., 2016, 2015; Deroulers et al., 2009; Painter and Hillen, 2015; Penington et al., 2011; Stevens, 2000). Combining such modelling approaches makes it possible to integrate computational simulations with analytical results, thus enabling a more extensive exploration of the model parameter space. This is a line of research that we will be pursuing in the near future by exploiting the formal methods that we have recently presented (Chaplain et al., 2018; Stace et al., 2019).

Acknowledgments

FRM is funded by the Engineering and Physical Sciences Research Council (EPSRC).