Developing an Agent-Based Mathematical Model for Simulating Post-Irradiation Cellular Response: A Crucial Component of a Digital Twin Framework for Personalized Radiation Treatment

We define three major cell states: Healthy, Arrested, and Dead. A Healthy cell maintains its basic proliferation potential with no or very light damage. An Arrested cell has its cycle halted in a specific cell-cycle phase. A Dead cell has suffered irreparable damage and suspends material exchange with the extracellular matrix (ECM). We assume that the dead state is an average phenomenon of all the possible cell death routes, such as apoptosis, necrosis, etc. Each of these cell states differs depending on the cell’s phase in the cell cycle: G₁, S, G₂, or M.

The allowed state transitions are:

• •

  From Healthy to Arrested or Dead.

• •

  From Arrested to Dead or Healthy.

Dead cells stay dead. Transitions depend only on a cell’s current state. These state transition rules are based on radiation biology experiments. For instance, 1) A high dose causing direct cell death corresponds to the state transition from Healthy to Dead; 2) A moderate dose causing cell-cycle arrest corresponds to the state transition from Healthy to Arrested; 3) Cell apoptosis after failed cell damage repair during cell-cycle arrest, corresponds to the state transition from Arrested to Dead; 4) A dead cell does not have capacity to repair damage, so theDead state is persistent. Because the dead state is an average phenomenon of different types of death, this formalism could be extended to include other definitions of cellular death types, e.g., cell apoptosis and necrosis. To model cell state transitions after radiation exposure, we propose state energy to quantify the radiation-induced damage to each cell. The state energy includes two parts: the direct state energy $E_{d}$, which quantifies direct radiation-induced damage:

$\displaystyle E_{d}=\alpha N$

(1)

where $\alpha$ is a constant and $N$ is the number of double-strand DNA breaks (DSBs) produced by direct radiation hitting the cell. The DSB is the most lethal DNA damage type to the cell, and here we aggregate all direct damage into an effective DSB number.

The indirect state energy $E_{i}$ quantifies the net damage from bystander signaling. The bystander signals emitted by a signaling cell can migrate and interact with other cells. In this work, we assume that the diffusion-reaction is the major process that contributes to the cell communication by bystander signals, and the receptor-ligand kinetics was considered to model the bystander signal and cell reaction process [9]. We describe the interaction by a second-order reaction depending on the bystander signal concentration and cell receptor concentration:

$\displaystyle E_{i}=\beta\int^{t}_{0}\gamma(t^{\prime})C_{r}(t^{\prime})C_{b}(t^{\prime})dt^{\prime}$

(2)

where $\beta$ is a constant, $t$ is the time after cell irradiation, $\gamma(t^{\prime})$ is reaction rate coefficient at time $t^{\prime}$, $C_{r}(t^{\prime})$ is the transient concentration of cell receptor concentration at time $t^{\prime}$, and $C_{b}(t^{\prime})$ is the transient bystander signal concentration at $t^{\prime}$.

Equation (2) integrates the absorption of bystander signals by the cell. For simplicity of denotation, we denote the integral of reaction rate as $C$, i.e.

$\displaystyle C=\beta\int^{t}_{0}\gamma(t^{\prime})C_{r}(t^{\prime})C_{b}(t^{\prime})dt^{\prime}$

(3)

To be parallel to equation (1), we can write the indirect state energy as:

$\displaystyle E_{i}=\beta C$

(4)

where $C$ is the bystander signal concentration absorbed in the cell. Both targeted and non-targeted cellular effect can lead to different level of DNA damage [10]. In this work, we propose that state energy is a measure of DNA damage induced by radiation including target effect and the bystander effect. State energy is dimensionless, and two types of state energy could be added together, then the total state energy of a cell could be as:

$\displaystyle E_{d}=\alpha N+\beta C$

(5)

