T Preece

Y Mao

Sustainability of dioecious and hermaphrodite populations on a lattice

Keywords: sexual reproduction, lattice, Allee effect

Sexual reproduction may be divided into two main categories; hermaphroditism and dioecy (Botany)/gonochorism (Zoology). Simultaneous hermaphrodites can function in both male and female roles whereas a dioecious/gonochorist population consists of distinct male and female individuals. Mean-field calculations, which ignore spatial aspects, suggest that self-incompatible hermaphrodites should have a twofold advantage over dioecious population when reproduction is limited by mating encounters. By use of stochastic spatial simulations we demonstrate that hermaphroditism has an even greater advantage when local interactions are considered. This result provides further support for the observation that hermaphroditism is associated with sedentary species, such as plants and animals with poor mate search efficiency. We also investigate the finite size effects associated with the well-known quadratic contact process.

A c c e p t e d m a n u s c r i p t 1. Introduction

Dioecy describes sexually reproducing plant species in which there are two distinct sexes, male and female. One individual of each sex is required for reproduction. The equivalent mating system in animal species is referred to as gonochorism. For simplicity, we will use the botanical term dioecy in this article. In contrast to dioecy, simultaneous hermaphrodites have the ability to reproduce in both female and male roles simultaneously.

The Allee effect describes a positive relationship between individual fitness and the density of conspecifics [START_REF] Odum | A note on the stable point of populations[END_REF]. The most typical single-sex (hermaphrodite) model of the Allee effect [START_REF] Lewis | Allee dynamics and the spread of invading organisms[END_REF] is represented by,

dh dt = ah(h -b)(c -h), (1)
where h is the population density, and a, b and c are positive constants (b < c). Equation (1) has three equilibria; h = 0 and h = c are stable, h = b is unstable. Below the threshold h = b the population converges to h = 0, so becomes extinct. Above the threshold h = b the population converges to the stable equilibrium, h = c. The most cited cause of the Allee effect is the difficulty of finding mates in sexually reproducing species [START_REF] Boukal | Single-species models of the allee effect: Extinction boundaries, sex ratios and mate encounters[END_REF].

It may be argued that hermaphroditism, as a mating strategy, is more advantageous than dioecy when reproduction is limited by mating encounters.

Consider a female in a dioecious population with a one to one ratio of females to males, for each encounter with another individual, the female encounters a male with probability 1/2. In a hermaphrodite population, every encounter

A c c e p t e d m a n u s c r i p t

of a hermaphrodite would be with a compatible mating partner. All else being equal (encounter rate, number of off-spring produced per mating pair), the hermaphrodite population can reproduce at twice the rate of the dioecious population. This reasoning provides support for the observation that hermaphroditism is often associated with sedentary species, such as plants and animals with poor mate search efficiency [START_REF] Eppley | Moving to mate: the evolution of separate and combined sexes in multicellular organisms[END_REF].

However this mean-field argument neglects spatial detail, which becomes important in systems with short range interactions between constituents and limited mixing caused by poor locomotion.

Stochastic spatial simulations have been extensively applied in chemical, ecological and sociological systems [START_REF] Durrett | Stochastic spatial models[END_REF][START_REF] Durrett | Stochastic spatial models a users guide to ecological applications[END_REF].

The Standard Contact Process (SCP) can be used to model populations that reproduce asexually. In the SCP lattice sites can be either occupied or vacant.

Individuals die at a rate λ and give birth at a rate r. Offspring are placed on a neighbouring site with a probability given by some displacement kernel.

In the Quadratic Contact Process (QCP) [START_REF] Durrett | Stochastic spatial models[END_REF], an occupied site becomes vacant at a rate λ, the same as in the SCP. However, in the QCP it takes two particles to create a new one. For this reason the QCP is often referred to as the "sexual reproduction process". A vacant site becomes occupied at a rate j/4, where j is the number of diagonally adjacent pairs of occupied neighbours.

The SCP and QCP are often identified as Schlögl's 1st and 2nd models for autocatalytic kinetics, respectively (Schlogl, 1972). For the 1st model, One artifact of the QCP is that if a population is contained inside a rectangle, it can never spread outside [START_REF] Durrett | Stochastic spatial models[END_REF]. This is due to the fact that an empty site needs at least two diagonally adjacent neighbours to become occupied. For lattices of finite size, this feature also implies that complete columns or rows of empty sites in the lattice cannot become occupied.

X
To assess the impact of these finite size effects, another spatial model similar to the QCP described above is introduced here and will be referred to as the Modified Quadratic Contact Process (MQCP). For the MQCP a hermaphrodite (female) produces offspring with rate proportional to the number of occupied (male) neighbouring sites. An offspring is placed on a randomly chosen neighbouring site. If the site is empty, the offspring takes residence. If the site is occupied the offspring is deleted. This approach is based on the sexual reproduction model of Stewart-Cox et al. (2005) and has the advantage that a population can grow outside a rectangle to which it is confined. A similar model, which also overcomes this finite size effect specific to the QCP, has been presented by Windus and Jensen (2007).

In the following section, the QCP is extended to accommodate a dioecious population and the MQCP is explained in more detail for both hermaphrodite and dioecious mating systems. Then, in section 3, mean-field analyses are presented for each model and section 4 presents results of our lattice simulations.

A c c e p t e d m a n u s c r i p t

Models

For the stochastic spatial simulations, all populations reside on a twodimensional square lattice of side length L, with periodic boundary conditions. For a hermaphrodite population, each lattice site can be occupied by a hermaphrodite (H) or be empty ( ). For a dioecious population each lattice site can be occupied by a male (M ), a female (F ) or be empty ( ). For both mating systems, an occupied state becomes vacant at a rate λ (deathrate).

In the hermaphrodite QCP [START_REF] Durrett | Stochastic spatial models[END_REF], a vacant site becomes occupied at rate equal to j 1 r 1 /4, where j 1 is the number of diagonally adjacent pairs of occupied neighbours and r 1 is the birth rate. Thus, j 1 = 0,0,2 and 4 for 0, 1, 3 and 4 occupied nearest neighbour sites respectively. j 1 = 1 (0) for two occupied nearest neighbours which are diagonally adjacent (on opposite sides of the vacant site).

In the dioecious QCP, a vacant site becomes occupied at rate equal to j 2 r 2 /4, where j 2 is the number of compatible diagonally adjacent pairs of occupied neighbours and r 2 is the birthrate. Thus, an (M, M ) or (F, F ) diagonally adjacent pair does not increment j 2 . The gender of the offspring is set to male with probability α and female with probability (1α).

In the hermaphrodite MQCP, an occupied site produces an offspring at a rate equal to k 1 r 3 /4, where k 1 is the number of occupied nearest neighbours and r 3 is the birth rate. A nearest neighbour site is chosen randomly to place the offspring on. If the chosen site is empty the offspring is placed on the site. If the site is occupied the offspring is deleted.

In the dioecious MQCP, a female site produces an offspring at a rate equal to k 2 r 4 /4, where k 2 is the number of male nearest neighbours and r 4

A c c e p t e d m a n u s c r i p t

is the birth rate. A nearest neighbour site is chosen randomly to place the offspring on. If the chosen site is empty the offspring is placed on the site.

If the site is occupied the offspring is deleted. The gender of the offspring is set to male with with probability α and female with probability (1α).

There exists a subtle difference between the MQCP and QCP. If

