# On a generalization of Jacobi’s elegantissima

Luc Haine  Luc Haine, Institut de Recherche en Mathématique et Physique,
Université catholique de Louvain, Chemin du Cyclotron 2, 1348 Louvain-la-
Neuve, Belgium<EMAIL_ADDRESS>To the memory of Hermann Flaschka

(Date: May 11, 2023; revised August 15, 2023)

###### Abstract.

We establish a generalization of Jacobi’s elegantissima, which solves the
pendulum equation. This amazing formula appears in lectures by the famous
cosmologist Georges Lemaître, during the academic years 1955-1956 and
1956-1957. Our approach uses the full power of Jacobi’s elliptic functions, in
particular imagi-nary time is crucial for obtaining the result.

###### Key words and phrases:

Elliptic functions, Landen’s transformation

###### 2020 Mathematics Subject Classification:

Primary: 33E05; Secondary: 97I80

###### Contents

  1. 1 Introduction
  2. 2 Jacobi’s elliptic functions and the pendulum
  3. 3 Lemaître’s generalization of Jacobi’s elegantissima
  4. 4 The circulatory motion
  5. 5 A sketch of Lemaître’s proof

##  1\. Introduction

In _”Fundamenta Nova Theoriae Functionum Ellipticarum”_ , Jacobi obtained the
following formula for the ”modular angle”

$\frac{\arcsin
k}{4}=\arctan\sqrt{q}-\arctan\sqrt{q^{3}}+\arctan\sqrt{q^{5}}-\ldots,$ (1.1)

which is an inversion formula for

$\displaystyle K$
$\displaystyle=\int_{0}^{1}\frac{dx}{\sqrt{(1-x^{2})(1-k^{2}x^{2})}},\;0<k<1,$
(1.2) $\displaystyle K^{\prime}$
$\displaystyle=\int_{0}^{1}\frac{dx}{\sqrt{(1-x^{2})(1-k^{\prime
2}x^{2})}},\;k^{\prime}=\sqrt{1-k^{2}},$ (1.3)

in terms of the so-called _nome_

$q=e^{-\pi\frac{K^{\prime}}{K}}.$

Immediately following his formula (see [6] 40, formula (47)), Jacobi wrote
_”quae inter formulas elegantissima censeri debet”_. The only modern reference
mentioning (1.1) seems to be [13], Example 22.5.2.

During the two academic years 1955-1956 and 1956-1957, Georges Lemaître, the
famous Belgian cosmologist, professor at the University of Louvain and father
of the big bang theory, taught his _”Leçons de mécanique. Le pendule”_ [9],
[10]. For a beautiful account of Georges Lemaître’s life and work, we refer
the reader to [7].

In the context of the oscillatory motion of a pendulum of length $l$ under the
influence of the force of gravity, as explained in Section 2, Jacobi’s formula
can be seen as a formula giving the maximum deviation angle in terms of the
period $4K\sqrt{l/g}$ of the oscillation and the _imaginary period_
$4K^{\prime}\sqrt{l/g}$ of the _complementary motion_ obtained by inverting
the direction of gravity, when the pendulum reaches the maximum angle. In his
lectures [9], [10], Georges Lemaître derived a formula of a similar type,
giving the angle of the oscillation at _any_ time, of which Jacobi’s formula
(1.1) becomes a particular case. According to his own words _”La théorie que
nous avons exposée peut être considérée comme une généralisation de
l’élégantissime de Jacobi”_.

Lemaître’s derivation of his generalization of Jacobi’s elegantissima is based
on a very ingenious infinite iteration of Landen’s transformation [8], a well
known transformation of elliptic functions, see [11], [12] for background. He
started his lectures with elementary mechanical and geometric proofs of
Jacobi’s interpretation [5] of Poncelet theorem and Landen’s transformation.
Though not mentioned in his lectures, some of the material can also be found
in Greenhill’s book ”The Applications of Elliptic Functions”; the French
translation [4] we refer to, with a preface by Paul Appell, is a revised
version of the original work.

To the best of our knowledge, the infinite iteration of Landen’s
transformation to obtain the solution of the pendulum, does not seem to appear
anywhere else in the literature. Though very intuitive, the passage to the
limit is not carefully justified in [9], [10]. The purpose of this note is to
give a direct proof of Lemaître’s formula, based on Jacobi’s theory of
elliptic functions, emphasizing _the role played by ima-ginary time_. This is
the content of Sections 3 and 4. In Section 5, we provide a sketch of
Lemaître’s proof. We don’t follow his geometric arguments, using instead a
standard version of Landen’s transformation, to be found in [13]. Our hope is
to attract the attention on Lemaître’s elementary and geometric approach of
elliptic functions.

##  2\. Jacobi’s elliptic functions and the pendulum

The motion of a pendulum of mass $m$ and length $l$ under the influence of the
force of gravity is given by

$ml\ddot{\varphi}(t)=-mg\sin\varphi(t)\Leftrightarrow\ddot{\varphi}(t)=-\frac{\sin\varphi(t)}{R},\;R=\frac{l}{g},$
(2.1)

with $g$ the acceleration of the gravity, and $\varphi(t)$ the angle made with
the descending vertical. Equation (2.1) has the following elementary
solutions, the two equilibria given by

$\varphi(t)=0,\;\forall\;t\in\mathbb{R},\quad\varphi(t)=\pi,\;\forall\;t\in\mathbb{R},$

and, assuming $\varphi(0)=0$, the two doubly asymptotic solutions

$\varphi(t)=4\arctan e^{\pm\frac{t}{\sqrt{R}}}-\pi.$ (2.2)

All other solutions can be expressed in terms of Jacobi elliptic functions.
Assuming $\varphi(0)=0$ and $\dot{\varphi}(0)>0$, as long as
$\dot{\varphi}(t)>0$, one has

$\sqrt{\frac{m}{2}}l\int_{0}^{\varphi(t)}\frac{d\varphi}{\sqrt{E-2mgl\sin^{2}\frac{\varphi}{2}}}=t,$
(2.3)

which follows from the conservation of energy

$\frac{ml^{2}\dot{\varphi}^{2}}{2}+2mgl\sin^{2}\frac{\varphi}{2}=E.$

_Oscillatory motions_ occur when

$0<E<2mgl.$

By the change of variable

$\sin\frac{\varphi}{2}=kx,$

with

$0<k=\sin\frac{\alpha}{2}<1,$ (2.4)

where $0<\alpha<\pi$, gives the maximum angle of oscillation, i.e.
$\dot{\varphi}(\alpha)=0$, equation (2.3) becomes

$\int_{0}^{x(t)}\frac{dx}{\sqrt{(1-x^{2})(1-k^{2}x^{2})}}=\frac{t}{\sqrt{R}}.$

Jacobi’s elliptic sine function ”sn” solves this inversion problem. Thus

$\sin\frac{\varphi(t)}{2}=k\;\mbox{sn}\Big{(}\frac{t}{\sqrt{R}},k\Big{)},\;\cos\frac{\varphi(t)}{2}=\mbox{dn}\Big{(}\frac{t}{\sqrt{R}},k\Big{)},\;\forall\;t\in\mathbb{R},$
(2.5)

where the second formula follows from the identity
$\mbox{dn}^{2}(z,k)=1-k^{2}\;\mbox{sn}^{2}(z,k)$, satisfied by Jacobi’s
elliptic sine ”sn” and delta amplitude ”dn” functions. Denoting by $T$ the
period of the motion,

$\frac{T}{4\sqrt{R}}=\int_{0}^{1}\frac{dx}{\sqrt{(1-x^{2})(1-k^{2}x^{2})}}\Leftrightarrow
T=4K\sqrt{R},$ (2.6)

with $K$ defined as in (1.2).

For an _oscillatory motion_ , the _complementary motion_ is obtained by
changing the sign of gravity, when the angle of oscillation reaches its
maximum value $\alpha$. This motion follows the arc of circle complementary to
the arc followed during the oscillation. It is characterized by

$\displaystyle\ddot{\phi}(t)=-\frac{\sin\phi(t)}{R},\;\phi(0)=\alpha-\pi,\;\dot{\phi}(0)=0.$

One immediately checks that the solution is given by

$\phi(t)=\varphi(it+K\sqrt{R})-\pi.$

Since

$\sin^{2}\frac{\phi(0)}{2}=\cos^{2}\frac{\alpha}{2}=1-k^{2},$

from (2.6), it follows that the period of the complementary motion is

$T^{\prime}=4K^{\prime}\sqrt{R},$

with $K^{\prime}$ as in (1.3). The mechanical interpretation of $T^{\prime}$,
using _imaginary time_ , is due to Paul Appell, see [1].

_Circulatory motions_ occur when

$E>2mgl\Leftrightarrow\frac{\dot{\varphi}^{2}(0)}{4}>\frac{1}{R}.$

By the change of variable $x=\sin\frac{\varphi}{2}$, equation (2.3) becomes

$\int_{0}^{x(t)}\frac{dx}{\sqrt{(1-x^{2})(1-k^{2}x^{2})}}=\frac{t}{k\sqrt{R}},\quad
0<k=\frac{2}{\sqrt{R}\;\dot{\varphi}(0)}<1,$

leading to

$\sin\frac{\varphi(t)}{2}=\mbox{sn}\Big{(}\frac{t}{k\sqrt{R}},k\Big{)},\;\cos\frac{\varphi(t)}{2}=\mbox{cn}\Big{(}\frac{t}{k\sqrt{R}},k\Big{)},\;\forall\;t\in\mathbb{R},$
(2.7)

where the second formula follows from the identity
$\mbox{cn}^{2}(z,k)+\mbox{sn}^{2}(z,k)=1$, satisfied by Jacobi’s elliptic sine
and cosine functions denoted ”sn” and ”cn”. Using that
$\mbox{sn}(z+2K,k)=-\mbox{sn}(z,k)$, one finds the period $T$ of the
circulatory motion

$T=2kK\sqrt{R}.$ (2.8)

##  3\. Lemaître’s generalization of Jacobi’s elegantissima

For the oscillatory motion, defining

$\theta(t)=\pi+\varphi(t),$

i.e. measuring the angle of oscillation from the ascending vertical, from
(2.5) one gets

$\displaystyle\cos\frac{\theta(t)}{2}$
$\displaystyle=-\sin\frac{\varphi(t)}{2}=-k\;\mbox{sn}\Big{(}\frac{t}{\sqrt{R}},k\Big{)},$
(3.1) $\displaystyle\sin\frac{\theta(t)}{2}$
$\displaystyle=\cos\frac{\varphi(t)}{2}=\mbox{dn}\Big{(}\frac{t}{\sqrt{R}},k\Big{)}.$
(3.2)

