# Optimal reaction coordinates and kinetic rates from the projected dynamics
of transition paths

Line Mouaffac  Sorbonne Université, Musée National d’Histoire Naturelle, UMR
CNRS 7590, Institut de Minéralogie, de Physique des Materiaux et de
Cosmochimie, IMPMC, F-75005 Paris, France  Karen Palacio-Rodriguez  Sorbonne
Université, Musée National d’Histoire Naturelle, UMR CNRS 7590, Institut de
Minéralogie, de Physique des Materiaux et de Cosmochimie, IMPMC, F-75005
Paris, France  Fabio<EMAIL_ADDRESS>Sorbonne Université, Musée National d’Histoire Naturelle, UMR CNRS 7590,
Institut de Minéralogie, de Physique des Materiaux et de Cosmochimie, IMPMC,
F-75005 Paris, France

###### Abstract

Finding optimal reaction coordinates and predicting accurate kinetic rates for
activated processes are two of the foremost challenges of molecular
simulations. We introduce an algorithm that tackles the two problems at once:
starting from a limited number of reactive molecular dynamics trajectories
(transition paths), we automatically generate with a Monte Carlo approach a
sequence of different reaction coordinates that progressively reduce the
kinetic rate of their projected effective dynamics. Based on a variational
principle, the minimal rate accurately approximates the exact one, and it
corresponds to the optimal reaction coordinate. After benchmarking the method
on an analytic double-well system, we apply it to complex atomistic systems:
the interaction of carbon nanoparticles of different sizes in water.

##  1 Introduction

Physico-chemical transformations such as phase transitions, chemical
reactions, and biomolecular conformational changes are characterized by
metastable states separated by free-energy barriers, so that transitions
between states are rare events. Atomistic computer simulations of rare events
– especially molecular dynamics (MD) – can play an important role,
complementary to experiments, by predicting mechanisms, free-energy
landscapes, and kinetic rates.

However, MD simulations face two prominent challenges. First, a gap – often of
many orders of magnitude – between the short time scale of atomic motion
(femtoseconds) and the long time scale of rare events hampers direct, brute-
force simulations. Second, the high-dimensional nature of configuration space
makes the analysis of the atomistic trajectories and the extraction of
relevant information intrinsically difficult. [1, 2]

To overcome or at least alleviate these challenges and gain insights into the
transformation processes, including quantitative thermodynamic and kinetic
properties, it is necessary to reduce the dimensionality of the problem by
projecting on an appropriate low-dimensional space. To this aim, collective
variables (CVs) – generic functions of atomic coordinates able to track
interesting structural changes – are often introduced, heuristically or by
machine learning. [3]

Among all possible CVs describing a transition, the optimal CV, or “reaction
coordinate" (RC), is widely considered the so-called committor function: in MD
simulations, this function is defined as the probability $p_{B}(x)$ that a
system will reach state $B$ before state $A$ starting from atomic positions
$x$ in $3N$-dimensional configuration space, with initial momenta randomly
drawn from the equilibrium distribution [4, 5, 6, 7, 8, 9, 10].

The committor is considered optimal since i) it allows predicting the fate of
atomic configurations towards reactants or products, and ii) it preserves the
kinetics when employed to build a reduced model of the dynamics (vide infra).
The concept dates back to the 1930’s, when Onsager studied the recombination
of a pair of ions in the presence of a uniform electric field [11].

The committor can be used to design tests for the “quality" of a RC: for
instance, a good RC $q$ is expected to display a distribution of values
$p_{B}(x)$ sharply peaked at 0.5 when considering a set of different
configurations $x$ corresponding to a same location $q(x)=q*$ in RC space,
i.e., precisely identifying the transition state configurations. [4, 12, 13]

Some remarks are in order. Given a generic CV $q$, the corresponding free
energy landscape is always well-defined mathematically, via the equilibrium
density $F(q)=-kT\log\rho_{\mathrm{e}q}(q)$: if the CV is able to resolve (at
least partially) reactants and products, such $F(q)$ will display a barrier.
However, different CVs can contain different amounts of information about the
transition mechanism. Moreover, the height of the barrier will depend on the
particular choice of the CV [14, 15] (see Figure 1).

The barrier along the optimal RC is supposed to contain more useful
information from the viewpoint of the calculation of the rate, keeping in mind
that any invertible transformation of the optimal RC results again in an
optimal RC, since such transformation cannot produce or destroy microscopic
information: in the case of a non-linear transformation, the free-energy
barrier and the diffusion coefficient are different for the two RCs [16, 14],
the differences compensating each other in such a way that the projected
dynamics keeps the same kinetic properties.

It is therefore important to distinguish the problem of converging the
calculation of the free energy barrier along any given CV from a statistical
viewpoint (a sampling problem) from the problem of identifying among all CVs
the optimal RC and computing the corresponding barrier, for the purpose of
estimating the rate.

Several methods, based on different theoretical principles, have been put
forward to estimate optimal RCs from MD data. Given the customary
identification of the optimal RC with the committor function, a possibility
consists in optimizing a parametrization of the committor function based on
$p_{B}(x)$ values estimated from transition path sampling (TPS), exploiting,
e.g., the histogram test [13], maximum-likelihood approaches [17, 9], or
artificial intelligence [18, 19].

We remark, however, that precisely estimating committor values very close to
zero or to one by shooting requires an unfeasible large amount of
trajectories: this hampers the accurate reconstruction of the optimal RC in
presence of large free-energy barriers elsewhere than in the vicinity of the
barrier top. This problem is addressed in the reweighted path ensemble
approach [20].

Several approaches have been developed to generate flexible RC representations
able to capture slow degrees of freedom connected to rare events [21, 22, 23,
24, 25, 26, 27]. Starting from the seminal work of Ma and Dinner [18], where a
neural network was trained with committor data, recent years saw a flourishing
of machine-learning approaches to RC optimization. [28, 29, 30, 31, 32]

In this work we address with a single methodology several prominent tasks of
molecular simulations: identifying an optimal RC, estimating the free energy
landscape and the diffusion along such RC, and estimating the rate of the
process.

Results in the literature indicate that i) a model of the projected dynamics
along the optimal RC gives more direct access to accurate kinetic rates [13,
8, 9], and ii) the accurate rate is the minimal one – attained for the optimal
RC – with respect to the rates computed from the projected dynamics along each
of all possible CVs [10]. The method we present exploits these principles:
putative CVs are generated in a Monte Carlo scheme, and the kinetic rates of
associated effective low-dimensional models are estimated and systematically
minimized, leading to optimal RCs and accurate kinetic predictions. An
important ingredient of our approach is the automatic construction of
stochastic (Langevin) models of CV evolution: for this purpose, we adopt the
maximum likelihood algorithm proposed in Ref. 33, based on short TPS-like MD
trajectories. We remark that such task is non-trivial from a numerical
viewpoint [34], and alternative approaches exist [35, 36, 37, 38].

##  2 Theoretical methods

###  2.1 Variational principle: the optimal coordinate provides the minimal
rate

We start by considering ergodic diffusive dynamics (as described by overdamped
Langevin equations) in a high-dimensional space, for a system with two
metastable states $A$ and $B$. The transition rate from state $A$ to $B$ is
denoted by $k_{A\rightarrow B}$. Next, we consider the effective dynamics
resulting from the projection of the full high-dimensional dynamics onto a
low-dimensional manifold of collective variables (CVs): such effective
dynamics is an approximation of the original one, formulated again in terms of
diffusive equations [39, 10]. the reaction rate of the effective dynamics
between states $\tilde{A}$ and $\tilde{B}$, defined in the low-dimensional
space, is $\tilde{k}_{\tilde{A}\rightarrow\tilde{B}}$. Zhang, et al.[10]
proved that the transition rate of the full dynamics is always smaller or
equal to the one computed using effective dynamics. In other words, the
optimal reaction coordinate yields a minimal rate,

$k_{A\rightarrow B}\leq\tilde{k}_{\tilde{A}\rightarrow\tilde{B}}$ (1)

The optimal projection, preserving the value of the rate constant
$k_{A\rightarrow B}=\tilde{k}_{\tilde{A}\rightarrow\tilde{B}}$, is achieved
when the CV corresponds to the committor function or any invertible function
of the latter.

As a complementary result, Ref. [40] formally proved that it is possible to
obtain accurate rates from an overdamped Langevin model when projecting high-
dimensional underdamped Langevin dynamics on a single CV, the committor,
irrespective of the existence of timescale separation and metastability in the
system. Considering that typical MD simulations, where Hamilton’s equations
are coupled with a thermostat, are akin to high-dimensional underdamped
Langevin equations, the previous result should hold also in conventional MD.

Our approach consists in optimizing a RC defined as a function of several
trial CVs by minimizing the kinetic rate of the projected dynamics with
respect to all possible definitions of the RC. Here we consider the projection
on a one-dimensional CV, assuming that the variational principle connecting
the optimal CV to the minimal rate holds also when the input data is
represented by MD simulations with a thermostat. The projected dynamics is
approximated by an overdamped Langevin equation, consistently with Refs. [10,
40]:

$\dot{q}=-\beta D(q)\frac{\partial F(q)}{\partial q}+\frac{\partial
D(q)}{\partial q}+\sqrt{2D(q)}\,\eta(t)\leavevmode\nobreak\ ,$ (2)

with $D(q)$ the diffusion coefficient, $F(q)$ the free energy profile and
$\eta(t)$ a Gaussian white noise of zero mean and unit variance. The use of
this stochastic model has several appealing features: the only parameters of
the equation consist in the functions $F(q)$ and $D(q)$, the mass and
velocities are not explicitly needed, the Markovian character allows for
simple likelihood expressions, and the mean first passage time (MFPT) of the
$A\rightarrow B$ transition, inverse of the kinetic rate $k$, can be directly
obtained by numerical integration [41, 14]:

$k^{-1}=\int_{q_{0}}^{b}dx\,\frac{e^{\beta
F(x)}}{D(x)}\int_{a}^{x}dy\,e^{-\beta F(y)}$ (3)

where $q_{0}$ is the starting point in $A$ for the MFPT computation, while $a$
and $b$ are respectively the reflecting and absorbing boundaries.

###  2.2 Algorithm for reaction coordinate optimization

The phase-space information used for RC optimization is represented by MD
trajectories spanning transitions between the reactants and products regions
of configuration space. We adopt TPS as a source of input data for several
reasons: [12, 13, 9] i) ergodic trajectories (spontaneously spanning the
transitions) are unfeasible in presence of free-energy barriers $\gg k_{B}T$,
while TPS has an affordable cost for many systems; ii) TPS trajectories are
free from biasing forces and, at statistical convergence, faithfully reproduce
ergodic trajectories; iii) recent numerical evidence indicate that $\sim 100$
TPS-like MD trajectories projected on a suitable CV are sufficient to
reconstruct accurate free-energies and rates by means of overdamped Langevin
models [33].

The proposed RC optimization algorithm starts from the projection of TPS
trajectories on a basis set of $N$ potentially-relevant CVs ${\bf
Q}(t)=(Q_{1}(t),Q_{2}(t)...Q_{N}(t))$. All these coordinates are put on the
same footing by shifting and scaling so that for each of them the range of
variation on the actual MD trajectories is between 0 and 1 (oriented so that
$A\rightarrow B$ for growing values). Note that such a linear transformation
of the coordinate does not deform the corresponding free energy landscape
besides an irrelevant additive constant: passing from $x$ to $y=ax+b$ the
probability densities transform according to
$e^{-\beta[F(x)-\tilde{F}(y)]}=\rho(x)/\tilde{\rho}(y)=dy/dx=a$. Moreover, the
new diffusion coefficient is scaled, $\tilde{D}_{y}/D_{x}=(dy/dx)^{2}=a^{2}$
[42], as can be easily seen in the simple case of a driftless constant-$D$
diffusion process: $\langle[x(t)-x(0)]^{2}\rangle=2D_{x}t$ implies that
$\langle[y(t)-y(0)]^{2}\rangle=2(a^{2}D_{x})t=2\tilde{D}_{y}t$.

A trial RC $q=\sum_{1}^{N}w_{i}\,Q_{i}\equiv{\bf w}\cdot{\bf Q}$ as a linear
combination of the basis-set CVs is generated from random weights normalized
as ${\bf w}^{2}=1$. A Monte Carlo loop is then started, at each step proposing
a new RC $q_{\text{new}}$ obtained by modifying the weights of the last
accepted step through small random variations ${\bf w}\rightarrow{\bf
w}+\delta{\bf w}$ (each $\delta w_{i}$ being drawn from a uniform distribution
between $[-0.05,0.05]$).

