### Cell culture

All experiments were performed on 6C2 cells, a chicken erythroblast cell line transformed by the avian erythroblastosis virus (*AEV*) [79, 80]. Cells were maintained in alpha minimal essential medium (Gibco-BRL,Gaithersburg, MD, USA) supplemented with 10% (v/v) fetal bovine serum, 1% (v/v) normal chicken serum, 100 µmol/l β-mercaptoethanol (Sigma-Aldrich, St Louis, MO, USA), 100 units/ml penicillin and 100 μg/ml streptomycin (Gibco-BRL), at a maximum density of 1 × 10^{6} cells per ml.

### Generation of stably transfected clones

Stably transfected clones, expressing a fluorescent reporter, were obtained as previously described [81]. Briefly, 6C2 cells were nucleofected in a transfection apparatus (Nucleofector™ II; Amaxa Nucleofector™ Technology) (T-16 program) using a commercial kit (Cell Line Nucleofector® Kit V; Lonza GmBH, Cologne, Germany) and a pT2.CMV-*mCherry*/pCAGGS-T2TP plasmid mix (ratio 5/1). The pT2.CMV-*mCherry* plasmid was constructed using the same strategy as described for the pT2.CMV-*hKO* plasmid [81], except that the *hKO* reporter gene was replaced by *mCherry*, extracted from the pRSET-B plasmid (kindly provided by Dr Roger Tsien, University of California, San Diego, CA, USA). mRNA birth/death fluctuations constitute a major source of stochasticity in gene expression because many mRNA species are present at very low molecular counts within cells [58, 70, 82, 83], thus we reduced this source of intrinsic noise by using the cytomegalovirus (CMV) promoter. Obtaining a strong signal also allowed us to overcome bias caused by autofluorescence in the flow-cytometry data. The integration into genomic DNA of the reporter is allowed by the Tol2 transposon system [84]; the CMV-*mCherry* sequence, flanked by Tol2 motifs, is recognized by a transposase (pCAGGS-T2TP), and randomly inserted into 6C2 genomic DNA. Seven days after transfection, stably transfected cells expressing the reporter gene were sorted and individually cloned in U-shaped 96-well microplates (Cellstar Greiner Bio-One GmBH, Frickenhausen, Germany) using a cytometer (FACSVantage SE; Becton-Dickinson, Franklin Lakes, NJ, USA).

### Molecular and cellular characterization of clones

For each clone, the genomic reporter insertion sites were identified using a splinkerette PCR method as previously described [81], in order to select only clones with a single insertion site. Briefly, genomic DNA isolated from clones expressing the gene reporter was purified by phenol extraction and ethanol precipitation, before being digested for 16 hours at 65°C with *Tai*I, a restriction enzyme with a 4 bp recognition site. The digested DNA was then ligated to a splinkerette adaptor for 1 hour at 22°C. After purification of the ligated product, two rounds of PCR (PCR1 and nested PCR2) were performed using primers specific for the reporter transgene *mCherry* and for the annealed splinkerette adaptor, and a commercial polymerase (AccuPrime™ Taq DNA Polymerase High Fidelity; Invitrogen Inc., Carlsbad, CA, USA). The PCR products were then purified and sequenced. Finally, the genomic reporter insertion sites were identified by similarity searches using the sequence analysis tool iMapper [85]. The identification of the insertion sites of the selected clones was confirmed using a high-throughput splinkerette-PCR method [86], allowing the analyses of hundreds of clones. This work will be described in details elsewhere.

For characterization of clones and analysis of treatment effects (see below), flow-cytometry analyses were performed (FACSCanto II; Becton-Dickinson) on cells extemporaneously pelleted and resuspended in Dulbecco's phosphate-buffered saline 1× solution (Gibco-BRL). Each sample was analyzed using an acquisition of 50,000 events (gated on living cells), and the positive fluorescence threshold was fixed using non-transfected cells. Possible variability resulting from flow-cytometer calibration was taken into account by systematically analyzing flow-calibration particles (SPHERO™ Rainbow; Spherotech Inc., Lake Forest, IL, USA), as a calibration reference.

Non-transfected cells were used to measure 6C2 native autofluorescence, and the difference between the fluorescence of transfected and non-transfected ones was used as an indicator of the transgene activity (note that autofluorescence was also systematically added to the model's output to compute the distribution distance scores).

For each clone, two indicators were systematically used: MFI (mean fluorescence intensity) and NV (the variance divided by the square mean).

For a given cell, the measured fluorescence *f* (from the flow cytometer) is *f* = *f*_{
t
} + *f*_{
a
}; that is, the sum of the true fluorescence *f*_{
t
} (coming from the reporter proteins) and the autofluorescence *f*_{
a
} (coming from the rest of the cell). The autofluorescence is not a constant, but has a distribution that is obtained using non-transfected cells. The two first moments of *f* read simply as

\u27e8f\u27e9=\u27e8{f}_{t}\u27e9+\u27e8{f}_{a}\u27e9

and

{\sigma}^{2}\left(f\right)={\sigma}^{2}\left({f}_{t}\right)+{\sigma}^{2}\left({f}_{a}\right).

Hence, with MFI and NV being the mean and normalized variance of the true fluorescence, we get:

MFI=\u27e8f\u27e9-\u27e8{f}_{a}\u27e9

and

NV=\frac{{\sigma}^{2}\left(f\right)-{\sigma}^{2}\left({f}_{a}\right)}{{\left(\u27e8f\u27e9-\u27e8{f}_{a}\u27e9\right)}^{2}}

Finally, to compare the theoretical distributions obtained from simulations (which only included the reporter fluorescence) with those obtained from experiments (which also included the autofluorescence), the model's output was first combined with the experimental autofluorescence. This was carried out by summing each simulation result with the value of a randomly selected cell from the autofluorescence distribution. The resulting distribution was the convolution between the theoretical and the autofluorescence distributions, and was then compared with the experimental distributions using a Kolmogorov-Smirnov test.

### Determination of *mCherry*mRNA and protein degradation rates

To determine the *mCherry* mRNA degradation rate, the mRNA concentration was estimated using quantitative reverse transcription (qRT)-PCR after transcription inactivation was achieved using actinomycin D treatment. Two clones (C5 and C11) were treated, in duplicate, for 0, 60, 124, 244 and 488 minutes with a final concentration of 10 µg/ml actinomycin D (A9415; Sigma-Aldrich), before extracting the mRNA after the instructions of RNeasy® Plus Mini Kit (Qiagen Inc., Valencia, CA, USA). To prepare the real-time PCR assay, 1 µg of total RNA from each sample was reversed transcribed using the SuperScript^{™} III First-Strand Synthesis System for RT-PCR (Invitrogen Inc.) in the presence of random hexamers. Quantification of mRNA levels by real-time PCR was performed in 96-well plates using a real-time PCR system (LightCycler 480; Roche Diagnostics, Basel, Switzerland). The measurement was performed in a final volume of 10 µl of reaction mixture (containing 2.5 µl of cDNA template diluted 1 in 5), prepared using a commercial kit (Light Cycler 480 SYBR Green I Kit; Roche Diagnostics) in accordance with the manufacturer's instructions, and with the primer set at a final concentration of 0.5 µmol/l (mCher-For: CCACCTACAAGGCCAAGAA, mCher-Rev: ACTTGTACAGCTCGTCCATG). An internal standard curve was generated using serial dilutions (from 2000 to 0.02 fg/µl) of purified PCR product. The reactions were initiated by activation of *Taq* DNA polymerase at 95°C for 5 minutes, followed by 45 three-step amplification cycles consisting of denaturation at 95°C for 15 seconds, annealing at 55°C for 15 seconds, and extension at 72°C for 15 seconds. The fluorescence signal was measured at the end of each extension step. After the amplification, a dissociation stage was run to generate a melting curve for verification of amplification-product specificity. The crossing point (CP) was determined by the second derivative maximum method in the LightCycler® 480 software (version 1.5.0). After normalization, taking into account cellular viability and mRNA quantity used for the retrotranscription step, the mRNA half-life was determined by fitting mRNA quantity evolution by a decreased exponential (least square) method.

To determine the mCherry protein degradation rate, we used flow cytometry to measure the protein half-life after translation inactivation using cycloheximide treatment. C5 and C11 clones were treated in duplicate for 0, 16, and 24 hours with a final concentration of 100 µg/µl cycloheximide (C4859; Sigma-Aldrich), and for each time point, the fluorescence of the treated cells was measured by flow cytometry. The autofluorescence component was removed as explained earlier. The protein half-life was determined using exponential fit of the fluorescence mean decrease curve, similarly to the procedure used for determining the mRNA half-life.

### Treatments with chromatin-modifying agents

To analyze the effect of chromatin state on the stochasticity of gene expression, clones were treated with TSA, a histone deacetylase inhibitor (P5026; Sigma-Aldrich) and 5-AzaC, an inhibitor of DNA methylation (A2385; Sigma-Aldrich). For each clone, kinetic treatment experiments were performed; clones were treated with 500 nmol/l TSA or 500 µmol/l 5-AzaC at five time points (0, 8, 24, 32, and 48 hours). For each time point, 1 × 10^{6} cells (for 0, 8, and 24 hours) or 5 × 10^{5} cells (for 32 and 48 hours) were treated with the relevant drug and characterized by flow cytometry.

### Model description

The two-state model of gene expression represents the chromatin activity as an 'on-off' process specified through the transition rates *k*_{
on
} and *k*_{
off
} (respectively representing the 'off-on' transition and the 'on-off' transition). To enable comparison with the experimental data, a simple model of mRNA and protein dynamics based on two production/degradation models completed the model. The production of mRNA was allowed only in the 'on' state (open chromatin) but completely forbidden in the 'off' state (closed chromatin). The model thus corresponds to the following equations:

\left\{\begin{array}{l}{k}_{T}=\frac{{k}_{on}}{{k}_{on}+{k}_{off}}\\ \frac{dR}{dt}=\rho {k}_{T}-\stackrel{\u0303}{\rho}R\\ \frac{dP}{dt}=\gamma R-\stackrel{~}{\gamma}P\\ f=\alpha P\end{array}\right.

(1)

where*,k*_{
on
} is the closed-to-open transition rate, *k*_{
off
} is the open-to-closed transition rate, and *k*_{
T
} is the resulting proportion of the 'on' state; *R* is the number of mRNAs, *ρ* is the mRNA production rate (when chromatin is open), and \stackrel{\u0303}{\rho} is the mRNA degradation rate; *P* is the number of mCherry proteins, *γ* is the mCherry production rate (per mRNA) and \stackrel{\u0303}{\gamma} is the mCherry degradation rate; *f* is the fluorescence intensity of the cell (after subtraction of the autofluorescence) and *α*, is a linear proportionality coefficient to convert the number of proteins into arbitrary fluorescence measures.

This model can be simulated with the SSA (see below) to ascertain the behavior of single cells and eventually to compute the fluorescence distributions. It can also be analytically derived to compute the MFI and NV of large cell populations at steady state.

### Analytical derivation of the model

Paulsson proposed an analytic expression of the mean quantity and NV of protein in the two-state model, as a function of chromatin-dynamics parameters and transcription-translation parameters [45]. In the case of a single gene and taking into account the parameter *α*, Paulsson's equation gives:

\left\{\begin{array}{l}MFI=\alpha \phantom{\rule{0.3em}{0ex}}\frac{\rho \phantom{\rule{0.3em}{0ex}}\gamma}{\stackrel{~}{\rho}\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\gamma}}\frac{{k}_{on}}{{k}_{on}+{k}_{off}}\\ NV=\left(\frac{\stackrel{~}{\rho}\phantom{\rule{0.3em}{0ex}}\stackrel{~}{\gamma}}{\rho \phantom{\rule{0.3em}{0ex}}\gamma}\phantom{\rule{0.3em}{0ex}}\frac{{k}_{on}+{k}_{off}}{{k}_{on}}\right)+\left(\frac{\stackrel{~}{\rho}}{\rho}\phantom{\rule{0.3em}{0ex}}\frac{{k}_{on}+{k}_{off}}{{k}_{on}}\phantom{\rule{0.3em}{0ex}}\frac{\stackrel{~}{\gamma}}{\stackrel{~}{\rho}+\stackrel{~}{\gamma}}\right)+\left(\frac{{k}_{off}}{{k}_{on}}\phantom{\rule{0.3em}{0ex}}\frac{\stackrel{~}{\gamma}}{\stackrel{~}{\rho}+\stackrel{~}{\gamma}}\phantom{\rule{0.3em}{0ex}}\frac{1+\frac{\stackrel{~}{\rho}}{\stackrel{~}{\gamma}+{k}_{on}+{k}_{off}}}{1+\frac{{k}_{on}+{k}_{off}}{\stackrel{~}{\rho}}}\right)\end{array}\right.

(2)

This equation can be used to express *k*_{
on
} and *k*_{
off
} as a function of MFI, NV, and the transcription-translation parameter sets. Rewriting the equation gives:

\left\{\begin{array}{c}A=\alpha \frac{\rho \gamma}{MFI\phantom{\rule{0.5em}{0ex}}\stackrel{~}{\rho}\stackrel{~}{\gamma}}\hfill \\ B=\frac{A\frac{\stackrel{~}{\rho}}{\rho}\left(\frac{\stackrel{~}{\gamma}}{\gamma}+\frac{\stackrel{~}{\gamma}}{\stackrel{~}{\rho}+\stackrel{~}{\gamma}}\right)-NV}{\left(A-1\right)\phantom{\rule{0.3em}{0ex}}\frac{\stackrel{~}{\gamma}}{\stackrel{~}{\rho}+\stackrel{~}{\gamma}}}\hfill \\ {\left(A\left(B\left(\stackrel{~}{\gamma}+\stackrel{~}{\rho}\right)+\stackrel{~}{\rho}\right)\right)}^{2}-4{A}^{2}B\stackrel{~}{\rho}\left(\stackrel{~}{\rho}+\stackrel{~}{\gamma}\left(B+1\right)\right)\phantom{\rule{0.3em}{0ex}}>0\hfill \\ {k}_{on}=\frac{-A\left(B\left(\stackrel{~}{\gamma}+\stackrel{~}{\rho}\right)+\stackrel{~}{\rho}\right)-\sqrt{{\left(A\left(B\left(\stackrel{~}{\gamma}+\stackrel{~}{\rho}\right)+\stackrel{~}{\rho}\right)\right)}^{2}-4{A}^{2}B\stackrel{~}{\rho}\left(\stackrel{~}{\rho}+\stackrel{~}{\gamma}\left(B+1\right)\right)\phantom{\rule{0.3em}{0ex}}}}{2{A}^{2}B}\hfill \\ {k}_{on}>0\hfill \\ {k}_{off}=\left(A-1\right){k}_{on}\hfill \end{array}\right.

(3)

### Parametric exploration of the analytical model

Because the clonal populations differed only in their insertion points (that is, their chromatin-dynamics parameters), equation 3 enabled us to find the clone-specific parameters from MFI and NV (measured by flow cytometry) and the transcription-translation parameters *ρ*, \stackrel{\u0303}{\rho}, *γ*, \stackrel{\u0303}{\gamma} and *α*.\stackrel{\u0303}{\rho} and \stackrel{\u0303}{\gamma} can be determined experimentally (see above) but *ρ*, *γ* and *α* remained unknown. We explored a wide range of these parameters (large enough to include all biologically relevant values): *ρ* = 6.0, 1.0, 0.5, 0.333, 0.25, 0.200, 0.1666, 0.14286, 0.125, 0.111, 0.100, 0.0500, 0.0333, 0.0250, 0.02, 0.01666, 0.01333, 0.0111, 0.00952, and 0.00833 mRNA/min, corresponding to one mRNA produced each 1/*ρ* = 10 seconds, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, and 50 minutes, and 1, 1.25, 1.5, 1.75, and 2 hours, when chromatin is in the open state; *γ* = 1.0, 0.200, 0.100, 0.0333, 0.01667, 0.006667, 0.00333, 0.00222, 0.001666, 0.001333, 0.00111, 0.0008333, 0.00069444, 0.00046296, and 0.0003472 protein/min/mRNA, corresponding to a protein produced each 1/*γ* = 1, 5, 10, and 30 minutes, 1, 2.5, 5, 7.5, 10, 12.5, 15, and 20 hours, and 1, 1.5 and 2 days per mRNA molecule; and *α* = 0.10, 0.15, 0.50, 1.0, 1.5, 5.0, 10.0, 15.0, 50.0, 100.0, 150.0, and 200.0 arbitrary units.

Exploring all values of *ρ*, *γ* and *α* gave us 3,600 couples (*k*_{
on
}; *k*_{
off
}), of which only 1,047 respected the condition mentioned in equation 3 (*k*_{
on
}>0).

### Comparison between the analytical model and the trichostatin A-treated clones

Equation 1 enabled us to compute the mean mRNA number (*R*) and the mean protein number (*P*) at steady state, from the values of the chromatin-dynamics and transcription-translation parameters:

\left\{\begin{array}{c}R=\frac{\rho}{\stackrel{~}{\rho}}\left(\frac{{k}_{on}}{{k}_{on}+{k}_{off}}\right)\\ P=\left(\frac{\rho}{\stackrel{~}{\rho}}\phantom{\rule{0.3em}{0ex}}\frac{\gamma}{\stackrel{~}{\gamma}}\right)\left(\frac{{k}_{on}}{{k}_{on}+{k}_{off}}\right)\end{array}\right.

(4)

Then, assuming that at *t* = 0, the cell switches to a new chromatin dynamics (because of the TSA treatment), compute *R*(*t*), (the evolution of mRNA number), and *P*(*t*), (the evolution of protein number following the TSA treatment) can be computed. If {k}_{on}^{TSA} and {k}_{off}^{TSA} are the new chromatin-dynamics parameters induced by the treatment, the equation is:

\left\{\begin{array}{c}R\left(t\right)=\frac{\rho}{\stackrel{~}{\rho}}\left(\frac{{k}_{on}}{{k}_{on}+{k}_{off}}-\frac{{k}_{on}^{TSA}}{{k}_{on}^{TSA}+{k}_{off}^{TSA}}\right){e}^{-\stackrel{~}{\rho}t}+\frac{{k}_{on}^{TSA}}{{k}_{on}^{TSA}+{k}_{off}^{TSA}}\frac{\rho}{\stackrel{~}{\rho}}\\ P\left(t\right)=\left(\frac{\rho}{\stackrel{~}{\rho}}\frac{\gamma}{\stackrel{~}{\gamma}}-\frac{\rho}{\stackrel{~}{\rho}}\phantom{\rule{0.3em}{0ex}}\frac{\gamma}{\left(\stackrel{~}{\gamma}-\stackrel{~}{\rho}\right)}\right)\left(\frac{{k}_{on}}{{k}_{on}+{k}_{off}}-\frac{{k}_{on}^{TSA}}{{k}_{on}^{TSA}+{k}_{off}^{TSA}}\right){e}^{-\stackrel{~}{\gamma}t}\\ \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}+\frac{\rho}{\stackrel{~}{\rho}}\frac{\gamma}{\left(\stackrel{~}{\gamma}-\stackrel{~}{\rho}\right)}\left(\frac{{k}_{on}}{{k}_{on}+{k}_{off}}-\frac{{k}_{on}^{TSA}}{{k}_{on}^{TSA}+{k}_{off}^{TSA}}\right){e}^{-\stackrel{~}{\rho}t}+\frac{\rho}{\stackrel{~}{\rho}}\frac{\gamma}{\stackrel{~}{\gamma}}\phantom{\rule{0.3em}{0ex}}\frac{{k}_{on}^{TSA}}{{k}_{on}^{TSA}+{k}_{off}^{TSA}}\hfill \end{array}\right.

(5)

The exact values of {k}_{on}^{TSA} and {k}_{off}^{TSA} remained unknown at this stage, but we could simulate the extreme situation by assuming that, under TSA treatment, the chromatin is fully open. Analytically, this gives:

\frac{{k}_{on}^{TSA}}{{k}_{on}^{TSA}+{k}_{off}^{TSA}}=1

(6)

Note that this equation represents an extreme situation, not the exact TSA influence on chromatin.

Introducing equation 6 into the dynamics of equation 5, we were able to compute, for a given transcription-translation parameter set, the maximum rate of protein concentration increase, and thus the maximum increase of reporter fluorescence. For each parameter set, we compared the predicted fluorescence increase under the extreme condition of a fully open chromatin. We then rejected all parameter sets for which the protein number did not increase sufficiently rapidly to account for the fluorescence increase measured experimentally during TSA treatment.

### Simulation of the model

The model can be simulated using an SSA, which is an exact continuous-time algorithm that enables simulation of chemical-reaction systems [54]. Each simulation represents one of the possible realizations of the system from a specified initial state and for a given kinetic parameter set (these parameters being here considered as probabilities). Each realization depends on a pseudo-random generator, and different realizations (that is, simulations of different cells issued from the same clone) can be computed by simply initializing this random generator with different seeds. The implementation of the two-state model (equation 1) in the SSA enables simulation of the entire system dynamics and visualization of the course of chromatin state, gene transcription, number of mRNAs, mRNA translation, number of proteins and, ultimately fluorescence, in a virtual single cell. By simulating a large number of such 'artificial cells', we were able to simulate 'virtual flow-cytometry experiments' and to compute MFI, NV, and full distribution for a given parameter set. We simulated 50,000 virtual cells for 30,000 minutes (a sufficiently long period to ensure that all cells were at a steady state, the concentration values being initialized at the theoretical values given by the analytical model). The fluorescence of each cell was then computed, and the simulated distribution generated through convolution with the autofluorescence of the 6C2 cells measured experimentally (see above). Simulated distributions were then compared with the experimental distribution using the Kolmogorov-Smirnov test. The quality of each parameter set was then evaluated (the score of a given parameter set being the mean Kolmogorov-Smirnov score of each clone). The best parameter set was thus the one that gave the best fit for all six clonal populations.

### Simulation of trichostatin A treatment in the model

Using the best parameter set, we simulated 50,000 cells of the two TSA-treated clones C5 and C11 for 30,000 minutes. The chromatin-dynamics parameters were then modified to account for the TSA treatment, and the two clones were simulated for a further 1,152 minutes (48 hours). For each clones, the simulated distributions were computed after 8, 24, 32 and 48 hours, and compared with the experimental distributions using a Kolmogorov-Smirnov test. The best chromatin-dynamics parameters ({k}_{on}^{TSA}; {k}_{off}^{TSA}) were those that gave the best mean score at the four time points. In total, 11 different chromatin-dynamics values were tested for each clone. Note that, knowing the MFI value of the treated clones, we could analytically compute the {k}_{T}^{TSA} value (using equation 5). Taking

{k}_{T}^{TSA}=\frac{{k}_{on}^{TSA}}{{k}_{on}^{TSA}+{k}_{off}^{TSA}}

(from equation 1), we can use this analytical value to simplify the parametric exploration.