Mathematical modeling reveals the factors involved in the phenomena of cancer stem cells stabilization

Summary

Cancer Stem Cells (CSC), a subset of cancer cells resembling normal stem cells with self-renewal and asymmetric division capabilities, are present at various but low proportions in many tumors and are thought to be responsible for tumor relapses following conventional cancer therapies. In vitro, most intriguingly, isolated CSCs rapidly regenerate the original population of stem and non-stem cells (non-CSCs) as shown by various investigators. This phenomenon still remains to be explained. We propose a mathematical model of cancer cell population dynamics, based on the main parameters of cell population growth, including the proliferation rates, the rates of cell death and the frequency of symmetric and asymmetric cell divisions both in CSCs and non-CSCs sub-populations, and taking into account the stabilization phenomenon. The analysis of the model allows determination of time-varying corridors of probabilities for different cell fates, given the particular dynamics of cancer cells populations; and determination of a cell-cell communication factors influencing these time-varying probabilities of cell behavior (division, transition) scenarios. Though the results of the model have to be experimentally confirmed, we can anticipate the development of several fundamental and practical applications based on the theoretical results of the model.

Keywords:

Cell cycle and cell division – Population dynamics – Stem cells – Cell death – Mathematical models – Stem cell therapy – Cancer stem cells – Tumor stem cells

Introduction

Stem cells are undifferentiated cells present in very low numbers in most tissues. Stem cells are responsible for tissue renewal and homeostasis, by giving rise to non-stem cells that proliferate and further differentiate in specialized cells. Stem cells show very specific features, notably regarding cell division: they are able to undergo asymmetrical division, dividing into a stem cell and non-stem cell; moreover, the rate of stem cells division is very low as compared to that of non-stem cells [1–3].

It has been demonstrated that in most malignant tumors, cancer cell populations appear to include a rare stem cell-like subpopulation suspected to be responsible for the initiation and maintenance of tumors in animals [4–14]. This subpopulation can be detected and purified using specific cellular probes or cell surface markers. In vitro, purified cancer stem cells (CSCs) are able to reconstitute the population heterogeneity whereas, in contrast, purified non-stem cells cannot. Also, CSCs were shown to be highly tumorigenic in xenografts experiments, and to be responsible for cancer metastasis, i. e. colonization of various tissues by the primary tumor. Because of these features, cancer stem cells are also called tumor-initiating cells [7]. However, not all cancer stem cells appear to have all the features of normal stem cells. For example, CSCs may have a diminished capacity to undergo asymmetrical division compared to normal stem cells [6,12,15–17].

CSCs have been demonstrated in most solid and hematologic tumors [8,10,16,18–28], with very well described common functional features, i. e. an indefinite self-renewal capability and the ability to undergo some asymmetric divisions. Importantly, CSCs are resistant to chemo- and radio- therapies [4,5,9,10,29–31], suggesting that they may be responsible for tumor relapses following chemo- or radio-therapy. This has important implications in therapeutic, as most of the current treatments target a regression of the tumors mass without accounting for the tumor functional heterogeneity.

Various mathematical models have been proposed for describing the dynamics of both normal [32–36] and cancer [30,31,37–47] stem cell populations behavior. These works suggest two different concepts for description of stem cells population behavior. One concept is based on the principle that stem cells act according to their intrinsic program, which may be deterministic or stochastic [30,31,35,39–41,45,48]. Alternatively, a concept of self-organization of stem cells suggests modeling of the entire system of cell-cell and cell-environment interactions, for which some authors also consider a stochastic behavior [35,36,49–51]. Modeling of cancer cell population behavior provides very important inferences as for understanding the nature of cancer growth so for clinical prognosis and treatment strategies. In many cases it allows to have evidence about factors that cannot be measured directly in clinical or in experimental investigations. For example, due to the mathematical model it was shown that the role of leukemia stem-like cells population on the course of disease is much greater than the one of leukemia blast cells [42,43,47]. The model-based estimation of prognostic factors in clinical data may help in designing treatment strategies [42–47]. Moreover, the quantitative information about the cancer stem cells (CSC) or tumor-initiating cells (TIC) fraction dynamics can be inferred by methods based on mathematical models. The method described in [44] identifies two characteristic equilibrium TIC regimes during expansion and regression of chronic myeloid leukemia.

However, the possibility of dedifferentiation of non-stem cancer cell to stem cell is one of the most important and yet unsolved question about CSC behavior, which has been, so far, addressed in only few of these works [30,38,52]. These (stochastic) models base their principle on the assumption that all possible transitions in subpopulations can occur spontaneously, with some probabilities.

Also, only few modeling approaches were proposed to gain insight into the intriguing phenomena of cancer stem cells population stability [38,44,53,54]. This phenomena detected in many cancer cell lines harboring measurable levels of cells with CSC features, is that over several years of cell passage the relative number of cancer stem cells fluctuates around a basal level, characteristic for each specific cell line (as illustrated in Fig 1, dotted red curve). Moreover, it has been shown that isolated cancer stem cells can rapidly regenerate in culture the heterogeneity of the parental cell line with the characteristic relative percentage of cancer stem cells (as illustrated in Fig 1, dark blue curve).

One work discussing this phenomenon models the CSC behavior as a Markov process [38]. The model is based on stochasticity of single-cell behaviors and does not consider cell-to-cell communications.

In our previous work [53,54] we constructed and analyzed a mathematical model that takes into account this intriguing characteristic of CSC population behavior. We suggested an instructive role of cell-to-cell signaling influencing the cell parameters and leading to CSC population equilibrium. The mathematical model accounts for all possible cancer stem and non-stem cell behaviors, i. e. type of division (symmetric or asymmetric), direct transition (differentiation or dedifferentiation) and cell death. The analysis of the model helped to elucidate some important characteristics of cancer stem cells evolution, in particular, a set of parameters of cell growth implying the necessity of non-stem to stem cell transition.

In this work we expand this mathematical model and address the question of “instructive signal(s)” underlying the phenomena of cancer cell population stability, aiming to provide meaningful predictions on its dynamics and nature. In the presented work we continue analysis of the model aiming to solve the following problems:

- determination of time-varying corridors of probabilities of different cell fates, given the dynamics of cancer cells populations;

- determination of a cell-to-cell communication factors, influencing time-varying probabilities of cell behavior (division, direct transition) scenarios.

We demonstrate that using data measured in the context of CSC population stabilization, our model is able to infer corridors of time-varying probabilities of cancer cell fates that provide significant insights into the cellular dynamics of heterogeneous tumors. Next we show how the set of curves of probabilities can help identifying a set and kinetics of secreted factors responsible for cell population behavior.

Methods

• Algorithm for the solution of the system of Eqs (14–17).

The system can be rewritten in the form:

A system (A1-A6) was approximated by a second order finite difference scheme

We consider that the initial conditions:

Then, the initial values for pi0,qj0,(i=1…4,j=1,2) are set by (A13) and the calculations from (A3-A5), using (A13). By that, from numerical solution of the system (A7-A12) we can obtain pin+1,qjn+1(i=1…4,j=1,2) for any time step n = 1,2,…. Accuracy of numerical simulations is controlled by decreasing time step.

Thus a triple of parameters p₁₀, p₂₀, q₂₀ allows obtaining pin,qjn for any time step. If at some time step the conditions

A user interface is developed allowing to input the parameters of the model and to output the results in the graphical form and in the form of numerical files.

The corresponding program code can be found at:

https://github.com/nickbessonov/CSC-article/blob/master/CSC%20(17.05.18).rar)

• 2

  Algorithm for the solution of the system of Eq (25).

The system of Eq (25) was solved numerically by well known trick when the system of steady equations was replaced by a system of simplest non-steady equations and the solution of a steady problem is sought by solving the non-steady problem and finding its steady solution.

In our case, for determining b_ik,a_k,c_k, k = 1…K, i = 1…6, the system Eq (25) was replaced by a system of the next non-steady equations

The system (B1) was solved with the help of iterations. Trial initial values aq0,cq0, and biq0 q = 1…K, i = 1…6 were set. Substituting these values into the right-hand side of the Eq (B1), we determined the values of unknown parameters in the next pseudo time step aq1,cq1, and biq1. Then this process was repeated until the required accuracy of the solution was reached. The only one difficulty arises at the beginning of application of this trick when it is necessary to determine the signs of time steps. This problem was solved easy by several numerical experiments. The sign minus in right part of (B1) is a result of these computational experiments. The corresponding program code can be found at:

https://github.com/nickbessonov/CSC-article/blob/master/CSC%20(17.05.18).rar)

• 3

  Experimental part

The measurements of cell division rates were performed on several human cancer cell lines, namely five NSCLC (non-small cell lung cancer) cell lines (A549,HCC827, H1299, H1650, H1975), four ovarian cancer cell lines (A2780, HEY, OV2008, TOV-112D cells) and colon carcinoma HCT 116 cell line. The detected proliferation rates of non-stem cells were similar or considerably close in all investigated cell lines, thus its averaged value was taken for the numerical simulations.

The percentage of cell death rates was measured by Fluorescence Activated Cell Sorting (FACS) analysis in A549 cell line at every day during 10 days after passage. The dead cells were stained by 7 AAD (7-AminoActinomycin D) and sorted accordingly. The analysis have shown the rate of cell death around 0.1 (10% of all cells) for all days during the period of cultivation. The second variant of death rate (0.5) was added to the numerical simulations to explore another possibility of death rate, close to a critical one, which was noticed in another investigated cell lines (see above), but without precise measurements.

Results and discussion

Mathematical model

Following our previous work [53] we suggest a model accounting for 4 scenarios of cell behaviors (different types of cell divisions and direct transition) for stem and non-stem cells, and assumed that each scenario can occur with some probability (Table 1):

Though we include the scenarios D→S+S and D→S for daughter cells in Table 1 as theoretically possible, for a time being we will consider in a model only two scenarios for daughter cells, which discriminates two principal modes of their behavior: the regular one (D→D+D) and the de-differentiation of daughter cells into S cells, which requires the activation of different genetic program. For the representation of the last case we consider only one scenario D→S+D (for which we assume the transition D→S as a part of it, meaning that D cell should undergo a transition D→S before the next asymmetric division). Thus, for the current model the probabilities q₃ and q₄ are taken to be equal to 0.

The probabilities should satisfy the usual restrictions on probabilities:

We considered that cell divisions/transition occur with the rate λ₁ in stem cells (S) and with the rate λ₂ in daughter cells (D), and that death rates are γ₁ and γ₂ in S and D cells respectively. We need to comment that though it seems more plausible to consider that each mode of cell divisions or transition occurs with its own rate (e.g., like in models [42,43], we will consider in our model one constant rate of cell divisions (or direct transition) for stem cells, and another constant rate of cell divisions for non-stem cells, independently from a scenario. Also, though it was shown that in general cell division rates may change in time, for example, due to regulation by signaling factors [55–57], we assume the possibility to consider them to be constant during the short period of stabilization of cell population. The same assumption is made in many mathematical models addressing the cancer cells population behavior [42–47].

The proposition (Table 1) gives a system of differential equations for the dynamics of S and D cancer cells:

Using (5), the system (4) can be re-written as:

In our previous work, using the system (6), we have analyzed the time-dependent evolution and the asymptotic behavior of the percentage of cancer stem cells s(t) in a cancer cell population with the initial conditions:

The fraction (percentage) of stem cells is calculated as:

Corridors of probabilities of cell fate scenarios

Here we continue analysis of the model aiming to determine the time-varying corridors of probabilities of different cell fates, given the dynamics of cancer cells populations. We assume that a curve s(t) of the percentage of CSC over time displaying a CSC stabilization after perturbation (as shown at Fig 1), and proliferation and death rates of cells are given as a result of experimental measurements in a particular cell line.

First, our numerical simulations show that a good fitting of model (6) to any reference curve s(t) (determined by the least squares method, with an average deviation equal ~ 5% or less) can be achieved in admission that the set of parameters λ_i,γ_i,p_i,q_i changes over time. As discussed in Section 1, we may admit for the following analysis that the parameters changing with time should be the probabilities of scenarios p_i(t),q_i(t), while parameters λ_i, γ_i, may be considered as constant.

This led us to Question 1:

Given the dynamics of percentage of CSC s(t), is it possible to find functions p_i(t), i = 1,2,3,4, q_i(t), i = 1,2 for a given constant set of λ_i and γ_i?

In order to solve this problem we first considered the following hypothesis: all changes in cell behavior scenarios up to stabilization should be minimized (the sum of all changes of p_i(t) and q_i(t) should be the minimal possible ones). The underlying biological logic of this hypothesis is that each type of scenario is implemented by different sets of biological mechanisms, which require specific proteins and other compounds to be prepared in a cell. This means that each time-step of cell population behavior towards stabilization should be more plausibly achieved by minimal change in each scenario existing at the previous time-step than by intermittent pattern of each of them.

This means that the function v:

Next we rewrite the system (6) and the conditions (2) and (3) as:

This gives us a system:

Thus, for six variables p_i(t) and q_i(t) we have a system of 6 Eqs (14–17). For the solution of that system we have to add three initial conditions:

All solutions for each p_i(t), q_i(t) corresponding to all possible sets of initial conditions (18), which satisfy the condition (1) for all time period T will provide a set of six corridors (ranges) of probabilities of scenarios p_i(t), q_i(t) of cell behavior in a given experimental case.

This means that, given the measured s(t) and a constant set of λ_i,γ_i, we can determine corridors of all possible probabilities of scenarios p_i(t), q_i(t) varying with time.

The results of the calculations of all possible corridors of probabilities for the reference curve s(t) (taken as the one in Fig 1), and different sets of biologically relevant parameters λ₁,λ₂,γ₁,γ₂ are presented in Fig 2; the corresponding algorithm and the program code are provided in Methods section. The considered rates of non-stem cell division (λ₂) were taken from experimental measurements in different cancer cell lines, while the rate of stem cell division (λ₁), which is difficult to measure, was chosen in order to vary the ratio of stem/non-stem division rates (1:1; 1:5; 1:10), corresponding to a statement about slower rate of stem cells divisions. The considered variants of cell death rates (0.1; 0.5) were taken from experimental measurements.

It is important to note, that the boundaries of each corridor can be (and usually are) comprised from the parts of different possible functions p_i(t) (or q_i(t)) existing inside a particular corridor. However, these boundaries can help to evaluate the possibility or impossibility of a scenario(s) in particular experimental systems.

For example, this simulation gives an answer to the question of the necessity of non-stem to stem cell transition scenario in the course of CSC stabilization toward equilibrium. As it can be seen from the Fig 2, for all considered biologically relevant sets of cell division and cell death rates λ₁,λ₂,γ₁,γ₂, the lowest boundary of the corridor of probability q₂(t) (corresponding to non-stem to stem cell transition) appears to be higher than zero, at all or at least at some time points.

Another important question is the existence of other types of stem cell behavior scenarios in addition to asymmetrical division, which is postulated to be predominant at least in normal stem cells. Results showed that for biologically relevant sets of parameters λ₁,λ₂,γ₁,γ₂, the highest level of possible probability p_i(t) (the upper corridor boundary) is around 60%, being much lower for some cases (Fig 2).

For many practical applications it is very important to be able to find a unique solution for each probability p_i(t), q_i(t). To this end, we addressed Question 2:

Given a measured curve s(t), and parameters λ_i,γ_i, what additional data are necessary and sufficient to get a unique solution for each p_i(t) and q_i(t)?

From the analysis done for the Question 1 it is clear that these additional data should be a set of three initial conditions (18).

We have explored this possibility and present two examples of the determination of such a set of unique functions p_i(t) and q_i(t) in Figs 3 and 4, considering two important biological questions, mentioned above. In the first simulation (Fig 3) we chose the initial conditions (18) such that q₂₀ be the minimal one, thus exploring the possibility of an absence of D to S cell transition the scenario for each particular case (six different cases (sets of parameters λ_i,γ_i,as in Fig 2) are considered). As we can see from our results, this possibility does not occur in any of these cases. In the second simulation (Fig 4) we search for the maximal probability of asymmetric division of stem cells in each particular case, by setting the initial conditions to maximize p₁₀. The results show that for all six cases the probability p₁ can not exceed 63%.

Determination of underlying field factors involved in cancer cell population dynamics

In our previous work [53] we have suggested that the coordinated dynamical change of the parameters of cell behavior, resulting in cancer stem cells subpopulation stabilization, occurs in response to a multiparametric biochemical signal produced in a system. In the simplest case, it may be a set of secreted factors influencing cell behavior; although it can also be induced by cell-to-cell contact, a hypothesis that is not taken into account here for the sake of simplification.

We noted this signal as an underlying field u(t), where u(t) is a set of biochemical compounds (factors) in a media, capable to influence cell fate. Generally, in can be presented as a matrix of factors u_ij(t), where each factor i may be produced by S or D cell at some time point t, and have an influence j on the behavior of S or/and D cell, possibly with some time delay t+r.

Here we further explore this idea, assuming that a field of secreted factors u(t) influences probabilities of cell fate scenarios, and thus that the structure of corridors of probabilities p_i and q_i depends on a changing set of concentrations of secreted biochemical compounds.

This means that in the Eq (10) all probabilities of division scenarios of S and D cells (p_i(t) and q_i(t) are the functions of an underling field u(t), changing with time:

The task of identification of molecular factors involved in the underlying field formation raises the Question 3:

Given a set of unique functions of probabilities p_i and q_i, is it possible to find a set of factors u(t) responsible for their dynamics?

First, for simplification, we will consider underlying field u(t) merely be a set of K biochemical compounds (factors) presenting in the media, i.e., as a set of functions u_k(t), k = 1,…,K. Later, we may consider the possible dependence of each factor u_k(t) on the amount of S and/or D cells, as these factors may be produced by one and/or the other type of cells.

Next, in order to determine the factors u_k(t) which influence the evolution of probabilities p_i(t) and q_i(t), we will perform decomposition of given functions p_i(t) and q_i(t) over the functions u_k(t). We will use the following form for the function u(t):

This form, by varying the coefficients a,b,c,m, allows modeling various kinetics of factors, e.g., linear with different parameters, constant, exponential, trapezoid, etc.

According to formula (20), the dependence of probabilities p_i,q_i on factors u_k(t) can be written as:

Therefore we may define function u_k(t) as

The test calculations and numerical simulation have shown that it is possible to set the value of the coefficient m = 6 for all factors, while coefficients b_ik,a_k,c_k should be selected by least-squares method. Thus, we will find such coefficients b_ik,a_k,c_k in the expression (21) that provides a minimum of the total quadratic deviation:

For that it is necessary to ensure equality to zero of the expressions:

By this we obtain a system of Eq (25) for determining b_ik,a_k,c_k,k = 1…K, i = 1…6, which was solved by the gradient descent method, starting with test initial values and obtaining the minimum of the function (24). The number of factors k is determined from the condition that the mean deviation y_i(t_n) from y_in should not exceed a given value of permissible variation ε_p:

The search starts from suggested initial number of factors K, and in the case when condition (26) is not fulfilled for a given ε_p, should be repeated for K+1, and further on until the satisfaction of (26).

Next a corresponding computer program, automatically performing all necessary computations for a chosen ε_p, was developed and applied for the particular set of curves p₁(t),p₂(t),p₃(t),p₄(t),q₁(t),q₂(t) presented on the Fig 3 (the corresponding algorithm and the program are provided in Methods section). The results of computations, determining the number and pattern of u_k(t) factors, providing the satisfaction of the condition ε ≤ 1% for each experimental case (the same sets of λ_i,γ_i considered before), are presented on Fig 6, where the height of each factor b_k.

The search started from suggested initial number of factors K = 4, and has elucidated the most plausible number of major u_k(t) factors for each experimental case (which varied from 3 to 6 factors depending on the case), together with their specific patterns (shape, height and timing).

It was found that in 3 cases out of 6 explored, the starting number of factors K = 4 was enough to determine a set of essential influencing field factors u_k(t) with sufficiently good value of variation ε.

Moreover, in two cases, namely, for the set λ₁ = 1,λ₂ = 1,γ₁ = 0,1,γ₂ = 0,1 (ε = 0,49%; Fig 6A) and for the set λ₁ = 1,λ₂ = 10,γ₁ = 0,1,γ₂ = 0,1 (ε = 1%; not shown), the program has detected only 3 factors as essential, thus showing the capability to find the minimal possible set of factors u_k(t) for the best fitting, independent of the starting number of factors K. In one case 4 factors were determined: λ₁ = 1,λ₂ = 1,γ₁ = 0,1,γ₂ = 0,5 (ε = 0,47%, Fig 6B). For the other 3 cases, for which with K = 4 unsatisfactory ε value was obtained, the search was continued with greater number of factors K up to optimal of fitting. Finally, for two cases the final number of essential influencing factors was 5: λ₁ = 1,λ₂ = 5,γ₁ = 0,1,γ₂ = 0,1 (ε = 0, 82%; Fig 6C) and λ₁ = 1,λ₂ = 5,γ₁ = 0,1,γ₂ = 0,5 (ε = 1,0%; Fig 6D), while in one case, for which K = 5 gave ε = 1,2, the final number of factors K appeared to be equal to 6: λ₁ = 1,λ₂ = 10,γ₁ = 0,1,γ₂ = 0,5 (ε = 0, 75%; not shown).

We expect that this model can give insight into a search for biochemical factors involved in controlling cell behavior. For example, if in a course of measurement of CSC population kinetics, also the Metabolome (Secretome) profiles of the culture media can be obtained, the secreted factors which appear, increase, decrease or disappear at each time point can be identified (Fig 7).

Next predictions on u_k(t) factors obtained by mathematical model for a given case (s(t), λ_i, γ_i) (Fig 6) should be compared with the results of the corresponding Secretome profile (Fig 7). The biochemical factors having the kinetics and patterns coinciding with the predicted underlying field functions u_k(t) could be the good candidates for being involved in controlling cell behavior.

It is probable that out of possibly several hundreds metabolites in Secretome, several will fit to predicted u_k(t). In this case two variants may be envisioned: either a cumulative effect of all (or several) indicated metabolites creates the u_k(t) signal, or only one among them has this special effect.

Important Remark. Though the system of Eq (25) cannot be solved uniquely, our computational experiments have shown that by requesting ε to be considerably small and by setting the approximation time steps in the program to be 10⁻⁴ arbitrary units, we get minimizing of the area of all u_k(t) in such a way that the difference between the various solutions is negligible, as it is below the sensitivity of experimental equipment used for monitoring the factors. For example, for ε≈1%, the maximal variation of coefficient b reflecting a concentration of a factor in a media, is 3%, and maximal variations of coefficients a and c reflecting the time of existing of the factor in media are 3%.

This means that the model suggests an approach to determine the secreted factors influencing the time-dependent cells behaviors leading to the stabilization of CSC population.

Predictions of putative production of secreted factors u_k(t) by S or D cells

The identification of the cell type, stem or non-stem, producing the factor(s) influencing CSC behavior might be important for basic research as well as for potential medical applications. In order to provide a tool for such a prediction, we assumed that some factors u_k(t) can be produced (completely or mostly) by S cells, others by D cells, while some factors may be without explicit belonging to any type of these cells.

To account for this possibility, we have considered a set of functions U_kc(t) which provides discriminatory dependence of a factor u_k(t) on S or D cells kinetics:

This means that having an optimal set of factors u_k(t), found at the previous step (Fig 6) for a given set of parameters λ_i,γ_i and given kinetics p_i(t), q_i(t),we will consider all variants of multiplication of each function u_k(t) on each of the functions s(t),d(t), and on the constant equal to 1, with determining ε value for each U_kc(t) according to formula (26).

We can assume that the result of such a multiplication that gives the smallest ε value represents the possible dependence of the factors upon the types of cells (S or D cells). The variant of multiplication on 1 means independence of a factor on S or D kinetic, which can reflect the production of these factors evenly by both types of cells or their production by the environment.

The results of the computation U_kc(t) for two selected cases are presented in the Table 2.

Two chosen experimental cases (A) and (B) for different sets of parameters λ₁,λ₂,γ₁,γ₂ are presented, with four secreted factors u_k for each of them. In Table 2, “S” or “D” means considering the dependence of a factor on S or D cells kinetics, while “1” means considering an independence of a factor on any specific type of cells.

The variant (1,1,1,1) marked in bold corresponds to the independence of all factors upon any specific type of cells. All variants with better results (ε is smaller than in (1,1,1,1)) and some chosen variants for the worse results (ε is larger than in (1,1,1,1)) are shown.

It can be seen that for the case (A): λ₁ = 1,λ₂ = 5,γ₁ = 0,1,γ₂ = 0,1, for K = 4 and without differentiating multiplication (corresponding to the variant 1,1,1,1) ε = 1.373% (Table 2A). The computation has shown that all other variants of multiplication give worse results (some selected ones are shown), except for only one variant (S,D,D,S), though decreasing the value of ε not considerably (ε = 1.298%). However, it still can be an indication of differential production of the factors (t) by S or D cells.

Much stronger dependence on S and D kinetics is demonstrated in the other case (B): (A) λ₁ = 1,λ₂ = 10,γ₁ = 0,1,γ₂ = 0,5, where for K = 4 and without differentiating multiplication the lowest ε value was found as ε = 1.435%, while in the best variant with differentiating multiplication (S,D,1,S), it decreases up to 1.013 (Table 2B). Also, the fact that in this case 25 variants with differentiating multiplication give better results than the variant without considering the dependence on S or D kinetics (1,1,1,1), can indirectly point on the necessity of concluding this additional calculation in this concrete case. On the other hand, this fact can be an indicator of instability of the found best variant of multiplication (for this case, all variants with better results (ε is smaller than in (1,1,1,1)) and some chosen variants for the worse results (ε is larger than in (1,1,1,1)) are shown in Table 2B).

It is important to note, that in order to show a larger difference in multiplicative and non-multiplicative outcome, we present in Table 2 the results of a differentiating multiplication, which was started from the optimal sets of non-multiplied factors u_k(t) obtained with the requirement ε≤1.5% corresponding to K = 4. The program started from the best non-multiplicative set obtained in stricter requirement ε≤1% (presented on Fig 6, and corresponding to K = 5 factors in the case (A) and K = 6 in the case (B)), gives the same qualitative results, i.e., the same dependence of particular factors on S and D kinetics, but with less significant decrease of ε value for the resultant best multiplication.

This means that in order to obtain reliable results about the dependence of the underlying field factors on S or D cells, one should apply the suggested computational program starting with the optimal sets of non-multiplied factors u_k(t) obtained for possibly high level of ε_p.

Conclusions

Using a set of theoretical assumptions and basic knowledge of cancer cells population behavior, we suggest a mathematical approach which may help experimentalists gain insight into a broad array of cancer cell fates using a limited amount of experimental measurements. The computational program based on our model allows determining several important characteristics of CSC behavior in a cancer cell population.

First, on the basis of a minimal set of experimentally measured values (rates of cell division, of cell death and CSC population kinetics s(t)), our model enables the prediction of probabilities of cell-line specific modes of S and D cells divisions, including the unusual D→S transition. Second, the model explores the dynamics of cell-cell interaction factors u_k(t) influencing cancer cells behavior (namely, the time-dependent probability of cell division modes p_i(t) and q_i (t)). Finally, the model allows to concern the likelihood of production of these factors u_k(t) by S (stem) or D (non-stem) cells in the cancer cell population.

A potential drawbacks of the proposed model is that the key assumptions of the model are pure theoretical and that the model is rather simple. For example, we do not take into the account the effects of deterministic-stochastic interplay which was first considered for CSC modeling in [37] and next developed in [58]. Another obvious limitation of the model is that the suggested applications of it still have to be confirmed by biological experiments.

However, we can anticipate the development of several epistemological and practical applications regarding cancer cell behavior, which cannot be accomplished by using biological methods only. As an example, the proposed work may give insight into evaluation of the risks of cancer to relapse upon radio- or chemotherapy, based on the experimental measurement of its CSC kinetics. Also, it may help to find the treatment achieving the most efficient suppression of cancer subpopulations, as well as the schedule achieving the best therapeutic results by considering the predicted underlying field behavior. An essential part of this schedule may depend on identifying the treatment time points at which elimination of specific underlying field factors will be predicted as crucial for cancer population abolishment.