Comparing Stochastic Differential Equations and Agent-Based Modelling and Simulation for Early-Stage Cancer

There is great potential to be explored regarding the use of agent-based modelling and simulation as an alternative paradigm to investigate early-stage cancer interactions with the immune system. It does not suffer from some limitations of ordinary differential equation models, such as the lack of stochasticity, representation of individual behaviours rather than aggregates and individual memory. In this paper we investigate the potential contribution of agent-based modelling and simulation when contrasted with stochastic versions of ODE models using early-stage cancer examples. We seek answers to the following questions: (1) Does this new stochastic formulation produce similar results to the agent-based version? (2) Can these methods be used interchangeably? (3) Do agent-based models outcomes reveal any benefit when compared to the Gillespie results? To answer these research questions we investigate three well-established mathematical models describing interactions between tumour cells and immune elements. These case studies were re-conceptualised under an agent-based perspective and also converted to the Gillespie algorithm formulation. Our interest in this work, therefore, is to establish a methodological discussion regarding the usability of different simulation approaches, rather than provide further biological insights into the investigated case studies. Our results show that it is possible to obtain equivalent models that implement the same mechanisms; however, the incapacity of the Gillespie algorithm to retain individual memory of past events affects the similarity of some results. Furthermore, the emergent behaviour of ABMS produces extra patters of behaviour in the system, which was not obtained by the Gillespie algorithm.


Introduction
In previous work, three case studies using established mathematical models of immune interactions with early-stage cancer were considered in order to investigate the additional contribution of ABMS to ODE models simulation [1]. These case studies were re-conceptualised under an agent-based perspective and the simulation results were compared with those from the ODE models. Our results showed that, apart from the well known differences between these approaches (as those outlined for example, in Schieritz and Milling [2]), further insight from using ABMS was obtained, such as extra population patterns of behaviour.
In this work we apply the Gillespie algorithm [3,4], which is a variation of the Monte Carlo method, to create stochastic versions of the original ODE models investigated in [1]. We aim to reproduce the variability embedded in the ABMS systems to the mathematical formulation and verify whether results resemble each other. In addition, due to the fact that the Gillespie algorithm also regards integer quantities for their elements, we hope that this method overcomes some differences observed when comparing atomic agents represented in the ABMS with possible fractions of elements observable in the ODE results.
To the best of our knowledge, current literature regarding the direct comparison of the Gillespie algorithm and ABMS is scarce.
We want therefore to answer research questions such as: (1) Does this new stochastic formulation produce similar results to ABMS? (2) Can these methods be used interchangeably for our case studies? (3) Does the stochastic model implemented using the Gillespie algorithm also find the extra patterns revealed by the ABMS? We aim to establish a methodological discussion regarding the benefits of each approach for biological simulation, rather than provide further insights into the biological aspects of the problems studied. We therefore intend to compare the dynamics of each approach and observe the outcomes produced over time. We hope that this study provides further insights into the potential usability and contribution of ABMS to systems simulation.

Case Studies
The case studies used for our comparison regard models with different population sizes, varying modelling effort and model complexity. We hope that by tackling problems with different characteristics, a more robust analysis of our experiments is performed. The features of each case study are shown in Table 1.
The first case study considered is based on an ODE model involving interactions between generic tumour and effector cells. The second case study adds to the previous model the influence of the IL-2 cytokine molecules in the immune responses. The third case study comprises a model of interactions between effector cells, tumour cells, and IL-2 and TGF-b molecules. For all case studies, three approaches are presented: the original mathematical model, its conversion into the Gillespie algorithm model and the ABMS model. To answer our research questions, the Gillespie and ABMS approaches outcomes are compared. Our results show that for most cases the Gillespie algorithm does not produce outcomes statistically similar to the ABMS. In addition, Gillespie is incapable to reproduce the extreme patterns observed in the ABMS outcomes for the last case study.
The remainder of this paper is organized as follows. The next section introduces the literature review comparing stochastic ODE models and ABMS for different simulation domains. First, we show general work that has been carried out in areas such as economics and operations research, and then we focus on research concerned with the comparison for immunological problems. Finally, we discuss gaps in the literature regarding cancer research. In the following section we introduce our Gillespie and agentbased modelling development processes and the methods used for conducting the experimentation. Subsequently we present our case studies, comparison results and discussions. In the final section we draw our overall conclusions and outline future research opportunities.

