Skip to main content

Module-based multiscale simulation of angiogenesis in skeletal muscle



Mathematical modeling of angiogenesis has been gaining momentum as a means to shed new light on the biological complexity underlying blood vessel growth. A variety of computational models have been developed, each focusing on different aspects of the angiogenesis process and occurring at different biological scales, ranging from the molecular to the tissue levels. Integration of models at different scales is a challenging and currently unsolved problem.


We present an object-oriented module-based computational integration strategy to build a multiscale model of angiogenesis that links currently available models. As an example case, we use this approach to integrate modules representing microvascular blood flow, oxygen transport, vascular endothelial growth factor transport and endothelial cell behavior (sensing, migration and proliferation). Modeling methodologies in these modules include algebraic equations, partial differential equations and agent-based models with complex logical rules. We apply this integrated model to simulate exercise-induced angiogenesis in skeletal muscle. The simulation results compare capillary growth patterns between different exercise conditions for a single bout of exercise. Results demonstrate how the computational infrastructure can effectively integrate multiple modules by coordinating their connectivity and data exchange. Model parameterization offers simulation flexibility and a platform for performing sensitivity analysis.


This systems biology strategy can be applied to larger scale integration of computational models of angiogenesis in skeletal muscle, or other complex processes in other tissues under physiological and pathological conditions.


Angiogenesis is a complex process whereby new capillaries are formed from pre-existing microvasculature. It plays important roles in many physiological processes including embryonic development, wound healing and exercise-induced vascular adaptation. In such processes, robust control of capillary growth leads to new healthy pattern of physiological vessel network that matches the metabolic demands of development, wound repair, or exercise [1]. In contrast, excessive or insufficient growth of blood vessels is associated with an array of pathophysiological processes and diseases, among which are malignant tumor growth, peripheral artery disease, diabetic retinopathy, and rheumatoid arthritis [1].

Systems-level studies of angiogenesis in physiological and pathophysiological conditions improve our quantitative understanding of the process and hence aid in therapeutic design. Extensive experimental studies of angiogenesis over the past two decades have revealed that the angiogenesis process is comprised of a series of events at multiple biological organization levels from molecules to cells, tissues, and organs. For example, as a first approximation, exercise-induced angiogenesis can be described as a sequence of the following events: i) Exercise increases oxygen consumption in tissue, followed by increased blood flow in the vasculature, thus affecting convection-diffusion oxygen transport processes [2]; ii) As exercise continues, insufficient oxygen delivery to the tissue leads to tissue cellular hypoxia, which results in activation of the transcription factor hypoxia-inducible factor 1α (HIF1α) [3] and the transcription coactivator peroxisome-proliferator-activated-receptor-gamma coactivator 1α (PGC1α) [4]; iii) These factors induce the upregulation of vascular endothelial growth factor (VEGF) expression [5]. VEGF is secreted from myocytes (and possibly stromal cells), diffuses through the interstitial space, and binds to VEGF receptors (VEGFRs) on microvascular endothelial cells; concomitantly, endothelial cell expression of VEGFRs is also altered [6]; iv) The increase in VEGF and VEGFR concentration and possibly VEGF gradients results in activation of endothelial cells and cause capillary sprouting. Thus new capillaries and anastomoses form and new capillary network patterns develop [7]; v) After exercise, VEGF and VEGFR expression remain elevated for a limited time and thereafter return to basal levels [6]. The signaling set in motion causes blood vessel remodeling to continue after exercise. Thus, the time scales of individual events range from seconds in oxygen convection-diffusion processes to hours in VEGF reaction-diffusion processes, to days or weeks in capillary sprouting processes. Spatial scales vary from nanometers at the molecular level to microns at the cellular level, to millimetres or centimetres at the tissue level.

The complexity of angiogenesis is a function not only of the multiscale characteristics in temporal and spatial domains, but also of the combinatorial interactions between key biological components across organizational levels. At the molecular level, multiple HIF-associated molecules and hundreds of genes activated by HIF form a complex transcriptional regulatory network [8]. Six isoforms of VEGF-A (VEGF121,145,165,183,189,206), three VEGFRs (VEGFR-1, -2, -3) and two coreceptors (neuropilin-1 and -2) constitute a complex ligand-receptor interaction network, regulating intracellular signaling and determining cellular response [9]. In addition, other VEGF proteins, like placental growth factor (PlGF) and VEGF-B, -C, -D, compete with VEGF-A for some of the same receptors. Matrix metalloproteinases (MMPs) also form a key molecular family with approximately 30 members; MMPs are capable of proteolyzing components of the extracellular matrix (ECM) thus decreasing the physical barriers encountered by a tip endothelial cell leading a nascent capillary sprout [10]. At the cellular level, endothelial cell activation, migration and proliferation are driven by local growth factor concentrations and gradients. Capillary sprouting is also governed by the interaction between a tip cell and its following stalk cells, and by cell adhesion to the ECM [11]. In addition, parenchymal cells, precursor cells and stromal cells, as well as the ECM, constitute the nascent sprout microenvironment, influencing endothelial cell signaling, adhesion, proliferation and migration.

Mathematical and computational models of angiogenesis have become useful tools to represent this level of biological complexity and shed new light on key control mechanisms. In particular, computational modeling of tumor-induced angiogenesis has been an active area of research over the past two decades and has also been extensively reviewed [1215]. Here we give a brief overview of the angiogenesis models relevant to building a multiscale model of angiogenesis in skeletal muscle using different modeling methodologies. The models can be classified into continuous, discrete and hybrid categories. Continuum models of growth factor activity often applied molecular-detailed reaction and reaction-diffusion differential equations. These models have been used to describe many aspects of angiogenesis, e.g., host tissue distribution of a chemotactic factor following its secretion from a tumor [16], VEGF-VEGFR interactions [17], a fibroblast growth factor-binding network [18], whole-body compartmental distribution of VEGF under exercise and peripheral artery disease conditions [19, 20], the contribution of endothelial progenitor cells to circulation of VEGF in organs and their effects on tumor growth and angiogenesis [21], and a VEGF reaction-transport model in skeletal muscle [22]. Models of other angiogenesis-associated proteins such as MMP2 and MMP9 also have been developed [23, 24]. By describing capillary networks in terms of endothelial cell densities, continuum models have also been developed to represent tumor-induced capillary growth [2527] and the wound healing process [28]. Discrete models such as cellular automata [29], cellular-Potts model [30], and agent-based models [3133] have been developed to describe tissue behavior stemming from the interaction between cells, extracellular proteins and the microenvironment. These cell-based models offer unique capability of representing and interpreting blood vessel growth pattern as an emergent property of the interactions of many individual cells and their local microenvironment. By combining the continuum approach with the cell-based modeling approach, hybrid modeling can be used to describe the in vivo vascular structure along with detailed molecular distributions [3437], providing appropriate computational resolution across various scales.

With a number of computational models currently available to describe different aspects of angiogenesis, integration of existing models along with new biological information is a promising strategy to build a complex multiscale model [38, 39]. While current advances mainly focus on the representation format of molecular interaction models (e.g., XML-based notation) and dynamic integration of these models (e.g., Cytosolve [40]), few strategies exist to combine existing models at multiple scales with mixed methodologies. Here we describe our development of a novel computational infrastructure to coordinate and integrate modules of angiogenesis across various scales of biological organization and spatial resolution. Using this approach can significantly reduce model development time and avoid repetitive development efforts. These modules can be adapted from previously-developed mathematical models. Our laboratory has developed a number of angiogenesis models including: oxygen transport [41]; VEGF reaction-diffusion [22]; capillary sprouting [33]; FGF-FGFR ligand-receptor binding kinetics [18]; MMP proteolysis [23, 24, 42]; and MMP-mediated VEGF release from the ECM [43]. We show results of a test case, which integrated a blood flow model, an oxygen transport model, a VEGF transport model and a cell-based capillary sprouting model. With the use of Java Native Interface functions, previously-developed angiogenesis models were redesigned as "pluggable modules" and integrated into the angiogenesis modeling environment. Another advantage of this simulation infrastructure is its flexibility, allowing integration of models written using different simulation techniques and different programming languages. Note that the primary aim of this study is building methodology for multiscale modeling, rather than obtaining novel physiological results; detailed simulations of skeletal muscle angiogenesis and comparison to experimental data will be presented elsewhere.

The computational scheme presented here fits into the Physiome Project defined as a computational framework allowing the integration of models and databases that intends to enhance the descriptive, integrative and quantitative understanding of the functions of cells, tissues and organs in human body [4446]. Integral parts of the Physiome Project are the Cardiac Physiome [47], the Microcirculaton Physiome [48, 49], and the EuHeart project, which are aimed at specific organs or physiological systems. The Virtual Physiological Human project is also aimed at a quantitative description of the entire human [5052]. To achieve the goals of these projects, it is essential to share computational models between a variety of modeling methodologies, computational platforms, and computer languages and incorporate them into integrative models. One approach in the past decade is to develop XML-based markup languages to facilitate model representation and exchange. The two most-accepted formats, SBML [53] and CellML [54], are designed to describe biochemical reaction networks in compartmental systems expressed by ordinary differential equations (i.e., they have no spatial description). FieldML [55] allowing for spatial description is under development. Alternatively, the object-oriented modeling methodology provides a strategy to describe the biological organizations and flexible solution to integrate currently available models. For example, universal modeling language (UML) [56, 57] and other meta-languages such as E-cell [58] have been proposed. However, the robustness of the integration of external models is dependent on the interface of these meta-languages. In the current study, we propose to use a natural object-oriented language, Java, as a modeling language to design the integration controller and link currently available modules at different scales.

