a.1. Cellular automata modeling.

In brief, a cellular automaton (CA) contains a lattice of any finite number in dimensions of cells. Each CA cell has a state. The number of state possibilities is typically finite. Each CA cell has a neighborhood. This can be defined in any number of ways, but it is typically a list of adjacent cells. The two most common options in a 3D grid of squares are the von Neumann neighborhood, where each cell interacts only with its six horizontal and vertical adjacent mates, and the 26 Moore neighborhood, comprising all the immediately adjacent cells.

The main benefit of using cellular automata in cancer modeling is the ability to observe emerging population level dynamics without a-priori knowledge of tumor behavior [7]. Summary of some important published approaches based on CA is shown in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of some important published CA models.

https://doi.org/10.1371/journal.pone.0183810.t001

a.2. Agent-based modeling.

Agent-based modeling (ABM) is another approach to modeling complex systems composed of interacting, autonomous agents. Agents have behaviors, frequently represented by simple rules, and interactions with other agents, which in turn influence their behaviors. Agent-based modeling suggests a technique to model social systems that are composed of agents who learn from their experiences, and adjust their behaviors so they are better suited to their environment [19]. In an agent-based model, each cell is often represented as an agent. Summary of some important published approaches based on ABM is shown in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of some important published agent-based models (M = Migration, P = Proliferation, Q = Quiescence, Mi = Microscopic, Me = Mesoscopic, Ma = Macroscopic).

https://doi.org/10.1371/journal.pone.0183810.t002

ABM is a more realistic modeling approach for many problems, especially problems in which there are multiple types of actors that interact in different ways [5]. However, it can become very intricate when they incorporate a lot of detail. Therefore, it can become computationally expensive, entailing exceedingly long computer run times for simulations and it asks for ancillary developmental resources. Another disadvantage of ABM is that model developers would have to work on the bottom level of abstraction pursue great endeavor to graphical display, memory management and synchronization mechanism [20].

a.1. Microscopic scale sub-models.

The microscopic scale refers to those phenomena that take place at the subcellular level. Signaling pathway is an important element to accede our proposed models to reality because the fate of a cell is determined by signals received from its environment. For this purpose, we used two important signaling pathways (TNF and EGFR) in our model [45][46]. The concentration of each element in the TNF and EGFR signaling are described by ODEs equations. Time step on this scale is second and the concentration of the materials in every point of the grid is updated.

EGFR signaling pathways: EGFR Signaling pathway starts with binding TGFα to EGFR and affects cell decision to migrate or proliferate. Each agent evaluates the concentration of EGFR at its current location. In our model, the EGFR signaling pathway of [47] is used and its equations are computed using numerical methods. EGFR Equations related to this signaling pathway are given in S1 Appendix.

The input to EGFR equations is TGFα concentration and their output is PLCγ concentration. After solving signaling pathway equations, the value of PLCγ is calculated. In the case that the PLCγ value of a cell is more than the average of PLCγ concentration of the other cells, it can divide, otherwise it can migrate.

TNF signaling pathways: TNF signaling pathway triggers by binding TNFα to TNFR and takes decisions of cell survival and cell death. TNFR is responsible for a diverse range of signaling events within cells, leading to necrosis or apoptosis and inhibiting tumorigenesis. TNF Equations related to this signaling pathway are given in S2 Appendix.

a.2. Mesoscopic scale sub-models.

The mesoscopic scale refers to the cellular level. The term mesoscopic refers to intermediate between short and long and applies to microstructures to be seen in between the atomic and the macroscopic length scales. In this scale, each cell (or agent) is characterized by a site or position (p(x,y,z)) and phenotype or action. The p(x, y, z) is used to express each site (or point) in the lattice, where x, y and z indicate the integer location in Euclidean terms.