When a cell experiences a radiation dose directly or absorbs bystander signals, its state energy increases. When DNA damage is repaired, the state energy decreases. When state energy increases, a cell will have more potential to jump to higher energy cell states. We introduce another term, external perturbation energy, $\Delta E$, to quantify the increase of the state energy. To quantify the energy state transition based on the transition rules, we introduce several basic propositions:

• •

  At a specific time, a cell can and only can stay in one specific energy state, and for each energy state, there is a corresponding state energy distribution. The details of the distribution will be explained below.

• •

  The cell state energy distribution will change if external perturbation energy is nonzero, and correspondingly, a cell will have certain probability of jumping to higher energy states if its state energy increases. A cell will have a certain probability of jumping to lower energy states if its state energy decreases, which corresponds to the cell damage reparation.

• •

  If a cell stays at the highest energy state (i.e. death state), there will be no further state transition. Also, if a cell stays at the lowest energy state (i.e. health state), its state energy cannot decrease.

For a given level of dose to cells, some cells of a given type that are apparently identical will stay healthy, some will have cell cycle arrest and still others will die. In order to quantify the distribution of cell states, we introduce a distribution of cell state energies for each cell state. In this study, we have three cell states: $S_{1}$, $S_{2}$, and $S_{3}$, representing the healthy state, arrested state and death state, respectively. Here we propose that for each cell state $S_{i}$ $(i=1,2,3)$ its state energy E follows a normal distribution $N(E_{i},\sigma^{2}_{i})$ . The state energy $\phi_{i}(E)$ distribution of state $S_{i}$ could be written as:

$\displaystyle\phi_{i}(E)=\frac{1}{\sqrt{2\pi}\sigma_{i}}exp\left(-\frac{(E-E_{i})^{2}}{2\sigma^{2}_{i}}\right)$

(6)

where $E_{i}$ is the mean state energy of state $S_{i}$, and is the standard deviation of state energy of state $S_{i}$. It is noteworthy that $E_{i}$ could be taken as the most feasible state energy for state $S_{i}$, and is a measure of the fluctuation of state energy distribution. One example of a state energy distribution is shown in Figure1. The biological implication of state energy distribution is that we often cannot determine a definite damage level for a certain cell phenotype after irradiation. Taking radiation-induced cell death for example, when we determine the lethal dose which induces cell death, we often observe that the lethal dose is within a dose interval instead of a single lethal dose. The same amount of radiation or bystander signal can cause stochastically different amounts of damage to identical cells.

When radiation-induced DNA damage is incurred, the cell state energy will increase, which will induce a change in the cell state energy distribution correspondingly. We assume that the shape of the state energy distribution will not change, but the mean state energy will shift. For instance, suppose the current energy state of the cell is $S_{i}$, then if cell absorbs perturbation energy $\Delta E$, its state energy distribution will be as:

$\displaystyle\phi_{i}(E)=\frac{1}{\sqrt{2\pi}\sigma_{i}}exp\left(-\frac{(E-(E_{i}+\Delta E))^{2}}{2\sigma^{2}_{i}}\right)$

(7)

After defining the state energy distribution, it is plausible to quantify the probability of state transitions. The cell state transition probability is quantified by calculating the overlapping integral of the state energy distribution.

The overlapping integral of two state energy distributions is written as:

$\displaystyle\left\langle\phi_{i}|\phi_{j}\right\rangle=\int_{D}\min[\phi_{i}(x),\phi_{j}(x)]dx$

(8)

where $\phi_{i}(x)$ and $\phi_{j}(x)$ are the state energy distribution functions of two cell energy states.

The computation of the overlapping integral for two state energy distributions, $N(E_{i},\sigma_{i}^{2}),N(E_{j},\sigma_{j}^{2})$, depends on whether the two variances are equal or not. Here we start with the simpler case where we have $\sigma_{1}^{2}=\sigma_{2}^{2}=\sigma^{2}$. We then have

$\displaystyle\left\langle\phi_{i}|\phi_{j}\right\rangle=2\Phi\left(-\frac{\lvert E_{i}-E_{j}\rvert}{2\sigma}\right)$

