Parameter Sensitivity and Qualitative Analysis of Dynamics of Ovarian Tumor Growth Model with Treatment Strategy

Md. Shah Alam¹, Md. Kamrujjaman^2,3, Md. Shafiqul Islam⁴

¹Department of Mathematics, University of Louisiana at Lafayette, Lafayette, LA, USA.
²Department of Mathematics, University of Dhaka, Dhaka, Bangladesh.
³Department of Mathematics and Statistics, University of Calgary, Calgary, AB, Canada.
⁴School of Mathematical and Computational Sciences, University of Prince Edward Island, PE, Canada.

DOI: 10.4236/jamp.2020.86073                

Abstract

In this paper, we are interested to find the most sensitive parameter, local and global stability of ovarian tumor growth model. For sensitivity analysis, we use Latin Hypercube Sampling (LHS) method to generate sample points and Partial Rank Correlation Coefficient (PRCC) method, uses those sample points to find out which parameters are important for the model. Based on our findings, we suggest some treatment strategies. We investigate the sensitivity of the parameters for tumor volume, y, cell nutrient density, Q and maximum tumor size, ymax. We also use Scatter Plot method using LHS samples to show the consistency of the results obtained by using PRCC. Moreover, we discuss the qualitative analysis of ovarian tumor growth model investigating the local and global stability.

Keywords

Parameter Sensitivity, Latin Hypercube Sampling, Partial Rank Correlation Coefficient, Scatter Plot, Monotonicity, Stability Analysis

1. Introduction

