Combination of singularly perturbed vector field method and method of directly defining the inverse mapping applied to complex ODE system prostate cancer model

ABSTRACT We propose a new method to solve a system of complex ordinary differential equations (ODEs) with hidden hierarchy. Given a complex system of the ODE, the hierarchy of the system is generally hidden. Once we reveal the hierarchy of the system, the system can be reduced into subsystems called slow and fast subsystems. This division of slow and fast subsystems reduces the analysis and hence reduces the computation time, which can be expensive. In our new method, we first apply the singularly perturbed vector field method that is the global quasi-linearization method. This method exposes the hierarchy of a given complex system. Subsequently, we apply a version of the homotopy analysis method called the method of directly defining the inverse mapping. We applied our new method to the immunotherapy of advanced prostate cancer.


Introduction
Prostate cancer is typical among men over 50, and diagnosis among young people is rare. Unlike other cancers, prostate cancer cells are common. Prostate cancer can remain unchanged for a long time. Thus approximately one-third of men over the age of 50, and the vast majority of men over the age of 80, have cancer cells in the prostate. Prostate cancer grows extremely slowly; among the elderly especially, it often does not cause any symptoms. However, prostate cancer can grow rapidly and even spread to other organs in the body, especially to the bones [2]. In the early stages, prostate cancer does not cause symptoms and does not exhibit any external signs. Symptoms appear when the tumour grows and places pressure on the urethral tube. Since prostate cancer grows slowly, the onset of symptoms can occur years after the tumour has appeared. Typical symptoms are difficulty in urinating, increased urinary frequency, especially at night, pain during urination, blood in urine or sexual contact (rare). If cancer has metastasized to the bones, bone pain is experienced. Sometimes, the first symptoms felt are muscle pain in the back, hips or pelvis. The initial diagnosis is a visit to a family doctor who performs a rectal finger test, PSA test, trans-rectal ultrasound and biopsy.
When the prostate cancer is at the advanced stage, it is more aggressive and spreads further throughout the body. In this case, the treatment required is chemotherapy, androgen deprivation therapy (ADT) or androgen suppression therapy, which is a hormonal treatment that reduces the effect of androgen. Androgens are male hormones that can stimulate cancer growth. ADT can decelerate and even stop cancer growth by decreasing androgen levels. ADT therapy reduces androgen-dependent (AD) cancer cells by preventing growth and inducing apoptosis [18]. At this hormone refractory stage, the AD cells are replaced by androgen-independent (AI) cells. These cells do not require the normal levels of androgen to sustain growth and are also resistant to the apoptotic effect of a lowandrogen environment. In some cases, an alternative therapy treatment called intermittent androgen deprivation (IAD) is used [1]. The advantage of this treatment is the availability of an off-treatment period to the patients that provides some relief from treatment-related side effects and reduces treatment costs. In addition, IAD may also delay the progression to castration-resistant prostate cancer. Immunotherapy dendritic vaccines (dendritic cell (DC) vaccines)-based vaccination is a treatment for harnessing the potential of a patients own immune system to eliminate tumour cells in metastatic hormone-refractory cancer. Patients that are treated by immunotherapy dendritic vaccines develop an immune response following vaccination. Some patients experience a significant decline in their PSA level, and others experience stabilization of their previous progressing disease [19,21,22].
A mathematical model of prostate cancer during androgen suppression was first proposed by [9,10]. Later, the first mathematical model of prostate cancer under IAS was proposed by [8]. Their model was extended to a model with partial differential equations [7,24,25], and a model with competition between different classes of prostate tumour cells [23]. However, these models share a typical drawback: they cannot describe the biphasic decline in PSA during the on-treatment periods. The authors of the paper [5] evaluated the probability of death by prostate cancer under various treatment options using a cell kinetic model. In the paper [13], the authors proposed a mathematical model for the immunotherapy of prostate cancer.
The model in the paper [8] has been further extended to allow many key parameters to be functions of the cell androgen content [20]. Other groups of researchers have considered more detailed biochemical dynamics of androgen and androgen receptors, and the related effects [6,11].
All these mathematical models are described by a system of linear or nonlinear ordinary differential equations (ODEs) or partial differential equations (PDE). In general, a set of differential equations describing a complex realistic phenomenon contains a number of different time scales (i.e. different rates of change) that correspond to subprocesses. Given such systems, it is highly difficult to reveal the hidden hierarchy and the implicit multitime scale of the original systems that governs the equations. Hence, one cannot apply asymptotic methods. On the one hand, discovering the hierarchical structure of systems requires considerable complicated numerical techniques; on the other hand, this known hierarchical structure allows a number of asymptotic approaches to be applied for the analysis of their behaviours. A different method for revealing the hierarchy of a given system is called global quasi-linearization or its other version called singular perturbed vector field (SPVF) [3,4]. The primary idea of this method is to transfer the original system of the governing equations with the hidden hierarchy (hidden multiscale structure) to its explicit form as a singularly perturbed system (SPS). When one finds this transformation from a hidden hierarchy model into a model with the standard SPS, the analysis of the original system can be treated by the highly powerful machinery of the standard SPS theory for model reduction and decomposition. An asymptotic method that can be applied to the ODE and PDE is the homotopy analysis method (HAM) that is also a semi-analytical technique, and was first devised in the paper [14,16]. In the current study, we combine the SPVF method and a new version of the HAM, called the method of directly defining the inverse mapping (MDDIM). It was presented previously [17] to the mathematical model that describes the advanced prostate cancer treatment by immunotherapy and ADT.