A newly proposed RC is accepted or rejected based on a Metropolis criterion
aimed at minimizing the kinetic rate estimated from a maximum-likelihood
Langevin model. The latter is optimized following the method in Ref. [33]: for
a sufficiently small time interval $\tau$, the transition probability density
$p$ (propagator) between points $q$ and $q^{\prime}$ in CV-space can be
approximated as [43]

$p(q^{\prime},t+\tau|q,t)\approx\frac{1}{\sqrt{2\pi\mu}}e^{-(q^{\prime}-q-\phi)^{2}/2\mu}$
(4) $\phi=a\tau+\frac{1}{2}(aa^{\prime}+Da^{\prime\prime})\tau^{2}\ \ ,\
\mu=2D\tau+(aD^{\prime}+2a^{\prime}D+DD^{\prime\prime})\tau^{2}$ (5)

where the prime indicates $d/dq$, $a=-\beta DF^{\prime}+D^{\prime}$ is the
drift in Eq. 2, and the approximation includes terms up to order $\tau^{2}$.
We verified that resorting to the less accurate first-order propagator
$\phi=a\tau=(-\beta DF^{\prime}+D^{\prime})\tau$, $\mu=2D\tau$ gives
unsatisfactory results for the systems considered in this work.

Based on Eq. 5, the likelihood of the MD trajectory, sampled with a time
resolution $\tau$ and projected on a CV $q$ is given by

$-\log\mathcal{L}(\theta)=\sum_{k=1}^{M-1}\Bigl{\\{}\frac{1}{2}\log[2\pi\mu_{k}(\tau)]+\frac{[q_{k+1}-q_{k}-\phi_{k}(\tau)]^{2}}{2\mu_{k}(\tau)}\Bigr{\\}}$
(6)

where $\theta$ represents the set of all parameters of the Langevin equation
(i.e., the profiles $F(q)$ and $D(q)$), $M$ is the number of trajectory
configurations, and $\phi_{k}(\tau)\equiv\phi(q_{k},\tau)$,
$\mu_{k}(\tau)\equiv\mu(q_{k},\tau)$. The optimal Langevin model for the given
$q$ and $\tau$, obtained by minimizing $-\log\mathcal{L}(\theta)$ as a
function of the parameters with an iterative stochastic algorithm (see Ref.
[33] for details), directly yields the free energy and diffusion profiles.

The kinetic rate is computed from the Langevin model using Eq. 3, with
integral boundaries defined in the following way: $a$ is the smallest observed
value of $q$ (for each CV the transition $A\rightarrow B$ is in the direction
of increasing $q$ values), while $q_{0}$ and $b$ are the average final values
of $q$ for shooting trajectories ending in $A$ and $B$, respectively, i.e.,
the bottom of the corresponding free-energy minima.

Having obtained the rate corresponding to a given proposed CV, a Metropolis-
like test is employed to accept or reject the CV based on the following
expression for the probability:

$P=\min\left[1\,,\,\Big{(}\frac{k_{\text{old}}}{k_{\text{new}}}\Big{)}^{\alpha}\times\frac{\tau^{\text{noise}}_{\text{old}}}{\tau^{\text{noise}}_{\text{new}}}\right]$
(7)

this expression tends to minimize the rate (with $\alpha$ playing the role of
an adjustable inverse Monte Carlo temperature, as discussed in the Results)
and simultaneously – coherently with the hypothesis of a Markovian Langevin-
like behavior – the memory of the stochastic process. The latter is quantified
by $\tau^{\text{noise}}$, the auto-correlation time of the “observed" noise
trajectory $G_{k}$ estimated from the $q_{k}$ trajectory via the optimal
Langevin model:

$G_{k}=\frac{q_{k+1}-q_{k}+[\beta
D(q_{k})F^{\prime}(q_{k})-D^{\prime}(q_{k})]\tau}{\sqrt{2D(q_{k})\tau}}$ (8)

(see Ref. [33] for details).

To summarize, the algorithm performs the following steps:

  1. 1.

Project TPS trajectories relaxing from a barrier top onto a basis set of CVs
$\\{Q_{i}\\}_{i=1,...,N}$

  2. 2.

Construct a first trial RC $q$ as a sum of the basis CVs with random weights

  3. 3.

Generate a new CV $q_{\text{new}}$ by randomly modifying the weights

  4. 4.

Estimate the free energy and diffusion profiles as well as the $A\rightarrow
B$ rate from an optimal Langevin model of the observed trajectory
$q_{\text{new}}(t)$ using likelihood maximization

  5. 5.

Accept or reject the new CV with the probability in Eq. 7, aimed at minimizing
the rate and enforcing Markovianity

  6. 6.

Go to step 3 (iterate until convergence of the rate)

At convergence, the algorithm provides an optimal RC, as well as the
corresponding free-energy and diffusion profiles and kinetic rate, all
starting solely from a set of $\leavevmode\nobreak\ 100$ TPS-like short MD
trajectories. If the basis set contains all the relevant degrees of freedom
for the transition process, the optimal rate should approach the exact rate.

###  2.3 Analytic double-well potential

A two-dimensional double-well free-energy surface is defined, for benchmark
purposes, as the sum of two Gaussian-shaped wells:

$F(x,y)=-C\left[e^{-\frac{(x-x_{0})^{2}}{2\sigma_{x}^{2}}}e^{-\frac{(y-y_{0})^{2}}{2\sigma_{y}^{2}}}+e^{-\frac{(x-x_{1})^{2}}{2\sigma_{x}^{2}}}e^{-\frac{(y-y_{1})^{2}}{2\sigma_{y}^{2}}}\right]$
(9)

with $C=20\,k_{B}T$, Gaussian centers $(x_{0},y_{0})$=$(0.1723,0.5058)$,
$(x_{1},y_{1})$=$(0.8060,0.5058)$, and widths $\sigma_{x}^{2}=0.03921$,
$\sigma_{y}^{2}=0.2519$.

We consider the region $x\in[0,1]$ and $y\in[0,1]$, see Fig. 1a. The diffusion
matrix was set to $\mathbf{D}=0.015\cdot\big{(}\begin{smallmatrix}1&0\\\
0&1\end{smallmatrix}\big{)}$ ps-1.

The reference MFPT was estimated directly from brute-force trajectories. For
this purpose, 10 long overdamped Langevin simulations with a time step of 0.1
fs were performed. For the aggregate simulation time (10 $\mathrm{\mu}$s), 16
039 jumps from state $A$ to $B$ were observed, yielding a MFPT of $5600\pm 40$
ps and a transition rate $k_{A\rightarrow B}=1.79\pm 0.01\cdot 10^{-4}$ ps-1.

Input trajectories for RC optimization were obtained by shooting 100
overdamped Langevin trajectories from the barrier top. From these
trajectories, 51 relaxed to state $A$ (left well) and 49 relaxed to state $B$
(right well). The basis set CVs, in this case, are simply $x$ and $y$,
yielding two one-dimensional projections $x(t)$ and $y(t)$ of the 100 relaxing
trajectories.

###  2.4 Fullerene dimers in water

The association and dissociation of C60 and C240 fullerene dimers (OPLS-AA
force-field [44]) in water solution (SPC force-field [45]) have been simulated
by MD using GROMACS v2019.4 [46, 47] patched with PLUMED 2.5.3 [48]. The
simulations are similar to those in Ref. 33: for more computational details
than what is summarized here, we refer to the latter article.

In the first system, two C60 molecules were solvated by 2398 water molecules
in a simulation box of $3.607^{3}$ nm3 with periodic boundary conditions. In
the second system, two C240 molecules were solvated by 5375 water molecules in
a simulation box of $5.22^{3}$ nm3. MD simulations were performed with a time
step of 1 fs in the $NPT$ ensemble at 298 K [49] (thermostat time constant = 1
ps) and 1 atm [50] (barostat time constant = 4 ps).

The reference free-energy profiles of the association/dissociation of the C60
fullerenes in water as a function of the 8 basis CVs were computed from
unbiased simulations of 1 $\mathrm{\mu}$s total aggregated time, from the
population histogram: $F(Q_{i})=-k_{B}T\log\rho_{\mathrm{eq}}(Q_{i})$. Error
bars were estimated as the standard deviation of the mean over 5 independent
replicas.

We used as input data for RC optimization a set of 100 aimless shooting [51,
17] trajectories of 20 ps, generated using the script publicly available at
https://github.com/physix-repo/aimless-shooting, employing a separation of 0.1
ps between successive shooting points.

For the C60 dimer we define the dissociated state based on the distance
between centers of mass as $d\geq 1.34$ nm and the associated state as $d\leq
1.17$ nm; for the C240 dimer we define the dissociated and associated states
as $d\geq 2.01$ nm and $d\leq 1.9$ nm, respectively. The reference rate
constants for the dissociation of the C60 and C240 dimers from the associated
state were estimated using the reactive flux formalism over 1000 aimless
shooting trajectories (see Ref. 33 for details), obtaining for C60 a MFPT of
$6.1\pm 1.2$ ns and for C240 a MFPT of $9.4\pm 1.1$ $\mathrm{\mu}$s.

The basis set CVs employed for projecting the fullerene dimers trajectories
are the following:

  1. 1.

$d$: the distance between the centers of mass of each fullerene molecule.

  2. 2.

$cc$: the number of carbon-carbon contacts, estimated summing continuous
coordination functions between any atom of the first fullerene (set $S_{1}$)
and any atom of the second fullerene (set $S_{2}$)

$cc=\sum_{i\in S_{1}}\sum_{j\in S_{2}}C_{ij}\ ,\ \ \
C_{ij}=\frac{1-\Big{(}\frac{r_{ij}}{r_{0}}\Big{)}^{n}}{1-\Big{(}\frac{r_{ij}}{r_{0}}\Big{)}^{m}}$
(10)

where $r_{ij}$ is the distance between atoms $i$ and $j$, with parameters
$r_{0}=0.35$ nm, $n=6$ and $m=10$.

  3. 3.

$c2w$: the number of carbon-water contacts, defined similarly to $cc$ in Eq.
10, in this case with $r_{0}=0.6$ nm, set $S_{1}$ including all carbon atoms
of the two fullerene molecules, and set $S_{2}$ including all oxygen atoms of
the water molecules.

  4. 4.

$c1w$: the water-carbon contacts for a single fullerene molecule, defined as
$c2w$ except for the inclusion of a single fullerene in set $S_{1}$.

  5. 5.

$sc$: the approximate carbon pair entropy, estimated using the expression [52]

$sc=-2\pi\rho k_{B}\int_{0}^{r_{\mathrm{max}}}\left[g(r)\ln
g(r)-g(r)+1\right]r^{2}\mathrm{d}r\leavevmode\nobreak\ ,$ (11)

where $\rho$ is the density, $r_{\mathrm{max}}$ is a cutoff set to 0.65 nm,
and $g(r)$ is the pair distribution function of carbon atoms, estimated via
Gaussian kernels as

$g(r)=\frac{1}{4\pi N\rho r^{2}}\sum_{i\neq
j}\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-(r-r_{ij})^{2}/(2\sigma^{2})}\leavevmode\nobreak\
,$ (12)

where $N$ is the number of carbon atoms and $\sigma=0.01$ nm.

  6. 6.

$sw$: the approximate water pair entropy, estimated with the same equations
and parameters as $sc$ applied to oxygen atoms.

  7. 7.

$ucc$: the Van der Waals carbon-carbon potential energy, computed over all
carbon pairs.

  8. 8.

$ucw$: the Van der Waals carbon-water potential energy, computed between all
the carbon atoms and all the solvent atoms.

In the Results section, for each putative RC in the optimization algorithm,
the likelihood of a Langevin model was maximized employing $3\cdot 10^{6}$ and
$2\cdot 10^{6}$ iteration steps for the $C_{60}$ and $C_{240}$ fullerenes,
respectively, adopting a time resolution $\tau=1$ ps for the projected MD
trajectory. The latter was tested as sufficient to yield time-decorrelated
noise when using the dimer center-of-mass distance $d$ as CV.

##  3 Results and discussion

###  3.1 Two-dimensional double well

To test the validity of the algorithm developed, we first benchmark it on an
analytic model: the double-well system detailed in the Methods section. By
construction the dynamics is of overdamped-Langevin nature, the trajectories
being obtained by integration of such stochastic differential equation over
the two-dimensional $F(x,y)$ landscape in Figure 1 with a constant isotropic
diffusion matrix.

After projecting the trajectories on a one-dimensional CV, several non-trivial
effects can be anticipated from theory. Figure 1 illustrates for instance
three simple choices of the CV, corresponding to the $x$ axis, the $y$ axis,
or the 45∘ diagonal.

Langevin trajectories shown in Figure 1(b) clearly resolve two distinct end-
states when adopting $q\equiv x$, while the separation becomes less clear upon
deviation of $q$ from the $x$ axis until vanishing for $q\equiv y$.