To derive from (3.1), (3.2) a nice formula for $e^{i\frac{\theta(t)}{2}}$, one
is tempted to use the following infinite product formula (see [13] 22.5, for a
proof)

$\mbox{dn}(z,k)+ik\;\mbox{sn}(z,k)=\\\
\prod_{n=1}^{\infty}\frac{\big{(}1+(-1)^{n}q^{n-\frac{1}{2}}e^{-\frac{i\pi
z}{2K}}\big{)}\big{(}1-(-1)^{n}q^{n-\frac{1}{2}}e^{\frac{i\pi
z}{2K}}\big{)}}{\big{(}1-(-1)^{n}q^{n-\frac{1}{2}}e^{-\frac{i\pi
z}{2K}}\big{)}\big{(}1+(-1)^{n}q^{n-\frac{1}{2}}e^{\frac{i\pi
z}{2K}}\big{)}},$ (3.3)

with

$q=e^{-\pi\frac{K^{\prime}}{K}},$ (3.4)

then to take the logarithm to get $\theta(t)$. However, it is not immediately
appropriate because of the imaginary exponentials which would appear in the
infinite product. Thus we first rewrite the solution using Jacobi’s imaginary
transformation (see [13] 22.41) combined with standard formulas, which can all
be obtained from the addition theorems for ”sn”, ”cn” and ”dn”. Since

$\mbox{sn}(z-K,k)=-\frac{\mbox{cn}(z,k)}{\mbox{dn}(z,k)},\;\mbox{dn}(iz+K^{\prime},k^{\prime})=\frac{k}{\mbox{dn}(iz,k^{\prime})}=k\;\frac{\mbox{cn}(z,k)}{\mbox{dn}(z,k)},$

we deduce that

$\mbox{sn}(z-K,k)=-\frac{1}{k}\;\mbox{dn}(iz+K^{\prime},k^{\prime}).$ (3.5)

Similarly, since

$\mbox{dn}(z-K,k)=\frac{k^{\prime}}{\mbox{dn}(z,k)},\;\mbox{sn}(iz+K^{\prime},k^{\prime})=\frac{\mbox{cn}(iz,k^{\prime})}{\mbox{dn}(iz,k^{\prime})}=\frac{1}{\mbox{dn}(z,k)},$

one has

$\mbox{dn}(z-K,k)=k^{\prime}\;\mbox{sn}(iz+K^{\prime},k^{\prime}).$ (3.6)

Combining (3.1), (3.2), (3.5) and (3.6), from (3.3), with $k^{\prime}$ instead
of $k$, we obtain

$\displaystyle e^{i\frac{\theta(t-K\sqrt{R})}{2}}$
$\displaystyle=\mbox{dn}\Big{(}i\;\frac{t}{\sqrt{R}}+K^{\prime},k^{\prime}\Big{)}+ik^{\prime}\;\mbox{sn}\Big{(}i\;\frac{t}{\sqrt{R}}+K^{\prime},k^{\prime}\Big{)},$
$\displaystyle=\prod_{n=1}^{\infty}\frac{\Big{(}1-i(-1)^{n}q^{\prime
n-\frac{1}{2}}e^{\frac{\pi
t}{2K^{\prime}\sqrt{R}}}\Big{)}\Big{(}1-i(-1)^{n}q^{\prime
n-\frac{1}{2}}e^{-\frac{\pi
t}{2K^{\prime}\sqrt{R}}}\Big{)}}{\Big{(}1+i(-1)^{n}q^{\prime
n-\frac{1}{2}}e^{\frac{\pi
t}{2K^{\prime}\sqrt{R}}}\Big{)}\Big{(}1+i(-1)^{n}q^{\prime
n-\frac{1}{2}}e^{-\frac{\pi t}{2K^{\prime}\sqrt{R}}}\Big{)}},$

with

$q^{\prime}=e^{-\pi\frac{K}{K^{\prime}}}.$ (3.7)

Since for $a\in\mathbb{R}$,

$e^{(-1)^{n}2i\;arctan\;a}=\frac{1+(-1)^{n}ia}{1-(-1)^{n}ia},$

taking the logarithm of the last formula, gives

$\theta\Big{(}t-\frac{T}{4}\Big{)}=\\\
\quad\sum_{n=1}^{\infty}(-1)^{n-1}4\Big{\\{}\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{\frac{\pi
t}{2K^{\prime}\sqrt{R}}}\Big{)}+\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{-\frac{\pi t}{2K^{\prime}\sqrt{R}}}\Big{)}\Big{\\}},$ (3.8)

with $T=4K\sqrt{R}$ as in (2.6), the period of the oscillatory motion.
Evaluated at $t=0$, remembering that $\varphi(0)=0$,

$\theta\Big{(}-\frac{T}{4}\Big{)}=\pi+\varphi\Big{(}-\frac{T}{4}\Big{)}=\pi-\alpha,$

with $\alpha$ the maximum angle of the oscillation, (3.8) reduces to

$\pi-\alpha=8\;\sum_{n=1}^{\infty}(-1)^{n-1}\arctan\sqrt{q^{\prime 2n-1}}.$

Since the maximum angle of oscillation of the complementary motion is
$\pi-\alpha$, permuting $K$ et $K^{\prime}$ and remembering (2.4) and (3.7),
we obtain

$2\arcsin
k=\alpha=\pi-(\pi-\alpha)=8\sum_{n=1}^{\infty}(-1)^{n-1}\arctan\sqrt{q^{2n-1}},$

with $q$ as in (3.4). This is exactly _Jacobi’s elegantissima_ (1.1). Thus,
formula (3.8), to be found in Lemaître [9], [10], deserves to be called _the
generalized elegantissima_.

##  4\. The circulatory motion

Putting again

$\theta(t)=\pi+\varphi(t),$

from (2.7), we now obtain

$\displaystyle\cos\frac{\theta(t)}{2}$
$\displaystyle=-\sin\frac{\varphi(t)}{2}=-\mbox{sn}\Big{(}\frac{t}{k\sqrt{R}},k\Big{)},$
(4.1) $\displaystyle\sin\frac{\theta(t)}{2}$
$\displaystyle=\cos\frac{\varphi(t)}{2}=\mbox{cn}\Big{(}\frac{t}{k\sqrt{R}},k\Big{)}.$
(4.2)

Once again, we need to pass to imaginary time. Now we have

$\displaystyle\mbox{cn}(z-K,k)$
$\displaystyle=k^{\prime}\;\frac{\mbox{sn}(z,k)}{\mbox{dn}(z,k)},$
$\displaystyle\mbox{cn}(iz+K^{\prime},k^{\prime})$
$\displaystyle=-k\;\frac{\mbox{sn}(iz,k^{\prime})}{\mbox{dn}(iz,k^{\prime})}=-ik\;\frac{\mbox{sn}(z,k)}{\mbox{dn}(z,k)},$

and thus

$\mbox{cn}(z-K,k)=i\frac{k^{\prime}}{k}\;\mbox{cn}(iz+K^{\prime},k^{\prime}).$
(4.3)

Combining (3.5), (4.1), (4.2) and (4.3), we obtain

$e^{i\frac{\theta(t-kK\sqrt{R})}{2}}=\frac{\mbox{dn}\Big{(}i\frac{t}{k\sqrt{R}}+K^{\prime},k^{\prime}\Big{)}-k^{\prime}\;\mbox{cn}\Big{(}i\frac{t}{k\sqrt{R}}+K^{\prime},k^{\prime}\Big{)}}{k}.$

Therefore, we are led to use the following infinite product formula (see [4]
266, formula (39), but it can also be proven along the same lines as (3.3)
following [13] 22.5, using an appropriate duplication formula)

$\frac{\mbox{dn}(z,k^{\prime})-k^{\prime}\;\mbox{cn}(z,k^{\prime})}{k}=\prod_{n=1}^{\infty}\frac{\big{(}1-q^{\prime
n-\frac{1}{2}}e^{-\frac{i\pi z}{2K^{\prime}}}\big{)}\big{(}1-q^{\prime
n-\frac{1}{2}}e^{\frac{i\pi z}{2K^{\prime}}}\big{)}}{\big{(}1+q^{\prime
n-\frac{1}{2}}e^{-\frac{i\pi z}{2K^{\prime}}}\big{)}\big{(}1+q^{\prime
n-\frac{1}{2}}e^{\frac{i\pi z}{2K^{\prime}}}\big{)}},$

with $q^{\prime}=e^{-\pi\frac{K}{K^{\prime}}}$. It gives

$e^{i\frac{\theta(t-kK\sqrt{R})}{2}}=\prod_{n=1}^{\infty}\frac{\Big{(}1+iq^{\prime
n-\frac{1}{2}}e^{\frac{\pi t}{2kK^{\prime}\sqrt{R}}}\Big{)}\Big{(}1-iq^{\prime
n-\frac{1}{2}}e^{-\frac{\pi
t}{2kK^{\prime}\sqrt{R}}}\Big{)}}{\Big{(}1-iq^{\prime
n-\frac{1}{2}}e^{\frac{\pi t}{2kK^{\prime}\sqrt{R}}}\Big{)}\Big{(}1+iq^{\prime
n-\frac{1}{2}}e^{-\frac{\pi t}{2kK^{\prime}\sqrt{R}}}\Big{)}}.$

Taking the logarithm, we now obtain the following formula

$\theta\Big{(}t-\frac{T}{2}\Big{)}=\\\
\sum_{n=1}^{\infty}4\Big{\\{}\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{\frac{\pi
t}{2kK^{\prime}\sqrt{R}}}\Big{)}-\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{-\frac{\pi t}{2kK^{\prime}\sqrt{R}}}\Big{)}\Big{\\}},$
(4.4)

with $T=2kK\sqrt{R}$ as in (2.8), the period of the circulatory motion.

##  5\. A sketch of Lemaître’s proof

In this Section, we don’t follow Lemaître’s arguments in detail, for which we
refer the reader to [10], the second version of his lectures during the
academic year 1956-1957. Our aim is to sketch the spirit of his approach, but
we use freely a standard form of Landen’s transformation to be found in [13].
We present his arguments in four steps. At each step, we refer to the sections
in the table of contents of [10], where elementary proofs using plane geometry
are given.

