New parameters estimation in an ordinary differential equation model for the multistep transformation of cancer cells using particle swarm optimization

2.1 Construction of the kinetic model of multistep transformation into cancer cells

Tumor growth and metastasis are influenced by the important role of tumor angiogenesis (A). Tumor angiogenesis is the proliferation of vascular endothelial cells. Then, these cells migrate following existing capillaries or veins, so that native blood vessels form a new vascular system. Tumor angiogenesis is a complex, dynamic process involving various types of cells and complex signaling networks. During its growth and metastasis, tumors are affected by the tumor microenvironment (TME), including the internal and external environment of tumor cells. The high rate of angiogenesis depends on vascular growth factors and the regulation of various angiogenesis-related signaling pathways. Therefore, the results of this study will be useful to further reveal that to highlight the importance of angiogenesis-related signaling pathways in antitumor therapeutic strategies [17].

Apoptosis or programmed cell death (D) is naturally responsible for the removal of old cells from the body. However, unregulated apoptotic signals, particularly the activation of the antiapoptotic system, allow cancer cells to escape this program. This condition will lead to uncontrolled proliferation, which results in tumor survival, therapeutic resistance, and cancer recurrence. Therefore, most anti-cancer therapies trigger the induction of apoptosis and associated cell death networks to eliminate malignant cells. This resistance is a complex phenomenon that originates from the interaction of various molecules and signaling pathways [18].

The potential for cancer development has been identified as a characteristic of genomic instability and mutations (G). In the past few decades, significant progress has been made by researchers in understanding the highly complex and varied causes of genomic instability in cancer. It has been proposed that aborted mitotic checkpoints, hypoxia, chromosome segregation errors, defective DNA damage repair, telomere shortening, reactive oxygen species, and oncogene-induced DNA replication stress cause this instability. Instability can be studied by looking at the defects associated with various components of the DNA maintenance machinery (the “babysitters” of the genome) caused by genetic mechanisms. Genetic instability is an important characteristic of cancer. Therefore, most cancers develop genetic instability at some stage in their development. Transient increases in instability are sometimes followed by a return to a relatively stable genome. However, based on dynamics, the role of instability in tumor development is not well understood [19].

Cancer development is an evolutionary process. Mutations act as a selection mechanism for cells during cancer development. Understanding the mechanistic basis of mutation accumulation is essential to elucidate the processes that shape tumor evolution. During the cell cycle, DNA replication (R) is a critical biological process. It involves careful duplication of the entire genome to ensure efficient and accurate replication, and also to limit the potential for acquisition of somatic changes. Each chromosome is divided into segments that are replicated in a defined and highly organized temporal sequence. This process is commonly called the replication timing (RT) program. In nonmalignant cells, the RT program is highly conserved in 50–70 % of the genome, while the remaining 30–50 % can vary dynamically during normal development, contributing to tissue specificity. Changes in the RT program during normal lineage differentiation are associated with differences in gene transcription levels [20].

Spencer and coworkers have presented a kinetic model of the multistep transformation of normal cells into cancer cells [13]. The equations of the kinetic model take the form of coupled ordinary differential equations (ODEs) and contain several important parameters. Each equation represents one of the four fundamental concepts of the mutation process. Some equations represent the processes of angiogenesis (A), immortality (D) (including evasion of cell death), genetic instability (G) (a function of mutation rates), and increased replication rate (R). The next step, namely, invasion and metastasis (M), is the final step that allows local tumor spread to occur.

Spencer and coworkers have reported that the number of normal cells (N) was 10⁸ [13], which is based on data from the literature [21]. These normal cells can induce angiogenesis (A), and these cells can also mutate to avoid death (D). These cells can mutate to cause genetic instability (G), and these cells can also mutate to increase the replication rate (R). In general, this kinetic model can explain the oncogenesis process. All parameters are based on breast cancer data.

Populations of cells that carry out two or three mutation processes will be given a symbol by listing the mutation processes together in alphabetical order. For example, N (normal cells) mutates into A, then mutates D to become AD (cells that have undergone two mutation processes), and then mutates R to ADR (the cell has carried out three mutation processes). A population of normal cells (N) that carried out four mutation processes was used as primary tumor cells (T). Although this kinetic model only addresses the genetic growth of primary tumor cells, it is possible to use primary tumor (T) cells because this model can invade and spread into metastatic cells (M).

According to previous research [22], human cells’ spontaneous mutation rate (k ₁) ranges from 10⁻⁷ to 10⁻⁶ mutations/gene/cell division. According to previous research [23], the loss of genes for repairing DNA can increase the genetic mutation instability factor (k ₂) rate from 10¹ to 10⁴. According to previous research [24], the estimated rate of mutation success from primary tumor cells (T) to metastatic cells (M) is in the range of 10⁻⁹–10⁻⁷ per cell division.

Without angiogenesis, which supplies blood to the tumor, the tumor cannot grow beyond 10⁶ cells [25]. The tumor size is limited to only approximately 10⁶ cells so that more than 10 % of the non-normal and non-metastatic cell population will have mutations in the A gene. Based on this fact, only a small portion of tumor cells must send angiogenesis signals to supply enough healthy blood for tumor cells. Furthermore, the load limit of a deadly tumor of approximately 10¹³ cells will restrict the number of abnormal cells that do not metastasize [26].

There are approximately 400 genes involved in primary tumor development, consisting of four categories [13]. However, it is assumed that of these 400 genes, there are only approximately 100 genes involved in four mutation categories, namely, A, D, G, and R. Based on this factual assumption, some genes can function in more than one category, thus allowing for double- and triple-state transition processes. However, it would reduce the number of genes involved from the order of 10 to 1, thus reflecting the likelihood that a single mutation would affect more than one category.

The coupled ordinary differential equations can be used to represent the kinetic model that replicates a normal cell population into a mutant cell population with the probability of the mutation process. This equation assumes that the normal cell population is constant. This is based on the assumption that normal cells leaving compartment N for another compartment are assumed to be few to avoid affecting the initial population,

The symbols of P _A , P _D , P _G , and P _R represent the number of non-metastatic cell populations with mutations in genes of categories A, D, G, and R. The symbol of P _N represents a constant normal cell population. It is assumed that the cell population leaving state N to become the population of other states is small enough, so that it does not affect the cell population of N. The number of N states is 10¹³ cells. The symbol of P N M ̄ represents the number of non-normal and non-metastatic cell populations, which is limited to a size of 10⁶ cells based on the logistic equation.

Equations (2)–(5) are the kinetic model equation format for each mutation process in a population. The probability diagram of the mutation process is used to derive the kinetic model equations (see Figure 1). For example, the process of Equation (3) can be described as follows: The number of cells performing the mutation process in compartment D increases. This finding is consistent with the growth of normal cells harboring mutations in one of the 100 genes in compartment D at a rate of k ₁ every b day. Because the cells in compartment D replicate every day without changing, the population of that compartment also grows. Deaths cause the population of compartment D to decline every d _D  day. A decrease in the population of cell exiting compartment D also. As a result, one of the other three groups (AD, DG, or DR), each with 100 genes, had a single mutation. Triple mutations are subsequently generated to create an ADGR with one gene implicated in the transition after double mutations are first obtained in one of three ways (DGR, ADG, or ADR), with ten genes involved in each transition. As a result, one of the other three groups (AD, DG, or DR), each with 100 genes, had a single mutation. Triple mutations are subsequently generated to create an ADGR with one gene implicated in the transition after double mutations are first obtained in one of three ways (DGR, ADG, or ADR), with ten genes involved in each transition. According to logistical regulations, the aggregate count of nonmetastatic and non-normal cells cannot exceed 10⁶ cells. Only when 10 % or fewer of non-normal and nonmetastatic cells contain category A mutations may angiogenesis be logistically restricted. Finally, a lethal tumor burden limit of 10¹³ cells limit the population of non-normal and nonmetastatic cells. Figure 1 describes all possible integrated mutation processes that can occur, while the parameters that may appear in the kinetic model are given in Table 1.

Figure 1:

The probability of mutation process proposed by Spencer and coworkers [13].

In this research, a simulation of the behavior of the mutation process of normal cells into cancer cells is presented in curve form. Therefore, the ordinary differential equation that describes the rate of the mutated cell population will be constructed from the ordinary differential equations from Equations (2)–(5). This research has developed 12 equations, consisting of Equations (6)–(17). The cells that enter a compartment from the preceding compartment have an impact on the population of cells in that compartment, which is the framework for creating a population equation for a mutant cell. Cells that replicate will remain in that compartment, while cells that leave that compartment will go to a new compartment. Cells that leave a mutation compartment to enter a new compartment are produced by obtaining genes from a type of mutation and becoming a new type of mutation. Cells that enter a mutation compartment are caused by acquiring one of the genes of that mutation type. Equations (1)–(17) will be solved numerically using the fourth-order Runge Kutta method, while the parameters of the kinetic model will be estimated using the PSO algorithm and programmed with Octave software.

The rate of cell division (1/b) and cell death (1/d) can occur every 1/10 day⁻¹. The rate of cell division accompanied by the R mutation process (1/b _R ) will increase, namely approximately 1/9.89 day⁻¹. The rate of cell death accompanied by D mutation (1/d _D ) decreased to approximately 1/10.08 day⁻¹ (see Appendix),

2.2 Modified particle swarm optimization (PSO) algorithm

In 1995, Kennedy and Eberhart discovered the particle swarm optimization (PSO) algorithm [31]. This PSO algorithm is easy to implement, and its convergence can be controlled through a small number of parameters [27], so PSO programming can be done using software and will achieve a fast level of convergence [32]. Due to the simple characteristics of the PSO algorithm, this algorithm has been widely used in various optimization problems. The PSO algorithm is a population-based metaheuristic stochastic algorithm that is very suitable for solving complex and multidimensional problems. However, the standard PSO method has a major weakness, namely that it cannot reach the global optimal region when facing more complicated inverse and optimal problems. This weakness occurs at the final stage of the search process, so a mismatch between local and global search will result in particles being attracted to the local optima. Often, this weakness makes convergence difficult to obtain [33]. The local best particle swarm optimization (MLBPSO) algorithm is a new modification of the PSO algorithm that aims to maximize the estimated probability of innovative parameter values. The MLBPSO algorithm shows through extensive simulations that this algorithm excels in terms of robustness to the accuracy of the estimated parameters while maintaining comparable computational costs [34].

The standard PSO method treats the g _best value of each individual as a volume-less particle in D-dimensional space. Therefore, the position and velocity of the ith particle are represented as X _i = (X _i1, X _i2, …, X _iD ) and V _i = (V _i1, V _i2, …, V _iD ), so that the moving particles can be expressed by the following equation:

In equations (17) and (18), iteration is expressed by the index t to express the position and speed that change with time. The learning parameters to determine the stochastic acceleration weight for each particle to be able to move towards the best particle individually or personally (p _best) and globally (g _best) in the evolutionary process are defined by the cognitive component C ₁ and the social component C ₂. Meanwhile, the other parameters, rand ₁() and rand ₂() are two random numbers that are generated uniformly in the range 0–1. The vector P _i = (P _i1, P _i2, …, P _iD ) is the previous best position of the i-th particle, this is usually called p _best. The vector P _g = (P _g1, P _g2, …, P _gD ) is the best particle position in the population, this is usually called g _best. The personal best position is the best position ever visited by an individual particle from the initial position to the position at time t. Meanwhile, the global best position is the best position found by each particle in the entire cluster since the first step. X _iD , V _iD , and P _iD are the D dimensions of the vectors X _i , V _i , and P _i . The parameter ω ^t is the inertial weight introduced to accelerate the PSO convergence speed based on the following formula:

where ω _max and ω _min are the maximum and minimum values of ω ^t , the values of ω _max and ω _min can be taken from the range of parameters to be optimized. Meanwhile, I and t are evolution iterations and the number of iterations. The maximum speed limit of these particles can be defined as the maximum speed, V _max. This maximum speed value can control the global position-finding ability of the particle. Meanwhile, population size is defined as a flock’s number of particles (P). Determining the size of a population (swarm) can be related to the optimization problem to be solved. Therefore, this population size is sensitive to the optimization problem, because each particle is a possible solution to the optimization problem. In some cases, the size of the initial population is an essential factor, because a larger population size can be used to meet specific needs. However, the large initial population size can result in large randomness as well.

In the standard PSO algorithm for the entire optimization process, the parameter values are set constant. In this research, three main parameters that play an important role in PSO performance based on the formulation in the description above were modified. This is to adjust the optimization process of the parameters of a model from time to time. This aims to control initial convergence, increase cluster diversity, and produce a more reliable PSO algorithm. The modified PSO algorithm must change these parameters [35],36]. In equation (17), the term ω t × V i D t has a larger inertia weight value. This will allow the PSO algorithm to take longer steps in the search domain, so this will benefit global search (exploration). This PSO algorithm can also explore finer regions of the search space when the inertia weight has a smaller value, so this benefits local search (exploitation). These two characteristics of inertial weights must find a balance in a good optimization algorithm.

In general, parameters C ₁ and C ₂ are static and their optimal values are determined empirically. Incorrect determination of the values of these initial parameters can result in divergent behavior [37]. In this research, the values of these parameters will be dynamically updated from the inertial weight on each particle by inserting several random parameters into ωi, so that the latest equation becomes as follows:

where index i represents the i-th particle. This modification of the inertial weight plays a role in accelerating the rate of convergence. This modification also ensures the balance between local and global searches. The algorithm of this modification is to prevent the particle speed from becoming stagnant. This provides different search capabilities because the inertial weights can have different values.

Two important parameters for determining PSO performance as a learning factor are cognitive constants, C ₁, and social constants, C ₂. When the C ₂ value is higher than the C ₁ value, the particle cannot escape the local optima. However, when the C ₁ value is greater than the C ₂ value, this will facilitate the particle’s ability to wander excessively in the optimization process [28]. In addition, the values of cognitive and social constants determine whether the PSO converges to local optima [38]. In this study, these two learning factors are updated dynamically using an equation that includes a random factor to prevent the particle from leaving the local optimum point and to ensure that the particle always converges to the global optimal solution:

The PSO algorithm has the biggest weakness because there are a large number of objective function evaluations. This will require longer computing time than traditional methods. However, this high objective function can be used for careful statistical analysis of the final results. Therefore, this method can determine reliable estimation results from the estimated parameters. The PSO algorithm is not very sensitive to initial guesses of a model’s parameters, making it attractive to use when there are a large number of unknown model parameters. In this research, the modified PSO algorithm is presented in detail as follows:

PSO algorithm scheme:

1. Determining the initial values of search parameters:

   Number of iterations, N _iter;

   Number of particles, N _Pt;

   Number of dimensions to be searched, N _D ;

   Vectors of length N _D with searching limit, x ^min and x ^max;

   PSO search parameters, C ₁, C ₂, ω _min, ω _max;

   set t = 0 (iteration counter).

2. Calculation of maximum particle velocity in each direction with dimension D:

   V D max = X D max − X D min / 2

3. Calculation of the position and velocity of the initial particle:

   X p , D t = X D min + r a n d X D max − X D min

   V p , D t = V D max 2 × r a n d − 1

4. Evaluate the objective function for each particle.

5. The particle positions and the particle objective function are written in the file that will be used to construct the confidence region.

6. Update X i D t , a vector with dimension N _D that contains the best position found by the whole particle swarm.

7. When the maximum number of iterations is achieved (t = N _iter), the search is terminated.

8. Calculate the inertial weight value:

   ω t = ω max − t I × ω max − ω min

   ω i = ω min + ω max − ω min × r a n d

9. Update the particle velocities for P = 1,…, N _pt; D = 1,…, N _D :

   C 1 = 1 − r a n d , C 2 = 1 + r a n d

   V i D t + 1 = ω t × V i D t + C 1 × r a n d 1 × P i D t − X i D t + C 2 × r a n d 2 × P g − X i D t

10. Update the particle positions:

    X i D t + 1 = X i D t + V i D t + 1

11. If the particle position is not inside the searching limits, the particle is placed at the violated searching limit.

12. The iteration is added to the iteration counter t = t + 1 , and the process returns to step 4.

In this research, three parameters (k ₁, k ₂, k ₃) contained in the kinetic model will be estimated using the PSO algorithm. The input three parameters of the kinetic rate constant are presented in Table 2. In the initial input, search parameters C ₁ and C₂ are made equal to 2.0, while the inertial weight ω is made equal to 1.0. The number of iterations is made equal to one thousand iterations carried out, while the number of particles used is equal to 20. The parameter search space is limited to the intervals listed in the Table 2 for the three parameters.

The PSO algorithm is an iterative optimization method for computation that minimizes errors by optimizing the problem. The PSO algorithm is a statistical process for optimizing parameter values. Particles cooperate, share information, and obey rules directly to determine the best solution to the optimization problem. This optimization process is a new approach to obtaining numerical solutions to coupled ordinary differential equations. The PSO algorithm is a technique to improve global search results and has a variety of parameter space characteristics. The application of the PSO algorithm is carried out by searching for approaching local optimal solutions through random initial conditions (Equation (21)), the optimal solution for one particle (individual particle optimal) and the optimal solution for all particles (group/global optimal) will converge to the local optimal solution. This allows to guarantee of the ability to find the global optimal solution, but the fast convergence ability will not be effective. Although the convergence of the PSO algorithm is fast, it may be trapped in a local optimum due to premature convergence. Therefore, a convergence speed controller (Equation (17)) is designed for the ultimate goal of solving this problem. Two adaptive approaches are proposed to adjust the convergence speed. The convergence speed is adaptively controlled to empower the PSO algorithm in solving large-scale numerical optimization problems. Currently, designing a proper convergence speed controller framework for the PSO algorithm has become a major challenge.

3.1 Sensitivity analysis of parameters variation of the kinetic model on cancer cell growth

This research estimates the parameter values contained in the kinetic model and then carries out a sensitivity analysis of these parameter values for the fastest mutation process from normal cells to cancer cells. The estimated parameters are listed in Table 2. The various values of the estimated parameters are selected based on the range of values given in Table 1. Except for the presentation of A cells, which are required to eliminate the angiogenesis limit, these cells will be tested in the range of 0 %–100 %, and genes involved in the transition will be tested with estimated values. The D:R contribution ratio parameter is analyzed with values in the range of Table 1 and is also analyzed with the inverse ratio value.

In this study, the contribution of each parameter to the growth and development of cancer cells was reviewed by estimating several parameters. Variations in parameter values are determined as listed in Table 2. When parameter values are being estimated, the values of other parameters are maintained at their default conditions. This study will review the values of optimized parameters that can contribute to the effect on the fastest path to cancer and also the fastest time to reach 10¹² M cells.

The calculation results obtained from the kinetic model regarding variations in parameter values are shown in Table 3. Table 3 shows that the results of estimating parameters from the kinetic model on cancer development are very sensitive to changes in parameter values. In Table 3 it can be seen that the values of the parameters tested are within the literature range for the kinetic model. The relative contribution parameter values of the D:R mutation, if the value of these parameters is changed to 8:2, do not result in a change in the time to reach 10¹² M cells from the fastest path to cancer.

Then, the calculation results that examine the variation in the relative contribution value of the D:R mutation are reversed to 3:7. The results of this calculation show the role of the D and R mutations. This result causes the fastest path mutation series to become the DRAG mutation series. These calculation results have also been proven by Spencer and coworkers in previous research [13]. If we look at the time to reach 10¹² M cells, the reversal of the value of the relative contribution of the D:R mutation makes cancer cell growth slower with an increase of 0.23 years (84 days). This proves that the role of the R mutation in the kinetic model of cancer development is greater than that of the D mutation.

The results of calculations by varying the estimated values for the gene parameters involved in the transition do not change the fastest path to cancer. The genes involved in the 500, 100, and 10 transitions only reduce the time by 1.02 years (372 days), and the genes involved in the 200, 20, and 2 transitions only reduce the time by 0.32 years (117 days). Therefore, these gene parameters do not have a major influence on the kinetic model of cancer progression. The results are the same as the results of previous research conducted by Spencer and coworkers [13].

A 10-fold increase in the value of k ₁ (cell mutation rate) reduces the time to reach 10¹² M cells by 0.86 years (314 days). Additionally, a 100-fold decrease in the value of k ₂ (the mutation rate of a cell after acquiring the G mutation) increases the time to reach 10¹² M cells by 0.15 years (10 days), and a 100-fold increase in time decreases the time to reach 10¹² M cells by 0.36 years (131 days). This shows that changes in the value of k ₁ have a greater influence than changes in the value of k ₂. This is because the G mutation occurs at the end of the fastest path to cancer, so the effect of k ₂ occurs at the end while k ₁ has an effect at the beginning of the path when it has not yet received the G mutation. So, the value of k ₁ will change to k ₂ when the cell gets the G mutation. Before getting the G mutation, the cell mutation rate had a value of k ₁. As for the fastest path to cancer, it can be seen that changing the value of k ₁ does not result in a change in the path while changing the value of k ₂ results in a change in the path when the value of k ₂ is reduced.

The value of cell division and death rate changed to 1 time every 5 days will reduce the time to reach 10¹² M cells by 2.2 years. Meanwhile, changing the value to 1 time every 20 days will increase the time to reach 10¹² M cells by 0.7 years (256 days). This means that increasing the value of the rate of cell division and death will accelerate the growth of cancer cells while decreasing the value of the rate of cell division and death will slow down the growth of cancer cells. However, this parameter does not have a greater influence than the tumor volume doubling time parameter. Concerning the fastest path to cancer, it can be concluded that changes in the parameters of cell division and death will result in changes in the fastest path from mutation of normal cells to cancer cells.

Overall, changes in parameter values only had a slight effect on the time to reach 10¹² M cells, except for the tumor volume doubling time parameter. This parameter has quite a big influence on the growth of cancer cells. Changing the parameter value to 200 will reduce the time to reach 10¹² M cells by 16.19 years. This parameter value also changes the fastest path to cancer to RDGA. In addition, changing the parameter value to 523 will increase the time to reach 10¹² M cells by 2.23 years, but the fastest path to cancer does not change. This shows that the change in the given parameter value is proportional to the resulting influence. This tumor volume doubling time parameter influences the term 1 / b − 1 / d . When a cell experiences a mutation in R and/or D, the terms in the equation will become 1/b _R and/or 1/d _D . The calculation to obtain the values 1/b _R and 1/d _D is listed in the attachment. Based on the results of these calculations, the tumor volume doubling time value is directly proportional to the b _R value and inversely proportional to the d _D value. It can be seen that when the tumor volume doubling time decreases, the rate of cell division with the R mutation will become faster. The rate of cell death with the D mutation will be longer so that the growth of cancer cells becomes much faster. This means that the greater the tumor volume doubling time, the slower the growth of cancer cells will be.

3.2 Kinetics of various pathways from normal cells to cancer cells

Figures 2–5 show the simulation results of the rate of cell population mutation within 100 years. Figure 2 shows that the order of speed of the mutated cell population begins with mutations R, D, A, and G. The population growth of the cells mutated with R is faster than that of cells mutated with D because the change in the rate of cell division (from b to b _R ) is greater than the change in the cell death rate (from d to d _D ). Therefore, the effect of adding a population of cells mutated with R from the term P R 1 b R − 1 d is greater than the effect of adding a population of cells mutated with D from the term P D 1 b − 1 d D . The cell population that mutates with G has the lowest speed because the increase in population is influenced only by P N M ̄ (normal cells), and the mutation rate when the population decreases increases by 1000 times (from k ₁ to k ₂). The increase in the population of cells mutated with A is also affected only by P N M ̄ (normal cells), but the mutation rate when reducing the population remains at k ₁.

Figure 2:

The population of cells with one type of mutation for the fastest path to cancer cells is (1) the population of cells with mutation R (blue solid line), (2) the population of cells with mutation D (green solid line), (3) the population of cells with mutation A (red solid line), and (4) the population of cells with mutation G (black solid line).

Figure 3:

The population of cells with two types of mutations for the fastest path to cancer cells is (1) the cell population with the DR mutation (red dotted line), (2) the cell population with the AD mutation (red solid line), (3) the cell population with the AR mutation (black solid line), (4) the cell population with the DG mutation (green solid line), (5) the cell population with the DR mutation (blue dotted line), and (6) the cell population with the AG mutation (blue solid line).

Figure 4:

The population of cells with three types of mutations for the fastest path to cancer are cell populations with ADR mutations (blue solid line), cell populations with DGR mutations (green solid line), cell populations with ADG mutations (red solid line), and cell populations with the AGR mutation (green solid line).

Figure 5:

The population of cells with four types of mutations for the fastest path to cancer cells is the cell population that has four types of mutations (T) (red solid line) and the cell population that has metastasized (M) (blue solid line).

Figure 2 also shows that the population of cells with R, D, and A mutations sharply increased until the 30th year. However, the increase in the population of mutated cells was not significant until the 62nd year for cells with R mutations, until the 57th year for cells with D mutations, or until the 45th year for cells with A mutations. After that, the population of mutated cells reaches a state of saturation. Moreover, the number of cells with G mutation increased sharply until the 2nd year. However, the increase was insignificant until the 10th year, when the cell population reached saturation. The saturation state means that the growth of the mutated cell population has reached its maximum limit.

Figure 3 shows the growth rate of the population of cells with mutations in DR, AD, and AR, which increased sharply until the 30th year. After that, an insignificant increase occurred until the 70th year for the cell population with DR mutations, until the 56th year for the population with AD mutations, and until the 57th year for the population with AR mutations. The population growth of cells mutated with DG and GR increased sharply until the 3rd year, and then, an insignificant increase occurred until the 16th year. After that, the increase in the cell population returned sharply until the 30th year. Then, before the saturation state was reached, there was an insignificant increase in the cell population until the 46th year for cells mutated with DG and until the 51st year for cells mutated with GR. The population growth of cells that were mutated with AG increased sharply until the third year. After that, an insignificant increase occurred until the 35th year, when the cell population reached a saturation state.

Figure 4 shows that the growth rate of the cell population changed with ADG, ADR, AGR, and DGR, which together increased until the 30th year. After that, the population of mutated cells reached a saturation state in the 60th year for mutations with ADR, in the 57th year for mutations with DGR, in the 52nd year for mutations with AGR, and in the 47th year for mutations with ADG.

Figure 4 also shows that the growth curve of the cell population with mutations in DGR, AGR, and ADG was greater than that with mutations in ADR, but the saturation point was lower than that with mutations in ADR. This is because cells with ADR mutations do not contain mutations in G, so their movement speed looks different from that of the other three cell populations. The three types of cell populations containing cells mutated with G increase the cell population of mutated cells containing mutations with G so that the initial speed increases rapidly. This is because the mutation rate with the G mutation increased 1000 times (from 10⁻⁷ to 10⁻⁴).

Figure 5 shows the growth rate of the T cell and M cell populations. T cells are tumor cells, namely, cells that have four mutations (RDAG), but these cells are still local. Tumor cells transform into M cells (metastatic cells) at a mutation rate of 10⁻⁹. These M cells can spread to other organs of the body and are also known to be cancer cells. According to the literature, tumor cells are cells that grow only in the same location and do not spread to other organs, while cancer cells can metastasize or spread to other organs [38].

Figure 5 also shows that the growth of the T-cell population was faster than that of the M cell population. This is because the T-cell population is dominated by a population of ADG-mutated cells that have R mutations, ADR-mutated cells that have G mutations, DGR-mutated cells that have A mutations, and AGR-mutated cells that have D mutations. The rate of mutations is 10⁻⁴ for each cell division, except for ADR mutation cells, which have a rate of G mutations of 10⁻⁷ for every cell division. This explanation is shown by Equation (16), while the cells that enter the M cell population are only T cells at a rate of 10⁻⁹ for each cell division, according to Equation (17).

Simulating cancer growth from this kinetic model produces the fastest path from normal cells to cancer cells that acquire the mutation sequence RDAG. This can be seen from Figure 2, which shows that mutation R has the highest speed. Figure 3 shows that the DR mutation has the highest speed. Figure 4 shows that ADR mutations have the highest speed. Based on these simulations, this study concludes that the fastest path to cancer cells is obtained by obtaining mutations R, D, A, and G. The sequence of mutations on the fastest path to cancer in this study is not the same as that in previous research conducted by Spencer and coworkers [13], namely, the fastest path to cancer cells in the results of their research is through the mutation pathway sequence DRAG. This difference is due to differences in the calculation of the parameter values b _R (cell division rate with R mutation) and d _D (cell death rate with D mutation) from the research of Spencer and coworkers [13]. These two parameter values influence the speed of the mutation populations D and R. The b _R parameter value set by Spencer and coworkers [13] is greater than the actual calculation results, so the speed of the R mutation cell population in the previous study was slower than that in this study. In addition, the d _D parameter value set by Spencer and coworker [13] is greater than the actual calculation results, so the speed of the D mutation population in the previous study was faster than that in this study (Figure 6).

Figure 6:

The fastest path to cancer in this study.