A mathematical model of the sterol regulatory element binding protein 2 cholesterol biosynthesis pathway

(cid:1) We formulate and analyse a nonlinear ODE model of the SREBP2 pathway. (cid:1) The mathematical model exhibits stable limit cycles under certain parameter conditions. (cid:1) Negative feedbacks in the SREBP2 pathway may help regulate cholesterol homeostasis. (cid:1) Our model provides a more accurate formulation of genetic regulation using nonlinear ODEs. Cholesterol is one of the key constituents for maintaining the cellular membrane and thus the integrity of the cell itself. In contrast high levels of cholesterol in the blood are known to be a major risk factor in the development of cardiovascular disease. We formulate a deterministic nonlinear ordinary differential equation model of the sterol regulatory element binding protein 2 (SREBP-2) cholesterol genetic regulatory pathway in a hepatocyte. The mathematical model includes a description of genetic transcription by SREBP-2 which is subsequently translated to mRNA leading to the formation of 3-hydroxy-3-methylglutaryl coenzyme A reductase (HMGCR), a main regulator of cholesterol synthesis. Cholesterol synthesis subsequently leads to the regulation of SREBP-2 via a negative feedback formulation. Parameterised with data from the literature, the model is used to understand how SREBP-2 transcription and regulation affects cellular cholesterol concentration. Model stability analysis shows that the only positive steady-state of the system exhibits purely oscillatory, damped oscillatory or monotic behaviour under certain parameter conditions. In light of our ﬁ ndings we postulate how cholesterol homeostasis is maintained within the cell and the advantages of our model formulation are discussed with respect to other models of genetic regulation within the literature.


Introduction and motivation
As an essential constituent of the plasma membrane of mammalian cells, cholesterol is used for the maintenance of both membrane structural integrity and selective permeability (Simons and Iknonen, 2000). However, superfluous cholesterol levels result in cellular toxicity (Yeagle, 1991;Tabas, 1997;Tangirala et al., 1994). Insufficient cholesterol causes cytotoxicity via compromised membrane structure.
Furthermore cellular cholesterol metabolism is a key modulator of plasma cholesterol, with the management of plasma hypercholesterolaemia at the cornerstone of population cardiovascular disease management (Grundy et al., 2004). It is therefore crucial that intracellular cholesterol levels are strictly regulated. Cellular cholesterol homeostasis, the property to maintain cholesterol concentration to within narrow ranges, results from a balance of three mechanisms: efflux, influx and biosynthesis.
Understanding the mechanisms which regulate cellular cholesterol content is vital to understanding pathology associated with sub-and supra-optimal cell and blood cholesterol concentrations. These levels are dependent on both the balance between dietary cholesterol intake and de novo synthesis of cholesterol within cells.
The low density lipoprotein receptor (LDLR) protein forms part of the lipoprotein metabolic pathway responsible for the clearance of cholesterol from the circulation (Brown and Goldstein, 1979;Goldstein et al., 1985). Biosynthesis of cholesterol is a multistep reaction in which the rate-limiting step is the reduction of 3hydroxy-3-methylglutaryl coenzyme A (HMG-CoA) in the reaction catalysed by the enzyme HMG-CoA reductase (HMGCR).
Over accumulation or excessive depletion of free cholesterol within the cell is prevented by a negative feedback loop that responds to elevations or depressions in intracellular cholesterol. This feedback loop exerts the majority of its control by regulating the synthesis of the two key proteins: HMGCR and LDLR. In brief, when the intracellular cholesterol level is low, both LDLR and HMGCR synthesis are activated, thereby increasing the influx of cholesterol via the LDLR pathway, and the biosynthesis of cholesterol in the cell. If conversely there are high cholesterol levels in the cell, synthesis of LDLR and HMGCR declines.
There has been much research conducted into the response of cell cholesterol to dietary intake, with the dietary fatty acid composition rather than cholesterol intake reported to have a greater impact on circulating cholesterol concentrations. In particular, partial replacement of saturated fat with either monounsaturated (found in olive oil) or n-6 polyunsaturated (found in vegetable oils such as sunflower oil) fatty acids have been associated with significant reductions in both total and LDLcholesterol concentrations (Mensink et al., 2003;Micha and Mozaffarian, 2010). Dietary fat composition is considered to influence circulating cholesterol concentrations via effects on hepatic cholesterol synthesis and the expression of genes involved in circulating LDL-cholesterol metabolism (Xu et al., 1999).
Previous mathematical modelling has included compartmental models of the lipoprotein metabolic pathway (Knoblauch et al., 2000;Packard et al., 2000;Adiels et al., 2005) and dynamic models of lipoprotein metabolism in conjunction with the LDLR pathway (August et al., 2007;Wattis et al., 2008). Of particular note in these dynamic models is the lack of explicit representation of the cholesterol biosynthesis reaction and as a consequence, the interplay between cholesterol biosynthesis, the LDLR uptake of lipoprotein cholesterol and cholesterol mediated negative feedback is not fully appreciated.
The cholesterol biosynthetic pathway is already the basis of the most common form of pharmaceutical treatment for high plasma cholesterol levels. HMGCR inhibitors, more commonly known as statins, act as competitive inhibitors of the HMGCR enzyme. By inhibiting the biosynthesis of cholesterol, statins deplete intracellular cholesterol concentration and promote the synthesis of both HMGCR and the LDLR, thereby increasing the uptake of lipoproteins (and plasma cholesterol) via the LDLR. It is recognised that individual response to statin treatment varies widely. Genetic variation in HMGCR has been associated with a diminished lipid lowering response (Chasman et al., 2004;Krauss et al., 2008), suggesting that the cholesterol biosynthetic pathway plays an important role in the control of plasma cholesterol levels.
However, relatively little modelling has been conducted to investigate the qualitative behaviour of the processes which govern de novo cholesterol synthesis at a genetic level, which may provide a better understanding of such phenomena. The mathematical model presented in this paper will examine the underlying genetic mechanisms governing cholesterol biosynthesis as a first step towards elucidating the dynamics of this pathway.
The paper is organised as follows. In Section 2 the biological processes which describe the genetic regulation of cholesterol biosynthesis are reviewed. Following this, the mathematical model is derived in Section 3 and details of model parameter values obtained from the literature are summarised in Section 4. Model analysis is undertaken in Sections 5-7 and the results are summarised and discussed in Section 8.

Regulated expression of cholesterol biosynthetic genes
A major point of control of the cholesterol biosynthetic pathway occurs at the level of gene expression in response to cellular cholesterol levels, as shown in Fig. 1. The insolubility of cholesterol dictates that it cannot directly influence a genetic response. The critical role in controlling the expression of a range of genes involved in the regulation of cellular lipid homeostasis falls to the three isoforms of the SREBP family of transcription factors, SREBP-1a, SREBP-1c and SREBP-2. In particular, the SREBP-2 isoform is relatively specific to regulating the expression of many enzymes involved in cholesterol biosynthesis (Brown and Goldstein, 1997).
SREBPs exist normally in a tight complex with the SREBP cleavage activating protein (SCAP) within the endoplasmic reticulum of cells. SCAP consists of two domains, one of which is responsible for the association with SREBP. The other domain contains a region known as the sterol sensing domain (SSD). When the cellular cholesterol concentration becomes depleted, SCAP escorts SREBP to the Golgi apparatus of the cell, where it undergoes sequential cleavage by proteases. The net effect of this is to liberate the transcription factor, nuclear SREBP which can then enter the cell nucleus (Eberlé et al., 2004). Here it binds to a regulatory binding site (a short sequence of DNA) on the promoter region of the target gene known as the sterol regulatory element (SRE) and activates its transcription (Soutar and Knight, 1990).
In the presence of replete cellular sterol concentrations, cholesterol binds directly to the SSD of SCAP. This causes a conformational change in SCAP which results eventually in the anchoring of the SCAP-SREBP complex to the endoplasmic reticulum (ER) membrane (Yang et al., 2002). This process is responsible for the retention of the SCAP-SREBP complex within the ER. Transcription of the target genes declines.
In the context of the HMGCR gene, when a cell 0 s cholesterol levels are low, the SCAP-SREBP complex is active and free to move. In such a state SREBP is formed and is able to reach the nucleus and activate HMGCR mRNA transcription and thus HMGCR synthesis, increasing the cholesterol concentration in the cell by upregulating its synthesis. If, conversely, there are high cellular cholesterol levels, then SCAP-SREBP is unable to move and effectively inactive. Consequently both HMGCR mRNA transcription and HMGCR translation decrease, and cholesterol synthesis is reduced.
In a simplified model of the gene expression response to cellular cholesterol concentration, the system can be seen as an end product negative feedback loop system, in the manner of the mathematical models of expression developed by, for example, Goodwin (1963Goodwin ( , 1965 and Griffith (1968). In such models, the response of the gene is directly dependent upon the concentration of cholesterol. A very low level of cholesterol will provoke a large response in the synthesis of HMGCR enzyme, and vice-a-versa. Theoretically, this results in a considerable range over which the model allows cholesterol concentration to vary. This is, however, uncharacteristic of the homeostatic property which the physiological system possesses, and which ensures that cellular cholesterol can only fluctuate within a narrow range of values, to avoid the cytotoxicity associated with extreme values.
The addition of the SREBP transcription factor function models the underlying biological mechanism, and also introduces complexity to the negative feedback loop in the form of an activator function which is suppressed by accumulation of an end product. In the following section a model of this interaction between SREBP and cholesterol, and the effect on gene expression are presented.

The model
The interactions characterising cellular cholesterol homeostasis and its regulation by transcription factors are many, and a full model of all variables and reactions is not necessarily pragmatic. Furthermore, the number of parameter values required will increase with complexity. Previous models have shown that excessive simplification can fail to reproduce dynamics which have been observed in experimental settings.
As an example, the work by Wattis et al. (2008) models nonlipoprotein cholesterol influx to the cell as proportional to the difference between cell cholesterol concentration and a predetermined ideal equilibrium value; this produces the correct dynamics for cell cholesterol response. An interesting consequence of this formalism, though, is that intracellular cholesterol concentration in the model reaches equilibrium rapidly (on a timescale of the order of minutes) after an influx of lipoprotein cholesterol to the cell. However, experimental results suggest that this may not be the case, with changes in intracellular cholesterol concentration occurring on timescales of 12-24 h (Liscum and Faust, 1987;Liscum et al., 1989). This suggests that not enough complexity is included here to capture the longer term dynamics of cholesterol synthesis at the level of the HMGCR gene.
A further requirement is that the system must respond naturally in the absence or presence of cholesterol as opposed to only acting reasonably under certain circumstances. For example, in the work of August et al. (2007), all cholesterol in the cell is assumed to be derived from lipoprotein sources. Whilst this reproduces the required qualitative behaviour under the conditions whereby extracellular lipoprotein is present, in the case where this is zero, the intracellular cholesterol level falls to zero, which is physiologically fatal for the cell.
The work presented in this paper is focused on formulating and analysing a nonlinear ordinary differential equation (ODE) model of the SREBP-2 cholesterol biosynthesis pathway. The goal of the work is to understand cholesterol regulation via the negative feedback between SREBP-2 transcription and cholesterol and to what extent this affects the steady-state cholesterol levels of the cell. In doing so we hope to more accurately capture cellular regulation of cholesterol and be able to understand it in the wider context of dietary cholesterol intake.

Mathematical model formulation
In this section we derive a system of nonlinear ODEs to describe the genetic regulation of cholesterol biosynthesis by SREBP-2 as summarised in Fig. 2.
The binding of SREBP-2 to the gene, subsequent transcription and translation to HMG-CoA mRNA and production of HMGCR and Fig. 2. The genetic regulation of cholesterol production by SREBP-2. The HMGCR gene G is transcribed at a rate μ m to produce HMGCR mRNA M. This is translated at a rate μ h into the HMGCR enzyme H. HMGCR then goes on to catalyse the reaction creating the metabolite cholesterol C at a rate μ c . This process is under the control of the transcription factor SREBP S which acts as a transcriptional activator for the pathway. Under conditions where cholesterol C is in excess S forms an inactive complex with C and transcription of the target gene declines. HMGCR mRNA, HMGCR and cholesterol are degraded at rates δm, δ h and δc, respectively. cholesterol can be described by the reaction equation Here x is the number of molecules of S required to bind to G to produce a functional effect. This binding reaction has an association rate κ 1 and a dissociation rate κ À 1 . M is transcribed at a rate μ d and H is translated at a rate μ h . The creation of C occurs at a rate μ c . δ m , δ h and δ c are respectively the degradation rates of M, H and C.
Similarly the binding of cholesterol to active SREBP-2 to form an inactive complex which down-regulates the transcription of cholesterol (negative feedback) is given by where y is the number of molecules of C required to bind to S to cause inactivation. This binding reaction has an association rate κ 2 and a dissociation rate κ À 2 .
We note two important biological concepts arising from the physiological mechanism of gene expression or protein synthesis, which will affect the form of the ODEs (Alberts et al., 2008) describing Eqs. (1) and (2).
(i) ½G : xS represents the concentration of DNA in an active state, which is able to undergo transcription. During transcription, activated DNA is copied by the action of an enzyme to produce mRNA. This process does not deplete ½G : xS. (ii) Protein is synthesised from mRNA via the action of ribosomes.
Following protein synthesis, mRNA detaches from the ribosome and the mRNA is free to participate in further synthesis reactions until it is degraded according to its half-life. Therefore, the synthesis of the enzyme, H, does not affect the concentration of M. That is, synthesis of H will not deplete M.
The governing ODEs equations are derived by application of the law of mass action to the biochemical reactions (1) and (2) which gives dg dt with initial conditions where in the above system of equations, we use the following notation in which square brackets denote concentration: g ¼½G, The coefficient x in the first term of Eq. (4) reflects that the dissociation of one active DNA complex releases x molecules of unbound transcription factor. The coefficient x in the second term of Eq. (4) states that the creation of one active DNA complex requires up to x DNA binding sites.
The number of genes within a cell is constant so adding Eqs. (3) and (5) leads to on using the initial conditions (10). We now assume that Eq. (5) reaches equilibrium rapidly (quasi-steady-state approximation) such that and using Eq. (11) we have which upon rearranging gives Here κ d is the dissociation constant for the reaction between S and

G.
We further observe that adding Eqs. (4), (5) and (9) gives Under the quasi-steady state assumption of Eq. (12) together with the initial conditions (10) we find that Also under the approximation (12) we see that s b C s b ð0Þ C 0. This is a valid assumption if we consider that the concentration of binding sites for a particular transcription factor on one particular gene is extremely small compared to the concentration of free transcription factor available in the cell, i.e. s b oos. We then obtain the following equation from (19): Finally we assume that the binding reaction between S and C reaches equilibrium rapidly such that κ 2 c y s Àκ À 2 ðs 0 À sÞ C0: Rearranging this result gives in which we define the constant κ c such that where κ s is the dissociation constant for the reaction between S and C. Using Eqs. (14), (20) and (22) to eliminate Eqs.
Non-dimensionalisation: Before proceeding to a complete analysis of the model, Eqs. (24)-(26) are non-dimensionalised. Time is scaled with respect to the synthesis rate of m such that where τ represents the non-dimensional time. The remaining variables are rescaled with respect to the concentration of total transcription factor, s 0 , such that This non-dimensionalisation leads to with the initial conditions where the non-dimensional parameters are given by The non-dimensional parameter values are summarised in Table 2.

Parameter estimation
A summary of the model parameter values is provided in Table 1 with details on how each was derived from the experimental literature given in Appendix A. Wherever possible data elicited from human liver cells (Hep G2) have been used. However, it has not been possible to determine all required parameters in this manner. In some cases the model parameters do not have a direct physiological counterpart since the biological processes occurring have been simplified in the mathematical modelling to reduce complexity; in others, the parameter value is not customarily measured in the required units, not least because of the difficulty in isolating the biosynthesis pathway. In these instances underlying biological principles have been used to estimate a realistic value, and to ensure that the model operates within a plausible physiologic domain.

Model analysis
In this section and continuing in Sections 6 and 7 we discuss the existence of steady-states of Eqs (30)-(32) and their stability.

Fixed point analysis
The steady states of equations (30)-(32) are given by the solution of where m ss , h ss and c ss are the steady state values of m, h and c respectively. Using Eqs. (36) and (37), Eq. (35) can be written as Expanding, we find that the steady states are given by the solution of the polynomial equation of degree ðxy þ1Þ, As all parameters are positive, we may apply the results of Descartes 0 Rule of Signs which states that the number of positive real roots of the system is either equal to the number of variations of signs in the coefficients of Eq. (40) or less than this by an even integer (Murray, 2002). As there is only one sign change in the sequence of coefficients Eq. (40), the system has only one positive real root, and therefore only one physiologically valid fixed point.

Fixed point stability
We now consider the stability of this fixed point by investigation of the eigenvalues of the linearised Jacobian matrix J of the system equations (30)-(32). The Jacobian is given by We note that ψ Z0 as all parameter values are positive and that c ss Z0 for physiologically valid parameter ranges. Calculation of the eigenvalues of J requires the solution of the equation where λ are the eigenvalues to be found and I is the identity matrix. Evaluation of Eq. (43) leads to the characteristic equation, which has three roots, the eigenvalues λ 1 , λ 2 and λ 3 . Firstly we note that ψ Z 0 ensures all coefficients of Eq. (44) are positive and thus by Descartes' Rule of Signs there can be no purely positive real eigenvalues. There are then two cases for the roots of (44), either three negative real eigenvalues or one negative real eigenvalue and a pair of complex conjugate eigenvalues.
The fixed point is stable if and only if the real parts of λ 1 , λ 2 and λ 3 are negative. To determine for which conditions this occurs, we apply the Routh-Hurwitz Stability criteria to Eq. (44) (Murray, 2002). Routh-Hurwitz 0 s criteria applied to a cubic equation of the form λ 3 þa 2 λ 2 þ a 1 λþa 0 ¼ 0 are satisfied if and only if a 0 4 0, a 1 40, a 2 4 0 and a 1 a 2 À a 0 4 0. That is, the necessary and sufficient condition for the roots of Eq. (44) to have negative real part requires Since all parameters are positive and real, conditions (45)- (47) hold. Thus the stability of the roots is dependent on condition (48). The possible dynamic behaviour of the system can be summarised as follows.
Case I: ρðδ m ; δ h ; δ c Þ 4 0. Here all real parts of all eigenvalues are negative. In this case the steady state is stable. This stable steady state may arise in one of two ways: (i) Case Ia: where all eigenvalues are real and negative. This will result in a stable node, where the concentrations of mRNA, protein and cholesterol will tend monotonically to a steady state; and (ii) Case Ib: where one eigenvalue is real and negative and two eigenvalues are complex conjugates with negative real part. In this case the fixed point is a stable spiral and the concentrations of mRNA, protein and cholesterol will demonstrate oscillatory convergence to a steady state.
Case II: ρðδ m ; δ h ; δ c Þ¼0. By substituting this value of ðδ m δ h δ c þ μ c ψÞ into Eq. (44), we now have the characteristic equation given by λ 3 þ γ 1 λ 2 þ γ 2 λþγ 1 γ 2 ¼ 0; ðλ þγ 1 Þðλ 2 þγ 2 Þ¼0; where we have γ 2 ¼ðδ m δ h þ δ m δ c þ δ h δ c Þ and γ 1 ¼ðδ m þ δ h þ δ c Þ. Therefore the characteristic equation has two conjugate roots λ 1;2 on the imaginary axis and one negative real eigenvalue λ 3 given by one negative real eigenvalue and two pure imaginary eigenvalues. The existence of two conjugate eigenvalues on the imaginary axis means that the stability of the equilibrium cannot be directly determined; this case is discussed in detail in Section 6.1.
Case III: ρðδ m ; δ h ; δ c Þ o0. Here one eigenvalue is real and positive and two eigenvalues are complex conjugates with positive real part. In this case the steady state is unstable, implying that end product concentration would grow unboundedly. This case is biologically infeasible and hence we ignore it for the remainder of this paper.
6. Fixed point stabilityvariation of μ c The eigenvalues of Eq. (44) can move between each case under the variation of system parameters. As an example we consider the effect of varying μ c on the system dynamics. This parameter may be varied so that a pair of complex conjugate eigenvalues can either move into the negative real half plane (a stable focus equilibrium) or into the positive real half plane (an unstable focus equilibrium). The point where the eigenvalues cross the imaginary axis (Case II) occurs at a critical value of μ c denoted by μ n c . At this point a unique, closed periodic orbit may bifurcate locally from the equilibrium as it changes stability. The isolated, closed trajectory is known as a limit cycle and causes oscillatory behaviour. This phenomenon is called a Hopf bifurcation (Guckenheimer and Holmes, 1983) and its existence dictates that the concentrations of m, h and c will oscillate.

Hopf bifurcation existence
According to the Hopf bifurcation theorem (Guckenheimer and Holmes, 1983), a bifurcation occurs for a critical value μ c ¼ μ n c , for which the following two conditions are fulfilled, at the equilibrium point ðm ss ; h ss ; c ss Þ: 1. The matrix J has two complex eigenvalues in some neighbourhood of μ n c and for μ c ¼ μ n c these eigenvalues are purely imaginary, that is, θðμ n c Þ¼0 and ωðμ n c Þ a 0: This non-hyperbolicity condition is a necessary condition for the Hopf bifurcation. 2. The relation described by holds in some neighbourhood of μ n c . This is a sufficient condition for the Hopf bifurcation and is also known as the transversality or Hopf crossing condition. It ensures that the eigenvalues cross the imaginary axis with non-zero speed and thus ensures that the crossing of the complex conjugate pair at the imaginary axis is not tangent to the imaginary axis. If this is not the case we may observe, for example, the occurrence in which the eigenvalues move up to the imaginary axis and then reverse direction without crossing.
We notice that the first condition has already been shown to hold at the critical value μ n c , given by the solution of where ψ is given by Eq. (42), together with the equation determining the equilibrium value of c ss for μ n c : From the results of Case II we know that at this value of μ n c the characteristic polynomial Eq. (44) has two purely imaginary roots where this represents an implicit function of two variables μ c and λ. The Implicit Function Theorem tells us that we may define μ c as a function of λ near the point ðμ n c ; λðμ n c ÞÞ, and the derivative of this function is given by providing ∂k ∂λ a0: We begin by computing the derivative of the function where, we have previously noted that ψ Z 0. However, in the case ψ¼0, the Jacobian matrix becomes which is lower triangular and hence has three negative real eigenvalues given by the entries of the leading diagonal, specifically Àδ m ; Àδ h and À δ c . This violates the requirement of condition 1 that the matrix J has two complex eigenvalues. Therefore we can conclude that in the case of a Hopf bifurcation, ψ a 0 and we need only be concerned with the strict inequality ψ 4 0.
Eqs. (54) and (55) and since ψ 4 0, and the second condition of the Hopf theorem is fulfilled. Thus we have proved the existence of a Hopf bifurcation at the critical value μ c ¼ μ n c .

Hopf bifurcation stability
Just as the steady states of a system may be stable or unstable, the limit cycle which branches from the fixed point in a Hopf bifurcation may also be stable or unstable. A stable limit cycle occurs if the Hopf bifurcation is supercritical whereas an unstable limit cycle is the product of a subcritical bifurcation.
At a subcritical bifurcation a unique and unstable limit cycle, which exists for μ c o μ Â c , is absorbed by a stable spiral equilibrium.
The equilibrium becomes unstable for μ c 4μ Â c ; in this case diverging oscillations and therefore unbounded growth in the evolution of variables is seen. In a supercritical bifurcation the equilibrium point prior to the Hopf bifurcation is a stable spiral, and concentrations of mRNA, protein and cholesterol display oscillatory decay before reaching a steady state value. At the bifurcation point μ c ¼ μ Â c , a limit cycle is born. At this point the equilibrium changes stability and becomes unstable. For μ c 4 μ Â c , this becomes a unique and stable small amplitude limit cycle; here the concentrations of mRNA, protein and cholesterol exhibit stable oscillations.
As the limit cycle is stable, any small perturbation from the closed trajectory causes the system to return to the limit cycle resulting in self sustained oscillations in concentrations of mRNA, protein and cholesterol within the region of some equilibrium value. Thus, as the occurrence of a supercritical Hopf bifurcation will result in behaviour which is analogous to the physiological process of homeostasis, it is necessary to determine the stability of the Hopf bifurcation. This analysis was undertaken as follows.
Numerical solutions to Eqs. (30)-(32) were obtained using MATLAB (MATLAB, 8.0.0.783, The MathWorks Inc., Natick, MA, 2012) and the characteristics of system bifurcations and limit cycles were explored using the MATLAB numerical continuation toolbox Matcont (Dhooge et al., 2003). The basic principle of this toolbox is to consider a system of ODEs with an equilibrium point at ðx 0 ; μ 0 Þ. The principle of numerical continuation requires finding a solution curve s of f ðx 0 ; μ 0 Þ¼0 with sð0Þ¼ðx 0 ; μ 0 Þ which describes how the equilibrium point varies. The curve s is traced by means of a predictor-corrector algorithm and bifurcations along s are detected using a suitable test function which changes sign at the bifurcation point.
Once the Hopf bifurcation has been detected, Matcont calculates the stability of the subsequent limit cycle by calculating the first Lyapunov coefficient or Lyapunov characteristic exponent, l 1 ð0Þ, of the dynamical system near the bifurcation point. This coefficient describes the average rate at which neighbouring trajectories in the phase space converge or diverge. Specifically, l 1 ð0Þ o 0 implies that the system is attracted to a stable periodic orbit and l 1 ð0Þ 4 0 implies that the system is attracted to an unstable periodic orbit.
In the case of Eqs. (30)-(32) with μ c as the bifurcation parameter, we find that a Hopf bifurcation is predicted to occur at the point (μ Â c , c n ) ¼(1.809, 0.011), whose values are the solution of the simultaneous equations (40) and (48). This bifurcation has a negative first Lyapunov coefficient which indicates that a stable limit cycle is produced and the bifurcation is supercritical.
The results of the Hopf bifurication existence and stability analysis of the governing system of Eqs. (30)-(32) can now be discussed in the context of cellular cholesterol homeostasis. Homeostasis is the tendency of a system to regulate its internal environment by maintaining a stable condition. All homeostatic mechanisms use feedback inhibition to facilitate a constant level. Essentially this involves controlling the concentration of a particular variable within a narrow range in the region of an optimal value. If this concentration alters, the feedback inhibition pathway automatically initiates a corrective mechanism which reverses this change and brings it back towards an equilibrium. In a system controlled by feedback inhibition, the equilibrium is never perfectly maintained, but constantly oscillates about an optimal level. Thus the existence of the Hopf bifurcation and the consequent appearance of small amplitude oscillations in the concentrations of mRNA, protein and cholesterol, are significant in its similarity to the behaviour of a homeostatic system.

Illustration of system behaviour
In this section we present numerical solutions to Eqs. (30)-(32) using the MATLAB stiff differential equation solver ODE15s (MATLAB, 8.0.0.783, The MathWorks Inc., Natick, MA, 2012) for various values of μ c to illustrate the system behaviour elucidated in the previous sections. All remaining parameters were held constant as detailed in Table 1. The parameter μ c was varied between 1.53 Â 10 À 2 s À 1 and 6.46 Â 10 À 2 s À 1 (physiologically valid limits) to demonstrate the variation of system behaviour through Cases I to II.
Simulation of Eqs. (30)-(32) starting with the initial value of 1.53 Â 10 À 2 s À 1 shows monotonic non-oscillatory convergence to a steady state, equivalent to Case Ia as illustrated in Fig. 3.
Continually increasing μ c shows the system shifting to Case Ib.
Thus we see oscillatory convergence to a steady-state as shown in   Table 2 except nondimensional μ c ¼0.462. Nondimensional initial conditions: m(0)¼ 3.65 Â 10 À 8 , h(0)¼ 1.10 Â 10 À 5 and c(0) ¼2.30 Â 10 À 2 . Note that the evolution of HMGCR has been rescaled to allow for easier comparison. Fig. 4. Stable spiral equilibrium (corresponding to Case Ib) where there is one negative real eigenvalue and a pair of complex conjugate eigenvalues with negative real part; variable concentrations exhibit oscillatory convergence towards a steady state value. Initial conditions are as in Fig. 3. All parameters are as in Table 2 except for μ c which has been increased 2 fold to μ c ¼ 0.923.
Following the bifurcation, the evolution of the concentrations of mRNA, HMGCR and cholesterol are purely periodic, with small amplitude stable oscillations. The period of the oscillations in Fig. 5 is approximately 16.9 h; further numerical investigations have revealed that on manipulation of system parameters, the oscillatory period can vary between approximately 12 and 24 h.
We also find that after μ c has passed through its critical Hopf bifurcation value, no further changes in dynamical behaviour occur. That is, once the system becomes oscillatory, it remains in this manner for all μ c 4 μ Â c .

Remaining parameter analysis and system behaviour
Further numerical investigation of the governing system of equations has shown that each of the system parameters, if varied whilst all other parameters are kept constant, are capable of inducing a Hopf bifurcation. In the case of synthesis rates, μ m and μ c , only one Hopf bifurcation occurs and is supercritical. Any oscillatory behaviour the system demonstrates is always stable and these oscillations persist for any parameter value greater than its critical bifurcation value.
We note, however, that if either the degradation rates, ðδ m ; δ h ; δ c Þ, or binding affinities ðκ m ; κ c Þ, are taken to be bifurcation parameters, qualitatively different behaviour from the case discussed above is seen. Calculation of the critical values for these parameters indicates that there are two physiologically valid points where a Hopf bifurcation may occur.
Examining the case of δ c we see that the critical value δ n c for which a Hopf bifurcation may occur is given by the solution of the equation which is quadratic in δ n c and hence, for the case of two positive real roots, gives rise to the possibility that there are two Hopf bifurcation points associated with this parameter. This result in turn affects the steady-states of c ss which are determined from The eigenvalues at this critical point are given by and so the first Hopf bifurcation condition holds. Proceeding in the manner of the calculation for μ n and therefore the second condition of the Hopf theorem also holds.
Here, the unique equilibrium value undergoes two Hopf bifurcation points. Before the first of these points, the equilibrium is a stable spiral. At the first bifurcation point a supercritical Hopf bifurcation leads to the appearance of a stable periodic orbit (as the eigenvalues of the system cross the imaginary axis from left to right). The amplitude of this limit cycle increases initially as δ c increases whilst later decreasing until the second bifurcation point where the limit cycle disappears (as the eigenvalues of the system cross the imaginary axis from right to left) and the equilibrium point becomes stable again. Numerical analysis demonstrates a negative first Lyapunov coefficient for Hopf bifurcations confirming their supercriticality. For any value of δ c falling between the two bifurcation values purely oscillatory behaviour is generated, whilst outside this region only stable non-oscillatory solutions exist as illustrated in Fig. 6.

Discussion
We have formulated and solved a deterministic ODE model of cholesterol biosynthetic regulation by SREBP-2 in a hepatocyte. The model predicts the existence of oscillatory behaviour within the system which we believe is important in understanding cholesterol homeostasis. In the HMGCR system, such a mechanism means that perturbations may be made to certain system variables without losing the required concentration within which cholesterol is allowed to vary to guard against cytotoxicity. Other advantages to this dynamic mechanism include limiting the time during which cholesterol concentration is necessarily elevated Fig. 6. The response of mRNA concentration to variation of δ c . A clearly demarcated region of purely stable oscillatory behaviour is visible between stable steady state solutions. All parameters are as in Table 2 except for δ c which is successively varied with nondimensional values as indicated. Initial conditions are as in Fig. 3.   Fig. 5. Following the occurrence of a Hopf bifurcation the (now unstable) equilibrium is attracted to a stable limit cycle. Variable concentrations exhibit purely oscillatory behaviour; the oscillations are stable. Initial conditions are as in Fig. 3. All parameters are as in Table 2 except for μ c which has been increased approximately 4 fold to μ c ¼1.946.
within the cell in response to increased demand. Furthermore, controlling cellular cholesterol levels in this manner may incur less demand on cellular energy supplies than sustained elevation. Dynamic oscillatory steady-state behaviour allows the system to vary between upper and lower bounds consistent with an oscillatory homeostatic mechanism (Ghosh and Chance, 1964;Waxman et al., 1966).
Following the work of Goodwin (1965) and Griffith (1968) negative feedback regulation of mRNA levels, which are modulated by end product concentration, are often modelled using a Hill type function such that the where m ¼ mðtÞ is the mRNA concentration, b ¼ bðtÞ is the concentration of the end product, μ is the rate of mRNA transcription, K is the Hill constant, n is the Hill coefficient and α is the rate of mRNA degradation. Goodwin (1965) and Griffith (1968) showed that such a system will exhibit oscillations should n Z8. Values of n greater than approximately 4 are, however, deemed biologically implausible. In comparison, our mathematical model formulation has explicitly accounted for the interaction between not only the transcriber (in this case SREBP-2), but the negative regulation of the transcriber by the end product (cholesterol). The form of Eq.
(24) accurately accounts for these interactions and allows biologically realistic values for them to be included in the model formulation whilst accurately accounting for the system dynamics. Whilst our mathematical model has been formulated in the context of cholesterol biosynthesis this specific process of transcription factor regulated gene expression could be applicable to other pathways regulated by SREBP-2, in addition to the modulation of other lipid regulating proteins by the SREBP-1a and SREBP-1c isoforms. Further experimental research is necessary to evaluate these mathematical results and to clarify the system behaviour. However, this model and its analysis may serve as a basis for the investigation of transcription factor mediated gene expression dynamics, and furthermore constitute an important component of the synthetic engineering of biological circuits (Zhang and Jiang, 2010).

Acknowledgments
BSB acknowledges the support of a Biotechnology and Biological Sciences Research Council (BBSRC) UK CASE studentship in collaboration with Unilever Research. MJT is grateful for the support of a Research Council UK Fellowship during the period in which this work was undertaken.

A.1. Determination of synthesis parameters
Information on the rates of transcription and translation are rarely quantified in terms of mass per unit time, instead these rates are often measured relative to a control process. Therefore in order to estimate practical values for these parameters, we consider the relevant biological mechanisms.

A.1.1. Rate of HMGCR mRNA transcriptionμ m
In order to produce an estimate for the transcription rate the assumption that any time delays involved in the initiation of transcription and promoter clearance are negligible, is made. It is also assumed that no abortive transcripts are produced. Liver cells are somatic and therefore the majority are diploid, meaning they contain two chromosomes and thus normal amounts of DNA. Ignoring the existence of both tetraploidy and double nuclei that can be present within some liver cells, we assume all liver cells to be diploid.
We estimate the rate of transcription in terms of nucleotides per unit time in a typical mammalian gene. It is possible for a cell to transcribe 14,000 base pairs in 20 min giving a rate of 12 bases per second (Darzacq et al., 2007). This value is used to provide a rough estimate of the rate of transcription, equivalent to the number of mRNA molecules produced per cell per unit time.
The human HMG-CoA reductase gene is 24,826 bases long (Goldstein and Brown, 1984). To transcribe one molecule of primary HMG-CoA reductase mRNA, from one gene, assuming a rate of 12 bases per second, takes 2068.83 s. We then take into account that the post transcriptional processing steps of mRNA cleavage, 5 0 capping and polyadenylation are rate limiting. A delay of approximately 30 min has been added to account for this. Therefore an approximation of the time it takes to transcribe one molecule of primary HMG-CoA reductase mRNA from one gene is 3868.83 s. Per gene, this equates to 2:59 Â 10 À 4 molecules HMGCR mRNA s À 1 : ðA:1Þ Therefore one gene can synthesise 2.59 Â 10 À 4 molecules of HMG-CoA reductase mRNA per second. Taking diploidy into account there are 5.17 Â 10 À 4 molecules HMG-CoA reductase mRNA being synthesised per cell per second.

A.1.2. Rate of HMGCR protein translationμ h
To calculate an estimate of the rate of translation the following assumptions are made. Firstly, any effects caused by the transport of mRNA from the nucleus to its localisation in the cytoplasm are ignored. Also ignored are the effects of protein folding on transcriptional regulation as well as biochemical interactions amongst proteins. Finally any time delays in the elongation phase of protein synthesis are considered to be negligible.
The in vivo estimate from Slobin (1991) suggests that the translation rate for eukaryotic cells is 6 amino acids per second, where one amino acid is encoded by 3 nucleotides or bases. Many ribosomes can translate the same mRNA molecule simultaneously (Granner and Weil, 2006). Because of their large size, ribosomes cannot attach to mRNA any closer than 35 nucleotides apart. This detail allows estimation of the rate of transcription, equivalent to proteins synthesised per unit time from mRNA.
A human HMG-CoA reductase mRNA transcript contains 4475 bases (Goldstein and Brown, 1984 that is, there are 3.32 Â 10 À 2 molecules of HMG-CoA reductase protein being synthesised per second. We thus take μ h ¼ 3:32 Â 10 À 2 s À 1 as an approximation of the rate of HMGR translation.

A.1.3. HMGCR activityμ c
The value of 52 pmol mevalonate formed per minute per mg protein for HMGCR activity in human liver microsomes is taken directly from Tanaka et al. (1982). This must now be converted into a turnover number for HMGCR which is the maximum number of substrate molecules that an enzyme can convert to product per mole of catalytic site of the enzyme per unit time. The activity of an enzyme is the moles of substrate converted per unit time and is a measure of the quantity of active enzyme present. The specific activity is the activity of an enzyme per mg of total protein, i.e. a measure of enzyme efficiency.
HMG-CoA reductase is a tetrameric protein, composed of monomers arranged in two dimers. Each dimer has two active sites (Istvan et al., 2000) and has a molecular mass ¼99,906 Da. (The Dalton, (Da), (alternatively atomic mass unit) is a unit of mass used to express atomic and molecular masses. It is defined so that one mole of a substance with atomic or molecular mass 1 Da will have a mass of precisely 1 g.) Therefore the full HMG-CoA reductase molecule is an enzyme containing 4 active sites with a molecular mass M r ¼2 Â 99,906¼ 199,812 Da, i.e. 1 mole HMG-CoA reductase has a mass of 199,812 g. The reaction catalysed by HMG-CoA reductase is given by which indicates a stoichiometry of one mole mevalonate being produced from one mole HMG-CoA. (Note all other substrates of this reaction are assumed to be present in excess). Segel (1993) suggests that there are approximately 1000 different enzymes in a cell. Thus the moles of enzyme in 1 mg of protein are equivalent to 1 Â 10 À 3 g 199; 812 g mol À 1 Â 1000 ! ¼ 5:00 Â 10 À 12 mol; ðA:7Þ and given 4 active sites per HMG-CoA reductase enzyme, there are 2.00 Â 10 À 11 moles of enzyme active sites in 1 mg of protein. Given the specific activity of an enzyme, 52 Â 10 À 12 mol min À 1 mg À 1 , we have 52 Â 10 À 12 mol min À 1 mg À 1 2:00 Â 10 À 11 mol mg À 1 ¼ 2:60 min À 1 : ðA:8Þ And so, an approximation to the turnover number for HMG-CoA reductase is given by μ c ¼ 4:33 Â 10 À 2 s À 1 :

ðA:9Þ
A.2. Determination of degradation ratesδ m , δ h and δ c The parameters δ m and δ h determine the degradation rate of the HMGCR mRNA, and HMGCR enzyme molecule respectively. The calculation of these values is based on the half-lives, t 1=2 ,o f the molecules as derived from an exponential decay model. For a decay constant δ, the rate of degradation of the variable is given by : ðA:10Þ Data from Wilson and Deeley (1995) give a half-life of t 1=2 ¼ 4:3h for HMG-CoA reductase mRNA, measured in Hep G2 cells and so δ m ¼ ln 2 15; 480 s ¼ 4:48 Â 10 À 5 s À 1 : ðA:11Þ Data from Brown et al. (1974) give a half-life of t 1=2 ¼ 3 h for the HMG-CoA reductase protein, measured in human fibroblast cells and so δ h ¼ ln 2 10; 800 s ¼ 6:42 Â 10 À 5 s À 1 : ðA:12Þ The parameter δ c determines the effective loss of cholesterol from the system as cholesterol does not degrade in a similar manner to the other species in the system. The value of this parameter may exhibit significant variation dependent on cellular conditions. For an initial value to be used in simulation, this parameter is assumed to be slightly faster than the degradation rates of protein and mRNA and δ c ¼ 1:20 Â 10 À 4 s À 1 is used.

A.3. Determination of binding affinity parametersκ m and κ c
The parameter defining the binding reaction between SREBP and the HMGCR gene is κ m . As no specific experimental data is available for the binding affinity between SREBP and the SRE of the HMGCR gene, we use the results of Yang and Swartz (2011) who quantify the DNA binding affinity for other regulatory transcription factors. Their results provide an estimate of DNA transcription factor binding affinity in the region of 54.2 nmol; these estimates are in keeping with the hypothesis that the interaction between SREBPs and SREs is characterised by relatively high affinity and low saturation kinetics (Amemiya-Kudo et al., 2002). To convert these affinities into units of molecules per ml the following calculation is used. A binding affinity of 54.2 nmol is equivalent to 54:2 Â 10 À 9 1000 moles=ml Â N A ðA:13Þ ¼ 3:26 Â 10 13 molecules=ml; ðA:14Þ and so κ m is taken to be Oð10 13 Þ.
The parameter which describes the binding of cholesterol to SREBP is κ c . Experimental investigations with cholesterol on the SSD contained in SCAP have revealed difficulty in calculating kinetic constants for this interaction owing to the insolubility of both cholesterol and SCAP, and the necessity of performing kinetic experiments in micelles (Radhakrishnan et al., 2004). Data from Radhakrishnan et al. (2004) indicates that the binding reaction between cholesterol and SCAP is saturable and half-maximal binding occurs at approximately 100 nmol. This value is converted into units of molecules per ml using the following calculation.

A.4. Binding coefficientsx and y
The work of Vallett et al. (1996) suggests that three binding locations on the HMGCR gene are available to SREBP and hence x ¼3. Experimental investigations on the SSD contained in SCAP have revealed its tetrameric nature (Radhakrishnan et al., 2004(Radhakrishnan et al., , 2008. That is, four molecules of cholesterol may bind to a SCAP-SREBP complex in order to promote inactivation and y¼4. A.5. Concentration of total transcription factors 0 In general, transcription factors have low gene expression, and are therefore present in relatively low concentrations within the cell (Sanguinetti et al., 2006). Since the SREBP signalling pathway has been simplified it is necessary to make some assumptions to obtain an estimate for this value. The majority of cellular cholesterol is located within the plasma membrane of the cell. Its concentration is governed by multiple regulatory proteins regulated in the ER which are themselves under the control of a small regulatory pool of ER cholesterol (Lange et al., 1999). The concentration of transcription factor, s 0 , is assumed to be of the order of this concentration, which allows small changes in the cellular cholesterol pool to be magnified in the transcription factor pathway. Lange et al. (1999) give an average of 0.45 nmol ER cholesterol per mg cell protein.