10.1 Competitive Release Introduction

To demonstrate how the aforementioned principles and components of HAL are applied, we now look at a simple but complete example of a hybrid modeling experiment. We will be designing a model of pulsed therapy based on a publication from our lab [1]. We also showcase the flexibility that the modular component approach brings by displaying 3 different parameterizations of the same model side by side in a “multiwell experiment”.

As in [1], the presented model assumes two competeing tumor-cell phenotypes: a rapidly dividing, drug-sensitive phenotype and a slower dividing, drug-resistant phenotype. There is also a drug diffusible that enters the system through the model edges and is consumed over time by the tumor cells.

Every tick (timestep), each cell has a probability of death and a probability of division. The division probability is affected by phenotype and the availability of space. The death probability is affected by phenotype and the local drug concentration.

An interesting outcome of the experiment is that pulsed therapy is better at managing the tumor than constant therapy. The modular design of HAL allows us to test 3 different treatment conditions, each with an identical starting tumor (No drug, constant drug, and pulsed drug) Under pulsed therapy the sensitive population is kept in check, while still competing spatially with the resistant phenotype and preventing its expansion. The rest of the section describes in detail how this abstract model is generated.

Figure 2 provides a high level look at the structure of the code. Table 1 provides a quick reference for the built-in HAL functions in this example. Any functions that are used by the example but do not exist in the table are defined within the example itself and explained in detail below the code. Those fluent in Java may be able to understand the example just by reading the code and using the reference table.

Built-in framework functions and classes used in the code are highlighted in green to make identifying framework components easier


Figure 2: Example program flowchart. Yellow font indicates where functions are first called.

Table 1: HAL functions used in the example. Each function is a method of a particular object, meaning that when the function is called it can readily access properties that pertain to the object that it is called from.

10.2 Main Function

1    public static void main(String[] args) { 
2        //setting up starting constants and data collection 
3        int x = 100, y = 100, visScale = 5, tumorRad = 10, msPause = 0; 
4        double resistantProb = 0.5; 
5        GridWindow win = new GridWindow("Competitive Release", x * 3, y, visScale); 
6        FileIO popsOut = new FileIO("populations.csv", "w"); 
7        //setting up models 
8        ExampleModel[] models = new ExampleModel[3]; 
9        for (int i = 0; i < models.length; i++) { 
10            models[i] = new ExampleModel(x, y, new Rand()); 
11            models[i].InitTumor(tumorRad, resistantProb); 
12        } 
13        models[0].DRUG_DURATION = 0;//no drug 
14        models[1].DRUG_DURATION = 200;//constant drug 
15        //Main run loop 
16        for (int tick = 0; tick < 10000; tick++) { 
17            win.TickPause(msPause); 
18            for (int i = 0; i < models.length; i++) { 
19                models[i].ModelStep(tick); 
20                models[i].DrawModel(win, i); 
21            } 
22            //data recording 
23            popsOut.Write(models[0].Pop() + "," + models[1].Pop() + "," + models[2].Pop() + "\n"); 
24            if ((tick) % 100 == 0) { 
25                win.ToPNG("ModelsTick" + tick + ".png"); 
26            } 
27        } 
28        //closing data collection 
29        popsOut.Close(); 
30        win.Close(); 
31    }

We first look at the main function for a bird’s-eye view of how the program is structured.

Note: Source code elements highlighted in green are already built into HAL.

