Exosomal miRs in Lung Cancer: A Mathematical Model

Lung cancer, primarily non-small-cell lung cancer (NSCLC), is the leading cause of cancer deaths in the United States and worldwide. While early detection significantly improves five-year survival, there are no reliable diagnostic tools for early detection. Several exosomal microRNAs (miRs) are overexpressed in NSCLC, and have been suggested as potential biomarkers for early detection. The present paper develops a mathematical model for early stage of NSCLC with emphasis on the role of the three highest overexpressed miRs, namely miR-21, miR-205 and miR-155. Simulations of the model provide quantitative relationships between the tumor volume and the total mass of each of the above miRs in the tumor. Because of the positive correlation between these miRs in the tumor tissue and in the blood, the results of the paper may be viewed as a first step toward establishing a combination of miRs 21, 205, 155 and possibly other miRs as serum biomarkers for early detection of NSCLC.


Introduction
Lung cancer is the leading cause of cancer-related deaths in the United States and worldwide, and non-small cell lung cancer (NSCLC) constitutes 85% of lung cancer deaths [1,2]. Five years survival rate for NSCLC is significantly higher for those diagnosed at early stage [3], but there are no reliable tools for early detection of lung cancer. Most lung cancers are first diagnosed on symptoms. Approximately 10% of patients present brain metastasis at the time of initial diagnosis and their mean survival is 4 months [4]. Hence, there is a need for novel noninvasive biomarkers for early lung cancer diagnosis [5].
Exosomes are nano-vesicles of size 30-100 nm in diameter, surrounded by a lipid bilayer, and containing fuctional proteins, mRNAs and microRNAs (miRs). Exosomes are released by various cells, including cancer cells [6]. A growing body of evidence suggests that exosomal miRs may be used as serum biomarkers for prognosis of malignant tumors [5,7]. Furthermore, exosomal miRs inhibitors have been evaluatedas anti-tumor drugs in experimental and clinical work for several types of cancer, including lung cancer [8,9]. and caspase 9, forming an apoptosome, which leads to apoptosis through activation of caspase 3 [26,27].
In NSCLC, the most expressed exosomal miRs are miR-21, miR-155 and miR-205 [10]. In  [8,25], and thus promotes activation of the MAPK and AKT pathways. Also, miR-21 and miR-205 block PTEN [28,29] and thus promote the activation of the AKT-mTOR pathway. Fig 2(b) shows the effect of overexpression of miR-155. MiR-155 blocks Apaf-1 expression [30] and thus also cellular apoptosis when DNA damage occurs. Hence overexpressions of miR-21 and miR-205 give rise to uncontrolled proliferation, while overexpression of miR-155 leads to reduced apoptosis. MiR-21 and miR-205 have also other targets; in particular it was suggested that miR-21 targets tumor suppressors involved in apopotsis, including Apaf-1, Pdcd4, RhoB and Faslg [31,32]; hence overexpression of miR-21 reduces apoptosis. However, for simplicity, we focus in this paper on what seems to be the main targets of miR-21 and miR-205 in NSCLC, as shown in Fig 2. In this paper we consider growth and invasion of lung tumor associated with mutations in EGFR, MAPK and AKT, and its treatment by anti exosomal miRNAs (miR-21, miR-205 and miR-155). We use the mathematical model to determine the efficacy of these drugs under these mutations.
We consider two aspects of tumor progression: (i) Invasion, in which a tumor planar front progresses, in time, away from the main body of the tumor; (ii) Proliferation, in which a small spherical tumor grows in time. In order to focus on the role of the exosomal-miRs, we do not include in the model the immune responses and angiogenesis; thus the model represents an early stage of lung cancer.

Mathematical model
The mathematical model is based on the network shown in Fig 2. For simplicity we use just one variable, MAPK, to represent the Ras-Raf-MEK-ERK signaling pathway, and AKT to represent the PI3K-AKT signaling pathway. We also combine miR-205 with miR-21 in modeling their effect on blocking PTEN. Table 1 lists the variables used in the mathematical model in unit of g/cm 3 .
Equations for proteins As in [33], the dynamics dỹ dt ¼f ðỹÞ of the proteins within cancer cells will appear in the form where C 0 is the steady state density of cancer cells.

Equation for EGF-EGFR (E). The equation for EGF-EGFR is given by
The coefficient λ E is the production rate of the EGF-EGFR complex and the factor 1/(1 + M/ K ME ) is the inhibition by ERK [14]; d E is the degradation rate of E.

Equation for MAPK (M).
The MAPK pathway is activated by the EGF-EGFR [13,14], a process resisted by TKI [8,13]. Hence where d M is the degradation rate of M.

Equation for AKT (A).
The activation of the AKT pathway initiates with the activation of PI3K (P 3 ) by EGF-EGFR directly and also through Ras (which is activated by EGF-EGFR) [15, 16,19,20]. In view of the TKI inhibition of EFG-EGFR [24,25], the equation for P 3 takes the form where S is the concentration of Ras. AKT is activated by PI3K which is negatively regulated by PTEN [15,25], so that We assume that the turnover of PI3K is very fast (the half-life of PI3K is very short [34]) and deduce from the steady state equation of P 3 that Substituting this into the equation for AKT, we obtain Assuming also that the concentration of Ras is proportional to that of MAPK, i.e. S = μM, we obtain the following equation for AKT:

Equation for TKI (T).
The production of TKI is inhibited by miR-21 [8,25]. Recalling that miR-21 is only a fraction of m 1 , we write the equation for TKI in the form Equation for PTEN (P). The expression of PTEN is inhibited by both miR-21 and miR-205 [28,29]. Hence Equation for Apaf-1-caspase 9 apoptosome (Ap). The expression of Apaf-1 is down-regulated by miR-155 [30]. Hence Equation for exosome (E C ). Cancer cells shed exosomes at a rate λ Ec C. We assume that exosomes are degraded, releasing their miRs, when merging with cancer cells. Taking the rate of this degradation to be d Ec E C Á C K C þC , the equation for the concentration of exosomes is given by where the term D Ec ΔE C represents dispersion (or diffusion) of exosomes. Equations for exosomal miR-21 (m 1 ) and exosomal miR-155 (m 2 ). MiR-21 and miR-155 are released from exosomes when exosomes merge with cancer cells. We take the exosomal production rate of miR-21 to be λ m 1 E C Á C/(K C + C), and obtain the equation Similarly, the equation for miR-155 is given by Equations for miR-21 (m i 1 ) and miR-155 (m i 2 ) in cancer cells. Since m i 1 and m i 2 lie in cancer cells, they diffuse with same coefficient as cancer cells. Hence, the equations for m i 1 and m i 2 are given by and @m i respectively. We will apply the mathematical model to consider two phases of tumor progression of lung cancer, invasion and proliferation.

Equation for cancer cell (C) in tumor invasion.
The equation for cancer cell is the following: Invasion of cancer cells is driven by competition for space and resources [35,36]. At the early stage of tumor invasion resources are not limited, hence cells undergo migration in the direction of decreased gradient of cancer cells density. On the left-hand side of Eq (12) the term χr Á (CrC) represents the directed migration of cancer cells in response to the competition for space; χ is the 'directed migration coefficient'. We assume a logistic growth with rate which depends on both MAPK and AKT, since both pathways lead to cell replication; this accounts for the first term on the right-hand side of Eq (12). In addition to natural apoptosis, at rate d C C, damage to cancer cells (at rate proportional to d D C) leads to apoptosis by formation of the Apaf-1-caspase 9 apoptosome [26,27]; this is accounted for by the second term on the right-hand side of Eq (12).

Boundary and initial conditions for tumor invasion model.
We assume that a solid tumor lies in the the half plane x < 0, and model the progression of the tumor front in the direction of increasing x. We assume that the tumor front is planar, and that it moves in the interval 0 x 2 from x = 0 cm towards the end-point x = 2 cm. We impose the following boundary conditions for the cancer cells, exosomes and miRs: Since initially the cancer is confined to x < 0, we take zero initial conditions: Model for tumor proliferation

Equations for cancer cells (C) and normal healthy cells (N) in tumor proliferation.
In the tumor invasion model, the directed migration coefficient χ represents the directed movement of the invading tumor cells. The precise range of the parameter χ is unknown. In order to visualize significant advance of the migrating tumor front we took χ to be the 'relatively' large. However, in the tumor proliferation model, proliferating cells grow faster than migrating cells [37][38][39], and the competition for space is primarily a competition with normal healthy cells (with density N) [40]. We therefore assume that the directed movement of cells is determined by the condition that the total cell density, C + N, is constant at each point in the tumor. The term −χr Á (CrC) in Eq (12) is then neglected and replaced by the term r Á ðuCÞ, where u represents the velocity of cells. The equation for cancer cells is given by where the competition for space with the normal healthy cells is represented by the term εN in the logistic growth. We assume that most exosomes shed by cancer cells release their content when they make contact with nearby cancer cells, and therefore keep Eq (7) unchanged. The equation for normal healthy cells, N, is given by The competition for space with cancer cells is represented by the term εC in the logistic growth term [40].
To simplify the computations, we assume that the tumor is spherical and denote its moving boundary (i.e. its radius) by r = R(t). We also assume that all the densities and concentrations are radially symmetric, that is, functions of (r, t), where 0 r R(t). In particular, u ¼ uðr; tÞe r , where e r is the unit radial vector.

Results for tumor invasion
In simulating the invasion of cancer cells we use the model Eqs (1)-(12), with boundary Conditions (13) and initial Conditions (14), and with the parameters of Tables 2 and 3. We explore how specific mutations affect the invasion of the tumor front, and how anti-miR drugs slow the invasion. We consider four cases: (i) the control case (with unspecified mutations), (ii) new mutation in EGFR, (iii) new mutation in MAPK, (iv) new mutation in AKT. In the control case all the parameters are taken to be the same as Tables 2 and 3, and χ is taken to be 3 × 10 −2 cm 5 g −1 day −1 . In the case of mutations in EGFR, MAPK or AKT, χ is unchanged but the production rates λ E , λ M and λ A are increased by some factor. The first row of Fig 3 shows the spatial profile of cancer cell density C(x, t) in the control case, and in the  three cases of mutations (in EGFR, MAPK and AKT) at different time points t = 5, 15, 30, 60 days. We see that under each of the three mutations the tumor advanced at day 60 by at least 10% more than in the control case. We note however that although the sizes of the invasion under the three mutations are nearly the same, we accounted for the three mutations differently, increasing the production rates of EGFR by a factor 1.3, of MAPK by a factor 1.6 and of AKT by a factor 1.8. The ratios between these factors suggest that a mutation of EGFR increases tumor invasion more than a mutation of MAPK, and a mutation of MAPK increases tumor invasion more than a mutation of AKT. These suggestions, however, need to be verified experimentally.
The second row of Fig 3 shows the effect of anti-miR-21 drug in the control case and in the cases of EGFR, MAPK and AKT mutations. We note that anti-miR-21 reduces the rate of invasion by approximately 17%. When both anti-miR-21 and anti-miR-155 are combined, the reduction is by 40%, as seen in the third row of     Cancer invasion depends on the directed migration coefficient χ. In [33] the range of the parameter χ was taken to be (3 × 10 −4 , 3 × 10 −2 )cm 5 g −1 day −1 . In the simulations of Figs 3-5, we took the largest value χ = 3 × 10 −2 cm 5 g −1 day −1 in order to visualize the invasion of the tumor front over a relatively short period of time. It is reasonable to expect that both tumor invasion and total mass will decrease if χ is decreased. This is illustrated in Fig 6 in the case of a tumor with the same EGFR mutation as in Figs 3 and 4. We denote by R χ the distance traveled by the tumor front by day 60, and by M χ the total linear mass of the cancer cells by day 60. Fig 6 shows the growth of R χ and M χ (at day 60) as χ increases from 3 × 10 −4 to 3 × 10 −2 cm 5 g −1 day −1 : R χ increases by a factor 8 and M χ increases by a factor 11 approximately.
We next apply anti-miR-21 and anti-miR-155 drugs to the tumor (by reducing λ m 1 and λ m 2 to λ m 1 /2 and λ m 2 /2, as in Figs 3 and 4) and denote the corresponding R χ and M χ by R m w and M m w . We represent the efficacy of the anti-miRs drugs by 0 R :¼ ðR w À R m w Þ=R w and 0 M :¼ ðM w À M m w Þ=M w , that is, by the percentage of reduction in R χ and M χ . Fig 7 shows that the efficacy of the drug increases as the directed migration coefficient χ increases. The efficacy f R is approximately 33% when χ = 3 × 10 −4 g −1 day −1 , and it increases to 40% when χ = 3 × 10 −4 g −1 day −1 . The efficacy f M is 57% when χ = 3 × 10 −4 g −1 day −1 , and it increases to 68% when χ = 3 × 10 −4 g −1 day −1 .
From Figs 6 and 7 we conclude that as χ increases, the tumor invasion and total mass increase, while at the same time the efficacy of anti-miRs drug also increases. The same results (not shown here) hold for other mutations as well as for the control case.

Results for tumor proliferation
The simulations of proliferation of cancer cells are based on the model Eqs (1)-(8), (16)- (18) with boundary Conditions (19) and initial Conditions (20). We increase both λ C1 and λ C2 by a factor 1.4 compared to the values in tumor invasion model in order to account for the fact that proliferating cells grow faster than migrating cells [37][38][39]. We also increase the steady state C 0 from 0.4 g/cm 3 to 0.46 g/cm 3 to reflect the fact that invading cancer cells have sparser density than proliferating cells. We take the steady state density of healthy cells, N 0 , to be 0.14 g/cm 3 so that Eq (17) holds. All the other parameter values are the same as in Tables 2 and 3. We take the initial tumor radius to be R(0) = 0.01 cm. Fig 9 shows the average concentrations of all the variables over a period of 60 days. Most of the concentrations are either monotone increasing or monotone decreasing in time: the cell growth inhibitors TKI, PTEN and Apaf-1 are decreasing, while the cell growth promoters are increasing. The only exception is the average density of E. It is initially increasing since MAPK density is small. But MAPK continues to increase (as T keeps decreasing), and after a few days the inhibition by MAPK (or actually ERK, see Fig 1) forces E to decrease, and it does so until it reaches a steady state.
We note that in estimating some of the parameters of the model equations we assumed steady-state of the various variables (cells, proteins, miRs). The steady state of the variables in Fig 9 agree approximately with those steady state values, and this establishes consistency of our assumed steady-state values. In particular, the average density of cancer cells stabilize at 0.4631 g/cm 3 , and the average density of normal healthy cells stabilize at 0.1337 g/cm 3 , while C + N remains approximately equal to 0.6 g/cm 3 at the entire time.

Treatment
It is well known that cancer cells in NSCLC lose sensitivity to anti-tumor drugs, for example to paclitaxel, gefitinib and cisplatin, and that some anti-miRs can restore some of this sensitivity. We use our model to explore the effect of anti-miR combined with paclitaxel, gefitinib and cisplatin.
Paclitaxel drugs (PTX) block progression of mitosis by protecting microtubules against disassembly and preventing chromosomes from achieving metaphase spindle configuration [41]. Researches have observed that paclitaxel-treated cells have defects in mitotic spindle assembly, chromosome segregation and cell division [41]. Experiments in vivo by Yung et al. [42] show that anti-miR-21 reduces tumor volume in NSCLC, and the combination of paclitaxel and anti-miR-21 demonstrated greater ability to reduce cancer cell proliferation than either agent administered alone. The simulations in Fig 12(a) mimic this experiment; the effect of PTX is accounted for by reducing 1.4λ Ci to 1.3λ Ci (i = 1, 2), and the effect of anti-miR-21 is accounted for by reducing λ m 1 to λ m 1 /2. We note however that in our model the cancer is at an earlier stage and its volume is much smaller compared to the volume of 0.8 cm 3 in [42].
Gefitinib is a drug used in the treatment of NSCLC. It blocks the production of EGF-EGFR and thus obstructs the MAPK and AKT pathways [43]. Tumor cells that are initially sensitive to gefitinib may eventually lose sensitivity due to the emergence of acquired resistance [44]. Alternative mechanisms are currently being explored aimed to overcome the development of gefitinib resistance in the patients of NSCLC [44,45]. Recent studies [8,25] show that miR-21  modulates gefitinib sensitivity. In particular, Shen et al. [8] demonstrated in vivo that reduction in miR-21 significantly restored gefitinib sensitivity by up-regulation of PTEN expression and the inactivation of AKT and MAPK pathways. We can use our model to represent the experimental results of Shen et al. [8]. We account for the effect of gefitinib by reducing 3λ E , in the case of EGFR mutation, to 1.5λ E , and the effect of anti-miR-21 by reducing λ m 1 to λ m 1 /2. Fig 12  (b) shows that anti-miR-21 alone reduces the growth of tumor volume, but in combination with gefitinib the reduction is significantly larger. This is in qualitative agreement with the results in [8], although here again our model considers an early stage of a tumor whereas, in [8], the tumor volume is already 0.5 cm 3 .
Cisplatin induces cancer cell apoptosis by inhibition of DNA synthesis and repair in cell cycle [46]. The efficacy of cisplatin is initially high, but in the majority of cancer patients it eventually drops due to cisplatin resistance. Many mechanisms of cisplatin resistance have been described, including changes in cellular uptake, drug efflux, increased detoxification of the drug, inhibition of apoptosis, and increased DNA repair [46]. Experiments by Zang et al. [30] show that down-regulation of miR-155 can enhance the sensitivity of lung cancer cells to cisplatin treatment through the induction of DNA damage and apoptosis via the restoration of the mitochondrial apoptotic pathway. The simulation in Fig 12(c) shows that anti-miR-155 alone reduces the tumor volume, but in combination with cisplatin the reduction is significantly higher.

Discussion and conclusion
Worldwide, lung cancer is the leading cause of cancer deaths, and approximately 85% of lung cancer cases are NSCLC [1,2]. Five years survival rate for NSCLC is significantly higher for those diagnosed at early stage [3]. Unlike mammography for breast cancer or colonoscopy for colon cancer there are no reliable tools for early detection of lung cancer; most lung cancers are first diagnosed on symptoms. Hence, there is increased focus on identifying biomarkers for detection of NSCLC at early stage [5].  = 1, 2); anti-miR-21 is accounted for by reducing λ m 1 to λ m 1 /2. (b) Both anti-miR-21 and the gefitinib, a drug that inhibits the EGFR, reduce the growth of tumor; anti-miR-21 is accounted for by reducing λ m 1 to λ m 1 /2, and EGFR-inhibitor gefitinib is accounted for by decreasing 3λ E to 1.5λ E . (c) Anti-miR-155 is accounted for by reducing λ m 2 to λ m 2 /1.2 and cisplatin, a drug that induces cancer cell apoptosis, is accounted for by increasing d C and d D by factor 1.1. A growing body of evidence suggests that exosomal miRs may be used as serum biomarkers for prognosis of malignant tumors [5,7]. In NSCLC the highest overexpressed miRs are miR-21, miR-205 and miR-155 [6,10]. Since exosomal miRs concentration in the blood are positively correlated to their concentrations in tissue [10][11][12], it is important to understand how the concentrations of miRs 21, 205 and 155 in NSCLC tissue are related to the progression of the cancer, both in terms of tumor growth and tumor-front invasion.
In the present paper we developed a mathematical model that relates the role of the above exosomal miRs in tissue to cancer cells proliferation and invasion. MiRs 21 and 205 regulate cell proliferation through MAPK and PI3K-AKT pathways, while miR-155 regulates apoptosis through the Apaf-1-Caspase 9 complex. The mathematical model includes, separately, invasion and proliferation phases of NSCLC. In invasion, the 'directed migration coefficient' χ plays a critical role: the tumor front increases as χ increases. In order to visualize this monotonic relation, we used a 'relatively' large χ; in vivo this parameter may be much smaller.
In Simulations of the model of tumor proliferation establish quantitative relations between the total mass of over-expressed miRs (21,205,155) and tumor volume. Because of the positive correlation between miRs in cancer tissue and serum [10][11][12], the present model may be viewed as a first step toward establishing a combination of miRs 21, 205, 155 and possibly additional miRs as serum biomarkers for early detection of NSCLC. As more experimental and clinical data become available, the model could then be refined by estimating more precisely some of the parameters, by expanding the genetic network of Fig 2, and by precisely relating the concentrations of miRs in serum to miRs in lung tissue.

Parameter values
In the sequal we shall use the following conversion of units: 1Da = 1g/mol, so that 1 mol of a protein with molecular weight mkDa has a mass of m Â 10 3 g: Also, 1 Molar = 1 mol/L = 10 −3 mol/cm 3 . Hence 1nM = 10 −12 mol/cm 3 , and Concentration of 1nM of a protein with molecular weight mkDa ¼ m Â 10 À 9 g=cm 3 : ð22Þ

Steady state concentrations
In some estimating parameters, we use steady state equations; we denote the steady state concentrations of species X by X 0 . EGF-EGFR Cancer cells express 2 − 3 × 10 6 EGFR proteins per cell [47]. We take the average to be 2.5 × 10 6 EGFR per cell, or 2:5Â10 6 N A fraction of a mole, per cell, where N A = 6.022 × 10 23 is Avogadro's number. The molecular weight of EGFR is 170kDa [47]. By Eq (21), the mass of EGFR in one cell is 2:5Â10 6 N A Â 170 Â 10 3 g ¼ 7:0573 Â 10 À 13 g. Assuming that one cell has a volume of 10 −9 cm 3 , we find that the concentration of EGFR is 7.0573 × 10 −4 g/cm 3 . We assume that the concentration of EGF-EGFR is not limited by the availability of EGF, hence E 0 = 7.0573 × 10 −4 g/cm 3 .
Exosome In breast cancer, 10 6 cancer cells release 5 × 10 8 exosomes in 24 hours [58]. Assuming that the number of cancer cells in lung cancer is 4 × 10 8 per cm 3 , and taking the average diameter of exosomes to be 70nm, we estimate the mass density of E C by E C0 = 3.6 × 10 −10 g/cm 3 .
MiR-21, miR-205 and miR-155 For simplicity, we assume steady state of Eqs (10) and (11), so that m i We also assume that the cellular concentration of miR-21 is proportional to the exosomal concentration of miR-21, in the sense that m i 10 ¼ gm 10 , where m 10 is the steady state of m 1 . Similarly, we assume m i 20 ¼ gm 20 , where m 20 is the steady state of m 2 . In the simulations, we take γ = 10, but the simulations do not change qualitatively if we use different values of γ of the same order of magnitude.

