Serum uPAR as Biomarker in Breast Cancer Recurrence: A Mathematical Model

There are currently over 2.5 million breast cancer survivors in the United States and, according to the American Cancer Society, 10 to 20 percent of these women will develop recurrent breast cancer. Early detection of recurrence can avoid unnecessary radical treatment. However, self-examination or mammography screening may not discover a recurring cancer if the number of surviving cancer cells is small, while biopsy is too invasive and cannot be frequently repeated. It is therefore important to identify non-invasive biomarkers that can detect early recurrence. The present paper develops a mathematical model of cancer recurrence. The model, based on a system of partial differential equations, focuses on tissue biomarkers that include the plasminogen system. Among them, only uPAR is known to have significant correlation to its concentration in serum and could therefore be a good candidate for serum biomarker. The model includes uPAR and other associated cytokines and cells. It is assumed that the residual cancer cells that survived primary cancer therapy are concentrated in the same location within a region with a very small diameter. Model simulations establish a quantitative relation between the diameter of the growing cancer and the total uPAR mass in the cancer. This relation is used to identify uPAR as a potential serum biomarker for breast cancer recurrence.


Introduction
Human breast cancer is a major cause of death in the United States and worldwide [1]. It is estimated that 230,000 women in the United States are diagnosed annually with invasive breast cancer, and more than 40,000 die from the disease [2]. A major factor that contributes to poor prognosis is the fact that diagnosis is often delayed due to limitation in mammography screening [3]. Poor prognosis occurs also in assessing the risk of recurrence in patients of low grade breast cancer; improving this assessment will help avoid unnecessary chemotherapy [4].
Risk factors associated with gene mutations such as BRCA1 and BRCA2, and with family history and aging have long been recognized [5]. More recent work is also looking for risk assessment that can be associated with serum biomarkers [6][7][8]. Three tissue biomarkers have been identified: urokinase plasminogen activator (uPA), plasminogen-activator-inhibitor (PAI-1), and tissue factor (TF) [3,4,9,10]. For uPA to become active it must bind to its receptor uPAR [11]. Active uPA is extracellular matrix-degrading protease that promotes tumor progression and metastasis. It binds to plasminogen and converts it to its activated form, plasmin, a process inhibited by PAI-1 [12][13][14][15][16]. Plasmin mediates the activation of matrix metaloproteinase (MMP) which enables cancer cells' migration [12,15,17]. TF promotes tumor by enhancing VEGF production [18]. Harbeck et. al [19] reported on an extensive 6-year study to assess the risk associated with node-negative breast cancer recurrence in terms of the levels of uPA and PAI-1. Based on this report and other studies it was concluded that tissue (uPA, PAI-1) provide predictive information about early breast cancer [4,20]. The American Society of Clinical Oncology also recommends uPA and PAI-1 as prognostic tumor markers for breast cancer [21].
Although uPA and PAI-1 levels are elevated in breast cancer tissue, these high levels are not detected in the blood. Indeed, as reported in Rha et al. [22], the blood level of uPA and PAI-1 of the plasminogen activation system correlated with that of breast tissue in order of R 2 = 0.35 for uPA and R 2 = 0.11 for PAI-1. So uPA and PAI-1 are not reliable serum biomarkers. On the other hand, it was reported by Rha et al. [22], that the correlation of the level of uPAR in the blood with that of tissue is significant, with R 2 = 0.61 (P = 0.001). Recently Soydine et al. [23] found that uPAR in serum and in urine of breast cancer patients (n = 180) was significantly higher than in healthy control (n = 60). Serum uPAR was also shown to be a prognostic biomarker in endometrial cancer [24].
The present paper is concerned with prognosis of breast cancer recurrence. Most commonly, recurrence occurs within 3-5 years, although the statistics of recurrence is not clear [25,26]. When breast cancer recurs, it most often recurs in the same location as the primary cancer [27]. In the present paper we address with a mathematical model the following question: Can uPAR be used as biomarker to recognize breast cancer recurrence? We assume that after treatment of the original cancer, some cancer cells survived in the same location, occupying a small spherical region of radius R 0 and that the cancer begins to grow while maintaining a spherical shape. Simulations of our mathematical model profile the cancer radius and the expression of uPAR as functions of R 0 and time t: R(R 0 , t) and uPAR(R 0 , t). If a patient's uPAR in tissue (surrounding the primary tumor) is measured at time t, say at t = 100 days, we can then use this measurement to determine R 0 and hence also the tumor radius R(R 0 , t) for any time t after 100 days. The articles of Rha et al. [22] and Soydine et al. [23] suggest that the uPAR level in serum is siginificantly correlated and hence proportional to the level of uPAR in the tissue, hence serum uPAR could serve as a potential biomarker. When clinical data become available to more reliably confirm this proportionality coefficient, the uPAR could then actually be used as serum biomarker for breast cancer recurrence.