Step 1. The basic idea (see [10] 10, 11) is to use Landen’s transformation to
express the circulatory motion as a sum of two _new_ circulatory motions
_identical up to a shift by half a period_. Precisely, if $\theta(t)$ is a
circulatory motion with period $T=2kK\sqrt{R}$ as in (2.8), there exists a
circulatory motion $\theta_{1}(t)$ with period $T_{1}=2k_{1}K_{1}\sqrt{R_{1}}$
such that

$\theta\Big{(}\frac{R}{R_{1}}t\Big{)}=\theta_{1}(t)+\theta_{1}\Big{(}t-\frac{T_{1}}{2}\Big{)},$
(5.1)

and the period $T_{1}$ of $\theta_{1}(t)$ is related to the period $T$ of
$\theta(t)$ as follows

$T_{1}=2\frac{R_{1}}{R}T.$ (5.2)

To establish this result, using (4.1), (4.2) and the addition theorem for ”sn”
and ”cn”, we rewrite (5.1) as follows

$\displaystyle\cos\frac{\theta\Big{(}\frac{R}{R_{1}}t\Big{)}}{2}$
$\displaystyle=\cos\frac{\theta_{1}(t)}{2}\cos\frac{\theta_{1}(t-k_{1}K_{1}\sqrt{R_{1}})}{2}$
$\displaystyle\quad-\sin\frac{\theta_{1}(t)}{2}\sin\frac{\theta_{1}(t-k_{1}K_{1}\sqrt{R_{1}})}{2},$
$\displaystyle=\mbox{sn}\Big{(}\frac{t}{k_{1}\sqrt{R_{1}}},k_{1}\Big{)}\mbox{sn}\Big{(}\frac{t}{k_{1}\sqrt{R_{1}}}-K_{1},k_{1}\Big{)}$
$\displaystyle\quad-\mbox{cn}\Big{(}\frac{t}{k_{1}\sqrt{R_{1}}},k_{1}\Big{)}\mbox{cn}\Big{(}\frac{t}{k_{1}\sqrt{R_{1}}}-K_{1},k_{1}\Big{)},$
$\displaystyle=-(1+k^{\prime}_{1})\frac{\mbox{sn}\Big{(}\frac{t}{k_{1}\sqrt{R_{1}}},k_{1}\Big{)}\mbox{cn}\Big{(}\frac{t}{k_{1}\sqrt{R_{1}}},k_{1}\Big{)}}{\mbox{dn}\Big{(}\frac{t}{k_{1}\sqrt{R_{1}}},k_{1}\Big{)}}.$
(5.3)

Now Landen’s transformation (see [13] 22.42) tells us that (5) can be written
as follows

$\cos\frac{\theta\Big{(}\frac{R}{R_{1}}t\Big{)}}{2}=-\mbox{sn}\Big{(}\frac{(1+k^{\prime}_{1})t}{k_{1}\sqrt{R_{1}}},k\Big{)},$
(5.4)

with

$k=\frac{1-k^{\prime}_{1}}{1+k^{\prime}_{1}}=\Big{(}\frac{k_{1}}{1+k^{\prime}_{1}}\Big{)}^{2}.$
(5.5)

Thus assuming, as prescribed by (4.1), that

$\cos\frac{\theta(t)}{2}=-\mbox{sn}\Big{(}\frac{t}{k\sqrt{R}},k\Big{)},$

to satisfy (5.4), we must pick

$R_{1}=\frac{R}{k},$ (5.6)

and choose $k_{1}$ as function of $k$ as in (5.5). It is a standard result
(see [13] 22.42) that

$K=\frac{(1+k^{\prime}_{1})K_{1}}{2},$ (5.7)

with $K_{1}$ as in (1.2) by replacing $k$ with $k_{1}$, hence

$T_{1}=2k_{1}K_{1}\sqrt{R_{1}}=4K\sqrt{R}=2\frac{R_{1}}{R}T,$

which establishes (5.2).

Iterating $j$ times the previous construction, we obtain a circulatory motion
$\theta_{j}(t)$ with period

$T_{j}=2k_{j}K_{j}\sqrt{R_{j}}=2^{j}\frac{R_{j}}{R}T,$ (5.8)

such that the circulatory motion $\theta(t)$ with period $T=2kK\sqrt{R}$ we
started with, can be expressed as follows

$\theta(t)=\sum_{n=-2^{j-1}}^{2^{j-1}-1}\theta_{j}\Big{(}\frac{R_{j}}{R}(t+nT)\Big{)},$
(5.9)

with $R_{j},k_{j}$ inductively defined by

$R_{j}=\frac{R_{j-1}}{k_{j-1}},\;k^{\prime}_{j}=\sqrt{1-k_{j}^{2}}=\frac{1-k_{j-1}}{1+k_{j-1}},j\geq
1,$ (5.10)

with $R_{0}=R,k_{0}=k$.

Lemaître also observes that connecting the $2^{j}$ successive circulatory
motions in (5.9) by segments, one gets a closed Poncelet polygon (see [10] 9,
15, 22). The next step is to establish that
$\lim_{j\to\infty}R_{j}=R_{\infty}<\infty$.

Step 2. We follow [10] (13, 14, 15). From (5.10), one easily establishes that

$\frac{R_{n}}{R_{n+1}}=\frac{\sqrt{R_{n-1}R_{n}}}{\frac{R_{n-1}+R_{n}}{2}},\;n\geq
1.$

Putting

$b_{0}=R<a_{0}=R_{1},$ (5.11)

as shown by Gauss [3], the two sequences

$a_{n}=\frac{a_{n-1}+b_{n-1}}{2},\;b_{n}=\sqrt{a_{n-1}b_{n-1}},\;n\geq 1,$

converge to a common limit denoted by $M(a_{0},b_{0})$, the so-called
arithmetic-geometric mean

$M(a_{0},b_{0})=\lim_{n\to\infty}a_{n}=\lim_{n\to\infty}b_{n},$

and moreover

$\frac{\pi}{2M(a_{0},b_{0})}=\int_{0}^{\frac{\pi}{2}}\frac{d\phi}{\sqrt{a_{0}^{2}\cos^{2}\phi+b_{0}^{2}\sin^{2}\phi}}.$
(5.12)

We refer the reader to [12] 2.3, for a beautiful account of the arithmetic-
geometric mean, and to [2] and [11] for more historical background and updated
research on the subject.

We have

$\frac{R_{1}}{R_{2}}=\frac{\sqrt{RR_{1}}}{\frac{R+R_{1}}{2}}=\frac{b_{1}}{a_{1}},$

and, by induction

$\frac{R_{n}}{R_{n+1}}=\frac{\sqrt{\frac{R_{n-1}}{R_{n}}}}{\frac{1}{2}(1+\frac{R_{n-1}}{R_{n}})}=\frac{\sqrt{a_{n-1}b_{n-1}}}{\frac{a_{n-1}+b_{n-1}}{2}}=\frac{b_{n}}{a_{n}}.$

Hence, using that

$b_{n}^{2}b_{n+1}^{2}\ldots b_{n+j}^{2}=a_{n-1}a_{n}\ldots
a_{n+j-1}b_{n-1}b_{n}\ldots b_{n+j-1},$

we get

$\displaystyle\frac{R_{n}}{R_{n+j}}$
$\displaystyle=\frac{R_{n}}{R_{n+1}}\frac{R_{n+1}}{R_{n+2}}\ldots\frac{R_{n+j-1}}{R_{n+j}}=\frac{b_{n}b_{n+1}\ldots
b_{n+j-1}}{a_{n}a_{n+1}\ldots a_{n+j-1}}$
$\displaystyle=\frac{a_{n-1}b_{n-1}}{b_{n+j}^{2}}=\frac{b_{n}^{2}}{b_{n+j}^{2}}.$

In particular, for $n=0$, we have

$\frac{R}{R_{j}}=\frac{b_{0}^{2}}{b_{j}^{2}},$

which gives

$R_{\infty}=\lim_{j\to\infty}R_{j}=\frac{R}{b_{0}^{2}}\lim_{j\to\infty}b_{j}^{2}=R\Big{(}\frac{M(a_{0},b_{0})}{b_{0}}\Big{)}^{2}.$

Remembering (5.6) and (5.11), $\frac{b_{0}}{a_{0}}=\frac{R}{R_{1}}=k$, we
compute from (5.12) that

$\frac{\pi}{2M(a_{0},b_{0})}=\frac{k}{R}\int_{0}^{\frac{\pi}{2}}\frac{d\phi}{\sqrt{1-k^{\prime
2}\sin^{2}\phi}}=\frac{k}{R}K^{\prime},$

with $K^{\prime}$ as in (1.3), leading to

$R_{\infty}=\Big{(}\frac{\pi}{2kK^{\prime}}\Big{)}^{2}R.$ (5.13)

The next step, the passage to the limit in (5.9), is not carefully justified
in Lemaître’s syllabi [9], [10]. We won’t attempt to justify it, since we have
obtained a different proof of his result in Section 4.

Step 3. This is developed in [10] (16, 17). Remembering (5.8), during the
iteration, the period essentially doubles at each step. Since the limit radius
$R_{\infty}$ is finite, at the limit we obtain the doubly asymptotic motion
(2.2) of the pendulum

$\lim_{j\to\infty}\theta_{j}(t)=\theta_{\infty}(t)=4\arctan
e^{\frac{t}{\sqrt{R_{\infty}}}},$

with $\theta_{\infty}(0)=\pi,\dot{\theta}_{\infty}(0)>0$. Taking _formally_
the limit of (5.9), we get

$\theta(t)=\sum_{n=-\infty}^{\infty}\theta_{\infty}\Big{(}\frac{R_{\infty}}{R}(t+nT)\Big{)},$

i.e., using (5.13), we obtain

$\theta(t)=\sum_{n=-\infty}^{\infty}4\arctan
e^{\frac{\pi}{2kK^{\prime}\sqrt{R}}(t+nT)}.$ (5.14)

However, this series is not convergent. Putting

$q^{\prime}=e^{-\pi\frac{K}{K^{\prime}}},$

and remembering (2.8), we can rewrite (5.14) as follows

$\displaystyle\theta\Big{(}t-\frac{T}{2}\Big{)}=\sum_{n=1}^{\infty}4\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{\frac{\pi t}{2kK^{\prime}\sqrt{R}}}\Big{)}$
$\displaystyle+\sum_{n=1}^{\infty}4\arctan\Big{(}q^{\prime-(n-\frac{1}{2})}\;e^{\frac{\pi
t}{2kK^{\prime}\sqrt{R}}}\Big{)},$
$\displaystyle=\sum_{n=1}^{\infty}4\Big{\\{}\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{\frac{\pi
t}{2kK^{\prime}\sqrt{R}}}\Big{)}-\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{-\frac{\pi t}{2kK^{\prime}\sqrt{R}}}\Big{)}\Big{\\}},$

