Revealing cell cycle control by combining model-based detection of periodic expression with novel cis-regulatory descriptors

Background We address the issue of explaining the presence or absence of phase-specific transcription in budding yeast cultures under different conditions. To this end we use a model-based detector of gene expression periodicity to divide genes into classes depending on their behavior in experiments using different synchronization methods. While computational inference of gene regulatory circuits typically relies on expression similarity (clustering) in order to find classes of potentially co-regulated genes, this method instead takes advantage of known time profile signatures related to the studied process. Results We explain the regulatory mechanisms of the inferred periodic classes with cis-regulatory descriptors that combine upstream sequence motifs with experimentally determined binding of transcription factors. By systematic statistical analysis we show that periodic classes are best explained by combinations of descriptors rather than single descriptors, and that different combinations correspond to periodic expression in different classes. We also find evidence for additive regulation in that the combinations of cis-regulatory descriptors associated with genes periodically expressed in fewer conditions are frequently subsets of combinations associated with genes periodically expression in more conditions. Finally, we demonstrate that our approach retrieves combinations that are more specific towards known cell-cycle related regulators than the frequently used clustering approach. Conclusion The results illustrate how a model-based approach to expression analysis may be particularly well suited to detect biologically relevant mechanisms. Our new approach makes it possible to provide more refined hypotheses about regulatory mechanisms of the cell cycle and it can easily be adjusted to reveal regulation of other, non-periodic, cellular processes.