Ovarian cancer is the fifth leading cause of death from non-skin cancers among women around the globe. “Silent killer” is another name of ovarian cancer, causes more deaths than any other gynecological malignancies. The American Cancer Society estimates 22,240 new cases of ovarian cancer and 14,070 deaths due to ovarian cancer only in United States in 2018. Almost 300,000 new patients have been diagnosed with ovarian cancer in 2018 (https://ocrahope.org/patients/about-ovarian-cancer/statistics/). Ovarian cancers were previously believed to begin only in the ovaries, but recent evidence suggests that many ovarian cancers may actually start in the cells in the far (distal) end of the fallopian tubes [1]. Only 20 percent ovarian cancers are detected at early age. Despite the advancement of last few decades, ovarian cancer still remains a major problem [2]. A mathematical expression called Droop’s cell quota model [3] governs the tumor growth, where cell quota represents the intracellular concentration of necessary nutrients provided through blood supply [4] [5].

For every mathematical model, input factors such as parameters are not always known with a sufficient degree of certainty because of natural variation, error in measurements or even simply a lack of current techniques to measure them. Our goal is to analyze the uncertainty of parameters of our model. Because uncertainty in parameter values chosen, introduces variability to the model’s prediction of resulting dynamics. So, the more uncertain parameters there are, the more significant the variability introduced. LHS-PRCC sensitivity analysis is an efficient tool often employed in uncertainty analysis to explore the entire parameter space of a model [6]. Scatter plot is an alternate of PRCC which works in a different way but able to give a meaningful understanding of the parameters behavior for a particular mathematical model [7].

The mathematical methods used in modeling biological systems vary according to different steps of the process. We focus on the mathematical representation of the system. However, other important steps in the modeling processes are parameters fitting and model selection [8]. Mathematical models of complex biological systems are central to systems biology [9] [10]. They can be used as an exploratory tool to complement and guide experimental work. Model simulations can be used to predict the system-wide effects of molecular targets, such as, determine the effects of molecular target(s) inhibition in specific populations [11] [12]. They can also serve as an important clinical tool, for example, classify benign and malignant tumors, predict disease prognosis for individual patients, and predict outcomes of treatments [13] [14] [15] [16]. All these studies employed to fit both on-treatment and off-treatment preclinical data using the same biologically relevant parameters. Using mathematical tools, similar study was considered in a proposed model that consisted of healthy cells, tumor cells, and mature vascular endothelial cells in the tumor [17]. The growth of the cancerous cells can also be limited by the lack of blood vessels, which carry important nutrients and supplies.

Scientists have been using ordinary and partial differential equations to model biological systems for a long time. As these models are utilized as a part of an attempt to better understanding of more and more complex phenomena, it is becoming obvious that the simple models cannot capture the complexity of dynamics observed in natural systems [18]. Different types of approaches can be taken to deal with these complexities. One way of dealing with this issue is constructing large system of ODE, which can be quite good at approximating observed behavior [19] [20]. Such models have the benefit of merging a simple, intuitive derivation with an extensive variety of possible behavior regimes for a single system. Also, these models hide lots of detailed workings of complex biological systems, where sometimes precise details are important for this system.

Ordinary differential equations (ODE) and delay differential equations (DDE) are useful in framing many biological phenomena [20]. Though delay differential equations and ordinary differential equations have many similarities, DDE have several features which make their analysis more difficult. ODE hold derivative which depend on the solution at current value of the independent variable (time), while DDE additionally contain derivatives, which are dependent on solutions at earlier times. One of the most significant difference between the ODE and DDE is the initial data. The solution of an ODE is resolved by its value at the initial point $t=a$. However, for the interval $a\le t\le b$, the term like $y\left(t-\tau \right)$ requires initial point to solve the system [20].

Basically, ovarian tumor growth model is DDE which has two phases namely on-treatment and off-treatment. Here we are investigating for only on-treatment case and for this reason, the model turns into ODE. We have introduced the Runge-Kutta method of order 5 to solve the system of non-linear differential Equation (1), as prescribed in the next Section 2.

The rest of the paper is organized as follows. The mathematical model is described in Section 2 with parameter estimations. The solution methodology is prescribed in Section 3. Also the Scatter plotting idea is articulated in this section. The results of numerical illustrations are presented in Section 4. LHS performance, Monotonicity plots analysis and PRCC studies are also investigated in this section. The contents of Section 5 are analyzing the treatment strategy to reduce the Ovarian cancer. Theoretical results such as local and global stability analysis are presented in Section 6. Finally, Section 7 concludes summary and discussion of the results.

The dynamics of mathematical model are integrated in the following section.

2. Mathematical Model of Ovarian Cancer Dynamics

Ovarian Tumor Growth Model is a simple vascularized model; a type of tumor that forms from cells that make blood vessels or lymph vessels. Vascular tumors may form on the skin, in the tissues below the skin, and/or in an organ. There are many types of vascular tumors. The most common type of vascular tumor is hemangioma, which is a benign tumor that usually occurs in infants and goes away on its own. In this model, the idea of nutrient limited induced angiogenesis has been used [21]. As already mentioned, the model is a DDE, but for on-treatment case it becomes an ODE where the delay part $\bar{y}=y\left(t_0-\tau \right)$ has established values ranging 200 - 10,300 [22]. We are investigating the case from the beginning of the on-treatment case for the time delay, $\tau =10$ and initial time $t_0=29$ which tends to give approximately $\bar{y}=y\left(t_0-\tau \right)=y\left(29-10\right)=y\left(19\right)=375$ (see Figure 7 in [22]).

Following is the ovarian tumor growth model [22],

$(\begin{array}{l}\frac{\text{d}y}{\text{d}t}={\mu }_m\left(1-\frac{q}{Q}\right)y-yd,\\ \frac{\text{d}Q}{\text{d}t}={\alpha }_1{p}_1\frac{y\left(t_0-\tau \right)}{y}-{\mu }_m\left(Q-q\right).\end{array}$ (1)

For convenience and parameter estimations, the variables and parameters of system (1) are described in Table 1.

We simulate the model for 100 days of on-treatment case to get both Tumor volume, y and Cell Quota Q, see Figure 1.

Now it’s time to describe the solution methods.

3. Methodology

3.1. LHS-PRCC

To explore the uncertainty of parameters, one of the most useful sensitivity analysis method is Latin Hypercube Sampling-Partial Rank Correlation Coefficient (LHS-PRCC). It determines the full parameter space of a model with an optimal number of computer simulations [6]. This method uses the combination of two statistical procedures, Latin Hypercube Sampling (LHS), which was first presented by McKay [23] in 1979 and Partial Rank Correlation Coefficient (PRCC) analysis.

Figure 1. Solution curves for Tumor Volume (y) and Cell Quota (Q).

Table 1. Variable and parameter list for ovarian tumor growth model.

Within a given range of parameters value, LHS samples them to generate different values at each simulation and PRCC uses those value to describe the relation of parameters with the output of a particular mathematical model [24]. PRCC is a sample based method to examine the correlation between a model output variable and parameters using sample points generated by LHS.

The goal of LHS-PRCC sensitivity analysis is to identify significant parameters which have great impact for model prediction and to rank these parameters depending on their contribution for a precise model prediction [25]. This section presents details on the steps of LHS. The method uses the following procedure:

1) Make a list of the parameters for the model with their consistent values.

2) We have to predict the uncertain parameters from the parameter lists. For some of these, it might not be difficult to find the possible range where the exact values might fall.