Model
The mathematical model is based on the diagram shown in Fig 1. The model includes, in addition to uPA, uPAR and PAI-1, also TF, VEGF, M-CSF, MMP and MCP-1. It also includes the cells that produce these proteins, or activated by them, namely cancer cells, fibroblasts, macrophages and endothelial cells. The variables of the model are listed in Table 1. The model is described by a system of partial differential equations (PDEs) in a radially symmetric tumor, with evolving radius R(t). The initial radius, R 0 = R(0), is a parameter which is patient-dependent. We shall assume that the total density of cells at each point in the tissue is constant; since tumor cells proliferate, the radius R(t) is increasing with time and cells are moving with velocity u which is time dependent.

Equation for tissue factor (T)
The tissue factor equation is given by where the second term is the production by cancer cells [28,29].

Equation for VEGF (V)
The evolution of VEGF concentration is modeled by The first term on the right-hand side accounts for production of VEGF by cancer cells [28,30,31], a process enhanced by tissue factor [18], and the second term accounts for VEGF produced by macrophages [30,31].

Equation for plasmin (P)
Plasminogen is the inactive precursor of trypsin-like serine plasmin. When it becomes activated, it is converted to plasmin. The generation of plasmin requires the binding of uPA to plasminogen, after uPA was released from the complex uPA-uPAR and became active [32], and this binding is inhibited by PAI-1 [12][13][14][15][16]. For simplicity we take the concentration of plasminogen to be constant. Hence the concentration of plasmin concentration satisfies the equation Equation for uPAR (u PR ) The uPA receptor uPAR is expressed by macrophages [11,33] and breast cancer cells [34]. Hence the equation of uPAR can be written as follows: Equation for uPA (u i P and u a P ) Inactive uPA is produced by fibroblasts [35]. Hence, uPA is activated when inactive uPA combines with uPAR. We take the equation for u a P to be @u a P @t À D u P Du a P ¼ l u u i P u PR K u PR þ u PR |fflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflffl ffl} production Àd u a P u a P |fflfflffl{zfflfflffl} degradation : ð6Þ The equation of the PAI-1 concentration is given by The first three terms on the right-hand side account for production of PAI-1 by cancer cells [11], fibroblasts [11,36] and macrophages [37].

Equation for M-CSF (q)
The M-CSF concentration satisfies the equation where the first term of the right-hand side account for production by cancer cells [38].

Equation for MCP-1 (p)
The equation of the MCP-1 concentration is given by The first term of the right-hand side accounts for production of MCP-1 by macrophages activated by M-CSF [31,39]. Here l p ðwÞ ¼ l p w w h if w < w h and λ p (w) = λ p if w > w h , where w h is an appropriate hypoxic level.

Equation for macrophages (M)
We assume that all cells are moving with common velocity u, and are subject to dispersion. The equation of macrophages density is given by where D M is the dispersion coefficient. The second term of the right-hand side accounts for chemotaxis [28,31,39]. Monocytes from the vascular system, with density M 0 , are attracted to the tissue by MCP-1, and they differentiate into macrophages. Macrophages are terminally differentiated cells. On the boundary of each blood vessel within the tumor there holds a flux condition, which we take to be As in [40] we can use a homogenization method to replace these fluxes by an average source of macrophages within the tissue, and this is the first term on the right-hand side of Eq (11). We assume that the macrophages are primarily tumor-associated-macrophages (TAM) of M2 phenotype.

Equation for endothelial cells (E)
Endothelial cells are chemoattracted by VEGF and their proliferation is increased by VEGF [31]. The equation for endothelial density is given by where we used the notation: The second term on the righthand side assumes a threshold V 0 below which proliferation of E does not occur [41,42].

