Cell-Oriented Modeling of Angiogenesis

Due to its significant involvement in various physiological and pathological conditions, angiogenesis (the development of new blood vessels from an existing vasculature) represents an important area of the actual biological research and a field in which mathematical modeling proved particularly useful in supporting the experimental work. In this paper, we focus on a specific modeling strategy, known as “cell-centered” approach. This type of mathematical models work at a “mesoscopic scale,” assuming the cell as the natural level of abstraction for computational modeling of development. They treat cells phenomenologically, considering their essential behaviors to study how tissue structure and organization emerge from the collective dynamics of multiple cells. The main contributions of the cell-oriented approach to the study of the angiogenic process will be described. From one side, they have generated “basic science understanding” about the process of capillary assembly during development, growth, and pathology. On the other side, models were also developed supporting “applied biomedical research” for the purpose of identifying new therapeutic targets and clinically relevant approaches for either inhibiting or stimulating angiogenesis.


INTRODUCTION
Computational models and simulations play many roles in science [1]. They are used to make precise and accurate predictions and to summarize data. They are used as heuristic approaches for designing experiments or to demonstrate surprising and counterintuitive consequences of particular forms of systematic organization. As far as biological systems are concerned, attempts have been made at modeling many different processes. In this respect morphogenesis represented a quite common target of modeling efforts. They addressed situations ranging from the formation of bacterial [2] and mesenchymal cell [3] aggregation patterns and Dictyostelium morphogenesis (see [4]) to tumor growth [5][6][7], limb patterning [8], and gastrulation [9].
Due to its significant involvement in various physiological and pathological conditions [10], angiogenesis (the development of new blood vessels from an existing vasculature) represents an important area of the actual biological research and a field in which mathematical modeling proved particularly useful. Perhaps the first mathematical analyses of vascular networks can be found in the seminal work of Wilhelm Roux (see [11]) and in the classic work of Thompson [12] where he studies "...a number of interesting points in connection with the form and structure of blood vessels." However, is within the past two decades that the application of mathematical and computational models has significantly supplemented experimental approaches in this field and enhanced our understanding of the main factors regulating the vascular pattern formation. One way to categorize the existing set of published models is according to the spatial scale they were developed to encompass [13,14].
Some computational studies focused on the "molecular level," building models of the intracellular dynamics (see [15,16]), such as signaling phenomena and gene expression. The coupling of many detailed single-cell models was suggested by some authors [17] as a possible modeling strategy to reproduce multicellular phenomena. However, very accurate models of a single-cell (see [18]) can, at best, treat clusters formed by a quite low number of cells.
On the other side several models have reproduced vessel-like patterns consistent with those observed in vitro [19][20][21][22][23][24] or in vivo (see for instance [25][26][27]) by following a "tissue level" approach (see [28]), in which the system is treated as a continuous substance, and the involved cells are described in terms of densities (using partial differential equations). Continuum models of this type average out the behavior of the individual elements and are capable of efficiently capturing features of angiogenesis at a "macroscale" (such as average sprout density, network expansion rates, etc.). They, however, are unable to provide detailed information at a "microscale" concerning the actual structure and morphology of the capillary network. In fact, the self-organization of the endothelial cells (EC) leading to the formation of new capillary branches is mainly the result of several intimately linked single-cell behaviors [29].
Thus, working at too coarse or fine a level of detail makes quite hard an accurate modeling of the complex process of angiogenesis. For this reason, "cell-centered" approaches, working at a "mesoscopic scale" and treating the cell as the fundamental module of development, have been devised [30]. They also proved quite useful to build multiscale models of the process, providing a sort of natural interface between "molecular level" and "tissue level" modeling.
This specific modeling strategy and the role it can play in the study of the angiogenic process are the focus of the present paper.