Related Work
Current in-silico approaches used in early-stage cancer research include computational simulation of compartmental models, individual-based models and rule-based models. Compartmental models adopt an aggregate representation of the elements in the system. They include deterministic methods such as ordinary differential equation (ODE) models, system dynamics (SD) models and partial differential equation (PDE) models. These models have been largely employed in the study of dynamics between cancer cells and tumour cells [5,6], therapies for cancer [7], tumour responses to low levels of nutrients [8][9][10][11] and tumour vascularization [12,13]. Although these models have been very useful to understand and uncover various phenomena, they present several limitations. For instance, they do not encompass emergent behaviour and stochasticity. In addition, it is difficult to keep a record of individual behaviour and memory over the simulation course [14,15]. Stochastic compartmental models include Monte Carlo simulation models, which are computational algorithms that perform random sampling to obtain numerical results [16]. Amongst others, they are useful for simulating biological systems, such as cellular interactions and the dynamics of infectious diseases [17]. As these methods rely on stochastic process to produce their outputs, they overcome some of the limitations of the deterministic compartmental models, as they allow for variability of outcomes. The individuals in these models, however, do not have any sort of memory of past events. Rule-based models are a relatively new research area mostly focused on modelling and simulating biochemical reactions, molecular interactions and cellular signalling. The literature regarding the application of rule-based models to interactions between the immune system and cancer cells, however, is scarce. Individual-based models, or agent-based modelling and simulation (ABMS), relax the aggregation assumptions present in compartmental approaches and allow for the observation of the behavior of the single cells or molecules involved in the system. This approach has also been applied to early-stage cancer research [1,18].
The differences between deterministic compartmental models and individual-based models are well known in operations research [2,[19][20][21] and have also been studied in epidemiology [22,23] and system's biology [24][25][26]. Deterministic compartmental models assume continuous values for the individuals in the system, whereas in ABMS individual agents are represented. This peculiarity of each approach highly impacts the simulation results similarity depending on the size of the populations [1]. There is still however the need for further investigations between the interchangeable use of some Monte Carlo methods and ABMS.

Approaches Comparison
As mentioned previously, there are few studies that compare the Gillespie algorithm with ABMS. Most of these studies regard research in economic models and immunology. To the best of our knowledge there is no literature regarding the direct comparison of these methods to early-stage interactions between the immune system and tumour cells. This section describes relevant researches in several areas, which provided further insights into the gaps in the current literature and the research questions addressed in this paper.
There are a few attempts of re-conceptualizing agent-based models into simpler stochastic models of complex systems in economics. For instance, Daniunas et al. [27] start from simple models with established agent-based versions (which they named ''the model's microscopic version'') and try to obtain an equivalent macroscopic behavior. They consider microscopic and macroscopic versions of the herding model proposed by Kirman [28] and the diffusion of new products, proposed by Bass in [29]. They conclude that such simple models are easily replicated in a stochastic environment. In addition, the authors state that for the economics field, only very general models, such as those studied in their article, have well established agent-based versions and can be described by stochastic or ordinary differential equations. However, as the complexity of the microscopic environment increases, it becomes challenging to obtain resembling results with stochastic simulations and further developments need to be pursuit. Furthermore, the authors debate that the ambiguity present in the microscopic description in complex systems is an objective obstacle for quantitative modeling and needs further studying.
Stracquadanio et al. [30] investigate the contributions of ABMS and the Gillespie method for immune modelling. The authors, however, do not apply both methods to the same problem. Instead, for the first approach, they chose to investigate a largescale model involving interactions of immune cells and molecules. This model's objective was to simulate the immune elements interplay over time. For the Gillespie approach, the authors investigate a stochastic viral infection model. The authors point out three factors that play a major role in the modeling outcome when comparing ABMS and Gillespie: simulation time, model precision and accuracy, and model applicability. Regarding time, the authors state that stochastic models implemented with the Gillespie algorithm are preferred. On the other hand, ABMS permits more control over simulation runtime as it keeps record of the behavior of each single entity involved in the system. Regarding applicability, the authors argue that traditional Gillespie methods do not account for spatial information, which can be detrimental to the model accuracy given the fact that many immune interactions occur within specific spatial regions of the simulation environment.
Karkutla [31] compares two biological simulators: GridCell, which is a stochastic tool based on Gillespie's, and his new developed ABMSim, which is a simulation tool based on ABMS. GridCell was developed to overcome the issues in traditional Gillespie's, as pointed out by Stracquadanio et al. [30]. It is a stochastic tool able to tackle non-homogeneity effectively by addressing issues of crowding and localization. In GridCell, however, the problem of tracking individual behaviour and determining particular characteristics to each element still exists. ABMSim was therefore developed to overcome these issues. GridCell has been compared with ABMSim qualitatively and quantitatively and the two tools have produced similar results in the author experiments.
In our work we further study the differences between the approaches outcomes and investigate whether under different problem characteristics for early-stage cancer we still obtain

Methods
This section introduces the research methodology used for the development of our simulation models and for the experimentation performed in the following sections. As mentioned previously, our investigations concern the use of three cases to answer research questions regarding the application of the Gillespie algorithm and ABMS interchangeably for early-stage cancer models. For each case study, there is a well established ODE model from the literature and its correspondent agent-based model, that we previously developed in [1]. The ODE model simulations are implemented using the ODE solver module from MATLAB (2011) The Gillespie algorithm is implemented by the direct conversion of the original mathematical equations into reactions and simulating them under the COPASI 4:8 (Build35) simulator environment. The method used for the stochastic simulations is the Gillespie algorithm adapted using the next reaction method [32], with interval sizes of 0:1, integration interval between 0 and 1 and maximum internal steps of 10 6 .
ABMS is a modelling and simulation technique that employs autonomous agents that interact with each other. The agents' behaviour is described by rules that determines how they learn, interact and adapt. The overall system behaviour is given by the agents individual dynamics as well as their interactions. Our agentbased models are implemented using the AnyLogic TM 6.5 educational version (XJ Technologies 2010) [33]. This approach is developed by using state charts and tables containing each agent description. The state charts show the different possible states of an entity and define the events that cause a transition from one state to another. In order to facilitate the understanding of the agentbased model, we reproduce here the models developments, which were based on [1] (The ABMS and Gillespie models are available for download in http://anytips.cs.nott.ac.uk/wiki/index.php/ Resources).

Methodology for Results Comparison
As the Gillespie and ABMS are both stochastic simulation methods, we ran five hundred replications for each case study and calculated the mean values for the outputs. For all approaches, the rates (for cellular death, birth, etc.) employed were the same as those established by the mathematical model.
In addition, in order to investigate any statically significant differences between the ABMS and Gillespie techniques for the case studies, we implement a mixed effect model. This is a type of regression that considers both fixed and random effects. This method accounts for correlation caused by repeating the measure over time (i.e., the tumour cell count is correlated over time for each simulation run). The mixed effect analysis was implemented in the programming language R using the package NLME [34].
As a mixed effect model requires finding parameters for a regression model, it is not suitable when considering the whole time period. This is because in cases 2 and 3 the tumour dynamics has a damping oscillation and the function describing this dynamics is unknown (see pages 11 and 14). Instead, the sequence of local maxima and minima are used. It can been seen that these are converging and any statistical deviation between these sequences for the different simulation techniques indicate differences between the output of the techniques. If the simulations from the ABMS and the Gillespie technique come from the same distribution, then there would be no statistical difference between the maxima and minima over time. Therefore, we investigate two  null hypotheses. The first is that the function of local maxima is the same for the ABMS and the Gillespie algorithm simulations. And the second is that the function of local minima is the same for the ABMS and the Gillespie algorithm simulations. We use a 1% significance level.
There is no standard technique for estimating the required sample size for non-linear mixed effect models for a defined power when the measure of effect is known [35]. Therefore, the simulations are run 500 times as this will increase the statistical power and increase the probability of a true positive in the statistical analysis. A false negative is still possible if there is only a small effect size, but if the effect is small, it is of less interest.

Case 1: Interactions between Tumour Cells and Generic Effector Cells
The first case considers tumour cells growth and their interactions with general immune effector cells, as defined in [8]. According to the model, effector cells search and kill the tumour cells inside the organism. They proliferate proportionally to the number of existing tumour cells. As the quantities of effector cells increase, their capacity of eliminating tumour cells is augmented. Immune cells proliferate and die per apoptosis, which is a programmed cellular death. In the model, cancer treatment is also considered and it consists of injections of new effector cells in the organism.
Mathematically, the interactions between tumour cells and immune effector cells are defined as follows [8]:   Table 2 shows the mathematical equations converted into reactions and their respective rate laws per cell.
In the agent-based model there are two classes of agents, the tumor cells and the effector cells, as described in [1]. Table 3 shows the parameters and behaviours corresponding to each agent state. For our agents, state charts are used to represent the different states each entity is in. In addition, transitions are used to indicate how the agents move from one state to another. Events are also employed and they indicate that certain actions are scheduled to occur in the course of the simulation, such as injection of treatment. The state chart representing the tumour cells is shown in Figure 1(a), in which an agent proliferates, dies with age or is killed by effector cells. In addition, tumour cells  contribute to damage to effector cells, according to the same rate as defined by the mathematical model (Table 4). Figure 1(b) shows the state chart for the effector cells. In the figure, the cell is either alive or dead by age or apoptosis. While the cell is alive, it is also able to kill tumour cells and proliferate. In the transition rate calculations, the variable TotalTumourCells corresponds to the total number of tumour cell agents; and the variable TotalEffectorCells is the total number of effector cell agents. In the simulation model, apart from the agents, there is also an event -namely, treatment -which produces new effector cells with a rate defined by the parameter s. Experimental design for the simulations. Similarly to the experiments from [1], four scenarios are investigated. The scenarios have different rates for the death of tumour cells (defined by parameter b), effector cells apoptosis (defined by parameter d) and different treatments (parameter s). The values for these parameters are obtained from [6] (Table 5). In the first three scenarios, cancer treatment is considered, while the fourth case does not consider any treatment. The simulations for the ABMS and the Gillespie algorithm are run five hundred times and the mean values are displayed as results.
Results and discussion. Figure 2 shows the results of our experiments. In the first column we display the results from the ODE model for guidance. The second column shows the results from the Gillespie algorithm and the third column presents the ABMS results. Each row of the figure represents a different scenario.
Results for Scenario 1 appear similar for the three approaches, although the effector cells curve from the ABMS show more variability. To evaluate whether the results are significantly different for the two simulation methods, we apply a mixed effect model. The null hypothesis is that there is no significant difference between the methods (and therefore there will be no significant fixed effect for the method type). We use a 1% significance level. We are testing the similarity for the population of effector cells.
The effector cells follow a dynamic similar to 1/x for the time between 1 and 100: We apply a mixed-effect model where the simulation run is considered to have a random effect on the parameter a and the simulation method has a fixed effect on a and b. The results are presented in Table 6. At a 1% significance level the results of the two techniques are significantly different as the p-values for the fixed effect of the method on the parameters are less than 0:1.
We believe that the variability observed in the ABMS and Gillespie curves, given their stochasticity, also influenced the statistical test results. The number of effector cells for all simulations follow a similar pattern, although the similarity hypothesis was rejected. This variability of Gillespie and ABMS is very evident with regards to the effector cells population as the size of the populations involved in the first scenario is relatively small, which increases the impacts of stochasticity in the outcomes.
Results for scenario 2 are shown in the second row of Figure 2. The outcomes seem fairly different. By observing the ODE results, during about the first ten days, the tumour cells decrease and then grow up to a value of about 240 cells, subsequently reaching a  steady-state. This initial decrease is also observed in both Gillespie and ABMS curves. However, only the Gillespie method shows a similar increase in the numbers of tumour cells when compared to the ODEs. Similarly to the previous scenario, the Gillespie and ABMS simulation curves present an erratic behaviour throughout the simulation days. There is, however, an unexpected decay of tumour cells over time in the ABMS simulation, which does not happen in the Gillespie outcome. We believe the difference observed in the ABMS is due to the the individual characteristics of the agents and their growth/death rates attributed to their instantiation. While both ODEs and Gillespie are compartmental models and therefore they apply the model rates to the cells population, ABMS on the contrary, employs these rates in an individual basis. As the death rates of the tumour cells agents are defined according to the mathematical model, when the tumour cell population grows, the newborn tumour cells have higher death probabilities, which leads to a considerable number of cells dying out. This indicates that the individual behaviour of cells can lead to a more chaotic behaviour when compared to the aggregate view observed in the compartmental simulation. For scenarios 3 and 4, shown in third and fourth rows of Figure 2, respectively, the results for the three approaches differ completely. The differences are even more evident for the tumour cells outcomes. The ODEs results for scenario 3 reveal that tumour cells decreased as effector cells increased, following a predator-prey trend curve. For the ABMS, however, the number of effector cells decreased until a value close to zero was reached, while the tumour cells numbers were very different from those in the ODEs results. The ODE pattern noticed was possible given its continuous character. In the ODE simulation outcome curve for the effector cells it is therefore possible to observe, for instance, that after sixty days the number of effector cells ranges between one and two. These values could not be reflected in the ABMS simulation, as it deals with integer values. Similarly, the Gillespie approach outcomes did not resemble those from the ODE model. There is more variability in the tumour cells curve than in the ABMS outcomes, although the number of tumour cells also reaches zero after around sixty days.
In the fourth scenario, although effector cells appear to decay in a similar trend for both approaches, the results for tumour cells vary largely. In the ODE simulation, the numbers of effector cells reached a value close to zero after twenty days and then increased to a value smaller than one. For the ABMS simulation, however, these cells reached zero and never increased again. For the Gillespie model results, a similar pattern as that from the ABMS model occurs, although there seems to be less variance in the outcome curve. In addition, the mean numbers for tumour cells for the Gillespie approach seem smaller that those observed in the ABMS.
Summary. An outcome comparison between an ABMS and a Gillespie algorithm model was performed for case study 1. We considered an ODE model of tumour cells growth and their interactions with general immune effector cells as the baseline for results validation. Four scenarios considering small population numbers were investigated and results from ABMS and Gillespie were different for both populations for all scenarios. Furthermore, both approaches differed largely from the original mathematical outcomes. These results indicate that, for this case study, the stochasticity applied to the population as a whole when compared to that applied to the individual has a higher impact given the small population sizes. The result analysis also reveals that conceptualizing the stochastic approaches from the mathematical equations does not always produce statistically similar outputs.

Case 2: Interactions between Tumour Cells, Effector Cells and Cytokines IL-2
Case two regards the interactions between tumour cells, effector cells and the cytokine IL-2. It extends the previous study as it considers IL-2 as molecules mediating the immune response towards tumour cells. These molecules interfere in the proliferation of effector cells, which occurs proportionally to the number of tumour cells in the system. For this case, there are two types of treatment, the injection of effector cells or the addition of cytokines.
The mathematical model used in case 2 is obtained from [9]. The model's equations described bellow illustrate the non-spatial  dynamics between effector cells (E), tumour cells (T) and the cytokine IL-2 (I L ): Equation 10 describes the rate of change for the effector cell population E [9]. Effector cells grow based on recruitment (cT) and proliferation ( ). The parameter c represents the antigenicity of the tumour cells (T) [5,9]. m 2 is the death rate of the effector cells. p 1 and g 1 are parameters used to calibrate the recruitment of effector cells and s1 is the treatment that will boost the number of effector cells.
Equation 11 describes the changes that occur in the tumour cell population T over time. The term a(1{bT) represents the logistic growth of T (a and b are parameters that define how the tumour cells will grow) and aaET g 2 zT is the number of tumour cells killed by effector cells. a a and g 2 are parameters to adjust the model.
The IL-2 population dynamics is described by Equation 12.
determines IL-2 production using parameters p 2 and g 3 . m 3 is the IL-2 loss. s2 also represents treatment. The treatment is the injection of IL-2 in the system. Table 7 shows the mathematical model converted into reactions for the Gillespie algorithm model. The first column of the table displays the original mathematical equation, followed by the equivalent reactions and rate laws in the subsequent columns.
As described in [1], the agents represent the effector cells, tumour cells and IL-2. Their behaviours are shown in Table 8. The state charts for each agent type are shown in Figure 3. The ABMS model rates are the same as those defined in the mathematical model and are given in Table 9. In the transition rate calculations, the variable TotalTumour corresponds to the total number of tumour cell agents, the variable TotalEffector is the total number of effector cell agents and TotalIL 2 is the total number of IL-2 agents. In the simulation model, apart from the agents, there are also two events: the first event adds effector cell agents according to the parameter s1 and the second one adds IL-2 agents according to the parameter s2.
Experimental design for the simulation. The experiment is conducted assuming the same parameters as those of the mathematical model (Table 10). For the ABMS and the Gillespie algorithm model, the simulation is run five hundred times and the average outcome value for these runs is collected. Each run simulates a period equivalent to six hundred days, following the same time span used for the numerical simulation of the mathematical model.
Results and discussion. The results obtained are shown in Figure 4 for tumour cells (left), effector cells (middle) and IL-2 (right), respectively. As the figure reveals, the results for all populations are analogous; the growth and decrease of all populations occur at similar times for all approaches. ABMS has a little more variability in the results, specially regarding IL-2. We believe that for this case, the large population sizes (around 10 4 )  produce a lower variability in the outcomes of the stochastic approaches. For statistical comparison, results were contrasted by applying a non-linear mixed effect model, as shown next. The local maxima sequence follows a second order polynomial function of the form: The local minima sequence follows a second order polynomial function of the form: For the mixed-effect model, we consider a and b to have fixed effects based on the type of simulation (e.g., ABMS or Gillespie algorithm) and a and b to have random effects based on the individual simulation run. The results of the mixed effect model  are presented in Tables 11 and 12. It can been seen that there is a significant difference between the a and b parameter values for the two different techniques. We therefore reject the null hypotheses and accept that there is a significance difference between the two techniques in terms of the sequence of maxima and sequence of minima, at a 1% significance level. Furthermore, the results show that the ABMS simulations tend to have larger local maxima and smaller local minima, which is clear in Figure 5. Summary. Interactions between tumour cells, effector cells and the cytokine IL-2 were considered to investigate the potential differences and similarities of ABMS and Gillespie algorithm outcomes. Statistical comparison between the Gillespie and the ABMS results show a significant difference in the outcomes. Compared to the original ODE model used as validation, ABMS displayed a little more variability in the results, whereas the Gillespie algorithm followed mostly the same patterns as those produced by the ODEs for all populations in the simulation. As for these simulations a bigger number of individuals was required, it was also observed that, regarding the use of computational resources, ABMS was far more time-and memory-consuming than the Gillespie approach.

Case 3: Interactions between Tumour Cells, Effector Cells, IL-2 and TGF-b
Case study three comprises interactions between tumour cells and immune effector cells, as well as the immune-stimulatory and suppressive cytokines IL-2 and TGF-b [5]. According to the ODE model developed by Arciero et al. in [5], TGF-b stimulates tumour growth and suppresses the immune system by inhibiting the activation of effector cells and reducing tumour antigen expression.
The mathematical model is described by the differential equations below: Equation 15 describes the rate of change for the effector cell population E. According to [5], effector cells are assumed to be recruited to a tumour site as a direct result of the presence of tumour cells. The parameter c in cT 1zcS represents the antigenicity of the tumour, which measures the ability of the immune system to recognize  on the presence of the cytokine IL-2 and is decreased when the cytokine TGF-b is present. p 1 is the maximum rate of effector cell proliferation in the absence of TGF-b, g 1 and q 2 are halfsaturation constants, and q 1 is the maximum rate of antiproliferative effect of TGF-b.
Equation 16 describes the dynamics of the tumour cell population. The term aT 1{ T K À Á represents a logistic growth dynamics with intrinsic growth rate a and carrying capacity K in the absence of effector cells and TGF-b. The term aaET g 2 zT is the number of tumour cells killed by effector cells. The parameter a a measures the strength of the immune response to tumour cells. The third term p 2 ST g 3 zS accounts for the increased growth of tumour cells in the presence of TGF-b. p 2 is the maximum rate of increased proliferation and g 3 is the half-saturation constant, which indicates a limited response of tumour cells to this growth-stimulatory cytokine [5].   The kinetics of IL-2 are described in equation 17. The first term p 3 ET (g 4 zT)(1zaS) represents IL-2 production which reaches a maximal rate of p 3 in the presence of effector cells stimulated by their interaction with the tumour cells. In the absence of TGF-b, this is a self-limiting process with half-saturation constant g 4 S [5]. The presence of TGF-b inhibits IL-2 production, where the parameter a is a measure of inhibition. Finally, m 2 I represents the loss of IL-2.
dS dt~p Equation 18 describes the rate of change of the suppressor cytokine, TGF-b. According to [5], experimental evidence suggests that TGF-b is produced in very small amounts when tumours are small enough to receive ample nutrient from the surrounding tissue. However, as the tumour population grows sufficiently large, tumour cells suffer from a lack of oxygen and begin to produce TGF-b in order to stimulate angiogenesis and to evade the immune response once tumour growth resumes. This switch in TGF-b production is modelled by the term p 4 T 2 h 2 zT 2 , where p 4 is the maximum rate of TGF-b production and t is the critical tumour cell population in which the switch occurs. The decay rate of TGF-b is represented by the term m 3 S. Table 13 presents the Gillespie algorithm model used for our simulations. The model was obtained by converting the ODEs into reaction equations.
The agents established for the ABMS represent the effector cells, tumour cells, IL-2 and TGF-b populations, as described in [1]. The agents' behaviour is defined in Table 14. The state charts for each agent type are illustrated in Figure 6.
The ABMS model rates corresponding to the mathematical model are given in Table 15. In the transition rate calculations, the variable TotalTumour corresponds to the total number of tumour cell agents; the variable TotalEffector is the total number of effector cell agents, TotalIL 2 is the total number of IL-2 agents and TotalTGFBeta is the total TGF-b agents. This model does not include events.

Experimental Design for the Simulation
The experiment is conducted assuming the same parameters as those defined for the mathematical model (Table 16). Similarly to the previous case studies, for the ABMS and Gillespie models the simulation is run five hundred times and the average outcome value for these runs is displayed as result. Each run simulates a   Table 16.
Results and discussion. The mean results of 500 runs for the Gillespie algorithm and the ABMS contrasted with the ODE model are shown in Figure 7. The left graph in the figure presents the outcomes for tumour cells; the graph in the middle shows the outputs for effector cells; the graph on the right shows the mean IL-2 outcomes (the TGF-b results have some particularities and therefore are discussed next). The figure shows that both Gillespie and ABMS do not match properly the original results from the mathematical model. Additionally, ABMS is far more dissimilar than what was anticipated. In order to understand why the mean values were that much different from what was expected, we plotted fifty individual runs for each approach, as shown in figures 8, 9, 10 and 11. These runs illustrate the variations observed in both ABMS (left side of the figures) and Gillespie (right side of the figures) approaches, due to its stochastic character. In the figures, the ODE model results were also plotted (dashed black line) in order to highlight the range of variation produced by the stochastic approaches. As it can be observed in the figures, both Gillespie and ABMS outcomes produce various slightly distinct starting times for the growth of populations. In addition to these variations, for a few runs the populations in ABMS decreased to zero, as previously reported in [1]. This behaviour was not reflected in the Gillespie algorithm results. This indicates that it is not always possible to replicate similar results within both approaches.
The use of ABS modelling has therefore led to the discovery of additional ''rare'' patterns, which we would have not been able to derive by using analytical methods or the dynamic Monte Carlo  method, i.e. the Gillespie algorithm. These ''extreme cases'' found by ABMS suggest that there might be circumstances where the tumour cells are completely eliminated by the immune system, without the need of any cancer therapies.
We believe that ABMS when compared to Gillespie produces extra patterns because of the agents individual behaviour and their interactions. While ODEs and the Gillespie algorithm always use the same values for the parameters over the entire population aggregate, ABMS rates vary with time and number of individuals. Each agent is likely to have distinct numbers for their probabilities and therefore have its own memory of past events (Gilespie, however, does not encompass individual memory for its elements). The agents individual interactions, which give raise to the overall behaviour of the system, are also influenced by the scenario determined by the random numbers used. By running the ABMS multiple times with different sets of random numbers, the outcomes vary according to these sets and the emerging interactions of the agents also produce the rare outcome patterns.
For further statistical comparison of the results that follow the same pattern of behaviour for ABMS and Gillespie, a mixed-effect model is used. In 236 of the 500 ABMS simulations the tumour cell population dies out early in time. The remaining 264 ABMS simulations (where the tumour cell population does not die out over the [0,600] time period) is compared with the 500 Gillespie simulations by a non-linear mixed effect model.
The local maxima sequence follows a second order polynomial function of the form: The local minima sequence follows a second order polynomial function of the form: For the mixed-effect model we considered a and b to have fixed effects based on the type of simulation (e.g., ABMS or Gillespie algorithm) and a and b to have random effects based on the individual simulation run. The results of the mixed effect model are presented in Tables 17 and 18. It can been seen that there is a significant difference between the a and b parameter values for the two different techniques. We therefore reject the null hypotheses and accept that there is a significance difference between the two techniques in terms of the sequence of maxima and sequence of minima, at a 1% significance level.
Tables 17-18 and Figure 12 show that the sequence of local maxima of the ABMS diverge from that of the Gillespie algorithm over time. In the ABMS the tumour cells tend to increase to a larger count than the Gillespie algorithm simulations causing the function of local maxima for the ABMS to be significantly greater than the Gillespie. A possible explanation for this, as mentioned previously, is the fact that agents have memory and therefore the rates (death, proliferation, etc) for a certain cell are determined by the cells (and their proportions) present in the system at the moment the cell was created. For the Gillespie algorithm, instead, the rates are applied globally to the entire population and remain constant over the simulation course. Consequently, for Gillespie, the individuals do not keep a record of the previous population dynamics. This explanation is supported by the observation that the function describing the sequence of local minima of the ABMS is significantly lower than the Gillespie algorithm over time, as the same argument would account for the ABMS simulations reaching lower levels.
Regarding the TGF-b outcomes, the ODEs results reveal numbers smaller than one ( Figure 11 on the right), which is not possible to achieve with the ABMS and the Gillespie algorithm. The simulation results regarding these molecules are therefore completely different for both stochastic approaches. By observing the multiple runs graph of ABMS, however, results indicate that the TGF-b grows at around 100 and 200 days, which resembles what occurs in the ODE simulation for the first two peaks of TGFb concentration. This suggests that ABMS, as opposite to Gillespie, is capable of capturing some of the behaviours of the analytical results even when the outputs are different. We believe that these observations need to be further investigated in order to determine whether this happens in other case studies. In addition, it is necessary to investigate in what circumstances and range of values ABMS is still capable to reflect behaviours of numbers smaller than one agent present in the ODE model.

Summary
The third case study simulations investigated interactions between effector cells, tumour cells and two types of cytokines, namely IL-2 and TGF-b. When compared to the original ODE results, both Gillespie algorithm and ABMS produced more variability in the outcomes. For each of the five hundred runs, a slightly different start of population growth was observed. In addition, ABMS produced extra patterns not observed in the original mathematical model and in the Gillespie results. These extra patterns have been reported previously in [1] and with the present work we wanted to find out whether the Gillespie algorithm simulation results would be as informative. This indicates that, for this case study, both methods should not be employed interchangeably, as some extra possible population patterns of behaviour might not be uncovered without ABMS. We believe that these emergent examples occur due to the individual interactions of the agents and their chaotic character. With these results, we answer our third research question that it is not possible in this case to obtain extreme patterns using the Gillespie algorithm.

Conclusions
In this work, we employed three case studies to investigate circumstances where we can use ABMS and the Gillespie algorithm interchangeably. We aimed at reproducing the variability embedded in the ABMS systems to the mathematical formulation and verify whether results resemble. Current literature  regarding the comparison of the Gillespie algorithm and ABMS is scarce and we wanted therefore to answer the questions: (1) Does the Gillespie algorithm produce similar results to ABMS? (2) Can these two methods be used interchangeably for our case studies? (3) Does the Gillespie algorithm also find the extra patterns revealed by the ABMS in the third case study? The case studies investigated regarded models with different characteristics, such as population sizes, modelling effort demanded and model complexity.
The first case study involved interactions with general immune effector cells and tumour cells. Four different scenarios regarding distinct sets of parameters were investigated and in the first three scenarios treatment was included. ABMS and Gillespie produced different results for all scenarios. It appears that two major characteristics of this model influenced the differences obtained: (1) The small quantities of individuals considered in the simulations (especially regarding the effector population size, which was always smaller than ten) that significantly increased the variability of both stochastic approaches; and (2) the stochasticity of the Gillespie algorithm is applied to the aggregates, while in the ABMS there is individual variability.
Case study 2 referred to the investigation of a scenario containing interactions between effector cells, cytokines IL-2 and tumour cells. For this case the Gillespie and ABMS approaches produced similar outcome curves, which also matched the pattern of behaviour of the mathematical model used for validation. As populations sizes had a magnitude of 104 individuals, the erratic behaviour of both stochastic approaches was no longer evident in the outcomes. However, although results seemed similar, further statistical tests reject their similarity hypothesis. It was observed that, for case 2 in general, ABMS simulation outcome curves tend to have larger local maxima and smaller local minima.
Case study 3 includes the influence of the cytokine TGF-b in the interactions between effector cells, cytokines IL-2 and tumour cells from the previous case. The simulation outcomes for the ABMS were mostly following the same pattern as those produced by the Gillespie algorithm, although the results were statistically different. In addition, Gillespie failed to replicate the alternative outcomes found by the ABMS. This indicates that for this case study the ABMS results are more informative, as they illustrate another set of possible dynamics to be validated in real-world. Furthermore, ABMS was also able to indicate two peaks where TGF-b concentrations have grown, although the corresponding values in the mathematical model were smaller than one.
In response to our research questions, we conclude that regarding the interchangeable use of Gillespie and ABMS, population size has a positive impact on result similarity. This means that bigger populations tend to result in close simulation output patterns. However, the stochasticity of both approaches and the memory present in the ABMS produce outcome differences which are statistically significant, although visually the outcomes look similar. Finally, the emergent behaviour of ABMS can contribute additional insight (extra patterns), which was not obtained by the aggregate stochasticity present in Gillespie, given its incapacity of retaining memory of past events for their elements.
Currently, it is acknowledged that the necessary levels of detail present a system's conceptual model might require the adoption of multiple simulation approaches, under different scales of abstraction, in order to produce satisfactory outcomes. For certain cases, multi-scale and multi-modelling solutions represent the only possible manner to properly implement a set of requirements, even though these practices are generally more resource-intensive.
Several models, especially those regarding biological systems, encompass elements with distinct time and length scales coupled together. Multiple layers involving different levels of granularity need therefore to be considered. Phenomena such as cellular gradient-based chemotaxis and nutrient-dependent cell cycles are better implemented by employing numerical simulation (such as partial and ordinary differential equations). On the other hand, cellular random movement, memory, adaptability, explicit representation of space, heterogeneous populations, etc. represent suitable concepts to be represented within an agent-based approach. There are also subsets of characteristics in a system that can be implemented under more than one approach.
Despite of the well-known domains of applicability of each simulation technique, other factors such as the real-world data availability, the computing resources, the required granularity of the model, the time available to develop the simulation, the research questions and the expertise of the simulation developer also influence the decision on the paradigms to be adopted. By understanding the strengths and weaknesses of each approach, it is possible decide on which route to take and also the implications of this decision to the final outcome. We hope that our study results can assist simulation researchers on this debate.