## Abstract

Angiogenesis is crucial during many physiological processes and is influenced by various biochemical and biomechanical factors. Models have proved useful in understanding the mechanisms of angiogenesis and also the characteristics of the capillaries formed as part of the process. We have developed a three-dimensional hybrid, agent–field model where individual cells are modelled as sprout-forming agents in a matrix field. Cell independence, cell–cell communication and stochastic cell response are integral parts of the model. The model simulations incorporate probabilities of an individual cell to transition into one of four stages—quiescence, proliferation, migration and apoptosis. We demonstrate that several features, such as continuous sprouts, cell clustering and branching, that are observed in microfluidic experiments conducted under controlled conditions using few angiogenic factors can be reproduced by this model. We also identify the transition probabilities that result in specific sprout characteristics such as long continuous sprouts and specific branching patterns. Thus, this model can be used to cluster sprout morphology as a function of various influencing factors.

## 1. Introduction

Angiogenesis is the formation of new blood vessels from a monolayer of cells or by the reorganization of capillaries via morphogenesis. It is crucial during development, and for physiological processes like menstrual cycles and reproduction. Two other cases are of particular interest. First, angiogenesis occurs during invasive tumour growth because of the additional nourishment required by the tumour. Second, it is essential for tissue engineering purposes because it is necessary to be able to predict and control capillary development in scaffolds during *in vitro* tissue development. The primary stages of angiogenesis can be categorized as (i) endothelial cell (EC) activation by chemotactic agents, (ii) secretion of proteases into the matrix, (iii) selection of ECs for sprouting, (iv) capillary growth, (v) tubulogenesis (vi) non-angiogenic growth of blood vessels via capillary loop formation, and (vii) formation of second-generation capillaries (Adams & Alitalo 2007).

The process of angiogenesis is regulated by a balance between several pro-angiogenic and anti-angiogenic factors. Over the years, numerous factors—VEGF, PLGF, PDGF, TNF-α, TGF-β, α-FGF, β-FGF, ENA/VASP (Folkman & Shing 1992), angiopoietin-1 (Ang-1), angiopoietin-2 (Ang-2) (Davis *et al.* 1996) and chemokines like PF4 (Slungaard 2005)—have been shown to influence the process. Additionally, it has also been demonstrated that mechanical factors such as flow, extra-cellular matrix (ECM) stiffness and surface shear stress each affects the extent and directionality of capillary formation. Helm *et al.* (2005) showed that interstitial flow of the order of 1 μm s^{−1} in combination with VEGF induced directionality in capillary structures (in the direction of flow) and caused fibrin-bound VEGF to be released via proteolysis. Yamamura *et al.* (2007) studied the effect of substrate stiffness on bovine pulmonary microvascular endothelial cells (BPMECs) and demonstrated that BPMECs cultured on flexible collagen gels form networks in 3 days and show extensive sprouting and formation of complex networks in 5 days, whereas the cells cultured on stiffer gels tend to form aggregates and thicker networks.

Angiogenesis modelling is a useful tool for understanding the interplay between all the factors that affect it and for the design of experiments of a predictive nature. Over the years, various models spanning different scales and focusing on different aspects of angiogenesis have been developed. These can be classified as continuum models (Dallon & Othmer 1997; Anderson & Chaplain 1998; Levine *et al.* 2000; Chaplain & Anderson 2004; Chaplain *et al.* 2006; Chaturvedi *et al.* 2005) and discrete models (Stokes & Lauffenburger 1991; Anderson & Chaplain 1998; Mantzaris *et al.* 2004; Chaplain *et al.* 2006).

The continuum models are based on conservation equations for chemotactic and haptotactic gradients. In one such continuum model, the Chaplain group models a ‘tissue response unit’, which includes an EC, tumour angiogenic factor and a generic matrix molecule. The numerical solution is obtained from a finite-difference approximation subject to no flux boundary conditions and a specified initial condition. They have also developed a discretized version of the continuum model where the motion of the capillary in response to a tumour is governed by the EC at the tip (Anderson & Chaplain 1998).

The cellular Potts model (CPM) is a lattice-based model developed to describe the behaviour of cellular structures and their interactions. It could be an agent-based model, which is a computational model that is based on one (or more) specific component(s) and its effect on the individual cells (agents) being modelled. A two-dimensional agent-based model of angiogenesis based on CPM has been developed by Peirce *et al.* (2004), where they identify multiple cell types and growth factors. Their cell-level rule-based model of network growth in mesenteric tissue predicts new vessel formation, vessel length extensions and recruitment of contractile perivascular cells in response to localized pressure, circumferential strain and focal application of growth factors. The Sherratt group has used an extension of the Potts model to simulate malignant cells and quantify invasion morphology as a function of cell–cell adhesion (Turner & Sherratt 2002). In a different approach, the Popel group has developed a multi-scale integrative model with specific modules for various growth factor receptor pairs and ECM proteolysis (Qutub *et al.* 2009). Their model considers oxygen delivery by haemoglobin-based oxygen carriers (Tsoukias *et al.* 2007), the cellular response to oxygen in skeletal muscles (Ji *et al.* 2005) and a cell-based model that results in angiogenesis via reorganization of existing capillaries (Qutub & Popel 2009). Other models include a random walk model (Plank & Sleeman 2004*a*,*b*) which is distinguished by the fact that it places no restrictions on the direction of capillary growth, an individual cell-based two-dimensional mathematical model of tumour angiogenesis in response to a diffusible angiogenic factor (Plank & Sleeman 2004*a*,*b*) and a fractal-based model in which the smaller pieces of the system show ‘statistical self-similarity’ to the whole and the anatomical entities are given a fractal dimension. Random walk models that incorporate biochemical cascades when VEGF binding occurs have also been developed (Levin *et al.* 2002). Physiological models like one of corneal angiogenesis have also been developed. Jackson & Zheng (2010) have developed one such model that integrates a mechanical model of elongation with a biochemical model of cell phenotype variation. Despite the wide variety in the modelling approach, most of the models focus on tumour angiogenesis instead of *in vitro* capillary growth, and may lack one of the following: stochasticity, a three-dimensional framework or simplified binding kinetics, and are therefore difficult to apply in practice for tissue engineering applications. A combination of these characteristics in a model used for tissue engineering applications will be very useful. Finally, model validation in many of the existing models is a challenge owing to the difficulty in controlling all important factors *in vivo* combined with limited capability of most *in vitro* systems to replicate angiogenesis.

The complex biological processes leading to capillary morphogenesis are a consequence of cell-level decisions that are based on global broadcast signals, limited near-neighbour communication and stochastic decision-making with feedback control. Integrating these factors, a cell becomes programmed to follow one of several state trajectories that could be characterized as quiescence, division, apoptosis or migration. We have developed a model to address the needs for greater understanding of the process and for a practical tool with predictive capabilities.

This is a three-dimensional coarse-grained multi-scale hybrid model where each cell is modelled as an individual decision-making entity and cell–cell interactions are incorporated via the combined effect of cells on the matrix and the effect of the surrounding matrix in the individual cell decision-making process. Thus, this model demonstrates the phenomenon that when individual cells are modelled independently according to a set of rules and when cell–cell communication is embedded, the cell ensemble results in capillaries with features that can be attained experimentally in bioreactors for controlled tissue engineering purposes. It also provides a platform for bracketing these cell-ensemble results into clusters with different sprout characteristics and identifying the factors that affect them. Most importantly, this paper presents a model framework designed alongside experimental constraints and one that can simulate capillaries like the ones generated in an *in vitro* microfluidic system. This model–experiment cross-talk is crucial in identifying the effect of individual influencing factors on angiogenesis.

A major driving force in this model is its usefulness in predicting the angiogenic response in a closely regulated experimental platform with the objective of providing validation and of implementing feedback control over a prototypical biological process. One such experimental platform is the microfluidic system recently developed by our group (Vickerman *et al.* 2008). The first step in analysing model predictions is to identify a set of relevant sprout characteristics that would be useful in measuring the effect of different conditions and could be easily measured in real time during an experiment. These features are then compared with the simulation results.

## 2. Methods

This section outlines the basic steps in model development and the experimental system used for validation.

### (a) Experimental system

A two-channel microfluidic device (figure 1) fabricated in PDMS was used for all experiments. Human microvascular endothelial cells (Lonza, CC-0207, Batch no. 0000097428) were seeded as a monolayer on collagen I, which acted as the matrix. The device has four ports for media, two for each channel, and two gel filling ports. It is similar to one described in Vickerman *et al.* (2008). Cells were stained with 0.01 μM Hoestch nuclear dye and 20 μM 5-chloromethylfluorescein diacetate (CMFDA; CellTracker Green CMFDA; Molecular Probes, Invitrogen, C7025) cytoplasmic dye. Phase-contrast and confocal microscopies were used to monitor capillary development under different treatment conditions.

### (b) Model development

The model is discrete and can be divided into two distinct parts or modules: one is stochastic and the other is deterministic. These modules communicate with each other to predict capillary formation as a function of the local microenvironment including the effects of growth factors and matrix properties. This interplay is depicted in figure 2. This is a three-dimensional lattice model based on Markov processes, where at any time point each lattice point can be occupied by a cell, the matrix or remain empty. A cell can be in one of three states: quiescent, migrating or proliferating. It can also undergo apoptosis, at which point the lattice point becomes empty. Every lattice point that is occupied by the matrix is referred to as the ‘field’ and has associated with it concentrations of the various growth factors and matrix properties. In its current form, the model only includes VEGF, Ang-1 and MMP concentrations and matrix stiffness. However, this can be easily extended to include any number of growth factors and is modular in nature. Capillary features successfully depicted in this three-dimensional model are continuous sprouts, secondary and tertiary branches and cell clusters—all of which are observed in experiments. Three-dimensional modelling enables capillary surface area and protrusion volume calculations in the simulations and their comparison to experiments conducted in microfluidic devices. It also enables the proper representation of the factors in the surrounding matrix and their effects on individual cell response.