Introduction to SPVF and MDDIM methods
In this section, we present in details of the SPVF. In addition, we present the algorithm for the SPVF and MDDIM methods, as well as the combination of these two methods.

SPVF method
Given a large and complex scientific model with nonlinear governing equations, the purpose of this research is to reduce the number of equations and to discover the hierarchy of the dynamical variables of the system, i.e. to decompose the system into subsystems with the fast and slow rates of the dynamical variables. Hence, one should seek new coordinates and a representation of the governing equations (of the original model) in the form of an SPS. Once these coordinates are found, the original system can be decomposed into slow and fast subsystems that allow for asymptotic and analytical methods to be applied. Shown below is the formalized general framework theory of the SPVF.
. . n. A general ODE system of order n has the form of˙ x = F( x).
We call a general ODE system, a fast-slow ODE (FS-ODE) if a small parameter exists such that the system exhibits the following form: .
We denote f n f and s n s as the dimensions of the slow and fast ODE subsystems, respectively.

Definition 2.3:
Given the descending order of the set of eigenvalues = {λ 1 , λ 2 , . . . , λ n }, the eigenvalue maximum gap is defined as follows: max i {|(λ i+1 /λ i )|} (we can also define the gap as min i {(λ i /λ i+1 )}). Let us denote by n s the index for which we obtain the maximum gap. We call s = {λ 1 , λ 2 , . . . , λ n s } the slow eigenvalues, and the remaining eigenvalues can again be re-indexed to have f = {λ 1 , λ 2 , . . . , λ n f } fast eigenvalues, where n f + n s = n. Definition 2.4: Let A n×m 1 , B n×m 2 denote the horizontal concatenation of the two matrices (A,B); therefore, the matrix of dimension n × m is obtained, where m = m 1 + m 2 .

Method explaining
Let˙ x = F( x) be a general system ODE of order n that describes a physical phenomenon with a hidden hierarchy of the rate change of variables. The rate of change is highly dependent on the choice of the coordinate system. Generally, two domains of the coordinate system exist: the domain where the system change is slow and the domain where the change is fast. If all the system variables in the coordinate system contain fast or slow time scales (but not both), then no noticeable difference is found between the rates of change. Thus the initial problem is not reduced. Our purpose is to obtain such systems of coordinates in which the ODE system will decompose into slow and fast subsystems, i.e. an FS-ODE system. Hence, we choose n arbitrary vectors { x 1 , . . . , x n } ∈ R n . We define the following n × n matrix T to be the images of these vectors under the vector field F: Let = {λ 1 , . . . , λ n } be the ascending eigenvalues by the absolute value of matrix T, and let V = { v 1 , . . . , v n } be the eigenvectors. Because the rate of change of each vector is determined by its eigenvalue, T will be decomposed according to its large and small absolute eigenvalues. The maximal gap of the eigenvalues can be determined by max i {|λ i+1 /λ i |}; we denote this index for which this expression is determined by n s . By this index, we can classify the eigenvalue/eigenvectors into two categories: fast and slow rates. If the index is equal to or smaller than n s , it belongs to the slow rate; otherwise, it is classified as fast rate.
where D slow , D fast are block diagonal matrices with the corresponding eigenvalues.
The matrix T contains eigenvectors that can be divided into two sets: fast eigenvectors and slow eigenvectors, corresponding to the large and small absolute eigenvalues (the larger the eigenvalue, the greater is the rate of change in the direction of matching eigenvector).
According to SPVF decomposition, if the vectors { x 1 , . . . , x n } are composed of a mixed rate of change, the matrix T can be decomposed to the sum Hence, we can decompose the matrix T to fast and slow subsystems.

