Cellular Potts Model: Applications to Vasculogenesis and Angiogenesis

The cellular Potts model (CPM, a.k.a. Glazier–Graner–Hogeweg or GGH model) is a somewhat liberal extension of probabilistic cellular automata. The model is derived from the Ising and Potts models and represents biological cells as domains of CA-sites of the same state. A Hamiltonian energy is used to describe the balance of forces that the biological cells apply onto one another and their local environment. A Metropolis algorithm iteratively copies the state from one site into one of the adjacent sites, thus shifting the domain interfaces and moving the biological cells along the lattice. The approach is commonly used in applications of developmental biology, where the CPM often interacts with systems of ordinary-differential equations that model the intracellular chemical kinetics and partial-differential equations that model the extracellular chemical signal dynamics to constitute a hybrid and multiscale description of the biological system. In this chapter we will introduce the cellular Potts model and discuss its use in developmental biology, focusing on the development of blood vessels, a process called vascular morphogenesis.Wewill start by introducing a range of models focusing on uncovering the basic mechanisms of vascular morphogenesis: network formation and sprouting and then show how these models are extended with models of intracellular regulation and with interactions with the extracellular micro-environment. We then briefly review the integration of models of vascular morphogenesis in several examples of organ development in S.E.M. Boas · R.M.H. Merks (B) · E.G. Rens Centrum Wiskunde & Informatica, Science Park 123, 1098 XG Amsterdam, The Netherlands e-mail: merks@cwi.nl S.E.M. Boas · R.M.H. Merks · E.G. Rens Leiden University, Mathematical Institute, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands S.A. Prokopiou Department of Integrated Mathematical Oncology, H. Lee Moffitt Cancer Center and Research Institute, 12902 Magnolia Drive, Tampa, FL 33612, USA Y. Jiang Department of Mathematics and Statistics, Georgia State University, Atlanta, GA 30303, USA e-mail: yjiang12@gsu.edu © Springer International Publishing AG 2018 P.-Y. Louis and F.R. Nardi (eds.), Probabilistic Cellular Automata, Emergence, Complexity and Computation 27, https://doi.org/10.1007/978-3-319-65558-1_18 279 280 S.E.M. Boas et al. health and disease, including development, cancer, and age-related macular degeneration. We end by discussing the computational efficiency of the CPM and the available strategies for the validation of CPM-based simulation models.

health and disease, including development, cancer, and age-related macular degeneration. We end by discussing the computational efficiency of the CPM and the available strategies for the validation of CPM-based simulation models.

Introduction
Probabilistic cellular automata (PCA) are widely applied as a modeling framework for biological and biomedical research. In particular, they are used in the study of biological pattern formation, to help to understand how biological structures can form from biological elements that follow simple rules. In this way, PCAs have been used in diverse applications, ranging from spatial structuring of ecosystems [37] to a range of biomedical problems, including the self-organization of the autonomous nervous system in the gut (Chap. 17 of this book), the growth and plasticity of tumors [3,27] or the formation of blood vessels [17]. The central question in these applications of PCA is how cells can self-organize into tissues and organs. The states of the PCA represent the types of biological entities, while the probabilistic nature of the PCA reflects the "noisiness" inherent to most physical and biological systems. The PCA here helps unravel how biological patterns can persist in the presence of homogenizing noise, or, perhaps more interestingly, PCAs demonstrated that noise can become a driving force of pattern formation, i.e., no patterns would form in a deterministic model [50,81]. The "PCAs" in the above applications typically deviate from the strict definition of PCA, as a system of locally coupled, homogenous system of Markov chains with synchronous updates. The updates can be asynchronous (one by one in random order), additional rules are applied (e.g., rules for mass-conserved random walks or diffusion), or the systems are hybridized with systems of partial-differential equations, e..g., to model diffusing molecular signals.