Equation for fibroblasts (f)
There is a mutual enhancement in the interaction between cancer cells and fibroblasts. Cancer cells secrete TGF-β which increases the activation and proliferation of fibroblasts, while EGF secreted by fibroblasts increases the proliferation of cancer cells [43][44][45][46][47][48][49][50]. We write simplified equations for TGF-β (T β ) and EGF(G): The growth of cancer cells occurs on a time scale of days, whereas the secretion and decay of cytokines occur on a time scale of minutes to hours [51]. In order to understand the growth of cancer we simplify the model by using quasi-steady-state approximation for the equations of T β and E, so that T β = c 1 C, G = c 2 f for some constants c 1 and c 2 . Then the fibroblasts-enhanced growth rate of fibroblast density (f) through T β may be replaced by λ fC C, and the tumor-cell-enhanced proliferation of cancer density (C) through G may be replaced by λ Cf f. We use the Michaelis-Menten law to express the enhanced proliferation of fibroblasts by T β (i.e., by c 1 C) because TGF-β activation of fibroblast and EGF enhancement of cancer cells may be limited due to the limited rate of receptors recycling associated with this process.
The equation of fibroblast density is then given by

Equation for cancer cells (C)
The equation of the cancer cells density is given by @C @t There are three terms in the bracket: The first term is for tumor growth at rate which is oxygen dependent; the second term represents enhancement by EGF produced by fibroblasts; the third term accounts for proliferation of cancer cells by uPA as it binds to its receptor uPAR on cancer cells [11,12,15,36], a process resisted by PAI-1 [14,16]. The Michaelis-Menten law is used to represent the limited rate of receptor recycling associated with the enhancements by f and by u PR . We take

Equation for oxygen (w)
The oxygen concentration evolves according to the equation The first term of the right-hand side accounts for infusion of oxygen through the blood, which is represented by the density of endothelial cells. The last three terms represent oxygen taken up by macrophages, tumor cells and fibroblasts.

Equations for MMP (Q) and TIMP (Q r )
We have the following sets of reaction-diffusion equations for MMP and TIMP: MMP and TIMP are activated by macrophages [52,53] and MMP activation is enhanced by plasmin P [17], and MMP is lost by binding with TIMP, while TIMP is depleted as it blocks MMP [38,54,55].

Equation for ECM (ρ)
Extracelluar matrix (ECM) is produced by fibroblasts [40], and is degraded by MMP [38,54]. The equation for the density of ECM is given by Equation for u We assume that the total density of all the cells plus the density of ρ is constant: We also assume that all cells are approximately of the same volume and surface area, so that the dispersion coefficients of the all cells have the same coefficient, D. By adding Eqs (11)-(14) and (17), we get an equation for rÁu: We can, conversely, derive Eq (18) from Eqs (11)- (14) and (19). We assume that if a breast cancer recurs, it is because some cancer cells survived in the initial location. We also assume, to simplify the computations, that these cells are contained in a sphere of small radius R 0 , and that as time increases the region containing the cancer cells is a sphere with increasing radius r = R(t), and all the model variables are radially symmetric, that is, they are functions of (r, t) where 0 r R(t), t > 0. In particular, jxj is the radial unit vector.

The free boundary equation
We assume that the tumor is spherical with radius r = R(t), so that where e r ¼ x jxj is the radial unit vector.

Boundary conditions
We assume flux boundary conditions due to oxygen transport of the form for X = w and E, where X 0 is the density of oxygen or of endothelial cells within the vascular system at the boundary of the tumor. The coefficients α X on the free boundary are chosen as follows: We assume no-flux boundary condition for all other variables except ρ. Since ρ satisfies the hyperbolic Eq (17) with ue r as the velocity of the free boundary, we do not need to prescribe a boundary condition for ρ on the free boundary. The radial velocity is determined from Eq (19).

Initial conditions
We choose the following initial conditions for ECM, oxygen and cells concentration in unit of g/cm 3 : ρ = 0.02, w = w h , M = 0.07, E = 0.05, f = 0.06, and C = 0.4. All the cytokines are taken to be initially zero. The initial conditions do not affect the simulation results after a few weeks. However, as will be shown, the very small value of R(0) plays a critical role in the tumor growth profile.