(9)

where $\Phi(x)$ is the cumulative distribution function of the normal distribution and it is as

$\displaystyle\Phi(x)=\int^{x}_{-\infty}\frac{1}{\sqrt{2\pi}}exp\left(-\frac{t^{2}}{2}\right)dt$

It is obvious that $0\leq\left<\phi_{i}|\phi_{j}\right>\leq 1$, and $\left<\phi_{i}|\phi_{j}\right>=1$ if and only if those two normal distributions are identical. We use the overlapping integral as a measure of cell state transition probability between two states, and the transition probability from state $S_{i}$ to $S_{j}$ is:

$\displaystyle P(S_{i}\rightarrow S_{j})=\left<\phi_{i}|\phi_{j}\right>$

(10)

Suppose initially at time $t$, a cell stays at state $S_{i}$. If it absorbs an external perturbation energy $\Delta E$ and its state energy distribution will shift to a possible higher energy state $S_{j}$, then the transition probability from $S_{i}$ to $S_{j}$ is:

$\displaystyle P(S_{i}\rightarrow S_{j})=2\Phi\left(-\frac{|E_{i}+\Delta E-E_{j}|}{2\sigma}\right)$

(11)

Based on equation (11), we can calculate the probabilities for all the possible transition routes. The possible cell state transition route is shown in Figure 2. The final transition route is selected based on a rejection algorithm based on Monte Carlo sampling [11]. The selection process is as follows:

• •

  if $\displaystyle\xi\leq\frac{P(S_{1}\rightarrow S_{1})}{P(S_{1}\rightarrow S_{1})+P(S_{1}\rightarrow S_{2})+P(S_{1}\rightarrow S_{3})}$, where $0<\xi<1$, the cell remains in $S_{1}$

• •

  if $\displaystyle\frac{P(S_{1}\rightarrow S_{1})}{P(S_{1}\rightarrow S_{1})+P(S_{1}\rightarrow S_{2})+P(S_{1}\rightarrow S_{3})}<\xi\leq\frac{P(S_{1}\rightarrow S_{2})}{P(S_{1}\rightarrow S_{1})+P(S_{1}\rightarrow S_{2})+P(S_{1}\rightarrow S_{3})}$, the cell transitions into $S_{2}$

• •

  Otherwise, cell transitions into $S_{3}$

Cells are known to have different radiation sensitivity in different cell phases. Cells in the G₂ and M phases usually have higher radiation sensitivity, lower radiation sensitivity in the G₁ phase, and the lowest radiation sensitivity during the latter part of the S phase. We define a radiation sensitivity factor $f$ for each cell phase. Here, we introduce a method to calculate the radiation sensitivity factor using the cell state transition rule we just introduced above.

For nomenclature simplicity, each cell phase is assigned an index i corresponding to the four cell phases where $i\in\{1,2,3,4\}$ and $i=1$ corresponds to G₁, $i=2$ corresponds to S, $i=3$ corresponds to G₂, and $i=4$ corresponds to M. We denote $E_{i,j}$ as the mean state energy of state $S_{j}$ at $i^{th}$ cell phase, i.e., $S_{i,j}$ where $i\in\{1,2,3,4\}$ and $j\in\{1,2,3\}$.

For cell phase G₁, the transition probability from $S_{1}$ to $S_{3}$ after cell absorbs external perturbation energy $\Delta E$ is

$\displaystyle P_{1}(S_{1}\rightarrow S_{3})=2\Phi\left(-\frac{|E_{1,1}+\Delta E-E_{1,3}|}{2\sigma}\right)$

(12)

Now we consider a less radiosensitive cell phase, S, the transition probability from $S_{1}$ to $S_{3}$ after absorbs energy $\Delta E$ is

$\displaystyle P_{2}(S_{1}\rightarrow S_{3})=2\Phi\left(-\frac{|E_{2,1}+\Delta E-E_{2,3}|}{2\sigma}\right)$