Background
Epstein-Barr virus (EBV) primary infection usually occurs early in childhood until teens, and then persists throughout life as latent infection in a fraction of B-lymphocytes in more than 90% of adults.One adverse consequence of the infection is an increased tumour risk [1].There are some 170.000new EBV-positive tumours occurring annually in the global human population, of which half derive from the hematopoietic compartment and half from epithelial precursors.The tumour risk is thought to be intrinsic to the viral strategy for survival and spread.Indeed, the ability of the virus to transiently induce proliferation of latently infected B-lymphocytes results in an increased pool of infected cells.This induction of proliferation depends on the switch between viral latent programs in the cell, which can bring the cell from resting state, into active cell cycle and back to resting.If this switch gets out of balance, more proliferation leads to a higher load of virus-infected cells and hypothetically increases the risk for lymphoma development.The mechanisms controlling induction of proliferation are not understood.The most upstream event that can be identified until now is the switch between two viral promoters, the Q promoter (Qp) and the C promoter (Cp), which in turn determines expression of key regulatory viral proteins.Here we present an in silico study of this viral switch, involving one viral and one human protein.The aim is to investigate whether this simplified model can explain the availabel experimental data, as well as opening the field of EBV research to modeling.
The EBV genome is a 172 kbp long double-stranded DNA which during latent infections is maintained in the nucleus of the host cells as an extra-chromosomal circular episome.Twelve viral genes can be expressed in different combinations during latent viral infection, while the remaining 70 major open reading frames are expressed during the replicative, lytic cycle [2].During latency EBV genes are expressed in four programs, denoted latency 0, I, II and III.In latency III all 12 latency genes are expressed, including six nuclear proteins, EBNA-1-6, three membrane proteins (LMP-1, LMP-2A and LMP-2B), BART and two small non-translated RNAs (EBER 1 & 2).In latency II, the viral genes for EBNA-1, the three membrane proteins and the EBERs are expressed while in latency state 0/I, only LMP2a and variably EBNA-1 are expressed [2,3].B-lymphocytes in latency III are proliferating, driven by the viral gene products, while the remaining latency forms are in non-proliferating, resting cells.An essential viral protein is EBNA-1; responsible for viral replication, episome partitioning as well as functioning as regulator of gene transcription [4].EBNA-1 is therefore expressed in all latent programs, but the mRNA transcripts results from different promoters; Cp in latency III and Qp in latency I [5][6][7][8].It is believed that the two promoters are mutually exclusive due to different EBNA-1 levels between the latency programs [6,9].
The Cp transcripts produce all six EBNA-proteins from bicistronic mRNAs modulated by alternative splicing [10].Its activity is directed by several regulatory sequences of which the enhancer 'family of repeats' (FR) region is the major regulatory site.FR consists of multiple EBNA-1 and octamer binding sites, in an alternating pattern [11], and EBNA-1 binding to FR is required for transcriptional activation [6,12,13].The enhancement of Cp by EBNA-1 bound at FR follows a complicated pattern, where at least eight of the 20 present sites need to be occupied by EBNA-1 for full transcriptional activation [14,15].The human transcription factors Oct-2 and Oct-1 in association with Bob.1, have been shown to activate Cp, while Oct-2 in association with members of the Groucho (Grg/TLE)family of proteins acts as an inhibitor of Cp [11,16].Oct-1 is however not the prime candidate for operating the switch since it is not B-cell specific, like Oct-2, and moreover there is little or no difference in Oct-1 levels between latency I and III cells (Almqvist J et al.: Repression of Epstein-Barr virus enhancer Family of Repeats mediated transcription by Oct and Grg/TLE transcriptional regulators, suggest an involvement in switching between latency states, submitted).Oct-2 levels, on the contrary, differ between latency I and III (as shown in this study), indicating that Oct-2 is the dynamic inhibitory protein, even though it requires to bind FR in complex with Grg/TLE to be inhibitory [16].
Between FR and the initiation site there are other regulatory sequences such as an EBNA-2-dependent enhancer [17,18], GRE, a glucocorticoid-responsive enhancer [19], and sites for Egr, Sp1 and NF-Y transcription factors [20].Deletion analysis however indicated that the GRE is not necessary for transcription initiation, and the EBNA-2 dependent enhancer is not sufficient to activate Cp [21].
Moreover, Sp1 and NF-Y are two ubiquitously expressed proteins known to interact with each other [22], and NF-Y is suggested to play a role in basal transcriptional regulation [22,23], with binding sites in 30% of all eukaryotic promoters [24].Taken together, these experimental data indicate that EBNA-1 is the major activator of the C promoter.
In contrast to Cp, Qp governs expression of a monocistronic EBNA-1 transcript and is thought to be a housekeeping promoter of latency I since it lacks TATA sequence upstream of the initiation site, like many cellular housekeeping genes [25].Experimental studies of Qp regulation have revealed three types of transcription factors binding in the promoter region; EBNA-1, E2F and IRF.IRF factors have been shown to constituitively activate Qp, with the role of directing general transcription factors in the absence of a TATA box [9,26].The information concerning E2F regulation is less consistent, since E2F has been shown to both activate and inhibit Qp [27,28].There is however no doubt that Qp is negatively auto-regulated through binding of EBNA-1 to the two binding sites downstream of the initiation site [6,8,27,29], and that one occupied site is enough to block Qp activity [6].
Transcriptional control of genes is not trivial to model, especially not for eukaryotic systems with complex regulatory machinery.However, a main dynamic regulatory mechanism is the binding of inhibitory or stimulatory transcription factors to operator sites on the DNA.Assuming that the rates for binding and unbinding of transcription factors are much faster than the rates for open complex formation and elongation, these bindings are in homeostatic equilibrium with the instantaneous transcription factor concentrations [30,31].This is the basis for modelling gene transcription probabilities with thermodynamic models, as has successfully been used to describe genetic switches in coliphage systems [30][31][32][33].
In this paper we construct a physico-chemical model of the switch between EBV latency I and latency III.The model is based on the hypothesis that promoter control by EBNA-1 and Oct-2 is sufficient to switch between Cp and Qp usage in the two latency programs.Even though this is a simplification, the experimental data summarized above point to EBNA-1 as the major regulatory factor in controlling both promoter activities, and to Oct-2 in complex with Grg/TLE as a candidate for inhibitory control of Cp.In our model, Cp is thus switched on by EBNA-1 and switched off by Oct-2+Grg/TLE, while Qp is switched off by EBNA-1 and is otherwise on.Figure 1 gives a schematic overview of the model.
While this system belongs to the same general class of systems studied previously, a distinguishing feature of our system is the relatively large number of number of binding sites, and the complicated observed dependence of activation on the number of EBNA-1 bound.This theoretical approach to understand EBV-induced cell proliferation is to be seen as an extension of experimental studies and a first step towards a system level description of EBV infections.

