Mathematical Modelling of Tumour Angiogenesis and the Action of Angiostatin as a Protease Inhibitor *

Tumour angiogenesis is the process whereby a capillary network is formed from a pre-existing vasculature in response to tumour secreted growth factors (TAF). The capillary network is largely composed of migrating endothelial cells which organise themselves into dendritic structures. In this paper we model angiogenesis via the theory of reinforced random walks, whereby the chemotactic response of the endothelial cells to TAF and their haptotactic response to the matrix macromolecule fibronectin is accomplished through transition probability rate functions. Tumour secreted growth and inhibitory factors are modelled on the basis of Michaelis–Menten kinetics. Particular attention is focussed on the action of anti-angiogenic agents (i.e. angiostatins). That is as a mechanism whereby angiostatin acts as a protease inhibitor. Numerical simulations yield results, which are in good agreement with the experimental observations obtained by Folkman and his co-workers in their classical rabbit eye cornea experiments. The model offers a theoretical understanding of how some angiostatins work to inhibit tumour growth.


INTRODUCTION
Angiogenesis is a morphogenic process whereby new blood vessels are induced to grow out of a pre-existing vasculature. It is fundamental to the formation of new blood vessels during embryonic development and contributes to the maintenance of tissue functionality in the adult. An important example of the latter is placental growth. Angiogenesis is also an important feature of various pathological processes such as wound healing and cancer progression.
Capillaries, the main micro-vessels involved in tumour angiogenesis, are composed of three components: (a) the basement membrane, which is a complex extra cellular matrix (ECM) encircling and supporting the cellular components, (b) the endothelial cells (EC), which form a mono-layer of flattened and extended cells lining the lumen and resting on the basement membrane and (c) pericytes which form a periendothelial cellular network embedded within the basement membrane.
The first event of tumour-induced angiogenesis is triggered by the secretion of a number of chemicals, collectively called Tumour Angiogenesis Factors (TAF) (Folkman and Haudenschild, 1980;Folkman and Klagsbrun, 1987;Folkman, 1985;1995) from a colony of cancerous cells of a solid tumour. These factors diffuse through the tissue creating a chemical gradient which eventually reaches neighbouring capillaries and other blood vessels. In response to TAF the EC in nearby capillaries appear to thicken and produce a proteolytic enzyme (protease), which in turn degrades the basement membrane. In further response to the TAF the normally smooth endothelial cell surface begins to develop pseudo podia which penetrate the weakened basement membrane. Capillary sprouts are formed by the accumulation of EC from the parent vessel. The sprouts grow in length, proliferate, form loops leading to micro-circulation of blood (i.e. anastomoses) and branch successively. The resulting capillary network continues to progress through the tissue ECM forming a micro-vasculature, and eventually invades the tumour colony leading to the rapid growth and metastasis. The means of progression through the tissue ECM is via chemotaxis and haptotaxis. Chemotaxis is the response of EC to chemical gradients set up by the TAF. A major component of the tissue ECM is fibronectin. It has been verified experimentally (Bowersox and Sorgente, 1982) that fibronectin stimulates directional migration of EC by establishing an adhesive gradient via haptotaxis.
Mathematical modelling of the complex processes involved in tumour angiogenesis has proceeded along two main fronts. Firstly, there is the continuum approach based on differential equations developed from the principle of mass conservation and chemical flux considerations. This allows for a description of macroscopic events such as evolution of EC density and cell kinetic and migration characteristics as a result of chemotaxis and haptotaxis. The second avenue of investigation has been via discrete or cellular automata models in order to gain insight into the microscopic world of capillary network formation and to capture observations obtained from experiments in vivo and in vitro. For recent surveys of both these approaches we refer to Stokes and Lauffenburger (1991), Sleeman (1996) and Chaplain and Anderson (1999).
In the case of continuum models we cite the recent work of Levine et al. (2000;2001a) in which a new approach to the initiation of angiogenesis is based on the idea of reinforced random walks, (Davis, 1990;Levine and Sleeman, 1997;Othmer and Stevens, 1997) a key ingredient of the modelling to be pursued in this paper. At the microscopic level the behaviour of individual cells is fundamental to the understanding of capillary formation; particularly that of branching and anastomoses (Stokes and Lauffenburger, 1991;Anderson and Chaplain, 1998 and the references given there). Therefore, there is the challenging problem of understanding the "bridge" between micro and macro-cellular events. In the recent work of Sleeman and Wallis (2000), it has been demonstrated that a possible key to establishing this connection is via the modelling of "transition rates" as exemplified in the work of Davis, (1990), Othmer -Stevens (1997) and Levine et al. (2000;2001a,b).
Angiostatin a proteolytic fragment of the protein plasminogen is known to be a potent inhibitor of tumour angiogenesis. Indeed several angiogenesis inhibitor proteins exist as internal fragments of larger proteins (e.g. Prolactin) suggesting that normal angiogenesis inhibitors may be stored within larger proteins. When the body needs to stop normal angiogenesis, after wounding healing or ovulation for example, these natural inhibitors may be available for immediate use by breaking down the larger proteins (Hanahan and Folkman, 1996;Jendraschak and Sage 1996;Falcone et al., 1998;Gorski et al., 1998;O'Reilly, et al., 1999). It has also been discovered that these naturally occurring angiogenesis inhibitors are normally controlled by the tumour suppressor gene p. 53, which has been implicated in various cancer types.
These important studies has lead to trials of antiangiogenic drugs as TNP-470 as a powerful new anticancer therapy. A primary goal of this paper is to model anti-angiogenic mechanisms and to attempt to understand how tumour angiogenesis may be inhibited.
The paper is organized as follows: In the first section, we describe the Boyden Chamber assay which is the basis for the modelling procedure. The chemistry of capillary sprout formation is considered in the second section while in the third section, we discuss chemical models for the action of angiostatin. Using the laws of mass action and Michaelis -Menten kinetics (third and fourth section) the chemistry is reformulated as mathematical models of chemical transport both within the capillary as well as the ECM.
Endothelial cell transport, using the idea of reinforced random walks, is considered in fifth and sixth section. This completes the mathematical modelling.
In the seventh section, we discuss the results of a number of numerical experiments.
We conclude with a brief analysis of the results of the paper together with a discussion of possible therapeutic strategies which could be pursued in achieving angiogenic inhibition.

BOYDEN CHAMBER ASSAY
Throughout this paper we base our modelling procedure on the simplified configuration shown in Fig. 1. That is in the x -y plane we envisage a capillary segment to be located along the x-axis on the interval ½0; L with a tumour colony, or source of TAF, located l units above the x-axis at y ¼ l: The ECM is modelled as a porous medium. Alternatively, as in the Boyden Chamber Assay, the ECM may be a matragel bed in which case the ECM can be modelled as a viscoelastic medium. Such a configuration has been considered by several authors. See the review article of Tranqui and Tracqui (2000) and Holmes and Sleeman (2000).
Time dependent quantities in the ECM will be represented by capital letters, e.g. Gðx; y; tÞ while quantities in the capillary will be denoted by lower case letters e.g. gðx; tÞ: It is important to note that in general Gðx; 0; tÞgðx; tÞ:

BIOCHEMICAL KINETICS OF SPROUT INITIATION
We begin by considering the initiation of elementary sprouts from the capillary wall. To better understand how tumour generated angiogenesis factors act on the EC we consider each EC has a certain number of receptors to which the angiogenic factor (ligand) binds. The bound molecules (intermediates) then stimulate the cell to produce a proteolytic enzyme and form new receptors.
Let V T be the angiogenic factor (TAF) released by the tumour colony. Then we model the process as follows where R E denotes some receptor site on the EC surface membrane, [R E V T ] is the intermediate and C T is a proteolytic enzyme produced as a consequence of this reaction. The proteolytic enzyme C T now acts to degrade fibronectin (F ) via the Michaelis -Menten catalytic reaction.
The Michealis-Menten condition is presumed to be in force here. This amounts to the statement that This reaction reflects the fact that the enzyme degrades fibronectin, converting it into products F 0 be means of catalysis. It is to be emphasized that while the chemotactic agent V T diffuses through the ECM from the tumour, it is converted almost immediately into receptor complexes upon arrival at the capillary wall via the above reactions. Therefore we expect that the diffusion of these proteins along the capillary lumen is negligible in comparison to their conversion into receptor complex (Levine et al., 2001b). However, the diffusion of these proteins cannot be TUMOUR ANGIOGENESIS AND ANGIOSTATIN neglected in the ECM. It is a major component of the transport mechanism for these proteins to and from the tumour.

MODELS FOR THE ACTION OF ANGIOSTATINS
There are several possible ways in which angiostatic agents could inhibit angiogenesis. Here we discuss one possibility which has been verified in the experimental literature (Stack et al., 1999). That is we consider angiostatin A as a direct inhibitor of protease viz.
where C A refers only to those molecules of proteolytic enzyme C T that are available for the degradation of fibronectin. C I denotes the proteolytic enzyme molecules which are inhibited by the angiostatin A from functioning as a catalyst for fibronectin degradation. In terms of where n e is the equilibrium constant for this step.
In general the reaction (3) is essentially complete. i.e. the equilibrium constant, n e , is quite large. However, the one example we have for this mechanism is to be found in the experimental work of Stack et al. (1999), where plasminogen activation is shown to be inhibited by angiostatin. (The tissue plasminogen activator (tPA) is produced in response to an angiogenic growth factor such as vascular endothelial cell growth factor (VEGF). Then tPA binds to plasminogen with the resultant product being plasmin (Pm). In the work of Stack et al., (1999), it is proposed that angiostatin binds directly to the intermediate [tPA -P g ] complex to inhibit the production of active protease. The angiostatin here is a fragment of plasminogen with a molecular weight of about 38 kDa. The literature value in Stack et al. (1999) given for n 21 e is of the order of 1 mM which is not very large.
(We are currently preparing an article in which this entire process is modelled using the ideas elucidated below. Here, we just employ the known equilibrium constant and the mechanism for angiostatin action.) In the mechanism (3), the critical equation for the concentration of active enzyme is Time course for EC propagation in the mother capillary.

CHEMICAL TRANSPORT EQUATIONS IN THE CAPILLARY
From the system of chemical reaction rate equations in "Biochemical kinetics of sprout initiation" and "Models for the action of angiostatins" section, we construct a mathematical model governing chemical transport in the capillary using the laws of mass-action and a careful use of Michaelis -Menten kinetics. We do not detail the derivations here but refer the reader to Levine et al. (2000;2001a,b) for a full account of the development.
Introduce the notation.
The chemical transport equations are given by This system models chemical transport when angiostatin is considered to be a direct inhibitor of protease as governed by Eq. (3). Here T rel is the relaxation time for the angiostatin decay.
In the system (5), n e is the equilibrium constant for Eq. (3) while l 1 ¼ k 1 d e h 0 ; l 2 ¼ k 3 ; n 1 ¼ k 1 =ðk 21 þ k 2 Þ and n 2 ¼ k 3 =k 4 : d e is the average number of receptors on an endothelial cell which bind to protease and h 0 is the EC density in a normal capillary. In addition the forcing term v r ðx; tÞ allows for the rate at which TAF  is produced by the tumour colony. Natural decay of the enzyme c is modelled by 2mc. It is known that EC produce fibronectin and this is accounted for by the introduction of the logistic growth term in Eq. (5) where f 0 is the maximum value of f. For therapeutic purposes any angiostatic agent must be supplied at a rate a r ðx; tÞ sufficient to "neutralise" the proteolytic enzyme until the tumour has been rendered inactive. A reasonable model for a r might be a constant function, representing a constant rate of supply of antiangiogenic drug.

CHEMICAL TRANSPORT EQUATIONS IN THE ECM
Within the ECM much the same chemistry applies as in the capillary. However, there are significant modifications necessary.
First of all we assume that the background fibronectin production is in considerable excess of that produced by the endothelial cells, so that the logistic growth term in the analogous fibronectin dynamics Eq. (5) is independent of Nðx; y; tÞ: Secondly, we assume that in so far as growth factor and angiostatin are concerned, the ECM is a porous medium through which these chemicals diffuse.
Furthermore we do not assume that the diffusivities for these species, namely D V , D A are constant. Thirdly, we allow for inhomogeneous diffusion of growth factor and for angiostatin. That is we assume that molecules of either type visualise the ECM as a porous medium but with variable diffusivities D V , D A . Furthermore, we assume that the porosity index "m" is the same for both species. This is based on the fact that these two proteins are about the same size. This is not the only way in which the ECM can be modelled. For instance it may be considered to be a visco-elastic medium as considered in Tranqui and Tracqui (2000) and Holmes and Sleeman (2000).
It is known that proteases must be bound to cell surfaces in order to effect cell migrations in the ECM (Alberts et al., 1994). Hence no corresponding diffusion terms need to be included in the rate law for the proteolytic enzyme.
Finally we need to include the diffusivity of fibronectin within the ECM. Generally fibronectin diffuses slowly. However, it is reasonable to assume that its diffusion is curvature dependent. In our two dimensional configuration we model this via mean curvature. ‡ ‡ If z ¼ fðx; yÞ then the curvature, for fixed z, of the level line is given by k ¼ 7 7f

j7fj
: The inclusion of this term in the fibronectin equation is motivated by an idea taken from the dendritic crystal growth. There, growth of dendrites occurs only at the tip where the local ratio of surface area to volume is greatest. From the above discussion, and keeping in mind the notational convention described in "Boyden chamber assay section", the kinetics within the ECM are taken to be If k , 0, the growth rate for fibronectin (F t ) is diminished whereas if k . 0 the growth rate is increased.
In the first equation of (6) we have included a source term V r ðx; y; tÞ to account for the situation in which growth factor is generated within the ECM. Furthermore we have included the source term a r ðx; tÞ ¼ ð1 2 F=F 0 Þ in the fourth equation of (6). This term allows us to introduce angiostatin in every region of the capillary network for which fibronectin density is below the background value F 0 in the ECM, at a rate which is proportional to the fibronectin deficit in the ECM.

CELL TRANSPORT IN THE CAPILLARY AND ECM
To model the transport of EC both within the capillary in response to the angiogenic growth factor and in the ECM we use, as in Levine et al. (2000;2001a,b) a continuum form of the idea of reinforced random walks developed by Othmer and Stevens (Levine and Sleeman, 1997;Othmer and Stevens, 1997). That is, if h denotes the EC in the capillary wall then In this equation t is called the transition probability function and is modelled to reflect the fact that the EC respond both chemotactically and haptotactically, respectively, to gradients in proteolytic enzyme c and fibronectin f. The interpretation of Eq. (7) is that the rate of change of EC is balanced by random diffusion and chemical sensitivity to c and f. The modelling of t is a challenging problem. However from biological considerations (Levine et al., 2000;2001a,b) we take in our simulations where a i ; b i ; i ¼ 1; 2 are empirical constants with 0 , a 1 ! 1 , a 2 ; b 1 . 1 @ b 2 . 0: These choices of parameter ranges reflect the fact the EC aggregate towards high concentrations of c a and to regions where f is most degraded. For full justification of the modelling of t we refer to Levine et al. (2000;2001a,b).
To model the transport of EC within the ECM we need to take account of both the proliferation and apoptosis or In this equation the first two terms on the right hand side are the two dimensional analogues of the corresponding terms on the right of Eq. (9). The third term incorporates several important issues. First we have included a "logistic" proliferation rate of EC and a linear "death rate". It is known that proliferation of EC occurs only near the tip of growing capillaries and to account for this we have included a factor Q(k ) which is small when tip curvature k is small and large when the tip curvature is large. Furthermore it is also known that cell proliferation respond to active protease C A in a biphasic manner and this has been modelled in Eq. (9) by the inclusion of the function G(C A ) shown in Fig. 2.
Finally proliferation only occurs if the active protease C A is above a certain threshold C A,0 . This is modelled in Eq. (9) by the inclusion of the "Switch" HðC A 2 C A;0 Þ where H is the Heaviside Step function defined by ( To complete this rather extensive mathematical model based on the chemical transport equations of sections "Models for the action of angiostatins" and "Chemical transport equations in the capillary" together with the EC equations of section "Cell transport in the capillary and ECM" we must close the system with appropriate boundary and initial conditions. This is detailed in Levine et al., (2001b) and will be referred to appropriately in the numerical experiments discussed in the "Numerical experiments section".

NUMERICAL EXPERIMENTS
To begin with we adopt the phrase "onset of angiogenesis" to mean that a capillary sprout has been initiated from the pre-existing vasculature. By "onset of vascularization" we mean that a newly formed capillary has just reached the tumour colony. In vivo one expects that tumour stimulated growth rate of new capillaries will depend on several variables. These would include chemical properties of growth factors and the proteases themselves as well as the protein structure of the ECM. The model we have developed above reflects this dependency.
To connect the EC transport equations with the capillary wall and the ECM we assume that the EC in the capillary wall move into the ECM when the fibronectin density f in the capillary wall falls below a certain threshold level f 1 . That is we take where 0 , c 1 # 1 is interpreted as a measure of the fraction of EC lining the lumen that are able to penetrate the ECM. We next assume that the tumour colony is supplying a prescribed flux of VEGF. In our experiments this flux is given by where s is a fixed constant normalized so that the mean flux is v 0 . By lowering the threshold value of f 1 or by decreasing the dosage rate n 0 will have the effect of slowing the rate of onset of angiogenesis. Also by changing the percentage of EC's c 1 will effect capillary tip growth rate. The coefficient Q(k ) in Eq. (9) is a curvature sensitivity factor and is taken to be This choice is made not only because we want the curvature sensitivity to be dimensionless but also because we need to control the sensitivity to proliferation. Note maximum sensitivity is Q max ¼ 21=1: Finding good estimates for the constants used in the model is not an easy task. Fortunately with the help of Mr John E. Hinrichsen we were able to locate some constants and obtain order estimates for some others. The entries used in Table A1 of the appendix include those found as well as others used in our computations. In Levine et al. (2001b) the reader may find extensive references for the entries in Table A1.
In the first numerical experiment we implanted a "tumour" 25 microns from a pre-existing capillary. In Fig. 3, we see the advance of EC across the ECM while in Fig. 4 we observe the EC aggregate along the capillary wall. Notice the peak of EC density near the tip in Fig. 3 and the bimodal distribution of EC density in Fig. 4. In Fig. 5 the degradation of fibronectin with time forming a channel through the ECM is seen. In Fig. 6, one sees the degradation of the capillary wall.
In Table I we give the travel times to various points in the ECM after tumour implantation. The tip speeds in the third column are quite close to the tip speeds observed in Folkman, et al. The extrapolated times of one and two millimeter implant distances are also consistent with the travel times they observed. Furthermore, they observed that the tip speed increases with tip approach to the tumour. This is also the case for our stimulated capillary tips.
Next we introduce angiostatin at the rate a r ðx; tÞ ¼ A r Hðt 2 T iv Þ: Here A r is the rate at which angiostatin is being supplied and T iv is elapsed time since the tumour began to leak growth factor into the ECM. In our computations T iv ¼ 4:1 h: We chose this time because for our choice n e , the reaction (3) is not complete. This means that one either has to swamp the mother capillary with very high concentrations of inhibitor to bind the active protease in order to prevent tumour vascularisation or else introduce the inhibitor relatively early in the growth of the daughter capillary. For purposes of illustration, we have opted for the latter. Figure 7 shows the time course of EC propagation in the ECM before and after inhibitor introduction. Figure 8 gives the corresponding time courses for fibronectin in the ECM. Notice that there has been a marked drop in EC density although the closing of the fibronectin channel is somewhat more muted. However, the closing of the capillary at the point where it begins to enter the ECM from the mother capillary can clearly be seen in Fig. 8d. In Fig. 9, we plot time courses within the mother capillary, of EC, firbronectin, protease and active protease in the presence of angiostatin. The dashed lines correspond to a time of 5.1 h while the solid lines correspond to the time of 4.1 h after the experiment has begun. Notice that the concentration of protease has fallen by about fifty percent while there is almost no active protease in the mother capillary. The opening (Fig. 9b) in the fibronectin is also closing.
In Figs. 10 and 11, we see the effect of the introduction of angiostatin on the time courses of VEGF, inactive protease and active protease.
We next turn to an analysis of the dependency of endothelial cell movement across the ECM to capillary tip curvature k based on the curvature sensitivity factor, e. To begin with we denote by T a for 0 # a # 1 the time, in hours, from the onset of "tumour" implantation to the time the capillary tip has moved a fraction a of the distance l from the capillary to the tumour colony. That is T a is the time taken for some EC's to have moved a distance al from the capillary located along the x-axis. The sensitivity of the travel times to the capillary tip sensitivity parameter e is recorded in Table II.
Notice that as e increases, cell travel time across the ECM slows down. This is to be expected since curvature sensitivity decreases with increasing e. Furthermore this will effect the rate of enzyme degradation of fibronectin due to the fact that the rate of growth of protease is proportional to EC concentration. Consequently, if protease concentration does not increase, the rate of fibronectin degradation will decrease. For example, when 1 ¼ 5; V ¼ 0:24 mm=d: To investigate how the travel time across the "ECM" to the capillary and back to the tumour depends on the distance from the tumour to the ECM, the number l, the sensitivity is fixed ð1 ¼ 1:0Þ and l is varied.
The results are shown in Table III.
The table clearly shows that as the distance to the tumour increases the travel time, as expected, increases but the mean tip speed decreases.

CONCLUSIONS AND DISCUSSION
In this paper a mathematical model of tumour angiogenesis has been described on the basis of the chemical kinetics governing the chemistry of angiogenic factors, protease activity as well as of fibronectin degradation and angiostatin interference in protease activity. The motion of endothelial cells both within the ECM and in the capillary wall is governed by a continuum limit of a reinforced random walk process.
Of particular importance is the angiostatic agent considered here as a direct inhibitor of protease as indicated in the experimental literature.
Using known literature values to the maximum extent possible, our numerical simulations give results that are in reasonable qualitative and quantitative agreement with the results of the rabbit eye cornea experiments of Folkman et al. In particular travel times of EC across the ECM towards the tumour are of the same order of magnitude as found in in vitro studies. The diameter of the growing (numerical) capillary is of the order of 6 -10 microns. It is also shown that capillary growth is sensitively dependent on tip curvature and hence on cell proliferation. The role of angiostatin as a direct protease inhibitor is included and shows that if angiostatin is administered soon after the onset of angiogenesis then growth of capillaries is greatly reduced and fibronectin degradation severely restricted Table A1. Mean tip speed is computed as n ¼ 0:017=ðT 1 2 T 0 Þ in mm/d. ðv 0 ¼ 3:0 mM mm=h; f 1 ¼ 0:40; l ¼ 17 mmmÞ: