Evidence that PUBO outperforms QUBO when solving continuous optimization problems with the QAOA

Jonas Stein

LMU Munich
Oettingenstr. 67

Farbod Chamanian

LMU Munich
Oettingenstr. 67

Maximilian Zorn

LMU Munich
Oettingenstr. 67

Jonas Nüßlein

LMU Munich
Oettingenstr. 67

Sebastian Zielinski

LMU Munich
Oettingenstr. 67

Michael Kölle

LMU Munich
Oettingenstr. 67

Claudia Linnhoff-Popien

LMU Munich
Oettingenstr. 67

Quantum computing provides powerful algorithmic tools that have been shown to outperform established classical solvers in specific optimization tasks. A core step in solving optimization problems with known quantum algorithms such as the Quantum Approximate Optimization Algorithm (QAOA) is the problem formulation. While quantum optimization has historically centered around Quadratic Unconstrained Optimization (QUBO) problems, recent studies show, that many combinatorial problems such as the TSP can be solved more efficiently in their native Polynomial Unconstrained Optimization (PUBO) forms. As many optimization problems in practice also contain continuous variables, our contribution investigates the performance of the QAOA in solving continuous optimization problems when using PUBO and QUBO formulations. Our extensive evaluation on suitable benchmark functions, shows that PUBO formulations generally yield better results, while requiring less qubits. As the multi-qubit interactions needed for the PUBO variant have to be decomposed using the hardware gates available, i.e., currently single- and two-qubit gates, the circuit depth of the PUBO approach outscales its QUBO alternative roughly linearly in the order of the objective function. However, incorporating the planned addition of native multi-qubit gates such as the global Mølmer-Sørenson gate, our experiments indicate that PUBO outperforms QUBO for higher order continuous optimization problems in general.

<concept_desc>Hardware Quantum computation</concept_desc>
<concept_desc>Mathematics of computing Continuous functions</concept_desc>

[500]Hardware Quantum computation
[500]Mathematics of computing Continuous functions

§ INTRODUCTION

Solving optimization problems is a central task in industries involving domains like production and logistics. Many of these problems concern scheduling, routing, packing and others, which are often NP-hard and thus demand for heuristic solvers. A particularly promising approach to solving such optimization problems is quantum computing, which has already shown results comparable to classical state-of-the-art methods for small problem sizes <cit.> despite current quantum hardware limitations. For a significant period of time, quantum optimization was driven by D-Wave System's Quantum Annealing devices, which are technically limited to solving problems written in Quadratic Unconstrained Binary Optimization (QUBO) form. This restriction was subsequently lifted in the Quantum Approximate Optimization Algorithm (QAOA) by Farhi et al., which essentially simulates the process of Quantum Annealing on a quantum gate computer and allows for additional generalization using the larger capabilities of a universal quantum computer <cit.>.

One particularly powerful generalization of the QAOA is its ability to solve higher order polynomial problems, i.e., it can natively work with Polynomial Unconstrained Binary Optimization (PUBO) problems. Instead of having to quadratize the a PUBO problem into QUBO form using ancillary qubits as is necessary for D-Wave's Quantum Annealers, needed multi-qubit interactions can be modelled using quantum gates <cit.>. While current quantum computers generally only support single- and two-qubit gates, e.g., trapped ion quantum computers are expected to implement multi-qubit gates such as the (global) Mølmer-Sørenson gate in the future[<https://ionq.com/docs/getting-started-with-native-gates>]. Such gates will allow the execution of the qubit interactions necessary to model PUBO problems in constant time without the currently needed decomposition in two- and single-qubit gates <cit.>, which scales linearly in the number of qubits involved.

While some binary, combinatorial optimization problems like Max-Cut or Number Partitioning are formulated in terms of QUBO natively, modelling intrinsically non-binary problems like the TSP for QUBO requires special encoding techniques like the one-hot encoding, which increase the search space beyond exigence <cit.>. For problems like these, it has been shown that their PUBO versions generally outperform their QUBO analogues in terms of solution quality as well as the required number of optimization steps and QAOA iterations <cit.>.

As many NP-hard problems such as scheduling or packing also involve continuous variables in higher order terms frequently in application <cit.>, we set out to compare the performance of PUBO and QUBO formulations for the QAOA on continuous optimization problems. Our two core contributions to this investigation are:

* an implementation of the QAOA capable of solving arbitrary polynomial optimization problems, that allows control over the used bit depth and the domains of the input variables, and

* an in-depth case-study evaluating the performance of PUBO and QUBO problem formulations on two established, continuous optimization benchmark functions.

This paper is structured into five sections. Following this introduction, we visit fundamental background knowledge necessary to comprehend our methodology in section <ref>. Section <ref> subsequently contains a detailed description of the concept used to solve higher order continuous optimization problems with the QAOA. Finally, the established approach is applied to conduct the aspired evaluation in section <ref> while concluding with a contextualization of the acquired results in section <ref>.

§ BACKGROUND

In this section, we describe the overall functionality of the QAOA and its initial motivation to get an overview of all its components possibly influencing the evaluation results.

The QAOA is inspired by Adiabatic Quantum Computing (AQC), which is an alternative paradigm of quantum computing besides the omnipresent Quantum Gate Model (QGM). The main difference of AQC to the QGM resides in its time evolution being inherently continuous instead of iteratively applying discrete gates, as done in the QGM. Drawing upon the adiabatic theorem, which essentially states that a physical system stays in its instantaneous eigenstate whenever the time evolution applied to it happens slowly enough and if there is a gap between the corresponding eigenvalue and the rest of the Hamiltonian's spectrum <cit.>, an optimization algorithm can be formulated as:

* Prepare an initial state $\ket{\psi}$ that is the ground state of a known Hamiltonian $\hat{H}_M$.

* Identify a Hamiltonian $\hat{H}_C$ modelling the objective function $f:\left\lbrace 0,1\right\rbrace^{n}\rightarrow \mathbb{R}$ where the eigenstates represent possible solutions to the input problem. The eigenvalues that correspond to the eigenstates embody the objective values of the respective solution.

*  Gradually evolve the initial state to the ground state of $\hat{H}_C$ corresponding to the global optimum of $f$ by applying the Hamiltonian $\hat{H}(t)=\left(1-t\right)\hat{H}_M + t\hat{H}_C$.

The standard choice for the Hamiltonian $\hat{H}_M$ is $\hat{H}_M\coloneqq - \sum_{i=1}^{n} \sigma^x_i$ which inherits the easy to prepare ground state $\ket{+}^{\otimes n}$, where $\sigma^x_i$ denotes the tensor product of $n-1$ identity matrices $I$ with the Pauli operator $\sigma_x$ at the $i$-th position. For $\hat{H}_C$, a possible definition is $\hat{H}_C\coloneqq \sum_{x\in\left\lbrace 0,1\right\rbrace^{n}} f(x)\ket{x}\bra{x}$ as this trivially matches its requirements stated above.

While Quantum Annealers are built to execute the procedure described in item <ref> for any given Ising Hamiltonian[Ising Hamiltonians represent the energy spectrum in a specific physical system. This system is described by an Ising model, which is a mathematical model of ferromagnetism in statistical mechanics. This Hamiltonian has the convenient property of being isomorphic to the NP-hard quadratic programming problem and hence naturally allows to model many interesting optimization problems with it.] $\hat{H}_C=\sum_i h_i \sigma^z_i + \sum_{i<j}J_{ij}\sigma^z_i\sigma^z_j$, discretization and Hamiltonian simulation techniques must be used to implement this time evolution in the QGM, which is the fundamental idea of the QAOA. The continuous time evolution of $\hat{H}(t)$ is discretized by iteratively simulating the time evolution of the Hamiltonians $\hat{H}(t_k)$ with equidistant $t_k\in\left[0, 1\right]$ strictly increasing from $0$ to $1$ and $k\in\left\lbrace 1, ..., P\right\rbrace$.

To perfectly approximate the continuous time evolution in the limit for $P\rightarrow \infty$, each Hamiltonian $\hat{H}(t_k)$ is chosen to act for time $\nicefrac{1}{p}$. However, especially for small $P$, it is typically unclear how quickly the time evolution should progress at each intermediate Hamiltonian. In this context, it has proven useful to introduce parameters associated with the duration of their time evolution. These parameters can then be used to, i.a., satisfy the conditions of the adiabatic theorem, given that $P$ is big enough. Notably the concrete implementations proposed for this parameterization use independent parameters for both Hamiltonians: $\gamma_k\in\mathbb{R}$ for the Hamiltonian $\hat{H}_{C}$ and $\beta_k\in\mathbb{R}$ for the Hamiltonian $\hat{H}_{M}$. This allows for increased flexibility, especially in the regime of low $P$. For the optimization of these parameters, many different approaches have been explored, foremost gradient based techniques like the parameter shift rule in combination with gradient descent <cit.>, but also other heuristic approaches focused on yielding results very quickly, such as the COBYLA optimizer <cit.>.

The QAOA algorithm can thus be understood as an algorithm, that simulates the time evolution of the Hamiltonian $\hat{H}(t)$ on gate based quantum computers. It does so using parameters guiding the time evolution speed as displayed in figure <ref>.

$\ket{0}^{\otimes n}$     H^⊗n     e^-ip (γ_1 (1- t_1)H_M + β_1 t_1Ĥ_C)      …     e^-ip (γ_p (1- t_p)H_M + β_p t_pĤ_C)     $\ket{\psi\left(\beta , \gamma\right)}$
The general form of the QAOA circuit.

§ CONCEPT

In this section, we show how the QAOA can be used to solve higher order continuous polynomial optimization problems. More specifically, we employ the following procedures:

* Discretization of the objective function

* Translating the objective function into a Hamiltonian

* Implementing the Hamiltonian using quantum gates

§.§ Discretization of the objective function