Results
The two states latency I and latency III can be studied in human EBV-positive B-cell lines either derived from tumors or transformed in vitro by the virus.Cells in latency I without exception show lower levels of EBNA-1 and higher levels of Oct-2 compared cells in latency III ( [11]; Almqvist J et al.: Repression of Epstein-Barr virus enhancer Family of Repeats mediated transcription by Oct and Grg/ TLE transcriptional regulators, suggest an involvement in switching between latency programs, submitted).An example of a qualitative comparison of EBNA-1 and Oct-2, as detected in Western blot in the two types of cells, is shown in Figure 2a.We have estimated the quantitative differences between the intensities of the bands in Figure 2a, and in several different pairs of latency I and III cells (data not shown), using the ImageJ densitometry program (NIH).The EBNA-1 signal was between 2.8-4.7 times higher in latency III compared to latency I cells, while Oct-2 was 2.1-3.3 times higher in latency I cells compared to latency III.In vivo binding of both EBNA-1 and Oct-2 to FR in these two cell types can be demonstrated by a CHIPassay (Figure 2b).

Latency states in the model
Our model displays regions of mono-stability, with either latency I or III as the only stable state, and a region of bistability.For low Oct-2+Grg/TLE levels, only latency III is stable, while, for higher levels, there is a bistable region.This eventually gives way to the monostable latency I, at even higher Oct-2+Grg/TLE levels.The net effects of production from the two promoters together with dilution and decay can conveniently be visualized by a production potential.The negative derivative of this potential, with respect to the number of EBNA-1, correspond to the net production rate of EBNA-1.Figure 3 displays this production potential for three different levels of Oct-2+Grg/TLE.At intermediate Oct-2+Grg/TLE levels, the system is bistable, i.e. the production potential has two local minima (dashed line).Increasing the Oct-2+Grg/TLE level eventually eliminates latency III (dotted line), while decreasing Oct-2+Grg/TLE level eliminates latency I (solid line).
In our model, the level of EBNA-1 in the two stable states is dependent on the Oct-2+Grg/TLE levels, as can be seen from the positions of the minima in Figure 3.The model has a latency III level of EBNA-1 around 34,000 proteins (this value is calibrated to published data, see Methods) when no Oct-2+Grg/TLE complexes are present.This level decreases to around 14,000 after which point this latency state vanishes.
The EBNA-1 level in latency I varies less with Oct-2+Grg/ TLE levels, remaining around 700-900 proteins (all values given for parameters presented in Table 1).

Promoter activity
Figure 4 shows computed Cp and Qp activity in our model, as functions of EBNA-1 and Oct-2+Grg/TLE levels, demonstrating the remarkably different activity levels of the two promoters.The Cp activity (blue surface) is 60-100% of maximum for all Oct-2+Grg/TLE and EBNA-1 levels shown.The Qp activity (red surface) is mostly around 1%, except at very low EBNA-1 levels.Since the half-life of EBNA-1 is long (see Methods), we expect Qp to only be active at a low level in latency I. Our model indeed shows very low Qp activity in latency I, both in the mono and bistable regions.It estimates the Qp activity to be about 1% of its full capacity in latency I, and close to zero in latency III.Cp on the contrary is active at 40-100% of maximum transcription rate in latency III, with the lower activity in the bistable region.In monostable latency I, the Cp activity is essentially zero and below 1% active for latency I in the bistable region.The Cp transcription is strongly dependent on Oct-2+Grg/TLE levels.However, at the system volumes tested here, Oct-2+Grg/TLE levels do not greatly affect Qp activity.To summarize, the two promoters are mutually exclusive in the monostable regions, although their activity levels differ, and also for the latency III state in the bistable region.In the latency I state in the bistable region, the Cp and Qp activities are some-what comparable.Table 2 summarizes the Cp and Qp transcription activities presented above.

Stability of the latency programs
The system volume was in our study estimated to be 2 * 10 -13 l (see Methods), but was increased and decreased ten-fold in sensitivity tests.For each volume size, the stable latency I and III levels of EBNA-1 was computed for different dimerization dissociation constants for EBNA-1; 10 -8 M, 10 -9 M and 10 -10 M, and varying levels of Oct-2.Stable steady state levels were also computed for a fivefold lower Oct-2 affinity to FR.The stability of the two latency states and their robustness to parameter changes can be quantified by the externally imposed change on EBNA-1 levels that induces the system to transit from one state to the other.This measure is of course appropriate only in the bi-stable region, where both states exist.As The genetic switch Figure 1 The genetic switch.This illustration describes the regulatory mechanisms of the modeled switch between the C promoter (Cp) and the Q promoter (Qp) in Epstein-Barr virus.Transcripts from Cp are spliced into bicistronic mRNAs coding either for EBNA-1 (red) or one of the other EBNA proteins (blue), while transcripts from Qp only codes for EBNA-1.EBNA-1 itself acts as a negative regulator of Qp, binding downstream of the transcription start, and positively regulates transcription from Cp through binding upstream to the Family of Repeats sequence (FR).FR bind also Oct-2 molecules (green), in complex with Grg/ TLE, acting as inhibitory regulators.shown in Figure 5, for most parameters latency III is more stable than latency I, i.e., a larger change in EBNA-1 levels is needed to induce a transition from latency III to latency I than in the opposite direction.

[see Additional file 1]
As the left plot in Figure 5 demonstrates, the model is certainly sensitive to changes in affinity of Oct-2 complex to FR.With a fivefold increase in affinity of Oct-2 the position of the bistable region is shifted towards five times higher Oct-2 levels, but the relative properties of the latency states do not change significantly.The right plot in Figure 5 illustrates the sensitivity of the model to changes in volume and dimerization constant for EBNA-1.The impact of changes in the dimerization constant depends on the system volume, where the larger volume is less robust.Comparing the red and blue lines in Figure 5, the effect of increasing the dimerization 10-fold results in a shift of the bistable region towards higher Oct-2 levels.This shift is substantially larger for the same dimerization parameter increase in the larger volume (black and green lines).The model's response to changes in the volume parameter is primarily that the bistable region shrinks with increasing volume, and latency I EBNA-1 levels rise.

Inducing transitions
In Figure 6 we illustrate how our model behaves when responding to changes in Oct-2+Grg/TLE.Initially, the cell is in latency I, with 850 EBNA-1 proteins molecules per cell and 15000 Oct-2+Grg/TLE complexes per cell.Decreasing the level of the Oct-2+Grg/TLE complex by 30% quickly leads to an increase in Cp activity and hence an increase in EBNA-1 production.Within four to five days EBNA-1 levels have reached the levels of stable latency III.Transitions from latency III to latency I on the other hand, requires quite high levels of the Oct-2+Grg/ TLE complex for an extended period of time.The affinity of Oct-2 to FR being about 300 times weaker than the affinity of EBNA-1, Oct-2+Grg/TLE needs to be in large excess over EBNA-1 in order to fairly compete for binding to FR.Given a level of about 34,000 EBNA-1 in latency III, and the non-linear dependence of Cp activity on EBNA-1 concentration, one needs an Oct-2+Grg/TLE level of about 100,000 molecules, at least transiently, to lower EBNA-1 levels sufficiently in 5 days (solid lines in Figure 6).However, the switch can proceed slower with a lower level of Oct-2 for a longer time (dashed lines in Figure 6).Each Paramater is presented with a name, also used in the Methods section, and its value and corresponding reference puplication.*The three parameters that were altered in robustness tests.The range is given in Parenthesis.

Potential landscape
Figure 3 Potential landscape.The potential landscape of the EBV kinetic system for three different levels of Oct-2+Grg/TLE.The two local minima in the potential landscape corresponds to latency I, at low EBNA-1 levels, and latency III, at high EBNA-1 levels.Increasing the Oct-2 levels in the system shift the latency III minima towards latency I and eventually it disappear completely, leaving the system in a monostable latency I state.Both transitions have their lower limits in time, due to maximum production rate and the slow degradation rate of EBNA-1.