r 1 = r 2
in the QCP, a hermaphrodite pair of diagonally adjacent neighbours are capable of producing offspring at the same rate as a male-female pair of diagonally adjacent neighbours. If r 3 = r 4 in the MQCP, a hermaphrodite pair of nearest neighbours can produce offspring at twice the rate of a malefemale pair of nearest neighbours. For example, consider an isolated pair of nearest neighbours, both members of a hermaphrodite pair can produce offspring at a rate r 3 /4 where as only the female can produce offspring in a dioecious pair, she does this at a rate r 4 /4. Thus, when the probability of encountering a compatible mate is included, the MQCP proposes a fourfold advantage of hermaphrodite over dioecious populations. With r 3 = r 4 , the MQCP proposes that a female and a hermaphrodite can produce offspring at equal rates provided a mating partner is in the vicinity. Depending on the species in question, either assumption may be valid. Fortunately, scaling the reproductive values appropriately can overcome this difference. For example, setting r 1 = 2r 2 in the QCP, proposes that females and hermaphrodites produce offspring at the same rate provided a mating partner is available.

Setting r 4 = 2r 3 , proposes that a hermaphrodite pair produce offspring at a rate equal to that of a male-female pair in the MQCP.

A c c e p t e d m a n u s c r i p t

3. Mean-field analysis

Hermaphrodite population

The mean-field kinetics for the hermaphrodite QCP population is given by

dh dt = -λh + r 1 h 2 (1 -h), ( 2
)
where h is the density of hermaphrodites residing on the lattice. Equilibrium points for equation ( 2) are plotted in figure 1 as a function of r 1 /λ. By appropriate scaling of time (τ = λt), the number of parameters in equation ( 2)

can be reduced to one (r 1 /λ) . Thus, figure 1 is unique for a given value of r 1 /λ. Two non-trivial equilibria (one node, one repellor) exist for r 1 /λ > 4.

We will refer to this as the survival/extinction phase. If the initial population density is below the unstable equilibrium, growth is negative and the population converges on the trivial equilibrium, h = 0. If the initial population density is above the unstable equilibrium the population converges on the stable non-trivial equilibrium. Only the trivial equilibrium exists for r 1 /λ < 4, thus extinction of the population always occurs.

Dioecious population

The mean-field kinetics of the dioecious QCP are given by,

dm dt = -λm + 2r 2 mf (1 -m -f )α (3) df dt = -λf + 2r 2 mf (1 -m -f )(1 -α) ( 4 )
where m and f are the densities of males and females respectively. Substituting r 2 = 3r 4 /8 into equations ( 3) and ( 4

A c c e p t e d m a n u s c r i p t

Results of lattice simulation

Previous simulations [START_REF] Liu | Quadratic contact process: Phase separation with interface-orientation-dependent equistability[END_REF][START_REF] Guo | Generic two-phase coexistence, relaxation kinetics, and interface propagation in the quadratic contact process: Simulation studies[END_REF] have shown that the QCP exhibits a discontinuous phase transition between an active state (finite population) and an absorbing state (extinction). It was also shown that coexistence of stable active and absorbing phases occurs over a finite parameter space. For any parameter set (λ, r) in this coexistence region, "droplets" of the absorbing state embedded in the active cannot grow indefinitely but rather die out, even though the the absorbing state is stable.

Similarly, "droplets" of the active state embedded in the absorbing state are short lived.

Figure 3 plots population density in the hermaphrodite QCP and MQCP models as a function of r/λ. For each point in the plot, the population was initialised with a full lattice and evolved for 10 3 generations. The densities plotted are averaged over the last 10 3 Monte Carlo Steps (MCSs) of the simulation, where 1 MCS = L 2 site selections. A transition from a stable active phase to a stable absorbing phase occurs at r 1 /λ ≈ 11. The location of the QCP transition is in agreement with previous simulations [START_REF] Liu | Quadratic contact process: Phase separation with interface-orientation-dependent equistability[END_REF][START_REF] Guo | Generic two-phase coexistence, relaxation kinetics, and interface propagation in the quadratic contact process: Simulation studies[END_REF] that have shown the coexistence region between r 1 /λ ≈ 10.5 and r 1 /λ ≈ 11.5 for the QCP. As all simulations in figure 3 were initialised with a fully populated lattice only the stable active state of the coexistence region was observed.

For the MQCP a transition from a stable active phase to a stable absorbing phase occurs at r 3 /λ ≈ 6 (figure 3) . This is in contrast to the mean-field approximation, which predicts that the transition in the MQCP model should occur at a larger value of r/λ larger than in the QCP (see

A c c e p t e d m a n u s c r i p t

figure 1). These results demonstrate that the finite size effects specific to the QCP do have a significant impact on stability.

Investigations of the dynamics of the interface between active and absorbing phases in the QCP [START_REF] Liu | Quadratic contact process: Phase separation with interface-orientation-dependent equistability[END_REF][START_REF] Guo | Generic two-phase coexistence, relaxation kinetics, and interface propagation in the quadratic contact process: Simulation studies[END_REF] showed that the active phase grows fastest (recedes slowest) for an interface orientation of 45

•
to the rows of the square lattice. This rate of growth (reduction) decreases (increases) as the angle of orientation is changed from 45 • . For an interface orientation of 0 • or 90 • to the lattice rows, the active phase cannot invade the absorbing phase. This is due to the fact that an empty site requires at least two diagonally adjacent nearest neighbours to become occupied. This results in the feature that if a population is contained inside a rectangle it can never give birth outside the rectangle [START_REF] Durrett | Stochastic spatial models[END_REF] and for lattices of finite size, complete columns or rows of empty sites cannot become occupied [START_REF] Liu | Quadratic contact process: Phase separation with interface-orientation-dependent equistability[END_REF]. For an interface orientation of 0 • or 90 • , in the MQCP, the active state can invade the absorbing state. This is expected to be the reason for the contrast between simulation and mean-field behaviour in the QCP and MQCP.

Due to finite size effects, analysis of size dependent behaviour is required to accurately determine the stability regions [START_REF] Liu | Quadratic contact process: Phase separation with interface-orientation-dependent equistability[END_REF][START_REF] Guo | Generic two-phase coexistence, relaxation kinetics, and interface propagation in the quadratic contact process: Simulation studies[END_REF]. The exact values of these parameters are not calculated here but the stability of the two phases is clarified by plotting waiting time to extinction of the population as a function of lattice size L (figures 4 and 5). Each point in figures 4 and 5 is the mean waiting time to extinction from 100 simulation runs; each run was initialised with a fully occupied lattice and lattice sites were set to male with probability 0.5 for dioecious populations. Figure 4 A c c e p t e d m a n u s c r i p t

shows that in the stable active phase extinction time increases exponentially with L. The solid lines in figure 4 are least squares fits for the data, their equations are given by;

τ = exp[2.98 + 0.471L],
for QCP hermaphrodite, (5)

τ = exp[2.02 + 0.263L],
for QCP dioecious, ( 6)

τ = exp[-192 + 1.00L],
for MQCP hermaphrodite, ( 7)

τ = exp[0.802 + 0.450L], for MQCP dioecious, (8)
where τ = λt is referred to as the time in generations (t is in MCSs).

In the stable absorbing phase, figure 5 shows that the extinction time increases logarithmically with L. The solid lines in figure 4 are least squares fits for the data, their equations are given by; 2006) and may be due to a metastable region, similar to the one observed for the QCP [START_REF] Liu | Quadratic contact process: Phase separation with interface-orientation-dependent equistability[END_REF][START_REF] Guo | Generic two-phase coexistence, relaxation kinetics, and interface propagation in the quadratic contact process: Simulation studies[END_REF]. Again, the stability of the two phases were clarified by plotting waiting time to These results demonstrate that, as expected, finite size effects specific to the QCP do have a significant impact on the stability of dioecious populations on a lattice.