We have four agents: healthy cell, cancerous cell, stalk vessel and tip vessel. Transitions between states (or selection a new action) are modeled by an online learning approach (see subsection ‎b). The general idea is to specify such model by means of different states in which it can be, actions which an agent can perform to affect the future behavior, and rewards or costs that depend on the state and decision. After an action has been chosen (based on a learning approach), our model changes its state depending on the action and the current state.

Cell state (agent state): The state of each cell is characterized by number of variables (Table 4). These variables are based on the simulation environment factors such as the value of signaling pathways output/input (Plcγ, TNFα), nutrients concentration (Oxygen, Glucose, VEGF), number of cell neighbors, etc. The value of each variable is continuous; hence we have an infinite set of states. For each of these variables, we define some rules.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Variables of state for each agent (for example the variables of the tip vessel state are VEGF, Doubling Age and Number of Capillary Cells).

https://doi.org/10.1371/journal.pone.0183810.t004

Cell phenotype (agent action): In the cellular scale, we concentrate on cell behaviors. Each cell (at time t) is assumed to have one phenotype (presented below) according to its type (Cancerous, Healthy, Tip, Stalk). Table 5 shows the phenotypes of each agent (See [48] for a discussion of the biological theory of this modeling construct).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Cells can have the following actions based on their types in each site (for example the Tip vessels undergo the actions of branch and expansion).

https://doi.org/10.1371/journal.pone.0183810.t005

• Apoptosis

Apoptotic healthy cells undergo programmed cell death in response to signaling events [48]. If healthy cells gain apoptotic phenotype, they become inactive and remove from simulation environment.

• Necrosis

Necrosis is a form of cancerous cell injury which results in the premature death of cells in living tissue by autolysis. On the contrary, apoptosis is a naturally occurring programmed and targeted cause of cellular death. The only effect of this phenotype will be a grid point occupation, and therefore other cells cannot enter at this point of the grid. By occurring necrosis, necrotic cell not only remains in simulation environment but also preserves from any changes.

• Migration

Cell always searches for a place with more nutrition to migrate or to delivers its offspring there. The candidate position (p(x,y,z)) is selected according to the following condition (Eq 1) [14]. (1) Where:

1. ∨ is OR operator.
2. = maximum Concentration of nutrient among all Moore neighbors of candidate position (p(x,y,z))
3. = Average nutrient concentration for candidate position among all its Moore neighbors
4. = Standard Deviation of nutrient concentration for candidate position among all its Moore neighbors

If the candidate site doesn’t have the above condition or two or more candidates have the above condition, all candidate locations will be ranked through Eq 2 and according to their ProbabilityOfSelection value, one candidate will be selected randomly. (2) Where:

1. 
2. = Average nutrient concentration for candidate position among all its Moore neighbors
3. 

Each lattice site can be occupied by one cancerous/healthy cell or four endothelial (vessel) cells. The cancer cells compete with healthy cells for the empty sites and the cancerous cell outperforms its normal counterpart. Therefore, if a cancerous cell and healthy one simultaneously select an empty site, the cancerous cell will gain the site. If the type of competitors is the same as each other, one of them will be selected randomly. If no movement at all is possible, the cell will remain stationary.

• Proliferation

Cell proliferation is the process that results in an increase of the number of cells. Cell proliferation is increased in cancer. Cancerous or healthy cells can enter the proliferation phenotype. In this phenotype, cell divides into two identical offspring cells. In our model, one of the offspring cells is located in their mother grid site and the other one is located in one of the free grid sites in a Moore neighborhood of its mother based on Eqs 1 and 2. However, as soon as its neighboring spaces are occupied by other cells, the cell moves to resting phase (quiescence phenotype) [11,49].

• Quiescence

The quiescent phenotype is the default phenotype for a cell. It represents G0 in terms of the cancerous or healthy cells cycle. For endothelial cells, this phenotype means they don’t do anything. It needs to be mentioned here that there is no reverse change of cell from quiescent back to other phenotypes in this model.

• Hypoxia

