\section{Introduction}
\label{sec_introduction}

Counting statistics is a scheme to calculate all statistics related to specific transitions
in a stochastic system.
In the counting statistics, a master equation with discrete states is used
to derive time-evolution equations for generating functions related to the specific transitions.
The scheme has been used to investigate F{\"o}rster resonance energy transfer,
and many successful results have been obtained \cite{Gopich2003,Gopich2005,Gopich2006}.
Although the scheme is basically formulated for a system with a finite number of states,
it is possible to use the scheme to investigate a system with an infinite number of states.
However, as exemplified later, we have non-closed equations in general,
so that it would be needed to develop approximation schemes suitable for specific systems.
As a first step, it is important to check
whether an approximation scheme for the counting statistics is available for
the system with an infinite number of states or not.

In the present paper, we focus on dynamics in genetic switches.
It has been shown that stochastic behavior
plays an important role in gene regulatory systems \cite{Elowitz2002,Rao2002,Kaern2005},
and there are many studies
for the stochasticity in the gene regulatory systems
from experimental points of view (e.g., see \cite{Gardner2000,Okano2008})
and theoretical ones (e.g., see
\cite{Hasty2000,Sasai2003,Hornos2005,Xu2006,Schultz2007,Shahrezaei2008,Walczak2009,Venegas-Ortiz2011}).
Not only studies by numerical simulations,
but also those by analytical calculations have been performed.
Some analytical expressions for
the static properties, i.e., stationary distributions
for the number of proteins or mRNAs, have already been obtained.
In addition, in order to investigate the role of the stochasticity in genetic switches,
dynamical properties, i.e., switching behavior between active and inactive gene states,
have also been studied.
Basically, such dynamical properties have been investigated
by numerical simulations (e.g., see \cite{Feng2011});
only for a simple system, analytical expressions for the first-passage time distribution
have been obtained \cite{Visco2009}.
The genetic switch is described by a master equation with an infinite number of states.
Hence, if we can use the scheme of the counting statistics
in order to investigate the dynamical properties in the genetic switches,
it will be helpful to obtain deeper understanding and intuitive pictures
for the genetic switches.

The aim of the present paper is
to seek the applicability of the counting statistics
in order to investigate the dynamical property in the genetic switches.
It immediately becomes clear that
a straightforward application of the counting statistics derives intractable non-closed equations.
In order to obtain simple closed forms,
we here employ an effective interaction approximation \cite{Ohkubo2011}.
As a result, we will show that the switching problem
can be treated as a simple two-state model approximately.
This result immediately gives us intuitive understanding for the switching behavior
and the non-Poissonian property.

The present paper is constructed as follows.
In Sec.~\ref{sec_model}, we give a brief explanation of a stochastic model for the genetic switch.
In Sec.~\ref{sec_counting_statistics},
the counting statistics is employed in order to count
the number of transitions in the genetic switch,
and, as a result, a simple two-state model is derived approximately.
The derived approximated results
are compared with those of Monte Carlo simulations in Sec.~\ref{sec_results}.
Section \ref{sec_conclusions} gives concluding remarks.

\section{Model}
\label{sec_model}

\begin{figure}
\begin{center}
\includegraphics[width=60mm]{fig1.eps}
\end{center}
\caption{
A schematic illustration of the self-regulating gene with repressed binding interaction.
When the regulatory proteins are combining the gene,
the gene is repressed and there is no production of proteins.
If the regulatory proteins are released from the gene,
the gene becomes active and it can produce the proteins.
We consider the transition between the active and inactive states
as a switch.
}
\label{fig_model}
\end{figure}

A gene regulatory system consists of many components, such as genes, RNAs, and proteins.
Here, a simplified model is used; mRNAs are neglected for simplicity,
and an activated gene assumes to directly increase the number of proteins.
In addition, in the simplified model, a repressed gene cannot produce any proteins.
The above model has been used to investigate the switching behavior
in previous works,
and, for example, see \cite{Schultz2007} for details of the model.

We summarize the model studied in the present paper in Fig.~\ref{fig_model}.
The binding interaction is assumed to be a repressed one,
and the gene is activated only when the regulatory proteins are not binding the gene.
The proteins are produced from the gene in the active state with rate $g$,
and proteins are degraded spontaneously with rate $d$.
The regulatory proteins bind the gene with a rate function $\mathcal{H}(n)$,
where $n$ is the number of free proteins.
For example, $\mathcal{H}(n) = h n$ for a monomer interaction case,
and $\mathcal{H}(n) = h n(n-1)/2$ for a dimer interaction case,
where $h$ is a rate constant for the binding.
$f$ is a rate constant with which the regulatory proteins are released
from the repressor site of the gene.

We here give short comments for the model from the viewpoint of experiments.
Using this simplified model,
we can discuss the connection among the model parameters,
the number of proteins, and the switching behaviors.
While the number of proteins $n$ can be observed or estimated experimentally,
as far as we know, there has not been an experimental technique
to observe the attachment and detachment of the regulatory proteins directly.
We hope that developments of single-molecule observations in future
would enable us to give information about the switching dynamics.

\section{Counting statistics for the number of transitions}
\label{sec_counting_statistics}

\subsection{Master equation for the number of proteins}
\label{subsec_master_equation}

Analytical treatments for the self-regulating gene system have been developed,
and an exact solution is known for the monomer interaction case, i.e.,
$\mathcal{H}(n) = h n$ \cite{Hornos2005,Visco2009}.
In order to simplify the analytical treatments,
an additional assumption has been used in some previous works \cite{Schultz2007,Ohkubo2011};
i.e., some of proteins are assumed to be inert when the gene state is active.
The inert proteins cannot repress the gene, and it is not degraded.
For the monomer interaction case, there is only one inert protein;
the number of inert protein for the dimer interaction case is two, and so on.
Note that the assumption of the inert proteins does not have physical meanings;
this only simplify the analytical treatments (for details, see \cite{Schultz2007}).
However,
it has been shown that this assumption has little influence of the gene system,
and then we employ the assumption in the present paper.

Let $\alpha_n$ and $\beta_n$ be states
in which there are $n$ free proteins for the active and inactive states, respectively.
The probabilities for $\alpha_n$ and $\beta_n$ at time $t$ satisfy
the following master equations;
\begin{align}
\frac{d P(\alpha_n,t)}{dt} =
& g [ P(\alpha_{n-1},t) - P(\alpha_n,t) ] \nonumber \\
& + d [(n+1) P(\alpha_{n+1},t) - n P(\alpha_n,t) ] \nonumber \\
&- h n P(\alpha_n,t) + f P(\beta_n,t),
\label{eq_master_monomer_1_exact}\\
\frac{d P(\beta_n,t)}{dt} =
& d [(n+1)P(\beta_{n+1},t) - n P(\beta_n,t)]  \nonumber \\
&+ h n P(\alpha_n,t) - f P(\beta_n,t),
\label{eq_master_monomer_2_exact}
\end{align}
where $P(\alpha_n,t)$ and $P(\beta_n,t)$
are probabilities for $n$ free proteins for the active and inactive states,
respectively.

As stated in Sec.~\ref{sec_introduction},
the exact solutions for stationary distributions
of the number of proteins have been derived,
and those are expressed
using the Kummer confluent hypergeometric functions.
For details, see \cite{Hornos2005,Schultz2007}.

\subsection{Counting statistics}
\label{subsec_counting_statistics}

Using the concept of the counting statistics \cite{Gopich2003,Gopich2005,Gopich2006},
it is possible to investigate dynamical properties, i.e.,
all statistics for the switching behavior between the active and inactive states.
In the present paper, as an example, we calculate
the number of transitions from the inactive state to the active state.
The generating functions for the transitions
are immediately obtained from the master equations \eqref{eq_master_monomer_1_exact}
and \eqref{eq_master_monomer_2_exact}.
A brief explanation of the counting statistics is given in the Appendix,
and we here give consequences of the counting statistics.

A probability, with which there are $k$ transitions from the inactive state to the active state
during time $t$, is denoted by $P(k|t)$.
The generating function for $P(k|t)$ is defined as
\begin{align}
F(\lambda,t) = \sum_{k=0}^{\infty} P(k|t) \lambda^{k},
\end{align}
where $\lambda$ is a counting variable.
The generating function gives all information related to ``inactive $\to$ active'' transitions.
According to the scheme of counting statistics,
we split $F(\lambda,t)$ into restricted generating functions
$\{\phi(\alpha_n,\lambda,t)\}$ and $\{\phi(\beta_n,\lambda,t)\}$,
where $\phi(\alpha_n,\lambda,t)$ and $\phi(\beta_n,\lambda,t)$ are
the generating functions for the system
in states $\alpha_n$ and $\beta_n$ at time $t$, respectively.
Using the scheme of the counting statistics,
we obtain the following time-evolution equations for the restricted generating functions
$\{\phi(\alpha_n,\lambda,t)\}$ and $\{\phi(\beta_n,\lambda,t)\}$:
\begin{align}
\frac{d \phi(\alpha_n,\lambda,t)}{dt} =
& g [ \phi(\alpha_{n-1},\lambda,t) - \phi(\alpha_n,\lambda,t) ] \nonumber \\
& + d [(n+1) \phi(\alpha_{n+1},\lambda,t) - n \phi(\alpha_n,\lambda,t) ] \nonumber \\
&- h n \phi(\alpha_n,\lambda,t) + \lambda f \phi(\beta_n,\lambda,t),
\label{eq_cs_1_exact}\\
\frac{d \phi(\beta_n,\lambda,t)}{dt} =
& d [(n+1)\phi(\beta_{n+1},\lambda,t) - n \phi(\beta_n,\lambda,t)]  \nonumber \\
&+  h n \phi(\alpha_n,\lambda,t) - f \phi(\beta_n,\lambda,t).
\label{eq_cs_2_exact}
\end{align}
Although Eqs.~\eqref{eq_cs_1_exact} and \eqref{eq_cs_2_exact}
are similar to Eqs.~\eqref{eq_master_monomer_1_exact} and \eqref{eq_master_monomer_2_exact},
note that the final term in the right hand side of Eq.~\eqref{eq_cs_1_exact} has
a factor $\lambda$.
The factor $\lambda$ is introduced in order to count the number of transitions,
and we can count the number of transitions related to this term (for details, see Appendix).
Using the above restricted generating functions,
the generating function $F(\lambda,t)$ is calculated as
\begin{align}
F(\lambda,t) = \sum_{n=0}^{\infty}
\left\{ \phi(\alpha_n,\lambda,t) + \phi(\beta_n,\lambda,t) \right\}.
\label{eq_total_generating_function}
\end{align}

Next, we introduce the following generating functions
for $\phi(\alpha_n,\lambda,t)$ and $\phi(\beta_n,\lambda,t)$:
\begin{align}
&\alpha(\lambda,z,t) \equiv \sum_{n=0}^{\infty} \phi(\alpha_n,\lambda,t) z^n,\\
&\beta(\lambda,z,t) \equiv \sum_{n=0}^{\infty} \phi(\beta_n,\lambda,t) z^n.
\end{align}
It is straightforward to derive the time-evolution equations
for the new generating functions $\alpha(\lambda,z,t)$ and $\beta(\lambda,z,t)$
from Eqs.~\eqref{eq_cs_1_exact} and ~\eqref{eq_cs_2_exact};
\begin{align}
\frac{d \alpha(\lambda,z,t)}{dt} =
& (z-1) \left[ g \alpha(\lambda,z,t)
- d \frac{\partial \alpha(\lambda,z,t)}{\partial z} \right] \nonumber
\label{eq_cs_1_exact_modified} \\
& - hz \frac{\partial \alpha(\lambda,z,t)}{\partial z} + \lambda f \beta(\lambda,z,t), \\
\frac{d \beta(\lambda,z,t)}{dt} =
& - (z-1) d \frac{\partial \beta(\lambda,z,t)}{\partial z} \nonumber \\
& +  hz \frac{\partial \alpha(\lambda,z,t)}{\partial z}
- f \beta(\lambda,z,t). \label{eq_cs_2_exact_modified}
\end{align}
Using the generating function $\alpha(\lambda,z,t)$ and $\beta(\lambda,z,t)$,
 the generating function $F(\lambda,t)$ is given by
\begin{align}
F(\lambda,t) = \alpha(\lambda,z=1,t) + \beta(\lambda,z=1,t),
\end{align}
and therefore it is enough to solve the following time-evolution equations
in order to calculate the generating function $F(\lambda,t)$:
\begin{align}
\frac{d \alpha(\lambda,t)}{dt} =&
- h \left. \frac{\partial \alpha(\lambda,z,t)}{\partial z} \right|_{z=1}
+ \lambda f \beta(\lambda,t), \label{eq_cs_1_reduced_exact}\\
\frac{d \beta(\lambda,t)}{dt} =&
 h \left. \frac{\partial \alpha(\lambda,z,t)}{\partial z} \right|_{z=1}
- f \beta(\lambda,t), \label{eq_cs_2_reduced_exact}
\end{align}
where we define $\alpha(\lambda,t) \equiv \alpha(\lambda,z=1,t)$
and $\beta(\lambda,t) \equiv \beta(\lambda,z=1,t)$.

Note that Eqs.~\eqref{eq_cs_1_reduced_exact} and \eqref{eq_cs_2_reduced_exact}
contain the derivative of $\alpha(\lambda,z,t)$ with respect to $z$.
Hence, the equations are not closed.
If these terms are expressed simply using $\alpha(\lambda,t)$,
we will have simultaneous differential equations written only by
the generating functions $\alpha(\lambda,t)$ and $\beta(\lambda,t)$;
i.e., we have closed equations and hence
the obtained equations may be solved analytically.
In the following analysis,
an effective interaction approximation is employed,
and we will show that the above statistics
can be approximated by a simple two-state model.

\subsection{Approximation for the interaction}
\label{subsec_approximation}

In the effective interaction approximation,
the interaction function $\mathcal{H}(n)$ is replaced as a constant value.
As shown in \cite{Ohkubo2011}, the dependence of $\mathcal{H}(n)$ on $n$
makes it difficult to obtain analytical results,
and it has been shown that the approximation gives qualitatively good results.

Replacing the interaction function $\mathcal{H}(n)$ as
\begin{align}
\mathcal{H}(n) = \tilde{h},
\end{align}
where $\tilde{h}$ is a constant,
we obtain the following equations instead of
Eqs.~\eqref{eq_cs_1_reduced_exact} and \eqref{eq_cs_2_reduced_exact}:
\begin{align}
\frac{d \alpha(\lambda,t)}{dt} =&
- \tilde{h} \alpha(\lambda,t) + \lambda f \beta(\lambda,t),
\label{eq_cs_1_reduced_approximate}\\
\frac{d \beta(\lambda,t)}{dt} =&
 \tilde{h} \alpha(\lambda,t) - f \beta(\lambda,t).
\label{eq_cs_2_reduced_approximate}
\end{align}
Note that Eqs.~\eqref{eq_cs_1_reduced_approximate} and \eqref{eq_cs_2_reduced_approximate}
are written only by $\alpha(\lambda,t)$ and $\beta(\lambda,t)$.
It means that the switching problem can be approximated as
a simple two-state model
\textit{if} the effective interaction $\tilde{h}$ is chosen adequately.

We here briefly explain the choice of the effective interaction $\tilde{h}$
using a simple example, i.e., the monomer binding interaction case.
For the monomer binding interaction,
the interaction function is calculated as follows \cite{Ohkubo2011}.
In this case, the interaction function is $hn$.
In order to obtain the effective interaction $\tilde{h}$,
the number of proteins $n$ is replaced
as the average number of proteins, i.e.,
\begin{align}
\tilde{h} = h \langle n \rangle_{\alpha},
\label{eq_effective_int_monomer}
\end{align}
where $\langle n \rangle_{\alpha}$ is the expectation of the number of free regulatory proteins
under a condition that the gene is in the active state (conditional expectation).

The conditional expectation can be calculated from
the stationary distribution of the number of proteins.
Note that the generating functions
$\alpha(\lambda,z,t)$ and $\beta(\lambda,z,t)$
are reduced to generating functions for the stationary distribution
of the number of proteins when $\lambda = 1$.
Hence, as shown in \cite{Ohkubo2011},
they are written as follows.
\begin{align}
\alpha(z) \equiv& \lim_{t\to \infty}\alpha(\lambda=1,z,t) = A F[a,b,N(z-1)], \\
\beta(z) \equiv& \lim_{t\to \infty}\beta(\lambda=1,z,t) \nonumber \\
= &
\left( 1+ \frac{\tilde{h}}{f} \right) A F[a-1,b-1,N(z-1)] - \alpha(z)
\end{align}
where $A = f/(f+\tilde{h})$ and
\begin{align*}
N = \frac{g}{d}, \quad a = 1 + \frac{f}{d}, \quad b = 1 + \frac{f+\tilde{h}}{d}.
\end{align*}
$F(p,q,r)$ is the Kummer confluent hypergeometric function,
\begin{align}
F(p,q,r) \equiv \sum_{n=0}^{\infty} \frac{(p)_n}{(q)_n} \frac{r^n}{n!},
\end{align}
where $(p)_n = p(p+1)(p+2)\cdots (p+n-1)$.
We, therefore, obtain
\begin{align}
\langle n \rangle_{\alpha}
\equiv \frac{1}{\alpha(1)}
\left. \frac{\partial}{\partial z} \alpha(z) \right|_{z=1}
= \frac{g(d+f)}{d(d+f+\tilde{h})}.
\label{eq_for_effective}
\end{align}
By inserting Eq.~\eqref{eq_for_effective} into Eq.~\eqref{eq_effective_int_monomer},
the following self-consistent equation is derived:
\begin{align}
\tilde{h}
= h \frac{g(d+f)}{d(d+f+ \tilde{h})}.
\label{eq_monomer_self_consistent}
\end{align}
Solving Eq.~\eqref{eq_monomer_self_consistent},
we obtain
\begin{align}
\tilde{h} =
\frac{-(d^2 + fd) + \sqrt{(d^2+fd)^2 + 3 h g d(d+f)}}{2d}.
\end{align}

We finally comment on a solution of the simple two-state model
(Eqs.~\eqref{eq_cs_1_reduced_approximate} and \eqref{eq_cs_2_reduced_approximate}).
The simple two-state model can be solved exactly \cite{Gopich2003,Gopich2006},
and the probability distribution $P(k|t)$ for the number of ``inactive $\to$ active'' transitions
during time $t$ is explicitly written as follows:
\begin{align}
P(k|t) =& \left( \frac{(1-\gamma^2) T}{2\gamma}\right)^k \frac{e^{-T}}{k! \sqrt{8 \gamma T / \pi}} \nonumber\\
&\times \left\{ 2 \gamma (k+T) I_{k-1/2}(\gamma T)
 + (1+\gamma^2) T I_{k+1/2}(\gamma T )  \right\},
\label{eq_prob_distribution}
\end{align}
where $T = (f + \tilde{h}) t /2$, $\gamma^2 = 1 - 4 f \beta(1)/(f+\tilde{h})$,
and $I_n(z)$ are modified Bessel functions of the first kind.
This expression~\eqref{eq_prob_distribution} immediately gives us
the non-Poissonian picture of the phenomenon.

\section{Numerical results}
\label{sec_results}

\begin{figure}
\begin{center}
\includegraphics[width=70mm]{fig2a.eps}\\
\includegraphics[width=70mm]{fig2b.eps}
\end{center}
\caption{Probability distributions for the number of ``inactive $\to$ active'' transitions.
(a) Monomer binding interaction case. (b) Dimer binding interaction case.
In each figure, filled circles and filled boxes are Monte Carlo results
for time $t = 10$ and $t = 100$, respectively.
Solid and dashed lines corresponds to
approximated analytical results of Eq.~\eqref{eq_prob_distribution}
for time $t = 10$ and $t = 100$, respectively.
}
\label{fig_results}
\end{figure}

In order to check the validity of the analytical treatments
and the approximations,
we here compare the analytical results with
those of Monte Carlo simulations.
The original genetic switch explained in Sec.~\ref{sec_model}
was simulated using a standard Gillespie algorithm \cite{Gillespie1977}.
The parameters used in the simulation are as follows:
$d = 1, \quad g = 50.0, \quad h = 0.004, \quad f = 0.1$.
Note that these parameters were selected as one of the typical values used
in the previous works \cite{Schultz2007,Ohkubo2011}.

Firstly, we consider the monomer binding interaction case.
According to the discussions in Sec.~3.3,
the value of the effective interaction $\tilde{h}$ is calculated as $\tilde{h} = 0.173$.
Figure~\ref{fig_results}(a) shows
the results of the analytical calculations (Eq.~\eqref{eq_prob_distribution})
and those of the Monte Carlo simulations.
Although there are quantitative differences,
the results shows that the approximated two-state model
captures the essential features of the phenomenon.

Next, we consider a dimer binding interaction case,
i.e., $\mathcal{H}(n) = h n(n-1)/2$.
In this case, the effective interaction is calculated as follows:
\begin{align}
\tilde{h} = h \frac{\langle n (n-1) \rangle_{\alpha}}{2}.
\end{align}
As shown in \cite{Ohkubo2011},
the effective interaction $\tilde{h}$ is obtained
by solving the following self-consistent equation:
\begin{align}
\tilde{h}
= \frac{h}{2} \frac{1}{\alpha(1)}
\left. \frac{\partial^2}{\partial z^2} \alpha(z) \right|_{z=0}.
\label{eq_dimer_self_consistent}
\end{align}
We here numerically solved the self-consistent equation
(Eq.~\eqref{eq_dimer_self_consistent}),
and the calculated value of the effective interaction is $\tilde{h} = 1.358$.
Using the calculated value,
we depict the analytical results and the corresponding Monte Carlo results
in Fig.~\ref{fig_results}(b).
From the comparison,
we confirmed that
the approximated two-state model is available even in the dimer binding interaction case.
Although results are not shown,
we performed numerical simulations for other some parameters,
and checked the validity of the analytical treatments.
For example, even for parameter regions in which
the probability distribution of the number of proteins has bistability,
the approximation scheme works well.

\section{Conclusions}
\label{sec_conclusions}

In the present paper, we studied an analytical scheme
to extract information related to the dynamical behavior in genetic switches.
Using an effective interaction approximation,
a simple two-state model is obtained,
and we confirmed that the two-state model captures the features of the phenomenon.
Note that in the analytical treatments,
we did not neglect the stochastic properties of the system
(except for the effective interaction approximation);
i.e., we can calculate all statistics for transitions approximately, including higher order moments.
It could be possible
to apply the above effective expression for the transitions between the active and inactive states
to more complicated gene regulatory networks
without loss of the stochasticity;
this would give us deeper understanding
for the switching behavior of the gene regulatory systems
including static, dynamical, and stochastic behaviors.
In addition, the idea of the effective interaction may be similar
to the mean-field approximation in statistical physics;
the interaction is replaced with the average.
It may be possible to develop higher-order approximations
using the analogy with the conventional approximation schemes in statistical physics;
this is an important future work.

We discussed properties only in the stationary states,
because the effective interaction approximation
has been applied only for the stationary states at the moment;
the average number of proteins (or higher moments) should be estimated adequately,
and it was calculated by using the analytical solutions
for the \textit{stationary} distributions of the number of proteins.
Recently, exact time-dependent solutions for a self-regulating gene have been derived
\cite{Ramos2011}.
Hence, it may be possible to extend the effective interaction approximation
to non-stationary states.
If so, the effective interaction $\tilde{h}$ would be time-dependent,
and, at least numerically,
it is possible to calculate various moments for the counting statistics
for time-dependent systems \cite{Ohkubo2010}.
We expect that the simple description developed in the present paper
is available for various cases,
such as complicated regulatory systems and time-dependent systems,
and that the description gives new insights
for the regulation mechanisms and stochastic behaviors.

\section*{ACKNOWLEDGMENTS}

This work was supported in part by grant-in-aid for scientific research
(Nos. 20115009 and 21740283)
from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan.


\section{Introduction}

Since the start of the first Large Hadron Collider (LHC) physics run in 2009,
the Compact Muon Solenoid (CMS) experiment~\cite{cmsdetpaper} has been
collecting data successfully. During Phase 1 of the LHC which will continue
until
about 2022, the instantaneous
luminosity delivered to the experiments is foreseen to increase gradually up
to more than twice its nominal value of $10^{34}$~cms$^{-2}$s$^{-1}$. At
present, two
shutdown periods are scheduled to give the machine and the experiments
the necessary time to
anticipate these luminosity increases: Long Shutdown 1 (LS1) in 2013/2014
should prepare the accelerator to run at its nominal luminosity and
Long Shutdown 2 (LS2) in
2018/2019 should take it to 2.2 times this value.

\begin{figure}[h]
\begin{center}
\includegraphics[width=.65\textwidth]{muonsystem.eps}\hspace*{.03\textwidth}
\includegraphics[width=.3\textwidth]{photoendcap.eps}
\caption{\label{fig:cmsmuonsystem}(left) One quarter of the CMS muon
  system. The colored detectors are the RPCs in the endcap system. RPCs
  enclosed in boxes have yet to be installed. (right) Photograph of endcap RPC
  detectors taken during installation of the CMS detector.}
\end{center}
\end{figure}

During these long shutdown periods the CMS Collaboration
intends to upgrade several subsystems of its
detector~\cite{cmsupgradeTP}. In particular, the instrumentation of
the muon system as depicted
in Fig.~\ref{fig:cmsmuonsystem} will be
extended in both endcaps
to ensure efficient muon triggering and reconstruction in that region
at high luminosities. In the endcaps, CMS is using Cathode Strip Chambers
(CSCs) as muon
tracking and trigger detectors, while Resistive Plate Chambers (RPCs) serve
as dedicated trigger detectors and improve the muon reconstruction.
At present, the 4th endcap disks remain largely uninstrumented:
CSCs are only installed in the very forward region
($1.8<|\eta|<2.4$), and
RPCs are missing completely and cover only the first three endcap disks
up to $|\eta|=1.6$.

During LS1
these 4th endcap stations will be instrumented further with
new CSCs below $|\eta|=1.8$ and new RPCs up to $|\eta|=1.6$.
In the following, we
describe the CMS plans for this first step in the upgrade of the RPC endcap
system.

For the latter subsystem the very forward region beyond
$|\eta|=1.6$ will still remain empty after LS1
and could in principle in a second step
be instrumented
up to $|\eta|=2.4$ matching the CSC system. However, the present
design of the endcap RPCs, made of a double High Pressure Laminate (HPL)
gas gap and operating in
avalanche mode, is not expected to be suitable for the particle rates
amounting to several tens of kHz/cm$^2$ in the scenario of an LHC luminosity
going up to $10^{34-35}$~cm$^{-2}$s$^{-1}$. Dedicated R\&D activities
to identify suitable
technologies to instrument that particular endcap region during LS2 is
ongoing, see for instance~\cite{ksleerpc2012, ieee2011paper}.

\section{Motivation for the RPC Upgrade}

Muons with high transverse momenta are one of the key objects in the
detection of possible new physics phenomena.
By 2015, the LHC luminosity should reach $10^{34}$~cm$^{-2}$s$^{-1}$. The
in-time pileup will be right at the edge of the CMS design envelope and will
present special challenges for the muon system to trigger on
high-$p_T$ muons.
The RPC upgrade is essentially driven by the impact of
the instantaneous peak luminosity on the trigger system.
In the endcaps, the RPC system
provides excellent timing with a somewhat worse momentum resolution compared
to the CSC system. To be
effective, the muon trigger
must achieve good enough resolution to identify high-$p_T$ tracks.
With the RPC trigger
requiring segments in at least three stations, the endcap system does not have
the necessary redundancy to control the
trigger rate at the increased luminosity while preserving high trigger
efficiency. The problem stems from mis-measurements of low-$p_T$ muons, that
are promoted to high-$p_T$ muons and contribute to the trigger rate.
With the much higher flux of
low-momentum muons at increased luminosity, such poorly measured muons
will dominate the trigger rate and make it unacceptably high, leading in turn
to unacceptably high muon-$p_T$ trigger thresholds.

\begin{figure}[h]
\begin{center}
\includegraphics[width=.6\textwidth]{RPCL1TriggerEff.eps}
\caption{\label{fig:RPCL1TriggerEff}Simulated RPC Level-1 trigger efficiency
  for the present system with three endcap layers compared to the upgraded
  situation with four endcap layers.}
\end{center}
\end{figure}

In Fig.~\ref{fig:RPCL1TriggerEff} the difference in the RPC Level-1 Trigger
performance is simulated between the present system with three stations and the
four station situation after the upgrade.
The present RPC trigger logic requires hits in
at least three layers, which causes the observed drop in efficiency for the
endcaps with only three stations. Adding the
4th layer in the endcaps, enabling a 3-out-of-4 trigger logic in
those regions, will bring the RPC endcap performance to a similar
level as the barrel system.

\section{Design of the New CMS Endcap RPCs}

The layout of the CMS RPC endcap stations and chambers is depicted in
Fig.~\ref{fig:layouts}. Each station
consists of three concentric rings,
called REx/1-3 (station x=1,2,3),
with chambers mounted in a staggered way. For the
instrumentation of the
4th station up to $|\eta|=1.6$, 144 new RPCs are required
to be mounted in two such concentric rings (RE4/2 and RE4/3) per
station, with 36 chambers per ring. As the design of the new chambers is
quite similar to the one of the existing RPCs, a total of 200 new chambers
will be built with 56 spare chambers for the RE2-4 stations.

\begin{figure}[h]
\begin{center}
\includegraphics[width=.4\textwidth]{presentRE2.eps}\hspace*{.03\textwidth}
\includegraphics[width=.55\textwidth]{ChamberLayout.eps}
\caption{\label{fig:layouts} (left) The second RPC endcap
  station during assembly of the CMS detector. (right)
  Schematic layout of an endcap chamber.}
\end{center}
\end{figure}

The endcap RPCs contain a double HPL gas gap with a
copper strip readout
panel placed in between. In fact, the design uses one large gas gap at the
bottom and two smaller gaps on top of the chamber.
The chamber is embedded in a honeycomb
box with the chamber services (readout electronics, gas in and outlets and
water cooling circuit) mounted on the outside.

The readout strip panel is divided in 3 $\eta$-partitions, with 32 strips each,
yielding a total of 96 strips per chambers. Strips are connected through coaxial
cables to Adaptor Boards linked to 3 Front End Boards (FEBs) per
chamber. Every Chamber has 1 Distribution Board for the electronics
control. The off-detector electronics consist mainly of Link Boards (LBs) that
receive signals in LVDS (Low-Voltage Differential Signaling) standard
from the FEBs and perform the synchronization with
the LHC clock and the transmission to the Trigger Electronics in the control
room. Each LB crate contains a Control Board (CB) that drives the crate,
provides inter-crate communication and takes care of the connection to the
readout and trigger systems. The LBs and CBs will include minor design
improvements related to low voltage regulation and protection circuits, to
cure problems observed in the past with the present system.

The working point of the present system depends strongly on the
environmental temperature, and the latter should therefore be kept inside
21-23$^{\circ}$C. To monitor the temperature,
each endcap chamber will be equipped with one Fibre Bragg
Grating
sensor mounted inside a heat conducting housing;
these sensors are radiation hard,
insensitive to magnetic fields, have a high precision, produce no
electrical noise, are easy to install,
require a minimal amount of cabling and finally cost about the same as
conventional
sensors. The latter cannot be used in this case due space limitations for
cabling.
As the new RE4 RPCs will be mounted inside the CMS detector
facing the CSC electronics, there is some danger of additional external
heating to the RPC gas gaps. The water cooling circuit on the chambers
consists of Cu pipes
soldered onto three copper plates where the FEBs are mounted. An aluminum
screen box is covering the cooling circuit and the FEBs.
To better stabilize the operating temperature, the layout and size of the
copper plates was optimized using
FloEFD/pro\footnote{FloEFD$^{\mbox{\tiny TM}}$ software package
by Mentor Graphics}
and the final configuration was
successfully tested in the laboratory with a detector inside a heated box. In
addition to that, every 10$^{\circ}$ sector of the endcap disk will get a
connection to the
cooling water manifolds in the cavern, i.e. only
2 chambers per connection.

\begin{figure}[h]
\begin{center}
\includegraphics[width=.47\textwidth]{RE4Mounting.eps}\hspace*{.03\textwidth}
\includegraphics[width=.47\textwidth]{RE4SuperModule.eps}
\caption{\label{fig:RE4Installation}
(left) Drawing of the new RPC modules installed
  on the CSC mounting posts on the backside of the 3rd endcap yoke. (right)
  Photograph of an RE4 Super Module test installation at CERN.}
\end{center}
\end{figure}

As shown in Fig.~\ref{fig:RE4Installation}, the new RE4 RPCs will be attached
on CSC mounting posts
on the back of the 3rd endcap yoke of the CMS detector. Before installation,
the corresponding
RE4/2 and RE4/3 chambers will be pre-mounted together
on Super Modules that cover a
10$^{\circ}$ sector of an endcap. This will not only reduce the amount of
cabling work to be done in the CMS cavern, but will also speed up the
installation of the chambers.

\section{RPC Production and Quality Control}

The fabrication of the RPC gas gaps starts with the
production of HPL sheets. The raw material is coming from the Puricelli
firm (Milano), where about 600 sheets of $1620\times 3200$~mm$^2$ are being
produced. The quality control (QC) of the HPL sheets is handled by
INFN Pavia. Next, the sheets are cut by Riva (Milano) and finally surface
cleaned by General Tecnica (Frosinone).

\begin{figure}[h]
\begin{center}
\includegraphics[width=.6\textwidth]{bakeliteresist_new.eps}
\caption{\label{fig:bakelite}Measured resistivity for recent production
  batches of HPL sheets.}
\end{center}
\end{figure}

CMS has put forward detailed specifications for the HPL sheets; main
requirements are a high uniformity in resistivity across the surface of
the sheets and also between the sheets, with a resistivity value between
$1-6\cdot10^{10}~\Omega$cm (measured at $20^{\circ}$C).
Fig.~\ref{fig:bakelite} displays the resistivity values measured for two recent
production batches. The data represent the average resistivity value and spread
measured at nine different positions across the sheet surface; 50 panels are
produced per batch and the obtained
resistivity values are clearly uniform and well
inside the CMS specifications.

\begin{figure}[h]
\begin{center}
\includegraphics[width=.4\textwidth]{gasgap.eps}
\caption{\label{fig:gasgap}Photograph of a completed HPL gas gap at Kodel.}
\end{center}
\end{figure}

The HPL sheets that pass the QC are sent to the KOrea DEtector Laboratory
(KODEL) at Korea University, where
about 700 gas gaps will be produced~\cite{sparkrpc2012}. The HPLs receive a
graphite coating, protected with a polyethylene terephthalate (PET) film.
Assembled gaps are treated with
linseed oil mixed with heptane.
Again very detailed technical
specifications and QC protocols were prepared by
CMS to ensure high quality gaps. Fig.~\ref{fig:gasgap} shows an example of a
HPL gas gap assembled at
KODEL.

The chamber mechanics (honeycomb boxes, screen boxes, read-out strip
panels, etc.) were produced by
two Chinese companies: Beijing Axicomb Technology Co., Ltd and Beijing
Gaonengkedi SGT Co., Ltd. In the past, those companies have already produced
the same components for the existing system.

The detector electronics are procured in Pakistan. For this project a total of
1350 Adaptor Boards, 650 Front End Boards and 250 Distribution Boards will be
produced. The new off-detector electronics will be produced by the INFN, while
the CMS Warsaw (Poland) group will take on the responsibility to integrate the
new electronics in the trigger system.

The final detector assembly is shared by three different institutions:
Barc (Mumbai, India), Ghent University (Ghent, Belgium) and CERN. Each of
these assembly sites has built up a new RPC lab, including setups for detailed
quality control of the HPL gas gaps and completed chambers.

A common quality control protocol for the chamber production has been carefully
prepared for each level of the construction. Individual components need to
satisfy detailed mechanical specifications. The HPL gas gaps will undergo visual
inspections, leak tests and dark current measurements. Assembled chambers will
be thoroughly tested (efficiency, cluster size, long term current monitoring,
noise) in cosmic test benches at each of the assembly sites. The
tests at chamber level will be repeated once the chambers have been
transported to CERN and yet again when the chambers have been mounted in the
Super Modules. The final checkout will occur after the installation of
the Super Modules in the CMS detector.

The full history of each chamber and every component is stored in an
Oracle based Construction Database. The database includes all measurements
performed before, during and after the detector assembly. This should enable
CMS to follow the evolution of each chamber in time, especially in case
problems should appear during operation later on.

\section{Summary}

As the instantaneous luminosity of the LHC continues to increase over time
during Phase 1,
the CMS experiment has to ensure its detector can operate in a stable and
efficient manner for increasingly high particle rates. In particular, the
muon system should be able to keep its trigger rate and efficiency under
control at a low enough muon $p_T$-threshold for physics studies.
To this end, CMS will among other things
extend its RPC system with a 4th layer in both endcaps.
The construction of 200 new RPCs of the standard CMS type
recently started. Components for the new chambers are being produced at
various places in the world. The assembly and quality control of the
chambers will be handled at three
different locations in India, Belgium and CERN.
The new chambers are to be installed inside the CMS detector
during the first LHC Long Shutdown in 2013-2014.

\acknowledgments

We acknowledge support from FWO (Belgium); MEYS (Bulgaria); CERN; INFN (Italy); NRF and WCU (Korea); Swiss Funding Agencies (Switzerland).


\section{Next-to-leading order correction to the decay amplitude}

In order to establish when the coupling constant $g$ can be regarded as small
such that the resummed one-loop approximation is justified, we calculate the
(upper limit of) the next-to-leading order diagram of the amplitude of the
decay process, see Fig. \ref{diagrams}. The decay amplitude results as the sum
of two terms
\begin{equation}
-i\mathcal{M}=-i\mathcal{M}_{1}-i\mathcal{M}_{2}\text{ ,
\end{equation}
where the tree-level and the triangle contributions depicted in Fig.
\ref{diagrams} read
\begin{equation}
-i\mathcal{M}_{1}=2ig\text{ ; }-i\mathcal{M}_{2}=4g^{3}I\text{ .}
\label{amplitudes
\end{equation}
The coefficients $2$ and $4$ in (\ref{amplitudes}) are the appropriate
symmetry factors. The quantity $I$ is the integra
\begin{equation}
I=I(m,M_{0})=-\int\frac{d^{4}q}{(2\pi)^{4}}\frac{\phi_{1}\phi_{2}\phi_{3
}{\left[  \left(  \frac{p}{2}-q\right)  ^{2}-m^{2}\right]  \left[  \left(
\frac{p}{2}+q\right)  ^{2}-m^{2}\right]  \left[  \left(  \frac{p}{2
+q-k_{1}\right)  ^{2}-M_{0}^{2}\right]  }\text{ .} \ ,\label{integral
\end{equation}
where $\phi_{i}$ is the cutoff function of the $i$-th vertex

\begin{figure}
[ptb]
\begin{center}
\includegraphics[
height=2.1819in,
width=5.3065in

{diagrams.eps
\caption{Tree-leve (left) and next-to-leading order (traingle diagram, right)
contributions to the decay amplitude.
\label{diagrams
\end{center}
\end{figure}

We have evaluated $I$ by using the free propagator of the $S$ virtual particle
exchanged by the two $\varphi$ fields. Thus, our results represent an upper
limit for $\left\vert I\right\vert $. Namely, $\left\vert I\right\vert $ would
be actually smaller by including the full propagator of $S$, i.e. taking into
account its finite width.

We define the coupling constant $g_{\max}$ as the value of $g$ for which
$\left\vert \mathcal{M}_{1}\right\vert =\left\vert \mathcal{M}_{2}\right\vert
.$ Namely, for $g=$ $g_{\max}=1/\sqrt{2\left\vert I\right\vert }$ the
next-to-leading order contribution equals the tree-level one. Thus, for
$g>g_{\max}$ one cannot neglect the role of higher order diagrams.

In particular, we obtain the following numerical results: For $M_{0}=1$ GeV,
$g_{\max}=19.503$ GeV; the critical value $g_{\ast},$ at which the new pole in
the I-Riemann sheet emerges, reads $3.379$ and the value $g_{E},$ for which
the stable state has a mass of $38$ MeV, reads $4.877$ GeV. Both values are
safely smaller than $g_{\max}$: the amplitude ratios read $\left[  \left\vert
\mathcal{M}_{2}\right\vert /\left\vert \mathcal{M}_{1}\right\vert \right]
_{g=g_{\ast}}=0.03$ and $\left[  \left\vert \mathcal{M}_{2}\right\vert
/\left\vert \mathcal{M}_{1}\right\vert \right]  _{g=g_{\ast}}=0.06,$ thus the
tree-level diagram is by far the dominant contribution for all the range of
coupling constants relevant for our study. For $M_{0}=0.45$ GeV one gets
$g_{\max}=3.666$ GeV, which should be compared with $g_{\ast}=1.399$ GeV and
$g_{E}=2.189$ GeV. The ratios read $\left[  \left\vert \mathcal{M
_{2}\right\vert /\left\vert \mathcal{M}_{1}\right\vert \right]  _{g=g_{\ast
}=0.14$ and $\left[  \left\vert \mathcal{M}_{2}\right\vert /\left\vert
\mathcal{M}_{1}\right\vert \right]  _{g=g_{\ast}}=0.36,$ which are still
safely small, although larger than in the previous example. It should once
more be stressed that our evaluation of $\left\vert I\right\vert $ represents
only an upper limit, and the value of $g_{\max}$ is indeed larger than the
estimate given here.

All these considerations enforce our point of view that the emergence of a
pole below threshold such as the one described in this manuscript is possible
in quantum field theory. Future evaluation of higher-loop contributions is planned.

\bigskip


\section{Introduction}
Industrial embedded systems are hardware systems that have associated software to operate them and which perform specific functionality. Such systems are often highly concurrent with several parallel components operating together \cite{AS11}. These components have synchronizations and dependencies between them that make such systems highly complex. The formal approaches are used to design models of completely new systems as well as to evaluate the existing systems. The correctness of these models becomes very crucial since they form input for the design of the actual software that controls the operations of the system hardware. Simulation is one approach that allows the analysis of certain behavioral aspects of the system. However, it is also essential to ensure that the system behaves accurately under all possible circumstances.
\par Embedded systems with hard real time constraints are very time-specific and require the verification of time-specific logistic rules. In such cases, model checking approaches that also allow time based verification needs to be considered. The main purpose of this work is to perform time based model checking on such an industrial system.

\par However the accurate time-based verification requires the translation of the current design models of the system to formal semantics of existing modeling languages and analysis of the entire state space of the system by means of state-of-the-art model checkers. Several model checkers are suitable for this kind of verification such as UPPAAL \cite{GDK04}\cite{JKF95}, mCRL2 \cite{JAMY06} and other temporal logic model checkers \cite{WIKI12}. These model checkers can be used to model highly concurrent systems and perform verification of required properties. In this paper we have used UPPAAL as  model checking tool which allows  verification of the required properties on the system for all its dynamic behaviors. This paper is the extended version of our earlier work \cite{REF13}.
\par
The remaining paper is organized as follows. Section \ref{1} explains the concept of timed automata. Section \ref{2} briefly describe about the various modeling approaches. Section \ref{3} explains various time models which affects the deadlock condition. Section \ref{4} describes the model checking tool UPPAAL. Section \ref{5} examines the limitation of the existing tiny two way buffer system. Section \ref{6} discusses about the proposed two way buffer system.  Section \ref{7} explains about the time analysis of proposed system. Section \ref{8} describes about the UPPAAL model of proposed buffer system in no time. Section \ref{9} describes about the UPPAAL model of proposed buffer system in deterministic time. Section \ref{10} describes about the UPPAAL model of proposed buffer system in non-deterministic time for no deadlock. Section \ref{11} discusses about buffer sizing of the proposed system which affects the efficiency of the overall system. Section \ref{12} describes the scope of the work and the future enhancement. Finally, section \ref{13} concludes this paper.

\section{Timed Automata}
\label{1}
A timed automaton (TA) is a finite-state machine extended with clock variables. It uses a dense-time model where a clock variable evaluates to a real number \cite{PUK12}. All the clocks progress synchronously. A timed automaton manipulate clocks, evolving continuously and synchronously with absolute time.\\
 The elements of TA include:

\begin{itemize}
	\item Each transition of such an automaton is labeled by a constraint over clock values (also called guard), which indicates when the transition can be fired
	\item A set of clocks to be reset when the transition is fired
	\item Each location is constrained by an invariant, which restricts the possible values of the clocks for being in the state, which can then enforce transition to be taken
	\item The time domain can be $\mathbb{N}$, the set of non-negative integers, or $\mathbb{Q}$, the set of non-negative rationals, or even $\mathbb{R}$, the set of non-negative real numbers
\end{itemize}

\textbf{Definition: Timed Automaton}: A timed automaton is a tuple (L, $l_0$, C, A, E, I) where,
\begin{itemize}
	\item L is a set of locations,
	\item $l_0 \in L$ is the initial location,
	\item $C$ is the set of clocks,
	\item $A$ is a set of actions, co-actions and the internal $\tau$-action,
	\item $E \subseteq L \times A \times B(C) \times 2^{C} \times L$ is a set of edges between locations with an action,
	\item $B(C)$ is a set of conjunctions over simple conditions of the form $x \bowtie c$ or $x - y \bowtie c$, where $x,y \in C$, $c \in \mathbb{N}$ and $\bowtie \in \{<, \geq, =, \leq, >\}$
	\item a guard and a set of clocks to be reset and
	\item $I : L \rightarrow B(C) $ assigns invariants locations.
\end{itemize}

\section{Various Modeling approaches}
\label{2}
These days there are many model checking tools used to model and adopt different aspects towards verification. SMV is the first model checker to use Binary Decision Diagrams (BDDs) \cite{REF10}, which are data structures used to represent boolean functions. mCRL2 \cite{REF8} is a model checker based on a specification language for describing concurrent discrete event systems and which uses mu-calculus \cite{REF9} for its specifications. Spin is used to perform verification on distributed systems and uses Linear Temporal Logic (LTL) to specify its properties \cite{REF11}. HyTech is designed for reasoning about temporal requirements in hybrid systems \cite{REF12}.
 \par Temporal logic model checkers are most reliable to perform model checking on the machines where the properties to be verified are time specific. The properties stating that an event can eventually happen or always happen after a given period of time is possible through temporal logic. Using temporal logic it is possible to represent time specific properties crucial for real time system we have designed.
\par There are various model checking tools that are used to perform the kind of verification for the two way buffer system which is a  real time system. Uppaal is the model checker for real time systems that uses the concept of timed automata for modeling purposes and a variant of Computation Tree Logic (CTL) for its specifications \cite{REF7}.

\section{Various Time Models}
\label{3}
The verification of the system requires the evaluation of various multiple time models as follows:
\begin{itemize}
	\item No time: It is a model of the system that does not consider the time taken by different actions and mainly considers the various possible sequences of interactions between the components. Verification on a system with no time ensures that all possible scenarios are considered and hence a system which satisfies a property with no time always satisfies the property.
	\item Deterministic time: In this model, specific time durations are assigned to all actions and the system is analyzed for this timing behavior. With deterministic time only a particular time sequence or behavior is considered. If a property is violated in the untimed version of the system, the system is not necessarily incorrect. This is because the property might be satisfied in the deterministic model of the system for all realistic time sequences which can occur in the actual machine.
	\item Non-deterministic time: This model allows actions to occur non deterministically within time intervals. This allows verification to be performed on a system with variable timings within given realistic ranges. Since the industrial systems have timing involved in them and the time duration of activities can vary, the non-deterministic models with time give the most realistic models of the industrial
machines. However, since the timing variations are within very small limits, the deterministic models can be efficiently used to analyze the behavior of the system.
\end{itemize}

\section{UPPAAL - Model Checking Tool}
\label{4}
In TA, notion of time is introduced by clock variables, which are used in clock constraints to model time-dependent behavior. Systems comprising multiple concurrent processes are modeled by networks of timed automata, which are executed with interleaving semantics and synchronization on channels. UPPAAL  is a tool set for the modeling, simulation, animation and verification of networks of timed automata. The UPPAAL model checker enables the verification of temporal properties, including safety, liveness and reachability properties. The simulator can be used to visualize counterexamples produced by the model checker in case the given specification is not satisfied by the model.
\par The UPPAAL modeling language extends timed automata by introducing bounded integer variables, binary and broadcast channels, and urgent and committed location. Timed automata are modeled as a set of locations, connected by edges. The initial location is denoted by $\circledcirc$. Invariants can be assigned to locations which enforce synchronization between them. Edges may be labeled with \textit{guards}, \textit{synchronizations} and \textit{updates}. Updates are used to reset clocks and to manipulate the data space. Processes synchronize by sending and receiving events through channels. Sending and receiving via a channel 'c' is denoted by c! and c?, respectively. Binary channels are used to synchronize one sender with a single receiver. A synchronization pair is chosen non-deterministically if more than one channel is enabled. Broadcast channels are used to synchronize one sender with an arbitrary number of receivers. Any receiver that can synchronize must do so. Urgent and committed locations are used to model locations where no time may pass. Committed locations enforce that the synchronization is atomic and are denoted by the symbol $\copyright$. Leaving a committed location has priority over leaving non-committed locations.
\par A UPPAAL model comprises three parts: \textit{global declarations}, \textit{parameterized timed automata} and a \textit{system declaration}. In the global declarations section, global variables, constants, channels and clocks are declared. In the system declaration, TA templates are instantiated and the system to be composed is given as a list of TA.

\section{Existing Tiny Two Way Buffer System}
\label{5}
The existing system consists of two generators, two buffers, and two exits as shown in Figure 1 \cite{AS11}. The first generator G1 repeatedly produces packets. This packet is received over channel 'a' by buffer M1 which sends it over channel 'c2'. The packet is received over this channel by a second buffer M2 which again sends it over channel 'e' to be received by the exit E2. Similarly, the second generator G2 repeatedly generates packets and sends them over channel 'b'. This is received by buffer M2 which sends it over channel 'c1' to buffer M1. M1 again forwards the packet over channel 'd' to exit E1. This system, obviously, has a deadlock which occurs when both G1 and G2 produce packets at the same instant of time.
\par In that case, the packet is received by the buffers M1 and M2 respectively at the same time. As a result, both M1 and M2 wait to send packets to each other and the system deadlocks.
\par Verification of the system has been done for various time models as for no time, deterministic time and non deterministic time. and it is giving significance results as for no time, deterministic time and nondeterministic time deadlock property not satisfied in the Uppaal user interface ie deadlock occurs with no time, deterministic time and nondeterministic time.
\par In order to study a system with no deadlock, a slightly modified version of the tiny two way buffer system was considered \cite{AS11}. In this system, both generator G1 and G2 produce a packet after a delay of 10 seconds. However, G1 has an initial shift of 5 seconds after which it begins sending a packet every 10 seconds.
Hence, initial delay for G1 is 15 seconds and for G2 it is 10 seconds. Due to this initial shift of 15 seconds at G1 and 10 seconds delay at G2 the time consumption in waiting to maintain synchronization in transfer of packets is more.
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.25]{fig1}
	\caption{Existing Model of Two Way Buffer System}
	\label{Fig1}
\end{figure}

\section{Proposed Two Way Buffer System}
\label{6}
In this paper we have proposed a two way buffer system which consists of two generators, three buffers and two exits. The generators generate the packets to be transferred until it reaches the exits. Figure 2 depicts the proposed model.

\begin{figure}[bpht!]
\centering
	\includegraphics[scale=1.4]{fig2}
	\caption{Proposed Model of Two Way Buffer System}
	\label{Fig2}
\end{figure}

First, the generator G1 repeatedly produces packets. This packet is received over the channel 'a' by buffer M1. The packet is received by buffer M2 over channel 'e' and sends it to exit E1 via the channel 'c'. Similarly, the second generator G2 repeatedly generates packets and sends them over channel 'b'. This is received by buffer M2 which sends it over an intermediate dummy buffer Md through channel 'g'. The dummy buffer sends the packet to buffer M1 via channel 'f'. M1 again forwards the packet over channel 'd' to exit E1. This process is repeated infinitely until it is stopped. We are incorporating the delay $(\theta)$ in transmission of packet at buffer M1 and waiting for M2 to be free. After waiting, buffer M1 will release the packet and send it to M2. On the other direction the packet released from M2 will reach at buffer Md, now Md will wait until M1 is free. Then Md will release the packet and send it M1. Based on the clocks being synchronized the data is transferred concurrently. The two way data packet transfer through the proposed system is concurrent. The deadlock is avoided as there is a continuous synchronization of time between each of the entities modeled. \\
The Timed Automata representation of the proposed two way concurrent buffer system is shown in Fig 3.

\begin{figure}[bpht!]
\centering
	\includegraphics[scale=1.4]{fig3}
	\caption{Timed Automata Representation of Proposed Two Way Buffer System}
\end{figure}

\section{Timing Analysis}
\label{7}
In this section, the timing analysis of individual buffers, entry and exit points has been done.\\
Notations used:
\begin{itemize}
	\item $\zeta$ : Time required over the channel for transfer of packets from one location to next consecutive location
	\item $\theta$ : Delay in transmission time of packet at buffer M1 and Md
	\item $\alpha$ : Iteration time for every transfer of packet from generator to exit in both the directions
\end{itemize}

At G1, time required to move the packet from G1 to M1 is $\zeta$ time units. At M1 it will wait for $\theta$ time units, because it has to wait for M2 to be free, after waiting $\theta$ time units M2 will be  free. So, total time at M1 is $\zeta+\theta$ time units. Time to reach from M1 to M2 will take $\zeta$ time units, so the total time at M2 is 2$\zeta+\theta$ time units.  The time to reach from M2 to E1 is $\zeta$ time units, so the total time at exit point E1 is 3$\zeta+\theta$ time units. Thus the time taken for transfer of packet from G1 to E1 is 3$\zeta+\theta$ time units. Additional delay of $\zeta$ time units is added to perform synchronization with G2. For next iteration G1 will release the packet at 4$\zeta+\theta$ time units.
\par On the other end at entry point G2 the packet is ready at $0^{th}$ time unit. Time taken by a packet to move from G2 to M2 requires $\zeta$ time units, total time required at M2 is $\zeta$ time unit. From M2 to Md , the time required is $\zeta$ time units, and there is a delay in the buffer Md, so total time at Md is 2$\zeta+\theta$ time units. At Md it will wait for $\theta$ time units. Because it has to wait for M1 to be free, after waiting $\theta$ time units M1 will be  free. The time taken to move the packet from Md to M1 is $\zeta$ time unit, now total time at M1 is 3$\zeta+\theta$. On the other side M1 will move the packet to exit E2 and the time taken is $\zeta$ time unit, now the total time at E2 is 4$\zeta+\theta$ time units. At this point on completion of one cycle G2 is ready at 4$\zeta+\theta$ time units. In this way the synchronization is maintained because G1 and G2, releasing the packets at 4$\zeta+\theta$ time units in either direction. Also the concurrency is maintained as packet transfer is concurrent in both the directions throughout the system.\\
Below is the mathematical form of timing analysis and certain assumptions:\\
Assumptions: Let $\alpha$ be a iteration number, $\theta$ be the delay at buffer and $\zeta$ be the delay duration of channel.

\begin{itemize}
	\item Time for subsequent packet transfer at generators G1 or G2 is
\\ For $\alpha$ = 0, the packet is ready at 0 time units;
\\ For $\alpha$ = 1, time consumed is 4$\zeta+\theta$ time units;
\\ For $\alpha$ = 2, time consumed is 8$\zeta+2\theta$ time units;
\\ In general $\displaystyle\sum\limits_{\alpha=0}^{\infty} [\alpha(4\zeta+\theta)]$.\\

At G2, first packet will transfer without any delay at $3\zeta+\theta$, but one additional delay is provided to synchronize with G2.

\item When M1 receives a packet from G1, the timing analysis is done as follows:
\\ For $\alpha$ = 0, time consumed is $\zeta+\theta$ ;
\\For $\alpha$ = 1, time consumed is $5\zeta+2\theta$ units;
\\For $\alpha$ = 2, time consumed is $9\zeta+3\theta$ units;
\\ In general $\displaystyle\sum\limits_{\alpha=0}^{\infty} [\zeta+\theta+\alpha(4\zeta+\theta)]$

\item When M1 receives a packet from Md, the timing analysis is done as follows:
\\ For $\alpha$ = 0, time consumed is $3\zeta+\theta$ ;
\\For $\alpha$ = 1, time consumed is $7\zeta+2\theta$ units;
\\For $\alpha$ = 2, time consumed is $11\zeta+3\theta$ units;
\\ In general $\displaystyle\sum\limits_{\alpha=0}^{\infty} [3\zeta+\theta+\alpha(4\zeta+\theta)]$

\item When M2 receives a packet from M1, the timing analysis is done as follows:
\\ For $\alpha$ = 0, time consumed is $2\zeta+\theta$ ;
\\For $\alpha$ = 1, time consumed is $6\zeta+2\theta$ units;
\\For $\alpha$ = 2, time consumed is $10\zeta+3\theta$ units;
\\ In general $\displaystyle\sum\limits_{\alpha=0}^{\infty} [2\zeta+\theta+\alpha(4\zeta+\theta)]$

\item When M2 receives a packet from G2, the timing analysis is done as follows:
\\ For $\alpha$ = 0, time consumed is $\zeta$ ;
\\For $\alpha$ = 1, time consumed is $5\zeta+\theta$ units;
\\For $\alpha$ = 2, time consumed is $9\zeta+2\theta$ units;
\\ In general $\displaystyle\sum\limits_{\alpha=0}^{\infty} [\zeta+\theta+\alpha(4\zeta+\theta)]$

\item When Md receives a packet from M2, the timing analysis is done as follows:
\\ For $\alpha$ = 0, time consumed is $2\zeta+\theta$ ;
\\For $\alpha$ = 1, time consumed is $6\zeta+2\theta$ units;
\\For $\alpha$ = 2, time consumed is $10\zeta+3\theta$ units;
\\ In general $\displaystyle\sum\limits_{\alpha=0}^{\infty} [2\zeta+\theta+\alpha(4\zeta+\theta)]$

\item When E1 receives a packet from M2, the timing analysis is done as follows:
\\ For $\alpha$ = 0, time consumed is $3\zeta+\theta$ ;
\\For $\alpha$ = 1, time consumed is $7\zeta+2\theta$ units;
\\For $\alpha$ = 2, time consumed is $11\zeta+3\theta$ units;
\\ In general $\displaystyle\sum\limits_{\alpha=0}^{\infty} [3\zeta+\theta+\alpha(4\zeta+\theta)]$
\item When E2 receives a packet from M1, the timing analysis is done as follows:
\\ For $\alpha$ = 0, time consumed is $4\zeta+\theta$ ;
\\For $\alpha$ = 1, time consumed is $8\zeta+2\theta$ units;
\\For $\alpha$ = 2, time consumed is $12\zeta+3\theta$ units;
\\ In general $\displaystyle\sum\limits_{\alpha=0}^{\infty} [4\zeta+\theta+\alpha(4\zeta+\theta)]$
\end{itemize}

The pseudocode for the system proposed:

\begin{verbatim}
// Initial state is current state
while(new_packet=TRUE)
  check(delay)
    if(buffer.next = TRUE)
      delay = 0  //Assigning zero delay
      if(delay = = 0)
        transition(l, l^i)
        l=l'
        l'=next_location
        if(buffer.next = free)
          decrement_delay
          check(delay)

\end{verbatim}

Here $l$, $l^{i}$ are the locations as per the timed automata formalism. $\textit{Check()}$ function verifies if there is any packet on the opposite way being transferred then it will delay the packet in the corresponding buffer.

\section{Uppaal Model of Proposed  Buffer System without Time}
\label{8}
This kind of model of the system has no delay incorporated as there is no time. Hence, the delay variables are all null in this system. This system has a deadlock when both generators produce a packet at the same time which are received by M1, M2 and Md. The buffers wait for each other to be free so they can forward their packets and the system deadlocks.  The location marked with red indicates where the control is at the moment.

\subsection{} Generators G1 and G2 has no delay so they do not wait. as the new packet generates it transfers that to next channel. as shown in the Fig 4.
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.75]{fig4}
		\caption{Generators G1 nad G2 in no time}
	\label{Fig 4}
\end{figure}

\subsection{} Similar to generators, buffers will transfer the packets to next channel as they received it without any delay. as shown in the Fig 5.
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.75]{fig5}
		\caption{Buffers M1, M2 and Md in no time}
	\label{Fig5}
\end{figure}

\subsection{} Similar to Generators and buffers Exits also forward the packets as they receive it. as shown in the Fig 6.
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.75]{fig6}
		\caption{Exit E1 and E2 in no time}
	\label{Fig 6}
\end{figure}

\subsection{} Verifier gives the result for deadlock, because "no time" approach violates the Proposed two way concurrent buffer system timings where we require slight delay at various positions to make the system concurrent. as shown in the Fig 7.
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.70]{fig7}
		\caption{verification result for no time}
	\label{Fig 7}
\end{figure}

\section{Uppaal Model of Proposed Buffer System in Deterministic Time}
\label{9}
In this section, some deterministic time is added to the system. The generator, G1, produces a packet every $\alpha$ seconds and the generator G2 produces packets every $\beta$ seconds. Also, after receiving a packet the intermediate buffers M1,M2 and Md delay for $\gamma$ seconds each before forwarding the packets to the respective exits. This process repeats over time. One possibility is that $\alpha$ is delay 10, $\beta$ is delay 1 and $\gamma$ is delay 2. It can be observed easily that there exists a deadlock in the system at time 10 seconds when both generators produce a packet at the same time.

\subsection{} Generators G1 and G2 will transfer the packets only when clock hits particular time not before and not after. so giving delay with this time system is logically difficult. in Fig 8 below transfer of packets will occur only at 0th and 5th second.
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.75]{fig8}
		\caption{Generators G1 and G2 in deterministic time}
	\label{Fig 8}
\end{figure}

\subsection{} Similar to generators the time is deterministic. the buffers will release the packets only at 0th, 1st, and second second. this is shown in Fig 9
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.75]{fig9}
		\caption{Buffers M1, M2 and Md in deterministic time}
	\label{Fig 9}
\end{figure}

\subsection{} Similar to generators and buffers the exits E1 and E2, release the packets at particular time as shown in the Fig 10 below.
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.75]{fig10}
		\caption{Exit E1 and E2 in deterministic time}
	\label{Fig10}
\end{figure}

\subsection{} The verification results for proposed model is shown in the Fig 11 below. The model is in Deadlock, because deterministic time approach violates the Two way concurrent buffer system timings where we require slight delay on or before a particular time, at various positions to make the system concurrent.
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.70]{fig11}
		\caption{verification result for deterministic time}
	\label{Fig11}
\end{figure}

\section{Uppaal Model of Proposed Buffer System in Non-deterministic Time with no Deadlock}
\label{10}
Timing model on various states in the UPPAAL model for no deadlock is non deterministic time. in the figure below the model is incorporated with the time which is non deterministic. For the proposed two way buffer system, each component of the system is denoted by a separate automaton on the model.

\subsection{} G1 and G2 has been modeled as an automaton with two locations as shown in Fig 12. In G1 on transition, the packet is sent over the channel 'a' by the send operation a!, with a guard condition of $x \geq 0$ and invariant $x\leq 5$ when it moves to next location. So that next packet will not be produced until that time and once the packet is transferred through all the components, the guard condition is checked for $x\geq 5$. In G2 on transition, the packet is sent over the channel 'd' by the send operation d!, with a update of x=0, when it moves to next location and checks for invariant $x\leq 5$ at that location.  This is depicted below in Fig 12.

\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.90]{fig12}
		\caption{Generators G1 and G2  for non deterministic time}
	\label{Fig12}
\end{figure}

\subsection{} The model of buffer M1 consists of three locations. Firstly, when it receives a packet 	from G1 over the channel 'a' by the receive operation a?, with a guard condition $x\geq 0$ and invariant $x\leq 3$ it reaches to the next location. The packet received and is ready to be sent over the channel 'b' by the send operation b! to buffer M2. This is done with a guard condition $x>2$ and it updates the clock time to $x=3$. Similarly when it receives packet from M2 through the channel 'f' by receive operation f?. The packet is then sent to E2 through the channel 'g' through send operation g!. Fig 13 depicts the model of Buffer M1.

\begin{figure}[bpht!]
\centering
	\includegraphics[scale=1.0]{fig13}
		\caption{Buffer M1 for non deterministic time}
	\label{Fig13}
\end{figure}

\subsection{} Buffer M2 is also modeled on similar lines but with different invariant constraints which is shown in Fig 14.

\begin{figure}[bpht!]
\centering
	\includegraphics[scale=1.0]{fig14}
		\caption{Buffer M2 for non deterministic time}
	\label{Fig14}
	\end{figure}

\subsection{} The dummy buffer $M_{d}$ , Clock variable, E1 and E2 are modeled as shown in Fig 15.

\begin{figure}[bpht!]
\centering
	\includegraphics[scale=1]{fig15}
		\caption{Clock, Exits E1 and E2, Dummy buffer Md for non deterministic time}
	\label{Fig15}
\end{figure}

\subsection{} The simulation results of the proposed two way concurrent system given by UPPAAL are shown in the Fig 16 where the packet goes in the same sequence as describes in section VII. In UPPAAL, the simulation has two ways to start, in first way G1 initiate the packet after that G2 which is shown in Fig 16 (a), in second way G2 initiate the packet after that G1, which is shown in Fig 16 (b).
\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.65]{fig16}
		\caption{Simulation results for Two way concurrent buffer system for non deterministic time}
	\label{Fig16}
\end{figure}

We have modeled the proposed two way buffer system in UPPAAL. The main purpose of this work is to verify the model with respect to a requirement specification. Like the model, the requirement specification is expressed in a formally well-defined language \cite{MHMR04}. UPPAAL tool has been used to test three properties \textit{reachability}, \textit{safety (Deadlock)} and \textit{liveness}.

\begin{enumerate}[(a).]
	\item \textbf{Reachability:} The reachability property is described as 'Does there exist a path starting at the initial state, such that $\varphi$ is eventually satisfied along that path' \cite{BDL06}. This specification in UPPAAL is expressed as \textbf{\textit{E$<>$ M1.G2\_receive}}. This property checks for a path from M1 to G2 and the proposed model satisfies this property. The verification result of reachability property is shown in Fig 17.

\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.75]{fig17}
	\caption{Reachability Property satisfied for non deterministic time}
	\label{Fig17}