Algorithm for SPVF method
The procedure above of the SPVF method depends on the choice of linear independent vectors { x 1 , . . . , x n }. The choice of the domain R n and these points is crucial for the algorithm and should be adapted to every particular model. The following steps are implemented (2) Compute the mean value of the vector field over the point from step 1: (3) Define the so-called control set as follows: cs = { x i ∈ : F(x i ) > F }. For simplicity, we reindex cs = { x 1 , . . . , x N cs }. (4) Build the ordered basis sets: (6) For each i = 1, 2, . . . , k, we compute the eigenvalues of the following matrix T i that correspond to the matching basis B i , . . , λ i n } be in the ascending order eigenvalues of T i . For each T i , we compute the maximum gap as follows: and obtain the desired hierarchy of the system in the new coordinates.

MDDIM method
In this section, we present the algorithm of MDDIM introduced by Liao [17].
Given a system of nth-order nonlinear differential equation, subject to the k linear boundary initial conditions where u(x) is an unknown function, x is an independent variable, I is an interval of x, N is the nonlinear operator, B i , x 0i , α i are linear operators and a constant, respectively, for 0 ≤ i ≤ k.
(2) Obtain the discrete data (the coordinates) of the plot graph from step (1) using the command ListPlot and Table from Mathematica.
Obtain the linear combination of functions from the set F ∞ , i.e. obtain the following (11) Obtain the linear combination of functions from the set F k , i.e. obtain the follow- . , k} such that the coefficients will be determined to satisfy the boundary conditions. (13) Obtain the linear combination of the functions from the setF ∞ , i.e. obtain the following space: Define the set of linearly independent functions that constitutes the basis functions for the expressions in the codomain of the operator N (given in (2)) as G cod Obtain the linear combination of functions from the set G cod ∞ , i.e. obtain the following space: cod is an initial guess that satisfies all linear boundary conditions given in (3). (17) Use the HAM [14] to write the deformation for the functions As proven by Liao in [15] the theorem of convergence, if the convergence-control parameter c 0 and the directly defined inverse mapping are chosen appropriately such that the series in step (9) is absolutely convergent, it is then the solution of the original system (2).

Mathematical model of immunotherapy prostate cancer
In this section, we present the mathematical model (governing equations) that describes advanced prostate cancer treatment by immunotherapy and ADT by a set of nonlinear ODEs. The primary assumptions of the model are as follows. Because T cells are activated and stimulated through interactions between proteins (antigens and cytokines) and receptors, the Michaelis-Menten kinetics are used for all immune response terms in the model [12]. Interleukin-2 (IL-2), a key cytokine with pleiotropic effects on the immune system, is included in the model to provide the clonal expansion dynamics of T helper cells. To simplify the model (the mathematical equations), we assume that the ratio of cytotoxic and helper T cells is constant; in addition, the rate of interaction between cytotoxic T cells and the tumour cells, which is based on antigen stimulation, is assumed to be the same for both AD and AI cells. The antigen-loaded DCs, denoted by D assume that the DCs undergo apoptosis at a constant rate and are not replenished by any mechanism other than further vaccinations. The vaccination contains D 1 antigen-loaded DCs, and the DCs are assumed to activate naive T cells based on the Michaelis-Menten kinetics model. In addition, we use the assumption that the ratio of cytotoxic and helper T cells is constant. The following mathematical model describes the interaction between androgen-dependent cancer cells ( where the function χ(t), that controlled by monitoring the serum PSA level, is given by where the function N(t) is the serum PSA concentration. The function χ(t) indicates that the androgen deprivation is switched on if the serum PSA concentration exceeds some level L 1 and is switched off when the serum PSA concentration drops below a certain level is the dynamical variable vector of the system. The initial conditions are as follows: The interpretation of the model (4) -(9) from the biological point of view: Equation (4) describes the rate of change of the number of androgen-dependent cancer cells that is a function of expressions of a difference and sums of proliferation and death, mutation to AI (androgen-independent), and killed by T cells, separately. Equation (5) describes the rate of change of the number of androgen-independent cancer cells that is a function of the expressions of a difference and sums of proliferation and death, mutation to AD (androgendependent), and killed by T cells, separately. Equation (6) describes the rate of change of a number of activated T cells that is a function of an expression of a difference and sums of activation by dendritic cells, natural death, and clonal expansion, separately. Equation (7) describes the rate of change of concentration of cytokines that is a function of an expression of a difference and sums of production by stimulated T cells, and clearance, separately. Equation (8) describes the rate of change of the concentration of androgen, which is a function of an expression of a difference and sums of homeostasis, and deprivation therapy, respectively. Equation (9) describes the rate of change of the number of DCs that is a function of an expression of natural death. The expressions for the androgen-dependent functions for AD cell growth and mutation that appear in Equations (4)-(5) are defined as follows: and respectively. When the value of A is zero, subsequently transfer mutation occurs from the AD cells to AI cells; for A = a 0 , no mutation occurs.

Analysis and results
In this section, we apply the combination of the SPVF and MDDIM methods to the mathematical model (4)-(9) subject to the initial conditions (10). The intention of this combination is to first apply the SPVF method and expose the hierarchy of the system, and then split the system into fast and slow subsystems followed by applying the MDDIM method, which is a semi-analytical method, and obtain the desired solutions of the system.

Application of the SPVF method
In this section, we apply the SPV field to the considered model with χ(t) = 1, i.e. the model is the ' on-treatment' portion of the therapy. In this case, the eigenvalues and the eigenvectors received from the system (4)-(9) by applying the SPVF algorithm are as follows: The eigenvalues and eigenvector above correspond to the matrix T i with gap max i = max n (|λ i n+1 (T i )|/|λ i n (T i |)) of all matrices calculated for the given system. According to the SPVF theory, the hierarchy of the system is λ This implies that the eigenvector u D that corresponds to the largest eigenvalue λ D indicates the fast direction of the system in the new coordinates. Subsequently, the other variables of the system change much more slowly compared to the direction of u D .
This result is satisfactory from the biological point of view because, given a treatment (χ(t) = 1), the cancer cells decay in time and the DCs (D) grow significantly. In the absence of treatment (χ(t) = 0), the DCs, androgen and cytokines grow slowly compared to the cancer cells that grow faster using the quasi-steady state approximation. The next step is to rewrite the mathematical model of immunotherapy prostate cancer in a new coordinate using the eigenvectors above. Hence, we write the following: ⎛ where t indicates the transpose operation. After the multiplication above, the next step is performed to write the dynamical variables of the original system as a function of the new coordinates, i.e. let˜ V = (Ñ 1 ,Ñ 2 ,T,Ĩ L ,Ã,D) subsequently N 1 = N 1 (˜ V), N 2 = N 2 (˜ V), T =

Application of the MDDIM method
The next step in our combination method is to apply the MDDIM method to the immunotherapy prostate cancer model in its new coordinates, i.e. apply the MDDIM method to the model (16). Hence, we recall again the vector definition of the dynamical variables of the system in the new coordinates:

),T(t),Ĩ L (t),Ã(t),D(t))
and define the following six nonlinear operators, as a vector form, for each dynamical variable of the system: such that (19) yields the original coupled equation (16). According to the MDDIM algorithm, to define the base functions set F ∞ that constitutes the base of the functions to the solution of the model, we solve the model numerically in the new coordinates, subsequently apply the least-squares approximation, and define the linearly independent set of base functions as follows: where α j is a free parameter to be determined later such that the approximation series of the solution converges to the real solutions of the model. Correspondingly, we define the space of functions that is their linear combination as This space of functions will approximate the vector solutions˜ V(t). We define a set of six functions (as we have six initial conditions) that will approximate the initial conditions as follows: Correspondingly, we define the space of functions that is their linear combination as follows: * We define the remaining set of F ∞ presented in (20) aŝ and defineˆ ∞ asˆ such that ( ∞ =ˆ ∞ ∪ 6 ). The next step in the algorithm is to define the set of basis functions for the expressions in the codomain of the nonlinear operators: g 1 (x), . . .}, and correspondingly the Span of this set as cod Following the HAM, we define the homotopy operators in the vector form, as follows: where L i are linear operators, q ∈ [0, 1], and β k (k = 1, 2, . . . , 6) is the well-known convergence control parameter. We assume that the physical variables of the model (16) can be expanded to a power series of the homotopy small parameter q in the vector form as for each dynamical variable of the system in the new coordinates. Additionally,˜ V 0 = (Ñ 1,0 ,Ñ 2,0 ,T 0 ,Ĩ L,0 ,Ã 0 ,D 0 ) ∈ 6 are the initial conditions that satisfy the boundary conditions of the considered model. According to the theory of the HAM when q varies from 0 to 1, the approximation solution varies from the initial approximation to the solution of the original problem (the real solution). To obtain the components of the series above, one should define the deformation equation as follows: for s = 0, in vector form where i is defined in (26) and and for s = 1, 2, 3, . . . in vector form where k = 1, 2, . . . , 6 and i s (t) =Ñ 1,s ,Ñ 2,s ,T s ,Ĩ L,s ,Ã s ,D s , where χ(·) is the unit step function and δ are the following expressions: (33) To obtain the homotopy series shown in (27), one should compute the func-tionsÑ 1,s ,Ñ 2,s ,T s ,Ĩ L,s ,Ã s ,D s iteratively beginning from the initial approximatioñ N 10 ,Ñ 20 ,T 0 ,Ĩ L,0 ,Ã 0 ,D 0 . To obtain these functions, using the HAM, one should apply the inverse linear operator to both sides of the deformation equations that, in scientific computing, can be time-consuming. The algorithm presented herein this problem and proposes a new solution that does not require the computation of the inverse linear operator, but rather directly defines the inverse mapping in the linear deformation equations [17]. Hence, we express the deformations equations in (31), as a vector form alternately as When the summation in the expression above is calculated differently for the various variables of the system, and J : cod ∞ −→ ∞ is the inverse mapping that has been defined as the same for all deformation equations. In our analysis, we define the inverse mapping J as where Aα 2 + Bα + C 0 is used to weigh the terms with larger values of α, and the constants A,B,C are parameters.