The exact one-dimensional $F(q)$ profile obtained by projecting on a generic
CV $q(x,y)$ can be computed from probability marginalization as

$e^{-\beta F(q)}=\int dx\int dy\,e^{-\beta F(x,y)}\,\delta[q-q(x,y)]$ (13)

Figure 1(c) shows that CV changes result in significant differences in the
one-dimensional free-energy profiles. In particular the barrier is maximal for
$q\equiv x$, it is strongly reduced for the diagonal, and it vanishes for
$q\equiv y$; analogously, the distance between the two minima is largest for
$q\equiv x$ and it vanishes for $q\equiv y$. These observations support the
intuitive idea that $x$ is the optimal RC for this problem, fully capturing
the progress of the transition mechanism. Moreover, the reduced barriers for
alternative CVs – simply the effect of overlapping contributions of the two
wells – suggest faster rates for the corresponding effective dynamics. The
effective one-dimensional diffusion coefficient $D(q)$ (Figure 1(c)) is also
affected by projection, increasing when $q$ deviates from $x$, however its
effect on the rate is linear, less important than barrier variations, which
have an exponential effect.

Figure 1: (a) Analytic double-well free-energy landscape, with arrows
indicating one-dimensional projections along the $x$ axis, the diagonal, and
the $y$ axis, respectively. (b) 100 Langevin bi-dimensional trajectories,
projected on a single CV, relax from the barrier-top into $A$ (blue) or $B$
(orange), and are used to optimize effective one-dimensional Langevin models
via likelihood maximization. (c) For each projection, the free-energy and
diffusion profiles inferred from the latter models are shown (continuous lines
with standard error over 10 reconstructions) alongside the exact $F(q)$
profiles and an alternative estimate of $D(q)$ (dashed lines), from the
Kramers-Moyal formula [42] $D(q)\approx(\langle\Delta
q^{2}\rangle_{q}-\langle\Delta q\rangle_{q}^{2})/2\tau$ where $\Delta q$ is
the CV displacement in a time interval $\tau$, with binned averages along $q$.

Starting from $x$ and $y$ as basis CVs, 10 independent runs of the RC
optimization algorithm are performed, aimed at minimizing the rate of the
effective projected model (at time resolution $\tau=0.01$ ps), each starting
with a different random linear combination as initial CV (see Figure 2).

The parameter $\alpha$ (effective Monte Carlo inverse temperature in Eq. 7) is
varied from 2 to 4: in all tested cases the initial rate decreases almost
monotonously, often by several orders of magnitude, until reaching an
equilibrium behavior with small fluctuations in a small interval. Low $\alpha$
values increase the acceptance of putative CVs in presence of an increase of
the rate, whereas high $\alpha$ values decrease such acceptance causing the
stochastic optimization to converge more rapidly towards a minimal rate, with
smaller equilibrium fluctuations. $\alpha$ plays therefore the role of an
adjustable parameter providing some control over the convergence of the rate
minimization algorithm.

Irrespective of the initial putative CV the algorithm successfully and
consistently retrieves an optimal combination including at least $95\%$ of
$x$. Given that after the initial relaxation a stationary distribution of
equilibrated $k$ values is reached, with a finite width that depends on the
parameter $\alpha$, introducing a degree of fuzziness, different criteria
could be envisaged to identify the optimal rate values and the corresponding
optimal RCs. In the case of the double well system, identifying the predicted
rate with the average value of the stationary distribution and its uncertainty
with the standard deviation yields $2.6\pm 1.0\cdot 10^{-4}$ ps-1 for
$\alpha=3$, in good agreement with the reference brute-force rate of $1.8\cdot
10^{-4}$ ps-1.

However, the results for the solvated fullerene dimers (see next section), a
realistic complex system, indicate the possibility of a few spurious outliers
at the lowest $k$ values, resulting from imperfect optimization of the
Langevin model for a few putative CVs. This suggests an alternative criterion
to identify the optimal rate predicted by the algorithm, as the 5th percentile
in the stationary distribution, i.e., the rate below which 5% of the rates are
found in the distribution. This criterion predicts (for $\alpha=3$) a 5th
percentile value of $1.6\cdot 10^{-4}$ ps-1: taking as optimal prediction the
average over the 5 rates closest to such value for each of the 10 independent
optimization runs, we get a predicted rate of $1.60\pm 0.04\cdot 10^{-4}$ ps-1
again in good agreement with the reference.

Figure 2: RC optimization by rate minimization applied to the double-well
benchmark: (a) the rate evolution is shown for 5 independent optimizations
with $\alpha=2,3,4$, respectively; (b) the histogram of the rate values in the
stationary regime, i.e., discarding the initial relaxation; (c) number of
correlated time steps (at O.O1 ps resolution) of the effective noise estimated
from the original trajectories using the optimal Langevin model for each CV.

We remark that, even though the original dynamics in $(x,y)$-space is
Markovian, non-Markovian effects could appear on the effective dynamics when
projecting on CVs different from the committor, requiring in principle an
underdamped (or generalized) Langevin description [40]. Therefore, the term
penalizing deviations from Markovianity in Eq. 7 plays, in general, an
important role, helping to reject proposed CVs that would introduce
significant memory effects rendering inappropriate the overdamped model.
Moreover, as shown in Ref. 33, from the viewpoint of the inference of Langevin
models non-Markovian behavior appears to be correlated with barrier
overestimation and thus rate underestimation, a spurious effect that must be
avoided to exploit correctly the variational principle connecting the minimal
rate to the optimal RC.

Our simulations show that for all values of $\alpha$ the vast majority of the
CVs accepted during the optimization process have an effective noise free from
time-correlation (Figure 2), with occasional exceptions, confirming the
ability of the algorithm to enforce Markovianity.

###  3.2 Interaction between fullerene nanoparticles in water

We applied the RC optimization algorithm to the case of the dissociation
process of fullerene dimers in water solution, respectively C60 and C240.
These systems are rather complex, including thousands of atoms, and feature an
associated state, corresponding to a free-energy minimum, characterized by the
carbon nanoparticles in close contact without the mediation of water
molecules, and a dissociated state where each fullerene is fully-solvated. The
transition state region features a water molecule bridging between the
fullerenes (Figure 3), corresponding to the top of a sizable free-energy
barrier when a CV able to resolve the metastable states, like the distance $d$
between centers of mass, is employed [33].

Figure 3: Characteristic atomic configurations for the C60 and C240 dimers in
bulk water solution. Water molecules distant $<4$ Å from carbon atoms are
shown. The progress from the associated state (left), passing through the
transition state (center), until the dissociated state (right), is followed by
CVs normalized to span between zero and one. A water molecule in the
transition state, bridging between fullerenes and breaking the direct
hydrophobic interactions between them, is highlighted.

We start by proposing a pool of eight basis CVs to be used as building blocks
for RC optimization: the set ranges from geometry-inspired CVs like the
distance $d$ and the total number of carbon-carbon or carbon-water contacts
($cc$, $c2w$, $c1w$) to physics-inspired CVs like the approximate two-body
carbon or water entropy ($sc$, $sw$) and the carbon-carbon or carbon-water
interaction energy ($ucc$, $ucw$). It is not obvious a priori which CVs or
combination thereof could be the best approximations of an optimal RC.

The TPS trajectories relaxing from the barrier top, together with the free-
energy and diffusion profiles inferred from likelihood maximization of
Langevin models for each of these CVs are reported for the C60 dimer in Figure
4, and display, as expected, strong variations depending on the CV choice (see
the Supporting Information for C240). Clearly, CVs able to resolve well the
associated and dissociated states feature a barrier, contrary to CVs lacking
such resolving power. Overall, the $F(q)$ landscapes compare well with the
reference brute-force ones. A time resolution $\tau=1$ ps is adopted for
Langevin modelling, as it yielded Markovian behavior in combination with the
CV $d$ for both fullerene sizes in Ref. 33.

The basis CVs that results in the lowest kinetic rates, of the order of
$10^{-3}$ ps-1, are $d$, $cc$ and $ucc$; the CVs providing the highest (hence,
less accurate) rates are instead $sw$ and $c1w$, of the order of $10^{-1}$
ps-1. This allows a first ranking of the quality of the basis CVs. The
position-dependent diffusion coefficient $D(q)$ is roughly between $0.01$ and
$0.03$ ps-1 for all CVs except $sw$, that reaches almost $0.08$ ps-1 (Figure
4).

Figure 4:  MD trajectories from TPS of the C60 fullerene dimer in water
relaxing to the associated (blue) or dissociated (orange) state (left), with
the corresponding $F(q)$, $D(q)$ profiles and kinetic rates (right), for each
CV of the basis set ((a) to (h), definitions in the Methods section). The
profiles $F(q)$ and $D(q)$ with continuous lines are inferred from likelihood
maximization (average and standard error over 10 independent runs). The
kinetic rate $k$ is calculated using Eq. 3. The dashed $F(q)$ are obtained
from equilibrium MD histograms.

We applied the RC optimization algorithm (500 accepted steps, with $10^{6}$
Langevin optimization iterations for each putative CV) starting from different
initial random combinations of the basis CVs and testing the parameter
$\alpha=2,3,4$. As shown in Figure 5, for both fullerene sizes the rate
relaxes to a stationary distribution after about 200 steps: rate fluctuations
are wider in these systems than for the double well, however the lower part of
the distribution has the same order of magnitude of the reference rates.

Figure 5: RC optimization results for the C60 and C240 fullerene dimers in
water. (a,e) 5 independent rate minimization runs are displayed, the grey bar
corresponding to the reference rate (thickness = error bar) estimated via
reactive flux. (b,f) Histogram of the rate values from 10 optimizations in the
stationary regime (i.e., discarding the first 300 iterations) for different
$\alpha$ values. (c,g) Number of correlated time steps (at 1 ps resolution) of
the effective noise. (d,h) Weights distribution for the optimal RCs from 10
independent runs: for each run, the 5 RCs closest to the 5th percentile rate
are included.

Employing the 5th percentile of the stationary $k$ distribution as criterion
(see previous section) at $\alpha=3$, we obtain as optimal rates the values
$1.26\pm\leavevmode\nobreak\ 0.15\cdot 10^{-4}$ ps-1 and
$1.8\pm\leavevmode\nobreak\ 1.1\cdot 10^{-8}$ ps-1 for the dissociation C60
and C240 dimers, respectively, that compare well with the references reactive-
flux rates $1.6\pm\leavevmode\nobreak\ 0.3\cdot 10^{-4}$ ps-1 and
$1.1\pm\leavevmode\nobreak\ 0.1\cdot 10^{-7}$ ps-1, respectively (see Methods
section).

The Markovianity appears enforced in an effective way, with the vast majority
of the accepted CVs having uncorrelated effective noise and a few ones
displaying $<10$ correlated steps (see Figure 5(c, g)).

Finally, the composition of the optimal RCs for the solvated C60 and C240
dimers is reported in Figure 5(d, h): since the basis CVs are not uncorrelated
between them, the optimal RC is not uniquely defined and different weights in
the combination $q=\sum_{1}^{N}w_{i}\,Q_{i}$ can lead to RCs with similar
quality and similar rates.

In fact, performing a principal component analysis on the members of the set
of optimal RCs (i.e., the 50 RCs with rate closest to the 5th percentile in
Figure 5(b,f)) provides explained variance ratios for the first four principal
components $[0.40,0.24,0.14,0.11]$ for C60 and $[0.56,0.17,0.16,0.04]$ for
C240 (see Supporting Information). This shows that it is not possible to
easily reduce the dimensionality in the optimal space of RCs by means of one
or two dominating components.

Nevertheless, some basis CVs contribute more than others, on average, to the
optimal RCs: in particular, Figure 5(d,h) suggest that $d$, $cc$, $sc$, and
$ucc$ play a prominent role for both fullerene sizes. The latter CVs do not
contain water degrees of freedom, pointing, somehow counter-intuitively, to a
limited role of the solvent in the optimal RC for these systems.

##  4 Concluding remarks

We presented an algorithm for the automatic optimization of the RC for
activated processes starting from $\approx 100$ short TPS MD trajectories and
a basis set of putative CVs. The approach is based on the variational
principle in Ref. 10: the rate predicted by a Langevin model of the effective
projected dynamics is minimized by the optimal RC, coinciding in this limit
with the true MD rate.

As an essential ingredient, we adopt the method in Ref. 33 to infer an optimal
overdamped Langevin model for each putative RC, thus reliably recovering the
corresponding free-energy and diffusion profiles as well as the kinetic rate.