A CELL-CENTERED APPROACH TO MODEL MORPHOGENESIS
The underlying principles of the "cell-centered" approach to modeling have been extensively discussed by Merks and Glazier [30], and its main characteristics will be only briefly recalled below.
The key concept on which cell-centered models are based is to assume the cell as the natural level of abstraction for mathematical and computational modeling of development. Thus, to a first approximation, the cell's internal properties (i.e., the details of the intracellular processes) are not explicitly taken into account and only its essential behaviors (such as movement, division, death, differentiation, adhesion, and secretion of chemicals) are considered. A significant advantage of this strategy is the relative simplicity of the models it generates. Systems composed by a quite large number of cells (up to 10 5 -10 6 cells) can be simulated, opening a concrete possibility to study how tissue-level processes could emerge from the collective dynamics of multiple interacting cells. It follows that cell-centered methods appear particularly suitable to investigate morphogenesis as also illustrated by very recent studies [31,32] showing how cell shape, most likely sensed by the mitotic spindle, serves as a major determinant of future cell and tissue development.
To achieve this goal, some methodological steps are required, in which cell-centered simulations are compared with experimental observations to identify the minimal set of single cell behaviors needed to produce certain tissue-level patterns. A typical flow-chart for this protocol of computational prediction and experimental validation is provided in Figure 1.
As far as the mathematical modeling and simulation techniques are concerned, several cell-centered computational strategies have been proposed to study morphogenesis.
Some of them were focused on tissue processes in which cells keep a fixed position with respect to each other [33]; others considered mobile cells and the physics of the adhesive forces between cells and the extracellular matrix (ECM) to simulate aggregates of thousand of cells (see [34]). In the Lagrangian Monte-Carlo method proposed by Drasdo et al. [35], for instance, attraction, compression, and bending energies determine movements of spheroidal cells to simulate cleavage and gastrulation [9] and avascular tumor growth [5], while a statistical mechanics-based approach was developed by Newman and Grima [36] to model chemotactic cell-cell interactions and to study cell ensembles analytically.
However, a quite popular strategy to model the self-organization of mobile cells is using cellular automata [37][38][39]. For this reason this computational technique will be here described in more detail. Cellular automata (CA) consist of discrete particles that occupy some or all sites of a regular lattice [40]. Each particle is characterized by one or more internal state variables and a set of rules describing the evolution of their state and position. Both the movement and change of state depend on the current state of the particle and those of neighboring particles. Again, the evolution rules may either be discrete or continuous, deterministic or probabilistic, and usually they are applied in time steps, following a synchronous or stochastic updating scheme. Philosophically, CA are attractive because they show some analogy with biological systems. In fact, their large-scale behaviors are completely self-organized [41,42]. An individual cell does not carry a road map, it can only respond to signals in its local environment. Furthermore, they need not privilege any single cell as pacemaker or director: all cells are fundamentally equivalent.
A relatively simple type of CA models is the so-called lattice-gas-based CA (LGCA) [39]. In LGCA individual particles on a discrete grid represent cells, each characterized by a velocity determined by the local interactions the cell experiments. At each time step, each cell will move to a neighboring site according to the velocity that cell had. Thus, in their biological applications LGCA treat cells as point-like objects with an internal state but no spatial structure. As a consequence LGCA models can be convenient and efficient for reproducing qualitative patterning in colonies where cells retain simple shapes during migration. Eukaryotic cells, however, often move by remodeling their cytoskeleton and changing their shapes. Since in some cases shape change significantly influences patterning, a modeling approach that takes into account cell shape is required. In this respect, a more efficient and complex CA is the Cellular Potts Model (CPM), in which a cell consists of a domain of lattice sites, thus describing cell volume and shape more realistically. Originally it was developed by Graner and Glazier [43] to simulate the cell rearrangement resulting from cell adhesion, in order to quantitatively simulate cell-sorting experiments. However, a number of cell behaviors can be quite easily implemented in this computational framework, and improvements to the CPM included the possibility to model cell growth, cell division, apoptosis and cell differentiation, chemotaxis, extracellular materials, and cell polarity (see [30]). The basic characteristic of the CPM is to represent the cell behaviors of interest in the form of terms within a generalized energy function which also includes the interactions with the ECM and parameters constraining individual cell behavior. As an example, a simple form of CPM is illustrated in Figure 2. Cells are represented on a rectangular numerical grid as patches of lattice sites, x, with identical nonzero indices σ x , while an index value 0 identifies the sites corresponding to the extracellular space. Grid points at patch interfaces represent cell surfaces, and the interaction between cell surfaces is modeled by defining coupling constants J σ x ,σ x representing the adhesion energy involved in the specified interaction. Each cell also has a set of attributes, including a "target" area and elongation, which poses some constraint on the possible cell shape changes. Thus, the following "energy function" can be defined for this system: where x represents the eight neighbors of x , λ A and λ L represent resistances to changes in size and elongation; respectively, A σ and L σ are the "target" values for cell area and length, a σ and σ are the actual cell area and length values, and the Kronecker delta is δ . The parameters involved in (2.1) can be estimated by specific experiments or based on biological considerations, and the dynamics of such an energy formalism can be solved using a variety of minimization methods, as, for instance, the well-known Metropolis or Kawasaki algorithms (see [30,39] for reviews).
LGCA-based models proved useful to describe the ripple formation in myxomycetes [47] and : Schematic representation of the Cellular Potts Model [43]. Cells are represented on a numerical grid as domains of pixels with identical index σ i (shown as a specific shade of gray), while the extracellular matrix is the set of the remaining pixels (white pixels) having σ = 0 by convention. Connections between neighboring lattice sites of unlike index (some of them are indicated with double-headed arrows) represent membrane bonds, with a characteristic bond energy J , which depends on the pair of objects in contact and determine the strength of their adhesion. Furthermore, because biological cells generally have a fixed range of sizes and shapes, additional elastic energy terms are considered whenever deviations from a target volume or elongation occur. All these contributions lead to (2.1), representing the "energy" of the system at each time instant. Cytoskeletally driven membrane fluctuations can be mimicked by randomly choosing a lattice site x and copying its index into a randomly chosen neighboring lattice site x (single-head arrow).
germinal center dynamics [48], and a quite wide range of biological problems (see [49][50][51][52][53]) were addressed with the CPM. In particular, as detailed in the next section, this type of cell-oriented strategy to modeling played a significant role in studies focused on angiogenesis.

CELL-CENTERED MODELS OF ANGIOGENESIS
Cell-oriented computational models of angiogenesis can be categorized around some key questions they have been developed to answer. From one side, they have generated "basic science understanding" about the process of capillary assembly during development, growth, and pathology. On the other side, models were also developed with the intention of supporting "applied biomedical research" for the purpose of identifying new therapeutic targets and clinically relevant approaches for either inhibiting or stimulating angiogenesis. These two targets motivating the development of the models will be the central thread of the present section.
One of the first, and most frequently cited, cell-centered models of angiogenesis has been developed by Stokes and Lauffenburger [54], who simulated individual cell movements by considering cell motility and chemotaxis as partially stochastic events. They used the model to assess microvascular endothelial cells migration in the presence or absence of acidic fibroblast growth factor (aFGF), and realistic capillary network structures were generated by incorporating rules for sprout branching and anastomosis. Although the model incorporated random motility and chemotaxis as mechanisms for cell migration, no account was taken of the interactions between the endothelial cells and the ECM. To account for this specific point, more recent models of cell migration improved the accuracy of the simulations by using force-based dynamics approaches to simulate internally generated forces and external traction forces, as well as matrix compliance and ECM stiffness [55].
Key morphological events involved in new vessel formation can be experimentally investigated by in vitro studies analyzing the endothelial cell self-organization in vitro [29,56]. In this context an important supporting tool for the interpretation of the observed patterns is represented by the CPM-based two-dimensional model by Merks et al. [57] simulating the process of in vitro vasculogenesis or the assembly of human umbilical vein endothelial cells (HUVEC) into networks of connected cells in a Matrigel environment (see [58]). The model considered a set of single cell behavior, including cell adhesion (between cells and with the ECM), chemotaxis, and cytoskeleton rearrangement. In this CPM cells are represented on a rectangular numerical grid as shown in Figure 2. By repeatedly replacing a value at cell interface by a neighboring grid point's value, it is possible to mimic active, random extension of filopodia and lamellipodia. If the resulting variation in effective energy is negative, then the cell change will be accepted, otherwise it will be accepted with the Boltzmann-weighted probability. The preferential extension of filopodia in the direction of chemoattractant gradients was implemented by including in (2.1) an extra reduction of energy whenever the cell protrudes into an area with a higher concentration of chemoattractant: where x is the neighbor into which site x moves (i.e., copies its value), χ is the strength of the chemotactic response, and c(x) is the local concentration of the chemoattractant. At each time instant, the concentrations c(x) were estimated from the following diffusion partial differential equation (PDE): where α is the rate at which the cells release chemoattractant, ε is the clearance rate of the chemoattractant, and D its diffusion coefficient. The Kronecker delta simply indicates that the release occurs at the cell locations, while the factor is cleared in the extracellular space. As shown in Figure 3, with a proper choice of the parameters, the model generates cell patterns in close agreement (from both a qualitative and quantitative point of view) with those generated in vitro by unstimulated HUVEC, suggesting that the three considered single-cell behaviors are essential for correct spatiotemporal vasculogenesis in vitro.
The same modeling approach was used by our group to analyze and interpret the results of in vitro angiogenesis experiments in conditions involving cell stimulation with proangiogenic factors or performed with nonendothelial cells potentially able to differentiate towards an endothelial phenotype. In the first study [59], a cell-centered mathematical modeling approach was used to determine essential cellular behaviors for pattern formation when human saphenous vein endothelial cells are stimulated by the pro-angiogenic factor adrenomedullin (AM) [60]. Cell culture measurements provided the key parameters to customize the model. When put to the test, the simulated pattern and morphometric parameters closely matched that of untreated EC, confirming that cell elongation, in conjunction with autocrine secretion of a chemoattractant, results in a cell-shape-dependent motility representing the key factor driving the formation of vascularlike morphologies by EC in vitro. However, the model failed to predict patterns of EC cultured with AM, revealing that it lacked input from an important cell behavior. Hypothesizing that the missing ingredient was cell proliferation, the model was extended to include it and called upon to estimate the percentage increase in cell number that would yield observed patterns. Remarkably, the proliferation rates predicted by the model showed consistency with bromodeoxyuridine incorporation experiments performed to verify such a prediction. In another study [61], aimed at investigating the vasculogenic potential of bone marrow macrophages from patients with multiple myeloma (MMMA), when seeded on Matrigel, a number of MMMA rapidly changed their morphology and developed an elongated shape. After 18 hours, the formation  [57,59]) for in vitro and in silico capillary-like patterns illustrating that the time course of pattern formation by HUVEC is well captured by the model. of a pattern consisting of cord-and tubular-like structures was observed, sometimes arranged to form closed meshes. When biophysical parameters consistent with the available experimental evidence were used to customize the model, it was quite accurate in quantitatively reproducing the observed in vitro patterns, provided that the possibility for cell differentiation was included in the model. In particular, it indicated that about 30% of the seeded MMMA were differentiated towards an endothelial phenotype, suggesting that in multiple myeloma a quite high number of MMMA could become involved in the process of capillarization by converting into a cell type at least similar to the endothelial one. Following the above-mentioned CPM-based modeling, Merks et al. [62] also showed that including VE-cadherin-mediated contact inhibition of chemotaxis in the simulation causes randomly distributed cells to organize into networks and cell aggregates to sprout, hence, reproducing aspects of both de novo and : Schematic summary of cell-based simulations of tumor angiogenesis (see [65,66]). Tumor cells secrete VEGF that stimulates EC, and the model distinguishes between tip and stalk EC phenotypes. The tip cells respond by moving chemotactically towards higher concentrations of VEGF using the matrix fibers for support. They are also capable of degrading the ECM. VEGF-stimulated stalk cells can proliferate and/or move behind a tip cell. Thus, vessels emerge that follow the chemotactic path of the endothelial tip cells.
As anastomoses occur, a network of vessels is formed. In the model by Mahoney et al. [66], oxygen is also secreted from the endothelial cells that belong to closed loops. It diffuses through the medium and reaches the tumor, where it is consumed by the tumor cells. Cellular dynamics (grey boxes and solid line arrows) are modeled by using CPM, the VEGF, and O 2 profiles by continuous models based on the diffusion equation (white boxes and dashed line arrows).
sprouting blood-vessel growth. This study, therefore, further confirmed the CPM as a potentially very helpful tool to investigate the whole spectrum of patterns formed during angiogenesis. As far as in vivo studies are concerned, a popular experimental setup for studying the sprouting of new vessels is the corneal pocket model, in which exogenous growth factors can be supplied in a controlled manner to induce reproducible angiogenic sprouts from the limbic vessels. This experimental model has recently been the subject of a cell-based mathematical model [63] allowing for a detailed study of the relative roles of EC migration, proliferation, and maturation in sprouts development. It showed that cell elasticity and cell-to-cell adhesion allow only limited sprout extension in the absence of proliferation, and the maturation process combined with bioavailability of VEGF can explain the localization of proliferation to the leading edge, consistently with experimental observations. The vascularized phase of tumor growth has dominated as the most common context in which to develop mathematical and computational models of angiogenesis. In this respect, the above-mentioned cellcentered model by Stokes and Lauffenburger [54] predicted, for the first time, that chemotaxis is needed to orient vascular growth toward the tumor. More complex cell-based models including a number of key events of angiogenesis (such as the migratory response of endothelial cells to cytokines secreted by a solid tumor, endothelial cell proliferation, endothelial cell interactions with ECM macromolecules, such as fibronectin, and capillary sprout branching and anastomosis) were proposed (see [13,64]). They provided capillary networks with a very realistic structure and morphology, capturing the early formation of loops, the essential dendritic structure of a capillary network, and the formation of the experimentally observed "brush border." A sophisticated cell-centered model of tumor angiogenesis was developed by Bauer et al. [65]. It describes diffusion, uptake, and decay of proangiogenic factors secreted by tumor cells and was used to understand the role of cell-cell and cell-matrix dynamics in regulating tumor angiogenesis. The model incorporates both discrete and continuous approaches (see Figure 4): a PDE describes diffusion and decay of tumor-secreted VEGF, while a CPM is used to describe EC migration, growth, division, and adhesion, as well as ECM degradation. Notably, this model is the first to capture anastomosis and branching without needing to predefine rules for these events: these properties emerged from the independent behaviors of the individual simulated EC. Furthermore, the model proved useful to address several questions, including the impact of ECM-binding affinity of VEGF on capillary morphology, the rate of capillary sprout elongation, and to what extent the composition of the stroma (ECM density and anisotropy) influences angiogenesis.
In almost all the above-mentioned models, the intracellular events are not modeled explicitly, and the information concerning the intracellular dynamics is embedded in the model parameters. Cell-centered approaches, however, can be extended to generate more realistic multiscale models of the complex process of angiogenesis. A recent example was provided by Scianna et al. [67]. The model spans three fundamental biological levels: at the extracellular level a continuous model describes secretion, diffusion, uptake, and decay of the autocrine VEGF, at the cellular level a CPM reproduces cell dynamics (migration, adhesion, chemotaxis), and at the subcellular level a set of reaction-diffusion equations describes a simplified VEGFinduced calcium-dependent intracellular pathway. The results agree with the known interplay between calcium signals and VEGF dynamics and with their role in malignant vasculogenesis.
Moving from the basic science of angiogenesis to the applied biomedical research, a number of celloriented models were developed to support the search for therapies and/or technologies aimed at favouring or inhibiting angiogenesis.
The development of tissue-engineered constructs greater than about 1 mm 3 is limited by the necessity to overcome oxygen diffusion limitations. Thus, the development of novel approaches for engineering microvascular networks ex vivo or inducing their ingrowth upon implantation of the construct is imperative. Jabbarzadeh and Abrams [68] developed a model of VEGF-mediated EC chemotaxis through a porous membrane in response to three different VEGF presentation strategies in order to assess which one could lead to the most extensive vascular coverage of the construct.
As far as the inhibition of angiogenesis is concerned, the early stages of tumor angiogenesis, in which EC escape the parent vessel and invade the ECM, were the focus of a cell-based mathematical model by Plank and Sleeman [69] with the aim to study the antiangiogenic potential of pharmacological strategies based on angiostatin.
A high fidelity simulation model of angiogenesis induced by solid tumors was developed by Mahoney et al. [66] as an evolution of the above-mentioned model by Bauer et al. [65]. The aim was to identify specific medically relevant intervention targets. The simulation system integrates (see Figure 4) the following: (a) a CPM that captures mechanisms of endothelial cell growth, cell adhesion, ECM fiber adhesion and degradation, and tip cell chemotaxis and haptotaxis, (b) a continuous model of VEGF secretion from the tumor, diffusion through the stroma (host tissue), and endothelial cell uptake and activation, (c) a flow model that estimates blood flow through the irregular network of vessels that emerge during angiogenesis, and (d) a continuous model of oxygen secretion from vessel loops, diffusion through the stroma, and uptake by the tumor. This model captures behaviors such as vessel branching, loop formation (anastomosis), progression and termination of tip movement, and activation and growth of new vessels. All these complex behaviors emerge from interactions among the simpler, biologically relevant component mechanisms of the model. The results of the simulations showed the effectiveness of this computational method in finding interventions that are currently in use (such as those aimed at disrupting VEGF) and, more interestingly, suggesting some new approaches that are counterintuitive yet potentially effective (such as those targeting the ECM to decrease the probability of the growing vessels forming loops).

CONCLUDING REMARKS
Developmental biology classically aims to understand how gene regulation leads to the development and morphogenesis of multicellular organisms. In this respect it has to be observed that genetic information influences the morphology and physiology of multicellular systems only indirectly. In fact, the gene network steers the biophysical properties of the cell by tuning the production of regulatory RNA sequences and proteins which in turn determine the behavior of the cell and its response to signals from its environment. In many aspects of biological development; therefore, what really matters are just these properties at the cell level (the type of signals released and the response to extracellular stimuli), and the cell's inner workings can be neglected. Thus, as proposed by some authors (see [30] for instance), two separate questions can be considered: the first one concerns how genetics drives cell behavior and the second one how cell behavior drives morphogenesis.
As far as the second question is concerned, the cell-centered modeling approach is certainly a significant tool to generate and test hypotheses in developmental biology, helping to understand which cell behaviors are essential to structure tissues. The studies on the angiogenic process represent a significant example of this concept. In fact, cell-oriented computational screening of the parameter dependence of patterning significantly helped to identify regulators of vascular development and suggest new hypotheses. Consistent with continuum models [20,22], the cell-centered approach confirmed the key role of chemotaxis in driving vascular formation both in vitro and in vivo [62,65]. Furthermore, it suggested that EC adhesion is essential to form stable vascular networks and that cell extension also strongly affects the multicellular patterns [57,59]. The understanding of the role of EC proliferation and of their interaction with the ECM on the formation of capillary sprouts in vivo was also greatly enhanced by cell-centered modeling efforts [65,66]. In this respect, continuous models have great difficulty assessing the role of these parameters. Thus, cell-oriented strategy to modeling angiogenesis represented a better tool to direct specific experiments, and recently a number of experimental validations of proposed models have been obtained [57,59]. As demonstrated by some of the studies here reviewed, this modeling approach can also assist in the identification of cell properties representing potential targets to improve tissue engineering constructs [68] or therapies [66,69] against angiogenesis-dependent pathologies.
To assist developmental biologists in the investigation of the question concerning which molecular processes are responsible for the essential cell behaviors leading to specific tissue-level or organism-level phenotypes, the cell-centered approximation will require extensions to both larger-and smaller-length scales [70]. The integration of models of the intracellular activity with cell-centered models of cell behavior seems to be possible in two fashions. The simpler strategy is likely to use microscopic models to provide estimates of the parameters controlling the cell-centered model. Alternatively, true hybrid models could be devised in which simulations of the inner cell processes function as components within cell-centered models (as in [71] and in the example provided by Scianna et al. [67]). Similarly, the cell-centered models can be interfaced with macroscale models of tissue or organ behavior either by providing estimates of tissue properties for continuum models or interacting directly with them in a hybrid model (an example is provided in [72]). As pointed out by Merks and Glazier [30], in this effort the key advantage of starting from a mesoscopic standpoint, such as a cell-oriented approach, is that we often need to introduce only a minimal additional algorithmic complexity and computation time to achieve results consistent with existing experimental data. Thus, this modeling strategy could also be a convenient tool to devise more complex models aimed to reach a better insight on the links between different levels of biological organization. Such a research effort may represent an important target for future work in the field of modeling angiogenic processes. In fact, it has to be observed that the development of multiscale models with capabilities to integrate processes spanning different spatial and temporal scales appears of particular relevance to complement experimental studies and address important questions concerning the vascular system, whose developmental and remodeling aspects seem intimately characterized by a complex multiscale logic [73].