Systems and Methods

We describe a computational platform capable of linking any number of modules. In the particular example of skeletal muscle angiogenesis, we integrate four modules: microvascular blood flow; oxygen transport; VEGF ligand-receptor interactions and transport; and a cell module describing capillary sprout formation. These four modules use diverse modeling methodologies: algebraic equations (blood flow), partial differential equations (PDEs, oxygen and VEGF transport) and agent-based modeling (ABM, cell model). An overview of the simulation scheme is shown in Figure 1A. Briefly, the model initiates with the input of a three-dimensional (3D) muscle tissue geometry that includes muscle fibers and a microvascular network; rat extensor digitorum longus (EDL) muscle is used as a prototype, as previously described [22]. This tissue geometry is first used to calculate blood flow in the vascular network, and then in the computation of oxygen distribution in the vascular and extravascular space, followed by simulation of VEGF distribution in the interstitial space and on the endothelial surface, and finally simulation of capillary sprouting and remodeling of the vascular network. Blood flow and hematocrit are simulated using the two-phase continuum model proposed by Pries et al[59]. Oxygen transport model [60] is used to calculate the spatial distribution of oxygen tension throughout the tissue. VEGF secretion from myocytes through an oxygen-dependent pathway is described by an experiment-based oxygen-dependent transfer function dependent on the factors HIF1α and PGC1α (details are below). A modified VEGF reaction-diffusion model [22, 61] is used to predict the spatial VEGF distribution in tissue interstitial space and at the surface of the endothelium. Using our agent-based model with this VEGF concentration profile defined as input [33], we further compute elongation, proliferation and migration of endothelial cells forming capillary sprouts. The result is a new capillary network. In turn, this new structure feeds back into the integrated model as an updated vascular geometry, and starts a new cycle with the flow model, oxygen transport model and VEGF reaction-diffusion model, thus simulating the dynamics of the angiogenesis process. Governing equations and a brief description of each individual module are given below.

Figure 1
figure 1

Schematics of Module-based Mulitscale Angiogenesis Modeling Methodology. A) Skeletal muscle angiogenesis is modeled as a multi-step process. It starts with a blood flow simulation followed by a simulation of oxygen convection-transport process. Using O2 tissue distribution, VEGF secretion by myocytes is computed as a function of oxygen-dependent transcription factors HIF1α and PGC1α; then a VEGF reaction-transport process is computed. Lastly, capillary formation is simulated based on VEGF concentration and gradients. Feedback loops increase the complexity of the model since a new geometry with nascent vessels will affect blood flow conditions, tissue hypoxia, and VEGF secretion and distributions. All four processes are simulated using a variety of modeling techniques and languages. We use Java as the language for modeling the controller, and apply JNI plugins to link these modules together. The controller is composed of four sub-packages, including Process, Biosystems, IO and Exceptions. B) Communications between different modules and Java codes in core package are implemented by transferring each module into a shared object library (SO file in Linux). Upper panel shows that two wrapper files (includes Java-to-C and C-to-Fortran wrapper) are written to communicate between the flow Java class defined in the controller and the Fortran flow module, to call the flow module in Fortran. Lower panel shows that a JNI C wrapper is required to transfer the data between the modeling controller (in Java) and the Oxygen/VEGF module (in C/C++).

Skeletal Muscle Tissue Geometry

A 3D representation of muscle tissue structure is constructed using a previously-described algorithm [22]; it includes cylindrical fibers arranged in regular arrays and a network of capillaries, small precapillary arterioles and postcapillary venules. The dimensions of the tissue studied are 200 μm width (x-axis), 208 μm height (y-axis) and 800 μm length (z-axis). The fiber and vascular geometry can be specified using different methods, including tissue-specific geometries with irregular-shaped fibers obtained from in vivo imaging; tissue dimensions can also be extended.

Flow Module

The in vivo hemorheological model [59, 60] is applied to calculate the distribution of blood flow rate (Q) and discharge hematocrit (H D ), among all the capillary segments under steady state conditions during exercise. The governing equations are derived from the mass conservation law for volumetric blood flow rate and red blood cell flow rate at the jth node (vascular bifurcation), as follows:


Here Q ij is the volumetric flow rate in vessel segment ij, a cylinder whose ends are the ith and jth nodes of the network. Flow rate is calculated as Q ij = πR4(p i -p j )/(8ηL), where p j is the hydrodynamic pressure at node j, R and L are the radius and length of the segment, and η is the apparent viscosity which is a function of R and H D (Fahraeus-Lindqvist effect). These equations are supplemented by the empirical equations governing red blood cell-plasma separation at vascular bifurcations. The system of nonlinear algebraic equations for all N segments is solved with respect to pressure and discharge hematocrit, from which flow in each segment is calculated.

Oxygen Module

Oxygen delivery from the microvasculature to skeletal muscle myocytes is one of the key functions of microcirculation. During exercise, oxygen consumption may increase many folds compared to resting state, affecting both extravascular and intravascular oxygen transport. The oxygen model consists of two partial differential equations, Eqns. 3 and 4, governing extravascular and intravascular oxygen transport, respectively, assuming muscle fibers and interstitial space are a single tissue phase [41, 60].

Oxygen tension in the tissue, P O2 = P(x,y,z,t), is governed by free oxygen diffusion, myoglobin-facilitated diffusion, and oxygen consumption by tissue cells:


Here and D Mb are the diffusivities of oxygen and myoglobin in tissue respectively; S Mb is the oxygen-myoglobin saturation; α tis is the oxygen solubility in tissue; is the binding capacity of myoglobin with oxygen; M c is the oxygen consumption rate coefficient for Michaelis-Menten kinetics; P crit is the critical P O2 at which oxygen consumption equals to 50% of M c ; and S Mb is defined as P/(P+P 50,Mb ) assuming the local binding equilibrium between oxygen and myoglobin, where P 50,Mb is the P O2 necessary for 50% myoglobin oxygen saturation.

Oxygen transport in the blood vessels is governed by:


Here is the oxygen-hemoglobin saturation in blood vessel; α RBC and α pl are oxygen solubility in red blood cell and plasma, respectively; P b is the oxygen tension in blood plasma; ν b is the mean blood velocity (ν b = Q/(πR2)); H T and H D are the tube and discharge hematocrit, calculated from blood flow model; is binding capacity of hemoglobin with oxygen; ξ is the distance along a vessel's longitudinal axis; J wall is the capillary wall flux; and is defined as assuming the binding equilibrium between oxygen and hemoglobin, where P 50,Hb is the P O2 necessary for 50% hemoglobin oxygen saturation.

In addition, continuity of oxygen flux at the interface between blood vessels and tissue yields:


where n is the unit normal vector, P wall is the local P O2 at the vessel wall and k 0 is the mass transfer coefficient estimated from an empirical equation k 0 = 3.15+3.26H T -9.71+9.74(H T )2 + 8.54()2. The system of nonlinear partial differential equations was solved using the finite difference method, with a grid size of 1 micron as described in [60].

VEGF module

VEGF is the most-studied molecular factor involved in angiogenesis, including exercise-induced angiogenesis. Among several splice isoforms in the VEGF family, VEGF120 and VEGF164 (in rodents; human isoforms are VEGF121 and VEGF165) are considered to be the major pro-angiogenic cytokines that induce proliferation and migration of endothelial cells. The molecular weights of VEGF164 and VEGF120 are 45 and 36 kDa respectively and thus their diffusion coefficients are slightly different; in addition, VEGF164 binds the heparan sulfate proteoglycans (HSPGs) while VEGF120 does not and thus the shorter isoform diffuses more freely through the ECM.

A reaction-diffusion model [22] is used to predict molecular distribution in the interstitial space and on the endothelial surface. The governing equations for VEGF164 and VEGF120 are:


where C V164 , C V120 , C H and CV 164•Hare the concentrations of VEGF164, VEGF120, HSPG and VEGF164·HSPG complex; D V164 and D V120 are the diffusion coefficients of VEGF164 and VEGF120; and k on,V164,H and k off,V164,H are the association and dissociation rate constants between VEGF164 and HSPG. The boundary conditions for VEGF164 and VEGF120 at the surfaces of muscle fibers and endothelial cells, and the complete details of ligand-receptor interactions, were described in [61].

The model describes the secretion of two VEGF isoforms from the muscle fibers, molecular transport of each isoform in the interstitial space, binding of VEGF164 to HSPG in the ECM, VEGF164/120 binding to VEGFR2 at the endothelial cell surface and internalization of these ligand-receptor complexes. The model also considers VEGFR1 and neuropilin-1 (NRP1) coreceptor binding with VEGF ligands.