Cellular Potts Model
A generalization of PCA, which is particularly widely applied to biomedical problems, is the cellular Potts model (CPM), also known as the Glazier-Graner-Hogeweg or GGH model [33]. The CPM is used to model structures of biological cells and extracellular materials. It is a generalization of the large q-Potts model, which derives from the Potts model. Glazier et al. [32] reviewed the derivation of the cellular Potts model from its predecessors in detail; a brief recap is useful in the present context to better understand the structure and notation of the CPM. The Potts model studies the interactions between domains on a lattice, e.g., during the solidification of a fluid. It is defined on a regular lattice Λ ⊂ Z 2 or Λ ⊂ Z 3 , with x ∈ Λ the coordinates in the lattice. The clusters of like spins, σ(x) ∈ {0, . . . , q}, represent individual domains, where the same spin can identify multiple domains if they are well separated spatially. Assuming that (without external fields) the spins follow Boltzmann statistics, the relative probability of each configuration of spins {σ(x)} is, , (18.1) where the Hamiltonian, H ({σ(x)}), describes the energy of the configuration, k is the Boltzmann constant, and T is the absolute temperature. In the Potts model, the interfaces between two domains contribute to the configuration energy: (1 − δ(σ(x), σ(x ′ )). (18.2) Here J (typically J ≥ 0, but see also Ref. [61]) is the energy associated with a unit length of the domain interfaces, (x, x ′ ) is a pair of adjacent lattice sites, and the Kronecker delta function (δ(a, b) = 1, if a = b; 0, if a ̸ = b, a ∧ b ∈ Z) selects adjacent lattice sites of unequal spin. Minimizing the Hamitonian energy using Monte Carlo methods tends to minimize the number of interfaces in the system by forming domains of identical spin, which would coarsen to form fewer and fewer domains.
The key innovation of the CPM was to associate each of the domains formed in the Potts model with a biological cell. To ensure that the volume of the cells (or area in 2D) is approximately conserved, the CPM adds a volume energy to the Hamiltonian, H volume = λ volume (v(σ) − V (σ)) 2 , where λ volume is a Lagrange multiplier to the volume constraint; similar terms with Lagrange multipliers can be added to represent additional optimization conditions. The volume energy of a domain with spin s is zero if its actual volume, v(s) = | {x|x ∈ Λ ∧ σ(x) = s} | (i.e., the number of sites in the lattice of spin equal to s), is equal to a target, or resting volume V (s). Any deviation of the actual volume to the target volume contributes elastically to the total energy. For the medium that surrounds the cells (locations x with σ(x) = 0), no volume constraint is applied. Note that the volume constraint adds a non-local, cellular scale to the system. The energy of a configuration in the CPM thus depends both on local, nearest-neighbor interactions, as well as on non-local properties of all sites in the lattice of equal spin value.
A further innovation of the CPM is a differentiation between interfacial energies, such that one type of interface may be favored over another. Each cell, σ(x), also has a type, τ (σ(x)) ∈ N, with each value of τ classifying the domain as a particular biological cell type (e.g., neuron, muscle cell), or a cell state (e.g., proliferating vs quiescent), or non-cellular material (e.g., fluid, substrate, and so forth). The interfacial energy, J , then becomes a function of the pair of types at the interface, (τ (σ(x)), τ (σ(x ′ )).
The full Hamiltonian of the CPM becomes, 3) The term H ′ identifies any additional constraints on the cell behavior, particularly those involving interactions with external fields [73,91] or additional constraints at cell level (e.g., constraints on the length of the cellular interfaces [40], or constraints on cell shape [56,99]).
Minimization of the Hamiltonian energy function corresponds to solving the balance of forces applied to the cell. The Hamiltonian in the CPM is usually minimized using Metropolis dynamics, which compares configurations differing by one spin at a time. In the CPM, the Metropolis algorithm is modified such that it mimics natural fluctuations of the membranes driven by the activity of the cytoskeleton. This feature introduces a temporal ordering into the energy minimization procedure, such that both the equilibrium configuration and the transition to equilibrium are of physical and biological interest.
More specifically, the CPM-version of the Metropolis algorithm selects a pair of adjacent lattice sites, (x, x ′ ), at random from the lattice; i.e., first a target site x ∈ Λ is selected at random, then a lattice site x ′ is selected at random from NB(x), the set of neighbors of x. On a square lattice, typical choices for NB(x) include the Moore neighborhood (the eight nearest neighbors); larger neighborhoods, e.g., the twenty neighbors of order 1-3 are also used to reduce lattice effects [52]. Next the algorithm attempts to change the spin of x (the target site), into the spin of its neighbor σ(x ′ ) (the source site). If the attempted update will reduce the energy, i.e., ∆H = H after − H before < 0, the change occurs with probability 1. If the attempt increases the energy (∆H > 0), the change will be accepted with Boltzmann probability: In contrast to the Potts model, in the CPM the temperature T is a cellular temperature, reflecting the amplitude of active cell membrane fluctuations. For lack of measurements of the distribution of the energy that is mechanically dissipated during these fluctuations, the CPM follows the Potts model by assuming Boltzmann probability. Time in the cellular Potts model is measured in Monte Carlo Steps (MCS), where one MCS corresponds with |Λ| copy attempts, i.e., the number of sites in the lattice. To identify the real time corresponding with one MCS, the kinetics of the CPM itself [33,95] or the kinetics of coupled models [59,92] is matched with the kinetics of experiments. Similar approaches are used to parameterize the real volume or area corresponding to one lattice site.

Generic Behavior of the CPM
The cellular Potts model (CPM) is widely applied to biomedical problems involving cell shape changes and cell-cell adhesion. It was introduced [33,34] in the early 1990s as a model for differential-adhesion-driven cell rearrangement: a proposed mechanism for spontaneous rearrangement of cell types in mixtures of embryonic cells [38]. The differential adhesion hypothesis suggests that these cellular rearrangements, also known as cell sorting, are driven by relative adhesion surface energies of different cell types [84]. This hypothesis has been tested with the CPM, using Typical time course of a cellular Potts simulation of binary cell sorting through differential adhesion. Cell type r (red) engulfs cell type b (blue) due to differences in adhesion energies; in this example, J dd < J ld < J ll < {J lM , J dM } a 0 MCS; b 100 MCS; c 500 MCS; d 1000 MCS; e 5000 MCS; f 10000 MCS only a volume constraint and adhesion energies [33,34]. The model was initialized with two cell types, light (l) and pigmented, dark (d) colored cells, mixed in a random aggregate surrounded by medium (M). If the adhesion energies are set as J dd < J ld < J ll < {J lM , J dM }, the two cell types would segregate into clusters of l and d, and the light aggregate would eventually engulf dark ones d ( Fig. 18.1). By changing the relative adhesion energies, other patterns can be generated, including checkerboard-patterns. The differential adhesion hypothesis has been tested experimentally by Krieg et al. [45]. They measured cell adhesion at the single cell level by using an atomic force microscopy and set the relative values of J in the CPM accordingly. As a result, the CPM sorted the cells the wrong side out: consistently the so called ectodermal cells ended up surrounded by mesodermal cells, whereas based on the measured cell adhesion values the CPM predicted the mesodermal cells should end up in the middle. Based on experimental observations they predicted that interfacial tensions generated by muscle-like actomyosin structures near the cell surface must also be considered to correctly predict the outcome of cell sorting. Krieg et al. [45] have incorporated the additional source of interfacial energy in the values of J . An alternative approach for modeling cortical tensions is to constrain the perimeter or surface area of the cells.

Hybrid Modeling
Although explicit modeling of shapes and adhesion is essential in many research problems dealing with cells, other components in biomedical systems benefit from a continuum modeling approach. For example, the diffusion of a chemical is typically modeled with a partial-differential equation (PDE). Many research problems lie at the interface of these two modeling approaches, in which case we can use hybrid modeling. Hybrid models combine multiple types of modeling techniques, such as discrete and continuum modeling.
A widely used approach is to couple a field or a set of fields representing, e.g., the distribution of a chemical signal, to the CPM. The CPM is then modified such that cells respond to the chemical field by moving to higher or lower concentrations, Chemotaxis toward an chemoattractant secreted by the cells themselves results in clustering of initially dispersed cells. Isolines (green lines) indicate ten chemoattractant levels relative to the maximum concentration in the simulation a mechanism called chemotaxis. The typical way to model chemotaxis is to modify ∆H during an attempted update, such that moves up the chemical gradient occur with higher (or lower) probability proportional to the gradient [73], with c the chemical field. The resulting bias lets cells gradually move up (or down) the gradient of the chemoattractant with a sensitivity λ chem . Figure 18.2 shows an example simulation of a hybrid CPM proposed by Merks et al. [56,59]), where the cells (colored in red) in medium (white background) aggregate because of a chemical signal. In this model, the cells secrete a chemical to attract surrounding cells. This chemoattractant diffuses and slowly decays in the medium, following the PDE, with α the secretion rate by cells, D the diffusion constant of VEGF and ϵ the decay of VEGF in the medium. After each MCS, the PDE (Eq. (18.6)) is solved numerically using a discretization matching the grid of the CPM, allowing the CPM to read out the chemoatractant concentration on each lattice site. The mutual attraction results in cell aggregation.  [57]; it is available at http://sourceforge.net/projects/tst. TST is a C++ library providing implementations of the CPM on two-dimensional, square lattices and functionality for hybrid CPM-PDE models. While limited in functionality relative to CompuCell3D, the simplicity of the TST allows more straightforward access to the CPM implementation. This makes it particularly well suited as a test bed for new CPM algorithms and extensions of the Hamiltonian. A recent tutorial [24] provides detailed instructions for how to adopt TST for one's own needs. Most simulations reported in Sects. 18.4 and 18.5 used the TST. The hybrid CPM and finite-element model reviewed in Sect. 18.4.4 was implemented in an independent C-code, released as supplements to its publication [92]; it may soon be merged with the TST.

