Functional characteristics of a double positive feedback loop coupled with autorepression
aa r X i v : . [ q - b i o . M N ] F e b Functional characteristics of a double positive feedbackloop coupled with autorepression
Subhasis Banerjee and Indrani BoseOctober 28, 2018
Department of Physics, Bose Institute, 93/1, A. P. C Road, Kolkata-700009, India
Abstract
We study the functional characteristics of a two-gene motif consisting of a double positive feed-back loop and an autoregulatory negative feedback loop. The motif appears in the gene regulatorynetwork controlling the functional activity of pancreatic β -cells. The model exhibits bistability andhysteresis in appropriate parameter regions. The two stable steady states correspond to low (OFFstate) and high (ON state) protein levels respectively. Using a deterministic approach, we show thatthe region of bistability increases in extent when the copy number of one of the genes is reduced fromtwo to one. The negative feedback loop has the effect of reducing the size of the bistable region. Lossof a gene copy, brought about by mutations, hampers the normal functioning of the β -cells giving riseto the genetic disorder, maturity-onset diabetes of the young (MODY). The diabetic phenotype makesits appearance when a sizable fraction of the β -cells is in the OFF state. Using stochastic simulationtechniques, we show that, on reduction of the gene copy number, there is a transition from the monos-table ON to the ON state in the bistable region of the parameter space. Fluctuations in the proteinlevels, arising due to the stochastic nature of gene expression, can give rise to transitions between theON and OFF states. We show that as the strength of autorepression increases, the ON → OFF statetransitions become less probable whereas the reverse transitions are more probable. The implicationsof the results in the context of the occurrence of MODY are pointed out.. .A.C.S. Nos.: 87.18Cf, 87.18Tt, 87.18Vf Introduction
Positive and negative feedback loops are frequently-occurring motifs in gene transcription regulatorynetworks and signaling pathways [1, 2]. The components of a feedback loop are genes, proteins andother molecules which are connected by regulatory interactions. Depending on the components and theirinteractions, feedback loops have distinct roles in diverse regulatory systems. A regulatory interactionis positive (negative) if an increase in the amount or activity of one component increases (decreases)the amount or activity of its interaction partner. A feedback loop is positive (negative) if the numberof repressing interactions is zero or even (odd). A large number of experiments and theoretical studieselucidate the major functional characteristics of feedback loops with simple structure [1, 2, 3, 4, 5, 6,7, 8]. Positive feedback in a gene transcription regulatory network (GTRN) tends to enhance proteinlevels whereas negative feedback favours homeostasis ,i.e., maintenance of proteins at a desired level.The simplest feedback loop has only one component which is thus self-regulating. For such a motif ina GTRN, a protein promotes / represses its own production via autoactivation / autorepression of theexpression of its gene. A positive feedback loop with two components and two regulatory interactionsis of two types: double negative and double positive. Again, considering a GTRN, the protein productsof the two genes in a double negative feedback loop repress each other’s synthesis. The construction ofa synthetic circuit, the genetic toggle, is based on this motif [9]. The double positive feedback loop isdefined by two genes, the protein products of which promote each other’s synthesis. There are severalexamples of two-component positive feedback loops in natural cellular networks [1, 2], a prominentexample being the cell division cycle, the regulatory network of which contains both double positiveand double negative feedback loops [10]. In this case, the loops control enzymatic activity. The doublenegative feedback loop, because of its more common occurrence, has been extensively studied in contrastto the double positive feedback loop.The next stage of complexity in feedback loops involves linked positive and negative feedback loops[2, 11, 12, 13]. The key variables in the dynamics of feedback loops are the concentrations of the com-ponent molecules. In the case of a GTRN, these may be the protein concentrations. In a deterministicdescription, the time evolution of the concentrations is determined by solving a set of coupled differentialequations, the number of equations being equal to the number of variables. In reality, the biochemical3 ig1.eps
Figure 1: The two-gene network model. The protein products of the GP and GS genes activate eachother’s synthesis. There is also an autorepressive loop in which the proteins of the GS gene repress theirown synthesis.events associated with gene expression are probabilistic in nature and this is reflected in the presence offluctuations (termed noise) around mean protein levels [14]. A stochastic description of time evolutionis thus more appropriate. A single positive feedback loop has a tendency to amplify noise, also the timetaken to reach the steady state protein level is longer than that in the case of an unregulated gene [1, 2].Interlinking of two positive feedback loops with slow and fast dynamics results in a switch with rapid ac-tivation and slow deactivation times and a marked resistance to noise in the upstream signaling pathway[11]. Addition of a single negative feedback loop leads to rapid deactivation in the absence of the signalwhich activates the switch [12]. The combination of positive and negative feedback loops may give riseto excitability with transient activation of protein levels. Recent experiments suggest that competencedevelopment in B. subtilis is achieved via excitability [15].In this paper, we study the functional characteristics of a two-gene double positive feedback loopcoupled with autorepression of the expression of one of the genes. The major motivation for studying thisspecific motif is its presence in the GTRN controlling the pancreatic β -cell function [16]. The hormoneinsulin is a small protein that is synthesized in the β -cells and secreted when an increase in the bloodglucose level is sensed. Glucose metabolism releases energy needed by cells to do useful work. Insulinis necessary to metabolize glucose and thereby control its level in the blood. Diabetes occurs due to anexcessive accumulation of glucose in the blood brought about by an insufficient production of or reducedsensitivity to insulin. The core of the β -cell transcriptional network consists of a double positive feedbackloop in which the transcription factors HN F − α and HN F − α , belonging to the nuclear hepatocytefamily, activate each other’s synthesis. There is also some evidence that HN F − α autorepresses itsown synthesis [16]. Mutations in the transcription factors HN F − α and HN F − α give rise to atype of diabetes known as maturity-onset diabetes of the young (MODY) which has an early onset withage less than usually years. There are six different forms of MODY of which MODY 1 and MODY4 are caused by mutations in the genes hnf − α and hnf − α respectively [17]. The structure of theregulatory network, of which the two genes hnf − α and hnf − α are integral components, is not fullyknown. A partial structure of the complex network is shown in [16, 17] involving the genes hnf − α,hnf − α, shp, hnf − β, hnf − β, hnf − γ, hnf − γ, and pdx − . The genes collectively controlthe transcription of a number of important genes involved in glucose metabolism in the β − cell. Theseinclude the glucose transporter ( Glut − gene, the glucokinase gene encoding the glycolytic enzymeglucokinase which acts as glucose sensor and also the insulin gene. Odom et al. [18] combined chromatinimmunoprecipitation assays with promoter microarrays to gain insight on the regulatory circuits formedby hnf − α and hnf − α. Both the proteins are found to control the activity of a large number of targetgenes in the β − cell. This recent finding as well as earlier experiments [16] indicate that the hnf − α and hnf − α genes play a prominent role in the pancreatic β − cell function. Mutations in the genes giverise to MODY resulting in the impairment of glucose-stimulated insulin secretion. Several experiments[16] provide clues on the possible molecular origins of MODY. The cross-regulatory interactions between HN F − α and HN F − α are switched on as pancreatic β -cells receive the signals to differentiate. Thedouble positive feedback loop has the potential for bistability, i.e., two stable steady states. The two statesare a basal state in which the two genes have low activity and an activated state which corresponds to highprotein levels. The states are analogous to the OFF and ON states of a switch. Normal functioning of thepancreatic β -cells requires the two-gene feedback loop to be in the ON state. The circuit operation is,however, vulnerable to decreased gene dosage caused by mutations (in a diploid organism each gene hastwo identical copies). Genetic disorders, termed haploinsufficiency, are known to occur due to reducedgene dosage resulting in decreased protein levels [19, 20, 21, 22]. Gene expression noise increases theprobability that a protein level falls below a threshold value so that the protein amount is insufficientfor meaningful activity. The loss of vital protein functions is responsible for the occurrence of geneticdisorders. MODY, brought about by reduced gene dosage, is thus an example of haploinsufficiency [16].We construct a mathematical model to study the dynamics of the core circuit consisting of a doublepositive feedback loop coupled with autorepression of the hnf − α gene. We use both deterministic andstochastic approaches to identify the functional features of the motif and discuss their possible relevancein the occurrence of MODY 5 ig2.eps Figure 2: The reaction kinetic scheme of the two-gene model. The meanings of the symbols are explainedin the text.
The circuit diagram of the motif to be studied is shown in the figure 1. GP and GS represent the genes hnf − α and hnf − α respectively. The arrow sign denotes activation by the appropriate proteinproduct and the hammerhead sign denotes repression. The chemical kinetic schemes corresponding tothe expression of genes GP and GS are shown in figures 2(a) and 2(b). The protein products of GP and GS are denoted by P and S . We assume that the regulation of gene expression is mediated by theprotein dimers S and P , K P and K S being the binding constants of dimerization. For each gene, thereare two rates of protein synthesis: a basal rate ( rate constants J P and J S ) and an activated rate ( rateconstants J P and J S ). In the second case, protein synthesis occurs in the activated state of the gene ( GP ∗ and GS ∗ ) attained via the binding of protein dimers S and P to the genes GP and GS respectively.The associated binding constants are K P P and K SS . The rate constants for protein degradation are γ P and γ s with φ denoting the degradation product. Dimer degradation is not taken into account as its rateis few-fold lower than the degradation rate of protein monomers. For the gene GS , there is an extrabiochemical event representing autorepression. The dimers S and P bind the promoter region of thegene GS competitively, i.e., the binding of one type of dimer excludes the binding of the other type.When the dimer S binds at GS , there is complete repression. The binding constant is denoted by K R .The protein concentrations S and P are the dynamical variables in the system. The time scale ofbinding events, in general, is much faster than that of protein synthesis and degradation. The boundcomplexes thus reach the steady state at an earlier time point. Taking this into account, the differentialrate equations describing the time evolution of the protein concentrations S and P are: dSdt = J S GS ∗ + J S GS − γ S S (1) dPdt = J P GP ∗ + J P GP − γ P P (2)6ith GS ∗ = n S T (cid:16) PM (cid:17) T (cid:16) PM (cid:17) , M = 1 K SS K P , T = 11 + K R K S S (3) GP ∗ = n P (cid:16) SN (cid:17) (cid:16) SN (cid:17) , N = 1 K P P K S (4)There are two conservation equations for the total concentrations n S and n P of the genes GS and GP . n S = GS + GS ∗ + GSS (5) n P = GP + GP ∗ (6)After an appropriate change in variables u = SJ S /γ S , v = PJ P /γ P , τ = γ S t (7)the differential rate equations (1) and (2) are transformed into dudτ = n S η β v (1 + µ u ) + β v − u (8) dvdτ = n P ξ α u α u − v (9)The different parameters are given by η = J S J S , ξ = J P J P , µ = J S γ S ! K S K R , α = J S γ S ! K P P K S , β = J P γ P ! K SS K P (10)The variable τ is dimensionless whereas the variables u and v have the dimensions of concentrationexpressed in units of [nm]. The parameters η and ξ are dimensionless while the parameters µ , α and β are expressed in units of nm ] . The dimensions of n s and n p are in units of [ nm ] with one gene copy7 ig3a.eps fig3b.eps Figure 3: u versus η curves showing bistability and hysteresis. The solid (dashed) lines represent stable(unstable) steady states for gene copy numbers (a) n P = 2 , n S = 2 and (b) n P = 2 , n S = 1 . Theparameter µ , a measure of the autorepression strength, is zero.corresponding to approximately nm ] . From now on, the units will not be explicitly mentioned. Table 1displays all the parameters and rescaled parameters of the two-gene model as well as their meanings anddefining formulae.We use the software package XPPAUT [23] to probe the dynamics of the double positive feedbackloop and the effect of autorepression of the S proteins on the dynamics . We focus on how the steadystate value of u (rescaled concentration of S proteins) changes as a function of the different parametersin equations (8) and (9). In the steady state, the rates of change dudτ and dvdτ are zero. Figure 3(a) shows aplot of u versus η when the autorepression strength given by µ is zero. The other parameters have values n S = n P = 2 , ξ = 30 . , J S = J P = 2 . , γ S = γ P = 1 . and α = β = 0 . . The plot showsthat a region of bistability separates two region of monostability. The two stable states in the bistableregion correspond to low and and high values of u . In this region and at a specific value of η , the choicebetween the stable steady states is history-dependent, i.e., depends on initial conditions [24]. If the valueof η is initially low, the system ends up in the low u state. As η increases, the system enters the regionof bistability but continues to be in the low expression state till a bifurcation point is reached. At thispoint, a discontinuous jump to the high u state occurs and the system becomes monostable. Bistability isaccompanied by hysteresis , i.e., the value of η at which the switch from the low to the high expressionstate occurs is greater than the value of η ( the lower bifurcation point ) at which the reverse transitiontakes place. The two stable branches are separated by a branch of unstable steady states (dash-dottedline) which are not experimentally accessible. There are now several known systems in which bistabilityand hysteresis have been observed experimentally [9, 13, 25, 26, 27, 28, 29]. Figure 3(b) shows the plotof u versus η for the same parameter values as in figure 3(a) except that the copy number of the GS geneis reduced from two to one, i.e., n S has the value . A comparison of figures 3(a) and (b) shows thatwith reduced copy number the extent of the region of bistability in considerably increased. The same8 arameter/Rescaled Parameter Meaning/Defining Formula J P , J S Rate constants for basal protein synthesis J P , J S Rate constants for activated protein synthesis K P , K S Binding constants of protein dimers S and P K P P , K SS Binding constants for the binding of proteindimers S and P at the genes GP and GSγ P , γ S Rate constants for protein degradation K R Binding constant for repressor dimer bindingat gene GS ; K R thus denotes the strength ofrepression η η = J S J S , ratio of activated and basal rateconstants for synthesis of S proteins ξ ξ = J P J P , ratio of activated and basal rateconstants for synthesis of P proteins µ µ = ( J S γ S ) K S K R ; with J S , γ S and K S keptfixed, µ can be varied by changing K R thusproviding a measure of repression strength α α = ( J S γ S ) K P P K S β β = ( J P γ P ) K SS K P T, M, N abbreviations defined in equations (3) and (4)Table 1: Parameters, rescaled parameters, their meanings and defining formulaeconclusion is reached when the steady state values of u are plotted versus the parameter β. The region ofbistability is lower in extent when the parameter µ, a measure of the autorepression strength, is increasedfrom zero. The value of µ is changed by modifying the value of K R (equation (10)), the binding constantfor repressor binding at the GS gene. Figure 4 shows the phase portrait corresponding to equations (8)and (9) with the parameter values ξ = 30 , η = 30 , α = 0 . , β = 0 . and µ = 0 . Thesystem is bistable for the parameter values quoted. The nullclines, obtained by putting dudτ = 0 , dvdτ = 0 , intersect at three points, the fixed points of the dynamics. The lower and upper fixed points are stablewhereas the intermediate fixed point is unstable, in fact, a saddle node [30]. The stable manifold of thesaddle node divides the uv − phase space into two basins of attraction. Trajectories starting in the lower(upper) basin of attraction end up at the lower (upper) stable fixed point as shown in figure 4. A trajectoryinitiated on the stable manifold stays on it and ends at the saddle node. A typical trajectory asymptoticallyapproaches the unstable manifold as t → ∞ . A trajectory is obtained by plotting the values of u and v at different time points, determined by solving equations (8) and (9). The arrow direction on a trajectorydenotes increasing time. 9 ig4.eps Figure 4: Phase portrait described by equations (8) and (9). The dark and light solid lines represent thenullclines intersecting at three fixed points. The stable and unstable fixed points are denoted by solidand empty circles respectively. The stable manifold divides the phase space into two basins of attraction.Some typical trajectories are shown with arrow directions denoting increasing time.Figures 5(a) shows the plot of ξ versus η exhibiting regions of monostability and bistability. Theparameter values are the same as before with α = β = 0 . and µ = 0 . The regions of bistability,enclosed within the red and black curves, correspond to n S = 1 and n S = 2 respectively. The differencein the locations of the two loops in the logarithmic plots clearly shows that the bistable region is of greaterextent when the gene copy number is reduced from two to one. The region of bistability is decreased inextent when autorepression is taken into account (Figure 5(b) with µ = 0 . . Figure 6 shows the µ − β plot with the regions of bistability falling within the red ( n S = 1 ) and black ( n S = 2 ) curves respectively.The value of µ is changed by varying K R (equation 10) with µ = 0 . K R .A major advantage of combining a double positive feedback loop operating between two genes withautorepression of the expression of one of the genes lies in dosage compensation [16]. This relates to thefact that the fall in steady state protein levels, brought about by a reduction in the gene copy number, isless when autorepression is included, compared to the case when there is no autorepression. A measureof dosage compensation is provided by the quantity G , termed percentage gain, defined as G ( µ ) = x ( µ ) − x ( µ = 0) x ( µ = 0) × (11)where x denotes the steady state concentration of S proteins when the copy number of the GS gene, n S , is one. The parameter µ is a measure of the repression strength. G is calculated by keeping themean level of S proteins to be the same in the two cases µ = 0 and µ = 0 when n S = 2 . This isachieved by adjusting the binding constant K P P contained in the parameter α in equation (10). Theother parameter values are η = ξ = 30 . , J S = J P = 2 . , γ S = γ P = 1 . and β = 0 . .10 ig5a.eps fig5b.eps Figure 5: Plots of ξ versus η showing regions of monostability and bistability when the parameter µ, ameasure of the autorepression strength, is zero (a) and 0.005 (b). The regions of bistability are enclosedwithin the red and black curves with gene copy numbers n P = 2 , n S = 1 and n P = 2 , n S = 2 respectively.Figure 7 shows the plot of G versus µ for the parameter values mentioned. As µ increases from zero,there is initially a sharp increase in the value of G followed by a slower growth which ultimately leadsto a near-saturation of G values. The results obtained in the deterministic approach provide insighton the advantages of autorepression in the non-occurrence of the genetic disorder MODY. The normalfunctioning of pancreatic β -cells requires the HN F − α and HN F − α protein levels to be high, i.e.,the two-gene system should be in the ON state. The genesis of MODY lies in a substantial fraction of the β -cells being in the OFF state. This is brought about by mutations in the hnf − α and hnf − α genesgiving rise to a fall in the steady state protein levels. In terms of the two-gene model studied by us, themonostable high state, in which the levels of the P and S proteins are both high, represents the ON stateof normal β -cells. The system may enter a region of bistability, in which both the ON and OFF states arepossible, due to the loss of a gene copy brought about by mutations. We will show in the next sectionthat fluctuations in the protein levels are responsible for transitions between the ON and OFF states. Inthe deterministic scenario, the major advantages of the autorepressive feedback loop appear to be dosagecompensation (figure 7) as well as a lesser possibility of the system being in the bistable region due to areduction in gene copy number. The continuance of the system in the monostable high state ensures thenormal functioning of cells. Similar conclusions are reached if the gene copy number n P is reduced fromtwo to one. There is, however, an asymmetry in the S and P protein levels as the expression of the gene GP is not autorepressed. 11 ig6.eps Figure 6: µ − β phase diagram with the regions of bistability falling within the red ( n P = 2 , n S = 1 )and black ( n P = 2 , n S = 2 ) curves respectively. fig7.eps Figure 7: Plot of percentage gain G (equation 11) versus µ , a measure of the autorepression strength. Consider the two-gene network to be originally in the monostable high state. In the deterministic for-malism, the system continues to be in the high, i.e., ON state even if it enters the region of bistabilitydue to the loss of a gene copy. This is due to history dependence, since the system is initially in theON state it continues to be in the ON state in the bistable region. The protein levels corresponding tothe ON state are, however, lower in magnitude in the bistable region. In the pancreatic β -cells, the oc-currence of MODY is possible only when a sizable fraction of cells is in the OFF state. The ON → OFFand OFF → ON transitions can be understood only when stochasticity in gene expression is taken intoaccount. We now give a simple physical picture of the origin of stochastic transitions [6]. In the caseunder consideration, the dynamical variables are the protein concentrations u and v . In the case of deter-ministic time evolution, trajectories starting in individual basins of attraction stay confined to the specificbasins with no possibility of a trajectory crossing from one basin to another. In the stochastic approach,the trajectories are no longer deterministic as the dynamical variables u ( t ) and v ( t ) are fluctuating. Inthe deterministic case, given the initial state defined by ( u ( t = 0) , v ( t = 0)) , the trajectory is fixed in the uv -phase space. In the stochastic case, different trajectories are generated in repeated trials. A transientfluctuation, if sufficiently strong, switches the system dynamics from one basin to the other brought aboutby the excursion of the trajectory across the boundary separating the two basins of attraction. In terms of12 ig8a.eps fig8b.epsfig8c.eps fig8d.eps Figure 8: Distribution of steady state GS protein levels, P ( u ) , in an ensemble of cells for repressorstrengths (a) µ = 0 . , (b) µ = 0 . , (c) µ = 0 . and (d) µ = 0 . respectively.the pancreatic β -cells, the switch to the OFF state hampers the normal functioning of the cells.For proper regulatory functions as transcription factors, the HN F − α and HN F − α proteinlevels are high with an optimal value as excessive protein amounts are known to be harmful rather thanbeneficial [16]. In this context, it is pertinent to undertake a comparison of the functional characteristicsof two-gene network models with and without the autorepressive loop and with the mean protein levelskept at the same high values in the two cases. The last condition ensures the normal functioning of thecells in both the cases. In section 2, we have identified certain advantages of the autorepressive loopas regards the system dynamics in a deterministic framework. Our goal is now to identify the desirablefeatures of the model incorporating both a double positive feedback loop and an autorepressive looptaking the stochastic aspects of the dynamics into consideration. This is done with the help of a detailedcomputer simulation based on the Gillespie algorithm [31]. The algorithm enables one to keep trackof the stochastic time evolution of the system. The different biochemical reactions considered in thesimulation are depicted in figures 2(a) and 2(b). The reactions are sixteen in number and are given by GS + P → GS ∗ (12) GS ∗ → S (13)13 S → S (14) S + S → S (15) GS ∗ → GS + P (16) S → S + S (17) S → Φ (18) GP + S → GP ∗ (19) GP ∗ → P (20) GP → P (21) P + P → P (22) GP ∗ → GP + S (23) P → P + P (24) P → Φ (25) GS + S → GSS (26) GSS → GS + S (27)The different symbols are as explained in section 2. The stochastic rate constants, associated with thereactions, are C(i), i=1,...,16, in the appropriate units. The results of the simulation are shown in figures8-9. Figures 8(a)-(d) show the distribution of GS protein levels, P ( u ) , in an ensemble of 4500 cellsfor repressor strengths µ = 0 . , . , . and . respectively after a simulation time of tmax = 2000 time units. The gene copy numbers are n P = 2 and n S = 1 so that the system is in theregion of bistability. The values of the stochastic rate constants are C (2) = 56 . , C (3) = 2 . , C (4) = 4 . ,C (5) = 280 . , C (6) = 100 . , C (7) = 1 . , C (8) = 10 . , C (9) = 50 . , C (10) = 2 . , C (11) =4 . , C (12) = 280 . , C (13) = 100 . , C (14) = 1 . , C (15) = 10 . . The value of µ is changed by varying14he stochastic rate constant C (16) . The value of the rate constant C (1) is changed to keep the meanprotein levels to be the same when n S = 2 , n P = 2 for all values of µ . The value µ = 0 implies thatonly the double positive feedback loop contributes to the dynamics. The distribution P ( u ) is found to bebimodal, i.e., has two distinct peaks corresponding to the OFF and ON states. In all the cases, the cellsare in the ON state at time t = 0 . One finds that the fraction of cells in the OFF state decreases as thevalue of µ increases. In fact, when µ = 0 , the number of cells which are in the OFF state is larger thanthat in the ON state. Since initially all the cells are in the ON state, a large number of ON → OFF statetransitions occur during the simulation time. For µ = 0 , the reverse transition is, however, much rarer.The role of the autorepressive loop thus appears to be to reduce the number of stochastic transitions fromthe ON to the OFF state. This makes the occurrence of MODY, brought about by a sizeable fractionof the cell population existing in the OFF state, less probable. There are two distinct time scales overwhich protein fluctuations occur. The probability distribution P ( u ) versus u has a two-peaked structure.Fluctuations on a short time scale confine the u values to lie predominantly within individual peaks.The long time scale corresponds to the time at which large fluctuations occur bringing about transitionsbetween states belonging to different peaks. The “escape time” is often very large and a quantitativemeasure is provided by the mean first passage time τ [32]. In the present case, the values of τ ON → OF F and τ OF F → ON are quite large for different values of µ. The maximum simulation time tmax is timeunits for all values of µ. For µ = 0 , τ ON → OF F is around time units whereas τ OF F → ON is even larger.As µ increases, τ ON → OF F increases whereas τ OF F → ON decreases. For µ = 0 . , τ ON → OF F is as largeas s. Because of large escape times, the probability distribution P ( u ) versus u is metastable on a largetime scale [32]. Over shorter periods of time, the shape of the distribution remains more or less invariant.The plots in figure 9 are obtained for an ensemble of 4500 cells. For gene copy numbers n P = 2 and n S = 2 , the mean protein levels are adjusted to be the same irrespective of the values of µ . Theparameter values are so chosen that the system is in the monostable high region. On reduction of n S to (one copy of the GS gene), the system enters the region of bistability and is in the ON state at time t = 0 . After a period T = 2000 time units of stochastic time evolution, the percentage of cells in theOFF state is determined. The red curve shows this percentage as a function of the repression strength µ .The drop in the percentage of cells in the OFF state is found to be exponential. The black curve shows15 ig9.eps Figure 9: For gene copy numbers n P = 2 , n S = 1 and after a time interval T = 2000 time unitsof stochastic time evolution, the percentage of cells in an ensemble of cells in the OFF state (redcurve) versus the repression strength µ with all the cells being in the ON state at time t = 0 . The blackcurve shows the percentage of cells in the ON state versus µ with all the cells being in the OFF state at t = 0 .the percentage of cells in the ON state after T = 2000 time units, with all the cells being initially in theOFF state. One finds that with increasing µ , the fraction of cells in the ON state becomes larger. Theautorepressive loop has the effect of making the ON state more stable and the OFF state more unstable.This feature enhances the probability of the nonoccurrence of MODY as there are infrequent transitionsfrom the ON to the OFF state. On the other hand, the system has a lesser probability of remaining stuckin the OFF state compared to the case when there is no autorepressive loop. In this paper, we have studied the functional characteristics of a motif consisting of a double positivefeedback loop operating between two genes and a negative feedback loop in which the protein productof one gene represses its own synthesis. The motif appears in the gene regulatory network controllingthe pancreatic β -cell function [16]. Loss of a gene copy due to mutations has been shown [16] to beresponsible for abnormal β -cell function resulting in MODY. We have studied the effect of reduced genecopy number on the dynamics of the model describing the two-gene motif. In a deterministic formalismbased on differential rate equations, we identified regions of bistability in appropriate parameter regions.The stable steady states, designated as the OFF and ON states, correspond respectively to low and highprotein levels. The normal β -cells are expected to be in the monostable ON state. The occurrence ofMODY is brought about by a fraction of β -cells being in the OFF state. The ON → OFF switch can occuronly in the bistable region. Negative feedback reduces the extent of the bistable region making it less16ikely that the cellular state falls in this part of the phase diagram. The region of bistability, however,increases in size on reduction of the gene copy number making the ON → OFF transitions more probable.Negative feedback also produces a mechanism of dosage compensation (figure 7). The results hold truefor a wide range of parameter values. Since switching to the OFF state is detrimental, one would havethought, from an evolutionary point of view, that two genes which are constitutively ON would be moreappropriate. In reality, the genes hnf − α and hnf − α form a positive feedback loop. Cross-regulationbetween the two genes is established when the pancreatic cells receive signals to differentiate [16]. Thepositive feedback loop provides a stable mechanism of gene expression since the two genes reset eachother’s activity to the functional state under physiological perturbations. This serves to self-perpetuatethe activity of the two genes and their targets in the pancreatic β − cells. Normal functioning of these cellsrequires both the protein levels to be high. The system of two genes that are constitutively ON are lessrobust under physiological perturbations since there is no resetting mechanism by which both the genesare in the functional ON state. The theoretical suggestions of bistability due to the existence of a positivefeedback loop [16, 18], backed up by the results of our mathematical model, should be tested in actualexperiments.The ON → OFF switch is brought about by protein fluctuations the origin of which lies in stochasticgene expression. Our major finding is that negative feedback makes the ON → OFF transitions less prob-able and the OFF → ON transitions more probable. Thus the function of the negative feedback appearsto be to protect the normal β -cell function since the cell is more likely to be in the ON state in this case.The asymmetric response to fluctuations prevents switching off and facilitates switching on of the highexpression state. In the deterministic scenario, one finds that the difference between the ON state andthe unstable steady state protein level increases as the autorepression strength is increased whereas thedifference between the unstable steady state and OFF state protein levels decreases on increasing the au-torepression strength. This may explain the asymmetry in the ON → OFF and OFF → ON switches whenstochasticity is taken into account. For moderate strengths of autorepression, the system is locked in theON state for extremely long times. In our simulations, we did not encounter ON → OFF switches for verylong trajectories ( ∼ seconds) with µ = 0 . . This translates into lifetimes measured in years andexplains the delayed onset of the diabetic phenotype [16]. The phenotype generally appears after several17ears indicating that the activation of the ON → OFF switch is rare. The average age at which MODY ismanifest could thus be dictated by the probability that a sufficient number of β -cells is locked in the OFFstate. We have considered the simplest form of negative autoregulation in our two-gene model. Thereare recent suggestion that negative autoregulation of the HNF-4 α gene in the pancreatic β -cells may bemore complex [33]. Also, the number of transcription factor binding sites of the two genes is not knownwith certainty. Cooperative binding at multiple sites is expected to promote the stability of the gene ex-pression states. Our two-gene motif constitutes a minimal model which seeks to explains the desirablefeatures of combining a double positive feedback loop with an autorepressive loop vis-á-vis the normalfunctional activity of β -cells. The insight gained from the model study is expected to provide a basis forthe investigation of more complex cases. 18 eferences [1] Alon U 2007 Network motifs: theory and experimental approaches Nat. Rev Genet. [9] Gardner T S, Cantor C R and Collins J J 2000 Construction of a genetic toggle switch in Escherichiacoli Nature β -cells: Implications for differentiation and haploin-sufficiency . 5
346 - 51[27] Ozbudak E M , Thattai M, Lim H N, Shraiman B I, Oudenaarden van A 2004 Multistability in thelactose utilization network of Escherichia coli Nature Vol. α gene expression by HNF-4 α GP and GS genes activate eachother’s synthesis. There is also an autorepressive loop in which the proteins of the GS gene repress theirown synthesis.Fig2. The reaction kinetic scheme of the two-gene model. The meanings of the symbols are explainedin the text.Fig3. u versus η curves showing bistability and hysteresis. The solid (dashed) lines represent stable(unstable) steady states for gene copy numbers (a) n P = 2 , n S = 2 and (b) n P = 2 , n S = 1 . Theparameter µ , a measure of the autorepression strength, is zero.Fig4. Phase portrait described by equations (8) and (9). The dark and light solid lines represent thenullclines intersecting at three fixed points. The stable and unstable fixed points are denoted by solidand empty circles respectively. The stable manifold divides the phase space into two basins of attraction.Some typical trajectories are shown with arrow directions denoting increasing time.Fig5. Plots of ξ versus η showing regions of monostability and bistability when the parameter µ, ameasure of the autorepression strength, is zero (a) and 0.005 (b). The regions of bistability are enclosedwithin the red and black curves with gene copy numbers n P = 2 , n S = 1 and n P = 2 , n S = 2 respectively.Fig6. µ − β phase diagram with the regions of bistability falling within the red ( n P = 2 , n S = 1 )and black ( n P = 2 , n S = 2 ) curves respectively.Fig7. Plot of percentage gain G (equation 11) versus µ , a measure of the autorepression strength.Fig8. Distribution of steady state GS protein levels, P ( u ) , in an ensemble of cells for repressorstrengths (a) µ = 0 . , (b) µ = 0 . , (c) µ = 0 . and (d) µ = 0 . respectively.Fig9. For gene copy numbers n P = 2 , n S = 1 and after a time interval T = 2000 time units ofstochastic time evolution, the percentage of cells in an ensemble of cells in the OFF state (redcurve) versus the repression strength µ with all the cells being in the ON state at time t = 0 . The blackcurve shows the percentage of cells in the ON state versus µ with all the cells being in the OFF state at t = 0 . 22 P GS S P+P P2K PP K P S0 Φ K S K SS S
J * S2J P *GPGP + S2 Φγ P J P0 S+S GSGS + P2 JGS + S2 GSS2K R u η u η
10 20 30 40 50 60 7001020304050
Unstable manifold u v Stable manifold η ξ η ξ β µ .000 0.005 0.010 0.015 0.020246810121416 G
10 20 30 4005001000150020002500 P ( u ) u
10 20 30 4005001000150020002500 P ( u ) u
10 20 30 4005001000150020002500 P ( u ) u
10 20 30 4005001000150020002500 P ( u ) u µ % o f c e llll