The deterministic component of the model includes the diffusion–convection–reaction equations that govern these growth factor concentrations. All cell–cell communication occurs via the field, i.e. the surrounding matrix. Direct cell–cell interactions were avoided to support the independence of each cell. Otherwise, the model would have to deal with joint probabilities of multiple cells, which are quite intricate. On the other hand, by making this assumption, we are still ensuring that the essence of cell–cell interactions is included while minimizing model complexity. This is done by a field variable termed the ‘connectivity factor’. Anything the cell ‘gives’ or ‘receives’ from a nearby cell is done via this field-associated factor.

A cell can exist in four main states: migrating (*M*), proliferating (*P*), quiescent (*Q*) and apoptotic (*A*). Its transition between states is depicted in figure 3. Once it undergoes transition to either a migratory or a proliferating state, it subsequently passes through a series of sub-states. The number of these sub-states is determined by the time of persistence of migration and the average cell cycle of an individual cell, thereby enabling us to account for differences in the time scales of these processes. Consequently, if a cell can exist in ‘*m*’ migratory sub-states and ‘*p*’ dividing sub-states, a transition probability matrix is established. Transition probabilities, e.g. *p*_{Q→Q}, *p*_{Q→M(1)}, *p*_{Q→P(1)}, *p*_{Q→A}, where *M*(1) and *P*(1) are the first steps of the respective processes, depend on the global (*g*) and local (*l*) conditions such as growth factor concentration and matrix properties in a manner described below.

The transition of each cell from one state at time *t* to another at time *t*+Δ*t* is stochastic and is a function of its current state, the condition of the surrounding matrix and external governing factors such as the presence or absence of flow, the growth factors present in neighbouring matrix and initial matrix stiffness. That is, transition from state *X*^{t} to *X*^{t+Δt} is described as *X*^{t+Δt}=*f*(*X*^{t},*U*^{t}(*g*),*U*^{t}(*l*)), where *U*^{t}(*g*) are the global or external factors and *U*^{t}(*l*) are the local factors or the characteristic of the surrounding matrix and measured by the output *Y* ^{t}. Certain variables like growth factor concentrations only affect the local environment while other factors like flow and pressure gradients have a global effect. The difference being, in the case of global influence, the same function is applied at all spatial locations, while in the case of local variables, the changes are calculated at each lattice point. This function is represented in the model as a transition probability.

#### (i) Field equations

The field equations in this model are presented below and can be simplified when required to provide feedback control. The field equations are written in this instance to include one growth factor, VEGF, and one protease, MMP. Similar equations for other factors can be easily incorporated. *C*_{vegf_s}, *C*_{vegf_b}, *C*_{vegf_r}, *C*_{MMP}, *M* and *M*_{cl} are the concentrations of soluble VEGF, VEGF bound to matrix, VEGF bound to receptor, MMP, matrix-binding sites available for binding to MMP and cleaved matrix, respectively.

The following governing equations describe *C*_{vegf_s}, *C*_{vegf_b}, *C*_{MMP} and *M*, where *D*_{vegf} is the diffusion coefficient of soluble VEGF, *D*_{MMP} is the diffusion coefficient of MMP and *v* is the interstitial flow velocity causing convective transport of VEGF. *k*_{on_m} and *k*_{off_m} are the binding constants for the reaction between soluble VEGF and binding motifs in the matrix. *k*_{on_c} and *k*_{off_c} are the binding constants for the reaction between soluble VEGF and VEGF receptors on the ECs. *ϕ* is the density of cells and *C*_{vegf_rec} is the number of VEGF receptors per unit volume, equal to the number of receptors on a single cell times *ϕ*. Pr_{vegf} and Pr_{MMP} are the rate of production of VEGF and MMP by a single cell, respectively. *ϕ* will be set equal to 1 for those grid points containing a cell (since we assume here a density of one cell per unit volume), and *ϕ*=0 for matrix. *k*_{c_m} is the rate at which the the matrix is cleaved by MMP.

The binding kinetics for the VEGF, MMP and matrix-binding sites (*M*) are:
(2.1)
rate constants *k*_{on_c} and *k*_{off_c}
(2.2)
production rate Pr_{MMP}
(2.3)
rate constant *k*_{c_m}.