Implementations of the CPM
Other implementations of the CPM are part of the free simulations environments Morpheus [83] and Chaste [67]; both these environments provide a range of biological modeling formalisms, including cellular automata, ordinary-differential equations and partial-differential equations, and off-lattice cell-based modeling techniques. Morpheus provides a high-level, XML-based declarative programming language to describe model rules and has an attractive user interface. The first releases of Morpheus were closed-source, limiting its applicability, but an open source release has been announced as of this writing (April 2015). Morpheus is available from http://imc.zih.tu-dresden.de/wiki/morpheus/doku.php. Chaste is a large C++ software library focusing on cardiac electrophysiology and cell-based modeling. The latter component contains cellular Potts functionality. Chaste is available from http:// www.cs.ox.ac.uk/chaste/.

Application of the Hybrid CPM to Blood Vessel Formation
Over the past decades, a number of mathematical and computational models have been developed to propose new models for the mechanisms of embryonic development. These mechanisms span all spatial and temporal scales encompassed by this complex process, ranging from the molecular level all the way to the organismal level and its environment, and ranging from microseconds (chemical reactions) up to years (homeostasis, aging, cancer). Although at present "computing" a human is beyond reach (although successful first steps have been taken for much simpler multicellular creatures [77]), cellular Potts modeling has been applied to the understanding of relatively more simple mechanisms, including the formation of blood vessels. During embryogenesis, vascular networks (blood vessels) are formed from initially dispersed endothelial cells (ECs), a process called vasculogenesis. Once the vasculature is established, capillary sprouts can branch off from the preexisting vasculature in response to externally supplied angiogenic stimuli, a process called angiogenesis. The new sprouts provide tissues and organs with oxygen and nutrients, and remove metabolic waste. Angiogenesis takes place in physiological situations, such as embryonic development, wound healing and reproduction [18]. The healthy body controls angiogenesis by balancing pro-and anti-angiogenic factors [19]. This balance, though, is sometimes disrupted and angiogenesis also appears in many pathologies, like diabetes [55], rheumatoid arthritis [43], cardiovascular ischemic complications [16], proliferative retinopathy [30], and cancer [29].
Sprouting angiogenesis typically starts from hypoxic tissues or cells (e.g., retinal astrocytes [76]) upregulating their production of pro-angiogenic factors such as vascular endothelial growth factor A (VEGFA) [28]. These angiogenic factors diffuse and bind to endothelial cell receptors on nearby blood vessels. Subsequently, the extracellular matrix (ECM) and basement membrane, surrounding the ECs, are degraded locally by activated proteases (e.g., matrix metalloproteinases, MMPs) produced by ECs.
Mathematical modeling is a useful tool for understanding the mechanisms of angiogenesis and to design experiments of a predictive nature. Since vessels often consist of only a few cells, explicitly considering individual cells is essential. In most modeling frameworks, the detailed investigation of cell-level properties, such as cell shape and cell adhesion, are mathematically difficult, if not impossible to consider. Therefore, the CPM with the advantage of representing cells as individual entities with a particular shape is an appropriate framework to study blood vessel formation. Over the past two decades, a plethora of mathematical and computational models have been developed to study aspects of angiogenesis. For a comprehensive review of mathematical and computational models in angiogenesis see [65], and references therein.
Typical modeling studies investigate how growth factors and receptors promote endothelial cell proliferation, how groups of endothelial cells assemble into individual vessels, and how tumors recruit the ingrowth of whole microvascular networks.
Here we briefly review cellular Potts approaches to the analysis of blood vessel formation, describing the required extensions to the CPM in technical detail.

Modeling Collective Cell Behavior During de Novo Vasculogenesis
Vasculogenesis is an embryogenic process during which endothelial cells organize into vascular networks. Computational modeling has been used in the search for generic cell behaviors that drive vasculogenesis. In this section, we will discuss four CPM-based models in which the individual behavior of initially dispersed cells collectively results in their organization into vascular networks.

Chemotactic Cell Aggregation
Serini et al. [78] showed in a continuum model that dispersed cells self-organize into polygonal patterns, when these cells secrete a chemoattractant to which all cells respond. However, Merks et al. [59] showed in a CPM-based model that with these assumptions on a longer time scale dispersed cells form rounded aggregates (see 18.2.2 and Fig. 18.2) rather than polygonal patterns. Merks et al. [59] suggested additional model assumptions to explain vascular network formation: (1) endothelial cells (ECs) adhere to one another with VE-cadherins and (2) VE-cadherin-binding inhibits VEGF signaling by interacting with the VEGF receptor 2, and (3) the ECs secrete a chemoattractant (e.g., VEGF [22]) that attracts other ECs. These assumption were simplified in the model as follows: the ECs only responded to the chemoattractant at regions of their membrane adjacent to the ECM, whereas at cell-cell interfaces the chemotaxis was inhibited, a mechanism called contact-inhibited chemotaxis. Indeed,  [59], b cell elongation with chemotaxis [56], c preferential adhesion to elongated cells [89], and d mechanical cell-matrix interactions [92]. Panel C reprinted from Ref. [89]; copyright (2008), with permission from the Biophysical Society the simulations showed self-organization of endothelial cells into vascular networks ( Fig. 18.3a). VE-cadherin mediated contact-inhibited chemotaxis was implemented in the cellular Potts model by only letting cells chemotact at cell-medium interfaces of their membrane, and not at cell-cell interfaces, resulting in λ chem = 0 at cell-cell interfaces and λ chem > 0 for cell-medium interfaces using chemotaxis as described in Sect. 18.2.2. The complete Hamiltonian of the model depends on a volume constraint for the cells, adhesion between the cells, and contact-inhibited chemotaxis. The concentration field of VEGF is described in Sect. 18.2.2.
Although VE-cadherin mediated contact-inhibited chemotaxis reproduces vasculogenesis, Köhn-Luque et al. [44] note that the diffusion speed assumed for VEGF by Merks et al. [59] is much lower than reported for most VEGF isoforms. Köhn-Luque et al. [44] propose an alternative CPM-based model for vascularization in which VEGF, containing ECM-binding domains, is secreted by the underlying endoderm. Endothelial cells scavenge VEGF by the secretion of ECM and subsequently chemotact more strongly to ECM-bound VEGF than to soluble VEGF, resulting in network formation.

Cell Elongation
Endothelial cells are often seen to elongate during network formation. Palm and Merks [63] showed with a CPM-based model that elongated, adhesive cells can self-organize into vascular structures. Cells aggregate into elongated structures that can only rotate very slowly, while connected in the branch points. If the model would run for infinity, the cells would form a spheroid, but this process is so slow that the cells dynamically arrest in a network-like pattern. Addition of chemotaxis to an auto-secreted chemoattractant to this cell elongation model [56] stabilizes network formation and speeds up the patterning process ( Fig. 18.3b). Cell elongation is modeled by the addition of a length constraint to the Hamiltonian: with L describing the target cell length, l the actual cell length and λ length the Lagrange multiplier. To preserve the integrity of the cells, a penalty is added to the Hamiltonian when a copy attempt would break up the cell. The length of the cell, l(σ), is usually estimated by taking the length of the long axis of an ellipse fitted to the cell by calculating the inertia tensor of the pixels belonging to the cell [56,99,100]. In cellular Potts implementations, the inertia tensors of the cells can be efficiently calculated by keeping track of the first and second order raw momenta of the cellular coordinates, in addition to the cellular area or volume [56]. The length constraint requires an additional connectivity constraint to prevent cells from splitting up into two disconnected patches; the connectivity constraint prevents updates that would split a cell in two patches [56]. Chemotaxis is implemented as described in Sect. 18.2.2.