3) Next step is to decide the sample size and to do this we need to determine the number of simulations we intend to run. Assume we decide to run N model simulations for analysis and we have K uncertain parameters, $v_i\mathrm{,1}\le i\le K$. Then the parameter space for the uncertain parameters would be defined by K dimensions.

4) K dimensions will correspond to uncertain parameters and N determines the length of dimensions. For every uncertain parameter, each of the N input values would be selected/determined by the LHS sampling scheme.

5) We need to specify a probability density or distribution function (pdf) for each uncertain parameter to implement this LHS sampling scheme. In this way, the variability in the pdf becomes a direct measure of the variability of the uncertain parameter.

6) Each probability density function is divided into N non-overlapping equiprobable intervals for sampling the values of each parameter.

7) Each equiprobable interval of each parameter is then randomly sampled once. The parameters are uncorrelated because each parameter is sampled independently.

8) Once step 7 is complete, each of the K uncertain parameters, $v_i\mathrm{,1}\le i\le K$, will have N values. Hence, we store the sampled values in an $N\times K$ table/matrix.

3.2. Scatter Plot

Scatter plots are occasionally used to examine the correlation between a model output variable and parameters visually [26]. It’s a variance based method to find a trend or pattern of the data obtained from LHS sampling. An output variable that is sensitive to the selected parameter will yield an obvious correlated pattern in the scatter plot. Generally, a Monte Carlo algorithm is used to sample the parameter space, and multiple scatter plots are drawn illustrating the relationship between each parameter and each output variable of interest. But in this paper, we use LHS instead of Monte Carlo algorithm for sampling the parameter space. Visual recognition of the correlation between parameters and model output values can be contingent on the choice of axis scales.

For the sake of comprehension and clarity, we state and discuss our illustrated key results at this point.

4. Results

4.1. Performing LHS on the Ovarian Tumor Growth Model

There are five (since $\bar{y}$ and $\tau$ have a fixed values) uncertain or Latin Hypercube Sampling (LHS) parameters in Table 1. To identify their various roles in the model predictions, we start performing LHS. In our analysis, 100 model simulations were performed with 1000 samples per run. Thus, the parameter space (LHS matrix) for the LHS parameters has dimension of length five with each dimension specifying an uncertain parameter vector of length 100. We assume the maximum and minimum values for each of the five LHS parameters (see Table 1). The baseline value for each LHS parameter has been set to a value at or near the middle of the range for a particular parameter. For every LHS parameter, each of the 100 input values are obtained by sampling a uniform probability density distribution. The 100 input values are then used to populate the LHS matrix from which we produced monotonic plots for each variable. We calculate tumor volume, y, cell quota or cell nutrient density, Q and maximum tumor size, ymax for each run after 100 days.

4.2. Analyzing the Monotonicity Plots

Analysis of Monotonicity plot is a precondition to apply PRCC on LHS generated samples. Three outcome measures mentioned in the last section are presented in the following figures. The subplots in Figure 2 (x-axis represents the parameter values and y-axis represents the respective outcome measure) have monotonic relation with all the outcome measures except reduction in nutrient uptake rate (p₁) with maximum tumor size, ymax.

To solve this issue one approach is to breakdown that graph into two monotonic regions. If instead of the small range of outcome measures observed for p₁ in ymax, the range had been several hundred or thousand units, we would have considered truncating the range and looking at each truncated half separately. However, the current effect of p₁ in ymax is minimal since the range is very small (i.e. 1000 - 1010) for number of orders of 10³. So, no action is needed. Hence, all LHS parameters of our model have a monotonic relationship with the outcome measures tumor Volume, y, cell quota, Q and maximum tumor size, ymax.

4.3. Analyzing the PRCC

In PRCC analysis, we consider the parameters with PRCC values > 0.5 (for direct relation) or <−0.5 (for inverse relation) and corresponding small P-values (<0.05) as the most influential parameters for the model.

Figure 2. Monotonicity Plots for Q and ymax.

In each PRCC plot (Figures 3-5), x-axis contains the parameters and the y-axis contains PRCC values. Also, for each PRCC plot there is a corresponding P-value plot where x-axis represents the parameters and the y-axis represents the corresponding P-values.

We observe that each PRCC and P-value plot show strong correlation of death rate, d and maximum growth rate, ${\mu }_m$ with all three outputs.

4.4. Analyzing the Scatter Plots

In Scatter plot analysis we try to find a pattern (relation or trend) for each output corresponding to each parameter of our model.

To get the result, we use 1000 sample in each run. The subplots of Figure 6 (x-axis represents outcome measures y, Q and ymax respectively in each subplot and y-axis represents range of parameters (Table 1)) show that the most important parameters are the death rate, d and the maximum growth rate, ${\mu }_m$ as we can clearly get a trend for these two parameters. Additionally, from the first and third subplot of Figure 6, we observe that nutrient uptake coefficient ( ${\alpha }_1$ ) shows a good trend for both tumor volume, y and maximum tumor size, ymax but fails to be considered as important as d and ${\mu }_m$ due to more spreading shape at the beginning. Moreover, it’s trending nature is far away from the trend of both d and ${\mu }_m$ for the output measure cell quota, Q (second subplot of Figure 6).

Hence death rate, d and maximum growth rate, ${\mu }_m$ are the two most sensitive parameters of our model.

Next, let us proceed to test the treatment strategies for sensitive parameters.

5. Treatment Strategy

In this section, we will suggest few treatment strategies depending on the results of model simulations for different values of the two most sensitive parameters (death rate, d and maximum growth rate, ${\mu }_m$ ).

· Figure 1 shows the primary solution curves of the model with base values (see Table 1) for d and ${\mu }_m$.

· Figure 7 and Figure 8 show the change in tumor volume and cell quota respectively due to the change in d while ${\mu }_m$ is fixed at 0.995.

· Figure 9 and Figure 10 show change in tumor volume and cell quota respectively due to the change in ${\mu }_m$ while d is fixed at 0.855.

Now we want to understand the relation of the two most sensitive parameters death rate, d and maximum growth rate, ${\mu }_m$ with tumor volume, y and and cell quota, Q.

· Figure 11 shows that the tumor volume, y decreases with the increase of death rate, d and increases as the maximum growth rate, ${\mu }_m$ increases. We also observe that for the value (approximately) of $d>0.95$ and ${\mu }_m<0.98$, tumor volume, y approaches zero.

· Figure 12 shows that the cell quota Q remains zero until death rate, d approaches to 1.4 (approximately) and then increases rapidly for a small increase in d. It also shows that Q decreases rapidly for a small increase in maximum growth rate, ${\mu }_m$ until it reaches approximately 0.45. After that, Q remains zero for any increasing value of ${\mu }_m$.

The following section shows the stability of equilibrium solutions.

Figure 3. PRCC and P-values Plot for Tumor Volume (y).

Figure 4. PRCC and P-values Plot for Cell Quota (Q).

Figure 5. PRCC and P-values Plot for Maximum Tumor Size (ymax).

Figure 6. Scatter plots for Tumor Volume, Cell Quota and Maximum Tumor Size, respectively.

Figure 7. Model Simulations for ${\mu }_m=0.995$ and $d=\mathrm{0.5,0.9936,1.2}$, respectively.

Figure 8. Model Simulations for ${\mu }_m=0.995$ and $d=\mathrm{0.5,0.9936,1.2}$, respectively.