where each term in the second sum is recomputed modulo $2\pi$ using

$4\arctan\frac{1}{x}=2\pi-4\arctan x,$

which agrees with formula (5) established in Section 4. Finally, and it is the
last step, Lemaître obtains his generalization of Jacobi’s elegantissima by
another application of Landen’s transformation.

Step 4. Following [10] (12, 18), we consider _two identical circulatory
motions_ $\theta_{1}(t)$, with period $T_{1}=2k_{1}K_{1}\sqrt{R_{1}}$, up to a
shift by half a period, going now in _opposite_ directions, i.e. we look at

$\theta(t)=\theta_{1}(t)-\theta_{1}(t-k_{1}K_{1}\sqrt{R_{1}}).$

By a similar computation as in (5) and (5.4), using (4.1) and (4.2), we get

$\cos\frac{\theta(t)}{2}=-\frac{1-k^{\prime}_{1}}{1+k^{\prime}_{1}}\;\mbox{sn}\Big{(}\frac{(1+k^{\prime}_{1})t}{k_{1}\sqrt{R_{1}}},\frac{1-k^{\prime}_{1}}{1+k^{\prime}_{1}}\Big{)},$

i.e. defining

$k=\frac{1-k^{\prime}_{1}}{1+k^{\prime}_{1}}=\Big{(}\frac{k_{1}}{1+k^{\prime}_{1}}\Big{)}^{2}\quad\mbox{and}\quad
R=kR_{1},$ (5.15)

we obtain

$\cos\frac{\theta(t)}{2}=-k\;\mbox{sn}\Big{(}\frac{t}{\sqrt{R}},k\Big{)}.$

Remembering (3.1), this shows that $\theta(t)$ is now an oscillatory motion
with period $T=4K\sqrt{R}$ as in (2.6). Using (5.7) and (5.15), one finds

$T_{1}=2k_{1}K_{1}\sqrt{R_{1}}=\frac{k_{1}}{1+k^{\prime}_{1}}4K\sqrt{R_{1}}=4K\sqrt{R}=T,$
(5.16)

i.e. the oscillatory motion $\theta(t)$ and the circulatory motion
$\theta_{1}(t)$ have now the same period. It is a standard result (see [13],
Example 22.4.5) that

$K^{\prime}=(1+k^{\prime}_{1})K^{\prime}_{1},$ (5.17)

with $K^{\prime}_{1}$ as in (1.3) by replacing $k^{\prime}$ with
$k^{\prime}_{1}$. Hence, using (5.15), we deduce

$2k_{1}K^{\prime}_{1}\sqrt{R_{1}}=\frac{k_{1}}{1+k^{\prime}_{1}}2K^{\prime}\sqrt{R_{1}}=2K^{\prime}\sqrt{R}.$
(5.18)

Using (5.7) and (5.17), we also have

$q^{\prime}_{1}=e^{-\pi\frac{K_{1}}{K^{\prime}_{1}}}=e^{-2\pi\frac{K}{K^{\prime}}}=(q^{\prime})^{2}.$
(5.19)

From (5), using (5.16), (5.18) and (5.19), we obtain

$\theta_{1}\Big{(}t-\frac{T_{1}}{2}\Big{)}=\\\
\sum_{n=1}^{\infty}4\Big{\\{}\arctan\Big{(}q^{\prime 2n-1}\;e^{\frac{\pi
t}{2K^{\prime}\sqrt{R}}}\Big{)}-\arctan\Big{(}q^{\prime 2n-1}\;e^{-\frac{\pi
t}{2K^{\prime}\sqrt{R}}}\Big{)}\Big{\\}},$

hence, by an easy computation, we get

$\displaystyle\theta\Big{(}t-\frac{T}{4}\Big{)}=\theta_{1}\Big{(}t-\frac{T_{1}}{2}+\frac{T_{1}}{4}\Big{)}-\theta_{1}\Big{(}t-\frac{T_{1}}{2}-\frac{T_{1}}{4}\Big{)}$
$\displaystyle=\sum_{n=1}^{\infty}(-1)^{n-1}4\Big{\\{}\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{\frac{\pi
t}{2K^{\prime}\sqrt{R}}}\Big{)}+\arctan\Big{(}q^{\prime
n-\frac{1}{2}}\;e^{-\frac{\pi t}{2K^{\prime}\sqrt{R}}}\Big{)}\Big{\\}},$

which is _the generalized elegantissima_ as in (3.8).

######  Acknowledgements.

_We thank the referees for helpful comments and for pointing out the relevance
of references[2] and [11] to this work._

## References

  * [1] P. Appell, _Sur une interprétation des valeurs imaginaires du temps en Mécanique_ , Comptes Rendus Hebdomadaires des Séances de l’Académie des Sciences, t.87 (30 déc. 1878) 1074-1077.
  * [2] J.M. Borwein and P.B. Borwein, _Pi and the AGM: A Study in Analytic Number Theory and Computational Complexity_ , Canadian Mathematical Society Series of Monographs and Advanced Texts, John Wiley & Sons Inc., New York (1987).
  * [3] C.F. Gauss, _Arithmetisch Geometrisches Mittel_ , in Werke, vol.3 (1799) 361-432.
  * [4] A.G. Greenhill, _Les fonctions elliptiques et leurs applications_ , Editions Jacques Gabay (1895).
  * [5] C.G.J. Jacobi, _Ueber die anwendung der elliptischen transcendenten auf ein bekanntes problem der elementargeometrie_ , Crelle Journal für die reine und angewandte Mathematik 3 (1828) 376-389.
  * [6] C.G.J. Jacobi, _Fundamenta Nova Theoriae Functionum Ellipticarum_ , Bornträger, Königsberg (1829).
  * [7] D. Lambert, _Un atome d’univers_ , La vie et l’oeuvre de Georges Lemaître, Editions Lessius, Bruxelles et Editions Racines, Bruxelles (2000).
  * [8] J. Landen, _An investigation of a general theorem for finding the length of any arc of any conic hyperbola, by means of two elliptic arcs, with some other new and useful theorems deduced therefrom_ , Philos. Trans. Royal Soc. London 65 (1775) 283-289.
  * [9] G. Lemaître, _Leçons de mécanique. Le pendule_. Notes de cours (1955-1956) https://archives.uclouvain.be/ark:/33176/dli000000aGtr#?c=0&m=0&s=0&cv=67&xywh=-479%2C-1%2C5309%2C2650.
  * [10] G. Lemaître, _Leçons de mécanique. Le pendule_. Notes de cours (1956-1957) https://archives.uclouvain.be/ark:/33176/dli000000aGNf#?c=0&m=0&s=0&cv=2&xywh=-605%2C0%2C5279%2C2634
  * [11] D.V. Manna and V.H. Moll, _Landen survey_ , Probability, Geometry and Integrable Systems, MSRI Publications, For Henry McKean’s Seventy-Fifth Birthday (Edited by M. Pinsky and B. Birnir) vol. 55 (2007) 287-319.
  * [12] H. McKean and V. Moll, _Elliptic curves_ , Function theory, geometry, arithmetic, Cambrige University Press (1997).
  * [13] E.T. Whittaker and G.N. Watson, _A course of modern analysis_ , Fifth edition edited and prepared for publication by Victor H. Moll, Cambridge University Press (2021).


# Counting unique molecular identifiers in sequencing using a multitype
branching process with immigration

Serik Sagitov  Mathematical Sciences, Chalmers University of Technology and
University of Gothenburg<EMAIL_ADDRESS>Anders Ståhlberg  Sahlgrenska
Center for Cancer Research, Department of Laboratory Medicine, Institute of
Biomedicine, University of Gothenburg<EMAIL_ADDRESS>Wallenberg
Centre for Molecular and Translational Medicine, University of Gothenburg
Region Västra Götaland, Sahlgrenska University Hospital, Department of
Clinical Genetics and Genomics, Gothenburg, Sweden

###### Abstract

Detection of extremely rare variant alleles, such as tumour DNA, within a
complex mixture of DNA molecules is experimentally challenging due to
sequencing errors. Barcoding of target DNA molecules in library construction
for next-generation sequencing provides a way to identify and
bioinformatically remove polymerase induced errors. During the barcoding
procedure involving $t$ consecutive PCR cycles, the DNA molecules become
barcoded by unique molecular identifiers (UMI). Different library construction
protocols utilise different values of $t$. The effect of a larger $t$ and
imperfect PCR amplifications is poorly described.

This paper proposes a branching process with growing immigration as a model
describing the random outcome of $t$ cycles of PCR barcoding. Our model
discriminates between five different amplification rates $r_{1}$, $r_{2}$,
$r_{3}$, $r_{4}$, $r$ for different types of molecules associated with the PCR
barcoding procedure. We study this model by focussing on $C_{t}$, the number
of clusters of molecules sharing the same UMI, as well as $C_{t}(m)$, the
number of UMI clusters of size $m$. Our main finding is a remarkable
asymptotic pattern valid for moderately large $t$. It turns out that
$\mathrm{E}(C_{t}(m))/\mathrm{E}(C_{t})\approx 2^{-m}$ for $m=1,2,\ldots$,
regardless of the underlying parameters $(r_{1},r_{2},r_{3},r_{4},r)$. The
knowledge of the quantities $C_{t}$ and $C_{t}(m)$ as functions of the
experimental parameters $t$ and $(r_{1},r_{2},r_{3},r_{4},r)$ will help the
users to draw more adequate conclusions from the outcomes of different
sequencing protocols.

Keywords: PCR barcoding, PCR amplification rate, UMI cluster, tree-
bookkeeping, PCR branching process, decomposable multitype Galton-Watson
process, growing immigration.

##  1 Introduction

Figure 1: Overview of an ultrasensitive sequencing protocol (the figure is
created using BioRender.com). The experimental workflow consists of two PCR
steps. In the first, barcoding PCR step, the UMIs are attached to target DNA
molecules during $t$ cycles of amplification. In the second, adapter PCR step,
sequencing adapters are attached to the barcoded DNA sequences during $x$
cycles of amplification. Typical numbers for the parameters $t$ and $x$ are
shown. Finally, libraries are purified and sequenced.