The governing conservation equations are those for soluble VEGF: (2.4) where for bound VEGF: (2.5) where for MMP: (2.6) where and for the concentration of matrix-binding sites: (2.7) where Assuming local chemical equilibrium and given the fact that VEGF (VEGF-165 Isoform) binding to the matrix (collagen I) at pH 7.4 is negligible (Goerges & Nugent 2004), the equations simplify further. Also, there is no convective transport owing to interstitial flow, but this can be added to the model easily.

The first two equations are simplified as follows: (2.8) where and (2.9) These equations are discretized and used in instances where all the parameters are known. While most have been experimentally derived, some of the parameter values would still be subject to considerable uncertainty. In those cases, an approximation was made.

The diffusion term in equations (2.8) is discretized spatially as follows:
(2.10)
where Δ*x*, Δ*y* and Δ*z* are the spacings of the computational grid in *x*, *y* and *z*, respectively, and are equal to 10 μm (figure 4). The diffusion term in equation (2.6) is discretized similarly.

For the reaction terms, the following methodology is used. The concentration of the biomolecules being modelled will locally change in a single time step by an amount that is dictated by local conditions. This change in the concentration is a sum of the local consumption and production of that molecule and is expressed as
(2.11)
and where
and
These could potentially be experimentally determined by measuring changes in the concentration of signalling molecules, MMP and matrix around a single cell. This can be done by using fluorescent collagen and measuring the matrix properties around individual cells. However, as such experimental data for individual cells are not available, these Δ*C* values can be empirically established relative to the rates of production and degradation of these substances by quantifying the qualitative effects on individual cells. When those are being ascertained, Δ*C*_{i,C,j} and Δ*C*_{i,P,j} represent step functions associated with this change at the *j*th level in the concentration of the *i*th biomolecule included in the model owing to consumption and production by the cells, respectively. In the simulation results included in this paper, we have used a two-level change for the biomolecules being considered, i.e. high production/consumption or low production/consumption. These changes depend on the cell state and the surrounding matrix characteristics. One can imagine that instead of having only a high and low value (e.g. where *j*=1 or 2) for each of these reactions, we can have finer increments where each element is dependent on the different permutations and combinations of other factors in the field. The existing model framework allows for such expansion depending on the different pro- and anti-angiogenic factors being considered.

The equations are discretized over time and the spatial derivatives are calculated at time *t*.

Equation (2.8) is discretized as
(2.12)
where
and where
and
At a time *t*, the amount of VEGF produced and consumed via binding to the cells is a function of the cell state and the state of the surrounding matrix concentration of soluble VEGF, i.e. *C*_{vegf_s} in the lattice. Thus, the reaction terms are evaluated based on the cell state. *k*_{on_c} and *k*_{off_c} values were obtained from the literature (Mac Gabhann & Popel 2007). Hence, the amount of VEGF consumed (Δ*C*_{vegf_s,C,i}) is determined by the rate equation. However, as it is difficult to obtain the amount of VEGF produced by an individual cell, it was approximated by the following method. If there is a migrating cell at a given lattice point, we assume it produces a certain amount of VEGF (Δ*C*_{vegf_s,P,1}) if the surrounding VEGF concentration is below a certain threshold and a different amount (Δ*C*_{vegf_s,P,2}) when it is above a threshold. This threshold of the surrounding VEGF is 12 units (or 24 ng ml^{−1} normalized to the minimum concentration of VEGF, i.e. 2 ng ml^{−1}). This was chosen because it has been shown that VEGF induces its own expression in microvascular ECs in a STAT3-dependent fashion when the cells are treated with 25 ng ml^{−1} VEGF (Bartoli *et al.* 2003).

Assuming negligible interstitial flow, equation (2.6) is discretized as
(2.13)
The change in MMP concentration owing to production in equation (2.13) was approximated by assuming that the change around a migrating cell (Δ*C*_{MMP,P,1}) and that around a quiescent or proliferating cell (Δ*C*_{MMP,P,2}) are different.

Equation (2.7) is discretized as
(2.14)
The value of *k*_{c_m} is obtained from the literature (Karagiannis & Popel 2004). The amount of matrix consumed (Δ*C*_{M,C,i}) is determined by the rate equation and the state of the cell as a migrating cell tends to release more MMP, cleaving the matrix.

#### (ii) Cell transitions

The transition probability of a cell from one state to another is dependent on the initial conditions, such as growth factor concentrations in the medium and matrix stiffness. Cell decisions are made once every hour or every ‘cellular’ time step. Thus, in a period of every hour, a certain number of cells migrate, proliferate, become quiescent and die. These probabilities are derived empirically. We have explored the range of probabilities and identified combinations that depict sprout characteristics seen under certain experimental conditions. These combinations are shown in figure 5. The transitional probability to a particular state at time *t*+1 can be written as
(2.15)
(2.16)
(2.17)
and
(2.18)
where *p*_{M,t+1}, *p*_{Q,t+1}, *p*_{P,t+1} and *p*_{A,t+1} are the probabilities that the cell state at time *t*+1 becomes *M*, *Q*, *P* or *A*, respectively.

