{'fraction_non_alphanumeric': 0.02934778971321487, 'fraction_numerical': 0.008265569052638625, 'mean_word_length': 4.657705414615676, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Edge usage, motifs and regulatory logic for cell cycling genetic networks

The cell cycle is a tightly controlled process, yet its underlying genetic network shows marked differences across species. Which of the associated structural features follow solely from the ability to impose the appropriate gene expression patterns? We tackle this question in silico by examining the ensemble of all regulatory networks which satisfy the constraint of producing a given sequence of gene expressions. We focus on three cell cycle profiles coming from baker's yeast, fission yeast and mammals. First, we show that the networks in each of the ensembles use just a few interactions that are repeatedly reused as building blocks. Second, we find an enrichment in network motifs that is similar in the two yeast cell cycle systems investigated. These motifs do not have autonomous functions, but nevertheless they reveal a regulatory logic for cell cycling based on a feed-forward cascade of activating interactions.


I. INTRODUCTION
The cell cycle -biomass accumulation, DNA replication and cell division -is at the heart of life. Disruptions in this cycling, due to environmental changes or genetic defects, can lead to cell death or to uncontrolled proliferation such as in tumorigenesis [1]. It is thus not surprising that the cell cycle is tightly controlled in all organisms, allowing a stereotyped sequence of events to progress in an orderly and timely fashion. Various molecular types participate in this orderly progress, but each species has its own specific regulatory machinery [2,3].
Detailed investigations along with reconstructions of network interactions have allowed the building of quantitative cell cycle models based on ordinary differential equations (ODEs), in particular for budding yeast [4], fission yeast [5] and mammals [6]. The mathematical complexity of these systems is very high. To render them more comprehensible, simplifications have been proposed while preserving most of the qualitative aspects of their dynamics [7].
Reducing complexity still further, a number of authors have provided simplified models by discretising time and by replacing the continuous concentrations of molecular species by presence/absence (binary) variables. Such a Boolean approach has a long history going back to the pioneering work of Kauffman [8] (see also the book [9] and references therein) and since then it has been used in quite a number of systems. In particular it has been applied to the cell cycle in the three previously mentioned species [10][11][12][13][14] with a fair amount of success. This Boolean approach is complementary to that based on ODEs: it is much less powerful quantitatively but it allows for simpler model construction and interpretation, leading to qualitative insights into the generic aspects of regulatory rules. The present work is of a similar vein, but uses a more realistic model of molecular interactions and focuses on the topology of regulatory networks. Adjusting a dynamical model to reproduce observed cell cycle dynamics can be challenging but it is also an under-specified problem: because of the plethora of parameters, there are generally many different solutions [15]. Certainly, nature takes advantage of this freedom, but a consequence is that related species can have quite different regulatory circuits, impeding attempts to extract common regulatory principles. To overcome this difficulty, we have taken an in silico approach that characterizes all possible regulatory circuits that are compatible with given cell cycle expression dynamics. This route allows us to determine which features are always or mostly present in these in silico circuits and it also makes possible more meaningful comparative studies of different biological networks.
The present work builds on regulatory models developed in refs. [16,17] where we revealed several consequences of toy expression constraints on network structure. Here, we use the actual expression patterns found in three biological systems. Our modeling of the regulatory networks uses thermodynamic considerations to describe the underlying molecular interactions. Interestingly, the actual number of effective inter-molecular interactions becomes an emergent feature, driven by biophysical constraints coupled to an appropriately defined fitness function mimicking the selection pressure operating in nature. Based on this choice of genetic regulatory network (GRN) modeling, we sample by Markov Chain Monte Carlo (MCMC) the space of all GRN that provide the proper gene expression patterns. We do so for the three ensembles of GRN where the gene expression profiles are taken from budding yeast, fission yeast and mammals derived in refs. [10,12,13]. For each, we extract the structural properties of the sampled GRN, from edge usage to the presence of network motifs. Finally, by comparing to the case of an idealized cell cycle, we reveal an underlying common regulatory logic in which feed-forward activating cascades play a central role. This logic is apparent when the associated motifs are externally driven, but the motifs on their own do not implement autonomous dynamics.
In the next section we briefly sketch the main aspects of our model, stressing the changes from our previous work [16,17] to deal with new aspects of the problem. In section III we present our results. We conclude in section IV.

A. Biophysical modeling of molecular interactions
In previous work [16,17] we developed a biophysical framework for modeling the regulatory interactions between a transcription factor (TF) and its DNA binding site. For cell cycling systems, a number of the actors, such as cyclins, are not TFs, so our previous framework must be generalized. The starting point is a network whose nodes and edges refer to molecules and activation/repression processes, be-they transcriptional, posttranscriptional, translational or post-translational. All interactions involve specific binding sites or surfaces that allow two molecules to physically establish contact and bind through the addition of small forces. Almost all of the facing elements (atoms, bases, amino acids, ...) have to "match" for the two molecules to bind. Following ref. [18], we consider that each of the facing regions is specified by an ordered lists of symbols that can be thought of as a string or a table characterizing that type of molecule and its site or surface dedicated to that binding partner. Then the mutual binding energy of two molecules is taken to be additive in the number of mismatches between their strings, each mismatch contributing a penalty ε.
Such a framework is easily justified for TF-DNA binding energies and leads to a thermodynamic formula for the probability of occupation of a given DNA binding site in the presence of n j TFs of type j [19]: Here d ij is the number of TF-DNA mismatches. The derivation of Eq.(1) in ref. [20] shows that it is of rather general validity, so we shall apply it to all of the molecular species in our cell cycle system. Note that P ij depends strongly on d ij and is appreciable only when d ij is small. With a realistic choice of model parameters, a small d ij is a priori very unprobable, so that a functional molecular contact is also very unprobable.

B. Time dynamics of expression levels
Refs. [16,17] considered transcriptional regulation; we have to generalize the framework to any kind of regulation because of the different molecular types driving the cell cycle. Let S j (t) ∈ [0, 1] (j = 1, . . . N ) denote the average level of expression of the molecular species j at time t, normalized by its maximum value. The N expression levels {S j (t)} will be referred to as the phenotype at time t. We are interested in processes of the type that must reproduce as well as possible a target cell cycle sequence, i.e., the sequence deduced from experimental measurements. For example, in the case of the fission yeast S. Pombe, this sequence pattern is given in Table I as specified in ref. [12].  [12]). Time runs from the top to the bottom. Successive, idealized expression levels, are either 0 or 1, at each time step. We keep the terminology from [12]: the initial state's expression is given by the first line of the table (referred to as a "START" state) and the last state is a fixed point associated with the cell size check point. When a critical cell size is reached, a signal (implemented by having a so-called "Start" node turn on) triggers the cycle again. The Cdc2 (Cdc2*) stands for Cdc2/Cdc13 (Cdc2/Cdc13*), where the * sign indicates the highly activated form of the complex [12]. Similarly Wee1 stands for the Wee1/Mik1 kinases.
Following refs. [16,17], the discrete time dynamics for the expression levels are based on a mean field approximation which neglects cooperative binding effects. Explicitly, we have where the P ij were defined in Eq.(1). In the present expression, j runs over activating and j ′ over inhibitory interactions respectively. It is easy to see that Eq.(3) embodies a simple logic for each S i : (i) each activator or repressor acts independently; (ii) the binding of at least one activator is required for activation of the target; (iii) the binding of even a single repressor is sufficient to veto the activation and so will turn off the expression of the target. For the sake of simplicity we set n j = nS j (t) in Eq.(1), with a value of n common to all molecular species. For our simulations, we set n = 1000, but we have checked that the results presented depend little on the magnitude of this parameter. Thus the fluctuations in the number of molecules of a given type will not matter much, further justifying our mean field approach.
Notice that Eq.(1) holds provided the system is in equilibrium. Thus, in writing Eq.(3) it is assumed that the binding-unbinding reaction is fast compared to that for full activation/repression. However, the rate of a binding-unbinding reaction is proportional to the concentration of reactants and the equilibration time becomes very large for small concentrations. To avoid this limitation, attempts to generate very small concentrations should be ignored. Furthermore, the unphysical assumption that there is always enough time to reach equilibrium sometimes leads to pathologies. In particular, the system suffers from an instability: even a relatively weak interaction has an a priori capacity of generating a spurious progressive amplification of physically unsignificant (because very low) expression levels. To avoid this bad behavior and improve the robustness of the model, we have recourse to a phenomenological correction, introducing a threshold H to modify Eq.(3) whenever the expression level is too small: In effect, S min can be thought of imposing a basal rate or some level of leakiness on the transcription while the threshold H ensures that expression levels stay low unless the input signal is sufficiently activatory. In practice we set S min = 1/n = 0.001 as if there were just one molecule of that species and H is set to H = 0.01. Then, S i (t) can leave the minimal level only when at least one d ij is small enough to have a dynamical significance [17].

C. Computational implementation
For each i and each j, (i, j), the interaction strength W ij is given by Eq.(1) where d ij can be considered to be the mismatch [18] between two strings of length L. In the simulation, rather than storing the 2 character strings for each such interaction, we simply store the binary difference string (also of length L) specifying which entries match (1) and which mismatch (0). There is also a sign associated with each interaction: + for an activator, -for a repressor. Mutations can change the sign of an interaction or they can transform a match into a mismatch and vice versa. Following the procedure used for TFs [16,17], and to keep the modeling simple, we consider that the original strings use a four letter code. Then the probability that a mutation replaces a match by a mismatch is 3/4, while the probability that a mismatch is converted into a match is 1/4, embodying the fact that it is easier to have a mismatch than a match. (Note that even for protein-protein interactions, this rule can be motivated by the fact that molecular changes are often the result of a point mutation at the DNA level.) We will refer to the matrix d ij and to the associated set of signs as the genotype. Typically d ij = d ji since these two numbers of mismatches correspond to different processes: the former to the activation of i by j, while the latter to the activation of j by i. In these two processes different active sites have to match.
For a given genotype, the dynamical system modeling the cell cycle is initialized in the "START" state and then the trajectory is generated by iterating Eq.(4). At each time step we compute the distance D(t) between the vectors of expression levels associated with the actual sequence of the dynamical system at hand and the target one (for fission yeast, the target sequence is given in Table I). The total distance D T = T t=0 D(t) is then used to define the fitness F of the genotype via where the parameter f gives the strength of the selection pressure to maintain "good" expression profiles. Indeed, deviations from the target expression profile are likely to be deleterious for the proper functioning of a cell and its progeny. Note that this fitness defines a weight for each genotype and thus an ensemble where fitness plays the role of a Boltzmann factor. Finally, to sample genotypes according to their fitness, we use Markov Chain Monte Carlo (MCMC). The fitness enters the associated Metropolis test which accepts or rejects a trial change (mutation) to the genotype. We always start with a completely random network and thermalize it until it has a good fitness. Then, we produce a long sequence of several thousand GRN genotypes, each being separated by 200 sweeps (a sweep is a series of N 2 L attempted mutations and sign flips). To deal with possible fragmentation of the search space [14], we checked that the statistical properties of the resulting genotypes were the same for several independent simulations, i.e. each simulation being initialized with a different random genotype.
The parameters values used in our computations are L = 12, ε = 1.75, f = 20, H = 0.01 and n = 1000. Compared to ref. [17], H is a new parameter, ε is slightly smaller while the remaining are the same. Interestingly, it is essential that L and ε be such that the a priori probability of a small mismatch is very low. Once that constraint is taken into account, the model results are robust with respect to the variation of the parameters.
The Metropolis test described above forces D T to be relatively small if f is large, but the MCMC is inefficient in that regime, forcing us to work with not too large values of f . Then in the T × N table of numbers S i (t), there may appear a few that are "ambiguous", that is that are away from the target value which in this study is either 0 or 1. For example an expression level may have the value 0.3 in a place where the target value is 0. Forcing a better agreement by increasing f is not computationally feasible. To overcome this difficulty, from our large MCMC sample, we have selected genotypes where such ambiguities are absent, imposing that expressed (respectively unexpressed) molecular types have in fact levels above 0.6 (respectively below 0.2). This selection is done mostly for esthetic reasons, the properties of the full sample are essentially the same because of the low frequency of these ambiguities.

A. Essential interactions
The networks generated by the MCMC are characterized at a microscopic level by weights W ij quantifying the regulatory interaction strength of a molecule of type j on the expression level of molecules of type i. Just as had been found in our previous work [16,17], most W ij remain negligible, but a small fraction have values that are far above the background level. These interactions are the ones important for reproducing the target sequence of expression states. We formalize this property as follows: if setting W ij = 0 causes the Metropolis test to reject the change 5 times in a row, we say that the interaction W ij is "essential". The set of essential interactions of a network allows for a summary representation that is well suited for visualising the influences of each molecular type. A genotype's set of essential interactions defines a directed network, the essential network, from which one may extract informative regulatory patterns.
It is well known that biological GRNs are sparse. As shown in refs. [16,17], our model very naturally generates sparse essential networks because relevanti.e. smallmismatches are a priori unprobable and arise only under strong selection pressure. Then simply because of en- tropy, the model strongly favors low numbers of essential interactions. However, this argument is qualitative only. Let us thus compare quantitatively the numbers of edges produced within the model to the numbers in the biological data as inferred by the different authors [10,12,13].
Hereafter, we will refer to those reconstructed networks as the "wild-type" networks. In Table II we give the average number of essential interactions along with the value in the wild-type network for baker's yeast [10], fission yeast [12] and mammals [13]. (Note: we do not count the "self-degradation" interactions added by hand by some of these authors which are specific to their modeling of the discrete time dynamics of the expression levels.) The mammalian wild-type network appears to use significantly more interactions than our model but the two yeast species seem to be in good agreement with our model. To quantify this, we have measured the fluctuations in each ensemble. As shown in Fig. 1, for the two yeast ensembles the distributions of the number of edges are bell shaped and narrow, so in fact the fluctuations are very mild. Furthermore, the difference between the wildtype numbers and the expected value in these ensembles are not much more than one standard deviation, a result that is quite non-trivial. Most of results presented hereafter concern the two species of yeast. The fact that our model does much better for yeasts than for mammals may not be a coincidence. Indeed, although baker's yeast and fission yeast are highly divergent evolutionarily, much more so than any mammals amongst themselves, they are both uni-cellular organisms. It is well known that the regulatory control of gene expression is intrinsically more complex in multi-cellular organisms, so a posteriori one may draw some satisfaction from the fact that our simple model seems to do well when regulation is simple but does much less well when regulation involves very complex mechanisms. In the remaining subsections we will omit results referring to the mammalian cell cycle. It is sufficient to say that they are systematically off the data.

B. Edge usage
Now consider the edge usage in the essential networks. The frequency with which there is an essential interaction from node j to node i can be represented by a square matrix with entries between 0 and 1. We have computed these matrices for the three cell cycles, and in Fig. 2 we display the results separately for activating and inhibitory interactions in baker's yeast. Particularly, there are 9 activators (7 common with the wildtype) and 7 repressors (4 common with the wild-type) almost always present in the ensemble for baker's yeast (cf. dark grey entries in Fig. 2). Further, the results can be compared with frequencies of activatory/inhibitory interactions from [11], where 6 activators and 4 repressors (all common with the wild-type) are found to be absolutely required for a network to produce the cell-cycle process. Specifically, among these absolutely required interactions, 8 edges are almost always present and two edges are present in more than 80% of networks in our ensemble.  Fig. 2 but here refer to the wildtype network of fission yeast as specified in ref. [12] (see also Fig. 6).
On a more general level there is a manifest dichotomy: certain interactions are very rarely if ever present, while others are often if not always present. Furthermore, there is a clear structure to the matrix: in the line just below the main diagonal, activating interactions are very frequent. A different pattern arises for the inhibitory interactions: there, all the entries near the diagonal are always absent, and elements 4 or 5 steps shifted from the diagonal are frequently or always present.
Qualitatively the same pattern arises in fission yeast but it is noisier (Fig. 3), and in the case of the mammalian cell cycle which involves only 7 genes, there are simply remanants of these patterns.

C. Overlaps with the wild-type networks
If an essential interaction of an in silico GRN is present in the wild-type network, we shall shall call this interaction a "hit". What fraction of the essential interactions of a GRN are hits? If the fraction is close to 1, there is a very strong overlap between that GRN and the wild-type network. As can be seen in Fig. 4, in practice the fraction varies from network to network but typically takes values in the range from 50% to 70% in both yeast species. The mean fractions are 65% and 57% for fission and baker's yeast so clearly there is a strong overlap between the in Histograms for the fraction of essential interactions that are "hits" in the baker's (left) and fission (right) yeast ensemble. We see that networks in our ensembles have roughly 50% to 70% of interactions in common with the wild-type networks. Notice, that the wild-type networks are not part of our ensemble.
silico and wild-type networks. The highest fraction of hits found in our MCMC is 73% for baker's yeast and 76% for fission yeast. Notice that this fraction is close to the overlap of hits with respect to all links in the wildtype network, since in our model E wt and E ess are similar for both yeast species (see Tab. II). To compare visually the associated networks to the corresponding wild-type networks, we display them in Fig. 5 for baker's yeast and in Fig. 6 for fission yeast. The main differences arise for the Sic1 and Clb1,2 genes which in the wild-type have higher degree and more mutually inhibitory interactions than in the GRN generated by MCMC. Such interaction pairs are quite natural in the context of protein-protein interactions and thus probably reveal a shortcoming of our modeling framework.

D. A multitude of different essential networks
The MCMC generates a very large number of genotypes from which we construct essential networks. In effect, each essential network corresponds to a different way interactions can be specified so as to reproduce the target gene expression dynamics. Many genotypes give rise to the same essential network, but are there many different essential networks? The answer of course is yes, we find hundreds of different such networks. One may then ask whether some arise much more frequently than others. On the right hand side of Fig. 5 (respectively Fig. 6) we show the most frequently found essential network for baker's yeast (respectively fission yeast); these arise with frequencies between 2 and 3%. Their characteristics are similar to those of the essential networks having the highest overlaps with the wild-type networks. These most frequent essential networks arise hundreds of times more frequently than the rarest ones. In Fig. 7 we show the associated rank histogram for the concurrence frequencies on a log-log plot. The fat tail in this distribution is roughly compatible with a power law. Such a shape for a rank histogram arises in a number of other systems, in particular in neutral networks where many genotypes give rise to the same phenotype.