Figure 9. Model Simulations for $d=0.855$ and ${\mu }_m=\mathrm{0.5,0.85604,1.2}$, respectively.

Figure 10. Model Simulations for $d=0.855$ and ${\mu }_m=\mathrm{0.5,0.85604,1.2}$, respectively.

Figure 11. Treatment Graph of Tumor Volume (y) for d and ${\mu }_m$.

Figure 12. Treatment Graph of Cell Quota (Q) for d and ${\mu }_m$.

6. Qualitative Analysis

In this section we are going to discuss the stability of our model. Since d and ${\mu }_m$ are the two most sensitive parameters then they should play vital roles in the stability analysis of this model. Stability are of two types namely local stability and global stability.

6.1. Local Stability

Local stability of a system of ODE occurs only surrounding by a small neighborhood of equilibrium points. If we move far away from the neighborhood, local stability can be altered. Using the system (1), we consider

$f\left(y\mathrm{,}Q\right)={\mu }_m\left(1-\frac{q}{Q}\right)y-yd$ (2)

$g\left(y\mathrm{,}Q\right)=\frac{375{\alpha }_1{p}_1}{y}-{\mu }_m\left(Q-q\right)$ (3)

We have used $\bar{y}=375$ [22] in the second equation of (1). Now for Equilibrium points

$\begin{array}{l}f\left(y,Q\right)=0\\ \Rightarrow {\mu }_m\left(1-\frac{q}{Q}\right)y-yd=0\\ \therefore y=0,Q=\frac{q{\mu }_m}{{\mu }_m-d}\end{array}$ $\begin{array}{l}g\left(y,Q\right)=0\\ \Rightarrow {\alpha }_1{p}_1\frac{\bar{y}}{y}-{\mu }_m\left(Q-q\right)=0\\ \therefore y=\frac{375{\alpha }_1{p}_1\left({\mu }_m-d\right)}{qd{\mu }_m},Q=\frac{q{\mu }_m}{{\mu }_m-d}\end{array}$

In this mathematical model, y being the tumor volume, cannot be zero (also, $y=0$ makes $g\left(y,Q\right)$ undefined). So, the only biological meaningful equilibrium point of the system is

$E\left(y\mathrm{,}Q\right)\equiv E\left(\frac{375{\alpha }_1{p}_1\left({\mu }_m-d\right)}{qd{\mu }_m}\mathrm{,}\frac{q{\mu }_m}{{\mu }_m-d}\right)$

and for the choice of ${\mu }_m>d$, E is always positive.

Now the Jacobian Matrix is defined as follows:

$J=\left[\begin{array}{cc}{\mu }_m\left(1-\frac{q}{Q}\right)-d& \frac{{\mu }_{m}qy}{Q^2}\\ \frac{-375{\alpha }_1{p}_1}{y^2}& -{\mu }_m\end{array}\right]$

at, E $J=\left[\begin{array}{cc}0& \frac{375{\alpha }_1{p}_1\left({\mu }_m-d\right)}{{\left(q{\mu }_m\right)}^{2}d}\\ \frac{-{\left(q{\mu }_{m}d\right)}^2}{375{\alpha }_1{p}_1{\left({\mu }_m-d\right)}^2}& -{\mu }_m\end{array}\right]$

So, the trace and determinant of the Jacobian Matrix are

$Tr\left(J\right)=-{\mu }_m<0\left(\text{foranychoiceof}{\mu }_m\right),Det\left(J\right)=\frac{d}{{\mu }_m-d}>0,$

for the choice of ${\mu }_m>d$, which concludes that E is locally asymptotically stable as long as ${\mu }_m>d$.

Let us now define the characteristic equation of the Jacobian matrix

${\lambda }^2-Tr\left(J\right)\lambda +Det\left(J\right)=0$

$\Rightarrow {\lambda }^2+\lambda {\mu }_m+\frac{d}{{\mu }_m-d}=0$

$\Rightarrow \lambda =\frac{-{\mu }_m\pm \sqrt{{\mu }_m^2-\frac{4d}{{\mu }_m-d}}}{2}$