(13)

From equation (12) and equation (13) we can know that transition probability is a function of . Here, we calculate the derivatives of equation (12) and equation (13) with respect to $\Delta E$. We can have

$\displaystyle\frac{dP_{2}(S_{1}\rightarrow S_{3})}{d\Delta E}=-\frac{1}{\sigma}\,g\,\left(-\frac{|E_{1,1}+\Delta E-E_{1,3}|}{2\sigma}\right)\frac{E_{1,1}+\Delta E-E_{1,3}}{|E_{1,1}+\Delta E-E_{1,3}|}$

(14)

where $g=\frac{1}{\sqrt{2\pi}}exp(-\frac{t^{2}}{2})$ which is the normal distribution density function. Similarly, we can get

$\displaystyle\frac{dP_{1}(S_{2}\rightarrow S_{3})}{d\Delta E}=-\frac{1}{\sigma}\,g\,\left(-\frac{|E_{2,1}+\Delta E-E_{2,3}|}{2\sigma}\right)\frac{E_{2,1}+\Delta E-E_{2,3}}{|E_{2,1}+\Delta E-E_{2,3}|}$

(15)

We define a radiation sensitivity factor f for each cell phase. We take $f_{1}=1$ for G₁ phase, and the radiation sensitivity factors in other cell phases could be normalized according to $f_{1}$. We propose that the radiation sensitivity factor $f_{2}$ for the cell phase G₂ as:

$\displaystyle f_{2}=\lim_{\Delta E\to 0}\dfrac{\dfrac{dP_{2}(S_{1}\rightarrow S_{3})}{d\Delta E}}{\dfrac{dP_{1}(S_{1}\rightarrow S_{3})}{d\Delta E}}$

(16)

Then according to equation (14) and equation (15), we can get

$\displaystyle f_{2}=exp\left(-\frac{(E_{2,1}-E_{2,3})^{2}}{8\sigma^{2}}+\frac{(E_{1,1}-E_{1,3})^{2}}{8\sigma^{2}}\right)$

(17)

Then we can have

$\displaystyle(E_{1,1}-E_{1,3})^{2}-(E_{2,1}-E_{2,3})^{2}=8\sigma^{2}\ln f_{2}$

(18)

Equation (18) shows the relationship between radiation sensitivity factor and the mean state energies.

It is known that the cell cycle could be temporarily or permanently arrested for a while after cell irradiation. Here, we propose a method to calculate the cell arrest duration based on the cell state transition model. As well-studied experimentally [12, 13], the kinetics of DSB rejoining exhibits a fast initial rate, which then decreases with repair time. The most widely accepted description of this kinetic behavior uses two first-order components (fast and slow). This general two-repair-components model is used widely [14, 15, 16]. The main common characteristic of the model is assuming the fast repair and slow repair follow a first-order exponential decay scheme. Then the total DSB number after irradiation could be described as:

$\displaystyle N(t)=f_{fast}\,e^{-\lambda_{1}t}+f_{slow}\,e^{-\lambda_{2}t}$

(19)

where $N(t)$ is the DSB number in cell at time $t$, $f_{fast}$ and $f_{slow}$ correspond to the weights of the two repair components respectively, $\lambda_{1}$ and $\lambda_{2}$ corresponds to the decay constants of repair components respectively. And we have $f_{fast}=1-f_{slow}$.

We introduce two definitions for mathematical formalization. The transient arrest state is the state which will be released after a certain amount of time. The transient arrest duration is the total time for a cell staying at transient arrest state.

The cell state energy proposed in our model is a measure of DNA damage level. In this study, we propose that state energy also follows the same kinetics as radiation-induced DSB. The state energy with respect to time could be expressed as:

$\displaystyle E(t)=E_{0}(f_{1}\,e^{-\lambda_{1}t}+f_{2}\,e^{-\lambda_{2}t})$

(20)

