Stability analysis of a multiscale model including cell-cycle dynamics and populations of quiescent and proliferating cells

1.   Introduction

Mammalian cell division patterns are critical to understanding human tumor progression. It has drawn the attention of many researchers, and it has been the topic of intense research for a long time. Several research works utilize age-structured frameworks to investigate the cell cycle. Some examples of age-structured growth models include epidemic ^[1,2,3], microscopic virus ^[4,5] and cell population ^[6,7,8,9] models. The concealed molecular complexity of a tissue, on the other hand, demand a more thorough modeling framework that includes special cellular and molecular interactions.

Cells that divide in living tissues can be divided into two categories: proliferating and quiescent. The proliferating cells go through different phases in the cell cycle ($G_1, S, G_2, M$) while dividing. Quiescent cells, however, do not grow or divide, but rather move to the $G_0$ phase, where they remain until they differentiate or undergo apoptosis. To preserve tissue homeostasis, cells must be capable to transition between proliferative and quiescent phases. The transition between proliferating and quiescent compartments, however, is dependent on signaling molecules called anti-growth or growth factors ^[10]. In a population of tumor cells, proliferating cells multiply till the tumor becomes active and aggressive. Additionally, to maintain homeostasis, the total number of cells in all sub-populations stays at equilibrium; hence, in a healthy and mutated cell populations, the proliferative compartment is constrained in size. Schematics in Figure 1 depicts multiscale-modeling framework used in this paper.

Mutations in genes that regulate cell growth and division can disturb the finely-tuned equilibrium that governs cell proliferation, which in turn can result in the emergence of cancer. Although other factors like environmental exposures and lifestyle choices may also play a role in cancer development, mutations in genes are a major contributor to this disease. Mutations can be of several types, those that are particularly significant for cancer involve increased potential for proliferation, decreased apoptosis, genetic instability, and reduced tumor suppression ^[11]. Furthermore, studies have shown that the transformation of a normal cell into a cancerous one usually requires the accumulation of one to ten mutations ^[11,12]. The mutated cell population also consists of its own progeny including different differentiation stages and quiescent cells.

This research focuses on the model development of the cell population in all healthy and mutated cell populations. We investigate the coupling dynamics of tissue cell density and cell cycle proteins. Previous studies have used age-structured models to investigate cell populations in quiescent phase ^[6] only, the proliferating phase ^[13,14] only, or both phases together ^{[7,15,16,17,18,19]}. Despite this, the impact of molecular interactions on the interplay between proliferative and quiescent phases at the subcellular level has yet to be examined. Thus, the primary aim of this paper is to develop a multiscale model utilizing mathematical techniques that can capture the intricacies of a complex system existing at the sub-cellular level.

The emphasis is laid upon the interaction between the macro- (population dynamics) and micro-scale (cell cycle dynamics) at two distinct levels. The concept of age is used to refer to the time since the last division, as described in previous studies ^[14,20]. In age-structured models, an additional variable, referred to as "age" ($a$), is introduced, which has a physiological rather than physical interpretation. The idea of cell age reflects the variability (biological) within a population of dividing cells and is distinct from the physical time variable ($t$). To model the behavior of cell populations in both the proliferative and quiescent stages at the macroscale, we use partial differential equations (PDEs). For predicting sub-cellular protein interactions related to cell cycle dynamics, we employ ordinary differential equations (ODEs). The two scales are connected via the feedback incorporated in both directions. Within a cycle of cell division ($G_1$, $S$, $G_2$, $M$), cells in the early proliferating phase ($\mathrm{G}_1$) can move to the quiescent-phase until the restriction point ($\mathrm{R}$).

It is clear that the restriction point ($\mathrm{R}$) cut the $\mathrm{G}_1$ phase into two portions, where before $\mathrm{R}$, the cells become quiescent but once $\mathrm{R}$ is passed, cells can no longer avoid division ^[21,22]. During the quiescent phase, cells do not grow or divide, but they still complete other cellular functions. The bidirectional transition of cells between quiescent and proliferation phases in both healthy and mutated cell populations plays a crucial role in maintaining tissue homeostasis, and it is controlled by the conditions of extracellular environmental ^[10]. In tumors, the transitional balance between two compartments is disturbed, and cells can grow uncontrollably ^[23]. Recent experimental results support that cyclin proteins are the most important to regulate any change in the cell-cycle phase ^[24]. Therefore, we utilize a vital aspect of cell cycle dynamics (i.e., the $G_1-S$ phase-transition) to predict how the transitional balance between proliferating and quiescent subpopulations evolve, which is essential for maintaining homeostasis.

Several proteins are involved in regulating the transitions between different phases of the cell cycle at the microscale. The interactions between these proteins have been mathematically modeled and simulated using ODEs by various researchers, as documented in ^{[25,26,27,28,29,30]} and related references. However, to keep things simple, we focus only on four proteins: Cyclin $\mathsf{D-CDK 4/6}$, $\mathsf{p21}$, $\mathsf{E2F}$, and $\mathsf{Rb}$. These proteins are mainly involved in Cyclin $\mathsf{D}$ activity and the transition of cells from $G_1$ phase to $S$ phase. Experimental data supports this choice, as Cyclin $\mathsf{D}$ has been shown to regulate the transition between the $G_0$ and $G_1$ phases ^[31,32,33]. In addition, over-expression of Cyclin $\mathsf{D}$ leads to cell division commitment in the proliferative phase, while under-expression leads to the quiescent phase. It is essential to note that these molecular-interactions are considered to occur in a rapidly growing cell population rather than a single cell. Furthermore, we assume the average protein concentrations in both quiescent and proliferating cell subpopulations, thus omitting the variability between individual cells. The proteins involved in the cell cycle play a vital role in regulating its progression. Cyclin proteins, which are structural proteins, along with their inhibitors cyclin-dependent kinases ($\mathsf{CDK}$), control the various phases of the cell cycle. Each phase of the cycle is regulated by a specific Cyclin/$\mathsf{CDK}$ complex.