$\therefore \lambda =\frac{-{\mu }_m}{2}\pm \frac{\sqrt{{\mu }_m^3-{\mu }_m^{2}d-4d}}{2\sqrt{{\mu }_m-d}}$

The non-zero (see Table 1 for the value of ${\mu }_m$ ) real part of the eigenvalues $\lambda$ shows that E is hyperbolic. So, the stability of E is robust i.e. for small perturbation on E, the stability may distort but does not change the phase portrait near the equilibrium qualitatively as shown in Figure 13.

Figure 13. Phase plane with solution trajectories of ovarian tumor growth model (1).

Also, the discriminant of the characteristic equation is

$D={\mu }_m^2-\frac{4d}{{\mu }_m-d}$

Being locally asymptotically stable, E will be a stable node if $D>0$ and will be a stable spiral if $D<0$ [27].

6.2. Global Stability

If equilibrium points are stable everywhere (beyond of the small neighborhood) then they will be globally stable. We can check the global stability of our model using Lyapunov stability. Finding a Lyapunov function for our model will be cumbersome due to the number of parameters and the non-linearity of the model. So, we will start with the usual Lyapunov function for a two dimensional system of ODE

$V\left(y\mathrm{,}Q\right)=\frac{1}{2}\left(y^2+Q^2\right)$ (4)

then $\stackrel{\dot{}}{V}\left(y\mathrm{,}Q\right)=y^2{\mu }_m-\frac{y^2{\mu }_{m}q}{Q}-y^{2}d+\frac{375{\alpha }_1{p}_{1}Q}{y}-{\mu }_m{Q}^2+{\mu }_{m}qQ$ (5)

$\Rightarrow \stackrel{\dot{}}{V}\left(y,Q\right)=\underset{a}{\underset{︸}{-\left(\frac{y^2{\mu }_{m}q}{Q}+y^{2}d+{\mu }_m{Q}^2\right)}}+\underset{b}{\underset{︸}{\left(y^2{\mu }_m+\frac{375{\alpha }_1{p}_{1}Q}{y}+{\mu }_{m}qQ\right)}}$ (6)

For global stability,

$\stackrel{\dot{}}{V}\left(y,Q\right)<0$

$\Rightarrow a+b<0$

$\Rightarrow y^3\left({\mu }_{m}Q-{\mu }_{m}q-Qd\right)+Q^2\left(375{\alpha }_1{p}_1+y{\mu }_{m}q\right)-y{\mu }_m{Q}^3<0$

which is vauge in terms of biological meaning. So, we want to adopt a different approach to cope with this situation.

Clearly, a will always dominate over b i.e. $a>b$ which yields $\stackrel{\dot{}}{V}\left(y,Q\right)<0$. To verify our statement numerically, we calculate $\stackrel{\dot{}}{V}\left(y,Q\right)$ for several (starting value, base value, ending value etc.) values of all parameters (using Table 1) and variables (y and Q). Every time, we get $\stackrel{\dot{}}{V}\left(y,Q\right)<0$ and by Lyapunov stability theorem, E is globally asymptotically stable.

We can also proceed with Poincaré-Bendixon theorem to verify that E is globally asymptotically stable. Using Equations (2) and (3), we get

$\frac{\partial f}{\partial y}+\frac{\partial g}{\partial Q}=-\frac{{\mu }_{m}q}{Q}-d<0$

So, by Bendixon’s negative criterion, we can find a simply connected and positively invariant set containing no closed orbits. Note that, both y and Q will be bounded (for a quick guess see Figure 1) on that connected positively invariant set.

Then by Poincaré-Bendixon theorem, every solution starting in that connected, positively invariant set will approach to E. Hence E is globally asymptotically stable.

7. Conclusions

Both PRCC and Scatter plot method techniques are very useful to identify the parameters that have significant impacts on a mathematical model. In our investigation, we used LHS to sample points for both these methods. Since both death rate, d and maximum growth rate, ${\mu }_m$ have significant effects on tumor volume, y, cell quota Q and maximum tumor size ymax, identifying them as the most sensitive parameters can help us to introduce new treatment strategies in this field.

