Fractional order modeling and analysis of dynamics of stem cell differentiation in complex network

1.   Introduction

Stem cells are central to regulate the body mechanism which governs development and tissue regeneration. These cells are self-improving with the capacity to both differentiate and multiply into two hundred cell types that make up the human body. Stem cells are primitive that have the potential to differentiate or develop into a variety of specific cell types. The two main varieties of stem cells in the biological system are embryonic and adult and have a unique capability to develop into well-known cell types in the human body. Moreover, these cells have an important role in the synchronization of the development of human body, tissue homeostasis as well as tissue regeneration. Their potential to differentiate and self-improvement distinguishes them. These two procedures are complementary in several ways in biologically dynamical processes. The self-regeneration ability refers to the stem cells that can create a new cell with similar characteristics to their progenitor, guaranteeing the survival of the basic stem cell population. However, differentiable measures the process by which these offspring become additional specialized and lose stem cell properties, in particular, the potential for self-regeneration ^[14].

Likewise, a stem cell, a progenitor cell is also a biological cell. Progenitor cells tend to differentiate into a fixed type of cell. But these cells are already more fixed than stem cells and are impelled to differentiate into their presence ^[7]. The main important difference between progenitor and stem cells is that progenitor cells can differentiate into more than two types of cells based upon their origin and ability but stem cells can reproduce indefinitely ^[21]. So, the conflict about the accurate definition remains and the abstraction is still in progress ^[22]. The progenitor cells are early descendents of stem cells that can differentiate to form one or more kinds of cells but they cannot be divided and reproduce indefinitely. There has been very small proof of this oligopotent concept of cells ^[8]. But progenitor cells speak of being in an extra stage of cell differentiation. Sometimes these cells lie in between stem and fully differentiable cells. This type of ability depends on the differentiation of their big stem cells and likewise on their function. But some progenitor cells could move through the body and migrate towards the tissue where they are needed. So many properties are captured by adult cells as well as progenitor cells.

Marciniak-Czochra and some other mathematical modeller proposed compartmental mathematical models of a discrete group of cell sub-populations to look into possible procedures for stabilization as well as regulation of red blood cell production after perturbations like bone narrow transplantation, which we used as a starting point for our research. Different forms of non-linearities in the mathematical model equations result from various conceivable regulatory feedback processes ^[18,6]. It was studied that to maintain such a stabilization cell system, two concurrent estimation populations are required ^[19].

The model provided with six compartments is an isolated negative response curve that describes how mature cells regulate the entire process ^[3]. In particular, the theories about how external signaling regulates hematopoiesis in feedback to a scarcity of mature cells are being investigated. The first hypothesis argues that the differentiation of cells is controlled solely by increasing the rate of accumulation, whereas the second hypothesis assumes that the prosecution of self-regeneration until differentiation is influenced by outside signals. A unique non-positive regulatory population (matching to hypotheses first or second) is potentially enough to manage the organization. The mathematical modeling in biology has been performed in various age structures ^[23,24]. The modelling of the COVID-19 pandemic has been observed in the environment ^[25,26]. We hope to validate this discovery in this work by performing a qualitative analysis of four compartmental mathematical models under the hypothesis first. The structure under Hypothesis Two was examined as well as requirements for the existence of partially free equilibrium points and positively for a four-compartment mathematical model, which have been established in the work of the authors ^[27]. The following properties of a stem cell population were discovered through mathematical and computational analysis of the biological system of equilibrium points: $(i)$ The mortality rate is lower than the basic reproduction rate for some cytokine levels, and $(ii)$ the signal intensity level of cytokine that is required to maintain the population density is lower than for all supplementary cell populations ^[16].