Discussion
EBV regulates fundamental properties of infected cells, and survives for long times in the host using different latency states, where we study the switch between the two extremes resulting in resting and proliferative states [1,3].
The molecular cell biology of the infection and EBVs association with various cancer types has been extensively studied.However, in order to fully understand how the virus survives in, and uses, its host cells we believe it is necessary to incorporate all information into a full, system level description.Our model of the switch between latency III and I is a first step towards such a model.The aim, besides opening the EBV research field to modeling, was to investigate whether EBNA-1 and Oct-2 regulation of Cp and Qp is sufficient to explain experimental results, and to identify important questions for further experimental studies.
A key effect in our model is the competitive binding of EBNA-1 and Oct-2+Grg/TLE to the Cp enhancer FR.Essential parameter values are hence binding affinities of EBNA-1 and Oct-2 to cognate sites.These binding affinities together determine the behaviour of the two involved promoters, and thereby also the switch between latency programs.Although it is well established that EBV can activate resting B-cells into proliferating blasts, the question how precisely this process is induced is still open.In our model, in order to keep the cell resting, using only Qp to produce EBNA-1, a high enough level of Oct-2+Grg/ TLE is necessary in the cell.When Oct-2+Grg/TLE levels drop, EBNA-1 can access FR and trigger Cp transcription.Since Cp is positively auto-regulated, working on a high production rate, EBNA-1 levels then quickly increase once Cp has been turned on.
We show that our model of the EBV latency switch displays expected features of the switch.At different Oct-2 concentrations, the switch is either monostable, and then exhibits states of either latency I or latency III, or is bistable.This behavior is robust to several parameter changes, although the boundaries and size of the stability domains vary.That the bistable region extends up to a higher level of Oct-2 in the system, when lowering Oct-2 affinity to FR, is to be expected, since a lower affinity requires a higher Oct-2 level to inhibit Cp transcription.Also, the fact that the model's sensitivity to changes in EBNA-1 dimerization is correlated with the system volume, can be explained.With a smaller system volume, the concentration of EBNA-1 is higher.With a concentration of EBNA-1 higher than, or of same order as, the dimerization dissociation constant, there will be no significant impact in changing this parameter.Besides correlating with dimerization of EBNA-1, a noticeable effect of volume increase is the shrinkage of the bistable region, with respect to Oct-2 levels.This correlates with the larger number of EBNA-  [35].This is compatible with our model as Grg/TLE can act as a chromatin remodelling protein by de-acetylating histones [36].We believe it likely that the transcription factor recruitment and regulation is upstream of the chromatin regulating events in the switch, i.e. creates the signal that leads to chromatin remodelling.
There is as yet no experimental system to directly operate this switch in B cells in vitro.However there are models in  [37].According to our model, for this transition to occur in 5 days, Oct-2 levels need to be elevated up to 10-fold.The higher the Oct-2 level already in latency III, the lower fold increase of Oct-2 is necessary.The increase is however only required temporarily, since stable latency I state can persist even for quite low Oct-2 levels.Although a qualitative comparison between our experimental and theoretical molecular levels is possible, it is at this stage not enough to determine whether the switch is carried out between monostable states or in the bistable region.Further experiments with silencing of Oct-2 RNA followed by time series measurements on EBNA-1 and Oct-2 levels, should shed light on this matter.
Noise plays an important role in gene regulation and in models of gene regulation when the levels of regulatory factors are low, and in primarily dynamic phenomena, when there is a choice between different pathways.Noise also decreases stability and robustness of steady states.For most parameter values in the bistable region, the latency I state is more easily perturbed than latency III.A consequence is that a switch from latency I to latency III can be expected to occur in response to smaller external stimuli.
However, since there is no compelling reason for the EBV switch to operate in the bistable region, the influence of noise is likely to be small.Furthermore, the protein levels are relatively high in both latency states, reducing the impact of stochastic effects.
Although the structure of the FR invites consideration of cooperative binding of EBNA-1 (and/or Oct-2+Grg/TLE) to these binding sites, no experimental data is at present

Promoter activity
Figure 4 Promoter activity.The probability of transcription from both promoters, Cp (blue) and Qp (red) as function of EBNA-1 and Oct-2 proteins in the system.This plot shows the case when Oct-2+Grg/TLE has the maximum binding affinity to FR, K dOFR = 2.5 nM.The Qp activity is remarkably low for most EBNA-1 levels compared with the Cp activity.available.Stress testing the model involving these additional effects is in progress, and will be reported elsewhere.The addition of cooperative binding would not affect the qualitative functions of the switch, while it might influence the quantitative measurements, such as protein levels and the switching times.Another issue in modelling gene regulatory systems is whether important entire elements are left out.While this cannot be deter-mined within a model, the propensity that it is the case can be adjudged from the model properties.We argue that our model is consistent with the knowledge present today, while it is clear that regulation of Oct-2 or Grg/TLE would logically complete the present model, as would information about other genetic feedback loops between the virus and its host cell.This then defines important tasks for the future.red lines are at the reference parameter values, in particular Oct-2 complex binding with affinity K dOFR = 2.5 nM and an EBNA-1 dimerization dissociation constant of 1 nM.Blue lines show a five-fold weaker Oct-2 affinity (K dOFR = 12.5 nM).Similar behaviour is then displayed at approximately five-fold higher Oct-2 level.Right figure: red solid and dotted lines at reference parameter values.Blue solid and dotted lines at a tenfold stronger EBNA-1 dimerization, and green solid and dotted lines at a tenfold greater volume.The latency III state is relatively robust towards either of these changes, while the latency I state changes more.Black solid and dotted lines show both a tenfold stronger EBNA-1 dimerization and a tenfold greater volume.This influences the latency III state more, essentially because EBNA-1 concentration in the latency III state is then comparable to the EBNA-1 dimerization dissociation constant.Lat I, V=2 −13 L,K doFR =2.5nM,K dEBNA =1nM

Conclusion
We have here shown that fundamental properties of the EBV infection, with high relevance for the virus as a risk factor in human cancer, can be modelled and understood mathematically.The model involves mechanisms which have, to our knowledge, only been used in models of simpler genetic systems before.However, also for this more complex system, these mechanisms nicely describe the dynamics and the properties of the viral genetic control.
A first conclusion is that EBNA-1 and Oct-2+Grg/TLE regulation of Cp and Qp is enough to establish mutually exclusive expression patterns, correlating with experimental data.Secondly, the switch manifests both mono-and bistable behavior, with the latency III state being markedly more robust.This in turn results in very different transition time scales between the latency states.
Our results and predictions moreover stress the importance, of the largely unknown regulation of the human transcription factor Oct-2 to the life-cycle of EBV.

Model base and main assumptions
As described in Background, our model is a simplification of the in vivo viral system, in the sense that it only involves promoter control by two transcriptional factors, EBNA-1 and Oct-2, with co-factors.The premise is that other regulatory factors documented in the literature, are less important as concerns the basic properties of the switch.
Regarding DNA binding of transcription factors to DNA, particularly to the FR region, the model assumes independent binding of both EBNA-1 and Oct-2.Moreover, bound proteins exclude only other bindings at that specific site, but do not affect neighbouring sites The octamer sites found in FR suggests a monomer binding of Oct-2, which we assume in our model, although the possibility of dimer bindings cannot be ruled out.The system is modelled with one kinetic equation, describing the dynamics of EBNA-1, coupled with thermodynamic equilibrium probabilities for transcription factor binding, to estimate transcription rates from Cp and Qp.As no regulation of Oct-2 levels by viral factors is known, Oct-2 level are essentially an external parameter to the model.The dynamics is modelled deterministically, disregarding noise in the transcription and translation (see Discussion).Details concerning parameter estimates are discussed in the following section.