Since we make the assumption that an overdamped stochastic equation be able to
faithfully model the projected dynamics, we explicitly disfavoured the
appearance of memory effects in the acceptance probability expression applied
to newly proposed CVs, thus keeping the Markovianity under control.

Numerical tests on an analytic double-well system as well as on MD simulations
of carbon nanoparticles in explicit water solution indicate that optimal RCs
are systematically identified in a robust and efficient way: this allowed us
to predict kinetic rates in the microsecond time scale from $\sim 10$ ps-long
trajectories.

A foremost advantage of the present algorithm is that RC optimization is done
as post-processing of a MD data set: a large number of putative CVs can be
screened with limited computer resources without the need to re-run expensive
MD calculations. Moreover, compared to alternative approaches based on machine
learning the committor, there is no need to estimate committor values over a
large sample of atomic configurations, a computationally-expensive endeavour
that, by construction, is limited to the region close to the barrier top. On
the contrary, the proposed approach only requires a small number of reactive
trajectories and, by construction, it learns the optimal RC across the full
span of the transition paths joining metastable states, even where the
committor assume values very close to zero or to one.

Here the basis set of CVs is built from heuristic considerations, while in
future works it could be advantageous to devise an unsupervised building
procedure. Likewise, from the simplistic formulation of optimal RCs as linear
combinations of basis CVs, it would be beneficial for complex processes to
adopt more flexible, nonlinear formulations, possibly including neural
networks. As a final word of caution, the effective description adopted here
in terms of overdamped Langevin equations is appropriate for a subset of all
possible rare-events systems and processes: underdamped or generalized (non-
Markovian) stochastic models can be necessary in other cases.

We expect that the approach introduced here could help discovering optimal RCs
and predicting accurate rates for challenging open problems like crystal
nucleation and protein-protein or protein-ligand dissociation, whose
complexity, so far, eluded systematic and predictive computational studies.

##  5 Acknowledgement

We gratefully acknowledge very insightful discussions with Christoph Dellago,
Gerhard Hummer, Hadrien Vroylandt, Arthur France-Lanord, Alessandro Barducci,
Bettina Keller, Tony Lelièvre and Edina Rosta. This work was performed with
the support of the Institut des Sciences du Calcul et des Données (ISCD) of
Sorbonne University (IDEX SUPER 11-IDEX-0004). Calculations were performed on
the GENCI-IDRIS French national supercomputing facility, under grant numbers
A0110811069, A0120901387, A0130811069.

##  6 Supporting information

Figure 6: Optimal RC composition for the benchmark double-well system with
$\alpha=3$ Figure 7:  MD trajectories from TPS of the C240 fullerene dimer in
water relaxing to the associated (blue) or dissociated (orange) state (left),
with the corresponding $F(q)$, $D(q)$ profiles and kinetic rates (right), for
each CV of the basis set ((a) to (h), definitions in the Methods section). The
profiles $F(q)$ and $D(q)$ with continuous lines are inferred from likelihood
maximization (average and standard error over 10 independent runs). The
kinetic rate is calculated using Eq. 3. Figure 8: Time auto-correlation of the
effective noise of the CV basis set for C60, evaluated on 100 TPS MD
trajectories. The solid lines are the average correlation of 10 independent
optimal models (see main text) while the shaded line represents the standard
deviation. Figure 9: Time auto-correlation of the effective noise of the CV
basis set for C240, evaluated on 100 TPS MD trajectories. The solid lines are
the average correlation of 10 independent optimal models (see main text) while
the shaded line represents the standard deviation. Figure 10: Principal
component analysis on the 50 members of the optimal set of RCs (a) for C60 and
(b) for C240.

## References

  * Pietrucci [2017] Fabio Pietrucci.  Strategies for the exploration of free energy landscapes: Unity in diversity and challenges ahead.  _Rev. Phys._ , 2:32–45, 2017.
  * Glielmo et al. [2021] Aldo Glielmo, Brooke E Husic, Alex Rodriguez, Cecilia Clementi, Frank Noé, and Alessandro Laio.  Unsupervised learning methods for molecular simulation data.  _Chem. Rev._ , 121(16):9722–9758, 2021.
  * Gkeka et al. [2020] Paraskevi Gkeka, Gabriel Stoltz, Amir Barati Farimani, Zineb Belkacemi, Michele Ceriotti, John D Chodera, Aaron R Dinner, Andrew L Ferguson, Jean-Bernard Maillet, Hervé Minoux, et al.  Machine learning force fields and coarse-grained variables in molecular dynamics: application to materials and biological systems.  _J. Chem. Theory Comput._ , 16(8):4757–4775, 2020\.
  * Bolhuis et al. [2002] Peter G Bolhuis, David Chandler, Christoph Dellago, and Phillip L Geissler.  Transition path sampling.  _Annu. Rev. Phys. Chem_ , 53:291–318, 2002.
  * Rhee and Pande [2005] Young Min Rhee and Vijay S Pande.  One-dimensional reaction coordinate and the corresponding potential of mean force from commitment probability distribution.  _J. Phys. Chem. B_ , 109(14):6780–6786, 2005\.
  * Metzner et al. [2006] Philipp Metzner, Christof Schütte, and Eric Vanden-Eijnden.  Illustration of transition path theory on a collection of simple examples.  _J. Chem. Phys_ , 125(8):084110, 2006.  doi: 10.1063/1.2335447.
  * E and Vanden-Eijnden [2010] Weinan E and Eric Vanden-Eijnden.  Transition-path theory and path-finding algorithms for the study of rare events.  _Annu. Rev. Phys. Chem._ , 61:391–420, 2010.
  * Banushkina and Krivov [2016] Polina V Banushkina and Sergei V Krivov.  Optimal reaction coordinates.  _Wiley Interdiscip. Rev. Comput. Mol. Sci._ , 6(6):748–763, 2016.
  * Peters [2016] Baron Peters.  Reaction coordinates and mechanistic hypothesis tests.  _Annu. Rev. Phys. Chem._ , 67:669–690, 2016.
  * Zhang et al. [2016] Wei Zhang, Carsten Hartmann, and Christof Schütte.  Effective dynamics along given reaction coordinates, and reaction rate theory.  _Faraday Discuss._ , 195:365–394, 2016.
  * Onsager [1938] L Onsager.  Initial recombination of ions.  _Phys. Rev._ , 54:554–557, 1938.
  * Hummer [2004] Gerhard Hummer.  From transition paths to transition states and rate coefficients.  _J. Chem. Phys_ , 120(2):516–523, 2004.
  * Best and Hummer [2005] Robert B Best and Gerhard Hummer.  Reaction coordinates and rates from transition paths.  _Proc. Natl. Acad. Sci_ , 102(19):6732–6737, 2005\.
  * Jungblut and Dellago [2016] Swetlana Jungblut and Christoph Dellago.  Pathways to self-organization: Crystallization via nucleation and growth.  _Eur. Phys. J. E_ , 39(8):1–38, 2016.
  * Dietschreit et al. [2022] Johannes CB Dietschreit, Dennis J Diestler, Andreas Hulm, Christian Ochsenfeld, and Rafael Gómez-Bombarelli.  From free-energy profiles to activation free energies.  _J. Chem. Phys_ , 157(8):084113, 2022.
  * Krivov and Karplus [2008] Sergei V Krivov and Martin Karplus.  Diffusive reaction dynamics on invariant free energy profiles.  _Proc. Natl. Acad. Sci_ , 105(37):13841–13846, 2008.
  * Peters et al. [2007] Baron Peters, Gregg T Beckham, and Bernhardt L Trout.  Extensions to the likelihood maximization approach for finding reaction coordinates.  _J. Chem. Phys_ , 127(3):034109, 2007.
  * Ma and Dinner [2005] Ao Ma and Aaron R Dinner.  Automatic method for identifying reaction coordinates in complex systems.  _J. Phys. Chem. B_ , 109(14):6769–6779, 2005\.
  * Jung et al. [2019] Hendrik Jung, Roberto Covino, and Gerhard Hummer.  Artificial intelligence assists discovery of reaction coordinates and mechanisms from molecular dynamics simulations.  _arXiv preprint arXiv:1901.04595_ , 2019.
  * Rogal et al. [2010] Jutta Rogal, Wolfgang Lechner, Jarek Juraszek, Bernd Ensing, and Peter G Bolhuis.  The reweighted path ensemble.  _J. Chem. Phys_ , 133(17):174109, 2010.
  * Pérez-Hernández et al. [2013] Guillermo Pérez-Hernández, Fabian Paul, Toni Giorgino, Gianni De Fabritiis, and Frank Noé.  Identification of slow molecular order parameters for markov model construction.  _J. Chem. Phys_ , 139(1):07B604_1, 2013.
  * Tiwary and Berne [2016] Pratyush Tiwary and BJ Berne.  Spectral gap optimization of order parameters for sampling complex molecular systems.  _Proc. Natl. Acad. Sci_ , 113(11):2839–2844, 2016\.
  * Noé and Clementi [2017] Frank Noé and Cecilia Clementi.  Collective variables for the study of long-time kinetics from molecular trajectories: theory and methods.  _Curr. Opin. Struct. Biol_ , 43:141–147, 2017.
  * M. Sultan and Pande [2017] Mohammad M. Sultan and Vijay S Pande.  tica-metadynamics: accelerating metadynamics by using kinetically selected collective variables.  _J. Chem. Theory Comput._ , 13(6):2440–2447, 2017\.
  * Pérez de Alba Ortiz et al. [2018] A Pérez de Alba Ortiz, A Tiwari, RC Puthenkalathil, and B Ensing.  Advances in enhanced sampling along adaptive paths of collective variables.  _J. Chem. Phys_ , 149(7):072320, 2018.
  * Wu and Noé [2020] Hao Wu and Frank Noé.  Variational approach for learning markov processes from time series data.  _J Nonlinear Sci_ , 30(1):23–66, 2020.
  * He et al. [2022] Ziwei He, Christophe Chipot, and Benoît Roux.  Committor-consistent variational string method.  _J. Phys. Chem. Lett._ , 13(40):9263–9271, 2022\.
  * Ribeiro et al. [2018] João Marcelo Lamim Ribeiro, Pablo Bravo, Yihang Wang, and Pratyush Tiwary.  Reweighted autoencoded variational bayes for enhanced sampling (rave).  _J. Chem. Phys_ , 149(7):072301, 2018.
  * Chen and Ferguson [2018] Wei Chen and Andrew L Ferguson.  Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration.  _J. Comput. Chem._ , 39(25):2079–2102, 2018.
  * Sidky et al. [2020] Hythem Sidky, Wei Chen, and Andrew L Ferguson.  Machine learning for collective variable discovery and enhanced sampling in biomolecular simulation.  _Mol. Phys._ , 118(5):e1737742, 2020.
  * Hooft et al. [2021] Ferry Hooft, Alberto Pérez de Alba Ortíz, and Bernd Ensing.  Discovering collective variables of molecular transitions via genetic algorithms and neural networks.  _J. Chem. Theory Comput._ , 17(4):2294–2306, 2021\.
  * Sun et al. [2022] Lixin Sun, Jonathan Vandermause, Simon Batzner, Yu Xie, David Clark, Wei Chen, and Boris Kozinsky.  Multitask machine learning of collective variables for enhanced sampling of rare events.  _J. Chem. Theory Comput._ , 18(4):2341–2353, 2022\.
  * Palacio-Rodriguez and Pietrucci [2022] Karen Palacio-Rodriguez and Fabio Pietrucci.  Free energy landscapes, diffusion coefficients, and kinetic rates from transition paths.  _J. Chem. Theory Comput._ , 18(8):4639–4648, 2022\.
  * Ferretti et al. [2020] Federica Ferretti, Victor Chardes, Thierry Mora, Aleksandra M Walczak, and Irene Giardina.  Building general langevin models from discrete datasets.  _Phys. Rev. X_ , 10(3):031018, 2020.
  * Hummer [2005] Gerhard Hummer.  Position-dependent diffusion coefficients and free energies from bayesian analysis of equilibrium and replica molecular dynamics simulations.  _New J. Phys._ , 7(1):34, 2005.
  * Zhang et al. [2011] Qi Zhang, Jasna Brujić, and Eric Vanden-Eijnden.  Reconstructing free energy profiles from nonequilibrium relaxation trajectories.  _J. Stat. Phys._ , 144(2):344–366, 2011.
  * Boninsegna et al. [2018] Lorenzo Boninsegna, Feliks Nüske, and Cecilia Clementi.  Sparse learning of stochastic dynamical equations.  _J. Chem. Phys_ , 148(24):241723, 2018.
  * Sicard et al. [2021] Francois Sicard, Vladimir Koskin, Alessia Annibale, and Edina Rosta.  Position-dependent diffusion from biased simulations and markov state model analysis.  _J. Chem. Theory Comput._ , 17(4):2022–2033, 2021\.
  * Legoll and Lelievre [2010] Frédéric Legoll and Tony Lelievre.  Effective dynamics using conditional expectations.  _Nonlinearity_ , 23(9):2131, 2010.
  * Lu and Vanden-Eijnden [2014] Jianfeng Lu and Eric Vanden-Eijnden.  Exact dynamical coarse-graining without time-scale separation.  _J. Chem. Phys_ , 141(4):07B6191, 2014.
  * Zwanzig [2001] Robert Zwanzig.  _Nonequilibrium statistical mechanics_.  Oxford university press, 2001.
  * Risken [1996] Hannes Risken.  _Fokker-planck equation_.  Springer, 1996.
  * Drozdov [1997] Alexander N Drozdov.  High-accuracy discrete path integral solutions for stochastic processes with noninvertible diffusion matrices.  _Phys. Rev. E_ , 55(3):2496, 1997.
  * Jorgensen et al. [1996] William L Jorgensen, David S Maxwell, and Julian Tirado-Rives.  Development and testing of the opls all-atom force field on conformational energetics and properties of organic liquids.  _J. Am. Chem. Soc._ , 118(45):11225–11236, 1996\.
  * Berendsen et al. [1981] Herman JC Berendsen, James PM Postma, Wilfred F van Gunsteren, and Jan Hermans.  Interaction models for water in relation to protein hydration.  In _Intermolecular forces: proceedings of the fourteenth Jerusalem symposium on quantum chemistry and biochemistry held in jerusalem, israel, april 13–16, 1981_ , pages 331–342. Springer, 1981.
  * Berendsen et al. [1995] Herman JC Berendsen, David van der Spoel, and Rudi van Drunen.  Gromacs: a message-passing parallel molecular dynamics implementation.  _Comput. Phys. Commun._ , 91(1-3):43–56, 1995\.
  * Abraham et al. [2015] Mark James Abraham, Teemu Murtola, Roland Schulz, Szilárd Páll, Jeremy C Smith, Berk Hess, and Erik Lindahl.  Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers.  _SoftwareX_ , 1:19–25, 2015.
  * Tribello et al. [2014] Gareth A Tribello, Massimiliano Bonomi, Davide Branduardi, Carlo Camilloni, and Giovanni Bussi.  Plumed 2: New feathers for an old bird.  _Comput. Phys. Commun._ , 185(2):604–613, 2014\.
  * Bussi et al. [2007] Giovanni Bussi, Davide Donadio, and Michele Parrinello.  Canonical sampling through velocity rescaling.  _J. Chem. Phys._ , 126(1):014101, 2007.  doi: 10.1063/1.2408420.  URL https://doi.org/10.1063/1.2408420.
  * Parrinello and Rahman [1981] Michele Parrinello and Aneesur Rahman.  Polymorphic transitions in single crystals: A new molecular dynamics method.  _J. Appl. Phys._ , 52(12):7182–7190, 1981.
  * Peters and Trout [2006] Baron Peters and Bernhardt L. Trout.  Obtaining reaction coordinates by likelihood maximization.  _J. Chem. Phys._ , 125(5):054108, 2006.  doi: http://dx.doi.org/10.1063/1.2234477.  URL http://scitation.aip.org/content/aip/journal/jcp/125/5/10.1063/1.2234477.
  * Piaggi et al. [2017] Pablo M Piaggi, Omar Valsson, and Michele Parrinello.  Enhancing entropy and enthalpy fluctuations to drive crystallization in atomistic simulations.  _Phys. Rev. Lett._ , 119(1):015701, 2017.


