Stochastic Analysis of Nonlinear Cancer Disease Model through Virotherapy and Computational Methods

: Cancer is a common term for many diseases that can affect anybody. A worldwide leading cause of death is cancer, according to the World Health Organization (WHO) report. In 2020, ten million people died from cancer. This model identiﬁes the interaction of cancer cells, viral therapy, and immune response. In this model, the cell population has four parts, namely uninfected cells (x), infected cells (y), virus-free cells (v), and immune cells (z). This study presents the analysis of the stochastic cancer virotherapy model in the cell population dynamics. The model results have restored the properties of the biological problem, such as dynamical consistency, positivity, and boundedness, which are the considerable requirements of the models in these ﬁelds. The existing computational methods, such as the Euler Maruyama, Stochastic Euler, and Stochastic Runge Kutta, fail to restore the abovementioned properties. The proposed stochastic nonstandard ﬁnite difference method is efﬁcient, cost-effective, and accommodates all the desired feasible properties. The existing standard stochastic methods converge conditionally or diverge in the long run. The solution by the nonstandard ﬁnite difference method is stable and convergent over all time steps.


Introduction
Cancer is a family of diseases associated with developing abnormal cells to seize or transmit to other parts of one's body. Cancer is the rapid emergence of abnormal cells which arise outside their normal limits, and it may occupy the linked parts of the body and transfer to the tissues afterward. Virotherapy is the treatment of cancer that detects or destroys cancer cells during the process, and healthy cells are not harmed. Tuwairqi et al. presented the qualitative analysis of cancer cells in the cell population [1]. In 2011, Crivelli et al. presented the cell dynamics with the recommendation of the control strategies [2]. In 2020, Nouni et al. analyzed the tumor cells' dynamics for the immune response of a virological model for cancer therapy [3]. Storey et al. developed a deterministic model for treating a tumor via oncolytic treatment [4]. Abernathy et al. investigated the dynamics of the cell population, including interactions between infected and uninfected brain tumor cells [5]. Matos et al. studied the new approaches for treating cancer-like infections [6]. Makaryan et al. analyzed the awareness strategies within immune cell actions and biological therapy techniques [7]. Malinzi et al. presented the wave propagation model for dynamics of chemo/virotherapy cancer [8]. Bajzer et al. developed a co-infection dynamic of the tumor with measles in the human body [9]. Timalsina et al. developed computational techniques to model tumor virotherapy in the cell population [10]. Rommeifanger et al. developed a melanoma tumor model in the cell dynamics and performed its qualitative analysis [11]. Wares et al. established a mathematical model for cell-cycle-specific cancer virotherapy [2]. Liu et al. launched a comparison analysis of a deterministic and stochastic model for tumor-immune responses to chemotherapy [12]. Eftimie et al. investigated the complex dynamics of cell populations with well-known epidemiology techniques [13]. Santiago et al. presented the optimal control interventions [14]. Kim et al. established the hybrid analysis of cancer in the cell population [15]. Berg et al. introduced multidimensional modeling of oncolytic tumor virotherapy [16]. Some notable models related to cervical cancer and many more are presented in References [17][18][19][20][21][22][23]. The well-known methods in the sense of stochastic are presented in References [24,25]. This work aimed to understand the complex interplay among tumor cells, oncolytic viruses, and immune response. Thus, a dynamical analysis was employed to investigate the optimal therapeutic strategies for cancer remission. Stochastic analysis of the cancer disease is more realistic, practical, accurate, and close to nature. The stochastic differential equations have no analytic solutions, due to the non-differentiable term of Brownian motion. Thus, there is a need for computational methods to solve the said problems; furthermore, our focus is on those methods that restore the model's dynamical properties. That is why we moved to construct the nonstandard finite difference method in the sense of stochastic. The rest of the paper is organized based on the following sections: In Section 2, the deterministic cancer model's formulation has fundamental properties. Sections 3 and 4 deal with the stochastic model's transition probabilities, positivity, boundedness and implementation methods, convergence, and comparative analysis. Finally, the conclusion is presented in Section 5.

