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

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 findings 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 3-hydroxy-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 LDL-cholesterol 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).

Fig. 1

Genetic regulation of cholesterol biosynthesis by SREBP-2. Hepatocytes synthesise HMGCR mRNA which in turn is translated into the enzyme HMGCR. HMGCR catalyses the synthesis of cholesterol which in turn influences its own transcription rate by interacting with the transcription factor SREBP; the transcription rate increases when cholesterol is low in the cell and declines when cholesterol is high. (SRE – sterol regulatory element; MH – HMGCR mRNA; C – cholesterol).

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'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 (1963, 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 non-lipoprotein 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.

Fig. 2

The genetic regulation of cholesterol production by SREBP-2. The HMGCR gene is transcribed at a rate to produce HMGCR mRNA . This is translated at a rate into the HMGCR enzyme . HMGCR then goes on to catalyse the reaction creating the metabolite cholesterol at a rate . This process is under the control of the transcription factor SREBP which acts as a transcriptional activator for the pathway. Under conditions where cholesterol is in excess forms an inactive complex with and transcription of the target gene declines. HMGCR mRNA, HMGCR and cholesterol are degraded at rates , and , respectively.

The binding of SREBP-2 to the gene, subsequent transcription and translation to HMG-CoA mRNA and production of HMGCR and cholesterol can be described by the reaction equation

Here x is the number of molecules of required to bind to to produce a functional effect. This binding reaction has an association rate and a dissociation rate . is transcribed at a rate and is translated at a rate . The creation of occurs at a rate . , and are respectively the degradation rates of , and .

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 bywhere y is the number of molecules of required to bind to to cause inactivation. This binding reaction has an association rate and a dissociation rate .

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).

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 .

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, , does not affect the concentration of . That is, synthesis of will not deplete .

The governing ODEs equations are derived by application of the law of mass action to the biochemical reactions (1) and (2) which giveswith initial conditionswhere in the above system of equations, we use the following notation in which square brackets denote concentration: , , , , , and .

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 toon using the initial conditions (10). We now assume that Eq. (5) reaches equilibrium rapidly (quasi-steady-state approximation) such thatand using Eq. (11) we havewhich upon rearranging giveswhereHere is the dissociation constant for the reaction between and .

We further observe that adding Eqs. (4), (5) and (9) givesUnder the quasi-steady state assumption of Eq. (12) together with the initial conditions (10) we find thatAlso under the approximation (12) we see that . 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. . We then obtain the following equation from (19):

Finally we assume that the binding reaction between and reaches equilibrium rapidly such thatRearranging this result givesin which we define the constant such thatwhere is the dissociation constant for the reaction between and .

Using Eqs. (14), (20) and (22) to eliminate Eqs. (3)–(5) and (9) from the system equations (3)–(9) we obtain the reduced systemwith the initial conditionsHere where is the maximal rate of transcription.

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 such thatwhere τ represents the non-dimensional time. The remaining variables are rescaled with respect to the concentration of total transcription factor, , such thatThis non-dimensionalisation leads towith the initial conditionswhere the non-dimensional parameters are given byThe non-dimensional parameter values are summarised in Table 2.

Table 2

Non-dimensional parameter values.

Parameter	Description	Nondimensional value
μc	Cholesterol synthesis rate	1.30
μm	HMGCR transcription rate	1.90×10−10
δm	HMGCR mRNA degradation rate	1.35×10−3
δh	HMGCR degradation rate	1.93×10−3
δc	Cholesterol utilisation rate	3.60×10−3
κm	SREBP-2-HMGCR gene dissociation constant	1.00×10−4
κc	SREBP-2-cholesterol dissociation constant	1.00×10−3

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.

Table 1

Summary of model parameter values. Details of parameter derivations are given in Appendix A.

Parameter	Description	Dimensional value
μ¯_m	HMGCR transcription rate	5.17×105 molecules ml−1 s−1
μ¯_h	HMGCR translation rate	3.33×10−2 s−1
μ¯_c	Cholesterol synthesis rate	4.33×10−2 s−1
δ¯_m	HMGCR mRNA degradation rate	4.48×10−5 s−1
δ¯_h	HMGCR degradation rate	6.42×10−5 s−1
δ¯_c	Cholesterol utilisation rate	1.20×10−4 s−1
κ¯_m	SREBP-2-HMGCR gene dissociation constant	O(10^13) molecules ml−1
κ¯_c	SREBP-2-Cholesterol dissociation constant	O(10^14) molecules ml−1
x	Molecules of SREBP-2 binding to gene	3
y	Molecules of cholesterol binding to SREBP-2	4

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 ofwhere m, h and c are the steady state values of m, h and c respectively. Using Eqs. (36) and (37), Eq. (35) can be written aswhereExpanding, we find that the steady states are given by the solution of the polynomial equation of degree , where and . As all parameters are positive, we may apply the results of Descartes' 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 bywithWe note that as all parameter values are positive and that for physiologically valid parameter ranges.