Stability
In this section, we analyse the stability of the model with and without immunotherapy.
In general, given a biological/chemical phenomenon described by a mathematical model, the most important tools for investigating the model are obtaining the equilibrium points and the stability of these equilibrium points. In this section, we first obtain the equilibrium points and subsequently investigate their stability as follows: (1) Determine the fixed point vector, V * by solving F 1 ( V * ), . . . , F 6 ( V * ) = 0.
(2) Construct the Jacobian matrix, The procedure above is applied to the model (4)- (9). We set all derivatives of the dynamical variables of the system to zero, and solve first for χ(t) = 0 followed by for According to the theory of stability for an ODE system, an equilibrium point V * is stable if all the eigenvalues of the Jacobian matrix, i.e. the Jacobian evaluated at V * , contain negative real parts. The equilibrium point is unstable if at least one of the eigenvalues contains a positive real part. The Jacobian matrix of the considered system is as follows:

Stability without immunotherapy
The model (4)-(9) without immunotherapy reduces to that presented in [8]. For χ(t) = 0, we obtained the trivial state equilibrium points E 0 = (N * 1 , N * 2 , T * , I * L , A * , D * ) = 0. In this case, the eigenvalues of the Jacobian matrix are as follows: Therefore, the equilibrium point E 0 is always a locally unstable saddle point.
The parameter values from Tables (1)-(2) are substituted into the equilibrium point (39) and subsequently into Equation (36) to yield the following eigenvalues: In our results, the fixed point vector V * is the locally stable point (as expected) because all the real parts of the eigenvalues are negative. 0.008/day AD cell death rate k 1 2 ng/mL AD cell proliferation rate dependence on androgen k 2 8 Effect of low androgen level on AD cell death rate k 3 0.5 ng/mL AD cell death rate dependence on androgen em 4 1 0.00005/day maximum mutation rate a 0 30 ng/mL Normal androgen concentration γ 0.08/day Androgen clearance and production rate ω 10/day Cytokine clearance rate μ 0.03/day T cell death rate c 0.14/day Dendritic cell death rate  4 5 · 10 −6 ng/mL/cell/day Maximum rate T cells produce IL-2 g 1 10 · 10 9 cells Cancer cell saturation level for T-cell kill rate g 2 400 · 10 6 cells DC saturation level for T-cell activation g 3 1000 ng/mL IL-2 saturation level for T-cell clonal expansion g 4 10 · 10 9 cells Cancer cell saturation level for T-cell stimulation c 1 1 · 10 −9 ng/mL/cell AD cell PSA-level correlation c 2 1 · 10 −9 ng/mL/cell AI cell PSA-level correlation Note: The biological meaning of e 1 = e 2 is that T cells are able to kill AD and AI cancer cells equivalently.
From the biological point of view, this implies that for χ(t) = 1, i.e. the PSA level increases above a threshold level that predetermined value (L 1 = 15[ng/mL] ), the ontreatment (vaccination) therapy starts and the system becomes stable. However, this is a momentary stable state, until the cancer cells (AI and AD cells) begin to grow when the vaccine begins to be removed from the patient's body, the system goes into an asymptotic unstable state, the level of PSA increases above the threshold and the vaccine is given again; subsequently, the system returns to the equilibrium state. Our aim in this study is to obtain the most optimal combination of the dynamic variables of the system with the appropriate weights to prevent the recurrence of cancer cells. Hence, we implemented the algorithm of the SPVF method. This method obtains the required optimal combination. When applying this algorithm, it obtains the fast direction of the system, i.e. the fast combination of the dynamic variables of the system and the other slow variables can be neglected. In practice, this method reduces the number of equations of the system, without losing the relational information of the original system as well as the dynamics of the original system. According to our results, the fast direction of the system in the new coordinates is the direction of the eigenvalues u D because its eigenvalues are the largest. Furthermore, the optimal combination of the dynamical variables of the system is the following: −0.000824768847795 · N 1 − 0.000416458998049 · N 2 − 0.003280160258353 · T +0.999956972186126·I L +0.004721050436208 · A+0.007221668270455 · D. As shown, the dominant variable of the new system in the coordinate system is the cytokines I L . This combination of treatment should be applied until the patient reaches a certain threshold of prostate-specific antigen (PSA), which is a biomarker of the disease. Upon reaching this threshold, the patient stops receiving the therapy until their PSA levels increase above a second threshold. Upon reaching this threshold, the patient stops receiving the therapy until his PSA levels increase above a second threshold. Subsequently, the treatment is applied again by the same combination of variables or the SPVF algorithm is reapplied with new parameters (that are specific to the patient), and a new combination of variables with different fast directions as well as a different coefficient to these variables (different weights) are obtained.
Choosing the fast direction of the system is contrary to our biological intuition because the cancer cells will grow quickly in the fast direction. However, according to the coefficients of the cancer cells that appear in the eigenvectors (i.e. the coordinates), they are small relative to the others.
To compare our analysis, we solve the system by applying three different methods (see Figures 1-6): 1: we solved the full system numerically, 2: MDDIM, 3: QSS, quasisteady state approximation. Subsequently, we transfer the model to the new coordinates according to the SPVF method and compare the results by applying the following methods (see : 1: QSS, quasi-steady state approximation, 2: MDDIM, 3: SPVF method, 4: the new method, the combination between MDDIN and SPVF methods, 5: solve the full model in its new coordinates. According to the results described by the figures presented herein, it is clear that the androgen-dependent cancer cells at the beginning of the treatment starting with the initial condition point increase up to ≈ 18 [cells], and subsequently decrease to zero ( Figure 1). Furthermore, the T cells exhibit the same mathematical behaviour (Figure 3).          Since we transfer the model using eigenvectors, we obtain a negative value ofÃ which is biologically insignificant. And hence, for sake of biological understanding, we must ignore these negative values. The androgen-independent cancer cells begin at the initial condition of the system and decreases directly to zero without increasing ( Figure 2). The concentration of cytokine increases exponentially and quickly. This result is practical from the biological point of view, because cytokines are secreted by the immune system, and immunotherapy actually stimulates the immune system as well as by other cells (Figure 4).
The graph showing the concentration of the androgens ( Figure 5) shows a decrease in androgens and is a positive aspect, because the aim is to reduce the androgen that stimulates cancerous growths. The usual aggressive treatments for the reduction in androgen are orchiectomy (also named orchidectomy, and sometimes shortened as orchi), androgen blocker receptors, and androgen blockers in the thyroid gland (the adrenal glands (also known as suprarenal glands)), blocks the production of gonadotropin-releasing hormones (GnRH or LHRH) such as LH from the hypothalamus. Most cancer cells lose their sensitivity to hormone therapy after several years. In this study, the androgens decrease to zero through immunotherapy for the graph obtained by numerical simulations and stabilize for the system solution using the MDDIM and QSS methods. The number of DCs ( Figure 6) decreases, coherent to the androgen decrease.
The solution profiles of the model presented in the new coordinates are unnecessary. The graphs do not provide additional information and biological insights beyond what that presented and analysed above. However, the presentation of these graphs presents mathematical importance, in that when the system is presented in the new coordinates, the resulting graphs are indeed a linear combination of the model solutions in the old coordinates with the corresponding coefficients (Figures 7-12) Figure 13. PSA serum level, with an injection rate of 0.078 billion cells for various values of e 1 , the maximum rate T cells kill cancer cells. A wider range of e 1 is able to suppress the growth of cancer and elongate the cycles of IAD. Figure 13 presents the PSA concentration level for the different cases of the controlled parameter e 1 . When the value of this parameter is small, the level of PSA increases exponentially; when the value of this parameter increases, the level of PSA decreases. This implies that the treatment lengths vary between patients depending on the value of e 1 . If we increase e 1 , subsequently the length of ' of-treatment' is much longer, resulting in more comfort for the patient.

Conclusion
In this research, we presented the mathematical model of evolution and treatment of prostate cancer. The treatment was performed using the immune system to attack the cancer cells. Our analysis included the effect of DC vaccines on the long-term behaviour of prostate cancer. This treatment improved the life quality of the patients and may delay the resistance toward treatment. Immunotherapy altered the body immune system to help fight this type of cancer. The mathematical model included a system of nonlinear ODEs of the first order and described the interaction between androgen-dependent cancer cells, androgen-independent cancer cells, activated T cells, concentration of cytokine, concentration of androgen and number of DCs. Because the system was complex and comprised a hidden hierarchy, we applied the SPVF method to expose this hierarchy. After implementing the SPVF algorithm, we applied different asymptotical methods because the system decomposed into fast and slow subsystems. In this study, we combined two semi-asymptotical-analytical methods called the SPVF and MDDIM methods to investigate the progression of the cancer cells. In addition, our analysis included a comparison between different methods for solving the system (described the interaction between the treatment and the cancer cells). The comparison involved the implementation of our new method (the combination of MDDIM and SPVF), MDDIM, SPVF, solving the full model numerically using the Runge-Kutta simulation and QSS approximation. In addition, we investigated the stability of the system around its equilibrium points.
In general, after applying the SPVF method, i.e. transferring the system of the governing equations to new coordinates and exposing the small parameters of the system, the system exhibited the same form as that assumed in the QSS approximation. Therefore, the question emerged as to why the SPVF method must be applied, To apply the QSS approximation, the system should be in the form of an SPS, i.e. a system with explicit separation to fast and slow subsystems, thereby implying that the system should contain explicitly small parameters. Otherwise, one should assume, based on intuition or the experimental results, which dynamical variables of the system change fast and which change slow. In the presented model, we assumed that the androgen, cytokines and DCs change slowly compared to the other variables of the system. This assumption was reasonable because the time scales at which these processes occurred were much shorter than that of the populations of growing cancer cells. In contrast to the QSS method, the SPVF method could be applied to any system without the initial conditions regarding the variables of the system, as well as without transferring the system to the SPS form. The SPVF method revealed that the fast direction of the system was a combination of the 'old' variables of the system such that the cancer cells changed (grew) slowly and decreased to zero. After we applied the SPVF method that only transferred the model to new coordinates, we applied the MDDIM method that is a semianalytical method to obtain the solution profiles of the system. Our results, in comparison with others, indicated that the optimal interaction of immunotherapy was the interaction described by the fast direction of the system with the appropriate coefficients. Nevertheless, the treatment should be adapted to the patient and other parameters should be considered for different patients. For example, the parameter e 1 could be a patient-specific parameter, as it is the maximum rate that T cells kill cancer cells per day. These parameters were controlled by the therapist; for patients with weaker immune systems, more frequent injections (appropriate combination between these parameters) could be much more effective. Similarly, the amount of DCs could be controlled by changing Equation (9), for example, to the form of (dD/dt) = δ − cD where δ is the injection rate.

Disclosure statement
No potential conflict of interest was reported by the authors.