τ = 11.2 + 18.4 ln (L) , for QCP hermaphrodite, (9) τ
The transition from the active to absorbing phase in the QCP occurs at a value of r/λ 30 times larger in the dioecious population than the hermaphrodite population (figures (3) and ( 6)). The transition from the active to absorbing phase in the MQCP occurs at a value of r/λ 15 times larger in the dioecious population than the hermaphrodite population (figures

(3) and ( 6)). These results demonstrate more drastic differences between the hermaphrodite and dioecious mating systems than predicted by mean-field analysis (see figure 1). If birth rates were scaled such that a hermaphrodite pair produce offspring at an equal rate to a male-female pair (r 4 = 2r 3 ), all else being equal, this would correspond to a transition from absorbing to active state at a value of r/λ ≈ 7.5 times larger in the dioecious MQCP model than the hermaphrodite MQCP.

Though the models predict different values for the advantage of single sex over dioecious populations on a lattice, both models suggest that a hermaphrodite population is stable over a large range of low birth rates for which a dioecious population is not sustainable. The relatively large birth rates required to sustain a dioecious population are likely to be due to local variation in the concentration of (M -M ) or (F -F ) pairs of cells, which hinder the birth process. As offspring are placed on nearest neighbour sites, areas with an optimum sex ratio (1/2) cannot compensate for areas with a poor sex ratio (close to 0 or 1). This can aid the propagation of the A c c e p t e d m a n u s c r i p t absorbing state. The observed advantage of hermaphroditism at low birth rates could contribute to the success of the hermaphrodite mating system in sedentary species. This provides additional support for the observations that hermaphroditism should be associated with sedentary species, such as plants and animals with poor mate search efficiency [START_REF] Eppley | Moving to mate: the evolution of separate and combined sexes in multicellular organisms[END_REF].

Conclusions

By means of stochastic spatial simulations we have shown hermaphrodite populations can be sustainable at birth rates much smaller than the birth rates required for a sustainable dioecious population. This implies that a hermaphrodite population has a much greater reproductive advantage over a dioecious population than the the twofold difference predicted by mean-field analysis. This result provides even greater theoretical support for the observations that hermaphroditism should be associated with sedentary species, such as plants and animals with poor mate search efficiency.

Our simulation results also show that finite size effects specific to the QCP introduce significant artifacts which can be overcome by simple modifications with the MQCP.

Acknowledgements

We would like to thank the Royal Society, the EPSRC and the University of Nottingham for financial support. For helpful discussions we would like to thank J. P. Garrahan and A. Davison. : Sum of male and female densities, m + f , versus r/λ after 10 3 generations for QCP ( ) and MQCP (+). Each value of r/λ is initialised with a fully populated lattice where each lattice site was set to male with a probability 0.5. The densities plotted were averaged over the last 10 3 MCS of the simulation. L = 100.

→

  2X and X → ; for the 2nd model, 2X → 3X and X → , where X → denotes the annihilation of a particle. The mean-field kinetics are ) for the 1st (2nd) model suggesting a continuous (discontinuous) transition to the absorbing state, . The mean-field kinetics of the QCP can lead to equation (1).

  figure 1 overlap for the QCP. It was also explained in section 2 that the hermaphrodite population can grow at four times the rate of the dioecious population in the MQCP model. If we let r 4 = 4r 3 , hermaphrodite and dioecious curves in figure 1 overlap for the MQCP. As shown by figure 2, deviations away from the optimum sex ratio of 1/2 can hinder the growth of the population. This can have a great effect when considering discrete spatial populations with local interactions.

Figure 6

 6 Figure6plots the population density (m + f ) as a function of r/λ for the dioecious population. A discontinuous phase transition between the active and absorbing states is observed at r 2 /λ ≈ 300 for the QCP and r 4 /λ ≈ 100 for the MQCP. Given more time, the two QCP points at densities close to 0.3 fall to zero. These slow relaxations close to the transition were also observed byTainaka et al. (2006) and may be due to a metastable region, similar to

  population as a function of lattice size L (figures 4 and 5).

Figure 1 :

 1 Figure1: Mean-field approximation equilibrium densities for; QCP hermaphrodite ( ), QCP dioecious ( ), MQCP hermaphrodite (×) and MQCP dioecious (+). The gender ratio of new born offspring α = 1/2 for the dioecious populations. Solid (dashed) lines represent stable (unstable) equilibria.

Figure 2 :

 2 Figure 2: Mean-field approximation for QCP: typical dynamics for a population of males and females. Plot shown for r 2 = 1.0, λ = 0.1, α = 0.5. Showing the stable equilibrium (filled circle) and unstable equilibrium (empty circle). The dashed lines are a projection of the stable (m + f ≈ 0.85) and unstable (m + f ≈ 0.1) equilibria for the hermaphrodite population with the same birth and death rates.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Hermaphrodite population density after 10 3 generations for lattice simulation of QCP ( ) and MQCP (×). L = 100.

  Figure 6: Sum of male and female densities, m + f , versus r/λ after 10 3 generations for


Xiaofeng Dai

Shannon Healy

Olli Yli-Harja

Andre S Ribeiro

Tuning cell differentiation patterns and single cell dynamics by regulating proteins' functionalities in a toggle switch

Keywords: Functionality regulation, Toggle Switch, Functional protein, Non-functional protein, Cell differentiation

We investigate how the regulation of protein multi-functionalities affect the dynamics of a stochastic model of a toggle switch and the differentiation pattern of cell population regulated by the switch. We study the effects of loss of functionality in DNA-binding and repression and the involvement in differentiation pathway choice. First is shown how the patterns of cell differentiation differ, when each of these functionalities is fully non-functional. Next, tuning the fraction of non-functional proteins regarding the ability to bind DNA is shown to allow fine tuning of the switch and cell differentiation pattern dynamics. Finally, biasing the probability of functionality of the two proteins biases the dynamics of the switch and cell differentiation patterns, especially when transcription factors retain the ability to bind DNA but have lost the ability to repress gene expression. Our results suggest that, besides transcriptional and translational levels of regulation, activation of functionalities in multi-functional proteins are an important regulator of gene networks.

Introduction

Protein expression and activation is highly regulated in cells. At several stages in a cell's lifetime, a protein's function and expression patterns can be altered by post-translational modifications, such as reversible phosphory-

A c c e p t e d m a n u s c r i p t

lation [START_REF] Ashcroft | Regulation of p53 Function and Stability by Phosphorylation[END_REF]. For example, p53 contains multiple phosphorylation sites which activates the protein, leading to cell cycle arrest and apoptosis [START_REF] Bates | p53 in signalling checkpoint arrest or apoptosis[END_REF].

A protein can have multiple functions, such as degrading a substance and repressing a gene's expression [START_REF] Korpela | A Recombinant Escherichia coli sensor strain for the detection of tetracyclines[END_REF]. These multiple functions typically arise from separable domains within the protein, which in some cases can act and be regulated independently of each other [START_REF] Rachfal | Structural and functional properties of CCN proteins[END_REF][START_REF] Thirumalai | Kinetics of folding of protein and RNA[END_REF]. As such, one domain may have reduced activity [START_REF] Kim | A quantitative model of transcription factoractivated gene expression[END_REF], or be nonfunctional, while the other is active, causing the protein to be "partially nonfunctional". Such partially nonfunctional proteins may have unique roles in cellular processes and growth, by manipulating the dynamics of the cell's gene regulatory network (GRN) [START_REF] Andrews | Selective disadvantage of non-Functional protein synthesis in Escherichia coli[END_REF].

Reversible phosphorylation is a critical mechanism in eukaryotic cells by which several of the functionalities of transcription factors can precisely tuned in a rapid fashion, in response to changes in the cellular environment [START_REF] Holmberg | Multisite phosphorylation provides sophisticated regulation of transcription factors[END_REF]. Phosphorylation regulation of the activity of transcription factors occurs via at least six mechanisms [START_REF] Holmberg | Multisite phosphorylation provides sophisticated regulation of transcription factors[END_REF]: (i) it can alter the degradation rate of the transcription factor [START_REF] Appella | Signaling to p53: breaking the posttranslational modification code[END_REF], (ii) it affects the cellular localization of a transcription factor [START_REF] Macian | Partners in transcription: NFAT and AP-1[END_REF]. (iii) it affects protein-protein interaction [START_REF] Whitmarsh | Regulation of transcription factor function by phosphorylation[END_REF], (iv) it affects the DNA-binding activity [START_REF] Whitmarsh | Regulation of transcription factor function by phosphorylation[END_REF], (v) regulates the transcriptional activity of factors such as HSF1 [START_REF] Kline | Repression of the heat shock factor 1 transcriptional activation domain is modulated by constitutive phosphorylation[END_REF], and (vi) can indirectly trigger the inactivation of transcription factors [START_REF] Chi | Negative regulation of Gcn4 and Msn2 transcription factors by Srb10 cyclin-dependent kinase[END_REF].

A well characterized example demonstrating these principles of regulation of partial functionality of a gene regulatory protein is the repressor Even-skipped (Eve) [START_REF] Um | The Transcriptional Repressor Even-skipped Interacts Directly with TATA-Binding Protein[END_REF]. Eve is a sequence-specific homeodomain repressor that is involved in pattern formation during early embryogenesis [START_REF] Smith | A spatially dynamic cohort of regulatory genes in the endomesodermal gene network of the sea urchin embryo[END_REF][START_REF] Um | The Transcriptional Repressor Even-skipped Interacts Directly with TATA-Binding Protein[END_REF]. It is comprised of two separable domains: a DNA-binding domain and a repressor domain, which directly interacts with the TATA-binding protein (TBP) [START_REF] Um | The Transcriptional Repressor Even-skipped Interacts Directly with TATA-Binding Protein[END_REF]. Phosphorylation of the repressor domain by glycogen-synthase kinase-3 was found to reduce Eve's interaction with TBP, thereby allowing transcription. Eve's DNA binding ability was not affected. While regulation of Eve's temporal and spatial expression through an extensive GRN has been characterized, it has not yet been considered how a partially nonfunctional Eve influences downstream events in embryogenesis [START_REF] Smith | A spatially dynamic cohort of regulatory genes in the endomesodermal gene network of the sea urchin embryo[END_REF].

It is therefore of importance to consider the dynamics of protein functionality when studying the dynamics of GRN, because the degree of functionality may be subject to deliberate modifications as a means to control the propensity of gene expression. The objective of this report is therefore to examine the pattern of outcomes generated from differential regulation of protein's functionality. To better understand the potential role of the regulation of

A c c e p t e d m a n u s c r i p t

proteins' functionality on GRN's dynamics, one can use computational models. Given the small number of most proteins in cells [START_REF] Bon | Many expressed genes in bacteria and yeast are transcribed only once per cell cycle[END_REF], the regulation of proteins functionality is likely to be a stochastic process. The dynamics of gene expression is noisy [START_REF] Kaern | Stochasticity in gene expression: from theories to phenotypes[END_REF][START_REF] Elowitz | Stochastic gene expression in a single cell[END_REF][START_REF] Golding | Real-time kinetics of gene activity in individual bacteria[END_REF][START_REF] Raser | Noise in gene expression: origins, consequences, and control[END_REF], affecting the dynamics of gene networks and cell differentiation pathway selection [START_REF] Arkin | Stochastic kinetic analysis of developmental pathway bifurcation in phage λ infected Escherichia coli cells[END_REF][START_REF] Mcadams | Stochastic mechanisms in gene expression[END_REF][START_REF] Ozbudak | Regulation of noise in the expression of a single gene[END_REF]. Temporal delays that inevitably exist in gene expression are also important regulators of GRN's dynamics [START_REF] Bratsun | Delayinduced stochastic oscillations in gene regulation[END_REF][START_REF] Lewis | Autoinhibition with transcriptional delay: a simple mechanism for the zebra fish somitogenesis oscillator[END_REF]. For example, the time needed to form the promoter complex [START_REF] Mcclure | Rate-limiting steps in RNA chain initiation[END_REF] affects the toggling frequency of genetic toggle switches (TS) [START_REF] Ribeiro | Dynamics of a two-dimensional model of cell tissues with coupled stochastic gene networks[END_REF]. The delayed Stochastic Simulation Algorithm (delayed SSA) [START_REF] Roussel | Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression[END_REF] (a modified version of the original SSA [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]), able to match gene expression at the single RNA and protein level [START_REF] Zhu | Studying genetic regulatory networks at the molecular level: delayed reaction stochastic models[END_REF][START_REF] Yu | Probing gene expression in live cells, one protein molecule at a time[END_REF], allows modeling multiple time delayed reactions [START_REF] Ribeiro | A general modelling strategy for gene regulatory networks with stochastic dynamics[END_REF], and thus is used here to simulate the dynamics of the GRN of each cell.

It has long been hypothesized that cell differentiation is based on bistable genetic sub-circuits that control many downstream genes [START_REF] Monod | Teleonomic mechanisms in cellular metabolism, growth, and differentiation[END_REF], and it has been shown that the TS, a network of two mutually repressing genes, can be used by cells to adopt different phenotypes [START_REF] Arkin | Stochastic kinetic analysis of developmental pathway bifurcation in phage λ infected Escherichia coli cells[END_REF] and as decision circuits of differentiation pathway selection [START_REF] Huang | Cell fates as a high-dimensional attractor states of a complex gene regulatory network[END_REF]. Importantly, the dynamical behavior of toggle switches driving differentiation have been shown to be regulatable [START_REF] Chang | Transcriptome-wide noise controls lineage choice in mammalian progenitor cells[END_REF] and their internal noise has been shown to play a crucial role in mouse hemapoetic cell line.

Recent analysis of human myeloid progenitor cells has provided experimental evidence for the existence of multiple metastable states [START_REF] Chang | Multistable and multistep dynamics in neutrophil differentiation[END_REF], consistent with the assumption of being driven by coupled toggle switches that can switch between noisy attractors [START_REF] Ribeiro | Noisy attractors and ergodic sets in models of gene regulatory networks[END_REF]. Also, the probably to switch of the TS has been shown to vary by increasing the dose of a stimulus [START_REF] Becskei | Positive feedback in eukaryotic gene networks: cell differentiation by graded to binary response conversion[END_REF][START_REF] Gardner | Construction of a genetic toggle switch in Escherichia coli[END_REF]. In [START_REF] Chang | Transcriptome-wide noise controls lineage choice in mammalian progenitor cells[END_REF] it was concluded that clonal heterogeneity of gene expression level is not only due to independent noise in the expression of individual genes, but reflects metastable states of a fluctuating transcriptome that may govern the reversible, stochastic priming of multipotent progenitor cells in cell fate decision.

Therefore, to study the effects of the regulation of proteins' functionalities in patterns of cell differentiation we model, in each cell of a cell population, a delayed stochastic TS [START_REF] Ribeiro | Dynamics and evolution of stochastic bistable gene networks with sensing in fluctuating environments[END_REF]. The model here used was built from the model first proposed and confronted with real data in [START_REF] Ribeiro | A general modelling strategy for gene regulatory networks with stochastic dynamics[END_REF].

We introduce in the model the possibility for each protein to be nonfunctional in some of its functionalities. The fractions of functional and nonfunctional proteins are controlled by the rates of competing stochastic chemi-

A c c e p t e d m a n u s c r i p t

cal reactions. To observe how this fraction affects cell differentiation patterns we introduce a differentiation mechanism in each cell that is stochastically activated at a random moment (drawn from an exponential distribution) of the cell's lifetime. We analyze the effects of varying: i) the 'transformation rate', i.e., the rate at which premature proteins are converted into functional or non-functional proteins, ii) the 'non-functional fraction', i.e., the fraction of premature proteins that become non-functional, and iii) the 'sub-system biases', i.e., the bias in the protein production of the sub-systems of the TS (each gene and its products are named a 'sub-system' of the TS).

In this work we conduct a preliminary analysis of the potential effect on gene network dynamics of the regulation of the combinatorial activation of proteins multi-functionalities. In each of the three scenarios, we study four cases to investigate how tuning the 'ratio' of functionality (ratio between functional and non-functional proteins) affects cells' dynamics and their differentiation patterns. The four cases are the possible combinations of preserved functionalities in partially non-functional proteins (for simplicity, partially non-functional proteins are here referred to as non-functional proteins). Namely, we consider the cases where non-functional proteins are unable to alter the gene's expression rate but can, or not, preserve the ability to bind to the correct binding site of the other gene, and to affect cells' differentiation choices.

We find that the regulation of proteins multi-functionalities is likely to be an important regulatory mechanism of the dynamics of GRN. This is further enhanced by the fact that such regulation occurs as a fast respond to changes in the cellular environment.

Methods

Model

To study the dynamics of single cells and differentiation patterns of cell populations, we simulate populations of identical cells at the single cell level. Each cell's dynamics is driven by the delayed SSA [START_REF] Roussel | Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression[END_REF], that has a waitlist to store delayed output events. The waitlist is a list of elements (e.g., proteins and RNA being produced), each to be released after a certain time interval (also stored in the waiting list). Simulations are done in SGN Sim [START_REF] Ribeiro | SGN Sim, a stochastic genetic networks simulator[END_REF].

The cell model has three sets of reactions. The first set defines the TS where proteins need to dimerize before binding to promoters, the second set

A c c e p t e d m a n u s c r i p t

regulates the transformation of premature proteins into functional and nonfunctional proteins and whether non-functional protein dimers can bind to their transcription factor binding site (TFBS) or not, and, finally, the third set of reactions regulates the choice of differentiation pathway and whether non-functional proteins are involved in this process or not. The first set of reactions is present in all cells, and the different combinations of reaction sets 2 and 3 form the four tested cases.

Toggle switch

The TS model [START_REF] Ribeiro | Dynamics and evolution of stochastic bistable gene networks with sensing in fluctuating environments[END_REF] consists of reactions 1 to 8, where i = 1, 2 if only index i is present, or i, j = 1, 2 (i = j) if both indices are present.

P ro i + RNAp kt -→ P ro i (τ 1 ) + RNAp(τ 2 ) + R i (τ 1 ) (1)
R i + Rib ktr -→ R i (τ 3 ) + Rib(τ 4 ) + P i (τ 5 , τ 5std ) (2)
P f,j + P f,j k dimer -→ D f,j (3)
D f,j k undimer -→ P f,j + P f,j (4)
P ro i + D f,j krep -- kunrep P ro i D f,j (5)
P ro i D f,j dp -→ P ro i (6) R i d rbs -→ ∅ (7)
P i dp -→ ∅ (8)
Gene expression is modeled by two multi-delayed reactions, one for transcription (reaction 1) and one for translation (reaction 2), to account for the time duration of the multiple steps involved in these processes [START_REF] Bratsun | Delayinduced stochastic oscillations in gene regulation[END_REF][START_REF] Mcclure | Rate-limiting steps in RNA chain initiation[END_REF][START_REF] Ribeiro | A general modelling strategy for gene regulatory networks with stochastic dynamics[END_REF][START_REF] Zhu | Studying genetic regulatory networks at the molecular level: delayed reaction stochastic models[END_REF]. P ro i is the promoter of gene i, RNAp is an RNA polymerase, Rib is a ribosome, R i is the ribosome binding site of the transcribed RNA, P i is a translated premature protein, P f,i is a functional protein and D f,i is a functional protein dimer. The delays (τ 1 to τ 5 ) account for the duration of the processes in transcription and translation [START_REF] Ribeiro | A general modelling strategy for gene regulatory networks with stochastic dynamics[END_REF]. When a product X has a delay τ , represented by X(τ ), it implies that when the reaction occurs, it takes τ seconds after that for X to appear in the cell.

In our model, mature functional proteins have to form dimers (reaction 3) before binding to the other gene's promoter region and repress its expres-

A c c e p t e d m a n u s c r i p t

sion (reaction 5). Dimmerization makes more likely that proteins are multifunctional and that the activation/inactivation of its functionalities can be done in an independent fashion. Unrepression can occur by the backward reaction of 5 or by the decay of the dimer bound to the promoter, via reaction 6. Degradation of RNAs, proteins and their dimers occur via reactions 7, 8 and 4, respectively.

Values of delays (in s) were extracted from measurements of the gene tsrvenus in E. coli [START_REF] Yu | Probing gene expression in live cells, one protein molecule at a time[END_REF]. τ 1 accounts for the promoter open complex formation (completion of transcription initiation). It follows a normal distribution with a mean of 40 s and standard deviation of 4s according to measurements in a lac promoter [START_REF] Mcclure | Rate-limiting steps in RNA chain initiation[END_REF][START_REF] Lutz | Dissecting the functional program of Escherichia coli promoters: the combined mode of action of Lac repressor and AraC activator[END_REF], used to drive the expression of the tsr-venus. Since the average rate of transcriptional elongation in E. coli is ∼ 50 nt/s, τ 2 is, on average, 90 s (τ 2 = τ 1 + 2500/50), given 2500 nt is the gene length [START_REF] Yu | Probing gene expression in live cells, one protein molecule at a time[END_REF]. Note that the delay associated with the time needed to form the ribosome binding site R i is also τ 1 (and not τ 2 as for a complete RNA molecule), since the time to form R i is almost identical to the promoter clearance time given that we model gene expression in prokaryotes, where translation is initiated when the ribosome binding site is formed [START_REF] Ribeiro | A general modelling strategy for gene regulatory networks with stochastic dynamics[END_REF][START_REF] Zhu | Studying genetic regulatory networks at the molecular level: delayed reaction stochastic models[END_REF]. In translation, the time needed for R i clearance, τ 3 , is set to 2 s [START_REF] Draper | Translation initiation, In Escherichia coli and Salmonella[END_REF]. Since the average translation rate is 15 amino acids per second, τ 4 is computed as τ 4 = τ 3 + 2500/45 = 58 s. The post-translational protein assembly process was observed to take 420 ± 140 s in [START_REF] Yu | Probing gene expression in live cells, one protein molecule at a time[END_REF], which are set, in accordance, as the values for τ 5 and τ 5std .