E. Network motifs
In the terminology of Alon [21], a network motif in a (biological) graph is a subgraph that is present at a higher frequency than expected at random. The term "at random" is generally defined using an ensemble of graphs where each node is constrained to have the same in and out degree as in the biological graph of interest. In practice this ensemble is generated using a randomization of the edges of the biological graph [22,23] and that is what we have done too.
Given the genotypes in each of our ensembles, produced by the MCMC, we have extracted the motifs present in the associated essential networks. The different motifs we find are represented graphically in Fig. 8. For two nodes, we find no significantly over expressed motifs. For three nodes, we find that for both yeast cell cycle ensembles several of the coherent feed-forward loops are over represented. In particular for the baker's yeast ensemble, these over-representation factors are close to 10. Finally, for four nodes, we find the presence of numerous motifs: two kinds of incoherent diamonds, two kinds of frustrated four-point loops and incoherent bifans. The detailed frequencies of each of the different motifs in these ensembles are given in Tab. III.
Finally, instead of working within the in silico GRNs produced by our ensembles, one can check for motifs in the wild-type networks. We find that the motifs are not all the same as in the associated ensemble. Among the strongest differences, the RR motif consisting of two mutually inhibitory nodes is strongly over-represented in the wild-type networks but not in the in silico networks. Differences exist also for a number of other motifs. For example, in the wild-type network of baker's yeast, the numbers are: motif RR = 3; A (C1-FFL) = 2, C (C4-FFL) = 1, D (neg 3-loop) = 2. In the case of fission yeast, the numbers are motif RR = 4; diamonds: G = 2, H = 2; bi-fan: I = 1. These numbers can be compared to those in the GRN ensemble using Tab. III.

F. Regulatory logic
The pattern of edge usage obtained in Fig. 2 depends on the order of the genes. The order we chose makes the successive expression states resemble as much as possible a left to right shifting block of 1's. As is visible in Table I, even with this "best" choice, irregularities remain, suggesting one consider an idealized case. We thus replaced the irregular shifting block by a regular one and then repeated our methodology on these idealized cell cycles. Interestingly, the associated essential networks have edge usage and motifs similar to the ones previously described for the wild-type cell cycles. Furthermore, the activating interactions typically form one long feed forward cascade. To illustrate this phenomenon, we show in Figure 9 the most frequently found essential network when idealizing the target sequence of Table I. All activating (non-self) interactions follow each other from "START" to the last expression state. The function of such a cascade is clear but if one extracts an associated motif, say three consecutive activating interactions, it has no autonomous dynamics. To understand the functioning of such a motif, one has to externally drive it, providing inputs and initial conditions that will initiate a "block" of on genes. The motif will then propagate as a cascade of falling dominoes [24] until the block of 1's is pushed out and the whole region consists of 0's.

IV. DISCUSSION AND CONCLUSION
It is quite remarkable that a model, based on rather elementary thermodynamic and probabilistic principles, succeeds in reproducing much of the cell cycle network topology found in both baker's yeast and fission yeast, two highly divergent organisms. Our in silico approach is based on considering all possible networks subject to the constraint that they give rise to the same gene expression patterns as the experimentally observed systems. In effect, this approach allows one to infer how network regulatory architectures are constrained by the network "function". Within this framework, our modeling was able to give a good indication of the number of interac- tions and of the actual edges that are used in the two yeast species. We also considered higher order features such as network motifs; this allowed us to reveal an ubiquitous design principle in the form of an activatory cascade that is realized perfectly within an idealized cell cycle and is manifest in the less regular cell cycles studied here. However in the case of the mammalian cell cycle network our results are less impressive to say the least. As already mentioned, this is not really a surprise. In view of the simplicity of the model, its failures are not devoid of interest: they indicate where a more sophisticated modeling is necessary. It is reassuring that the model does better for uni-cellular organisms than multi-cellular organisms which have far more levels of regulatory control; the contrary would be suspect! We were seeking to determine which aspects of GRN structure are expected from simple modeling and argu-  ments of a generic nature. Hence, we have assumed that all interactions are a priori possible and only allowed the expression phenotype to constrain the networks. But it is clear that ignoring chemistry in a biological process is somewhat limiting and may be responsible for some of the features arising in the biological networks but not in our in silico networks. A likely candidate for this is the motif consisting of 2 mutually repressing molecular species: it is very frequent in the two yeast cell cycle net-works but not in our ensembles of GRN. This absence occurs in spite of the fact that for other phenotypic constraints, namely having fixed point expression patterns, our model does produce such motifs. Thus it would be of interest to see the impact on our results of forbidding interactions known to be unlikely on biochemical grounds, but such a study is beyond the scope of the present work.