Massive parallel sequencing is implemented in a wide range of applications
within basic research and for clinical practice. Numerous protocols and
technologies are developed to accurately detect and quantify differences in
molecular sequences and detect variants. In cancer management, sequencing is
applied in diagnostics to identify mutations that can be targeted with
specific therapies. The standard massive parallel sequencing techniques can
detect variants with frequencies down to the range of 1–5% [21, 23]. However,
this sensitivity is not sufficient for several emerging applications. For
example, detection of circulating tumor-DNA in liquid biopsies requires
technologies that have the ability to detect variants with frequencies lower
than $0.1$% in clinically relevant samples [1, 7, 11]. The main source of
sequencing noise is due to the polymerase induced errors that occur during
library construction and sequencing [5].

To reduce the sequencing noise, unique molecular identifiers (UMIs), also
known as DNA barcodes, can be used to enable ultrasensitive sequencing [13].
The UMIs typically consist of 8-12 randomised nucleotides that are
experimentally attached to each target DNA molecule. The UMIs are introduced
through a limited number of PCR cycles followed by a general amplification
step (Figure 1). After sequencing, all reads with the same UMI can be tracked
back to the original DNA molecule, allowing to control the polymerase induced
errors and quantification biases.

Experimentally it is challenging to introduce UMIs, since randomised sequences
easily produce non-specific PCR products. To address this challenge, several
barcoding PCR cycles can be applied, which simplifies the experimental
protocol [4, 20]. However, the number of different UMIs and their distribution
is not easily estimated with the increasing number of barcoding PCR cycles,
limiting the use of UMIs in different applications.

Figure 2 depicts the outcome of three perfectly successful barcoding PCR
cycles. For each double-stranded molecule, the upper segment represents a
single-stranded molecule in the direction $5\textprime\to 3\textprime$ and the
lower segment represents a single-stranded molecule in the direction
$3\textprime\to 5\textprime$. The target double-stranded DNA molecule is
placed at the top level, $t=0$. According to Figure 2, the first PCR cycle
produces two double-stranded molecules shown at the level $t=1$: the left
pair, consisting of the target sense molecule plus the antisense molecule with
a reverse primer, and the right pair, consisting of the target antisense
molecule plus the sense molecule with a forward UMI primer. All four single
stranded molecules at $t=1$ are incomplete with one generated UMI. The
complete molecules start appearing at the level $t=2$. With perfect PCR
amplifications, the sizes of UMI clusters grow geometrically as shown by Table
1. For example at $t=3$, as illustrated on Figure 2, there are six UMI
clusters labelled by A, B, C, D, E, F, with the clusters B and D having size
two, and the clusters A, C, E, F being singletons.

Figure 2: Schematic representation for three cycles of barcoding PCR assuming that all fourteen PCR amplifications were successful. The two target DNA strands are shown in blue and red (sense and antisense, respectively). The forward primers with UMI and the reverse primers are shown as the solid and dashed black segments. A single-stranded molecule needs primer sequences at both ends to be complete for the use in downstream PCR. Capital letters mark complete molecules with different UMIs, while non-capitalized letters indicate incomplete molecules with different UMIs. The amplification rates $r_{1},r_{2},r_{3},r_{4},r$, in general, may differ from each other. Table 1: The UMI cluster numbers for the perfect PCR amplifications.  Cluster size | 1 | 2 | 3 | 4 | 5 | Total
---|---|---|---|---|---|---
$t=2$ | 2 | 0 | 0 | 0 | 0 | 2
$t=3$ | 4 | 2 | 0 | 0 | 0 | 6
$t=4$ | 8 | 4 | 2 | 0 | 0 | 14
$t=5$ | 16 | 8 | 4 | 2 | 0 | 30
$t=6$ | 32 | 16 | 8 | 4 | 2 | 62

In practice, PCR amplifications are imperfect and a single-stranded molecule
gets successfully amplified with a probability $r$, which we call the
amplification rate. Despite the attempts to optimise the primers and reaction
conditions, the PCR amplification rates are rarely close to 1. With imperfect
amplifications, Figure 2 and Table 1 should be modified to reflect the fact
that the outcome of each PCR cycle is random. Figure 3 illustrates a possible
realisation of six barcoding PCR cycles, summarised by Table 2. In particular,
at $t=6$, as shown on the bottom rows of Figure 3 and Table 2, A and D form
UMI clusters of size 3, B is a UMI cluster of size 2, and the clusters C, E,
F, G, H, I are singletons.

Figure 3: A possible output of the first six cycles for the barcoding PCR step with imperfect amplifications. The failed PCR amplifications are marked by the dashed arrows and dashed double-arrows.  Table 2: The UMI cluster numbers of imperfect PCR cycles depicted in Figure 3.  Cluster size | 1 | 2 | 3 | 4 | 5 | Total
---|---|---|---|---|---|---
$t=2$ | 0 | 0 | 0 | 0 | 0 | 0
$t=3$ | 1 | 0 | 0 | 0 | 0 | 1
$t=4$ | 4 | 0 | 0 | 0 | 0 | 4
$t=5$ | 4 | 2 | 0 | 0 | 0 | 6
$t=6$ | 6 | 1 | 2 | 0 | 0 | 9

The PCR amplification rate is dependent on both sequence context and sample
quality. In particular, the target DNA molecules are often long, containing
thousands or even millions base pairs. Moreover, some molecules may have
inhibitors attached to the DNA. As a rule, the amplification rates of the
original molecules are smaller compared to the later formed DNA molecules.
Addressing these features, our model discriminates between five amplification
rates $(r_{1},r_{2},r_{3},r_{4},r)$ as indicated by Figures 2 and 3, assuming
that

$0<r_{1},r_{2}\leq r_{3},r_{4}\leq r\leq 1.$ (1)

Here, $r_{1}$ and $r_{2}$ represent the amplification rates of the original
sense and antisense DNA molecules, respectively. The rates $r_{1}$ and $r_{2}$
may differ, since their sequence contexts are different, in so far as two
complementary sequences have different nucleotide sequences. The parameters
$r_{3}$ and $r_{4}$ refer to the amplification rates of the DNA molecules that
have been amplified once. Again, $r_{3}$ and $r_{4}$ may be different since
they have different sequence contexts. We assume that both $r_{3}$ and $r_{4}$
are larger than each of the rates $r_{1}$ and $r_{2}$, since the DNA molecules
become shorter after being amplified once. Finally, we assume that the
shortest amplicon that is amplified exponentially is amplified with the
highest rate $r$ disregarding the difference between the sense and antisense
strands.

This paper introduces and studies a mathematical model for the outcome of the
barcoding PCR experiment starting from a single double-stranded DNA molecule.
Our stochastic model is defined in terms of the five parameters
$(r_{1},r_{2},r_{3},r_{4},r)$ using the framework of the multitype Galton-
Watson processes [6]. The use of branching processes as a stochastic model for
counting the molecules in the repeated PCR amplification cycles is well
established in the literature, see [6, 9, 12, 14, 15] and references therein.
In [17] the branching process approach is applied to the second, adapter PCR
step mentioned in Figure 1. However, to our knowledge, the use of branching
processes is novel for modelling the barcoding PCR step. The multitype
branching process with immigration of this paper is a special example of the
multitype Galton-Watson processes with neutral mutations examined in [3].

##  2 Results

We study the possible outcomes of the barcoding PCR step, starting from a
single double-stranded DNA molecule. Our stochastic model for counting the
unique UMIs is built upon an efficient bookkeeping system for the DNA
barcoding procedure, presented in Section 2.1. We discriminate between six
different types of single-stranded molecules emerging during the barcoding PCR
procedure and introduce a multitype Galton-Watson process [6] with immigration
describing the random process of reproduction of the single-stranded
molecules, see Section 2.2.

We study this model by focussing on $C_{t}$, the number of clusters of
molecules sharing the same UMI, as well as $C_{t}(m)$, the number of UMI
clusters of size $m$. The underlying branching properties of the model yield
recursive formulas for the expected values $\mathrm{E}(C_{t})$ and
$\mathrm{E}(C_{t}(m))$, see Section 2.3. Our main finding, based on the
analysis of the proposed multitype Galton-Watson process, is the following
asymptotic result

$\mathrm{E}(C_{t}(m))/\mathrm{E}(C_{t})\to 2^{-m},\quad m\geq 1,\quad
t\to\infty.$ (2)

The remarkable feature of this relationship is that the limits are the same
irrespectively of the parameter values $(r_{1},r_{2},r_{3},r_{4},r)$. This
approximation may work well already at $t=10$, as illustrated by Figure 6
below.

###  2.1 Tree-bookkeeping system for barcoding PCR

In this section we introduce a convenient bookkeeping system for the outcome
of multiple barcoding PCR cycles. We start by considering the simple case of
perfect PCR amplifications. Compared to the schematic representation of Figure
2, the tree-graph approach of Figure 4 allows us to neatly depict more than
three PCR cycles using the same space. At any tree level $t$, the vertical
branches (lineages) of the tree on Figure 4 are pairwise connected designating
the double-stranded molecules of Figure 2.

Out bookkeeping system distinguishes between six different types of the tree
lineages. At any level $t$ there is exactly one target sense lineage labeled
by 0, and one target antisense lineage labeled by 1. The other four type of
lineages $T_{1},T_{2},T_{3},T_{4}$ are defined by the following lineage
generation rules:

$\displaystyle 0\to T_{1},\quad 1\to T_{2},\quad T_{1}\to T_{3},\quad T_{2}\to
T_{4},\quad T_{3}\to T_{4},\quad T_{4}\to T_{3},$ (3)

where

-

the lineages $0,1,T_{1},T_{2}$ represent incomplete molecules,

-

the lineages $T_{3},T_{4}$ represent complete molecules,

-

newly generated lineages $T_{2}$ and $T_{3}$ represent molecules with a novel
UMI,

-

the lineages $T_{4}$ represent molecules that inherit the UMI of the parental
molecule.

Figure 4: A tree-graph summary of the four cycles of perfect barcoding PCR.
Each vertical lineage represents a molecule that, once appeared, persists over
the succeeding PCR cycles. The dashed and solid branches discriminate between
the incomplete and complete molecules. Labels $0,1,T_{1},T_{2},T_{3},T_{4}$
indicate different types of the molecules with and without UMIs, where the
different levels $t$ of the tree represent the consecutive PCR cycle numbers.

Let $Z_{t}^{(i)}$ stand for the number of the $T_{i}$-lineages at the level
$t$, for $i=1,2,3,4$ and $t\geq 0$. Given