\end{figure}

	\item \textbf{Safety:} Safety property is of the form 'For the path formula $\varphi$ covering all the states (locations) there is no deadlock in the system'. The safety specification in UPPAAL is expressed as \textbf{\textit{A[] no deadlock}}. This property checks for deadlock and the proposed model satisfy this property. The Fig 18 depicts the  satisfiability of the safety property.

\begin{figure}[bpht!]
\centering
	\includegraphics[scale=0.75]{fig18}
		\caption{Safety Property satisfied for non deterministic time}
	\label{Fig18}
\end{figure}

	\item \textbf{Liveness:} In the model when a packet is sent it should eventually be received. This specification in UPPAAL is expressed as \textbf{\textit{G1.send $--$$>$ E1.receive}}. It checks whether E1 is eventually receiving the packet sent from G1 and in the result which is shown in Fig 19 this property is satisfied.

\begin{figure}[bpht!]
\centering
	\includegraphics[height=6cm, width=8cm]{fig19}
		\caption{Liveness Property satisfied}
	\label{Fig19}
\end{figure}
\end{enumerate}

\section{Buffer Sizing}
\label{11}
Buffer system widely used in routers, Buffer sizing is very important aspect of the buffer system. As the proposed Two Way Buffer System may be over buffered in a high speed network system and be unusable and can lead to many complaints of congestion \cite{REF1}. But keeping aspects of buffer sizing standards can reduce the chances of congestion or any other flaws in high speed networks too.  The paragraph below is giving light on the buffer sizing constraints.

