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

n class="Disease">Lung cancer is the leading cne">ause of an class="Disease">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 n class="Chemical">lipid bilne">ayer, and contne">aining fuctionne">al proteins, mRNan class="Gene">As 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].

In the present pn class="Gene">aper we develop a mathematical model thne">at relne">ates the role of the exosomal an class="Gene">miRs in lung cancer tissue to cancer cells proliferation and invasion. Since there is a positive correlation between exosomal miRs in serum and tissue in lung cancer [10-12], the model may serve as a first step toward establishing miRs as reliable serum biomarkers for early detection.

n class="Gene">A simple schematic of a cell proliferation in the context of an class="Disease">cancer is shown in Fig 1. When epidermal growth factor (EGF) ligands to its receptor EGFR, it initiates activation of the Ras-Raf-MEK-ERK pathway [13, 14] and the PI3K-AKT pathway [15-17]. Both pathways lead to cell proliferation [18-21] through activation of mTOR [16, 20, 22]. EGF-EGFR is negatively regulated by ERK [14, 23].

Fig 1

Two pathways, Ras-Raf-MEK-ERK and PI3K-AKT lead to the cell proliferation.

Sharp arrows indicate activation/enhancement and blocked arrows indicate inhibition.

n class="Gene">EGFR tyrosine kinne">ase inhibitor (TKI) modulne">ates the activne">ation of both Ran class="Gene">AS and PI3K, thereby inhibiting the activation of the downstream pathways of EGF-EGFR [13, 24, 25]. PTEN modulates the activation of AKT through converting PIP3, generated by PI3K, to PIP2 by dephosphorylation [15]. When DNA damage occurs, a signaling pathway activates Apaf-1 and caspase 9, forming an apoptosome, which leads to apoptosis through activation of caspase 3 [26, 27].