Gleb Dubosarskii et al

Gleb Dubosarskii, Department of Electrical and Computer Engineering, Western
University, Ontario, Canada.

# Stochastic Geometry of Network
of Randomly Distributed Moving Vehicles
on a Highway

Gleb Dubosarskii  Serguei Primak  Xianbin Wang<EMAIL_ADDRESS>Dubosarskii G
S. Primak  X. Wang

(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Summary]Vehicular ad-hoc networks (VANETs) have become an extensively studied
topic in contemporary research. One of the fundamental problems that has
arisen in such research is understanding the network statistical properties,
such as the cluster number distribution and the cluster size distribution. In
this paper, we analyze these characteristics in the case in which vehicles are
located on a straight road. Assuming the Rayleigh fading model and a
probabilistic model of intervehicle distance, we derive probabilistic
distributions of the aforementioned connectivity characteristics, as well as
distributions of the biggest cluster and the number of disconnected vehicles.
All of the results are confirmed by simulations carried out for the realistic
values of parameters.

, , and (2019), Stochastic Geometry of Network of Randomly Distributed Moving
Vehicles on a Highway, Intern. J. Commun. Systems, 2019;00:1–6.

###### keywords:

vehicle network, clustering, network evolution

††articletype: Article Type

##  1 Introduction

Vehicular ad-hoc networks (VANETs) are developed to provide intelligent
transportation that improves road safety, reduces transportation time, and
decreases road congestion. One of the possible signal propagation scenarios is
vehicle-to-vehicle (V2V) communication, which means a direct connection
between two nearby vehicles. Another considered scenario is vehicle-to-
infrastructure (V2I) communication, in which a vehicle sends and receives a
signal from Road Side Units (RSUs) at the side of the road. RSUs form
infrastructure that improves the performance of VANETs and enhances
information distribution.

Multihop information propagation in VANET is limited by high vehicular speed
and frequent disconnections. To address these issues, significant progress was
made in routing protocols 1, 2, 3, 4, 5 developed to reduce latency and
stabilize the connection between vehicles, allowing vehicles to receive real-
time travel-related information.

The high mobility of vehicles and the changing network topology lead to
difficulties in the distribution of content, including images and video. To
mitigate these shortcomings, cloud architectures are proposed that allow for
the sharing of abundant resources, and for the regulating of network
connectivity specifications. In 6, the scheme, in which parked vehicles can
play the role of RSUs, is proposed. In 7, RSU cloud is introduced that has the
ability to dynamically reconfigure data forwarding rules in the network to
meet frequently changing service demands. Article 8 integrates both RSU and
vehicles into one computational network that provides services, and proposes
an optimal strategy for resource allocation.

VANET is a highly dynamic network, making it difficult to analyze from a
deterministic point of view. This is the reason that statistical methods are
used to make predictions under the assumption of particular probabilistic
distributions of vehicles on the road and models of signal transmission.
Considering this problem on a real map seems to be an impossible task, so, in
most articles, the authors limit themselves to the cases of a highway or a
crossroad. Articles 9, 10, 11, 12, 13 are devoted to the investigation of
statistical characteristics of the network, such as the number of clusters,
cluster size, and the number of disconnected vehicles. Let us discuss the
content of these articles in more detail. In 9, it is assumed that the
probabilistic distribution of distance between cars is known, and as soon as
the distance between neighboring vehicles does not exceed transmission range
$L$, they are always able to establish a connection. Under these assumptions,
probability distribution of information propagation distance is found, as well
as its expected value and variation. Articles 10, 13 discuss more realistic
communication channel models, such as Rayleigh, Rician, and Weibull fading
channels. For each of these models, the probability of connection between
neighboring vehicles and the probability of full connectivity of the network
are derived. The authors of 11 consider a more complex case in which the cars
are located between two RSUs, but use a simple the model of the communication
channel: they assume (as in 9) that the vehicles are always connected within
the constant transmission range and disconnected otherwise. Such advanced
network characteristics as expected number of clusters, network connectivity
probability, expected number of vehicles required for a VANET to be fully
connected, and expected network capacity are found within the framework of
this model. In 12, the authors suggest that vehicles in the cluster do not
communicate with each other, but, instead, all vehicles communicate with the
main car. Under presumption of the Rayleigh fading channel model, the Markov
model of packet transmission, and the PHY decoding failure model, the average
packet loss probability is derived. In 14, 15, 16, the evolution of the
network over time is considered. In 14, a case in which a signal is
transmitted to vehicles moving in the opposite direction and then sent back to
the initial side is investigated. The probability of such a successful two-hop
connection is derived under the assumption that vehicles move at a constant
speed. Article 15 focuses on the study of link duration in the specific case,
in which the speed of the vehicles grows linearly, reaches the limit, and then
remains constant. In our article 16, the network evolution is investigated
under the presumption that the connection between consecutive cars is
preserved according to the Markov probabilistic model; its parameters are
explicitly expressed through the network macro parameters. Under these
conditions, explicit formulas are obtained for the probability distribution of
link duration between consequent vehicles, the cluster existence, and other
fundamental network characteristics.

In this article, we consider the case of a vehicle network on a highway. Our
goal is to give an extended description of the network, not only in terms of
average values and variance, as is done in the previous articles (such as 11),
but to find the whole distribution of the number of clusters, cluster size,
and the number of disconnected vehicles (here we call them idle vehicles). We
use a more realistic model of the communication channel (Rayleigh fading
channel) than in 9, 11. We advance further than the authors of 10, 12, 13 and
find the distribution of the number of clusters and cluster size. The
difference between this article and 14, 15, 16 is that we do not consider the
evolution of the network over time, but concentrate in detail on more nuanced
connectivity characteristics at a fixed moment in time. Moreover, to the best
of our knowledge, we are the first to investigate the distribution of the
biggest cluster size. We assume that every car communicates only with the
nearest car, in front, and behind it. This assumption is made in the previous
statistical research in this area, in order to make the model simple enough to
analyze. Under the presumption of Rayleigh fading model and known density
function of intervehicle distance, we express the above mentioned
distributions in terms of known parameters of the network. Our model can be
applied to other scenarios if the probability of connection between every pair
of consecutive vehicles has the same value $p$. We derive several interesting
results describing network connectivity, namely, that the average size of
cluster is a constant, approximately equaling $1/(1-p)$; the average number of
clusters is proportional to the number of vehicles $n$ equaling
$1+(n-1)(1-p)$; and in the network, on average there is a constant fraction
$\approx 1/(1-p)^{2}$ of the idle vehicles. The average largest cluster of the
network, however, is not a constant and grows as $\log_{1/p}n$. We confirm
these and other theoretical results by simulations carried out in the cases of
different car densities and number of cars. The averages and variances of the
studied properties are summarized in the following table (where
$\gamma=0.577\ldots$ is an Euler constant):

Property | Average | Variance
---|---|---
Number of clusters | $1+(n-1)(1-p)$ | $(n-1)p(1-p)$
Size of clusters | $(1-p^{n})/(1-p)$ | $(2p^{n+1}n-2p^{n}n-p^{2n}-p^{n+1}+p^{n}+p)/(1-p)^{2}$
Size of the biggest cluster | $\approx\log_{1/p}\\{(n-1)(1-p)\\}+\gamma/\ln(1/p)+0.5$ | $\approx\pi^{2}/\ln^{2}(1/p)+\frac{1}{12}$
Number of idle cars | $2(1-p)+(n-2)(1-p)^{2}$ |  $-3np^{4}+10np^{3}-11np^{2}+4np+$ $8p^{4}-22p^{3}+18p^{2}-4p$

The studied statistical characteristics have several practical applications.
Cluster size distribution and maximum cluster size distribution provide us
with tools for network load prediction. Also, this statistical information can
be used for security purposes. The cluster size estimation, for example,
allows for the prediction of the number of infected vehicles in the event of
an attack on a cluster. Estimation of the number of disconnected vehicles
makes it possible to evaluate the quality of connection, in order to make the
percentage of disconnected vehicles acceptably low.

The paper is organized as follows. In section 2 we describe the connectivity
model and derive a formula for the probability of connection between two
consecutive cars. In 3 we derive all the probabilistic distributions mentioned
above. Finally, in section 4 the simulation results are presented and compared
to the calculations done by the formulas from section 3. The appendix is
devoted to a quick introduction to the theory of generating functions needed
in section 3. We use the following abbreviations:

We use the following abbreviations for the connectivity characteristics of the
network:

ClustNum | number of clusters in the network
---|---
ClustSize | size of the cluster
BiggestClust | size of the biggest cluster
IdleCars | number of disconnected vehicles

The following variables and functions are used:

$p$ | probability of connection between two consecutive vehicles
---|---
$G_{T}$ | transmit antenna gain
$G_{R}$ | receive antenna gain
$P_{tx}$ | transmit power
$\alpha$ | path loss exponent
$K$ | constant associated with the path loss model
$d$ | distance between cars
$C$ | speed of light
$W$ | thermal noise power
$k$ | Boltzmann constant, $k=1.38\times 10^{-23}J/K$
$f_{c}$ | carrier frequency
$T_{0}$ | room temperature
$B$ | transmission bandwidth
$f_{\gamma}(x)$ | Signal-to-noise ratio probability density function
$f_{d}(x)$ | intervehicle distance probability density function
$\binom{k}{s}$ | binomial coefficient
$coeff_{x^{n}}f(x)$ | coefficient of the term $x^{n}$ in series $f(x)$
$clust(N)$ | number of clusters of vehicle network $N$
$count(r;N)$ | number of clusters of vehicle network $N$ having size $r$
$Num(s,k)$ | number of networks having $k$ clusters, under condition that $s$ of them have size $r$.
$\gamma$ | Euler constant, $\gamma=0.577\ldots$

##  2 Network model

We consider the network of randomly distributed cars on a highway and suppose
that every vehicle of the network can establish a connection only with the
closest front and back neighbours. We denote by $d$ the distance between two
consecutive vehicles. We assume that the distance between vehicles has random
distribution with density function $f_{d}(x)$. Typically it is assumed that
distance between cars has exponential distribution, normal distribution, gamma
distribution or log-normal distribution. These distributions correspond to
different flow conditions. For instance, exponential distribution corresponds
to low traffic flow conditions, while high traffic flow conditions are
described by normal distribution.

Let us suppose that there are $n$ cars on a road and distribution of distance
between cars has a mean value $M$, therefore, the average distance between the
first and the last car is $(n-1)M$, consequently, despite the fact that the
road is long, it is improbable that the distance between first and last
vehicle is very large.

We assume that the average channel Signal-to-noise ratio (SNR) between two
consecutive vehicles $\overline{\gamma}$ is calculated using the following
formula from 13:

$\overline{\gamma}=\frac{P_{tx}K}{d^{\alpha}W},$ (1)

where $P_{tx}$ is the transmit power, $\alpha$ is the path loss exponent and
$K$ is a constant associated with the path loss model. The parameter $K$ is
given by the formula

$K=\frac{G_{T}G_{R}C^{2}}{(4\pi f_{c})^{2}},$ (2)

where $G_{T}$ and $G_{R}$ are the transmit and receive antenna gains, $C$ is
the speed of light and $f_{c}$ is the carrier frequency. We assume that the
antennas are omni directional

$G_{T}=G_{R}=1$ (3)

and $f_{c}=5.9GHz$. The thermal noise power can be determined by the formula

$W=kT_{0}B,$ (4)

where $k=1.38\times 10^{-23}J/K$ is the Boltzmann constant, $T_{0}$
($T_{0}=300^{\circ}K$) is the room temperature, and $B$ ($B=10$ MHz) is the
transmission bandwidth.

The SNR density function $f_{\gamma}(x)$ (is different from $f_{d}$) is given
by the formula

$f_{\gamma}(x)=\frac{1}{\overline{\gamma}}e^{-x/\overline{\gamma}}.$ (5)

It is assumed that there is a connection between consecutive cars if $\gamma$
is greater than the threshold $\Psi$. Therefore,

$P(\gamma>\Psi)=\int_{\Psi}^{\infty}f_{\gamma}(x)dx=e^{-\Psi/\overline{\gamma}}.$
(6)

Let us denote the probability that two consecutive vehicles can establish a
connection by $p$. By using (1), (6) and following the logic of 11 we can
derive formula expressing value of the probability $p$ through intervehicle
density function $f_{d}(x)$:

$p=\int_{0}^{\infty}P(\gamma(x)>\Psi)f_{d}(x)dx=\int_{0}^{\infty}e^{-\frac{\Psi
x^{\alpha}W}{P_{tx}K}}f_{d}(x)dx.$ (7)

The value of the parameter $p$ lies within an interval $(0,1)$ and plays a
prominent role in the further investigations. It determines the quality of
communication. A larger value of parameter $p$ leads to a better quality of
connection and consequently a bigger size of clusters.

##  3 Connectivity of the network

###  3.1 Distribution of number of clusters

The main achievement of this article is that we drive not only the average
values of various connectivity characteristics, as in the previous articles
(see the Introduction), but also their probability distributions and
variations. In this section we investigate such a significant characteristic
of the network as distribution of number of clusters. In this section we
establish the following formula for the probability of number of clusters
$ClustNum$ equaling $r$:

$P(ClustNum=r)=\binom{n-1}{r-1}p^{n-r}(1-p)^{r-1},$ (8)

where parameter $p$ is determined in (7). In other words, number of clusters
has binomial distribution.

To prove formula (8) we need to better understand the structure of the
network. The idea behind the derivation is to consider connections between
vehicles as Bernoulli trials. There are $n$ vehicles in the network
establishing no more than $n-1$ connections between each other, thus, we can
think about them as $n-1$ independent trials, in which each of them is
successful if consecutive vehicles can connect, and unsuccessful otherwise.
Clusters are assumed to be formed by a series of several consecutive vehicles,
therefore, if there are $r$ clusters in the network, then exactly $r-1$
vehicles cannot establish a connection. From the above we could conclude that
the probability $P(ClustNum=r)$ equals the probability that in $n-1$ trials
exactly $r-1$ are unsuccessful with the probability of success $p$. This
probability is given by well-known formula for binomial distribution, which in
our case is exactly formula (8).

By using known formulas for mean and variance of binomial distribution one can
derive that the average number of clusters and their variance are given by the
formulas

$\textbf{E}[ClustNum]=1+(n-1)(1-p),$ (9) $\textbf{Var}[ClustNum]=(n-1)p(1-p).$
(10)

Here we consider a disconnected vehicle as a separate cluster. The average
number of clusters formed by at least two cars can be calculated by
substracting from (9) the average number of idle cars derived later in section
3.4 (see the formula (40)). Therefore, it is given by the following formula:

$1+(n-1)(1-p)-2(1-p)-(n-2)(1-p)^{2}.$ (11)

Remark.  We see from (9) that the average number of clusters in the system
grows linearly with the number of vehicles with a growth factor $1-p$.

###  3.2 Distribution of cluster size

In this subsection we concentrate on finding probability
$P(ClustSize~{}=~{}r)$ that size of arbitrary cluster in the network is $r$.
We prove the following formula:

$P(ClustSize=r)=\begin{cases}p^{r-1}(1-p),&\mbox{if }r<n,\\\ p^{n-1},&\mbox{if
}r=n.\end{cases}$ (12)

Unfortunately, we do not know a simple proof. So, we use mathematical
apparatus of generating functions, which allows for reformulating the problem
in terms of series and use different analytical methods to operate with them.

We denote number of clusters of some vehicle network $N$ by $clust(N)$ and
number of clusters having size $r$ by $count(r;N)$. Let us establish the
formula

$P(ClustSize=r)=\sum_{k=1}^{n}\sum_{N:clust(N)=k}\frac{count(r;N)}{k}P(N),$
(13)

where $P(N)$ is the probability of the network $N$ and the second summation is
carried out over all networks $N$ with the restriction $clust(N)=k$. The
network could have from $1$ to $n$ clusters, so $k$ varies from $1$ to $n$,
and for each value of the parameter $k$, we consider all networks with $k$
clusters. Finally, for each such case, we sum up the probabilities that a
randomly selected cluster has size $r$. Since the network $N$ contains
$count(r;N)$ clusters of size $r$, then the probability of choosing a cluster
of size $r$ among $k$ clusters is given as $\frac{count(r;N)}{k}$. We should
multiply the last probability by $P(N)$ in order to calculate the probability
that the cluster of size $k$ is chosen in the network $N$.

We denote networks having $k$ clusters, under condition that $s$ of them have
size $r$ by $N_{s,k}$, therefore, $count(r;N_{s,k})=s$. Since the network has
$k$ clusters, $s$ varies from $1$ to $k$. We rewrite the summation over all
networks $N$ as the summation over the networks $N_{s,k}$:

$P(ClustSize=r)=\sum_{k=1}^{n}\sum_{s=1}^{k}\sum_{N_{s,k}}\frac{count(r;N_{s,k})}{k}P(N_{s,k})=\sum_{k=1}^{n}\sum_{s=1}^{k}\sum_{N_{s,k}}\frac{s}{k}P(N_{s,k}).$
(14)

Repeating the steps of derivation of (8) we conclude that

$P(N_{s,k})=(1-p)^{k-1}p^{n-k}.$ (15)

From (14) and (15) we derive that

$P(ClustSize=r)=\sum_{k=1}^{n}\sum_{s=1}^{k}\sum_{N_{s,k}}\frac{s}{k}(1-p)^{k-1}p^{n-k}=\sum_{k=1}^{n}(1-p)^{k-1}p^{n-k}\sum_{s=1}^{k}\frac{sNum(s,k)}{k}.$
(16)

Let us establish the formula

$Num(s,k)=coeff_{x^{n}}\binom{k}{s}x^{rs}(x+x^{2}+x^{3}+\ldots-x^{r})^{k-s},$
(17)

where $\binom{k}{s}$ is a binomial coefficient and $coeff_{x^{n}}P(x)$ is a
coefficient of the term $x^{n}$ in polynomial $P(x)$. The important step is to
understand that the value of $Num(s,k)$ equals the number of representations
of the integer $n$ as a sum of $k$ summands under condition that exactly $s$
of them equal $r$. This is true, because the clusters are formed by group of
consecutive vehicles. Therefore, we can use the formula (52) in the appendix
identical to the formula (17) for calculation of $Num(s,k)$.

From (17) we obtain

$\sum_{s=1}^{k}\frac{sNum(s,k)}{k}=\frac{1}{k}coeff_{x^{n}}\sum_{s=1}^{k}\binom{k}{s}sx^{rs}(x+x^{2}+x^{3}+\ldots-x^{r})^{k-s}=\frac{1}{k}coeff_{x^{n}}\sum_{s=1}^{k}\binom{k}{s}sx^{rs}\Big{(}\frac{x}{1-x}-x^{r}\Big{)}^{k-s}.$
(18)

The following identity

$\sum_{s=1}^{k}\binom{k}{s}sx^{rs}y^{k-s}=k(x^{r}+y)^{k-1}x^{r}$ (19)

holds. We substitute $y=\frac{x}{1-x}-x^{r}$ into (19) and use (18) to derive
the following:

$\sum_{s=1}^{k}\frac{sNum(s,k)}{k}\\\
=coeff_{x^{n}}\Big{(}x^{r}+\frac{x}{1-x}-x^{r}\Big{)}^{k-1}x^{r}\\\
=coeff_{x^{n}}\Big{(}\frac{x}{1-x}\Big{)}^{k-1}x^{r}.$ (20)

Combining equations (16) and (20) we obtain

$P(ClustSize=r)=\sum_{k=1}^{n}coeff_{x^{n}}\Big{(}\frac{x}{1-x}\Big{)}^{k-1}x^{r}(1-p)^{k-1}p^{n-k}\\\
=p^{n-1}\sum_{k=1}^{n}coeff_{x^{n}}\Big{(}\frac{x}{1-x}\Big{)}^{k-1}x^{r}\Big{(}\frac{1-p}{p}\Big{)}^{k-1}=p^{n-1}\sum_{k=1}^{\infty}coeff_{x^{n-r}}\Big{(}\frac{x}{(1-x)}\frac{1-p}{p}\Big{)}^{k-1}.$
(21)

The last equality in (21) is true because of the identity

$p^{n-1}\sum_{k=n+1}^{\infty}coeff_{x^{n-r}}\Big{(}\frac{x}{(1-x)}\frac{1-p}{p}\Big{)}^{k-1}=0.$

The last identity is valid, because the series
$\sum_{k=n+1}^{\infty}\Big{(}\frac{x}{(1-x)}\frac{1-p}{p}\Big{)}^{k-1}$ has
nonzero coefficients only of $x^{i},i\geq n$. Therefore, the coefficient of
$x^{n-r}$ equals zero. Applying the formula of the sum of geometric
progression to (21) we finally derive

$P(ClustSize=r)=p^{n-1}coeff_{x^{n-r}}\sum_{k=1}^{\infty}\Big{(}\frac{x}{(1-x)}\frac{1-p}{p}\Big{)}^{k-1}=p^{n-1}coeff_{x^{n-r}}\frac{1}{1-\frac{x}{(1-x)}\frac{1-p}{p}}=p^{n-1}coeff_{x^{n-r}}\frac{1-x}{1-x(1+\frac{1-p}{p})}\\\
=p^{n-1}coeff_{x^{n-r}}(1-x)\sum_{k=0}^{\infty}x^{k}p^{-k}=\begin{cases}p^{n-1}\Big{(}p^{-n+r}-p^{-n+r+1}\Big{)},&\mbox{if
}r<n,\\\ p^{n-1},&\mbox{if }r=n\end{cases}=\begin{cases}p^{r-1}(1-p),&\mbox{if
}r<n,\\\ p^{n-1},&\mbox{if }r=n.\end{cases}$ (22)

Formula (12) is proven.

Now we are ready to calculate the average value and the variance of the
distribution (12). Firstly, we find the first two moments by the formulas

$M_{1}=np^{n-1}+(1-p)\sum_{\alpha=1}^{n-1}\alpha p^{\alpha-1},$ (23)
$M_{2}=n^{2}p^{n-1}+(1-p)\sum_{\alpha=1}^{n-1}\alpha^{2}p^{\alpha-1}.$ (24)

By using the identities

$\sum_{\alpha=1}^{n-1}\alpha
p^{\alpha-1}=\frac{1+((n-1)p-n)p^{n-1}}{(1-p)^{2}},$ (25)
$\sum_{\alpha=1}^{n-1}\alpha^{2}p^{\alpha-1}=\frac{1+p-((n-1)^{2}p^{2}+(-2n^{2}+2n+1)p+n^{2})p^{n-1}}{(1-p)^{3}}$
(26)

we derive

$M_{1}=\frac{1-p^{n}}{1-p},$ (27)
$M_{2}=\frac{2p^{n+1}n-2p^{n}n-p^{n+1}-p^{n}+p+1}{(1-p)^{2}}.$ (28)

From formulas (27) and (28) we finally obtain

$\textbf{E}[ClustSize]=M_{1}=\frac{1-p^{n}}{1-p},$ (29)
$\textbf{Var}[ClustSize]=M_{2}-M_{1}^{2}\\\
=\frac{2p^{n+1}n-2p^{n}n-p^{2n}-p^{n+1}+p^{n}+p}{(1-p)^{2}}.$ (30)

Remark.  It follows from the formulas (29) and (30) that the average cluster
size approximately equals ${1/(1-p)}$ (constant!) and its variance
approximately equals ${p/(1-p)^{2}}$, since $p^{n}\rightarrow
0,n\rightarrow\infty$.

###  3.3 Distribution of size of the biggest cluster

In the previous section we explored the question about distribution of cluster
size. It was established that the average cluster size is approximately
$1/(1-p)$, therefore, being a constant. However, it is natural to investigate
this problem further and understand how the biggest cluster size deviates from
the average.

We denote the probability that the biggest cluster of the network is formed by
exactly $r$ vehicles by $P(BiggestClust=r)$. Let us introduce the function
$g(k)$ by the formula

$g(k)=\sum_{m=1}^{\lfloor
n/(k+1)\rfloor}(-1)^{m}p^{mk}(1-p)^{m-1}\binom{n-1-mk}{m-1}+\sum_{m=0}^{\lfloor(n-1)/(k+1)\rfloor}(-1)^{m}p^{mk}(1-p)^{m}\binom{n-1-mk}{m},$
(31)

where $\lfloor x\rfloor$ means rounding down. We establish the following
formula:

$P(BiggestClust=r)=\begin{cases}p^{n-1},&\mbox{if }r=n,\\\
g(r)-g(r-1),&\mbox{if }1<r<n,\\\ (1-p)^{n-1},&\mbox{if }r=1.\end{cases}$ (32)

The cases $r=n$ and $r=1$ are trivial, because under these assumptions all
vehicles are connected or disconnected respectively. Thus, we need to prove
the formula (32) only in the case $1<r<n$. We denote the length of the longest
run in the sequence of the $n-1$ Bernoulli trials by $L_{n-1}$. The following
equation is proven in 17:

$P(L_{n-1}\leq r-1)=g(r),$ (33)

where $g(k)$ is determined in (31). The number of vehicles in the network is
$n$, consequently, there are $n-1$ connections between them and we can
consider them as a sequence of Bernoulli trials as before. The probability
$P(BiggestClust~{}=~{}r)$ equals the probability that the network has the
longest run equaling $r-1$ (number of connections is less than the number of
cars by 1). From the above we derive

$P(BiggestClust=r)=P(L_{n-1}=r-1)=P(L_{n-1}\leq r-1)-P(L_{n-1}\leq
r-2)=g(r)-g(r-1).$ (34)

Formula (32) is proven.

Despite the fact that the formula (32) can be used for the calculation of the
probabilities $P(BiggestClust=r)$, it is overly complicated to predict
behaviour of the network with growing number of vehicles. Fortunately, the
asymptotic of the longest successful run with growing $n$ is already found
(see, for instance, 18, 19), in which the expected value and variance of the
longest run are derived for sufficiently large values of $n$. Taking into
account that there are $n-1$ connections in our network, the fact that the
number of connections is less that number of vehicles in a cluster by 1, and
the formulas from above mentioned articles we conclude that the expected value
and variation of the biggest cluster of our network are

$\textbf{E}[BiggestClust]\approx\log_{1/p}\\{(n-1)(1-p)\\}+\gamma/\ln(1/p)+0.5,$
(35) $\textbf{Var}[BiggestClust]\approx\pi^{2}/\ln^{2}(1/p)+\frac{1}{12},$
(36)

where $\gamma=0.577\ldots$ is Euler constant.

Remark.  From the previous derivations we know that the size of average
cluster is $\approx 1/(1-p)$, however, from (35) we conclude that the size of
the biggest cluster grows logarithmically with the number of vehicles in the
network.

###  3.4 Distribution of idle cars

We call car idle if it fails to establish a connection with any of its
neighbours. Unfortunately, there is no simple formula for calculating the
probability $P(IdleCars=r)$ that the network has exactly $r$ idle vehicles.
One can prove by using method of generating functions that the probability
$P(IdleCars=r)$ can be calculated by the formula

$P(IdleCars=r)=\frac{p^{n+1}}{(1-p)^{2}(n-r)!}\frac{d^{n-r}}{dx^{n-r}}\Big{\\{}\frac{(1-x)^{r+1}}{((1-x)\frac{p}{1-p}-x^{2})^{r+1}}\Big{\\}}\Big{|}_{x=0},$
(37)

where $k!$ is a factorial and $\frac{d^{k}}{dx^{k}}f(x)\Big{|}_{x=0}$ is
$k$-th derivative of function $f(x)$ calculated at the point $x=0$.

Despite the fact that the formula of distribution of idle vehicles is quite
complicated, we can derive simple expression for the expected number of idle
vehicles. Let us denote the random variable that equals 1 it $k$-th vehicle is
idle and 0 otherwise by $\xi_{k}$. Let us prove that the expected value of
$\xi_{k}$ can be determined by this formula

$\textbf{E}[\xi_{k}]=\begin{cases}(1-p)^{2},&\mbox{if }1<k<n,\\\ 1-p,&\mbox{if
}k=1\mbox{ or }k=n.\end{cases}$ (38)

It is enough to prove the formula (38) in the case $1<k<n$, since derivations
in other cases are similar. If $1<k<n$, then the vehicle has exactly two
neighbours. Since the probability of disconnection equals $(1-p)$, the
probability of car being idle equals $(1-p)^{2}$. Finally,

$\textbf{E}[\xi_{k}]=1\times(1-p)^{2}+0\times(1-(1-p)^{2})=(1-p)^{2}.$ (39)

Formula (38) is proven.

Therefore, the average number of idle cars is

$\textbf{E}[IdleCars]=\textbf{E}\Big{[}\sum_{k=1}^{n}\xi_{k}\Big{]}=\sum_{k=1}^{n}\textbf{E}[\xi_{k}]=2(1-p)+\sum_{k=2}^{n-1}(1-p)^{2}=2(1-p)+(n-2)(1-p)^{2}.$
(40)

Let us calculate the variance of distribution of number of idle vehicles. It
is given as follows:

$\textbf{Var}[IdleCars]=\textbf{E}\Big{[}\Big{(}\sum_{k=1}^{n}\xi_{k}\Big{)}^{2}\Big{]}-\textbf{E}[IdleCars]^{2}=\sum_{k=1}^{n}\textbf{E}[\xi_{k}^{2}]+2\sum_{k=1}^{n-1}\sum_{j=k+1}^{n}\textbf{E}[\xi_{k}\xi_{j}]-\textbf{E}[IdleCars]^{2}.$
(41)

If two vehicles are not neighbours then the variables $\xi_{k}$ and $\xi_{j}$
are independent and, therefore,
$\textbf{E}[\xi_{k}\xi_{j}]=\textbf{E}[\xi_{k}]\textbf{E}[\xi_{j}]$, otherwise
$\textbf{E}[\xi_{k}\xi_{j}]$ can be derived by analogy with the derivation of
(38). Thus, the following formula

$\textbf{E}[\xi_{k}\xi_{j}]=\begin{cases}\textbf{E}[\xi_{k}]\textbf{E}[\xi_{j}],&\mbox{if
}j>k+1,\\\ (1-p)^{3},&\mbox{if }j=k+1\mbox{ and }1<k<n-1,\\\
(1-p)^{2},&\mbox{if }j=k+1,k=1\mbox{ or }k=n-1.\end{cases}$ (42)

is valid for $j>k$. Using (42), we derive

$\sum_{k=1}^{n-1}\sum_{j=k+1}^{n}\textbf{E}[\xi_{k}\xi_{j}]=\sum_{k=1}^{n-1}\sum_{j=k+2}^{n}\textbf{E}[\xi_{k}]\textbf{E}[\xi_{j}]+\sum_{k=2}^{n-2}\textbf{E}[\xi_{k}\xi_{k+1}]+\textbf{E}[\xi_{1}\xi_{2}]+\textbf{E}[\xi_{n-1}\xi_{n}]\\\
=\sum_{k=1}^{n-1}\sum_{j=k+2}^{n}\textbf{E}[\xi_{k}]\textbf{E}[\xi_{j}]+(n-3)(1-p)^{3}+2(1-p)^{2}.$
(43)

Consequently, by using (41) and (43) and taking into account that
$\xi_{k}^{2}=\xi_{k}$ we obtain

$\textbf{Var}[IdleCars]=\sum_{k=1}^{n}\textbf{E}[\xi_{k}^{2}]+2\sum_{k=1}^{n-1}\sum_{j=k+1}^{n}\textbf{E}[\xi_{k}\xi_{j}]-\Big{(}\textbf{E}\Big{[}\sum_{k=1}^{n}\xi_{k}\Big{]}\Big{)}^{2}=\sum_{k=1}^{n}\textbf{E}[\xi_{k}^{2}]+2\sum_{k=1}^{n-1}\sum_{j=k+2}^{n}\textbf{E}[\xi_{k}]\textbf{E}[\xi_{j}]\\\
+2(n-3)(1-p)^{3}+4(1-p)^{2}-\sum_{k=1}^{n}\textbf{E}[\xi_{k}]^{2}-2\sum_{k=1}^{n-1}\sum_{j=k+1}^{n}\textbf{E}[\xi_{k}]\textbf{E}[\xi_{j}]\\\
=\sum_{k=1}^{n}\textbf{E}[\xi_{k}]-\sum_{k=1}^{n}\textbf{E}[\xi_{k}]^{2}-2\sum_{k=1}^{n-1}\textbf{E}[\xi_{k}]\textbf{E}[\xi_{k+1}]+2(n-3)(1-p)^{3}+4(1-p)^{2}.$
(44)

We derive the explicit formulas for every sum in (44). Equality (40) gives us

$\sum_{k=1}^{n}\textbf{E}[\xi_{k}]=2(1-p)+(n-2)(1-p)^{2}.$ (45)

By analogy with derivation of (40) we can deduce that

$\sum_{k=1}^{n}\textbf{E}[\xi_{k}]^{2}=2(1-p)^{2}+(n-2)(1-p)^{4}.$ (46)

The last sum in (44) can be derived using (38) as follows:

$\sum_{k=1}^{n-1}\textbf{E}[\xi_{k}]\textbf{E}[\xi_{k+1}]=\sum_{k=2}^{n-2}\textbf{E}[\xi_{k}]\textbf{E}[\xi_{k+1}]+\textbf{E}[\xi_{1}]\textbf{E}[\xi_{2}]+\textbf{E}[\xi_{n-1}]\textbf{E}[\xi_{n}]=(n-3)(1-p)^{4}+2(1-p)^{3}.$
(47)

Combining formulas (44)–(47) we obtain

$\textbf{Var}[IdleCars]=2(1-p)+(n-2)(1-p)^{2}-2(1-p)^{2}-(n-2)(1-p)^{4}-2(n-3)(1-p)^{4}-4(1-p)^{3}\\\
+2(n-3)(1-p)^{3}+4(1-p)^{2}=-3np^{4}+10np^{3}-11np^{2}+4np+8p^{4}-22p^{3}+18p^{2}-4p.$
(48)

Remark.  Formula (40) leads us to an important conclusion that in average
there is a constant fraction $\approx(1-p)^{2}$ of idle vehicles in the
network.