Results
In Figs 2 and 3 we take R(0) = 10 −2 cm = 100μm; at this initial size the tumor may contain a few thousand cells, including cancer cells. We define the average density/concentration (Ave) and total mass (TM) as follows  hypoxic conditions as the density of endothelial cells is still small (VEGF is still small). The average of oxygen is decreasing as cancer cells are increasing and consuming more oxygen, while the consumption of oxygen by macrophages, fibroblasts, and epithelial cells, combined, is also slightly increasing. uPA is growing after 80 days when fibroblast density begins to increase; PAI-1 mimics the profile of active uPA, while the concentrations of uPAR is relatively stable.
In estimating some of the parameters we assumed steady state of averages of densities of cells and concentrations of proteins, except for p and w. Harbeck et al. [20] measured the concentrations of uPA and PAI-1 in breast cancer survivors by the protein nanogram of protein antigen per milligram of tissue protein of breast tissue. They determined that when the risk of recurrence is high, PAI-1%5×uPA. This proportion between PAI-1 and active uPA is seen also in Fig 2. For the purpose of developing serum biomarkers we are interested in the total mass of uPAR in the growing tumor, rather than in the average density. Fig 3 shows that the total mass of cells and cytokines are increasing in time, but at different rates.
So far we have taken R(0) fixed at R 0 = 10 −2 cm. We next want to use the model for diagnostic purposes. Our goal is to determine from one measurement of the total mass of uPAR at time t 0 after the primary breast cancer treatment the radius of the tumor at the same time t 0 and at any subsequent time t 1 , t 1 > t 0 . To do this we take R 0 2 (10 −2 , 5 × 10 −2 ) cm as a potential  Tables 2 and 3. initial radius. For each R 0 , we compute the tumor radius at time t, R = R(R 0 , t) and the total mass of uPAR at time t, uPAR(R 0 , t). Fig 4 shows R and Fig 5 shows the total mass of uPAR, at any (R 0 , t) in the range R 0 2 (10 −2 , 5 × 10 −2 ) cm and 0 < t < 1000 days. From these two figures we can generate a mapping from total uPAR(t) to R(t) which is independent of R 0 , as follows: From a value of uPAR at some time t after primary breast cancer treatment, we estimate, by using Fig 4, the corresponding parameter R 0 . We can then use Fig 5 to predict the value of R corresponding to this R 0 and the time t. The mapping is uPAR(t)!R(t) is shown in Fig 6, where the color bar determines the value of R(t) for a given pair of uPAR and t.
The color map in Fig 6 is a prognostic map for recurrent breast cancer: when a patient's uPAR is measured t days after the primary breast cancer treatment, the color bar in Fig 6 predicts the size of the recurrent tumor.