When a cell jumps to $S_{2}$ from $S_{1}$, the state energy of $S_{2}$ will decrease with time because of the DNA damage repair. The probability of cell state transition from $S_{2}$ to$S_{1}$ at time $t$ will be as:

$\displaystyle P(t)=2\Phi\left(-\frac{\lvert E_{2}(t)-E_{1}\rvert}{2\sigma}\right)\,$

(21)

According to equation (20), we can know that

$\displaystyle E_{2}(t)=E_{2}(0)(f_{1}\,e^{-\lambda_{1}t}+f_{2}\,e^{-\lambda_{2}t})$

(22)

where $E_{2}(0)$ is the initial cell state energy at $S_{2}$. We take $E_{1}=0$, so equation (21) could be written as:

$\displaystyle P(t)=2\int^{-\frac{E_{2}{t}}{2\sigma}}_{-\infty}\frac{1}{\sqrt{2\pi}}\,exp\left(-\frac{x^{2}}{2}\right)dx\,,$

(23)

Then we can know the jumping probability rate is as:

$\displaystyle\frac{dP(t)}{dt}=-\frac{1}{\sqrt{2\pi}\sigma}\,exp\left(-\frac{E^{2}_{2}(t)}{8\sigma^{2}}\right)\frac{dE_{2}(t)}{dt}\,,$

(24)

Based on equation (24), we can obtain the mean life time of $S_{2}$ which is the mean duration of transient arrest state. The mean lifetime of $S_{2}$ could be obtained by

$\displaystyle T_{m}=\int^{\infty}_{0}t\,\frac{P(t)}{dt}dt$

(25)

Substituting $E_{2}(t)$ back into equation (23) then the integration in equation (25) could be changed to as:

$\displaystyle T_{m}=\frac{E_{2}(0)}{\sqrt{2\pi}\sigma}\int^{\infty}_{0}t\,exp\left(-\frac{E^{2}_{2}(t)}{8\sigma^{2}}\right)(f_{1}\,\lambda_{1}\,e^{-\lambda_{1}t}+f_{2}\,\lambda_{2}\,e^{-\lambda_{2}t})dt$

(26)

We can approximate equation (26) with the expansion of the term $exp(-\frac{E^{2}_{2}(0)}{8\sigma^{2}})$ neglecting the higher order terms in the integration and we get the mean lifetime of $S_{2}$ that could be calculated as:

$\displaystyle T_{m}\approx\frac{E_{2}(0)}{\sqrt{2\pi}\sigma}\left(\frac{f_{1}}{\lambda_{1}}+\frac{f_{2}}{\lambda_{2}}\right)$

(27)

We see that the cell arrest duration could be determined by DNA repair kinetics and state energy distribution. From equation (27) we can know that the slow repair phase of DSB, which has the smaller half-life will dominate the transition duration. This is in line with one speculative view published in [17] that the slow homologous recombination repair (HRR) of DSB determines the length of cell arrest duration. Transient arrest duration is a measure of the time needed for cell state recovery.

As soon as the duration of a given cell phase has passed, the transition to the next phase of the cell cycle occurs. The time at which the transition takes place varies in a random manner according to a distribution of durations of the cell cycle phases. A probability density function $f_{i}(t)$ can be assigned to each phase such that the probability of the duration time for the $i^{th}$ phase having a value in the interval around t is given by $f_{i}(t)dt$. In this work, $f(t)$ is defined as a Gaussian distribution [19]. In each phase, the phase duration time, $\bar{t}_{i}$, and a variance $\sigma_{i}^{2}$, which are identified with the mean and variance of $f_{i}(t)$. The probability density function of duration time is a half Gaussian distribution, and it is written as:

$$f_{i}(t)=\begin{cases}0,&\text{$t<T_{i}$}\\ \dfrac{1}{\sqrt{2\pi}\sigma}\,exp(-\dfrac{(t-\bar{t}_{i})^{2}}{2\sigma^{2}}),&\text{$t\geq T_{i}$}\end{cases}$$

(29)