##  4 Simulations

We use the model described by the equations (1)–(7) and make simulations for
$n=10$ and $n=20$ vehicles on the road assuming exponential distribution of
the intervehicle distance with probability density function $\rho e^{-\rho
x}$, where $\rho$ is a vehicle density. We consider two scenarios. In the
first case, the density is low $\rho=0.01$ (10 cars per kilometer) which is
the case for rural traffic, in the second case $\rho=0.05$ (50 cars per
kilometer), which corresponds to urban traffic. We use the following values of
the parameters taken from 13 (see their descriptions in section 2):

$G_{T}=G_{R}=1$, $f_{c}=5.9$ GHz, $\alpha$ = 2.5,

$T_{0}=300^{\circ}$ K, $B=10$ MHz, $\Psi=10$ dB.

We use the transmit power value $P_{tx}=4$ dBm. This value corresponds to the
message transmission over short distances up to 99 meters 20.

Graphs of all network properties are plotted in the cases a) $n=10$,
$\rho=0.01$, b) $n=10,\rho=0.05$, c) $n=20,\rho=0.01$, and d)
$n=20,\rho=0.05$. In the case $\rho=0.01$, the probability of connection
between two consecutive cars is $p=0.5576$, in the case $\rho=0.01$ the
probability of connection is $p=0.9525$ (it does not depend on the number of
cars within the model). To verify the obtained formulas we generate sequence
of vehicles randomly 100000 times and then calculate network connectivity
distributions from these sample values. The graphics of the simulations are
identical to the results obtained by the formulas, which confirms their
correctness.

On Figure 1 the graphics of distribution (8) of number of clusters are
presented as well as simulation results for scenarios a, b, c, and d. As one
can see from Figure 1 in cases a and c, the connection is less stable
($p=0.5576$), therefore, the network on average has sufficiently large number
of clusters 4.9812 and 9.4048, respectively. In cases b and d, the connection
occurs with a much higher probability $p=0.9525$, therefore, in these cases
the network tends to be fully connected, which is expressed in a high
probability of having 1–3 clusters. In these cases average cluster sizes are
significantly lower: 1.4278 and 1.9031, respectively. The fact that the
distributions on Figure 1 are close to normal is not surprising, since, as
well known, the binomial distribution in the case $n\to\infty$ tends to
normal.

((a)) $n=10$,$\rho=0.01$

((b)) $n=10$,$\rho=0.05$

((c)) $n=20$,$\rho=0.01$

((d)) $n=20$,$\rho=0.05$

Figure 1: Probability distribution of number of clusters in the network for
different $n$ and $\rho$

On Figure 2 the graphs of the cluster size distribution are presented. The
results are obtained by the formula (12) and compared to the simulation
results. The same logic as in the previous paragraph leads us to the
conclusion that the cluster size in cases a and c should be smaller than in
cases b and d. This hypothesis is confirmed, because the average cluster sizes
in these cases are 2.2541 and 2.2606, respectively, versus 8.1109 and 13.0948
in cases b and d. The decrease in graphs in cases a and c is explained by the
fact that in these cases the connection is less stable, therefore, the
probability of a cluster having a size $i$ decreases with increasing of the
parameter $i$. Graphs b and d have a pronounced peak, since in the case of
stable connection, the network tends to be fully connected, therefore, a
cluster having size $n$ has a maximum probability.

((a)) $n=10$,$\rho=0.01$

((b)) $n=10$,$\rho=0.05$

((c)) $n=20$,$\rho=0.01$

((d)) $n=20$,$\rho=0.05$

Figure 2: Probability distribution of size of a cluster for different $n$ and
$\rho$

The graphs of the biggest cluster size distribution obtained by the formula
(32) and simulation results are presented on Figure 3. By analogy with the
explanation of the behavior of cluster size distribution graphs, the behavior
of the maximum cluster size distribution can be explained. The average values
in cases a, b, c, and d equal 4.0425, 8.9048, 5.2253, and 16.0203,
respectively.

((a)) $n=10$,$\rho=0.01$

((b)) $n=10$,$\rho=0.05$

((c)) $n=20$,$\rho=0.01$

((d)) $n=20$,$\rho=0.05$

Figure 3: Probability distribution of size of the biggest cluster for
different $n$ and $\rho$

Finally, the graphs of the distribution of the number of idle vehicles are
depicted on Figure 4. The simulation results are compared with the
probabilities calculated by the formula (37). One can explain the graphs of
distribution of idle cars by the fact that in the case of a more stable
connection (graphs b and d), the probability that idle cars are not present in
the network is close to 1, and then quickly decreases with increasing their
number. On graphs a and c, on the contrary, the probability of the network
having several idle vehicles is quite high due to the less stable connection.
This can also be explained in terms of average numbers of idle vehicle in
cases a, b, c, and d equaling 2.4501, 0.1131, 4.4069, and 0.1357 respectively.

((a)) $n=10$,$\rho=0.01$

((b)) $n=10$,$\rho=0.05$

((c)) $n=20$,$\rho=0.01$

((d)) $n=20$,$\rho=0.05$

Figure 4: Probability distribution of the number of idle vehicles for
different $n$ and $\rho$

##  5 Conclusion

This article aims to present versatile investigation of VANET connectivity
properties in the case, in which vehicles are distributed on the road. It
expresses results in terms of parameters of known probability distributions of
intervehicle distance and connectivity model. We derive distributions of
number of clusters, cluster size, size of the biggest cluster and number of
idle vehicles as well as calculate expected value and variance of every of
these characteristics. The results are confirmed by the simulations in the
cases of urban and rural traffic flow.

##  Appendix A Method of generating functions

In this subsection we give a brief exposition of mathematical method of
generating functions. This elegant and effective method is used in this
article to obtain distribution of cluster size. The essence of the method is
that it treats infinite sequence of numbers $a_{k}$ as the coefficients of a
power series $\sum_{k=0}^{\infty}a_{k}x^{k}.$ We use this method to establish
a relation between the number of different representations of integer number
as a sum of integer numbers with some restrictions. We touch only one aspect
of this theory, for other applications we recommend to read the book 21. Here
we formulate several problems in ascending order of complexity. We need the
result of the problem 3, but in order to obtain it we solve problems 1 and 2.

Problem 1. What is the number of representations of positive integer number
$n$ as a sum of $k$ positive integer summands?

We denote coefficient of $x^{n}$ in power series $F(x)$ by
$coeff_{x^{n}}F(x)$. It is stated that the following coefficient:

$coeff_{x^{n}}\Big{(}x+x^{2}+x^{3}+\ldots\Big{)}^{k}$ (49)

equals the number of representations of $n$ as $k$ summands. Let us prove it.
We can rewrite the sum $\Big{(}x+x^{2}+x^{3}+\ldots\Big{)}^{k}$ as

$\sum_{\alpha_{1}}\sum_{\alpha_{2}}\ldots\sum_{\alpha_{k}}x^{\alpha_{1}}x^{\alpha_{2}}\ldots
x^{\alpha_{k}}=\sum_{\alpha_{1}}\sum_{\alpha_{2}}\ldots\sum_{\alpha_{k}}x^{\alpha_{1}+\alpha_{2}+\ldots+\alpha_{k}}=\sum_{l=1}^{\infty}x^{l}\sum_{\alpha_{1}+\alpha_{2}+\ldots+\alpha_{k}=l}1.$
(50)

Therefore, the coefficient (49) indeed equals the number of representations of
$n$ as a sum of $k$ summands $\alpha_{1},\alpha_{2},\ldots,\alpha_{k}$.

Problem 2. What is the number of representations of positive integer number
$n$ as a sum of $k$ positive integer summands not equal $r$?

The following formula gives the solution to this problem:

$coeff_{x^{n}}\Big{(}x+x^{2}+x^{3}+\ldots-x^{r}\Big{)}^{k}.$ (51)

The only difference between formulas (49) and (51) is that we subtract $x^{r}$
because integer summands do not equal $\alpha$. It can be proven similar to
the proof of formula (49).

Problem 3. What is the number of representations of positive integer number
$n$ as a sum of $k$ positive integer summands with the following restriction:
among them exactly $s$ numbers equal $r$?

The answer is given by the following equality:

$coeff_{x^{n}}\Big{\\{}\binom{k}{s}x^{rs}\Big{(}x+x^{2}+x^{3}+\ldots-x^{r}\Big{)}^{k-s}\Big{\\}}.$
(52)

because we can choose $s$ summands equaling $r$ in $\binom{k}{s}$ ways and
other $k-s$ summands should not be equal $r$ (see problem 2).

## References

  * 1 Hashem Eiza M, Owens T, Ni Q. Secure and Robust Multi-Constrained QoS Aware Routing Algorithm for VANETs. IEEE Trans. Dependable Secure Comput. 2016; 13(1): 32–45.
  * 2 Rak J. LLA: A New Anypath Routing Scheme Providing Long Path Lifetime in VANETs. IEEE Commun. Lett. 2014; 18(2): 281–284.
  * 3 Sahu PK, Wu EH, Sahoo J, Gerla M. BAHG: Back-Bone-Assisted Hop Greedy Routing for VANET’s City Environments. IEEE Trans. Intell. Transp. Syst. 2013; 14(1): 199–213.
  * 4 Wu Q, Liu Q, Zhang L, Zhang Z. A trusted routing protocol based on GeoDTN+Nav in VANET. China Commun. 2014; 11(14): 166–174.
  * 5 Eiza MH, Ni Q. An Evolving Graph-Based Reliable Routing Scheme for VANETs. IEEE Trans. Veh. Technol. 2013; 62(4): 1493–1504.
  * 6 Gong H, Yu L, Liu N, Zhang X. Mobile content distribution with vehicular cloud in urban VANETs. China Commun. 2016; 13(8): 84–96.
  * 7 Salahuddin MA, Al-Fuqaha A, Guizani M. Software-Defined Networking for RSU Clouds in Support of the Internet of Vehicles. IEEE Internet Things J. 2015; 2(2): 133–144.
  * 8 Lin C, Deng D, Yao C. Resource Allocation in Vehicular Cloud Computing Systems With Heterogeneous Vehicles and Roadside Units. IEEE Internet Things J. 2018; 5(5): 3692–3700.
  * 9 Wang BX, Adams TM, Jin W, Meng Q. The process of information propagation in a traffic stream with a general vehicle headway: A revisit. Transp. Research Emerg. Technol. C 2010; 18(3): 367 – 375.  11th IFAC Symposium: The Role of Control.
  * 10 Babu A, Muhammed Ajeer V. Analytical model for connectivity of vehicular ad hoc networks in the presence of channel randomness. Int. J. Commun. Syst. 2013; 26(7): 927–946.
  * 11 Kwon S, Kim Y, Shroff NB. Analysis of Connectivity and Capacity in 1-D Vehicle-to-Vehicle Networks. IEEE Trans. Wireless Commun. 2016; 15(12): 8182–8194.
  * 12 Wang H, Liu RP, Ni W, Chen W, Collings IB. VANET Modeling and Clustering Design Under Practical Traffic, Channel and Mobility Conditions. IEEE Trans. Wireless Commun. 2015; 63(3): 870–881.
  * 13 Chandrasekharamenon NP, AnchareV B. Connectivity analysis of one-dimensional vehicular ad hoc networks in fading channels. EURASIP Journal on Wireless Communications and Networking 2012; 2012(1).
  * 14 Kesting A, Treiber M, Helbing D. Connectivity Statistics of Store-and-Forward Intervehicle Communication. IEEE Trans. Intell. Transp. Syst. 2010; 11(1): 172–181.
  * 15 Yan G, Olariu S. A Probabilistic Analysis of Link Duration in Vehicular Ad Hoc Networks. IEEE Trans. Intell. Transp. Syst. 2011; 12(4): 1227–1236.
  * 16 Dubosarskii G, Primak SL, Wang X. Evolution of Vehicle Network on a Highway. IEEE Transactions on Vehicular Technology 2019; 68(9): 9088–9097.
  * 17 Pekoz EA, Ross SM. A simple derivation of exact reliability formulas for linear and circular consecutive-k-of-n: F systems. Journal of Applied Probability 1995; 32(2): 554–557.
  * 18 Schilling MF. The Longest Run of Heads. College Math. J. 1990; 21(3): 196–207.
  * 19 Gordon L, Schilling MF, Waterman MS. An extreme value theory for long head runs. Probab. Theory Relat. Fields 1986; 72(2): 279–287.
  * 20 Rezha FP, Siadari TS, Shin SY. Adaptive transmission power in cluster-based routing VANET. In: Jeju Island, Korea. ; 2012: 539–543.
  * 21 Wilf HS. Generatingfunctionology.  CRC Press, 3rd edition .  2005\.