Effects of Matrix Metalloproteinase on Tumour Growth and Morphology via Haptotaxis

A tumour invasion model has been developed to examine the influence of diffusible matrix metalloproteinases (MMPs) and the mediated proteolysis of non-diffusible, membrane type MMPs (MT-MMPs) on tumour growth and morphology via haptotaxis. Our results are the first to explore the influence of localized degradation of extracellular matrix (ECM) by MT-MMPs on the morphology of tumours growing in nutrient-rich and nutrient-poor microenvironments. Two-dimensional numerical simulation, using a level-set-based tumour host interface-capturing method, reveals that haptotaxis due to ECM degradation by MMP causes greater instability than with ECM degradation by MT-MMP in low-nutrient environments, even at low proliferation rates; whereas the resulting morphologies are similar for high apoptosis rates. Our simulation results show that while haptotaxis leads to completely different tumour growth rates and morphologies depending on proliferation and apoptosis rates in low-nutrient environment, there are no significant variations when we compare the haptotaxis due to ECM degradation by MMP and MT-MMP, except for low proliferation rates. Focusing on the differences between MMP and MT-MMP mediated effects; our study has important implications in MMP-target validation and MMP-inhibitor-drug development for anti-cancer clinical trials. Effects of Matrix Metalloproteinase on Tumour Growth and Morphology via Haptotaxis