In this study, we observed and listed out the main findings:

1) When there is no treatment therapy for maximum growth rate ${\mu }_m=1.58$ (see Table 1), then limiting the supply of nutrient affects the growth of cancer of cells, which results in higher death rate, d.

2) On the other hand, when there is a no treatment therapy for death rate, $d=0.28$ (see Table 1), then reducing the maximum growth rate, ${\mu }_m$ will reduce the number of cancer cells.

3) Controlling nutrient supply for cancer cells with some level of treatment can have remarkable effects on cancer treatment strategies.

4) It is concluded about the treatment by referring Figure 11 that, if we can increase the death rate ( $d>0.95$ ) (by restricting nutrient supply or by using medicine) or control the maximum growth rate ( ${\mu }_m<0.85$ ) of tumor cells then the tumor volume will eventually get smaller and approach to zero (die out).

5) Both ${\mu }_m$ and d controlled the stability of the model. Based on their relation ( ${\mu }_m>d$ ), the system gets local asymptotic stability around the equilibrium point, see Figure 13.

6) Finally, we identified that the equilibrium point of our system obtains global asymptotic stability.

Acknowledgements

The authors are acknowledged to the anonymous reviewer for his constructive suggestions to improve the manuscript. The author M. Kamrujjaman research was partially supported by TWAS, Grant: 2019_19-169 RG/MATHS/AS_I.

Author’s Contributions