Hypoxia is a phenotype for a cancerous cell in which the cancerous cell is deprived of adequate oxygen supply. Decrease oxygen availability (hypoxia) stimulates cancerous cell to produce VEGF much more. Cancerous cells enter hypoxia when oxygen levels drop below a defined threshold and will enter necrosis if they remain at that level for too long.

• Branch

The endothelial cells spearheading the vascular sprouts are known as the endothelial tip cells. Tip cells take the decisions of vessel branch, thereby defining the route in which the new sprout grows [30].

In our model, tip cells deliver two new offspring in two free sites of the grid in its Moore neighborhood of its mother based on Eqs 1 and 2.

• Expansion

Tip cells are the leading cells of the sprouts; they guide following endothelial cells (stalk cells) and sense their environment for guidance cues. As the tip cells move out from the existing blood vessel, the stalk cells follow to sprout. As tip cell moves forward, a new stalk cell is created in back of it.

• Sprout

Following the tip cells are the endothelial stalk cells, which are highly proliferative, establish adherent and tight junctions to ensure the stability of the new sprout. Stalk cell always searches for a place with more VEGF to deliver its daughter cell there [50]. The candidate place is selected according to Eqs 1 and 2. Type of new vessel generated from the existing stalk cells is the same as its parent. Arteriole and Venule are separated from artery and vein respectively.

Material diffusion: We considered five materials in the environment (Oxygen, Glucose, TGFα, VEGF and TNFα). The concentration (C) of each material is identified within each of environment grid points. Environment diffuses the material by PDE diffusion equation based on its source (S), uptake (U) and waste (W). In this model, duty of endothelial cells are the production of Oxygen, Glucose (which is consumed by cancerous/healthy cells) and consumption (or uptake) of VEGF. The duty of cancerous/healthy cells are the production of VEGF (which is consumed by endothelial cells), TGFα and TNFα (which is consumed at mesoscopic scale based on equations in S1 Appendix and S2 Appendix).

• Oxygen, glucose TGFα and TNFα diffusion

The concentration of each material is identified within each of environment grid points. Environment diffuses the material by the diffusion equation. We rewrite Diffusion equation for discrete environments in Eq 3 for oxygen, glucose, TGFα and TNFα based on [18] and [51]. For TGFα and TNFα, uptake is computed based on their ODEs equations (see S1 Appendix and S2 Appendix). Grid sites occupied by healthy and cancer cells are sinks of oxygen, glucose, TGFα and TNFα. (3) Where:

1. D_nutrient is nutrient diffusion coefficient,
2. 
3. 
4. 
5. 
6. 
7. r_art and r_cap are radius of arteriole and radius of capillary respectively
8. 
9. 
10. nutrient = {oxygen,glucose,TGFα,TNFα}
11. β_c is constant variable (Table A in S3 Appendix)
12. ΔS = the size of grid point

• Vascular endothelial growth factor diffusion

Vascular endothelial growth factor (VEGF) stimulates the ingrowth of a new blood supply from the host vasculature via angiogenesis. Cancerous cells release growth factors (VEGF), in the surrounding normal tissue. We rewrite Diffusion equation for discrete environments in Eq 4 for VEGF. (4) Where:

1. D_VEGF is VEGF diffusion coefficient,
2. 
3. 
4. 
5. r_art and r_cap are radius of arteriole and radius of capillary respectively
6. 
7. 
8. 
9. ω_VEGF and ϕ_c are constant variables (Table A in S3 Appendix)
10. ΔS = the size of grid point

a.3. Macroscopic scale sub-models.

The macroscopic scale pertains to processes happening at the tissue level. At the tissue level, we consider vessels and tumor directly.