Deterministic Formulation
For any time, the states of the model are described as follows: x(t) represents the uninfected cancer cells, y(t) denotes the infected cancer cells, v(t) gives the virus-free cells, and z(t) shows the immune cells. Furthermore, the incoming and outgoing ratios are defined as λ, which is the growth rate of cells; C, which is the carrying capacity; d, which represents the death rate of cells that die due to natural causes or with infection; β, which is the rate of the oncolytic virus on cancer cells; δ, which is the rate of infected cells; b, which is the rate at which new version particles are released with burst size; and γ, which is the rate of decay of free virus-cell. Other ratios are defined as follows: α is the rate of uninfected cells due to immunity of the body, µ is the rate of infected cells due to weak immunity, k is the rate of elimination at which cancer cells become immune, h 2 is the rate of stimulation of uninfected cells by immunity, h 1 is the response rate at which infected cells become immuned, and ρ is the rate of decay of immune cells. The systematic flow of cancer disease is presented in Figure 1. The nonlinear ordinary differential equations by using the law of mass action are as follows:

Analysis of Model
In this section, we define the feasible region of the system (1)-(4) as follows:   . Otherwise, the system is unbounded.
Proof. Consider the cell population function as follows: By using the Gronwall's inequality, we obtain the following: as required.
Consider the x = 0, y = 0, v = 0, and z = 0 in the given system, and then we obtain the disease-free equilibrium (DFE-D 1 ) = r 1 d 1 , 0, 0, 0 . Then, by entering the x = 0, y = 0, v = 0, and z = 0 into the system and solving them simultaneously, the endemic equilibrium is obtained as follows:

Reproduction Number
To determine the reproduction number, let us apply the next-generation matrix method to find the transmission and transition matrices of the system (1)-(4) after substituting the disease-free equilibrium as follows: Hence, the dominant eigenvalue of the AB −1 is called the reproduction number and denoted as R 0 = abr 1

Local Stability
In this section, we present two well-known theorems to determine the local stability of the model. Again, consider Equations (1)-(4) as a function of F, G, H, and K as follows: The Jacobian matrix for the system (5)-(8) becomes the following: Theorem 1. If the reproduction number R 0 < 1, then the disease-free equilibrium (DFE), , 0, 0, 0 , is locally and asymptotically stable. Otherwise, the system is unstable.
Proof. The Jacobian matrix (9) at the disease-free equilibrium, , 0, 0, 0 , is evaluated as follows: All coefficients of the polynomial are positive. Therefore, the Routh-Hurwitz stability criterion for a 2nd-degree polynomial is satisfied. Therefore, the disease-free equilibrium is stable.
, is locally and asymptotically stable.
Proof. The Jacobian matrix (9) at the endemic equilibrium, E 1 = (x 1 , y 1 , v 1 , z 1 ), is as follows: By considering the coefficients of the characteristic equation above as a 4th-order polynomial, we have the following: It satisfies the Routh-Hurwitz stability criterion for a 4th-degree polynomial. Therefore, the endemic equilibrium is also stable.

Stochastic Cancer Virotherapy Model
Let us consider the vector M = [x, y, v, z] T for the system (1)-(4). Transition probabilities are presented in Table 1 for the expectations and variance. Table 1. Transition probabilities of cancer epidemic model.

Euler Maruyama Method
In this section, we utilize the Euler Maruyama scheme to determine the numerical solution of differential Equation (11), and the scientific parameters of the model are presented in Table 2 as follows: where ∆t is the time-step size, and ∆B n is the discretization parameter independent paths. By using the values of the parameters presented in Table 2, we plotted the graphs of the Euler Maruyama scheme for disease-free equilibrium (DFE) and endemic equilibrium (EE) (see Figures 2 and 3) with MATLAB software.