This work was carried out in collaboration among all authors. All authors read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Leunig, M., et al. (1992) Angiogenesis, Microvascular Architecture, Microhemodynamics, and Interstitial Fluid Pressure during Early Growth of Human Adenocarcinoma LS174T in SCID Mice. Cancer Research, 52, 6553-6560.
[2]	Siegel, R., Ward, E., Brawley, O. and Jemal, A. (2011) The Impact of Eliminating Socioeconomic and Racial Disparities on Premature Cancer Deaths. CA: A Cancer Journal for Clinicians, 61, 212-236. https://doi.org/10.3322/caac.20121
[3]	Droop, M.R. (1966) Vitamin B12 and Marine Ecology. Journal of Marine Biological Association of the United Kingdom, 46, 659-671. https://doi.org/10.1017/S0025315400033403
[4]	Kuang, Y., Nagy, J.D. and Eikenberry, S.E. (2016) Introduction to Mathematical Oncology (Chapman & Hall/CRC Mathematical Biology Series). Chapman & Hall, London.
[5]	Liu, J.G. (2019) Analysis and Computation of Some Tumor Growth Models with Nutrient: From Cell Density Models to Free Boundary Dynamics. Discrete & Continuous Dynamical Systems B, 24, 3011-3035. https://doi.org/10.3934/dcdsb.2018297
[6]	Gomero, B. (2012) Latin Hypercube Sampling and Partial Rank Correlation Coefficient Analysis Applied to an Optimal Control Problem. Master’s Thesis, University of Tennessee, Tennessee.
[7]	Friendly, M. and Denis, D. (2005) The Early Origins and Development of the Scatter Plot. The Journal of the History of the Behavioral Sciences, 41, 103-130. https://doi.org/10.1002/jhbs.20078
[8]	Bhadauria, A.S. (2015) Mathematical Epidemiology Special Attention to Cancer Growth and Its Treatment. University of Lucknow, Lucknow.
[9]	Mesiano, S., Ferrara, N. and Jaffe, R.B. (1998) Role of Vascular Endothelial Growth Factor in Ovarian Cancer: Inhibition of Ascites Formation by Immunoneutralization. The American Journal of Pathology, 153, 1249-1256. https://doi.org/10.1016/S0002-9440(10)65669-6
[10]	Hong, W.-S. and Zhang, G.-Q. (2019) Simulation Analysis for Tumor Radiotherapy Based on Three-Component Mathematical Models. Journal of Applied Clinical Medical Physics, 20, 22-26. https://doi.org/10.1002/acm2.12516
[11]	Lord, C.J. and Alan, A. (2010) Biology-Driven Cancer Drug Development: Back to the Future. BMC Biology, 8, Article No. 38. https://doi.org/10.1186/1741-7007-8-38
[12]	Hidrovo, I., et al. (2017) Experimental Method and Statistical Analysis to Fit Tumor Growth Model Using SPECT/CT Imaging: A Preclinical Study. Quantitative Imaging in Medicine and Surgery, 7, 299-309. https://doi.org/10.21037/qims.2017.06.05
[13]	Zhukov, N.V. and Tjulandin, S.A. (2008) Targeted Therapy in the Treatment of Solid Tumors: Practice Contradicts Theory. Biochemistry (Moscow), 73, 605-618. https://doi.org/10.1134/S000629790805012X
[14]	Vander, H.M.G. (2011) Targeting Cancer Metabolism: A Therapeutic Window Opens. Nature Reviews Drug Discovery, 10, 671-684. https://doi.org/10.1038/nrd3504
[15]	Motta, S. and Pappalardo, F. (2012) Mathematical Modeling of Biological Systems. Briefings in Bioinformatics, 14, 411-422. https://doi.org/10.1093/bib/bbs061
[16]	Hossine, Z., Meghla, A.A. and Kamrujjaman, M. (2019) A Short Review and the Prediction of Tumor Growth based on Numerical Analysis. Advances in Research, 19, 1-10. https://doi.org/10.9734/air/2019/v19i130115
[17]	Kuang, Y., Nagy, J.D. and Elser, J.J. (2004) Biological Stoichiometry of Tumor Dynamics: Mathematical Models and Analysis. Discrete & Continuous Dynamical Systems, 4, 221-240. https://doi.org/10.3934/dcdsb.2004.4.221
[18]	Panetta, J.C. (1997) A Mathematical Model of Breast and Ovarian Cancer Treated with Paclitaxel. Mathematical Biosciences, 146, 89-113. https://doi.org/10.1016/S0025-5564(97)00077-1
[19]	Erguler, K. and Stumpf, M.P.H. (2011) Practical Limits for Reverse Engineering of Dynamical Systems: A Statistical Analysis of Sensitivity and Parameter Inferability in Systems Biology Models. Molecular BioSystems, 7, 1593-1602. https://doi.org/10.1039/c0mb00107d
[20]	Shampine, L. and Thompson, S. (2001) Solving ODEs with MATLAB. Applied Numerical Mathematics, 37, 441-458. https://doi.org/10.1016/S0168-9274(00)00055-6
[21]	Kuang, Y., Nagy, J.D. and Elser, J. (2004) Biological Stoichiometry of Tumor Dynamics: Mathematical Models and Analysis. Discrete and Continuous Dynamical Systems Series B, 4, 221-240. https://doi.org/10.3934/dcdsb.2004.4.221
[22]	Everett, R.A., Nagy, J.D. and Kuang, Y. (2016) Dynamics of a Data Based Ovarian Cancer Growth and Treatment Model with Time Delay. Journal of Dynamics and Differential Equations, 28, 1393-1414. https://doi.org/10.1007/s10884-015-9498-y
[23]	McKay, M.D., Beckman, R.J. and Conover, W.J. (1979) A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code. Technometrics, 21, 239-245. https://doi.org/10.1080/00401706.1979.10489755
[24]	Islam, M.R. and Peace, A.L. (2018) Parameter Sensitivity Analysis and Control Strategies of a Multi-Stage Epidemic Model. Poster, Conference: Multiscale Dynamics of Infections, Ohio State University, Columbus, OH, 23-27 April 2018.
[25]	Blower, S.M. and Dowlatabadi, H. (1994) Sensitivity and Uncertainty Analysis of Complex Models of Disease Transmission: An HIV Model, as an Example. International Statistical Review, 62, 229-243. https://doi.org/10.2307/1403510
[26]	Wu, J., Dhingra, R., Gambhir, M. and Remais, J.V. (2013) Sensitivity Analysis of Infectious Disease Models: Methods, Advances and Their Application. Journal of the Royal Society Interface, 10, 2012-1018. https://doi.org/10.1098/rsif.2012.1018
[27]	Allen, L.J.S. (2007) Introduction to Mathematical Biology. Pearson/Prentice Hall, Upper Saddle River.