We previously applied an empirical equation to describe the relationship between P O2 and VEGF secretion rate to estimate local fiber VEGF secretion [22]. The empirical relationship was derived by combining experimentally-based relationships between intracellular HIF1α concentration and P O2 in vitro, and HIF1α concentration and VEGF secretion in vivo in skeletal muscle. However, PGC1α has recently been found to be another important regulator of VEGF in exercise [4, 62]. It was first discovered as a cold-inducible transcriptional coactivator for nuclear hormone receptors in brown fat and an enhancer of mitochondrial metabolism and function [63, 64]. Recently, a series of experimental studies [6567] have shown that VEGF gene expression and protein levels are highly dependent on the presence and concentration of PGC1α, through HIF-independent and HIF-dependent pathways. Thus, we have modified the equation by incorporating the effect of PGC1α:


Here S VEGF is the VEGF secretion rate, S B,VEGF is basal secretion rate at normoxic levels of [HIF1α]. It is defined as a function of [PGC1α], written as a sigmoidal form, , where [PGC1α] is the PGC1α concentration normalized relative to normoxic expression for wild type skeletal muscle, n p is the Hill constant, and k h , A, and B are empirical constants. S 0,VEGF is defined as basal VEGF secretion rate at normoxic levels of [HIF1α] and [PGC1α] for wild type skeletal muscle. Equations for oxygen-dependent [PGC1α] under wild type, knockout and over-expression conditions are shown in additional file 1 (Eqns. S4-S6). Data fitting based on an array of experimental data [4, 66, 67] results in the following parameters: A = 2.3167, B = 0.35, k h = 2.5641, n p = 1.086, = 3.

Cell Module

The Cell module is adapted from our 3D agent-based in vitro model [33] to describe how capillary endothelial cells respond to stimuli, specifically VEGF concentration and gradients, during the time course of sprout formation. The model applies logical rules to define cell activation, elongation, migration and proliferation events, based on extensive published experimental data [33]. The model makes predictions of how single-cell events contribute to vessel formation and patterns through the interaction of various cell types and their microenvironment.

The primary rules used in the in vitro model [33] are specified as follows: i) Endothelial cells are activated at an initial time point and the number of activated cells is constrained by a specified maximum number per unit capillary length. This activation initiates development of a tip cell segment of a sprout, later followed by the formation of a first stalk cell segment. ii) Following invasion of new sprouts into the tissue by extending the leading tip and stalk cell segments, the tip cell continues to migrate in the interstitial space following VEGF gradients and moves towards higher VEGF concentration. In addition, the tip cell can also proliferate with a certain probability, and the stalk cell can elongate and proliferate in a specified fashion; note that the probability of tip cell proliferation is much smaller than that of the stalk cells. The combined effect of these two cell phenotypes can be simulated as a biophysical push-pull system. iii) Branching occurs with a specified probability after a designated time threshold has elapsed at either a stalk or tip cell. The branching angle is selected stochastically and is less than 120 degrees. In the original model [33] the frequency of branching events during the spouting process was a function of the expression of ligand Dll4 and receptor Notch on the endothelial cells. Details of other rules and the parameters were described in [33].

To simulate in vivo conditions in the skeletal muscle vessel network, we modified some of the previously-defined rules and introduced additional rules to the model. Since muscle fibers and vasculature occupy respectively 79.7% and 2.5% of the tissue volume, the interstitial space totals 17.8%. Hence the freedom of tip endothelial cell migration during sprout formation is constrained to occur in a small volume of interstitial space. Note that in the model the endothelial cells consist of cylindrical cell segments (10 μm length and 6 μm diameter per capillary segment; 4 segments per cell defined in this study); the rules are formulated for these segments rather than for whole cells. This part of the model can be readily modified. The additional rules imposed in this study are as follows: i) Elongation or migration of cell segments follows the original rules as developed and defined in [33], except when the cell may encounter a fiber by following the growth factor gradient, we assume that the tip cell filopodia will sense the fiber and instead the cell follows the second largest VEGF gradient direction alternatively to elongate or migrate. ii) Anastomoses are formed when the tip cell senses an existing capillary or a sprout within 5 microns. iii) Since the function of Dll4-Notch is not clearly defined in skeletal muscle, their effects are not taken into account in the current simulations, but this effect can be readily added. For model simplification and demonstration purpose, tip cell elongation and the branching are not allowed in the present study.

Integration of computational modules

The development of an anatomically-, biophysically- and molecular-detailed spatio-temporal model by integration of different modules is a novel and challenging task. One of the main objectives of this study is to create a platform for integration of different modules written in different programming languages and using mixed modeling methodologies. The component modules may be created in the same or different laboratories, and could also be selected from a public model database. The difficulty of this task stems from the fact that few standards and open-source software/libraries for PDE solvers and ABMs exist. As a result, modules are dependent on their native languages and on differential equation solvers, making the integration difficult. Another problem facing the integration of modules is how to define and implement the connectivity between them, i.e., the exchange of data between the modules. Here we solve these two problems using a novel computational infrastructure and object-oriented design as described below.

Computational Infrastructure

To overcome the language barrier between the four modules selected in this study (Fortran for the Flow module, C/C++ for the Oxygen and VEGF modules and Java for Cell Module, Figure 1A), we choose Java to design the controller, which provides a flexible high-level interface and object-oriented facilities. Instead of rewriting the codes in each module in Java, we use a mixed-language programming environment to link the modules and save repetitive effort. The native codes in Fortran and C run faster than Java, and this compromise solution can also inherit advantages of these two languages. Another important technical aspect that renders this hybrid system feasible is the existence of Java Native Interface (JNI) API to convert functions and data type from native codes (Fortran and C) automatically to Java. To fulfil this purpose, we redesigned the native codes for the Flow, Oxygen, and VEGF modules to the format of functions and subroutines, and compiled these codes into the Java-readable libraries, turning all four modules into "pluggable" libraries that can be called by the controller coded in Java (Figure 1B). Furthermore, these libraries can be dynamically linked, making the simulation of dynamic angiogenesis processes feasible. Thus, using the controller as a bridge between each module, communication between different modules is relatively easy to implement. Last, it is easy to use Java to implement the connection between core codes and a new parameter database file used by the four selected modules.

To achieve high performance of native codes, parallel computing is implemented in the Oxygen and VEGF modules, as they require extensive computing resources. The current version of the modules adds OPENMP (open multi-processing) support, an industry standard for memory-shared parallel systems, to shorten the simulation time when the PDE solver is called by the controller. Numerical simulation time for the Oxygen and VEGF PDEs has been sped up ~5 fold using an 8 quad-core processor.

Object-oriented design

As a starting point to integrate all four modules, we focused on robust design of the controller, providing the connectivity between the modules, rather than providing solvers/software for mathematical models (i.e., rather than focusing on the capability of solving equations or agent-based models numerically). This will not impose constraints on the modeling methodology used in the modules, and it will allow a wide choice of modules to be integrated, providing more flexibility. More detail on the procedure for integration of modules at multiple time scales using object-oriented classes is given below in the section "Integration of modules at multiple time scales".

Regarding design at the upper level (i.e., package design), we devised four sub-systems: Biosystem, Process, IO (Input/Output) and Exception as shown in lower panel in Figure 1A. The Biosystem is a repository of hierarchical structure of biological and biophysical information in the tissue as shown in Figure 2. The Process is composed of angiogenesis subprocesses occurring at various stages in growth process. The IO is used to provide an interface for the interaction between the system and the user, for example, to read parameters used for the simulation of the Oxygen and VEGF modules. The Exception provides the capability to debug and handle runtime errors.

Figure 2
figure 2

Object-oriented design for the angiogenesis modeling package. Major classes across tissue and cell scales in the modeling controller are shown. They include SkeletalMuscle, Myofiber, Vessel, Grid, Segment and Node classes in the Biosystems subpackage, and BloodFlow, O2Diffusion, VEGFRxnDiffusion, CellSprouting, and StartAngio classes in the Process subpackage. The hierarchical structure of relationships between the classes is represented by arrows.

At the lower level (i.e., class design), we applied object-oriented concepts to describe a hierarchical structure of skeletal muscle (EDL) and events in the angiogenesis process. The class diagram shown in Figure 2 depicts the relationships among the major classes proposed. In the Biosystem sub-package, six classes are defined to describe the entities composing the tissue of interest: SkeletalMuscle, Myofiber, Vessel, Grid, Node and Segment. In addition, angiogenesis event classes defined in the sub-package Process include BloodFlow, O2Diffusion, VEGFRxnDiffusion, CellSprouting and StartAngio.

The Node and Segment classes in the Biosystem sub-package are defined below the cell scale. The Node class represents the circular surface of cylindrical segments at their ends. It contains spatial information, including the circle center position and the circle radius. It may also have biophysical information such as blood pressure, flow velocity and hematocrit if the node is contained in a blood vessel. The Segment class represents a cylinder in 3D space, corresponding to either a fraction of a blood vessel or muscle fiber (assuming the fiber and blood vessel are cylindrical in shape). A Segment object contains two Node objects at its ends and the length of the cylinder. Each Segment also contains biophysical information such as blood flow, pressure and hematocrit, and VEGF receptor density if the segment type is a blood vessel.