Calculation of the eigenvalues of J requires the solution of the equationwhere λ 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 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's criteria applied to a cubic equation of the form are satisfied if and only if , , and . That is, the necessary and sufficient condition for the roots of Eq. (44) to have negative real part requiresSince 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: . 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: . By substituting this value of into Eq. (44), we now have the characteristic equation given bywhere we have and . Therefore the characteristic equation has two conjugate roots on the imaginary axis and one negative real eigenvalue λ3 given byone 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: . 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.

Fixed point stability – variation of μ

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 μ 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 μ denoted by . 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 , for which the following two conditions are fulfilled, at the equilibrium point :

The matrix J has two complex eigenvaluesin some neighbourhood of and for these eigenvalues are purely imaginary, that is,This non-hyperbolicity condition is a necessary condition for the Hopf bifurcation.

The relation described byholds in some neighbourhood of . 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 , given by the solution of where ψ is given by Eq. (42), together with the equation determining the equilibrium value of c for :

From the results of Case II we know that at this value of the characteristic polynomial Eq. (44) has two purely imaginary roots , given in Eqs. (49) and (50), where

To show that the second condition holds we use the Implicit Function Theorem. For each and the corresponding system, Eqs. (30)–(32), define as a function of two variables μ and λ, where is the characteristic polynomial of the system equations (30)–(32) defined by Eq. (44).

Let the complex eigenvalues be the roots of the characteristic polynomial. Hence, for these eigenvalues we havewhere this represents an implicit function of two variables μ and λ. The Implicit Function Theorem tells us that we may define μ as a function of λ near the point , and the derivative of this function is given byproviding

We begin by computing the derivative of the function with respect to λ, and evaluating this at the critical point . Thus we have,

Simplifying in conjunction with the fact that from Eq. (51), we obtain

Furthermore, from the characteristic polynomial Eq. (44), we havewhere, we have previously noted that . 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 and . 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, and we need only be concerned with the strict inequality .

Eqs. (54) and (55) together with Eq. (53) yield

Upon the rationalisation of the denominator of this complex fraction we obtainand since , and the second condition of the Hopf theorem is fulfilled. Thus we have proved the existence of a Hopf bifurcation at the critical value .

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 , is absorbed by a stable spiral equilibrium. The equilibrium becomes unstable for ; 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 , a limit cycle is born. At this point the equilibrium changes stability and becomes unstable. For , 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 ODEswith an equilibrium point at . The principle of numerical continuation requires finding a solution curve σ of with which describes how the equilibrium point varies. The curve σ is traced by means of a predictor-corrector algorithm and bifurcations along σ 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, , 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,In the case of Eqs. (30)–(32) with μ as the bifurcation parameter, we find that a Hopf bifurcation is predicted to occur at the point (, c⁎)=(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.

implies that the system is attracted to a stable periodic orbit and

implies that the system is attracted to an unstable periodic orbit.

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 μ 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 μ shows the system shifting to Case Ib. Thus we see oscillatory convergence to a steady-state as shown in Fig. 4. Still further increases in μ see the system reaching Case II, where we have pure oscillations in mRNA, HMGCR and cholesterol; this is illustrated in Fig. 5. The Hopf bifurcation occurs at the transition from Case Ib to Case II.

Fig. 3

Stable node equilibrium (corresponding to Case Ia) where there are three negative real eigenvalues; variable concentrations exhibit monotonic convergence towards a steady state value. All parameters are as in Table 2 except nondimensional μ=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 μ which has been increased 2 fold to μ=0.923.

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 μ which has been increased approximately 4 fold to μ=1.946.

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 μ 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 .

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, μ and μ, 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, , or binding affinities , 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 δ we see that the critical value for which a Hopf bifurcation may occur is given by the solution of the equationwhich is quadratic in 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 which are determined from

The eigenvalues at this critical point are given byand so the first Hopf bifurcation condition holds. Proceeding in the manner of the calculation for , we findand 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 δ 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 δ 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.

Fig. 6

The response of mRNA concentration to variation of δ. 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 δ which is successively varied with nondimensional values as indicated. Initial conditions are as in Fig. 3.

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 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 is the mRNA concentration, 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 . 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).