On a Linear Polymerisation Process
Sergei Zuyev
Chalmers University of Technology and University of Gothenburg, Department of
Mathematical Sciences, 412 96 Gothenburg, Sweden. sergei.zuyev@chalmers.se
http://www.math.chalmers.se/~sergei
Abstract. We consider a Markov fusion-disaggregation process of ‘mole-
cules’ which are linear arrangements of ‘atoms’. Molecules may fuse at
the ending atoms by creating a bond in between at a given rate. Each
bond joining two atoms in a molecule can break so that a molecule disag-
gregates into two shorter molecules at another rate. We give en explicit
expression for the stationary distribution of the number of molecules and
the conditional distribution of their lengths given their number. We show
an asymptotic normality of the distribution of the number of molecules
when the number of atoms in the system grows.
Keywords: Fusion-Disaggregation Process · Reversibility · Stationary
Distribution · Markov Chain
1 Introduction
The model we consider in this article has origins in colloidal chemistry. Colloid
is a mixture of molecules which stay evenly dispersed in the volume neither
settling to the botton nor floating to the surface of the container. Although the
geometry of molecules can be rather complex, some materials, like lacquers, have
long chains of atoms in a linear molecules which combine to form a net when the
lacquer sets. The molecule bonds may also break causing emergence of shorter
molecules, this is the main mechanism of degradation of the material. The two
competing mechanisms: fusion and disaggregation compete with each other and
should be in balance in a stable system.
To mathematically model a colloid with linear molecules, we consider a sys-
tem consisting of N elements called atoms. n ≥ 2 atoms may form a molecule
of size n (or an n-molecule) which is a linearly arranged set of atoms linked by
bonds. We will also call a singular atom a 1-molecule.
Time evolution of the system is described by a homogeneous continuous time
Markov chain (MC) determined by two competing processes: fusion (or cluster-
ing) and disaggregation (or severance). Two molecules close to each other may
combine by creating a bond between terminal atoms. Ignoring spatial considera-
tions, each n-molecule fuses to an r molecule at rate λ to form an n + r molecule.
In other words, each pair of molecules fuses at the rate 2λ > 0 independently of
their sizes. Each bond breaks with the rate 1, so that an n-molecule disaggregates
to two shorter molecules at rate n − 1.
Copyright © 2020 for this paper by its authors. Use permitted under
Creative Commons License Attribution 4.0 International (CC BY 4.0).
It is clear that the longer molecule is, the more bonds it contains, the more
is its rate of disaggregation. On the other hand, the more molecules are in the
system, the more is the number of their terminal atoms which may fuse, so the
two competing mechanisms should balance each other in the long run leading to
a stationary distribution. Since the system is finite, the stationary distribution
is also its asymptotic distribution which we aim to characterise.
The similar (and even more general) model has been considered by Frank
Kelly who obtained a product form of the stationary distribution [3, Th.8.1].
Unfortunately, no details on how the expression was obtained is given there, the
proof consists in just checking that the solution satisfies the so-called Detailed
Balance Equations which is a necessary and sufficient condition for the reversibil-
ity of the corresponding Markov chain. In this short article, we fill this gap by
explicitly describing a procedure which can be used to establish reversibility
and complement Kelly’s result with finding exact asymptotics of the stationary
distribution when the number of atoms in the system grows.
In the next section we give the necessary definitions and describe the model
we are working with. In Section 3 we find an explicit expression for the stationary
distribution of the model and then study its asymptotic properties when the
number of atoms grows in Section 4.
2 Markov Model and Reversibility
The states of the system can be described by a vector n = (n1 , n2 , . . . , nN ),
PN
k=1 knk = N with nk being the number of k-molecules in the system. Let
X(t), t ≥ 0 denotes the corresponding MC in this state space and ek stands
for the k-th coordinate vector of RN . If τ is the time of the first change of the
system state, then X(τ ) − X(0) with positive probability may assume only the
following values:
1. vij = −ei − ej + ei+j , i 6= j, when two unequal size i- and j- molecules fuse.
The intensity of this jump is 2λni nj .
2. vii = −2ei +e2i , when two i-molecules fuse. The intensity is then λni (ni −1).
3. −vij , i 6= j, when an i + j-molecule breaks into two uneven parts. There are
two bonds in each i+j-molecule which breaking results in i- and j-molecules,
so the intensity of this transition is 2ni+j .
4. −vii , when a 2i-molecule breaks in the middle into two i-molecules, the
intensity of this is n2i .
Not all these transitions are possible from all the states: the ones which would
lead to negative values of any coordinate are forbidden.
Obviously, any state of the system is attainable from any other state by a
sequence of positive-intensity jumps: for instance, by disaggregating to singular
atoms and then fusing to the other configuration. This means the associated MC
is irreducible and hence geometrically ergodic, i.e. there exists a unique station-
ary distribution to which convergence of the distribution at a time t happens
exponentially fast whatever is the starting state X(0), see, e.g., [2].
An important property some MCs possess is time reversibility. Recall that
a continuous time irreducible MC X(t) with states {1, 2, 3, . . . } and transition
rates q(i, j) is called reversible if (X(t1 ), . . . , X(tn )) has the same distribution as
(X(t0 −t1 ), . . . , X(t0 −tn )) for all t0 , t1 , . . . , tn . Then X(t0 −t) is a distributional
copy of X(t), hence the name.
A MC is reversible if and only if there is a probability distribution (π(i))
such that the following Detailed Balance Equations are satisfied for all pairs of
states i, j:
π(i)q(i, j) = π(j)q(j, i), (1)
see, e.g., [3, Th. 1.2]. Such (π(i)) is then the stationary distribution of the MC.
The Kolmogorov’s criterion of reversibility states that MC is reversible if
and only if for any cyclic sequence of states i, i1 , i2 , . . . , in , i, the product of the
transition intensities on one and opposite directions coincide:
q(i, i1 )q(i1 , i2 ) . . . q(in , i) = q(i, in )q(in , in−1 ) . . . q(i1 , i), (2)
see, e.g., [3, Th. 1.8]. Both (1) and (2) allow to express the stationary distribution
of a reversible MC in the form:
q(i0 , i1 )q(ii , i2 ) . . . q(in , j)
π(j) = π(i0 ) (3)
q(j, in )q(in , in−1 ) . . . q(i1 , i0 )
for any chosen state i0 and a sequence of states (a path) i0 , i1 , . . . , in , j for which
the product of the intensities
P above is non-zero. The value of π(i0 ) itself is then
found from normalisation: j π(j) = 1.
3 Stationary Distribution
In this section we give a closed form expression for the stationary distribution of
the MC modelling the molecular system described in the previous section. For
this, we first establish that the MC is reversible. Property (2) is hard to check.
Instead, we evaluate expressions (3) and see if the obtained sequence (π(j))
satisfies the detailed balance equations (1). Then the MC is reversible and such
(π(j)) is its stationary distribution. This the way we prove the main result of
this section.
Theorem 1. The Markov chain describing evolution of the system of N atoms
is reversible. Its stationary distribution satisfies, for n0 = (N, 0, . . . , 0),
N!
π(n) = π(n0 )λL(n) Q , (4)
k≥1 nk !
P
where L(n) = k≥2 (k − 1)nk is the number of links in the state n.
Proof. Take the state n0 = (N, 0, . . . , 0), where all atoms are separated, to ex-
press the probability of all other states in the form of (3). We choose the fol-
lowing path to reach a state n = (n1 , . . . , nN ) 6= n0 from n0 . Let m, m ≥ 2 be
the length of the longest molecule in n, i.e. all nm+1 = · · · = nN = 0. At first,
the atoms join one by one to create an m-molecule, then, if nm > 1, the second
m-molecule, the third, etc., until all m-molecules are created. Then it comes the
turn of m − 1-molecules (or the next smaller size molecules if nm−1 = 0 until
all the molecules of that size are created. Then, the next smaller size molecules,
etc., until all k-molecules, k ≥ 2, are created. n1 = n0 − 2e1 + e2 . The next
one is n2 = n0 − e1 − e2 + e3 and so on. The state with one m-molecule is
obtained by adding −e1 − em−1 + em to the previous state. This creation process
is exemplified on Figure 1. Each state can be represented uniquely as an array
of segments representing molecules spanning the grid points {1, 2, . . . , N }. The
k-molecules are represented by disjoint segments [j, j + m − 1] arranged by their
size, the longest on the left.
Fig. 1. The path between the starting and the terminal states in a system with 9
atoms. The molecules are represented by segments ordered by their size. The transition
intensities are shown on the sides.
N
The rate at which the first link is established is λN (N − 1) since there are
2 pairs of atoms in configuration (N, 0, . . . , 0) and each pair joins at rate 2λ.
The next atom joins the newly formed 2-molecule at the rate 2λ(N − 2), the
next one at the rate 2λ(N − 3) and so on, the last atom to complete the first
m-molecule joins at the rate 2λ(N − m + 1). Now there are N − m separate
atoms in the system and one m-molecule. The product of the intensities of all
the above steps is then
N!
2m−2 λm−1 N (N − 1) . . . (N − m + 1) = 2m−2 λm−1 .
(N − m)!
The intensities of construction of the next molecule just repeat the above ex-
pressions with N replaced by N − m. For instance, if nm > 1, the product of
intensities of all the steps until construction of the second m-molecule is given
by
N! (N − m)! N!
2m−2 λm−1 × 2m−2 λm−1 = 22(m−2) λ2(m−1)
(N − m)! (N − 2m)! (N − 2m)!
After all nm molecules are created, the product is
N!
2(m−2)nm λ(m−1)nm .
(N − mnm )!
Now there are N −mnm atoms in the system and the construction of (m−1)-
molecules proceed P similarly. Once all k-molecules for k ≥ 3 are created, there
are N 0 = N − k≥3 knk atoms. The first 2-molecule is created with intensity
λN 0 (N 0 − 1), the second – λ(N 0 − 2)(N 0 − 3) etc. until the last n2 -nd one with
intensity λ(N 0 − 2(n2 − 1))(N 0 − 2(n2 − 1) − 1). The product of intensities of the
transitions involving creation of all 2-molecules is thus
P P
n2 N 0! n2
(N − k≥3 knk )! n2
(N − k≥3 knk )!
λ =λ =λ .
(N 0 − 2n2 )!
P
(N − k≥2 knk )! n1 !
This gives the following product of the intensities of all the transitions from n0
to n:
P P N! N!
2 k≥3 (k−2)nk λ k≥2 (k−1)nk = 2I(n) λL(n) , (5)
n1 ! n1 !
P
where I(n) = k≥3 (k − 2)nk is the number ofPinner atoms (i.e. connected to
two other atoms in the molecules) and L(n) = k≥2 (k − 1)nk is the number of
links.
The reversed moving along the same path consists, first, in one of n2 2-
molecules to break if n2 ≥ 1. This happens at the rate n2 , because there are
that many links in all 2-molecules. The next 2-molecule breaks at the rate n2 −1,
etc. The last one breaks at the rate 1, so that the product of intensities of all
these breaks is n2 !. If n3 ≥ 1, the next step is a link of a 3-molecule to break,
this happens at rate 2n3 , and then the second link of the same, now 2-molecule,
breaks at rate 1. Now the number of 3-molecules in the system is decreased by
one, so if there are still 3-molecules, the next one disintegrates into single atoms
at rate 2(n3 − 1). Thus, the product of intensities of the steps from n to the state
where all 2- and 3-molecules are disintegrated into atoms equals n2 ! 2n3 n3 !. After
all the molecules of the sizes smaller than k have disintegrated into single atoms,
the first link in a k-molecule breaks at rate (k −1)nk , the second link in the same
molecule breaks at rate 2 (since there are two ends) and so on until the last one
breaks at rate 1. So the first disintegrated k-molecule contributes 2k−2 nk to the
product of intensities. Similarly, the second contributes 2k−2 (nk − 1) and so on.
Hence, the product of the intensities of all the steps from n back to n0 is given
by Y Y
2nk (k−2) nk ! = 2I(n) nk ! . (6)
k≥2 k≥2
This allows us to define the quantities
2I(n) λL(n) nN1!! N!
π(n) = π(n0 ) I(n) Q = π(n0 )λL(n) Q . (7)
2 k≥2 kn ! k nk !
The next step consists in verifying that such defined quantities satisfy the
detailed balance equations (1) for the four possible type of jumps described at
the beginning of Section 2.
Indeed, for the transition n ↔ n + vij , i 6= j, the detailed balance equation
becomes
π(n) 2λni nj = π(n + vij ) 2(ni+j + 1).
Since L(n + vij ) = L(n) + 1, this reduces to
N! ni nj N !
λL(n) Q 2λni nj = λL(n)+1 Q 2(ni+j + 1)
n
k k ! (n i+j + 1) k nk !
which is true. Similarly one can easily verify that for transition n ↔ n + vii , the
detailed balance equation
π(n) λni (ni − 1) = π(n + vii ) (n2i + 1)
also holds. Thus by [3, Th. 1.2], the Markov chain is reversible and its stationary
distribution is given by (4).
P
Corollary 1 Denote by M (n) = k nk the total number of molecules (including
atoms) in the state n. Under the stationary regime, the probability mass function
(p.m.f.) of M is given by
λ−m N − 1
−1
π(M (n) = m) = Z (λ, N ) , m = 1, . . . , N, (8)
m! m − 1
the normalising constant being
N
λ−m N − 1
X 1 1
Z(λ, N ) = = L (−1/λ), (9)
m=1
m! m−1 λN N −1
where Lαn (x) is the generalised Laguerre polynomial of degree n, see, e.g., [4].
The distribution is unimodal with the mode at the point
& p '
∗ −(1 + λ) + (1 + λ)2 + 4λN p
m = ∼ O( N/λ) as N → ∞. (10)
2λ
Its Laplace transform is
Z(et λ, N ) e−t L1N −1 (−e−t /λ))
LM (t) = E e−tM = = , (11)
et Z(λ, N ) L1N −1 (−1/λ)
the first two moments are given by
l2 l2 l1 l3 − l22
EM = 1 + ; var M = + , (12)
λl1 λl1 λ2 l12
where lk = LkN −k (−1/λ), k = 1, 2, 3.
The p.m.f. of the number M of molecules under the stationary regime are shown
on Figure 2.
0.15
0.4
●
λ=20
●
●
●
0.3
●
●
λ=20 N=10 000
0.10
●
● N=100 ●
λ=5 ●
●
●
●
●
pmf
pmf
0.2
●
●
●
λ=5
●● ●
●
● ●
λ=1 ●
● ● ●
●
● ●●
0.05
● ●●
●
● ●
● ●
●
λ=1
● ● ●
●
● ● ● ●
●
0.1
λ=0.1
● ● ●
● ● ●
● ● ●
● ● ● ● ● ● ●●
●
●●
● ● ●
λ=0.01
● ● ●
●
● ● ● ●
●
λ=0.1
● ● ● ● ● ●
● ● ● ●
● ●
● ● ● ● ●
● ● ●
● ● ● ● ●
● ● ● ● ● ●
● ● ● ● ●
● ● ● ● ●
● ● ● ● ●
● ● ● ● ● ● ●
● ● ● ● ●
● ● ● ● ●
● ● ● ● ● ● ●
● ● ● ● ● ●●
● ● ●● ●
● ● ●● ● ● ● ●
● ● ● ● ● ● ●
0.00
● ● ● ● ●● ● ● ● ●
● ●●
● ●
0.0
● ● ●●●● ● ●●● ● ● ● ●● ●
●
● ● ●
●●●
●● ● ●● ● ● ●
● ●
● ●●
●
0 20 40 60 80 100 0 100 200 300 400
m m
Fig. 2. Top prot: the p.m.f. of the number M of molecules under the stationary regime
for the cases N = 100 and N = 10 000. Bottom plot: the histograms of M in a computer
simulation.
P
Proof. Taking into account that k knk = N , expression (4) gives that
X N!
π(M (n) = m) = π(n0 )λN −m Q . (13)
n:
P k nk !
P k nk =m
k knk =N
Assume the system contains exactly m molecules. Number them all somehow
from 1 to m and let (l1 , . . . , lm ) be their sizes. If we represent each k-molecule
by a closed segment of length k − 1, then they can be put from left to right
starting from the first to the last segment disjointly on [1, N ] so that each k-
molecule spans exactly k integer points, P like on Figure 1, but not ordered by
their lengths, but by their number. Since k knk = N , all the integer points of
[1, N ] are covered and there are exactly m−1 open segments of the form (i, i+1)
which separate the molecules. Thus all the sequences (l1 , . . . , lm ) of the sizes of
ordered molecules can be obtained by throwing away
N −1
m − 1 open unit segments
out of total N − 1 of [1, N ] \ N, so there are m−1 of these.
Pm
Alternatively, the sequences (l1 , . . . , lm ), li ≥ P satisfying P i=1 li = N can
1
be produced from n = (n1 , n2 , . . . ) satisfying k nk = m, k knk = N by
assigning n1 1’s, n2 2’s, etc. Q to places 1, . . . , m, so that each such n generates
the multinomial number m!/ k nk ! of such distinct sequences. Therefore,
X m! N −1
Q =
P k nk ! m−1
n:
P k nk =m
k knk =N
yielding
N! N − 1
π(M (n) = m) = π(n0 )λN −m .
m! m − 1
PN
The normalising condition m=1 π(M (n) = m) = 1 allows to express
" N #−1
X N ! N − 1
π(n0 ) = λN −m (14)
m=1
m! m − 1
that gives for M (n) distribution (8) after cancelling out N !λN .
Notice that
N −m
π(M (n) = m + 1) = π(M (n) = m).
λm(m + 1)
The multiplicative coefficient given by the fraction is monotonely decreasing and
becomes less than 1 for m = m∗ as in (10). Thus the distribution is unimodal.
The expression for the Laplace transform is immediate from (9) and hence
the moments follow from its Taylor expansion at t = 0.
The obtained expression (14) for π(n0 ) plugged into (4) together with the
distribution (8) leads to the main result of this section:
Theorem 2. The stationary distribution π is given by
λ−M (n)
π(n) = Z −1 (λ, N ) Q (15)
k nk !
with Z(λ, N ) as in (9). The conditional distribution of the molecule lengths given
their total number is
m
m! n1 ,...,nN
π(n M (n) = m) = N −1 Q = N −1
(16)
m−1 k nk ! m−1
which is the symmetricPMultinomial distribution Mult(m; N −1 , . . . , N −1 ) condi-
tioned on the set {n : k knk = N }.
4 Asymptotic Distribution for Molecule Number
In this section we consider the limit distributions for the system when the number
N of molecules grow. In the previous section we expressed the Laplace transform
of the stationary distribution of the total number M = MN of molecules in the
system with N atoms in terms of the generalised Laguerre polynomials. In [1],
the following asymptotics in n for fixed parameters α and z was obtained:
√ !
α e−z/2 e2 (n+1)z 1
Ln (−z) = √ 1+O √ . (17)
2 π z 1/4+α/2 (n + 1)1/4−α/2 n+1
Substituting this expression into (11) with n = N − 1 and z = 1/λ, we get
n t z √ o 1
LMN (t) = exp − − e−t − 1 + 2 N z e−t/2 − 1
1+O √ . (18)
4 2 N
√ p
Introduce aN = N z = N/λ and denote
M N − aN
µN = p .
aN /2
For the Laplace transform of µN we have the identity
√ t
LµN (t) = et 2aN LMN p
aN /2
and it is straightforward to check using (18) that limN →∞ LµN (t) = t2 /2. Thus
we have proved the following Central limit theorem:
Theorem 3. The stationary distributions of the number MN of molecules in
the system with N atoms when properly centred and scaled weakly converge to
the standard Normal law:
MN − aN p
p =⇒ N (0, 1), where aN = N/λ. (19)
aN /2
Note that from (12) pand (17), both E MN and var Mn under the stationary
regime have order of N/λ, so the centering and scaling in (19) is as could be
expected.
Using the standard methods, a local limit theorem can also be established
which we give without a proof, see illustrative Figure 3:
0.4
0.4
● ● ●●● ●
●
● ●
● ● ●
● ● ●
● ●
●
●
●
●
● ●
●
●
●
●
●
0.3
0.3
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
● ●
0.2
0.2
●
●
● ●
●
●
● ●
●
●
●
● ●
●
● ●
●
●
λ=5
0.1
0.1
●
●
λ=5
●
●
● ●
● ●
● ●
●
N=10^4 ● ● ●
●
● N=10^6 ●
●
● ● ● ●
● ●
● ●
● ● ● ●
● ●
● ●
● ● ●●
● ●● ●●
● ●● ●●
● ●● ●●
0.0
0.0
● ● ●●
−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3
x x
Fig. 3. The scaled p.m.f. (the dots) of the number M of molecules under the stationary
regime for the case N = 104 (left plot) and N = 106 (right plot) and the density of the
standard normal distribution (the curve).
√
Theorem 4. The value π(MN = m) for m = O( N ) is approximated when
p → ∞ by the density of the Normal N (aN , aN /2) distribution, where aN =
N
N/λ, i.e.
1 n (m − a )2 o
N
π(MN = m) = √ exp − + o(N −1/4 ). (20)
πaN aN
5 Discussion
The model we considered in this paper allows for generalisation in various direc-
tions. First of all, the linear structure of the molecules we impose can be dropped,
instead, terminal atoms to where the molecules can attach (the molecule end-
points as now) may be defined, possibly of varying accommodation capacities, see
[3, Sec. 8]. It is interesting if explicit expressions and the asymptotic behaviour
can also be established for the stationary distribution in these cases.
Another direction is to allow some molecules to disappear and to come into
the system. This could model the diffusion, its parameters, for instance, may
depend on the molecule size. One may expect a Poisson distribution to appear
in the stationary regime in such an open system, as in [3, Th. 8.2].
Acknowledgement
The author thanks his Master students: Lydia Andersson, Karin Furufors, Jes-
sica Gilmsjoe and Emelie Lööf, for the simulation studies which confirmed the
theoretical results reported here.
References
1. Borwein, D., Borwein, M., Crandall, R.: Effective Laguerre asymptotics. SIAM J.
Number. Anal. 46(6), 3285–3312 (2008)
2. Diaconis, P., Miclo, L.: On quantitative convergence to quasi-stationarity. Ann. Fac.
Sci. Toulouse Math. 24(4), 973–1016 (2015)
3. Kelly, F.: Reversibility and stochastic networks. Cambridge University Press (1979)
4. Rusev, P.: Classical orthogonal polynomials and their associated functions in com-
plex domain. Mann Drinov Academic Publishing House, Sofia (2005)