\par Internet routers are packet switches, and therefore buffer packets during times of congestion.  Arguably, router buffers are the single biggest contributor to uncertainty in the Inter-net. Buffers cause queuing delay and delay-variance; when they overflow they cause packet loss, and when they under-flow they can degrade throughput. Given the significance of their role, we might reasonably expect the dynamics and sizing of router buffers to be well understood, based on a well-grounded theory, and supported by extensive simulation and experimentation \cite{REF1}.

\par Router buffers are sized today based on a rule-of-thumb commonly attributed to as discussed in \cite{REF3} . Using experimental measurements of at most eight TCP flows on a 40 Mb/s link, they concluded that because of the dynamics of TCP's congestion control algorithms  a router needs an amount of buffering equal to the average round trip time of a flow that passes through the router, multiplied by the capacity of the router's net-work interfaces.
\par This is the well-known \cite{REF1}
\begin{equation}
B = RTT * C_{rule}
\end{equation}

\par Network operators follow the rule-of-thumb and require that router manufacturers provide 250ms (or more) of buffering \cite{REF4}. The rule is found in architectural guidelines \cite{REF5}, too. Requiring such large buffers complicates router design, and is an impediment to building routers with larger capacity. For example, a 10Gb/s router line card needs approximately 50ms�10Gb/s= 2.5Gbits of buffers, and the amount of buffering grows linearly with the line-rate \cite{REF2}.