Naturally, without any perturbation, a cell will go through the normal cell cycle, and we can simulate the cell cycle progression based on the probability density distribution in equation [29]. We can use Monte Carlo sampling to determine the duration of one cell phase. For phase G_i, we can sample a duration length for it as

$\displaystyle T_{i}=t_{i}+\sigma_{i}\,\sqrt{-2\ln\xi_{1}}\cos(\frac{\pi}{2}\xi_{2})$

(30)

where $\xi_{1}$ and $\xi_{2}$ are random numbers.

We can determine the transition rule determined by time evolution as:

$$F_{1}(S^{t}_{i},\eta^{t}_{i})=\begin{cases}\text{cell goes to next phase},&\text{$t>T_{i}$}\\ \text{cell stays in current phase},&\text{$t\leq T_{i}$}\end{cases}$$

(31)

where $F_{1}(S_{i}^{t},\eta_{i}^{t})$ is used to describe the cell phase transition result.

Cell cycle progression is dependent on external conditions, specifically the local nutrient availability (such as glucose concentration in the medium) and interaction with neighboring cells (integral pressure) [20]. In this work, we assume that all the cells are in good nutrient condition. The interaction with neighboring cells is considered by checking the contact inhibition condition of the cell. The cell will stay at quiescent state if there is no space for the daughter cells [21]. Regarding where to position the two daughter cells, one of the common rules is to position them randomly at the adjacent vacant sites. Taking an example as shown in Figure 3, if the grid lattice $L(i,j)$ is filled by a cell, then we define a variable $P(i,j)=1$, otherwise $P(i,j)=0$. Then we can determine whether there is space for new cells in the system by checking the value of $\eta_{i}^{t}$.

$\displaystyle\eta^{t}_{i}=P(i-1,j)+P(i+1,j)+P(i,j-1)+P(i,j+1)$

(32)

and apparently $\eta^{t}_{i}\in{0,1,2,3,4}$.

The transition rule determined by neighbor condition is as:

$$F_{2}(S^{t}_{i},\eta^{t}_{i})=\begin{cases}\text{cell goes to next phase},&\text{$\eta^{t}_{i}>0$}\\ \text{cell stays in current phase},&\text{$\eta^{t}_{i}=0$}\end{cases}$$

(33)

Then we can write the transition rule for cell phase transition as:

$\displaystyle F(S^{t}_{i},\eta^{t}_{i})=F_{1}(S^{t}_{i},\eta^{t}_{i})\cap F_{2}(S^{t}_{i},\eta^{t}_{i})$

(34)

It is worth noting that the transition rule could be easily extended by incorporating other conditions which determine the cell phase transition process.

In this study, we used the cell phase ratio to describe the cell phase distribution, and the cell phase ratio of phase $P_{i}$ is defined as:

$\displaystyle f_{P_{i}}=\frac{N_{P_{i}}}{N_{P}}$

(35)

where $f_{P_{i}}$ is the cell phase ratio of $i^{th}$ cell phase, $N_{P_{i}}$ is the number of cells staying at $i^{th}$ phase and $N_{P}$ is the total number of cells.

The cell state transition is a continuous stochastic process evolving with time after cell irradiation. The cell state transition probability obtained by equation (11) is the total probability during the time until actual cell state transition (phenotype) is observed. In the agent-based model, the time is discretized into time steps and the probability of cell state transition should be determined for each time step. So here we introduce a model to quantify the probability of cell state transition in each time step. Firstly, the external perturbation energy $\Delta E$ for each cell is determined, then the probability of corresponding cell state transition is calculated based on $\Delta E$ in the time step. Secondly, the cell state transition decision is made based on the calculated transition probability.

The cell state transition is dependent on sources which contribute to $\Delta E$. For instance, DSBs and bystander signals all can induce $\Delta E$. The direct DSB is a more effective lethal damage compared to bystander signals, so it can induce the fast cell state transition. Relatively, the bystander signal is a less damaging agent which is subject to relatively constant and cyclic stresses to cells since the bystander signal will perpetuate in the cell culture until it decays away. In this work, considering the difference between those two external perturbation factors, we propose two types of cell state transition modes, i.e., delayed transition and instantaneous transition.