Each 'cell' is initialized without proteins or RNAs, one promoter of each gene, 40 RNA polymerases [START_REF] Record | Escherichia coli and Salmonella[END_REF], and 100 ribosomes [START_REF] Zhu | Studying genetic regulatory networks at the molecular level: delayed reaction stochastic models[END_REF]. The basic transcription rate is set to k t = 0.01 causing the expected time for an RNAp to bind to the promoter to be 1 kt×RN Ap = 2.5 [START_REF] Kennell | Transcription and translation initiation frequencies of the Escherichia coli lac operon[END_REF]. Reaction rates k tr = 0.00042 and d rbs = 0.01 are set to match the average observed number of proteins produced per RNA E. coli, i.e., 4.2 [START_REF] Yu | Probing gene expression in live cells, one protein molecule at a time[END_REF]. Proteins decay is set to d p = 0.001, ten times slower than mRNA [START_REF] Creighton | Proteins[END_REF][START_REF] Wang | Precision and functional specificity in mRNA decay[END_REF]. We additionally set heuristically k rep = 1, k unrep = 0.1, k dimer = 5 • 10 -4 and k undimer = 5 • 10 -5 to attain a toggling frequency within the order of magnitude of experimental observations [START_REF] Gardner | Construction of a genetic toggle switch in Escherichia coli[END_REF]. The transformation rate k trans is tunable.

Transformation of premature proteins into functional and non-functional

Transformation of premature proteins into functional and non-functional proteins are modeled by reactions 9 to 12:

P i ktrans -→ P f,j (9)
A c c e p t e d m a n u s c r i p t

P i ktrans×f i -→ P nf,i (10)
P f,i k d -→ ∅ (11)
P nf,i k d -→ ∅ (12)
where P i , P f,i and P nf,i denote premature, functional and non-functional proteins, respectively. A premature protein P i becomes either P f,i or P nf,i via, respectively, reactions 9 and 10. Reactions 11 and 12 model the degradation of functional and non-functional proteins. Tuning f i allows adjusting the expected fraction of premature proteins that become non-functional and functional. Note that in premature proteins none of its functionalities has been activated, while non-functional proteins are, as mentioned, only partially non-functional. The conversion of premature proteins into functional or non-functional is assumed here to occur at the post-translational modification stage. Both non-functional and functional proteins are able to be part of a dimer. The loss of the two basic functionalities of dimers, i.e., molecular recognition and enzymatic activity [START_REF] Mathews | Biochemistry[END_REF][START_REF] Merkel | Functional protein microarrays: just how functional are they?[END_REF], are modeled by 'case N ' and 'case B ', respectively. In 'case N ', non-functional dimers cannot bind to promoters while in 'case B ', non-functional dimers, although able to bind to the genes' TFBS, do not repress gene expression. This is modeled by adding reactions for binding and unbinding of D nf,i to P ro j , and reactions for the transcription and the decay of the complex P ro j D nf,i in reactions 15, 16, and 17, respectively. 16)

P f,j + P nf,j k dimer -→ D nf,j (13)
D nf,j k undimer -→ P f,j + P nf,j (14)
P ro i + D nf,j krep -- kunrep P ro i D nf,j (15)
P ro i D nf,j + RNAp kt -→ P ro i D nf,j (τ 1 ) + RNAP (τ 2 ) + R i (τ 1 ) (
P ro i D nf,j k d -→ P ro i (17)
Reaction 16 allows the gene to express while bound to non-functional dimers. This is expected to affect the dynamics of the TS, if significant amounts of non-functional dimers are present in the cell, since they decrease the probability that the promoter is repressed by functional dimers.

A c c e p t e d m a n u s c r i p t

The following assumptions in the model should be justified. First, we assume that k rep and k unrep have equal values for both the functional and non-functional dimers in the case where non-functional dimers can bind to DNA but do not repress transcription. This occurs, e.g., in the example described of the repressor Even-skipped (Eve) where phosphorylation of the repressor domain by glycogen-synthase kinase-3 allows transcription without affecting Eve's DNA binding ability [START_REF] Um | The Transcriptional Repressor Even-skipped Interacts Directly with TATA-Binding Protein[END_REF]. It is further noted that the domain where the repressor binds to DNA can be upstream to the elongation sequence and thus, a non-functional dimer bound to DNA does not necessarily affect transcription (e.g. by blocking the passage of the RNA polymerase) even though it competes with the functional dimer for the same promoter binding site (e.g., shadow enhancer sites [START_REF] Hong | Shadow enhancers as a source of evolutionary novelty[END_REF]). Also, we do not assume the formation of dimers from two non-functional proteins. These dimers would not be able to bind to DNA and so the only effect they would have on the dynamics would be reducing the number of functional-nonfunctional dimers. In this scenario, to reach the same ratio between functional and non-functional dimers (composed by one functional and one nonfunctional protein) as when not assuming dimers of two nonfunctional proteins, one would just vary f further, to qualitatively have the same variations in the differentiation patterns reported using our model. Due to the fact that we vary f in a very wide range, we allowed this simplification in the model.

Finally, we assume in one scenario that the kinetic parameters of transcription (k t ) are not significantly affected by the binding of a non-functional dimer, which can nevertheless competitively inhibit the binding of functional dimers. Eve is an example of such situations since, as described, has two separable domains and the reduction of its ability to repress transcription does not affect Eve's DNA binding ability [START_REF] Um | The Transcriptional Repressor Even-skipped Interacts Directly with TATA-Binding Protein[END_REF].

Cell differentiation mechanism

Our simplified model of pathway choice of cell differentiation assumes that the two dimers expressed by the TS are also transcription factors of another gene (or genes), responsible for determining the differentiation pathway [START_REF] Ribeiro | Variability of the distribution of differentiation pathway choices regulated by a multipotent delayed stochastic switch[END_REF]. We assume two cases, one where non-functional dimers affect the differentiation pathway choice, and one where they do not. The cell differentiation, where only functional dimers are involved, is modeled by reactions 18 to 22.

∅ kx -→ X, (18)
A c c e p t e d m a n u s c r i p t

X + D f,1 + Cell 0 k dif -→ Cell 1 (19)
X + D f,2 + Cell 0 k dif -→ Cell 2 (20)
X + D f,1 + D f,2 + Cell 0 0.5×k dif -→ Cell 3 (21)
X + Cell 0 0.5×k dif -→ Cell 4 (22)
In the model of cell differentiation, differentiation is assumed irreversible, and both the moments at which a cell commits to differentiate (reaction 18) and the choice of differentiation pathway (reactions 19 to 22) are stochastic. The cell commits to differentiate after a protein 'X' becomes present in the cell via reaction 18 (internally produced or external signal). Its rate, k x , is set to 10 -6 s -1 , implying that, on average, X will appear in the cell at approximately t = 0.5 × 10 6 s, which is half of a cell's lifetime (here set to 10 6 s). This unrealistically long lifetime provides better statistics, but one could instead follow the dynamics of cells over many generations (as done experimentally in [START_REF] Neubauerz | Immunity phase shift in defective lysogens: nonmutational hereditary change in early regulation of λphage[END_REF]). Reaction 18 occurs at most once during a cell's lifetime. Once it occurs and X is present in the cell, reactions 19 to 22 can then occur. The choice of differentiation pathway depends only on the amounts of dimers (D f,i or both D f,i and D nf,i ) present at the moment of differentiation (reactions 19 to 22), where k dif is set to 10 -5 s -1 ). These reactions require 'X' and 'cell 0 ' as substrates, thus are competing reactions and only one of them, at most, occurs in a cell, since there is only one molecule X. The differences between the rate constants of these reactions approximately compensate for the differences in propensities caused by the different amounts of dimers needed in different differentiation pathways.

cell 0 represents a stem cell and cell i (where i ∈ {1, . . . , 4}) are markers of different cell fates. Reactions 19 to 22 define how D 1 and D 2 , at the moment of differentiating, stochastically determine the cell's differentiation pathway. For example, when these reactions become active, the cell will most likely differentiate into cell type 1 (reaction 19) if only D 1 is present (it is possible to also differentiate into type 4 but with lower propensity). If both D 1 and D 2 are present in abundance, the cell probably differentiates into cell type 3 (reaction 21) and, if both dimers are absent, into cell type 4 (reaction 22).

Reactions 19 to 22 assume that only functional dimers are involved in determining cell differentiation pathway. However, one can not exclude the possibility that non-functional dimers, while not affecting the dynamics of TS, are functional from the point of view of affecting differentiation pathway

A c c e p t e d m a n u s c r i p t

choice. We model this case, assuming that functional and non-functional dimers are indistinguishable in what concerns the selection of cell differentiation pathway, by adding reactions 23 to 27 to reactions [START_REF] Hong | Shadow enhancers as a source of evolutionary novelty[END_REF][START_REF] Holmberg | Multisite phosphorylation provides sophisticated regulation of transcription factors[END_REF][START_REF] Huang | Cell fates as a high-dimensional attractor states of a complex gene regulatory network[END_REF][START_REF] Kaern | Stochasticity in gene expression: from theories to phenotypes[END_REF].

X + D nf,1 + Cell 0 k dif -→ Cell 1 (23) X + D nf,2 + Cell 0 k dif -→ Cell 2 (24)
X + D f,1 + D nf,2 + Cell 0 0.5×k dif -→ Cell 3 (25)
X + D nf,1 + D f,2 + Cell 0 0.5×k dif -→ Cell 3 (26) X + D nf,2 + D nf,2 + Cell 0 0.5×k dif -→ Cell 3 ( 27
)
We described four hypothetical models which we next study in detail. First, non-functional dimers can both bind to the genes and affect differentiation pathway choice ('case B D'). Second, non-functional dimers can bind to the genes but not affect the differentiation pathway choice ('case B N'). Third, non-functional dimers cannot bind to the genes but affect the differentiation pathway choice ('case N D'). Fourth, non-functional dimers neither bind to the genes nor affect differentiation pathway choice ('case N N'). Note that in all cases non-functional dimers do not repress gene expression.

Results

In each of the four cases we analyze the effects of varying the 'transformation rate', the 'non-functional fraction', and the 'sub-system biases' in the cell differentiation pattern. The results in each case are then confronted. Each cell population consists of 1000 cells to allow sufficient sampling.

We also measure the noise levels, computed as the ratio between standard deviation and mean level of the protein dimers' temporal levels [START_REF] Ribeiro | Noisy attractors and ergodic sets in models of gene regulatory networks[END_REF]. The dimers' total mean levels and standard deviations are calculated from the time series of both functional and non-functional dimers of each cell population (the time series of each population are conjuncted). We discard the initial transient (1000 s) since cells are initialized without proteins or RNA.

Effects of the transformation rate

The time needed for a premature protein to become functional (or nonfunctional) is determined by the 'transformation rate', k trans . In this section,

A c c e p t e d m a n u s c r i p t

we study the effects of varying this time length (which follows, for each protein, an exponential distribution) in the patterns of cell differentiation. We fix f i to 0.1 and assume no sub-system bias (i.e., f 1 = f 2 where both proteins are such that their propensity to be non-functional are identical given the cellular environment). We generate five cell populations differing in the value of k trans (reactions 9 and 10). Since k trans cannot be lower than 0.001 (below which most premature proteins decay before maturing) and since the decay of premature proteins could be neglected when k trans is high enough (k trans > 100), we set k trans = 0.001, 0.01, 0.1, 1, 100 in cells of population 1 to 5, respectively.

Varying k trans within these boundaries does not affect the differentiation patterns nor the stochasticity of the TS in any of the four cases, since k trans does not affect the expected fractions of functional and non-functional proteins present in the cell, but only the time for a premature protein to become functional or not. This case therefore allows us to determine how the dynamics of the TS differs when non-functional dimers are able, or not, to bind to the promoter, given that they do not repress it when bound.

The differentiation pattern for k trans = 1 in the four cases are shown in Fig. 1. The loss of the ability of non-functional dimers to bind to promoters causes significant differences in the patterns, namely, less cells differentiate into type 3 and instead differentiate to types 1 and 2 (equally likely) because the binding of non-functional dimers to promoters alters the TS dynamics in comparison with the non-binding case. Specifically, when non-functional dimers compete for promoter binding sites, they protect the target genes of repression, allowing both genes to express simultaneous, thus disrupting the bistability of the TS. When non-functional dimers are involved in cell fate determination, the probability of observing the co-existence of dimers from both sub-systems further increases. Thereby, by preserving either of these functionalities, non-functional dimers, although do not perform their repression function, can affect the dynamics of TS and thereby cell differentiation.

Fig. 2 plots cell population's noise level against different 'log 10' transformed transformation rate ('trans') for both cases B and N. It reveals that the noise level of 'case B ' is much lower than 'case N ' for all values of k trans for the reasons mentioned. When non-functional dimers bind to promoters (case B ) and protect them from repression, more proteins are accumulated while the proteins fluctuations are dampened due to the loss of toggling behavior.

Also, from Fig. 2, both noise levels remain almost constant as the transformation rate varies (except for a slight decrease at k trans = 0.001). A

A c c e p t e d m a n u s c r i p t

premature protein can either become functional, or non-functional, or be degraded. The relative propensity of these choices (conversion vs. decay) can be compared by computing R td , which is the ratio between the transformation rate and the decay rate of premature proteins, R td = ktrans dp+2ktrans . Since d p ≤ k trans , and both d p , k trans > 0, 1 3 ≤ R td ≤ 1 2 . When R td is high enough (R td > 10 21 ) so that the decay of premature proteins is negligible, the amounts of functional and non-functional proteins in the cell are stable over time after a transient, resulting in stable noise levels across different cell populations. However, when R td = 1 3 (the extreme scenario depicted at k trans = 0.001), where decay is an equally competitive pathway, the mean and standard deviation of functional and non-functional proteins are significantly decreased. These two variations 'compensate' one another so that the noise level of the TS and the differentiation patterns remain mostly invariant as R td varies.

Effects of the non-functional fraction