Preferential Attraction to Elongated Structures
Similarly to the work of Palm and Merks [63], a mechanism for vascular patterning excluding external factors has been proposed by Szabó et al. [88,89]. The reason for a mechanism not involving chemical and mechanical forces originates from experiments showing vascular patterning under normal tissue conditions on a solid substrate [88]. Based on experimental results of [89] a new hypothesis for vascular network formation was proposed. These experiments showed elevated cell motility within the presence of elongated structures and cells were observed to migrate faster within narrow sprouts, while cells in wider sprouts have a decreased motility. Furthermore, the width of more elongated sprouts increases with a faster rate. The led to authors to propose that cells are highly attracted to elongated structures. Szabó et al. [88] implemented this attraction to elongated structures in a particle based method in previous work [88], where cell shape was not resolved.
In their follow-up work, Szabó and coworkers [89] added their mechanism to the CPM by adding a bias at the time of copying (cf. Eq. (18.5)), as the two eigenvalues of the cellular inertia tensor, representing the long and short axis of the cell. Thus, θ(σ) is a measure of the eccentricity of a cell with spin equal to σ. The summation in Eq. (18.8) only goes over the neighboring sites of x ′ that belong to cells other than σ(x), σ(x ′ ). The term 1 − δ(x, 0)) − (1 − δ(x ′ , 0)) ensures that the medium is not influenced by this preference, and that copies at cell-cell interfaces are independent of this preference, such that no cell has more advantage than an other cell. If Eq. (18.8) is rewritten, see [89], it becomes clear that it can be considered as an asymmetric extension of the adhesive energy J . Thus, the attraction to elongated structures is actually a preferential adhesion to elongated structures.
Simulations of this model with initially dispersed patches of connected cells show network formation (Fig. 18.3c). Subsequent to an initial budding of one cell from a connected patch, other cells are attracted to the base of the sprout and follow the leading elongated cell. Migration of the sprout continues until a branch is established and stabilized. Due to surface tension, branches can break up, whereas new branches form continuously. The resulting networks are thus quasi-stable, the networks change continuously, while overall the statistical properties (number of branches, wavelength, and so forth) of the pattern are stable. A biological basis for the preferential adhesion to elongated structures is not yet established. Szabó et al. [89] propose that the preference can arise from mechanical tension of elongated structures, which cells can respond to by VE-cadherin based mechanosensing.

Mechanical Cell-ECM Interactions
The models described above explained vascular network formation based on chemical interactions between cells. The models by Palm and Merks [63] and Szabó et al. [89] suggested that cells are able to form network-like structures in absence of a substrate to transmit the chemical signal by elongating or by preferentially adhering to elongated structures, respectively. Another explanation for network formation can be found when considering the mechanical environment of cells. The extracellular matrix (ECM), a network of extracellular proteins that surround most cells in tissues, dictates the mechanical environment. The rigidity of the ECM influences cell behavior; cells have been observed to migrate in the direction of higher stiffness [48], and orient to the direction of stretch [36]. Further, focal adhesions, macromolecular assemblies by which the cytoskeleton connects to the ECM, stabilize under mechanical force [71] or on rigid substrates [66]. Cells do not only respond to the mechanical properties of the ECM, but also actively deform it [35,97]. By applying traction forces, induced by stress fibers within the cell, cells can locally orient [93] and stiffen [97] the substrate they adhere to. This allows for cells to mechanically communicate with each other [70,97]. Califano et al. [15] have shown that on polyacrylamide gels of sufficient compliance, cells self-organize into vascular-like networks, while they are unable to do so on very rigid substrates [14].
Some previous cell-based modeling has already been dedicated to mechanical cell-ECM interactions [7,20], where cells contract the matrix and in response align to each other. In [51], a continuum model where cell and ECM density dynamics are regulated by chemical and mechanical forces is presented that leads to network formation. However, in this model, strains in the matrix did not significantly influence network formation. To further investigate the influence of mechanical cell-ECM interactions via strains in the ECM, van Oers et al. [92] have developed a hybrid CPM and Finite Element Model to study network formation. The traction forces that the cells apply to the ECM are described with a model proposed by Lemmon and Romer [46]. This experimentally validated model treats the cytoskeleton as a single cohesive unit, as a result of which the cell forces that cells generate at each point depend on the local cell shape. The strains that are generated in the ECM are calculated using finite elements, where the finite elements correspond to the lattice sites of the CPM. Subsequently, the cells respond to the strains in the matrix. It is assumed that cells preferentially protrude in the direction of higher strain. This was implemented by adding the following bias to the Hamiltonian at the time of copying, with g(x, x ′ ) = 1 for extensions and g(x, x ′ ) = −1 for retractions, λ durotaxis is a parameter that describes the mechanical sensitivity of cells. v m = x − x ′ , is the direction of copying, and ϵ 1 and ϵ 2 , and v 1 and v 2 are the eigenvalues and eigenvectors of ϵ that represent the principal strains and strain orientation. Thus, extension in the direction of higher strain are promoted and likewise retractions are inhibited. The sigmoid function h(E) = 1/(1 + exp(−β(E − E θ ))), starts at zero, goes up when there is sufficient stiffness, and eventually reaches a maximum. This means that a certain level of stiffness, due to strain stiffening, E(ϵ) = E 0 (1 + (ϵ/ϵ st )1 ϵ≥0 ) is needed to cause a cell to spread.
On a single cell level, this model predicts that a single cell elongates due to a positive feedback loop of increasing traction and strain stiffening, as previously suggested by Winer et al [97]. Further, two cells in each others vicinity locally align. On a collective scale, these cell-level dynamics lead to vascular-like network formation (Fig. 18.3d); Cells are seen to elongate and locally align to form connected patches of aligned cells. Notably, network formation only occurs on substrates of intermediate stiffness and the simulated networks continuously remodel. Bridging events occur, where two groups of cells penetrate an existing lacuna, forming two lacunae. The paths that cells follow to divide a lacuna is directed by strain lines. Such bridging events have been observed in experimental conditions as well [92].

Modeling Sprouting During Angiogenesis
So far we have seen four independent mechanism that can lead to vascular-like network formation and thus give different explanations of the mechanisms of de novo vasculogenesis. A natural question then is whether these mechanism can also give rise to sprouting angiogenesis? In this section, we will explain how these four CPMbased models can drive sprouting from spheroids. In the next section, we will discuss how other CPM studies have contributed to investigating sprouting angiogenesis.