Angiogenesis: Angiogenesis is termed as the growth of new blood capillaries on the basis of pre-existing vessels, which is a critical step in cancer growth and metastasis. Tumorigenesis is a prime case of an emergent tissue patterning phenomenon reliant on many cell types, both normal and aberrant cell behaviors. Models of angiogenesis fall into three classes: models of vasculature growth without tumor or other tissue cells, models of vasculature growth with only tumor cells, and models of vasculature growth with both tumor and healthy tissue cells [31]. Vasculatures in the system give cells (agents) materials (oxygen, glucose, TGFα and TNFα), which they require to survive. Once vascularized, cancer has access to a huge nutrient source and rapid growth ensues. At the other hand, VEGF stimulates the ingrowth of a new blood supply from the host vasculature via angiogenesis. In this section, we present a model of angiogenesis and vascular tumor growth. We assume that there are eight main parallel parent vessels (four veins and four arteries located on the four vertical edges and four vertical surfaces of cube respectively) as the major part of the blood supply in our 3D simulation environment as shown in Fig 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Four veins and four arteries located on the middle of four surfaces and four edges of cube.

https://doi.org/10.1371/journal.pone.0183810.g002

In our model, the general order of blood flow is as follow: Arteries → arterioles → venules → veins. When we refer to an “active” vessel, we are referring to an arteriole in which it can connect to a venule. As arteriole connects to venule, the blood flow circulated throughout it.

b.1. Regression sub-model (SVR-NSGA-II).

It is well known that SVR generalization performance (estimation accuracy) depends on a good setting of C, Ԑ and kernel function parameter (e.g. gamma). Existing software implementations of SVM regression usually treat these parameters as user-defined inputs. Selecting a particular kernel type and kernel function parameters are usually based on application-domain knowledge and also should reflect the distribution of input (x) values of the training data.

Parameter C decides the tradeoff between the model complexity (flatness) and the degree to which deviations larger than Ԑ are tolerated in optimization formulation. Parameter Ԑ controls the width of the Ԑ-insensitive zone, used to fit the training data. The value of Ԑ can affect the number of support vectors used to construct the regression function. Parameter gamma (γ) is the free parameter of some kernels (e.g. RBF, polynomial and sigmoid kernels).

Kernel function maps the data into higher dimensional spaces in the hope that in this higher-dimensional space, the data could become more easily separated. Each kernel function can extract a specific type of information from a given dataset. Now the question arises how to get an appropriate kernel function. Since each kernel has some degree of variability in practice, there is nothing else for it but to experiment with different kernels and regulate their parameters.

In a previous study, we proposed a method (SVR-NSGA-II) [52] that is an improved algorithm based on SVM. SVR-NSGA-II optimizes the above parameters (C, Ԑ and γ) by implementing the evolutionary process (NSGA-II). The Non-dominated Sorting Genetic Algorithm II (NSGA-II) is a Multiple Objective Optimization (MOO) algorithm and is an instance of an Evolutionary Algorithm from the field of Evolutionary Computation.

For single objective optimization problems, the best single design is the goal. But for multi-objective problems, with several and possibly conflicting objectives, there is usually no single optimal solution. A suitable solution should provide for acceptable performance over all objectives.

We convert the RL problem into regression problem, which the observed states and actions are viewed as input variables and the Q values as output variables. To collect enough training samples, simulation iterations must be high enough for accurate results. In order to improve learning performance, the training samples are slide with the manner of time-window and the newest samples are considered as train dataset.

Assume that we are given a training set D, consisting of pairs (x_i,y_i), for i = 1,…,m. Each sample x_i is a vector that describes the cell states, phenotypes (action) and Q values. In other words, the phenotypes, the variables of cell states, constitute the features of the dataset. The label y_i associated with x_i is a discrete value between -1 and 1 and describe the Q value.

An initial population of chromosomes is generated randomly. Then the population is evolved for a number of generations. The chromosome is a fixed length list of genes or parameters.

We have implemented three kernels Sigmoid, Polynomial and RBF for our method. SVR-NSGA-II aggregate these three kernels which are performed through a simple voting based on allocating weight to each of them. This weight is based on F-Measure of its kernel for train dataset (Eq 5).