In this study, we mainly concentrate on circulation relating to the cell differentiation procedure and its stability analysis. In a few stabilizations cell models, like the cascade of circumstances originating in stem cells earlier than in thoroughly differentiable cells, the hematopoietic system is well outlined ^[2]. On the other hand, in some blood tissues like the erythrocytes, thrombocytes, and leukocytes, the condition is almost transparent as the total upcoming stages are unknown ^[13]. As a result, we use mathematical models with as well as without progenitor cells to investigate the impact of intermediary stages among stem, progenitor, and mature cells on stem cell dynamics. For the development of a new concept in differentiation in cell dynamics, incorporated the concept of a myeloid cell. This is the fresh work in cell division modeling through fractional calculus. In hematopoiesis, myelogenous or myeloid cells arise from the progenitor cells after their differentiation ^[9,10]. They were differentiated from descendants of common progenitors that derived from hematopoietic cells. Also, these cells(common progenitors) are the well-known defenders against infection in the human body.

So, based on cell differentiation facts, we concentrate briefly on the stability analysis of the characterization of stem cells. Here are two important questions still open: $(i)$ Does the stem cell collection remain constant throughout the life of a living entity$?$ $(ii)$ Does it decay the age of a living entity passes$?$ So some moderation of the model from to explore the dynamical behavior of the system less than the satisfaction of all cells, together with stem, progenitor, myeloid and mature cells, explode into replicative senescence has been given in ^[4,11]. Based on this understanding of differentiation in cell dynamics class in the proposed model. For the implementation of this class, the mathematical model is constructed by using the four compartmental model. The numerical simulations of the proposed model abcd utilizing the estimated data were presented.

This pipeline is described with the assistance of differential as well as delayed differential equations. But in this work, we apply fractional differential equations due to their good explanations in the real world. The mathematical modeling approaches through the fractional differential equation have gained interest as well as another service in the last few years, in the researcher's opinion. This is due to the concept that fractional-order mathematical models have the capability of gaining the crossover dynamical behavior and fading memory that is shown by the biological process. Most of the absorbing facts are the data fitting that can be described as non-integer order rather than ordinary-order ^[5,20]. Moreover, some fractional derivatives that are identified in their kernels, such as singular and non-singular, are utilized very efficiently to address real-world problems ^[15]. The fractional differential operator mathematical model generated through the Atangana-Baleanu derivative in the Caputo sense is more important and explains the dynamics systems more effectively than the Caputo-Fabrizio because of its derivative introduced at the basic point of departure in the Mittag-Leffler kernel, which is non-local and non-singular ^[1]. Thus, the mathematical models created through the Atangana-Baleanu derivative in the Caputo sense are more important and better survey the dynamics of an epidemic. In the ongoing research work, many well-known researchers have established mathematical models with fractional-order derivatives, which are proposed as a good tool for the current life problem response. We extended the model proposed by ^[28] by using the Atangana-Baleanu derivative in the Caputo sense into 'abcd' cell model by incorporating the myeloid cell class which is an important class of stem cell.

The aim of this mathematical model is to simulate the concentration of stem and progenitor cells. The remaining work on the paper is arranged as follows. Section 2 is devoted to the basic preliminaries related to non-integer fractional derivatives which are utilized in the formulation of mathematical model considerations. The mathematical model description in the Atangana-Baleanu derivative in Caputo sense is given in Section 3. Some basic properties of the mathematical model are given in Section 4. The analysis work is done in Section 5. The numerical solution is carried out in Section 6 by using the Adams Bashforth Moulton method. The numerical results are facilitated in Section 7. Finally, the conclusion is drawn in Section 8.

2.   Preliminaries on fractional calculus

In this section, we present some valuable concepts and basic definitions of fractional operators.

Definition 2.1. For $\varphi(t) \in C^{v}$, the Caputo derivative is defined as:

Obviously, $_{0}^{C}D^{\alpha}_{t} \varphi(t) \rightarrow D^{\alpha}_{t} \varphi(t)$ whenever $\alpha \rightarrow 1.$

Definition 2.2. The Atangana-Baleanu derivative in the Caputo sense is defined as:\\ $_{a}^{ABC}D^{\alpha}_{t} \varphi(t) = \frac{AB(\alpha)}{1-\alpha} \int_{a}^{t} \varphi^{\prime}(x) \times E_{\alpha}\bigg[-\alpha \frac{(t - x)^{\alpha}}{1- \alpha} \bigg]dx,$ where $0 \leq \alpha \leq 1$, $\varphi$ is in $C[a, b]$ space and $AB(\alpha)$ is the normalization function, such that $AB(0) = AB(1) = 1.$

Definition 2.3. The new fractional derivative associated to integral with the Mittag-Leffler kernel is defined as:

3.   Mathematical model

In this section, the integer order mathematical model is formulated. This mathematical model is formulated by the sub-division of the total cell population of stem cell by $\Delta(t)\; as\; t \rightarrow \infty.$ The proposed model is described by the sub-division of the total cell population given by stem cells $a(t)$, progenitor cells $b(t)$, myeloid cells $c(t)$ and mature cells $d(t)$. The other biological parameters involved in the formulation of the model are $u_{a}, u_{b}, u_{c}$. These three biological parametric cells are the rate of stem, progenitor, and myeloid cells. The division of stem, progenitor, and myeloid cells that regenerate are represented by $v_{a}, v_{b}, v_{c}$. Also, $v_{a}(d), v_{b}(d), v_{c}(d)$ represents the divisional quantities of the stem, progenitor, and myeloid cell dynamical populations, and additionally, $\mu_{a}, \mu_{b}, \mu_{c}, \mu_{d}$ represents the mortality rate of the corresponding cell populations. On division, the flow of stem cells population is

The transition diagram of stem cells dynamics is depicted in Figure 1.

The population of progenitor cells is given by

due to the differentiation of stem cells. Also, the population of progenitor cells is given by

In the same fashion, the biological system class of myeloid cells due to differentiation of progenitor cells is

Similarly, the flow populations of myeloid cells are

Also, in the same fashion, the flow in the class of mature cells due to differentiation of myeloid cells is

The model is given as follow:

with initial conditions:

The comprehensive cells division quantities measure the rule of the two procedures by the mature cell production. Take this for $u_{a}, u_{b}, u_{c} \in [0, 1)$ and $\mu_{a}, \mu_{b}, \mu_{c}, \mu_{d} > 0.$ The division amounts of the stem, progenitor, and myeloid cells are examined by extracellular biological signaling particles. Specifically, the division amounts of stem, progenitor, and myeloid cells are given by

where $e_{a}, e_{b}$ and $e_{c}$ are the uncontrolled divisions of stem, progenitor, myeloid cells, and $f(t)$ is signal strength at time $t$. Suppose that cytokines are discharged by comprehensive cells and that this production is synchronized by mechanisms sensitive to the rate of mature cells,

To estimate the sensitivity to the number of mature cells. This explicit dependency can be rationalized by using a pseudo equilibrium point approximation of the plausible dynamical behavior of the cytokine particles. This estimation throws back the supposition that signal strength attains its efficiency at a lower threshold than the significance of mature cells and mitigates the asymptotic towards zero if the quantity of mature cells increases. The state variables and the biological parameters are shown in Table 1 and Table 2 respectively as below.

Thus, the spitting amounts of stem, progenitor, and myeloid cells depend on the size of mature cells: $v_{a}(d) = \frac{e_{a}}{1 + ld}, v_{b}(d) = \frac{e_{b}}{1 + ld}, v_{c}(d) = \frac{e_{c}}{1 + ld},$ where $e_{a}, e_{b}, e_{c} > 0$.

We explain the total volume of self-regeneration as the chance that a progeny cell remains in the same position of differentiable as the parent cells. We are additionally supposing that mature cells do not increase and thus really proceed into the mature cell mathematical compartment complementary to differentiable stem cells (two-dimensional mathematical compartmental model), progenitor cells (three-dimensional mathematical compartmental model), or myeloid cells (four-dimensional mathematical compartment model).

We present mathematical model equations in dimensionless expressions divided by $\mu_{a},$ and establish the dimensionless time $\tilde{t} : = \mu_{a}t$ :

For the sake of convenience, we are dropping the tilde, Eq $(3.1)$ suit

Fractional-order model

The fractional-order derivative of Atangana-Baleanu in the Caputo sense is extra comprehensive as the kernel elaborated in this operator is both non-singular as well as non-local. Also, this aforesaid operator tends to catch the memory effects as well as the crossover dynamical behavior that persists in the system of biological cells. On taking these points into mind, in this subsection, we again formulate the abcd mathematical model $(3.3)$ with the assistance of the fractional-order operator of Atangana-Baleanu in the Caputo sense, with the order i.e., $\alpha \in [0, 1)$. Therefore, the resulting fractional-order mathematical model demonstrated below is acquired with the substitution of the ordinary derivative by the fractional-order Atangana-Baleanu in the Caputo sense operator:

Thus, in the above proposed mathematical model $(3.4)$, the fractional operator is the Atangana-Baleanu in Caputo sense.

4.   Basic properties of the model

In this section, the well-posedness of the proposed model is discussed below:

4.1.   Existence and uniqueness

In this subsection, we have the powerful feature of the abcd model, said to be existence and uniqueness. For existence and uniqueness purpose we apply the fixed point theory. Initially, the fractional-order model $(3.4)$ is reformulated through the following one-point boundary-value problem:

From the Eq $(4.1)$, the vector demonstrated by $\Theta = (a(t), b(t), c(t), d(t))$ represents the corresponding state variables, and simply the vector $\Theta_{0} = (a(t_{0}), b(t_{0}), c(t_{0}), d(t_{0}))$ demonstrates the corresponding one-point boundary-value conditions. Furthermore, $H$ a continuous vector function as described below:

Also, the function $H$ satisfies the Lipschitz-criteria demonstrated as below:

Moreover, the below theorem gives out the required existence and uniqueness result.

Theorem 4.1. Let $\psi_{1}(\alpha) = \frac{1 - \alpha}{ABC(\alpha)},$ and $\psi_{2}(\alpha) = \frac{\alpha}{ABC(\alpha)\Gamma(\alpha)}.$

If

then the unique solution of the abcd model will exist.

Proof. To take the integral in ABC sense of the Eq $(4.1)$, we have the non-linear Volterra integral equation as follows:

To going next, we interpret $\aleph = (0, T)$ and the correlated function $H : C(\aleph, \mathbb{R}_{+}^{4}) \rightarrow C(\aleph, \mathbb{R}_{+}^{4})$ defined as:

So, the Eq $(4.4)$ takes the following form:

We represent the supremum norm on $\aleph$ by $\|.\|_{\aleph}$ states as follows:

The Banach space is constructed by continuous space $C(\aleph, \mathbb{R}_{+}^{4})$ along with the norm $\|.\|_{\aleph}$. The following inequality can be easily demonstrated:

with $\Theta \in C(\aleph, \mathbb{R}_{+}^{4})$ and $H(t, x) \in C(\aleph, \mathbb{R})$ such that

By the use of Eq $(4.6)$, we have the following result

By using Eqs $(4.2)$ and $(4.8)$ as well as triangular inequality, the Eq $(4.9)$ the following form:

Thus, we obtained

where, $P = \psi_{1}(\alpha)E + \psi_{2}(\alpha)T^{\alpha}_{max}E.$ Hence, if the function $H$ holds the Eq $(4.3)$, this function will be a contradiction. This implies that the Eq $(4.1)$ has a unique solution.

4.2.   Invariant region and attractivity

In this subsection, we show the feasible region in which our solution lies. The feasible region for cell population is $\Theta \subset \mathbb{R}^{4}_{+},$ such that $\Theta = \bigg\{\big(a(t), b(t), c(t), d(t)\big) \in \mathbb{R}^{4}_{+}, a(t) > 0, b(t) \geq 0, c(t) \geq 0, d(t) \geq 0 : a(t) + b(t) + c(t) + d(t) \leq \frac{\Theta}{1 + l d(t)} \bigg\}$.

Theorem 4.2. The biological feasible region $\Theta \subset \mathbb{R}^{4}_{+}$ is positive invariant and attractive for the fractional-order model $(3.4)$ with respect to one-point boundary-value conditions in $\mathbb{R}^{4}_{+}$.

Proof. Let $\Theta = a(t) + b(t) + c(t) + d(t).$ By adding all the four equations on the left hand side of model $(3.4)$, we have

This implies,

As a result of the laplace transformation, we have the following:

Therefore, $\Theta(t) \rightarrow \frac{s_{d}}{1 + ld(t)}$ whenever $t \rightarrow \infty$ for all positive values of time, all the terms given by the closed period with its exact one-point boundary values of the fractional-order abcd model $(3.4)$ factors in $\Theta$ persist in $\Theta$. Hence, the biological feasible region $\Theta$ is said to be positively invariant and will fascinate all the solutions in $\mathbb{R}^{4}_{+}.$

4.3.   Positivity

In this subsection, we explore the positivity of solution of the model $(3.4)$.

Theorem 4.3. For $t > 0$, the solution to the fractional-order abcd model $(3.4)$ is positive as well as bounded for $(a(t_{0}), b(t_{0}), c(t_{0}), d(t_{0})) \in \mathbb{R}^{4}_{+}.$

Proof. We have:

Thus, we have the following feasible region: $\Theta = \{(a(t_{0}), b(t_{0}), c(t_{0}), d(t_{0})) \in \mathbb{R}^{4}_{+} : a(t_{0}), b(t_{0}), c(t_{0}), d(t_{0}) \geq 0\}$. Hence, the sum of all four terms is positive.

5.   The analysis

In this section, we explore the stability of equilibrium points by using the reproduction number of stem cell, progenitor cell and myeloid cell. The reproduction number of stem cell entering in the first compartment mixing the number of stem cell leaving out from the same compartment after division in the life time of stem cell which is denoted by $\mathcal{R}_{a} = (2u_{a} - 1)s_{a}.$ Similarly for progenitor cell, the reproduction number is $\mathcal{R}_{b} = \frac{(2u_{b} - 1)s_{b}}{g_{b}}$ for myeloid cell, the reproduction number is $\mathcal{R}_{c} = \frac{(2u_{c} - 1)s_{c}}{g_{c}}.$

The combination of the basic reproduction numbers $\mathcal{R}_{a}, \mathcal{R}_{b}$ and $\mathcal{R}_{c}$ characterizes the existence and stability analysis of the equilibrium point. For the sake of convenience, we take the following

Existence and stability of equilibrium points

In this subsection, we discuss the existence of equilibrium points, namely free equilibrium points, the absence of stem and progenitor cell equilibrium and endemic equilibrium points.

Theorem 5.1. $(a)$ There exists a free equilibrium $E^{*}_{0} = (0, 0, 0, 0)$ of model $(3.4)$.

$(b)$ There exists a no stem and progenitor cell equilibrium $E^{*}_{1} = (0, 0, c_{1}, d_{1})$ of model $(3.4)$ given by

if

$(c)$ There exists an endemic equilibrium $E_{2} = (a_{2}, b_{2}, c_{2}, d_{2})$ of model $(3.4)$ given by

if

and

Theorem 5.2. The free equilibrium point $E^{*}_{0}$ of the model $(3.4)$ has four real eigenvalues $\mathcal{R}_{a} - 1, g_{b}(\mathcal{R}_{b} - 1), g_{c}(\mathcal{R}_{c} - 1)$ and $-g_{d}$ and is locally asymptotically stable if

and unstable if

Proof. At $E^{*}_{0} = (0, 0, 0, 0)$, the variational matrix of the model $(3.4)$ is described as below:

From Eq $(5.9)$ we observed that the eigenvalues of the variational matrix at free equilibrium points are $\lambda_{1} = \mathcal{R}_{a} - 1, \lambda_{2} = g_{b}(\mathcal{R}_{b} - 1), \lambda_{3} = g_{c}(\mathcal{R}_{c} - 1)$ and $\lambda_{4} = -g_{d}$ and is locally asymptotically stable if

and unstable if

Now, we will discuss the stability analysis of the absence of stem and progenitor cell equilibrium points. Suppose the condition $(5.4)$ holds. We introduced two positive constants for the sake of eigenvalues connected with the absence of stem and progenitor cell equilibrium points.

Theorem 5.3. The absence of stem and progenitor cell equilibrium $E^{*}_{1}$ of model $(3.4)$ has four eigen values $\lambda_{i}, i = 1, 2, 3, 4$ where

and $\lambda_{3, 4}$ lies on the left side of the half plane. The absence of stem and progenitor cell steady point $E^{*}_{1}$ is locally asymptotically stable if

and unstable if

Proof. Since, we have $f_{c, 1}(d_{1}) - g_{c} = 0\; and\; c_{1} = \frac{g_{d}d_{1}}{f_{c, 2}(d_{1})}.$ By dropping the sign from $c_{1}$ and $d_{1}$, we have the following characteristic equation:

We obtain roots,

from Eq $(5.2)$ and the identity

It follows $E^{*}_{1}$ is locally asymptotically stable if $\mathcal{R}_{b} > \mathcal{R}_{a}\; and\; \mathcal{R}_{c} > \mathcal{R}_{b}$ and unstable if $\mathcal{R}_{b} < \mathcal{R}_{a}\; and\; \mathcal{R}_{c} < \mathcal{R}_{b}.$

Now finally we will look into the stability of the endemic equilibrium points. Suppose that the condition $(5.5)$ holds, and we define

Then,

Theorem 5.4. If $\mathcal{R}_{a} > 1$ holds, then $E_{2}$ is locally asymptotically stable if

and is unstable if

Proof. We have $f_{a, 1}(d_{2}) = 1, b_{2} = \frac{g_{d}d_{2}}{f_{b, 2}(d_{2})}\; and\; a_{2} = \frac{(g_{b} - f_{b, 1}(d_{2}))g_{d}d_{2}}{f_{a, 2}(d_{2})f_{b, 2}(d_{2})},$ by omitting the index from $a_{2}, b_{2}$ and $d_{2}$, we have the characteristic equation:

We have the following expressions:

Therefore, Eq $(5.21)$ becomes

By using the $g_{b}(\mathcal{R}_{b}/\mathcal{R}_{a}) = \rho(\epsilon)$, in Eq $(5.22)$, we have following characteristic equation

where,

By using the Routh-Hurwitz criteria ^[12] for Eq $(5.23)$, all the roots have negative real parts if and only if

6.   Numerical solution

In this section, we employed the iterative approach of the fractional-order cell division model $(3.4)$ with Atangana-Baleanu in the Caputo sense to obtain the graphical results of abundant biological parameters and memory index on cell division. The well-organized numerical technique ^[17] is utilized for this determination. First of all, a brief derivation of the approximate approach is introduced for Eq $(4.1)$ and then implemented on the system for $(3.4)$ in compensation.

Using definition $2.3$ on both sides of the initial condition $(4.1)$. Then we obtain the integral equation as below:

where, $\psi_{1}(\alpha) = \frac{1 - \alpha}{ABC(\alpha)}$ and $\psi_{2}(\alpha) = \frac{\alpha}{ABC(\alpha)\Gamma(\alpha)}.$

Moreover, at $t = t_{i+1} = (i+1)h$, we have

In addition to the support of a polynomial interpolation scheme, the operator $H(y, z(y))$ is approximated by employing a two-step Lagrange polynomial on the closed interval $[t_{q}, t_{q+1}]$. So, we started as follows:

The Eq $(6.2)$ becomes:

Furthermore, for the simplification of the Eq $(6.4)$, the approximate solution is considered as follows:

where,

$p_{1} = i - q, p_{2} = i + 1 - q, p_{3} = i + 2 - q + \alpha, p_{4} = i + 2 - q + 2 \alpha,$

$z_{1} = H_{1}(t_{i-1}, z(t_{i-1})), z_{2} = H_{2}(t_{i-1}, z(t_{i-1})), z_{3} = H_{3}(t_{i-1}, z(t_{i-1})), z_{4} = H_{4}(t_{i-1}, z(t_{i-1})).$

7.   Numerical results

In this section, we simulate the mathematical model $(3.4)$ to validate the analytic results. The default parameter are given in Table 2. Figures 2–4 show the variation of all three types of cell population versus time for different of their self-regeneration rates. Figure 2 displays the variation of stem cell population $a(t)$ with respect to time for different values of self-regeneration rate of stem cell. It is seen that as we increase the value of $u_{a}$ from $0.52, 0.54, 0.56$, the population of stem cell increases but it decreases with time. Figure 3 also displays the variation of progenitor cell $b(t)$ with respect to time for different values of self-regeneration rate of progenitor cell. It is seen that as we increases the value of $u_{b}$ from $0.90, 0.92, 0.94$, the population of progenitor cell increases but it also decreases with time. Figure 4 shows the effects of self-generation rate on the dynamics of myeloid cell $c(t)$. It is clear from the figure that the value of self-generation significantly affects the myeloid cell population.

Figure 5 demonstrates the effects of division rate $v_{a} = 0.1, 0.2, 0.3$ on the dynamics of stem cell. It is seen that the population of stem cells decreases with the increase in the values of the division of stem cells. The effects of the self-regeneration rate of the stem cell on the population are seen. Figure 6 shows the effects of the division rate on the population of progenitor cells. It is noticed that as we increase the value of division rate $v_{b}$ from $0.5, 0.7, 0.9$, the population of progenitor cells decreases. It suddenly attains maximum value at $t = 2$ and afterward, it remains almost constant. Figure 7 depicts the variation of myeloid cell with respect to time for different values of division rate of myeloid cell.

Figures 8–11 exhibit the effect of $\alpha$ on all four cell populations. In Figure 8, it is seen that by varying the $\alpha$ from $0.3, 0.5, 0.7$ the stem cell $a(t)$ is affected significantly and it increases with time but decreases with increases in $\alpha$. We can also find similar trends in Figures 9 and 10, but in Figure 11, the population of $d(t)$ sharply falls and decreasing with time and become almost zero.

8.   Conclusions

In this study, a fractional-order mathematical model is analyzed to understand the dynamical behavior of the stem cell maturation that is responsible for self-renewal, cell death, and regulation of cell division. We considered four compartmental models for a different kind of cell population differentiated for stem cells. The concept of reproduction number is introduced that facilitates our analysis. This reproduction number quality is a threshold parameter that is used to characterize the existence of equilibria and decide the stability of the equilibrium points. It is found that the equilibrium points are stable when $max\{\mathcal{R}_{a}, \mathcal{R}_{b}, \mathcal{R}_{c}\} < 1$ and unstable when $max\{\mathcal{R}_{a}, \mathcal{R}_{b}, \mathcal{R}_{c}\} > 1$. Our results indicate that the rate of stem cell division and stem cell self-renewal directly depends upon the environment heterogeneity. In further, these results could be verified by using real data from realistic sources to validate the proposed model.

Acknowledgement

The authors would like to acknowledge the support from Taif University researchers supporting project (no. TURSP-2020/10), Taif University, Taif, Saudi Arabia.

Conflict of interest

The authors declare that there are no conflict of interest.