Parameters *p*_{Q,t},*p*_{Mm,t} and *p*_{Pp,t} are the probabilities defined as
(2.19)
(2.20)
and
(2.21)
and *p*_{Q→M,t}, *p*_{Q→Q,t}, *p*_{Q→P,t} and *p*_{Q→A,t} are the transition probabilities to migration, quiescence, proliferation and apoptosis, respectively, under the instantaneous global and local conditions. These are expressed as a probability distribution of percentages of transition (equations (2.22)–(2.25)). These probabilities are to be determined experimentally. Alternatively, they could be specified functions of the activation of relevant intracellular signalling pathways.
(2.22)
(2.23)
(2.24)
and
(2.25)
For the current simulation, these transitional probabilities are arbitrarily varied to depict the presence of VEGF, Ang-1 and differences in matrix stiffness. Certain ‘rules’ for cell transition are outlined below.

— Initially, a uniform monolayer of cells is assumed to exist on the surface of the matrix. New sprouts are allowed to be initiated for the first 4 h. This is implemented to prevent multiple sprout formation and depicts the biological state when only certain cells in the monolayer sprout and act as ‘inhibitors’ on nearby cells. This lateral inhibition is enforced via the field.

— Each cell undergoes a ‘decision-making’ process where it can decide to migrate, divide, die or stay quiescent.

— If a cell has ‘decided’ to migrate, the direction of migration is stochastic, though it is biased towards the lattice position occupied by the matrix that is associated with the highest concentration of chemoattractants and MMP by giving the cell a higher probability to migrate towards that lattice point. It has a lower probability to migrate into any of the other matrix-occupied lattice positions.

— If a cell has decided to divide, the new cell occupies the matrix-occupied lattice position that is associated with the highest concentration of MMP as it causes local degradation of the matrix.

— If a cell dies, its position is occupied by empty space.

— A cell can migrate only into a lattice position occupied by matrix but can divide into either empty space or matrix.

— If two cells choose the same lattice point into which to divide or migrate, a tie-breaking rule is applied.

— A migrating or dividing cell releases MMP into each of the adjacent matrix elements, thus influencing the migration of itself and neighbouring cells. The amount released depends on the state of the cell and can be classified as ‘low’ or ‘high’—the values of which are recorded in table 1.

— VEGF is both released and consumed by the cells according to the state of the cell and the surrounding VEGF concentration. As explained earlier and recorded in table 1, the surrounding local concentration above which VEGF production is high is 12 and that above which VEGF consumption by a cell is high is 10 units.

— Cells are allowed to divide both at the monolayer and in the stalks.

— In the current simulations, the number of migration sub-states (

*m*) is 3 and the number of proliferation sub-states (*p*) is 20, where each ‘cellular’ time step corresponds to approximately 1 h.— VEGF influences endothelial decision making via both paracrine and autocrine signalling. So, the deterministic model accounts for the diffusion, consumption and production of VEGF and MMP by the cells.

— While the sample simulations included in this paper account for change in concentration of VEGF and MMP, and initial addition of VEGF (40 or 20 ng ml

^{−1}) and/or Ang-1 (500 ng ml^{−1}) to the field, the module-based algorithm and signalling molecules interacting with the cells in a modular array format ensure that additional molecules can be included easily.— In the simulations, the three-dimensional lattice is normalized to the characteristic length of an EC of approximately 10 μm. The concentrations of the different growth factors and MMP are normalized to the maximum concentration in every simulation.

— All simulations are recorded 50 h after the monolayer is established. All experimental images are documented 2 days after cell seeding.

## 3. Results

### (a) Sprout morphology

In all the simulations, cells were seeded as a monolayer and a VEGF gradient was established from *z*=20 before the simulations were started. Simulations were run with different transition probabilities for *p*_{Q→Q}, *p*_{Q→M(1)}, *p*_{Q→P(1)} and *p*_{Q→A} transitions.

Figure 5*a*–*e* is generated for five such sets of transition probabilities and includes a pictorial representation of sprout formation under different conditions and phase-contrast images of experimental conditions that give rise to similar sprout formation characteristics. The experiments were conducted in the microfluidic devices referred to above, with or without two important growth factors: VEGF and Ang-1. Various experiments conducted in the microfluidic device in our laboratory have demonstrated that capillary stabilization is achieved when Ang-1 is used in the media along with VEGF. Since the experiments in our laboratory have been limited to VEGF and Ang-1, it was decided to incorporate these two pro-angiogenic factors in the first stage of the modelling. It has also been shown that Ang-1 has an impact on VEGF signalling irrespective of the presence or absence of Ang-2 (Zhu *et al.* 2002). Ang-2 is a natural antagonist and disrupts blood vessel formation. As the primary goal of the paper is to determine conditions that will result in the evolution of well-developed capillaries, it was decided to first focus on such factors. In our experience, high transition probability of migration and proliferation corresponds to high concentrations of growth factors and matrix of optimal stiffness, while low probabilities correspond to a lower concentration of growth factors and a softer or stiffer matrix. High concentration of Ang-1 implies a lower rate of regression owing to capillary stabilization. These figures are included here to demonstrate that certain experimental conditions can be reproduced by some combination of transition probabilities. The exact matching of these probabilities to specific conditions will require further experimentation.