At the tissue scale, SkeletalMuscle, Myofiber, Vessel and Grid are defined in the Biosystem sub-package. The Vessel class is used for representation of microvessels including capillaries, venules and arterioles. The Myofiber class is used for representation of skeletal muscle myocytes. These two classes are each composed of a series of segments. Interstitial space information including VEGF and HSPG concentrations is described in the Grid class. The Grid class also contains the local P O2 value. Ensemble components of fiberTissue, capillary, venule, arteriole, voxel, tissueSize and gridUnitSize render the SkeletalMuscle class, which describes detailed skeletal muscle structure.

In the Process sub-package, four classes including BloodFlow, O2Diffusion, VEGFRxnDiffusion and CellSprouting provide connections with Java-readable libraries compiled from their corresponding modules. The calling of these Java classes involves three steps: i) The skeletal muscle object (realization of SkeletalMuscle class) will be initialized and then transferred as the input to each module. ii) Each specific module will compute using their intrinsic numerical solvers and then the results will be transferred to the Java interface class. iii) Finally, the skeletal muscle object will be updated with the solutions. The StartAngio class contains methods defined to specifically simulate exercise-induced angiogenesis.

Developing the computational environment

The simulation experiments were run on a computer with 64 bit Linux Ubuntu system, 8 quad-CPU and 128 Gbyte memory. Eclipse is used as an Integrated Development Environment for coding purposes. JDK (Java Development Kit) 1.6.16 (Oracle, Redwood Shores, CA) is used as Java compiler, and Intel Fortran/C++ compiler suite (v.11.1) (Intel, Santa Clara, CA) is used as Fortran/C++ compiler. We also incorporate Java 3D™ API (Oracle, Redwood Shores, CA) for 3D programming purposes and the Log4j package to log all the runtime messages for the purpose of debugging. We use Bazaar as our source code version control system since as an industry standard for software development it provides support for a large scale project development by a team of programmers. It has advantages in terms of branching, merging and keeping revision versions. The GNU Make tool is chosen for automation of building executable programs and libraries from source codes, and running programs from binary codes. In particular, it is useful for a hybrid system (i.e., mixed programming environment) in a single program. The MASON package is used as the agent-based modeling library. Unit tests are performed using the JUNIT 4.0 package

Integration of modules at multiple time scales

We performed the integration of the modules using a sequential method, that is, the modules are run sequentially rather than in parallel. This is based on a time scale analysis of each process integrated into the multiscale angiogenesis model. We computed the module at the fastest time scale first. The outcome of that module was considered as a pseudo-steady state and used as an input feeding into the modules at slower time scales. This continued sequentially until we computed the module at the slowest time scale. Specifically, the blood flow regulation and oxygen distribution modules reach equilibrium within seconds to minutes; VEGF gradients at time scales of minutes to hours; and capillary sprouting from hours to weeks. Thus we computed the steady-state flow and oxygen module first, then used their simulation results as an input to VEGF reaction-diffusion model and run PDE solver to compute VEGF profile, and then finally run the agent-based model to simulate angiogenesis patterns for single-bout exercise. When endurance exercise for days or weeks is simulated, the updated model geometry will be used as the new input to run flow, oxygen and cell modules sequentially.

Example run of a simulation of single-bout exercise

Using the object-oriented design concept, we constructed a controller which is capable of interacting with each individual module defined in our multiscale model. A run of the simulation starts with input of geometry files and parameter database file (written in database file format). These inputs initialize an object of "SkeletalMuscle" class with the 3D coordinate information of segments and nodes for skeletal muscle fiber and blood vessel network. The parameter database file assigns values to the biochemical and physiological parameters defined in the model. Following this, the controller calls the flow module as a dynamic linked library (dll) as described above, using the geometric and biochemical parameters as input. The flow module has its own built-in solver that returns (to the controller) results including blood pressure, hematocrit, viscosity and flow rates in the various vessel segments. The information is stored in the controller "SkeletalMuscle" object. This object, now with the updated flow information and the increased oxygen consumption rate, is passed by the controller to the oxygen module. This module computes the oxygen tension at each grid point defined within the skeletal muscle and passes the values back to the controller to be stored in the voxel field (Grid class) of the "SkeletalMuscle" object (element P O2 is defined in Grid Class). The controller passes the O2-updated object to the VEGF module, which includes O2-PGC1-HIF-VEGF empirical equations. Upon completing its simulation, the VEGF module passes the VEGF, HSPG, and other concentrations throughout the tissue back to the controller. This spatial concentration profile information is then passed to the cell module to simulate capillary growth. During the 8-hour post-exercise period, due to the cessation of exercise, blood flow rate and oxygen consumption rate return to basal level. We thus assume that capillary sprouting is the only active module during this period; blood flow rate and oxygen distributions are assumed to be at basal levels.


Using the integration strategy and simulation package described above, we performed a series of computational experiments to simulate activity-induced angiogenesis in EDL during a single bout of exercise. First, a 3D simulation was conducted to sequentially integrate modules from blood flow to oxygen convection-diffusion, to VEGF diffusion-reaction, to capillary sprouting. Second, sensitivity analysis was used to examine sensitivities of key biochemical, biophysical and physiological parameters to capillary growth. Third, angiogenesis patterns were compared between different exercise intensities, and high- or low-inspired oxygen conditions.

Sequential integration of modules from blood flow to capillary sprouting

Single-bout exercise was simulated by dynamic integration of four modules: Flow, Oxygen, VEGF and Cell. The simulations are based primarily on physiological parameters used in our previously developed models [33, 60, 61], as shown in Table 1. Most of the parameters were taken from experimental data for rat skeletal muscle while others were adopted from measurements in other tissues or theoretical estimates when parameters are unavailable from literature (e.g., diffusion coefficients for VEGF164/120, secretion rates for VEGF164/120). Here we describe several key parameters used in the simulation.

Table 1 Parameters of the multiscale model*

Oxygen consumption rate, M c , increases to 2- to 50-fold the basal rate during exercise [68]. Previous models simulated the system under 6- and 12-fold oxygen consumption [22, 61]. Here we used oxygen consumption of 9-fold the basal rate, corresponding to moderate exercise intensity conditions. These moderate exercise conditions are typical of experimental rat (mouse) aerobic treadmill exercise [69, 70]. Another important parameter of the exercise environment is oxygen saturation of arterioles feeding into capillaries (S O2A ). The model uses values of 0.6 and 0.3 for S O2A , corresponding to a normoxic environment and a low oxygen environment (hypoxic hypoxia). Hydrodynamic pressure drop between small arteries and venules (ΔP) is specified at 10 mmHg and inlet discharge hematocrit is 0.4, consistent with physiological observations [41, 71]. For exercise-induced angiogenesis, experimental observations show that VEGF mRNA expression is elevated significantly within 1 hour post-exercise, and it remains elevated until returning to basal levels after 8 to 24 hours [69, 70, 7274]. Due to the lack of experimental measurement of VEGF protein secretion rate from the myofibers in vivo, we assume that VEGF protein expression follows the similar time course as mRNA. Thus, we consider VEGF secretion as a step function starting at high levels at the onset of the post-exercise period, remaining constant for an 8-hour interval, and finally ending after 8 hours. We simulate the system for eight hours as the period of capillary growth; beyond this time elevated VEGF secretion ceases and the stimulus to sprout declines.

Using skeletal muscle geometry (Figure 3A) the blood flow simulation results (Figure 3B) show heterogeneous distribution of blood flow velocity in the microcirculatory network, with variation in velocity between blood vessels across the x-y plane. The mean flow velocity of all capillaries in the network is 0.11 cm/s and total blood flow is 227.82 ml·(100 g tissue)-1 min-1, consistent with experimental measurements [7577]. EDL tissue displays oxygen spatial heterogeneity not only in the x-y plane, but also in the z-direction, as shown in the 3D oxygen distribution plot (Figure 3C), generated from the Oxygen module simulation results. Most of the hypoxic tissue space has few nearby capillaries, confirming that the oxygen tension gradient and distribution is dependent on capillary spatial positions in the skeletal muscle space. Calculating all oxygen tensions in the tissue, we obtained an average tissue P O2 of 13.23 mmHg with a coefficient of variation 0.28.

Figure 3
figure 3

3D simulation of blood flow, oxygen, and VEGF distribution during the single-bout exercise: A) 2D cross section of skeletal muscle: gray circles represent fibers and red circles represent capillaries; B) Blood flow velocity distribution in skeletal muscle microvascular network; C) Oxygen distribution throughout the tissue; D) VEGF secretion level along the muscle fibers; E) total VEGFR-bound VEGF distribution (including both 120 and 164 isoforms) on vascular surface; and F) free VEGF concentration distribution in interstitial space.