In n class="Disease">NSCLC, the most expressed exosomal an class="Gene">miRs are miR-21, miR-155 and miR-205 [10]. In Fig 2(a), we simplified the network of Fig 1 by using MAPK and AKT to represent the Ras-Raf-MEK-ERK and PI3K-AKT pathways. Fig 2(a) also shows the effect of overexpression of miR-21 and miR-205 on NSCLC proliferation. Recent studies established that miR-21 blocks TKI [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.

Fig 2

Abbreviated version of Fig 1 depicting the roles of miR-21, miR-205 and miR-155.

MAPK represents the Ras-Raf-MEK-ERK signaling pathway and AKT represents the PI3K-AKT signaling pathway. (a) MiR-21 blocks TKI; miR-21 and miR-205 block PTEN. (b) MiR-155 blocks Apaf-1+caspase 9. Sharp arrows indicate activation/enhancement and blocked arrows indicate inhibition.

Mn class="Gene">APK represents the Ras-Raf-MEK-ERK signaling pathway and AKT represents the PI3K-AKT signaling pathway. (a) MiR-21 blocks TKI; miR-21 and miR-205 block PTEN. (b) MiR-155 blocks Apaf-1+caspase 9. Sharp arrows indicate activation/enhancement and blocked arrows indicate inhibition.

In this pn class="Gene">aper we consider growth and invne">asion of an class="Disease">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 aspn class="Gene">ects of an class="Disease">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, Mn class="Gene">APK, to represent the Ras-an class="Gene">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/cm3.

Table 1

List of variables in unit of g/cm3.

Notation	Description
E	EGF-EGFR concentration
M	active MAPK concentration
P3	PI3K concentration
A	active AKT concentration
T	TKI concentration
P	PTEN concentration
Ap	Apaf-1-caspase 9 apoptosome concnetration
EC	cancer-shed exosome concentration
m1	exosomal miR-21 concentration (inlcuding miR-205)
m1i	cellular miR-21 concentration (inlcuding miR-205)
m2	exosomal miR-155 concentration
m2i	cellular miR-155 concentration
C	cancer cells density
N	normal healthy cells density

n class="Chemical">Equne">ations for proteins an class="Gene">As in [33], the dynamics of the proteins within cancer cells will appear in the form where C0 is the steady state density of cancer cells.

n class="Chemical">Equne">ation for an class="Gene">EGF-EGFR ( The equation for EGF-EGFR is given by The coefficient λ is the production rate of the EGF-EGFR complex and the factor 1/(1 + M/K) is the inhibition by ERK [14]; d is the degradation rate of E.

n class="Chemical">Equne">ation for Man class="Gene">APK ( The MAPK pathway is activated by the EGF-EGFR [13, 14], a process resisted by TKI [8, 13]. Hence where d is the degradation rate of M.

n class="Chemical">Equne">ation for an class="Gene">AKT ( The activation of the AKT pathway initiates with the activation of PI3K (P3) 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 P3 takes the form where S is the concentration of Ras.

n class="Gene">AKT is activne">ated by PI3K which is negne">atively regulne">ated by an class="Gene">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 P3 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: where λ = λ λ/d, K = K/μ.

n class="Chemical">Equne">ation for TKI ( The production of TKI is inhibited by an class="Gene">miR-21 [8, 25]. Recalling that miR-21 is only a fraction of m1, we write the equation for TKI in the form

n class="Chemical">Equne">ation for an class="Gene">PTEN ( The expression of PTEN is inhibited by both miR-21 and miR-205 [28, 29]. Hence

n class="Chemical">Equne">ation for an class="Gene">Apaf-1-caspase 9 apoptosome ( The expression of Apaf-1 is down-regulated by miR-155 [30]. Hence

n class="Chemical">Equne">ation for exosome ( an class="Disease">Cancer cells shed exosomes at a rate λC. We assume that exosomes are degraded, releasing their miRs, when merging with cancer cells. Taking the rate of this degradation to be , the equation for the concentration of exosomes is given by where the term DΔE represents dispersion (or diffusion) of exosomes.

n class="Chemical">Equne">ations for exosomal an class="Gene">miR-21 ( 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 λE ⋅ C/(K + C), and obtain the equation Similarly, the equation for miR-155 is given by

n class="Chemical">Equne">ations for an class="Gene">miR-21 () and miR-155 () in cancer cells. Since and lie in cancer cells, they diffuse with same coefficient as cancer cells. Hence, the equations for and are given by and respectively.

We will n class="Gene">apply the mathematical model to consider two phases of tumor progression of lung cancer, invasion and proliferation.

Model for tumor invasion

n class="Chemical">Equne">ation for an class="Disease">cancer cell ( 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 χ∇ ⋅ (C∇C) 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 dC, damage to cancer cells (at rate proportional to dC) 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).

Bounn class="Gene">dary and initial conditions for an class="Disease">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

n class="Chemical">Equne">ations for an class="Disease">cancer cells ( 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-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 −χ∇ ⋅ (C∇C) in Eq (12) is then neglected and replaced by the term , where 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 n class="Disease">tumor is spherical and denote its moving bounan class="Gene">dary (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, , where is the unit radial vector.

n class="Chemical">Equation for : We assume thne">at the combined densities of healthy and an class="Disease">cancer cells in the tumor is constant (θ), and take We also assume that D = D. Adding Eqs (16) and (15), we obtain so that We assume that the free boundary r = R(t) moves with the velocity of cells:

Bounn class="Gene">dary and initial conditions for an class="Disease">tumor proliferation model. We impose the boundary conditions: We assume that the concentrations of proteins TKI, PTEN and Apaf-1 which inhibit tumor growth are ‘relatively’ high, i.e., above the steady state, so that they initially decrease as the tumor begins to increase. We also assume that the remaining proteins which are involved in promoting tumor growth, are initially below their steady state. One choice of initial conditions is given below:

Results

Results for tumor invasion

In sin class="Gene">mulating the invne">asion 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.

Table 2

Summary of parameter values.

Notation	Description	Value used	References
DEc	diffusion coefficient of exosomes	1.23 × 10−4 cm2day−1	[82] & estimated
Dm1	diffusion coefficient of miR-21	0.13028 cm2day−1	[61, 83] & estimated
Dm2	diffusion coefficient of miR-155	0.13028 cm2day−1	[61, 83] & estimated
DC	diffusion coefficient of cancer cells	8.64 × 10−7 cm2day−1	[77]
DN	diffusion coefficient of normal healthy cells	8.64 × 10−7 cm2day−1	[77]
χ	directed migration coefficient	3 × 10−4 − 3 × 10−2 cm5g−1day−1	[33] & estimated
λE	production rate of EGF-EGFR	1.1741 × 10−3 day−1 g/cm3	[47] & estimated
λM	production rate of MAPK	1.6499 × 10−2 day−1	[47, 48] & estimated
λA	production rate of AKT	2.9422 × 10−3 day−1	[67, 68] & estimated
λMA	activation rate of AKT by MAPK (Ras)	1/2	[75] & estimated
λT	production rate of TKI	8.4013 × 10−4 day−1 g/cm3	estimated
λP	production rate of PTEN	2.3352 × 10−4 day−1 g/cm3	[47, 70] & estimated
λAp	production rate of Apaf-1	4.4095 × 10−3 day−1 g/cm3	[73] & estimated
λEc	production rate of exosome by NSCLC cells	9.81 × 10−9 day−1	[58] & estimated
λm1	production rate of miR-21 & miR-205 by EC	0.8626 × 10−3 day−1	[59, 74] & estimated
λm2	production rate of miR-155 by EC	0.4313 × 10−3 day−1	[59, 74] & estimated
λC1	growth rate of NSCLC cells due to MAPK	0.6133 day−1	[77] & estimated
λC2	growth rate of NSCLC cells due to AKT	0.3067 day−1	[77] & estimated
λN	growth rate of normal healthy cells	0.092 day−1	[77] & estimated
ε	competition for space coefficient	0.1	[40]
dE	degradation rate of EGF-EGFR	0.8318 day−1	[47] & estimated
dM	degradation rate of MAPK	0.6931 day−1	[63–66]
dA	degradation rate of AKT	0.6931 day−1	[67, 68]
dT	degradation rate of TKI	0.3466 day−1	[69] & estimated
dP	degradation rate of PTEN	22.1807 day−1	[70, 71]
dAp	degradation rate of Apaf-1	2.7726 day−1	[73]
dEc	degradation rate of exosome	21.8 day−1	fitted
dm1	degradation rate of miR-21 and miR-205	0.5545 day−1	[74]
dm2	degradation rate of miR-155	0.5545 day−1	[74]
dC	natural death rate of cancer cells	0.023 day−1	[77] & estimated
dD	death rate of cancer cells due to DNA damage	0.414 day−1	[77] & estimated
dN	natural death rate of normal healthy cells	0.023 day−1	[77] & estimated

Table 3

Summary of parameter values (continued).

Notation	Description	Value used	References
E0	S.S. 1concentration of EGF-EGFR	7.0573 × 10−4 g/cm3	[47]
M0	S.S. concentration of MAPK	8.4 × 10−6 g/cm3	[48, 50]
P30	S.S. concentration of PI3K	1.56 × 10−6 g/cm3	[53, 54]
A0	S.S. concentration of AKT	9.362 × 10−7 g/cm3	[53]
T0	S.S. concentration of TKI	8.5050 × 10−5 g/cm3	estimated
P0	S.S. concentration of PTEN	1.88 × 10−7 g/cm3	[47]
Ap0	S.S. concentration of Apaf-1-caspase 9	2.84 × 10−5 g/cm3	[57]
EC0	S.S. concentration of C-shed exosome	3.6 × 10−10 g/cm3	[58]
m10	S.S. concentration of miR-21 and miR-205	2.8 × 10−13 g/cm3	[59, 60]
m20	S.S. concentration of miR-155	1.4 × 10−13 g/cm3	[10, 59]
C0	S.S. density of cancer cell	0.4 g/cm3	[62]
N0	S.S. density of normal healthy cells	0.14 g/cm3	[62] & estimated
CM	carrying capacity of NSCLC cells	0.8 g/cm3	[62]
KME	inhibition of EGF-EGFR by ERK	3.936 × 10−5 g/cm3	[53, 54] & estimated
KTM	inhibition of MAPK by TKI	8.5050 × 10−5 g/cm3	estimated
KTA	inhibition of AKT by TKI	8.5050 × 10−5 g/cm3	estimated
KPA	inhibition of AKT by PTEN	1.88 × 10−7 g/cm3	[47]
KmT	inhibition of TKI by miR-21	0.56 × 10−13 g/cm3	[59, 60] & estimated
KmP	inhibition of PTEN by miR-21 and miR-205	0.56 × 10−13 g/cm3	[59, 60] & estimated
Km2	inhibition of Apaf-1 by miR-155	0.28 × 10−13 g/cm3	[10, 59] & estimated
KMA	half-saturation of MAPK (Ras) in AKT activation	8.4 × 10−6 g/cm3	[48, 50] & estimated
KC	half-saturation of cancer cell on merging with exosome	0.4 g/cm3	[62] & estimated
KM	half-saturation of MAPK on cancer cell proliferation	8.4 × 10−6 g/cm3	[48, 50] & estimated
KA	half-saturation of AKT on cancer cell proliferation	9.362 × 10−7 g/cm3	[47] & estimated
KAp	half-saturation of Apaf-1 on NSCLC cell apoptosis	2.84 × 10−5 g/cm3	[57] & estimated

1S.S. refers to steady state.

We explore how spn class="Gene">ecific mutne">ations affan class="Gene">ect 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 cm5g−1day−1. In the case of mutations in EGFR, MAPK or AKT, χ is unchanged but the production rates λ, λ and λ 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.

Fig 3

Invasion of cancer cells with density C(x, t).

The parameters are as in Tables 2 and 3, and χ = 3 × 10−2 cm5g−1day−1. (a) The first row: The control case; a mutation in EGFR where λ is increased by 1.3-fold; a mutation in MAPK where λ is increased by 1.6-fold; a mutation in AKT where λ is increased by 1.8-fold. (b) The second row: Using anti-miR-21, where λ is reduced by a factor 2 compared to the first row. (c) The third row: Using both anti-miR-21 and anti-miR-155, where λ and λ are reduced by a factor 2 compared to the first row. The time is in unit of day, and x is in unit of cm.

The parameters are as in Tables 2 and 3, and χ = 3 × 10−2 cm5g−1n class="Gene">day−1. (a) The first row: The control case; a mutne">ation in an class="Gene">EGFR where λ is increased by 1.3-fold; a mutation in MAPK where λ is increased by 1.6-fold; a mutation in AKT where λ is increased by 1.8-fold. (b) The second row: Using anti-miR-21, where λ is reduced by a factor 2 compared to the first row. (c) The third row: Using both anti-miR-21 and anti-miR-155, where λ and λ are reduced by a factor 2 compared to the first row. The time is in unit of day, and x is in unit of cm.

The sn class="Gene">econd 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 Fig 3.

Fig 4 sin class="Gene">mulates the dynne">amics of the totne">al linear mass of the an class="Disease">cancer cells in the control case and the cases of the three mutations when (i) no drug is applied, (ii) anti-miR-21 is applied, and (iii) both anti-miR-21 and anti-miR-155 are applied. We see that, by day 60: (i) each mutation increased the total cancer linear mass by approximately 25% compared to the control case; (ii) anti-miR-21 reduced the total tumor linear mass by approximately 40%; and (iii) in combination with anti-miR-155 the reduction was 65%, and the total linear tumor mass did not grow faster than in the untreated control case.

Fig 4

Total linear mass of cancer cells.

The mutations and anti-miRs are the same as in Fig 3, and all parameters are the same as in Fig 3. In the legend, anti-m1 indicates anti-miR-21, and anti-m1&m2 indicates the combination of anti-miR-21 and anti-miR-155. The time is in unit of day, and the mass is in unit of g.

The mutations and anti-n class="Gene">miRs are the sne">ame as in Fig 3, and all parameters are the sne">ame as in Fig 3. In the legend, anti-m1 indicates anti-an class="Gene">miR-21, and anti-m1&m2 indicates the combination of anti-miR-21 and anti-miR-155. The time is in unit of day, and the mass is in unit of g.

Similar results can be obtained in the case of n class="Gene">multiple mutne">ations. Fig 5 illustrates the case of two mutne">ations (the first row), with reduction in invne">asion front by an class="Gene">approximately 40% at day 60 when the cancer is treated with both anti-miR-21 and anti-miR-155 drugs (the second row).

Fig 5

Invasion of cancer cells with density C(x, t).

The parameters are as in Tables 2 and 3, and χ = 3 × 10−2 cm5g−1day−1. (a) Mutations in both EGFR and MAPK, where λ and λ are increased by 1.15-fold and 1.3-fold, respectively; mutations in both EGFR and AKT, where λ and λ are increased by 1.15-fold and 1.4-fold, respectively; mutations in both MAPK and AKT, where λ and λ are increased by 1.3-fold and 1.4-fold, respectively. (b) Using both anti-miR-21 and anti-miR-155, where both λ and λ are reduced by a factor 2 compared to the panels in (a). The time is in unit of days, and x is in unit of cm.

The parameters are as in Tables 2 and 3, and χ = 3 × 10−2 cm5g−1n class="Gene">day−1. (a) Mutne">ations in both an class="Gene">EGFR and MAPK, where λ and λ are increased by 1.15-fold and 1.3-fold, respectively; mutations in both EGFR and AKT, where λ and λ are increased by 1.15-fold and 1.4-fold, respectively; mutations in both MAPK and AKT, where λ and λ are increased by 1.3-fold and 1.4-fold, respectively. (b) Using both anti-miR-21 and anti-miR-155, where both λ and λ are reduced by a factor 2 compared to the panels in (a). The time is in unit of days, and x is in unit of cm.

n class="Disease">Cancer invne">asion depends on the diran class="Gene">ected migration coefficient χ. In [33] the range of the parameter χ was taken to be (3 × 10−4, 3 × 10−2) cm5g−1day−1. In the simulations of Figs 3–5, we took the largest value χ = 3 × 10−2 cm5g−1day−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 cm5g−1day−1: R increases by a factor 8 and M increases by a factor 11 approximately.

Fig 6

The growth of R and M.

(a) The distance that the invasion front traveled by day 60 (R, in unit of cm). (b) The total linear mass of cancer cells at day 60 (M, in unit of g). χ ranges from 3 × 10−4 to 3 × 10−2 cm5g−1day−1. All other parameter values are the same as in the EGFR-mutation case of Fig 3.

(a) The distance that the invasion front traveled by n class="Gene">day 60 (R, in unit of cm). (b) The totne">al linear mass of an class="Disease">cancer cells at day 60 (M, in unit of g). χ ranges from 3 × 10−4 to 3 × 10−2 cm5g−1day−1. All other parameter values are the same as in the EGFR-mutation case of Fig 3.

We next n class="Gene">apply anti-an class="Gene">miR-21 and anti-miR-155 drugs to the tumor (by reducing λ and λ to λ/2 and λ/2, as in Figs 3 and 4) and denote the corresponding R and M by and . We represent the efficacy of the anti-miRs drugs by and , 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 ϕ is approximately 33% when χ = 3 × 10−4g−1day−1, and it increases to 40% when χ = 3 × 10−4g−1day−1. The efficacy ϕ is 57% when χ = 3 × 10−4g−1day−1, and it increases to 68% when χ = 3 × 10−4g−1day−1.

Fig 7

The growth of ϕ and ϕ.

(a) The efficacy (ϕ) of anti-miR drugs in reducing the distance traveled by the tumor front at day 60. (b) The efficacy (ϕ) of anti-miR drugs in reducing the total linear mass of cancer cells at day 60. χ varies from 3 × 10−4 to 3 × 10−2 cm5g−1day−1. The anti-miR drugs reduce both λ and λ by a factor 2. All other parameter values are the same as Fig 6.

(a) The efficacy (ϕ) of anti-n class="Gene">miR drugs in reducing the distne">ance traveled by the an class="Disease">tumor front at day 60. (b) The efficacy (ϕ) of anti-miR drugs in reducing the total linear mass of cancer cells at day 60. χ varies from 3 × 10−4 to 3 × 10−2 cm5g−1day−1. The anti-miR drugs reduce both λ and λ by a factor 2. All other parameter values are the same as Fig 6.

From Figs 6 and 7 we conclude that as χ increases, the an class="Disease">tumor invne">asion and totne">al mass increase, while at the sne">ame time the efficacy of anti-an class="Gene">miRs drug also increases. The same results (not shown here) hold for other mutations as well as for the control case.

Fig 8 shows the relationship between the invasion distance in the control case to the total mass of n class="Gene">miR-21 and the totne">al mass of an class="Gene">miR-155 after the first 60 days. Since the concentration of miRs in serum are positively correlated to their concentration in lung cancer tissue [10-12], Fig 8 suggests that miR-21 and miR-155 could potentially be used as serum biomarkers for NSCLC, in line with suggestions made in [5, 7].

Fig 8

Growth of invasion distance as a function of the total mass of m1 and total mass of m2.

(a) The invasion distance of cancer cells in the control case as a function of the total mass of miR-21 and miR-205. (b) The invasion distance of cancer cells as a function of the total mass of miR-155. All the parameter values are the same as the control case in Fig 3.

(a) The invasion distance of n class="Disease">cancer cells in the control case as a function of the totne">al mass of an class="Gene">miR-21 and miR-205. (b) The invasion distance of cancer cells as a function of the total mass of miR-155. All the parameter values are the same as the control case in Fig 3.

Results for tumor proliferation

The sin class="Gene">mulations of proliferation of an class="Disease">cancer cells are based on the model Eqs (1)–(8), (16)–(18) with boundary Conditions (19) and initial Conditions (20). We increase both λ and λ 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-39]. We also increase the steady state C0 from 0.4 g/cm3 to 0.46 g/cm3 to reflect the fact that invading cancer cells have sparser density than proliferating cells. We take the steady state density of healthy cells, N0, to be 0.14 g/cm3 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 n class="Gene">days. Most of the concentrations are either an class="Chemical">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.