When the transition probabilities represent extreme states such as 85 per cent migration/5 per cent proliferation/5 per cent quiescence/5 per cent proliferation, many migratory cells are seen that protrude far into the matrix. Under a different extreme state, e.g. 85 per cent proliferation/5 per cent migration/5 per cent quiescence/5 per cent proliferation, the entire monolayer starts to proliferate and forms large cell clusters. These extreme cases are not representative of any physiological condition, but were conducted to demonstrate the range of applicability of the model. As observed in figure 5, longer sprouts are observed with a higher migration probability and almost no sprouts are observed when the cells have a higher probability to remain quiescent.

Sprout formation under conditions of VEGF (20 ng ml^{−1}) and high Ang-1 (500 ng ml^{−1}), which act as sprout-stimulating and sprout-stabilizing agents, respectively, in a collagen I gel made with concentration 2 mg ml^{−1} and simulation generated with transition probabilities that result in a similar morphology to the above mentioned condition are shown in figure 5*a*. The branching network is well developed under such conditions. Simply by removing Ang-1 and increasing VEGF concentration, so that there is a high concentration of the sprout-stimulating agent (VEGF 40 ng ml^{−1}) but no sprout-stabilizing agent, cell clusters and short branches are prevalent (figure 5*b*). With only Ang-1 (500 ng ml^{−1}) present and no VEGF, the absence of a sprout-stimulating agent leads to relatively less sprout formation. Very few stable capillaries are observed relative to small sprouts and single migrating cells (figure 5*c*). With only VEGF (40 ng ml^{−1}) and no Ang-1 in a stiffer collagen I gel (2.5 mg ml^{−1}), there is a sprout-stimulating agent but no sprout-stabilizing agent. But owing to the stiffer gels, the cells find it harder to migrate; hence they cannot form many branches (figure 5*d*). When both VEGF (20 ng ml^{−1}) and Ang-1 (500 ng ml^{−1}) are present in a stiffer gel (2.5 mg ml^{−1} collagen I), some capillary-like structures are observed, but fewer and longer branches and more cell clusters are observed at the monolayer than would be seen in a less stiff matrix under the same conditions (figure 5*e*). In this case, the branches have fewer cells per branch compared with the same growth factor combination in a softer matrix. These sample results demonstrate that by proper adjustments of the parameters, the model can produce a variety of morphologies. In addition, the model can also be used to explore how specific features depend on the imposed conditions, as in the following examples.

Consider first the average protrusion length of cells into the matrix (figure 6). This is determined by calculating the depth to which cells penetrate at each lattice point in the monolayer in the *z* direction and assuming unit volume occupied by each cell. Corresponding sets of transitional probabilities are listed in the figure. In the first three transitional sets, the probabilities to migrate and die are held constant, the probability to divide decreases and that to remain quiescent increases. In the next three transitional probability sets, the probabilities to divide and die are held constant, the probability to migrate decreases and that to remain quiescent increases. As expected, certain transition probabilities result in higher cell protrusion into the matrix. As the probability of the cell to remain in the quiescent state increases in the first dataset, the protrusion distance increases because of the lower proliferation rate. Thus, fewer sprouts are initiated and the ones that are, grow long. This implies that while proliferation occurs both in the monolayer and in the sprouts, addition of any proliferating agent affects the proliferation rate in the monolayer more significantly. Hence, the addition of growth factors that simulate proliferation should result in a larger network of sprouts and/or more cell clusters. In the second data cluster, with a higher migration probability, longer sprouts are formed. While the difference in sprout length is not significant with small changes in migration probability, the trends observed are very significant. It provides an insight into the conditions that would give rise to network-like versus single-sprout-like capillary structures.

### (b) Continuous sprouts versus cell clustering

In order to understand the effect of each of the different transition probabilities, the number of different morphologies are recorded when two of the four transitions are kept constant and the other two are systematically varied. Figure 7 shows an average of 30 simulations for each of these conditions.