The instantaneous transition happens immediately after the cell absorbs external perturbation energy. The delayed transition is a relatively slow transition that happens after a certain time from the time when the cell absorbs external perturbation energy. For the delayed transition, the time-to-transition, $t$, varies. In this work, we introduce a model to quantify the state transition as a failure process. First, we introduce a term, observation time $T_{o}$, within which the radiation-induced cellular phenotype could be observed. The exact time of state transition is not known, and we can assume the time-to-transition as a random variable, $T$, and it has a density function $f(t)$, then we can have

$\displaystyle F(t)=P(T\leq t),\,\,\,\,\,\,\,\,\,t\geq T_{o}$

(36)

where $F(t)$ is the cell state transition distribution function. The density function can be described in terms of $T$ as:

$\displaystyle f(t)=\lim_{\Delta t\to 0}P(t<T\leq t+\Delta t)$

(37)

This can be interpreted as the probability that the cell state transition time $T$ will occur between the time $t$ and the next time $t+\Delta t$.

Within the observation time, the total probability of cell state transition could be calculated by

$\displaystyle p_{o}=\int^{T_{o}}_{0}f(t)dt$

(38)

Considering the cell state transition from $S_{i}$ to $S_{j}$, then $p_{o}$ is determined as:

$\displaystyle p_{o}=2\Phi\left(-\frac{|E_{i}+\Delta E-E_{j}|}{2\sigma}\right)$

(39)

Here we propose that $f(t)$ follows exponential distribution, and it is

$\displaystyle f(t)=\lambda e^{-\lambda t}$

(40)

Then we can know that

$\displaystyle\int^{T_{o}}_{0}\lambda e^{-\lambda t}dt=p_{o}$

(41)

and from equation (41) we can solve out as

$\displaystyle\lambda=-\frac{\ln(1-p_{o})}{T_{o}}$

(42)

Without losing the generality, considering a time step $\Delta t$ starting at time $t$, then we can know the probability of state transition in the time step will be as

$\displaystyle p(t<t^{\prime}\leq t+\Delta t\,|\,T>t)=\frac{p(t<T\leq t+\Delta t)}{p(T>t)}$

(43)

where $t^{\prime}$ is the time when cell completes cell state transition within the time step, then we can obtain

$\displaystyle p(t<t^{\prime}\leq t+\Delta t\,|\,T>t)=\frac{\int^{t+\Delta t}_{t}\lambda e^{-\lambda x}dx}{1-\int^{t}_{0}\lambda e^{-\lambda x}dx}=1-e^{-\lambda\Delta t}$

(44)

which is the cell state transition probability in each time step.

The time step $\Delta t$ is chosen reasonably small considering the calculation cost and the accuracy. In this work, $\Delta t$ is chosen as 1 minute. The cell state transition probability $p$ in each time step will be used to determine whether the cell will go through cell state transition or not. In this work, the rejection sampling rule is used to determine the transition results. The transition rule is defined as:

$$F(\lambda,\Delta t)=\begin{cases}\text{cell undergoes state transition},&\text{if $p<\xi$}\\ \text{cell stays in current phase},&\text{if $p\geq\xi$}\end{cases}$$

(45)

where $\xi$ is a random number, and $F(\lambda,\Delta tt)$ is used to describe the cell state transition result.

The cell state distribution of the cell culture could be obtained according to the cell state transition rule. In this study, we used the cell state ratio to describe the cell phase distribution, and the cell state ratio of state $S_{i}$ is defined as

$\displaystyle f_{S_{i}}=\frac{N_{S_{i}}}{N_{S}}$

(46)

where $N_{S_{i}}$ is the number of cells staying in cell state $S_{i}$, and $N_{S}$ is the number of total cells.