F-measure is a harmonic mean between precision (or Confidence) and recall (or Sensitivity). Recall is the percentage of real positive cases that are correctly predicted positive. Precision indicates the percentage of predicted positive cases that are correctly real positives. Also, NSGA-II optimizes Recall, F-Measure and precision in a multi-objective way. (5) where:

1. 
2. TP is the number of correct predictions in which an instance is positive;
3. FN is the number of incorrect predictions in which an instance is negative
4. FP is the number of incorrect predictions in which an instance is positive
5. TN is the number of correct predictions in which an instance is negative
6. k = {Sigmoid,RBF,Polynomial}

b.2. Q-learning sub-model.

Q-learning is a simple incremental algorithm developed from the theory of dynamic programming for delayed reinforcement learning. The objective of Q-learning is to quantify the Q value for an optimal value when the system dynamics features are unexplored. Q-learning approach is fulfilled as follows: a Q-learning system observes the current state s_t and executes an action a_t at each time step t. Next, perceive the subsequent state s_t+1 and receives an immediate reward r_t. The value of Q can be adapted based on Eq 6 where η is a learning rate that controls the learning speed and γ is a discount factor used to determine the proportion of delay to the future rewards.

(6)

Policy (cellular phenotype decision): To enable healthy/cancerous cells to evolve in the 3D domain, we need to define some rules. These rules clarify the relationship between cell states (see Table 4) and cell actions (see Table 5). For each rule, if the probability of selection for some phenotypes increases, it decreases for other phenotypes. In Table 6, policies are determined on a scientific basis using evidences derived from various studies.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. The relationship between state variables and actions according to policies.

https://doi.org/10.1371/journal.pone.0183810.t006

Reward: The objective of the learner is to choose actions maximizing discounted cumulative rewards over time. The value of reward at time t (r_t) is calculated based on the following reward function (RF). In Eq 7, for each policy (p) the value of the reward function is equal to 0 when the value of x is equal to the threshold and also its sign changes at this point. It is obvious that the value of RF is between -1 and +1. A gradient is a constant number that describes both the direction and the steepness of the line. RF is increasing/decreasing if it goes up/down from left to right if the gradient is positive/negative respectively (see Fig 4). (7) where:

1. x equal to value of state variables

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Four reward functions examples.

Threshold point (or inflection point) is a point on a curve at which the curve changes from being concave (concave downward) to convex (concave upward), or vice versa. Gradient (or slope) specifies the rate of RF change.

https://doi.org/10.1371/journal.pone.0183810.g004

In Table B in S3 Appendix, we specify a threshold for all rules that the sign of its probability changes at this point. We use the average of above function for all policies to calculate the total value of the reward.

Train phase.

1. Step 1. Check the current state s_t
2. Step 2. Construct test sample (s_t,a_k) comprised of each action a_k in action set (A) and current state s_t
3. Step 3. Constitute dataset for training (select the d newest samples)
4. Step 4. For all available k actions, predict Q_k corresponding to (s_t,a_k) by SVR-NSGA-II
5. Step 5. Randomly select action a_t based on its Q_t value ()
6. Step 6. Run proposed hybrid multiscale model (run diffusion and signaling pathways) and perform action a_t, obtain r_t and the successor state s_t+1
7. Step 7. Update Q_t according to Eq 6
8. Step 7. Add the new observed data into dataset (s_t, a_t, Q_t) if r_t>θ (θ is a predefined threshold value)

Test phase.

1. Step 1. Check the current state s_t
2. Step 2. Construct test sample (s_t,a_k) comprised of each action a_k in action set (A) and current state s_t
3. Step 3. Constitute dataset for training
4. Step 4. For all available k actions, predict Q_k corresponding to (s_t, a_k) by SVR-NSGA-II
5. Step 5. Select action with highest Q_t value
6. Step 6. Run proposed hybrid multiscale model (run diffusion and signaling pathways) and perform action