Defines all of the constants that will be needed to setup the model and display.
Creates a GridWindow of RGB pixels for visualization and for generating timestep PNG images. x*3, y define the dimensions of the pixel grid. X is multiplied by 3 so that 3 models can be visualized side by side in the same window. The last argument is a scaling factor that specifies that each pixel on the grid will be viewed as a 5x5 square of pixels on the screen.
Creates a file output object to write to a file called populations.csv
Creates an array with 3 entries to fill in with Models.
Fills the model list with models that are initialized identically. Each model will hold and update its own cells and diffusible drug. See the Grid Definition and Constructor section and the InitTumor Function section for more details.
Setting the DRUG_DURATION constant creates the only difference in the 3 models being compared. In models[0] no drug is administered (the default value of DRUG_DURATION is 0). In models[1] drug is administered constantly (DRUG_DURATION is set to equal DRUG_CYCLE). In models[2] drug will be administered periodically. See the ExampleModel Constructor and Properties section for the default values.
Executes the main loop for 10000 timesteps. See the ModelStep Function for where the Model tick is incremented.
Requires every iteration of the loop to take a minimum number of milliseconds. This slows down the execution and display of the model and makes it easier for the viewer to follow.
Loops over all models to update them.
Advances the state of the agents and diffusibles in each model by one timestep. See the Model Step Function for more details.
Draws the current state of each model to the window. See the Draw Model Function for more details.
Writes the population sizes of each model every timestep to allow the models to be compared.
Every 100 ticks, writes the state of the model as captured by the GridWindow to a PNG file.
After the main for loop has finished, the FileIO object and the visualization window are closed, and the program ends.

10.3 ExampleModel Constructor and Properties

1  public class ExampleModel extends AgentGrid2D { 
2    //model constants 
3    public final static int RESISTANT = RGB(0, 1, 0), SENSITIVE = RGB(0, 0, 1); 
4    public double DIV_PROB_SEN = 0.025, DIV_PROB_RES = 0.01, DEATH_PROB = 0.001, 
5            DRUG_DIFF_RATE = 2, DRUG_UPTAKE = 0.91, DRUG_TOXICITY = 0.2, DRUG_BOUNDARY_VAL = 1.0; 
6    public int DRUG_START = 400, DRUG_CYCLE = 200, DRUG_DURATION = 40; 
7    //internal model objects 
8    public PDEGrid2D drug; 
9    public Rand rng; 
10    public int[] divHood = MooreHood(false); 
12    public ExampleModel(int x, int y, Rand generator) { 
13        super(x, y, ExampleCell.class); 
14        rng = generator; 
15        drug = new PDEGrid2D(x, y); 
16    }	 		
17 }


This section covers how the grid is defined and instantiated.

The ExampleModel class, which is user defined and specific to this example, is built by extending the generic AgentGrid2D class. The extended grid class requires an agent type parameter, which is the type of agent that will live on the grid. To meet this requirement, the <ExampleCell> type parameter is added to the declaration.
Defines RESISTANT and SENSITIVE constants, which are created by the Util RGB function. These constants serve as both colors for drawing and as labels for the different cell types.
Defines all constants that will be needed during the model run. These values can be reassigned after model creation to facilitate testing different parameter settings. In the main function, the DRUG_DURATION variable is modified for the Constant-Drug, and Pulsed Therapy experiment cases.
Declares that the model will contain a PDEGrid2D, which will hold the drug concentrations. The PDEGrid2D can only be initialized when the x and y dimensions of the model are known, which is why we do not define them until the constructor function.
Declares that the Grid will contain a Random number generator, but take it in as a constructor argument to allow the modeler to seed it if desired.
Defines an array that will store the coordinates of a neighborhood generated by the MooreHood function. The MooreHood function generates a set of coordinates that define the Moore Neighborhood, centered around the (0,0) origin. The neighborhood is stored in the format [0102,...,0n,x1,y1,x2,y2,...,xn,yn] . The leading zeros are written to when MapHood is called, and will store the indices that the neighborhood maps to. See the CellStep function for more information.
Defines the model constructor, which takes as arguments the x and y dimensions of the world and a Random number generator.
Calls the AgentGrid2D constructor with super, passing it the x and y dimensions of the world, and the ExampleCell Class. This Class is used by the Grid to generate a new cell when the NewAgentSQ function is called.
The random number generator argument is assigned and the drug PDEGrid2D is defined with matching dimensions.

10.4 InitTumor Function

