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

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.


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, in-tuitive 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 t b ≤ ≤ , the term like ( ) y t τ − 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 de- The dynamics of mathematical model are integrated in the following section.

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  Figure 7 in [22]). Following is the ovarian tumor growth model [22], Journal of Applied Mathematics and Physics 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.    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, ,1

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.

Performing LHS on the Ovarian Tumor Growth Model
There are five (since y and τ 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.

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  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 1 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 1 in ymax is minimal since the range is very small (i.e. 1000 -1010) for number of orders of 10 3 . 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.

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. Journal of Applied Mathematics and Physics 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, m µ with all three outputs.

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, 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 ( 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 m µ due to more spreading shape at the beginning. Moreover, it's trending nature is far away from the trend of both d and m µ for the output measure cell quota, Q (second subplot of Figure 6).
Hence death rate, d and maximum growth rate, m µ are the two most sensitive parameters of our model.
Next, let us proceed to test the treatment strategies for sensitive parameters.

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, m µ ).
• Figure 1 shows the primary solution curves of the model with base values (see Table 1) for d and m µ .
• Figure 7 and Figure 8 show the change in tumor volume and cell quota respectively due to the change in d while m µ is fixed at 0.995. Journal of Applied Mathematics and Physics • Figure 9 and Figure 10 show change in tumor volume and cell quota respectively due to the change in 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, 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, m µ increases. We also observe that for the value (approximately) of 0.95 d > and 0.98 m µ < , 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, m µ until it reaches approximately 0.45. After that, Q remains zero for any increasing value of m µ .
The following section shows the stability of equilibrium solutions.

Qualitative Analysis
In this section we are going to discuss the stability of our model. Since d and 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.

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 We have used 375 y = [22] in the second equation of (1). Now for Equili-Journal of Applied Mathematics and Physics

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 For global stability, We can also proceed with Poincaré-Bendixon theorem to verify that E is globally asymptotically stable. Using Equations (2) and (3), we get 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.

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, 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 1.58 m µ = (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, 0.28 d = (see Table 1), then reducing the maximum growth rate, 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 ( 0.95 d > ) (by restricting nutrient supply or by using medicine) or control the maximum growth rate ( 0.85 m µ < ) of tumor cells then the tumor volume will eventually get smaller and approach to zero (die out). 5) Both m µ and d controlled the stability of the model. Based on their relation ( 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.