Data Curation
In 2020, the total population of both sexes was approximately 220,892,332, which includes the male sex (113,672,007) and the female sex (107,220,324) [20]. The number of new cases of both sexes and all ages is shown in Figure 4. Using the least-square-curve method technique, the desired fitting is presented in Figure 5. The estimated values from the data are presented in Table 2, and the value of the reproduction number is R 0 = 1.4167. Furthermore, Figure 6 depicts the residual of the actual data. Hence, the desired values of the transmission rates are helpful to study the system (1)-(4) graphically.

Non-Parametric Perturbation of Model
This section introduces the stochasticity in each compartment of the system (1)-(4). Then Equations (1)-(4) become as follows: where, for σ i, i = 1, 2, 3, 4 is the randomness of the model; and B(t) is the Brownian motion.

Positivity and Boundedness of Stochastic Model
Consider U(t) = (x(t), y(t), v(t), z(t)), and the norm: In addition, denote C 2,1 1 R 4 × (0, ∞) : R + as the families of all positive function V(U, t) defined on R 4 × (0, ∞) respectively. Let the function be twice differentiable in U and once in t, and then we have the following: If "L" acts on a function U * C 2,1 R 4 × (0, ∞) : R 4 + , then where t means transpose.
where we set inf ϕ = ∞(ϕ is an empty set).
Since τ m is increasing as m → ∞ , we have the following: Then τ ∞ ≤ τ e . Now we wish to show that τ ∞ = ∞, as desired.

Computational Methods
This section deals with well-known methods, such as the stochastic Euler, the stochastic Runge Kutta, and the proposed stochastic nonstandard finite difference method with the given non-negative initial conditions as follows:

Stochastic Euler
The stochastic Euler method could be defined on the system (13)- (16). (See Appendix A).

Stochastic Runge Kutta
The stochastic Runge Kutta method could be developed on the system (13)- (16). (See Appendix B).

Stochastic NSFD
The stochastic nonstandard finite difference could be developed for the system (13)- (16).

Stability Analysis
This section determines the model's stability by considering the equilibrium state and the theorem. The model is linearized for the stochastic nonstandard finite difference. (See Appendix D).

Comparison Section
The stochastic nonstandard finite difference method was compared with other stochastic numerical methods. It is easy to see that other stochastic numerical methods conditionally converge or diverge with larger time-step values by looking at the numerical solutions. On the other hand, the stochastic nonstandard finite-difference scheme remains convergent for all time-step sizes. Figures 7-12 show these results.      Here, notice that, when we increase a time-step size, the Stochastic Runge-Kutta method fails to restore the dynamical properties of the model. However, the stochastic nonstandard finite difference method is still convergent.

Conclusions
Stochastic modeling is a reliable and efficient technique for handling highly nonlinear problems' natural phenomena. The non-parametric perturbation technique was used in establishing the stochastic cancer virotherapy model. The nonstandard finite difference method gives dynamically consistent, positive, and bounded solutions. It is believed that existing algorithms did not restore the dynamical properties of the model, such as positivity, boundedness, and dynamical consistency, when we take a considerable time-step size. However, the nonstandard finite-difference scheme solved the problem competitively and out-performed the standard methods that were compared in this work. This is visibly noticed in the unconditional convergence of the solution given by the nonstandard finitedifference scheme against the conditional convergence or divergence offered by the existing methods. It further shows that the approach by the nonstandard finite difference is novel and can be applied to solving other nonlinear stochastic problems. The following steps are helpful to reduce the risk of getting cancer by making healthy choices, such as maintaining a healthy weight, avoiding tobacco, limiting the amount of alcohol, and protecting the skin. In the future, we will extend the idea used in this work to different types of stochastic modeling; moreover, its stabilities are presented in References [26][27][28][29][30][31][32][33][34].  Similarly, the decomposition of these equations is shown in (A2)-(A4) y n+1 = y n + h[ax n v n − cy n z n − d 1 y n − by n + σ 2 y n ∆B n ], (A2) v n+1 = v n + h[by n − h 2 y n z n − d 1 v n − m 1 v n + σ 3 v n ∆B n ], (A3) z n+1 = z n + h[cy n z n + h 2 y n z n − d 1 z n + m 1 v n + σ 4 z n ∆B n ],