Fig 9

Average densities/concentrations of all the variables in the model.

λ and λ are increased by a factor 1.4. θ = 0.6 g/cm3, C0 = 0.46 g/cm3 and N0 = 0.14 g/cm3. All other parameter values are the same as in Tables 2 and 3.

λ and λ are increased by a factor 1.4. θ = 0.6 g/cm3, n class="Chemical">C0 = 0.46 g/cm3 and N0 = 0.14 g/cm3. an class="Gene">All other parameter values are the same as in Tables 2 and 3.

We note that in estimating some of the parameters of the model n class="Chemical">equations we assumed steady-stne">ate of the vne">arious vne">ariables (cells, proteins, an class="Gene">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/cm3, and the average density of normal healthy cells stabilize at 0.1337 g/cm3, while C + N remains approximately equal to 0.6 g/cm3 at the entire time.

Fig 10 shows the growth of the n class="Disease">tumor radius and volume, and of the totne">al mass of m1 and m2 for the first 60 days. From these profiles we can deduce relations between the total volume of the tumor and the total mass of m1 and of m2, at day 60. These relations are shown in Fig 11.

Fig 10

The growth of tumor radius R(t), tumor volume, total mass of m1 and total mass of m2 for the first 60 days.

λ and λ are increased by a factor 1.4. θ is taken to be 0.6 g/cm3. All other parameter values are the same as in Tables 2 and 3.

Fig 11

Volume of tumor as a function of the total mass of miR-21 and that of miR-155.

(a) The volume of tumor as a function of the total mass of miR-21 after the first 60 days. (b) The volume of tumor as a function of the total mass of miR-155 after the first 60 days. All the parameter values are the same as Fig 10.

λ and λ are increased by a factor 1.4. θ is taken to be 0.6 g/cm3. an class="Gene">All other parameter values are the sne">ame as in Tables 2 and 3.

(a) The volume of n class="Disease">tumor as a function of the totne">al mass of an class="Gene">miR-21 after the first 60 days. (b) The volume of tumor as a function of the total mass of miR-155 after the first 60 days. All the parameter values are the same as Fig 10.

Fig 11 may suggest that n class="Gene">miR-21 and miR-155 could be used as biomarkers for determining the volume of NSCLC when this volume is still extremely small; however this suggestion is highly speculative at this time, since other miRs shed from both the cancer cells and immune cells are also circulating in the blood.

We note that the corresponding Fig 8 sian class="Gene">mulates a different situation, where the tumor is already well established in {x < 0} and its front begins to invade into {x > 0}.

Treatment

It is well known that n class="Disease">cancer cells in an class="Disease">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.

n class="Chemical">Paclitaxel drugs (an class="Chemical">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λ to 1.3λ (i = 1, 2), and the effect of anti-miR-21 is accounted for by reducing λ to λ/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 cm3 in [42].

Fig 12

Tumor growth under different drugs. EGFR mutation is accounted for by increasing λ by factor 3.

(a) Paclitaxel inhibits the division of cells, where λ are reduced from 1.4λ to 1.3λ (i = 1, 2); anti-miR-21 is accounted for by reducing λ to λ/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 λ to λ/2, and EGFR-inhibitor gefitinib is accounted for by decreasing 3λ to 1.5λ. (c) Anti-miR-155 is accounted for by reducing λ to λ/1.2 and cisplatin, a drug that induces cancer cell apoptosis, is accounted for by increasing d and d by factor 1.1.

(a) n class="Chemical">Paclitaxel inhibits the division of cells, where λ are reduced from 1.4λ to 1.3λ (i = 1, 2); anti-miR-21 is accounted for by reducing λ to λ/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 λ to λ/2, and EGFR-inhibitor gefitinib is accounted for by decreasing 3λ to 1.5λ. (c) Anti-miR-155 is accounted for by reducing λ to λ/1.2 and cisplatin, a drug that induces cancer cell apoptosis, is accounted for by increasing d and d by factor 1.1.

n class="Chemical">Gefitinib is a drug used in the trene">atment of an class="Disease">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λ, in the case of EGFR mutation, to 1.5λ, and the effect of anti-miR-21 by reducing λ to λ/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 cm3.

n class="Chemical">Cisplatin induces an class="Disease">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, n class="Disease">lung cancer is the leading cne">ause of an class="Disease">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].

n class="Gene">A growing body of evidence suggests thne">at exosomal an class="Gene">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-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 pn class="Gene">aper we developed a mathematical model thne">at relne">ates the role of the above exosomal an class="Gene">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 the model of n class="Disease">tumor proliferation, whne">at makes the an class="Disease">tumor volume grow is the fact that the combined densities of cancer cells (C) and healthy cells (N) is limited, i.e, C + N = constant. Simulations of the proliferating model show qualitative agreement with experimental results in the treatment of NSCLC. In these experiments, miR-21 and miR-155 are used to re-establish sensitivity of cancer cells to specific chemotherapeutic drugs.

Sin class="Gene">mulations of the model of an class="Disease">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-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 sn class="Chemical">equne">al we shall use the following conversion of units: 1an class="Gene">Da = 1g/mol, so that Also, 1 Molar = 1 mol/L = 10−3 mol/cm3. Hence 1nM = 10−12 mol/cm3, and

Steady state concentrations

In some estimating parameters, we use steady state an class="Chemical">equne">ations; we denote the steady stne">ate concentrations of span class="Gene">ecies X by X0.

n class="Gene">EGF-an class="Gene">EGFR Cancer cells express 2 − 3 × 106 EGFR proteins per cell [47]. We take the average to be 2.5 × 106 EGFR per cell, or fraction of a mole, per cell, where N = 6.022 × 1023 is Avogadro’s number. The molecular weight of EGFR is 170kDa [47]. By Eq (21), the mass of EGFR in one cell is . Assuming that one cell has a volume of 10−9 cm3, we find that the concentration of EGFR is 7.0573 × 10−4 g/cm3. We assume that the concentration of EGF-EGFR is not limited by the availability of EGF, hence E0 = 7.0573 × 10−4 g/cm3.

Mn class="Gene">APK The molne">ar concentration of Ras is 0.4μM [48] and its molan class="Gene">ecular weight is 21kDa [49]. Hence, by Eq (22), the concentration of Ras is M01 = (0.4 × 103) × 21 × 10−9 g/cm3 = 8.4 × 10−6 g/cm3. The molar concentration of Raf is 0.013μM [48, 50], and its molecular weight is 72kDa [48]. Hence the concentration of Raf is M02 = (0.013 × 103) × 72 × 10−9 g/cm3 = 0.936 × 10−6 g/cm3. The molar concentration of MEK is 1.4μM [48, 49] and its molecular weight is 43kDa [51]. Hence the concentration of MEK is M03 = (1.4 × 103) × 43 × 10−9 g/cm3 = 6.02 × 10−5 g/cm3. The molar concentration of ERK is 0.96μM [48, 49] and its molecular weight is 41kDa [52]. Hence the concentration of MEK is M04 = (0.96 × 103) × 41 × 10−9 g/cm3 = 3.936 × 10−5 g/cm3. All the M0 are of the same order of magnitude. In estimating parameters, we shall use the steady state concentration M0 of the MAPK. We take it to be that of Ras, that is M0 = M01, since steady state of Ras leads to steady state in the MAPK pathway.

PI3K The molar concentration of PI3K is 8nM [53, 54], and its moln class="Gene">ecular weight is 195kDa [55]. Hence, by Eq (22), the concentration of PI3K is P30 = 8 × 195 × 10−9 g/cm3 = 1.56 × 10−6 g/cm3.

n class="Gene">AKT The molar concentration of AKT is 15.1nM [47], and its molecular weight is 62kDa [47]. Hence, by Eq (22), the concentration of AKT is given by A0 = 15.1 × 62 × 10−9 g/cm3 = 9.362 × 10−7 g/cm3.

n class="Gene">PTEN The molar concentration of PTEN is 4nM [56], and its molecular weight is 47kDa [56]. Hence, by Eq (22), the concentration of PTEN is P0 = 4 × 47 × 10−9 g/cm3 = 1.88 × 10−7 g/cm3.

TKI TKI inhibits the activation of PI3K, and n class="Gene">PTEN inhibits the activne">ation of an class="Gene">AKT. We assume that these two inhibitions are proportional, that is T0/E0 = P0/P30. Hence T0 = P0E0/P30 = 8.5050 × 10−5 g/cm3.

n class="Gene">Apaf-1 The molne">ar concentration of an class="Gene">Apaf-1 ranges from 0.1μM to 0.5μM [57]. We take it to be 0.2μM. The molecular weight of Apaf-1 is 142kDa [57]. Hence, by Eq (22), the concentration of Apaf-1 is (0.2 × 103) × 142 × 10−9 g/cm3 = 2.84 × 10−5 g/cm3. We assume that the concentration of Apaf-1-caspase-9 apoptosome is not limited by the availability of caspase 9, hence the concentration is A = 2.84 × 10−5 g/cm3.

Exosome In n class="Disease">breast cancer, 106 an class="Disease">cancer cells release 5 × 108 exosomes in 24 hours [58]. Assuming that the number of cancer cells in lung cancer is 4 × 108 per cm3, and taking the average diameter of exosomes to be 70nm, we estimate the mass density of E by E = 3.6 × 10−10 g/cm3.

n class="Gene">MiR-21, an class="Gene">miR-205 and miR-155 For simplicity, we assume steady state of Eqs (10) and (11), so that . We also assume that the cellular concentration of miR-21 is proportional to the exosomal concentration of miR-21, in the sense that , where m10 is the steady state of m1. Similarly, we assume , where m20 is the steady state of m2. 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.

The range of molar concentration of n class="Gene">miR-21 in healthy individuals is 0.1-0.326 amol/μL [59, 60], and we tne">ake it to be 0.2 amol/μL, thne">at is . The molan class="Gene">ecular weight of miR-21 is 7kDa [61]. Hence, by Eq (22), the cellular concentration of miR-21 is (0.2 × 10−3) × 7 × 10−9 g/cm3 = 1.4 × 10−12 g/cm3. The concentrations of miR-205 and miR-155 are approximately the same as concentration of miR-21 [10]. Since in our model we combine miR-21 and miR-205, we take , while . Hence, m20 = 1.4 × 10−13g/cm3 and m10 = 2.8 × 10−13 g/cm3.

an class="Disease">NSCLC We take, for the invne">asion model, C = 0.8 g/cm3 and an class="Chemical">C0 = 0.4 g/cm3 [62].

Parameter estimation

In the sn class="Chemical">equne">al, in expressions of enhne">ancement of the form or inhibition of the form , the parameter K, the ‘hne">alf-sne">aturation’ of X, will be tne">aken to be the steady stne">ate of X. Thus K = K = M0, K = an class="Gene">A0, K = A and K = C0, and K = M0, K = K = T0, K = P0. 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 = m10/5, K = m10/5 and K = m20/5.

For a spn class="Gene">ecies with concentration X and half-life t1/2, the dynne">amics of its degradation or death is given by

Parameter estimation for n class="Chemical">Eq (1): The half-life of an class="Gene">EGFR ranges from 8 to 24 hours [47]. We take it to be 20 hours, i.e. t1/2 = 5/6 days, so that d = ln2/t1/2 = 0.8318 days−1. From the steady state of Eq (1) with K = M0, we deduce that λ = 2dE0 = 1.1741 × 10−3 day−1 ⋅ g/cm3.

Parameter estimation for n class="Chemical">Eq (2): The hne">alf-life of an class="Gene">KRAS is 12 hours [63]; the half-life of Raf is 30 hours [64]; the half-life of MEK is 8 hours [65], and the half-life of ERK is 24 hours [66]. We take the half-life of the combined pathway MAPK to be 24 hours, that is t1/2 = 1 day, then d = 0.6931 day−1. From the the steady state of Eq (2) with K = T0, we find that λ = 2dM0/E0 = 1.6499 × 10−2 day−1.

Parameter estimation for n class="Chemical">Eq (3): The hne">alf-life of an class="Gene">AKT ranges from 12 to 36 hours [67, 68]. We take it to be 24 hours, that is, t1/2 = 1 day, so that d = 0.6931 day−1. We assume that the activation of AKT pathway by Ras is weaker than the activation by EGF-EGFR, and take λ = 1/2. From the steady state of Eq (3) with K = T0, K = P0 and K = M0, we find that λ = (16/5) ⋅ dA0/E0 = 2.9422 × 10−3 day−1.

Parameter estimation for n class="Chemical">Eq (4): The hne">alf-life of the TKI drugs an class="Chemical">erlotinib, ASD9291 and sunitinib are 36, 50 and 40-60 hours [69]. We take the half-life of TKI to be 48 hours, that is, t1/2 = 2 days. Hence d = 0.3466 day−1. From the steady state equation of Eq (4) with K = m10/5, we find that λ = 28.5dT0 = 8.4013 × 10−4 day−1 ⋅ g/cm3.

Parameter estimation for n class="Chemical">Eq (5): The half-life of an class="Gene">PTEN is 45 minutes [70, 71], that is t1/2 = 0.03125 days, so that d = 22.1807 day−1. From the steady state of Eq (5) with K = m10/5, we find that λ = 56dP0 = 2.3352 × 10−4 day−1 ⋅ g/cm−3.

Parameter estimation for n class="Chemical">Eq (6): The hne">alf-life of an class="Gene">Apaf-1 is 1.81h [72]; the half-life of caspase-9 is 6.6h [73]. We assume that Apaf-1-caspase-9 apoptosome is as stable as caspase-9, and take the half-life of the apoptosome to be 6h, that is t1/2 = 0.25 days. Then d = 2.7726 day−1. From the steady state of Eq (6) with K = m20/5, we get λ = 56dAp0 = 4.4095 × 10−3 day−1 ⋅ g/cm3.

Parameter estimation for n class="Chemical">Eq (7): The rate of breakdown of exosomes upon contne">act with an class="Disease">cancer cells is unknown. We take this rate to be d = 21.8 day−1. From the steady state of Eq (7) with K = C0, we get λ = dE/(2C0) = 9.81 × 10−9 day−1.

Parameter estimation for n class="Chemical">Eq (8): The hne">alf-life of an class="Gene">miRs is greater than 24 hours [74]; we take half-life of miR-21 to be 30 hours, i.e. t1/2 = 1.25 day. Hence d = 0.5545 day−1. From the steady state of Eq (8) with K = C0, we get λ = 2dm10/E = 0.8626 × 10−3 day−1.

Parameter estimation for n class="Chemical">Eq (9): We tne">ake the half-life of an class="Gene">miR-155 to be 30 hours [74], and then d = 0.5545 day−1. From the steady state of Eq (9) with K = C0, we obtain λ = 2dm20/E = 0.4313 × 10−3 day−1.

Parameter estimation for n class="Chemical">Eq (12): The most common mutne">ations in an class="Disease">NSCLC occur in tumor suppressors TP53 and ALK, and in oncogenes PTEN, EGFR, KRAS, LKB1 and BRAF, and mutations seem to occur more frequently in MAPK (KRAS, BRAF) than in PI3K-AKT [21, 75, 76]. We accordingly assume that the proliferation rate of cancer cells through the MAPK pathway is higher than the proliferation rate through the AKT pathway, and take λ = 2λ. We also assume that in steady state and take λ = 0.46 day−1 [77]. Hence λ = (4/3)λ = 0.6133 day−1 and λ = (2/3)λ = 0.3067 day−1.

From the steady state n class="Chemical">equne">ation of an class="Chemical">Eq (12) with K = M0, K = A0, K = Ap0 and C0 = C/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 to be larger than d, so that d = 0.9λ and 2d = 0.1λ. Hence d = 0.414 and d = 0.023.

Diffusion coefficients: We take D = 8.64 × 10−7 cm2n class="Gene">day−1 [77]. Diffusion of a sphere is inversely proportionne">al to its diameter. We assume thne">at the average diameter of exosome and cells are 70nm and 10μm respan class="Gene">ectively. Then we get . By [79], we have the relation , where D and M are the diffusion coefficient and molecular weight of vascular endothelial growth factor (VEGF), respectively, and D = 8.64 × 10−2 [77], and M = 24kDa [80]. The molecular weight of miR-21 is M = 7 kDa, hence D = 0.13028 cm2day−1, and similarly D = 0.13028 cm2day−1. The directed migration coefficient χ is taken to be in the range 3 × 10−4 − 3 × 10−2 cm5g−1day−1 [33]. In cell invasion, χ should be much larger than in cell growth, so we take χ = 2 × 10−2 cm5g−1day−1 for the model of tumor invasion.

Parameter estimation for n class="Chemical">Eq (16): an class="Disease">Cancer cells can survive in hostile environment better than normal cells, so the apoptosis rate d should be somewhat larger than d; we take d = 1.1d = 0.0253 day−1 and λ = 0.8λ = 0.368 day−1. Since cancer cells replication is less susceptible to damage, d should be larger than d; we take d = 1.1d = 0.4554 day−1. We choose ε = 0.1, as in [40]. We assume that C0 in the proliferation phase to be somewhat larger than the average density 0.4 g/cm3 in the invasion phase; we take C0 = 0.46 g/cm3 and N0 = 0.14 g/cm3.

Sensitivity analysis

We performed sensitivity analysis on some of the production parameters of the system n class="Chemical">Eqs (1)–(9), (15)–(18); we also included the importne">ant parameter d which wne">as only fitted. Following the method of [81] we performed Lne">atin hypercube sne">ampling and generated 1000 sne">amples to calculne">ate the partial rank correlne">ation coefficients (PRCC) and the p-vne">alues with respan class="Gene">ect 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.

Fig 13

Statistically significant PRCC values (p-value < 0.01) for R(t) at day 60.

We see that the production rates that increase proliferation through n class="Gene">EGF-an class="Gene">EGFR → MAPK and EGF-EGFR → AKT pathways, namely, (λ, λ) and (λ, λ, λ) are positively correlated to tumor radius. On the other hand the production rates of cell-replication inhibitors, λ and λ, and the production rate of apoptosis-promotor apoptosome, λ, are negatively correlated. Since miR-21 blocks the inhibitors T and P, so if λE grows the tumor volume will increase. Hence λ is positively correlated and d is negatively correlated.

Computational method

We employ moving mesh method to numerically solve the free bounn class="Gene">dary problem for the an class="Disease">tumor proliferation model. To illustrate this method, we take Eq (16) as example and rewrite it as the following form: where F represents the term in the right hand side of Eq (16). Let and denote numerical approximations of i-th grid point and , respectively, where τ is the size of time-step. The discretization of Eq (23) is derived by the fully implicit finite difference scheme: where , , , and . The mesh moves by , where is solved by the velocity equation.