$\displaystyle(r_{1},r_{2},r_{3},r_{4},r)=(1,1,1,1,1),$ (4)

the evolution of the vector
$(Z_{t}^{(1)},Z_{t}^{(2)},Z_{t}^{(3)},Z_{t}^{(4)})$ is deterministic and in
accordance with rule (3) satisfies the following recursions

$\displaystyle Z_{t}^{(1)}$ $\displaystyle=Z_{t-1}^{(1)}+1,$ $\displaystyle
Z_{t}^{(2)}$ $\displaystyle=Z_{t-1}^{(2)}+1,$ $\displaystyle Z_{t}^{(3)}$
$\displaystyle=Z_{t-1}^{(3)}+Z_{t-1}^{(1)}+Z_{t-1}^{(4)},$ $\displaystyle
Z_{t}^{(4)}$ $\displaystyle=Z_{t-1}^{(4)}+Z_{t-1}^{(2)}+Z_{t-1}^{(3)},$

valid for $t\geq 1$ under the initial condition

$\displaystyle Z_{0}^{(1)}=Z_{0}^{(2)}=Z_{0}^{(3)}=Z_{0}^{(4)}=0.$ (5)

From these recursions it is easy to see that $Z_{t}^{(1)}=Z_{t}^{(2)}=t-1$,
implying $Z_{t}^{(3)}=Z_{t}^{(4)}$. Moreover,

$Z_{t}^{(3)}=2Z_{t-1}^{(3)}+t-2,$

so that

$Z_{t}^{(2)}+Z_{t}^{(3)}=2(Z_{t-1}^{(2)}+Z_{t-1}^{(3)})+1,$

yielding

$Z_{t}^{(2)}+Z_{t}^{(3)}=2^{t}-1,\quad t\geq 0.$ (6)

At any given level $t$, we split the set of complete lineages into the
clusters of lineages representing the UMI clusters of molecules sharing the
same UMI. According to our bookkeeping system, there are two different types
of lineage clusters: $T_{2}$-clusters stemming from the $T_{2}$-lineages and
$T_{3}$-clusters stemming from the $T_{3}$-lineages. A $T_{2}$-cluster
consists of the $T_{4}$-lineages, which are the daughter lineages of the stem
$T_{2}$-lineage, with the stem $T_{2}$-lineage not being part of the cluster
as it represents an incomplete molecule. A $T_{3}$-cluster consists of the
stem $T_{3}$-lineage and its daughter $T_{4}$-lineages.

If $C_{t}(m)$ is the number of $T_{2}$ and $T_{3}$-clusters of size $m$
observed at the level $t$ of the lineage tree, then

$C_{t}=C_{t}(1)+\ldots+C_{t}(t-1)$ (7)

gives the total number of lineage clusters at the level $t$. Observe that in
the special case (4), we have

$C_{t}=Z_{t}^{(2)}+Z_{t}^{(3)}-1,$

so that by (6),

$C_{t}=2^{t}-2,\quad t\geq 1.$ (8)

Furthermore, under (4), we get

$C_{t+1}(1)=C_{t}+2,\quad C_{t+1}(m+1)=C_{t}(m),\quad m\geq 2,\quad t\geq 1.$
(9)

The first equality in (9) says that each cluster at the level $t$ produces at
the next level $t+1$ one novel $T_{3}$-cluster, in addition to a new
$T_{2}$-singleton and a new $T_{3}$-cluster generated by a $T_{1}$-lineage.
(By a $T_{2}$-singleton we mean a $T_{2}$-lineage that has not yet produced a
daughter lineage.) The second equality in (9) says that each cluster of size
$m$ at the level $t$ turns into a cluster of size $m+1$ at the next level
$t+1$. By (8) and (9),

$C_{t}(m)=2^{t-m},\quad m\geq 1,\quad t\geq m+1,$ (10)

cf Table 1. Relations (8) and (10) imply an interesting rule of thumb saying
that given (4), at any given $t$, the increase of the cluster size by 1
reduces the number of clusters by half. More precisely,

$C_{t}(m)/C_{t}\to 2^{-m},\quad m\geq 1,\quad t\to\infty.$ (11)

Figure 5: The tree view of the Figure 3. Depending on the corresponding type
of the underlying DNA molecule, the five amplification rates
$r_{1},r_{2},r_{3},r_{4},r$, are assigned to the corresponding horizontal
arrows.

If (4) does not hold, so that some PCR amplifications may fail, the cluster
numbers $C_{t}(m)$ and $C_{t}$ become random. In the next section, we obtain
recursion relations for the corresponding mean values, from which we will
derive (2), a far-reaching extension of the deterministic relation (11).

###  2.2 Multitype branching process with immigration

We are going to write $\xi\sim\text{Ber}(r)$ to say that the random variable
$\xi$ has a Bernoulli distribution with

$\mathrm{P}(\xi=1)=r,\quad\mathrm{P}(\xi=0)=1-r.$

Assuming that the amplification rates $(r_{1},r_{2},r_{3},r_{4},r)$ satisfy
(1), we restate (3) in the form

$\displaystyle 0\to 0+\xi_{0}T_{1},\quad T_{1}\to T_{1}+\xi^{(1)}T_{3},\quad
T_{3}\to T_{3}+\xi^{(3)}T_{4},$ (12) $\displaystyle 1\to 1+\xi_{1}T_{2},\quad
T_{2}\to T_{2}+\xi^{(2)}T_{4},\quad T_{4}\to T_{4}+\xi^{(4)}T_{3},$ (13)

involving six random inputs

$\xi_{0}\sim\text{Ber}(r_{1}),\quad\xi_{1}\sim\text{Ber}(r_{2}),\quad\xi^{(1)}\sim\text{Ber}(r_{3}),\quad\xi^{(2)}\sim\text{Ber}(r_{4}),\quad\xi^{(3)},\xi^{(4)}\sim\text{Ber}(r).$

We clarify these relations by referring to $T_{3}\to T_{3}+\xi^{(3)}T_{4}$,
which says that a $T_{3}$-lineage existing at any given level $t$, at the next
level $t+1$ infallibly reproduces itself and produces a new lineage of type
$T_{4}$ with probability $r$.

The sequence of random vectors

$\displaystyle\\{(Z_{t}^{(1)},Z_{t}^{(2)},Z_{t}^{(3)},Z_{t}^{(4)}):\quad
t=0,1,\ldots\\}$ (14)

forms a Markov chain with the initial state (5). Treating relations (12) and
(13) as the reproduction rules involving four types of individuals along the
dicrete time $t$, we may view this Markov chain as a multitype branching
process with immigration [16]. This is an example of a decomposable multitype
branching process, since the types $T_{3}$ and $T_{4}$ may generate each other
but not the types $T_{1}$ and $T_{2}$. There are two sources of immigration
for this four-type branching process: the $0$-lineage generates the
$T_{1}$-individuals at the rate $r_{1}$, and the $1$-lineage generates the
$T_{2}$-individuals at the rate $r_{2}$. The types $T_{1}$ and $T_{2}$ without
directly communicating with each other, give rise to the types $T_{3}$ and
$T_{4}$ respectively.

The reproduction rules (12) and (13) yield the following recursive relations

$\displaystyle Z_{t}^{(1)}$ $\displaystyle=Z_{t-1}^{(1)}+\xi_{t,0},\quad
Z_{t}^{(2)}=Z_{t-1}^{(2)}+\xi_{t,1},$ (15) $\displaystyle Z_{t}^{(3)}$
$\displaystyle=Z_{t-1}^{(3)}+\sum_{j=1}^{Z_{t-1}^{(1)}}\xi_{t,j}^{(1)}+\sum_{j=1}^{Z_{t-1}^{(4)}}\xi_{t,j}^{(4)},\qquad
Z_{t}^{(4)}=Z_{t-1}^{(4)}+\sum_{j=1}^{Z_{t-1}^{(2)}}\xi_{t,j}^{(2)}+\sum_{j=1}^{Z_{t-1}^{(3)}}\xi_{t,j}^{(3)},$
(16)

involving independent Bernoulli random values

$\xi_{t,0}\sim\text{Ber}(r_{1}),\quad\xi_{t,1}\sim\text{Ber}(r_{2}),\quad\xi_{t,j}^{(1)}\sim\text{Ber}(r_{3}),\quad\xi_{t,j}^{(2)}\sim\text{Ber}(r_{4}),\quad\xi_{t,j}^{(3)},\xi_{t,j}^{(4)}\sim\text{Ber}(r),$

each indicating whether the underlying PCR amplification is successful or not.

In the current setting, the supercritical four-type branching process (14)
could be described in terms of a single type branching process with a growing
immigration. The specific reproduction rules (12) and (13) allow the types
$T_{3}$ and $T_{4}$ to be treated as a single type, say $T$, such that the
type $T$ individuals produce $(1+r)$ offspring on average. In terms of the
process (14), the number $Z_{t}$ of $T$-individuals at time $t$ can be
expressed as the sum

$Z_{t}=Z_{t}^{(3)}+Z_{t}^{(4)},$

and the Markov chain $(Z_{t})_{t\geq 0}$ can be treated as a branching process
with growing immigration. By (16), the number of type $T$ immigrants at time
$t$ is given by the sum $I_{t}=I_{t}^{(1)}+I_{t}^{(2)}$ of two independent
random variables

$\displaystyle I_{t}^{(1)}=\sum_{j=1}^{Z_{t-1}^{(1)}}\xi_{t,j}^{(1)},\qquad
I_{t}^{(2)}=\sum_{j=1}^{Z_{t-1}^{(2)}}\xi_{t,j}^{(2)}$

having binomial distributions $I_{t}^{(1)}\sim\text{Bin}(t,r_{1}r_{3})$ and
$I_{t}^{(2)}\sim\text{Bin}(t,r_{2}r_{4})$. Observe that the generating
function for the number of immigrants
$h_{t}(s)=\mathrm{E}(s^{I_{t}})=\mathrm{E}(s^{I_{t}^{(1)}})\mathrm{E}(s^{I_{t}^{(2)}})$
is computed explicitly

$\displaystyle
h_{t}(s)=(1-r_{1}r_{3}+r_{1}r_{3}s)^{t}(1-r_{2}r_{4}+r_{2}r_{4}s)^{t},\quad
0\leq s\leq 1,\quad t\geq 0.$ (17)

According to Theorem 2b from Section 4 of [18], the long term population size
growth of the supercritical branching process with growing immigration is
regulated by its reproduction rate $(1+r)$:

$\displaystyle(1+r)^{-t}Z_{t}\to W\text{ in }L_{2},\quad t\to\infty.$ (18)

Here, the limit $W$ is a strictly positive random variable, whose Laplace
transform

$\mathrm{E}(e^{-\lambda W})=\prod_{k=1}^{\infty}h_{k}(\lambda(1+r)^{-k})$

is determined by the five amplification rates $(r_{1},r_{2},r_{3},r_{4},r)$ in
terms of the generating functions (17) and the limiting Laplace transform
$\phi(\lambda)$ for the branching process without immigration satisfying the
functional quadratic equation

$\phi((1+r)\lambda)=(1-r)\phi(\lambda)+r\phi^{2}(\lambda).$

The main concern of this paper is not the decomposable multitype branching
process (14) per se, but certain functionals thereof, especially the number of
clusters of size $m$,

$C_{t}(m)=X_{t}(m)+Y_{t}(m),\quad 1\leq m\leq t-1,$

where $X_{t}(m)$ is the number of $T_{2}$-clusters and $Y_{t}(m)$ is the
number of $T_{3}$-clusters of size $m$ at the level $t$. By (16), we have

$\displaystyle Y_{t}(1)$
$\displaystyle=\sum_{j=1}^{Z_{t-1}^{(1)}}\xi_{t,j}^{(1)}+\sum_{j=1}^{Z_{t-1}^{(4)}}\xi_{t,j}^{(4)}+\sum_{j\in\mathbb{Y}_{t-1}(1)}(1-\xi_{t,j}^{(3)}).$
(19)

Here, $\mathbb{Y}_{t}(m)$ is the set of $T_{3}$-lineages which have exactly
$(m-1)$ daughter lineages of the type $T_{4}$ at the level $t$. Furthermore,
again by (16), for $1\leq m\leq t-1$,

$\displaystyle X_{t}(m)$
$\displaystyle=\sum_{j\in\mathbb{X}_{t-1}(m-1)}\xi_{t,j}^{(2)}+\sum_{j\in\mathbb{X}_{t-1}(m)}(1-\xi_{t,j}^{(2)}),$
(20) $\displaystyle Y_{t}(m)$
$\displaystyle=\sum_{j\in\mathbb{Y}_{t-1}(m-1)}\xi_{t,j}^{(3)}+\sum_{j\in\mathbb{Y}_{t-1}(m)}(1-\xi_{t,j}^{(3)}),$
(21)

where $\mathbb{X}_{t}(m)$ is the set of $T_{2}$-lineages which have exactly
$m$ daughter lineages of the type $T_{4}$ at the level $t$, provided $m\geq
1$, while $\mathbb{X}_{t}(0)$ is the set of $T_{2}$-singletons at the level
$t$.

Let $X_{t}(0)$ be the number of $T_{2}$-singletons at the level $t$. (To
illustrate, the example of Figure 5 gives $X_{0}(0)=X_{1}(0)=X_{5}(0)=0$,
$X_{2}(0)=X_{3}(0)=X_{4}(0)=X_{6}(0)=1$.) Then, due to (15) and (16),

$\displaystyle X_{t}(0)$
$\displaystyle=\xi_{t,1}+\sum_{j\in\mathbb{X}_{t-1}(0)}(1-\xi_{t,j}^{(2)}).$
(22)

Since

$Z_{t}^{(2)}=X_{t}(0)+X_{t}(1)+\ldots+X_{t}(t-1),$

the total number of clusters at the level $t$ equals

$\displaystyle
C_{t}=C_{t}(1)+\ldots+C_{t}(t-1)=X_{t}(1)+Y_{t}(1)+\ldots+X_{t}(t-1)+Y_{t}(t-1)=Z_{t}^{(2)}-X_{t}(0)+Z_{t}^{(3)}.$
(23)

In the expression (23) for the total number of clusters $C_{t}$, the
dominating term is $Z_{t}^{(3)}$, which according to (18) is of order
$(1+r)^{t}$. Observe that relation (18) implies that both the mean number of
clusters $\mathrm{E}(C_{t})$ and the standard deviation
$\sqrt{\mathrm{Var\hskip 0.56905pt}(C_{t})}$ are growing proportionally to
$(1+r)^{t}$ as $t\to\infty$.

###  2.3 The expected values

In this section, we denote

$c_{t}=\mathrm{E}(C_{t}),\quad c_{t}(m)=\mathrm{E}(C_{t}(m)),\quad 1\leq m\leq
t-1.$

and show first that

$\displaystyle c_{t}$
$\displaystyle=\alpha(1+r)^{t}+\alpha_{1}t-\alpha_{2}+\alpha_{3}(1-r_{4})^{t}+\alpha_{4}(1-r)^{t},$
(24)

where

$\alpha:=\frac{r_{1}r_{3}+r_{2}r_{4}}{2r^{2}},\quad\alpha_{1}:=r_{2}(1-r_{4}r^{-1}),\quad\alpha_{2}:=r_{1}r_{3}r^{-2}+r_{2}r_{4}^{-1},\quad\alpha_{3}:=r_{2}r_{4}^{-1},\quad\alpha_{4}:=\frac{r_{1}r_{3}-r_{2}r_{4}}{2r^{2}},$

and then derive the main result (2) of this paper. Notice that $\alpha_{4}=0$
if $r_{1}r_{3}=r_{2}r_{4}$, and in the deterministic case (4), relation (24)
turns into (8). Our results concerning the expected values are illustrated by
Figure 6. On the right panel of Figure 6, the four lines, representing
different combinations of the parameter values, almost coincide demonstrating
that asymptotic relation (2) works well already for $t=10$.

Figure 6: Left panel: the expected values of the total number of clusters
$\mathrm{E}(C_{t})$. Middle ($t=5$) and right ($t=10$) panels show the plots
of $\mathrm{E}(C_{t}(m))/\mathrm{E}(C_{t})$ over the cluster sizes
$m=1,\ldots,t-1$. Different colours represent different sets of the parametrs
$(r_{1},r_{2},r_{3},r_{4},r)$: black $(1,1,1,1,1)$, green
$(0.6,0.6,0.7,0.8,0.9)$, red $(0.4,0.4,0.6,0.7,0.9)$, and blue
$(0.3,0.3,0.4,0.5,0.9)$.

Put

$z_{t}^{(i)}=\mathrm{E}(Z_{t}^{(i)}),\quad x_{t}(m)=\mathrm{E}(X_{t}(m)),\quad
y_{t}(m)=\mathrm{E}(Y_{t}(m)).$

The proofs of (24) and (2) rely on the recursive relations

$\displaystyle x_{t}(0)$ $\displaystyle=r_{2}+x_{t-1}(0)(1-r_{4}),$ (25)
$\displaystyle y_{t}(1)$
$\displaystyle=r_{3}z_{t-1}^{(1)}+rz_{t-1}^{(4)}+y_{t-1}(1)(1-r),$ (26)
$\displaystyle x_{t}(m)$
$\displaystyle=x_{t-1}(m-1)r_{4}+x_{t-1}(m)(1-r_{4}),\quad 0\leq m\leq t-1,$
(27) $\displaystyle y_{t}(m)$
$\displaystyle=y_{t-1}(m-1)r+y_{t-1}(m)(1-r),\quad 0\leq m\leq t-1,$ (28)

following from (22), (20), (19), and (21).

Proof of (24). From (15) we obtain

$z_{t}^{(1)}=r_{1}+z_{t-1}^{(1)},\quad z_{t}^{(2)}=r_{2}+z_{t-1}^{(1)},\quad
z_{0}^{(1)}=z_{0}^{(2)}=0,$

so that $z_{t}^{(1)}=r_{1}t$, $z_{t}^{(2)}=r_{2}t$, and therefore by (16),

$\displaystyle z_{t}^{(3)}$
$\displaystyle=z_{t-1}^{(3)}+r_{1}r_{3}(t-1)+rz_{t-1}^{(4)},$ $\displaystyle
z_{t}^{(4)}$ $\displaystyle=z_{t-1}^{(4)}+r_{2}r_{4}(t-1)+rz_{t-1}^{(3)}.$

This yields

$\displaystyle z_{t}^{(3)}+z_{t}^{(4)}$
$\displaystyle=(r_{1}r_{3}+r_{2}r_{4})(t-1)+(1+r)(z_{t-1}^{(3)}+z_{t-1}^{(4)})$
$\displaystyle=(r_{1}r_{3}+r_{2}r_{4})\sum_{j=1}^{t-1}(t-j)(1+r)^{j-1}=2\alpha((1+r)^{t}-rt-1),$

where we used the relation

$\sum_{j=1}^{t-1}(t-j)(1+r)^{j-1}=(1+r)^{t}\sum_{j=2}^{t}(j-1)(1+r)^{-j}=r^{-2}((1+r)^{t}-rt-1).$

On the other hand, we have

$\displaystyle z_{t}^{(3)}-z_{t}^{(4)}$
$\displaystyle=(r_{1}r_{3}-r_{2}r_{4})(t-1)+(1-r)(z_{t-1}^{(3)}-z_{t-1}^{(4)})$
$\displaystyle=(r_{1}r_{3}-r_{2}r_{4})\sum_{j=1}^{t-1}(t-j)(1-r)^{j-1}=2\alpha_{4}((1-r)^{t}+rt-1),$

so that

$\displaystyle z_{t}^{(3)}$
$\displaystyle=\alpha(1+r)^{t}-r_{2}r_{4}r^{-1}t-r_{1}r_{3}r^{-2}+\alpha_{4}(1-r)^{t},$
$\displaystyle z_{t}^{(4)}$
$\displaystyle=\alpha(1+r)^{t}-r_{1}r_{3}r^{-1}t-r_{2}r_{4}r^{-2}-\alpha_{4}(1-r)^{t}.$

By (23), we have

$c_{t}=z_{t}^{(2)}+z_{t}^{(3)}-x_{t}(0),$

and the stated formula (24) follows from the obtained expression for
$z_{t}^{(3)}$ and the next consequence of (25):

$\displaystyle x_{t}(0)=r_{2}r_{4}^{-1}(1-(1-r_{4})^{t}).$ (29)

Proof of (2). Relation (26) implies

$\displaystyle y_{t}(1)$ $\displaystyle=\alpha
r(1+r)^{t-1}-r_{2}r_{4}r^{-1}-\alpha_{4}r(1-r)^{t-1}+y_{t-1}(1)(1-r),$

which entails

