This section provides a detailed description of the SPVF method. The mathematical model has the dimension n = 10. The following steps are implemented

1: Select N vectors, Γ = { x → 1 , … , x → N } , x → i ∈ R n , where N ≫ n .

2: Compute the mean value of the vector filed over the point from step 1 : F ¯ = 1 N ∑ i = 1 N F → ( x → i ) ,

3: Define the following set: Γ c s = { x → i ∈ Γ : ‖ F ¯ ( x i ) ‖ > ‖ F ¯ ‖ } for simplicity let reindex Γ c s = { x → 1 , … , x → N c s } .

4: Build the ordered basis sets:

B i = { x → ( i − 1 ) n + 1 , … , x → i n } with the corresponding matrix

and let B = { B 1 , B 2 , … , B m } , A = { A 1 , A 2 , … , A m } , where m = ⌊ N c s n ⌋ .

5: Select only the reference basis set from step 4 which have | D e t ( A i ) | above the average level over all determinate basis i.e., let Ω = 1 m ∑ i = 1 m | D e t ( A i ) | , then the reference basis is B r b = { B i : | D e t ( A i ) | ≥ Ω , i = 1 , … , m } . Again let us reindex, B r b = { B 1 , B 2 , … , B k } with the matching reindex of vectors x → .

6: For each i = 1 , 2 , … , k compute the eigenvalues of following matrix T i that correspond to the matching basis B i ,

i.e., compute the determinant of the following matrix | T i − λ I | where I is the unit matrix, and solve the equation:

7: let { λ 1 i , λ 2 i , … , λ n i } be in ascending ordered eigenvalues of T i . For each T i the maximum gap is computed as:

8: Denote by i m a x the index for which g a p m a x i is maximal. Compute the eigenvectors of T i m a x , i.e, solve the system of equations:

and obtain the eigenvectors: { w → 1 i m a x , w → 2 i m a x , … , w → n i m a x } , that correspond to { λ 1 i m a x , λ 2 i m a x , … , λ n i m a x } consist of the desired coordinate system. Let n s - be the index for which ( | λ n + 1 i m a x ( T i m a x ) | / | λ n i m a x ( T i m a x | ) ) is maximal. Then the vectors { w → 1 i m a x , w → 2 i m a x , … , w → n s i m a x } and, { w → n s + 1 i m a x , w → n s + 2 i m a x , … , w → n i m a x } are the new slow and fast vectors of the slow and fast system correspondingly.

9: Rewrite the original system in the new coordinate using the eigenvectors { w → 1 i m a x , w → 2 i m a x , … , w → n i m a x } .

In this section, the mathematical model is transformed to new coordinates using the eigenvectors of the vector field.

By applying the SPVF method to the system of Equations 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, The following eigenvalues and eigenvectors were obtained:

According to the algorithm of the S P V F , the maximum gap is g a p m a x i = λ 1 λ 2 = 353.840 . The corresponding eigenvectors are as follows:

where T denote the transpose operator. This means that the original system of equations can be decomposed into fast and slow subsystems, where the fast direction of the system is in the direction of the eigenvector w → 1 corresponding to the eigenvalue λ 1 , and the slow direction of the system is in the direction of the eigenvectors w → 2 − w → 10 corresponding to the eigenvalue λ 2 − λ 10 .

The next step of the S P V F method is to transform model (1-10) using the above eigenvectors; hence, let W → be a vector of the dynamical variables of the mathematical model:

W → = ( C C , N K , W B C , C T L , A Z D n c , A Z D c , P a n c , P a c , ℱ , ℋ ) and, V → = ( x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 , x 9 , x 10 ) be the variables of the model in the new coordinates. Hence, the system can be rewritten as

where the matrix A contains the eigenvectors obtained by applying the SPVF method.

The next step is to express the old system variables as functions of the new variables. To achieve this, multiply the set of Equations 14 by the inverse matrix of A

Take the derivative of the system (14) with respect to time:

Then substitute the expressions of the R H S (right-hand side) from the system (1) a ^ (10) instead of d W → d t in (16); that is,

Where

Finally, substitute Equation 15 into Equation 17 to obtain the original mathematical model in the new coordinates with the initial conditions as follows:

The system of differential equations obtained (19) is a system of equations that only describes a mathematical model and has no biological or physical meaning because the new variables are a combination of the old variables without any expression meaning. However, the great advantage of this system is that in these new coordinates, the hierarchy of the system is precisely known; therefore, the fast and slow variables are exactly known. This procedure enables division of the new ODE system into fast and slow subsystems. The procedure for splitting into fast and slow subsystems is as follows:

Where x → = ( x 1 , … , x 10 ) is the vector of dynamical variables of the model in the new coordinates. As shown, variable x 1 is a fast variable of the new system that corresponds to the direction of the eigenvector corresponding to the largest eigenvalue λ ₁ obtained by applying the SPVF algorithm. This procedure allows us to reduce the system (19) written in new coordinates to only one significant subsystem (in this case, one equation). According to the theory of asymptotic analysis for a given system in the singular perturbed system (SPS) form, a fast subsystem can be studied while the slow system is frozen. When investigating a fast system, no important or relevant information regarding the entire system is lost. The main aim of this study was to determine the equilibrium points of the system and their stability. The advantage of the SPVF method is that the eigenvalues do not change in size under a linear transformation. Therefore, by determining the equilibrium points of the new system and analyzing their stability, an inverse transformation can be performed (using the inverse matrix of the eigenvectors) to find the equilibrium points of the original system and guarantee that these will be stable equilibrium points of the original mathematical model. The equilibrium point of the new system is determined by solving the following equation:

for d x 1 d t = 0 , i.e.,

While the other variables of the system remain constant (frozen), they can be considered values of the initial conditions in the new system. As stated before, the mathematical model is transferred from the coordinates of the dynamical variables presented by the vector to new coordinates presented by the vector using the eigenvectors w → i . The transformation process involved expressing the new dynamic variables of the system as functions of the original model’s dynamical variables, i.e., the new variables are combinations of the old variables. The results are presented in Figures 1 and 2 for different parameter values. In these figures, the black line represents the solution profile of the combination of the original variables, which cause the cancer cells to achieve stability. The red line represents the solution profile of cancer cells that achieve stability at the equilibrium point.

The parameters data that used in this research are presented in the relevant Figure. Figure 1 : This combination of variables and parameters behaves in a roughly cyclic manner, meaning that it rises and falls; however, the general trend is downward. The intervals are approximately constant; however, the values on the y-axis vary. For this combination, the cancer cells stabilize at a relatively fast rate, meaning that there is a very sharp decrease at the beginning, which then reaches an equilibrium state very quickly. The sharp decrease at the beginning of treatment is attributed to the high combination of variables and parameters at the beginning of treatment, which causes a sharp decrease in cancer cells.

Figure 3 : For these combinations of parameters and variables, it can be observed that, initially, cancer cells increase relatively sharply and then gradually fall, but not as in the previous case, where they fell cyclically. In this case, cancer cells do not decrease quickly and stabilize, but rather take relatively more time to stabilize.

Figure 2 : For this combination of parameters, the values constantly increase in the dynamic variables of the original system. This indicates the aggressiveness of the treatment, that is, the variables increase over time with the treatment. With this combination, the cancer cells stabilized very quickly, indicating that they initially decreased very quickly and then stabilized in a straight step.

In all cases, a correlation is observed between the variables and parameters to the state where the equilibrium points stabilize.

The equilibrium points of the mathematical model are examined in this section (written in the new coordinates) and their stability. In general, given a mathematical model presented by nonlinear ODE system, where the hierarchy of the variables is hidden, it is very hard and even impossible to study the stability of the equilibrium points analytically. Therefore, the great advantage of transforming the mathematical model to new coordinates is first of all exposing the hierarchy of the dynamic variables of the system, that is, of the mathematical model, subsequently, the model is split into a fast subsystem and a slow subsystem. by orders of magnitude of the eigenvalues ˆaaˆof the matrix that represents the vector field of the mathematical model. After splitting the mathematical model into subsystems, the fast subsystem can be analyzed, the equilibrium points can be found analytically in most cases as a function of the system parameters, while the other variables remain constant, and can be taken as the values of the initial conditions of the model. After finding the equilibrium points, their stability was analyzed, followed by an inverse transformation to the equilibrium points, based on principles from linear algebra, the equilibrium points and their stability are preserved under the linear transformation.

According to the results of the eigenvalues presented above the fast subsystem contain only the first variable i.e., x 1 while x 2 , … , x 10 are the slow variables.

The following steps are implements for finding the equilibrium points of the mathematical model in the new coordinates and determining their stability.

1. Substitute the slow variables as a constant into the fast subsystem (one can take the initial condition of the slow variables as the constants).

2. Setting the fast variable derivative to zero i.e., solve the fast subsystem (with slow variables as constants) and find the equilibria points of the fast variables x 1 *

(where the star notation denoted the fast equilibrium points).

3. Substitute the equilibrium points from step 2( x 1 * ) into the slow subsystem, solve the slow subsystem and find the equilibrium points of the slow variables ( x 2 * , … , x 10 * )

Here, the system consists of eight equations and eight unknown variables ( x 3 * , … , x 10 * ) .

4. Substitute the equilibrium points from steps 2 and 3 at the Jacobian matrix of the full system (the model at the new coordinates).

5. Compute the eigenvalues of the Jacobian matrix of the system (the model at the new coordinates) for each set of equilibrium point (the stable points are those with a negative real part of the eigenvalues).

6. Transform only the equilibrium points that are stable from steps 2 and 3 to the original coordinates using the inverse matrix of the eigenvectors, i.e., compute

where i indicates for different stable equilibrium points. In Figure 4 The solution profiles of the fast variable, depending on time t, are presented. As one can see from the plot graph the stability point reached after t ≈ 30 days. It is very important to note at this stage that the mathematical model in the new coordinates has no biological meaning since the new variables are a linear combination of the old variables. The only variables that have biological significance are the old variables. The important parameter from the stability analysis that can be extracted from the graph shown in Figure 4 is time. An inverse linear transformation to the equilibrium points shows that after the same time parameter, both the cancer tumor and the whole system, as described by the mathematical model, stabilize.