Introduction
Tumour invasion is driven by proliferation and migration into the surrounding tissue. Invasive cells have different characteristics than those of normal cells in that they are less adhesive and more mobile, mitotic, and metabolically active than normal cells [1][2][3][4]. Among these characteristics, cell motility is a crucial aspect of tumour invasion as cells with motile capabilities can access new nutrient sources and infiltrate the surrounding tissue which is essential for metastasis. Cell motion is commonly described by both random motion of the cells and directed motion stimulated by chemical (nutrient or other growth factors) or structural (ECM) gradients. The gradients are sensed by the migrating cells and influence the direction of movement. Chemo taxis is the phenomenon by which the movement of cells is directed in response to a soluble extracellular chemical gradient; i.e., of nutrient or oxygen or other growth factors. Cells produce matrix degrading enzymes (MDE); e.g., matrix metalloproteinase (MMPs) and urokinase plasminogen activators that degrade the ECM locally to make room for migration as cells cannot move into regions of the tissue which are too dense. The directed movement of cells in response to gradients of fixed or bound (non-diffusible) chemicals, such as the ECM, is known as haptotaxis.
Most MMPs are secreted as soluble enzymes but six of them are membrane-type MMPs (MT-MMPs) that are associated with the cell membrane [5]. Soluble MMPs, such as MMP-2, are secreted as latent pro-enzymes (proMMPs) and the activation of proMMPs often involves MT-MMPs and tissue inhibitors of matrix metalloproteinases (TIMPs). In addition to activating other MMPs, MT-MMPs degrade a number of ECM macromolecules [6][7][8]. MT1-MMP (MMP14) has emerged as an important collagenase that cancer cells use to degrade and invade in a collagen-rich environment [9,10]. While poorly invasive breast adenocarcinoma cells undergo apoptosis when confronted with a collagen-rich environment, the production of MT1-MMP endows these cells with the capacity to escape from collagen-induced apoptosis [11]. The localized activity of MT1-MMP helps tumour cells to overcome higher collagen density found in some peritumoural stroma [12]. Sabeh, et al. [13] have demonstrated that when cancer cells are faced with structural barriers created in reconstituted gels by covalently cross-linked fibrils of type I collagen, or that exist in the stromal environment of the mammary gland, invasion is dependent on MT1-MMP-mediated proteolysis.
Experimental evidence suggests that hypoxia (shortage of oxygen/ nutrient) can contribute importantly to instability and increase the invasive behaviour of a tumour. Instability (i.e. non-circular shape) provides a mechanism for tumour invasion that does not require development of a neo-vasculature to supply essential nutrients. Poplawski et al. and Poplawski et al. [14,15] demonstrated that the development of instability depends primarily on the diffusional limitation parameter (ratio of tumour growth rate to diffusion rate of nutrients) while the morphological details depend on cell-cell adhesion. Their works show that the lack of competition for nutrients (i.e. high nutrient environment) promotes spherical, non-invasive tumours and low concentrations of nutrients in the environment (which cause tumour-cell competition), or cells with a very high substrateconsumption rate generate a fingering instability and irregular, invasive tumours. Their results agree with the in vitro and in vivo experiments conducted in and with other tumour model predictions showing invasive fingered morphology if the nutrient supply is too small. Similar behaviour has been reported when haptotaxis is considered along with the impact of nutrient availability on tumour morphology [16][17][18][19][20][21][22][23][24][25][26][27].
Apoptosis (programmed cell death) serves as a natural barrier to cancer development and therefore, the ability of cancer cells to evade apoptosis is an essential hallmark of cancer [2,28,29]. Tumour cells evolve a variety of strategies to limit or circumvent apoptosis and often show weaker response to hypoxia-induced apoptosis [30,31]. Research has revealed how apoptosis is attenuated in some tumours that succeed in progressing to high-grade malignancy states and resistance to therapy [28,29]. Anderson et al. [32] demonstrated that a high necrosis threshold (below which cell death occurs) leads to more fingers in low nutrient environment. Previous studies (reviewed in) also suggests that defects along apoptotic pathways play a pivotal role in cancer development and many novel treatment strategies targeting apoptosis may be used [28].
Deakin et al. [12] presented a continuum model of cancer cell invasion to examine the differences between the invasion mediated by the localized degradation by MT1-MMP and the more extensive degradation by the soluble and diffusible MMP-2 or MMP-9. The model included cancer cell proliferation (with competition for space with the ECM), random diffusion, a hypotactic response to ECM gradients and the production (by cancer cells) and the natural decay of the MDEs. The results showed that invasion through soluble MMPs produces a greater invasion depth and greater amount of degraded ECM than those through MT-MMPs. Watanabe et al. [33] studied the dynamic nature of MT1-MMP at invadopodia and demonstrated the importance of its transient peak (and subsequent degradation of the ECM) in the activity of MT1-MMP followed by steady state activity. This transient activity was due to the inhibition by TIMP-2, and the steady state activity of MT1-MMP. When the transient activity was forcibly suppressed in computer simulations, the ECM degradation was heavily suppressed, indicating the essential role of this transient activity of MT1-MMP in the ECM degradation.
Other mathematical models that have examined the activation of MMP-2 by MT1-MMP include [34][35][36]. The experiments and simulations conducted by Hoshino et al. [36] established the role of the rapid turnover of MT1-MMP in ECM degradation at invadopodia (protrusions in the cell membrane extended into the ECM). The degree of the reduction in ECM degradation depended on the degree of the reduction in the MT1-MMP turnover. Furthermore, their simulations suggested synergetic contributions of proteolysis activity and the MT1-MMP turnover to ECM degradation because there was a nonlinear and significant reduction in ECM degradation if both factors were reduced simultaneously. But none of these works have considered the effect of haptotaxis due to ECM gradients created by MT-MMP on tumour growth and morphology with the variations of nutrient availability in the tumour microenvironment. In this study we investigate and compare the influences of haptotaxis due to MT-MMP and MMPmediated effects for different proliferation and apoptosis rates in low and high nutrient environments.
Primarily, tumour growth depends on the balance between expansive forces (caused by cell proliferation) and cell-cell adhesion forces (which maintain the tumour's compactness) [37]. Lownutrient tumour growth and morphology are found to be strongly affected by adhesion and haptotaxis. Increasing cell-cell and cell-ECM adhesion of the microenvironment can help to limit the rate of tumour fragmentation and the extent of tissue invasion [24]. There is a combination of the strength of adhesion and haptotaxis in which finger like shapes are observed [38]. Change in the adhesion properties of cells and the disruption of the hypotactic movement result in less extensively disperses invasion fronts, as confirmed in vitro studies [39,40]. Cell proliferation, which depends on the availability of nutrient in the tumour microenvironment, is the predominate factor for tumour growth. Therefore, the effect of MT-MMP on tumour growth and morphology in nutrient-poor and nutrient-rich tissues with varying extents of adhesion must be understood for the design of drugs targeting haptotaxis for controlling malignant tumour growth and invasion.
In the present work, a continuum model of tumour cell invasion is proposed, based on previous continuum models of solid tumour growth to investigate the effect of haptotaxis on tumour growth in nutrient-poor and nutrient-rich tissues [41,42]. Our model accounts for tumour cell migration due to proliferation, cell-cell and cell-ECM adhesion and increasing gradients of nutrient and ECM. We capture the response of adjacent healthy tissue to tumour growth by solving for velocity and pressure both inside and outside the tumour. Our model enables qualitative predictions regarding the invasive ability of tumour cells. Using two-dimensional numerical simulations, we will present the results of a parameter study of quantitative aspects of the tumour progression, such as the size and shape of tumour with a range of biophysical and taxis parameters, and investigate the causal link between these parameters on tumour growth and morphology. Our model focuses solely on the interactions between the tumour cells and the surrounding tissue, so we have not included angiogenesis. The objectives of the present work are to differentiate the effect of haptotaxis for localized and extensive ECM degradation by MT-MMP and MMP, respectively, on tumour growth and morphology in high and lownutrient environments with varying extents of adhesion.

The Mathematical Model
A one millimetre radius spheroidal tumour contains roughly one million cells. In large scale systems where the cancer cell population is on the order of one million or more, continuum models are more suitable than discrete cell models. In a continuum model, the cell species velocity is obtained from the inertia less momentum conservation equation based on Darcy's law, representing the instantaneous equilibrium among forces associated with pressure, cell adhesion, elastic forces, forces exchanged with the ECM leading to haptotaxis and chemo taxis due to gradients of nutrient and other growth factors, and other mechanical effects. The continuum models incorporate the adhesion among cancer cells by a surface tension at the tumour surface [20][21][22]24,37,41,[43][44][45][46][47][48][49][50]. Surface tension also models cell-ECM adhesion and the presence of a membrane, or capsule, formed of ECM macromolecules that may encapsulate tumours in vitro and in vivo. This representation of adhesion is relatively simple and indirect, involving no explicit modelling of cell-cell or cell-ECM contact. In these approaches cell to cell and cell to ECM adhesion are modelled as tumour shape stabilizing mechanical forces. This choice is supported by the experimentally observed presence of surface tension at tissue boundaries [51]. The tumour and healthy tissues are modelled as constant-density fluids. A sharp interface separates the tumour and healthy tissue. Cell-cell and cell-ECM adhesion is modelled using a surface-tension force at the tumour interface. The cell velocity is determined by Darcy's law, and growth occurs due to pressure gradients induced by mitosis (cell proliferation) and haptotaxis. Here we model non-necrotic tumours while accounting for apoptosis (programmed cell death).
We consider a rectangular region as shown in Figure 1 Containing the tumour mass Ω T with boundary Σ T and the non-cancerous tissue Ω H consisting of the ECM, healthy cells, and any other material immediately surrounding the tumour. We choose to focus on the key variables involved in tumour cell growth and invasion; nutrient concentration ρ N , MDE concentration ρ M , ECM density ρ E , cellular velocity field u and solid pressure p. We derive a system of coupled partial differential equations to model tumour growth in the surrounding tissue and present numerical solutions in two spatial dimensions describing the macroscopic dynamics of invasion in different tumour microenvironments due to haptotaxis. The surrounding tissue contains a variety of host cells like fibroblasts, macrophages and blood vessels, all of which have been shown to be important factors in tumour growth [52]. However, these factors are neglected in order to focus on the impact of nutrient concentration and the haptotaxis on tumour growth.

Nutrient transport
The net effect of nutrients and growth-promoting factors are characterized by a single vital nutrient (oxygen) that is required for cell survival and mitosis. The concentration of nutrient ρ N , which diffuses through the ECM and is consumed by the tumour cells, satisfies the following reaction-diffusion equation.
Here, D N is the nutrient diffusivity, δ N is the rate of nutrient uptake by the tumour cells (assumed constant) and I ΩT is a characteristic function equal to 1(0) inside (outside) of the tumour region Ω T . Because nutrient diffusion and uptake occur much more quickly than tumour growth, the quasi-steady assumption applies. We assume that the nutrient and nutrient flux are continuous across the tumour boundary. The two cases that we consider are (i) a low-nutrient environment (e.g., no pre-existing blood vasculature in the computation domain) and (ii) a high-nutrient environment (e.g., the presence of a pre-existing network of nutrient-supplying blood vessels in the domain). For case (ii), nutrient delivery by the blood vasculature and uptake by noncancerous cells are assumed to be in balance outside of Ω H so that the blood vasculature sufficiently delivers nutrient so that ρ N is a constant (ρ N∞ ) at the boundary of the computational domain. So for case (i), we solve the Eq. (1) in the entire domain and for case (ii) we solve it only in the tumour region. We assume that tumour cells uptake nutrient at a greater rate than noncancerous cells and that nutrient uptake is negligible in Ω H [53][54][55]. Furthermore, we assume that there is no nutrient decay in the entire domain.

The MDE and ECM density
MMPs are produced (or activated) by the tumour cells at a constant rate λ M and diffuse freely throughout the tissue. Removal of these enzymes takes place due to natural decay and by deactivation of the enzymes at a constant rate λ DM . The equation governing the evolution of MMP concentration is therefore given by, Because the diffusion coefficient for MMP (D M ) is much smaller than that for oxygen, the full time-dependent diffusion equation is used [30]. As MT-MMP is bound to tumour cells, it moves with the tumour cells with velocity u. So following, (but not considering random cell movement and ECM as physical growth restraint as modelled in) the MT-MMP equation reads [12].
MT-MMPs are produced by the tumour cells at a constant rate λ MT and decreased due to natural decay and by deactivation of the enzymes at a constant rate λ DMT . The reference value ρ M0 /ρ MT0 is the maximum possible MMP concentration, assuming that tumour cells can sense the amount of MDE produced and the production of MDE is a function of ρ M0 /ρ MT0 [42]. The ECM is a complex mixture of macromolecules, but for the sake of simplicity we represent it with a single concentration. Similarly, we characterize MMP and MT-MMP with a single concentration. Although MMPs have multiple roles either promoting or inhibiting tumour progression or metastasis, we only focuses on pro-invasive roles in this paper [56]. In our model, we do not consider MMP-activation process and the effects of TIMP. As cell-ECM interaction is very complex, it is beneficial to investigate each of these effects separately.
Degradation of the ECM by MDEs increases the ability of the tumour to push into the surrounding tissue, both by reducing the mechanical rigidity of the surrounding tissue and by creating extra space for the growing tumour [57]. The ECM is considered as non-motile matter and changes in its distribution are due to its local degradation by MDEs upon contact at a rate λ DE . The ECM density ρ E is given by, For degradation by MT-MMP, ρ M and ρ M0 will be replaced by ρ MT and ρ MT0 , respectively. In the far-field (at the boundary of the computational domain), we use the Dirichlet boundary conditions ρ M =0, ρ MT =0, and ρ E =ρ E∞ (ECM density at far field). The ECM may deform and gets remodelled in response to pressure and insoluble ECM macro-molecules released by the cells, but we have not considered this.

Cellular velocity field
Tumour cells are assumed to move with a single cellular velocity, and cellular motion through the ECM is assumed to be similar to fluid flow in a porous medium. The proliferating tumour cells generate an internal (oncotic) mechanical pressure that also exerts force on the surrounding noncancerous tissue in Ω H . We assume that the motion of the cancer cells is governed by the proliferation pressure and response to nutrient and ECM gradients by chemo taxis and haptotaxis, respectively. Cells thus move with a mass-averaged velocity arising from a generalized Darcy-type constitutive law for velocity from excess forces due to chemo taxis and haptotaxis.
Here p is the proliferation pressure and µ, χ N and χ E are the mobility, chemotaxis and haptotaxis coefficients, respectively. The cellular mobility µ characterizes the overall ability of the tumour tissue to respond to pressure gradients as well as the permeability of the tumour cells. The growth rate of the tumour is characterized by its rate of volume change and is accounted for with the continuity equation. Cell birth and death are in balance in Ω H and so there is no change in the volume of that region. The diffusion of tumour cells is an order of magnitude smaller than advection Zheng et al. [41], so we neglect diffusion of tumour cells. Although the ECM functions directly as a growth restraint for tumour cells, we have not included this role in our model in order to focus on the effect of directional migration (via haptotaxis) on tumour growth and morphology. The haptotaxis coefficients are assumed to be non-negativeand homogeneous in Ω T , but they can be functions of time, space and ECM density. Experiment Zaman et al. [58] has shown that cancer cell motility depends not only on ECM gradients (haptotaxis) but also on ECM density (haptokinesis). However, since our aim here is to focus on gradient driven migration, we effectively choose to ignore the dependence of the cancer cell motility on the ECM density in the absence of ECM gradients. We have also neglected cell movement along ECM fibers (collagen), which is known as contact guidance [59,60].
Using a mass-conservation argument, we get the following equation for the tumour cell density ρ T : Where λ p and λ A are the rates of volume gain and loss due to cellular mitosis and apoptosis respectively. The tumour mass increases through cell proliferation and decreases through apoptosis for a non-necrotic tumour. As cells prefer to stay bound to each other, we assume that the total tumour cell density ρ T is constant, defining what is known as the close-packing density. Due to the sharp interface between the tumour and host tissue, we assume that ρ T = 0 in Ω H in and ρ T = 1 in Ω T so that equation (5) We also assume that in the proliferating region, cell-mitosis is proportional to the amount of nutrient present and that apoptosis may occur. There is no proliferation in the host microenvironment. We assume that the normal velocity is continuous across the tumour boundary.

Mechanical pressure
As the tumour phase, encompassing cells, interstitial fluid, and extracellular matrix, is treated as porous media, no distinction between interstitial fluid hydrostatic pressure and mechanical pressure due to cell-cell and cell ECM interactions is made. When there is no taxis term, continuity of the normal velocity across the tumour boundary dictates that there is no jump in the normal derivative of p across Σ T , i.e., (∂p/∂n) = 0. Assuming a uniform cell-cell and cell-matrix adhesion throughout the tumour, the adhesive pressure can be incorporated as a jump boundary condition at the tumour-host interface Σ T of the Laplace-Young surface-tension type; i.e., [p]= p outer -p inner = κΓ, where κ is the mean curvature of the interface and Γ is a constant adhesion parameter analogous to the surface-tension coefficient. Thus, the adhesive forces are modelled by a curvature boundary condition on the interface Σ T . The non-cancerous tissue in Ω H is assumed to be close enough to the tumour to be affected by the pressure changes within the computational domain. The pressure is assumed to satisfy the homogeneous boundary condition, p=0. Experiment has found that the activity of enzymes is cantered almost exclusively in a narrow region adjacent to the tumour surface where it is balanced by the opposing anti-proteases secreted by nearby healthy tissue [4,61]. Therefore, we take taxis coefficients as piecewise constant terms i.e., in the tumour region they are positive numbers and outside the tumour, they are set to zero. Thus, pressure satisfies the following normal-gradient jump condition. [ Defining the net pressure responsible for tumour growth as p net , such that p net =p-(χ N ρ N +χ E ρ E )/μ Eq. (5) becomes, . net u p µ =− ∇  (9) Combining Eqs. (7) and (8), we obtain With the following jump conditions on pressure and its normal gradient across the tumour boundary.
The cell velocity u can be determined from Eq. (9) once the pressure is known.

Non-dimensionalization
As the nutrient-concentration equation reveals intrinsic diffusionlength and relaxation-time scales, we non-dimensionalize space and time with L and the cell proliferation (mitosis) time scale τ, where The length scale L corresponds to the diffusion-penetration length of oxygen in tissue and therefore characterizes the maximum invasion distance at the early stage of tumour invasion. The characteristic tumour pressure p T =L 2 λ p,max/ is that which results in maximum cell proliferation with a cell speed of unity. The field variables are made dimensionless as specified in Table 1, while non-dimensional parameters are defined in (Tables 2-4). As reference densities, we use the far-field nutrient concentration ρ N∞ , the maximum sustainable MMP and MT-MMP concentration ρ M0 and ρ MT0 , respectively and the far-field ECM density ρ E∞ . The dimensionless governing equations in terms of the nondimensional variables and parameters are given below. We use the same symbols for the non-dimensional variables as those used for the respective dimensional variables defined above. , The following non-dimensional initial and boundary conditions are applied.

Numerical Method
The set of governing equations, Eqs (13 and 14), are solved We solve the Eq. (16) with the semi Lagrangian scheme developed by Aldredge [63] using upwind transient interpolation modelling [64,65]. The scheme provides third-order spatial accuracy and shape preservation. We use the immersed interface method [66] to solve for pressure, which is discontinuous across the tumour-host interface, in accordance with the boundary conditions specified in Eq (15). This requires calculation of normal vectors and curvature at the interface using the level-set function. For accurate normal φat and every curvature time step calculation by we reinitialize by solving the following equation [67,68].
Where φ 0 is the level set function prior to re-initialization and τ is the re-initialization time, reset to and 0 at the beginning of each reinitialization process. We discretize the temporal derivative using a third-order total-variation-diminishing Runge-Kutta method (TVD-RK) and approximate sign(φ 0 ) |grad(φ)| with the fifth-order WENO scheme. However when two interfaces are in close contact, level set functions develop singularities that yield inaccurate normal vectors and curvatures when applying traditional discretization methods. In these cases, we construct the level-set function on a local sub grid where we can accurately calculate the normal vectors and curvature using standard discretization methods [69,70]. As the velocity field near the interface is very sensitive to variations in the curvature, we use a gaussian smoothing filter to damp high frequency perturbations in the speed and interface position, in accordance with [71].

Results and Discussion
We employ a computational Ω of domain dimension [0,20] × [0,20] and grid size 256 × 256. The initial tumour shape is slightly elliptical and cantered at the domain and is represented by the following equation. x y − − + = (17) To minimize the influences of the computational-domain boundary, we stop the simulation when the tumour boundary is within four distance units of the boundary of the computational domain. We assume that the ECM is initially not degraded (i.e., ECM=1 and MDE=0 in the entire domain). To assess tumour growth, we calculate the tumour area and surface perimeter at regular intervals. Moreover, following [24], we define a shape parameter S = (Perimeter) 2 / π*Ar(4 ea) and length scale LS=2*Area/Perimeter to characterize the morphological characteristics of the tumour. The shape parameter S is a measure of the non-circularity the initial tumour shape; when S=1, the initial tumour shape is a circle. The length scale LS is a measure of the smallest dimension of the tumour. For example, for a rectangular tumour with width W and length L, LS = LW/ (L + W) and LS ≈ W if W<<L.

Effect of haptotaxis due to ECM gradients from MT-MMP mediated proteolysis
In this section we examine the influences of haptotaxis due to localized ECM degradation by MT-MMP with variations in nutrient availability in the tumour microenvironment subject to the variations in surface tension, proliferation and apoptosis rates. We also compare these effects with MMP-mediated effects.

Nutrient concentration
Pressure P net p tot The extent of haptotaxis of intact-ECM gradients [43] 0.2, 1, 5  throughout the entire computational domain, on a uniform 2D cartesian mesh, except for the Poisson equation for ρ N , which for case (ii) is solved for the tumour region Ω T only. Because we anticipate frequent and complex morphological tumour-surface changes we use the levelset method, based on a reformulation of a model proposed in [63]. Let φ be an auxiliary level-set function whose zero level set denotes the boundary Ω T of an avascular tumour growing into a surrounding, noncancerous tissue, satisfying φ < 0inside Ω T , φ> 0 outside on the tumour boundary Σ T . We update the position of the interface by solving the PDE,

Effect of MT-MMP in a nutrient-poor microenvironment:
Here, we investigate the influence MT-MMP mediated proteolysis on tumour growth and morphology in a nutrient-poor environment, which has not yet been reported in the literature. Figure 2 shows the evolution of the tumour morphology with time for haptotaxis coefficient χ E =1 for both MMP and MT-MMP mediated proteolysis. Relevant parameter values are listed in (Tables 2 and 3). During early times the growth is similar for both (t=48) and forms a neck by creating bulbs on each side of the y-axis that grow in time. The influence of MMP-mediated proteolysis on the growth rate becomes pronounced, as corrugation of the tumour surface continues and bulges are formed on the tumour surface symmetric about both axes (at t=80), which grow in time while MT-MMP creates no bulge along x-axis. The contours of MDEs and ECM are shown in Figure 3 where we can see the diffusible and non-diffusible MDEs and the resulting ECM degradations. The ECM gradients near the boundary of the tumour are responsible for these differences in tumour morphologies via haptotaxis, as shown in (Figure 4). We notice that gradients are largest on the opposing sides of the bulge (regions A and C in Figure 4a for MMP, indicating further hapotaxis induced growth of the bulge which are smooth for MT-MMP despite having larger gradients than that for MMP (Figure 4c). Similarly the growth of top bulges for both MMP and MT-MMP can be explained by Figure  4b and 4d.
Qualitatively, similar results are found for five-fold increase of the haptotaxis coefficient for MT-MMP whereas the tumour develops four large, symmetric bulbs with very narrow necks, which continue to spread in the computational domain and eventually spilt into pieces for MMP. The effect of an increase in the haptotaxis coefficient on tumour growth and morphology is shown in ( Figure 5). As seen in the plot of tumour area versus time (Figure 5a), the higher the rate of haptotaxis, the greater variations for MMP and MT-MMP. This trend is also reflected in the plot of tumour perimeters with time, presented in (Figure 5b). Growth is faster χ for MMP and more distinguishable as E increases.  It is not known whether decrease in adhesion can cause similar growth and morphologies due to MMP and MT-MMP mediated proteolysis in nutrient-poor microenvironments. Therefore, we decrease the surface tension by 10-fold for the simulation results shown in (Figure 6). We get more fingers for MMP compared to that for MT-MMP (c.f., Figure 6). The difference between the growth rates and morphological details also increase as surface tension decrease as shown in Figure 7 comparing with Figure 5. Therefore, the results of our simulations show that MT-MMP mediated proteolysis has less effect on tumour growth and morphology compared to that of MMP mediated proteolysis in a hypoxic microenvironment. For similar cases, the tumour develops more fingers for MMP than MT-MMP. So MMP mediated proteolysis leads to more tumour invasion.    are given in (Tables 2 and 3). The evolution of the tumour morphology with time for both proteases is displayed in Figure 8 for χ E =1. The tumour develops an H-shaped structure in this case, also observed in Zheng, et al. [41] for no taxis movement although the tumour grows faster due to the combined effect of cell proliferation and haptotaxis. The tumour morphologies are very similar for both MMP and MT-MMP even after an increase in the haptotaxis coefficient by a factor of 5 (not shown here). The only difference being that the buds tend to bend inward more for MT-MMP. Although the buds approach each other, they never reconnect. Rather, they continue to grow toward the boundary of the computational domain. Figure 9 shows the effect of haptotaxis on tumour growth and morphology for both proteases in high nutrient environment. We note that the growth is faster for MT-MMP for χ E =0.2 and χ E =1 unlike for the case of tumour growth in a nutrient-poor environment. However, for χ E =5, both proteases gives comparable rate of growth. This indicates lower haptotactic sensitivity for nutrient-rich tumour growth. Still, qualitatively, the tumour morphologies for both proteases in nutrientrich environment are very similar to that found for a tumour with no taxis even with increase in haptotaxis coefficients. Now, we investigate the effect of decrease in cell adhesion on tumour growth dynamics and morphologies in nutrient-rich environment. If we decrease the surface tension by ten-fold, we get no change in morphologies for MMP and MT-MMP. However, we notice that MMP gives faster growth as the coefficient of haptotaxis χ E =5, whereas for χ E =0.2 and 1, their growth rates are very similar (c.f. Figure10). Growth is faster for MT-MMP when χ E is small, but for larger χ E , comparable growth. So for nutrient rich environment, the surface tension and haptotaxis coefficient mainly determine the growth rates for MMP and MT-MMP. The results of our simulations, therefore, clearly show that the effects of MMP and MT-MMP on tumour growth and morphologies are not significantly different for a tumour growing in a nutrient-rich    (Tables 2 and 4). We investigate the effect of haptotaxis under the same conditions described in section 4.1.1, except with a haptotaxis coefficient, χ E =10 and proliferation rate, λ p = 0.1 and with no apoptosis. For the case of the larger proliferation rate, λ p = 0.5, haptotaxis results in stable λ tumour growth. However, for p = 0.1, Figure 11 shows that the tumour grows into four symmetric buds for MMP in constrast to the MT-MMP case where two symmetric buds are formed. The environment. effects of haptotaxis on tumour growth and morphology are presented in (Figure 12). When χ E =5, the tumour area, perimeter and length scale all increase linearly with time; very similar to behaviour found for the case without taxis. However, the shape parameter is constant, indicating stable circular growth. When taxis coefficient is increased twice the tumour grows linearly at early times but much more rapidly at later times. The LS first increases, then decreases with time, indicating narrowing of the fingers. Similarly, the shape parameters are nearly constant at early times but increase rapidly with decreasing length scale, an indication of tumour surface instability. Thus at low proliferation rates, in the absence of apoptosis, MMP gives faster growth and more fingers than MT-MMP. Similar results we reported in section 4.1.1, where apoptosis and high proliferation rates are considered. So in the absence of necrosis, MMP gives faster growth and more fingers in lownutrient environment irrespective of proliferation rate. Similar tumour morphologies for both MDEs are obtained for higher apoptosis rates. Thus, the difference between the effects of MMP and MT-MMP on tumour morphology is not significant for high apoptosis rates.

The effect of MT-MMP at low proliferation rate and high apoptosis
The results of our simulation clearly show that there is no significant variation in tumour growth dynamics and shapes when we consider the haptotaxis due to diffusible MMP and non-diffusible MT-MMPmediated ECM degradation, except for low proliferation rate in lownutrient environment, whereas there are significant variations due to haptotaxis in nutrient-poor microenvironment depending on cell proliferation and apoptosis rates.

Summary and Conclusion
A better understanding of how the tumour environment affects cancer progression should provide new targets for the isolation and destruction of cancer cells via interference with the complex crosstalk established between cancer cells, host cells, and their surrounding extracellular matrix [5]. So we have proposed a two-dimensional continuum model for the growth and invasion of tumour cells into healthy tissue that focuses on four key components implicated in the invasion process: tumour cells, ECM, MDE (MMPs or MT-MMPs) and nutrient which are initially homogeneous in the tumour microenvironment. We consider both diffusible MMPs and non-diffusible MT-MMP mediated proteolysis and their separate effects on tumour growth via haptotaxis. Our model is the first to examine the effect of haptotaxis for localized ECM degradation by MT-MMP on tumour growth and morphology in high and low-nutrient environments and differentiate the influences of haptotaxis for various cell proliferation and apoptosis rates in hypoxic environment for soluble and non-soluble MMPs.
With the inclusion of haptotaxis, tumour growth and morphology are substantially characterized by whether the tumour is growing in a nutrient-rich or nutrient-poor microenvironment. The invasive fingering morphologies are found in nutrient-poor microenvironments, as a result of haptotaxis, consistent with the results of previous hybrid models whereas in nutrient-rich environments there is no significant influence of haptotaxis on fingered tumour morphology. These results show that ECM degradation by MMP leads to greater instability than that for MT-MMP in low-nutrient environments, even at low proliferation rates; whereas for high apoptosis rates, these both lead to similar unstable morphologies. In summary, our study shows that while haptotaxis lead to completely different tumour growth rates and morphologies, depending on proliferation and apoptosis rates, in lownutrient environment as demonstrated earlier, there are no significant variations when we compare haptotaxis due to ECM degradation by diffusible MMP and surface bound MT-MMP except for low proliferation rate [71][72][73].