\section*{Supplemental Material}

In this document we provide a detailed analysis of various properties
of the HAL module.

\section{Analysis of full model}

We first consider the full model, which is given by:
\begin{eqnarray}
\left\{ 
\begin{array}{ccl}
\vspace{2mm}
\frac{d[G]}{dt}   &=&\omega (1-[G])   -\alpha [G] [A]        \\
\vspace{2mm}
\frac{d[M]}{dt}   &=&\mu_M [G]         + \mu_M^A(1-[G])       -\delta_M [M]          \\
\vspace{2mm}
\frac{d[A]}{dt}   &=&\mu_A [M]        -\delta_A [A]          -\gamma_{A\! B} [A] [B] 
\\ \vspace{2mm}
&& +\lambda_{A\! B} [AB]  +\omega (1-[G]) -\alpha [G] [A]\\
\vspace{2mm}
\frac{d[B]}{dt}   &=&\mu_B            -\delta_B [B]          -\gamma_{A\! B} [A] [B]
\vspace{2mm}
+\lambda_{A\! B} [AB]\\
\vspace{2mm}
\frac{d[AB]}{dt}  &=&\gamma_{A\! B} [A] [B] -\lambda_{A\! B} [AB] -\delta_{A\! B} [AB]
\end{array}
\right.
\label{suppl:full_model}
\end{eqnarray}
Table~\ref{BiochemicalReactions} lists all the reactions of the HAL
module shown in Fig. 1 of the main text and the corresponding mass action
terms. There are $12$ rate constants.

The first equation governs the time evolution of a variable $[G]$
which represents an average gene activity. Even in the cases where
gene activity is considered as a stochastic variable alternating
between two values (active and inactive), such equations
can be derived from moment expansions of the chemical master
equation~\cite{Wang14:_stoch_oscil_selfrepressed_gene}. The
form used here is valid when the variances of the
stochastic variables can be neglected. If the results of
~\cite{Wang14:_stoch_oscil_selfrepressed_gene} can be transposed here,
this would be the case when $\omega/\delta_M \geq 1$. However, we have
checked with stochastic simulations that this heuristic bound is too
pessimistic, because the main discrepancy observed for lower values of
$\omega/\delta_M$ is only a slightly higher variability in interpulse
time intervals. Thus, Equations~\eqref{suppl:full_model} are adequate
for most parameter sets considered in our analysis.


In the limit of fast gene dynamics, the quasi-steady-state 
approximation $d[G]/dt=0$ yields $[G]= (1+\alpha[A]/\omega)^{-1}$. Substituting this in
\eqref{suppl:full_model}, we get a system of four equations:
\begin{eqnarray}
\left\{ 
\begin{array}{ccl}
\vspace{2mm}
\frac{d[M]}{dt}   &=&\frac{\omega \mu_M + \alpha \mu_M^A [A]}{\omega + \alpha [A]}       
-\delta_M [M] \\
\vspace{2mm}
\frac{d[A]}{dt}   &=&\mu_A [M] -\delta_A [A] -\gamma_{A\! B} [A] [B] 
+\lambda_{A\! B} [AB]\\
\vspace{2mm}
\frac{d[B]}{dt}   &=&\mu_B            -\delta_B [B]          -\gamma_{A\! B} [A] [B]
+\lambda_{A\! B} [AB]\\
\vspace{2mm}
\frac{d[AB]}{dt}  &=&\gamma_{A\! B} [A] [B] -\lambda_{A\! B} [AB] -\delta_{A\! B} [AB]
\end{array}
\right.
\label{full_model_MM}
\end{eqnarray}
which recovers the standard Michaelis-Menten form for the mRNA synthesis.

\begin{figure}[t]
\begin{tikzpicture}
\clip node (m) [matrix,matrix of nodes,
fill=black!20,inner sep=0pt,
nodes in empty cells,
nodes={minimum height=0.6cm,minimum width=0 .5cm,anchor=center,outer sep=0,font=\sffamily},
row 1/.style={align=center,nodes={fill=black,text=white}},
column 1/.style={text width=2.7cm,align=left,every even row/.style={nodes={fill=white}}},
column 2/.style={text width=1.7cm,align=center,every even row/.style={nodes={fill=white}}},
column 3/.style={text width=0.9cm,align=center,every even row/.style={nodes={fill=white}}},
column 4/.style={text width=1.8cm,align=center,every even row/.style={nodes={fill=white}}},
column 5/.style={text width=1.7cm,align=center,every even row/.style={nodes={fill=white}}},
prefix after command={[rounded corners=4mm] (m.north east) rectangle (m.south west)}
] {
	Reaction        & Reactants     & $\xrightarrow{\mbox{rate}}$     & Products      & Term\\
	$\ G$        transcription: & $G$           & $\xrightarrow{\mu_M}$           & $G + M$       &
	$\mu_M [G]$\\
	$\ M$        degradation:   & $M$           & $\xrightarrow{\delta_M}$        & $\varnothing$ &
	$\delta_M [M]$\\
	$\ M$        translation:   & $M$           & $\xrightarrow{\mu_A}$           & $M + A$       &
	$\mu_A [M]$\\
	$\ A$        degradation:   & $A$           & $\xrightarrow{\delta_A}$        & $\varnothing$ &
	$\delta_A [A]$\\
	$\ B$        production:    & $\varnothing$ & $\xrightarrow{\mu_B}$           & $B$           &
	$\mu_B$ \\
	$\ B$        degradation:   & $B$           & $\xrightarrow{\delta_B}$        & $\varnothing$ &
	$\delta_B [B]$\\
	$\ G$        repression:    & $G + A$       & $\xrightarrow{\alpha}$          & $G_{\! A}$    &
	$\alpha [G][A]$\\
	$\ G_{\! A}$ deregulation:  & $G_{\! A}$    & $\xrightarrow{\omega}$          & $G + A$       &
	$\omega (1-[G])$\\
	$\ G_{\! A}$ transcription: & $G_{\! A}$    & $\xrightarrow{\mu_M^A}$         & $G_{\! A} + M$&
	$\mu_M^A (1-[G])$\\
	$\ AB$       complexation:  & $A + B$       & $\xrightarrow{\gamma_{A\! B}}$  & $AB$          &
	$\gamma_{A\! B}[A][B]$\\
	$\ AB$       dissociation:  & $AB$          & $\xrightarrow{\lambda_{A\! B}}$ & $A + B$       &
	$\lambda_{A\! B}[AB]$\\
	$\ AB$       degradation:   & $AB$          & $\xrightarrow{\delta_{A\! B}}$  & $\varnothing$ &
	$\delta_{A\! B}[AB]$\\
};
\end{tikzpicture}
\caption{List of all the biochemical reactions which define the HAL
module. By convention rates are denoted by: $\mu$ for production rates,
$\delta$ for degradation rates, $\alpha$ for binding rates, $\omega$
for unbinding rates, $\gamma$ for complexation rates and $\lambda$
for dissociation rates. The rightmost column gives the corresponding rates
in the differential equations as obtained from mass action kinetics.}
\label{BiochemicalReactions}
\end{figure}

Figure~\ref{fig_phased} shows a plot of ``phase diagrams" of the system.
To compute it we fixed the parameters to the following values $\mu_M^{-1}
= 0.5$, $\delta_M^{-1} =20$, $\mu_A^{-1} = 0.04$, $\mu_B^{-1} =
10^{-2}$, $\delta_A^{-1} =\delta_B^{-1} = 10^3$, $\delta_{A\! B}^{-1}
= 10$, $\gamma_{A\! B}^{-1}= 10$, $\omega^{-1} = 100$, $[A]_0 = 1$.
Six of these parameters $\alpha$, $\omega$, $\mu_M$, $\mu_A$,
$\mu_B$, $\delta_A$ and $\delta_B$ were varied two at a time while
keeping four of them fixed (recall that $\alpha = \omega/[A]_0$).
This procedure generates $6\cdot 5/2 =15$ two dimensional slices of
the phase diagram. The analysis consists in numerically integrating
Eqs.~\eqref{suppl:full_model} for every set of input rates while
identifying if the solution is pulsating or stationary.  The pulsating
domain is show as black in Fig.~\ref{fig_phased}. The axes in the phase
diagrams in Fig.~\ref{fig_phased} are in logarithmic scale and each
axes covers a variation of two orders of magnitude centered around the
values of rates given above. Hence the selected point is rather far from
the phase boundaries.  Apart from the pulsating solution we distinguish
two types of stationary solutions with high $A$ (high $B$) shown as red
(green) in Fig.~\ref{fig_phased}. In these phases one of the two proteins
has typically much higher concentration than the other.

\begin{figure*}[t]
\includegraphics[width=0.45\columnwidth]{theta_alpha.pdf}
\includegraphics[width=0.45\columnwidth]{alpha_dM.pdf}
\includegraphics[width=0.45\columnwidth]{theta_muM.pdf}

\includegraphics[width=0.45\columnwidth]{alpha_muM.pdf}
\includegraphics[width=0.45\columnwidth]{theta_dM.pdf}
\includegraphics[width=0.45\columnwidth]{theta_muA.pdf}

\includegraphics[width=0.45\columnwidth]{muM_dM.pdf}
\includegraphics[width=0.45\columnwidth]{muM_muA.pdf}
\includegraphics[width=0.45\columnwidth]{muM_muB.pdf}

\includegraphics[width=0.45\columnwidth]{muA_muB.pdf}
\includegraphics[width=0.45\columnwidth]{alpha_muB.pdf}
\includegraphics[width=0.45\columnwidth]{alpha_muA.pdf}

\includegraphics[width=0.45\columnwidth]{dM_muA.pdf}
\includegraphics[width=0.45\columnwidth]{dM_muB.pdf}
\includegraphics[width=0.45\columnwidth]{theta_muB.pdf}
\caption{Slices of phase diagrams of the model \eqref{suppl:full_model}
depicting various phases as a function of the six ``most relevant"
parameters $\alpha$, $\omega$, $\mu_M$, $\mu_A$, $\mu_B$, $\delta_A$
and $\delta_B$.  The black area corresponds to the pulsating domain,
while the red (resp. green) to steady states characterized by high
concentration of A (resp. B).}
\label{fig_phased}
\end{figure*}

\FloatBarrier

\clearpage

\begin{figure}[t]
\includegraphics[width=1.0\columnwidth]{suppl_fig.pdf}
\includegraphics[width=1.0\columnwidth]{suppl_fig2.pdf}
\vspace{-6mm}
\caption{Plots of protein concentrations vs. time (a$_1$-d$_1$) and of
mRNA concentrations vs. time (a$_2$-d$_2$) for four different sets of rate
constants. In the case (a) the rates are those reported in Fig.~2 of the
main paper. (b,c and d) are obtained by (a) by a change of a single rate
constant as follows: (b) $\mu_B$ is doubled ($\mu_B^{-1}=0.00075$),
(c) $\omega$ is halved ($\omega^{-1}=200$) (d) $\mu_A$ is halved
($\mu_A^{-1}=0.1333$).}
\label{fig_param_variation}
\end{figure}

Figure~\ref{fig_param_variation}(a-d) shows the effect on change of
rates on the concentrations of the proteins and mRNA in the pulsating
regime. Fig.~\ref{fig_param_variation}(a) reproduces the same rates as
in Fig.~2 of the main paper, while the cases (b), (c) and (d) correspond
to variations of a single rate with respect to the case (a). The duration
of the $A$ peaks is quite robust against the parameters variation, while
$T_b$ varies: (a) $T_a=2.01$, $T_b=2.13$ (b) $T_a=1.71$, $T_b=3.81$
(c) $T_a=1.93$, $T_b=3.48$ and (d) $T_a=1.76$, $T_b=3.76$.  In the case
(b) the $B$ protein synthesis rate $\mu_B$ is doubled with respect of
(a). This has a strong effect in the height and duration of the peaks of
$B$, but a milder effect on the duration of the peaks of $A$. Halvening
the value of (c) $\omega$ and of (d) $\mu_A$ is also affecting strongly
the peaks of $B$. While the duration of the peaks of $A$ is robust, their
height is not. This appears to be mostly affected by a change in $\mu_B$,
in agreement with the analysis of the reduced model in the next section.

\section{Analysis of the reduced model}

We present here the details of the analytical calculations for the
durations of the A and B phases for the reduced model:
\begin{subnumcases}{}
\frac{dg}{d\tau}\,   = \Omega (1-g)       -\sigma g a \label{equG}\\
\frac{dm}{d\tau}     = \quad g     \qquad -       m   \label{equM}\\
\frac{da}{d\tau}\,   = \quad k_a  m \quad -       a b \label{equA}\\
\frac{db}{d\tau}\,   = \quad k_b   \qquad -       a b \label{equB}
\end{subnumcases}

As shown in Fig.~2 of the main text two phases can be identified in
the pulsating domain: in one phase, $[A]$ is peaked and $[B]$ is small,
while in the other phase $[B]$ is peaked and $[A]$ small. We refer to
these as to the A-phase and to the B-phase, respectively. The two phases
are due to the mutual sequestration of $A$ and $B$. In the A-phase the
$G_A$ gene is strongly repressed, the mRNA synthesis is stopped and the
mRNA concentration $m(t)$ decreases due to degradation. As long as $k_a
m > k_b$ (rescaled variables and parameters, see Eqs.~(\ref{equA}) and
(\ref{equB})) the production of $a$ dominates over the production of
$b$. Once the mRNA concentration drops and $k_a m < k_b$, the production
of $B$ becomes dominant and the concentration of $A$ starts decreasing
till $A$ is completely sequestered out of the system and the transition to
the B-phase is made. In the B-phase the protein $A$ is released from its
gene $G_A$ promoter site. The absence of repression produces a rise in
the mRNA concentration; when the mRNA concentration reaches the threshold
value $k_am = k_b$, the concentration of $B$ starts decreasing and one
is back to the A-phase again.

We compute now the duration of the two phases.  Let us start from the
A-phase.  A first assumption is that the gene is constantly repressed
when the concentration of $A$ is high, hence $g=0$. We can thus eliminate
the variable $g$ from Eq.~(\ref{equM}) to obtain the solution:
\begin{equation}
m(\tau) = m_A\, e^{-\tau}
\label{MinAphase}
\end{equation}
where we set the origin of time $\tau=0$ at the beginning of
the A-phase. Using the same notation as the paper $m_A$ and $m_B$ 
indicate the mRNA concentrations at the beginning of the A- and of
the B-phases. A second assumption is that in the A-phase the
concentration of $B$ is stationary hence $db/dt=0$, which implies $k_b=ab$
from Eq.~(\ref{equB}). Substituting this into Eq.~(\ref{equA}) and using
(\ref{MinAphase}) we get the following equation for the evolution of $a$:
\begin{equation}
\frac{da}{d\tau} =	k_a m_A\, e^{-\tau} - k_b
\label{dadt_approx}
\end{equation}
The solution of the previous equation with initial condition $a(0)=0$ is
\begin{equation}
a(\tau) = k_a m_A \, \left( 1 - e^{-\tau} \right)  - k_b \tau
\label{eq:a_tau}
\end{equation}
which is a function with a single maximum beyond which it decreases
monotonically and it becomes negative at long times, which is obviously
an unphysical result.  We can estimate the duration of the
A-phase from the requirement $a(T_a)=0$, which gives:
\begin{equation}
\frac{T_a}{ 1 - e^{-T_a}} = \frac{k_a m_A }{k_b} 
\label{eqTa}
\end{equation}
For the B-phase we assume that the concentration of free $A$ in
solution is very small so that the binding rate to the gene promoter
site is negligible. We can approximate Eq.~(\ref{equG}) with $dg/d\tau
\approx \Omega (1-g)$, from which we get the following solution:
\begin{equation}
g(\tau) = 1 - e^{-\Omega (\tau-T_a)} 
\label{g_phaseB}
\end{equation}
where we used the initial condition $g(T_a)=0$ in the B-phase, $T_a
\leq \tau \leq T_a + T_b$.  We approximate further the previous expression
to the first order in the exponential:
\begin{equation}
g(\tau) \approx \Omega (\tau-T_a) 
\label{g_phaseB2}
\end{equation}
and which is valid for $\tau-T_a \ll \Omega^{-1}$. We now plug  in the
previous expression into Eq.~(\ref{equM}) and solve it to get for the
mRNA concentration in the B-phase ($T_a \leq \tau \leq T_a+T_b$):
\begin{eqnarray}
m(\tau) &=& \Omega \left(\tau - T_a - 1\right) +
e^{-\left(\tau-T_a \right)} \left(m(T_a)+ \Omega\right) \nonumber\\
        &=& \Omega \left(\tau - T_a - 1\right) +
e^{-\left(\tau-T_a \right)} \left(m_A  e^{-T_a} + \Omega\right)
\label{MinBphase}
\end{eqnarray}
where we have used Eq.~(\ref{MinAphase}): $m(T_a) = m_A  e^{-T_a}$.
We proceed as done for the A-phase. We assume that $a$ is stationary
in the B-phase, i.e. $da/d\tau =0$ which yields $k_a m = ab$
(Eq.~(\ref{equA})). Substituting this result in Eq.~(\ref{equB}) we get
the following Equation for the growth of $b$:
\begin{equation}
\frac{db}{d\tau} = k_b - k_a m
\end{equation}
with $m(\tau)$ given by Eq.~(\ref{MinBphase}). 
Using the initial condition $b(T_a) =0$ we get:
\begin{eqnarray}
b(\tau) &=& k_b \left(\tau-T_a \right) - k_a 
\Omega \left[ \frac{1}{2} \left( \tau - T_a \right)^2 
- \left( \tau - T_a \right) \right] \nonumber\\
&-& k_a \left(1 - e^{-(\tau-T_a)}\right)
\left(m_A  e^{-T_a} + \Omega \right)
\end{eqnarray}
We obtain the length of the B-phase from the requirement 
that $b(T_a + T_b) =0$, which leads to the following relation
\begin{equation}
\left(k_b+k_a\Omega\right)T_b = \frac{1}{2} k_a\Omega T_b^2 +
k_a \left(1 - e^{-T_b}\right)
\left(m_A  e^{-T_a} + \Omega \right)
\label{BinBphase}
\end{equation}
An additional relation is obtained by requiring that at the end of
B-phase: $m(T_a+T_b)=m_A $ which yields from Eq.~(\ref{MinBphase}):
\begin{equation}
m(0) = \Omega \frac{ -1+ T_b + e^{-T_b}}{1 - e^{-T_a - T_b }}
\label{pbc}
\end{equation}
Inserting the previous equation in Eq.~\eqref{eqTa} we get:
\begin{equation}
\frac{T_a}{e^{T_a}-1} =  \beta \frac{-1+T_b+e^{-T_b}}{e^{T_a}-e^{-T_b}}
\label{per1_final}
\end{equation}
where we defined $\beta \equiv k_a \Omega/k_b$.
We now use Eq.~\eqref{eqTa} to get an expression for $m(0)$ which we substitute
in \eqref{BinBphase} to get:
\begin{equation}
\frac{T_a}{e^{T_a}-1} = 
\frac{\beta \left(T_b - \frac{T_b^2}{2}\right) +  T_b}{1 - e^{-T_b}} - \beta
\label{per2_final}
\end{equation}
We also note that for $\beta=2$ the exact solution of
Eqs.~(\ref{per1_final}) and (\ref{per2_final}) is $T_a=T_b=2$, i.e. the
two phases have equal duration.  For $\beta > 2$ ($\beta < 2$) one has
$T_a > T_b$ ($T_a < T_b$). In terms of the original kinetic constants,
the parameter $\beta$ reads:
\begin{equation}
\beta = \frac{\omega \mu_M \mu_A }{\mu_B \delta_M^{\, 2}}
\end{equation}
It characterizes the relative importances of the A-phase and B-phase.
The A-phase dominates if $\omega$ (unbinding rate of the repressor A
from its gene), $\mu_M$ (mRNA synthesis rate) or $\mu_A$ (protein A
synthesis rate) are large. The B-phase is favored when $\mu_B$ (protein
B synthesis rate) or $\delta_M$ (mRNA degradation rate) are large.

\subsection{On the robustness of $T_a$}

We analyze now the dependence of $T_a$ and $T_b$ on $\beta$.
Eq.~\eqref{per2_final} is of the form
\begin{equation}
f_1(T_a) = f_2(T_b,\beta)
\label{func_eq}
\end{equation}
where $f_1$ and $f_2$ are the following functions:
\begin{equation}
f_1(x) = \frac{x}{e^x-1}
\end{equation}
and
\begin{equation}
f_2(x,\beta) = \frac{\beta \left(x - \frac{x^2}{2}\right) +  x}{1 - e^{-x}} - \beta
\end{equation}
For any $x > 0$, the function $f_1(x)$ satisfies $0 < f_1(x) < 1$.  This
implies that $T_b$, solution of \eqref{per1_final} and \eqref{per2_final}
must be such that $0<  f_2(T_b,\beta) < 1$.

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth]{fig_supp03.pdf}
\caption{Plot of the function $f_2(x,\beta)$ showing that 
the solutions for $T_b$ of Eq.~\eqref{func_eq} are in a limited
range of $x$.}
\label{fig_last}
\end{figure}

Figure~\ref{fig_last} shows a plot of $f_2(x,\beta)$ for two values of
$\beta$. For large values of $x$, the function becomes negative and
its value satisfies $0< f_2(T_b,\beta) < 1$ only for a limited
range of $x$. This range varies strongly with $\beta$, which implies a
variation of $T_b$ with $\beta$. The analysis of $f_2(x,\beta)$
shows that $\lim_{\beta \to 0} T_b = \infty$ and 
$\lim_{\beta \to \infty} T_b =0$.

To proceed further we combine \eqref{per1_final} and \eqref{per2_final}
to eliminate $\beta$. We obtain:
\begin{equation}
\frac{T_a}{e^{T_a}-1} = 
\frac{T_b}{1-e^{-T_b} - 
\frac{T_b-\frac{T_b^2}{2}-1+e^{-T_b}}{-1+T_b+e^{-T_b}} \left(e^{T_a}-e^{-T_b}\right)}
\label{Ta_vs_Tb}
\end{equation}
In the limit $T_b \to 0$ (large $\beta$) the previous relation becomes: 
\begin{equation}
\frac{T_a}{e^{T_a}-1} = \frac{3}{e^{T_a}+2}
\label{eqTa_beta-infty}
\end{equation}
which has as unique solution $T_a \approx 2.149$.
In the opposite limit $T_b \to \infty$ we get from \eqref{Ta_vs_Tb}:
\begin{equation}
\frac{T_a}{e^{T_a}-1} = 2 \, e^{-T_a}
\label{eqTa_beta0}
\end{equation}
which has as solution $T_a \approx 1.594$. Hence this analysis shows
that while $T_b$ is unbounded and assumes any positive values when
$\beta$ is varied, $T_a$ is bounded in the interval $[1.594,2.149]$.
As discussed in the paper, the changes in rate constants which could
potentially affect $T_a$ are compensated by a change in $m_A$, the mRNA
concentration at the beginning of the A-phase, such that the ratio
\begin{equation}
c \equiv \frac{k_a m_A}{k_b}
\label{eq:def_c}
\end{equation}
remains constant. Using Eq.~\eqref{eqTa}, we can compute a range for $c$
using the estimated range of values of $T_a$. The result is $2 \leq c
\leq 2.4$.

\subsection{On the amplitude of $A$}

One can get some insights on the amplitude of $A$ from the analysis 
of the simplified model. The maximum of $a$ is obtained from 
Eq.~\eqref{eq:a_tau}:
\begin{equation}
\max_\tau {a} = k_b \left( \frac{k_a m_A}{k_b} - 1 - \log \frac{k_a m_A}{k_b}
\right)
\end{equation}
which shows that this quantity is not robust. Indeed, we have shown that
the solution of Eqs.~\eqref{per1_final} and \eqref{per2_final} are
such that the ratio \eqref{eq:def_c} is robust. The maximum of $a$
depends on this ratio, but it is also is proportional to $k_b$.
Transforming back to the original concentration units we find for 
the peak of $A$
\begin{equation}
\max_t {[A]} = \frac{\mu_B}{\delta_M} \left( c - 1 - \log c \right)
\end{equation}
where $c$ is defined in \eqref{eq:def_c}. This suggests that,
besides fixing $\delta_M$ which determines the overall timescale,
to control the height of the peaks of $[A]$ one needs to control
$\mu_B$, the $B$ production rate.  This is consistent with the plots of
Fig.~\ref{fig_param_variation}: the height of the $A$ peaks is mostly
affected by a change in $\mu_B$ (case (b)).

\begin{figure}[t]
\includegraphics[width=0.95\columnwidth]{histogram_TaTb.pdf}
\caption{Histograms of $T_a$ and $T_b$, the duration of the two
phases obtained from stochastic simulations with the Gillespie
algorithm~\cite{gill77} using two different choices of rates corresponding
to two values of $\beta$. While $\beta$ varies of almost an order of
magnitude the distibution of $T_a$ is weakly modified, whereas $T_b$
is strongly affected.}
\label{fig:stoch_distribution}
\end{figure}

\section{Stochastic analysis}

We extended the analysis of the HAL to the stochastic regime, performing
simulations using the Gillespie algorithm. Typical outputs of these
simulations are given in Fig.~5 of the main text, which shows that the
protein concentrations evolve through peaks of variable duration and
height due to stochastic fluctuations. To quantify the variability
in the dimensionless durations of the two phases, we studied their
probability distribution for two different parameter sets, as shown in
Figure~\ref{fig:stoch_distribution}. The two top graphs are obtained
using the parameter values $\mu_M^{-1} = 1.11$, $\delta_M^{-1} =
16.67$, $\mu_A^{-1} = 0.59$, $\mu_B^{-1} = 0.05$, $\delta_A^{-1} =
\delta_B^{-1} = 10^3$, $\delta_{A\!B}^{-1} = 10$, $\gamma_{A\! B}^{-1}
= 0.02$, $\omega^{-1} = 100$, $[A]_0 = 1$ (and $\lambda_{A\! B}$ and
$\mu_M^A$ fixed as in Table I). These constants correspond to a value of
$\beta=0.21$. The two bottom graphs are generated using the same rates
except for $\delta_M^{-1} =50$, which corresponds to $\beta=1.91$. We
note that $\beta$ strongly influences the duration of the peaks of $B$;
conversely the distribution of $T_a$ is only weakly affected while $\beta$
is varied of almost an order of magnitude. This is consistent with the
deterministic analysis developed in the main text.

\begin{table}[t]
\centering
\begin{tabular}[c]{c|c|ccc}
      $\beta$ &Quantity & Stochastic & Deterministic & Analytical\\ \hline
      0.21 & $\langle T_a \rangle$ & 1.47 & 1.30 & 1.70 \\
           & $c_v(T_a)$ & 0.49 \\
           &$\langle T_b \rangle$ & 1.68 & 1.48& 10.00\\
                      & $c_v(T_b)$ & 0.60 \\ \hline
                      1.91 & $\langle T_a \rangle$ & 1.71& 1.44& 1.98\\
           & $c_v(T_a)$ & 0.44 \\                      
           &$\langle T_b \rangle$ & 0.51& 0.35& 0.50\\
           & $c_v(T_b)$ & 0.83 \\ \hline
    \end{tabular}
\caption{Average durations obtained from (i) stochastic simulations
of the chemical reactions in Fig.~\ref{BiochemicalReactions},
(ii) deterministic simulation of Eqs.~\eqref{suppl:full_model} and
(iii) Eqs.~\eqref{per1_final} and \eqref{per2_final} (analytical
approximation).}
\label{tab:indicat}
\end{table}

\FloatBarrier

The average values $\langle T_s \rangle$ and coefficients of
variation $c_v(T_s) = \sigma(T_s)/\langle T_s \rangle$, where
$\sigma(T_s)$ is the standard deviation of $T_s$, are given in
Table~\ref{tab:indicat}, which compares them to the values of $T_a$
and $T_b$ obtained from the deterministic simulations of the full
model [Eqs.~\eqref{suppl:full_model}], as well as from the analytical
approximation [Eqs.~\eqref{per1_final} and \eqref{per2_final}].
The values obtained confirm that $T_a$ is much less sensitive than
$T_b$ to $\beta$. For each value of $\beta$ and each average duration,
the agreement between the three estimates is reasonable except for
$\langle T_b \rangle$ at low $\beta$. This seems to indicates that
for some parameter sets with a low value of $\beta$, the analytical
approximation severerely overestimates $T_b$, perhaps because it
misses an ingredient leading to a faster dynamics.  However, this does
not affect our conclusion that $T_a$ is relatively constant, nor that
$T_b$ is largely tunable.  Table~\ref{tab:indicat} also shows that the
stochastic variability affects more the duration of the $B$ phase than
that of the $A$ phase.

Summarizing, the stochastic analysis of the system supports the
conclusions drawn from the study of the deterministic model:
while the A-phase appears to be rather robust, the duration of the 
B-phase is tunable and more subject to stochastic fluctuations.
One interesting issue to be left for future investigations is 
whether it is possible to find simple extensions of this genetic
module for which the stochastic fluctuations in $T_a$ can be 
further reduced.