Parameters
It is evident that all modelling results depend on the parameters used, resulting in a demand for exact parameters, or good estimates based on experimental studies.In Table 1 we list the parameters used in this study.These include the EBNA-1 dimerization constant, DNA binding dissociation constants, steady state level of EBNA-1 in latency III together with EBNA-1 half-life time, cell division time and the volume of the system.
Of central importance in this model are the binding affinities for our transcription factors to the FR region and Qp sites.The dissociation constant for EBNA-1 binding to FR and the relative dissociation constant for the Qp sites was determined by Ambinder et al., 1990.Their experimental results give a dissociation constant, K dEFR , for EBNA-1 from FR that is 15 pM and a K dQ from the Qp binding sites that is 14 times higher [39].Regarding EBNA-1 dimerization, EBNA-1 is known to be in dimer form in solution and bind DNA as a dimer [40], but the exact dimerization constant has however not been experimentally determined.We have used a K dE value of 1 nM as reference  value, but varied this parameter to examine its impact on the system.
Oct-2 belongs to the POU family of proteins, which share two distinct DNA-binding subdomains; a POU-specific domain and a POU homeodomain.The ATGC half of the octamer sites is recognized by the POU-specific domain, and the homeodomain binds to the AAAT half [41].Oct proteins therefore bind octamer sites as monomers, although they are able to form dimers and bind to similar palindromic sites known as MORE and PORE.In the FR region, the binding sites present are octamer-like; hence Oct-2 is assumed to bind as a monomer.The Oct-2 dissociation constant for the consensus octamer site, ATGCAAAT is known to be 2.5 nM [41].The octamer sites in FR are however not perfect, and most sites noticeably have an A6G variant in the last half of the sequence, and at least two base exchanges in the first half.From published mutational analyses of Oct-2 and Oct-1 octamer sites we can however deduce that the homeodomain generally has a higher binding affinity that the POU-specific domain, and that an A6G mutation does not dramatically affect the affinity [41,42].The dissociation constant of the Oct-2+Grg/TLE complex to sites within FR, K dOFR , is therefore here assumed to be 2.5 M but the effects of a higher dissociation constant was also tested.
The transcription rates for Cp and Qp and the translation rate for EBNA-1 transcripts are not experimentally known.As our model does not describe the transcription and translation steps individually but uses a total production rate, we incorporate both rates into one.Estimation of this production rate comes from experimentally determined steady state levels of EBNA-1 in latency III cell lines.Sternås et al. measured these EBNA-1 levels to range from 25,000 to 44,000 with an approximate mean value around 34,000 [43].This roughly corresponds to a production of 34,000 new EBNA-1 proteins every cell cycle, when production is balanced by dilution and degradation of EBNA-1.
The exact half-life of EBNA-1, τ 1/2 , is not known, but experiments indicate that it is at least 48 hours, probably due to the Gly-Ala repeat domain which inhibits proteasomal degradation [44,45].In our study we have hence chosen a half-life of 48 hours, and the cell division time is set to be 24 hours [46].The model assumes an homogeneous distribution of molecules in the system and does not include spatial movements between different cellular compartments.EBNA-1 is mostly located in the nucleus, due to its functions in transcriptional and translational control [47].Our system volume therefore corresponds to the nuclear volume of B-cells, a spherical volume with a diameter of 7 µm [48].

Physico-chemical model
The EBNA-1 dynamics is studied with a non-linear differential equation, (Eq. 1) describing production of EBNA-1 (E) from both promoters, Cp and Qp, protein degradation and also dilution, in the case of proliferating cells.The production from Cp is computed as the probability of gene transcription P c , times the rate of protein production r c (including transcription and translation), and likewise for Qp.Dilution is computed with a continuous dilution rate r dil , for cells that are in proliferating state.Degradation is computed with the continuous degradation rate r deg .The switch from resting to proliferating state is modelled instantaneously, meaning that the dilution rate is turned on as soon as Cp activity is higher than Qp activity, although in the real cell there is probably some delay.The total number of EBNA-1 molecules, E tot , in the system includes the free monomers, E, the dimers, E d , and the dimers bound specifically to DNA, E DNA , (Eq. 2).E DNA is computed as the mean value of bound EBNA-1 dimers and has to be calculated iteratively from the binding probabilities at each time step of a dynamic simulation, since it is dependent on the dimer concentration.The kinetic equation for EBNA-1 then reads; where E tot is defined as;

Transcriptional probability
The promoter transcription probability is computed from the probability of having a combination of transcription factors bound at the operator that allows for transcriptional activation.Full activation of Cp requires at least eight bound EBNA-1 dimers out of the 20 available sites at FR [14,15].The independent binding of EBNA- Q promoter activity, these are computed in the same manner, only dependent on EBNA-1 free energies, E eQ , and concentration [E] (Eq.4).The number of binding sites, N, is 20 for FR and 2 for Qp.