Sprouting-Like Behavior of Cells in de Novo models
Merks et al. [56,59] showed that cells that secrete and chemotact toward a chemoattractant sprout from a spheroid when either cells elongate or exhibit contact-inhibition chemotaxis. With just plain chemotaxis, sprouting was merely possible for a small range of diffusion constants or strong cell-cell adhesion. So, what gives these two mechanisms, contact-inhibited chemotaxis and cell elongation, the ability to drive chemotactic cells to sprout? The initiation of sprouts at the surface of the spheroid are thought to occur due to a buckling instability; cells in the core of the cluster are compressed due to the pressure the cells on the surface of the spheroid apply inward due to chemotaxis toward the increasing chemoattractant concentration inside the Fig. 18.4 Angiogenesis models. Simulation results of angiogenesis driven by a contact-inhibited chemotaxis [59], b cell elongation with chemotaxis [56], c preferential adhesion to elongated cells [87], and d by mechanical cell-matrix interactions [92] spheroid. Then, for the case of contact-inhibited chemotaxis, chemoattractant gradients at convex regions of the cell aggregate are more shallow than the gradients at concave regions. This makes it more likely for cells to protrude from convex regions of the surface, i.e., at the tips of sprouts (Fig. 18.4b). In spheroids of elongated cells, sprouts start to extend due to local alignment of cells (Fig. 18.4b).
Szabó and Czirók [87] also investigated under what conditions sprouting from a spheroid occurs. It turns out that the assumed mechanism of preferential adhesion to elongated structures suffices for the cells to sprout. Sprouting also occurs by adhesion only, but sprouts can quickly break down. Thus, this attraction stabilizes sprout extensions. Finally, it is argued that the inclusion of leader cells, that polarize and have a persistence in migratory direction, is required to obtain sprouting dynamics that are more similar to experimental results (Fig. 18.4c).
The mechanical model by Van Oers et al. [92] is also able to reproduce sprouting from a spheroid (Fig. 18.4d). Similar to vascularization, only sprouts are formed on substrates of intermediate stiffness. Due to random motility, one cell protruding from the spheroid increases the strain in front of it and subsequently follows it. This instigates a positive feedback loop of strain development and cells extending from the surface that are guided by the strain lines; forming the sprout.

CPM Models of Sprouting Angiogenesis
In the works by Bauer et al. [5,6] and Daub and Merks [23], VEGF and the ECM are incorporated to study how gradients of VEGF and properties of the ECM can influence sprout formation. The fibrous nature of the ECM influences cell migration in various ways. Fiber orientation directs cell migration, by contact guidance. Further, cells exhibit haptotaxis; migration toward higher ECM densities and haptokinesis; increased movement on intermediate ECM densities. In both models, sprout formation is investigated in the context of sprouting from a blood vessel toward a tumor secreting VEGF.
Bauer et al. [5,6] model the ECM geometry by including ECM fibers and interstitial fluids, as CPM pixels and frozen tissue-specific cells. The endothelial cells preferentially adhere to the ECM fibers and migrate toward the tumor. Endothelial cells interact with other tissue cells via adhesion. The vessel and tumor are located at opposite sides of the domain. The tumor secretes VEGF that diffuses and is degraded in the model domain and is taken up by the endothelial cells. Endothelial cells interact with the ECM by uptaking and degrading it, and chemotact toward VEGF. Haptotaxis is incorporated by high cell-fiber adhesion. In addition, a distinction between tip and stalk cells is made, by letting tip cells perform chemotaxis and degrade the ECM. Sprout migration is then made possible as tip cells degrade the matrix and stalk cells follow by means of haptotaxis. The model shows that speed, direction and branching of sprouting is dependent on ECM fiber density and composition.
Daub and Merks [23] investigated the effects of ECM densities on sprout morphology and branching by coupling the CPM with a PDE describing ECM density dynamics. The model set-up resembles the one used by Anderson and Chaplain [2], who used a stochastic discrete PCA-like model based on PDE discretization, to describe sprouting toward a VEGF secreting tumor. In [23], a VEGF gradient is presented to the CPM cells, to which cells respond by chemotaxis. Furthermore, VEGF induces the cellular secretion of proteolytic enzymes that degrade the ECM. Cells in turn respond to the ECM by haptokinesis and haptotaxis. Haptokinesis promotes the formation of branches and increases the sprout velocity on intermediate ECM densities. The degree of sprouting is most influenced by the haptotaxis parameter. Again, this work has showed the importance of cells interacting with the ECM properties to sprout formation.

Multiscale Models
The above CPM models of blood vessel formation asked how a single, stereotypic set of cell behaviors results in multicellular patterns. This can be an accurate representation of the situation in some in vitro cell cultures, but in realistic situations, i.e., in actual organisms, the situation is usually much more complex. Blood vessels typically consist of a number of cell types, including endothelial cells, pericytes, and smooth muscle cells, each of which require their own description using the CPM. Usually this is done by assigning each cell type a different value of τ (see Eq. (18.3)), and assigning different parameters to each cell type (e.g., by giving a different value of J to each combination of cell types; cf. Sect. 18.2.1).
The situation becomes more complicated if the cells change type depending on the signals they receive from adjacent cells, a common situation in biology. In angiogenesis two phenotypes of endothelial cells are distinguished, the tip cell that has many protrusions and is highly migratory but rarely divide, and the stalk cell that has few cellular protrusion and can proliferate. The differentiation into two types is mediated by the cell-cell contact signaling (e.g., through the Notch pathway [39]). To model this situation, the CPM is often extended with sets of coupled ODEs, where each cell (or spin σ) obtains its own set of ODEs. The ODEs can then also be coupled with the ODEs of adjacent cells, to model chemical signaling, or with sets of PDEs, e.g., to describe the diffusion and reaction dynamics of chemoattractants. To model cell differentiation, the CPM parameters of the cells (i.e., the target areas A(σ), the cell rigidity, λ(σ), the interfacial tension parameters J , and so forth) can be replaced for functions of the intracellular ODEs. As a result, the dynamics of the ODE can lead to changes in the behaviors and positions of the cells in the CPM, which can in turn affect the ODE, resulting in interesting multiscale dynamics. An example of this approach is given in Sect. 18.6.1.
In these examples, the cells in the CPM are still treated as homogenous structures, whereas in actual organisms the internal structure, e.g., the cytoskeleton, affect the behavior of cells in the tissue. The regulatory networks simply regulate the parameters of the cells. In a number of problems, it becomes important to describe the internal structure of the cells in more detail. To explain the dynamics of a type of highly motile skin cell, the keratocyte, Marée et al. [54] extended the CPM with an intracellular, dynamic model of the actin cytoskeleton, an approach that was later generalized to study the response speed to chemotactic cues in eukaryotic cells [53]. This work made use of an intracellular set of PDE's to describe the polymerization and orientation of actin filaments and of the enzymes regulating the polymerization rates. The polymerization model then biased the extensions and retractions in the CPM by modifying ∆H during the copy attempt (cf. Eqs. (18.5) and (18.9)). Another approach to include internal structure is the compartmental cellular Potts model. In this approach multiple spin domains are bundled together to form one biological cell. This approach was used by Boas and Merks [10] to model the formation of the lumen, i.e., the hollowing out of new blood vessels such that blood can flow through. Section 18.6.2 will briefly review this approach.