1    public void InitTumor(int radius, double resistantProb) { 
2        //get a list of indices that fill a circle at the center of the grid 
3        int[] tumorNeighborhood = CircleHood(true, radius); 
4        int hoodSize = MapHood(tumorNeighborhood, xDim / 2, yDim / 2); 
5        for (int i = 0; i < hoodSize; i++) { 
6            if (rng.Double() < resistantProb) { 
7                NewAgentSQ(tumorNeighborhood[i]).type = RESISTANT; 
8            } else { 
9                NewAgentSQ(tumorNeighborhood[i]).type = SENSITIVE; 
10            } 
11        } 
12    }

The next segment of code is a function from the ExampleModel class that defines how the tumor is first seeded after the ExampleModel is created.

The arguments passed to InitTumor function are the approximate radius of the circular tumor being created and the probability that each created cell will be of the resistant phenotype.
Sets the circleCoords array using the built-in CircleHood function, which stores coordinates in the form [01,02,...,0n,x1,y1,x2,y2,...xn,yn]. These coordinate pairs define a neighborhood of all squares whose centers are within the radius distance of the center (0,0) origin square. The leading 0s are used by the MapHood function to store the mapping indices. The boolean argument specifies that the origin will be included in this set of squares, thus making a completely filled circle of squares.
Uses the built-in MapHood function to map the neighborhood defined above to be centered around xDim/2,yDim/2 (the dimensions of the AgentGrid). The results of the mapping are written as indices to the beginning of the tumorNeighborhood array. MapHood returns the number of valid indices found, and this will be the size of the starting population.
Loops from 0 to hoodSize, allowing access to each mapped index in the tumorNeighborhood.
Samples a random number in the range (0 - 1] and compares to the resistantProb argument to set whether the cell should have the resistant phenotype or the sensitive phenotype.
Uses the built-in NewAgentSQ function to place a new cell at each tumorNeighborhood position. In the same line we also specify that the type should be either resistant or sensitive, depending on the result of the rng.Double() call.

10.5 ModelStep Function

1    public void ModelStep(int tick) { 
2        ShuffleAgents(rng); 
3        for (ExampleCell cell : this) { 
4            cell.CellStep(); 
5        } 
6        int periodTick = (tick - DRUG_START) % DRUG_CYCLE; 
7        if (periodTick > 0 && periodTick < DRUG_DURATION) { 
8            //drug will enter through boundaries 
9            drug.DiffusionADI(DRUG_DIFF_RATE, DRUG_BOUNDARY_VAL); 
10        } else { 
11            //drug will not enter through boundaries 
12            drug.DiffusionADI(DRUG_DIFF_RATE); 
13        } 
14    }

This section looks at the main step function which is executed once per tick by the Model.

The ShuffleAgents function randomizes the order of iteration so that the agents are always looped through in random order.
Iterates over every cell on the grid, and calls the CellStep function on every cell.
The GetTick function is a built-in function that returns the current Grid tick. The If statement logic checks if the tick is past the drug start and if the tick is in the right part of the drug cycle to apply drug. (See the Grid Definition and Constructor section for the values of the constants involved, the DRUG_DURATION variable is set differently for each model in the Main Function)
If it is time to add drug to the model, the built-in DiffusionADI function is called. The default Diffusion function uses the standard 2D Laplacian and is of the form: δCδt = D2C, where D in this case is the DRUG_DIFF_RATE. DiffusionADI uses the ADI method which is more stable and allows us to take larger steps than the 2D Laplacian can support. The additional argument to the DiffusionADI equation specifies the boundary condition value DRUG_BOUNDARY_VAL. This causes the drug to diffuse into the PDEGrid2D from the boundary.
Without the second argument the DiffusionADI function assumes a reflective boundary, meaning that drug concentration cannot escape or enter through the sides. Therefore the only way for the drug concentration to decrease is via consumption by the Cells. See the CellStep function section, line 6, for more information.

10.6 CellStep Function and Cell Properties

1  class ExampleCell extends AgentSQ2Dunstackable { 
2    public int type; 
4    public void CellStep() { 
5        //Consumption of Drug 
6        G.drug.Mul(Isq(), G.DRUG_UPTAKE); 
7        double deathProb, divProb; 
8        //Chance of Death, depends on resistance and drug concentration 
9        if (this.type == RESISTANT) { 
10            deathProb = G.DEATH_PROB; 
11        } else { 
12            deathProb = G.DEATH_PROB + G.drug.Get(Isq()) * G.DRUG_TOXICITY; 
13        } 
14        if (G.rng.Double() < deathProb) { 
15            Dispose(); 
16            return; 
17        } 
18        //Chance of Division, depends on resistance 
19        if (this.type == RESISTANT) { 
20            divProb = G.DIV_PROB_RES; 
21        } else { 
22            divProb = G.DIV_PROB_SEN; 
23        } 
24        if (G.rng.Double() < divProb) { 
25            int options = MapEmptyHood(G.divHood); 
26            if (options > 0) { 
27                G.NewAgentSQ(G.divHood[G.rng.Int(options)]).type = this.type; 
28            } 
29        } 
30    } 
31  }

We next look at how the Agent is defined and at the CellStep function that runs once per Cell per tick.

The ExampleCell class is built by extending the generic AgentSQ2Dunstackable class. The extended Agent class requires the ExampleModel class as a type argument, which is the type of Grid that the Agent will live on. To meet this requirement, we add the <EampleModel> type parameter to the extension.
Defines an internal int called type. each Cell holds a value for this field. If the value is RESISTANT, the Cell is of the resistant phenotype, if the value is SENSITIVE, the cell is of the sensitive phenotype. The RESISTANT and SENSITIVE constants are defined in the ExampleGrid as constants.
The G property is used to access the ExampleGrid object that the Cell lives on. G is used often with Agent functions as the AgentGrid is expected to contain any information that is not local to the Cell itself. Here it is used to get the drug PDEGrid2D. The drug concentration at the index that the Cell is currently occupying (Isq()) is then multiplied by the drug uptake constant, thus modeling local drug consumption by the Cell.
Defines deathProb and divProb variables, these will be assigned different values depending on whether the ExampleCell is RESISTANT or SENSITIVE.
If the cell is resistant the deathProb variable is set to the DEATH_PROB value alone, if the cell is sensitive, an additional term is added to account for the probability of the cell dying from drug exposure, using the concentration of drug at the cell’s position (Isq())
Samples a random number in the range (0-1] and compares to deathProb to determine whether the cell will die. If so, the built-in agent Dispose() function is called, which removes the agent from the grid, and then return is called so that the dead cell will not divide.
Sets the divProb variable to either DIV_PROB_RES for resistant cells, or DIV_PROB_SEN for sensitive cells.
Samples a random number in the range (0 - 1] and compare to divProb to determine whether the cell will divide.
If the cell will divide, the built-in MapEmptyHood function is used which displaces the divHood (the moore neighborhood as defined in the Grid Definition and Constructor section) to be centered around the x and y coordinates of the Cell, and writes the empty indices into the neighborhood. The MapEmptyHood function will only map indices in the neighborhood that are empty. MapEmptyHood returns the number of valid divison options found.
If there are one or more valid division options,a new daughter cell is created using the builitin NewAgentSQ function and its starting location is chosen by randomly sampling the divHood array to pull out one if its valid locations. Finally with the .type=this.type statement, the phenotype of the new daughter cell is set to the phenotype of the pre-existing daughter that remains in place, thus maintaining phenotypic inheritance.

10.7 DrawModel Function

1    public void DrawModel(GridWindow vis, int iModel) { 
2        for (int x = 0; x < xDim; x++) { 
3            for (int y = 0; y < yDim; y++) { 
4                ExampleCell drawMe = GetAgent(x, y); 
5                if (drawMe != null) { 
6                    vis.SetPix(x + iModel * xDim, y, drawMe.type); 
7                } else { 
8                    vis.SetPix(x + iModel * xDim, y, HeatMapRGB(drug.Get(x, y))); 
9                } 
10            } 
11        } 
12    }

We next look at the DrawModel Function, which is used to display a summary of the model state on a GridWindow object. DrawModel is called once for each model per timestep, see the Main Function section for more information.

Loops over every lattice position of the grid being drawn, xDim and yDim refer to the dimensions of the model.
Uses the built-in GetAgent function to get the Cell that is at the x,y position.
If a cell exists at the requested position, the corresponding pixel on the GridWindow is set to the cel’s phenotype color. to draw the models side by side, the pixel being drawn is displaced to the right by the model index.
If there is no cell to draw, then the pixel color is set based on the drug concentration at the same index, using the built-in heat colormap.

10.8 Imports

1  package Examples._6CompetitiveRelease; 
2  import Framework.GridsAndAgents.AgentGrid2D; 
3  import Framework.GridsAndAgents.PDEGrid2D; 
4  import Framework.Gui.GridWindow; 
5  import Framework.GridsAndAgents.AgentSQ2Dunstackable; 
6  import Framework.Gui.UIGrid; 
7  import Framework.Tools.FileIO; 
8  import Framework.Rand; 
9  import static Examples._6CompetitiveRelease.ExampleModel.*; 
10 import static Framework.Util.*;

The final code snippet looks at the imports that are needed. Any modern Java IDE should generate import statements automatically.

The package statement is always needed and specifies where the file exists in the larger project structure
Imports all of the classes that we will need for the program.
Imports the static fields of the model so that we can use the type names defined there in the Agent class.
Imports the static functions of the Util file, which adds all of the Util functions to the current namespace, so we can natively call them. Statically importing Util is recommended for every project.

10.9 Model Results

Table 2 displays the model visualization at tick 0, tick 400, tick 1100, tick 5500, and tick 10,000. The Figure caption explores the notable trends visible in each image. Figure 3 displays the population sizes as recorded by the FileIO object at the end of every timestep.

Table 2: Selected model visualization PNGs. Blue cells are drug sensitive, Green cells are drug resistant, background heatmap colors show drug concentration. At timestep 0 and timestep 400 (right before drug application starts), all 3 models are identical. At tick 1100 the differences in treatment application show different effects: when no drug is applied, the rapidly dividing sensitive cells quickly fill the domain, when drug is applied constantly, the resistant cells overtake the tumor. Pulsed drug kills some sensitive cells, but leaves enough alive ot prevent growth of the resistant cells. At tick 5500, the resistant cells have begun to emerge from the center of the pulsed drug model. At tick 10000, all domains are filled. Interestingly, the sensitive cells are able to survive in the center of the domain because drug is consumed by cells on the outside. This creates a drug-free zone in which the sensitive cells outcompete the resistant cells.


Figure 3: FileIO population output. This plot summarizes the changes in tumor burden over time for each model. This plot was constructed in R using data accumulated in the program output csv file. Displayed using GGPlot in R

This modeling example illustrates the power of HAL’s approach to model building. Writing relatively little complex code, we setup a 3 model experiment with nontrivial dynamics along with methods to collect data and visualize the models. We now briefly review the model results.

As can be seen in Figure 3 and Table 2, the pulsed therapy is the most effective at preventing tumor growth, however the resistant cells ultimately succeed in breaking out of the tumor center and outcompeting the sensitive cells on the fringes of the tumor. It may be possible to maintain a homeostatic population of sensitive and resistant cells for longer by using a different pulsing schedule or possibly by using adaptive therapy. As the presented model is primarily an example, we do not explore how to improve treatment further. For a more detailed exploration of the potential of adaptive therapy for prolonging competitive release, see [1].


[1]   Jill A Gallaher, Pedro M Enriquez-Navas, Kimberly A Luddy, Robert A Gatenby, and Alexander RA Anderson. Adaptive vs continuous cancer therapy: Exploiting space and trade-offs in drug scheduling.