The initiation of a cell cycle is triggered when the microenvironment of a cell receives sufficient growth signals ^[34], which induce the activity of Cyclin-$\mathsf{CDK}$ 4-6 complex. During the $G_1$ phase, Cyclin $\mathsf{D}$ is activated only by growth signals, and in case of low or no growth factors, its concentration decreases, and the cell doesn't enter the cycle. Growth factors bind to specific receptors which are situated on the cell's external cytoplasmic membrane, activating intracellular signaling pathways (such as Ras/Map/Rap kinase), eventually which directing to the production of Cyclin $\mathsf{D}$ (see ^[35,36,37], for more details). The active complex of Cyclin $\mathsf{D}$ and $\mathsf{CDK 4/6}$ is synthesized at maximum. Cyclin complex then triggers the transcription factor $\mathsf{E2F}$'s activation by phosphorylating retinoblastoma protein $\mathsf{Rb}$. Consequently, the transcription factor $\mathsf{E2F}$ accumulates and triggers some other important cyclins. It is important to note that these interactions occur in a fast growing population of cells rather than in a single cell.

In summary, our contribution is to develop a multiscale model that describes the coupling between two predominant scales and focuses on the challenges related to the disruption of cell transitioning between proliferating and quiescent states, leading to uncontrolled tumor growth. Specifically, we investigate whether Cyclin $\mathsf{D}$ complex is one of the important cause in creating a deregulation in cell transitioning between proliferating and quiescent cells.

The structure of rest of this paper is as follows. In Section 2, we delve into the details of multiscale mathematical modeling of proliferating and quiescent cells with regards to the dynamics of cell cycle. In Section 3, we investigate the existence and uniqueness of non-negative solutions utilizing spectral and semigroup theory. Moving on to Section 4, we start with deriving steady-state solutions and then establish spectral criteria for their local stability. Specifically, we show that if the linearised semigroup's growth bound is negative, the steady-state solution is locally asymptotically stable, while a positive growth bound implies instability of the steady-state solution. Lastly, Sections 5 and 6 offer results discussion and the concluding remarks of this paper, respectively.

2.   Mathematical modeling

2.1.   Age-structured model

The cell densities of healthy and mutated cells in the proliferative and quiescent compartments are described by nonlinear hyperbolic transport PDEs that relate the cell density distribution to both physiological age $a$ and time $t$. Specifically, the densities of healthy cells in the proliferative and quiescent phases are expressed as follows:

where the rate evolution of a cell cycle is denoted by $g_h(a)$ in the equation. The first term on the right side, $\gamma_h(N)q_h(a, t)$, represents the transition to proliferating from quiescent cells, while the term $\tau_h(a)p_h(a, t)$ represents the cell densities for completing cell division in some age of the proliferating phase. Cells which move to the quiescent phase without undergoing division are represented by $\alpha_h(a, x_1)p_h(a, t)$. The loss in the proliferating cells due to apoptosis/necrosis is represented by the death rate $\mu_{p_h}(a)$. The inflow from healthy proliferating cells, regulated by the microscale variable of Cyclin complex concentration $x_1$ for each age, is denoted by the first term in Eq (2.2): $\alpha_h(a, x_1)p_h(a, t)$. The next term depicts a loss in the quiescent cells due to either by returning to proliferating phase at the rate $\gamma_h(N)$ or by cell death due to apoptosis (or necrosis), as represented by the death rate $\mu_{q_h}(a)$.

Next, the cell density of mutated cells in mutated proliferating $(p_c)$ and quiescent $(q_c)$ phases, respectively, is presented:

where the terms used are similar to those in the case of healthier cells, as shown in Eqs (2.1) and (2.2). The total cell number in both healthy and mutated populations of cells in quiescent and proliferating phases is denoted by $N(t)$ and is defined in Eq (2.5). In the case of quiescent cells, aging does not occur (i.e., the cells stop aging), so the convection term related to physiological age $a$ is absent in Eqs (2.2) and (2.4). The total number of cells, represented by $N(t)$, represents the sum of all cells in the proliferating and quiescent phases throughout all ages, and can be expressed as

where maximum age of the cells is given by $a^\star$. The initial conditions are given below:

The boundary conditions are given as follows:

for $t > 0$, where the number 2 shows the two newborn cells initializing in the proliferating phase, and the parameter $m$ represents the mutation rate. Since healthy cell can acquire a mutation only during a division process, therefore, new born mutated cells start will start at age 0.

The function $\tau_i(a)$ represents the cell number that finish dividing at a particular age in both healthy and mutated proliferating phases. Here, the index $i$ indicates whether the compartment is healthy or cancerous, denoted by $h$ and $c$, respectively. The function $\tau_i(a)$ is regulated by the age of the cell, denoted by $a$, and is almost zero until a minimum cell age. Afterward, the function increases until it reaches the age of $a^*$.

The maximum proliferation rate is represented by $\rho_{1, i}$, while $\rho_{2, i}$ is the age t achieve half-maximum. The Hill coefficient is represented by the exponent $\gamma_{1, i}$.

Next, we establish the rate at which cells transition to the quiescent phase from the proliferating compartments, which depends on both the age of the cell ($a$) and the quantity of the Cyclin $\mathsf{D-CDK 4/6}$ complex ($x_1$).

The function $\alpha_i(a, x_1)$ depicts the number non-dividing cells due to anti-growth factors. The age-dependence of $\alpha_i$ is motivated by the fact that cells transition to the quiescent phase from the proliferating phase only until they reach a specific age that marks a restriction point ($\mathrm{R}$) in the cell cycle (which is also $G_1-S$ phase transition). However, before the restriction point, the Cyclin complex's concentration $x_1$ must be below a certain value to enable cells to exit the proliferating phase. In Eq (2.10), the Hill coefficients are represented by $\gamma_{2, i}$ and $\gamma_{3, i}$, while $\sigma_{2, i}$ and $\sigma_{3, i}$ denote the concentration of the Cyclin $\mathsf{D-CDK 4/6}$ complex $x_1$ and the age $a$, respectively. After $\gamma_{2, i}$ and $\gamma_{3, i}$, the rate function $\alpha$ decreases asymptotically to zero, preventing cells from transitioning to the quiescent phase. This implies that at age $\sigma_{3, i}$, cells are inevitably committed to entering the proliferation phase. Finally, $\sigma_{3, i}$ represents the threshold concentration of the Cyclin complex to determine the restriction point $\mathrm{R}$.

The function $\gamma_i(N)$, which determines the number of cells transitioning to the proliferating phase from the quiescent phase, is represented by a Hill function of $N$ that decreases monotonically:

where the Hill function are defined as follows: $\nu_i$ specifies the maximum rate at which cells transition to proliferating from quiescent population, when there are no cells, i.e., $N = 0$; $\kappa_i$ is the Hill coefficient, and $\theta_i$ represents the proportion of the total cell population that reaches half the maximum value of $\nu_i$. This implies that the number of quiescent cells transitioning to the proliferative compartment decreases to zero as the cell population increases, illustrating density inhibition.

Cell growth is controlled by proteins such as cytokines and other factors that regulate proliferation ^[38]. Cytokines bind to specific receptors, activating signaling pathways ^[39]. The cytokine signals that regulate cell numbers are reliant on the total population of cells, as demonstrated by various studies ^[40]. For a detailed explanation of cytokine signal dynamics, refer to ^[41,42]. Using the quasi-steady-state approximation, we can express the quantity of growth factors ($g_f$) produced by the entire cell population ($N$) as

2.2.   Cell cycle model

We restrict ourselves to only four microscale proteins (states of the model) in our cell cycle model, as mentioned earlier, that are sufficient to account for the reversible transitions between the quiescent and proliferating phases.We utilize Michaelis-Menten kinetics to depict the chemical reactions occurring between enzymes and substrates during the cell cycle. The basics of this approach are briefly outlined below. When adequate growth factors are present, Cyclin $\mathsf{D}$ protein combines with its catalytic partner CDK4-6, resulting in the formation of Cyclin $\mathsf{D-CDK 4/6}$ complex. As outlined in ^[31,43], this complex is responsible for phosphorylating various proteins involved in the cell cycle. These phosphorylated proteins are essential for advancing through the initial growth phase of the cell cycle and crossing the restriction point $\mathrm{R}$. To provide further details, the Cyclin $\mathsf{D-CDK 4/6}$ complex induces phosphorylation of the retinoblastoma protein $\mathsf{Rb}$, resulting in its deactivation and subsequent release of the transcription factor $\mathsf{E2F}$. This activates several signals that promote cell growth, facilitating the progression of the cell cycle. In contrast, $\mathsf{p21}$ regulates the cell cycle by inhibiting the functions of various $\mathsf{CDK}$ proteins. These proteins are described in Table 1.

In this study, we examine the behavior of a single cell to represent the dynamics of all cells within a population. Assuming that all cells behave similarly, we utilize a single ordinary differential equation (ODE) model with similar parameters for all cells in the population to describe the underlying cell cycle dynamics on a microscale. We also account for cells with shorter cycles on a macroscale using the function $\tau_i(a)$. Our model considers a specific age $a^\star$ at which the representative cell completes division. The cell cycle dynamics are described by the ODE system presented in ^[44].

A detailed description of all the terms and parameters involved in the model equations (2.13)–(2.16) can be found in ^[45]. Although we will not delve into the full derivation of these equations here, inquisitive readers can refer to ^[44] for a comprehensive explanation. To aid in understanding, we have included simulations of the four microscale states mentioned earlier in Figure 2.

3.   Existence and uniqueness of non-negative solution

This section presents the uniqueness of the solution to the initial-boundary value problem (2.1)–(2.7) and (2.13)–(2.16), which we will simplify by using the microscale model for the entire time $t$, instead only until age $a$. We introduce Banach spaces, $X = L^1(0, a^\star)\times L^1(0, a^\star)\times L^1(0, a^\star)\times L^1(0, a^\star)$ and $Y = L^1(0, a^\star)\times L^1(0, a^\star)\times L^1(0, a^\star)\times L^1(0, a^\star)$, with the norm $|\phi| = \sum_{i = 1}^4|\phi_i|1$ for $\phi(a) = (\phi_1(a), \phi_2(a), \phi_3(a), \phi_4(a))^\mathrm{T}\in X$ and $|\varphi| = \sum{i = 1}^4|\varphi_i|1$ for $\varphi(a) = (\varphi_1(a), \varphi_2(a), \varphi_3(a), \varphi_4(a))^\mathrm{T}\in Y$, where $|\cdot|1$ is the standard norm of $L^1(0, a^\star)$. We first treat the initial-boundary value problem of system (2.1)–(2.7) as an abstract Cauchy problem on Banach space $X$. We assume that $g{h_a}, g{h_{aa}}, g_{c_a}, g_{c_{aa}}\in L^{\infty}((0, a^\star)\times\mathbb{R}^+)$, and non-negative death rates, that is, $\mu_{p_h}(\cdot) = \mu_{q_h}(\cdot)\geq0$ and $\mu_{p_c}(\cdot) = \mu_{q_c}(\cdot)\geq0$, and are locally integrable on $[0, a^\star)$. The transition rate $\alpha_i(a, x_1)\in L^{\infty}((0, a^\star)\times(0, a^\star))$, and $\tau_i(a)\in L^{1}(0, a^\star)$. We start by defining a linear operator $\mathsf{A}_1$ as follows:

The symbol $\mathrm{T}$ denotes the transpose of the vector, and the domain $\mathrm{D}(\mathsf{A}_1)$ is given by the following:

The nonlinear operator $\mathsf{F}_1:X \times Y\to X$ is given by

where the linear operator $N$ on $L^1(0, a^\star)\times L^1(0, a^\star)\times L^1(0, a^\star)\times L^1(0, a^\star)$ is given by

Consider $\upsilon(t) = (p_h(\cdot, t), q_h(\cdot, t), p_c(\cdot, t), q_c(\cdot, t))^\mathrm{T}\in X$. We can define the initial-boundary value problem (2.1)–(2.7) as an abstract semilinear initial value problem (IVP) in $X$, as shown below:

where $\upsilon_0(a) = (p_{h_0}(a), q_{h_0}(a), p_{c_0}(a), q_{c_0}(a))$.

Next, we define IVP (2.13) as a Cauchy problem on the Banach space $Y$. Suppose $\mathsf{A}_2$ is a linear operator which reads

where the domain $\mathrm{D}(\mathsf{A}_2)$ is

We define the nonlinear operator $\mathsf{F}_2:X \times Y\to Y$ by

where $\phi\in X, \, \varphi\in Y$. Take $v(t) = (x_1(t), x_2(t), x_3(t), x_4(t))^\mathrm{T}\in Y$. Then (2.13)–(2.16) can be expressed as an abstract semilinear IVP in $Y$:

where $v_0(t) = (x_1^0, x_2^0, x_3^0, x_4^0)$, we can now establish a joint Cauchy problem for (3.1) and (3.2) as shown below:

where $\zeta = (\upsilon, v)$, $\zeta_0 = (\upsilon_0, v_0)$, $\mathsf{A} = \left(\begin{array}{cc} \mathsf{A}_1 & 0\\ 0 & \mathsf{A}_2 \end{array}\right)$, $\mathsf{F} = \left(\begin{array}{c} \mathsf{F}_1 \\ \mathsf{F}_2 \end{array}\right)$ and the Banach space is $Z = \{X, Y\}$. Assuming that $T(t)$ is a $C_0$-semigroup generated by $\mathsf{A}$ for $t\geq0$, and the operator $\mathsf{F}$ is continuously Fréchet differentiable on $Z$ (specifically, both $\mathsf{F}_1$ and $\mathsf{F}_2$ are Fréchet differentiable on $X$ and $Y$), a continuous mild solution $t\to \zeta(t, \zeta_0)$ exists and is unique for each $\zeta_0\in Z$ on a maximal interval $[0, t_1)$ in $Z$.

where $t_1$ can be either $+\infty$ or $\lim_{t\to t_1^-}|\zeta(t, \zeta_0)| = \infty$. Moreover, if $\zeta_0\in \mathrm{D}(\mathsf{A})$, then $\zeta(t, \zeta_0)\in \mathrm{D}(\mathsf{A})$ for $0\leq t < t_1$, and the function $\zeta\to \zeta(t, \zeta_0)$ is continuously differentiable and satisfies (3.3) on $[0, t_1)$. This has been established in Proposition 4.16 ^[46,47].

Remark 3.1. Let's take $p_{h, \max}$, $q_{h, \max}$, $p_{c, \max}$, $q_{c, \max}$, $x_{1, \max}$, $x_{2, \max}$, $x_{3, \max}$ and $x_{4, \max}$ to represent the maximum values of the solution variables. If we normalise the governing equations using $N(a) = p_h(a, t)+q_h(a, t)+p_c(a, t)+q_c(a, t)+x_1(a)+x_2(a)+x_3(a)+x_4(a)$, then an a-priori estimate on these would lead to $p_h(a, t)+q_h(a, t)+p_c(a, t)+q_c(a, t)+x_1(a)+x_2(a)+x_3(a)+x_4(a) = 1$.

Lemma 3.1. Let $\Omega = \{(p_h, q_h, p_c, q_c, x_1, x_2, x_3, x_4)\in Z|p_h\geq0, q_h\geq0, p_c\geq0, q_c\geq0, x_1\geq0, x_2\geq0, x_3\geq0, x_4\geq0\}$ and let $\Omega_0 = \{(p_h, q_h, p_c, q_c, x_1, x_2, x_3, x_4)\in$ $Z|0\leq p_h\leq p_{h, \max}, 0\leq q_h\leq q_{h, \max}, 0\leq p_c\leq p_{c, \max}, 0\leq q_c\leq q_{c, \max},$ $0\leq x_1\leq x_{1, \max}, 0\leq x_2\leq x_{2, \max}, 0\leq x_3\leq x_{3, \max}, 0\leq x_4\leq x_{4, \max}\}$. Then after a finite time, the mild solution $\zeta(t, \zeta_0)$ of (3.3), where $\zeta_0\in \Omega$, enters a positively invariant set $\Omega_0$.

Proof. To obtain the solution of Eq (2.1), we will begin by utilizing transformations $\tilde{p}_h(a, t) = g_h(a)p_h(a, t)$ and $\tilde{q}_h(a, t) = g_h(a)q_h(a, t)$ for $t\in[0, t_1]$ and $a\in[a_0, a^\star)$. Then, for $t\in(0, t_1)$ and $a\in(a_0, a^\star)$, we have from Eq (2.1)

Next, we apply the parameter transform to remove the term $g_h(a)$ and define a new age variable $\eta$ for both $p_h$ and $q_h$, Lemma 3.1 ^[41]. This yields the expression:

Therefore, from Eq (3.5), it follows that

To obtain the explicit relation of $\tilde{p}_h(a(\eta), t)$, we will utilize the method of characteristics (MOC). Specifically, we assume that $\tilde{p}_h(a(\eta), t)$ is governed by an ODE along the curve $(a(\psi_1(y)), \psi_2(y)) = \psi(y)$, and therefore, we have

where $c_1, c_2\in\mathbb{R}$ are constants. Then, it follows

We can now write $\tilde{p}_h$ using an ODE (3.7) so that

Now, we define the boundary set $\Gamma$ as ${[a_0, a^\star)\times{0}}\cup{{0}\times[0, t_1]}$, which enables us to use the boundary condition to determine $\tilde{p}_h(a(c_1), c_2)$ if a curve $(a(\psi_1(y)), \psi_2(y))$ begins in $\Gamma$. In order for $(a(y+c_1), y+c_2)$ to lie on $\Gamma$, either $c_1 = 0$ or $c_2 = 0$. Therefore, we have the following two scenarios:

In the first scenario, we can randomly choose $c_1 = 0$ and $c_2\in[0, t_1)$. In this case, we have

The solution in ${(a(\eta), t)|t\in[0, t_1], \eta\in[0, \min(\eta^*, t))}$ can be obtained using the characteristic solution as follows:

which implies

Using the above equation, we can obtain the expression for $g_h(a(\eta))p_h(a(\eta), t)$ when $\eta < t$. By choosing an arbitrary $c_1\in[0, \eta^*)$ and setting $c_2 = 0$, we obtain

Using the characteristic solution, we can obtain a solution in the set ${(a(\eta), t)|t\in[0, t_1], \eta\in[t, \eta^*)}$ as follows:

which results into

Hence, the relation for $g_h(a(\eta))p_h(a(\eta), t)$ is now established for $\eta > t$. As a result, the ultimate solution for $g_h(a(\eta))p_h(a(\eta), t)$ can be expressed as

where the boundary condition $\tilde{p}_h(a(0), t-\eta)$ is denoted as $h(t-\eta)$. Note that for positive initial data, the above expression is positive and for $g_h(a)q_h(a, t)\geq0$.