Parameter estimation
In the sequal, in expressions of enhancement of the form X KþX or inhibition of the form 1 1þX=K , the parameter K, the 'half-saturation' of X, will be taken to be the steady state of X. Thus K MA = K M = M 0 , K A = A 0 , K Ap = A p0 and K C = C 0 , and K ME = M 0 , K TM = K TA = T 0 , K PA = P 0 . However, we make an exception in the case of miRs; we assume that the inhibition of protein expressions by miRs is more significant than inhibition by signaling proteins, and take K mT = m 10 /5, K mP = m 10 /5 and K m 2 = m 20 /5. For a species with concentration X and half-life t 1/2 , the dynamics of its degradation or death is given by 2λ C2 . We also assume that in steady state and take λ C = 0.46 day −1 [77]. Hence λ C1 = (4/3)λ C = 0.6133 day −1 and λ C2 = (2/3)λ C = 0.3067 day −1 .
From the steady state equation of Eq (12) with K M = M 0 , K A = A 0 , K Ap = Ap 0 and C 0 = C M / 2, we get We assume that apoptosis rate through intrinsic apoptosis pathway (Apaf-1/Caspace-9) is higher than apoptosis through the extrinsic signaling pathway [78], and take d D to be larger than d C , so that d D = 0.9λ C and 2d C = 0.1λ C . Hence d D = 0.414 and d C = 0.023. Diffusion coefficients: We take D C = 8.64 × 10 −7 cm 2 day −1 [77]. Diffusion of a sphere is inversely proportional to its diameter. We assume that the average diameter of exosome and cells are 70nm and 10μm respectively. Then we get D E C ¼ 8:6Â10 À 7 7Â10 À 3 ¼ 1:23 Â 10 À 4 cm 2 day À 1 . By [79], we have the relation D m 1 where D V and M V are the diffusion coefficient and molecular weight of vascular endothelial growth factor (VEGF), respectively, and D V = 8.64 × 10 −2 [77], and M V = 24kDa [80]. The molecular weight of miR-21 is M m 1 = 7 kDa, hence D m 1 = 0.13028 cm 2 day −1 , and similarly D m 2 = 0.13028 cm 2 day −1 . The directed migration coefficient χ is taken to be in the range 3 × 10 −4 − 3 × 10 −2 cm 5 g −1 day −1 [33]. In cell invasion, χ should be much larger than in cell growth, so we take χ = 2 × 10 −2 cm 5 g −1 day −1 for the model of tumor invasion.
Parameter estimation for Eq (16): Cancer cells can survive in hostile environment better than normal cells, so the apoptosis rate d N should be somewhat larger than d C ; we take d N = 1.1d C = 0.0253 day −1 and λ N = 0.8λ C = 0.368 day −1 . Since cancer cells replication is less susceptible to damage, d DN should be larger than d D ; we take d DN = 1.1d D = 0.4554 day −1 . We choose ε = 0.1, as in [40]. We assume that C 0 in the proliferation phase to be somewhat larger than the average density 0.4 g/cm 3 in the invasion phase; we take C 0 = 0.46 g/cm 3 and N 0 = 0.14 g/cm 3 .

Sensitivity analysis
We performed sensitivity analysis on some of the production parameters of the system Eqs (1)-(9), (15)-(18); we also included the important parameter d Ec which was only fitted. Following the method of [81] we performed Latin hypercube sampling and generated 1000 samples to calculate the partial rank correlation coefficients (PRCC) and the p-values with respect to the tumor radius at day 60. We have taken the range of each parameter from 1/2 to twice its value in Table 2. The results are shown in Fig 13. We see that the production rates that increase proliferation through EGF-EGFR ! MAPK and EGF-EGFR ! AKT pathways, namely, (λ E , λ M ) and (λ E , λ M , λ MA ) are positively correlated to tumor radius. On the other hand the production rates of cell-replication inhibitors, λ T and λ P , and the production rate of apoptosis-promotor apoptosome, λ Ap , are negatively correlated. Since miR-21 blocks the inhibitors T and P, so if λ m 1 E C grows the tumor volume will increase. Hence λ m 1 is positively correlated and d Ec is negatively correlated.

Computational method
We employ moving mesh method to numerically solve the free boundary problem for the tumor proliferation model. To illustrate this method, we take Eq (16) as example and rewrite it as the following form: @Nðr; tÞ @t ¼ D N DNðr; tÞ À divðuNÞ þ F; ð23Þ where F represents the term in the right hand side of Eq (16). Let r k i and N k i denote numerical approximations of i-th grid point and Nðr k i ; ntÞ, respectively, where τ is the size of time-step. The discretization of Eq (23) is derived by the fully implicit finite difference scheme: where N r ¼ i . The mesh moves by where u kþ1 i is solved by the velocity equation.