To determine the effect of increasing specific transitional probabilities, we varied specific ratios (migration (*M*) to quiescent (*Q*); proliferation (*P*) to quiescence (*Q*); and migration (*M*) to proliferation (*P*)) systematically by keeping the other two constant in each case. The results observed for each of the above-mentioned variations are depicted in figures 7 and 8. As the *M*/*Q* ratio increases (figure 7*a*), the number of individual migrating cells increases and the number of cell clusters at the monolayer decreases. This is because at a certain fixed transition probability of proliferation, as the tendency of cells to migrate increases, cluster formation decreases. However, discontinuous sprouts or clusters formed away from the monolayer show a bimodal response. At the first peak, these clusters tend to be broader, and at the second peak (when the tendency to migrate is higher than that to proliferate), these clusters tend to be longer. When the ratio of migration to quiescence nears one, the number of such ‘floating clusters’ decreases. This is also the condition that causes an increase in the number of continuous cell sprouts. This shows that these two morphologies lie in two distinct areas in the space map of the different transition probabilities.

A crucial element in capillary morphology prediction would be in identifying the balance between the transition probabilities of migration and proliferation. Figure 7*b* shows the effect of varying these two transition probabilities keeping the transition to quiescence and transition to apoptosis fixed at 15 per cent and 5 per cent, respectively. The most relevant prediction is that the number of continuous probabilities peak when this ratio tends to one. Though the number of cell clusters in the matrix does not change much, these clusters, like the clusters at the monolayer, tend to become larger as this ratio increases.

The simulations in figure 7*c* were generated while keeping the transition probability of migration fixed at 15 per cent and the transition probability of apoptosis fixed at 5 per cent. While the number of continuous sprouts and discontinuous sprouts or cell clusters in the matrix remains unchanged as the ratio increases, it is also important to note that at this fixed transition rate for migration, very few continuous sprouts are formed. Since this lies in the first half of the curve that determines the effect of migration on sprout formation as shown in figure 7*a*, it is verified that sprout formation is optimal at a higher migration transition probability. The number of single migrating cells decreases and the number of clusters at the monolayer increases as this ratio increases.

The customary capillary morphology found *in vivo* is continuous sprouts, and its ratio to each of the other morphologies as a function of the ratio of the transition probabilities between migration and quiescence, proliferation and migration and proliferation and quiescence are shown in figure 8. As expected, the curves peak when there is a balance between migration and quiescence, which translates into effective capillary induction and stabilization (figure 8*a*). The model predicts that the ratio between proliferation and quiescence should be increased in order to increase the number of continuous sprouts in comparison to individual migrating cells (figure 8*c*). Again, as expected after observing figure 7*b*, the number of continuous sprouts is maximized in comparison to any other capillary morphology when the ratio of transitioning to proliferation to that of transitioning to migration tends towards one (figure 8*b*).

### (c) Branching

This model is also a useful tool for studying and analysing branching patterns in capillary formation. Very little is known about the characteristics of capillary networks and the factors that influence them. In order to use this model to study such patterns, we first need to identify certain significant characteristics. We limited these to number of sprouts, number of total branches, number of branches per sprout and number of anastomoses per simulation. These were analysed as a function of changing transition probabilities based on the expectation that most growth factors will have an impact on migration, proliferation, apoptosis and quiescence and hence will fall in the transition data space being analysed. The results shown in figure 9 can be used to narrow down the transition probabilities that provide the desired network characteristics.

### (d) Predicting sprout characteristics with and without VEGF

A very simple prediction is made using the described model, and experiments were conducted in the microfluidic devices to demonstrate the model–experiment overlap. The conditions used to demonstrate this are the presence and absence of VEGF—the most crucial pro-angiogenic factor.

The concentrations used are no VEGF and 40 ng ml^{−1} VEGF in a collagen I matrix. Absence of VEGF implies a set of conditions that increases the bias for the cells to remain quiescent and decreases their bias to migrate or proliferate. Hence, they fall in the following bracket of transition probabilities: [0.85, 0.05, 0.05, 0.05]. Presence of 40 ng ml^{−1} VEGF implies that there is comparatively a higher bias for cell migration and proliferation. Thus, they fall in the following bracket of transition probabilities: [0.7, 0.15, 0.1, 0.05]. The experiments were conducted in the microfluidic devices that are described earlier and stained with phalloidin (Invitrogen) after being fixed. The images were obtained using an Olympus confocal microscope and processed using ImageJ software (NIH, Bethesda, MD, USA). As recorded, most of the sprouts formed under no VEGF condition were located along a surface of the device and showed only two-dimensional invasion. The sprouts under the pro-angiogenic condition of 40 ng ml^{−1} VEGF showed three-dimensional invasion and well-developed branching. This was reflected in the model predictions as well.

Figure 10 summarizes the results obtained from experiments and simulations under the aforementioned conditions. The model results are an average of 30 simulations each and the experiments are an average of 16 experimental regions. The mean number of sprouts, number of branches, number of anastomoses and number of migrating cells predicted by the model are reflected in the experiments conducted in the devices.

## 4. Discussion

### (a) Intracellular signalling

An advantage of this model is that it allows for the incorporation of an intracellular signalling model. Though it has not been included in the sample simulations provided in this paper, the model can call intracellular modules during the different states: migration, division and quiescence. These functions will predict the signalling pathway that will be used by the different growth factors and the cellular responses to them in a probabilistic manner. The output of these functions (cellular responses) will act as the input for the next step of cell ensemble response. Currently, cell-signalling experiments under the same input conditions that provided the cell morphology responses are underway. The data-driven intracellular modelling will be done by applying either PLSR or decision-tree analysis to the obtained signalling data. Decision-tree analysis has been applied to many different areas of biomedical decision-making processes (Adam *et al.* 2002; Hua *et al.* 2006). This section of the model continues to be developed, and since the model is modular, it can be included without altering any other framework.

### (b) Controls framework

This entire model will be included in a controls framework that is currently being developed. The details of the procedure are outlined in Wood *et al.* (2009), where a less detailed capillary formation model is included in the controls framework presented. By applying this framework to the detailed stochastic angiogenesis model discussed in this paper, we will be able to predict capillary characteristics as a function of the transition probability space to an even greater detail. It will also enable us to study the temporal dynamics of the different transitions.

### (c) Limitations

The simplifications made in the algorithm of this model may make it difficult to identify factors that affect the sprout characteristics in small gradients. Additionally, the usage of the ‘field’ as the cell communication medium is primarily used for modulation of endothelial cell–cell contact, connectivity in the sprout and lateral inhibition owing to notch signalling. Since the cells cannot change shapes in the current format of the model, such communication restricts the ability to maintain the cell junction. This limitation of the model can be overcome by allowing the cells to change shape in future versions. Additionally, cell polarity and intercellular communication, which is a function of cell polarity, is ignored in the current form of the model. It can be foreseen that this can be incorporated when cells are allowed to change shape. In that event, the cell communication function will have to evolve to account for that complexity.

### (d) Advantages

The advantage of the three-dimensional agent-based coarse-grain stochastic angiogenesis model outlined here is that several angiogenic factors can be included and that it can function both as a stand-alone model and inside the control framework mentioned above. Most importantly, all simulations can be compared in-house with experiments conducted in the mircrofluidic devices described. By modelling cell populations via the application of global controls that are collectively modified by the aggregate cell population, cell–cell interactions are introduced into this model in an explicit manner. Several previous works have modelled capillary growth as a function of the leading cell or developed models that involve distinct zones of capillaries. For example, Balding & McElwain (1985) modelled two regions of EC density: the sprout proper and the tip; however, they neglected dynamics of the cells remaining on the monolayer. Such cell–cell interactions are crucial for developing a model that reflects the impact of cellular signalling molecules and the ECM because of the feedback loop that exists between the aggregate effect of cells on the matrix and the impact of the matrix on each individual cell. Since then, as mentioned in the introduction, several types of models that provide crucial insights into the process have been developed. However, the model we present is unique in its ability to cross-talk with experimental results, and thus provides a platform for predicting non-tumour angiogenesis in *in vitro* bioreactors.

By developing a model based on the influence of global and local conditions, the coarse-grain combined impact of various biochemical and biomechanical factors on angiogenesis can be determined. The model can be extended to include any number of biomolecules that are either being added to the bioreactor (affecting global conditions) or being secreted/consumed by migrating/dividing cells (affecting local conditions) because the ‘molecular matrix’, i.e. the components associated with each lattice in the field, can be of any size. Such a methodology is essential in order to develop a system-based, experimentally verifiable approach to the process of capillary formation.

This model enables us to simulate individual cell behaviour and compare cell population behaviour at the end of the simulation. Certain basic characteristics observed in these simulations that reflect experimental observations are the presence of continuous sprouts, their differential formation as a function of added growth factors, horizontal and vertical migration and cell clustering. Different transitional probabilities reflect these experimental observations, and future work will explore the predictive nature of such a model.

## 5. Conclusion

In conclusion, we demonstrate here the successful implementation of a three-dimensional agent–field model for capillary formation that can generate sprouts with characteristics similar to those observed in three-dimensional *in vitro* experiments. This is a hybrid stochastic–deterministic model where individual cells are modelled independently and their communication via cell–cell interactions is included in a mathematically simplified manner that captures the importance of that action but does not complicate the implementation. The simulation results obtained are characteristically similar to experimental observations as shown, indicating that this model can be used for various demonstrative and predictive applications.

## Acknowledgements

We would like to acknowledge NSF-EFRI grant no. 0735997 and the Singapore-MIT Alliance for Research and Technology for funding.

## Footnotes

One contribution of 13 to a Theme Issue ‘The virtual physiological human: computer simulation for integrative biomedicine II’.

- © 2010 The Royal Society