We then derive the solution expression from (2.2) as shown below:

As a direct consequence, we observe that $q_h(a, t)$ is non-negative for positive initial data and whenever $g_h(a)q_h(a, t)\geq0$. Similarly, we can obtain the solution expression for $p_c$ and $q_c$.

Next, to ensure the positivity of the coupled ODE model (2.13), we express the system of ODEs as follows:

where $\mathsf{F}_1$, $\mathsf{F}_2$, $\mathsf{F}_3$ and $\mathsf{F}_4$ correspond to the vector fields of the microscale states $x_1$–$x_4$. It is worth noting that in (3.9), $\mathsf{F}_1$ does not depend on $N$ (i.e., $p_h$, $q_h$, $p_c$ and $q_c$), as $N$ changes with time and is a fixed constant at each time step, which determines the growth factors entire age range.

To ensure that the solutions of all ODEs is positive, it is essential to verify that the vector fields $\mathsf{F}_1, \mathsf{F}_2, \mathsf{F}_3, \mathsf{F}_4$ are smoothly differentiable and oriented in a direction that points away from the negative regions in the state space. Starting with the ODE for $x_1$ from (3.9), we set $x_4 = 0$ in $\mathsf{F}_1(x_1, x_4)$ to obtain $\dot{x}_1 = \mathsf{F}_1(x_1)$. We can observe that $\mathsf{F}_1(x_1) = k_{1s}\bigg(\frac{g_f}{k_{gf}+g_f}\bigg)-k_{1d}\bigg(\frac{x_1}{k_{1}+x_1}\bigg) > 0$ for all $a > 0$ when $k_{1s}\bigg(\frac{g_f}{k_{gf}+g_f}\bigg) > k_{1d}\bigg(\frac{x_1}{k_{1}+x_1}\bigg)$. This implies that the concentration of $x_1$ consistently rises more than it falls over time, as the sole source of an increase in $x_1$ concentration is from growth factors. Consequently, during periods when growth factors are at their minimum, the concentration of $x_1$ is also at its minimum, meaning the amount of degradation or decrement cannot surpass the activation of the $x_1$ complex. Given that the solution to system (2.13)–(2.16) is unique for each initial condition, as can be observed from (3.3) and (3.4), we can infer that the solution remains in the first quadrant for any $x_4 > 0$. Therefore, the positivity of the solution for $x_1$ is guaranteed. To obtain an ODE for $\dot{x}_2$, we assume $x_1 = 0$ in $\mathsf{F}_2(x_1, x_2, x_3)$. This results in $\dot{x}_2 = \mathsf{F}_2(x_2, x_3)$, whose solution takes the form $x_2(a) = x_2^0e^{-(k{32}x_3(a)-k{2d})a}$. Hence, for any positive initial data, $x_2(a)$ remains positive for all ages and $x_3(a)$ values. Similarly, we substitute $x_3 = 0$ in $\mathsf{F}_2(x_1, x_2, x_3)$ to obtain a nonlinear ODE $\dot{x}_2 = \mathsf{F}_2(x_1, x_2)$ for $x_2$. Although an explicit solution cannot be computed, the phase portrait of $(x_1, x_2)$ reveals that the solution trajectories move away from the axis separating the positive and negative space for positive initial data. By following a similar procedure, we can establish sufficient conditions for the positivity of solutions for $x_3(a)$ and $x_4(a)$. Therefore, we conclude that if $\zeta_0\in\Omega$, then $\zeta(t, \zeta_0)\in\Omega$ for all $t > 0$.

The above analysis implies that the local solution $\zeta(t, \zeta_0)$ of (3.3) with initial conditions $\zeta_0\in \mathrm{D}(\mathsf{A})\cap \Omega$ has a well-defined and finite norm. Consequently, we obtain our final result.

Theorem 3.1. The abstract Cauchy problem (3.3) has a unique global classical solution on $Z$ with respect to the initial data $z_0\in\Omega\cap D(\mathsf{A})$.

As a result of having positive initial data, the IVP (2.1)–(2.4) possesses a singular positive solution.

4.   Existence and stability of steady-state

This section aims to determine the steady-state solution of the model and to present sufficient conditions for the existence of a positive steady-state. To this end, we specify some notation. Let $X$ be a real or complex Banach space, and let $X^\star$ denote its dual space. We denote the value of $F\in X^\star$ at $\psi\in X$ as $\langle F, \psi\rangle$. Additionally, we define a cone $X_+$ as a non-zero set that satisfies $X_+\cap(-X_+) = {0}$, $\lambda X_+\subset X_+$ for $\lambda\geq0$, and $X_++X_+\subset X_+$. Furthermore, we define the dual cone, denoted as $X^\star_+$, as the subset of the dual space.

4.1.   Existence of steady-states

The steady-states of the system (2.1)–(2.4) and (2.13)–(2.16) are denoted by $\bar{p}_h(a)$, $\bar{q}_h(a)$, $\bar{p}_c(a)$, $\bar{q}_c(a)$ and $\bar{x}_1-\bar{x}_4$. These steady-states must satisfy the following set of time-invariant ordinary differential equations:

where $\bar{\gamma}_i = \gamma_i(\bar{N})$, $\bar{g}_f = g_f(\bar{N})$ and $\bar{N} = \int_0^{a^\star}(\bar{p}_h(a)+\bar{q}_h(a)+\bar{p}_c(a)+\bar{q}_c(a)) \mathrm{d} a$. Since the cell cycle model's ODEs are age-dependent and the system is in a steady-state due to the input of growth factors, all cell cycle states attain a steady-state. As a result, we can determine the steady-states of the quiescent and proliferating cell populations, represented by $\bar{p}_h(a)$, $\bar{q}_h(a)$, $\bar{p}_c(a)$ and $\bar{q}_c(a)$, without explicitly solving the equations of microscale model. Solving the system (4.10) for $\bar{p}_h$, $\bar{q}_h$, $\bar{p}_c$ and $\bar{q}_c$ allows us to obtain the values of $\bar{q}_h$ and $\bar{q}_c$:

and substituting the aforementioned expressions for $\bar{q}_h$ and $\bar{q}_c$ into the equations for $\bar{p}_h$ and $\bar{p}_c$, respectively, results in the following expressions:

Solving Eq (4.12) for $\bar{p}_h(a)$ and $\bar{p}_c(a)$, yields steady-state solutions for $\bar{p}_h(a)$, $\bar{q}_h(a)$, $\bar{p}_c(a)$ and $\bar{q}_c(a)$ as follows:

It is evident that the system described in Eqs (2.1)–(2.4), (2.13) always has a trivial steady-state.

4.2.   Stability analysis of steady-state solutions

Our next objective is to obtain the stability criteria for a positive steady-state solution. Suppose $p_h(a, t) = \bar{p}_h$, $q_h(a, t) = \bar{q}_h$, $p_c(a, t) = \bar{p}_c$, $q_c(a, t) = \bar{q}_c$, $\forall t\geq0$ represent equilibrium solutions to the PDE model (2.1)–(2.4) and $p_h^*(a, t)$, $q_h^*(a, t)$, $p_c^*(a, t)$ and $q_c^*(a, t)$ represent the corresponding perturbation terms:

After substituting the aforementioned expressions into the PDE model (2.1)–(2.4), we obtain

where $n(t): = \int_0^{a^\star}\left(p_h^*(a, t)+q_h^*(a, t)+p_c^*(a, t)+q_c^*(a, t)\right) \mathrm{d} a$. Then, take the derivative with respect to epsilon $\epsilon$, leads to

which simplifies to

In the limit as $\epsilon$ approaches zero, we arrive at a linear system of partial differential equations:

where $\gamma_i(\bar{N}) = \nu_i\theta_i^{\kappa_i}/ (\theta_i^{\kappa_i}+\bar{N}^{\kappa_i})$, where $i = \{h, c\}$. Next, we formulate (4.14) as semilinear problem:

where the generator $\mathrm{C}$ is defined on the Banach space $X$ as follows:

where

where $\mathrm{D}(\mathrm{C})$ is defined below:

Next, we take the resolvent equation for the operator $\mathrm{C}$:

which leads to

and

By solving (4.18) and (4.20), we get

which after substituting in Eqs. (4.17) and (4.19) gives

Substituting $\phi_1(a)$ and $\phi_3(a)$ back in Eq (4.21) yields

Lemma 4.1. The resolvent of operator $\mathrm{C}$ is compact and its spectrum, denoted by $\sigma(\mathrm{C})$, satisfies the condition:

Here, $\sigma_P(\mathrm{C})$ refers to the point spectrum of $\mathrm{C}$, and $U_\lambda$ is an operator dependent on $\lambda$.

Proof. The expression of $\phi_1(a)$ and $\phi_3(a)$ can be re-written as

where the linear operator on Banach space, $U_{i, \lambda}$ is given as

where

Similarly, we rewrite $\phi_2(a)$ and $\phi_4(a)$ as

where the linear operator on Banach space, $V_{i, \lambda}$ is given as

Let $\Lambda = {\lambda\in\mathbb{C}, |1\in \sigma(U_\lambda)}$. For $\lambda\in\mathbb{C}\backslash \Lambda$, the operators $U_{i, \lambda}$ and $V_{i, \lambda}$ are compact operators from $X$ to $L^1(0, a^\star)$, implying that $\phi_1(a)$ and $\phi_3(a)$ are represented by compact operators, and similarly, $\phi_2(a)$ and $\phi_4(a)$ are also represented by compact operators. As a result, the operator $\mathrm{C}$ has a compact resolvent, which confirms that its spectrum $\sigma(\mathrm{C})$ constitutes only isolated eigenvalues, i.e., $\sigma(\mathrm{C}) = \sigma_P(\mathrm{C})$ (see Theorem 6.29 on page 187 in ^[48]). Hence, $\mathbb{C}\backslash\Lambda\subset\rho(\mathrm{C})$, where $\rho(\mathrm{C})$ is the resolvent of operator $\mathrm{C}$. Therefore, $\sigma_P(\mathrm{C}) = \sigma(\mathrm{C})\subset\Lambda$. Since $U_\lambda$ is a compact operator, we have $\sigma(U_\lambda)\backslash{0} = \sigma_P({U_\lambda)\backslash{0}}$. If $\lambda\in\Lambda$, there exists an eigenfunction $\psi_\lambda$ such that $U_\lambda\psi_\lambda = \psi_\lambda$. It is easy to see that $(\phi_1(a), \phi_2(a), \phi_3(a), \phi_4(a))^\mathrm{T}$ provides an eigenvector of $\mathrm{C}$ for the eigenvalue $\lambda$. Thus, we have $\Lambda\subset\sigma_P(\mathrm{C})$, and we can conclude that (4.22) is satisfied.

Lemma 4.2. Consider the operator ${\mathrm{C}}$ which generates $C_0$-semigroup for $t\geq0$. Then, ${\bf T}(t)$ is eventually norm continuous (ENC), and we have

where $s(\mathrm{C})$ represents the spectral bound of the operator $\mathrm{C}$, and $\omega_0(\mathrm{C})$ denotes the growth bound of the semigroup ${\bf T}(t)$.

Proof. To begin, we express the bounded operator $\mathrm{C}$ as

for $\phi\in X$, $k_h(a) = -\frac{\partial g_h(a)}{\partial a}-\tau_h(a)-\alpha_h(a, \bar{x}_1)-\mu_{p_h}(a)$ and $k_c(a) = -\frac{\partial g_c(a)}{\partial a}-\tau_c(a)-\alpha_c(a, \bar{x}_1)-\mu_{p_c}(a)$. To establish the compactness of $\mathrm{C}$, our strategy is to demonstrate that for any bounded sequence $(\phi^n){n\in \mathbb{N}}$ in $X$, the sequence $(\mathrm{C}\phi^n){n\in \mathbb{N}}$ contains a subsequence that converges uniformly. To accomplish this, we invoke the Arzelà-Ascoli Theorem, which requires us to verify that $(\mathrm{C}\phi^n){n\in \mathbb{N}}$ is uniformly bounded and uniformly equicontinuous. To prove boundedness, since we assume $(\phi^n){n\in \mathbb{N}}$ to be bounded, we get

which determines that $(\mathrm{C}\phi^n){n\in \mathbb{N}}$ is also bounded. For uniform equicontinuity, we consider