Sprouting Morphogenesis with Tip Cell Selection
A suitable example problem to illustrate the structure and dynamics of multiscale models that include a model of cell differentiation, is the selection of tip and stalk cells via the Delta-Notch molecular signaling pathway. Delta-Notch signaling acts as a lateral-inhibition mechanism, where a high expression of Delta activates the expression Notch in adjacent cells, which in turn suppresses the activity of Delta. Delta-Notch signaling is involved in a variety of processes in developmental biology, including the formation of body segments: somitogenesis [26], asymmetric cell division [4,72], neuronal plasticity [1,47]), and the initiation of angiogenesis [8]. Delta-Notch is mediated by interactions between Notch receptors and Delta/Serrate/LAG-2 (DSL) ligands [13]. In angiogenesis, extracellular VEGF has been shown to initiate the endothelial Delta-Notch signaling leading to the dynamic stalk-tip cell selection [8] (see Fig. 18.5).
Prokopiou and coworkers [68,69] have introduced a detailed, multiscale model of sprouting angiogenesis based on the CPM. In this model, each cell includes an ODE model of the Delta-Notch-VEGF signaling pathway, which acts to regulate task division between adjacent cells. We will review their model in detail here; similar approaches have been taken recently by Palm et al. [62] and by Boas and Merks [11].

Growth Factor and Extracellular Matrix Fields
Prokopiou's model considers the interaction between endothelial cells and the extracellular matrix (ECM). In the standard CPM (see, e.g., Sect. 18.2) the substrate, or medium, is represented as a homogenous cell covering the whole computational domain. Such a homogenous material can also be a suitable description of the ECM if the spatial inhomogeneity of ECM does not affect the problem under study. However in the case of EC migration, empirical evidence (e.g., in the developmental retina) showed that ECM form fiber bundle networks that are highly inhomogeneous and that ECs follow the tracks of fiber bundles. It is therefore necessary to model the ECM as a discrete field. In order to model such a non-homogenous ECM, the ECM fibers were distributed randomly. In particular, these fibers are modeled as a static field in the numerical domain. Each pixel in the numerical domain occupied by an ECM fiber is given a non zero (=1) value (and zero elsewhere). Cell tracking along the fibers is modeled as preferential adherence to the fibers. Thus, haptotaxis, the directional migration of cells up the ECM density field, is incorporated as an additional mecha-nism. In the CPM, haptotaxis can be implemented similar to chemotaxis (Eq. 18.5), with the main difference that the ECM field does not diffuse. During a copy attempt, the following term is added to ∆H ,

Subcellular Level: Modeling Lateral-Inhibition
At the subcellular level, tip cell differentiation is regulated via the Delta-Notch signaling pathway, which is activated by VEGF. The contact lateral-inhibition effect for the exchange of the endothelial (stalk-tip) phenotype is implemented using a modification of a well mathematical model proposed by Collier et al. [21], where a system of coupled ODEs describes the dynamic processes of Delta and Notch activation and inhibition between cells that are in contact with each other. Motivated by the experimental work of Lobov et al. [49], which showed that VEGF induces Delta in the retinal vasculature, Prokopiou and coworkers [68,69] extended the model of Collier et al. [21] to incorporate the contribution of VEGF (as defined in Eq. (18.11)), the non-dimensionalized form is: where D j , N j , represent the levels of Delta and Notch expression in cell j.
[VEGF j ] = (1/a j ) ω i VEGF ji is the average VEGF in a cell j; that is, the sum of VEGF at each pixel i inside cell j over the cell area, a j , where ω is the total number of pixels in cell j. V EG F h is the VEGF level at which the production rate of Delta is half maximal. The trans-Delta (D j ) is taken to be the sum over the immediate (contacting) neighbors i of cell j. P j is the perimeter of cell j, and P i j is the common area of cell j with its neighbor cells i, which is defined as The summation is over all pairs of adjacent sites in the lattice. Equation (18.12) describe (i) the activation of Notch production within each cell as a function of the levels of (trans-) Delta expressed by neighboring cells, (ii) the inhibition of Delta expression by Notch, and (iii) the activation of Delta production by extracellular VEGF. In the absence of VEGF signaling, there is no up-regulation of Delta and, therefore, no tip cell activation. Equation (18.12) was implemented using the Systems Biology Workbench and integrated within the CompuCell3D framework.

Coupling of ODE Model to CPM
To simulate the effect of the regulation by the signaling network on cell behavior, we let the level of Delta in the ODEs (Eq. (18.12)) determine the cell type τ ∈ {tip, stalk} (cf., Eq. (18.3)): if D(σ) > θ tip , the cell type becomes τ = tip, or else the cell type becomes τ = stalk. Each cell type is associated with a prescribed set of properties. The tip cells have a higher chemotactic coefficient than stalk cells (λ chem (tip) > λ chem (stalk)); the stalk cells if they are adjacent to tip cells, can grow by gradually increasing their target areas. The latter property is to implement the assumption that only stalk cells adjacent to tip cells proliferate, because the proliferation of all stalk cells would lead to a thick/swollen sprout and parent vessel.
To avoid any predefined or probabilistic rules of cell growth and division, we assign to each cell a clock φ(σ) that progresses only for stalk cells adjacent to tip cells. The clock progresses at a rate a = 0.01 h per MCS. In addition, considering that we want a cell to divide after doubling in size, the target volume of a cell (V ; cf. Eq. (18.3)) should grow by one initial target volume, V (t = 0) during one cell cycle of 17 hours, yielding a growth rate,

Simulation Results
The Notch-Delta pathway, the VEGF source, and the ECM heterogeneity all work in concert to sprouting angiogenesis. We first analyzed Eq. (18.12) to find the parameter values for α such that the homogeneous steady state become unstable. For a string of cells, the solution is a dynamic 'salt and pepper' pattern of cells with alternating high and low Delta values; for a 2D sheet of cells, the solution is a dynamic 'checker board' pattern. When we define a tip cell as a cell with a Delta-level above a threshold, the solutions then lead to a dynamic interchange of phenotypes between stalk and tip cells. Hence in our simulations, the phenotype distribution of ECs along the capillary sprout is determined by two main mechanisms: the astrocyte-derived VEGF that activates the Delta activity in each cell, and the Notch-Delta signaling pathway that yields the 'salt-pepper' pattern. Following the tip/stalk selection, the ECs then migrate chemotactically to VEGF distribution and haptotactically to ECM distribution. Figure 18.6 shows the development of multiple tip cells, each can potentially lead the formation of a sprout; when the head tip cell of two growing sprouts meet, the Notch-Delta pathway re-establishes the tip, resulting in the apparent fusion of two sprouts into one, in a process termed anastomosis. We summarize the effect of different VEGF and ECM profiles (Table 18.1) on the resulting morphology of the capillary sprouts. Figure 18.7 shows representative snapshots of sprout evolution in each scenario.

No VEGF Gradient (Scenarios 1 and 2)
In scenarios 1 and 2, there is no VEGF gradient. A sufficiently high level of VEGF activates the Notch-Delta pathway, and leads to selection of tip cells (yellow) and Homogeneous VEGF and heterogeneous ECM 3.
Static VEGF gradients and homogeneous ECM 4.
Static VEGF gradients and heterogeneous ECM 5.
Heterogeneous VEGF and heterogeneous ECM stalk cells (red). However the cells do not receive directional guidance from a gradient of VEGF, resulting in a much reduced migration. Figure 18.7 shows that cell proliferation and elongation are undirected and, therefore, stalk and tip cells evenly fill the space. This morphology was observed in experiments [31,60], where a spatial gradient in VEGF was removed in the retina, by increasing expression levels of VEGFA in transgenic mouse models. In scenario 2, the addition of a non-uniform ECM has a weak effect, because the ECM does not offer an overall gradient.

Static VEGF Gradient (Scenarios 3 and 4)
In these two scenarios, we incorporated static VEGF gradients, which eventually lead to either a swollen (scenario 3) or a thin (scenario 4) sprout formation. The results of scenarios 1-4 look quite similar up to approximately 12 h. The sprouts are dominated by single, elongated tip cells. However, differences become visible at later time points. Particularly, in scenario 4, cell proliferation is focused onto a single sprout as a result of the VEGF gradients and the heterogeneous ECM.

Dynamic VEGF from Single Source (Scenarios 5 and 6)
Here, a fixed astrocyte (VEGF source) is responsible for the VEGF gradients. Figure 18.7 (scenarios 5 and 6) demonstrates the model's ability to reproduce realistic capillary sprout morphologies (up to ∼38 h). Scenario 5 (with homogeneous ECM) can give a polarized sprout, but the emerged sprout in scenario 6 (with heterogeneous ECM) has the right extension speed (∼1.6 µm/h) as it was evaluated from our experimental (unpublished) data. Therefore, we suggest that scenario 6 provides a close approximation to a growing vascular sprout. However, since the astrocyte cannot move away, scenario 6 does not allow for the formation of longer sprouts, because at late time points (85-100 h) a mass of cells starts surrounding the astrocyte.

Lumen Formation
The work by Scianna [74] and the model by Boas and Merks [10] illustrate the use of the so called compartmental CPM [74,82,86] to treat subcellular structures during angiogenesis. In the compartmental CPM, the Potts domains (clusters of  spins, σ) represent parts of cells, rather than individual cells. The compartments are then bundled together to represent one biological cell using a cluster identifier, ξ.
All spins belonging to the same cell then have the same cluster identifier, ξ(σ). Additional constraints can be imposed on the whole cell (see, e.g., Ref. [10]), on individual components [75] or both [10]. Example applications include a model of the nuclei of endothelial cells [74], in particular the way how the nuclei can slow down the migration of the cells in the ECM if the nuclei are larger than the typical pore size [75], and the model of lumen formation that we will review in more detail here.
Once new blood vessels are formed, they must hollow out to allow the perfusion of blood. The mechanisms of hollowing or lumen formation have been debated for centuries. Experimental research has led to two main hypotheses: vacuolation [9,25,41,96] and cell-cell repulsion [85]. During vacuolation, vacuoles are suggested to form by the fusion of pinocytotic vesicles. Initially, lumens were thought to form intracellularly by spanning the cell with a large vacuole that then fuses to the cell membrane on both sides of the cell [25,41]. Later, lumens were also suggested to form extracellularly by the secretion of vacuoles between cells [9,96]. During cellcell repulsion, cell membranes of adjacent cells are suggested to repulse each other to form an extracellular lumen between the cells [85]. Both hypotheses are supported by strong experimental evidence, leaving the debate unresolved. To address this debate, Boas and Merks [10] developed a computational model of lumen formation that can represent both hypotheses.
The lumen formation model is initialized with twelve endothelial cells in a branched blood vessel, surrounded by immobile extracellular matrix (ECM). Each cell is modeled as a cluster of CPM compartments (σ) with the same cluster identifier (ξ) to allow for polarization of the cell membrane and for the formation of vesicles and vacuoles within the cell. The cell polarizes into two cell membrane compartments and a cytosol compartment upon contact with the ECM, representing cell membrane polarization by integrin signaling from the ECM. All membrane pixels that are in contact with the ECM form a membrane compartment of type τ (σ) = basolateral. The adjacent second neighbor order membrane pixels hereof are added to this compartment to represent tight junctions between cells, and the rest of the membrane becomes the second membrane compartment of type τ = apical. The membrane is repolarized every other time step. To mimic cell-cell repulsion, apical membranes of opposing cells are assigned a high adhesion energy.
During vacuolation, membrane pixels that internalize into the cytosol compartment have a probability to become cell compartments of type vesicle to represent pinocytosis. These single-pixel vesicles move through the cell following a biased random walk, by swapping the position of a vesicle with a neighboring pixel. Acceptance of a swap depends on a constant probability P A multiplied by a Boltzmann probability P Boltzmann (∆H ), with ∆H the change in effective energy resulting from changes in adhesion energy between compartments due to the swap. Vesicles prefer to adhere to other vesicles and vacuoles. Once vesicles meet, they fuse together into a single compartment of type vacuole, which moves by regular CPM dynamics. Vesicles and vacuoles are secreted when in contact with the apical membrane, forming a new extracellular compartment of type luminal fluid. Upon contact, luminal fluid compartments fuse into a single lumen.
Continuous lumens can be formed in the model through the branched blood vessel by vacuolation as well as by cell-cell repulsion (Fig. 18.8). However, lumen formation is far more robust to parameter values changes when the two hypotheses are combined, suggesting that the two hypotheses work synergistically. Vacuolation can help lumen formation by cell-cell repulsion by piercing cells and by enlarging the luminal space in-between cells. The cell-cell repulsion hypothesis assists One may question synergy of the two hypotheses as experimentalists mostly find evidence for one or the other hypothesis. It is important to realize that lumen formation by vacuolation is mostly studied in small intersegmental vessels (ISV) of zebra fish, while cell-cell repulsion is mostly studied in aortae of mice. Interestingly, when lumen formation by synergy of the two hypotheses is performed in the model initialized with a one-cell thick vessel, the resulting lumen formation visually resembles vacuolation. In contrast, when lumen formation by synergy of the two hypotheses is performed in the model initialized with a multicell thick vessel, the resulting lumen formation visually resembles cell-cell repulsion. In conclusion, the computational model of lumen formation suggests that vacuolation and cell-cell repulsion work synergistically and that the discrepancy between observations of different experimental groups might be explained by the vessel sizes they are studying.

Integrating Angiogenesis Models into CPM Models of Organogenesis
Angiogenesis is a key process in many developmental and pathological processes. For this reason, the simple models of endothelial cell-cell interactions have been integrated in larger models of organ development and tumor growth. Shirinifard et al. [79] have integrated an angiogenesis model similar to the one proposed by [59] with a tumor growth model, where the growth of tumor cells was made dependent on the availability of oxygen. Kleinstreuer et al. [42] integrated a cellular Potts model of in vitro angiogenesis with a large dataset of the US Environmental Protection Agency (EPA) of pesticides and their effects on vascular morphogenesis. By linking the adverse effects of pesticides to individual cell behaviors of the cell types involved in vasculogenesis, they could construct a first toxicological, predictive model based on the cellular Potts model. A further example of how cellular Potts models of angiogenesis can be integrated into larger models of tissue development is on age-related macular degeneration, by Shirinifard et al. [80]. Age-related macular degeneration (AMD) is the main source of vision loss in the elderly and a looming epidemic for our aging society. There are two basic forms of AMD, the "dry" form and the "wet" form. In dry AMD, the layer of retinal pigment epithelial cells (RPE) in the macula degenerate and die (atrophy). These RPE cells support the light sensitive photoreceptor cells that are critical to vision. Dry AMD can progress slowly and culminate with the more advanced stage called Geographic Atrophy, where a patch of photoreceptor cells die off. The wet AMD is due to the abnormal blood vessels (known as choroidal neovascularization or CNV) growing under the retina and macula. These new blood vessels may then bleed and leak fluid, causing the macula to bulge or lift up from its normally flat position, thus distorting or destroying central vision. Under these circumstances, vision loss may be rapid and severe.
In CNV, after capillaries initially penetrate basement membrane under the RPE (called the Bruch's membrane or BrM), invading vessels may either regress or expand. Clinically, during early and late CNV, the expanding vasculature usually spreads in one of three distinct patterns: in a layer between BrM and the retinal pigment epithelium (sub-RPE or Type 1 CNV), in a layer between the RPE and the photoreceptors (sub-retinal or Type 2 CNV) or in both loci simultaneously (combined pattern or Type 3 CNV). Most previous studies hypothesized that CNV primarily results from growth-factor effects or holes in BrM, but failed to explain the initiation nor progression patterns of CNV. Shirinifard et al. [80] used 3D CPM of the normal and pathological maculae to recapitulate these three growth patterns (Fig. 18.9). The key feature of these tissue models are the adhesions within and between different tissue layers: BrM, RPE, and photoreceptor outer segment (POS) (Fig. 18.9a), in addition to endothelial cell dynamics, VEGF dynamics and MMP degradation of ECM. These models aimed to test the hypothesis that CNV results from combinations of impairment of adhesion, in particular: RPE-RPE epithelial junctional adhesion, adhesion of the RPE basement membrane complex to BrM (RPE-BrM adhesion), and adhesion of the RPE to the photoreceptor outer segments (RPE-POS adhesion). Figure 18.9b shows a time sequence of snapshots from a typical simulation of Type 3 CNV, where new blood vessel invades both under and above the RPE layer. Results from all combinations of adhesion parameters were summarized into tables and risk maps.   [80] under the terms of the Creative Commons Attribution License apical RPE-POS or epithelial RPE-RPE adhesion (e.g., due to inflammation) results in Early sub-retinal CNV. (5) Simultaneous reduction in RPE-RPE epithelial binding and RPE-BrM adhesion results in either sub-RPE or sub-retinal CNV which often progresses to combined pattern CNV. These findings suggest that defects in adhesion dominate CNV initiation and progression. This conclusion is both novel and surprising, but coherently explain the heterogeneous range of CNV growth patterns and dynamics.

Conclusion
In this chapter, we have introduced the cellular Potts model and discussed how it can be seen as a special case of PCA. In contrast to the formal definition of PCA, the cellular Potts model is asynchronous, and its rules are not strictly local. The dynamics, as guided by the Hamiltonian (Eq. (18.3)), depend on the local neighborhood of the lattice sites, as well as on the properties of the whole biological "cell," i.e., the set of all lattice sites x that have the same state, or spin σ(x). Examples of such non-local dependencies include the volume constraint Eq. (18.3) and the length constraint (Eq. (18.7)). As an advantage of this approach relative to traditional PCA that represent biological 'cells' with individual lattice sites, cells in the CPM can assume arbitrary shapes, which can be given by the model (see e.g., Sect. 18.4.2) or change dynamically during the simulations. The resulting simulation images and movies are often perceived by biological researchers as "realistic," allowing for one on one visual and quantitative comparison with microscopic data.
Of course, such flexibility comes at a cost. With the current speed of serial processors, typical simulations of the CPM are fast enough that large scale parameter studies can be performed [64], although full three-dimensional simulations can be limited to at most several million cells for individual simulations. What is currently largely out of reach are formal mathematical analyses of the CPM similar to those performed for PCA, making it practically impossible to generalize or proof any insights obtained with the CPM beyond what was tested numerically for individual parameter sets. Fortunately first attempts to formalize the CPM have been made, as shown in Chap. 19 (see also Ref. [94]). Apart from its non-locality, another mathematical limitation of the CPM is its required asynchronicity. Apart from complicating formal treatment, it hinders its implementation on graphical processing units (GPUs). To reach optimal speedup, GPUs rely heavily on the synchronicity and the locality of the algorithm. Although GPU-implementations of the CPM have been proposed [90,98], a fully synchronous, local reformulation of the CPM would help dramatically speedup CPM-simulations.
After introducing hybrid CPMs, in which the CPM dynamics affects the kinetics of a PDE model and vice versa, we illustrated the applicability of the CPM to biomedical problems. Here we focused on the modeling and simulation of blood vessel growth: angiogenesis and vasculogenesis. After discussing in Sects. 18.4 and 18.5 how CPMs have been instrumental in proposing and analyzing new hypothesis for the cell behavior that is responsible for the formation of blood vessel like structures, Sect. 18.6 showed how such models can be incorporated into more complete, multiscale models. Such multiscale models typically contain more detailed models of the intracellular kinetics, implemented using ODEs or using a compartmental CPM. Finally, we have shown in Sect. 18.6.3 how these can be incorporated into larger scale models of organ development or disease progression. Mathematically, such models can become complicated: we have coupled systems of PDEs, ODEs and compartmental CPMs, where several of such models might operate at of the system.
In such cases, model validation might become a serious concern. First, the behavior of the system becomes difficult to determine. Parameter sweeps are key tools for determining the behavior simulation models [64], but they can only be performed starting from one or a few sets of nominal parameter values, as a result of which some interesting or false behavior of the model might be missed. To get better insight into the whole parameter space, useful methodology includes global sensitivity analyses, which have recently been tested on a simple CPM of vascular morphogenesis [12], Second, among the plethora of potential biological mechanisms represented by our models, the ones that best describe the actual mechanism must of course be selected. Sections 18.4 and 18.5 showed that a range of different mechanism is able to describe vasculogenesis and angiogenesis (see also Ref [58]). Since these models all result in sprout and branch formation reminiscing the experimental data, they are all plausible explanations for these phenomena. To this end, experiments should be designed in order to further validate these models in order to rule out some of the different hypotheses. We must however keep in mind that it is possible that different mechanisms operate in different tissues or time periods in development. Or, most likely, different mechanisms work together in order to effectively create and stabilize vessels. In order to gain a better understanding of angiogenesis, we must figure out how and when certain mechanisms play a role and how they influence each other. Computational modeling using the CPM serves as a good starting point to get insight into the roles and interactions of alternative mechanisms of vasculogenesis in a combined model, driving the development of new, testable hypotheses.