In this section we vary the non-functional fraction (fraction of premature proteins that become non-functional) and study its effect on the dynamics of the TS and on cell differentiation patterns. For this, we vary f i (reaction 10) equally for both genes of the TS (namely, f 1 = f 2 = f ), without altering the gene expression rates or the transformation rate (k trans = 1). As mentioned, that are several mechanisms by which cells can regulate the activity of multi-functional proteins [START_REF] Holmberg | Multisite phosphorylation provides sophisticated regulation of transcription factors[END_REF] and, as assumed here, the fraction of functional proteins, pertaining a given functionality.

We simulate six cell populations, with 'f' set to 0.0001, 0.001, 0.01, 0.1, 1, and 10, respectively, for cell populations 1 to 6. The case f = 0.0001, serves as a control model in which there are virtually no non-functional proteins, since there are less than 1000 proteins in total in a cell at any given time. This implies that on average, for f = 0.0001 there are no non-functional proteins in the cell or these exist in a quantity that does not significantly affect the cell's dynamics.

The differentiation pattern of each case is plotted against the 'log 10' transformation of 'f' in Fig. 3. More cells differentiate into type 3 in detriment of types 1 and 2 as the fraction of non-functional dimers increases, in all four cases, meaning that the intervals when both proteins are present at the same time during the cells lifetime increases with f .

A small amount of non-functional proteins (f < 0.0001) in the cells is sufficient to cause cells to opt more for types 1 and 2 instead of type 3 in

A c c e p t e d m a n u s c r i p t

the cases B N and B D, while the critical point where changes in differentiation pattern occur is at f = 0.1 and f = 1, respectively, in cases N D and N N. This agrees with the results that both functionalities render more cells to differentiate into type 3 but do so at different extents (see section 3.1). This also applies to the difference in the differentiation patterns of cases B D and B N. Since non-functional proteins compete with functional ones to form dimers, when most functional proteins are forming dimers with nonfunctional ones due to the overwhelming dominance of non-functional ones (f > 1), cells differentiate more to types 1 or 2 in case B N. However, the differentiation pattern remains unchanged in case B D since non-functional dimers are functional regarding differentiation. The noise measured in cells of 'case B ' and 'case N ' is shown in Fig. 4. The noise varies in case B but not in case N (for f within realistic intervals).

Effects of the sub-system bias

In this section, we study the effects of sub-system bias on cell differentiation. The sub-system bias, 'r = f 2 f 1 , is such that for r = 1 one of the proteins is more prone to be non-functional. We set r = 0.1 . . . 0.9 with a step size 0.2, and set k trans = 1 and f 1 = 0.1. Since r < 1 it implies that P 1 is more prone to be non-functional.

The fraction of cells differentiated into each cell type is plotted against r in Fig. 5. The sub-system bias does not affect the patterns when nonfunctional dimers can bind to promoters (cases B N and B D). Accordingly, the noise is invariant with sub-system bias in such cases (Fig. 6). This is again caused by the simultaneous expression of both proteins in these cases. Further notice that, the internal noise level of case B is much lower than that of case N, similarly to the case depicted in Fig. 2.

In the cases that r affects the differentiation patterns, the degree of the changes are proportional to the variation of r. Also, as seen from Fig. 6, the noise increases slightly in case N as r decreases (stronger bias). This is because as r increases, gene 1 is less subject to repression and thus, the TS regains bistability which, consequently, increases the noise level.

Finally, more cells differentiate into type 3 in case N D than in case N N. This is explained by the increased probability of the simultaneous presence of both dimers in case N D, as discussed in section 3.1.

A c c e p t e d m a n u s c r i p t

Discussion

We studied how the multi-functionalities of repressor proteins can be manipulated to regulate the dynamics of a toggle switch, leading to different patterns of cell differentiation. To this objective, we varied the transformation rate, the non-functional fraction and the sub-system bias. We observed that changing the transformation rate did not affect the patterns of cell differentiation, while altering the non-functional fraction or the sub-system bias did. Further, the magnitude of such regulation depends on the extent and which functionalities are preserved in partially non-functional dimers.

In our protein model, we assumed two regulatory protein domains, a DNA-binding domain and a repressor domain [START_REF] Monod | Teleonomic mechanisms in cellular metabolism, growth, and differentiation[END_REF]. Using this model, two functional cases were studied as non-functional protein dimers which can bind ('case B') or cannot bind ('case N') DNA. We found that, in both cases, nonfunctional protein affects the patterns of cell differentiation by interfering with the bistability of the TS. As well, we discovered that an increased proportion of nonfunctional protein resulted in a decreased in the noise level of the proteins' of the TS. These results have interesting implications for the involvement of nonfunctional dimers in cell differentiation.

Various extracellular signaling molecules or changes in the extracellular environment can lead to changes in gene expression. This occurs by mechanisms which target signaling pathways to the control of transcription factor activity. One of the more common methods of transcription factor regulation is by phosphorylation or dephosphorylation of the protein. This modification can regulate a transcription factor by manipulating its localization, its ability to interact with other regulatory proteins or DNA, or targeting it for proteolytic degradation [START_REF] Whitmarsh | Regulation of transcription factor function by phosphorylation[END_REF]. Importantly, transcription factors can have independent domain with specific functions, such as DNA-binding and proteinprotein interactions, and functional regulation of each domain can be directly localized without interfering with other domains. In this study, we examined, using computer simulations, the dynamics of regulating the functionality of a hypothetical multidomain protein which is involved in a TS. We found a interesting ability to vary the physiological outcome of cell differentiation by tuning this protein's functionality. It is therefore interesting to speculate that the cell has evolved similar mechanisms to vary protein functionality, using methods of activation such as phosphorylation, as a method of controlling gene expression and thereby differentiation.

Our results define a possible regulatory role of the control of multi-

A c c e p t e d m a n u s c r i p t

functionality in protein in small genetic circuits. They suggest that, in addition to transcriptional and translational levels of regulation, the activation of various functionalities in multi-functional proteins might be an important level of regulation in the dynamics of genetic networks. These findings can be tested, and prove valuable, in synthetic biology by engineering genetic circuits for specific purposes. Finally, one can, from these results, hypothesize that the mechanisms that regulate protein multi-functionality are subject to selection and thus, are likely to evolve. Also, while we assume the regulation of protein multifunctionality by regulatory molecules in the cell, such as in phosphorylationdephosphorylation, the probability that a proteins functionality may also be sequence dependent to some extent, thus, evolvable.

Figures Legends

Differentiation pattern of a cell population (with k trans = 1) for the cases 'BD', 'BN', 'ND' and 'NN'. x-label represent different cases, y-label is the number of cells of each cell type after differentiation. The number of cells of type 4 is not zero but is always below 1%.

Figure 2

 2 Figure 2Noise level in cell populations against different transformation rates (k trans ) for cases 'B ' and 'N '.

Figure 3

 3 Figure 3Differentiation patterns for various non-functional fractions (f = f 1 = f 2 ) for the cases 'BD', 'BN', 'ND' and 'NN'. x-label is the log 10 of f , y-label is the number of cell types after differentiation.

Figure 4

 4 Figure 4 Noise level in cell populations for various values of non-functional fractions (f = f 1 = f 2 ) for cases 'B ' and 'N '.

Figure 5

 5 Figure 5Differentiation patterns for various sub-system biases (f 2 /f 1 ) for the cases 'BD', 'BN', 'ND' and 'NN'. x-label is the bias towards sub-system 2, and

Figure 6

 6 Figure 6Noise level in cell populations for various sub-system biases (f 2 /f 1 ) for cases 'B ' and 'N '.