Hence, we have shown that $(\mathrm{C}\phi^n){n\in \mathbb{N}}$ is equicontinuous, and by the Arzelà-Ascoli Theorem, we can conclude that there exists a uniformly convergent subsequence of $(\mathrm{C}\phi^n){n\in \mathbb{N}}$. Hence, $\mathrm{C}$ is compact, which in turn implies that ${\bf T}$ is an ENC semigroup. Since the spectral mapping theorem can be applied to ENC semigroups, we have the spectral determined growth condition given by $\omega_0(\mathrm{C}) = s(\mathrm{C})$. Thus, we obtain (4.25).

The local exponential asymptotic stability of the steady-state solution $\omega = 0$ of (4.15) is established when $\omega_0(\mathrm{C}) < 0$. Specifically, there exist constants $\epsilon > 0$, $M\geq1$, and $\gamma < 0$ such that if $x\in X$ and $|x|\leq\epsilon$, then the solution $\omega(t, x)$ of (4.15) exists globally and satisfies $|\omega(t, x)|\leq M \exp(\gamma t)|x|$ for all $t > 0$. In order to examine the stability of steady states, it is necessary to identify the dominant singular point within the set $\Lambda$, which corresponds to the element with the highest real value. By utilizing (4.22) and (4.25), we can then determine the growth bound of the semigroup ${\bf T}$.

Lemma 4.3. For any $\lambda\in \mathbb{R}$, the operator $U_{i, \lambda}$ is nonsupporting with respect to $X_+$ and

holds.

Proof. By Eqs. (4.23) and (4.24), we can conclude that the operator $U_{i, \lambda}, \lambda\in \mathbb{R}$ is strictly positive. To prove that $U_\lambda, \lambda\in \mathbb{R}$ is non-supporting, we can easily demonstrate the inequality

where the linear function $f_{i, \lambda}$, is

Thereby, it leads us to $U_{i, \lambda}^{n+1}\psi\geq\langle f_{i, \lambda}, \psi\rangle\langle f_{i, \lambda}, c\rangle^n c, \, \forall n.$ which holds for all $\psi\in X_+$, where $f_{i, \lambda}$ is strictly positive and the constant function $c = 1$ is a quasi-interior point of $L^1(0, a^\star)$. This implies that $\langle F, U_{i, \lambda}^n\rangle > 0$ for every pair $\psi\in X_+\backslash{0}$, $F\in X_+^{*}\backslash{0}$, and therefore $U_{i, \lambda}, , \lambda\in \mathbb{R}$ is non-supporting. We then use inequality (4.27) and take the duality pairing with the eigenfunctional $F_{i, \lambda}$ of $U_{i, \lambda}$ corresponding to $r(U_{i, \lambda})$, yielding

Assuming $\psi = c$, we obtain the inequality:

where

It follows that

By using the positivity of $\gamma_i(\bar{N}), \, \mu_{p_i}, \, \mu_{q_i}, \, \alpha_i$ and $\tau_i$, we conclude the following:

Hence proved. The previous lemma implies that the function $\lambda\to r(U_{i, \lambda})$ is decreasing for all $\lambda\in \mathbb{R}$. Moreover, if there exists a $\lambda\in \mathbb{R}$ such that $r(U_{i, \lambda}) = 1$, then it follows that $\lambda\in \Lambda$, as $r(U_{i, \lambda})\in\sigma_P(U_{i, \lambda})$. Combining this with the monotonicity property of $r(U_{i, \lambda})$ and inequality (4.26), we obtain the following result.

Lemma 4.4. There exists a unique $\lambda_0\in\mathbb{R}\cap\Lambda$ such that $r(U_{i, \lambda_{0}}) = 1$, and $\lambda_0 > 0$ if $r(U_{0}) > 1$; $\lambda_0 = 0$ if $r(U_{0}) = 1$; $\lambda_0 < 0$ if $r(U_{0}) < 1$.

We will demonstrate that $\lambda_0$ is a dominant singular point, utilizing Theorem 6.13 in ^[49].

Lemma 4.5. If there exists a $\lambda\in\Lambda, \lambda\neq\lambda_0$, then $\mathit{\text{Re}}\lambda < \lambda_0$.

Proof. Let $\lambda\in \Lambda$ and $U_{i, \lambda}\psi = \psi$. Then $|\psi|(a) = |\psi(a)|$, and we have $|U_{i, \lambda}\psi| = |\psi|$. Therefore, we obtain $U_{i, \text{Re}\lambda}\psi\geq\psi$. By taking the duality pairing with $F_{\text{Re}\lambda}\in X^\star_+$, we get $r(U_{i, \text{Re}\lambda})\langle F_{\text{Re}\lambda}, |\psi|\rangle\geq\langle F_{\text{Re}\lambda}, |\psi|\rangle$. We have $r(U_{i, \text{Re}\lambda})\geq 1$, as $F_{\text{Re}\lambda}$ is strictly positive. Since $r(U_{i, \lambda}), \lambda\in\mathbb{R}$ is a declining function, we conclude that $\text{Re}\lambda\leq\lambda_0$. Suppose $\text{Re}\lambda = \lambda_0$. Then $U_{i, \lambda_0}|\psi| = |\psi|$. If we assume $U_{i, \lambda_0}|\psi| > |\psi|$, then taking the duality pairing with the eigenfunctional $F_0$ corresponding to $r(U_{i, \lambda_0}) = 1$ results in $\langle F_{0}, |\psi|\rangle > \langle F_{0}, |\psi|\rangle$, which is a contradiction. Therefore, we have $U_{i, \lambda_0}|\psi| = |\psi|$, and we can deduce that $|\psi| = c\psi_{0}$, where constant $c$ is assumed to be $1$, and $\psi_0$ is the eigenfunction relating to $r(U_{i, \lambda_0}) = 1$. Hence, we have $\psi(a) = \psi_{0}(a)\exp(iv(a))$ for a real-valued function $v(a)$. Substituting this into $U_{i, \lambda_0}\psi_0 = |U_{i, \lambda}\psi|$ yields

