- Email: [email protected]
Numerical interactions of random and directed motility during cancer invasion
MATHEMATICAL COMPUTER MODELLING PERGAMON
Mathematical
and Computer Modelling
30 (1999) 123-133 www.elsevier.nl/locate/mcm
Numerical Interactions Directed Motility During
of Random and Cancer Invasion
A. J. PERUMPANANI Department of Surgery, MassachusettsGeneral Hospital Harvard Medical School, Building 1400W., One Kendall Square Cambridge, MA 02139,U.S.A. [email protected]
J. NORBURY Mathematical Institute, University of Oxford Oxford OX3 9DU, U.K. (Received and accepted January
1999)
Abstract-The continuum modelling of cell migration during cancer invasion results in the coupling of parabolic and hyperbolic partial differential equations (PDEs) arising from the random motility of normal tissue and the directed movement up substrate gradients of malignant cells. The numerical solution of such systems of equations require different stability criteria being simultanc+ ously satisfied. We show that in such a coupled system, the origins of numerical instability can be identified by analyzing the fastest growing mode in a numerically unstable solution. In general, stability can be achieved by choosing an appropriate grid size representing the more stringent of the conditions for hyperbolic and parabolic stability. However, this induces variable degrees of numerical diffusion because of a changing CFL (Courant, Friedrichs, and Lewy) number. Solving the hyperbolic and parabolic PDEk on separate grids results in a better convergence of the solution. Finally, we discuss the use of higher-order schemes for the solution of such problems. Cancer modelling brings together directed and random motility in a unique way thereby presenting interesting new numerical problems. @ 1999 Elsevier Science Ltd. All rights reserved.
Keywords-Mixed
hyperbolic parabolic, Cancer invasion, Operator splitting.
1. INTRODUCTION The morbidity and mortality associated with cancer is consequent to the ability of cancer cells to invade neighbouring tissue, blood vessels, and lymphatics, and thereby spread to various organs away from the primary site. This sequence of steps, called metastasis, begins with the
directed movement, of invasive cancer cells into the surrounding extracellular matrix and normal cells [l]. From a simplistic viewpoint, cancer invasion may be viewed as the overwhelming of the minimal random motility of normal cells by the enhanced directed motility of mutated cancerous cells. Directed cell movement along substrate gradients is an early feature of the developing embryo which is subsequently manifested only in specialized situations such as wound healing, We thank J. Sherratt and I. Dawkins for several useful discussions on this paper and A. Tilles for a careful reading of the manuscript. 08957177/99/s - see front matter. @ 1999 Elsevier Science Ltd. All rights reserved. PII: SO8957177(99)00169-7
‘b-et
by -4M’W
124
A. J. PERUMPANANI AND J. NORBURY
inflammation, placentation, and cancer metastasis [2]. In contrast, random motility is associated with the normal physiological processes of growth and development. Examples of simple random diffusion in the body include the nutrition of tissues, distribution of growth factors, and the respiratory exchange of gases. There has been much recent interest in the modelling of cancer metastasis both to understand the basic processes involved and to generate hypothesis that can be tested in the laboratory. Marusic et al. [3] and Perumpanani et al. [4] have developed models for tumour growth while McElwain et al. [5] have shown that cancer cells undergo passive motion towards the center of a multicell spheroid. Sherratt [6], Gatenby et al. [7], and Perumpanani et al. [8] have modelled various aspects of cancer invasion while Sherratt [9] and Albert et al. [lo] have studied the immune response to cancer. A closely related phenomenon which shares many phenotypic features with invasion is the process of tumour angiogenesis by which new blood vessels grow into a tumour; this has been modelled by Byrne et al. [l l] and Chaplain et al. [12]. This paper discusses a model for the invasion of normal tissue by a growing mass of cancerous tissue with particular focus on the numerical behaviour of the model equations. An interesting feature of this model is the interaction of the directed movement of the invading cancer cells (which results in a hyperbolic PDE) and the random motility of the proteases (which results in a parabolic PDE). A detailed derivation of the model and its biological implications are discussed elsewhere [8]; in this paper we focus on the numerical behaviour of this system of PDEs.
2. THE
MATHEMATICAL
MODEL
Briefly, the interaction of invasive cancer cells (~(z,t)) with the surrounding extracellular matrix (ECM, c(z,t)) through the liberation of a protease (p(z,t)) is modelled using a system of conservation equations. Here x and t represent one-dimensional space and time, respectively. The local proliferation of the cancer cells is described by the logistic growth model [13] subject to a local maximum (the carrying capacity) which has been scaled to unity. Invasive cancer cells have been shown to move up gradients of connective tissue [l], which is modelled as an invasive flux proportional to the ECM gradient .g in equation (1) below. The interaction of invasive cells with the ECM has previously been shown to induce production of proteases [14] and is modelled by a protease production term proportional to UC. Besides being naturally degraded (-krp), the proteases are also neutralized by the production of anti-proteases (Tissue Inhibitors of Metalloproteases (TIMPs)). The TIMPs are produced by the cellular elements of the connective tissue and at a lower rate by the invasive cells. Hence, the neutralisation of MMPs by TIMPs are modelled by two separate terms -kspc and -kGpu. Since the cells are the source of the TIMPs their concentration is taken as proportional to the concentration of the cells. Being small biochemicals the proteases are capable of diffusing into the surrounding tissue which we modelled as simple Pickian diffusion in equation (3). The extracellular matrix molecules which are fixed to the surrounding tissue are degraded by the action of the proteases which we modelled by the degradation term proportional to pc in equation (2). The complete model is now written as proliferation dU -=aat
directed motility ,-a-,
k2g
[ 1 UE
,
ac - = -kapc} ECM degradation, at 8P at=
* protease production
-~pc-kspll-k~~+~. protease degradation
8% diffusion
CONDITIONS. In the context of malignant invasion, we are studying a domain which is much larger than the size of the tumour, since the tumour is local&d to one part of the tissue.
BOUNDARY
Numerical Interections
125
In these simulations, we used a very large spatial domain so that the effects of the boundary were minimal. Zero Aux boundary conditions were used because the interest here is in the behaviour of an isolated system free from any external inputs. In our simulations we solved system (J)-(3) using step function initial conditions after nondimensionalisation. The total tissue density (u+c) through the whole domain (0 < z < L) was uniform. The invasive cells (u) were initially confined near z = 0, (u(z,O) = 1 for 0 c z < L/10, u(z, 0) = 0 elsewhere) and the remainder of the domain had ECM (c(z, 0) = 1 - u(z, 0)). There was no protease to start with (p(z, 0) = 0). INITIAL
CONDITIONS.
3. NUMERICAL
SOLUTIONS
We investigated the model in Section 2 to study the rate of invasion of the ECM (c) by the invasive cells (u) as a function of the system parameters (the subscripted ks). A typical simulation is shown in Figure 1, in which an advancing wave of invasive cells replaces a receding wave of ECM which is digested by an intervening pulse wave of protease. These simulations can be used to determine the parametric changes which will produce maximal inhibition of invasion, which can then be tested in laboratory assays as potential therapeutic modalities. For example, a decrease in k2 would correspond to therapeutic blockade of the haptotactic receptor which causes directed migration up fixed gradients. Because of the nonlinear nature of these equations the speed of the invading wave of cancer cells could depend on the parameters in a nonlinear fashion. In such a case, therapeutic blockade of the receptor could produce an unintended but potentially dangerous augmentation of invasion which can be predicted by modelling. The numerical solution of (l)-(3) is complicated by the interaction of the diffusional behaviour of the protease and the directed motility of the invasive cells. While equation (3) is parabolic because of the ksp,, term, equation (1) is a nonlinear convection equation where the convection coefficient is kgc,,. This combination of different types of PDE behaviour imposes different stability criteria on the numerical solutions. Furthermore, to get an accurate measure of the speed of invasion numerical diffusion should be kept to a minimum in the solution of equation (1).
4. STABILITY In order to numerically solve system (l)-(3), with upwinding of the form “+I3
-‘y
CRITERIA we first used an explicit finite difference scheme
= klUJ? (1 _ Uj” ) - k,&CS-JJJ’
At q+l - C? ’ = -k3PjnC;, At
- k$Yj’S;Cjn,
where UT, CT, and PT are the values of u, c, and p at the jth spatial point and the nth time step. Here 6, is the central difference operator (( UT+l - Ujn_i)/2Az) and S-, is the backward difference operator ((UT -UT-,)/A z ) w h ic h are used in accordance with the upwind scheme [15]; Ax and At are the spatial and temporal step sizes, respectively. Since (l)-(3) is a mixed parabolichyperbolic system of equations, the criteria for stability impose different constraints on the mesh and its refinement path. To illustrate this, consider the two equations representing the dispersion terms alone in (l)-(3)
where a is equivalent
ut + au= = 0,
(7)
Pt - Pzz = 0,
(8)
to cz in the full model and kz = ks = 1. For explicit
finite difference
126
A. J. PERUMPANANI
J. NORBURY
AND
0.9
0.8 0.7 0.6
0.9 0.8
0.6
p:: 0.3 0.2 0.1 n
Figure 1. Finite difference solutions of the system of PDEs (l)-(3) using the Hyperbolic-Parabolic Separation (HPS) method described in Section 4. Here the u equation (1) and the p equation (3) are solved on separate grids with similar spatial grid points. The time step (At) at each iteration is chosen so as to maximise the CFL number for the u equation and to satisfy the stability condition for the p equation. The solutions are temporally matched periodically with intervening values for u de termined by linear interpolation. The parameter values used in this solution are kl = 1.72, kg = 4.0, ka = 0.345, k4 = 5.0, kg = kg = 0.346, k7 = 0, ]cs = 0.5, and L = 30.
solutions the condition for stability of (7) is aAt/Ax I 1 (commonly referred to as the Courant, Friedrichs, and Lewy condition or the CFL condition) and for (8) is At/(Ao)2 5 0.5 [15]. When the two equations are solved simultaneously on similar grids, the (AZ, At) plane can be divided into four regions based on the stability criteria as shown in Figure 2. They are Region Region Region Region
1: 2: 3: 4:
parabolically unstable, hyperbolically stable; parabolically stable, hyperbolically unstable; both unstable, and both stable.
Numerical solutions using the scheme (4)-(6) s h owed that this is also true for the model system (l)-(3). Though numerical instability started in one of the equations, the system eventually failed as a whole because of the coupling between the equations. In order to understand the origins of the instability, and in particular, to distinguish instabilities arising in Regions 1 and 2 in Figure 2 we studied the solutions of (7)-(g) using Fourier analysis and compared the results with the numerical solution of (l)-(3). Standard Fourier analysis shows that p(s, t) = exp-ICZt sin kx satisfies differential equation (8) for all values of the number k. It can be shown that a similar Fourier mode is an exact solution of the difference equation arising from the finite difference approximation of (8). Substituting
Numerical Interactions
127 2At =Axj
Region 3
2-
Both unstable 1.5.
At
I-
Both Stable
0.5'
0.2
0.4
0.6
0.8
1
1.2
1.4
AX Figure 2. Shows the four regions in the (AZ, At) plane where the set of hyperbolic and parabolic equations described in (7) may be stable or unstable either individuaby or together. In Section 4, we discuss how hyperbolic and parabolic instability arising in a system of PDEk can be distinguished. n&jAz)
PT = (8
I where i = &i’,
into the finite difference approximation
p;+l At and dividing
of (8), i.e.,
PT _ Pjn+1 - 2P7 +#-I
-
(Az)~
(9)
’
by py shows that the Fourier mode is a solution for all values of n and j provided x z X(k) = 1+ V (@AZ - 2 + .c-ikAZ) = 1 - 4vsin2 :kAx,
where v = At/(Az)2 numerical
and X(k) is called the amplification
approximation
can be written
factor for the mode. For k = rnn the
in the form
(11) The low frequency terms in this expansion give a good approximation to the exact solution of differential equation (8), because the series approximation for X(k) and exp(-k2At) match
reasonably well. Under the restriction v 5 l/2 the above scheme is convergent and the error is bounded [15]. For large values of k the modes in the exact solution are rapidly damped by the exponential terms exp( -k2t). But in the numerical solution the damping factor (X(k)1 will become greater than unity for large values of k, if v > l/2; in particular this will happen when kAx = K, for then X(k) = 1-4~. These Fourier modes will then grow unboundedly as n increases. The above analysis is valid for a solution that admits all modes. However, on a discrete mesh there are only a finite number of distinct modes; modes with wave numbers kl and kz are indistinguishable if (kl - kz)A z is a multiple of 27~. Therefore, 21: can be expanded as a linear combination of the distinct modes corresponding to k = mn, where m = -(J - 1)) -(J 2) )...) -l,O,l,..., J and J = l/Ax. The highest mode that can be carried by the mesh has k = Jr or kAx = n; this mode has values fl at alternate points on the mesh. Equation (10) shows that this is also the most unstable mode for the difference scheme and has the amplification factor X(J?r) = 1 - 4~ and S[X] = 0. It is the fastest growing mode when u > l/2 and hence, eventually dominates the solution. However, in solving the hyperbolic equation (7) using the upwind
technique,
the Fourier analysis gives for the constant X E X(k) = l-
(I_
,+kAr)
a > 0 the amplification
= 1 - 5 (1 - e-ikA=)
,
factor
128
A. J. PERUMPANANI
where t = uAt/Ax
is the CFL number. lx(k)]2
AND J. NORBURY
This leads to
= ((1 -Q
+
+ [
(13)
= 1 - 4<(1 - I) sin2 :kAx.
It follows that IX(k) ] 5 1 for all k provided that 0 < [ 5 1. However, when 6 > 1, the highest frequency mode, kAx = X, or uj 0; (-l)j grows, but in this case the fact that S[X] =
-140
5I
I
10I
15
20
25.
301
35
40I
-14b
1
2
X
X
(4
(b)
3
4
I 5
a
9
I 10
2 2
0 -2 -4 -6 II
-a -10 -12 -140
I
5
I
10
15
20
25
30
35
40
6
7
X
X
(4
(4
Figure 3. Numerical solutions of the system of equations (l)-(3) using a finite difference explicit scheme showing two instances of failure. (a) illustrates failure because of the violation of the stability condition for diffusion, that is, ksAt/(A~)~ 5 0.5, where ks is the diffusion coefficient of the protease p. (b) shows a magnified view of the growing oscillations that causes the numerical failure which corresponds to the highest mode that can be carried by the mesh. (c) shows an instance of failure because the CFL number is too high (1.1) which appears similar to (a) but the magnified view in (d) shows that this has a lower frequency than the failure mode in (b).
Numerical
Interactions
5. TRUNCATION
129
ERROR
A suitable way to obtain stable solutions would be to restrict the mesh to the region in the (AZ, At) plane shown in Figure 2, where both the equations are stable. To do this we chose a fixed At and restricted the choice of Ax such that both solutions were stable. In Table 1, we show a numerical estimate of the wave velocity using such a strategy. Even for relatively fine temporal meshes. the velocity did not converge rapidly enough before the conditions for stability were violated resulting in numerical instability. To understand this more clearly we studied the truncation error of the numerical solution of equation (7). Table 1. A comparison of the computed wave speeds for the numerical solutions of (l)-(3) by different numerical schemes. When the temporal grid size At is held constant and the spatial grid size AZ is refined, failure occurs before the wave velocity has converged adequately. By holding the CFL number constant, better convergence is obtained, though instabilities again restrict the refinement of the spatial mesh. However, by solving the three equations in (l)-(3) using the high resolution scheme described in Section 6 on separate grids with different time steps and periodically matching up the solutions better convergence of wave velocities was obtained as shown under HPS (Hyperbolic-Parabolic Separation). For the intermediate values of u solved on a larger temporal grid than the p equation linear interpolation was used. Using the first-order upwind method required Ax = 0.007 to obtain comparable convergence. The parameter values used in this simulation are Icl = 1.72, Ic2 = 4.0, k3 = 0.345, kq = 5.0, kg = ks = 0.346, k7 = 0, and ]cs = 0.1. For the HPS solution we set the CFL number to 0.6.
Ax
At = 0.005
I 0.2 I
1.2201
1
CFL=0.3
HPS
1.1426
1 1.015
0.1
1.1784
1.1390
1.015
0.06
1.1582
1.1361
1.021
0.05
1.1426
1.1327
1.039
0.04
1.1327
1.1275
1.043
tn) t)
The truncation error is defined as the difference between the two sides of the difference equation for (7) when the approximation UT is replaced throughout by the exact solution u(zj, of the differential equation. At any point away from the boundary the truncation error T(x, is defined as T(x,t)
=
6+tu(x, At
t>
-
&r’11(x, Ax
t)
(14)
*
n
The Taylor series expansion may be used to express the truncation error as an infinite series so that u, - ;Axu,, +a .a j (15) * = f (Atu,, - aAzu,,) + . . . . Now even in the case where a is constant, so that utt = a2uz2, T(x, t) = -f (1 - <)aAzu,,
+ ... .
(16)
A. J.
130
PERUMPANANI
AND
J. NORBURY
Hence, the method is only first-order accurate and the truncation error depends on the CFL number c. In the numerical solution, when we changed the mesh size without keeping < constant, the CFL number changed at each refinement with a consequent change in truncation error explaining the poor convergence of the solution. Because of the above observation, we used a second approach which iixed the CFL number fairly low and used this fixed CFL number to determine the time step as At = CFL x Ax/a. The result of this numerical scheme is shown in Table 1. We found that for any particular choice of CFL number the spatial mesh can be refined up to a particular point at which the stability condition for the diffusion equation is violated causing the scheme to fail. The larger the CFL number chosen, the coarser was the spatial mesh size before which failure resulted. In the case demonstrated in Table 1, diffusional instability occurred even before the wave speed had satisfactorily converged. Hence, a scheme with a fixed CFL number was of limited value in solving the model in Section 2. Furthermore, it is important to choose CFL numbers as close to 1 as possible because low CFL numbers introduce large degrees of numerical diffusion as shown in equation (16). Hence, we needed to find a numerical scheme which would keep the CFL number fixed and close to 1, but would not violate the stability condition for the diffusion equation. A satisfactory solution to this problem was found by solving the hyperbolic and the parabolic equations on separate grids and taking the time steps so as to maximise the CFL number for the u equation and satisfy the diffusion stability condition for the p equation. Similar methods have been used in the past for as many time steps as necessary. The advantage of such an approach is that either of the time steps can be broken down into smaller steps should it be necessary to satisfy the stability criteria for that stage without prejudice to the other stage. Wu [lS] described such a method to study constant and variable flow velocities. Valocchi and Malmstead [17] have described the application of operator-splitting to a problem involving advection, dispersion, and reaction terms where they treated the reaction terms in a separate stage. We used a similar approach to solve the model in Section 2. The finite difference equation of (1) was iterated so as to satisfy the CFL condition while keeping the CFL number as closeto 1 as possible. Equation (3) was then iterated choosing a time step so that the diffusion condition was satisfied. For the range of parameters used, this typically meant that several time steps were needed for the p equation to match each time step of the u equation. The intermediate values were determined by linear interpolation. The fourth column of Table 1 shows the convergence of the wave velocity in the numerical solutions using the above strategy under the heading HPS (Hyperbolic-Parabolic Separation). The convergence of this scheme compared favourably with other schemes using a fixed time step or a fixed CFL number. A numerical solution using this scheme is shown in Figure 1.
6. HIGH
RESOLUTION
SCHEMES
First-order accurate schemes for hyperbolic equations suffer from numerical diffusion, but higher-order schemes, such as the Lax-Wendroff scheme [18] while giving higher resolution to discontinuities of the solution, exhibit spurious oscillations around such points. Much effort has been directed at obtaining second-order schemes which give high resolution while keeping spurious oscillations to the minimum. For example, Van Leer (191 and Roe [20] have proposed such high resolution schemes which incorporate terms that dampen spurious oscillations. Such damping terms are called “flux limiters”. This flux is the difference between the flux of a higher-order scheme and that of a low order scheme which has been “limited” in such a way as to ensure the resulting scheme is oscillation free (a term frequently used in the literature is “Total Variation Diminishing”). A high resolution scheme such as the Lax-Wendroff scheme [18] or the Warming and Beam scheme [21] may be used to solve (1). A flux limiter is introduced to dampen the oscillations in the wake of the travelling wave. We solved equation (1) using the Lax-Wendroff scheme with a
Numerical Interactions flux limiter
131
of the form
where vi = k&C~At/Ax and & = kzAt/Ax2 and the $js are the flux limiters. Setting 4 z 1 in (17) yields the Lax-Wendroff scheme and setting 4 = 0 gives the first-order upwind scheme. In order to calculate the flux limiters a quantity rj is first defined such that rj = A-/A+, the ratio of the backward and forward difference operators at each mesh point. Various flux limiters have been described by Van Leer [19], Roe [20], and several other authors. We used the minmod scheme [22] which gave the fastest convergence and defines 4 as
4j =
0,
rj 5 0,
Tj,
0
1,
Tj
5 1,
(18)
> 1.
Use of the minmod scheme in conjunction with HPS described in Section 5 resulted in a faster convergence of the computed wavespeed, though the end results were comparable to the first-order upwind scheme.
7. DISCUSSION The modelling of a large number of biological problems involve the features of directed and random motility which naturally bring together hyperbolic and parabolic equations. For instance, wound healing involves the directed migration of fibroblasts and the random motility of growth factors [23]. During angiogenesis, endothelial cells show chemotactic migration (directed migration along soluble gradients) up gradients of Tumour Angiogenesis Factors (TAFs), while TAFs themselves show random motility [24]. A similar combination is also seen during the emigration of leucocytes [25] and in the formation of a haemochorial placenta [26]. The numerical features mentioned in this paper will have a bearing on the modelling of such problems which study the migration of cells. The model aims at determining the rate of invasion of the cancer cells as a function of the parameters of the model. The rate of invasion was defined as the speed of the cancer cell wave which was measured numerically as the rate of translation of a fixed concentration on the u wave profile. However, the numerical solutions showed a sensitive dependence on the initial conditions. The system supports travelling waves of varying speeds depending on the rate of decay of the initial conditions. A slower rate of decay gave a faster wave. The slowest waves obtained corresponded to initial conditions with compact support which is in agreement with similar observations made with the Fisher equation [27] and the wave solutions of other nonlinear PDEs [5]. Since invading tumours correspond to initial conditions with a compact support, we . m. um wavespeed omin as the index of tumour invasiveness. measured the mmi A biologically interesting prediction made by the model is the dependence of the minimum wave speed omin on the diffusion coefficient of the protease ks. For extremely low values of ks, the minimum wave speed computed is in good agreement with the minimum wave speed obtained from the numerical solution of the purely hyperbolic system obtained with ks = 0, which has been described elsewhere [28]. As ks is increased the computed wave speed initially increased and then decreased. This feature is illustrated in Figure 4; we investigated it for a wide range of parameter values and found it to be a consistent feature of the model. Attempts are now underway in our lab to verify this prediction using an experimental array of collagen gel invasion assays. If confirmed, this prediction could explain the propensity of invasive cells to metastasise to particular tissues in the body where the diffusivity of the proteases they produce could make them maximally
A. J. PERU~PANAN~ AND J. NORBURY
132 1.2
t.15
I.1 cI m&l 1.05
I
0.95
0.9
0
0.2
0.3
0.6
0.8 k8
I
I .2
I.4
1.6
i.8
2
Figure 4. Computed minimum wave velocities a min for a finite difference solution of (l)-(3) for various values of the protease diffusion coefficient .Ics. As 4 is increased, amin initially increases and then decreases. The therapeutic implications of this finding are discussed in Section 7. The parameter values used in this simulation are Icl = 1.72, Icz = 4.0, Ic3 = 0.345, k4 = 5.0, ks = ]cs = 0.346, and k7 = 0.
invasive. The modelling of cell migration not only presents new numerical challenges because of the interaction of different types of motility, but also helps in the formulation of biological hypothesis which can be used to design novel experiments.
REFERENCES 1. S. Aznavoorian, ML. Stracke, H. Krutzsch, E. Schiffman and L.A. Liotta, Signal transduction for chemotaxis and haptotaxis by matrix molecules in tumour cells, J. of CelE Biology 110, 1427-1438, (1990). 2. W.G. Stetler-Stevenson, S. Aznavoorian and L.A. Liotta, Tumor cell interactions with the extracellular matrix during invasion and metastasis, Annu. Rev. Cell. B’iol. 9, 541-73, (1993). 3. M. Mar&c, Z. Bajzer, S.V. Pavlovic and J.P. Freyer, Tumour growth in tivo and as multicellular spheroids compared by mathematical models, Bull. Math. Biol. 56 (4), 617-631, (1994). 4. A.J. Perumpanani, J.A. Sherratt and J. Norbury, Mathematics modelliig of capsule formation and multinodularity in benign tumour growth, Nonlinearity 10, 1599-1614, (1997). 5. D.L.S. McElwain and G.J. Pettet, Cell migration in multicell spheroids-swimming against the tide, Bull. Math. Biol. 55, 655-674, (1993). 6. J.A. Sherratt, Chemottiis and chemokinesis in eukaryotic cells: The keller-segal equations as an approximation to a detailed model, Bull. Math. Biol. 56, 129-146, (1994). 7. R.A. Gatenby and E.T. Gawlinski, A reaction-diffusion model of cancer invasion, Cancer Research 5%,5?455753, (1996). 8. A.J. Perumpanani, J.A. Sherratt, J. Norbury and H.M. Byrne, Biological inferences from a mathematical model for malignant invasion, Invasion Metastasis 28, 209-221, (1996.). 9. J.A. Sherratt and M.A. Nowak, Oncogenes, anti-oncogenes and the immune response to cancer: A mathematical model, PFOC. R. Sot. Land 3 248, 261-271, (1992). 10. A. Albert, M. Freedman and A.S. Perelson, amours and the immune system: The effects of a tumour growth modulator, Math. Biosci. 50, 25-58, (1980). 11. H.M. Byrne and M.A.J. Chaplain, Mathe~tic~ models for tumour progeny-num~cal simul&tio~ and non-linear wave solutions, Bz~ll. Math. Biol. 57, 461-486, (1995.). 12. M.A.J. Chaplain and A.M. Stuart, A model mechanism for the chemotactic response of endothelial cells to tumour angiogenesis factor, Iii&A J. Math. Appt. Med. Biol. 10 (3), 149-168, (1993.). 13. V.G. Vaidya and F.J. Alexandro, Jr., Evaluation of some mathematical models for tumour growth, Int. J. Biomed. Comput. 13, 19, (1982). 14. B. Xie, C.D. Bucana and I.J. Fidler, Density-dependent induction of 92-kD type type-IV coilagenase activity in cultures of A431 human epidermoid carcinoma cells, Am. J. Pathol. 144(5), 1958-67, (1994). 15. K.W. Morton and D.F. Mayers, Nplmericnl Solution of Purtial Diflerential Equations, Cambridge University Press, (1994).
Numerical
Interactions
133
16. J.K. Wu, Wave equation model for solving advection-diffusion equation, Znt. J. of Numerical Methods in Engineering 37, 2717-2733, (1994). 17. A.J. Valochi and M. Malmstead, Accuracy of operator splitting for s&&ion-dispersion-reaction problems, Water Resources Research 5, 1471-1476, (1992). 18. P.D. Lax and B. Wendroff, Systems of conservation laws, Pure Appl. Math. 13, 217-237, (1960). 19. B. Van Leer, Towards the ultimate conservative difference scheme ii. Monotonicity and conservation combined in a second-order scheme, J. Comp. Phys. 14, 361-370, (1974). 20. P.L. Roe, Numerical algorithms for the linear wave equation, Technical Report 81047, Royal Aircraft Estab lishment, (1981). 21. R.F. Warming and R.M. Beam, Upwind second-order difference schemes and applications in aerodynamics, AZAA J. 14, 1241-1249, (1976). 22. P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal 21, 995-1011, (1984). 23. J.A. Sherratt and J.D. Murray, Model of epidermal wound healing, Proc. Roy. Sot. Lond. B 241, 29-36, (1990). 24. M.A. Chaplain and A.R. Anderson, Mathematical modelling simulation and prediction of tumour induced angiogenesis, Invasion Metastasis 16 (4-5)) 222-234, (1996). 25. T.A. Springer, Adhesion receptors of the immune system, Nature 346, 425434, (1990). 26. C.H. Graham and P.K. Lala, Mechanisms of placental invasion of the uterus and their control, Biochem. Cell Biol. 70, 867-874, (1992). 27. J.D. Murray, Mathematical Biology, Springer Verlag, (1990). 28. A.J. Perumpanani, J. Norbury, J.A. Sherratt and H.M. Byrne, A two parameter family of travelling waves with a singular barrier arising from the modelling of matrix mediated malignant invasion, Physica D, (submitted).
Mathematical
and Computer Modelling
30 (1999) 123-133 www.elsevier.nl/locate/mcm
Numerical Interactions Directed Motility During
of Random and Cancer Invasion
A. J. PERUMPANANI Department of Surgery, MassachusettsGeneral Hospital Harvard Medical School, Building 1400W., One Kendall Square Cambridge, MA 02139,U.S.A. [email protected]
J. NORBURY Mathematical Institute, University of Oxford Oxford OX3 9DU, U.K. (Received and accepted January
1999)
Abstract-The continuum modelling of cell migration during cancer invasion results in the coupling of parabolic and hyperbolic partial differential equations (PDEs) arising from the random motility of normal tissue and the directed movement up substrate gradients of malignant cells. The numerical solution of such systems of equations require different stability criteria being simultanc+ ously satisfied. We show that in such a coupled system, the origins of numerical instability can be identified by analyzing the fastest growing mode in a numerically unstable solution. In general, stability can be achieved by choosing an appropriate grid size representing the more stringent of the conditions for hyperbolic and parabolic stability. However, this induces variable degrees of numerical diffusion because of a changing CFL (Courant, Friedrichs, and Lewy) number. Solving the hyperbolic and parabolic PDEk on separate grids results in a better convergence of the solution. Finally, we discuss the use of higher-order schemes for the solution of such problems. Cancer modelling brings together directed and random motility in a unique way thereby presenting interesting new numerical problems. @ 1999 Elsevier Science Ltd. All rights reserved.
Keywords-Mixed
hyperbolic parabolic, Cancer invasion, Operator splitting.
1. INTRODUCTION The morbidity and mortality associated with cancer is consequent to the ability of cancer cells to invade neighbouring tissue, blood vessels, and lymphatics, and thereby spread to various organs away from the primary site. This sequence of steps, called metastasis, begins with the
directed movement, of invasive cancer cells into the surrounding extracellular matrix and normal cells [l]. From a simplistic viewpoint, cancer invasion may be viewed as the overwhelming of the minimal random motility of normal cells by the enhanced directed motility of mutated cancerous cells. Directed cell movement along substrate gradients is an early feature of the developing embryo which is subsequently manifested only in specialized situations such as wound healing, We thank J. Sherratt and I. Dawkins for several useful discussions on this paper and A. Tilles for a careful reading of the manuscript. 08957177/99/s - see front matter. @ 1999 Elsevier Science Ltd. All rights reserved. PII: SO8957177(99)00169-7
‘b-et
by -4M’W
124
A. J. PERUMPANANI AND J. NORBURY
inflammation, placentation, and cancer metastasis [2]. In contrast, random motility is associated with the normal physiological processes of growth and development. Examples of simple random diffusion in the body include the nutrition of tissues, distribution of growth factors, and the respiratory exchange of gases. There has been much recent interest in the modelling of cancer metastasis both to understand the basic processes involved and to generate hypothesis that can be tested in the laboratory. Marusic et al. [3] and Perumpanani et al. [4] have developed models for tumour growth while McElwain et al. [5] have shown that cancer cells undergo passive motion towards the center of a multicell spheroid. Sherratt [6], Gatenby et al. [7], and Perumpanani et al. [8] have modelled various aspects of cancer invasion while Sherratt [9] and Albert et al. [lo] have studied the immune response to cancer. A closely related phenomenon which shares many phenotypic features with invasion is the process of tumour angiogenesis by which new blood vessels grow into a tumour; this has been modelled by Byrne et al. [l l] and Chaplain et al. [12]. This paper discusses a model for the invasion of normal tissue by a growing mass of cancerous tissue with particular focus on the numerical behaviour of the model equations. An interesting feature of this model is the interaction of the directed movement of the invading cancer cells (which results in a hyperbolic PDE) and the random motility of the proteases (which results in a parabolic PDE). A detailed derivation of the model and its biological implications are discussed elsewhere [8]; in this paper we focus on the numerical behaviour of this system of PDEs.
2. THE
MATHEMATICAL
MODEL
Briefly, the interaction of invasive cancer cells (~(z,t)) with the surrounding extracellular matrix (ECM, c(z,t)) through the liberation of a protease (p(z,t)) is modelled using a system of conservation equations. Here x and t represent one-dimensional space and time, respectively. The local proliferation of the cancer cells is described by the logistic growth model [13] subject to a local maximum (the carrying capacity) which has been scaled to unity. Invasive cancer cells have been shown to move up gradients of connective tissue [l], which is modelled as an invasive flux proportional to the ECM gradient .g in equation (1) below. The interaction of invasive cells with the ECM has previously been shown to induce production of proteases [14] and is modelled by a protease production term proportional to UC. Besides being naturally degraded (-krp), the proteases are also neutralized by the production of anti-proteases (Tissue Inhibitors of Metalloproteases (TIMPs)). The TIMPs are produced by the cellular elements of the connective tissue and at a lower rate by the invasive cells. Hence, the neutralisation of MMPs by TIMPs are modelled by two separate terms -kspc and -kGpu. Since the cells are the source of the TIMPs their concentration is taken as proportional to the concentration of the cells. Being small biochemicals the proteases are capable of diffusing into the surrounding tissue which we modelled as simple Pickian diffusion in equation (3). The extracellular matrix molecules which are fixed to the surrounding tissue are degraded by the action of the proteases which we modelled by the degradation term proportional to pc in equation (2). The complete model is now written as proliferation dU -=aat
directed motility ,-a-,
k2g
[ 1 UE
,
ac - = -kapc} ECM degradation, at 8P at=
* protease production
-~pc-kspll-k~~+~. protease degradation
8% diffusion
CONDITIONS. In the context of malignant invasion, we are studying a domain which is much larger than the size of the tumour, since the tumour is local&d to one part of the tissue.
BOUNDARY
Numerical Interections
125
In these simulations, we used a very large spatial domain so that the effects of the boundary were minimal. Zero Aux boundary conditions were used because the interest here is in the behaviour of an isolated system free from any external inputs. In our simulations we solved system (J)-(3) using step function initial conditions after nondimensionalisation. The total tissue density (u+c) through the whole domain (0 < z < L) was uniform. The invasive cells (u) were initially confined near z = 0, (u(z,O) = 1 for 0 c z < L/10, u(z, 0) = 0 elsewhere) and the remainder of the domain had ECM (c(z, 0) = 1 - u(z, 0)). There was no protease to start with (p(z, 0) = 0). INITIAL
CONDITIONS.
3. NUMERICAL
SOLUTIONS
We investigated the model in Section 2 to study the rate of invasion of the ECM (c) by the invasive cells (u) as a function of the system parameters (the subscripted ks). A typical simulation is shown in Figure 1, in which an advancing wave of invasive cells replaces a receding wave of ECM which is digested by an intervening pulse wave of protease. These simulations can be used to determine the parametric changes which will produce maximal inhibition of invasion, which can then be tested in laboratory assays as potential therapeutic modalities. For example, a decrease in k2 would correspond to therapeutic blockade of the haptotactic receptor which causes directed migration up fixed gradients. Because of the nonlinear nature of these equations the speed of the invading wave of cancer cells could depend on the parameters in a nonlinear fashion. In such a case, therapeutic blockade of the receptor could produce an unintended but potentially dangerous augmentation of invasion which can be predicted by modelling. The numerical solution of (l)-(3) is complicated by the interaction of the diffusional behaviour of the protease and the directed motility of the invasive cells. While equation (3) is parabolic because of the ksp,, term, equation (1) is a nonlinear convection equation where the convection coefficient is kgc,,. This combination of different types of PDE behaviour imposes different stability criteria on the numerical solutions. Furthermore, to get an accurate measure of the speed of invasion numerical diffusion should be kept to a minimum in the solution of equation (1).
4. STABILITY In order to numerically solve system (l)-(3), with upwinding of the form “+I3
-‘y
CRITERIA we first used an explicit finite difference scheme
= klUJ? (1 _ Uj” ) - k,&CS-JJJ’
At q+l - C? ’ = -k3PjnC;, At
- k$Yj’S;Cjn,
where UT, CT, and PT are the values of u, c, and p at the jth spatial point and the nth time step. Here 6, is the central difference operator (( UT+l - Ujn_i)/2Az) and S-, is the backward difference operator ((UT -UT-,)/A z ) w h ic h are used in accordance with the upwind scheme [15]; Ax and At are the spatial and temporal step sizes, respectively. Since (l)-(3) is a mixed parabolichyperbolic system of equations, the criteria for stability impose different constraints on the mesh and its refinement path. To illustrate this, consider the two equations representing the dispersion terms alone in (l)-(3)
where a is equivalent
ut + au= = 0,
(7)
Pt - Pzz = 0,
(8)
to cz in the full model and kz = ks = 1. For explicit
finite difference
126
A. J. PERUMPANANI
J. NORBURY
AND
0.9
0.8 0.7 0.6
0.9 0.8
0.6
p:: 0.3 0.2 0.1 n
Figure 1. Finite difference solutions of the system of PDEs (l)-(3) using the Hyperbolic-Parabolic Separation (HPS) method described in Section 4. Here the u equation (1) and the p equation (3) are solved on separate grids with similar spatial grid points. The time step (At) at each iteration is chosen so as to maximise the CFL number for the u equation and to satisfy the stability condition for the p equation. The solutions are temporally matched periodically with intervening values for u de termined by linear interpolation. The parameter values used in this solution are kl = 1.72, kg = 4.0, ka = 0.345, k4 = 5.0, kg = kg = 0.346, k7 = 0, ]cs = 0.5, and L = 30.
solutions the condition for stability of (7) is aAt/Ax I 1 (commonly referred to as the Courant, Friedrichs, and Lewy condition or the CFL condition) and for (8) is At/(Ao)2 5 0.5 [15]. When the two equations are solved simultaneously on similar grids, the (AZ, At) plane can be divided into four regions based on the stability criteria as shown in Figure 2. They are Region Region Region Region
1: 2: 3: 4:
parabolically unstable, hyperbolically stable; parabolically stable, hyperbolically unstable; both unstable, and both stable.
Numerical solutions using the scheme (4)-(6) s h owed that this is also true for the model system (l)-(3). Though numerical instability started in one of the equations, the system eventually failed as a whole because of the coupling between the equations. In order to understand the origins of the instability, and in particular, to distinguish instabilities arising in Regions 1 and 2 in Figure 2 we studied the solutions of (7)-(g) using Fourier analysis and compared the results with the numerical solution of (l)-(3). Standard Fourier analysis shows that p(s, t) = exp-ICZt sin kx satisfies differential equation (8) for all values of the number k. It can be shown that a similar Fourier mode is an exact solution of the difference equation arising from the finite difference approximation of (8). Substituting
Numerical Interactions
127 2At =Axj
Region 3
2-
Both unstable 1.5.
At
I-
Both Stable
0.5'
0.2
0.4
0.6
0.8
1
1.2
1.4
AX Figure 2. Shows the four regions in the (AZ, At) plane where the set of hyperbolic and parabolic equations described in (7) may be stable or unstable either individuaby or together. In Section 4, we discuss how hyperbolic and parabolic instability arising in a system of PDEk can be distinguished. n&jAz)
PT = (8
I where i = &i’,
into the finite difference approximation
p;+l At and dividing
of (8), i.e.,
PT _ Pjn+1 - 2P7 +#-I
-
(Az)~
(9)
’
by py shows that the Fourier mode is a solution for all values of n and j provided x z X(k) = 1+ V (@AZ - 2 + .c-ikAZ) = 1 - 4vsin2 :kAx,
where v = At/(Az)2 numerical
and X(k) is called the amplification
approximation
can be written
factor for the mode. For k = rnn the
in the form
(11) The low frequency terms in this expansion give a good approximation to the exact solution of differential equation (8), because the series approximation for X(k) and exp(-k2At) match
reasonably well. Under the restriction v 5 l/2 the above scheme is convergent and the error is bounded [15]. For large values of k the modes in the exact solution are rapidly damped by the exponential terms exp( -k2t). But in the numerical solution the damping factor (X(k)1 will become greater than unity for large values of k, if v > l/2; in particular this will happen when kAx = K, for then X(k) = 1-4~. These Fourier modes will then grow unboundedly as n increases. The above analysis is valid for a solution that admits all modes. However, on a discrete mesh there are only a finite number of distinct modes; modes with wave numbers kl and kz are indistinguishable if (kl - kz)A z is a multiple of 27~. Therefore, 21: can be expanded as a linear combination of the distinct modes corresponding to k = mn, where m = -(J - 1)) -(J 2) )...) -l,O,l,..., J and J = l/Ax. The highest mode that can be carried by the mesh has k = Jr or kAx = n; this mode has values fl at alternate points on the mesh. Equation (10) shows that this is also the most unstable mode for the difference scheme and has the amplification factor X(J?r) = 1 - 4~ and S[X] = 0. It is the fastest growing mode when u > l/2 and hence, eventually dominates the solution. However, in solving the hyperbolic equation (7) using the upwind
technique,
the Fourier analysis gives for the constant X E X(k) = l-
(I_
,+kAr)
a > 0 the amplification
= 1 - 5 (1 - e-ikA=)
,
factor
128
A. J. PERUMPANANI
where t = uAt/Ax
is the CFL number. lx(k)]2
AND J. NORBURY
This leads to
= ((1 -Q
+
+ [
(13)
= 1 - 4<(1 - I) sin2 :kAx.
It follows that IX(k) ] 5 1 for all k provided that 0 < [ 5 1. However, when 6 > 1, the highest frequency mode, kAx = X, or uj 0; (-l)j grows, but in this case the fact that S[X] =
-140
5I
I
10I
15
20
25.
301
35
40I
-14b
1
2
X
X
(4
(b)
3
4
I 5
a
9
I 10
2 2
0 -2 -4 -6 II
-a -10 -12 -140
I
5
I
10
15
20
25
30
35
40
6
7
X
X
(4
(4
Figure 3. Numerical solutions of the system of equations (l)-(3) using a finite difference explicit scheme showing two instances of failure. (a) illustrates failure because of the violation of the stability condition for diffusion, that is, ksAt/(A~)~ 5 0.5, where ks is the diffusion coefficient of the protease p. (b) shows a magnified view of the growing oscillations that causes the numerical failure which corresponds to the highest mode that can be carried by the mesh. (c) shows an instance of failure because the CFL number is too high (1.1) which appears similar to (a) but the magnified view in (d) shows that this has a lower frequency than the failure mode in (b).
Numerical
Interactions
5. TRUNCATION
129
ERROR
A suitable way to obtain stable solutions would be to restrict the mesh to the region in the (AZ, At) plane shown in Figure 2, where both the equations are stable. To do this we chose a fixed At and restricted the choice of Ax such that both solutions were stable. In Table 1, we show a numerical estimate of the wave velocity using such a strategy. Even for relatively fine temporal meshes. the velocity did not converge rapidly enough before the conditions for stability were violated resulting in numerical instability. To understand this more clearly we studied the truncation error of the numerical solution of equation (7). Table 1. A comparison of the computed wave speeds for the numerical solutions of (l)-(3) by different numerical schemes. When the temporal grid size At is held constant and the spatial grid size AZ is refined, failure occurs before the wave velocity has converged adequately. By holding the CFL number constant, better convergence is obtained, though instabilities again restrict the refinement of the spatial mesh. However, by solving the three equations in (l)-(3) using the high resolution scheme described in Section 6 on separate grids with different time steps and periodically matching up the solutions better convergence of wave velocities was obtained as shown under HPS (Hyperbolic-Parabolic Separation). For the intermediate values of u solved on a larger temporal grid than the p equation linear interpolation was used. Using the first-order upwind method required Ax = 0.007 to obtain comparable convergence. The parameter values used in this simulation are Icl = 1.72, Ic2 = 4.0, k3 = 0.345, kq = 5.0, kg = ks = 0.346, k7 = 0, and ]cs = 0.1. For the HPS solution we set the CFL number to 0.6.
Ax
At = 0.005
I 0.2 I
1.2201
1
CFL=0.3
HPS
1.1426
1 1.015
0.1
1.1784
1.1390
1.015
0.06
1.1582
1.1361
1.021
0.05
1.1426
1.1327
1.039
0.04
1.1327
1.1275
1.043
tn) t)
The truncation error is defined as the difference between the two sides of the difference equation for (7) when the approximation UT is replaced throughout by the exact solution u(zj, of the differential equation. At any point away from the boundary the truncation error T(x, is defined as T(x,t)
=
6+tu(x, At
t>
-
&r’11(x, Ax
t)
(14)
*
n
The Taylor series expansion may be used to express the truncation error as an infinite series so that u, - ;Axu,, +a .a j (15) * = f (Atu,, - aAzu,,) + . . . . Now even in the case where a is constant, so that utt = a2uz2, T(x, t) = -f (1 - <)aAzu,,
+ ... .
(16)
A. J.
130
PERUMPANANI
AND
J. NORBURY
Hence, the method is only first-order accurate and the truncation error depends on the CFL number c. In the numerical solution, when we changed the mesh size without keeping < constant, the CFL number changed at each refinement with a consequent change in truncation error explaining the poor convergence of the solution. Because of the above observation, we used a second approach which iixed the CFL number fairly low and used this fixed CFL number to determine the time step as At = CFL x Ax/a. The result of this numerical scheme is shown in Table 1. We found that for any particular choice of CFL number the spatial mesh can be refined up to a particular point at which the stability condition for the diffusion equation is violated causing the scheme to fail. The larger the CFL number chosen, the coarser was the spatial mesh size before which failure resulted. In the case demonstrated in Table 1, diffusional instability occurred even before the wave speed had satisfactorily converged. Hence, a scheme with a fixed CFL number was of limited value in solving the model in Section 2. Furthermore, it is important to choose CFL numbers as close to 1 as possible because low CFL numbers introduce large degrees of numerical diffusion as shown in equation (16). Hence, we needed to find a numerical scheme which would keep the CFL number fixed and close to 1, but would not violate the stability condition for the diffusion equation. A satisfactory solution to this problem was found by solving the hyperbolic and the parabolic equations on separate grids and taking the time steps so as to maximise the CFL number for the u equation and satisfy the diffusion stability condition for the p equation. Similar methods have been used in the past for as many time steps as necessary. The advantage of such an approach is that either of the time steps can be broken down into smaller steps should it be necessary to satisfy the stability criteria for that stage without prejudice to the other stage. Wu [lS] described such a method to study constant and variable flow velocities. Valocchi and Malmstead [17] have described the application of operator-splitting to a problem involving advection, dispersion, and reaction terms where they treated the reaction terms in a separate stage. We used a similar approach to solve the model in Section 2. The finite difference equation of (1) was iterated so as to satisfy the CFL condition while keeping the CFL number as closeto 1 as possible. Equation (3) was then iterated choosing a time step so that the diffusion condition was satisfied. For the range of parameters used, this typically meant that several time steps were needed for the p equation to match each time step of the u equation. The intermediate values were determined by linear interpolation. The fourth column of Table 1 shows the convergence of the wave velocity in the numerical solutions using the above strategy under the heading HPS (Hyperbolic-Parabolic Separation). The convergence of this scheme compared favourably with other schemes using a fixed time step or a fixed CFL number. A numerical solution using this scheme is shown in Figure 1.
6. HIGH
RESOLUTION
SCHEMES
First-order accurate schemes for hyperbolic equations suffer from numerical diffusion, but higher-order schemes, such as the Lax-Wendroff scheme [18] while giving higher resolution to discontinuities of the solution, exhibit spurious oscillations around such points. Much effort has been directed at obtaining second-order schemes which give high resolution while keeping spurious oscillations to the minimum. For example, Van Leer (191 and Roe [20] have proposed such high resolution schemes which incorporate terms that dampen spurious oscillations. Such damping terms are called “flux limiters”. This flux is the difference between the flux of a higher-order scheme and that of a low order scheme which has been “limited” in such a way as to ensure the resulting scheme is oscillation free (a term frequently used in the literature is “Total Variation Diminishing”). A high resolution scheme such as the Lax-Wendroff scheme [18] or the Warming and Beam scheme [21] may be used to solve (1). A flux limiter is introduced to dampen the oscillations in the wake of the travelling wave. We solved equation (1) using the Lax-Wendroff scheme with a
Numerical Interactions flux limiter
131
of the form
where vi = k&C~At/Ax and & = kzAt/Ax2 and the $js are the flux limiters. Setting 4 z 1 in (17) yields the Lax-Wendroff scheme and setting 4 = 0 gives the first-order upwind scheme. In order to calculate the flux limiters a quantity rj is first defined such that rj = A-/A+, the ratio of the backward and forward difference operators at each mesh point. Various flux limiters have been described by Van Leer [19], Roe [20], and several other authors. We used the minmod scheme [22] which gave the fastest convergence and defines 4 as
4j =
0,
rj 5 0,
Tj,
0
1,
Tj
5 1,
(18)
> 1.
Use of the minmod scheme in conjunction with HPS described in Section 5 resulted in a faster convergence of the computed wavespeed, though the end results were comparable to the first-order upwind scheme.
7. DISCUSSION The modelling of a large number of biological problems involve the features of directed and random motility which naturally bring together hyperbolic and parabolic equations. For instance, wound healing involves the directed migration of fibroblasts and the random motility of growth factors [23]. During angiogenesis, endothelial cells show chemotactic migration (directed migration along soluble gradients) up gradients of Tumour Angiogenesis Factors (TAFs), while TAFs themselves show random motility [24]. A similar combination is also seen during the emigration of leucocytes [25] and in the formation of a haemochorial placenta [26]. The numerical features mentioned in this paper will have a bearing on the modelling of such problems which study the migration of cells. The model aims at determining the rate of invasion of the cancer cells as a function of the parameters of the model. The rate of invasion was defined as the speed of the cancer cell wave which was measured numerically as the rate of translation of a fixed concentration on the u wave profile. However, the numerical solutions showed a sensitive dependence on the initial conditions. The system supports travelling waves of varying speeds depending on the rate of decay of the initial conditions. A slower rate of decay gave a faster wave. The slowest waves obtained corresponded to initial conditions with compact support which is in agreement with similar observations made with the Fisher equation [27] and the wave solutions of other nonlinear PDEs [5]. Since invading tumours correspond to initial conditions with a compact support, we . m. um wavespeed omin as the index of tumour invasiveness. measured the mmi A biologically interesting prediction made by the model is the dependence of the minimum wave speed omin on the diffusion coefficient of the protease ks. For extremely low values of ks, the minimum wave speed computed is in good agreement with the minimum wave speed obtained from the numerical solution of the purely hyperbolic system obtained with ks = 0, which has been described elsewhere [28]. As ks is increased the computed wave speed initially increased and then decreased. This feature is illustrated in Figure 4; we investigated it for a wide range of parameter values and found it to be a consistent feature of the model. Attempts are now underway in our lab to verify this prediction using an experimental array of collagen gel invasion assays. If confirmed, this prediction could explain the propensity of invasive cells to metastasise to particular tissues in the body where the diffusivity of the proteases they produce could make them maximally
A. J. PERU~PANAN~ AND J. NORBURY
132 1.2
t.15
I.1 cI m&l 1.05
I
0.95
0.9
0
0.2
0.3
0.6
0.8 k8
I
I .2
I.4
1.6
i.8
2
Figure 4. Computed minimum wave velocities a min for a finite difference solution of (l)-(3) for various values of the protease diffusion coefficient .Ics. As 4 is increased, amin initially increases and then decreases. The therapeutic implications of this finding are discussed in Section 7. The parameter values used in this simulation are Icl = 1.72, Icz = 4.0, Ic3 = 0.345, k4 = 5.0, ks = ]cs = 0.346, and k7 = 0.
invasive. The modelling of cell migration not only presents new numerical challenges because of the interaction of different types of motility, but also helps in the formulation of biological hypothesis which can be used to design novel experiments.
REFERENCES 1. S. Aznavoorian, ML. Stracke, H. Krutzsch, E. Schiffman and L.A. Liotta, Signal transduction for chemotaxis and haptotaxis by matrix molecules in tumour cells, J. of CelE Biology 110, 1427-1438, (1990). 2. W.G. Stetler-Stevenson, S. Aznavoorian and L.A. Liotta, Tumor cell interactions with the extracellular matrix during invasion and metastasis, Annu. Rev. Cell. B’iol. 9, 541-73, (1993). 3. M. Mar&c, Z. Bajzer, S.V. Pavlovic and J.P. Freyer, Tumour growth in tivo and as multicellular spheroids compared by mathematical models, Bull. Math. Biol. 56 (4), 617-631, (1994). 4. A.J. Perumpanani, J.A. Sherratt and J. Norbury, Mathematics modelliig of capsule formation and multinodularity in benign tumour growth, Nonlinearity 10, 1599-1614, (1997). 5. D.L.S. McElwain and G.J. Pettet, Cell migration in multicell spheroids-swimming against the tide, Bull. Math. Biol. 55, 655-674, (1993). 6. J.A. Sherratt, Chemottiis and chemokinesis in eukaryotic cells: The keller-segal equations as an approximation to a detailed model, Bull. Math. Biol. 56, 129-146, (1994). 7. R.A. Gatenby and E.T. Gawlinski, A reaction-diffusion model of cancer invasion, Cancer Research 5%,5?455753, (1996). 8. A.J. Perumpanani, J.A. Sherratt, J. Norbury and H.M. Byrne, Biological inferences from a mathematical model for malignant invasion, Invasion Metastasis 28, 209-221, (1996.). 9. J.A. Sherratt and M.A. Nowak, Oncogenes, anti-oncogenes and the immune response to cancer: A mathematical model, PFOC. R. Sot. Land 3 248, 261-271, (1992). 10. A. Albert, M. Freedman and A.S. Perelson, amours and the immune system: The effects of a tumour growth modulator, Math. Biosci. 50, 25-58, (1980). 11. H.M. Byrne and M.A.J. Chaplain, Mathe~tic~ models for tumour progeny-num~cal simul&tio~ and non-linear wave solutions, Bz~ll. Math. Biol. 57, 461-486, (1995.). 12. M.A.J. Chaplain and A.M. Stuart, A model mechanism for the chemotactic response of endothelial cells to tumour angiogenesis factor, Iii&A J. Math. Appt. Med. Biol. 10 (3), 149-168, (1993.). 13. V.G. Vaidya and F.J. Alexandro, Jr., Evaluation of some mathematical models for tumour growth, Int. J. Biomed. Comput. 13, 19, (1982). 14. B. Xie, C.D. Bucana and I.J. Fidler, Density-dependent induction of 92-kD type type-IV coilagenase activity in cultures of A431 human epidermoid carcinoma cells, Am. J. Pathol. 144(5), 1958-67, (1994). 15. K.W. Morton and D.F. Mayers, Nplmericnl Solution of Purtial Diflerential Equations, Cambridge University Press, (1994).
Numerical
Interactions
133
16. J.K. Wu, Wave equation model for solving advection-diffusion equation, Znt. J. of Numerical Methods in Engineering 37, 2717-2733, (1994). 17. A.J. Valochi and M. Malmstead, Accuracy of operator splitting for s&&ion-dispersion-reaction problems, Water Resources Research 5, 1471-1476, (1992). 18. P.D. Lax and B. Wendroff, Systems of conservation laws, Pure Appl. Math. 13, 217-237, (1960). 19. B. Van Leer, Towards the ultimate conservative difference scheme ii. Monotonicity and conservation combined in a second-order scheme, J. Comp. Phys. 14, 361-370, (1974). 20. P.L. Roe, Numerical algorithms for the linear wave equation, Technical Report 81047, Royal Aircraft Estab lishment, (1981). 21. R.F. Warming and R.M. Beam, Upwind second-order difference schemes and applications in aerodynamics, AZAA J. 14, 1241-1249, (1976). 22. P.K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal 21, 995-1011, (1984). 23. J.A. Sherratt and J.D. Murray, Model of epidermal wound healing, Proc. Roy. Sot. Lond. B 241, 29-36, (1990). 24. M.A. Chaplain and A.R. Anderson, Mathematical modelling simulation and prediction of tumour induced angiogenesis, Invasion Metastasis 16 (4-5)) 222-234, (1996). 25. T.A. Springer, Adhesion receptors of the immune system, Nature 346, 425434, (1990). 26. C.H. Graham and P.K. Lala, Mechanisms of placental invasion of the uterus and their control, Biochem. Cell Biol. 70, 867-874, (1992). 27. J.D. Murray, Mathematical Biology, Springer Verlag, (1990). 28. A.J. Perumpanani, J. Norbury, J.A. Sherratt and H.M. Byrne, A two parameter family of travelling waves with a singular barrier arising from the modelling of matrix mediated malignant invasion, Physica D, (submitted).