Discussion
There are several mathematical models of breast cancer focusing on different aspects of the disease [28,38,[56][57][58]. Our model is the first one to focus on biomarkers associated with the risk of breast cancer recurrence.
Cancer recurrence occurs in 10 to 20 percent of all breast cancer survivors. In order to avoid unnecessary radical treatment, it is important to diagnose a recurrent cancer as early as possible. Although several tissue biomarkers have been identified, biopsy cannot be frequently repeated. This motivated the present study of focusing on potential serum biomarkers which are non-invasive. Of all the plasmiogen system tissue biomarkers only uPAR concentration significantly correlates with uPAR concentration in blood [22,23]. For this reason we developed in this paper a mathematical model that quantifies the relation between tissue uPAR and the size of a recurrent cancer (Fig 6).
The mathematical model is represented by a system of partial differential equations in a growing tumor with radius R(t). We assume a very small initial radius R 0 2 (10 −2 , 5 × 10 −2 ) cm, corresponding to cancer cells that survived after primary breast cancer treatment. The radius R 0 may vary from one patient to another. Nevertheless, Fig 6 shows that a patient's uPAR measured at any time within 1000 days after primary breast cancer treatment, can be used to estimate the cancer size (i.e., its radius) by the color bar. Furthermore, the initial R 0 can then be determined from While the measurement of serum biomarkers for patient survival after primary breast cancer treatment is still being debated [8], we propose here uPAR as a potential serum biomarker. When clinical data become available to enable us to estimate the precise proportion of uPAR tissue concentration to plasma concentration, uPAR could then be used as plasma biomarker which informs the size of a recurrent tumor. The simulation results of Figs 4-6 can be extended to include smaller values of R 0 , e.g. 0.005 R 0 0.01. For such values, it takes a significantly longer time for the tumor radius to reach a size that can be detected by self examination. Color map for the total mass of uPAR. R 0 ranges from 0.01 to 0.05 cm and t ranges from t = 0 to t = 1000 days. Color represents the total mass of uPAR. All the parameters are as in Tables 2 and 3.
The simulation results rely heavily on parameters taken from the literature, sometimes in a different context, and on parameters estimated in this paper, sometimes rather crudely. For this reason we performed sensitivity analysis to determine how the radius of the tumor will vary when some of the parameters are increased or decreased.
In developing the mathematical model we made several simplifying assumptions: (i) the tumor is radially symmetric; (ii) the total density of cells at each point of the tissue is constant; and (iii) the cells are approximately of the same volume and surface area making the dispersion coefficient of cells equal. Furthermore, the mathematical model is minimal in the sense that includes just the plasminogen system (of uPA, uPAR, PAI-1), plasmin and MMP, tissue factor, VEGF, MCP-1, M-CSF, and the cells that produce these proteins or are activated by them, namely, macrophages, endothelial cells, fibroblasts and cancer cells. For these reasons, this work should be viewed as providing just a first step upon which a more elaborate and more comprehensive model could be developed in the future. When new data of uPAR expression level in patients of breast cancer recurrence become available, some of the parameters of the model will accordingly be adjusted to make the model simulations agree with patients data. Additional cytokines may then also be included, and the assumption of a spherical tumor may be changed to better reflects tumor histology.

Diffusion coefficients
The diffusion coefficients of proteins (Y) are proportional to the molecular surface area [54], which is proportional to M 2=3 Y , where M Y is the molecular weight [54]. Accordingly, we can have the following relation: where M V and D V are the molecular weight and diffusion coefficient of VEGF. Since D V = 8.64 × 10 −2 cm 2 day −1 [59], M V = 24 kDa [60] and M P = 92 kDa [61], M u P = 38 kDa [56], Eq (1) • d T : The tissue factor half life is 9 hours [64], hence d T = 1.85 day −1 .
• A T : The concentration of tissue factor in cancer is 3.5 × 10 −5 g/ml [65]. The ratio of tissue factor in healthy to disease in plasma is 1:2 [66]. We assume larger ratio in breast tissue, so that in the healthy case T = 2 × 10 −4 g/ml and take K T = 10 −4 g/ml. From Eq (1), in steady state for the healthy case, we get A T = d T T = 3.23 × 10 −5 g/ml/day.
• λ TC : We assume that most of the tumor is populated with cancer cells, and take C = 0.5 g/ml. As in the deviation of Eq (12), to simplify the model and to estimate some of the parameters, we assume that cytokine equations are in steady-state. From the steady-state of Eq (1) in disease state we have A T + λ C C = d T × 3.5 × 10 −5 , so that λ TC = 5.7 × 10 −5 day −1 .

Eq (2)
According to (31,39) λ VC varies in the range of 10 −21 to 10 −20 in units of g/s/cell, and the production of VEGF by tumor associated macrophages is far larger than the production of VEGF by cancer cells. Accordingly we take λ VC = 2 × 10 −8 /day and λ VM = 2 × 10 −6 /day. We also assume that T enhances cancer-cells production by less than 200% and take λ VT = 2.
• λ P : The molecular weight of plasminogen is 92 kDa [61] while molecular weight of uPA is 38 kDa [56]. The concentration of uPA in breast cancer is 1.8 × 10 −6 g/ml [70,71]. We assume that the number of uPA proteins in plasma is the same as the number of plasmin proteins, so that the concentration of plasmin the concentration of uPA ¼ 92 kDa 38 kDa : Therefore, the concentration of plasmin in breast cancer is approximately 4.4 × 10 −6 g/ml. By Eq (3), we have which implies λ P = 2.42 × 10 −6 g/ml/day.
• λ u PR M and λ u PR C : The concentration of uPAR in breast cancer is 1.8 × 10 −6 g/ml [70,71]. The concentration of uPAR in normal healthy breast is significantly smaller [73]; we take it to be approximately 10 times smaller. In healthy case, M = 0.04 g/ml [74], and from l u PR M M À d u PR u PR ¼ 0 we get, λ u PR M = 6.21 × 10 −6 day −1 .
Eqs (5) and (6) • d u i P and d u a P : The half life of activated uPA is 5 hours [76,77], therefore d u a P ¼ 3:2/day. We assume that inactivated uPA degrades slower, and take d u i P ¼ 2:4/day. • λ uf : The concentration of uPA in normal healthy tissue is estimated to u i P ¼ 6 Â 10 À6 g/ml [78], and fibroblast density is 0.07 g/ml [40]. From we get λ uf = 2.057 × 10 −4 day −1 .

Eq (12)
In Eq (12) production term l fC f C K C þC is due to cytokines secreted by cancer cells. We assume that this term is only a fraction of the death rate d f f of fibroblasts, where d f = 16.6 × 10 −3 /day, and take λ fC = 5 × 10 −3 /day.

Eq (13)
λ C (w) accounts for the proliferation rate minus the death rate by necrosis, while d C is the death rate by apoptosis. In transgenic mice λ C (w) is large [31,38], and cancer develops within a few days. Since breast tumor in human develops much slower, on the time scale of years, we take l C ðwÞ ¼ 0:6 w w h day −1 if w < w h and λ C (w) = 0.6 day −1 if w > w h , while d C = 0.5 day −1 [31]. We assume that the enhanced proliferation rate by fibroblast and by u a P binding to u PR are small, and take λ Cf = 0.06 day −1 and λ Cup = 0.05 day −1 .

Eq (14)
We assume that fibroblasts consume oxygen at the same rate as macrophages, so that d wf = d wM = 80 cm 3 /g/day by [31].

Eq (15)
From the steady state of Eq (3), with no activated uPA, we get P ¼ l P d P . Since d P = 1.39/day [67] and λ P was estimated (in Eq (3)) by 2.42 × 10 −6 g/cm 3 /day, we get P = 1.74 × 10 −6 g/cm 3 . With active uPA in Eq (3), this value of P should be increased by the factor 1 þ l P u u a P K P A þP A , so accordingly we take the half-saturation K P to be 4.4 × 10 −6 g/cm 3 .

Boundary conditions
Since most VEGF is produced by tumor associated macrophages [31,39]. the steady state of Eq (2) yields where ε < 1, d V = 12.6/day [31] and λ VM = 2 × 10 −6 /day. Taking M = 0.3 g/ml and ε = 1/2, we get the approximate value V = 7 × 10 −8 g/ml. We accordingly take K V = 7 × 10 −8 . We also take e a w and e a E to be of order 1, and for simplicity choose e a w ¼ e a E ¼ 1; however, other choices affect the model simulations only very little (not shown here).

Sensitivity analysis
We performed sensitivity analysis on the all production parameters of the system (1)- (17). Following the method in [81], we performed Latin hypercube sampling and generated 1000 samples to calculate the partial rank correlation (PRCC) and the p-values with respect to the radius of the tumor at day 600. The results are shown in Fig 7 (The p-value was <0.01).
The most positively correlated production parameters are λ VM (the production of VEGF by macrophages), λ u (the production uPA activator), λ Cf (λ Cf and λ fC represent the cross-talk between cancer cells and fibroblasts, which increases the number of cancer cells). The most negatively correlated production parameters are λ Pf and λ PC (which increase, together with λ PM , the production of PAI-1, thus increasing the blockage on uPA and the consequently proliferation of cancer cells), and λ P (which increases plasmin, and hence also PAI-1).
The remaining parameters are mildly correlated to tumor growth, and their correlation (+ or -) is agreement with the model description in Fig 1.   Fig 7. The sensitivity analysis for the cytokine production rates. The figure shows the partial rank correlation (PRCC) between the cytokine production rate and the radius of tumor. All the parameters are as in Tables 2 and 3. doi:10.1371/journal.pone.0153508.g007

Computational Method
In order to illustrate our numerical method, we consider the following convection-diffusion equation: where F X accounts for all the 'active' terms. Since the model we consider is a free boundary problem, we employ the moving mesh method to compute it. We write Eq (33) can be written in the total derivative form dXðrðtÞ; tÞ dt þ divðvÞX ¼ Dr 2 XðrðtÞ; tÞ þ F X : Let r n i , X n i denote numerical approximations of i-th grid point and Xðr n i ; tÞ, respectively, when t = nτ, where τ is the time stepsize. The discretization is derived by the explicit Euler finite difference scheme, i.e., where X r ¼ h 2 À1 X n iþ1 Àh 2 1 X n iÀ1 Àðh 2 1 Àh 2 À1 ÞX n i h 1 ðh 2 À1 Àh 1 h À1 Þ , X rr ¼ 2 h À1 X n iþ1 Àh 1 X n iÀ1 þðh 1 Àh À1 ÞX n i h 1 ðh 1 h À1 Àh 2 À1 Þ , and h À1 ¼ r n iÀ1 À r n i , h 1 ¼ r n iþ1 À r n i . Then the mesh is moving by r nþ1 where v n i is solved by the velocity equation. In order to make the Euler method stable, we take t minfh 1 ;h À1 g 2D .