$\displaystyle y_{t}(1)$ $\displaystyle=\sum_{j=1}^{t-1}\Bigg{(}\alpha
r(1+r)^{t-j}-r_{2}r_{4}r^{-1}-\alpha_{4}r(1-r)^{t-j}\Bigg{)}(1-r)^{j-1}$
$\displaystyle=\tfrac{1}{2}\alpha(1+r)((1+r)^{t-1}-(1-r)^{t-1})-r_{2}r_{4}r^{-2}(1-(1-r)^{t-1})-\alpha_{4}r(1-r)^{t-1}(t-1)$
$\displaystyle=\tfrac{1}{2}\alpha(1+r)^{t}-r_{2}r_{4}r^{-2}-\alpha_{4}r(1-r)^{t-1}(t-1)+\alpha_{5}(1-r)^{t-1},$

where $\alpha_{5}=r_{2}r_{4}r^{-2}-\tfrac{1}{2}\alpha(1+r).$ Thus,

$y_{t}(1)\sim\tfrac{1}{2}\alpha(1+r)^{t},\quad t\to\infty.$

Using this as the initiation step for the induction over $m$ based on
recursion (28), we find that

$\displaystyle y_{t}(m)\sim 2^{-m}\alpha(1+r)^{t},\quad m\geq 1,\quad
t\to\infty.$

This gives (2) in view of (24) and

$c_{t}(m)=x_{t}(m)+y_{t}(m),\quad 1\leq m\leq t-1,$

where $x_{t}(m)=o((1+r)^{t})$ in accordance with (27) and (29).

Figure 7: Simulation results for the proportions $C_{t}(m)/C_{t}$,
$m=1,\ldots,t-1$ with $t=10$ and the amplification rates
$(r_{1},r_{2},r_{3},r_{4},r)=(0.3,0.3,0.4,0.5,0.9)$. The left panel presents
ten individual simulation results. The right panel summarises 1 000 000
simulation results for the proportions $C_{t}(m)/C_{t}$: the red line connects
the simulation averages, the vertical intervals attached to the red line show
the means $\pm$ one standard deviations obtained from the simulations. Also on
the right panel, the blue line connects the theoretical values for
$\mathrm{E}(C_{t}(m))/\mathrm{E}(C_{t})$, this is the same blue line as on the
right panel of Figure 6.

##  3 Discussion

The number of approaches and applications that use UMIs in sequencing is
rapidly increasing. In cancer diagnostics, the use of UMIs is crucial since it
allows to correct for both polymerase induced errors and amplification biases
[10]. Many sample types and matrices are challenging to analyse due to limited
amounts of DNA and enzymatic inhibitors. In this paper we propose a convenient
bookkeeping system for annotating the emerging UMI clusters during $t$
consecutive barcoding PCR cycles.

The proposed tree based bookkeeping system leads to a branching process model
for the counts $C_{t}(m)$ of the UMI clusters of sizes $m=0,1,\ldots,m-1$. Our
model distinguishes between five PCR amplification rates
$(r_{1},r_{2},r_{3},r_{4},r)$. A key feature of interest for such a model is
the set of proportions $C_{t}(m)/C_{t}$, $m=1,\ldots,t-1$, where
$C_{t}=C_{t}(1)+\ldots+C_{t}(t-1)$ is the total number of the UMI clusters.

The main theoretical finding of this paper, convergence (2), claims that the
ratio between the expected counts $\mathrm{E}(C_{t}(m))/\mathrm{E}(C_{t})$ is
approximately $2^{-m}$ regardless of the underlying parameters
$(r_{1},r_{2},r_{3},r_{4},r)$. It was demonstrated that this approximation
formula works well even for moderately large values of $t$, see the right
panel of Figure 6. We hypothesise a biologically more relevant asymptotic
result

$\mathrm{E}(C_{t}(m)/C_{t})\to 2^{-m},\quad m\geq 1,\quad t\to\infty.$ (30)

To address this hypothesis, a simulation study based on our model was
performed by Hongui Zhan and Yizhe Gu, two master students at the Chalmers
University of Technology. Their simulation results summarised by Figure 7,
support the approximation formula (30) for the moderate value of $t=10$ and a
particular choice of the amplification rates
$(r_{1},r_{2},r_{3},r_{4},r)=(0.3,0.3,0.4,0.5,0.9)$.

Our model of the barcoding PCR step uses five different amplification rates as
the key model parameters $(r_{1},r_{2},r_{3},r_{4},r)$. In sequencing,
amplification rates are rarely assessed and there is no general method to
determine $(r_{1},r_{2},r_{3},r_{4})$. In the framework of quantitative PCR,
the overall amplification rate $r$ can be assessed by standard curves [22].
The amplification rate $r$ varies between assays due to different sequence
context as well as between samples due to sample inhibition [2, 19]. The
amplification rate also decreases during the last PCR cycles, when reagents
become sparse. Scientist working with PCR and sequencing are experimentally
used to the fact that some samples and sequences suffer from poor
amplification rates or that the original DNA molecules never become amplified.
In future studies, it will be important to verify our model with experimental
data to estimate the importance of different model parameters. Such a verified
model will be valuable in development of improved sequencing protocols and our
ability to detect and quantify individual DNA molecules with single nucleotide
resolution.

### Acknowledgements

We are grateful to professor Peter Jagers for bringing us together and for
fruitful discussions. This research was partially funded by Region Västra
Götaland, Sweden; Swedish Cancer Society (20-1098); Swedish Research Council
(2020-01008); Swedish Childhood Cancer Foundation (MTI2019-0008 and
2020-0007); the Swedish state under the agreement between the Swedish
government and the county councils, the ALF-agreement (ALFGBG-965065);
Sweden’s Innovation Agency and the Sjöberg Foundation.

## References

  * [1] Andersson, D., Kristiansson, H., Kubista, M., Ståhlberg, A. Ultrasensitive circulating tumor DNA analysis enables precision medicine: experimental workflow considerations. Expert Rev. Mol. Diagn. 21 (2021) 299–310.
  * [2] Bar T, Ståhlberg A, Muszta A, Kubista M. Nucleic Acids Res. Kinetic Outlier Detection (KOD) in real-time PCR. 31(2003 ) e105.
  * [3] Bertoin J. The structure of the allelic partition of the total population for Galton-Watson processes with neutral mutations. Ann. Probab. 37 (2009) 1502–1523.
  * [4] Cohen, J.D. et al. Detection and localization of surgically resectable cancers with a multi-analyte blood test. Science 359 (2018) 926–930.
  * [5] Filges, S., Yamada, E., Ståhlberg A, Godfrey T.E. Impact of polymerase fidelity on background error rates in next-generation sequencing with unique molecular identifiers/barcodes. Sci. Rep. 9 (2019) 3503.
  * [6] Haccou, P., Jagers, P., Vatutin, V. Branching processes: variation, growth, and extinction of populations. Cambridge Studies in Adaptive Dynamics, Cambridge University Press (2005).
  * [7] Ignatiadis, M., Sledge, G.W., Jeffrey, S.S. Liquid biopsy enters the clinic—Implementation issues and future challenges. Nat. Rev. Clin. Oncol. 18 (2021) 297–312.
  * [8] Jagers,P. Branching processes with biological applications. John Wiley (1975).
  * [9] Jagers, P., Klebaner, F.C. Random variation and concentration effects in PCR. J. Theoret. Biol. 224 (2003) 299–304
  * [10] Johansson G, Kaltak M, Rîmniceanu C, Singh AK, Lycke J, Malmeström C, Hühn M, Vaarala O, Cardell S, Ståhlberg A. Ultrasensitive DNA Immune Repertoire Sequencing Using Unique Molecular Identifiers. Clin Chem. 66 (2020) 1228–1237.
  * [11] Heitzer, E., Haque, I.S., Roberts, C.E., Speicher, M.R. Current and future perspectives of liquid biopsies in genomics-driven oncology. Nat. Rev. Genet. 20 (2019) 71–88.
  * [12] Kimmel, M., Alexrod, D. E. Branching processes in biology. Interdisciplinary Applied Mathematics Vol. 19 (2002).
  * [13] Kinde, I., Wu, J., Papadopoulos, N., Kinzler, K.W., Vogelstein, B. Detection and quantification of rare mutations with massively parallel sequencing. Proc. Natl. Acad. Sci. USA 108 (2011) 9530–9535.
  * [14] Krawczak, M., Reiss, J., Schmidtke, J., Rösler, U. Polymerase chain reaction: replication errors and reliability of gene diagnosis. Nucleic Acids Res. 17 (1989) 2197–2201.
  * [15] Lalam, N. Estimation of the reaction efficiency in polymerase chain reaction. J. Theor. Biol. 242 (2006) 947–953
  * [16] Mode, C.J. Multitype branching processes: theory and applications. American Elsevier (1970).
  * [17] Pflug, F.G., von Haeseler, A. TRUmiCount: correctly counting absolute numbers of molecules using unique molecular identifiers. Bioinformatics 34 (2018) 3137–3144.
  * [18] Rahimov, I. Homogeneous branching processes with non-homogeneous immigration. Stochastics and Quality Control 36 (2021) 165–183.
  * [19] Ruijter JM, Barnewall RJ, Marsh IB, Szentirmay AN, Quinn JC, van Houdt R, Gunst QD, van den Hoff MJB. Efficiency Correction Is Required for Accurate Quantitative PCR Analysis and Reporting. Clin. Chem. 67 (2021) 829–842.
  * [20] Ståhlberg, A., Krzyzanowski, P.M., Jackson, J.B., Egyud, M., Stein, L., Godfrey, T.E. Simple, multiplexed, PCR-based barcoding of DNA enables sensitive mutation detection in liquid biopsies using sequencing. Nucleic Acids Res, 44 (2016) 105.
  * [21] Stead, L.F., Sutton, K.M., Taylor, G.R., Quirke, P., Rabbitts, P. Accurately identifying low‐allelic fraction variants in single samples with next‐generation sequencing: Applications in tumor subclone resolution. Hum Mutat, 34 (2013) 1432–1438.
  * [22] Svec D, Tichopad A, Novosadova V, Pfaffl MW, Kubista M. How good is a PCR efficiency estimate: Recommendations for precise and robust qPCR efficiency assessments. Biomol. Detect. Quantif. 3 (2015) 9–16.
  * [23] Xu, H., DiCarlo, J., Satya, R.V., Peng, Q., Wang, Y. Comparison of somatic mutation calling methods in amplicon and whole exome sequence data. BMC Genomics 15 (2014) 244.