\par It is not well understood how much buffering is actually needed, or how buffer size affects net-work performance \cite{REF6}.

\par Overall  study of buffer sizing lights on the usability of small buffers in the high speed network.  significantly smaller buffers could be used in backbone routers  without a loss in network utilization \cite{REF1}.

\section{Scope for Future Enhancement}
\label{12}
The Two way buffer system can be used in the backbone of the network where two way transfer of the packet is required on a single line. The two way system will provide slight delay in the channel and support for concurrent packet transfer. The buffer system mainly used in routers. There are various kind of routers available in the market which fulfills different purposes. The Two way buffer system is most appropriate for core routers which are used in backbone of the network, subscriber edge routers this type of router belongs to an end user (enterprise) organization. The Two Way Buffer System is also appropriate for the system where there is heavy traffic of packets from one way and less traffic from other way.
\par This research is beyond the scope of Ultra high speed network (GBPS) where optical links are used to transfer the data. This can be done by making "Two Way Concurrent Buffer System" more specific to this kind of networks, by keeping constraints of such system in mind.

\section{Conclusion}
\label{13}
The existing two way tiny buffer system when modeled using the UPPAAL tool it satisfy only the safety property only. The system is not time efficient because of unusual delays in transmission of packets of some buffer locations. In this paper we have tried to overcome the drawback of the existing tiny two way buffer system by incorporating a dummy buffer Md. The entire model has been modeled and verified in UPPAAL which incorporates the concepts of Timed Automata. The simulation result has satisfied all the three properties namely Safety, Liveness and Reachability, also it is fully concurrent.

\bibliographystyle{IEEEtran}
\balance