Using oxygen tension distribution as input to the VEGF simulation parsed by the integration core package, we calculated VEGF secretion rate (Figure 3D). The average VEGF secretion level is 2.12-fold with a coefficient of variation of 0.32. The increase in VEGF secretion levels results in the increase of soluble VEGF concentration in tissue interstitial space (average free VEGF concentration is 1.9 pM). Simulation results from the VEGF module (Figure 3E) show that VEGF gradients in the longitudinal direction (along the fibers) are much smaller than in the cross-sectional planes, as indicated by their mean values (0.56% and 3.93% change in VEGF concentration over 10 microns, respectively). Similarly, VEGF binding to the endothelial surface receptor VEGFR2 is non-uniform across the vascular network (Figure 3F). The mean surface concentration of total receptor-bound VEGF is 1023 × 10-7 pmol/μm2 (coefficient of variation of 0.20), of which 8% is VEGF120 and 92% is VEGF164; NRP1-facilitated VEGF120-VEGFR1 and VEGF164-VEGFR2 complexes contribute 4.69% and 67.09% of the total receptor complexes, respectively.

We simulated the capillary activation process using our modified in vivo cell model. The precondition for an endothelial cell to be considered as activated is that any of its cell segments is exposed to a VEGF concentration greater than a specified VEGF activation threshold. A random selection algorithm is used to determine the sprouting position in the activated endothelial cells and the number of activated endothelial cells is also constrained by the allowed maximum number of activated tip cells per capillary (5 cells per capillary) and the minimum distance between tip cells (40 μm). These rules allow us to define activated tip cell and cell segments as the sprouting point for new capillary sprouts. Using a VEGF activation threshold of 2.0 pM, 2% percent of cells are activated and these cells initiate the sprouting process.

We next simulate cellular sprouting and predict dynamic development of nascent vessels during 8 hours post-exercise, as defined in the timeline scheme in Figure 4A. Figure 4B-F shows new vascular network formation at 1, 3, 5, 7 and 8 hours following exercise. One of the assumptions in the model is that new capillary formation does not affect muscle fiber oxygenation, VEGF secretion or VEGF internalization, due to the slow capillary maturation process. Thus, a feedback loop is not implemented for a single bout of exercise, and a static VEGF profile is used as input for the cell model. This is a condition that can be relaxed in future models. We also assume periodic conditions for capillary growth, i.e., when capillary vessel sprouts encounter the boundary, they continue to grow beyond the boundary (e.g. arrow 1 at Figure 4C) and are represented by a corresponding capillary growing from the opposite face (e.g. arrow 2 at Figure 4E). As shown in Figure 4D-F, capillary sprouts on the same capillary initially migrate in a similar direction, while at a later time, they may form anastomoses with adjacent vasculature while migrating (arrows 3, 4, 5). Using parameter values listed in Table 1, capillary proliferation and migration rules defined in the cell module result in a 12.6% increase in total capillary length compared to the initial length before exercise (Figure 4F). However, for the conditions of a single bout of exercise considered in this study most of these new sprouts are likely to retract; sustained exercise conditions with formation of multiple functional anastomoses will be considered in subsequent studies.

Figure 4
figure 4

3D simulation of capillary network growth during the single-bout exercise. A) Time course of capillary growth was simulated based on the timeline scheme for exercise-induced angiogenesis. The simulation leads to new capillary network at: B) 1 h post-exercise; C) 3 h; D) 5 h; E) 7 h; F) 8 h. When new capillary segments grow out of boundaries (Arrow 1), capillary will grow from corresponding periodic boundaries (Arrow 2). Arrows 3-5 refer to the capillary anastomoses during growth. Vessels colored in red represent the pre-existing microvascular vessels, ones in light blue are the newly-developed blind-ended capillaries, and ones in dark blue are new capillary anastomoses. G) Relative capillary length change vs time.

Sensitivity analysis of cellular parameters to define viable parameter space

We further explore the effect of parameters defined in the multiscale model of capillary growth. For this purpose, all biochemical, biophysical, and physiological parameters are extracted from each module and are deposited in a parameter database file. We show an example of how endothelial cell activation affects the capillary growth pattern and growth rate (two key variables related to capillary growth). We vary the threshold level of VEGF for endothelial cell activation within a range of 1.5 to 2.3 pM and this change significantly affects the capillary growth rate. Figure 5A shows that above a VEGF threshold of 2.3 pM, none of the endothelial cells are activated and the capillary network has no growth. Over the range of 2.3-1.8 pM, capillary growth rate increases slightly as the threshold value decreases. Below 1.8 pM, capillary growth increases dramatically since a large number of capillaries are activated. The response to these threshold levels depends on the mean concentration of VEGF in the tissue (here, 1.9 pM).

Figure 5
figure 5

Sensitivity of the capillary length and the number of anastomoses formed with respect to VEGF threshold for EC activation. A) Relative capillary length as a function of time and VEGF threshold; B) Number of anastomoses formed as a function of VEGF threshold at 4 h and 8 h. Simulation sample size is five for each VEGF threshold at a given time.

The number of anastomoses is another important variable that indicates vascular adaptation of the microvasculature since the increase of blood flow results from these newly-formed vessel loops (and not from blind-ended vessels). We show in Figure 5B that the number of anastomoses increases in a similar fashion when threshold values decrease from 2.3 pM to 1.5 pM. Interestingly, more than half of the anastomoses are predicted to form during the first 4 hours.

Simulation of exercise-induced angiogenesis under different conditions

Varying experiment-based parameters defined in the Oxygen modules, we can simulate vessel growth in skeletal muscle during different exercise conditions. For example, different oxygen consumption rates, M C, correspond to different exercise intensities and arterial oxygen saturation values S O2A are associated with different available oxygen levels. We compare capillary growth patterns between three exercise conditions: i) M C = 1 × 10-3 ml O2/ml tissue/s and S O2A = 0.6 (light intensity exercise in normoxia; results not shown); ii) M C = 1.5 × 10-3 and S O2A = 0.6 (moderate intensity exercise in normoxia; Figure 6B and 4F); iii) M C = 1.5 × 10-3 and S O2A = 0.3 (moderate intensity exercise in hypoxic hypoxia; Figure 6A). Simulation results show that during low-intensity exercise (6-fold basal level oxygen consumption rate) with normoxic conditions, the capillary network will not change (case i). Compared to exercise condition of a 9-fold basal oxygen consumption rate in normoxia (case ii), exercise under hypoxic hypoxia conditions at a 9-fold basal oxygen consumption rate further increases capillary growth to 125% (case iii).

Figure 6
figure 6

Angiogenesis pattern during different exercise conditions at 8 h post-exercise. A) Oxygen consumption at 9-fold the basal level under hypoxic hypoxia conditions (S O2A = 0.3); B) Oxygen consumption at 9-fold the basal level under normoxic conditions (S O2A = 0.6).

Discussion and Conclusions

Using a module-based computational platform, we connected four modules, each of which describes different aspects of the angiogenesis process (i.e, blood flow, oxygen distribution, VEGF distribution, and capillary growth). This architecture allowed the individual modules to remain in their native language and required rewriting only of the input/output routines to handle the passing of common data between the submodules and the integrator. This significantly reduced the time required to integrate the available modules. JNI plugins are used to link the codes written in different languages, and they allow ready communication between the core package and each module. As a validation of our simulation strategy and integration package, we have conducted a series of computational experiments to simulate dynamic activity-induced angiogenesis during a single bout of exercise, using rat EDL skeletal muscle. We illustrated how the computational design of a core controller module is beneficial to integrating different aspects of biological knowledge with a variety of numerical simulation schemes: i) sequential integration; ii) sensitivity analysis; and iii) simulation of various exercise conditions.

Object-oriented programming (OOP) is attracting more modelers' attention and has been increasingly used in biomedical computational projects [56, 57, 78]. We applied OOP concepts to describe a complex biological system which encompasses biological components and processes from the molecular level to the tissue level. With our model architecture, we also depicted the hierarchical relationships among the different classes abstracted from such systems. For the core package, Java was used since it has a natural object-oriented design and method inheritance structure that is intuitively appealing to describe the biological system.

Our approach is complementary to other ongoing projects in multiscale modeling fields, e.g. CompuCell [79], the Multiscale Systems Immunology (MSI) simulation framework [80] and Virtual Cell [81]. These projects share the same objective: the development of multiscale modeling software or libraries with a versatile user interface. CompuCell implements cellular Potts models to simulate cell dynamics, MSI uses agent-based models and Virtual Cell uses differential equations. While we did not attempt to provide standards or software for the angiogenesis models, we chose to focus on the extensibility and generality of an integration methodology that allows additional individual modules to be integrated at any time, and does not impose constraints on the modeling techniques employed.

As a starting point, current study shows an example of the simulation by integrating four modules. We used existing models for the first modules to be able to compare the results of the integrated model to published results of each submodule. This computational paradigm will allow us to incorporate additional modules developed in the laboratory, e.g., MMP [23, 24], HIF1α [82] and FGF2 [18], using an SBML/CellML-based open source simulation packages to represent the molecular interaction networks. In addition, this platform can incorporate other relevant published models from the literature or public model repositories, e.g. CellML Models Repository and BioModels Database It is also feasible to improve on or add to the individual modules currently comprising the integrated model. For example, the geometry used by the modules can be updated to that of a different muscle type without rewriting the other module codes; or more detailed information such as the fiber type of each muscle fiber can be included as an additional data notation, with corresponding changes in local oxygen consumption and O2 diffusion coefficient and solubility.