Lemma 6.12 ^[49] implies that $\text{Im}\lambda+v(\zeta)$ equals a constant $\Theta$. Using the fact that $U_{i, \lambda}\psi = \psi$, we obtain the equation $\exp(i\Theta)U_{i, \lambda_0}\psi_{\lambda_0} = \psi_{\lambda_0}\exp(iv(\zeta))$. This equation shows that if $\Theta = v(\zeta)$, then $\text{Im}\lambda = 0$. Therefore, the proof is complete.

Theorem 4.1. The equilibrium state $(\bar{p}_h(a), \bar{q}_h(a), \bar{p}_c(a), \bar{q}_c(a))^{\mathrm{T}},$ for (2.1)–(2.4), is locally asymptotically stable if $r(U_{0}) < 1$ and locally unstable if $r(U_{0}) > 1$.

Proof. Lemmas 4.4 and 4.5 suggests that $\sup{\text{Re}\lambda:1\in\sigma_P(U_{i, \lambda})} = \lambda_0$. This implies that if $r(U_{0}) < 1$, then $s(\mathrm{C}) = \sup{\text{Re}\lambda:1\in\sigma_P(U_{i, \lambda})} < 0$. Conversely, if $r(U_{0}) > 1$, then $s(\mathrm{C}) = \sup{\text{Re}\lambda:1\in\sigma_P(U_{i, \lambda})} > 0$. Therefore, the proof is complete.

5.   Results and discussion

In this section, we show simulations of the model to investigate the evolution of healthy and mutated sub-populations of quiescent and proliferative cells in relation to the cell cycle dynamics. These simulations are performed in the Matlab and we have used finite volume method with discretization central upwind scheme. Table 2 displays the model parameters that were used in the simulations. A maximum cell age of $50$ is assumed, and the time step and spatial step size are set to $\Delta t = 0.02$ and $\Delta a = 0.5$, respectively. Additionally, we use unit speed, such that $g_h(a) = g_c(a) = 1$.

Steady-state dynamics of healthy and mutated cell populations:

This case study aims to study the steady-state dynamics of both healthy and mutated cell populations, where the death rates are set to $\mu_{p_h} = \mu_{q_h} = 0.0014$ and $\mu_{p_c} = \mu_{q_c} = 0.0014$. Additionally, $\gamma(N) = 6.8964\times10^{-6}$ and $\rho_1 = 1.0$, while the mutation rate is set to $m = 0.1$. The initial conditions for the healthy cell populations are defined as normal distributions $p_h(a, 0) = q_h(a, 0) = \frac{k_0}{\sqrt{2\pi\sigma^2}}\exp\big(-\frac{(a-\mu)^2}{2\sigma^2}\big)$, where $k_0 = 10^6$, $\mu = 2$, and $\sigma^2 = 200$. Similar initial distributions are used for the mutated cell populations, but with $k_0 = 10^2$. A normal distribution is used because it provides a good approximation of the cell distribution within a population by incorporating heterogeneity with respect to cell age in a population.

Figure 3 illustrates the number density distribution of different cell populations, namely healthy (a) proliferating, (b) quiescent, and mutated (c) proliferating, (d) quiescent, as they evolve and eventually reach a steady-state. Here, the cell age is measured in time and cell density is measured as cells per cubic millimeter. The mutation rate is set to $m = 0.1$. Over time, the population of mutated proliferating and quiescent cells gradually declines as they grow faster and occupy the tissues completely. However, the total population of cells (consisting of both healthy and mutated proliferating and quiescent cells) denoted by N(t) increases rapidly, as shown in Figure 4(a), before eventually reaching a steady-state. Furthermore, Figure 4(b) demonstrates that the growth factors, which are influenced by the cell population N(t), initially increase due to the low cell count and subsequently decrease until reaching an equilibrium. Finally, Figure 4(c) portrays the transition rate of cells, denoted by $\gamma(N)$, from quiescent phase to proliferating phase.

Exponential growth of mutated cell populations:

For this case study, we selected a mutation rate of $m = 0.2$ and used death rates of $\mu_{p_h} = \mu_{q_h} = 0.0014$ and $\mu_{p_c} = \mu_{q_c} = 0.0010$. Additionally, we altered several other parameters, including $\nu_{1, c} = 0.045$, $\rho_{1, c} = 1.0$, $\rho_{2, c} = 30$, $\sigma_{1, c} = 0.040$ and $\sigma_{2, c} = 0.45$. Figure 5 shows the cell density distribution of healthy and mutated proliferating and quiescent cells, respectively. Both healthy subpopulations exhibit the trends of reaching a trivial steady-state with time. However, the mutated cell populations $p_c(a, t)$ and $q_c(a, t)$ increase exponentially, mimicking cancerous behavior. Figure 6 illustrate the total number of cells, growth factors, and the switching function from quiescent to proliferating phase, $\gamma_i(N)$, respectively. The total cell population, which comprises healthy and mutated proliferating and quiescent cell populations, exhibits an exponential increase in cell number over time. The growth factors are initially maximum due to the low cell count and gradually decline until reaching extremely low levels. Finally, the transition functions $\gamma_h$ and $\gamma_c$ also decline as the cell population grows.

The proposed model has some limitations that should be considered. Firstly, exclusion of cell heterogeneity, that is a essential aspect to account for cellular noise. Additionally, the feedback model that is comprised of growth factors is fairly simple straight forward, and to characterize the activation of the Cyclin $\mathsf{D}$ complex, all signaling pathways should be considered. Moreover, at the microscale, it would have been necessary to model cell cycle dynamics separately for healthy and mutated cell populations to investigate their respective compartments' more heterogeneous behavior. Moreover, while the Cyclin-CDK$\mathsf{CDK 4/6}$ are crucial for the $G_1$ to $S$ phase transition, the model overlooks the other restriction point that detects DNA damage in the $S$ phase.

6.   Conclusions

This reserch presents non-linear, multiscale modeling of physiologically-structured healthy and mutated quiescent and proliferating cells in relation to cell cycle dynamics. We modeled reversible transitioning from quiescent to proliferating cells and vice versa. We checked the wellposedness of the model, derive non-trivial equilibrium solutions and find spectral criteria for local stability. We also performed numerical simulations to study the impact of Cyclin $\mathsf{D-CDK 4/6}$ complex on the transition between two sub-populations. Furthermore, we predict that the Cyclin' complex plays an important role in the reversible transition, and any abnormality in this transition can result in cancer.

Conflict of interest

All authors declare no conflicts of interest in this paper.