\section{Introduction}
For some time now, fluorescent microscopy methods have provided quantitative measurements of gene expression cell by cell, see for instance \cite{elowitz_stochastic_2002,ozbudak_regulation_2002}.  The measurements have shown that the protein production is highly variable process, even for genetically identical cells in constant environmental conditions.  The fluctuations can affect  genetic expression and impact the behavior of the  cell, see \cite{losick_stochasticity_2008}, or, on the contrary be beneficial as they can participate in strategies to adapt to changing environment \citep{balaban_bacterial_2004,acar_stochastic_2008}.

More recently, \cite{taniguchi_quantifying_2010} performed an extensive quantification of the variability of the gene expression of around a thousand different genes in \emph{E. coli}. For each gene, an experiment measured in each cell the total fluorescence of the protein associated to the gene; this quantity was then normalized by the cell volume and the fluorescence of a single protein; then results in the protein concentration in each cell of the population.  Furthermore for a significant part of the considered genes, the messenger RNA productions are also quantified in each cell: mRNA-sequencing technique (RNA-seq), they were able to estimate the average production of 841 different types of mRNA; and using fluorescence \emph{in situ} hybridization, they could even measure of some types of mRNAs (137 types) with a single molecule precision in each cell of a population. Statistics over the population then give the average protein and RNA messenger production and the coefficient of variation (CV) for both quantities.

By comparing the behavior of the CV with the one associated with a classical two-stage model of the gene expression proposed in \cite{rigney_stochastic_1977,berg_model_1978} and reviewed in \cite{paulsson_models_2005} (see below), they show that the behavior of the messengers variability with respect to their abundance resembles the one expected on the basis of the two-stage model.  For protein variability, they identify two regimes of protein variability depending on the average protein concentration. For low expressed proteins, the protein CV is shown to be inversely proportional to the average production, in accordance with the one expected with the two-stage model (as it it shown in \cite{taniguchi_quantifying_2010}). By contrast, for highly produced proteins, the CV becomes independent of the average protein production and the quantified variability is then significantly larger than the one expected on the basis of the two-stage model. This effect does not seem specific to the type of protein as all the highly expressed proteins are similarly impacted; therefore a possible gene-specific phenomenon (such as the variability induced by the regulation of the gene) seems unlikely to explain the additional noise observed. Therefore authors then proposed in their study cell-scale phenomenons as possible explanations for this shift, in particular the fluctuations in the availability of RNA-polymerases and ribosomes in the production of different proteins. But other cell-scale mechanisms can also contribute significantly to the protein variability as for example the partitioning at division where each compounds (mRNA or protein) can go either to the two daughter cells \citep{huh_non-genetic_2011,huh_random_2011,swain_intrinsic_2002}, the gene replication where the transcription rate is doubled at some point in the cell cycle \citep{swain_intrinsic_2002}. Finally, some assumptions considered in the two-stage model are questionable when it is applied in the context of Taniguchi~et~al. as for example the fact that a death process acting on the proteins models the dilution effect due to the cell volume growth (see \cite{fromion_stochastic_2013}).

The extensive study of Taniguchi~et~al. gives an extensive set of measurements for a majority of genes of \emph{E. coli}. This turns out very useful to link theoretical models with experimental data. Not only to determine parameters of the models but  also to validate (or not) the predictions of these models to experimental results. This has been developed by Taniguchi~et~al. for the two-stage model, they have concluded that cell-scale mechanisms lack in this classical model to fully reflect the experimental protein variability observed. Yet they did not try with models that include cell-scale mechanisms to definitely corroborate their hypothesis.

To tackle this limitation, we propose in this paper a model that integrates several cell mechanisms which are not present in the two-stage model and which are usually considered as possible contributors to the variability of proteins. To fix the model parameters, we used the measures of 841 genes provided by Taniguchi~et al. (those for which both the average mRNA and protein production have been measured) to conduct a simulation derived from Gillespie algorithm. The results of those simulations reveal and quantify the contribution to the gene expression variability of each of the cell mechanisms considered.

In the next subsection, we present a review of the classical two-stage model as it is broadly used in the literature. We will explicitly present its limitations in regards to cell-scale mechanisms that are suspected of playing a determinant role in protein variability. We then present  aspects that need to be changed in order to represent these cell-scale mechanisms. The Section~\ref{results} presents the predictions and the results in three steps: we successively study the effect on the protein variance of random partition at division, of the gene replication and finally of fluctuations in the availability of RNA-polymerases and ribosomes. The Section~\ref{methods}  exhaustively describe our complete model, derived from the two-stage model and the procedure used to analyze it.


\subsection{\label{twostagemodel} Limitations of the Two-Stage Model}

We based our approach on the one of the simplest classical stochastic model of protein production, the two-stage models \citep{rigney_stochastic_1977,berg_model_1978,paulsson_models_2005,shahrezaei_analytical_2008}.  This model describes the production and degradation of mRNAs and proteins of one particular type. In contrast to a three-stage model \citep{paulsson_models_2005}, it considers a constitutive gene as it does not integrate a regulatory stage at the transcription initiation level and thus represents the expression of a constitutive gene. The additional noise observed in Taniguchi~et~al. uniformly impact highly expressed proteins, it therefore seems to be the result of cell-scale phenomena. It is therefore unlikely that the gene-specific variability induced by the gene regulation could explain the additional variability. We will come back to this aspect in the discussion section of the paper.

The two-stage model represents the evolution of the random variables $M$ and $P$ respectively representing the number of mRNAs and proteins associated with this gene inside a single cell (see Figure~\ref{fig:mod_paulsson_scheme}).  It  is important to note that the quantities considered in this model are numbers. In Figure~\ref{fig:mod_paulsson_scheme}, $M$ and $P$ respectively denotes the mRNA and protein number inside the cell; but their respective concentrations are not explicit in the model as there is no notion of cell volume in this model.

Both the transcription and the translation mechanisms are represented in the model, as well as mRNA degradation.  There is also a protein decay mechanism that can either represents the degradation of the protein, or the dilution effect induced by the growth of the cell volume, since the proteolysis often occurs in a timescale much longer than the cell cycle (see \cite{koch_protein_1955}). All the events, like the productions and degradations of mRNAs and proteins, are considered as occurring at random times exponentially distributed whose rates depend on the current state of the system.
\begin{figure}
  \begin{center}
\begin{tikzpicture}
\matrix [column sep=20mm, row sep=10mm,ampersand replacement=\&] {
  \node (g) [draw=none,fill=none] {}; \&
  \node (m) [draw, shape=rectangle] {$M$}; \&
  \node (p) [draw, shape=rectangle] {$P$}; \\
  \&
  \node (em) [draw=none,fill=none] {$\emptyset$}; \&
  \node (ep) [draw=none,fill=none] {$\emptyset$}; \\
  };
\draw[-> ] (g) -- (m) node[above,midway,align=left]{$\lambda_1$ };
\draw[-> ] (m) -- (em) node[right,midway] {$\sigma_1M$};
\draw[-> ] (m) -- (p) node[midway,above] {$\lambda_2M$};
\draw[-> ] (p) -- (ep) node[right,midway] {$\sigma_2P$};
\draw (m.north)+(0,0.5) node[above]{mRNA};
\draw (p.north)+(0,0.5) node[above]{Protein};

\draw (m.west) node[below left,inner sep=2] {\tiny +1};
\draw (m.south) node[below left,inner sep=2] {\tiny -1};
\draw (p.west) node[below left,inner sep=2] {\tiny +1};
\draw (p.south) node[below left,inner sep=2] {\tiny -1};
\end{tikzpicture}
\end{center}
\caption{\label{fig:mod_paulsson_scheme}Classical two-stage model of one constitutive gene. $M$ and $P$ represent the respective number of mRNAs and proteins in a given cell. \rmodif{The parameter $\lambda_1$ is the rate at which mRNAs are produced, $\lambda_2M$ will be the rate at which proteins are translated, $\sigma_1$ indicates the mRNA degradation rate and $\sigma_2$ is the rate at which protein leaves the cell (by the effect of volume growth and division). }}
\end{figure}
The simplicity of the two-stage model allows the identification of closed-form formulas for the mean and the variance of mRNAs and proteins (see \cite{paulsson_models_2005}); and as it only represents mechanisms that are specific to the protein production, the variance predicted can be naturally considered as representative of the intrinsic noise, as it has been done by Taniguchi~et~al. But the model lacks several aspects that are yet mandatory to make quantitative comparisons with experimental results and it neglects several mechanisms that might have an impact on the protein variance.

\begin{description}


\item [No representation of the volume] The fact that the cell volume is not explicitly represented in the two-stage model (and its derivative) does not seem to have been specifically highlighted in the literature. It appeared to us that it is an important limitation of classical models (such as the two-stage model) for the three intertwined following reasons.
    \begin{description}

    \item [Quantitative comparison with experimental concentrations] As previously said, the two-stage model describes the evolution of the \emph{number} of messengers and proteins in a single cells. It appears that the measures of Taniguchi~et~al. represent concentrations, as this quantity is often a more relevant than the numbers as it determines each reaction speed through the law of mass action. The direct quantitative comparison of mean and variances obtained through the two-stage model and those obtained experimentally are not completely licit. In order to have a model for growing bacteria than can be compared with the experimental data, it is necessary to explicitly describe the volume, in order to represent concentrations in the cytosol, and not only numbers.

    \item [Impact of volume evolution on concentration] In real cells, the simple growth of the cell volume mechanically has a direct impact on compound concentrations through dilution. This effect of dilution of proteins is approximated in the two-stage model by a degradation process at random times exponentially distributed whose rate $\sigma_2P$ depends on the volume growth rate and the protein number. But in reality, dilution is a different process and should be represented as an explicit volume increase rather than a protein disappearance rate. Furthermore, the effect of volume growth on messenger concentration is completely neglected in the two-stage model (due to their quick degradation).

    \item [Impact of volume evolution on reaction rates] The volume growth has also another impact on the protein production as it tends to slow down the reaction rates by diminishing the reactant concentrations in the law of mass action. It will have important consequences as we will explicitly represent RNA-polymerases and ribosomes (see below), and the volume growth will impact their respective concentrations and hence on the global protein production rate.

    \end{description}

\item [No partition at division] The two-stage model represents the expression of one gene in a single cell, so protein decay represents the tendency of proteins to disappear from the cell through the effect of division. In reality, division is a relatively sudden event that partition both mRNAs and proteins in the two daughter cells.

\item [No gene replication] The two-stage model does not consider genetic replication events, it implicitly assumes that the gene promoter number is the same, thus keeping the rate of transcription $\lambda_{1}$ constant. In reality, liked to the volume variation across the cell cycle, the gene activity is linked to its promoter concentration that tends to decrease with the cell cycle, until the promoter is replicated.  At this point the rate of transcription is doubled, thus possibly inducing a possible transcription burst.

\item [Constant availability of RNA-polymerases and ribosomes] RNA-polymerases and ribosomes are needed in protein production for respectively  transcription and translation. As these resources are shared among all the productions of all the different proteins, their availability fluctuates through time. Since the two-stage model is a gene-centered model (it represents the expression of only one gene), it is not able to represent the competitive interactions between the different productions for available RNA-polymerases and ribosomes. In particular, the mRNA production rate $\lambda_{1}$ and protein production per mRNA rate $\lambda_{2}$ are constant, as if the concentrations of available RNA-polymerases and ribosomes remain constant across the cell cycle.
\end{description}
Even if some of these external mechanisms have been theoretically studied independently (random partitioning at division in \citep{huh_random_2011,soltani_intercellular_2016}, gene replication \citep{soltani_intercellular_2016} or fluctuations in the availability of ribosomes \citep{fromion_stochastic_2015}), a complete analysis of their impact on the protein variance has not been performed. Moreover, their quantitative contributions to the global protein variability, with realistic biological parameters, remain largely unknown.




\section{Results}\label{results}

\rmodif{In the Materials and Methods section (Section~\ref{methods}), we describe in details our approach that considers models that integrates the features that are missing in the two-stage model previously described. In order to determine the relative impact of each of these aspects, we have proceed in successive steps of increasing complexity. Bellow is presented the results concerning three intermediate models that successively incorporate a specific feature; the the impact of each of these features on the protein variance will be studied one at a time. }The three intermediate models are the the following:
\begin{itemize}
\item The first intermediate model only consider the growth of the cell and study the impact of the partitioning at division; in particular, each gene concentration, as well as the free RNA-polymerase and ribosome concentrations, are considered as constant during the whole cell cycle. It is focusing on the effect of the partition at division on the protein variance. Its results are presented in Section~\ref{subsec:model1}).
\item The second intermediate model then study the effect of gene replication, while still considering cell volume growth and partitions at division but keeping free RNA-polymerase and ribosome concentrations as constant.  The corresponding results are presented in Section~\ref{subsec:model2}.
\item Then is considered the complete model, with fluctuating free RNA-polymerases and ribosomes, replication, volume growth and partitions at division. Its results are presented in Section~\ref{subsec:model3}.
\end{itemize}
For each of the models, the parameters are fitted in order to correspond to each of the genes of Taniguchi~et~al. and the protein concentration variance is determined using the procedure described in Subsection~\ref{method:analysis}.

\subsection{\label{subsec:model1}Impact of Random Partitioning}

\begin{figure}
\subfloat[\label{fig:mod1_principle}]{\begin{centering}
\includegraphics[width=6cm]{model1/Figure_20160721_BioModel_ultime11}
\par\end{centering}
\centering{}}~~~\subfloat[\label{fig:mod1_scheme}]{\begin{tikzpicture}
\linespread{1}
\matrix [column sep=13mm, row sep=10mm,ampersand replacement=\&] {
  \node (g) [draw=none,fill=none] {}; \&
  \node (m) [draw, shape=rectangle] {$M$}; \&
  \node (p) [draw, shape=rectangle] {$P$}; \\
  \&
  \node (em) [draw=none,fill=none] {$\emptyset$}; \&
  \node (ep) [draw=none,fill=none,color=red,align=center] {
  Exact/random  \\   partitioning every $\tau_{D}$}; \\
  };
\draw[-> ] (g) -- (m) node[above,midway,align=left]{$V(s)\cdot\lambda_1$ };
\draw[-> ] (m) -- (em) node[right,midway] {$\sigma_1M$};
\draw[-> ] (m) -- (p) node[midway,above] {$\lambda_2M$};
\draw[-> ,red] (m) -- (ep)
\draw[-> ,red] (p) -- (ep)
\draw (m.north)+(0,0.5) node[above]{mRNA};
\draw (p.north)+(0,0.5) node[above]{Protein};

\draw (m.west) node[below left,inner sep=2] {\tiny +1};
\draw (m.south) node[below left,inner sep=2] {\tiny -1};
\draw (p.west) node[below left,inner sep=2] {\tiny +1};


\end{tikzpicture}}
\begin{centering}
\subfloat[\label{fig:mod1_profile}]{\includegraphics[height=5.5cm,page=1]{model1/Figure_20160629_three_examples_mod1}

}~~\subfloat[\label{fig:mod1_var_ratio}]{\begin{tikzpicture}
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[height=5.5cm,page=1]{model1/Figure_20160713_general_results}};

\node[anchor=south west,inner sep=0] (im1) at (4,.9) {

\includegraphics[width=4.4cm]{model1/Figure_20160704_one_simulation}};

\end{tikzpicture}}
\par\end{centering}
\caption{\label{fig:mod1}Intermediate model with volume growth and partition at division.
\textbf{(A)}:  The model considers only transcription,  translation and  partitioning at division. \textbf{(B)}: Model
of production for one type of protein. The parameters $\lambda_{1}$,
$\sigma_{1}$ and $\lambda_{2}$ are specific to the type of protein.
The gene concentration is considered as constant, giving a rate of
transcription proportional to the cell volume. \textbf{(C)}: The evolution
of the normalized protein concentration during the cell cycle for
three proteins. The central thick line represents the average production
during the cell cycle and the colored area  the standard
deviation. \textbf{(D)} Main: For each type of protein, the variance
in the case of exact partition divided by the variance in the case
of random partition as a function of the Fano Factor. The proportion
of variance attributed to the effect of the random partitioning globally
follows the prediction of the simplified model in blue dashed line
(see Section~\ref{subsec:S_mod1_simple_model_partition}). Inset:
Comparison between the two partitions for the protein OmpC. The thick
line is the average protein concentration during the cell cycle and
the colored areas correspond to the standard deviation in the two
cases. }
\end{figure}
The first intermediate model focuses on the effect of random partitioning. Fluctuations in the availability of RNA-polymerases and ribosomes, as well as the gene replication is therefore not considered. As the global sharing of ribosomes and RNA-polymerases are not represented, it results in a model where there are no interactions between the productions of different proteins; each productions can be considered separately. It is a gene-centered model and one can focus on the production of one particular type of protein; $M$ and $P$ being respectively the number of mRNAs and proteins of this type.

Figure~\ref{fig:mod_paulsson_scheme} depicts this intermediate model. The parameter $\lambda_{1}$ of transcription (resp. $\lambda_{2}$ for the translation) implicitly includes some aspects specific to the gene (the promoter\textendash polymerase affinity, for instance) and the effective constant concentration of free RNA-polymerases (resp. free ribosomes). As the gene is in constant concentration, the rate of mRNA creation increases alongside the volume; it is therefore equal to $\lambda_{1}V(s)$, with $V(s)$ the volume of the cell at time $s$. Hence, the rate of mRNA production per volume unit remains constant.

For this intermediate model, a deterministic volume growth is considered.  Based on experimental results (see \cite{wang_robust_2010}), we represent the growth of the bacteria volume as exponential. For this model, if $s$ is the time spent since the last division, the volume is given by
\[
V(s)=V_{0}2^{s/\tau_{D}},
\]
with $V_{0}$ being the typical size of a cell at birth and $\tau_{D}$ the time of the cell cycle. The explicit description of the volume allows to consider the concentration of mRNA $M(s)/V(s)$ and proteins $P(s)/V(s)$ at any time $s$ of the cell cycle.

In this first model, we study the effect of  partitioning on the protein variance. Two mechanisms of segregating compounds at division are compared, either an exact or a random partitioning.
\begin{itemize}
\item For the exact partitioning,  at division the number of proteins and mRNAs are equally allocated between the two daughter cells. Clearly, this mechanism does not have an impact on the variance of the mRNA and protein concentrations.
\item For the random partitioning at division, each mRNA and protein has an equal chance to be in either one of the two daughter cells (so that with probability $1/2$, they are in the next cell of interest). An additional variability,  due to this random allocation, should be therefore  added in this case.
\end{itemize}
As depicted in Materials and Methods, we first perform a theoretical analysis of this model in order to predict its mean concentrations of mRNA and protein averaged over the cell cycle: respectively $\Ehat{M/V}$ and $\Ehat{P/V}$ as they are defined in Equations~\ref{eq:Ehat_P} of the Materials and Methods section. In this case, for any time $s$ of the cell cycle, the mean concentration of mRNAs and proteins remain constant during the cell cycle and are
given by,
\begin{equation}
\Ehat{M/V}=\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}\qquad\text{and}\qquad\Ehat{P/V}=\frac{\lambda_{2}\tau_{D}}{\log2}\cdot\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}\label{eq:mod1_EM_EP}
\end{equation}
Proofs of these formulas can be found in Section~\ref{subsec:S_mod1_mrna_dynamic} and \ref{subsec:S_mod1_prot_dynamic} of the Appendix. These formulas allow, for every gene considered in \cite{taniguchi_quantifying_2010}, to fix the parameters $\lambda_{1}$, $\sigma_{1}$ and $\lambda_{2}$ in order to have an average production that corresponds to those experimentally measured. This gives a series of parameters corresponding to a representative sample of real bacterial genes (more details in Section~\ref{subsec:S_mod1_param} of the Appendix). These parameters, for every genes, are summarized in the box-plots of Figure~\ref{tab:Quantitative-summary-Model0-1}. As described in Section~\ref{subsec:S_mod1_Gillespie-Algorithms} of the Appendix, simulations are performed using an algorithm derived from Gillespie method in order to determine the evolution of the protein production across the cell cycle and determine its variance averaged over the cell cycle $\Vhat{P/V}$ -- as it is defined in Equation~\ref{eq:Vhat_P} of Materials and Methods.  We then check that the average protein and mRNA concentrations correspond to those experimentally measured (see Figure~\ref{fig:S_mod1_mRNA-prot-check}).

Figure~\ref{fig:mod1_profile} shows the profile of protein concentration during the cell cycle for three representative genes (\emph{adk}, \emph{fabH} and \emph{yjiE}); which are respectively highly, moderately and lowly expressed. As predicted theoretically, the mean concentration $\E{P_{s}/V(s)}$ does not change across the cell cycle (it is due to the fact that the gene concentration remains constant). The figure shows the global tendency that the CV (defined as $\sqrt{\Vhat{P/V}}/\E{P/V}$) decreases as the average production increases.

Figure~\ref{fig:mod1_var_ratio} shows the proportion of variance that is added with the introduction of the random partitioning.  In appears that for all genes, their protein variance indeed increase with random partition, up to be doubled for some genes (like for the protein OmpC whose profile is shown in the inset, where the average production remains constant in both cases, but the random partition increases the protein variance at the beginning of the cell cycle.).

The x-axis of Figure~\ref{fig:mod1_var_ratio} is somewhat unusual as it is the protein Fano Factor, defined as $\Vhat{P/V}/\Ehat{P/V}$.  It is used because the proportion of variance added by the random partition shows a remarkable clear dependence to the Fano Factor: proteins with a low Fano Factor are particularly more impacted.  Note that rare  proteins also tend to have low Fano factor (see Figure~\ref{fig:S_mod1-FanoFactor}), the global tendency remains the same with having the average production as an x-axis, even if this dependence is less strong (see Figure~\ref{fig:S_mod1-VarRatioVsMean}).

In order to explain this clear dependence on the Fano factor, a simplified model, that only focuses on the partition effect without considering the volume growth, has been analyzed (details in Section~\ref{subsec:S_mod1_simple_model_partition} of the Appendix). Its predictions are shown in blue dotted line of the main Figure~\ref{fig:mod1_var_ratio}. It globally predicts the proportion of noise that can be attributed to the random partitioning. It confirms that this effect is only significant for proteins with very low Fano factor. It also shows that it is unlikely that the variance observed for a given gene  can mainly be attributed to random partitioning.



Analogously to the dual reporter technique, we can decompose, using
the environmental state decomposition, the variance of each protein
$\Vhat{P/V}$ into the two terms $\Vhatt{P/V}{int}$ and $\Vhatt{P/V}{ext}$
that would respectively represent the intrinsic and extrinsic contribution
to the protein variance (see Subsection~\ref{methods:Intrinsic extrinsic} of Materials and Methods). In this model, as the only external environment
considered is the cell cycle, the variance of each protein $\Vhat{P/V}$
can be decomposed into the sum of
\begin{align}
\Vhatt{P/V}{int} & =  \frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\V{P_{s}/V(s)}\mathop{}\mathopen{}\mathrm{d} s,\label{eq:var_1}\\
\Vhatt{P/V}{ext} & =  \frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\E{P_{s}/V(s)}^{2}\mathop{}\mathopen{}\mathrm{d} s-\Ehat{P/V}^{2}.\label{eq:var_2}
\end{align}
(see Section~\ref{subsec:S_mod1_ESD} of~the Appendix).  In this intermediate model, even with random partitioning, such a decomposition shows surprisingly no external contribution $\Vhatt{P/V}{ext}$ for every protein (since the protein concentration $P_{s}/V(s)$ remains constant across the cell cycle). It is therefore remarkable that this decomposition only captures a part of what is generally accepted as the extrinsic noise.

\subsection{\label{subsec:model2} Impact of Gene Replication}

\begin{figure}
\subfloat[\label{fig:mod2_principle}]{
\includegraphics[width=6.5cm]{model2/Figure_20160721_BioModel_ultime21}}~~~\subfloat[\label{fig:mod2_scheme}]{
\begin{tikzpicture}
\linespread{1}
\matrix [column sep=13mm, row sep=10mm,ampersand replacement=\&] {
  \node (g) [draw=none,fill=none] {}; \&
  \node (m) [draw, shape=rectangle] {$M$}; \&
  \node (p) [draw, shape=rectangle] {$P$}; \\
  \&
  \node (em) [draw=none,fill=none] {$\emptyset$}; \&
  \node (ep) [draw=none,fill=none,color=red,align=center] {
  Random partitioning \\  every $\tau_{D}$}; \\
  };
\draw[-> ] (g) -- (m) node[pos=-0.1,above,align=left]
		{Before repl. $\lambda_1$ };
\draw[-> ] (g) -- (m) node[pos=-0.1,below,align=left]
		{After repl. $2\lambda_1$ };
\draw[-> ] (m) -- (em) node[right,midway] {$\sigma_1M$};
\draw[-> ] (m) -- (p) node[midway,above] {$\lambda_2M$};
\draw[-> ,red] (m) -- (ep)
\draw[-> ,red] (p) -- (ep)
\draw (m.north)+(0,0.5) node[above]{mRNA};
\draw (p.north)+(0,0.5) node[above]{Protein};

\draw (m.west) node[below left,inner sep=2] {\tiny +1};
\draw (m.south) node[below left,inner sep=2] {\tiny -1};
\draw (p.west) node[below left,inner sep=2] {\tiny +1};


\end{tikzpicture}
}

\subfloat[\label{fig:mod2_profiles}]{\includegraphics[height=5.5cm]{model2/Figure_20160629_three_examples_mod2}
}~~\subfloat[\label{fig:mod2_va_ratio_and_profile_AdK}]{\begin{tikzpicture}
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[height=5.5cm,page=1]{model2/Figure_20160713_general_results}};

\node[anchor=south west,inner sep=0] (im1) at (1.3,.9) {

\includegraphics[height=3.3cm]{model2/Figure_20160819_variation_on_one_gene}};

\end{tikzpicture}}
\caption{\label{fig:mod2}Intermediate model with cell cycle and gene replication.
\textbf{(A)}: Features of the intermediate model. The model now takes
into account replication. \textbf{(B)}: The production of one particular
type of protein. The number of mRNAs and proteins are respectively
$M$ and $P$; the difference with the previous model is the introduction
of replication at a time $\tau_{R}$ after the cell birth where the
transcription rate is doubled. \textbf{(C)}: The evolution of the
normalized protein concentration during the cell cycle for three proteins.
The central thick line and the colored areas represent the same quantities
as in Figure~\ref{fig:mod1_profile}. \textbf{(D)} Main: for each
type of protein, the protein variance of this model divided by the
one in the previous model (gene in constant concentration and random
partitioning). The replication paradoxically tends to slightly diminish
the fluctuations of the concentration. Insight: on the left, the concentration
through the cell cycle for the protein Adk (a close-up of the one presented
in Figure~\ref{fig:mod2_profiles}); \rmodif{on the right, profile of a modified version of Adk with parameters chosen in order to minimize the variance.}}
\end{figure}

Taking back the model previously described with volume growth and random partitioning at division, we introduce the notion of gene replication. As in the case of the slowing growing bacteria of \cite{taniguchi_quantifying_2010}, we consider only one DNA replication per cell cycle. The gene is represented for now on as an entity that is replicated at some instant $\tau_{R}$ of the cell cycle, hence doubling the transcription rate (see Figure~\ref{fig:mod2_scheme}): before a time $\tau_{R}$ after the cell birth, the rate of transcription will be $\lambda_{1}$, after the time $\tau_{R}$, the rate of transcription is doubled. The time $\tau_{R}$ at which the gene is replicated is specific to each gene, and is determined according to the position of the gene on the chromosome.

The mean of mRNAs concentration can be determined at any moment of the cell cycle. At any time $s$ in the cell cycle (with $0\leq s<\tau_{D}$ ), it is given by,
\begin{equation}
\E{M_{s}/V(s)}=\frac{\lambda_{1}}{\sigma_{1}V(s)}\left[1-\frac{e^{-(s+\tau_{D}-\tau_{R})\sigma_{1}}}{2-e^{-\tau_{D}\sigma_{1}}}+\ind{s\geq\tau_{R}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right],\label{eq:mod2-EMs}
\end{equation}
with $\ind{x\geq y}=1$ if $x\geq y$ and $0$ otherwise. One can refer to Section~\ref{subsec:S_mod2_mRNA-number} for the proof.

Similarly, the mean of  protein concentration can also be explicitly determined. Before the replication, if $s$ indicates the time after the birth of the cell (i.e. $0\leq s<\tau_{R}$ ), we can determine the mean of $P_{s}/V(s)$ as a function of the mean of $M_{0}$ and $P_{0}$. Similarly, after the replication, the mean of $P_{s+\tau_{R}}/V(s+\tau_{R})$ (with $s$ such as $0\leq s<\tau_{D}-\tau_{R}$ ) is known as a function of the mean of $M_{\tau_{R}}$ and $P_{\tau_{R}}$.  They are given by the formula,
\begin{equation}
\E{P_{s+\tau}/V(s+\tau)}=\frac{1}{V(s+\tau)}\left[\E{P_{\tau}}+\frac{\eta\lambda_{1}}{\sigma_{1}}\lambda_{2}s+\left(\E{M_{\tau}}-\frac{\eta\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right],\label{eq:mod2-EPs}
\end{equation}
with $\tau=0$ and $\eta=1$ for the case before the replication
and $\tau=\tau_{R}$ and $\eta=2$ for the case after the replication.
At equilibrium we are able to determine explicit values for $\E{P_{0}}$
and $\E{P_{\tau_{R}}}$ as a function of the parameters (see Section~\ref{subsec:S_mod2_Protein-number}).

As for the previous model, Equations~(\ref{eq:mod2-EMs}) and~(\ref{eq:mod2-EPs}) allow, for every gene experimentally measured in \cite{taniguchi_quantifying_2010}, to give a set of parameters $\lambda_{1}$, $\sigma_{1}$ and $\lambda_{2}$ that corresponds to it (Section~\ref{subsec:S_mod2_param} and Figure~\ref{fig:S_mode2_TableParam}).  To determine the variance of protein concentration, simulations can be performed, but it is noticeable that we also managed to have  formulas for the mRNA and protein variances. These formulas greatly simplify the analysis of this intermediate model (see Section~\ref{subsec:S_mod2_mRNA-number} and \ref{subsec:S_mod2_Protein-number}).

In the left figure of the inset of Figure~\ref{fig:mod2_va_ratio_and_profile_AdK} is shown the profile of the protein Adk: its average concentration (thick line) and its standard deviation (blue area) during the cell cycle are shown. It appears that the mean concentration at any given time $s$ of the cell cycle $\E{P_{s}/V(s)}$ is not constant during the cell cycle, as it was the case in the model of Figure~\ref{fig:mod1}.  The curve of $\E{P_{s}/V(s)}$ fluctuates around $2\%$ of the global average protein production $\Ehat{P/V}$. This small fluctuation is barely perceptible for the normalized profiles presented in Figure~\ref{fig:mod2_profiles}: globally this figure is not significantly different from the profiles of the previous intermediate model presented in Figure~\ref{fig:mod1_profile}.

The main Figure~\ref{fig:mod2_va_ratio_and_profile_AdK} shows the effect of replication on the variance: it represents the ratio of protein variance between the current model with gene replication and the previous model with the gene in constant concentration and random partitioning.  For all the genes, the variances predicted show little difference from the previous intermediate model. The ratio is even surprisingly slightly above one for many genes, indicating that for these genes the replication tends to reduce the variance.

As for the previous intermediate model, we can use the environmental state decomposition to separate the part of variance $\Vhatt{P/V}{int}$ (defined in Equation~(\ref{eq:var_1})) specific to the gene expression and $\Vhatt{P/V}{ext}$ (defined in Equation~(\ref{eq:var_2})) attributed to cell cycle fluctuations. As the mean concentration, $\E{P_{s}/V(s)}$ is no longer constant during the cell cycle $\Vhatt{P/V}{ext}$ is no longer null. Yet it appears that $\Vhatt{P/V}{ext}$ only represents a very small part of the global variance $\Vhat{P/V}$ (for $99\%$ of the genes, it represents less than $2\%$). For  this intermediate model, the extrinsic contribution of DNA replication computed with this decomposition is small.

Using the analytical formula of the protein variance (Section~\ref{subsec:S_mod2_Protein-number}) and by performing variations on certain parameters, we can analyze the effect of several aspects on the protein variance in this model. By  keeping the average protein production constant, we modify the position of the gene on the chromosome, or the mRNA lifetime or change the average mRNA number. Results are shown in Figure~\ref{fig:S_mod2_variation-gene}.  Changing the gene position have almost no impact on the global protein variance $\Vhat{P/V}$. The effect of mRNA lifetime is more noticeable as a shorter mRNA lifetime can diminish the protein variance at most about $40\%$.

The mRNA number seems to have the most important effect on the protein production: for the same average protein production, having more mRNAs greatly diminish the protein variance; such effect has been experimentally observed \citep{blake_noise_2003,ozbudak_regulation_2002}. The right insight of Figure~\ref{fig:mod2_va_ratio_and_profile_AdK} shows an example of such a protein with reduced variance: this protein is based on Adk, the protein average production is the same but there are ten times more mRNAs, with a ten times storter lifetime. The variance in indeed reduced but with the cost of the production of additional mRNAs. Yet, even in this case, the protein expression is not strongly cycle-dependent (see inset of Figure~\ref{fig:S_mod2_recovering_distrib}); in particular, this profile is not precise enough to be used as a ``trigger'' for periodic cell events (such as DNA replication initiation, or partition at division): the evolution of the protein concentration across the cell cycle is not precise enough to robustly distinguish different phases of the cell cycle.

\subsection{Impact of the sharing of RNA-polymerases and ribosomes }\label{subsec:model3}

\begin{figure}
\subfloat[\label{fig:mod3_principle}]{\centering{}\includegraphics[height=4.cm]{model3/Figure_20160721_BioModel_ultime32}}~~~\subfloat[\label{fig:mod3_scheme}]{\begin{tikzpicture}

\matrix [column sep=6mm, row sep=3.mm,ampersand replacement=\&] {
  \&
  \node (fy) [draw, shape=rectangle] {$F_Y$}; \&
  \&
  \node (fr) [draw, shape=rectangle] {$F_R$}; \&
  \\

  \node (g) [draw, shape=rectangle] {$G_i$}; \&
  \node (ey) [draw, shape=rectangle] {$E_{Y,i}$}; \&
  \node (m) [draw, shape=rectangle] {$M_i$}; \&
  \node (er) [draw, shape=rectangle] {$E_{R,i}$}; \&
  \node (p) [draw, shape=rectangle] {$P_i$};\\

  \&
  \node (none1) [draw=none,fill=none] {}; \&
  \node (empty) [draw=none,fill=none] {$\emptyset$}; \&
  \&
  \node (none3) [draw=none,fill=none,align=left] {\phantom{P}\\\phantom{P}}; \\

  \node (gj) [draw, shape=rectangle] {$G_j$}; \&
  \node (eyj) [draw, shape=rectangle] {$E_{Y,j}$}; \&
  \node (mj) [draw, shape=rectangle] {$M_j$}; \&
  \node (erj) [draw, shape=rectangle] {$E_{R,j}$}; \&
  \node (pj) [draw, shape=rectangle] {$P_j$}; \\
  \&
  \node (none1j) [draw=none,fill=none] {}; \&
  \node (emptyj) [draw=none,fill=none] {$\emptyset$}; \&
  \&
  \node (none3j) [draw=none,fill=none,align=left] {\phantom{P}\\\phantom{P}}; \\
};
\draw[->] (g) -- (ey);
  \draw[->] (g.east) -- +(0.5,0) |- (none1.center)  -- +(0.6,0) |- (m.west);
  \draw (none1.center) node[above] {};
\draw[->] (fy.west)  to [left,bend right=45]  (ey.west);
\draw[->] (ey.east)  to [left,bend right=45]  (fy.east);
\draw[] (ey.east)+(0.7,1.2) node[right] {};
  \draw (fy.east)+(0.3,0)  node[right,above] {\tiny +1};
  \draw (fy.west)+(-0.3,0) node[left,above]  {\tiny -1};
\draw[-> ] (m) -- (empty) node[midway,right] {};
\draw[->] (m) -- (er) node[below,midway] {};
\draw[->] (er) -- (p) node[below,midway] {};
\draw[->] (fr.west)  to [left,bend right=45]  (er.west);
\draw[->] (er.east)  to [left,bend right=45]  (fr.east);
  \draw (fr.east)+(0.3,0)  node[right,above] {\tiny +1};
  \draw (fr.west)+(-0.3,0) node[left,above]  {\tiny -1};
\draw[->,color=red] (m)  --  (none3.west);
\draw[->,color=red] (p.south)  --  (none3.north);
\draw (none3)+(.5,0) node[draw=none,fill=none,color=red,align=left] {Random \\partitioning};

\draw[->] (gj) -- (eyj);
 \draw[->] (gj.east) -- +(0.5,0) |- (none1j.center)  -- +(0.6,0) |- (mj.west);
  \draw (none1j.center) node[above] {};
\draw[->] (eyj.east)  to [left,bend right=45]  (fy.east);
\draw[->] (fy.west)  to [left ,bend right=45]  (eyj.west);
\draw[] (eyj.east)+(0.7,1.2) node[right] {};
\draw[-> ] (mj) -- (emptyj) node[midway,right] {};
\draw[->] (mj) -- (erj) node[below,midway] {};
\draw[->] (erj) -- (pj) node[below,midway] {};
\draw[->] (fr.west)  to [left ,bend right=45]  (erj.west);
\draw[->] (erj.east)  to [left,bend right=45]  (fr.east);
\draw[->,color=red] (mj)  --  (none3j.west);
\draw[->,color=red] (pj.south)  --  (none3j.north);
\draw (none3j)+(.5,0) node[draw=none,fill=none,color=red,align=left] {Random \\partitioning};
\end{tikzpicture}}

\subfloat[\label{fig:mod3_profiles}]{\includegraphics[height=6cm]{model3/Figure_20160912_MCDivisionThreeExamples_vs_MSGeneReplication}

}~~\subfloat[\label{fig:mod3_var_ratio}]{\begin{tikzpicture}
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[height=6cm]{model3/Figure_20160912_MCDivision_vs_MSGeneReplicationSI}};

\node[anchor=center,inner sep=0] (im1) at (4.,2.4) {

\includegraphics[height=3cm,page=2]{model3/Figure_20160912_MCDivisionThreeExamples_vs_MSGeneReplication}};

\end{tikzpicture}}
\caption{\label{fig:mod3}Complete model. (A): The model now considers the sharing of RNA-polymerases and ribosomes between the different productions. (B): The model of production of all proteins. The $i$-th gene is associated with one particular type of mRNA (whose number is $M_{i}$) and proteins (whose number is $P_{i}$).  The number of free RNA-polymerases (resp. ribosomes) is $F_{Y}$ (resp.  $F_{R}$), the number of those sequestered on the $i$-th gene is $E_{Y,i}$ (resp. $E_{R,i}$). (C): The evolution of the normalized protein concentration during the cell cycle for three proteins. The central thick line and the colored areas represent the same quantities as in Figure~\ref{fig:mod1_profile}. (D) Main: variance in the model with sharing of RNA-polymerases and ribosomes divided with the model with replication. Insight: one the left, the concentration through the cell cycle for the protein FabH (a close-up of the one presented in Figure~\ref{fig:mod3_profiles})}
\end{figure}

We now consider the complete model as it is  described in details in Subsection~\ref{method:complete model} and that includes a limited amount of RNA-polymerases and ribosomes. Now, RNA-polymerases and ribosomes are explicitly represented in the model and each of these macromolecules is considered either allocated (i.e. sequestered on a gene if it is an RNA-polymerase, or on an mRNA if it is a ribosome), or free (i.e. either moving freely in the cytoplasm or, in the case of RNA-polymerases, non-specifically sliding on the DNA).

The previous intermediate models were ``gene-centered'', each class of proteins was considered independently from each other.  The common sharing of RNA-polymerases and ribosomes is an additional key feature that leads us to investigate a multi-protein model where all the genes of the bacteria are considered altogether. For each type of protein $i$, we denote by $G_{i}(s)$, $M_{i}(s)$ and $P_{i}(s)$ respectively the number of gene copies, of messengers and of proteins at time $s$. For each gene $i$, denote by $E_{Y,i}(s)$ the number of RNA-polymerases sequestered on the $i$-th gene for transcription and by $E_{R,i}(s)$, the number of ribosomes sequestered on a mRNA of type $i$ for translation. The non-allocated RNA-polymerases and ribosomes are respectively denoted by $F_{Y}(s)$ and $F_{R}(s)$.  In a first step, we have considered that the gene pool of \cite{taniguchi_quantifying_2010} ($841$ genes with their mRNA and protein expression measured) would represent the whole genome. We will later see that the addition of new genes does not change significantly  our results (see Section~\ref{subsec:S_Other-Additional-genes} of the Appendix).

New ribosomes and RNA-polymerases are added to the system as the cell volume increases: in a first step, these macro-molecules are added regularly such as their concentration in the cell remains constant during the cell cycle; we will later consider a more realistic way to represent RNA-polymerases and ribosomes production.

The previous intermediate models represented the production of a specific type of protein immersed into a ``background environment'' where the cell grows and divides, the current model includes simultaneously all the genes altogether. In this model, as we are on the scale of the whole cell, we would like to model the impact of the global protein production on the cell growth. We therefore can no longer consider that the production of each type of protein has no effect on the global performance of the cell.
The volume $V(s)$ depends now on the global production of proteins, and it is not an independent and deterministic feature anymore. As the density of cell components tends to be constant in real-life experiment \citep{marr_growth_1991} and as the dry mass of the cell is essentially proteins \citep{neidhardt_chemical_1996}, the model considers the volume as proportional to the current total mass of proteins in the cell. The mass of each protein is given by the length of its peptide chain (see Section~\ref{subsec:S_mod3_DeterModel-pres} for an exhaustive description of the model).

The processes of mRNA and protein productions are both separated in two parts: the binding and initiation on one side, and the elongation and termination on the other side. The rate at which an RNA-polymerase fixes on a gene of type $i$ at time $s$ depends on the copy number of the $i$-th gene $G_{i}(s)$, the free RNA-polymerase concentration $F_{Y}(s)/V(s)$ and a parameter $\lambda_{1,i}$ specific to the gene that takes into account the RNA-polymerase\textendash promoter affinity. The elongation rate of each mRNA only depends on the average transcription speed and the length of the gene. The mechanism for translation is similar.

As this model is more complex than the previous ones, the complete analytical description of mRNA and protein dynamics seems to be out of reach. For this reason we have used a simplified model to represent the average behavior of the stochastic and used its predictions to fix the parameters (see Section~\ref{subsec:S_mod3_DeterModel-pres}).



At the cell scale, simulations describe exponentially growing cells, see Figure~\ref{fig:S_mod3_exp-growth}. The amount of available ribosomes and RNA-polymerases also changes rapidly, of the order of the second for the RNA-polymerases, and of the order of one tenth of a second for the ribosomoes, Figures~\ref{fig:S_mod3_Frnap-profile} and~\ref{fig:S_mod3_Fribo-profile}.Taniguchi et al.


Figure~\ref{fig:mod3_var_ratio} compares the results of the simulations with the previous intermediate model with gene replication and random partitioning. It shows that, for $90\%$ of the genes, the interactions between protein productions only represent at most $10\%$ of variability.  In inset is shown the example of the protein FabH profile during the cell cycle, showing that sharing of RNA-polymerases and ribosomes introduces little change.

We also analyze the model using the environmental state decomposition.  Two genes of the dual reporter technique would undergo the same volume growth with the same evolution of free RNA-polymerases and ribosomes; as a consequence, the common cellular environment on which the decomposition is operated now includes the concentration of free RNA-polymerases and ribosomes (in addition to the cell cycle). Details can be found in Section~\ref{subsec:S_mod3_ESD} of the Appendix.  It appears that for all the genes, the extrinsic contribution of the variance $\Vhatt{P_{i}/V}{ext}$ represents only a very small portion of $\Vhat{P_{i}/V}$ (for $99\%$ of the cells, the ratio $\Vhatt{P_{i}/V}{ext}/\Vhat{P_{i}/V}$ represents less than $1\%$).

We compared these simulation results with a simplified theoretical model: Section~\ref{sec:S_mod3_Other-influences} presents a multi-protein model, that is inspired by the one described in \cite{fromion_stochastic_2015}.  Even if it is a multi-protein model as it represents the expression of a large number of genes altogether, it is a simpler model than the one presented here as it considers separately transcription and translation, and it does not consider neither volume growth, partitioning at division, nor DNA replication. We show that the predicted distributions of free RNA-polymerases (and ribosomes in the adapted model) fits well the one observed in our simulations (see Figure~\ref{fig:S_mod3_influence_fRNAP_fRibo}).  As we will see in the discussion, this good correspondence between the models would suggest that the mean-field mathematical properties proven for the simplified model could be applied to our complete model.

\subsection{Model and Parameter Sensitivity}

The complete model supposes a series of modeling and parameter choices that might legitimately influence the protein production. We have analyzed several of these aspects and have  shown that they do not appear to  significantly change the results previously presented.
\begin{description}
\item [{Quantity of free RNA-polymerases and ribosomes}] The average concentrations of free RNA-polymerases and ribosomes in a cell cannot be deduced from Taniguchi~et al. They are nevertheless needed to estimate the parameters of our model, see Section~\ref{subsec:S_mod3_param}. Globally, one can expect to have a low concentration of free ribosomes and a higher concentration of free RNA-polymerases, see Section~\ref{subsec:F-RNAP-and_ribo}. But precise numbers  seem to be difficult to obtain.  We therefore perform several simulations different values for these concentrations (for each macro-molecule, a concentration taking successively 1, 10, 100 and 1000~$\text{copies}/\text{\textmu m}^3$), without significant changes. See Section~\ref{sec:Impact-RNAP-ribo} for details.

\item [{Additional Genes}] As previously said, to perform our first simulations, only 841 genes from which the average mRNA and protein concentration have been measured in Tanuguchi~et~al. are considered. To have a pool of proteins that might represent a global gene expression in \emph{E. coli}, we studied the case of a simulation with a  set of genes that represent about 2000 genes, more in  accordance to the expected number of genes expressed in a growing condition. To propose realistic parameters for these fictional genes, we sample them according to different empirical distributions observed in the empirical data, and also by taking into account the possible correlations observed (the correlation that exists between the average mRNA and protein production for instance). Details can be found in Section~\ref{subsec:S_Other-Additional-genes}.

\item [{Non-specific binding of RNA-polymerases}] It has been proposed that many of the RNA-polymerases are non-specifically bound on the DNA (see \cite{klumpp_growth-rate-dependent_2008} for instance). We have done a simulation where RNA-polymerase can bind non-specifically on the DNA. When in this state, they are not available for the transcription. See Section~\ref{subsec:Other-Sequestered-polymerases}.

\item [{Production of RNA-polymerases and ribosomes}] The total amount of RNA-polymerases and ribosomes (whether free or not) were at first considered in constant concentration: the RNA-polymerases and ribosomes were added as the cell volume increases. We have done a simulation that considers a way to represent their production to have a more realistic representation: both RNA-polymerases and ribosomes are produced as if they were one of the proteins of the system (this goal of the simulation is just to have an insight of the effect of RNA-polymerase and ribosome production, not to represent precisely their production mechanisms). \rmodif{The introduction of such mechanism indeed changes some aspects of the simulation: in particular, the growth of the cell is more erratic at it then directly correlated on the total number of ribosomes (see Figure~\ref{fig:production-RNAP-ribo}). But as the production of proteins increase with a higher number of ribosomes, so does the volume of the cell. In terms, of concentration, the induced fluctuations in the number of ribosomes has little   impact in terms of protein concentration variance. See Section~\ref{subsec:S_Other-Production-of-RNAP-ribo} for more details.}


\item [{Precision in the division and DNA replication initiation timing}] The initial simulation triggers DNA replication initiation and division when the cell reaches some volume. We investigate approximate models of division and of DNA replication initiation by introducing a volume-dependent rate of division as it is commonly used in the literature, see \cite{tyson_sloppy_1986,wang_robust_2010,soifer_single-cell_2014,osella_concerted_2014}. See Section~\ref{subsec:Uncertainty-in-divi-init}.



\end{description}


\section{Discussion}

\subsection{Interpretation of the model predictions}

The experimental data of Taniguchi et al. gave us the opportunity to systematically and quantitatively inspect the impact on protein variance of many cell mechanisms that are not often considered in stochastic models of protein production. The broad variety of genes experimentally measured in Taniguchi et al. has been a good opportunity for us to realistically test our models for a wide number of different genes, with different mRNA and protein productions, different mRNA lifetimes or gene position on the chromosome or gene length.

From this analysis, it appears that among the different features included in the model, the random partitioning has the most significant effect on the variance of protein concentration, especially for the less expressed proteins.  The gene replication induces little difference (it even tends, to a small extent, to reduce the variance in some cases); the important fluctuations of free RNA-polymerases and ribosomes have little impact on the protein production, which does not fit the hypothesis made in \cite{taniguchi_quantifying_2010}. It is confirmed by the environmental state decomposition, that separates the intrinsic and extrinsic contribution on protein variance (in a analog way as it is done with the dual reporter technique): the extrinsic contributions represent at most a few percents of the total variance.

We interpret the surprising little impact of sharing of RNA-polymerases and ribosomes on the proteins variance by noticing the similarities of our model with the one described in \cite{fromion_stochastic_2015}. Indeed, as previously explained, the global behavior of free RNA-polymerases and ribosomes can be predicted by a simplified model derive from \cite{fromion_stochastic_2015}, where the RNA-polymerases and ribosomes are shared among the different productions. The main result of this article is a mean field theorem: as the number of genes increases, the production process of different types of proteins can be seen as independent production processes. The reason is that the dynamic of free RNA-polymerases and ribosomes is much faster than the production of mRNAs and proteins of one particular type. The rate at which an mRNA or protein is elongated only depends on the ``local equilibrium'' concentrations of free RNA-polymerases and ribosomes (a similar phenomenon can be found in \cite{dessalles_stochastic_2017b}). Our model seems to display such mechanism: with a global sharing of RNA-polymerases and ribosomes by a large number of protein productions, the dynamic of free RNA-polymerases and ribosomes is faster than the production of each mRNA and protein of each type. As a consequence, it is not surprising to see that this multi-protein model, which takes into account the production of all proteins display little difference with the intermediate ``gene-centered'' model due to a mean-field effect.

\subsection{Comparison with experimental measures}

In the end, we can compare the results of our models with the experiments. Firstly, one can remark that the profile of the mean production $\E{P_{s}/V(s)}$ (the plain line in Figure~\ref{fig:mod2_va_ratio_and_profile_AdK}) during the cell cycle corresponds to the one observed experimentally in \cite{walker_generation_2016}. Furthermore, the maximum deviation of the average $\E{P_{s}/V(s)}$ around the global average protein $\Ehat{P/V}$ (red dashed line in Figure~\ref{fig:mod2_va_ratio_and_profile_AdK}) is between $2\%$ and $4\%$ for all the proteins of our models, and  Walker et al. measure such fluctuations also around $2\%$ of the global average for genes at different positions on the chromosome (see Figure 1.d and Figure S6.b of \cite{walker_generation_2016}).

Secondly, we compared the global protein fluctuations predicted in our models with those measured in Taniguchi et al.  Figure~\ref{fig:final_pannel} shows, for every gene, the protein CV (defined as $\Vhat{P/V}/\Ehat{P/V}^{2}$) against the average protein concentration $\Ehat{P/V}$; it is compared with the same results obtained by \cite{taniguchi_quantifying_2010} (corresponding to the Figure~2.B of this article). It appears that the noise approximately scales inversely the average protein concentration like in the first ``intrinsic noise'' regime of \cite{taniguchi_quantifying_2010}.  But unlike in the \cite{taniguchi_quantifying_2010} experiment, there is no lower plateau for highly expressed proteins: for the highest produced proteins, the CV should be in the order of $10^3$ fold higher than the one predicted. It confirms that the features considered here cannot correctly explain the noise observed experimentally.

\begin{figure}
\subfloat[\label{fig:production-RNAP-ribo}]{\includegraphics[page=1,height=5.2cm]{model3/Figure_20160802_elong_example}}
\subfloat[\label{fig:final_pannel}]{\includegraphics[height=5.2cm]{model3/Figure_20160912_MCDivisionGeneral}}

\caption{(A): A variation on the complete model where the RNA-polymerases and ribosomes are produced by the cell. There is a clear correlation between the total number of ribosomes and the volume derivative (the time step $\Delta t = 30\,\text{sec}$). (B): The coefficient of variation of protein concentration
(defined as $\protect\Vhat{P/V}/\protect\Ehat{P/V}$) for each genes
predicted by the complete model. In yellow the point cloud predicted
by the experimental measures of \cite{taniguchi_quantifying_2010}.}
\end{figure}

We can propose different interpretations to explain the discrepancy between the predictions of the models and the experimental measures. For the  biological processes not included in our models and that can have an impact on the variability, one can first mention the gene regulation as in our models, all the genes are considered as constitutive. The introduction of a gene regulation can indeed induce the a large variability on the protein production \citep{paulsson_models_2005,shahrezaei_analytical_2008}. Nonetheless, the ``extrinsic noise plateau'' observed in Taniguchi et al. only concerns the proteins and not the mRNAs (the mRNAs CV is globally inversely proportional to the mRNA average production, for every types of mRNAs, there is no ``plateau'' for highly expressed mRNAs). As a consequence, one can expect that the mechanism explaining the extrinsic noise plateau takes place at the translation level and not at the transcription. Moreover for highly expressed proteins, the protein CV is independent of the protein expression; it is therefore not gene-dependent as it would be the case for gene regulation. Finally, we have considered a simple model with gene regulation (like the three-stage model of \cite{paulsson_models_2005}), and determined the regulations parameters in order to predict the protein variance observed in \cite{taniguchi_quantifying_2010}; we came with an activation/deactivation timescale has to be very high (several times the doubling time) in order to reproduce the variance experimentally observed, which is way above the typical biologically expected parameters.

One can also mention other possible mechanisms not represented in our models such as the fluctuations of availability of amino-acids or free RNA nucleotides in the medium, thus inducing additional fluctuations in the translation speed. Even if one can see here a clear analogy with the fluctuations in RNA-polymerases and ribosomes availability (which also impact the transcription and translation speeds), the different timescales of the dynamics of amino-acids or free RNA nucleotides might result in a different effect. One can also challenge the exponential hypothesis by proposing different times for some process: for instance, elongation time would be better represented as having an Erlang distribution, that has the sum of independent exponential random variables; even if some results incline to say that it has a limited impact \citep{fromion_stochastic_2013}. Also, in this model, the binding and initiation (either of RNA-polymerases or ribosomes) are considered as a single event. A more precise representation would be to describe them as two different processes (\cite{siwiak_transimulation_2013} gives for instance a median transcriptional initiation time of $15\,\text{s}$ which is of the same order of magnitude as the elongation time).

One can also consider that, as this effect mainly affects proteins with the highest fluorescence, it is possible that some saturation induces a bias in the estimation of variance of highly produced proteins. To our knowledge, exhaustive measures of \cite{taniguchi_quantifying_2010} for the mean and variance of protein and mRNA concentration have not been replicated at the same scale, so we have been able to confront our results to other measures.

\section{\label{methods}Materials and Methods}

\subsection{\label{method:complete model} The complete model}
In this section is presented in details the complete model of Section~\ref{subsec:model3} (the other models being simplifications of this complete model). It is an extension of the classical two-stage model and is presented in Figure~\ref{fig:Model-general}. It represents, for any gene, both the number of mRNAs and proteins associated with the gene inside a given cell. But, contrary to the two-stage model it also explicitly represents the changing volume of the cell $V(s)$ (changing across the time $s$), so that, if $M_s$ and $P_s$ represent the respective number of mRNAs and proteins of a given gene at any time $s$, one can now explicitly represent their concentration by
$$M_s/V(s)\qquad\text{and}\qquad P_s/V(s).$$
Furthermore, contrary to the two-stage model, all genes in the bacteria are represented altogether here (in order to represent the global sharing of RNA-polymerases and ribosomes in the different protein productions) and the division (to represent a partitioning of components at septation). When dividing, the model focuses on only one of the two daughter cells in order to follow one lineage of cells.

\begin{figure}
\begin{centering}
\subfloat[\label{fig:Model-general}]{\includegraphics[width=6cm]{other_figs/Figure_20160721_BioModel_ultime31}}~~~~
\subfloat[\label{fig:procedure}]{
\begin{tikzpicture}[scale=0.7, every node/.style={scale=0.7},ampersand replacement=\&]
\tikzset{
    mynode/.style={rectangle,rounded corners,draw=black, top color=white, bottom color=blue!30,very thick, inner sep=1em, minimum size=3em,text centered,align=center},
    myarrow/.style={->, >=latex', shorten >=1pt, thick},
    mylabel/.style={text width=7em, text centered}
}
\linespread{1}
\def -2.5 {-2.5}
\def 2.8 {2.8}

\node[mynode] (model) {Model};
\node[mynode] (means) at ($(model)+(0,-2.5)$) {Theoretical means\\$\Ehat{P/V}$ and $\Ehat{M/V}$};
\node[mynode] (param) at ($(means)+(0,-2.5)$) {Parameters\\$\lambda_1$, $\sigma_1$, $\lambda_2$, etc.};
\node[mynode] (sim_mean) at ($(param)+(2.8,-2)$){Simulated means\\$\Ehat{P/V}$ and $\Ehat{M/V}$};
\node[mynode,draw=red] (sim_var) at ($(param)+(-2.8,-2)$) {Simulated variance\\$\Vhat{P/V}$};
\node[mynode,bottom color=green!30] (exper) at ($(model)+(2.8,0)$) {Experimental \\measures};


\draw[->] (model) -- (means) node[midway,left,align=right]{theoretical \\analysis};
\draw[->] (means) -- (param) node[midway,left,align=right]{parameters \\estimation};
\draw[->] (param.south) |- (sim_mean.west);
\draw[->] (param.south) |- (sim_var.east) node[midway,below]{simulations};
\draw[->] (sim_mean.north)+(0.3,0) -- ($(exper.south)+(0.3,0)$)
	node[midway,above,right,align=left]{verifying \\ correspondence};
\draw[->] (exper.south)+(-0.3,0) -- +(-0.3,-3) -| (param.north);


\end{tikzpicture}}


\subfloat[\label{fig:mrna_prot_taniguchi}]{
\includegraphics[height=4.75cm]{other_figs/Figure_20160910_sumary_gene_creation}}~~
\subfloat[\label{tab:Quantitative-summary-Model0-1}]{\includegraphics[height=4.75cm]{model1/Figure_20171110_boxParameters}}
\end{centering}
\caption{Experimental data and modeling principles. (A): The complete model. (B): Scheme of analysis to determine the variance of protein concentration for every gene predicted in each intermediate models.  (C): \cite{taniguchi_quantifying_2010} measures of mRNAs and proteins for $841$ genes. (E): Box-plots presenting the parameters deduced from experiments to corresponds to the model of Figure~\ref{fig:mod1}; for the other models, these rates are in the same order of magnitude. }
\end{figure}

\begin{description}
\item [{Volume~growth~and~division}] The volume of the cell is represented in the model and increases alongside the growth of the cell. As the number of mRNAs and proteins of each type is represented, this volume also allows to explicitly represent their concentration inside the cell. When the cell approximately doubles its volume, division occurs and all the compounds (mRNAs and proteins) are randomly partitioned in the two daughter cells (this partition is considered as equally likely  as each compound has an equal chance to be in either one of the two daughter cells, we do not consider strong asymmetry in the cell volume division.
Then the model only follows one of the two daughter cells beginning a new cell cycle.

\item [{Units~of~Production}] Each type of protein has a specific type of mRNA and a unique gene associated with (in particular, there is no notion of operons). For each unit of production, the number of gene copies, mRNAs and proteins inside the cell is explicitly represented. Each copy of the gene can be transcribed into an mRNA.  The mRNA can then be translated into a protein until its degradation, the degradation rate is specific to the type of mRNA. We do not consider any rate of protein disappearance.  The proteolysis occurring in a timescale much longer than the cell cycle (see \cite{koch_protein_1955}) for most proteins, its decay is then dominated by protein partition that occurs at division.
\item [{DNA-Replication}] Each gene can be present in one or two copies in the cell (only one DNA replication is considered as in the slowly growing cells of Taniguchi et al.).  The instant of replication of each gene is simply determined by its position on the chromosome. When replicated, the rate of transcription of the gene is doubled.
\item [{RNA-polymerases~and~ribosomes}] The production of mRNAs and proteins respectively requires RNA-polymerases and ribosomes. The concentrations of non-allocated (or free) RNA-polymerases and ribosomes respectively determine the rates of transcriptions and translations. During the time of elongation, the RNA-polymerase (resp. ribosome) remains sequestered on the DNA (resp. the mRNA). As the cell grows, new RNA-polymerases and ribosomes are created and participate in the production of proteins.
\end{description}




\subsection{\label{method:analysis} Analysis of each intermediate models}

Each of the intermediate models is systematically analyzed with the same method, see Figure~\ref{fig:procedure}. The average behavior of the model is theoretically predicted (either with explicit formulas for the intermediate models, or approximately for the complete model); it allows, for each of the 841 genes of Taniguchi~et~al.  (for which both the average protein and mRNA production have been measured), to determine the set of  parameters of the model. An overview of the parameters hence determined is shown in Figure~\ref{tab:Quantitative-summary-Model0-1}.  Simulations are then performed \textemdash{} using methods derived from Gillespie algorithm \textemdash{} with these parameters in order to check the accuracy of the average concentration and to predict the variance of the concentration of proteins. For each model, we then have the variance of protein concentration predicted for the wide range of genes measured in \cite{taniguchi_quantifying_2010}.

\subsection{\label{methods:mean variance} Means and Variances}
Throughout this paper we use the notation $\E{X}$ and $\V{X}$ for the mean and variance of a random variable $X$. We introduce the concentrations for mRNAs and proteins. If $M_{s}$ and $P_{s}$ denote the number of mRNAs and proteins of a given type at time $s$ and  $V(s)$ is the volume at this instant, the corresponding concentrations are $M_{s}/V(s)$ and $P_{s}/V(s)$. One of the goals of this work is of studying the properties of the  mean, $\E{P_{s}/V(s)}$, and variance, $\V{P_{s}/V(s)}$,  of these concentrations.  These quantities correspond to the mean and the variance of a population of synchronized cells of volume $V(s)$. The measurements of \cite{taniguchi_quantifying_2010} consider a cell population in exponential growth; by consequence, we also have to define the notions of global mean and variance for a heterogeneous population (see \cite{collins_rate_1962,sharpe_bacillus_1998,robert_division_2014} for the population distribution in exponential growth). By denoting by $\nu$ the age distribution of the cell population, we can define the global average $\Ehat{P/V}$ and variance $\Vhat{P/V}$ averaged over the population,
\begin{align}
\Ehat{P/V} & =\int_{0}^{\tau_{D}}\E{P_{s}/V(s)}\nu(s)\,\mathop{}\mathopen{}\mathrm{d} s\label{eq:Ehat_P}\\
\Vhat{P/V} & =\int_{0}^{\tau_{D}}\left[\E{\left(P_{s}/V(s)\right)^{2}}-\Ehat{P/V}^{2}\right]\nu(s)\,\mathop{}\mathopen{}\mathrm{d} s.\label{eq:Vhat_P}
\end{align}
We observe that the choice of $\nu$  the age distribution
(either uniform or corresponding to an exponentially growing population)
does seem to have much impact, see Section~\ref{subsec:other_pop_distrib} for more details.

\subsection{\label{methods:Intrinsic extrinsic} Intrinsic and Extrinsic Effects}
As for the previous studies of \cite{elowitz_stochastic_2002,swain_intrinsic_2002}, we want to decompose the protein variance that can be attributed to the intrinsic expression from the one due to extrinsic mechanisms.  It appears that this decomposition can be computed using two different ways.

First, there is the method used in Taniguchi~et~al. As the intrinsic noise is usually attributed to protein production mechanism alone, a model that represents only transcription and translation,
as the classical two-stage model, are usually considered as predicting the intrinsic noise.

Yet, in our case, the consideration of a more realistic mechanism for protein disappearance (through segregation at division rather than regular decay) prevent us from using directly the two-stage models as a quantitative representation of the intrinsic variance. Our first intermediate model considers genes in constant concentration and no interaction between the protein productions; in addition to that, we will propose a particular version of it (considered as our baseline model), which also do not consider any effect of compound partition at division (see Section~\ref{subsec:model1}). This baseline model will be very close to the two-stage model in that sense that it contains no other features than those intrinsically linked to the protein production. Therefore, we will use the protein variance predicted by this baseline model as the intrinsic protein variance.

From this baseline model, any additional variance added by the introduction to the model of external mechanisms (random partition at division, gene replication, etc.) would be considered as extrinsic. For every type of protein, we will look how the global variance $\Vhat{P/V}$ changes with the subsequent introduction of the external mechanisms.

Secondly, it appears that the previous method of extrinsic noise estimation is not exactly the same from the first attempt made by \cite{elowitz_stochastic_2002} using the dual reporter technique. \cite{hilfinger_separating_2011} showed that decomposition using the dual reporter technique can be interpreted as an estimator of the environmental state decomposition (also called the law of total variance). It decomposes the protein variance between the effects specifically due to the stochastic nature of the instants of birth and death of mRNAs and proteins (intrinsic noise) and the external effect of the biological environment (extrinsic noise).  If $Z$ represents the state of the cell, the number of RNA-polymerases, the volume, etc\ldots, then the protein concentration $P/V$ can be decomposed such as,
\begin{align}
\Vhat{P/V}  =\underbrace{\vphantom{\left[\Ehat{P|}\right]}\Ehat{\Vhat{P/V|Z}}}_{\begin{array}{c}
\Vhatt{P/V}{int}\\
\text{unexplained by }Z
\end{array}}+\underbrace{\Vhat{\Ehat{P/V|Z}}}_{\begin{array}{c}
\Vhatt{P/V}{ext}\\
\text{explained by }Z
\end{array}}.\label{eq:ESD}
\end{align}

The variable $Z$ represents the common environment in which the two similar genes of the dual reporter technique evolve; yet mathematically, it is dependent on the model that we consider (for each model, it is what is considered as being part of the "environment" of the gene). We have therefore explicitly described for each intermediate models what it represents in this context \rmodif{(the section~\ref{subsec:model2} shows an illustrative example of this decomposition)}. Once the variable $Z$ fixed, the decomposition is explicit and separate the total variance in two parts: $\Vhatt{P/V}{int}$ that corresponds to the intrinsic contributions and  $\Vhatt{P/V}{ext}$ corresponds to the variance induced by external contribution represented by $Z$ (volume fluctuation, concentrations of free RNA-polymerases and ribosomes, etc.) The term $\Vhatt{P/V}{int}$ is indeed the variance that can be expected from a model without any external fluctuation. A model like the two-stage model does not consider any change in the environment of the protein production, so the second term of the decomposition would be null.

In the intermerdiate models, we quantify the external contributions by these two means, either by looking at the increase of the global variance $\Vhat{P/V}$, or by performing the environmental state decomposition and looking the portion of external variance $\Vhatt{P/V}{int}/\Vhat{P/V}$ predicted. We will see that these two ways to quantify the external contributions of the protein variance are not strictly equivalent.

\section{Acknowledgments}

\newpage
\section*{Appendix}

\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thesection}{S\arabic{section}}
\renewcommand\thethm{S\arabic{thm}}

\setcounter{figure}{0}
\setcounter{equation}{0}
\setcounter{section}{0}
\setcounter{thm}{0}

Each of the three different models are analyzed following the scheme presented in Figure~\ref{fig:procedure} of the main article. We begin with a theoretical analysis to derive the average production of mRNAs and proteins along the cell cycle. Using these results, we then fit the parameters of the model to make them correspond to the experimental data of \cite{taniguchi_quantifying_2010}.  Then simulations are performed to predict the protein variance for each gene considered.


\section{Impact of Random Partitioning}
We present now  results for the model with volume growth, constant gene concentration and partitioning at division which has been presented in Figure~\ref{fig:mod1} of the main article. The results of the theoretical analysis part are similar for both cases of exact or random partitioning at division.


\subsection{Theoretical Analysis}

\subsubsection{Messenger-RNA Dynamic}\label{subsec:S_mod1_mrna_dynamic}

Here is shown the Proposition~\ref{prop:EM_s} that describes the average number of mRNAs at any instants of the cell cycle.

For any time $t\in\mathbb{R}_{+}$, denote by $M_{t}$ the number of mRNAs at this instant. We suppose that the initial time $t=0$ is a time of division; in this case, at each time $i\cdot\tau_{D}$ with $i\in\mathbb{N}$ are moments of division. For any $i\in\mathbb{N}$ , $M_{i\tau_{D}}$ denotes the number of mRNAs at the beginning of $i$-th cell cycle and $M_{i\tau_{D}-}$ the number of mRNAs in the $(i-1)$-th cell cycle just \emph{before} division.

We suppose that a lot of cell divisions have already occurred even before time $t=0$, and hence the considered cell cycle takes place when the embedded Markov chain $\left(M_{i\tau_{D}}\right)_{i}$ has already reached its equilibrium: the distribution $M_{i\tau_{D}}$ is the same as the distribution of $M_{(i+1)\tau_{D}}$.  If the equilibrium is already reached at time $0$, it implies that the distribution of any $M_{i\tau_{D}+s}$ for any $i\in\mathbb{N}$ and $s\in[0,\tau_{D}[$ is equal to the distribution of $M_{s}$. As a consequence, we can only consider the first cell cycle $s\in[0,\tau_{D}[$ to fully characterize the behavior of $M_{s}$ at any time $s\in\mathbb{R}_{+}$.

We propose here to describe the evolution of $\left(M_{s}\right)$ between times $0$ and $\tau_{D}$ (during this period of time, the number of mRNA approximately doubles). We first divide mRNAs into two categories,
\begin{itemize}
\item First group: mRNAs which were present at the birth of the cell. Each mRNA $i$ of the first group is characterized by $E_{\sigma_{1}}^{i}$, its lifetime given by an exponential random variable of rate $\sigma_{1}$.  The $i$-th mRNA still exists at time $s$ if and only if $E_{\sigma_{1}}^{i}>s$.  As a consequence, the number of mRNAs of this group still existing at time $s$ is given by
\begin{equation}
\sum_{i=1}^{M_{0}}\ind{E_{\sigma_{1}}^{i}>s}.\label{eq:mrna-first-group}
\end{equation}
\item Second group: mRNAs which have been created since the birth of the cell. The description of the number of mRNA of this group is more complicated. It is necessary to resort to the framework of Marked Poisson Point Processes (MPPP). An MPPP is a two-dimensional process.  It is based on a Poisson process where each of its random point is ``marked'' with another random variable; each point of a MPPP is a couple $(x,y)$ where $x$ is part of a Poisson point process and $y$ is the mark distributed according to a certain distribution.  One can refer to the first Chapter of \cite{robert_stochastic_2010} or \cite{kingman_poisson_1993} for the main results concerning MPPP.

We use this tool to characterize the number of mRNAs of the second group. In our case, the first variable $x$ represents the time at which the mRNA is created and the second variable $y$ represents the mRNA lifetime. Define $\mathcal{N}$ an MPPP of intensity
\[
\nu(\mathop{}\mathopen{}\mathrm{d} x,\mathop{}\mathopen{}\mathrm{d} y)=\lambda_{1}V(x)\mathop{}\mathopen{}\mathrm{d} x\otimes\sigma_{1}e^{-\sigma_{1}y}\mathop{}\mathopen{}\mathrm{d} y.
\]

It is noticeable that the underlying Poisson Process of this MPPP is not homogeneous. If the $i$-th mRNA of this group is born at time $x_{i}$ and its lifetime is $y_{i}$, then it exists at time $s$ if and only if $(x_{i},y_{i})\in\Delta_{s}$ with
\[
\Delta_{s}=\left\{ (x,y)\in\mathbb{R}_{+}^{2},\,0<x<s,\,y>s-x\right\} .
\]
One can refer to Figure~\ref{fig:MPPP-1}. Therefore the number of mRNAs of this group still present at time $s$ is given by
\begin{equation}
\mathcal{N}\left(\Delta_{s}\right)=\iint_{\mathbb{R}_{+}^{2}}\ind{(x,y)\in\Delta_{s}}\mathcal{N}\left(\mathop{}\mathopen{}\mathrm{d} x,\mathop{}\mathopen{}\mathrm{d} y\right).\label{eq:mrna-second-group}
\end{equation}

\end{itemize}
\begin{figure}
\begin{centering}
\begin{tikzpicture}[
        scale=1.6,
        IS/.style={blue, thick},
        LM/.style={red, thick},
        axis/.style={very thick, ->, >=stealth', line join=miter},
        important line/.style={thick}, dashed line/.style={dashed, thin},
        every node/.style={color=black},
        dot/.style={circle,fill=black,minimum size=4pt,inner sep=0pt,
            outer sep=-1pt},
    ]
	\fill[black!10]  (1.8,0) --  (1.8,2) -- (0,2) --  (0,1.8)  -- cycle;
	\draw[line width=0.5mm,green] (2.8,2)--(2.8,0)
			node[pos=0.21,below,green,rotate=90]{Division};
	\draw[line width=0.5mm,green] (0.025,2)--(0.025,0)
			node[pos=0.15,below,green,rotate=90]{Birth};

	\node[black] at (1.2,1.3){$\Delta_s$};
    \draw[axis,<->] (3,0) node(xline)[right,align=left] {Time} -|
                    (0,2) node(yline)[above] {Lifetime};
	\coordinate (starpt) at (1.25,1);
	\foreach \x/\y/\z in {0.5/0.25/0, 1.4/0.6/0, 2/1.3/0.5, 0.2/0.4/0,
						2.3/0.7/0.2, 2.55/0.2/0}
		{
    	\node[dot,color=blue] at (\x,\y) {};
		\draw[dashed line,blue] (\x,\y)--($(\x+\y-\z,\z)$);
		}

	\node[dot,color=blue] at (starpt) {};
	\draw[dashed line,red] let \p1=(starpt) in
		(starpt)--($(\x1,0)$) node[below,red]{$t_i$}
		(starpt)--($(0,\y1)$)
		;
    \draw[dashed line,red] let \p1=(starpt) in
		(starpt)--($(\x1+\y1,0)$) node[below,red]{$t_i+\delta_i$};

	\draw[<->,red,line width=0.25mm] let \p1=(starpt) in
		(-0.1,0)--(-0.1,\y1) node[above,midway,rotate=90,red]{
						$\delta_i\sim{\cal E}(\sigma_{1})$};

	\draw[<->,black,line width=0.25mm] (0,-0.3)--(2.8,-0.3)
		node[below,midway]{Point Poisson Process ($\lambda_1$)};
	\draw[line width=0.5mm] (1.8,2)--(1.8,0) node[below]{$s$};

\end{tikzpicture}
\par\end{centering}
\caption{\label{fig:MPPP-1}Illustration of the Marked Point Poisson Process describing the dynamic of mRNAs: each mRNA is characterized by the point $(x_{i},y_{i})$ (with $(x_{i},y_{i})$ following the MPPP $\mathcal{N}$, whose distribution is of intensity $\nu$). The random variable $x_{i}$ represents the time at mRNA creation and $y_{i}$ its lifetime, hence this mRNA exists from volume $x_{i}$ up to volume $x_{i}+y_{i}$; that is to say that mRNA is still present at time $s$, if and only if the point $(x_{i},y_{i})$ is in the set with $\Delta_{s}=\left\{ (x,y)\in\protect\mathbb{R}_{+}^{2},\,0<x<s,\,y>s-x\right\} $. }
\end{figure}
By summing the number of mRNAs for each group (Equations~(\ref{eq:mrna-first-group})
and~(\ref{eq:mrna-second-group})), it follows the total number of
mRNAs present at time $s\in[0,\tau_{D}[$:
\begin{equation}
M_{s}=\sum_{i=1}^{M_{0}}\ind{E_{\sigma_{1}}^{i}>s}+\mathcal{N}\left(\Delta_{s}\right).\label{eq:M_s-1}
\end{equation}
This description of the dynamic of $M_{s}$, together with the equilibrium
hypothesis which implies that $M_{0}\overset{\mathcal{D}}{=}M_{\tau_{D}}$,
allows to prove the next proposition.
\begin{prop}
\label{prop:EM_s}At equilibrium, the concentration of mRNAs at time
$s\in[0,\tau_{D}[$ of the cell cycle is
\[
\E{M_{s}/V(s)}=\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}.
\]
\end{prop}
\begin{proof}
 By taking the mean of Equation~(\ref{eq:M_s-1}), it follows for
any time $s$ of the cell cycle:
\[
\E{M_{s}}=\E{\sum_{i=1}^{M_{0}}\ind{E_{\sigma_{1}}^{i}>s}}+\E{\mathcal{N}\left(\Delta_{s}\right)}.
\]

Since all $\left(E_{\sigma_{1}}^{i}\right)_{i}$ are i.i.d. and independent
of $M_{0}$, the first term is given by
\[
\E{\sum_{i=1}^{M_{0}}\ind{E_{\sigma_{1}}^{i}>s}}  =\E{M_{0}}e^{-s\sigma_{1}}.
\]
For the second term, one has to remark that as $\mathcal{N}$ is a
MPPP, $\mathcal{N}\left(\Delta_{\tau_{D}-}\right)$ is a Poisson random
variable (Proposition~1.13.a of \cite{robert_stochastic_2010}).
The parameter of this Poisson random variable is given by
\[
\nu\left(\Delta_{s}\right)  =  \iint_{\Delta_{s}}\nu\left(\mathop{}\mathopen{}\mathrm{d} x,\mathop{}\mathopen{}\mathrm{d} y\right)
  =  V_{0}\frac{\lambda_{1}\sigma_{1}}{\log2+\sigma_{1}\tau_{D}}\left(2^{s/\tau_{D}}-e^{-\sigma_{1}s}\right).
\]
As a consequence, one gets that, for any time $s$ in the cell cycle,
\[
\E{M_{s}}=\E{M_{0}}e^{-s\sigma_{1}}+V_{0}\frac{\lambda_{1}\tau_{D}}{\tau_{D}\sigma_{1}+\log2}\cdot\left(2^{s/\tau_{D}}-e^{-s\sigma_{1}}\right).
\]

We still have to specify the mean number of mRNAs at birth $\E{M_{0}}$.
At the end of the cell cycle, for $s=\tau_{D}-$, the average number
of mRNAs is given by
\[
\E{M_{\tau_{D}-}}=\E{M_{0}}e^{-\tau_{D}\sigma_{1}}+V_{0}\frac{\lambda_{1}\tau_{D}}{\tau_{D}\sigma_{1}+\log2}\cdot\left(2-e^{-\tau_{D}\sigma_{1}}\right),
\]
and since at equilibrium,
\[
\E{M_{\tau_{D}}}=\E{M_{0}}=\E{M_{\tau_{D}-}}/2.
\]
Hence
\[
\E{M_{0}}\left(2-e^{-\tau_{D}\sigma_{1}}\right)=V_{0}\frac{\lambda_{1}\tau_{D}}{\tau_{D}\sigma_{1}+\log2}\cdot\left(2-e^{-\tau_{D}\sigma_{1}}\right),
\]
which gives the result.
\end{proof}

In particular, as the mean $\E{M_{s}/V(s)}$ does not change across the cell cycle, the global mRNA average $\Ehat{M/V}$ (as it is defined in Equation~(\ref{eq:Ehat_P})), does not depend on the choice of the population distributions $\nu$ in this case. It is given by,
\begin{equation}
\Ehat{M/V}=\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}.\label{eq:EhatM/V}
\end{equation}


\subsubsection{Protein Dynamic}\label{subsec:S_mod1_prot_dynamic}
The mean number of mRNAs is now determined for any moment of the cell
cycle. Each of the mRNAs potentially produces proteins at rate $\lambda_{2}$.
As for the mRNAs, we describe the number of proteins at time $s$
by grouping them into two categories.
\begin{itemize}
\item The $P_{0}$ proteins that were present at birth and which remain
in the bacteria during all the cell cycle (as said in the main article
the proteolysis is not considered in this model).
\item The proteins that have been created during the current cell cycle.
The rate of production is depending on the current number of mRNAs.
We consider $\mathcal{N}_{\lambda_{2}}^{i}$ (for $i\in\mathbb{N}$ and $i\geq1$)
independent Poisson Point Processes of intensity $\lambda_{2}$. If
the $i$-th mRNA exists at time $s$ (that is to say if $i\leq M_{s}$),
then the number of proteins produced by this mRNA between $s$ and
$s+\mathop{}\mathopen{}\mathrm{d} s$ is $\mathcal{N}_{\lambda_{2}}^{i}(\mathop{}\mathopen{}\mathrm{d} s)$.
\end{itemize}
To sum up, the number of proteins at a time $s$ of the cell cycle is given by
\begin{equation}
P_{s}=P_{0}+\sum_{i=1}^{\infty}\int_{0}^{s}\ind{i\leq M_{u}}\mathcal{N}_{\lambda_{2}}^{i}\left(\mathop{}\mathopen{}\mathrm{d} u\right)\mbox{.}\label{eq:P_s-1}
\end{equation}
The first term is the number of proteins at birth, and the second
take into account all the proteins created between times $0$ and
$s$. One can then determine the mean number of proteins at any time
$s$ of the cell cycle.
\begin{prop}
\label{prop:EP_s}At equilibrium, the concentration of proteins at
any time $s\in[0,\tau_{D}[$ of the cell cycle is
\[
\E{P_{s}/V(s)}=\frac{\lambda_{2}\tau_{D}}{\log2}\cdot\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}.
\]
\end{prop}

\begin{proof}
By taking the average of Equation~\eqref{eq:P_s-1}, one gets
\begin{multline*}
\E{P_{s}}  =\E{P_{0}}+\sum_{i=1}^{\infty}\E{\int_{0}^{s}\ind{i\leq M_{u}}\mathcal{N}_{\lambda_{2}}^{i}\left(\mathop{}\mathopen{}\mathrm{d} u\right)}=\E{P_{0}}+\sum_{i=1}^{\infty}\E{\int_{0}^{s}\ind{i\leq M_{u}}\lambda_{2}\,\mathop{}\mathopen{}\mathrm{d} u}\\
  =\E{P_{0}}+\lambda_{2}\int_{0}^{s}\E{\sum_{i=1}^{\infty}\ind{i\leq M_{u}}}\,\mathop{}\mathopen{}\mathrm{d} u=\E{P_{0}}+\lambda_{2}\int_{0}^{s}\E{M_{u}}\,\mathop{}\mathopen{}\mathrm{d} u.
\end{multline*}
As we know the mean number of mRNAs $\E{M_{u}}$ at  time $u$
of the cell cycle with Proposition~\ref{prop:EM_s},
\[
\E{P_{s}}=\E{P_{0}}+\frac{\lambda_{2}\tau_{D}}{\log2}\cdot\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}\cdot\left(V(s)-V_{0}\right).
\]
Since the system is at equilibrium, we have for time $\tau_{D}-$,
$\E{P_{\tau_{D}-}}=2\E{P_{0}}$; so
\[
\E{P_{0}}  =\frac{\lambda_{2}\tau_{D}}{\log2}\cdot\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}\cdot\left(V(\tau_{D}-)-V_{0}\right)
  =\frac{\lambda_{2}\tau_{D}}{\log2}\cdot\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}\cdot V_{0}.
\]
Consequently, for any time $s$ of the cell cycle,
\[
\E{P_{s}}=\lambda_{2}\frac{\tau_{D}}{\log2}\cdot\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}\cdot V_{0}\left(1+2^{s/\tau_{D}}-1\right);
\]
hence the result.
\end{proof}

In particular, as the mean $\E{P_{s}/V(s)}$ does not change across
the cell cycle, the global protein average is given by,
\begin{equation}
\Ehat{P/V}=\frac{\lambda_{2}\tau_{D}}{\log2}\cdot\frac{\lambda_{1}\tau_{D}}{\sigma_{1}\tau_{D}+\log2}.\label{eq:EhatP/V}
\end{equation}
\subsection{Parameters Estimation}\label{subsec:S_mod1_param}
For each gene measured in \cite{taniguchi_quantifying_2010}, we want to identify the set of corresponding parameters $\lambda_{1}$, $\sigma_{1}$ and $\lambda_{2}$. We also need to determine the ``global'' quantities $\tau_{D}$ and $V_{0}$. We first determine the parameters common to all genes. The division time $\tau_{D}$ is set to  $150$min in the article and the volume at birth $V_{0}$ is taken equal to $1.3\,\text{\textmu m}^{3}$.\footnote{The value of $V_{0}$, even if it is not explicitly given in \cite{taniguchi_quantifying_2010} can be deduced from the typical width given in its supplementary materials.}

Then we have to determine for each gene the three gene-specific parameters $\lambda_{1}$, $\sigma_{1}$ and $\lambda_{2}$. We consider the genes of the article for which was measured the empirical mean of messengers $\mu_{m}$ and proteins $\mu_{p}$ concentrations, as well as the mRNA half-life time $\tau_{m}$.

First we determine the rate of mRNA degradation for each gene with the measured mRNA half-life time $\tau_{m}$.  A half-life $\tau_{m}$ indicates that a mRNA has a probability $1/2$ to disappear within a duration $\tau_{m}$, hence $e^{-\sigma_{1}\tau_{m}}=1/2$. From that, we can compute the rate $\sigma_{2}$ (specific for each type of mRNA),
\[
\sigma_{1}=\log2/\tau_{m}\mbox{.}
\]

Then we can identify the averages of mRNA and protein concentrations of the model (respectively $\Ehat{M/V}$ and $\Ehat{P/V}$) with the empirical averages $\mu_{m}$ of mRNA concentration and $\mu_{p}$ of protein concentration of the article. With Equations~(\ref{eq:EhatM/V}) and~(\ref{eq:EhatP/V}), the parameters $\lambda_{1}$ and $\lambda_{2}$ are
\[
\lambda_{1}  =\mu_{m}\cdot\frac{\sigma_{1}\tau_{D}+\log2}{\tau_{D}},\qquad
\lambda_{2}  =\mu_{p}\cdot\frac{\log2}{\tau_{D}}\cdot\frac{\sigma_{1}\tau_{D}+\log2}{\lambda_{1}\tau_{D}}.
\]

A summary of the different parameters can be seen in Figure~\ref{tab:Quantitative-summary-Model0-1}.  Having determined all the parameters allows to perform simulations of the model using stochastic algorithm in order to assess the variability of every protein and compare them with those experimentally obtained in \cite{taniguchi_quantifying_2010}.

\subsection{\label{subsec:S_mod1_Gillespie-Algorithms}Simulations}
When performing simulations, one needs to take care of the non-homogeneity of the Poisson processes describing mRNA creation times: the rate protein production $\lambda_{1}V(s)$ is not a homogeneous rate as it changes with time. That does not allow a direct application of Gillepsie method \citep{gillespie_exact_1977},  an extension for non-homogenous processes has to be used.

\cite{gillespie_exact_1977} describes an algorithm to simulate stochastic trajectories such as the quantities of different chemical species interacting together. The main idea is to consider the state of a system (for instance the number of each chemical compounds) and to compute the first reaction to occur, as well as the time when it happens. Once both pieces of information computed, one change the current state of the system accordingly with the reaction, and update the time.

One important hypothesis is that all reactions occur at exponential times (even if the rates of these exponential times may depend on the current state of the system).  In the current intermediate model, at any time $s$, the state is described by $\left(M_{s},P_{s}\right)$ (respectively, the number of mRNAs and proteins), and the rate of mRNA production is $\Lambda(s)=\lambda_{1}V(s)$ with $\lambda_{1}$ a parameter and $V(s)$ the non-constant volume of the cell. The parameter $\Lambda(s)$ does not depend on the state $\left(M_{s},P_{s}\right)$ but is time dependent through $V(s)$; for this reason, it is not an exponential time.

In this case, the duration of time $T$ until the next mRNA production is characterized by
\[
\P{T>x}=\exp\left(-\int_{0}^{x}\Lambda(s)\,\mathop{}\mathopen{}\mathrm{d} s\right),  \quad x>0.
\]
which is not an exponential distribution as $\Lambda$ is non-constant.  To compute $T$, we consider that $\Lambda(s)$ is strictly positive for any $s\in\mathbb{R}_{+}$, as a consequence $F(x):=\int_{0}^{x}\Lambda(s)\,\mathop{}\mathopen{}\mathrm{d} s$ is strictly increasing. Let $E$ be an exponential random variable with  parameter $1$. We have hence
\[
\forall y>0\quad\P{E>y}=\exp\left(-y\right).
\]
If we consider the case of $y=F(x)$, since $F$ is strictly increasing, hence
\[
\P{E>y}=\exp\left(-F(x)\right)\qquad\text{and}\qquad\P{E>y}=\P{F^{-1}(E)>x}.
\]
As a consequence the random variable $F^{-1}(E)$ has the same distribution as $T$.

Based on that we can propose a new version of the algorithm of Gillespie that can take into account non-exponential times such as $T$.
\begin{algo}
The equivalent of Gillespie algorithm that considers non-homogeneous events is
\begin{enumerate}
\item Initialization: Initialize time of molecules in the system and the
time.
\item Next exponential event: determine the next event that occurs at an exponential
time as in Gillespie algorithm.
\item Next non-homogeneous event: determine the next event that occurs at
non-homogeneous rates with the method previously described.
\item Update: choose between events of Step 2 or Step 3 that happen first.
Update the time and the molecule count accordingly.
\item Iterate: Consider again the Step 2 unless it has reached the end of
the simulation.
\end{enumerate}
\end{algo}

\begin{figure}
\subfloat[\label{fig:S_mod1_mRNA-prot-check}]{\begin{tikzpicture}[scale=0.75]
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[height=5.cm,page=1]{model1/Figure_20170130_comparision_mod0_exp}};

\node[anchor=south west,inner sep=0] (im1) at (3.,1.2) {

\includegraphics[width=4.cm,page=2]{model1/Figure_20170130_comparision_mod0_exp}};

\end{tikzpicture}

}~~\subfloat[\label{fig:S_mod1-FanoFactor}]{\includegraphics[height=5cm,page=2]{model1/Figure_20160713_general_results}
}

\subfloat[\label{fig:S_mod1-VarRatioVsMean}]{\includegraphics[height=5cm,page=3]{model1/Figure_20160713_general_results}
}

\caption{Intermediate model with volume growth and partition at division. (A) Comparison
of the average productions of proteins (main) and mRNAs (inset) obtained
in the simulations and those experimentally measured. (B) The Fano
factor as a function of the average protein production: the production
with the lowest Fano Factor tends to be the less expressed. (C) For each
type of protein, the variance in the case of exact partition divided by
the variance in the case of random partition as a function of the
protein average production.
}
\end{figure}


\subsection{\label{subsec:S_mod1_ESD}Environmental State Decomposition}

As explained in the main article, the dual reporter technique \citep{elowitz_stochastic_2002}
that compares the expression of two similar genes (with the same promoter
and RBS, at an equivalent position on the chromosome) in the same
cell can be interpreted as an estimator of the environmental state
decomposition \citep{hilfinger_separating_2011}. If $Z$ the cell
state (i.e., the common environment in which the two genes of the
dual reporter technique are expressed), then we can apply the law
of total variance on the protein number,
\[
\Vhat{P/V}=\underbrace{\vphantom{\left[\Ehat{P|}\right]}\Ehat{\Vhat{P/V|Z}}}_{\Vhatt{P/V}{int}}+\underbrace{\Vhat{\Ehat{P/V|Z}}}_{\Vhatt{P/V}{ext}}.
\]
In order to be applied here, we need to specify what does the gene
environment $Z$ refers to. In our model, two similar genes (with
the same parameters $\lambda_{1}$, $\sigma_{1}$,$\lambda_{2}$ )
in the same cell would undergo the same volume growth. But each gene
would undergo a specific partition at division (for instance, in the
case of random partitioning, the partitions of the proteins of one
of the gene, is uncorrelated with the partition of the proteins of
the other gene).

With only considering the volume as the gene environment, we end up
with the following decomposition,
\begin{align*}
\Vhatt{P/V}{int} & =  \frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\V{P_{s}/V(s)}\mathop{}\mathopen{}\mathrm{d} s,\\
\Vhatt{P/V}{ext} & =  \frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\E{P_{s}/V(s)}^{2}\mathop{}\mathopen{}\mathrm{d} s-\Ehat{P/V}^{2}.
\end{align*}


\subsection{\label{subsec:S_mod1_simple_model_partition}Simplified model for
the random partitioning}

Here we explain the simplified model used to make the prediction of
protein noise ratio in blue dash line of Figure~\ref{fig:mod1_var_ratio}.

For a  given  quantity $P$ associated to the cell, the partitioning can be performed in two ways, either exact or random. The result in each case  will  be denoted respectively by $P_{e}$ and $P_{r}$. During division, the volume is divided by two, changing from $2V_{0}$ to $V_{0}$. In order to be plotted in Figure~\ref{fig:mod1_var_ratio}, we need to consider the variance of protein concentration after division
\[
\eta:=\frac{\V{P_{e}/V_{0}}}{\V{P_{r}/V_{0}}}\quad\text{as a function of\quad}x:=\frac{\V{P/(2V_{0})}}{\E{P/(2V_{0})}}.
\]

\begin{prop}
The Variance ratio $\eta$ as a function of $x$ is given by
\[
\eta=\frac{2V_{0}x}{2V_{0}x+1}.
\]
\end{prop}


\begin{proof}
We have the quantity before division $P$. Since by definition, we have
that $P_{e}=P/2$ and
\[
\E{P_{e}/V_{0}}  =\E{P/(2V_{0})}, \quad  \V{P_{e}/V_{0}}  =\V{P/(2V_{0})}.
\]

For the effect of binomial division, see Lemma~\ref{lem:Binom_div}, it describes the effect of the binomial division on the means and on the variances of several quantities. By the volume in order to observe the concentrations, one gets that
\begin{align*}
\E{P_{r}/V_{0}} & =\frac{\E P}{2V_{0}}=\E{P/\left(2V_{0}\right)}
\end{align*}
and
\[
\V{P_{r}/V_{0}}  =\frac{\V{P_{r}}}{V_{0}^{2}}=\frac{\V P+2\E{P_{r}}}{4V_{0}^{2}}
  =\V{P/\left(2V_{0}\right)}+\frac{\E{P/(2V_{0})}}{2V_{0}}.
\]
As a consequence, this gives the relation
\[
\eta=\frac{\V{P/(2V_{0})}}{\V{P/\left(2V_{0}\right)}+\E{P/(2V_{0})}/(2V_{0})}=\frac{2V_{0}x}{2V_{0}x+1}
\]
\end{proof}

\section{Impact of Gene Replication}

We present here results relative to the model presented in Figure~\ref{fig:mod2} of the main article with volume growth, partition at division and gene replication.  For this intermediate model, we still follow the analysis scheme presented in Figure~\ref{fig:procedure}, but here, we are also able to produce analytic results for the variance of mRNAs and proteins.

\subsection{Theoretical analysis}

\subsubsection{\label{subsec:S_mod2_mRNA-number}Dynamics of mRNA number }

The aim of this section is to prove Theorem~\ref{thm:xs_rep} which states that at any time of the cell, the mRNA number follows a Poisson distribution. To do so, we first give a description of the number mRNAs at any time in the cell cycle using a Marked Poisson Point Process.  With this description, we will be able to show Proposition~\ref{prop:Distrib_M_0}, that the distribution of $M_{0}$ at the beginning of the cell cycle is a Poisson distribution. This proposition will allow to finally prove the main theorem of the subsection.

If time $s=0$ is the beginning of a new cell cycle and if the system is already at equilibrium in the same sense as the previous models (see Section~\ref{subsec:S_mod1_mrna_dynamic}).  We consider that $M_{0}$, the number of mRNAs at birth is known.  As in Section~\ref{subsec:S_mod1_mrna_dynamic}, we assort mRNAs in independent groups; here we consider three categories.
\begin{itemize}
\item mRNAs which were present at the birth of the cell. Each of them is characterized by its lifetime given by an exponential time of rate $\sigma_{1}$. The $i$-th mRNA is still present at time $s$ if and only if $E_{\sigma_{1}}^{i}>s$, with $\left(E_{\sigma_{1}}^{i}\right)$ being i.i.d. exponential random variables of parameter $\sigma_{1}$.
\item mRNAs created since the birth of the cell by the first copy of the gene. The $i$-th mRNA of this group is characterized by the time of creation $t_{i}$ given by a Poisson Process of rate $\lambda_{1}$ and its lifetime $\delta_{i}$ given by an exponential time of rate $\sigma_{1}$.
\item mRNAs created since the gene replication by the second copy of the gene. As in the previous group, the $i$-th mRNA is characterized by the time of creation $t_{i}$ given by a Poisson Process of rate $\lambda_{1}$ and its lifetime $\delta_{i}$ given by an exponential time of rate $\sigma_{1}$. But here, the Poisson Process of rate $\lambda_{1}$ begins at time $\tau_{R}$, the time of replication of the gene.
\end{itemize}
As in Section~\ref{subsec:S_mod1_mrna_dynamic}, one can represent the number of mRNAs of the second and the third group as two independent MPPPs $\mathcal{N}$ and $\mathcal{N}'$. The first variable $x$ of each of these MPPPs will represent the time. The intensity of each of the MPPP is the same,
\begin{align*}
\nu(\mathop{}\mathopen{}\mathrm{d} x,\mathop{}\mathopen{}\mathrm{d} y) & =\lambda_{1}\mathop{}\mathopen{}\mathrm{d} x\otimes\sigma_{1}e^{-\sigma_{1}y}\mathop{}\mathopen{}\mathrm{d} y.
\end{align*}
The only difference between $\mathcal{N}$ and $\mathcal{N}'$ is
the fact that they begin at time $0$ for $\mathcal{N}$ and at time
$\tau_{R}$ for $\mathcal{N}'$ (see Figure~\ref{fig:MPPP-2}). As
a consequence, if we consider an mRNA of either group, the conditions
of its existence at time $s\in[0,\tau_{D}[$ are respectively,
\begin{itemize}
\item if it is in the second group: $(t_{i},\delta_{i})\in\Delta_{s}$ with
$\Delta_{s}=\left\{ (x,y),\,0<x<s,\,y>s-x\right\} ,$
\item if it is in the third group: $(t_{i},\delta_{i})\in\Delta'_{s}$ with
$\Delta'_{s}=\left\{ (x,y),\,\tau_{R}<x<s,\,y>s-x\right\} .$
\end{itemize}
\begin{center}
\begin{figure}
\begin{centering}
\subfloat[\label{fig:MPPP-2_a}Representation of $\mathcal{N}$]{\begin{centering}
\begin{tikzpicture}[
        scale=1.6,
        IS/.style={blue, thick},
        LM/.style={red, thick},
        axis/.style={very thick, ->, >=stealth', line join=miter},
        important line/.style={thick}, dashed line/.style={dashed, thin},
        every node/.style={color=black},
        dot/.style={circle,fill=black,minimum size=4pt,inner sep=0pt,
            outer sep=-1pt},
    ]
	\fill[black!10]  (1.8,0) --  (1.8,2) -- (0,2) --  (0,1.8)  -- cycle;
	\draw[line width=0.5mm,green] (2.8,2)--(2.8,0)
			node[pos=0.21,below,green,rotate=90]{Division};
	\draw[line width=0.5mm,green] (0.025,2)--(0.025,0)
			node[pos=0.15,below,green,rotate=90]{Birth};
	\draw[line width=0.5mm,green] (1,2)--(1,0)
			node[pos=0.25,below,green,rotate=90]{Replication};

	\node[black] at (0.7,1.6){$\Delta_s$};
    \draw[axis,<->] (3,0) node(xline)[right,align=left] {Time} -|
                    (0,2) node(yline)[above] {Lifetime};
	\coordinate (starpt) at (1.25,1);
	\foreach \x/\y/\z in {0.5/0.25/0, 1.4/0.6/0, 2/1.3/0.5, 0.2/0.4/0,
						2.3/0.7/0.2, 2.55/0.2/0}
		{
    	\node[dot,color=blue] at (\x,\y) {};
		\draw[dashed line,blue] (\x,\y)--($(\x+\y-\z,\z)$);
		}

	\node[dot,color=blue] at (starpt) {};
	\draw[dashed line,red] let \p1=(starpt) in
		(starpt)--($(\x1,0)$) node[below,red]{$t_i$}
		(starpt)--($(0,\y1)$)
		;
    \draw[dashed line,red] let \p1=(starpt) in
		(starpt)--($(\x1+\y1,0)$) node[below,red]{$t_i+\delta_i$};

	\draw[<->,red,line width=0.25mm] let \p1=(starpt) in
		(-0.1,0)--(-0.1,\y1) node[midway,above,rotate=90,red]{
						$\delta_i\sim{\cal E}(\sigma_{1})$};

	\draw[<->,black,line width=0.25mm] (0,-0.3)--(2.8,-0.3)
		node[below,midway]{Point Poisson Process ($\lambda_1$)};
	\draw[line width=0.5mm] (1.8,2)--(1.8,0) node[below]{$s$};

\end{tikzpicture}
\par\end{centering}
}~~\subfloat[\label{fig:MPPP-2_b}Representation of $\mathcal{N}'$]{\begin{centering}
\begin{tikzpicture}[
        scale=1.6,
        IS/.style={blue, thick},
        LM/.style={red, thick},
        axis/.style={very thick, ->, >=stealth', line join=miter},
        important line/.style={thick}, dashed line/.style={dashed, thin},
        every node/.style={color=black},
        dot/.style={circle,fill=black,minimum size=4pt,inner sep=0pt,
            outer sep=-1pt},
    ]
	\fill[black!10]  (1.8,0) --  (1.8,2) -- (1,2) --  (1,0.8)  -- cycle;

	\draw[line width=0.5mm,green] (2.8,2)--(2.8,0)
			node[pos=0.21,below,green,rotate=90]{Division};
	\draw[line width=0.5mm,green] (0.025,2)--(0.025,0)
			node[pos=0.15,below,green,rotate=90]{Birth};
	\draw[line width=0.5mm,green] (1,2)--(1,0)
			node[pos=0.25,below,green,rotate=90]{Replication};

	\node[black] at (1.5,1.6){$\Delta_s'$};
    \draw[axis,<->] (3,0) node(xline)[right,align=left] {Time} -|
                    (0,2) node(yline)[above] {Lifetime};
	\coordinate (starpt) at (1.25,1);
	\foreach \x/\y/\z in { 1.1/0.6/0, 1.65/0.3/0,
						2/1.6/0.8, 2.3/0.4/0, 2.55/0.6/0.35}
		{
    	\node[dot,color=violet] at (\x,\y) {};
		\draw[dashed line,violet] (\x,\y)--($(\x+\y-\z,\z)$);
		}

	\node[dot,color=violet] at (starpt) {};
	\draw[dashed line,red] let \p1=(starpt) in
		(starpt)--($(\x1,0)$) node[below,red]{$t_i$}
		(starpt)--($(0,\y1)$)
		;
    \draw[dashed line,red] let \p1=(starpt) in
		(starpt)--($(\x1+\y1,0)$) node[below,red]{$t_i+\delta_i$};

	\draw[<->,red,line width=0.25mm] let \p1=(starpt) in
		(-0.1,0)--(-0.1,\y1) node[midway,above,rotate=90,red]{
						$\delta_i\sim{\cal E}(\sigma_{1})$};

	\draw[<->,black,line width=0.25mm] (1,-0.3)--(2.8,-0.3)
		node[below,midway]{Point Poisson Process ($\lambda_1$)};
	\draw[line width=0.5mm] (1.8,2)--(1.8,0) node[below]{$s$};

\end{tikzpicture}
\par\end{centering}
}
\par\end{centering}
\caption{\label{fig:MPPP-2}The illustration of the Marked Point Poisson Processes
that describe the dynamic of mRNAs: each mRNA is characterized by
the point $(t_{i},\delta_{i})$, with $(t_{i},\delta_{i})$ following
MPPP $\mathcal{N}$ or $\mathcal{N}'$. Both MPPPs are of intensity
$\nu$. The random variable $t_{i}$ represents its birth time and
$\delta_{i}$ its lifetime, hence this mRNA exists from time $t_{i}$
up to time $t_{i}+\delta_{i}$. The only difference for the two processes
is the starting time: the process $\mathcal{N}$ in \textbf{(A)} begins
at birth (in particular, an mRNA is still present at time $s$, if
and only if the point $(t_{i},\delta_{i})$ is in the set $\Delta_{s}=\left\{ (x,y),\,0<x<s,\,y>s-x\right\} $);
the process $\mathcal{N}'$ , in \textbf{(B)}, begins at replication
(in particular an mRNA is still present at time $s$, if and only
if the point $(t_{i},\delta_{i})$ is in the set $\Delta'_{s}=\left\{ (x,y),\,\tau_{R}<x<s,\,y>s-x\right\} $). }
\end{figure}
\par\end{center}

Hence, we can describe the number of mRNAs at any time $s\in[0,\tau_{D}[$ as follows,
\begin{equation}
M_{s}=\sum_{i=1}^{M_{0}}\ind{E_{\sigma_{1}}^{i}>s}+\mathcal{N}\left(\Delta_{s}\right)+\ind{s\geq\tau_{R}}\mathcal{N}'_{\lambda_{1}}\left(\Delta'_{s}\right).\label{eq:M_s-2}
\end{equation}
Each term corresponds to each group of mRNAs previously described.

At first we want to characterize the distribution of $M_{0}$, the
number of mRNAs at the birth of the cell. To do so, we use the equilibrium
hypothesis that implies that $M_{0}\overset{\mathcal{D}}{=}M_{\tau_{D}}$.
\begin{prop}
\label{prop:Distrib_M_0}At equilibrium, the number of mRNAs at birth
$M_{0}$ follows a Poisson distribution of parameter:
\[
x_{0}=\frac{\lambda_{1}}{\sigma_{1}}\left[1-\frac{e^{-(\tau_{D}-\tau_{R})\sigma_{1}}}{2-e^{-\tau_{D}\sigma_{1}}}\right]\mbox{.}
\]
\end{prop}
\begin{proof}
When $s{=}\tau_{D}-$,  by Relation~\eqref{eq:M_s-2},
\[
M_{\tau_{D}-}=\sum_{i=1}^{M_{0}}\ind{E_{\sigma_{1}}^{i}>\tau_{D}-}+\mathcal{N}_{\lambda_{1}}\left(\Delta_{\tau_{D}-}\right)+\mathcal{N}'_{\lambda_{1}}\left(\Delta'_{\tau_{D}-}\right)\mbox{.}
\]
The first term corresponds to initial messengers not degraded after the time $\tau_{D}$. Suppose that $M_{0}$ is distributed according to a Poisson distribution with parameter $x_{0}$, then the random variable
\[
\sum_{i=1}^{M_{0}}\ind{E_{\sigma_{1}}^{i}>\tau_{D}-}
\]
follows also a Poisson distribution with parameter $x_{0}e^{-\tau_{D}\sigma_{1}}$. Operation of thinning of Poisson processes, see \cite{kingman_poisson_1993} for example.

The second term corresponds to mRNAs that were created by the first
copy of the gene and which are still present at division. Since $\mathcal{N}$
is a MPPP, $\mathcal{N}\left(\Delta_{\tau_{D}-}\right)$ is a Poisson
random variable (Proposition~1.13 of \cite{robert_stochastic_2010})
with parameter
\[
\nu\left(\Delta_{\tau_{D}-}\right)=\int_{0}^{\tau_{D}}\int_{\tau_{D}-x}^{\infty}\lambda_{1}\sigma_{1}e^{-\sigma_{1}y}\,\mathop{}\mathopen{}\mathrm{d} y\,\mathop{}\mathopen{}\mathrm{d} x=\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-\tau_{D}\sigma_{1}}\right)\mbox{.}
\]

The third term corresponds to mRNAs that were created by the second copy of the gene (replicated at time $\tau_{R}$) and which are still present at division. As before, $\mathcal{N}'\left(\Delta'_{\tau_{D}-}\right)$ is a Poisson random variable with parameter
\[
\nu\left(\Delta'_{\tau_{D}-}\right)=\int_{\tau_{R}}^{\tau_{D}}\int_{\tau_{D}-x}^{\infty}\lambda_{1}\sigma_{1}e^{-\sigma_{1}y}\,\mathop{}\mathopen{}\mathrm{d} y\,\mathop{}\mathopen{}\mathrm{d} x=\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-\left(\tau_{D}-\tau_{R}\right)\sigma_{1}}\right).
\]

As $M_{\tau_{D}-}$ is the sum of three independent Poisson random
variables, one gets that
\begin{align*}
M_{\tau_{D}-} & \sim\mathcal{P}\left(x_{0}e^{-\sigma_{1}\tau_{D}}+\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-\tau_{D}\sigma_{1}}\right)+\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-\left(\tau_{D}-\tau_{R}\right)\sigma_{1}}\right)\right)\\
 & \sim\mathcal{P}\left(x_{0}e^{-\sigma_{1}\tau_{D}}+\frac{\lambda_{1}}{\sigma_{1}}\left(2-e^{-\tau_{D}\sigma_{1}}-e^{-\left(\tau_{D}-\tau_{R}\right)\sigma_{1}}\right)\right).
\end{align*}

Between $\tau_{D}-$ and $\tau_{D}$, with the random sampling, each
mRNA has an equal chance to stay or to disappears, therefore
\[
M_{\tau_{D}}=\sum_{1=0}^{M_{\tau_{D}-}}B_{1/2,i}
\]
with $\left(B_{1/2,i}\right)$ i.i.d. Bernoulli random variables with parameter $1/2$. The random variable $B_{1/2,i}$ determines if the $i$-th mRNA is in the next considered cell or not. The random variable $M_{\tau_{D}}$ hence follows a Poisson distribution such that
\[
M_{\tau_{D}}\sim\left.\mathcal{P}\left(\left[x_{0}e^{-\sigma_{1}\tau_{D}}+\frac{\lambda_{1}}{\sigma_{1}}\left(2-e^{-\tau_{D}\sigma_{1}}-e^{-\left(\tau_{D}-\tau_{R}\right)\sigma_{1}}\right)\right]\right/2\right).
\]
Since the system is at equilibrium, one has $M_{0}\overset{\mathcal{D}}{=}M_{\tau_{D}}$,
therefore
\[
x_{0}=\frac{1}{2}\left(x_{0}e^{-\sigma_{1}\tau_{D}}+\frac{\lambda_{1}}{\sigma_{1}}\left(2-e^{-\tau_{D}\sigma_{1}}-e^{-\left(\tau_{D}-\tau_{R}\right)\sigma_{1}}\right)\right),
\]
which gives
\[
x_{0}=\frac{\lambda_{1}}{\sigma_{1}}\left[1-\frac{e^{-(\tau_{D}-\tau_{R})\sigma_{1}}}{2-e^{-\tau_{D}\sigma_{1}}}\right]\mbox{.}
\]
Since the equilibrium distribution is unique, the number of mRNAs
at birth follows a Poisson distribution of parameter $x_{0}$ at equilibrium.
\end{proof}

We have determined the equilibrium distribution of the embedded Markov Chain $\left(M_{i\tau_{D}}\right)_{i\in\mathbb{N}}$. Now, we analyze the distribution of mRNA number at any instant $s$ of the cell cycle.
\begin{thm}
\label{thm:xs_rep}At equilibrium, at  time $s$ in the cell cycle,
the distribution of the mRNA number $M_{s}$ is Poisson with  parameter
\[
x_{s}=\frac{\lambda_{1}}{\sigma_{1}}\left[1-\frac{e^{-(s+\tau_{D}-\tau_{R})\sigma_{1}}}{2-e^{-\tau_{D}\sigma_{1}}}+\ind{s\geq\tau_{R}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right]\mbox{.}
\]
In particular, the mean and the variance of mRNA concentration are known at any time $s$ of the cell cycle,
\begin{align*}
\E{M_{s}/V(s)} & =  \frac{\lambda_{1}}{\sigma_{1}V(s)}\left[1-\frac{e^{-(s+\tau_{D}-\tau_{R})\sigma_{1}}}{2-e^{-\tau_{D}\sigma_{1}}}+\ind{s\geq\tau_{R}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right],\\
\V{M_{s}/V(s)} & =  \frac{\lambda_{1}}{\sigma_{1}V(s)^{2}}\left[1-\frac{e^{-(s+\tau_{D}-\tau_{R})\sigma_{1}}}{2-e^{-\tau_{D}\sigma_{1}}}+\ind{s\geq\tau_{R}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right].
\end{align*}

\end{thm}
\begin{proof}
At a moment $s$ of the cell cycle,  the moment-generating
function of $M_{s}$ at $\xi<0$ is given by
\[
\E{\exp\left(\xi M_{s}\right)}=\E{\exp\left(\xi\left(\sum_{i=1}^{M_{0}}\ind{E_{\sigma_{1}}^{i}>s}+\mathcal{N}\left(\Delta_{s}\right)+\ind{s\geq\tau_{R}}\mathcal{N}'_{\lambda_{1}}\left(\Delta'_{s}\right)\right)\right)}.
\]
Since $M_{0}$, $E_{\sigma_{1}}^{i}$, $\mathcal{N}_{\lambda_{1}}$
and $\mathcal{N}'_{\lambda_{1}}$ are all independent, it follows
that
\[
\E{\exp\left(\xi M_{s}\right)}=\E{\exp\left(\sum_{i=0}^{M_{0}}\xi\ind{E_{\sigma_{1}}^{i}>s}\right)}\cdot\E{\exp\left(\xi\mathcal{N}\left(\Delta_{s}\right)\right)}\cdot\E{\exp\left(\xi\ind{s\geq\tau_{R}}\mathcal{N}'\left(\Delta'_{s}\right)\right)}.
\]
For the first factor, since all the random variables $\ind{E_{\sigma_{1}}^{i}>s}$ are i.i.d. Bernoulli variables with parameter $e^{-s\sigma_{1}}$ and independent of $M_{0}$, one has
\[
\E{\exp\left(\sum_{i=0}^{M_{0}}\xi\ind{E_{\sigma_{1}}^{i}>s}\right)}=\E{\E{\exp\left(\xi\ind{E_{\sigma_{1}}^{1}>s}\right)|M_{0}}^{M_{0}}}=\E{\exp\left(1+e^{-s\sigma_{1}}(e^{\xi}-1)\right)^{M_{0}}}.
\]
With Proposition~\ref{prop:Distrib_M_0}, $M_{0}$ is known to be a Poisson random variable of parameter $x_{0}$, hence, with the probability generating function of a Poisson random variable,
\[
\E{\exp\left(\sum_{i=0}^{M_{0}}\xi\ind{E_{\sigma_{1}}^{i}>s}\right)}=\E{\exp\left(x_{0}e^{-s\sigma_{1}}\left(e^{\xi}-1\right)\right)}
\]
holds. For the second factor, one can recall that $\mathcal{N}\left(\Delta_{s}\right)$ is a Poisson random variable. As in Proposition~\ref{prop:Distrib_M_0}, its parameter can be calculated
\[
\nu\left(\Delta_{s}\right)=\int_{0}^{s}\int_{\tau_{D}-x}^{\infty}\lambda_{1}\sigma_{1}e^{-\sigma_{1}y}\,\mathop{}\mathopen{}\mathrm{d} y\,\mathop{}\mathopen{}\mathrm{d} x=\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-s\sigma_{1}}\right).
\]
Identically for the third factor, $\mathcal{N}'\left(\Delta'_{s}\right)$ is a Poisson random variable of parameter
\[
\nu\left(\Delta'_{s}\right)=\int_{\tau_{R}}^{s}\int_{\tau_{D}-x}^{\infty}\lambda_{1}\sigma_{1}e^{-\sigma_{1}y}\,\mathop{}\mathopen{}\mathrm{d} y\,\mathop{}\mathopen{}\mathrm{d} x=\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-\left(s-\tau_{R}\right)\sigma_{1}}\right).
\]

As a consequence, the moment generating function of $M_{s}$ is
\begin{align*}
\E{\exp\left(\xi M_{s}\right)} & =  \E{\left(x_{0}e^{-s\sigma_{1}}+\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-s\sigma_{1}}\right)+\ind{s>\tau_{R}}\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right)\left(e^{\xi}-1\right)}\\
 & =  \E{\left(x_{0}e^{-s\sigma_{1}}+\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-s\sigma_{1}}+\ind{s>\tau_{R}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right)\right)\left(e^{\xi}-1\right)}
\end{align*}
which is the moment-generating function of a Poisson random variable
of parameter
\[
x_{0}e^{-s\sigma_{1}}+\frac{\lambda_{1}}{\sigma_{1}}\left(1-e^{-s\sigma_{1}}+\ind{s>\tau_{R}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right).
\]
\end{proof}

\subsubsection{\label{subsec:S_mod2_Protein-number}Dynamics of Proteins }

\global\long\def\Es#1#2{\left\langle #1\right\rangle _{#2}}

As for the previous analysis of the mRNA number, we search an expression
for the protein production through the cell cycle. This case is more
complicated than the mRNA case and we will only calculate analytical
expressions only for the first two moments of $P_{s}$.

Propositions~\ref{prop:EPs_and_VPs} and~\ref{prop:EPs_and_VPs-2}
are the main theoretical results of this section: for any time $s$
of the cell cycle, it gives explicit expressions for the mean $\E{P_{s}}$
and the variance $\V{P_{s}}$ of the protein number. This result is
important as it will be used to directly calculate  the mean $\Ehat{P/V}$
and variance $\Vhat{P/V}$ of the protein concentration averaged across
the cell cycle without using simulations: only with the parameters
of the model ($\lambda_{1}$, $\sigma_{1}$, $\lambda_{2}$, $\tau_{R}$
and $\tau_{D}$), we will be able to know the behavior of the protein
concentration in terms of variance.

In order to prove the Propositions~\ref{prop:EPs_and_VPs} and~\ref{prop:EPs_and_VPs-2},
we will characterize $\E{P_{s}}$ and $\V{P_{s}}$ in the two following
cases:
\begin{enumerate}
\item First, we consider the case before replication ($s<\tau_{R}$ ). We begin by considering that the state of the cell at birth $(M_{0},P_{0})$ is known and we calculate the first two moments of $P_{s}$ for any time $s<\tau_{R}$ (Corollary~\ref{cor:EPs_EPs2}). Then, we integrate over all the possible initial states $(M_{0},P_{0})$ to determine expressions for $\E{P_{s}}$ and $\V{P_{s}}$ for any time $s<\tau_{R}$ (Proposition~\ref{prop:EPs_and_VPs}). These expressions are dependent of the first moments of $(M_{0},P_{0})$: they depend on $\E{M_{0}}$, $\E{P_{0}}$, $\V{M_{0}}$, $\V{P_{0}}$ and $\Cov{M_{0},P_{0}}$.
\item Then we consider the case after replication ($s\geq\tau_{R}$ ). Similarly the first case, we will consider that the state of the cell at replication $(M_{\tau_{R}},P_{\tau_{R}})$ is known and we calculate the first two moments of $P_{s}$ for any time $\tau_{R}\leq s<\tau_{D}$ (Proposition~ \ref{prop:EPs_and_VPs-2}). After integration, expressions for $\E{P_{s}}$ and $\V{P_{s}}$ for any time $s$ after replication are determined, these expressions depend on $\E{M_{\tau_{R}}}$, $\E{P_{\tau_{R}}}$, $\V{M_{\tau_{R}}}$, $\V{P_{\tau_{R}}}$ and $\Cov{M_{\tau_{R}},P_{\tau_{R}}}$ (Proposition~\ref{prop:EPs_and_VPs-2}).
\end{enumerate}
In the end, in Propositions~\ref{prop:EPs_and_VPs} and~\ref{prop:EPs_and_VPs-2}, are presented the mean and variance of protein number at any time $s$ of the cell cycle, only depending on the first moments of $\left(M_{0},P_{0}\right)$ and $\left(M_{\tau_{R}},P_{\tau_{R}}\right)$. Additional results then determine explicitly the first moments of $\left(M_{0},P_{0}\right)$ and $\left(M_{\tau_{R}},P_{\tau_{R}}\right)$ so that the mean and variance of protein number will be fully characterized.


\paragraph*{Description of the  Process of the Number of Proteins  }
Before beginning, we describe the number of proteins $P_{s}$ at any time $s$. We will use this description in the following proofs. Similarly to mRNA case (\ref{eq:M_s-2}), we group them into two categories.
\begin{itemize}
\item The $P_{0}$ proteins that were there at birth and which remain in the cell during all the cell cycle (as said in the main article the proteolysis is not considered in this model).
\item The proteins that were created during the cell cycle. The rate of production depends on the current number of mRNAs. For that we consider $\left(\mathcal{N}_{\lambda_{2}}^{i}\right)_{i\in\mathbb{N}}$, a sequence of i.i.d. Poisson Point Processes of intensity $\lambda_{2}$; if the $i$-th mRNA exists at time $s$ (that is to say if $i\leq M_{s}$), then the number of proteins produced by this mRNA between $s$ and $s+\mathop{}\mathopen{}\mathrm{d} s$ is $\mathcal{N}_{\lambda_{2}}^{i}(\mathop{}\mathopen{}\mathrm{d} s)$. Hence, the total number of proteins produced between $s$ and $s+\mathop{}\mathopen{}\mathrm{d} s$ is then $\sum_{i=1}^{\infty}\ind{i\leq M_{u}}\mathcal{N}_{\lambda_{2}}^{i}\left(\mathop{}\mathopen{}\mathrm{d} s\right)$.
\end{itemize}
To summarize, the number of proteins at a time $s$ of the cell cycle is
\begin{equation}
P_{s}=P_{0}+\sum_{i=1}^{\infty}\int_{0}^{s}\ind{i\leq M_{u}}\mathcal{N}_{\lambda_{2}}^{i}\left(\mathop{}\mathopen{}\mathrm{d} u\right)\mbox{.}\label{eq:P_s-2}
\end{equation}
The first term is the number of proteins at birth, and the second
takes into account all proteins created between times $0$ and $s$.

\paragraph*{Protein Number Before Replication}
We begin  with the case before replication,  $s<\tau_{R}$. We use the notation $\Es{\cdot}{M_{0},P_{0}}$ as the conditional expectation given $(M_0,P_0)$, i.e. $\Es{\cdot}{M_{0},P_{0}}=\E{\cdot|\left(M_{0},P_{0}\right)}$.  We first characterize the first two moments of $P_{s}$ conditionally on $\left(M_{0},P_{0}\right)$.  As for the mRNAs, we determine at first the moment-generating function of $P_{s}$.
\begin{prop}
\label{prop:MGF_Ps}For any $s\in[0,\tau_{R}[$, the conditional moment generating
function of $P_{s}$ can be expressed as
\[
\Es{\exp\left(\xi P_{s}\right)}{M_{0},P_{0}}=\exp\left(\xi P_{0}\right)\cdot h_{s}\left(\lambda_{2}\left(e^{\xi}-1\right)\right)
\]
for any $\xi<0$ and such as $h_{s}$ is the moment generating function
of $\int_{0}^{s}M_{u}\,\mathop{}\mathopen{}\mathrm{d} u$. The expression of $h_{s}$ is given by
\[
h_{s}(\xi):=\exp\left[M_{0}\log\left[\frac{\sigma_{1}-\xi e^{-(\sigma_{1}-\xi)s}}{\sigma_{1}-\xi}\right]+\lambda_{1}\frac{\xi}{\sigma_{1}-\xi}\left(s-\frac{1-e^{-(\sigma_{1}-\xi)s}}{\sigma_{1}-\xi}\right)\right].
\]
\end{prop}
\begin{proof}
With Equation~(\ref{eq:P_s-2}), it is easy to show that
\[
\Es{\exp\left(\xi P_{s}\right)}{M_{0},P_{0}}  =  \exp\left(\xi P_{0}\right)\cdot\Es{\prod_{i=1}^{\infty}\Es{\exp\left(\xi\int_{0}^{s}\ind{i\leq M_{u}}\mathcal{N}_{\lambda_{2}}^{i}\left(\mathop{}\mathopen{}\mathrm{d} u\right)\right)|\left(M_{u}\right)_{u\leq s}}{M_{0},P_{0}}}{M_{0},P_{0}}.
\]
We then consider the Laplace functional of the Poisson process $\mathcal{N}_{\lambda_{2}}^{i}$,
\begin{align*}
\left.\E{\exp\left(\xi\int_{0}^{s}\ind{i\leq M_{u}}\mathcal{N}_{\lambda_{2}}^{i}\left(\mathop{}\mathopen{}\mathrm{d} u\right)\right)\right|\left(M_{u}\right)_{u\leq s}, P_0}& =  \exp\left[\lambda_{2}\int_{0}^{\infty}\left(\exp\left(\xi\ind{i\leq M_{u}}\ind{u\leq s}\right)-1\right)du\right]\\
 & =  \exp\left[\lambda_{2}\left(e^{\xi}-1\right)\int_{0}^{s}\ind{i\leq M_{u}}du\right].
\end{align*}
By making the product for $i$ from $1$ to infinity, one gets
\[
\prod_{i=1}^{\infty}\Es{\exp\left(\xi\int_{0}^{s}\ind{i\leq M_{u}}\mathcal{N}_{\lambda_{2}}^{i}\left(\mathop{}\mathopen{}\mathrm{d} u\right)\right)|\left(M_{u}\right)_{u\leq s}}{M_{0},P_{0}}=  \exp\left[\lambda_{2}\left(e^{\xi}-1\right)\int_{0}^{s}M_{u}du\right].
\]
As a consequence, it indeed follows that
\[
\Es{\exp\left(\xi P_{s}\right)}{M_{0},P_{0}}=\exp\left(\xi P_{0}\right)\cdot h_{s}\left(\lambda_{2}\left(e^{\xi}-1\right)\right).
\]
Using the expression~\eqref{eq:M_s-2} of $M_{s}$, integrated between time $0$ and $s<\tau_{R}$ gives the result. For more details of the calculations, see Chapter~3 of \cite{dessalles_stochastic_2017}).
\end{proof}
As the moment generating function of $P_{s}$ has been characterized, it is possible to deduce, by derivation, the first two moments of $P_{s}$ knowing $(M_{0},P_{0})$ for any time $s$ before the gene replication.
\begin{cor}
\label{cor:EPs_EPs2}At equilibrium, for  $s\in[0,\tau_{R}[$,
 the first two conditional moments of $P_{s}$ are given by
\begin{align*}
\Es{P_{s}}{M_{0},P_{0}} & =  P_{0}+\lambda_{2}\left(\frac{\lambda_{1}}{\sigma_{1}}s+\left(M_{0}-\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right),\\
\Es{P_{s}^{2}}{M_{0},P_{0}} & =  \left(\Es{P_{s}}{M_{0},P_{0}}\right)^{2}+M_{0}\frac{\lambda_{2}}{\sigma_{1}}\left(1-e^{-\sigma_{1}s}+\frac{\lambda_{2}}{\sigma_{1}}\left[1-e^{-\sigma_{1}s}\left(e^{-\sigma_{1}s}+2s\sigma_{1}\right)\right]\right)\\
 &  +\frac{\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\left[s\sigma_{1}-1+e^{-\sigma_{1}s}+2\frac{\lambda_{2}}{\sigma_{1}}\left(\sigma_{1}s\left(1+e^{-\sigma_{1}s}\right)-2\left(1-e^{-\sigma_{1}s}\right)\right)\right]
\end{align*}
\end{cor}

\begin{proof}
The first two moments of $P_{s}$ can be obtained by derivation of the moment generating function of  Proposition~\ref{prop:MGF_Ps},
\[
\Es{P_{s}}{M_{0},P_{0}}  =  \lim_{\xi\rightarrow 0}\frac{d}{d\xi}\left[\exp\left(\xi P_{0}\right)h_{s}\left(\lambda_{2}\left(e^{\xi}-1\right)\right)\right]
  =  P_{0}+\lambda_{2}h_{s}'\left(0\right)
\]
and
\[
\Es{P_{s}^{2}}{M_{0},P_{0}}  =  \lim_{\xi\rightarrow0}\frac{d^{2}}{d\xi^{2}}\left[\exp\left(\xi P_{0}\right)h_{s}\left(\lambda_{2}\left(e^{\xi}-1\right)\right)\right]
  =  \left(\Es{P_{s}}{M_{0},P_{0}}\right)^{2}+\lambda_{2}h_{s}'\left(0\right)+\left(\lambda_{2}\right)^{2}\left(h_{s}''\left(0\right)-h_{s}'\left(0\right)^{2}\right)
\]

The calculations of $h_{s}'\left(0\right)^{2}$, $h_{s}''\left(0\right)^{2}$
allow to show the result (see Chapter~3 of \cite{dessalles_stochastic_2017}
for the details of the calculation).
\end{proof}
The previous corollary gives expressions for $\Es{P_{s}}{M_{0},P_{0}}$ and $\Es{P_{s}^{2}}{M_{0},P_{0}}$. In the next proposition, we integrate these expressions over all birth states $\left(M_{0},P_{0}\right)$ to find formulas for $\E{P_{s}}$ and $\V{P_{s}}$ for any time $s<\tau_{R}$ before replication. These expression depends on joint moments of $M_{0}$ and $P_{0}$.
\begin{prop}
\label{prop:EPs_and_VPs}At any time $s\in[0,\tau_{R}[$ before replication, the mean and the variance of $P_{s}$ are given by
\begin{align*}
\E{P_{s}} & =  \E{P_{0}}+\lambda_{2}\left(\frac{\lambda_{1}}{\sigma_{1}}s+\left(x_{0}-\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right),\\
\V{P_{s}} & =  \V{P_{0}}+2\lambda_{2}\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\Cov{P_{0},M_{0}}+\left(\lambda_{2}\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right)^{2}x_{0}\\
  & +x_{0}\frac{\lambda_{2}}{\sigma_{1}}\left(1-e^{-\sigma_{1}s}+\frac{\lambda_{2}}{\sigma_{1}}\left[1-e^{-\sigma_{1}s}\left(e^{-\sigma_{1}s}+2s\sigma_{1}\right)\right]\right)\\
  & +\frac{\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\left[s\sigma_{1}-1+e^{-\sigma_{1}s}+2\frac{\lambda_{2}}{\sigma_{1}}\left(\sigma_{1}s\left(1+e^{-\sigma_{1}s}\right)-2\left(1-e^{-\sigma_{1}s}\right)\right)\right]
\end{align*}
where $x_{0}$ is defined in Proposition~\ref{prop:Distrib_M_0}.
\end{prop}

\begin{proof}
By considering the mean of the random variable $\E{P_{s}|(M_{0},P_{0})}$ in Corollary~\ref{cor:EPs_EPs2}, the result for
$\E{P_{s}}$ is easy to get. For the variance, consider the expression of $\E{P_{s}^{2}|(M_{0},P_{0})}$
\begin{multline*}
\E{P_{s}^{2}} =  \E{\Es{P_{s}}{M_{0},P_{0}}^{2}}+\E{M_{0}}\frac{\lambda_{2}}{\sigma_{1}}\left(1-e^{-\sigma_{1}s}+\frac{\lambda_{2}}{\sigma_{1}}\left[1-e^{-\sigma_{1}s}\left(e^{-\sigma_{1}s}+2s\sigma_{1}\right)\right]\right)\\
 +\frac{\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\left[s\sigma_{1}-1+e^{-\sigma_{1}s}+2\frac{\lambda_{2}}{\sigma_{1}}\left(\sigma_{1}s\left(1+e^{-\sigma_{1}s}\right)-2\left(1-e^{-\sigma_{1}s}\right)\right)\right]
\end{multline*}
and
\begin{multline*}
\V{P_{s}}  =  \E{\Es{P_{s}}{M_{0},P_{0}}^{2}}-\E{P_{s}}^{2}+\E{M_{0}}\frac{\lambda_{2}}{\sigma_{1}}\left(1-e^{-\sigma_{1}s}+\frac{\lambda_{2}}{\sigma_{1}}\left[1-e^{-\sigma_{1}s}\left(e^{-\sigma_{1}s}+2s\sigma_{1}\right)\right]\right)\\
   +\frac{\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\left[s\sigma_{1}-1+e^{-\sigma_{1}s}+2\frac{\lambda_{2}}{\sigma_{1}}\left(\sigma_{1}s\left(1+e^{-\sigma_{1}s}\right)-2\left(1-e^{-\sigma_{1}s}\right)\right)\right].
\end{multline*}
Now, for the expression of $\E{\Es{P_{s}}{M_{0},P_{0}}^{2}}-\E{P_{s}}^{2}$,
\begin{align*}
\E{\Es{P_{s}}{M_{0},P_{0}}^{2}}&-\E{P_{s}}^{2}  =  \E{P_{0}^{2}}+\E{\left(\lambda_{2}\left(\frac{\lambda_{1}}{\sigma_{1}}s+\left(M_{0}-\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right)\right)^{2}}\\
  & +2\E{P_{0}\times\lambda_{2}\left(\frac{\lambda_{1}}{\sigma_{1}}s+\left(M_{0}-\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right)}
  -\E{P_{0}}^{2}+\E{\lambda_{2}\left(\frac{\lambda_{1}}{\sigma_{1}}s+\left(M_{0}-\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right)}^{2}\\
  & -2\E{P_{0}}\E{\lambda_{2}\left(\frac{\lambda_{1}}{\sigma_{1}}s+\left(M_{0}-\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right)}\\
&=  \V{P_{0}}+\V{\lambda_{2}\left(\frac{\lambda_{1}}{\sigma_{1}}s+\left(M_{0}-\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right)}\\
&\hspace{2cm}   +2\Cov{P_{0},\lambda_{2}\left(\frac{\lambda_{1}}{\sigma_{1}}s+\left(M_{0}-\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}s}}{\sigma_{1}}\right)}.
\end{align*}
Finally, one just has to remark that due to Proposition~\ref{prop:Distrib_M_0}
$\E{M_{0}}=\V{M_{0}}=x_{0}$.
\end{proof}


\paragraph*{Protein Number After Replication}
For a time $s$ such as $\tau_{R}\leq s<\tau_{D}$ . We adopt a similar approach as for the previous case, the state just after replication $\left(M_{\tau_{R}},P_{\tau_{R}}\right)$ is known, and we want to determine the first two moments of $P_{s}$ for any time $s$ after the replication.
\begin{prop}
\label{prop:EPs_EPs2-2}At equilibrium, for a time $s\in[\tau_{R},\tau_{D}[$,
conditionally on the state of the cell at replication $(M_{\tau_{R}},P_{\tau_{R}})$,
the first two moments of $P_{s}$ are given by
\[
\Es{P_{s}}{M_{\tau_{R}},P_{\tau_{R}}} =  P_{\tau_{R}}+\lambda_{2}\left(2\frac{\lambda_{1}}{\sigma_{1}}\left(s-\tau_{R}\right)+\left(M_{\tau_{R}}-2\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}}{\sigma_{1}}\right),
\]
\begin{multline*}
\Es{P_{s}^{2}}{M_{\tau_{R}},P_{\tau_{R}}}  =  \left(\Es{P_{s}}{M_{\tau_{R}},P_{\tau_{R}}}\right)^{2}
   +M_{\tau_{R}}\frac{\lambda_{2}}{\sigma_{1}}\left(1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}+\frac{\lambda_{2}}{\sigma_{1}}\left[1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}\left(e^{-\sigma_{1}\left(s-\tau_{R}\right)}+2\left(s-\tau_{R}\right)\sigma_{1}\right)\right]\right)\\
   +2\frac{\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\left[\left(s-\tau_{R}\right)\sigma_{1}-1+e^{-\sigma_{1}\left(s-\tau_{R}\right)}+2\frac{\lambda_{2}}{\sigma_{1}}\left(\sigma_{1}\left(s-\tau_{R}\right)\cdot\left(1+e^{-\sigma_{1}\left(s-\tau_{R}\right)}\right)-2\left(1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}\right)\right)\right].
\end{multline*}
\end{prop}
\begin{proof}
After the replication, the rate of mRNA production is doubled, but otherwise, the dynamic is identical as it was before the replication.  One can hence easily adapt the proofs of Proposition~\ref{prop:MGF_Ps} and Corollary~\ref{cor:EPs_EPs2}, by replacing the initial state by the state at replication $(M_{\tau_{R}},P_{\tau_{R}})$, by considering that the mRNA production rate is $2\lambda_{1}$, and that the time spent since the initial state is $s-\tau_{R}$.
\end{proof}
We can then integrate the previous expressions on all possible states
at replication $(M_{\tau_{R}},P_{\tau_{R}})$. It follows that
\begin{prop}
\label{prop:EPs_and_VPs-2}At any time $s\in[\tau_{R},\tau_{D}[$
after replication, depending on joint moments of $P_{\tau_{R}}$ and
$M_{\tau_{R}}$, the mean and the variance of $P_{s}$ are given by
\begin{align*}
\E{P_{s}} & =  \E{P_{\tau_{R}}}+\lambda_{2}\left(2\frac{\lambda_{1}}{\sigma_{1}}\left(s-\tau_{R}\right)+\left(x_{\tau_{R}}-2\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}}{\sigma_{1}}\right),\\
\V{P_{s}} & =  \V{P_{\tau_{R}}}+2\lambda_{2}\frac{1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}}{\sigma_{1}}\Cov{P_{\tau_{R}},M_{\tau_{R}}}+\left(\lambda_{2}\frac{1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}}{\sigma_{1}}\right)^{2}x_{\tau_{R}}\\
  & +x_{\tau_{R}}\frac{\lambda_{2}}{\sigma_{1}}\left(1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}+\frac{\lambda_{2}}{\sigma_{1}}\left[1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}\left(e^{-\sigma_{1}\left(s-\tau_{R}\right)}+2\left(s-\tau_{R}\right)\sigma_{1}\right)\right]\right)\\
  & +2\frac{\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\left[\phantom{2\frac{\lambda_{2}}{\sigma_{1}}}\left(s-\tau_{R}\right)\sigma_{1}-1+e^{-\sigma_{1}\left(s-\tau_{R}\right)}+2\frac{\lambda_{2}}{\sigma_{1}}\left(\sigma_{1}\left(s-\tau_{R}\right)\left(1+e^{-\sigma_{1}\left(s-\tau_{R}\right)}\right)-2\left(1-e^{-\sigma_{1}\left(s-\tau_{R}\right)}\right)\right)\right],
\end{align*}
with $x_{\tau_{R}}$ as defined in Theorem~\ref{thm:xs_rep}.
\end{prop}
\begin{proof}
It is similar to the proof of Proposition~\ref{prop:EPs_and_VPs}.
\end{proof}

\paragraph*{Protein Number in the Whole Cell Cycle}

In order to have an analytic expression for the mean $\E{P_{s}}$ and variance $\V{P_{s}}$ for any time $s$ of the cell cycle, we need to have expressions for the means $\E{P_{0}}$ and $\left\langle P_{\tau_{R}}\right\rangle $, the variances $\V{P_{0}}$ and $\V{P_{\tau_{R}}}$ as well as the covariances $\Cov{P_{0},M_{0}}$ and $\Cov{P_{\tau_{R}},M_{\tau_{R}}}$.  The general idea is to use the equilibrium properties that give a relation between the distributions at birth and at division. Indeed, it gives:
\[
P_{\tau_{D}}\overset{\cal D}{=}P_{0}\quad\text{and}\quad\left(M_{\tau_{D}},P_{\tau_{D}}\right)\overset{\cal D}{=}\left(M_{0},P_{0}\right).
\]

Indeed, between times $\tau_{D}-$ and $\tau_{D}$, the proteins undergo a random partitioning, and since the system is at equilibrium, the distribution of the number of proteins after division $P_{\tau_{D}}$ is the same as the distribution of proteins at birth $P_{0}$. As a consequence:
\[
\sum_{i=1}^{P_{\tau_{D}-}}B_{i,1/2}\overset{\cal D}{=}P_{0}
\]
with $\left(B_{i,1/2}\right)$ being independent Bernoulli random
variables of parameter $1/2$ and being all independent of $P_{\tau_{D}-}$.
\begin{lem}
\label{lem:Binom_div}The mean and the variance of $P_{0}$ depend on the mean and the variance of $P_{\tau_{D}-}$ in the following way
\[
\E{P_{\tau_{D}-}}  =  2\E{P_{0}}\qquad \V{P_{\tau_{D}-}}  =  4\V{P_{0}}-2\E{P_{0}}.
\]
\end{lem}
\begin{proof}
With the moment-generating function of $P_{0}$, one gets
\[
\E{\exp\left[\xi P_{0}\right]}=\E{\prod_{1=1}^{P_{\tau_{D}-}}\E{\exp\left[B_{i,1/2}\right]}}=\E{\left(\frac{1+e^{\xi}}{2}\right)^{P_{\tau_{D}-}}}=\E{\exp\left[\log\left(\frac{1+e^{\xi}}{2}\right)P_{\tau_{D}-}\right]}
\]
As a consequence, by denoting $\eta(\xi):=\E{\exp\left[\xi P_{\tau_{D}-}\right]}$
the moment generating function of $P_{\tau_{D}-}$, it follows:
\begin{align*}
\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d}\xi}\E{\exp\left[\xi P_{0}\right]} & =  \frac{e^{\xi}}{1+e^{\xi}}\cdot\eta'\left(\log\left(\frac{1+e^{\xi}}{2}\right)\right)\\
\frac{\mathop{}\mathopen{}\mathrm{d}^{2}}{\mathop{}\mathopen{}\mathrm{d}\xi^{2}}\E{\exp\left[\xi P_{0}\right]} & =  \frac{e^{\xi}}{\left(1+e^{\xi}\right)^{2}}\cdot\eta'\left(\log\left(\frac{1+e^{\xi}}{2}\right)\right)+\left(\frac{e^{\xi}}{1+e^{\xi}}\right)^{2}\cdot\eta''\left(\log\left(\frac{1+e^{\xi}}{2}\right)\right).
\end{align*}
As $\xi$ goes to $0$, one getss
\[
\E{P_{0}}  =  \frac{\E{P_{\tau_{D}-}}}{2} \text{ and }
\E{P_{0}^{2}}  =  \frac{1}{4}\cdot\E{P_{\tau_{D}-}}+\frac{1}{4}\cdot\E{P_{\tau_{D}-}^{2}}.
\]
The lemma is proved.
\end{proof}


We then use this Lemma to calculate the means $\E{P_{0}}$ and $\left\langle P_{\tau_{R}}\right\rangle $
and the variances $\V{P_{0}}$ and $\V{P_{\tau_{R}}}$.
\begin{prop}
\label{prop:EP0_EPtauR}For $\eta=1,2$, denote,
\[
f_{\eta}(s):=\eta\frac{\lambda_{1}}{\sigma_{1}}\left(s-\tau\right)+\left(x_{\tau}-\eta\frac{\lambda_{1}}{\sigma_{1}}\right)\frac{1-e^{-\sigma_{1}\left(s-\tau\right)}}{\sigma_{1}}
\]
with $\tau=0$ in the case of $\eta=1$ (before replication) and
$\tau=\tau_{R}$ for the case $\eta=2$ (after replication). In that
case, we have that:
\[
\E{P_{0}}  =\lambda_{2}\left(f_{1}\left(\tau_{R}\right)+f_{2}\left(\tau_{D}\right)\right) \text{ and }\E{P_{\tau_{R}}}  =\lambda_{2}\left(2f_{1}\left(\tau_{R}\right)+f_{2}\left(\tau_{D}\right)\right).
\]
\end{prop}

\begin{proof}
With Propositions~\ref{prop:EPs_and_VPs} and~\ref{prop:EPs_and_VPs-2},
one gets
\[
\E{P_{\tau_{D}}}  =\E{P_{\tau_{R}}}+\lambda_{2}\left(f_{2}\left(\tau_{D}\right)\right)
  =\E{P_{0}}+\lambda_{2}\left(f_{1}\left(\tau_{R}\right)+f_{2}\left(\tau_{D}\right)\right).
\]
We conclude with the Lemma~\ref{lem:Binom_div}.
\end{proof}

\begin{prop}
\label{prop:EP0_VP0-1}For $\eta=1,2$, define
\begin{multline*}
g_{\eta}(s)  :=\left(\lambda_{2}\frac{1-e^{-\sigma_{1}\left(s-\tau\right)}}{\sigma_{1}}\right)^{2}x_{\tau}\\
  +x_{\tau}\frac{\lambda_{2}}{\sigma_{1}}\left(1-e^{-\sigma_{1}\left(s-\tau\right)}+\frac{\lambda_{2}}{\sigma_{1}}\left[1-e^{-\sigma_{1}\left(s-\tau\right)}\left(e^{-\sigma_{1}\left(s-\tau\right)}+2\left(s-\tau\right)\sigma_{1}\right)\right]\right)\\
  +\eta\frac{\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\left[\left(s-\tau\right)\sigma_{1}-1+e^{-\sigma_{1}\left(s-\tau\right)} +2\frac{\lambda_{2}}{\sigma_{1}}\left(\sigma_{1}\left(s-\tau\right)\left(1+e^{-\sigma_{1}\left(s-\tau\right)}\right)-2\left(1-e^{-\sigma_{1}\left(s-\tau\right)}\right)\right)\right].
\end{multline*}
with $\tau=0$ in the case of $\eta=1$ (before replication) and
$\tau=\tau_{R}$ for the case $\eta=2$ (after replication). In that
case, we have that:
\[
\V{P_{0}}  =\frac{1}{3}\left\{ 2\E{P_{0}}+2\frac{\lambda_{2}}{\sigma_{1}}\left[\left(1-e^{-\sigma_{1}\tau_{R}}\right)\Cov{P_{0},M_{0}}+\left(1-e^{-\sigma_{1}\left(\tau_{D}-\tau_{R}\right)}\right)\Cov{P_{\tau_{R}},M_{\tau_{R}}}\right]+g_{1}\left(\tau_{R}\right)+g_{2}\left(\tau_{D}\right)\vphantom{2\E{P_{0}}+2\frac{\lambda_{2}}{\sigma_{1}}\left[\left(1-e^{-\sigma_{1}\tau_{R}}\right)\Cov{P_{0},M_{0}}+\left(1-e^{-\sigma_{1}\left(\tau_{D}-\tau_{R}\right)}\right)\Cov{P_{\tau_{R}},M_{\tau_{R}}}\right]}\right\} .
\]
\end{prop}
\begin{proof}
By considering the expressions of Proposition~\ref{prop:EPs_and_VPs-2}
for $s=\tau_{D}-$,
\[
\V{P_{\tau_{D}-}}=\V{P_{\tau_{R}}}+2\lambda_{2}\frac{1-e^{-\sigma_{1}\left(\tau_{D}-\tau_{R}\right)}}{\sigma_{1}}\Cov{P_{\tau_{R}},M_{\tau_{R}}}+g_{2}\left(\tau_{D}\right).
\]
Similarly, the expression of Proposition~\ref{prop:EPs_and_VPs}
for $s=\tau_{R}-$ gives the expression of $\V{P_{\tau_{R}}}$ by
continuity. We have
\[
\V{P_{\tau_{D}-}}  =  \V{P_{0}}+2\frac{\lambda_{2}}{\sigma_{1}}\left[\left(1-e^{-\sigma_{1}\tau_{R}}\right)\Cov{P_{0},M_{0}}+\left(1-e^{-\sigma_{1}\left(\tau_{D}-\tau_{R}\right)}\right)\Cov{P_{\tau_{R}},M_{\tau_{R}}}\right]\\
   +g_{1}\left(\tau_{R}\right)+g_{2}\left(\tau_{D}\right).
\]
Lemma~\ref{lem:Binom_div} describes the effect of the binomial sampling
between $\tau_{D}-$ and $\tau_{D}$ on the mean and the variance
of $P$. Since, we are at equilibrium of cell cycles, one has
\[
3\V{P_{0}}  =  2\E{P_{0}}+2\frac{\lambda_{2}}{\sigma_{1}}\left[\left(1-e^{-\sigma_{1}\tau_{R}}\right)\Cov{P_{0},M_{0}}+\left(1-e^{-\sigma_{1}\left(\tau_{D}-\tau_{R}\right)}\right)\Cov{P_{\tau_{R}},M_{\tau_{R}}}\right]
   +g_{1}\left(\tau_{R}\right)+g_{2}\left(\tau_{D}\right).
\]
\end{proof}
The expression of $\V{P_{\tau_{R}}}$ can then be deduced from Proposition~\ref{prop:EPs_and_VPs}.
\begin{prop}
\label{prop:Cov_P0_M0}For $\eta=1,2$, define
\begin{align*}
k_{\eta}\left(s\right) & :=  \frac{\eta\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\E{M_{\tau}}e^{-\left(s-\tau\right)\sigma_{1}}\left(\left(s-\tau\right)\sigma_{1}-\left(1-e^{-\sigma_{1}\left(s-\tau\right)}\right)\right)\\
  & +\frac{\eta\lambda_{1}}{\sigma_{1}}\E{P_{\tau}}\left(1-e^{-\left(s-\tau\right)\sigma_{1}}\right)
   +\frac{\eta\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\E{M_{\tau}}\left(1-e^{-\left(s-\tau\right)\sigma_{1}}\right)^{2}\\
  & +\frac{\lambda_{2}}{\sigma_{1}}e^{-\left(s-\tau\right)\sigma_{1}}\left(\left(\E{M_{\tau}^{2}}-\E{M_{\tau}}\right)\left(1-e^{-\sigma_{1}\left(s-\tau\right)}\right)+\sigma_{1}\left(s-\tau\right)\E{M_{\tau}}\right)\\
 & +\frac{\eta\lambda_{1}\lambda_{2}}{\sigma_{1}^{2}}\left[\frac{\eta\lambda_{1}}{\sigma_{1}}\left(1-e^{-\left(s-\tau\right)\sigma_{1}}\right)\left(\left(s-\tau\right)\sigma_{1}-\left(1-e^{-\sigma_{1}\left(s-\tau\right)}\right)\right)+\left(1-e^{-\sigma_{1}\left(s-\tau\right)}\left(\left(s-\tau\right)\sigma_{1}+1\right)\right)\right].
\end{align*}
with $\tau=0$ in the case of $\eta=1$ (before replication) and
$\tau=\tau_{R}$ for the case $\eta=2$ (after replication). In that
case,  the covariances can be expressed as
\[
\Cov{M_{0},P_{0}}=\frac{1}{\left(4-e^{-\tau_{D}\sigma_{1}}\right)}\left\{ k_{1}\left(\tau_{R}\right)e^{-\left(\tau_{D}-\tau_{R}\right)\sigma_{1}}+k_{2}\left(\tau_{D}\right)\right\} -\E{M_{0}}\E{P_{0}}
\]
and
\[
\Cov{M_{\tau_{R}},P_{\tau_{R}}}=\left(\Cov{M_{0},P_{0}}+\E{M_{0}}\E{P_{0}}\right)e^{-\tau_{R}\sigma_{1}}+k_{1}\left(\tau_{R}\right)-\E{M_{\tau_{R}}}\E{P_{\tau_{R}}}.
\]
\end{prop}
\begin{proof}
The proof follows the same arguments as in the proof of the previous propositions. Details of the calculations can be found in Chapter~3 of \cite{dessalles_stochastic_2017}.
\end{proof}


\subsection{\label{subsec:S_mod2_param}Parameter Estimation}

As in the previous intermediate model, we set the doubling time $\tau_{D}$ to $150\,\text{min}$ and the volume at birth $V_{0}=1.3\,\text{\textmu m}^{3}$.  For each gene, we have to determine four different parameters $\lambda_{1}$, $\sigma_{1}$, $\lambda_{2}$ and $\tau_{R}$. We have considered the genes of \cite{taniguchi_quantifying_2010} for which the empirical mean of messengers $\mu_{m}$ and proteins $\mu_{p}$ concentrations, as well as the mRNA half-life time $\tau_{m}$ have been measured.  We still deduce the mRNA degradation rate $\sigma_{1}$ with the mRNA half-life time $\tau_{m}$ (such that $\sigma_{1}=\log2/\tau_{m}$).

\begin{figure}
\subfloat[\label{fig:S_mode2_TableParam}]{\includegraphics[height=4.5cm]{model2/Figure_20171110_boxParameters}
}~~
\subfloat[\label{fig:S_mod2_mRNA-prot-check}]{\begin{tikzpicture}
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[page=1,height=5cm]{model2/Figure_20170828_comparision_mod2_exp}};

\node[anchor=south west,inner sep=0] (im1) at (2,1) {

\includegraphics[height=3cm,page=2]{model2/Figure_20170828_comparision_mod2_exp}};

\end{tikzpicture}

}

\subfloat[\label{fig:S_mod2_variation-gene}]{\includegraphics[height=2.5cm]{model2/Figure_20160721_variation_on_one_gene_SI}

} ~~\subfloat[\label{fig:S_mod2_recovering_distrib}]{\includegraphics[height=5cm,page=2]{model2/Figure_20160819_variation_on_one_gene}

}\caption{Intermediate model with cell cycle and gene replication. \textbf{(A)}: Quantitative summary of the parameters for this model. \textbf{(B)}: Comparison of the average productions of protein (main) and mRNAs (inset) obtained in the simulations and those experimentally measured .\textbf{(C)}: Evolution of $\protect\Vhat{P/V}$ while varying successively the gene position in the DNA, the mRNA number and the mRNA lifetime while keeping $\protect\Ehat{P/V}$ constant.\textbf{ (D)}: Main: Profile of a modified version of AdK with higher transcription rate (approximately ten times more) and a lesser mRNA lifetime (ten times less). The variance is reduced, but it is not enough to clearly separate between the distributions at birth (at time $s=0$) and at the replication of the gene (at time $s=\tau_{R}$) (Inset).}
\end{figure}


\subsubsection{\label{subsec:AppGene-replication-times-1}A Model for the Instants of Gene Replication}

In this model, the time at which each gene is replicated is estimated as follows: we first determine the time of DNA replication initiation (the time $\tau_{I}$ in the cell cycle); as we consider that the DNA-polymerase replicates DNA at constant speed, we can deduce the time of replication of each gene only by knowing its position in the DNA.

The article \cite{wallden_fluctuations_2015} investigates the replication initiation. It is shown that the initiation occurs at a fixed volume per replication origin, and thus independently from the time since the previous division. Furthermore, this  volume seems to be constant for different conditions. For slow growing bacteria (with only one DNA replication per cell cycle), such as those in \cite{taniguchi_quantifying_2010}, the volume at which DNA replication initiation occurs is $V_{I}=1.8\,\text{\textmu m}^{3}$.  As in our model, the volume is considered as growing exponentially, we define the time of replication initiation $\tau_{I}$ as
\[
\tau_{I}=\frac{\tau_{D}}{\log2}\log\frac{V_{I}}{V_{0}}.
\]
The initiation of DNA replication occurs at $\tau_{I}$, the remaining delay to gene replication of each gene is considered as deterministic (we consider the speed of DNA replication as  constant).  The whole chromosome is replicated in around $40\,\text{min}$ \citep{grant_dnaa_2011}, therefore the distance of the gene from the origin of replication is sufficient to determine the time it takes for the DNA-polymerase to replicate it. The position of each gene was determined with Ecogene database \citep{zhou_ecogene_2013}.

\subsubsection{Estimation of $\lambda_{1}$ and $\lambda_{2}$ in an  Homogeneous  Population}

We still have to determine the rates $\lambda_{1}$ and
$\lambda_{2}$. One can interpret the empirical
average mRNA and protein concentration of the experiment (respectively
$\mu_{m}$ and $\mu_{p}$) as the global average of mRNA and protein
concentrations of the model (respectively $\Ehat{M/V}$ and $\Ehat{P/V}$).

Contrary to the previous intermediate model, the mean concentrations
$\E{M_{s}/V(s)}$ and $\E{P_{s}/V(s)}$ change during the cell cycle.
As depicted in Equation~\ref{eq:Ehat_P} of the main article, in
order to consider the concentrations $\Ehat{M/V}$ and $\Ehat{P/V}$
averaged over the cell population, one have to explicit the age distribution
$\nu$ of the population. We consider at first that the distribution
is homogeneous between age $0$ and $\tau_{D}$ (see Section~\ref{subsec:other_pop_distrib}
for a more realistic distribution). Then, the global averages are
known through the integration over the cell cycle of the mean formulas
of Theorem~\ref{thm:xs_rep} and the Propositions~\ref{prop:EPs_and_VPs}
and~\ref{prop:EPs_and_VPs-2} we can write the global average of
mRNA and protein concentrations as
\begin{align*}
\Ehat{M/V} & =  \frac{\lambda_{1}}{\sigma_{1}}\frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\frac{1}{V_{0}2^{s/\tau_{D}}}\left(1-\frac{e^{-(s+\tau_{D}-\tau_{R})\sigma_{1}}}{2-e^{-\tau_{D}\sigma_{1}}}+\ind{s\geq\tau_{R}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right)\,\mathop{}\mathopen{}\mathrm{d} s,\\
\Ehat{P/V} & =  \lambda_{2}\frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\frac{1}{V_{0}2^{s/\tau_{D}}}\left(f_{1}\left(\tau_{R}\right)+f_{2}\left(\tau_{D}\right)+f_{1}\left(\tau_{R}\wedge s\right)+\ind{s\geq\tau_{R}}f_{2}(s)\right)\,\mathop{}\mathopen{}\mathrm{d} s.
\end{align*}
As a consequence, parameters $\lambda_{1}$ and $\lambda_{2}$ can
be expressed as follows:
\begin{align*}
\lambda_{1} & =  \sigma_{1}\tau_{D}\mu_{m}\left(\int_{0}^{\tau_{D}}\frac{1}{V_{0}2^{s/\tau_{D}}}\left(1-\frac{e^{-(s+\tau_{D}-\tau_{R})\sigma_{1}}}{2-e^{-\tau_{D}\sigma_{1}}}+\ind{s\geq\tau_{R}}\left(1-e^{-(s-\tau_{R})\sigma_{1}}\right)\right)\,\mathop{}\mathopen{}\mathrm{d} s\right)^{-1},\\
\lambda_{2} & =  \tau_{D}\mu_{p}\left(\int_{0}^{\tau_{D}}\frac{1}{V_{0}2^{s/\tau_{D}}}\left(f_{1}\left(\tau_{R}\right)+f_{2}\left(\tau_{D}\right)+f_{1}\left(\tau_{R}\wedge s\right)+\ind{s\geq\tau_{R}}f_{2}(s)\right)\,\mathop{}\mathopen{}\mathrm{d} s\right)^{-1}.
\end{align*}
For each gene, all parameters can be hence determined.

\subsubsection{\label{subsec:other_pop_distrib}Impact of the Distribution of the Population of Cells}

As previously noticed, the definitions of $\Ehat{M/V}$ and $\Ehat{P/V}$ depends on the population age distribution. In real experimental populations of cells (like in \cite{taniguchi_quantifying_2010}) the number of cells in the population is exponentially growing: any dividing cell gives birth to two daughter cells. The distribution of ages is therefore not uniform.

Using a classical age distribution $\nu$ in the definitions $\Ehat{M/V}$ and $\Ehat{P/V}$ (Equation~\ref{eq:Ehat_P}) for exponentially growing populations (see \cite{collins_rate_1962,sharpe_bacillus_1998,robert_division_2014} for instance), we have performed a parameter estimation that takes into account this effect. For any gene, the protein variance is estimated in both cases: either with an uniform population or an exponentially growing population. The variances in both cases are almost identical (the histogram ratio of both variances is centered around $1$ with a standard deviation of $8\cdot10^{-3}$).

The distribution considered does  not have a significant impact on the variance of the model. This is due to the fact that the mean concentration $\V{P_{s}/V(s)}$ of any protein remains approximately constant during the cell cycle, there is therefore no significant difference of protein concentration dosage at the beginning or at the end of the cell cycle. We observe the same effect in the case of the complete model of next section (with the sharing of RNA-polymerases and ribosomes).

\section{Impact of the Sharing of RNA-Polymerases and Ribosomes}

\subsection{\label{subsec:S_mod3_Detailed_descr}A Detailed Description of the Model}

The unit of production of one particular protein is presented in Figure~\ref{fig:S_mod3_Model-schema}.
We recall that, for any time $s$, the copy number of the $i$-th
gene is $G_{i}(s)$, the number of mRNA is $M_{i}(s)$ and the number
of proteins is $P_{i}(s)$, the number of RNA-polymerases sequestered
on the $i$-th gene is $E_{Y,i}(s)$ and the number of ribosomes sequestered
on an mRNA of type $i$ is $E_{R,i}(s)$. The number of non-sequestered
RNA-polymerases and ribosomes are respectively denoted as $F_{Y}(s)$
and $F_{R}(s)$.

\begin{figure}
\begin{centering}
\begin{tikzpicture}
\linespread{1}
\matrix [column sep=20mm, row sep=10mm,ampersand replacement=\&] {
  \&
  \node[gray] (fy) [draw, shape=rectangle] {$F_Y$}; \&
  \&
  \node[gray] (fr) [draw, shape=rectangle] {$F_R$}; \&
  \\

  \node (g) [draw, shape=rectangle,blue] {$G_i$}; \&
  \node (ey) [draw, shape=rectangle,blue] {$E_{Y,i}$}; \&
  \node (m) [draw, shape=rectangle,blue] {$M_i$}; \&
  \node (er) [draw, shape=rectangle,blue] {$E_{R,i}$}; \&
  \node (p) [draw, shape=rectangle,blue] {$P_i$};\\

  \&
  \node (none1) [draw=none,fill=none,blue] {}; \&
  \node (empty) [draw=none,fill=none,blue] {$\emptyset$}; \&
  \node (none2) [draw=none,fill=none] {}; \&
  \node (none3) [draw=none,fill=none,color=red,align=center] {
	Divisions at \\the volume $2 V_0$}; \&


  \\
   \node (gj) [draw, shape=rectangle,green] {$G_j$}; \&
  \node (eyj) [draw, shape=rectangle,green] {$E_{Y,j}$}; \&
  \node (mj) [draw, shape=rectangle,green] {$M_j$}; \&
  \node (erj) [draw, shape=rectangle,green] {$E_{R,j}$}; \&
  \node (pj) [draw, shape=rectangle,green] {$P_j$}; \\
  \&
  \node (none1j) [draw=none,fill=none,green] {}; \&
  \node (emptyj) [draw=none,fill=none,green] {$\emptyset$}; \&
  \node (none2j) [draw=none,fill=none,green] {}; \&
  \node (none3j) [draw=none,fill=none,color=red,align=center] {
	Divisions at \\the volume $2 V_0$}; \\
};
\draw[->,blue!50,thick] (fy.west)  to [left,bend right=70]  (ey.west);
\draw[->,blue!50,thick] (ey.east)  to [left,bend right=70]  (fy.east);
\draw[->,blue!50,thick] (fr.west)  to [left,bend right=70]  (er.west);
\draw[->,blue!50,thick] (er.east)  to [left,bend right=70]  (fr.east);
\draw[->,green!50,thick] (fy.west)  to [left,bend right=70]  (eyj.west);
\draw[->,green!50,thick] (eyj.east)  to [left,bend right=70]  (fy.east);
\draw[->,green!50,thick] (fr.west)  to [left,bend right=70]  (erj.west);
\draw[->,green!50,thick] (erj.east)  to [left,bend right=70]  (fr.east);


\draw[->,blue] (g) -- (ey);
  \draw[->,blue] (g.east) -- +(1,0) |- (none1.center)  -- +(1.7,0) |- (m.west);
  \draw[blue] (none1.center) node[above] {$\lambda_{1,i} G_i  F_Y/V$};
\draw[blue] (ey.east)+(-0.85,.9) node[right] {$\mu_{1,i} E_{Y,i}$};
  \draw[gray] (fy.east)+(0.3,0)  node[right,above] {\tiny +1};
  \draw[gray] (fy.west)+(-0.3,0) node[left,above]  {\tiny -1};
\draw[->,blue] (m) -- (empty) node[midway,right] {$\sigma_{1,i} M_i$};
\draw[->,blue] (m) -- (er) node[above,midway] {$\lambda_{2,i} M_i  F_R /V$};
\draw[->,blue] (er) -- (p) node[above,midway] {$\mu_{2,i} E_{R,i}$};
  \draw[gray] (fr.east)+(0.3,0)  node[right,above] {\tiny +1};
  \draw[gray] (fr.west)+(-0.3,0) node[left,above]  {\tiny -1};
  \draw[blue] (m.west)+(-0.2,0) node[left,above]  {\tiny +1};
  \draw[blue] (m.south)+(-0.2,0) node[left,below]  {\tiny -1};
  \draw[blue] (p.west)+(-0.2,0) node[left,above]  {\tiny +1};

\draw[-> ,red] (m) -- (none3);
\draw[-> ,red] (p) -- (none3);

\draw[->,green] (gj) -- (eyj);
  \draw[->,green] (gj.east) -- +(1,0) |- (none1j.center)  -- +(1.7,0) |- (mj.west);
  \draw[green] (none1j.center) node[above] {$\lambda_{1,j} G_j  F_Y/V$};
\draw[green] (eyj.east)+(-0.4,.9) node[right] {$\mu_{1,j} E_{Y,j}$};
 
 
\draw[->,green] (mj) -- (emptyj) node[midway,right] {$\sigma_{1,j} M_j$};
\draw[->,green] (mj) -- (erj) node[above,midway] {$\lambda_{2,j} M_j  F_R /V$};
\draw[->,green] (erj) -- (pj) node[above,midway] {$\mu_{2,j} E_{R,j}$};
 
 
  \draw[green] (mj.west)+(-0.2,0) node[left,above]  {\tiny +1};
  \draw[green] (mj.south)+(-0.2,0) node[left,below]  {\tiny -1};
  \draw[green] (pj.west)+(-0.2,0) node[left,above]  {\tiny +1};

\draw[-> ,red] (mj) -- (none3j);
\draw[-> ,red] (pj) -- (none3j);

\end{tikzpicture}
\par\end{centering}
\caption{\label{fig:S_mod3_Model-schema}Production of production unit of the
$i$-th and $j$-th protein with the common pools of free RNA-polymerases
and ribosomes.}
\end{figure}

\begin{description}
\item [{Transcription}] In the current model, the process of mRNA production is considered as taking part in two steps: first, the binding of the RNA-polymerase and initiation; and second, the elongation and termination of the mRNA. For the first step, \emph{inside a unit volume}, the rate at which an RNA-polymerase binds on the promoter of the $i$-th gene is given by the law of mass action
\[
\lambda_{1,i}\frac{G_{i}(s)}{V(s)}\frac{F_{Y}(s)}{V(s)}.
\]
with $\lambda_{1,i}$ accounts for the specificity of the promoter (its affinity for the RNA-polymerase, the chromosome conformation, etc.). As we are interested in the rate of reactions inside the whole cell of volume $V(s)$, the rate of reaction is then
\[
\lambda_{1,i}G_{i}(s)\frac{F_{Y}(s)}{V(s)}.
\]
The elongation time is given by an exponential random variable of rate $\mu_{1,i}$. Once the elongation terminates, the RNA-polymerase is released in the cytoplasm (increasing the number of free RNA-polymerases $F_{Y}$ by one unit). A messenger is considered created as soon as its elongation begins: the reason for it is that in bacteria (unlike eukaryotes), since transcriptions and translations happen in the same medium, a translation can begin on an mRNA on which the transcription is not finished. As for the previous models, each messenger of type $i$ has a lifetime given by an exponential random variable of rate $\sigma_{1,i}$.
\item [{Translation}] Similarly to the transcription, the rate at which a ribosome encounters an mRNA of type $i$ and initiate translation is $\lambda_{2,i}M_{i}(s)F_{R}(s)/V(s)$ where $\lambda_{2,i}$ will account for mRNA specific aspects (RBS affinity for ribosomes, etc.). The total number of ribosomes sequestered on messengers of type $i$ is $E_{Y,i}(s)$ and each elongation time follows an exponential distribution of rate $\mu_{2,i}$. Here we consider that the protein is created after the termination (since the protein is usually fully functional once its translation is completed); the number of proteins $P_{i}(s)$ is then increased by one unit.  As previously we do not consider protein proteolysis since it usually occurs at much longer timescale than cell cycle.
\item [{DNA~Replication~and~Division~of~the~Cell}] At a time $s$, each gene $i\in\{1,...,K\}$ is characterized by the gene copy number $G_{i}(s)$. As previously only one DNA replication per cell cycle is considered: as a consequence then, for each $i\in\{1,...,K\}$, $G_{i}(s)$ is constant and equal to $1$ (before replication) or to $2$ (after replication). There is two modeling choice for when the DNA replication is initiated: it can occur at a fixed time after the last division or when the cell reaches a certain volume $V_{I}$.  The first simulations are made by considering the volume-dependent initiation event, but as we will see in Section~\ref{subsec:Other-deterministic-times}, simulations with the other modeling choice show no noticeable difference.  The volume $V_{I}$ is fixed to $1.8\,\text{\textmu m}^{3}$ (see \cite{wallden_fluctuations_2015} and Section~\ref{subsec:AppGene-replication-times-1} about this choice). We consider the speed of DNA replication as constant; as a consequence, once known the replication time $\tau_{I}$, the delay until the replication of $i$-th gene is fixed, and is given by the gene position.

For the division, we considered at first that, like in the previous models, the division occurs when the cell reaches exactly the volume $2V_{0}$ (with $V_{0}=1.3\,\text{\textmu m}^{3}$ as it was the case for the  intermediate models considered earlier). We will consider in Section~\ref{subsec:Uncertainty-in-divi-init} the case where the division timing is not as precise. As before, the effect of septation is a random sampling of messengers and proteins: each of them has an equal chance to be in the next considered cell or not. Moreover, at division, all genes have only one copy.
\item [{Volume~Increase}] As said in the main article, the volume $V(s)$ is no longer deterministic as it was the case in the previous intermediate models and it is considered as proportional to the current total mass of proteins in the cell. We denote by $\beta_{P}$ represents ratio mass-volume and by $w_{i}$ the mass of a type $i$ protein. In that case, we have by definition
\begin{equation}
V(s)=\sum_{i=1}^{K}w_{i}P_{i}(s)/\beta_{P}.\label{eq:volume_mass}
\end{equation}
Thus each protein of type $i$ created increases the total
volume of the cell with respect to the factor $w_{i}/\beta_{P}$.
The mass $w_{i}$ of a protein is determined according to its gene
length.
\item [{Production~of~RNA-polymerases~and~ribosomes}] The total number
of RNA-polymerases and ribosomes (whether allocated or not) are respectively
denoted by $N_{Y}(s)$ and $N_{R}(s)$. In a first step, we consider
that the both these quantities are in constant concentration, that
is to say
\[
N_{Y}(s)=\left\lfloor \beta_{Y}V(s)\right\rfloor \quad\text{and}\quad N_{R}(s)=\left\lfloor \beta_{R}V(s)\right\rfloor ,
\]
with $\beta_{Y}$ and $\beta_{R}$ constant parameters and where $\left\lfloor \centerdot\right\rfloor $ is the notation for the floor function. As the cell grows, new RNA-polymerases and ribosomes are added to the system in the corresponding proportion. When division occurs, ribosomes and RNA-polymerases will be set accordingly to the new volume. In Section~\ref{sec:S_mod3_Other-influences} we will consider the more complex case where both RNA-polymerases and ribosomes are directly produced through a gene expression process.
\end{description}

\subsection{Theoretical Analysis}

This complete model is more complex than the previous ones. It is
due in part to the feedback loop that proteins have on their own production:
the more proteins, the more the volume increases, thereby increasing
the total amount of ribosomes and hence the translation rates. This
complicates the complete analytical description of mRNA and protein
mean productions. In this section, we propose a description that mimics
the average behavior of our stochastic model: the goal is to be able
to fit parameters to real measures and use them for stochastic simulations.

\subsubsection{\label{subsec:S_mod3_DeterModel-pres}Presentation of the Deterministic
Production Model}

The description chosen to reflect the average behavior of the stochastic
model previously described is a system of ordinary differential equations
(ODEs) that describes the kinetics of each compound concentration
of the system.

We consider $K$ genes, each of them has a corresponding type of mRNA and protein. For a gene of type $i$, the concentration of gene copies is given by $g_{i}(s)$; mRNAs and protein concentrations are denoted by $m_{i}(s)$ and $p_{i}(s)$. Similarly, $f_{Y}(s)$ and $f_{R}(s)$ respectively represent the concentrations of free RNA-polymerases and free ribosomes; while $e_{Y,i}(s)$ and $e_{R,i}(s)$ denote the concentrations of RNA-polymerases and ribosomes currently sequestered to produce type $i$ proteins. All these quantities correspond to concentrations and not numbers of entities (their stochastic counterparts would be the concentrations $G_{i}(s)/V(s)$, $M_{i}(s)/V(s)$, $P_{i}(s)/V(s)$, etc.).

The reactions between different compounds are given by the \emph{law of mass action}, that is to say that the rate of chemical reaction is proportional to the reactants abundance. We will study the evolution of $m_{i}$, the concentration of mRNAs of type $i$ . The creation of a type $i$ mRNA is the result of a reaction between a free RNA-polymerase (whose concentration is $f_{Y}(s)$) and the gene $i$ (whose concentration is $g_{i}(s)$); $\lambda_{1,i}$ is interpreted as the affinity constant of the reaction. The type $i$ mRNA degradation is the result of a reaction that occurs at rate $\sigma_{1,i}$.

As in the usual description of the cell (see \cite{goelzer_cell_2011} for instance), one also must consider the dilution: without any molecule creation, the concentration of the compound still decreases as the volume grows due to dilution. If we consider that the cell volume is growing exponentially, doubling of volume in a time $\tau_{D}$, then the rate of dilution is $\log2/\tau_{D}$. The exponential growth corresponds to the volume dynamics of real bacteria \citep{wang_robust_2010}, and we will see in Section~\ref{subsec:DeterModel-validation} that it is a good approximation of the growth of cells in stochastic simulations.

All these aspects considered altogether, the kinetics of the concentration
of mRNAs of type $i$ is given by the ODE:
\begin{equation}
\frac{\mathop{}\mathopen{}\mathrm{d} m_{i}}{\mathop{}\mathopen{}\mathrm{d} s}(s)=\lambda_{1,i}g_{i}(s)\cdot f_{Y}(s)-\sigma_{1,i}m_{i}(s)-\frac{\log2}{\tau_{D}}\cdot m_{i}(s).\label{eq:sys_mi}
\end{equation}
The first term represents the mRNA creation; the second, the mRNA
degradation; and the last, the dilution.

For the other reactions, for $i\in\{1,...,K\}$, one has
\begin{align}
\frac{\mathop{}\mathopen{}\mathrm{d} p_{i}}{\mathop{}\mathopen{}\mathrm{d} s}(s) & =\mu_{2,i}e_{R,i}(s)-\frac{\log2}{\tau_{D}}\cdot p_{i}(s),\label{eq:sys_pi}\\
\frac{\mathop{}\mathopen{}\mathrm{d} e_{Y,i}}{\mathop{}\mathopen{}\mathrm{d} s}(s) & =\lambda_{1,i}g_{i}(s)\cdot f_{Y}(s)-\mu_{1,i}e_{Y,i}(s)-\frac{\log2}{\tau_{D}}\cdot e_{Y,i}(s),\label{eq:sys_eyi}\\
\frac{\mathop{}\mathopen{}\mathrm{d} e_{R,i}}{\mathop{}\mathopen{}\mathrm{d} s}(s) & =\lambda_{2,i}m_{i}(s)\cdot f_{R}(s)-\mu_{2,i}e_{R,i}(s)-\frac{\log2}{\tau_{D}}\cdot e_{R,i}(s).\label{eq:sys_eri}
\end{align}

As for the stochastic model of the previous section, assume that the concentration of RNA-polymerases (allocated or not) is constant and equal to $\beta_{Y}$, i.e.
\begin{equation}
\beta_{Y}=f_{Y}(s)+\sum_{i=1}^{K}e_{Y,i}(s),\label{eq:sys_fy}
\end{equation}
since $\sum_{i}e_{Y,i}$ and $f_{Y}$ represent the concentrations
of respectively the allocated and non-allocated RNA-polymerases. It
is similar to the ribosomes as we have:

\begin{equation}
\beta_{R}=f_{R}(s)+\sum_{i=1}^{K}e_{R,i}(s).\label{eq:sys_fr}
\end{equation}

The classical strategy in literature to study such system (an analogous model is presented in \cite{borkowski_translation_2016}) is to consider the system in steady state growth: the gene concentration $g_{i}$ is considered as constantly equal to its average value during the cell cycle, and then one can calculate the concentrations of $m_{i}$, $p_{i}$, $e_{Y,i}$ and $e_{R,i}$ at steady state by writing the Equations~\ref{eq:sys_mi} to~\ref{eq:sys_fr} with the derivative term as null. Using such method to determine parameters are unfortunatly not precise enough: there is a clear shift between the stochastic protein concentration and the one that shoule be obtained.

We have described the cell during one cycle with a non-constant gene concentration. The instant  of replication of gene $i$ replication within the cycle  is denoted by $\tau_{R,i}$. in particular, at time $s$, the $i$-th gene copy number is known: $g_{i}(s)=\left(1+\ind{s\geq\tau_{R,i}}\right)/(V_{0}2^{s/\tau_{D}})$ (the factor $V_{0}2^{s/\tau_{D}}$ represents the volume). By analogy with the equilibrium condition presented in Section~\ref{subsec:S_mod1_mrna_dynamic}, it is likely that  a large number of cell cycles have already occurred, so that the concentration of any entities is the same at the beginning and at the end of the cell cycle. For each unit of production, the concentrations $m_{i}$, $p_{i}$, $e_{Y,i}$ and $e_{R,i}$ are such as
\begin{equation}
\forall i\in\{1,...,K\}\quad\begin{cases}
p_{i}(0)=p_{i}(\tau_{D}),\quad & m_{i}(0)=m_{i}(\tau_{D}),\\
e_{Y,i}(0)=e_{Y,i}(\tau_{D}),\quad & e_{R,i}(0)=e_{R,i}(\tau_{D}).
\end{cases}\label{eq:bondery-conditions}
\end{equation}
With these considerations, we have a system of ODEs to describe the average behavior of the main stochastic model during the cell cycle. In the next section, under some simplifications, we propose to give expressions for $m_{i}(s)$, $p_{i}(s)$, $e_{Y,i}(s)$, $e_{R,i}(s)$, $f_{Y}(s)$ and $f_{R}(s)$ as a function of all  parameters ($\lambda_{1,i}$, $\sigma{}_{1,i}$, etc.) and $g_{i}(s)$.

\subsubsection{\label{subsec:S_mod3_DeterModel-dyn}Dynamics of the Average Production Model}

In order to estimate the parameters, one needs to have expressions
for $m_{i}$, $e_{Y,i}$, $p_{i}$, $e_{R,i}$, $f_{R}$ and $f_{Y}$
of the previous ODEs for any time $s$ of the cell cycle. But the
interdependence between $e_{Y,i}$ and $f_{Y}$ on one hand and $e_{R,i}$
and $f_{R}$ on the hand raises difficulties when integrating these
equations. Explicit solution for the dynamics $m_{i}$, $e_{Y,i}$,
$p_{i}$, $e_{R,i}$, $f_{R}$ and $f_{Y}$ are therefore not easy
to obtain directly.

In order to have expressions for these quantities, we have chosen to   some simplifications. In the next sections, the stochastic simulations show a good correspondence between their average concentration of free RNA-polymerase and ribosomes and the ones predicted here; it will therefore justify \emph{a~posteriori }the simplifications that we make in this section.

For the RNA-polymerases, we denote by $\widetilde{\mu_{1}}:=\sum_{i}\mu_{1,i}/K$
the average elongation rates of transcription and the function $h$ such
as
\[
h(s):=\sum_{i=1}^{K}e_{Y,i}(s)\frac{\widetilde{\mu_{1}}}{\mu_{1,i}}.
\]
The dynamic of $h$ is given by summing the equations \ref{eq:sys_eyi}
for $i$ from $1$ to $K$, and by using Equation~(\ref{eq:sys_fy}):
\begin{equation}
\frac{d}{\mathop{}\mathopen{}\mathrm{d} s}h(s)=f_{Y}(s)\cdot\widetilde{\mu_{1}}\left(1+\sum_{i=1}^{K}\frac{\lambda_{1,i}}{\mu_{1,i}}g_{i}(s)\right)-\beta_{Y}\widetilde{\mu_{1}}-\frac{\log2}{\tau_{D}}\cdot h(s).\label{eq:d_ey_i}
\end{equation}

The $h$ is simply a weighted sum of the $e_{Y,i}$ allocated RNA-polymerases. We decided to consider that such weighting has little influence, and that $h$ does not greatly differ from the uniform sum $\sum_{i}e_{Y,i}$, that is to say:
\[
h(s)=\sum_{i=1}^{K}e_{Y,i}(s)\frac{\widetilde{\mu_{1}}}{\mu_{1,i}}\quad\simeq\quad\sum_{i=1}^{K}e_{Y,i}(s)=\beta_{Y}-f_{Y}(s).
\]
It would be in particular true if all elongation rates $\mu_{1,i}$
are identical for all genes (i.e. if $\mu_{1,i}\equiv\widetilde{\mu_{1}}$
for all $i$).

With this simplification, from Relation~\eqref{eq:d_ey_i}, one obtains a differential equation for  $f_{Y}$
\begin{equation}
\frac{\mathop{}\mathopen{}\mathrm{d}}{\mathop{}\mathopen{}\mathrm{d} s}f_{Y}(s)=\widetilde{\mu_{1}}\beta_{Y}\left(\frac{\log2}{\widetilde{\mu_{1}}\tau_{D}}+1\right)-\widetilde{\mu_{1}}\left(1+\frac{\log2}{\widetilde{\mu_{1}}\tau_{D}}+\sum_{i=1}^{K}\frac{\lambda_{1,i}}{\mu_{1,i}}g_{i}(s)\right)f_{Y}(s).\label{eq:d_fy}
\end{equation}
One can remark that the concentrations of free RNA-polymerases is on a quick timescale. Indeed, as there are of the order of  $1.4\times10^{3}$ mRNAs in the cell, see~\cite{neidhardt_chemical_1996} that last approximately $4$ minutes, see~\cite{taniguchi_quantifying_2010}, it gives of the order of  $6$ translations per second. As a consequence, one can expect that that $f_{Y}$ quickly reaches its equilibrium during the cell cycle. This consideration will be justified \emph{a~posteriori} by the agreement with  stochastic simulations.

With these considerations, we set  the derivative term of Equation~\eqref{eq:d_fy} to be null, hence
\[
f_{Y}(s)=\beta_{Y}\frac{{\displaystyle 1+\frac{\log2}{\widetilde{\mu_{1}}\tau_{D}}}}{{\displaystyle \sum_{i=1}^{K}\frac{\lambda_{1,i}}{\mu_{1,i}}g_{i}(s)+1+\frac{\log2}{\widetilde{\mu_{1}}\tau_{D}}}}.
\]

In the next section it will be shown that $\log2/(\widetilde{\mu_{1}}\times\tau_{D})\sim10^{-3}\ll1$, we will therefore neglect the contribution of this term.
 With a similar argument for free ribosomes, we get
\[
f_{Y}(s)=\beta_{Y}\frac{{\displaystyle 1}}{1+{\displaystyle \sum_{i=1}^{K}\frac{\lambda_{1,i}}{\mu_{1,i}}g_{i}(s)}}\quad\text{and}\quad f_{R}(s)=\beta_{R}\frac{{\displaystyle 1}}{1+{\displaystyle \sum_{i=1}^{K}\frac{\lambda_{2,i}}{\mu_{2,i}}m_{i}(s)}}.
\]

With global quantities $f_{Y}$ and $f_{R}$ known, we are able
to give expression for gene-specific variables. For each $i\in\{1,...,K\}$,
one can integrate Equation~(\ref{eq:sys_mi}) and find that:
\[
\frac{\mathop{}\mathopen{}\mathrm{d} m_{i}}{\mathop{}\mathopen{}\mathrm{d} s}(s)=\lambda_{1,i}g_{i}(s)\cdot f_{Y}(s)-\sigma_{1,i}m_{i}(s)-\frac{\log2}{\tau_{D}}\cdot m_{i}(s).
\]
With the boundary conditions of Equation~(\ref{eq:bondery-conditions}),
it is easy to deduce that:
\begin{equation}
m_{i}(s)=\lambda_{1,i}\frac{e^{-\sigma_{1,i}s}}{2^{s/\tau_{D}}}\left[\int_{0}^{s}2^{u/\tau_{D}}e^{\sigma_{1,i}u}g_{i}(u)f_{Y}(u)\,\mathop{}\mathopen{}\mathrm{d} u+\frac{\int_{0}^{\tau_{D}}2^{u/\tau_{D}}e^{\sigma_{1,i}u}g_{i}(u)f_{Y}(u)\,\mathop{}\mathopen{}\mathrm{d} u}{2e^{\sigma_{1,i}\tau_{D}}-1}\right].\label{eq:mi_final}
\end{equation}
Since the quantities $g_{i}$ , $f_{Y}$ are known, we have an explicit
solution for $m_{i}$.

Similarly for $e_{Y,i}(s)$ and $e_{R,i}(s)$,
\begin{align*}
e_{Y,i}(s) & =\lambda_{1,i}\frac{e^{-\mu_{1,i}s}}{2^{s/\tau_{D}}}\left[\int_{0}^{s}2^{u/\tau_{D}}e^{\mu_{1,i}u}g_{i}(u)f_{Y}(u)\,\mathop{}\mathopen{}\mathrm{d} u+\frac{\int_{0}^{\tau_{D}}2^{u/\tau_{D}}e^{\mu_{1,i}u}g_{i}(u)f_{Y}(u)\,\mathop{}\mathopen{}\mathrm{d} u}{2e^{\mu_{1,i}\tau_{D}}-1}\right],\\
e_{R,i}(s) & =\lambda_{2,i}\frac{e^{-\mu_{2,i}s}}{2^{s/\tau_{D}}}\left[\int_{0}^{s}2^{u/\tau_{D}}e^{\mu_{2,i}u}m_{i}(u)f_{R}(u)\,\mathop{}\mathopen{}\mathrm{d} u+\frac{\int_{0}^{\tau_{D}}2^{u/\tau_{D}}e^{\mu_{2,i}u}m_{i}(u)f_{R}(u)\,\mathop{}\mathopen{}\mathrm{d} u}{2e^{\mu_{2,i}\tau_{D}}-1}\right].
\end{align*}
Consider now the type $i$ protein concentration. By integrating
the Equation~(\ref{eq:sys_mi}), and by considering the boundary
condition of Equation~\ref{eq:bondery-conditions}, one gets the relation
\begin{equation}
p_{i}(s)=\frac{\mu_{2,i}}{2^{s/\tau_{D}}}\int_{0}^{\tau_{D}}\left(1+\ind{u<s}\right)2^{u/\tau_{D}}e_{R,i}(u)\,\mathop{}\mathopen{}\mathrm{d} u.\label{eq:pi_final}
\end{equation}

As in the previous models, we are interested in average concentrations
over the cell cycle. Since, in the system of ODEs, we define average
concentrations over the cell cycle of free RNA-polymerases and ribosomes
respectively as
\begin{equation}
\overline{f_{Y}}=\frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\beta_{Y}\frac{1}{{\displaystyle \sum_{i=1}^{K}\frac{\lambda_{1,i}}{\mu_{1,i}}g_{i}(s)}+1}\,\mathop{}\mathopen{}\mathrm{d} s\quad\text{and}\quad\overline{f_{R}}=\frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\beta_{R}\frac{1}{{\displaystyle \sum_{i=1}^{K}\frac{\lambda_{2,i}}{\mu_{2,i}}m_{i}(s)+}1}\,\mathop{}\mathopen{}\mathrm{d} s.\label{eq:fy-fr}
\end{equation}
We defined similarly the concentrations $\overline{m_{i}}$ and $\overline{p_{i}}$
averaged over the cell cycle. By integrating Equations~(\ref{eq:pi_final})
and~(\ref{eq:pi_final}), it follows:

\begin{equation}
\overline{m_{i}}=\frac{\lambda_{1,i}}{\sigma_{1,i}\tau_{D}+\log2}\int_{0}^{\tau_{D}}g_{i}(u)f_{Y}(u)\,\mathop{}\mathopen{}\mathrm{d} u\quad\text{and}\quad\overline{p_{i}}=\frac{\lambda_{2,i}\mu_{2,i}\tau_{D}}{\log2\left(\mu_{2,i}\tau_{D}+\log2\right)}\int_{0}^{\tau_{D}}m_{i}(u)f_{R}(u)\,\mathop{}\mathopen{}\mathrm{d} u.\label{eq:av-mi-pi}
\end{equation}

Now we have expressions of the average concentrations of $\overline{m_{i}}$,
$\overline{p_{i}}$, $\overline{f_{R}}$ and $\overline{f_{Y}}$ for
any time $s$ in the cell cycle that will be used in the next subsection
to determine the parameters.

\subsection{\label{subsec:S_mod3_param}Estimation of Parameters}
The stochastic model of this section are used to describe the production of all proteins of the cell. Recall that \cite{taniguchi_quantifying_2010} has only considered $1018$ genes, out of which only $841$ have their mRNA production measured. In a first step, we only take into account the $841$ genes with protein and mRNA production measured and consider that it would represent the whole genome; in Section~\ref{subsec:S_Other-Additional-genes} we will study the case of a simulation with a complete set of genes representing a full genome of about $2000$ genes.

The determination of the model parameters $\sigma_{1,i}$ of mRNA
degradation of type $i$ , of the doubling time $\tau_{D}$ and the
time $\tau_{R,i}$ of gene replication is the same as for the previous
intermediate models (see Section~\ref{subsec:S_mod1_param} and~\ref{subsec:S_mod2_param}).

We still need to determine all reaction rates for every protein type ($\lambda_{1,i}$, $\mu_{1,i}$, $\lambda_{2,i}$ and $\mu_{2,i}$ for $i\in\{1,...,K\}$) as well as concentration parameters of RNA-polymerases, and ribosomes (respectively $\beta_{Y}$ and $\beta_{R}$ ), the proportion between the volume and the protein mass $\beta_{P}$, the mass of each proteins $w_{i}$ and the copy number $g_{i}$ of any gene.

\cite{taniguchi_quantifying_2010} does not give the quantities of non-allocated RNA-polymerases or ribosomes. To determine the set of parameters, we fix the average concentration of free RNA-polymerases and ribosomes. Note that we can have multiple sets of parameters depending on this choice.  In the simulations, we will examine several simulations with different values for average free RNA-polymerase and ribosome concentrations to see their impact on the dynamic of the model (Section~\ref{sec:Impact-RNAP-ribo}).

The rates $\mu_{1,i}$, $\mu_{2,i}$ of mRNAs and protein elongation rates can be deduced from the gene length of the $i$-th gene. In the description of the model, we have considered that the length of the mRNA is characterized by its length; so a rate the parameter $\mu_{1,i}$ is given by the mRNA elongation speed ($39\,\text{Nucl}/\text{s}$ in \cite{bremer_modulation_1996} for slowly growing cells) divided by the length of the $i$-th gene. Similarly $\mu_{1,i}$ is given by the protein elongation speed ($12\,\text{aa}/\text{s}$ in \cite{bremer_modulation_1996} for slowly growing cells) divided by the number of amino-acids coded by the $i$-th gene divided . The mass of each protein $w_{i}$ is also deduced from the length of the gene as it determines the number of amino-acids of the protein.

What remains to determine are the concentration parameters of RNA-polymerases, and ribosomes ($\beta_{Y}$ and $\beta_{R}$), the proportion between the volume and the mass of proteins $\beta_{P}$, as well as the activities of the gene and the mRNA (respectively $\lambda_{1,i}$ and $\lambda_{2,i}$) in each unit of production $i\in\{1,...,K\}$. To do so, we interpret the mRNA and protein concentration of each type measured in \cite{taniguchi_quantifying_2010} as the average concentration of each mRNA and proteins over the cell cycle of this model (respectively $\overline{m_{i}}$ and $\overline{p_{i}}$).  Moreover, as previously said, the average concentrations of free RNA-polymerases $\overline{f_{Y}}$ and free ribosomes $\overline{f_{R}}$ are fixed.

We want now to compute $\beta_{P}$, $\beta_{Y}$, $\beta_{R}$, $\lambda_{1,i}$ and $\lambda_{2,i}$ based on known values for $\overline{f_{Y}}$, $\overline{f_{R}}$, $\overline{m_{i}}$ and $\overline{p_{i}}$.  We determine the parameter $\beta_{P}$. In the description of the stochastic model, Equation~(\ref{eq:volume_mass}) states that at any moment, the volume is considered to be proportional to the total mass of proteins. Interpreting $\overline{p_{i}}$ as the average concentration of the protein of type $i$ leads by integration of Equation~(\ref{eq:volume_mass}) to
\[
\beta_{P}=\sum_{i=1}^{K}w_{i}\overline{p_{i}}.
\]

We continue with the parameters relevant to the transcription:
$\lambda_{1,i}$ and $\beta_{Y}$. With Equations~\ref{eq:fy-fr}
and~\ref{eq:av-mi-pi}, $\beta_{Y},\lambda_{1,1},...,\lambda_{1,K}$
are solution of the system
\begin{equation}
\begin{cases}
\beta_{Y} & =\overline{f_{Y}}\left({\displaystyle \frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}{\displaystyle \left(\sum_{i=1}^{K}\frac{\lambda_{1,i}}{\mu_{1,i}}g_{i}(s)+1\right)^{-1}}\,\mathop{}\mathopen{}\mathrm{d} s}\right)^{-1}\\
\lambda_{1,i} & =\overline{m_{i}}\cdot{\displaystyle \left(\sigma_{1,i}\tau_{D}+\log2\right)\cdot\left({\displaystyle \int_{0}^{\tau_{D}}g_{i}(u)f_{Y}(u)\,\mathop{}\mathopen{}\mathrm{d} u}\right)^{-1}}\quad\forall i\in\{1,...,K\}.
\end{cases}\label{eq:fixe_point1}
\end{equation}
Since $\overline{f_{Y}}$, $\overline{m_{i}}$ and $g_{i}(s)$ have
already been settled, we can use a fixed point optimization procedure
to determine $\beta_{Y}$ and all $\lambda_{1,i}$. Then, as these
parameters are determined, we now have an explicit expression for
$f_{Y}(s)$ for any time $s$ of the cell cycle.

We have to determine the parameters relevant to translation, namely $\lambda_{2,i}$ and $\beta_{R}$. Here again, we use a fixed point optimization procedure to deliver the result. With Equation~(\ref{eq:fy-fr}) and the expression of $\overline{p_{i}}$ in Equation~(\ref{eq:av-mi-pi}), $\beta_{R},\lambda_{2,1},...,\lambda_{2,K}$ are solutions of the system
\begin{equation}
\begin{cases}
\beta_{R} & =\overline{f_{R}}\times\left({\displaystyle \frac{1}{\tau_{D}}\int_{0}^{\tau_{D}}\left(\sum_{i=1}^{K}\frac{\lambda_{2,i}}{\mu_{2,i}}m_{i}(s)+1\right)^{-1}\,\mathop{}\mathopen{}\mathrm{d} s}\right)^{-1}\\
\lambda_{2,i} & =\overline{p_{i}}\times\left({\displaystyle {\displaystyle \frac{\mu_{2,i}\tau_{D}}{\log2\left(\mu_{2,i}\tau_{D}+\log2\right)}\int_{0}^{\tau_{D}}m_{i}(u)f_{R}(u)}\,\mathop{}\mathopen{}\mathrm{d} u}\right)^{-1}\quad\forall i\in\{1,...,K\}.
\end{cases}\label{eq:fixe_point2}
\end{equation}
By fixing the average amount of free RNA-polymerases and ribosomes, it is possible, through this procedure to determine  parameters with  the experimental measures.
\subsection{\label{subsec:DeterModel-validation}Validation of the Average Production Model}
\begin{figure}
\subfloat[\label{tab:S_mod3_Quantitative-summary}]{\includegraphics[height=5cm]{model3/Figure_20171110_boxParameters}}

\subfloat[\label{fig:S_mod3_mRNA-prot-check}]{\begin{tikzpicture}
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[height=5cm,page=1]{model3/Figure_20170201_firstfig_alpha_fr_fy}};

\node[anchor=south west,inner sep=0] (im1) at (3,.75) {

\includegraphics[page=2,height=2.7cm]{model3/Figure_20170201_firstfig_alpha_fr_fy}};

\end{tikzpicture}

}~~\subfloat[\label{fig:S_mod3_exp-growth}]{\includegraphics[height=5cm,page=3]{model3/Figure_20170201_firstfig_alpha_fr_fy}}

\subfloat[\label{fig:S_mod3_Frnap-profile}]{\includegraphics[height=2.75cm,page=1]{model3/Figure_20160912_firstfig_globals_sim00}

}~~\subfloat[\label{fig:S_mod3_Fribo-profile}]{\includegraphics[height=2.75cm,page=2]{model3/Figure_20160912_firstfig_globals_sim00}}

\caption{Complete model.\label{fig:Correspondance}\textbf{(A)}: Quantitative summary of the parameters. (Different choice of $\overline{f_{Y}}$ and $\overline{f_{R}}$ when computing the parameters induce little changes for the rate of transcription per gene  $\lambda_{1,i}\overline{f_{Y}}$ and the rate of translation per mRNA  $\lambda_{2,i}\overline{f_{R}}$). \textbf{(B)}: Ratio between the average concentration for protein (main figure) and mRNA (inset) in simulation and in experiments. \textbf{(C)}: a simulation sample that shows that cell volume grows exponentially from around $V_{0}$ up to around $2V_{0}$ (inset); the growth rate distribution is centered around the expected growth rate $\log2/\tau_{D}$ (main figure). \textbf{(D)} and \textbf{(E): }the respective means of free RNA-polymerases and ribosomes at each moment of the cell cycle in the simulations (solid lines) and the ones predicted by the system of ODEs (dashed lines). Inset: an example of the dynamics of free RNA-polymerases and ribosomes for one simulation.}
\end{figure}

The description of the average production through the system of ODE (Section~\ref{subsec:S_mod3_DeterModel-pres}) makes the computation of parameters of the stochastic model possible. We need to check that the deterministic description globally corresponds to the average behavior of the stochastic model; for instance, one has to validate that stochastic simulations with the  parameters previously determined, are consistent with the number  of mRNAs and proteins observed.

Here, we present the results of a particular simulation, whose parameters are presented  in Figure~\ref{tab:S_mod3_Quantitative-summary}. Its average behavior will be compared with the expressions derived from the system of ODEs. The simulation presented here takes into account  the $841$ genes with protein and mRNA production described in \cite{taniguchi_quantifying_2010}, and we have fixed the number of free RNA-polymerases and ribosomes in order to compute the parameters.

The system of ODEs assumes that the volume growth is exponential with rate $\log2/\tau_{D}$. In Figure~\ref{fig:S_mod3_exp-growth}, the volume of the cell indeed seems to grow exponentially in the simulations; the growth rate corresponds to the expected a doubling time of $\tau_{D}$.

For each type of gene, Figure~\ref{fig:S_mod3_mRNA-prot-check} shows the ratio between the protein production observed in the simulations divided by the protein production expected (and similarly for the mRNAs in inset). It appears that the correspondence is correct, especially for the highly expressed proteins. It is less precise for the protein less expressed but, globally, the correspondence  seems good enough.

Computed from the stochastic simulations, the main Figures~\ref{fig:S_mod3_Frnap-profile} and~\ref{fig:S_mod3_Fribo-profile} present the mean number of free RNA-polymerases and ribosomes as a function of the cell volume. The mean of each free entity is not constant during the cell cycle. The dashed lines represent the expected value of free entities given by the model of ODEs (Equation~(\ref{eq:fy-fr})). It is indeed a good approximation for the behavior of free RNA-polymerases and ribosomes.  The stochastic simulation displays relative quick timescale for the evolution of free RNA-polymerases (of the order of the second) and even quicker for the free ribosomes (insets of Figures~\ref{fig:S_mod3_Frnap-profile} and~\ref{fig:S_mod3_Fribo-profile}).

All these results support the idea that the expressions derived from the system of ODEs are accurate to describe the average behavior of the stochastic model.
\subsection{\label{sec:Impact-RNAP-ribo}Impact of Free RNA-polymerases and Ribosomes}
As in Section~\ref{subsec:S_mod3_DeterModel-pres}, the parameter computation supposes that  the average concentrations of free RNA-polymerases $\overline{f_{Y}}$ and ribosomes $\overline{f_{R}}$ are fixed.  In Figure~\ref{fig:S_mod3_influence_fRNAP_fRibo}, are presented several simulations where the average concentrations of these free entities are changed.

\subsubsection{\label{subsec:F-RNAP-and_ribo}Few Free Ribosomes and Many Free RNA-polymerases}

The first simulation, corresponding to Figures~\ref{fig:S_mod3_sim_11} and~\ref{fig:S_mod3_sim_12}, considers a low concentration of free ribosomes and a high concentration of free RNA-polymerases. This situation seems to reasonable in this biological setting.   Consequently, free ribosomes and free polymerases are subject to a large competition between transcripts (see \cite{warner_economics_2001} in the case of the yeast). At the same time, the parameters are fixed  in such a way that most of  RNA-polymerases are non-allocated, i.e. not bound on the DNA, see~\cite{klumpp_growth-rate-dependent_2008}, see Figure~\ref{fig:mod3} of the article and  Figure~\ref{tab:S_mod3_Quantitative-summary}.  The Figure~\ref{fig:S_mod3_sim_11} compares the variance of the complete stochastic model, with the one predicted by the previous model (with gene replication and random partitioning); it shows that for $90\%$ of the genes, the protein variance ratio is above $0.9$ (the mean of the ratio is $0.96$).

In these simulations, we look at the distributions of free RNA-polymerases and ribosomes. In Figure~\ref{fig:S_mod3_sim_12}, we show these distributions at three different phases in the cell cycle: we have selected cells of a given volume, either $1.40\,\text{\textmu m}^{3}$,$1.95\,\text{\textmu m}^{3}$ or $2.50\,\text{\textmu m}^{3}$, which correspond to the beginning, middle and end of the cell cycle. These distributions change as the volume increase (so that the average follows the curves shown in Figures~\ref{fig:S_mod3_Frnap-profile} and~\ref{fig:S_mod3_Fribo-profile}).

In order to interpret the observed distributions of free RNA-polymerases and ribosomes at a certain extend, we can propose a simplified model of RNA-polymerase and ribosome allocation for which  translation and the translation are then considered separately, and that there is no notion of cell growth. The idea would be to approach the ``local'' equilibrium of RNA-polymerases and ribosomes before any significant change in the volume (for more details about the simplified model, see Section~\ref{subsec:S_mod3_simplified_models}).

This simplified description predicts that for a given volume $V$, the distribution of free RNA-polymerases and ribosomes would be both a binomial distribution. These predicted binomial distributions are plotted in Figure~\ref{fig:S_mod3_sim_12} in thick lines. In the RNA-polymerase case, the binomial distribution globally fit the histograms.  The ribosome distribution is  singular: the parameters of the binomial distribution $\left(N,\phi\right)$ are such that $\phi\ll N$. It is due to the low concentration of free ribosomes chosen for the parameters computation. But even this denatured case shows a good correspondence between the binomial distribution and the simulation histograms.

\begin{figure}
\subfloat[\label{fig:S_mod3_sim_11}]{
\includegraphics[height=5cm,page=2]{model3/Figure_20160912_MCDivision_vs_MSGeneReplicationSI}
}~~\subfloat[\label{fig:S_mod3_sim_12}]{\begin{tikzpicture}
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[page=4,height=2.4cm]{model3/Figure_20170201_firstfig_alpha_fr_fy}};

\node[anchor=south west,inner sep=0] (im1) at (0,-2.4) {\includegraphics[page=5,height=2.4cm]{model3/Figure_20170201_firstfig_alpha_fr_fy}};

\end{tikzpicture}}

\subfloat[\label{fig:S_mod3_sim_21}]{
\includegraphics[height=5cm,page=3]{model3/Figure_20160912_MCDivision_vs_MSGeneReplicationSI}}~~\subfloat[\label{fig:S_mod3_sim_22}]{\begin{tikzpicture}
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[page=3,height=2.4cm]{model3/Figure_20160912_firstfig_alpha_fr_fy_sim04}};

\node[anchor=south west,inner sep=0] (im1) at (0,-2.4) {\includegraphics[page=4,height=2.4cm]{model3/Figure_20160912_firstfig_alpha_fr_fy_sim04}};

\end{tikzpicture}}

\subfloat[\label{fig:S_mod3_sim_31}]{\begin{centering}
\includegraphics[height=5cm,page=4]{model3/Figure_20160912_MCDivision_vs_MSGeneReplicationSI}
\par\end{centering}
}~~\subfloat[\label{fig:S_mod3_sim_32}]{\begin{tikzpicture}
\node[anchor=south west,inner sep=0] (im) at (0,0) {\includegraphics[page=3,height=2.4cm]{model3/Figure_20160912_firstfig_alpha_fr_fy_sim11}};

\node[anchor=south west,inner sep=0] (im1) at (0,-2.4) {\includegraphics[page=4,height=2.4cm]{model3/Figure_20160912_firstfig_alpha_fr_fy_sim11}};

\end{tikzpicture}

}\caption{\label{fig:S_mod3_influence_fRNAP_fRibo}Three simulations (1st line: few free ribosomes, many free RNA-polymerases; 2nd line: few free ribosomes, few free RNA-polymerases; 3rd line: many free ribosomes, many free RNA-polymerases). \textbf{(A)}, \textbf{(C)} and \textbf{(E)}: Ratio between the protein variance of the gene-centered and the multi-protein models. Inset: the histogram of these variances. \textbf{(B)}, \textbf{(D)} and \textbf{(F)}: Free RNA-polymerase (above) and free ribosome (below) number distribution for cells each of the volumes. In thick lines the binomial distribution predicted for the simplified model (see Section~\ref{subsec:S_mod3_simplified_models}).}
\end{figure}


\subsubsection{\label{subsec:S_mod3_Impact-free-rnap}Influence of Free RNA-polymerase Concentration}

By keeping the low concentration of free ribosomes, we have produced a series of parameters where the average concentration of free RNA-polymerases was fixed successively to $1$, $10$, $100$ and $1000$ $\text{copies}/\text{\textmu m}^{3}$.  In each case, we have deduced a set of parameters, where the affinity constants $\lambda_{1,i}$, $\lambda_{2,i}$ are still calculated so that average mRNA and protein concentrations still correspond to the experimental measures. By performing simulations, we observe that the protein variability remains in the same order of magnitude.  Even more, as shown in Figure~\ref{fig:S_mod3_sim_21}, the gap between the multi-protein model and the previous gene centered model (with volume growth, random partition and gene replication), seems to be reduced for the lowest free RNA-polymerase concentration: the variance ratio between the two models is $0.96$ ($90\%$ of the genes have a variance ratio above $0.92$).

Even for extremely low RNA-polymerase concentrations, the distribution
of free RNA-polymerases and ribosomes are still well predicted by
the simplified model presented in Section~\ref{subsec:S_mod3_simplified_models}
(see Figure~\ref{fig:S_mod3_sim_22}).

\subsubsection{\label{subsec:S_mod3_Impact-free-ribo}Influence of Free Ribosome Concentration}

With a high concentration of free RNA-polymerases, we observe the influence of the quantity of free ribosomes on the protein variance. We have computed a set of parameters based on average concentrations of non-allocated ribosomes of $1$, $10$, $100$ and $1000$ $\text{copies}/\text{\textmu m}^{3}$.  It can be first remarked that for very high free concentrations, the binomial fit of the simplified model (described in Section~\ref{subsec:S_mod3_simplified_models}) is not relevant to describe the free ribosome distribution.

In this case again, changes to the average concentration of free ribosomes are negligible. As the average concentration of free ribosomes increases, the variance of each protein decreases. As shown in Figure~\ref{fig:S_mod3_sim_31}, for a concentration of $1000\,\text{copies}/\text{\textmu m}^{3}$, the variance of the multi-protein represents on average $0.98$ of the one predicted by the gene-centered model ($90\%$ of the genes have a variance ratio above $0.93$).

Fluctuations in the number of free ribosomes seem to be the main source of the small additional variability observed in the multi-protein model (compared with the previous intermediate model); and this effect seems less important as the number of free ribosomes is high. But in real bacteria, the number of free ribosomes usually seems quite low due to the high cost of ribosome production; then, a low number of free ribosomes (like in the simulation of Section~\ref{subsec:F-RNAP-and_ribo}) seems more plausible than this simulation.

To confirm the specific influence of fluctuations of ribosomes on the protein variability, we have performed  simulations with a modified version of the model. The multi-protein model has been changed in such a way that the concentration of non-allocated ribosomes is fixed  during the whole simulation (meanwhile the free RNA-polymerases are still fluctuating). Results about protein variability are similar to what is shown in Figure~\ref{fig:S_mod3_sim_31}: the variance of each protein concentration is equivalent to what was described by the gene-centered model.

\subsection{\label{sec:S_mod3_Other-influences}Other Possible Influences on Protein Variability }

\begin{figure}

\subfloat[\label{fig:Sampling-genes}]{\includegraphics[height=5.5cm,page=1]{model3/Figure_20160910_sumary_gene_creation}

}~~\subfloat[\label{fig:Model4-Elong-1}]{
\includegraphics[height=3.5cm,page=1]{model3/Figure_20160912_FabH_example_elong}
}

\subfloat[\label{fig:Model-seques-rnap}]{\begin{tikzpicture}[scale=4]
\node[draw,shape=rectangle] (fy) at (0,0){$F_Y$};
\node[draw,shape=rectangle] (sy) at (-1,0){$D_Y$};
\node[draw,shape=rectangle] (eyi) at (1,0){$E_{Y,i}$};

\node[draw,shape=rectangle] (ey1) at (1,0.4){$E_{Y,1}$};
\node[draw,shape=rectangle] (eyk) at (1,-0.4){$E_{Y,K}$};
\node[]  at (1,0.2){...};
\node[]  at (1,-0.2){...};

\draw[transform canvas={yshift=.3ex},->]
	(fy) -- (sy) node[above,midway]
	{$\lambda_{+}F_{Y}G/(K\cdot V)$};
\draw[transform canvas={yshift=-.3ex},->]
	(sy) -- (fy) node[below,midway] {$\lambda_{-}\cdot D_Y$};
\draw[transform canvas={yshift=.3ex},->]
	(eyi) -- (fy) node[above,midway] {$\lambda_{1,i} G_i  F_Y/V$};
\draw[transform canvas={yshift=-.3ex},->]
	(fy) -- (eyi) node[below,midway] {$\mu_{1,i} E_{Y,i}$};

\draw[->,black!60]
	[bend right=27] ($(ey1.west)+(0,0.015)$) to
	 ($(fy.east)+(0,0.055)$);
\draw[->,black!60]
	($(fy.east)+(0,0.025)$) to
	[bend right=-27] ($(ey1.west)+(0,-0.015)$);

\draw[->,black!60]
	[bend right=-27]	($(eyk.west)+(0,-0.015)$) to
	($(fy.east)+(0,-0.055)$);
\draw[->,black!60]
	($(fy.east)+(0,-0.025)$) to
	[bend right=27] ($(eyk.west)+(0,0.015)$);


\draw (sy.north)+(0,0.3) node[above]{Non-specifically bound};
\draw (fy.north)+(0,0.3) node[above]{Free};
\draw (ey1.north)+(0,0.05) node[above]{Specifically bound};

\draw (ey1.north west)+(-0.05,0.05)   node (a) {};
\draw (eyk.south east)+(0.05,-0.05)   node (b) {};
\path[draw]  (a) rectangle (b);


\end{tikzpicture}}

\caption{Other features tested. \textbf{(A)}: Sampling example of mRNA and
protein concentration in order to create additional genes. \textbf{(B)}:
Model variant with production of RNA-polymerases and ribosomes. Distributions
of the number of FabH protein for cells of volume $1.40\,\text{\textmu m}^{3}$,$1.95\,\text{\textmu m}^{3}$
or $2.50\,\text{\textmu m}^{3}$. For each volume, the distribution
in the model with production of RNA-polymerases and ribosomes is similar
to the model with RNA-polymerase and ribosomes in constant concentration.
\textbf{(C)}: Model variant with non-specific RNA-polymerase binding.
An RNA-polymerase can be either among the specifically bound on the
DNA (in $E_{Y,i}$ for $i\in\{1,...,K\}$) or free (among $F_{Y}$)
or non-specifically bound on the DNA (among $D_{Y}$).}
\end{figure}

In this section, based on the set of parameters of the simulation
Section~\ref{subsec:F-RNAP-and_ribo} (with few free ribosomes and
many free RNA-polymerases), we make variations on some modeling choices
for some cellular mechanisms: a larger set of genes, RNA-polymerases
and ribosomes as a result of gene expression, the introduction of
RNA-polymerase non-specific binding on the DNA, considering uncertainty
in the division and DNA replication processes, etc. We will show that
the protein variability is quite robust to any of these changes: as
for the results presented in Section~\ref{subsec:F-RNAP-and_ribo},
the protein variance is still increased by at most $10\%$ compared
to the gene-centered model.

\subsubsection{\label{subsec:S_Other-Additional-genes}Additional Genes }
The genome of \emph{E. coli } has approximately $2000$ expressed genes.  But the measures of \cite{taniguchi_quantifying_2010} take into account only a part of it.  A total $1018$ protein types were considered in this reference, and among them, only $841$ types have their mRNA production estimated. In order to better represent the complete genome of the bacteria, we have created a set of parameters with an extended pool of additional randomly created genes so that the total number of genes would be $2000$.

For each new gene, we have sampled its average protein and mRNA concentration, an mRNA lifetime and  gene position. By studying the data of \cite{taniguchi_quantifying_2010}, we have investigated the possible statistical correlations between these quantities; it appears that only the mRNA and protein concentration are correlated (as it is shown in the Figure~\ref{fig:mrna_prot_taniguchi} of the main article). We therefore have sampled the mRNAs lifetime and the gene position and length independently from the two other quantities.

As in the dataset, the genes are evenly distributed in the DNA, the gene position is assumed to be  uniformly distributed. The empirical mRNA lifetime distribution fitted a log-normal distribution; we have chosen the mRNA lifetime accordingly.

For the mRNA and the protein expression, we have taken into account their correlation. The first step of the procedure is to sample the protein production of the new gene according protein empirical distribution of the dataset. We then obtain a realistic protein production for the new gene. In a second step, the mRNA production is chosen depending on its protein production already determined. We have subdivised the dataset in 10 classes according to the protein production (see the different colors of Figure~\ref{fig:Sampling-genes}); we then consider the classes in which the protein production of the new gene falls in. We then sample its mRNA production according to the empirical distribution of the mRNA productions of this specific class only. Thus,  the mRNA production of the new gene is sampled with a specific distribution that depends on its protein production With this procedure, the newly created genes have a protein and mRNA productions that seem to be in adequation with the original dataset (see Figure~\ref{fig:Sampling-genes}).

Simulations with the completed genome show no significant difference in terms of protein variability. In particular, the variance ratio between protein concentration of the gene-centered model and the multi-protein model is not different as in Section~\ref{subsec:F-RNAP-and_ribo}.
\subsubsection{\label{subsec:S_Other-Production-of-RNAP-ribo}Production of RNA-polymerase and Ribosomes}

In the complete stochastic model as it was presented in Section~\ref{subsec:S_mod3_Detailed_descr}, all ribosomes and all RNA-polymerases are supposed to have constant concentrations (respectively $\beta_{R}$ and $\beta_{Y}$). In reality, both RNA-polymerases and ribosomes are composed of different subunits, each subunit is either a protein or, in the case of ribosomes, a functional RNA. The variability of the production of these subunits can have an overall impact on the global production.

We have performed a preliminary simulation that takes into account this aspect: the goal is not to have a precise description of mechanisms of RNA-polymerase and ribosome productions, but rather to have an insight in the magnitude of additional variability it can induce.  In this version of the model, the expression of one gene represents RNA-polymerase production and the expression of another gene represents the ribosome production. It refers to a case where the RNA-polymerases and ribosomes would be composed of only one proteic subunit.

We therefore created two genes, whose protein production was fixed to correspond to the measured concentration of RNA-polymerases and ribosomes.  The mRNA production and lifetime, the gene position and length have been chosen by the same procedure as described in the previous subsection.

This simulation brings an additional variability in the growth rate: the volume growth is more variable. These fluctuations are directly correlated with the number of ribosomes in the cell (Figure~\ref{fig:Model4-Elong-1}, above). But surprisingly, these additional variability has no significant impact in the protein variability. The Figure~\ref{fig:Model4-Elong-1} (below) also shows the distribution of the protein FabH for cells of different volumes. This case does not differ from the case where the total amount of RNA-polymerases and ribosomes are in constant concentration.

We can propose a possible interpretation of these results. The fluctuations in the total number of ribosomes seems to influence primarily the speed of growth as shown in Figure~\ref{fig:Model4-Elong-1}.  When ribosomes are produced, it accelerates the global production of all types of proteins thus increasing the volume. As a consequence, both the production of each type of protein and the volume are co-regulated.  Fluctuations in the total number of ribosomes affect the volume growth and the production of the $i$-th protein in the same way such as in a cell of a given volume, the $i$-th protein distribution is relatively unchanged.

\subsubsection{\label{subsec:Other-Sequestered-polymerases}Non-specifically Bound
Polymerases}

In the complete stochastic model as it was described in Section~\ref{subsec:S_mod3_Detailed_descr}, RNA-polymerases are either on the DNA involved in a transcription process, or is among the $F_{Y}$ free RNA-polymerases that freely evolve in the cytoplasm. But it has been experimentally shown that large portion of RNA-polymerases can bind non-specifically on the DNA, without initiating transcription. For instance, \cite{klumpp_growth-rate-dependent_2008} estimated that around $90\%$ of the RNA-polymerases are non-specifically bound to the DNA.

We have created an alternative version of the stochastic model to introduce a third possible class for RNA-polymerases. RNA-polymerases can also bound non-specifically on the DNA. The binding rate is modeled as follows: at any time $s$, a free RNA-polymerase bind on the DNA at a rate that depends on the number of free RNA-polymerases $F_{Y}(s)$ and on the DNA concentration $G(s)/(K\cdot V(s))$; the global rate is hence $\lambda_{+}F_{Y}(s)G(s)/(K\cdot V(s))$ where $\lambda_{+}$ is a parameter that represents the natural affinity of RNA-polymerases for the DNA. Once an RNA-polymerase is bound, it is released in a time represented by an exponential random variable of rate $\lambda_{-}$ (see Figure~(\ref{fig:Model-seques-rnap})).

We performed a simulation where the parameters $\lambda_{+}$ and $\lambda_{-}$ are chosen such that around $90\%$ of the RNA-polymerases are sequestered on the DNA at any time, as it was the case in \cite{klumpp_growth-rate-dependent_2008}.  The variability of protein  concentration does not seem to be impacted in this case.

\subsubsection{\label{subsec:Uncertainty-in-divi-init}Uncertainty in the
Replication Initiation and Division Timing}

In the complete stochastic model as it was initially described, replication initiation and division occur when the cell reaches the respective volumes of $V_{I}$ and $2V_{0}$. In practice it does not occur in this way. We propose here a modification of the stochastic model to take into account this aspect.

The way replication and division occur is still a disputed topic, see for example \cite{tyson_sloppy_1986,wang_robust_2010,soifer_single-cell_2014,osella_concerted_2014}.  For the division, one hypothesis (referred as ``sizer model'') is that the division decision depends on the current size of the cell (the size can refer to the mass or the volume, but as explained in Section~\ref{subsec:S_mod3_Detailed_descr}, the density constraint, see \cite{marr_growth_1991}, ensures a proportionality relation between these two quantities). With this hypothesis, at any instant, the instantaneous probability to divide depends only on the current cell size. In a first approximation, the cell size distributions observed experimentally can be explained by this ``sizer model'', see \cite{robert_division_2014,osella_concerted_2014}.  It is therefore this framework that we have considered to represent the cell division decision.

At time $s$ of the simulation, with a cell of volume $V(s)$, we introduce an instantaneous division rate $b_{D}(V(s))$ with $b_{D}$ is a positive function (the probability to divide between times $s$ and $s+\mathop{}\mathopen{}\mathrm{d} s$ is given by $b_{D}(V(s))\,\mathop{}\mathopen{}\mathrm{d} s$ ). The division decision is hence only volume dependent. The function $b_{D}$ is chosen so that the division occurs around the volume $2V_{0}$ with $V_{0}=1.3\,\text{\textmu m}^{3}$ and division precision can be fixed (for more information about the function $b_{D}$, see Chapter~4 of \cite{dessalles_stochastic_2017}).

Similarly for the replication initiation decision, the stochastic model initially described considers a fixed volume $V_{I}$ at which the DNA replication is initiated. We introduce variability in this cell decision, in the same way as we do for the division: at time $s$, we consider a replication initiation rate $b_{I}(V(s))$ such as the function $b_{I}$ is chosen in order to have a replication initiation that occurs around volume $V_{I}$.

We did several simulations with different functions $b_{D}$ and $b_{I}$ in order to have different precisions in the division and replication initiation decisions. The protein variability does seem to be changed significantly by any of these scenarios.


\subsubsection{\label{subsec:Other-deterministic-times}Deterministic Time for Replication }

When the stochastic model has been initially presented (in Section~\ref{subsec:S_mod3_Detailed_descr}), we have proposed two ways to model the time of DNA replication initiation $\tau_{I}$ . It can either be a deterministic time after the last division, or it can happen when the cell reaches the specific volume $V_{I}$. We have checked that this modeling choice has no significant influence on the global dynamic of the system, in particular in the protein noise.

\subsection{\label{subsec:S_mod3_ESD}Environmental State Decomposition}

In order to perform the environmental state decomposition, one has to specify the environment $Z$ with which the conditioning is made in Equation~\ref{eq:ESD}. In order to have a decomposition analog to the result of the dual reporter technique, one has to take into account the environmental aspects that would be similar in the expression of two identical genes in the same cell. In our model, two gene expressions in the same cell would undergo the same volume growth (as it was the case for the first two intermediate models), the same number of free RNA-polymerases, and the same number of free ribosomes. But on the contrary, the division would be specific to each expression. For a protein of type $i$, the decomposition would therefore be:
\[
\V{P_{i}/V}=\underbrace{\vphantom{\left[\Ehat{P_{i}|}\right]}\Ehat{\Vhat{P_{i}(t)/V(t)|\left(V(t),F_{Y}(t),F_{R}(t)\right)}}}_{\Vhatt{P_{i}/V}{int}}+\underbrace{\Vhat{\Ehat{P_{i}(t)/V(t)|\left(V(t),F_{Y}(t),F_{R}(t)\right)}}}_{\Vhatt{P_{i}/V}{ext}}.
\]
Estimations of $\Vhatt{P_{i}/V}{int}$ and $\Vhatt{P_{i}/V}{ext}$ are directly made on the results of the simulations. The result is that the extrinsic contribution of the variance $\Vhatt{P_{i}/V}{ext}$ represent only a very small portion of $\Vhat{P_{i}/V}$ for any gene $i$.

\subsection{\label{subsec:S_mod3_simplified_models}Simplified Models for Transcription
and Translation}

\cite{fromion_stochastic_2015} propose a multi-protein model for translation, with a  limited number of ribosomes, where each type of mRNA is supposed to be in constant quantities and where the maximum number of ribosomes on one single mRNA is limited.  The system also evolves in a fixed volume as it is the case for classical models.

In order to have a prediction for the number of respectively free
RNA-polymerases and free ribosomes; we consider two analog models
that are slightly simplified versions of the model of \cite{fromion_stochastic_2015}.
They respectively represent the transcription and translation part
and they are completely independent.

The goal is, for each of these models, to provide the equivalent of
the first results of \cite{fromion_stochastic_2015} and we will show
that the expected distribution of free RNA-polymerases (or free ribosomes)
is binomial in these simplified cases.

\subsubsection*{A Model for Transcription}

As explained in Section~\ref{subsec:F-RNAP-and_ribo}, one can interpret the model of \cite{fromion_stochastic_2015} as taking place in a fixed volume $V$ (it would correspond to a small portion of the cell cycle in the complete stochastic model, a portion where the volume of the cell does not change much). We also consider that the gene copy of each unit of production remains constant; as a consequence, the gene copy number of the $i$-th gene $G_{i}$ is constant and known. As in the complete stochastic model, and contrary to the model of \cite{fromion_stochastic_2015}, we consider that there is no limiting number of elongating RNA-polymerases on one gene.

In a pool of $K$ genes, denote by $N_{Y}$ the constant total number of polymerases. We consider the random variables $E_{Y,i}$ for $i\in\{1,...,K\}$ be the number of RNA-polymerases attached to the $i$-th gene. As a consequence, the random variable
\begin{equation}
F_{Y}:=N_{Y}-\sum_{i=1}^{K}E_{Y,i}\label{eq:rnap-decomp}
\end{equation}
is the number of free RNA-polymerases in the system.

The process $X(t)=\left(E_{Y,i}(t),i\in\{1,...K\}\right)$ takes place
in the state place $S$ the subset of $\mathbb N^{K}$ such as
\[
S:=\left\{ x\in\mathbb{N}^{K},\sum_{i=1}^{K}x_{i}\leq N_{Y}\right\} .
\]
There are at most $N_{Y}$ RNA-polymerases that can be attached to genes at the same time. We can describe the Markov process transition by the following \emph{Q}-matrix: by setting the vector $e_{i}=(\delta_{i'=i})_{i'\in\{1,...,K\}}$ ($\delta$ is used here as the Kronecker delta), for any $x,y\in S$,
\[
\begin{cases}
q(x,x+e_{i})=\lambda_{1,i}G_{i}\lambda_{1,i}f(x)/V & \text{for any }i\in\{1,...K\},\\
q(x,x-e_{i})=\mu_{1,i}x_{i} & \text{for any }i\in\{1,...K\}\text{, if \ensuremath{x_{i}>0}},\\
q(x,y)=0 & \text{if }\left\Vert x-y\right\Vert >1.
\end{cases}
\]
where
\[
f_{Y}(x):=N_{Y}-\sum_{i=1}^{K}x_{i}
\]
the number of free RNA-polymerases. Equation~\ref{eq:rnap-decomp}
leads in particular to $f(x-e_{i})=f(x)\boldsymbol{+}1$ for all $i$.

As in \cite{fromion_stochastic_2015}, we are looking for an invariant reversible probability measure $\pi$ of this Markov process, i.e.  for any $i\in\{1,\dots,K\}$,
\[
\pi(x)\mu_{1,i}x_{i}  = \pi(x-e_{i})\cdot\lambda_{1,i}G_{i}\left(f_{Y}(x)+1\right)/V.\label{eq:equilibre-point-fixe}
\]
\begin{prop}
The invariant measure $\pi$ of the number of RNA-polymerases in each gene has the following form
\[
\pi(x)=\frac{1}{Z}\cdot\frac{1}{f_{Y}(x)!}\prod_{i=1}^{K}\frac{\left(G_{i}\lambda_{1,i}/(V\mu_{1,i})\right)^{x_{i}}}{x_{i}!}
\]
for any $x\in S$ and with $Z>0$ the normalization constant.
\end{prop}
\begin{proof}
We only need to check that $\pi$ satisfies Equation~\ref{eq:equilibre-point-fixe}. For a gene $i\in\{1,...,K\}$, we take $x\in S$ such that $x-e_{i}\in S$ .
\begin{align*}
\pi(x)\mu_{1,i}x_{i} & =  \frac{1}{Z}\cdot\frac{1}{f_{Y}(x)!}\prod_{i'=1}^{K}\frac{\left(G_{i'}\lambda_{1,i'}/(V\mu_{1,i'})\right)^{x_{i'}}}{x_{i'}!}\mu_{1,i}x_{i}\\
 & =  \frac{1}{Z}\cdot\frac{1}{(f_{Y}(x)+1)!}\prod_{i'\ne i}^{K}\left(\frac{\left(G_{i'}\lambda_{1,i'}/(V\mu_{1,i'})\right)^{x_{i'}}}{x_{i'}!}\right)\cdot\frac{\left(G_{i}\lambda_{1,i}/(V\mu_{1,i})\right)^{(x_{i}-1)}}{(x_{i}-1)!}\cdot\frac{G_{i}\lambda_{1,i}}{V}(f_{Y}(x)+1)\\
 & =  \pi(x-e_{i})\cdot G_{i}\lambda_{1,i}(f_{Y}(x)+1)/V.
\end{align*}
So $\pi$ satisfies Equation~\ref{eq:equilibre-point-fixe}.
\end{proof}
We can now derive from the previous proposition the equilibrium distribution of $F_{Y}$, the number of free RNA-polymerases of the process.
\begin{prop}
\label{prop:app-FY}The distribution of the  number of free polymerases $F_{Y}$ is given by
\[
\mathbb P\left[F_{Y}=n\right]=\binom{N_{Y}}{n}\frac{\Lambda_{Y}^{N_{Y}-n}}{(1+\Lambda_{Y})^{N_{Y}}}\mbox{,}
\]
with $\Lambda$ defined such as
\[
\Lambda_{Y}:=\sum_{i=1}^{K}G_{i}\lambda_{1,i}/(V\mu_{1,i}).
\]
\end{prop}
$F_{Y}$ follows a binomial distribution ${\cal B}\left(\phi,N\right)$ with parameters $\phi=(1+\Lambda_{Y})^{-1}$ and $N=N_{Y}$.
\begin{proof}
From  Equation~\ref{eq:equilibre-point-fixe}, it follows that for  $n\in\{1,...,N_{Y}\}$,
\begin{align*}
\mathbb P\left[F_{Y}=n\right]
 & =  \sum_{x\in S}\frac{1}{Z}\cdot\frac{1}{n!}\prod_{i=1}^{K}\frac{\left(G_{i}\lambda_{1,i}/(V\mu_{1,i})\right)^{x_{i}}}{x_{i}!}\ind{\sum_{i}x_{i}=N_{Y}-n}\\
 & =  \frac{1}{Z}\cdot\frac{1}{n!}\sum_{x\in S}\prod_{i=1}^{K}\frac{\left(G_{i}\lambda_{1,i}/(V\mu_{1,i})\right)^{x_{i}}}{x_{i}!}\ind{\sum_{i}x_{i}=N_{Y}-n}
 =  \frac{1}{Z_{F_{Y}}}\cdot\frac{1}{n!}\mathbb P\left[\sum_{i=1}^{K}\sum_{k=1}^{G_{i}}C_{i,k}=N_{Y}-n\right]
\end{align*}
with $\forall i\in\{1,\dots,K\}$ and $\forall k\in\{1,\dots,G_{i}\}$,
$C_{i,k}\sim\cal P(\lambda_{1,i}/(V\mu_{1,i}))$.

Since the independent random variables $C_{p,k}$ are Poisson,  their sum is also Poisson  with  parameter $\Lambda:=\sum_{i=1}^{K}G_{i}\lambda_{1,i}/(V\mu_{1,i})$.
By summing up the previous relation, one gets that

\begin{align*}
1 & =  \frac{1}{Z_{F_{Y}}}\cdot\sum_{n=0}^{N_{Y}}\frac{1}{n!}\mathbb P\left[C_{1,1}=N_{Y}-n\right] =  \frac{1}{Z_{F_{Y}}}\cdot\sum_{n=0}^{N_{Y}}\frac{1}{n!}e^{-\Lambda}\frac{\Lambda^{N_{Y}-n}}{(N_{Y}-n)!} =  \frac{1}{Z_{F_{Y}}}\cdot\frac{1}{N_{Y}!}e^{-\Lambda}\cdot\sum_{n=0}^{N_{Y}}\binom{N_{Y}}{n}\Lambda^{N_{Y}-n},
\end{align*}
hence $Z_{F_{Y}}=e^{-\Lambda}(1+\Lambda)^{N_{Y}}/N_{Y}!$.
\end{proof}

\subsubsection*{A Model for Translation}
The model for translation considered here is analogous to the transcription case. We still consider that the volume $V$ is fixed and that for each gene, the number $M_{i}$ of mRNAs of type $i$ is known and constant (because of these, the process describes here is independent of transcription).  As in the complete stochastic model, and contrary to the model of \cite{fromion_stochastic_2015}, we consider that there is no limiting number of elongating ribosomes on one mRNA.

Similarly to the transcription, we can define $N_{R}$ (the total number of ribosomes), $E_{R,i}$ (the number of ribosomes elongating an mRNA of type $i$) and $F_{R}$ (the number of free ribosomes) such as
\[
F_{R}:=N_{R}-\sum_{i=1}^{K}E_{R,i}.
\]
The rate at which a ribosome is sequestered on a type $i$ mRNA is
therefore $M_{i}\lambda_{2,i}/V$, and the rate at which an elongation
terminates on a type $i$ mRNA is $\mu_{2,i}E_{R,i}$.

As this model is analog to the transcription case, we can also prove
that
\begin{prop}
\label{prop:app-FR}The number of free ribosomes $F_{R}$ follows
\[
\mathbb P\left[F_{R}=n\right]=\binom{N_{R}}{n}\frac{\Lambda_{R}^{N_{R}-n}}{(1+\Lambda_{R})^{N_{R}}}\mbox{,}
\]
with $\Lambda_{R}$ defined such as
\[
\Lambda_{R}:=\sum_{i=1}^{K}M_{i}\lambda_{2,i}/(V\mu_{2,i}),
\]
 $F_{R}$ follows a binomial distribution ${\cal B}\left(\phi,N\right)$
for which $\phi=(1+\Lambda_{R})^{-1}$ and $N=N_{R}$.
\end{prop}