For discretizing a given objective function $f:\left[a,b\right] \rightarrow \mathbb{R}$ with $a<b\in \mathbb{R}$, we need to select a suitable bit encoding. For the sake of simplicity, we choose the sign–magnitude representation which maps any integer to its native binary encoding while initially disregarding its sign, to then finally represent its sign using an extra bit at the start, e.g.: $3_{10} \mapsto{} 0\, 11_2$ and $-3_{10} \mapsto{} 1\, 11_2$. In addition to that simplification, we also restrict the possible domain spaces of each variable to be of the form $\left]-2^n, 2^n\right[$ where $n\in\mathbb{N}$, to alleviate needed precautions for intervals that are unbalanced or away from powers of two. This decision allows us to incorporate numbers beyond the whole numbers in a straightforward manner, i.e., by using standard floating point representation with a freely selectable bit resolution $m\in\mathbb{N}$. The complete binary encoding of a given $x\in\left]-2^n, 2^n\right[$ and bit resolution $m\in\mathbb{N}$ can thus be described by the following approximation:

\begin{align}
x \approx \left(2x_0-1\right)\left(\sum_{i=1}^{n} 2^{n-i} x_i + \sum_{i=1}^{m}x_{n+i}2^{-i}\right)
\end{align}

As desired, this discretization leads to the bit string representation $x\approx x_0 \, x_1 ... x_n,x_{n+1} ... x_{n+m}$, so that, e.g., $\left]-2^2, 2^2\right[\ni -2,75_{10} \mapsto{} 1\, 10,110_2$ for a bit resolution of $m=3$. Note however, that the borders of the domain space can only be approached when increasing the bit resolution $m$, while every additional bit contributes with advancement of $\nicefrac{1}{2^{m+1}}$. Using this bit encoding, we can also represent functions with higher dimensional input spaces by following the described substitution procedure for every dimension and then concatenating the resulting bit strings.

§.§ Translating the objective function into a Hamiltonian

As described in section <ref>, there is a native mapping between binary functions $f:\left\lbrace 0,1\right\rbrace^{n}\rightarrow \mathbb{R}$ and Hamiltonians, i.e., $\hat{H}_C\coloneqq \sum_{x\in\left\lbrace 0,1\right\rbrace^{n}} f(x)\ket{x}\bra{x}$. This method can be very inefficient however, if we only have access to $f$ as a black box function, because the Hamiltonian can be comprised of exponentially many non-zero terms. Given that we have access to $f$ in a white box manner, we can conduct this mapping much more efficiently, i.e., by substituting every $x_i\in\left\lbrace 0,1\right\rbrace$ with a $s_i\in\left\lbrace -1,1\right\rbrace$ as in $x_i\mapsto\nicefrac{\left(s_i+1\right)}{2}$. In the case of $f$ having higher degree interactions than two in its input bits (i.e., e.g., a term like $\alpha x_0x_1x_2$ with $\alpha\in\mathbb{R}$), inserting a suitable quadratization step is obligatory for the QUBO version. Typically this step is done before translating into the spin configuration domain $\left\lbrace -1,1\right\rbrace$ by adding ancillary bits to the input space and a penalty term to the function $f$, as exemplified in equation <ref>. For details on this quadratization step, we reference to the python package qubovert, which we used for this step in our implementation[<https://github.com/jtiosue/qubovert>]. Notably, finding the optimal quadratization in terms of minimizing the number of needed ancillary qubits is NP-hard, as pointed out in <cit.>.

\begin{align}
f(x_0,x_1,x_2)&=\alpha x_0x_1x_2 \nonumber\\
\mapsto\; f(x_0,x_1,x_2,z)&=\alpha zx_2+2\alpha \left(x_0x_1-2\left(x_0+x_1\right)z+3z\right)
\label{eq:quadratization}
\end{align}

In order to translate the resulting function of spin configurations $f':\left\lbrace -1,1\right\rbrace^{n}\rightarrow \mathbb{R}$ into a quantum mechanical Hamiltonian, we can simply substitute all spins $s_i$ with Pauli operators using the trivial map $s_i\mapsto\sigma_{i}^z$. <cit.>

§.§ Implementing the Hamiltonian using quantum gates

To implement the quantum circuit of the QAOA, we have to conduct Hamiltonian simulation of $\hat{H}_M$ and $\hat{H}_C$. While $\hat{H}_M$ can easily be simulated using parameterized $X$ gates, $\hat{H}_C$ involves higher order terms (as e.g., $\alpha \sigma_{i}^z \sigma_{j}^z \sigma_{k}^z $ where $\alpha \in \mathbb{R}$) for the PUBO variant. As pointed out in <cit.>, Hamiltonians of this form can be simulated using the generic architecture shown in figure <ref>, naturally expanding from the well-know quadratic case $\alpha \sigma_{i}^z \sigma_{j}^z$. When having access to a suitable multi-qubit gate such as the (global) Mølmer-Sørenson gate, combining the information presented in figure 4.19 in <cit.> and figure 5 from <cit.>, we can simulate arbitrary degrees of Pauli matrices using one extra ancillary qubit with an overhead of merely two extra circuit operations. As all terms in $\hat{H}_C$ commute, the Hamiltonian simulation simplifies into a concatenation of the gates used to implement all terms in the sum notation of $\hat{H}_C$ as exemplified in figure <ref>, concluding this section.

Hamiltonian simulation of $\dfrac{\theta}{2}\sigma_{0}^z\sigma_{1}^z$.

1         1

Hamiltonian simulation of $\dfrac{\theta}{2}\sigma_{0}^z\sigma_{1}^z$.

1                 1

1         1

Hamiltonian simulation of $\dfrac{\theta}{2}\sigma_{0}^z\sigma_{1}^z\sigma_{2}^z$.

Hamiltonian simulation of the components in the cost Hamlitonian $\hat{H}_C$.

§.§ Example

We now demonstrate how all described steps of transforming the objective function into the corresponding QAOA circuit can be done in practice using the following example:

\begin{align}
f:\left]-2^2,2^2\right[&\rightarrow \mathbb{R}\\
x&\mapsto x^2 + 2x
\end{align}

Choosing a zero bit resolution $m=0$ for simplicity, the bit encoding is displayed in the following map:

\begin{align}
x\mapsto \left(2x_0-1\right)\left(2^1 x_1 + 2^0 x_2\right).
\end{align}

Therefore, $f$ can now be written in discretized form as follows:

\begin{align}
f\left(x_0,x_1,x_2\right)= &\left(\left(2x_0-1\right)\left(2^1 x_1 + 2^0 x_2\right)\right)^2\nonumber \\
&+ 2\left(2x_0-1\right)\left(2^1 x_1 + 2^0 x_2\right) \nonumber \\
=& 4\left(4x_0x_1+x_0x_2+x_1x_2\right)
\end{align}

This then translates to the spin configuration function $f'$ as described in equation <ref> below.

\begin{align}
f'\left(s_0,s_1,s_2\right)=& 4\left(4\dfrac{s_0+1}{2}\dfrac{s_1+1}{2}+\dfrac{s_0+1}{2}\dfrac{s_2+1}{2}+\dfrac{s_1+1}{2}\dfrac{s_2+1}{2}\right) \nonumber \\
\label{eq:spin-conf-fun}
\end{align}

Using the mapping from a spin configuration function to a quantum Hamiltonian as described in section <ref>, we get:

\begin{align}
\hat{H}_C=4\left(\sigma_{0}^z\sigma_{1}^z+\sigma_{0}^z\sigma_{2}^z+\sigma_{1}^z\sigma_{2}^z+2\sigma_{0}^z+2\sigma_{1}^z+2\sigma_{2}^z+3I^{\otimes 3}\right)
\end{align}

Subsequently, we can use the combination of CNOT gates wrapping a parameterized rotation gate $R_z(\theta)$ applied on the target qubit to construct the circuit simulating the Hamiltonian $\hat{H}_C$, as indicated figure <ref>.

[column sep=12pt, row sep=20pt,between origins]
$\ket{0}$     H[3,steps=1,style=dashed, rounded corners,fill=blue!20, inner xsep=2pt, background,label style=label position=below,anchor= north,yshift=-0.2cm] State prep.     R_z(16γ_1)[3,steps=10,style=dashed, rounded corners,fill=blue!20, inner xsep=2pt, background,label style=label position=above,anchor= north,yshift=0.5cm]Parameterized Hamiltonian simulation of $\hat{H}_C$     1         1     2         2                 R_x(2β_1)[3,steps=1,style=dashed, rounded corners,fill=blue!20, inner xsep=2pt, background,label style=label position=below,anchor= north,yshift=-0.2cm]Parameterized Hamiltonian simulation of $\hat{H}_M$      …

$\ket{0}$     H     R_z(16γ_1)          R_z(8γ_1)                      1         1      R_x(2β_1)      …

$\ket{0}$     H     R_z(16γ_1)                      R_z(8γ_1)               R_z(8γ_1)         R_x(2β_1)      …
QAOA circuit implementation using single-qubit and CNOT-gates for the example in section <ref> showing $P = 1$ iterations.

§ EVALUATION

To compare the performance of the QAOA for PUBO and QUBO formulations of higher order continuous optimization functions, we run experiments on two established benchmark functions (see figures <ref> and <ref>): The 1-Dimensional Styblinski-Tang function $s(x)=\nicefrac{\left(x^4-16x^2+5x\right)}{2}$ (denoted as 1D-ST) <cit.>, and the 2-Dimensional Rosenbrock function $r(x,y)=100\left(y-x^2\right)^2 +\left(x-1\right)^2$ (denoted as 2D-Rb) <cit.>. These functions where chosen for their different requirements in terms of the number of needed qubits to model them (for details see figure <ref>) and their hardness[According to the results from Global Optimization Benchmarks and AMPGO by Andrea Gavana, see <http://infinity77.net/global_optimization/index.html>]. Having to specify input domain spaces in which the search for the optimal value is to be conducted, we choose the interval $\left]-4,4\right[$ for the 1D-ST function and $\left]-4,4\right[^2$ for the 2D-Rb function. These domain spaces allow us to find the  global optimum of each function and enable us to investigate many different bit resolutions while staying within reasonable simulation times of a couple of hours. More specifically, these input domains allow exploring bit resolutions of 0 to 3 for the 1D-ST function and 0 to 1 for the 2D-Rb function.

The 1D Styblinski-Tang function for a bit resolution of one.

The 2D Rosenbrock function for a bit resolution of zero.

Visualizations of benchmark functions used for the evaluation.

In the following, we explore the performance differences between the PUBO and QUBO approaches in terms of three criteria:

* The solution quality

* The parameter training

* The circuit width and depth

For all of the following experiments, we used Qiskit's qasm simulator, the COBYLA optimizer because of its short runtime, and 1024 shots as a standard for all circuit runs. In addition to that, we initialized all parameters using ramp initialization, as it consistently showed the best results in our experiments. Notably, the ramp initialization simply corresponds to the choosing equidistantly spaced intervals for the discretized Hamiltonian simulation described in section <ref>. Furthermore, we conducted our studies for a very high number of QAOA iterations compared to related work, i.e., $1\leq P \leq 40$, as this allows for a better performance estimation in terms of scaling.

§.§ Solution quality

To evaluate the solution quality of both approaches (PUBO and QUBO), we now examine their performance at different bit resolutions and varying QAOA iterations $P$ as exemplified in figure <ref>.

For figure <ref>, we chose to display a baseline result, i.e., the 1D-ST function at zero bit resolution, as this function is a QUBO problem by nature. With both plots showing very similar behavior, it becomes apparent that PUBO performs completely analogously to the QUBO for quadratic functions.

Examining figures <ref> and <ref>, we can see that the PUBO approach consistently outperforms the QUBO approach for higher order functions, as the expected value, the median and the overall variance are significantly lower for PUBO. This becomes increasingly apparent for the harder 2D-Rb function displayed in figure <ref>, as we can see that the QUBO approach essentially plateaus for increasing $P$, while the PUBO performance clearly benefits from higher $P$. These results are especially promising if this trend does continue for higher $P$, which is to be explored in future work.

While these plots merely display exemplified results, our full evaluation results clearly substantiate the trends visible in the selected plots.

Results for the 1D-ST function with a bit resolution of zero.

Results for the 1D-ST function with a bit resolution of one.

Results for the 2D-Rb with a bit resolution of zero.

Box plots showing the quality of the solutions found using the QAOA for the QUBO (blue) and PUBO (red) approaches for different numbers of QAOA iterations $P$. The seeked global minimum for the Rosenbrock function is $0$ and $-39.16599$ for the Styblinksi-Tang.

§.§ Parameter Training

Following the recommendation of <cit.>, we select the COBYLA optimizer to train the QAOA parameters. This optimizer has a built in stopping criterion, terminating the learning process when the last couple optimization iterations did not increase the objective value above a specific threshold (in our case $1e-4$). To prevent this procedure from exceeding a reasonable execution time, the user can also specify a number of maximum possible iterations. We use this functionality by capping the number of optimization steps at 1000, relying on results from preliminary experiments that showed, that almost no problem instances exceeded this number of optimization iterations. This allows us to compare the number of optimization steps between the PUBO and QUBO approaches unimpaired of this hyperparameter, as almost all parameter trainings run until completion.

Examining the number of optimization steps for different $P$ shown in figure <ref>, it becomes clear that both approaches need roughly the same number of optimization steps. In general, we can also observe that the number of optimization steps for the Styblinski-Tang function is generally higher compared to the Rosenbrock function. We suspect this being caused by the flatter landscape of the Rosenbrock function leading to below-threshold training improvements sooner. In addition to that, we can observe that the QUBO approach has a tendency to decrease its ascend in training time earlier when the solution quality is worse than the PUBO (which is the case for the 1D-ST function at a bit resolution of three, as this function has a very similar plot to the one displayed in figure <ref>).

1D-ST function at a bit resolution of zero.

1D-ST function at a bit resolution of three.

2D-Rb function at a bit resolution of one.

Number of parameter training iterations for different numbers of QAOA iterations $P$.

When simulating quantum circuits using classical hardware, execution times play an important role, as they limit what can be learned about their properties such as scaling behavior using non-quantum hardware. As displayed in figure <ref>, the training time does not differ significantly if the function only has a small amount of higher order terms involved (see figure <ref>), while execution time increases massively for the QUBO approach the more qubits are needed and the more higher order terms appear. Notably, that difference is mostly dominated by the number of qubits involved (13 for the QUBO formulation of 1D-ST at a bit resolution of three versus the 17 qubits needed for the QUBO formulation of the 2D-Rb function at a bit resolution of 1). This clearly demonstrates the performance of PUBO for simulation on classical hardware, also allowing for a deeper scaling analyses, which are very valuable in practice.

2D-Rb function at a bit resolution of zero.

1D-ST, at a bit resolution of three.

2D-Rb, at a bit resolution of one.

Wall-clock training time of the QAOA parameters for different numbers of QAOA iterations $P$.

§.§ Circuit width and depth

For the execution of the proposed approaches on real hardware, two criteria are essential: the circuit width (i.e., the number of qubits) and the circuit depth (i.e., the number of subsequent gate operations). Figure <ref> exemplifies the both using the 1D-ST function, as it allows for a bigger scaling analysis in terms of bit resolution.

Before comparing the number of needed qubits for both approaches, we recall that the number of required qubits is entirely determined by the bit depth and the dimensions of the input domain, according to our chosen discretization. For PUBO, we can easily calculate the number of required qubits by adding up the number of bits used to represent each dimension of the input domain. For QUBO we can calculate this number by determining the number of required ancillary qubits and adding it to the number of qubits required for the PUBO formulation as a result of the quadratization. The number of ancillary qubits however relies heavily on the exact function and the techniques used for quadratization. We used a combination of different techniques based on the python package qubovert[<https://github.com/jtiosue/qubovert>] and boolean algebra simplifications. The resulting number of qubits for the 1D-ST function are displayed in figure <ref>. Comparing the PUBO and QUBO approaches, we can clearly see a higher number of needed qubits in the QUBO variant, which gradually increases with bit resolution, as the number of qubic and quartic terms accumulate according to the chosen discretization, as described in section <ref>.

Continuing with the circuit depth (also displayed in figure <ref>), we can observe a clear disadvantage of the PUBO approach when executed on a device that doesn't inherit a suitable gate set: While the QUBO's overall circuit depth at the highest complexities caps at less than 1400, the PUBO's overall depth reaches around 4000. The substantially higher circuit depth for current hardware raises an important potential drawback when deciding on whether to incorporate PUBO on current NISQ devices, where gate-fidelity is a significant constraint. For future quantum computers implementing suitable multi-qubit gates however, the scaling in terms of circuit depth would roughly equal that of the QUBO approach. Possibly, even less gates might be needed, as no interactions with ancillary qubits are needed. Another promising observation is the easier use of classical circuit simulators in PUBO, as they can generally provide arbitrary gate sets and thus allow for even shorter circuits while also needing fewer qubits.

1D-ST, QUBO

1D-ST, PUBO

Number of required qubits and circuit depth for different numbers of QAOA iterations $P$ and bit resolutions. The circuit depth is calculated after a transpilation targeted towards a gate set without multi-qubit gates beyond CNOTs.

§ CONCLUSION

The conducted experiments clearly indicate that PUBO formulations achieve superior result quality over their quadratized QUBO analogues for continuous polynomial objective functions of a higher order. Until suitable multi-qubit gates become available, this manifests in a trade-off between the number of needed qubits (linearly higher for QUBO) and the circuit depth (linearly higher for PUBO). In terms of parameter training steps, both approaches performed equally. When using a quantum circuit simulator however, the wall-clock times for the PUBO formulations showed much better results, most probably because of the lower number of qubits that need to be simulated. For NISQ hardware, the performance difference is still mostly unclear and should be investigated in future work. We expect a strong dependence on the objective function, the input domain and bit resolution as well as their interplay with the error rates to be decisive. Finally, in the future we plan on exploring the combination of our findings with the existing positive results on using PUBO for combinatorial optimization problems to investigate the performance of PUBO formulations for NP-hard mixed integer problems.

This work was partially funded by the German BMWK project QCHALLenge (01MQ22008A). The authors want to thank Johannes Kolb for his contributions to this research.

#1     #1#1#1     #1     #1     #1     #1#1          #1#1

[Albash and LidarAlbash and Lidar2018]
authorpersonTameem Albash and
personDaniel A. Lidar. year2018.
Demonstration of a Scaling Advantage for a Quantum
Annealer over Simulated Annealing.
journalPhys. Rev. X  volume8
(dateJul year2018), pages031016.
Issue 3.

[Born and FockBorn and Fock1928]
authorpersonM Born and personV Fock.
Beweis des Adiabatensatzes.
journalZeitschrift für Phys.
volume51, number3 (year1928),

[Boros and HammerBoros and Hammer2002]
authorpersonEndre Boros and personPeter L.
Hammer. year2002.
Pseudo-Boolean optimization.
journalDiscrete Applied Mathematics
volume123, number1 (year2002),

[Denchev, Boixo, Isakov, Ding, Babbush,
Smelyanskiy, Martinis, and NevenDenchev et al2016]
authorpersonVasil S. Denchev, personSergio
Boixo, personSergei V. Isakov, personNan Ding,
personRyan Babbush, personVadim Smelyanskiy,
personJohn Martinis, and personHartmut Neven.
What is the Computational Value of Finite-Range
journalPhys. Rev. X  volume6
(dateAug year2016), pages031015.
Issue 3.

[Ebadi, Keesling, Cain, Wang, Levine,
Bluvstein, Semeghini, Omran, Liu, Samajdar, Luo, Nash, Gao, Barak, Farhi,
Sachdev, Gemelke, Zhou, Choi, Pichler, Wang, Greiner, Vuletić, and
LukinEbadi et al2022]
authorpersonS. Ebadi, personA. Keesling,
personM. Cain, personT. T. Wang, personH.
Levine, personD. Bluvstein, personG. Semeghini,
personA. Omran, personJ.-G. Liu, personR.
Samajdar, personX.-Z. Luo, personB. Nash,
personX. Gao, personB. Barak, personE.
Farhi, personS. Sachdev, personN. Gemelke,
personL. Zhou, personS. Choi, personH.
Pichler, personS.-T. Wang, personM. Greiner,
personV. Vuletić, and personM. D. Lukin.
Quantum optimization of maximum independent set
using Rydberg atom arrays.
journalScience volume376,
number6598 (year2022), pages1209–1215.

[Farhi, Goldstone, and GutmannFarhi
et al2014]
authorpersonEdward Farhi, personJeffrey
Goldstone, and personSam Gutmann.
titleA Quantum Approximate Optimization Algorithm.
[arxiv]1411.4028 [quant-ph]

[Floudas and LinFloudas and Lin2005]
authorpersonChristodoulos A Floudas and
personXiaoxia Lin. year2005.
Mixed Integer Linear Programming in Process
Scheduling: Modeling, Algorithms, and Applications.
journalAnn. Oper. Res. volume139,
number1 (year2005), pages131–162.

[Glos, Krawiec, and ZimborásGlos
et al2022]
authorpersonAdam Glos, personAleksandra
Krawiec, and personZoltán Zimborás.
Space-efficient binary optimization for
variational quantum computing.
journalnpj Quantum Inf. volume8,
number1 (year2022), pages39.

[Maslov and NamMaslov and Nam2018]
authorpersonDmitri Maslov and
personYunseong Nam. year2018.
Use of global interactions in efficient quantum
circuit constructions.
journalNew Journal of Physics
volume20, number3 (datemar
year2018), pages033018.

[Mitarai, Negoro, Kitagawa, and FujiiMitarai
et al2018]
authorpersonK. Mitarai, personM. Negoro,
personM. Kitagawa, and personK. Fujii.
Quantum circuit learning.
journalPhysical Review A volume98,
number3 (datesep year2018).

[Nielsen and ChuangNielsen and
authorpersonMichael A. Nielsen and
personIsaac L. Chuang. year2010.
booktitleQuantum Computation and Quantum
Information: 10th Anniversary Edition.
publisherCambridge University Press.

authorpersonM. J. D. Powell.
booktitleA Direct Search Optimization Method That
Models the Objective and Constraint Functions by Linear Interpolation.
publisherSpringer Netherlands,
addressDordrecht, pages51–67.

authorpersonH. H. Rosenbrock.
An Automatic Method for Finding the Greatest or
Least Value of a Function.
journalComput. J. volume3,
number3 (date01 year1960),

[Salehi, Glos, and MiszczakSalehi
et al2022]
authorpersonÖzlem Salehi, personAdam
Glos, and personJarosław Adam Miszczak.
Unconstrained Binary Models of the Travelling
Salesman Problem Variants for Quantum Optimization.
journalQuantum Information Processing
volume21, number2 (datefeb
year2022), numpages30 pages.

[Styblinski and TangStyblinski and
authorpersonM.A. Styblinski and
personT.-S. Tang. year1990.
Experiments in nonconvex optimization: Stochastic
approximation with function smoothing and simulated annealing.
journalNeural Networks volume3,
number4 (year1990), pages467–483.

[Tabi, El-Safty, Kallus, Hága, Kozsik, Glos,
and ZimborásTabi et al2020]
authorpersonZsolt Tabi, personKareem H.
El-Safty, personZsófia Kallus, personPéter Hága,
personTamás Kozsik, personAdam Glos, and
personZoltán Zimborás. year2020.
Quantum Optimization for the Graph Coloring Problem
with Space-Efficient Embedding. In booktitle2020 IEEE
International Conference on Quantum Computing and Engineering (QCE).

authorpersonThe Qiskit Team.
titleSimulating molecules using VQE.


2021

These authors contributed equally to this work.

These authors contributed equally to this work.

[1,4]Carmen Molina-París

[1]Martín López-García

1]Department of Applied Mathematics, School of Mathematics, University of
Leeds, Leeds, LS2 9JT, UK

2]Research $\&$ Development, Pharmaceuticals Systems Pharmacology $\&$
Medicine, Bayer AG, Leverkusen, 51368, Germany

3]Instituto de Investigación Tecnológica (IIT), Universidad Pontificia
Comillas, Madrid, Spain

4]T-6 Theoretical Biology and Biophysics, Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM, USA

5]Department of Mathematics, University of Benghazi, Benghazi, Libya

# A stochastic model of cell proliferation and death across a sequence of
compartments

Hanan Dreiwi<EMAIL_ADDRESS>Flavia Feliciangeli<EMAIL_ADDRESS>Mario
Castro<EMAIL_ADDRESS>Grant Lythe<EMAIL_ADDRESS>molina-
<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [  [  [  [

###### Abstract

Cells of the human body have nearly identical genome but exhibit very
different phenotypes that allow them to carry out specific functions and react
to changes in their surrounding environment. This division of labour is
achieved by cellular division and cellular differentiation, events which lead
to a population of cells with unique characteristics. In this paper, we model
the dynamics of cells over time across a sequence of compartments. Cells
within a compartment may represent being at the same spatial location or
sharing the same phenotype. In this sequence of compartments, cells can either
die, divide or enter an adjacent compartment. We analyse a set of ordinary
differential equations to describe the evolution of the average number of
cells in each compartment over time. We also focus on the progeny of a founder
cell in terms of a stochastic process and analyse several summary statistics
to bring insights into the lifespan of a single cell, the number of divisions
during its lifespan, and the probability to die in each compartment. Numerical
results inspired by cellular immune processes allow us to illustrate the
applicability of our techniques.

###### keywords:

Cellular dynamics, stochastic process, progeny, single-cell analysis,
symmetric and asymmetric division, self-renewal

##  1 Introduction

Cellular proliferation, death and differentiation are crucial and ongoing
processes during all stages of life. The most famous example of cell
differentiation is embryonic development, where unspecialised cells give rise
to properly functional and differentiated cells Martin1981 ; Evans1981 .
However, this progressive sequence of cell type changes takes place not only
in embryo cells but also in adult organisms. Cells gradually acquire a
different structure, and this process leads to the emergence of a distinct
cell population with specific functions. Some examples are skin cells
zhang2015modelling , T cell development in the thymus Sawicka2014 , or the T
cell exhaustion dysfunctional process Wherry2015 ; Chen2019 . The general idea
is that starting from a “precursor” or “stem-cell” population, cells can
differentiate, and thus, change their “state” (phenotype) to produce a stable
population of differentiated and functional specific cells.

Recent techniques in experimental cell biology have made it possible to track
individual cell states and the progeny of a specific cell Hodzic2016 ;
Perie2014 . However, we are still far from fully understanding the cellular
heterogeneity arising from the development and cell differentiation processes.
Thus, mathematical and computational approaches can help answer some open
questions regarding the underlying mechanisms of cellular fate. During the
last decades, the development of mathematical models in Biology has
significantly increased the ability to gain a better quantitative
understanding of the dynamics of cell populations over time. Deterministic
models, which do not incorporate randomness and are typically easier to
analyse, allow one to describe the dynamics over time of many entities, such
as animals or cells. However, when looking at a smaller scale, such as when
tracking single-cell behaviour, stochastic fluctuations arising from complex
cellular interactions and the molecular events which regulate cellular fate
must be considered.

Within the area of stochastic analysis, the theory of branching processes has
been widely used in cellular dynamics Kimmel2002 . For example, the classic
Galton-Watson model, which was originally developed to study the extinction of
family surnames watson1875probability  has been successfully applied to
cellular dynamics Kimmel2002 . The theory of branching processes can answer
questions related to the limiting behaviour (e.g., probability of extinction
versus unlimited cellular growth) of cellular populations. A natural
generalization is the multi-type branching process, where individuals are not
all of the same type. These models can be effectively used to represent cells
changing their spatial location de2019fate , or their phenotype Nordon20117
over time. These approaches model cellular dynamics across compartments, where
cells in the same compartment are assumed to be homogeneous and behave
identically. We can highlight here the seminal work by Matis Matis1971 ,
proposing a stochastic compartmental model for cellular dynamics.
Deterministic compartmental models are widely used in pharmacodynamics, as
well as in Mathematical Ecology Matis1994 , Immunology Eftimie2016  and
Epidemiology Capasso2020 .

Here we consider a sequence of compartments and model the process by which,
from a stem-cell pool (in compartment $C_{1}$), cells can undergo
differentiation and migration events across adjacent compartments, potentially
leading to the rise of a differentiated terminal progeny population (in
compartment $C_{N}$). Cells of each compartment can divide, die or transit to
adjacent compartments (e.g., representing potentially reversible phenotype
change). In our general stochastic model, differentiation can be either
considered to be irreversible (e.g., cancer cell mutations) or arising from a
mutation event that has a negligible probability of being reversed (e.g.,
embryonic cell program). On the other hand, cells are short-lived compared to
the host, so a continuous regeneration process is required to maintain cell
populations. Recent advances in genetic labelling Klein2011  have confirmed
the existence of different types of division, which play an essential role in
the ability and flexibility of a cell pool to expand (or contract) when
required. Different division events can drive different biological processes
according to whether the less differentiated pool where the division happens
expands to one new cell (self-renewal), stays of the same size (asymmetric
division, where one of the daughter cells changes phenotype), or shrinks
(symmetric division, where both daughter cells change phenotype) Barile2020 .
The inclusion of symmetric and asymmetric cell division allows one to analyse
how a population of cells needs to be flexible in expanding and how epigenetic
information is established and inherited in the cells of a complex and multi-
cellular organism.

In Section 2 we describe the dynamics of cells dividing, dying or exiting
across a sequence of compartments over time, in terms of a continuous-time
Markov chain. The mean behaviour of the system is analysed in Section 2.1. In
Section 2.2 we focus on studying the proliferative potential of the system at
the population level by quantifying the number of cells within the genealogy
of a single progenitor cell. In both sections, we consider either the
situation where differentiation events can be reversible, or a mathematically
more simple “irreversible model”, where differentiation to the next
compartment cannot be reversed. In Section 2.3 our focus is on investigating
summary statistics related to a single cell which we track over time. We
conclude, in Section 3, by proposing a set of numerical results inspired by
cellular immune processes in order to illustrate our approach.

##  2 Stochastic model

We propose a stochastic model of cell proliferation, death and differentiation
(or migration) across an ordered sequence of compartments. Cells in a given
compartment may represent co-location in the same spatial compartment, or
belong to a common phenotype, or share some common characteristics. We
consider a sequence of compartments $C_{i}$, $i\in\\{1,\dots,N\\}$, which
cells follow, behaving independently of each other as typically considered
within the theory of continuous-time branching processes Kimmel2002 .

We start with a general stochastic model in this section, and consider
cellular events, inspired by some recent mathematical models Barile2020 ;
Sawicka2014 ; zhang2015modelling ; de2019fate . We assume that each of these
events takes place at a certain per cell rate. Particular situations of
interest arise from setting some of these rates equal to zero, so that some
events are not allowed to happen, as we illustrate in our case studies in
Section 3. Each cell in a given compartment, $C_{i}$, can divide, die or exit
to one of the two adjacent compartments. When a division occurs, daughter
cells might both belong to the same compartment as the mother (typically
referred to as self-renewal), both daughter cells might both instantaneously
belong to the next compartment (symmetric division), or one daughter cell
might belong to the same compartment as the mother while the other belongs to
the next compartment (asymmetric division) Barile2020 ; zhang2015modelling .

Figure 1: Schematic representation of the general stochastic model for cell
proliferation, death and differentiation across an ordered sequence of
compartments.

Our stochastic model, depicted in Figure 1, can be described in terms of the
continuous-time Markov chain (CTMC) ${\cal
X}=\\{({\bm{C}}_{1}(t),{\bm{C}}_{2}(t),\dots,{\bm{C}}_{N}(t)):\ t\geq 0\\}$,
where ${\bm{C}}_{i}(t)$ represents the number of cells in compartment $C_{i}$
at time $t$, with state space given by ${\cal
S}=\\{0,1,2,\dots\\}^{N}=\mathbb{N}_{0}^{N}$. Cellular events, labelled E1 to
E5, represent transitions across states $(n_{1},\dots,n_{N})\in{\cal S}$ as
follows:

  * (E1)

Self-renewal (cellular division where both daughter cells remain in the same
compartment as the mother) can occur in any compartment $C_{i}$, with per cell
rate $\lambda_{i}$, for $i\in\\{1,\dots,N\\}$,

$\displaystyle(n_{1},\dots,n_{i-1},n_{i},n_{i+1},\dots,n_{N})$
$\displaystyle\xrightarrow{\lambda_{i}n_{i}}$
$\displaystyle(n_{1},\dots,n_{i-1},n_{i}+1,n_{i+1},\dots,n_{N}).$

  * (E2)

Asymmetric division (cellular division where one of the daughter cells remains
in the same compartment as the mother, while the other goes to the next
compartment) occurs in compartment $C_{i}$ with per cell rate $a_{i}$, for
$i\in\\{1,\dots,N\\}$,

$\displaystyle(n_{1},\dots,n_{i-1},n_{i},n_{i+1},\dots,n_{N})$
$\displaystyle\xrightarrow{a_{i}n_{i}}$
$\displaystyle(n_{1},\dots,n_{i-1},n_{i},n_{i+1}+1,\dots,n_{N}).$

  * (E3)

Symmetric division (cellular division where both daughter cells
instantaneously move to the next compartment) occurs in compartment $C_{i}$
with per cell rate $s_{i}$, for $i\in\\{1,\dots,N\\}$,

$\displaystyle(n_{1},\dots,n_{i-1},n_{i},n_{i+1},\dots,n_{N})$
$\displaystyle\xrightarrow{s_{i}n_{i}}$
$\displaystyle(n_{1},\dots,n_{i-1},n_{i}-1,n_{i+1}+2,\dots,n_{N}).$

  * (E4)

Differentiation (or migration) between adjacent compartments can occur with
rates $(\nu_{i},\xi_{i})$, for $i\in\\{1,\dots,N\\}$,

$\displaystyle(n_{1},\dots,n_{i-1},n_{i},n_{i+1},\dots,n_{N})$
$\displaystyle\xrightarrow{\nu_{i}n_{i}}$
$\displaystyle(n_{1},\dots,n_{i-1},n_{i}-1,n_{i+1}+1,\dots,n_{N}),$
$\displaystyle(n_{1},\dots,n_{i-1},n_{i},n_{i+1},\dots,n_{N})$
$\displaystyle\xrightarrow{\xi_{i}n_{i}}$
$\displaystyle(n_{1},\dots,n_{i-1}+1,n_{i}-1,n_{i+1},\dots,n_{N}),$

with $\nu_{N}=\xi_{1}=0$.

  * (E5)

Cells can die in any compartment $C_{i}$ with per cell rate $\mu_{i}$,
$i\in\\{1,\dots,N\\}$,

$\displaystyle(n_{1},\dots,n_{i-1},n_{i},n_{i+1},\dots,n_{N})$
$\displaystyle\xrightarrow{\mu_{i}n_{i}}$
$\displaystyle(n_{1},\dots,n_{i-1},n_{i}-1,n_{i+1},\dots,n_{N}).$

We note that cells in the last compartment $C_{N}$ cannot symmetrically or
asymmetrically divide, or differentiate.

###  2.1 Number of cells in each compartment

We first describe the dynamics of the process by focusing on the time course
of the average numbers of cells, $\mathbb{E}[{\bm{C}}_{i}(t)]$, that obey the
following system of differential equations Matis1971

$\displaystyle\frac{d\,\mathbb{E}[{\bm{C}}_{1}(t)]}{dt}$ $\displaystyle=$
$\displaystyle-(\mu_{1}+\nu_{1}+s_{1}-\lambda_{1})\mathbb{E}[{\bm{C}}_{1}(t)]+\xi_{2}\mathbb{E}[{\bm{C}}_{2}(t)],$
$\displaystyle\frac{d\,\mathbb{E}[{\bm{C}}_{i}(t)]}{dt}$ $\displaystyle=$
$\displaystyle(\nu_{i-1}+a_{i-1}+2s_{i-1})\,\mathbb{E}[{\bm{C}}_{i-1}(t)]-(\mu_{i}+\nu_{i}+\xi_{i}+s_{i}-\lambda_{i})\,\mathbb{E}[{\bm{C}}_{i}(t)]$
(1) $\displaystyle+$
$\displaystyle\xi_{i+1}\,\mathbb{E}[{\bm{C}}_{i+1}(t)],\,\,\,\,\
\quad\quad\quad\quad\quad\quad i\in\\{2,...,N-1\\},$
$\displaystyle\frac{d\,\mathbb{E}[{\bm{C}}_{N}(t)]}{dt}$ $\displaystyle=$
$\displaystyle(\nu_{N-1}+a_{N-1}+2s_{N-1})\,\mathbb{E}[{\bm{C}}_{N-1}(t)]-(\mu_{N}+\xi_{N}-\lambda_{N})\mathbb{E}[{\bm{C}}_{N}(t)],$

where $\mathbb{E}[{\bm{C}}_{i}(t)]$ represents the expectation of the random
variable ${\bm{C}}_{i}(t)$. These equations constitute a homogeneous first-
order linear system of ODEs with constant coefficients, which can be written
more succinctly in matrix form as follows

$\frac{d{\bm{\Phi}}}{dt}={\bf A}{\bm{\Phi}},$ (2)

where

${\bm{\Phi}}=\left(\begin{array}[]{c}\mathbb{E}[{\bm{C}}_{1}(t)]\\\
\mathbb{E}[{\bm{C}}_{2}(t)]\\\ \vdots\\\ \mathbb{E}[{\bm{C}}_{N-1}(t)]\\\
\mathbb{E}[{\bm{C}}_{N}(t)]\end{array}\right),\quad{\bf
A}=\left(\begin{array}[]{ccccc}-\Delta_{1}&\xi_{2}&0&\cdots&0\\\
\Lambda_{1}&-\Delta_{2}&\xi_{3}&\dots&\vdots\\\ 0&\ddots&\ddots&\ddots&0\\\
\vdots&\vdots&\Lambda_{N-2}&-\Delta_{N-1}&\xi_{N}\\\
0&\cdots&0&\Lambda_{N-1}&-\Delta_{N}\end{array}\right),$ (3)

with

$\displaystyle\Delta_{1}$ $\displaystyle=$
$\displaystyle\mu_{1}+\nu_{1}+s_{1}-\lambda_{1},$ $\displaystyle\Delta_{i}$
$\displaystyle=$
$\displaystyle\mu_{i}+\nu_{i}+s_{i}+\xi_{i}-\lambda_{i},\quad\Lambda_{i-1}\ =\
\nu_{i-1}+a_{i-1}+2s_{i-1},\ i\in\\{2,\dots,N-1\\},$ $\displaystyle\Delta_{N}$
$\displaystyle=$ $\displaystyle\mu_{N}+\xi_{N}-\lambda_{N},\quad\Lambda_{N-1}\
=\ \nu_{N-1}+a_{N-1}+2s_{N-1}.$

The initial value problem in (2) with ${\bf\Phi}_{0}={\bf\Phi}(0)$ has a
unique solution (allen2007introduction, , Theorem 4.1)

${\bm{\Phi}}(t)=e^{{\bf A}t}{\bm{\Phi}}_{0},$ (4)

where $e^{{\bf A}t}$ represents the matrix exponential,

$e^{{\bf A}t}={\bf I}+{\bf A}t+{\bf A}^{2}\frac{t^{2}}{2!}+{\bf
A}^{3}\frac{t^{3}}{3!}+...=\sum_{i=0}^{+\infty}\frac{({\bf A}t)^{i}}{i!}.$

The system of equations (2) admits
$\lim_{t\rightarrow+\infty}{\bm{\Phi}}(t)={\bf 0}_{N}$ (column vector of
zeros) as an asymptotic solution; this is exponential stable (all solutions of
the system from any initial conditions converge to 0) if and only if each
eigenvalue of ${\bf A}$ has a negative real part Callier1991linear .

Figure 2: Schematic diagram of the general compartmental model without reverse
differentiation events.

We note that, in many biological applications, some of the parameters in Fig.
1 will be zero, and thus, the analysis of such systems will be much simpler.
For instance, when differentiation is irreversible Sawicka2014 ; Barile2020 ,
so that $\xi_{i}=0$ for $i\in\\{2,\dots,N\\}$ as shown in Figure 2. We will
refer to this situation as the irreversible model. In this case, and if one
focuses on a single cell starting in compartment $C_{1}$ at time $t=0$,
${\bm{\Phi}}(0)=(1,0,\dots,0)^{T}$, it is possible to get the average
behaviour of the system in terms of the following explicit expressions

$\displaystyle\mathbb{E}[{\bm{C}}_{i}(t)]$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}e^{-\Delta_{i}t},&i=1,\\\
\left(\prod\limits_{l=1}^{i-1}\Lambda_{l}\right)\sum\limits_{j=1}^{i}e^{-\Delta_{j}t}\prod\limits_{\begin{subarray}{c}m=1\\\
m\neq
j\end{subarray}}^{i}(\Delta_{m}-\Delta_{j})^{-1},&i\in\\{2,\dots,N\\}.\end{array}\right.$
(7)

The expression above is only well-defined if $\Delta_{i}\neq\Delta_{j}$ for
all pairs $(i,j)$. If this is not the case, alternative similar explicit
solutions can be found. For example, when $\Delta_{i}=\Delta_{j}$ for all
$i,j\in\\{2,\dots,N\\}$, then the formula above simplifies to

$\mathbb{E}[{\bm{C}}_{i}(t)]=\left(\prod\limits_{l=1}^{i-1}\Lambda_{l}\right)\frac{t^{i}}{i!}e^{-\Delta_{i}t}.$
(8)

From the expressions above, it is clear that in the irreversible model
$\lim_{t\rightarrow+\infty}\mathbb{E}[{\bm{C}}_{i}(t)]=0$ when $\Delta_{i}>0$,
$\forall i\in\\{1,\dots,N\\}$. This is consistent since the parameters
$-\Delta_{i}$ are the eigenvalues of ${\bf A}$ in this case. For the
applications we have in mind Sawicka2014 ; zhang2015modelling , it is of
interest to quantify the cumulative number of cells, on average, that arrive
to the final compartment $C_{N}$ starting with a single cell in compartment
$C_{1}$. To answer this question, one can set $\lambda_{N}=\mu_{N}=0$, so that
cells which arrive to $C_{N}$ accumulate and can be counted, and then compute
in the irreversible model, from the expression above

$\lim_{t\rightarrow+\infty}\mathbb{E}[{\bm{C}}_{N}(t)]=\prod_{i=1}^{N-1}\frac{\Lambda_{i}}{\Delta_{i}}.$

A particular feature of this system is that cells behave independently from
each other. This means that the dynamics of the genealogy of a set of $M$
precursor (or progenitor cells) in compartment $C_{1}$ at time $t=0$, can be
analysed as $M$ independent stochastic processes. Thus, in Section 2.2 we
focus on analysing a number of summary statistics of interest which are
related to the genealogy of a single cell starting in a given compartment
(typically compartment $C_{1}$).

###  2.2 The genealogy of a single precursor cell

We consider a cell starting in a given compartment $C_{i}$, and we define
$G_{i}$, the random variable representing the total number of cells in the
genealogy of this cell. Cells in the genealogy are the daughters,
granddaughters, …, of the precursor cell, which arise from division events
(either self-renewal, asymmetric or symmetric) in any compartment over time,
not including the precursor cell. It is a summary statistic of the process
which quantifies the proliferative potential of a single cell. For example, in
Figure 3, we represent a particular realisation of the stochastic process,
where $G_{1}=8$.

The mean number of cells in the genealogy of a precursor cell,
$m_{i}=\mathbb{E}[G_{i}]$, for any initial compartment of interest
$i\in\\{1,\dots,N\\}$, can be obtained via first-step arguments by
conditioning on the next event that occurs in the stochastic process. This
approach leads to the following system of equations

$\displaystyle\Delta_{1}m_{1}$ $\displaystyle=$
$\displaystyle\Lambda_{1}m_{2}+2(\lambda_{1}+a_{1}+s_{1}),$
$\displaystyle\Delta_{i}m_{i}$ $\displaystyle=$
$\displaystyle\Lambda_{i}m_{i+1}+\xi_{i}m_{i-1}+2(\lambda_{i}+a_{i}+s_{i}),\quad
i\in\\{2,\dots,N-1\\},$ $\displaystyle\Delta_{N}m_{N}$ $\displaystyle=$
$\displaystyle\xi_{N}m_{N-1}+2\lambda_{N}.$

The system above can be expressed in matrix form via the column vectors ${\bf
m}=\left(m_{1},\dots,m_{N}\right)^{T}$ and ${\bf
b}=\left(2(\lambda_{1}+a_{1}+s_{1}),\dots,2(\lambda_{N-1}+a_{N-1}+s_{N-1}),\right.$
$\left.2\lambda_{N}\right)^{T}$ as follows

${\bf J}\,{\bf m}={\bf b},$ (9)

with a tri-diagonal coefficient matrix

${\bf J}=\left(\begin{array}[]{cccccc}\Delta_{1}&-\Lambda_{1}&0&0&\cdots&0\\\
-\xi_{2}&\Delta_{2}&-\Lambda_{2}&0&\cdots&0\\\
0&-\xi_{3}&\Delta_{3}&-\Lambda_{3}&\cdots&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&\cdots&0&-\xi_{N-1}&\Delta_{N-1}&-\Lambda_{N-1}\\\
0&\cdots&0&0&-\xi_{N}&\Delta_{N}\end{array}\right).$ (10)

Figure 3: A realisation of the stochastic process following the genealogy of a
single precursor cell which starts in compartment $C_{1}$. The cell tracked
(see Section 2.3) is depicted as striped, and where the colour indicates the
compartment where it is at any given time. Here, the tracked cell dies in
$C_{3}$ (brown), while its genealogy continues up to $C_{4}$. In this example,
$G_{1}=8=G_{1}(1)+G_{1}(2)+G_{1}(3)+G_{1}(4)=4+0+2+2$.

One can exploit the tri-diagonal structure of ${\bf J}$ to get an explicit or
recursive solution of this system. In particular, by following a Gaussian
forward-elimination backward-substitution approach, such as the Thomas
algorithm conte1972elementary ; thomas1949elliptic , one can obtain the
recursive equations

$\displaystyle m_{N}$ $\displaystyle=$ $\displaystyle\rho_{N},\quad m_{i}\ =\
\rho_{i}-\gamma_{i}m_{i+1},\quad i\in\\{1,\dots,N-1\\},$ (11)

where $\gamma_{1}=-\Delta^{-1}_{1}\Lambda_{1}$,
$\rho_{1}=2\Delta^{-1}_{1}(\lambda_{1}+s_{1}+a_{1})$, and

$\displaystyle\gamma_{i}$ $\displaystyle=$
$\displaystyle-\frac{\Lambda_{i}}{\Delta_{i}+\xi_{i}\gamma_{i-1}},\quad
i\in\\{2,\dots,N-1\\},$ $\displaystyle\rho_{i}$ $\displaystyle=$
$\displaystyle\frac{2(\lambda_{i}+s_{i}+a_{i})+\xi_{i}\rho_{i-1}}{\Delta_{i}+\xi_{i}\gamma_{i-1}},\quad
i\in\\{2,\dots,N\\}.$

This recursive scheme leads to the explicit solution

$\displaystyle m_{i}$ $\displaystyle=$
$\displaystyle\sum_{j=i}^{N}\,(-1)^{j-i}\rho_{j}\,\left(\prod_{l=i}^{j-1}\gamma_{l}\right),\quad
i\in\\{1,\dots,N\\},$ (12)

where $\prod_{l=i}^{i-1}\gamma_{l}=1$. A condition on the parameters arises
during the implementation of the recursive scheme,

$\displaystyle\Delta_{i}+\xi_{i}\gamma_{i-1}$ $\displaystyle>$ $\displaystyle
0,\quad i\in\\{2,\dots,N\\},$

for the average values $m_{1},\dots,m_{N}$ to be finite and non-negative, for
all $i\in\\{1,\dots,N\\}$. This ensures that the number of cells in the
genealogy of the precursor cell is finite with probability one,
$\mathbb{P}(G_{i}<+\infty)=1$, since
$m_{i}=\mathbb{E}[G_{i}]=\mathbb{E}[G_{i}|G_{i}<+\infty]\mathbb{P}(G_{i}<+\infty)+\mathbb{E}[G_{i}|G_{i}=+\infty]\mathbb{P}(G_{i}=+\infty)$.

In the irreversible case —where $\xi_{i}=0$, $i\in\\{2,\dots,N\\}$—, the
solution above simplifies to

$\displaystyle
m_{i}=2\sum^{N}_{j=i}\frac{\lambda_{j}+a_{j}+s_{j}}{\Delta_{j}}\left(\prod_{l=i}^{j-1}\frac{\Lambda_{l}}{\Delta_{l}}\right),\quad
i\in\\{1,\dots,N\\},$ (13)

where for $j=N$, we set $a_{N}=s_{N}=0$. For $i=N$, the empty product above is
equal to one, so that $m_{N}=\frac{2\lambda_{N}}{\mu_{N}-\lambda_{N}}$. In
this model the condition on the parameters for finite and non-negative
solutions becomes $\Delta_{i}>0$ for all $i\in\\{1,\dots,N\\}$. This is
consistent with direct inspection of Figure 2, where $\Delta_{i}>0$ avoids
unlimited accumulation of cells in compartment $C_{i}$.

We note that cells in the genealogy of the precursor cell belong to different
compartments, as depicted by the colours in Figure 3. For a precursor cell
starting in compartment $C_{i}$, those compartments $C_{j}$,
$j\in\\{1,2,\dots,N\\}$ with more proliferative potential will contribute more
to the value of $G_{i}$. This will mainly depend on the parameters
$(\lambda_{j},a_{j},s_{j},\mu_{j},\nu_{j},\xi_{j})$, but also on the number of
cells in the genealogy arriving into that compartment (and which can thus
divide). It is of interest then to split $G_{i}=\sum_{j=1}^{N}G_{i}(j)$ in
terms of which compartments the progeny belong to, where $G_{i}(j)$ represents
the number of cells that belong to compartment $C_{j}$ in the genealogy of the
precursor cell starting in compartment $i$. For example, for the stochastic
realization in Figure 3,
$G_{1}=G_{1}(1)+G_{1}(2)+G_{1}(3)+G_{1}(4)=4+2+0+2=8$.

One can follow similar arguments to the ones above to compute the mean
quantities $m_{i}(j)=\mathbb{E}[G_{i}(j)]$. In particular, for an initial
compartment $C_{i}$, a first-step argument yields the equations

$\displaystyle\Delta_{i}m_{i}(i)$ $\displaystyle=$ $\displaystyle
2\lambda_{i}+a_{i}+\Lambda_{i}\,m_{i+1}(i)+\xi_{i}m_{i-1}(i),$
$\displaystyle\Delta_{i}m_{i}(i+1)$ $\displaystyle=$ $\displaystyle
2s_{i}+a_{i}+\Lambda_{i}\,m_{i+1}(i+1)+\xi_{i}m_{i-1}(i+1),$
$\displaystyle\Delta_{i}m_{i}(j)$ $\displaystyle=$
$\displaystyle\Lambda_{i}m_{i+1}(j)+\xi_{i}m_{i-1}(j),\quad
j\notin\\{i,i+1\\},$

where we implicitly set $\xi_{1}=0$ and $m_{N+1}(j)=0$ for notational
convenience in the equations above. Following a similar recursive approach to
the one above, one can get for any $j\in\\{1,\dots,N\\}$,

$\displaystyle m_{N}(j)$ $\displaystyle=$ $\displaystyle\rho_{N}(j),\quad
m_{i}(j)\ =\ \rho_{i}(j)-\gamma_{i}m_{i+1}(j),\quad i\in\\{1,\dots,N-1\\},$

where $\gamma_{1}=-\Delta^{-1}_{1}\Lambda_{1}$,
$\rho_{1}(j)=\Delta^{-1}_{1}d_{(1,j)}$, and

$\displaystyle\gamma_{i}$ $\displaystyle=$
$\displaystyle-\frac{\Lambda_{i}}{\Delta_{i}+\xi_{i}\gamma_{i-1}},\quad\quad
i\in\\{2,\dots,N-1\\},$ $\displaystyle\rho_{i}(j)$ $\displaystyle=$
$\displaystyle\frac{d_{(i,j)}+\xi_{i}\rho_{i-1}(j)}{\Delta_{i}+\xi_{i}\gamma_{i-1}},\quad
i,j\in\\{2,\dots,N\\},$

with

$\displaystyle d_{(i,j)}$ $\displaystyle=$
$\displaystyle\begin{cases}2\lambda_{i}+a_{i}&\quad\text{if}\quad j=i,\\\
2s_{i}+a_{i}&\quad\text{if}\quad j=i+1,\\\
0&\quad\text{otherwise.}\end{cases}$

This recursive scheme leads to the explicit solution

$m_{i}(j)=\sum_{k=i}^{N}(-1)^{k-i}\rho_{k}(j)\,\left(\prod_{p=i}^{k-1}\gamma_{p}\right).$
(14)

This expression simplifies further in the irreversible model, where
$\xi_{i}=0$ for all $i\in\\{2,\dots,N\\}$. In this instance $m_{i}(j)=0$
whenever $i>j$, and the equations above simplify to

$\displaystyle m_{i}(i)$ $\displaystyle=$
$\displaystyle\Delta_{i}^{-1}(2\lambda_{i}+a_{i}),$ $\displaystyle m_{i}(j)$
$\displaystyle=$
$\displaystyle\Delta^{-1}_{j-1}\left(\prod_{p=i}^{j-2}\Delta^{-1}_{p}\Lambda_{p}\right)\left(d(j-1,j)+d(j,j)\Delta^{-1}_{j}\Lambda_{j-1}\right),\quad
j\in\\{i+1,...,N\\},$

for any $i\in\\{1,...,N-2\\}$. Here,
$m_{N}(N)=\frac{2\lambda_{N}}{\Delta_{N}}$.

###  2.3 Single-cell analysis

We consider here the case where no symmetric or asymmetric division occurs in
the irreversible model, $s_{i}=a_{i}=\xi_{i}=0$ for all $i\in\\{1,\dots,N\\}$.
In the previous sections, we have analysed several summary statistics related
to the population of cells over time (e.g., the genealogy of a precursor
cell). In this section, instead, we aim to investigate a number of summary
statistics related to the lifespan of a single cell, extending the single-cell
analysis originally proposed in Ref. de2019fate  which considered T cells in
the periphery.

We propose to track the dynamics of a single cell, starting in compartment
$C_{i}$, until the cell dies. In any given compartment, this cell can divide
(and, since both daughters are indistinguishable, we keep tracking one of
them), move to the next compartment or die. This dynamics can be represented
then in terms of the CTMC ${\cal Y}=\\{Y(t):\ t\geq 0\\}$ over the state space
${\cal S}=\\{C_{1},C_{2},\dots,C_{N},\phi\\}$, where $Y(t)$ represents the
state of the cell at time $t$: either being in some compartment $C_{j}$, or
having already died (state $\phi$). A schematic representation of the process
${\cal Y}$ is given in Figure 4, while a particular realisation of this
stochastic process is depicted in Figure 2, where the tracked cell is shown as
striped.

We propose to analyse the dynamics of this cellular system by focusing on the
following summary statistics:

  * •

the lifespan, $T_{i}$, of the cell, starting in compartment $C_{i}$,

$\displaystyle T_{i}$ $\displaystyle=$ $\displaystyle\inf\\{t\geq 0:\
Y(t)=\phi|Y(0)=i\\},$

which quantifies the survival potential of cells in the system depending on
their initial location;

  * •

the number of divisions carried out by the cell during its lifespan, $D_{i}$,
which quantifies the proliferation potential of cells in the system; and

  * •

the probability of this cell dying in each compartment $C_{j}$; that is,
$\beta_{i}(j)=\mathbb{P}(Y(T_{i}-\Delta t)=C_{j})$ for a small enough $\Delta
t$.

Figure 4: Schematic representation of the process ${\cal Y}$.

####  2.3.1 Lifespan of a single cell

Let $T_{i}$ be the lifespan of a single cell starting in compartment $C_{i}$,
$i\in\\{1,\dots,N\\}$,

$\displaystyle T_{i}$ $\displaystyle=$ $\displaystyle\inf\\{t\geq
0:Y(t)=\phi|Y(0)=i\\},$

and consider $\tau_{i}=\mathbb{E}[T_{i}]$, its average lifespan. By
conditioning on the next event that the cell can undergo in this process, one
can get the recursive relation

$\displaystyle(\mu_{i}+\nu_{i})\tau_{i}$ $\displaystyle=$
$\displaystyle\nu_{i}\tau_{i+1}+1,\quad i\in\\{1,\dots,N-1\\},$

with $\tau_{N}=\mu_{N}^{-1}$. These equations can be explicitly solved as
follows

$\displaystyle\tau_{i}=\sum^{N-i}_{j=0}\frac{1}{\mu_{N-j}+\nu_{N-j}}\left(\prod_{r=i}^{N-(j+1)}\frac{\nu_{r}}{\mu_{r}+\nu_{r}}\right),\quad
i\in\\{1,\dots,N-1\\},$

and $\tau_{N}=\mu_{{N}}^{-1}$, since $T_{N}\sim Exp(\mu_{N})$. A similar
approach allows one to compute the Laplace-Stieltjes transform of $T_{i}$, and
any of its order moments. For example, the second order moment of the lifespan
of a cell starting in compartment $C_{i}$ is given by

$\displaystyle\mathbb{E}[T_{i}^{2}]=\sum^{N}_{j=i}R_{j}\left(\prod_{r=i}^{j-1}\frac{\nu_{r}}{\mu_{r}+\nu_{r}}\right),\quad
i\in\\{1,\dots,N\\},$

where
$R_{i}=\frac{2(\lambda_{i}\tau_{i}+\nu_{i}\tau_{i+1}+1)}{(\mu_{i}+\nu_{i})(\mu_{i}+\nu_{i}+\lambda_{i})}$,
and $R_{N}=\frac{2}{\mu_{N}^{2}}$. We note that if the cell starts in the last
compartment $C_{N}$, $\mathbb{E}[T_{N}^{2}]=2\mu_{N}^{-2}$ since $T_{N}\sim
Exp(\mu_{N})$.

####  2.3.2 Number of divisions during the lifespan

We denote by $D_{i}$ the number of division events triggered by the tracked
cell during its lifespan, starting in compartment $C_{i}$. One can focus on
computing $\eta_{i}=\mathbb{E}[D_{i}]$, which by first-step arguments
satisfies

$\displaystyle(\lambda_{i}+\mu_{i}+\nu_{i})\eta_{i}$ $\displaystyle=$
$\displaystyle\lambda_{i}(\eta_{i}+1)+\nu_{i}\eta_{i+1},\quad
i\in\\{1,\dots,N-1\\},$ $\displaystyle(\lambda_{N}+\mu_{N})\eta_{N}$
$\displaystyle=$ $\displaystyle\lambda_{N}(\eta_{N}+1).$

These equations lead to the explicit solution

$\displaystyle\eta_{i}$ $\displaystyle=$
$\displaystyle\sum^{N-i}_{k=0}\frac{\lambda_{N-k}}{\mu_{N-k}+\nu_{N-k}}\left(\prod_{j=i}^{N-(k+1)}\frac{\nu_{j}}{\mu_{j}+\nu_{j}}\right),\quad
i\in\\{1,\dots,N-1\\}.$

We note that $\eta_{N}=\lambda_{N}\mu_{N}^{-1}$ since $D_{N}\sim
Geometric(\frac{\lambda_{N}}{\mu_{N}+\lambda_{N}})$. These proliferation
events can occur at any time during the lifespan of the cell, which is
visiting different compartments over time. Thus, one can quantify the
proliferation potential of the cell during its eventual visit to each
compartment by considering $D_{i}=\sum_{j=i}^{N}D_{i}(j)$, where $D_{i}(j)$ is
the number of divisions which occur exactly in compartment $C_{j}$ during the
lifespan of the cell, which started in compartment $C_{i}$. The average values
$\eta_{i}(j)=\mathbb{E}[D_{i}(j)]$ can be computed similarly, leading to

$\displaystyle\eta_{i}(i)$ $\displaystyle=$
$\displaystyle\frac{\lambda_{i}}{\mu_{i}+\nu_{i}},$ $\displaystyle\eta_{i}(j)$
$\displaystyle=$
$\displaystyle\frac{\lambda_{j}}{\mu_{j}+\nu_{j}}\left(\prod_{k=i}^{j-1}\frac{\nu_{k}}{\mu_{k}+\nu_{k}}\right),\quad
j\geq i+1.$

We note that the expression above is consistent with the interpretation that
$D_{i}(j)\sim Geometric(\frac{\lambda_{j}}{\lambda_{j}+\nu_{j}+\mu_{j}})$,
restricted to the arrival of the cell to compartment $C_{j}$. In general, we
note that

$\displaystyle\eta_{i}(j)$ $\displaystyle=$
$\displaystyle\mathbb{E}[D_{i}(j)]=\mathbb{E}[D_{i}(j)\ |\ \hbox{\it cell ever
visits $C_{j}$}]\mathbb{P}(\hbox{\it cell ever visits $C_{j}$})$
$\displaystyle+$ $\displaystyle\mathbb{E}[D_{i}(j)\ |\ \hbox{\it cell never
visits $C_{j}$}]\mathbb{P}(\hbox{\it cell never visits $C_{j}$}),$

and, since $\mathbb{E}[D_{i}(j)\ |\ \hbox{\it cell never visits $C_{j}$}]=0$,
the quantity $\eta_{i}(j)$ also accounts for the probability of the cell not
visiting this compartment.

####  2.3.3 Cellular fate

Finally, the cell will in at a given compartment, and one can compute the
probability of this cell dying in each compartment at the end of its lifespan,
which we denote by

$\displaystyle\beta_{i}(j)$ $\displaystyle=$
$\displaystyle\mathbb{P}(Y(T_{i}-\Delta t)=C_{j}),\quad j\geq i,$

for small enough $\Delta t$. These probabilities are given by

$\displaystyle\beta_{i}(i)$ $\displaystyle=$
$\displaystyle\frac{\mu_{i}}{\mu_{i}+\nu_{i}},\quad\beta_{i}(j)\ =\
\frac{\mu_{j}}{\mu_{j}+\nu_{j}}\prod_{k=i}^{j-1}\frac{\nu_{k}}{\mu_{k}+\nu_{k}}\quad
j\geq i+1.$

##  3 Results

We illustrate our approach with two applications. In Section 3.1, we implement
the methods from Section 2.1 and Section 2.2 to explore the impact of
asymmetric and symmetric division events on the cellular dynamics, for the
case $N=4$ compartments. We carry out some sensitivity analysis on the
probabilities of self-renewal, asymmetric and symmetric division events
happening. In Section 3.2 we focus on a particular application of our single-
cell analysis from Section 2.3 to an existing model of T cell development in
the thymus Sawicka2014 .

###  3.1 Asymmetric and symmetric division: the case of four compartments

We consider here $N=4$ compartments, where the last compartment $C_{4}$ does
not involve any death, division or differentiation events,
($\mu_{4}=\lambda_{4}=\xi_{4}=0$) to represent terminal accumulation of cells
in it. This allows us to quantify the number of cells that exit the system
formed by the first three compartments (as opposed to dying), which can be of
interest in processes such as thymic development Sawicka2014 . We select
$\mu_{i}=1.0$ for all $i\in\\{1,2,3\\}$, so that the unit of time for the
rates is the average lifetime of a cell in the system. We aim to analyse the
impact of asymmetric and symmetric division on the temporal dynamics, and
thus, set $\nu_{1}=\nu_{2}=\nu_{3}=0.5$, in the Irreversible model.

Cells can divide in each compartment at a rate $\omega$, and this division
represents self-renewal with probability $p_{SR}$, asymmetric division with
probability $p_{AD}$, and symmetric division with probability $p_{SD}$. This
is equivalent to setting, with the notation introduced in Section 2,
$\lambda_{i}=p_{SR}\cdot\omega$, $s_{i}=p_{SD}\cdot\omega$ and
$a_{i}=p_{AD}\cdot\omega$, for $i\in\\{1,2,3\\}$. We choose $\omega=0.9$, so
that the system has significant proliferative potential, and focus on the
following scenarios of interest:

  1. [leftmargin=3cm]

  2. Scenario 1.

$p_{SR}=1.0$, $p_{SD}=0$, $p_{AD}=0$.

  3. Scenario 2.

$p_{SR}=0.8$, $p_{SD}=0.1$, $p_{AD}=0.1$.

  4. Scenario 3.

$p_{SR}=0.1$, $p_{SD}=0.8$, $p_{AD}=0.1$.

  5. Scenario 4.

$p_{SR}=0.1$, $p_{SD}=0.1$, $p_{AD}=0.8$.

Our aim is to explore the impact that asymmetric or symmetric proliferation
has on the dynamics of the system (Scenarios 2, 3 and 4), compared to the
situation where only self-renewal takes place (Scenario 1). In Figure 5 we
plot the mean number of cells, $\mathbb{E}[{\bm{C}}_{i}(t)]$, in compartments
$i\in\\{1,2,3,4\\}$ in Scenarios 1-4. For compartments $C_{i}$ with
$i\in\\{1,2,3\\}$, and since $\Delta=\Delta_{1}=\Delta_{2}=\Delta_{3}$ and
$\Lambda=\Lambda_{1}=\Lambda_{2}=\Lambda_{3}$, one can directly use Eq. (8),
while for compartment $C_{N}$, $N=4$, one has

$\displaystyle\mathbb{E}[{\bm{C}}_{N}(t)]$ $\displaystyle=$
$\displaystyle\left(\frac{\Lambda}{\Delta}\right)^{N-1}-\sum_{j=0}^{N-1}\mathbb{E}({\bm{C}}_{j})\left(\frac{\Lambda}{\Delta}\right)^{N-j}.$

In Figure 5, we consider initial conditions ${\bm{C}}_{1}(0)=10^{2}$,
${\bm{C}}_{2}(0)={\bm{C}}_{3}(0)={\bm{C}}_{4}(0)=0$, representing $10^{2}$
initial cells in the first compartment and no cells in the other compartments.
We observe that an exponential decay in the number of cells in $C_{1}$ is
followed by sequential increases in the other two compartments until a steady-
state number of cells is achieved in the exiting compartment, $C_{4}$.
Interestingly, the dynamics of the system is faster once symmetric or
asymmetric division is considered (Scenarios 2-4 compared to 1), so that the
decay in compartment $C_{1}$ is quicker and the steady-state is reached
sooner.

The fastest dynamics is observed for Scenario 3, where symmetric division is
more likely, and the two daughters of a progenitor cell move directly to the
next compartment. It is also worth highlighting that symmetric or asymmetric
division does not only affect the speed of the dynamics but also the total
average number of cells exiting the system (i.e., reaching compartment
$C_{4}$), where $\lim_{t\rightarrow+\infty}\mathbb{E}[{\bm{C}}_{4}(t)]$ is
significantly larger when asymmetric and (especially) symmetric division can
occur. We note here that, importantly, the per cell division rate $\omega$ is
equal in all four scenarios. This suggests that, in this type of systems,
asymmetric or symmetric division (as opposed to self-renewal proliferation)
represents a mechanism which facilitates generating a larger population of
exiting cells with the same proliferative effort. That is, in Scenario 1 where
only self-renewal can occur, a larger number of divisions is required in each
compartment for enough cells to escape death and differentiate to the next
compartment, and to eventually reach $C_{4}$.

Figure 5: Mean number of cells, $\mathbb{E}[{\bm{C}}_{i}(t)]$, predicted over
time for compartments $i\in\\{1,2,3,4\\}$ in scenarios 1, 2, 3 and 4.

Our comments above are also consistent with the results of Figure 6, where we
plot the mean number of cells $m_{1}(j)$ in the genealogy of a single cell
starting in $C_{1}$, belonging to compartments $C_{j}$, $j\in\\{1,2,3,4\\}$,
for Scenarios 1-4. In Scenario 1, where only self-renewal can occur, the mean
number of cells in compartments $C_{1}$, $C_{2}$, $C_{3}$, and $C_{4}$ in the
genealogy of the progenitor cell decreases monotonically across the sequence
of compartments, $m_{1}(1)>m_{1}(2)>m_{1}(3)>m_{1}(4)$. We note that
$m_{1}(4)=0$ can be explained since it only accounts for progeny cells which
arrive into compartment $C_{4}$ as a direct result of cell proliferation, and
no symmetric or asymmetric division is considered in Scenario 1. We also
stress here that this monotonic decrease happens even though division and
differentiation rates are equal in all compartments $j\in\\{1,2,3\\}$. This is
related to the fact that some cells in the genealogy will die before reaching
compartments $C_{2}$ and/or $C_{3}$. On the other hand, and related to our
comments corresponding to Figure 5, Scenario 1 leads to the largest progeny,
suggesting that it is an inefficient system for reaching a desired population
of exiting cells in $C_{4}$. Asymmetric and (especially) symmetric division
events significantly reduce the number of descendants from the single
progenitor cell in all compartments, while maximising the number of exiting
cells as described by Figure 5. In particular, Scenario 3 is characterised by
the highest total average number of cells exiting the system,
$\mathbb{E}[{\bm{C}}_{4}(+\infty)]$, as well as the smallest progeny size,
$m_{1}$, while leading to the largest progeny, $m_{1}(4)$, in the final
compartment.

Figure 6: Mean number of cells $m_{1}(j)$ in the genealogy of a single cell
starting in $C_{1}$, belonging to compartments $C_{j}$, $j\in\\{1,2,3\\}$, for
Scenarios 1-4.

###  3.2 Tracking a thymocyte during thymic development

We consider now the model proposed in Ref. (Sawicka2014, , Model 2), and shown
in Figure 7. Double negative (DN) thymocytes differentiate to become pre-
selection DP thymocytes (pre-DP). In this model, pre-DP is the first
compartment, which will contain an initial number of cells (initial condition,
$C_{1}(0)$). Pre-DPs will undergo maturation in the thymus. These cells can
progress to the double positive stage (post-DP), where thymocytes express both
CD4 and CD8 co-receptors. Post-DP cells that are positively selected
transition to the single positive (SP) stage, where they can express either
the CD4 or CD8 co-receptor. Some of these cells will then reach the periphery
as CD4 or CD8 SP cells.

Figure 7: Schematic representation of the thymic development model proposed in
Ref. Sawicka2014 . Grey cells represent cellular death.

We exploit this particular model to illustrate the applicability of our
single-cell analysis proposed in Section 2.3. This also allows us to show how
our methods can be easily adapted even under slightly modifications to the
network topology of compartments (here, a bifurcation in the last
compartmental stage is observed, rather than a linear sequence of
compartments).

First, it is of interest to quantify cellular fate in this system: what is the
percentage of pre-DP thymocytes that are predicted to die in each of the
compartments during maturation, and what percentage will successfully reach
the periphery instead (either as a CD4 or CD8 SP cell). It is clear that our
arguments in Section 2.3 can be easily adapted to quantify this. In
particular, one can slightly redefine the probabilities $\beta_{i}(j)$ in
Section 2.3, focusing on $i=1$ (i.e., a single pre-DP thymocyte being
tracked), as

$\displaystyle\beta_{1}(1)$ $\displaystyle=$ $\displaystyle\hbox{\small
probability that the pre-DP thymocyte dies in the pre-DP compartment},$
$\displaystyle\beta_{1}(2)$ $\displaystyle=$ $\displaystyle\hbox{\small
probability that the pre-DP thymocyte dies in the post-DP compartment},$
$\displaystyle\beta_{1}(4)$ $\displaystyle=$ $\displaystyle\hbox{\small
probability that the pre-DP thymocyte dies in the CD4 SP compartment},$
$\displaystyle\beta_{1}(8)$ $\displaystyle=$ $\displaystyle\hbox{\small
probability that the pre-DP thymocyte dies in the CD8 SP compartment},$
$\displaystyle\beta_{1}(4P)$ $\displaystyle=$ $\displaystyle\hbox{\small
probability that the pre-DP thymocyte reaches the periphery as a CD4 SP
cell},$ $\displaystyle\beta_{1}(8P)$ $\displaystyle=$
$\displaystyle\hbox{\small probability that the pre-DP thymocyte reaches the
periphery as a CD8 SP cell}.$

Similar solutions to those in Section 2.3 can be obtained by incorporating the
compartmental bifurcation in the first-step analysis, leading to

$\displaystyle\beta_{1}(1)$ $\displaystyle=$
$\displaystyle\frac{\mu_{1}}{\mu_{1}+\nu_{1}},$ $\displaystyle\beta_{1}(2)$
$\displaystyle=$
$\displaystyle\frac{\nu_{1}}{\mu_{1}+\nu_{1}}\frac{\mu_{2}}{\mu_{2}+\nu_{24}+\nu_{28}},$
$\displaystyle\beta_{1}(4)$ $\displaystyle=$
$\displaystyle\frac{\nu_{1}}{\mu_{1}+\nu_{1}}\frac{\nu_{24}}{\mu_{2}+\nu_{24}+\nu_{28}}\frac{\mu_{4}}{\mu_{4}+\nu_{4}},$
$\displaystyle\beta_{1}(8)$ $\displaystyle=$
$\displaystyle\frac{\nu_{1}}{\mu_{1}+\nu_{1}}\frac{\nu_{28}}{\mu_{2}+\nu_{24}+\nu_{28}}\frac{\mu_{8}}{\mu_{8}+\nu_{8}},$
$\displaystyle\beta_{1}(4P)$ $\displaystyle=$
$\displaystyle\frac{\nu_{1}}{\mu_{1}+\nu_{1}}\frac{\nu_{24}}{\mu_{2}+\nu_{24}+\nu_{28}}\frac{\nu_{4}}{\mu_{4}+\nu_{4}},$
$\displaystyle\beta_{1}(8P)$ $\displaystyle=$
$\displaystyle\frac{\nu_{1}}{\mu_{1}+\nu_{1}}\frac{\nu_{28}}{\mu_{2}+\nu_{24}+\nu_{28}}\frac{\nu_{8}}{\mu_{8}+\nu_{8}}.$

The proliferation potential of thymocytes during thymic development directly
depends on them reaching the CD4 SP and CD8 SP compartments, where they are
able to divide, before they reach the periphery. Thus, the average number of
divisions initiated by a single pre-DP thymocyte during its thymic development
journey, $\eta_{1}=\eta_{1}(4)+\eta_{1}(8)$, can be obtained as

$\displaystyle\eta_{1}(4)$ $\displaystyle=$
$\displaystyle\frac{\nu_{1}\nu_{24}}{(\mu_{1}+\nu_{1})(\mu_{2}+\nu_{24}+\nu_{28})}\frac{\lambda_{4}}{\mu_{4}+\nu_{4}},$
$\displaystyle\eta_{1}(8)$ $\displaystyle=$
$\displaystyle\frac{\nu_{1}\nu_{28}}{(\mu_{1}+\nu_{1})(\mu_{2}+\nu_{24}+\nu_{28})}\frac{\lambda_{8}}{\mu_{8}+\nu_{8}}.$

Finally, the average lifespan of a pre-DP cell during thymic development
(i.e., the mean time until it dies or it reaches the periphery) is given by

$\displaystyle\tau_{1}=\frac{1}{\mu_{1}+\nu_{1}}\left[\frac{\nu_{1}}{\mu_{2}+\nu_{24}+\nu_{28}}\left(\frac{\nu_{24}}{\mu_{4}+\nu_{4}}+\frac{\nu_{28}}{\mu_{8}+\nu_{8}}+1\right)+1\right].$

We select parameter values from those reported in Ref. (Sawicka2014, , Section
3.2) (see Table 1). Our methodology allows one to predict that the average
lifespan of a pre-DP cell during the maturation process (until it dies or
reaches the periphery) is $2.84$ days. During this lifetime, the cell may
undergo differentiation and proliferation, before dying in one of the
compartments without ever reaching the periphery, or reaching the periphery
either as a CD4 or CD8 cell. In Figure 8, we plot the predicted cellular fate
probabilities. We observe that the most likely outcomes correspond to cell
death, especially during the early stages of maturation (pre-DP and post-DP
compartments). Once the cell reaches the CD4 SP compartment, it is more likely
for this cell to reach the periphery than dying in that compartment, while
these probabilities are comparable in the CD8 case. The average number of
proliferation events triggered by this single pre-DP cell is
$\eta_{1}=\eta_{1}(4)+\eta_{1}(8)=0.0139+0.0046=0.0185$, meaning that out of
$100$ pre-DP cells entering the thymic development process only around $18$
cells are expected to be produced via direct proliferation from the original
cells when visiting the CD4 SP or CD8 SP compartments. These small quantities
are directly related to the small probabilities of reaching these compartments
at all, so that the cell can actually proliferate. Our results here are in
agreement with conclusions and insights from Sawicka2014 , where the authors
state that thymic development is an stringent process characterized by an
extremely low success rate.

Rate | $\mu_{1}$ | $\nu_{1}$ | $\mu_{2}$ | $\nu_{24}$ | $\nu_{28}$ | $\lambda_{4}$ | $\lambda_{8}$ | $\mu_{4}$ | $\mu_{8}$ | $\nu_{4}$ | $\nu_{8}$
---|---|---|---|---|---|---|---|---|---|---|---
Value | $0.263$ | $0.137$ | $1.369$ | $0.07$ | $0.054$ | $0.216$ | $0.093$ | $0.04$ | $0.11$ | $0.21$ | $0.14$

Table 1: Parameter values as reported in Ref. (Sawicka2014, , Section 3.2), in
units $days^{-1}$. Figure 8: Probabilities of a single pre-DP cell to die in
each of the compartments (pre-DP, post-DP, CD4 SP, or CD8 SP) before reaching
the periphery, or reaching the periphery as a CD4 or CD8 SP cell.

##  4 Discussion

We have presented a general model for cell proliferation and differentiation
in terms of a stochastic process across a sequence of compartments. Cells can
divide, die or move to adjacent compartments over time. We have explicitly
quantified the average number of cells in each compartment over time under
different scenarios of interest (e.g., the irreversible model, where
differentiation cannot be reversed), and analysed the genealogy of a single
progenitor cell in terms of summary statistics appropriately defined. Single-
cell analysis allows one to track the journey of one cell during its lifetime
across the system of compartments, where we are able to quantify the cell
lifetime, its proliferation potential and the likelihood of the different
potential cellular fates. We have presented a number of numerical results to
illustrate the applicability of our techniques, and the impact that model
parameters can have on the corresponding predictions.

A particular limitation of our approach is that it relies on cells behaving
independently from each other, which is directly related to the fact that the
corresponding system of ODEs for the average behaviour of the process is
linear. This allows one to implement techniques from the theory of branching
processes, and makes the single-cell analysis implementation feasible (where
one can follow a single tracked cell while neglecting the dynamics of its
neighbours). On the other hand, cellular interactions might not be negligible
in some situations (e.g., under competition for resources, where logistic
growth-type models might be more appropriate). Relaxing this particular
assumption is thus, the aim of future work.

##  5 Acknowledgments

H. Dreiwi was a Daphne Jackson Fellow during the period 2019-2021 and would
like to thank the Daphne Jackson Trust and the University of Leeds, School of
Mathematics, for their financial and academic support. This work has been
supported by the European Commission through the Marie Sklodowska-Curie Action
(H2020-MSCA-ITN-2017) Innovative Training Network Quantitative T cell
Immunology and Immunotherapy (QuanTII), project number 764698 (FF, GL, and
CMP).

## References

  *   * (1) Martin, G.: Isolation of a pluripotent cell line from early mouse embryos cultured in medium conditioned by teratocarcinoma stem cells. Proc. Natl. Acad. Sci. (1981)
  * (2) Evans, M., Kaufman, M.: Establishment in culture of pluripotent cells from mouse embryos. Nature (1981)
  * (3) Zhang, H., Hou, W., Henrot, L., Schnebert, S., Dumas, M., Heusèle, C., Yang, J.: Modelling epidermis homoeostasis and psoriasis pathogenesis. Journal of The Royal Society Interface 12(103), 20141071 (2015)
  * (4) Sawicka, M., Stritesky, G., Reynolds, J., Abourashchi, N., Lythe, G., Molina-Paris, C., Hogquist, K.: From pre-DP, post-DP, SP4, and SP8 thymocyte cell counts to a dynamical model of cortical and medullary selection. Front Immunol. (2014). https://doi.org/%****␣sn-article.bbl␣Line␣100␣****10.3389/fimmu.2014.00019
  * (5) Wherry, E.J., M, M.K.: Molecular and cellular insights into T cell exhaustion. Nat Rev Immunol. (2015). https://doi.org/10.1038/nri3862
  * (6) Chen, Z., Ji, Z., Ngiow, S.F., Manne, S., Cai, Z., Huang, A.C., Johnson, J., Staupe, R.P., Bengsch, B., Xu, C., Yu, S., Kurachi, M., Herati, R.S., Vella, L.A., Baxter, A.E., Wu, J.E., Khan, O., Beltra, J.-C., Giles, J.R., Stelekati, E., McLane, L.M., Lau, C.W., Yang, X., Berger, S.L., Vahedi, G., Ji, H., , Wherry, E.J.: TCF-1-centered transcriptional network drives an effector versus exhausted CD8 T cell-fate decision. Immunity (2019). https://doi.org/10.1016/j.immuni.2019.09.013
  * (7) Hodzic, E.: Single-cell analysis: Advances and future perspectives. Bosn J Basic Med Sci. (2016)
  * (8) Perié, L., Hodgkin, H.P.D., Naik, S.H., Schumacher, T.N., de Boer, R.J., Duffy, K.R.: Determining lineage pathways from cellular barcoding experiments. Cell reports (2014)
  * (9) Kimmel, M., Axelrod, D.E.: Branching Processes in Biology, 9th edn. Springer, USA (2002)
  * (10) Watson, H.W., Galton, F.: On the probability of the extinction of families. The Journal of the Anthropological Institute of Great Britain and Ireland 4, 138–144 (1875)
  * (11) de la Higuera, L., López-García, M., Castro, M., Abourashchi, N., Lythe, G., Molina-París, C.: Fate of a naive t cell: a stochastic journey. Frontiers in Immunology 10, 194 (2019)
  * (12) Nordon, R.E., Ko, K.-H., Odell, R., Schroeder, T.: Multi-type branching models to describe cell differentiation programs. Journal of Theoretical Biology (2011). https://doi.org/10.1016/j.jtbi.2011.02.006
  * (13) Matis, J.H., Hartley, H.O.: Stochastic compartmental analysis: Model and least squares estimation from time series data. International Biometric Society (1971)
  * (14) Matis, J.H., Wehrly, T.E.: Compartmental models of ecological and environmental systems. Handbook of Statistics, vol. 12, pp. 583–613. Elsevier (1994). https://doi.org/10.1016/S0169-7161(05)80019-7
  * (15) Eftimie, R., Gillard, J.J., Cantrell, D.A.: Mathematical models for immunology: Current state of the art and future research directions. Bull Math Biol (2016). https://doi.org/10.1007/s11538-016-0214-9
  * (16) Capasso, V.: Mathematical structures of epidemic systems (2020)
  * (17) Klein, A.M., Simons, B.D.: Universal patterns of stem cell fate in cycling adult tissues. Development (2011). https://doi.org/10.1242/dev.060103
  * (18) Barile, M., Busch, K., Fanti, A.-K., Greco, A., Wang, X., Oguro, H., Zhang, Q., Morrison, S.J., Rodewald, H.-R., Höfer, T.: Hematopoietic stem cells self-renew symmetrically or gradually proceed to differentiation. bioRxiv (2020). https://doi.org/10.1101/2020.08.06.239186
  * (19) Allen, L.J.: An Introduction to Mathematical Biolog. Pearson/Prentice Hall, USA (2007)
  * (20) Callier, F.M., Desoer, C.A.: Linear System Theory. Springer, USA (1991)
  * (21) Conte, S., De Boor, C.: Elementary numerical analysis: An algorithmic appr mcgraw-hill. New York (1972)
  * (22) Thomas, L.H.: Elliptic problems in linear difference equations over a network. Watson Sci. Comput. Lab. Rept., Columbia University, New York 1, 71 (1949)