Experiments
Western Blot Protein concentration of nuclear extracts of latency I and III cells was determined by Bio-Rad Dc protein assay (Bio-Rad, Hercules, CA).The same amount of nuclear extract was loaded on every lane.Proteins were fractionated by 9% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to nitrocellulose membranes.After blocking for 1 h at room temperature with 5% milk, made up in PBS-0.1% Tween 20 (PBST), the membranes were probed overnight at 4°C with the following antibodies at indicated dilutions.Anti Beta-Actin (Sigma) was used as second protein control; Oct-2 (Santa Cruz Biotechnologies, Santa Cruz, Calif.) were used at 1:2000 and anti-EBNA1 used at 1:1000 (OTX-1, a kind gift from Jaap Middeldorp, Amsterdam Free university).The second antibody used was horseradish peroxidase-conjugated (HRP) and bound immunocomplexes were detected by enhanced chemiluminscence (ECL; Amersham Life Science, Little Chalfont, Buckinghamshire, United Kingdom).
Chromatin immunoprecipitation assay (ChIP) 2 * 10 6 EBV latency type I Rael cell or latency III CBMI-Ral-Sto cell were collected.Proteins were cross-linked to DNA by adding formaldehyde to 1% of culture medium.Cell pellets were collected after 10 min incubation at 37°C.The cell pellets were resuspended in SDS lysis buffer, DNA sheared to 500-1000 bp by sonicating the lysate.The sonicated samples were diluted 10 fold and 2 ug of immunoprecipitating antibody was added respectively: OTiX, mouse monoclonal antibody anti EBNA1, (kindly provided by Dr Jaap Middeldorp, Amsterdam), anti-OCT2 (C-20); rabbit polyclonal antibody, Santa cruz Biotechnology, California; Goat polyclonal antibody Santa cruz Biotechnology, California, and rotated at 4C overnight.100 ul of salmon sperm DNA/protein A agarose -50% slurry was added to collect antibody/protein/DNA complex.After washes and elution, the complexes were reverse crosslinked by adding sodium chloride to final concentration 0.2 M and kept at 65C for at least 6 hrs.Proteins were digested by proteinase K then DNA was purified by phenol/Chlorofome/ethanol precipitation.PCR primers were FRjz1S: TCCCTCTGGGAGAAGGGTAT and FRjz1A:TTT-TCGCTGCTTGTCCTTTT for family of repeats (FR).For control experiments to exclude non-specific binding of antibodies to DNA, the immunopreciptated DNA was analyzed with primers to beta-actin: AC-S: ATCATGTTTGAGACCTTCAA and AC-A: CATCTCTT-GCTCGAAGTCCA.DNA from latency III cell lysate (Mutu III) was immunopreciptated with 4 ug of anti-Oct2 or 4 ug of anti-EBNA1 and amplified by PCR with primers to beta-actin.Water and reaction without antibody was used as negative controls, while cell lysate DNA prior to immunopreciptation was used as positive control.

Figure 2
Protein levels and DNA binding.a) Western blot showing EBNA-1 and Oct-2 in latency I type cell Rael (I) and latency III type cell CBMI Ral-STO (III).Monoclonal antibodies were used to detect the proteins as described in Methods.b) Chromatin immunoprecipitation(CHiP)-assay of EBNA-1 and Oct-2 probed by PCR over the left hand of FR.The antibody-protein crosslinked and precipitated DNA was probed undiluted, diluted 1:10 and 1:100.The sources of the protein-FR complexes were the same cell lines as in a).Controls were immunpreciptation leaving out the first specific antibody, a cloned, purfied piece of FR (positive control) and PCR with no DNA (negative control) c) Control of non-specific binding of antibodies to DNA.DNA from latency III cell lysate was immunopreciptated with anti-Oct2 or anti-EBNA1, with negative and positive controls.Results show no significant background binding of the antibody.

Figure 5
Latency escapes and model robustness.Figure showing the minimum instantaneous change in EBNA-1 protein number necessary to switch from latency I to III and vice versa, as a function of the number of Oct-2+Grg/TLR proteins.For small numbers only the latency III state exists, while for large numbers only latency I state exists, compare Figure 3 and Figure 4. Left figure:

Table 2 : Computed C and Q promoter activities for the mono- and bi-stable states in both latencies.
Latency I and III levels of EBNA-1.This pdf-file includes tables with computed stable latency I and III levels of EBNA-1 for the different parameter sets tested.