Correspondingly, the modeling controller module will be further expanded to include the classes representing other components, which may be relevant to skeletal muscle angiogenesis when new biological information and models are available. In the long run, definition of unified XML-based notation to describe structure and spatial information of physiological systems, ordinary/partial differential equations, agent-based rules and system parameters will make each module more portable and more easily integrated. Progress has been initiated in this direction such as insilicoML [83]. In addition, development of the ontology for the field of angiogenesis modeling will also be beneficial to the sharing, usage and integration of different modules, similar to gene ontology [84], and the open biological and biomedical ontologies [85]. These ongoing and future developments should contribute to the overall efforts of the Systems Biology and Physiome Projects [46, 86], and could be applied to cardiovascular disorders such as coronary and peripheral artery diseases [87].

Along with this framework for modeling the system, there is now the opportunity to include other useful functionality such as local and global parameter sensitivity modules, or parameter optimization modules. While local parameter sensitivity calculations can be performed with the current model, the computational expense of the model (in particular the PDE solver) is such that parameter estimation would be very slow. We aim to use cluster computing in order to permit these activities in the future.

In summary, we have developed the computational methodology capable of the integration of biologically- and computationally-heterogeneous modules and this systems biology strategy can be applied to larger scale integration of computational models of angiogenesis in skeletal muscle.


  1. Figg WD, Folkman MJ, SpringerLink (Online service): Angiogenesis An Integrative Approach From Science to Medicine. 2008, Boston, MA: Springer Science+Business Media, LLC

    Google Scholar 

  2. Popel AS: Theory of oxygen transport to tissue. Crit Rev Biomed Eng. 1989, 17: 257-321.

    CAS  PubMed  Google Scholar 

  3. Hirota K, Semenza GL: Regulation of angiogenesis by hypoxia-inducible factor 1. Crit Rev Oncol Hematol. 2006, 59: 15-26. 10.1016/j.critrevonc.2005.12.003.

    Article  PubMed  Google Scholar 

  4. Arany Z: PGC-1 coactivators and skeletal muscle adaptations in health and disease. Curr Opin Genet Dev. 2008, 18: 426-434. 10.1016/j.gde.2008.07.018.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  5. Olfert IM, Howlett RA, Wagner PD, Breen EC: Myocyte vascular endothelial growth factor is required for exercise-induced skeletal muscle angiogenesis. Am J Physiol Regul Integr Comp Physiol. 2010, 299: R1059-1067. 10.1152/ajpregu.00347.2010.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  6. Egginton S: Invited review: activity-induced angiogenesis. Pflugers Arch. 2009, 457: 963-977. 10.1007/s00424-008-0563-9.

    Article  CAS  PubMed  Google Scholar 

  7. Bloor CM: Angiogenesis during exercise and training. Angiogenesis. 2005, 8: 263-271. 10.1007/s10456-005-9013-x.

    Article  PubMed  Google Scholar 

  8. Lundby C, Calbet JA, Robach P: The response of human skeletal muscle tissue to hypoxia. Cell Mol Life Sci. 2009, 66: 3615-3623. 10.1007/s00018-009-0146-8.

    Article  CAS  PubMed  Google Scholar 

  9. Mac Gabhann F, Popel AS: Systems biology of vascular endothelial growth factors. Microcirculation. 2008, 15: 715-738. 10.1080/10739680802095964.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  10. van Hinsbergh VW, Koolwijk P: Endothelial sprouting and angiogenesis: matrix metalloproteinases in the lead. Cardiovasc Res. 2008, 78: 203-212. 10.1093/cvr/cvm102.

    Article  CAS  PubMed  Google Scholar 

  11. Gerhardt H, Golding M, Fruttiger M, Ruhrberg C, Lundkvist A, Abramsson A, Jeltsch M, Mitchell C, Alitalo K, Shima D, Betsholtz C: VEGF guides angiogenic sprouting utilizing endothelial tip cell filopodia. J Cell Biol. 2003, 161: 1163-1177. 10.1083/jcb.200302047.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  12. Qutub AA, Mac Gabhann F, Karagiannis ED, Vempati P, Popel AS: Multiscale models of angiogenesis. IEEE Eng Med Biol Mag. 2009, 28: 14-31. 10.1109/MEMB.2009.931791.

    Article  PubMed Central  PubMed  Google Scholar 

  13. Peirce SM: Computational and mathematical modeling of angiogenesis. Microcirculation. 2008, 15: 739-751. 10.1080/10739680802220331.

    Article  PubMed Central  PubMed  Google Scholar 

  14. Byrne HM: Dissecting cancer through mathematics: from the cell to the animal model. Nat Rev Cancer. 2010, 10: 221-230. 10.1038/nrc2808.

    Article  CAS  PubMed  Google Scholar 

  15. Anderson AR, Chaplain MA: Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull Math Biol. 1998, 60: 857-899. 10.1006/bulm.1998.0042.

    Article  CAS  PubMed  Google Scholar 

  16. Chaplain MA, Stuart AM: A mathematical model for the diffusion of tumour angiogenesis factor into the surrounding host tissue. IMA J Math Appl Med Biol. 1991, 8: 191-220. 10.1093/imammb/8.3.191.

    Article  CAS  PubMed  Google Scholar 

  17. Mac Gabhann F, Popel AS: Interactions of VEGF isoforms with VEGFR-1, VEGFR-2, and neuropilin in vivo: a computational model of human skeletal muscle. Am J Physiol Heart Circ Physiol. 2007, 292: H459-474. 10.1152/ajpheart.00637.2006.

    Article  CAS  PubMed  Google Scholar 

  18. Filion RJ, Popel AS: A reaction-diffusion model of basic fibroblast growth factor interactions with cell surface receptors. Ann Biomed Eng. 2004, 32: 645-663. 10.1023/B:ABME.0000030231.88326.78.

    Article  PubMed  Google Scholar 

  19. Stefanini MO, Wu FT, Mac Gabhann F, Popel AS: A compartment model of VEGF distribution in blood, healthy and diseased tissues. BMC Syst Biol. 2008, 2: 77-10.1186/1752-0509-2-77.

    Article  PubMed Central  PubMed  Google Scholar 

  20. Wu FT, Stefanini MO, Mac Gabhann F, Kontos CD, Annex BH, Popel AS: Computational kinetic model of VEGF trapping by soluble VEGF receptor-1: effects of transendothelial and lymphatic macromolecular transport. Physiol Genomics. 2009, 38: 29-41. 10.1152/physiolgenomics.00031.2009.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  21. Stoll BR, Migliorini C, Kadambi A, Munn LL, Jain RK: A mathematical model of the contribution of endothelial progenitor cells to angiogenesis in tumors: implications for antiangiogenic therapy. Blood. 2003, 102: 2555-2561. 10.1182/blood-2003-02-0365.

    Article  CAS  PubMed  Google Scholar 

  22. Mac Gabhann F, Ji JW, Popel AS: VEGF gradients, receptor activation, and sprout guidance in resting and exercising skeletal muscle. J Appl Physiol. 2007, 102: 722-734. 10.1152/japplphysiol.00800.2006.

    Article  CAS  PubMed  Google Scholar 

  23. Karagiannis ED, Popel AS: A theoretical model of type I collagen proteolysis by matrix metalloproteinase (MMP) 2 and membrane type 1 MMP in the presence of tissue inhibitor of metalloproteinase 2. J Biol Chem. 2004, 279: 39105-39114. 10.1074/jbc.M403627200.

    Article  CAS  PubMed  Google Scholar 

  24. Vempati P, Karagiannis ED, Popel AS: A biochemical model of matrix metalloproteinase 9 activation and inhibition. J Biol Chem. 2007, 282: 37585-37596. 10.1074/jbc.M611500200.

    Article  CAS  PubMed  Google Scholar 

  25. McDougall SR, Anderson AR, Chaplain MA: Mathematical modelling of dynamic adaptive tumour-induced angiogenesis: clinical implications and therapeutic targeting strategies. J Theor Biol. 2006, 241: 564-589. 10.1016/j.jtbi.2005.12.022.

    Article  PubMed  Google Scholar 

  26. Byrne HM, Chaplain MA: Mathematical models for tumour angiogenesis: numerical simulations and nonlinear wave solutions. Bull Math Biol. 1995, 57: 461-486.

    Article  CAS  PubMed  Google Scholar 

  27. Chaplain MA, Ganesh M, Graham IG: Spatio-temporal pattern formation on spherical surfaces: numerical simulation and application to solid tumour growth. J Math Biol. 2001, 42: 387-423. 10.1007/s002850000067.

    Article  CAS  PubMed  Google Scholar 

  28. Schugart RC, Friedman A, Zhao R, Sen CK: Wound angiogenesis as a function of tissue oxygen tension: a mathematical model. Proc Natl Acad Sci USA. 2008, 105: 2628-2633. 10.1073/pnas.0711642105.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  29. Peirce SM, Van Gieson EJ, Skalak TC: Multicellular simulation predicts microvascular patterning and in silico tissue assembly. FASEB J. 2004, 18: 731-733.

    CAS  PubMed  Google Scholar 

  30. Bauer AL, Jackson TL, Jiang Y: A cell-based model exhibiting branching and anastomosis during tumor-induced angiogenesis. Biophys J. 2007, 92: 3105-3121. 10.1529/biophysj.106.101501.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  31. Bentley K, Gerhardt H, Bates PA: Agent-based simulation of notch-mediated tip cell selection in angiogenic sprout initialisation. J Theor Biol. 2008, 250: 25-36. 10.1016/j.jtbi.2007.09.015.

    Article  CAS  PubMed  Google Scholar 

  32. Bentley K, Mariggi G, Gerhardt H, Bates PA: Tipping the balance: robustness of tip cell selection, migration and fusion in angiogenesis. PLoS Comput Biol. 2009, 5: e1000549-10.1371/journal.pcbi.1000549.

    Article  PubMed Central  PubMed  Google Scholar 

  33. Qutub AA, Popel AS: Elongation, proliferation & migration differentiate endothelial cell phenotypes and determine capillary sprouting. BMC Syst Biol. 2009, 3: 13-10.1186/1752-0509-3-13.

    Article  PubMed Central  PubMed  Google Scholar 

  34. Milde F, Bergdorf M, Koumoutsakos P: A hybrid model for three-dimensional simulations of sprouting angiogenesis. Biophys J. 2008, 95: 3146-3160. 10.1529/biophysj.107.124511.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  35. Harrington HA, Maier M, Naidoo L, Whitaker N, Kevrekidis PG: A hybrid model for tumor-induced angiogenesis in the cornea in the presence of inhibitors. Math Comput Model. 2007, 46: 513-524. 10.1016/j.mcm.2006.11.034.

    Article  Google Scholar 

  36. Das A, Lauffenburger D, Asada H, Kamm RD: A hybrid continuum-discrete modelling approach to predict and control angiogenesis: analysis of combinatorial growth factor and matrix effects on vessel-sprouting morphology. Philos Transact A Math Phys Eng Sci. 2010, 368: 2937-2960.

    Article  CAS  Google Scholar 

  37. Shirinifard A, Gens JS, Zaitlen BL, Poplawski NJ, Swat M, Glazier JA: 3D Multi-Cell Simulation of Tumor Growth and Angiogenesis. PLoS One. 2009, 4: e7190-10.1371/journal.pone.0007190.

    Article  PubMed Central  PubMed  Google Scholar 

  38. Hetherington J, Bogle IDL, Saffrey P, Margoninski O, Li L, Rey MV, Yamaji S, Baigent S, Ashmore J, Page K: Addressing the challenges of multiscale model management in systems biology. Comput Chem Eng. 2007, 31: 962-979. 10.1016/j.compchemeng.2006.10.004.

    Article  CAS  Google Scholar 

  39. Meier-Schellersheim M, Fraser ID, Klauschen F: Multiscale modeling for biologists. Wiley Interdiscip Rev Syst Biol Med. 2009, 1: 4-14. 10.1002/wsbm.33.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  40. Ayyadurai VAS, Dewey CF: Cytosolve: A scalable computational methodology for dynamic integration of multiple molecular pathway models. Cell Mol Bioeng. 2010, 4: 28-45. 10.1007/s12195-010-0143-x.

    Article  PubMed Central  PubMed  Google Scholar 

  41. Goldman D, Popel AS: A computational study of the effect of capillary network anastomoses and tortuosity on oxygen transport. J Theor Biol. 2000, 206: 181-194. 10.1006/jtbi.2000.2113.

    Article  CAS  PubMed  Google Scholar 

  42. Karagiannis ED, Popel AS: Distinct modes of collagen type I proteolysis by matrix metalloproteinase (MMP) 2 and membrane type I MMP during the migration of a tip endothelial cell: insights from a computational model. J Theor Biol. 2006, 238: 124-145. 10.1016/j.jtbi.2005.05.020.

    Article  CAS  PubMed  Google Scholar 

  43. Vempati P, Mac Gabhann F, Popel AS: Quantifying the proteolytic release of extracellular matrix-sequestered VEGF with a computational model. PLoS One. 2010, 5: e11860-10.1371/journal.pone.0011860.

    Article  PubMed Central  PubMed  Google Scholar 

  44. Hunter PJ, Borg TK: Integration from proteins to organs: the Physiome Project. Nat Rev Mol Cell Biol. 2003, 4: 237-243. 10.1038/nrm1054.

    Article  CAS  PubMed  Google Scholar 

  45. Hunter PJ, Crampin EJ, Nielsen PM: Bioinformatics, multiscale modeling and the IUPS Physiome Project. Brief Bioinform. 2008, 9: 333-343. 10.1093/bib/bbn024.

    Article  CAS  PubMed  Google Scholar 

  46. Popel AS, Hunter PJ: Systems Biology and Physiome Projects. Wiley Interdiscip Rev Syst Biol Med. 2009, 1: 153-158. 10.1002/wsbm.67.

    Article  PubMed Central  PubMed  Google Scholar 

  47. Bassingthwaighte J, Hunter P, Noble D: The Cardiac Physiome: perspectives for the future. Exp Physiol. 2009, 94: 597-605. 10.1113/expphysiol.2008.044099.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  48. Bassingthwaighte JB: Microcirculation and the physiome projects. Microcirculation. 2008, 15: 835-839. 10.1080/10739680802388906.

    Article  PubMed Central  PubMed  Google Scholar 

  49. Popel AS, Pittman RN: The Microcirculation Physiome. The Biomedical Engineering Handbook. Edited by: Bronzino JD, Peterson DR. Taylor & Francis, 4,

  50. Hunter P, Coveney PV, de Bono B, Diaz V, Fenner J, Frangi AF, Harris P, Hose R, Kohl P, Lawford P: A vision and strategy for the virtual physiological human in 2010 and beyond. Philos Transact A Math Phys Eng Sci. 2010, 368: 2595-2614. 10.1098/rsta.2010.0048.

    Article  Google Scholar 

  51. Cooper J, Cervenansky F, De Fabritiis G, Fenner J, Friboulet D, Giorgino T, Manos S, Martelli Y, Villa-Freixa J, Zasada S: The Virtual Physiological Human ToolKit. Philos Transact A Math Phys Eng Sci. 2010, 368: 3925-3936. 10.1098/rsta.2010.0144.

    Article  Google Scholar 

  52. Gianni D, McKeever S, Yu T, Britten R, Delingette H, Frangi A, Hunter P, Smith N: Sharing and reusing cardiovascular anatomical models over the Web: a step towards the implementation of the virtual physiological human project. Philos Transact A Math Phys Eng Sci. 2010, 368: 3039-3056.

    Article  Google Scholar 

  53. Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, Cornish-Bowden A: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003, 19: 524-531. 10.1093/bioinformatics/btg015.

    Article  CAS  PubMed  Google Scholar 

  54. Beard DA, Britten R, Cooling MT, Garny A, Halstead MD, Hunter PJ, Lawson J, Lloyd CM, Marsh J, Miller A: CellML metadata standards, associated tools and repositories. Philos Transact A Math Phys Eng Sci. 2009, 367: 1845-1867. 10.1098/rsta.2008.0310.

    Article  CAS  Google Scholar 

  55. Christie GR, Nielsen PM, Blackett SA, Bradley CP, Hunter PJ: FieldML: concepts and implementation. Philos Transact A Math Phys Eng Sci. 2009, 367: 1869-1884. 10.1098/rsta.2009.0025.

    Article  Google Scholar 

  56. Webb K, White T: UML as a cell and biochemistry modeling language. Biosystems. 2005, 80: 283-302. 10.1016/j.biosystems.2004.12.003.

    Article  CAS  PubMed  Google Scholar 

  57. Roux-Rouquie M, Caritey N, Gaubert L, Rosenthal-Sabroux C: Using the Unified Modelling Language (UML) to guide the systemic description of biological processes and systems. Biosystems. 2004, 75: 3-14. 10.1016/j.biosystems.2004.03.009.

    Article  PubMed  Google Scholar 

  58. Takahashi K, Ishikawa N, Sadamoto Y, Sasamoto H, Ohta S, Shiozawa A, Miyoshi F, Naito Y, Nakayama Y, Tomita M: E-Cell 2: multi-platform E-Cell simulation system. Bioinformatics. 2003, 19: 1727-1729. 10.1093/bioinformatics/btg221.

    Article  CAS  PubMed  Google Scholar 

  59. Pries A, Secomb T, Gaehtgens P, Gross J: Blood flow in microvascular networks. Experiments and simulation. Circ Res. 1990, 67: 826-834.

    Article  CAS  PubMed  Google Scholar 

  60. Ji JW, Tsoukias NM, Goldman D, Popel AS: A computational model of oxygen transport in skeletal muscle for sprouting and splitting modes of angiogenesis. J Theor Biol. 2006, 241: 94-108. 10.1016/j.jtbi.2005.11.019.

    Article  PubMed  Google Scholar 

  61. Ji JW, Mac Gabhann F, Popel AS: Skeletal muscle VEGF gradients in peripheral arterial disease: simulations of rest and exercise. Am J Physiol Heart Circ Physiol. 2007, 293: H3740-3749. 10.1152/ajpheart.00009.2007.

    Article  CAS  PubMed  Google Scholar 

  62. Handschin C, Spiegelman BM: The role of exercise and PGC1alpha in inflammation and chronic disease. Nature. 2008, 454: 463-469. 10.1038/nature07206.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  63. Puigserver P, Wu Z, Park CW, Graves R, Wright M, Spiegelman BM: A cold-inducible coactivator of nuclear receptors linked to adaptive thermogenesis. Cell. 1998, 92: 829-839. 10.1016/S0092-8674(00)81410-5.

    Article  CAS  PubMed  Google Scholar 

  64. Knutti D, Kaul A, Kralli A: A tissue-specific coactivator of steroid receptors, identified in a functional genetic screen. Mol Cell Biol. 2000, 20: 2411-2422. 10.1128/MCB.20.7.2411-2422.2000.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  65. Arany Z, Foo SY, Ma Y, Ruas JL, Bommi-Reddy A, Girnun G, Cooper M, Laznik D, Chinsomboon J, Rangwala SM: HIF-independent regulation of VEGF and angiogenesis by the transcriptional coactivator PGC-1alpha. Nature. 2008, 451: 1008-1012. 10.1038/nature06613.

    Article  CAS  PubMed  Google Scholar 

  66. Leick L, Hellsten Y, Fentz J, Lyngby SS, Wojtaszewski JF, Hidalgo J, Pilegaard H: PGC-1alpha mediates exercise-induced skeletal muscle VEGF expression in mice. Am J Physiol Endocrinol Metab. 2009, 297: E92-103. 10.1152/ajpendo.00076.2009.

    Article  CAS  PubMed  Google Scholar 

  67. O'Hagan KA, Cocchiglia S, Zhdanov AV, Tambuwala MM, Cummins EP, Monfared M, Agbor TA, Garvey JF, Papkovsky DB, Taylor CT, Allan BB: PGC-1alpha is coupled to HIF-1alpha-dependent gene expression by increasing mitochondrial oxygen consumption in skeletal muscle cells. Proc Natl Acad Sci USA. 2009, 106: 2188-2193.

    Article  PubMed Central  PubMed  Google Scholar 

  68. Roy TK, Popel AS: Theoretical predictions of end-capillary PO2 in muscles of athletic and nonathletic animals at VO2max. Am J Physiol. 1996, 271: H721-737.

    CAS  PubMed  Google Scholar 

  69. Breen EC, Johnson EC, Wagner H, Tseng HM, Sung LA, Wagner PD: Angiogenic growth factor mRNA responses in muscle to a single bout of exercise. J Appl Physiol. 1996, 81: 355-361.

    CAS  PubMed  Google Scholar 

  70. Kivela R, Silvennoinen M, Lehti M, Jalava S, Vihko V, Kainulainen H: Exercise-induced expression of angiogenic growth factors in skeletal muscle and in capillaries of healthy and diabetic mice. Cardiovasc Diabetol. 2008, 7: 10.1186/1475-2840-7-13.

    Google Scholar 

  71. Ellsworth ML, Popel AS, Pittman RN: Assessment and impact of heterogeneities of convective oxygen transport parameters in capillaries of striated muscle: experimental and theoretical. Microvasc Res. 1988, 35: 341-362. 10.1016/0026-2862(88)90089-1.

    Article  CAS  PubMed  Google Scholar 

  72. Gustafsson T, Ameln H, Fischer H, Sundberg CJ, Timmons JA, Jansson E: VEGF-A splice variants and related receptor expression in human skeletal muscle following submaximal exercise. J Appl Physiol. 2005, 98: 2137-2146. 10.1152/japplphysiol.01402.2004.

    Article  CAS  PubMed  Google Scholar 

  73. Jensen L, Pilegaard H, Neufer PD, Hellsten Y: Effect of acute exercise and exercise training on VEGF splice variants in human skeletal muscle. Am J Physiol Regul Integr Comp Physiol. 2004, 287: R397-402. 10.1152/ajpregu.00071.2004.

    Article  CAS  PubMed  Google Scholar 

  74. Gavin TP, Drew JL, Kubik CJ, Pofahl WE, Hickner RC: Acute resistance exercise increases skeletal muscle angiogenic growth factor expression. Acta Physiol (Oxf). 2007, 191: 139-146. 10.1111/j.1748-1716.2007.01723.x.

    Article  CAS  Google Scholar 

  75. Egginton S, Zhou AL, Brown MD, Hudlicka O: Unorthodox angiogenesis in skeletal muscle. Cardiovasc Res. 2001, 49: 634-646. 10.1016/S0008-6363(00)00282-0.

    Article  CAS  PubMed  Google Scholar 

  76. Armstrong RB, Laughlin MH: Blood Flows within and among Rat Muscles as a Function of Time during High-Speed Treadmill Exercise. J Physiol Sci. 1983, 344: 189-208.

    Article  CAS  Google Scholar 

  77. Armstrong RB: Magnitude and distribution of muscle blood flow in conscious animals during locomotory exercise. Med Sci Sports Exerc. 1988, 20: S119-123. 10.1249/00005768-198810001-00005.

    Article  CAS  PubMed  Google Scholar 

  78. Liu G, Marathe DD, Matta KL, Neelamegham S: Systems-level modeling of cellular glycosylation reaction networks: O-linked glycan formation on natural selectin ligands. Bioinformatics. 2008, 24: 2740-2747. 10.1093/bioinformatics/btn515.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  79. Izaguirre JA, Chaturvedi R, Huang C, Cickovski T, Coffland J, Thomas G, Forgacs G, Alber M, Hentschel G, Newman SA, Glazier JA: CompuCell, a multi-model framework for simulation of morphogenesis. Bioinformatics. 2004, 20: 1129-1137. 10.1093/bioinformatics/bth050.

    Article  CAS  PubMed  Google Scholar 

  80. Mitha F, Lucas TA, Feng F, Kepler TB, Chan C: The Multiscale Systems Immunology project: software for cell-based immunological simulation. Source Code Biol Med. 2008, 3: 6-10.1186/1751-0473-3-6.

    Article  PubMed Central  PubMed  Google Scholar 

  81. Moraru II, Schaff JC, Slepchenko BM, Blinov ML, Morgan F, Lakshminarayana A, Gao F, Li Y, Loew LM: Virtual cell modelling and simulation software environment. Systems Biology, IET. 2008, 2: 352-362. 10.1049/iet-syb:20080102.

    Article  CAS  Google Scholar 

  82. Qutub AA, Popel AS: A computational model of intracellular oxygen sensing by hypoxia-inducible factor HIF1 alpha. J Cell Sci. 2006, 119: 3467-3480. 10.1242/jcs.03087.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  83. Asai Y, Suzuki Y, Kido Y, Oka H, Heien E, Nakanishi M, Urai T, Hagihara K, Kurachi Y, Nomura T: Specifications of insilicoML 1.0: a multilevel biophysical model description language. J Physiol Sci. 2008, 58: 447-458. 10.2170/physiolsci.RP013308.

    Article  PubMed  Google Scholar 

  84. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT: Gene Ontology: tool for the unification of biology. Nature Genetics. 2000, 25: 25-10.1038/75556. Nature Publishing Group; 25

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  85. Smith B, Ashburner M, Rosse C, Bard J, Bug W, Ceusters W, Goldberg LJ, Eilbeck K, Ireland A, Mungall CJ: The OBO Foundry: coordinated evolution of ontologies to support biomedical data integration. Nat Biotechnol. 2007, 25: 1251-1255. 10.1038/nbt1346.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  86. Garny A, Cooper J, Hunter PJ: Toward a VPH/Physiome ToolKit. Wiley Interdiscip Rev Sys Biol Med. 2010, 2: 134-147.

    Article  Google Scholar 

  87. Mac Gabhann F, Qutub AA, Annex BH, Popel AS: Systems biology of pro-angiogenic therapies targeting the VEGF system. Wiley Interdiscip Rev Syst Biol Med. 2010, 2: 694-707.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

Download references


The authors thank the Imaging Science Center at Johns Hopkins University for use of the computer cluster; Anthony Kolasny, IT Architect, for technical assistance; Nicholas Tan for his contribution to the PGC1α model; Dr. Corban Rivera for the discussion. This work was supported by the National Institutes of Health grants R01 HL101200, R33 HL087351, R01 HL079653 and R01 CA138264 (ASP), and R00 HL093219 (FMG).

Author information

Authors and Affiliations


Corresponding author

Correspondence to Gang Liu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

Conceived and designed the simulations: GL AAQ PV FMG ASP; Performed the simulations: GL; Analyzed the data: GL AAQ PV FMG ASP; Wrote the paper: GL AAQ PV FMG ASP. All authors read and approved the final manuscript.

Electronic supplementary material

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Liu, G., Qutub, A.A., Vempati, P. et al. Module-based multiscale simulation of angiogenesis in skeletal muscle. Theor Biol Med Model 8, 6 (2011).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: