Abstract
The goal of this paper is to construct a new algorithm for the numerical simulations of the evolution of tumour invasion and metastasis. By means of mathematical model equations and their numerical solutions we investigate how cancer cells can produce and secrete matrix degradative enzymes, degrade extracellular matrix, and invade due to diffusion and haptotactic migration. For the numerical simulations of the interactions between the tumour cells and the surrounding tissue, we apply numerical approximations, which are spectrally accurate and based on small amounts of grid-points. Our numerical experiments illustrate the metastatic ability of tumour cells.
1. Introduction
The analysis of data obtained from the World Health Organization (WHO) [1] and the UN [2] databases shows that, at present, cancer is and probably will remain to be among the leading causes of death worldwide [3–5] being surpassed only by cardiovascular diseases. According to the data provided by the WHO, cancer disease is the cause of the death of roughly six million people yearly [1]. This explains the major significance of the fight against the malignant conditions, which includes prevention [6], cure [7, 8], and cancer research.
Tumour development is a very complex multistep process involving many intracellular and extracellular phenomena which are strongly nonlinear and time varying [4, 9–11]. Genomic changes as well as microenvironmental factors such as the extracellular matrix (ECM), various growth factors, and substrate concentrations have been shown to play a major role in the process of carcinogenesis [12].
Generally, tumours can be classified as benign and malignant. The growth of benign tumours is self-limiting and their cells tend to stay in the same place. Malignant tumours may grow without limitations and their constituent cells are prone to migrate or metastasize to other parts of the organism [13–15]. The ability of malignant cancer to invade the local tissue and to spread throughout the organism is their most insidious and dangerous property. Metastasis is the predominant cause of most cancer deaths [14, 16, 17].
The process of metastasis includes angiogenesis and invasion. Tumour angiogenesis (rapid growth of blood vessels near the tumour cells) is induced by a secretion of various growth factors such as vascular endothelial growth factor (VEGF). These vessels facilitate the influx of oxygen and other nutrients needed for the development of the cancer [18]. The process of angiogenesis is followed by invasion and penetration of cancer cells into surrounding tissues and possibly by dissemination of cancer cells through blood vessels. Thus, tumour cells can be carried to a distant site of the body. There they can implant and initiate the development of a secondary tumour [14, 16, 19]. An important role in the process of cancer invasion is performed by matrix degradative enzymes (MDEs) such as metallo-proteases (MMPs). They are produced by tumour cells and digest the ECM, which enables the migration of cancer cells through the tissue [13, 14, 17].
In the last half century, many mathematical models describing the process of tumourigenesis have been the subject of active research. Mathematical and computational methods have contributed to clarifying the factors that are sufficient to explain experimental and clinical data, to defining these factors in precise terms and to suggesting experiments for calculation of these factors [20]. In addition, analyses and simulations of mathematical models have been used for the reduction of the amounts of costly experiments needed for the development of therapies [21, 22]. It is strongly believed that mathematical and computational methods will play a significant role in cancer research in the future. They may improve the understanding of some complicated features and details of tumour evolution as well as be effectively used in clinical laboratories, by means of appropriate model-based decision support systems [4]. We refer the readers to special issues [23–27] for more complete bibliography regarding the applications of mathematical and computational methods to cancer research.
Gatenby and Gawlinski present one of the first models of tumour invasion in the papers [28, 29]. Gatenby [28] considers the competition between healthy host cells and modified (tumour) cells and proposes and analyses several models formulated in terms of ordinary differential equations. Gatenby and Gawlinski [29] present a reaction-diffusion model for the investigation of the role of the alteration of the microenvironmental acidity induced by cancer cells for their invasion into the organism. Subsequently, the series of papers, among others [30–43], have appeared offering models and detailed analysis of diverse features of cancer invasion. In this paper, we study the continuum models of avascular tumour growth investigated by Chaplain et al. (cf., e.g., [31, 34–37]). The first model of this series is proposed in Anderson et al. [31]. The authors consider three major variables involved in the process of cancer invasion, namely, cancer cells, ECM, and MDEs. In order to study in detail mainly the influence of the surrounding tissue on the process of migration of tumour cells, the proliferation of the latter is not included in the continuum model. The authors analyse numerically in one and two dimensions the impact of ECM gradients resulting from the destruction of ECM by MDE and the role of haptotaxis on cancer invasion. An extension of this model is presented in Chaplain and Anderson [34] who consider the role of oxygen as a nutrient for the tumour cells. The authors propose also a new model equation for endogenous inhibitors, such as tissue inhibiting MMPs, that can neutralize MDEs. We include this equation in our model (8), see Section 2 below. The model of Chaplain and Anderson [34] has been further developed by Lolas [37] and Chaplain and Lolas [35, 36] who have considered terms describing chemotaxis, proliferation of cancer cells and reestablishment of the ECM. Lolas [37] examines a variety of continuum models, in particular incorporating the effects of just chemotaxis, and haptotaxis, and their combination, and so forth. One of the conclusions of the author is that the mechanism of chemotaxis without haptotaxis cannot lead to a successful cancer invasion if there is no proliferation of tumour cells and reestablishment of ECM. Further novel ordinary differential equations that describe the cancer cell proliferation and the remodeling of the extracellular matrix re-establishment function allowing the incorporation of the plasminogen activation cycle are included in the model of Chaplain and Lolas [36] that also investigates the role of the uPA system for the cancer invasion. uPA inhibitors and plasmin have also been investigated in the model by Chaplain and Lolas [35]. Clear and detailed description of the biological processes observed during the cancer invasion and metastasis is provided in [31, 34–37]. In particular, in these paper, the key stages of the metastatic cascade, the structure and functions of the major constituents of the ECM and the basic representatives of the MDEs participating in the interactions between the healthy and cancer cells are systematically presented on the basis of broad theoretical and experimental bibliography.
In our paper, we propose a different numerical approach than the approach used, for example, in [31, 34–37]. The goal of the paper is to obtain numerical results which are based on small amounts of spatial grid points applied to the model equations so that low-dimensional vectors of data are used to make the numerical computations fast. We construct a new algorithm for the systems [31, 34, 36, 37] by using spectrally accurate approximations to the terms that model the tumour cell random motility, the haptotaxis, the MDE diffusion, and the diffusion of the endogenous inhibitors. Since the algorithm computes the solutions with spectral accuracy, it is based on smaller amounts of spatial grid points than the amounts of grid points used for the less accurate finite difference approximations (strategy applied, e.g., in [31, 34–37]), which consequently saves computational time. The idea of using small amounts of spatial grid point and saving time for computing one solution for one set of parameters, which has to be repeated many times for many sets, is important, for example, for the numerical experiments carrying the goal of estimating parameter values from laboratory data. This idea is applied in [44] to estimate parameter values of one of the models presented in [36, 37] from the in vivo experimental data [45] developed by using transgenic mouse models. The numerical approach from [44] is based on a different approximation to the haptotactic term than the approximations used in this paper and our numerical schemes are constructed for systems which are various variants and generalizations of the model investigated in [44]. Furthermore, because of considering different variants of boundary conditions the schemes in this paper differ from that of the paper [44].
Additionally to the model presented in [36, 37] and applied in [44], in this paper, we investigate other models, which are presented in [34] or are combinations of the model equations from [34, 36, 37]. Moreover, in [44], the parameter values are evaluated quantitatively from the laboratory data [45] so that the solutions of the model equations correlate with the data. Contrarily to [44], in this paper, we choose the parameter values qualitatively in order to observe and compare solutions computed with different parameters. This comparison allows to analyse the influence of the parameters on the shape of the solutions and we conclude that complicated interactions between tumour cells, ECM, MDEs, and endogenous inhibitors can be directed by choosing the parameter values. Our sequence of numerical simulations is initiated from the solutions obtained with the parameter values chosen in [34] (for comparison) and next we gradually change the values and analyse their influence on the solutions. Animated graphical visualization of the solutions and how they change according to the parameters is helpful in observing the influence of the parameters on the shape of the solutions. The idea of using small amounts of spatial grid points and saving time for computing solutions of the model equations is crucial in the effective utilization of, for example, animated simulations of tumours, which can be used as a predictive and visualized tool in clinical applications. Decreasing the amounts of spatial grid points used for such visualizations saves not only the time of demonstrations but also the computer memory. It is not possible to demonstrate the animated simulations in papers and we only note that they are interesting and help in visualization of the complicated biological processes. Instead of the animated simulations we include snapshots at different stages in time.
The contents of this paper is as follows: the model equations are described in Section 2, the algorithm is introduced in Section 3, the results of numerical experiments and simulations are presented in Section 4, and Section 5 includes our concluding remarks and future research work.
2. Mathematical Model
In this section, we investigate various models of tissue invasion by cancerous cells. In Section 2.1, we investigate the Chaplain and Anderson model [34] focusing on interactions between ECM and cancer tumour and metastatic abilities of cancer cells. In Section 2.2, we investigate further expansions of the model and its different versions with additional terms connected with proliferation of tumour cells, ECM renewal, and different functions modelling the production of MDEs by the tumour cells. Section 2.3 deals with a more general model with an additional equation, which describes evolution of endogenous inhibitors.
2.1. Cell-Matrix Interactions and Cell Migration
In the next section, we construct a numerical scheme for the following model of tissue invasion: with the space variable belonging to the scaled domain of tissue, and time . The model equations (1) describe interactions between tumour cells, MDEs, and ECM. The interacting variables are -tumour cell density, -ECM density, and -MDEs concentration. The system (1) is derived in detail in [34] and is a part of a more general system consisting of (1) with an additional fourth equation for endogenous inhibitor concentration denoted by . In [34], it is assumed that the tumour cells, the MDEs, and the inhibitors remain within the space domain and zero-flux boundary conditions are imposed. The fourth equation for the endogenous inhibitor concentration is dropped under the additional assumption that negative effect of the endogenous inhibitors is overcame by the MDEs in an actively invading tumour. This assumption implies that and the general system of four equations is reduced to (1). In Section 2.3, we investigate the model with all four equations.
2.2. Migration and Proliferation of Cancer Cells, ECM Renewal, and MDE Production
We also investigate further expansions of the model (1), which are introduced, for example, in [36, 37]. After adding the proliferation term to the right-hand side of the equation governing tumour cell motion (the first equation in (1)) and the ECM renewal term to the right-hand side of the equation for the ECM (the second equation in (1)), we obtain the following model: where is the proliferation rate of the tumour cells and is the growth rate of the ECM.
We also make experiments with the following modification of (2): where the MDE production is modeled by . The motivation for choosing such form of the MDE production in [36, 37] follows from experimental observations of polarized expression of MDEs at the invading leading edge of tumour, see, for example, Estreicher et al. [46].
We investigate the model equations (1), (2), and (3) supplemented by the zero-flux boundary conditions at and either the Dirichlet conditions or the zero-flux boundary condition at . As in [34], we assume that the initial tumour is centered at , the initial MDE concentration is proportional to the initial tumour cell density with as the constant of the proportionality, and the MDE has already degraded the ECM, thus we consider the same initial conditions as in [34], which are the following: for . The parameter values for the model equations are specified in Section 4.
2.3. Production of Endogenous Inhibitors
We additionally consider the general model where the last equation describes evolution of endogenous inhibitors (concentration of which is denoted by ). This equation is the fourth equation in the model (10.5) proposed by Chaplain and Anderson in [34], where it is assumed that endogenous inhibitors are produced by ECM as a response to the MDEs and the function models the inhibitor production. The term models neutralization of the MDEs and models decay of the inhibitors. We assume that the initial inhibitor concentration is and impose the zero-flux boundary conditions
Our goal is to construct a new efficient algorithm for solving the models (1), (2), (3), and (8) and investigate the ability of cancer cells to produce and secrete the MDE, which then degrade the ECM, and allow the cells to start their migration towards healthy parts of the tissue.
3. Construction of Numerical Approximations to Tumour Cells, ECM, and MDEs
In this section, we construct numerical solutions to the model equations (1), (2), and (3) supplemented by the initial conditions (7) and the boundary conditions (4) and (5). For the numerical solutions, we consider the Chebyshev-Gauss-Lobatto points with , in the scaled domain of tissue. Our goal is to construct approximations to , , and , for ; the values of the solutions at are known from (5).
Let and we use similar notations for and . We shall replace the spatial derivatives in (1) by numerical approximations constructed for the vectors , , , in the first equation and for in the third equation. For , we apply the following spectrally accurate approximations with , where is the first-order differentiation matrix based on the points (11), see [47, 48]. We also apply the analogous spectrally accurate approximations for and , respectively.
Since the exact value of is given by (4), the approximation (13) is not needed at the first point , that is, for the first component of . Therefore, from (13), the first approximation in (15) with and , and from the boundary conditions (4) and (5) we obtain where From (13) and (16) we obtain the following approximation for the second-order derivatives where
We now construct approximations to , , , and . From the spectrally accurate approximations (15) and from (4) and (5) we obtain According to (20), we have the following approximation for the second-order derivative of with the notation From (21) we have with the similar notation for
We now replace the spatial derivatives in the model (1) by their corresponding approximations. We apply (18), (16), (20), and (22) to the first equation in (1) and obtain its discrete version written in the following form where stands for the component-wise multiplication between two vectors. The discrete version of the second equation in (1) is written in the form and from (24) we obtain the following discrete form of the third equation in (1) The resulting system (26)–(28) is composed of ordinary differential equations and is a semidiscrete version of (1). Note that since the spatial derivatives in (1) are approximated with the spectral accuracy, much smaller numbers of grid-points are needed for (26)–(28) than for finite difference schemes, and time integration of the smaller systems is more robust and more efficient than time integration of the finite difference systems.
For the models (2) and (3) supplemented with (4) and the right-hand side boundary condition (6), which is different than (5), we need to apply different approximations than (16), (18), (20), and (22) as they include (5) instead of (6). For this problem, from (13), instead of (16), we obtain where Further, from (29), we obtain As in (13) and (15), for the vector of approximations to the haptotactic term we obtain and instead of (24), from (6), we obtain
From (29), (31), (33), and (34) we obtain the following scheme for the problem (3), (4), and (6) where is a vector entries of which are all 1-s. For (2), the component needs to be replaced by in the last equation of (35). From the boundary conditions (10), we obtain the approximation for the diffusion of the endogenous inhibitors and the semi-discrete version for the last equation in (8) is written in the following form where we assume that the inhibitor production is modelled by . The semi-discrete equations have to be closed by initial conditions chosen according to (7) and (9).
4. Numerical Experiments
We apply the approximations introduced in Section 3 and begin our series of numerical simulations from (26)–(28), which correspond to model (1). Results of our numerical experiments are presented in Figures 1–6. We use the parameter values , , , , , and different values of and specified in the captions of the figures.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
We apply , that is 32 grid-points , for the numerical experiments presented in Figures 1–4. The time of integration of the system (26)–(28) based on 32 grid points is 0.22 sec to compute the numerical solutions presented in Figures 1 and 2 and 0.38 sec to compute the numerical solutions from Figures 3 and 4. For Figures 5 and 6, we apply , that is 45 grid points, and in this case, the time of integration of (26)–(28) is 0.35 sec.
Figures 1, 3, and 5 show snapshots in time and Figures 2, 4, and 6 show continuous evolution in time of tumour cells, ECM, and MDE, and their interactions for all in the space domain. The numerical results presented in Figures 1 and 2 were obtained with and , see also [34, Figure 10.2]. The results from Figures 3 and 4 were obtained with and .
Two distinct clusters of tumour cells are seen in Figures 1 and 3 at and . The numerical results show that the new clusters, which are not seen at and appear at and , are created at the leading edge of the tumour as a result of the diffusion and haptotactic migration modeled by the two components from the right-hand side of the first equation in (1): random motility and haptotaxis , respectively. Since is greater for Figure 3 than for Figure 1, because of larger haptotactic migration in Figure 3 than in Figure 1, the two clusters seen in Figure 3 are more separated from each other than the two clusters in Figure 1. The pictures show the effect of haptotaxis. The small clusters of cells, which break away from the main body of the tumour, illustrate the potential for the cancer cells to degrade the surrounding tissue, migrate, and start the metastatic cascade. The migrations of the small clusters may not be detected during the processes of medical treatments, and even after resections of the main tumours, the new small clusters may initiate recurrences of the disease. A new cluster of tumour cells broken away from the main body of the tumour is also observed in Figures 5 and 6, which present numerical data computed with and .
The next part of our numerical experiments concerns the models (2) and (3) supplemented by (4), (6), and (7). The results for model (2) are presented in Figures 7, 8, 11, and 12 and for model (3) in Figures 9, 10, 13, and 14. These experiments start from the initial condition (7) corresponding to the snapshot in time in Figure 1. We observe that the small clusters of cancer cells separated from the main tumours are better formed at than at and as time evolves the haptotactic migration together with the production of new cancer cells spread the shapes of the tumours over the -domain. We also observe that the snapshots in time in Figures 1, 7, and 8 look similar to each other and the models (1), (2), and (3) give similar results for although they are supplemented by the different boundary conditions, either (5) or (6), and solved with different parameters and . However, these similarities are observed only for and as time evolves the corresponding solutions of the models (1), (2), and (3) differ from each other. For example, already at , Figure 8 shows greater production of tumour cells than Figure 7. Moreover, at and , Figure 7 shows weaker MDE production and greater production of the tumour cells and the ECM than in Figure 1 due to the fact that , , and are greater in Figure 7 than in Figure 1. On the other hand, also at and , Figure 8 shows greater MDE production than Figure 1. Although for Figure 8 and for Figure 1, since the MDE production is greater in Figure 8 than in Figure 1, the MDE concentration is greater in Figure 8 than in Figure 1 and consequently the ECM degradation is more progressive in Figure 8 than in Figure 1. It can also be observed that although the parameters , , and in the third equation of (1) and (2) are the same for Figures 1, 7, and 8, the MDE curves show different MDE concentrations in all of these figures.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
In Figures 7–10, we observe differences due to the MDE production terms and in (2) and (3), respectively. The parameter values for Figure 9 are the same as for Figure 7 and the parameter values for Figure 10 are the same as for Figure 8 but the MDE production is weaker in Figure 9 than in Figure 7 and also weaker in Figure 10 than in Figure 8. The MDE concentrations in Figures 9 and 10 are more uniformly spread out across the -domain than in Figures 7 and 8. Furthermore, Figures 9 and 10 illustrate that the system (3) models a decreasing MDE production as the tumour invasion progresses and consequently less MDE concentrations in the regions where the high tumour cells densities are situated (possibly in these regions the tumour cells already finalized their invasion due to lack of space and move to other regions) than in the regions of less advanced stage of invasion. This shows that cancer cells may not need MDEs in the regions where they deal with lack of space. Similar features are observed in Figures 11–14, where the clusters separating from the main tumours are more profound than in Figures 7–10. The same parameter values were used for Figure 11 as for Figure 13 and the same parameter values were used for Figure 12 as for Figure 14. Figures 11 and 12 illustrate numerical solutions to (2), while Figures 13 and 14 illustrate numerical solutions to (3). Comparison of Figures 11–14 confirms that the term in (3) models lower MDE production than the term in (2) (with the same parameter values). We also observe that since MDE production is lower in Figures 13 and 14 than in Figures 11 and 12, the ECM degradation is smaller in Figures 13 and 14 than in Figures 11 and 12.
Figure 15 shows snapshots in time of the four solutions to the more general model (8) describing the interactions between the tumour cells, ECM, MDEs, and endogenous inhibitors. All four profiles show that, as time evolves, the ECM produces endogenous inhibitors, concentration of which increases in time and their higher concentration is located in the regions where ECM is not yet entirely degraded rather than in the regions where the degradation is already effectively developed. The inhibitor profile shows that the ECM responds to the MDEs by producing the endogenous inhibitors.
(a)
(b)
(c)
(d)
5. Concluding Remarks and Future Directions
We have constructed a new numerical algorithm for fast computations of the solutions of the mathematical models proposed by Chaplain et al. in [34, 36, 37], which consist of systems of nonlinear partial differential equations describing interactions between tumour cells, ECM, and MDEs. The algorithm is based on spectrally accurate approximations and small amounts of grid points, which results in ordinary differential systems of small dimensions and fast computations. We have applied the algorithm and presented and compared the numerical simulations with a variety of model equations. The simulations demonstrate that the models describe important features of the interactions between tumour cells and the surrounding tissue, and in particular the initiation of a new colony of cells and metastasis.
Our future research work will address the question for which parameter values and domains the model [34] and the kinetic type model proposed in [49] are equivalent. We will also address numerical methods with spectrally accurate approximations for the models with two-dimensional spatial domain and with different kinds of the function modelling the inhibitor production [34].
Acknowledgment
The authors wish to express their gratitude to anonymous referees for the useful remarks and comments which led to the improvements in the presentation of the paper.