A model of flux regulation in the cholesterol biosynthesis pathway: Immune mediated graduated flux reduction versus statin-like led stepped flux reduction.

The cholesterol biosynthesis pathway has recently been shown to play an important role in the innate immune response to viral infection with host protection occurring through a coordinate down regulation of the enzymes catalysing each metabolic step. In contrast, statin based drugs, which form the principle pharmaceutical agents for decreasing the activity of this pathway, target a single enzyme. Here, we build an ordinary differential equation model of the cholesterol biosynthesis pathway in order to investigate how the two regulatory strategies impact upon the behaviour of the pathway. We employ a modest set of assumptions: that the pathway operates away from saturation, that each metabolite is involved in multiple cellular interactions and that mRNA levels reflect enzyme concentrations. Using data taken from primary bone marrow derived macrophage cells infected with murine cytomegalovirus or treated with IFNγ, we show that, under these assumptions, coordinate down-regulation of enzyme activity imparts a graduated reduction in flux along the pathway. In contrast, modelling a statin-like treatment that achieves the same degree of down-regulation in cholesterol production, we show that this delivers a step change in flux along the pathway. The graduated reduction mediated by physiological coordinate regulation of multiple enzymes supports a mechanism that allows a greater level of specificity, altering cholesterol levels with less impact upon interactions branching from the pathway, than pharmacological step reductions. We argue that coordinate regulation is likely to show a long-term evolutionary advantage over single enzyme regulation. Finally, the results from our models have implications for future pharmaceutical therapies intended to target cholesterol production with greater specificity and fewer off target effects, suggesting that this can be achieved by mimicking the coordinated down-regulation observed in immunological responses.

Introduction

Cholesterol is central to a diverse range of cellular functions, including membrane development and maintenance [1], lipid raft formation and vesicular transport [2], steroid hormone synthesis [3], neurological development [4], and oxysterol [5] and vitamin D synthesis [6]. Recently, the cholesterol metabolism has been shown to have an important role in host–pathogen interactions. It has been documented to be perturbed in response to infection [7,8] and, conversely, cholesterol and its associated metabolites have been shown to alter inflammatory mediators [9,10].

Cholesterol synthesis is one step in a pathway of metabolic interactions that is subject to catalytic regulation [11] and evidence suggests that this pathway is critical to the optimal growth of a range of viruses and microbes including cytomegalovirus (CMV), Hepatitis C (HCV), HIV, Japanese Encephalitis (JEV), West Nile (WNV), Dengue (DENV), Measles viruses (MV), African Swine Fever Virus (ASFV), Mycobacteria and Salmonella [7,8,12-19].

The cholesterol biosynthesis pathway itself comprises a sequence of metabolic interactions that occur across several organelles, starting with the processing of Acetyl-Coenzyme A (henceforth denoted ACoA – see Supplementary section 1 for a list of all metabolite abbreviations) in the mitochondria and ending with cholesterol synthesis in the endoplasmic reticulum [20,21]. This pathway branches in the peroxisome and endoplasmic reticulum into the sterol arm and the non-sterol arms (prenylation and dolichylation), the latter arms carrying flux away from the main sterol arm.

Coordinate transcriptional control of the enzymes of the cholesterol biosynthesis pathway is mediated by SREBP2 and feedback control occurs through regulation of SREBP2 transport. The SCAP:SREBP2 complex is ordinarily chaperoned to the Golgi complex where SREBP2 is cleaved, before it migrates to the nucleus to activate the suite of enzymes associated with the pathway. However, in the presence of relatively high concentrations of intracellular cholesterol or side-chain hydroxylated cholesterol, in particular 25-hydroxycholesterol, SCAP:SREBP2 is retained instead in the endoplasmic reticulum. Retention of SCAP:SREBP2 acts to down-regulate transcription of the enzymes acting on the pathway until ordinary levels of cholesterol and its derivatives have been restored. Hence, the pathway undergoes transcriptionally mediated regulation through changes to enzyme concentrations [20].

Recently, we reported a modest, but statistically significant, decrement in the concentrations of enzymes associated with the cholesterol biosynthesis pathway, in response to both infection and interferon treatment in macrophages. This was observed at the transcriptional level and was shown to correlate with reduced protein concentrations [12]. This decrement was found to be part of the innate immune response, intended to suppress viral growth. However, the mechanism through which such changes decrease the activity of the cholesterol biosynthesis pathway is something that has yet to be fully elucidated in the published literature.

We have sought to investigate this mechanism of regulation, exploring the impact on flux that results from such enzyme decrements. This problem is experimentally challenging, but tractable with computational methods.

Flux is the natural quantity to consider when studying metabolic pathway function and flux studies have been employed both theoretically [22,23] and experimentally [24,25]. The flux through the pathway describes the stoichiometrically adjusted rate of production of each metabolite and so captures whether and how the production rate of the metabolites affect each other. Ultimately, the final flux value in the pathway describes the rate of cholesterol synthesis. Metabolic Control Analysis (MCA) and Flux Balance Analysis (FBA) are two typical approaches to studying flux in a pathway system. However, MCA approaches focus on the effects of individual infinitesimal changes in enzyme activity rather than the compound effects of multiple finite changes and FBA, in its standard form, does not relate flux changes to substrate concentration changes. As a result, they are inappropriate for our study in which we validate the pathway model at the level of the substrate, implement multiple finite enzyme decrements and model the effects of chemical inhibition. We build an ordinary differential equation (ODE), dynamical model of the sterol pathway using Michaelis–Menten and mass action kinetics that incorporates additional interactions to represent the consumption of metabolites in non-sterol related processes. We demonstrate that multiple small decreases in enzyme activity can suppress the flux through the main cholesterol biosynthesis pathway. This suppression presents itself as a graduated reduction when the profile of flux is considered.

Cholesterol levels have also been demonstrated to be an important risk factor in cardiovascular disease [26,27] and their control is an active area of research [28]. Current therapies involve the use of statins to competitively inhibit the enzyme HMGCR which is responsible for catalysis of the interaction transforming 3-hydroxy-3-methyl-glutaryl Coenzyme A (HCoA) to Mevalonate (M). However, the efficacy of such therapies is limited by drug toxicity and off target effects [29-31]. Here, we show that statin treatment regulates the flux through the pathway in a manner that is markedly different to that following infection. The metabolic interaction catalysed by HMGCR is significantly upstream of cholesterol biosynthesis and we show that the impact of a statin-like treatment is to suppress flux throughout most of the pathway. This impacts significantly upon many of the metabolites upstream of cholesterol and upon the non-sterol arms, thereby incurring off-target effects. In contrast, because coordinate enzyme regulation leads to a graduated reduction in flux along the pathway, it has a less dramatic impact upon the branches upstream of cholesterol production.

This manuscript is organized as follows. In Section 2, we describe the experimental and mathematical methods employed to determine enzyme and metabolite levels in response to infection and IFNγ treatment and to model the pathway. In Section 2.1, we describe the experimental method and in Sections 2.2 and 2.3, we describe how the model was built, how the initial conditions were defined and how the model was used to simulate pathway activity. In Section 3, we present the results of using the model to study the flux through the pathway, with Sections 3.1–3.4 describing the validation of the model and the impact on the flux of the response to IFNγ treatment, to CMV infection and to statin intervention, respectively. In Section 4, we discuss these results, their relationship, the off-target effects and their implications for specific, targeted regulatory strategies. In Section 5, we summarize our results. Supplementary material in support of the results presented here is available online.

Materials & methods

Experimental measurements

Enzyme levels were inferred from gene expression measurements of bone-marrow derived macrophage cells in two time course experiments, one in which cells were infected with murine cytomegalovirus (mCMV) and one in which cells were treated with IFNγ. Measurements were taken at half hour intervals for 12 h using Agilent microarray platforms and at 24 h for select members using QPCR. Agreement between mRNA expression and protein concentrations was validated by quantitative western blotting for selected members of the pathway [12].

Intracellular cholesterol concentration was determined enzymatically using the Amplex-Red cholesterol assay kit (Molecular Probes) according to manufacturer recommendations. Briefly, cells were washed with 1 ml ice cold PBS and then lyzed in 200 μl cold Lipid buffer containing 0.5 M of potassium phosphate, pH 7.4, 0.25 mM cholic acid, and 0.5% triton X-100. Cell lysates were sonicated on ice with three 10-s pulses at high intensity. 20 μl were then used to determine protein concentration using a standard BSA assay to normalize the protein concentration. For cholesterol measurement, 20 μl of each sample were added to 80 μl assay solution, which contained 300 μM Amplex Red reagent, 2 U per ml HRP and 2 U per ml cholesterol oxidase, 0.1 M of potassium phosphate, pH 7.4, 0.05 mM cholic acid, and 0.1% triton X-100. After preincubation for 30 min at 37 °C under light exclusion conditions, fluorescence was measured using excitation at 530 ± 2.5 nm and fluorescence detection at 590 ± 2.5 nm with a Polarstar Optima Multifunction Microplate Reader (BMG Labtech, UK). The values were corrected from the background. The relative amount of free cholesterol to the mock treated samples was calculated using the manufacturer's supplied standard curve.

For the measurement of metabolite concentration, an Agilent high performance liquid chromatography (HPLC) system coupled with an Applied Biosystem Triple Quadrupole/Ion Trap mass spectrometer (4000Qtrap) was used for quantification of individual polar lipids (phospholipids and sphingolipids). Electrospray ionization-based multiple reaction monitoring (MRM) transitions were set up for the quantitative analysis of various polar lipids. HPLC atmosphere chemical ionisation (APCI)/MS was carried out for analysis of sterols [12].

Model construction

Because regulation and feedback occur through transcriptional control of enzyme activity, we chose to model the impact of enzyme activity on the pathway flux. This obviated the need to explicitly consider SREBP2 mediated feedback in the pathway as any feedback would be accounted for in our measurements of enzyme activity. From the representations available in the KEGG pathway database [11], we assembled the description of the pathway shown in Fig. 1A (presented in SBGN notation [32]). From this description, we built a deterministic model of the pathway in which catalysed metabolic transitions were described with Michaelis–Menten kinetics and autocatalysed metabolic transitions were described with mass action kinetics. Metabolites play a role in a range of cellular processes and undergo degradation. To capture this, each metabolite was also considered to be consumed in an interaction competing with the main pathway, but at a rate lower than its consumption in the main pathway. The competing reactions were modelled with mass action kinetics (not shown in Fig. 1A).

Fig. 1

A) The cholesterol pathway represented in SBGN notation, starting with the metabolites Acetyl-Coenzyme A (ACoA) and ending in cholesterol synthesis. B)–D) The normalized concentrations of 14-lanosterol (B), zymosterol (C) and cholesterol (D) at 0 h (solid, black) and 12 h (diagonal stripes) after mCMV infection and after IFNγ treatment. We show results from experiment and simulation. Experimental measurements were normalized against measurements from a mock time course and simulated measurements were normalized against the concentration at 0 h.

Parameter values for the cholesterol biosynthesis pathway were taken from the Brenda enzyme database [33]. Where parameters were not known, they were approximated with the mean of the corresponding known parameters. The mean values were  hr−1 and  mM. Normalized enzyme levels were inferred from the microarray time courses described above. To obtain an absolute scale for these measurements, we assumed an average number of 5000 enzyme proteins [34] in a region of the cell (the endoplasmic reticulum) of volume 10−14 l [35-37]. This gave a concentration scale of  mM which we took to correspond to the mean normalized enzyme level measurement across both experiments at the start of each time course (mean normalized measurement = 1279.2). We subsequently transformed the normalized expression levels to concentrations using this equivalence and took expression levels to be commensurate with protein levels following the validation described above.

From Supplementary section 2, we can see that, in the limit of low substrate, a Michaelis–Menten interaction acts much like a mass action interaction with a rate constant k = kcatE/k, where E is the enzyme concentration. Using , and , this gave k = 156 h−1 and we used this value as the rate constant for all the autocatalysed, mass–action interactions in the main pathway.

In order for the behaviour of the cholesterol biosynthesis pathway activity to dominate over the other dynamical effects, the competing interactions were taken to have mass action rate constants two orders of magnitude smaller than the corresponding mass action constants formed from the main pathway parameters. Thus, for substrate i, consumed in a Michaelis–Menten interaction, the corresponding mass action constant would be and we assigned to the mass action constant of the competing interaction the value , where and are the turnover and Michaelis–Menten constants of the Michaelis–Menten pathway interaction consuming the metabolite and is the enzyme concentration at the start of the time course. Cholesterol was assumed to be consumed at the same rate that it was created in order to avoid accumulation.

A complete list of parameter values is shown in Supplementary sections 3 and 4 and a complete list of normalized enzyme measurements is shown in Supplementary sections 5 and 6.

Modelling strategy

Near saturation, small fluctuations in enzyme concentration lead to small changes of Vmax in Michaelis–Menten interactions (Vmax = kcat × enzyme concentration) and this, in turn, leads to large changes in metabolite concentration. Hence, we would expect the pathway to operate in a regime away from saturation where the dynamics would be more stable and robust. Here, we assumed that the pathway was operating away from saturation of the Michaelis–Menten interactions.

We assumed a level of flux into the pathway ∼2/3 rds of the lowest Vmax value obtained across both time courses. From the parameters and the enzyme concentrations at the start of each time course, we calculated the concentrations of each metabolite that would allow the pathway to continue in dynamic equilibrium if no changes were observed in the enzyme concentrations (see Supplementary section 7). The pathway was then numerically integrated, allowing the enzyme concentrations to change in accordance with the concentrations measured in the time course. In the interval between time points, enzyme concentrations were calculated through linear interpolation.

In each interaction of the pathway, the flux was defined and calculated as the rate at which each metabolite is produced. Because this pathway is stoichiometrically trivial, this allowed a direct comparison between all interactions in the pathway. Numerical integration of the pathway was conducted in a two step process, with all the fluxes being calculated from metabolite concentrations and then all the metabolite concentrations being updated from the net fluxes. For a general interaction in the interior of a sequential pathway the update rule took the form , where m is the concentration of metabolite i. Modifications to this rule were required at branch points in the pathway and at the start and end. In order to determine a size for Δt, simulations were run in iterations, with values for Δt progressively decreasing. Iterations continued until a value for Δt was reaching at which the results stabilized with no variation in the first four significant figures of the pathway output. This determined the Δt size.

Results

Model validation

Our first step was to assess the quality of the model by testing whether it behaved in a manner consistent with the observed underlying biology. We did this by comparing the concentrations of three metabolites that cover the pathway at 12 h following immune challenge.

In Fig. 1B–D, we see normalized metabolite concentrations at 0 and 12 h following mCMV infection or following IFNγ treatment, from both experiment and simulation. Experimentally determined concentrations were normalized against the mock treatment time course; computationally determined concentrations were normalized against the concentration calculated at the start of the time course, determined as part of the initial conditions (see Section 2.3).

From a comparison of the experimental and the computationally determined values, we can see that the behaviour of the model is in qualitative agreement with the experimentally observed response of the macrophage cells to mCMV infection and IFNγ treatment. This agreement becomes even clearer if we consider the values at 0, 12 and 24 h post infection or post treatment (Supplementary Fig. 1). This gives us confidence that the model can be used to address questions surrounding the relationship between changes in enzyme levels, metabolite concentrations and flux.

In response to IFNγ treatment, the flux through the cholesterol biosynthesis pathway is significantly suppressed in a graduated manner

It is not known how coordinate enzyme control impacts upon flux through the pathway and cholesterol biosynthesis. Thus, we first chose to assess the response of the pathway to the enzyme time courses measured following IFNγ treatment.

We simulated the pathway activity by taking the measurements from the microarray time course of macrophage cells following IFNγ treatment to represent enzyme concentrations (we have previously shown good correlation between mRNA concentrations and enzyme levels [12]). Fig. 2A shows the resulting profile of flux along the pathway and how this profile develops over the 12 h. For presentational purposes, we numbered the interactions with 1 representing the input flux and 17 representing cholesterol synthesis (for the full numbering, see Supplementary section 12). The pathway forks at Zymosterol with flux split down both forks. Cholesterol synthesis occurs on both forks and for presentational simplicity, we omitted the details of each fork, but retained the cholesterol synthesis rate, calculated as the sum of the synthesis rates from each fork.

Fig. 2

The flux through the cholesterol biosynthesis pathway following treatment with IFNγ. A) The development of the flux through the pathway in simulation is shown from 0 to 12 h following treatment. Interactions are numbered from 1 (the input flux) to 17 (cholesterol production). For the full numbering, see Supplementary section 12. At 0 h, the flux through the pathway is relatively constant. However, by 12 h the flux has been significantly suppressed along the pathway. B) The profile of flux through the pathway at 0 and 12 h following treatment. These profiles represent cross sections of the surface shown in A). The flux is dramatically reduced in the first 12 h. Interactions can be classified as dominant (Squa-23Ox) and non-dominant (the remainder) depending on their degree of impact on the pathway flux. C) The flux through the non-dominant interactions between ACoA–HCoA and FPP-Squa, normalized against the flux through the ACoA–HCoA interaction. The flux through these non-dominant interactions is suppressed more modestly than in the dominant interactions. D) The flux through the non-dominant interactions between Squa-23Ox and 4MZC-3K4M normalized against the flux through Squa-23Ox. The flux through these non-dominant interactions is suppressed more modestly than in the dominant interactions.

From Fig. 2A, we can see that the profile of flux is relatively constant across the pathway at the point of treatment, but that, as time advances, the flux along the pathway becomes suppressed, leading to a much reduced rate of cholesterol synthesis at 12 h post treatment.

In order to explore the down-regulation in pathway activity further, we examined the cross sections taken at 0 and 12 h post treatment. The resulting profiles are shown in Fig. 2B. At the point of treatment (0 h), we can see that the pathway undergoes a very modest reduction in flux along its length, attributable to the flux lost through the interactions competing for each metabolite. However, as a result of the coordinate down regulation of enzyme activity across the time course, this modest reduction is dramatically amplified.

This profile of flux reduction along the pathway can be considered in terms of (a) dominant interactions, in which the flux reduces significantly and (b) non-dominant interactions, in which the flux reduces more modestly. From Fig. 2B, we can see that the dominant interaction is Squalene-2,3Oxidosqualene (henceforth Squa-23Ox). In order to further explore the degree of suppression in pathway flux that takes place in the non-dominant interactions, we investigated the profile of flux leading up to Squa synthesis and from Squa to 3-keto-4-methyl-zymosterol (henceforth 3k4m), normalizing the flux values against the flux through the first interaction in the sequence. The flux profiles are shown in Fig. 2C and D, respectively. In both Fig. 2C and D, we can see that the flux through the sequence of interactions is significantly suppressed after 12 h.

In response to mCMV infection, the flux through the cholesterol biosynthesis pathway shows some suppression in a graduated manner

Ligand activation of the IFNγ gamma receptor by the IFNγ cytokine is involved in immune activation of macrophages. To compare modulation of the pathway in response to IFNγ with its modulation in response to infection, we next modelled the pathway activity using the time course recorded in response to mCMV infection. Fig. 3A shows the resulting profile of flux along the pathway and how this profile developed over time. As in Fig. 2, we numbered the interactions with 1 representing the input flux and 17 representing cholesterol synthesis (see Supplementary section 12). Because the pathway forks at Zymosterol, we omitted the details of the flux along each fork, but retained the cholesterol synthesis rate, calculated as the sum of the synthesis rates on each fork.

Fig. 3

The flux through the cholesterol biosynthesis pathway following infection with mCMV. A) The development of the flux through the pathway in simulation is shown from 0 to 12 h post infection. Interactions are numbered from 1 (the input flux) to 17 (cholesterol production). For the full numbering, see Supplementary section 12. At 0 h, the flux through the pathway is relatively constant. However, by 12 h the flux has been significantly suppressed along the pathway. B) The profile of flux through the pathway at 0 h and 12 h post infection. These profiles represent cross sections of the surface shown in A). The flux is dramatically reduced in the first 12 h following infection. Interactions can be classified as dominant (ACoA–HCoA and Squa-23Ox) and non-dominant (the remainder) depending on their degree of impact on the pathway flux. C) The flux through the non-dominant interactions between ACoA–HCoA and FPP-Squa, normalized against the flux through the ACoA–HCoA interaction. The flux through these non-dominant interactions shows a mild suppression towards the end of the sequence. D) The flux through the non-dominant interactions between Squa-23Ox and 4MZC-3K4M normalized against the flux through Squa-23Ox. The flux through these non-dominant interactions shows no suppression between 0 and 12 h.

The profile of flux in Fig. 3A is relatively constant across the pathway at the moment of infection, but as time increases, the flux along the pathway reduces, leading to a suppressed pathway and a much reduced rate of cholesterol synthesis. To analyse this response further, we again took cross sections from this surface at 0 and 12 h post infection. In Fig. 3B, we can see that there is a modest reduction in flux at the point of treatment (0 h) and that this is amplified as a result of the response to treatment.

As mentioned previously, the profile of flux reduction along the pathway can be considered in terms of interactions that make a dominant contribution to flux reduction and those that make a non-dominant contribution. Interestingly, we see that the distribution of dominant interactions is distinct from the distribution seen in the response to IFNγ treatment. Here, the dominant interactions are ACoA–HCoA and Squa-23Ox. In order to explore further the regulation of flux through the non-dominant interactions, we examined the flux profiles between dominant interactions, normalizing the flux through each sequence against the flux through the first interaction in the sequence. The results can be seen in Fig. 3C and D. In the first of these profiles, it is clear that the flux at 12 h post infection is similar to that at infection, but that some modest reduction does occur towards the end of the sequence as a result of the enzyme regulation. In the second profile, the flux profile clearly does not alter significantly between 0 h and 12 h post infection.

Under statin-like intervention, the pathway acquires a step reduction in flux

From Figs. 2 and 3, we can see that the pathway has the ability to respond to perturbation throughout its length. We chose to compare these flux profiles to that which is likely to be obtained when statin-like inhibitors that target the HMGCR enzyme are introduced into the pathway.

Fig. 4A shows the profile of flux obtained in the unperturbed pathway and the profiles at 12 h following IFNγ treatment and mCMV infection. Fig. 4A also shows the profile of flux obtained when the statin-like inhibitor is introduced with a concentration sufficient to suppress cholesterol production to the mean of the production rates for IFNγ treatment and mCMV infection at 12 h. Here we can see that the profile takes a very different form, impacting dramatically upon the interactions upstream of Squa. The profile corresponding to statin-like inhibition is less graduated than in the physiological responses both in terms of the dominant interactions and the non-dominant interactions, and such a flat profile arises because the single HMGCR inhibited interaction takes the entire regulatory burden of reducing the flux along the pathway with the other interactions contributing no more flux reduction than in the unperturbed case.

Fig. 4

A) The effect of a statin-like inhibitor on pathway activity. The profile of flux at 0 and 12 h following mCMV infection and IFNγ treatment together with the profile of flux in an unperturbed cell following the introduction of a statin-like inhibitor which targets the enzyme HMGCR. The effect of a statin-like inhibitor is to step down the flux through the interactions catalysed by HMGCR. This impacts upon the pathway significantly upstream of the point of cholesterol synthesis and creates a flux profile dramatically different to that which arises from the biological response to mCMV infection or IFNγ treatment. B) The profile of flux achieved along the pathway when inhibitor concentrations are chosen so that each interaction contributes equally to the regulation of flux (inhibitor levels listed in Supplementary section 9).

Discussion

Figs. 2–4 give us insight into how the cholesterol biosynthesis pathway operates and its various modes of regulation. Fig. 2 shows that, in response to IFNγ treatment, the pathway undergoes a reduction in flux that is mediated by both the dominant and non-dominant interactions. The distribution of flux at the point of IFNγ treatment shows a very mild reduction along the pathway. This is attributable to the flux lost through the off-pathway, competing interactions. However, at 12 h post treatment, we can see that the flux becomes significantly suppressed along the pathway due to the coordinate enzyme reduction. This occurs because a reduction in the enzyme concentration in an interaction has the effect of decreasing the rate of the interaction. This, in turn, increases the proportion of flux shunted through the off-pathway interaction competing for the same metabolite. At 12 h post treatment, we can see that most of the interactions on the pathway play a suppressive role. This implies that, for most of the metabolites, there is an increase in the proportion of flux shunted through the off-pathway interactions, if flux conservation is to be maintained.

In Fig. 3, we can see that the pathway also undergoes a reduction in flux in response to mCMV infection. In this case, the reduction is mediated mostly by the dominant interactions with very few non-dominant interactions contributing to the suppression of flux between 0 and 12 h post infection. Here, we infer that the coordinate enzyme control leads to an increased proportion of flux being shunted through the off-pathway interactions that compete for metabolites with the dominant interactions. However, we see little suppression of flux through the non-dominant interactions between 0 and 12 h post infection, implying that the proportion of flux being shunted through the off-pathway interactions, competing for the same metabolites, stays approximately the same as at the point of infection.

In Fig. 4A, we see that a statin-like pharmaceutical intervention comparable to the immune-led responses introduces a profile of flux with a significantly different shape to that of the immune led response. A statin-like intervention leads to a step down in flux towards the start of the pathway with the result that the flux passing through many of the interactions in the upper half of the pathway is significantly lower than the levels observed in immune-led responses. We can infer that a statin-like intervention shunts a significantly increased proportion of flux through the off-pathway interaction competing for HCoA. In the remaining interactions, no significant change in the proportion of flux passing through the main pathway interactions occurs.

If we consider the flux profiles shown in Fig. 4A, we see that the introduction of a statin-like inhibitor steps down the flux at the interaction HCoA-M and that this has consequences for many of the interactions downstream of HCoA-M. In the upper regions of the pathway, the level of flux after the introduction of a statin-like inhibitor is well below that achieved in immune-led responses. This is particularly significant for the interactions that consume and produce Isopentyl-PP (IsPP) and Farnesyl-PP (FPP) as these metabolites are also consumed in the prenylation and dolichylation non-sterol arms that fork from the main sterol pathway [20,21]. Hence, the difference in profiles between statin-like intervention and the immune-led regulation is likely to have significant consequences for these non-sterol arms due to the difference in flux passing through the branch point in the pathway. Off target effects are known to be a particular concern for statin based treatments [29-31] and the difference in flux profiles may provide a clue as to why this is the case.

If we consider the profiles of flux associated with coordinate enzyme regulation and with single enzyme regulation, it is clear that the coordinate regulation provides a more robust strategy for regulating the pathway. Single enzyme regulation offers no protection against metabolite surges occurring downstream of the regulated interaction. Therefore, distributing flux reduction over multiple interactions, rather than a single interaction, protects the regulatory mechanism from cellular perturbations. The regulation of the cholesterol biosynthesis pathway by IFNγ has been shown to be part of an antiviral response and a distributed flux reduction also increases the difficulty for a pathogen to manipulate the pathway to proviral effect. For these reasons, it is likely that coordinate control has enjoyed a selective advantage in the evolution of immune regulation.

From Figs. 2 and 3, we can see that all the interactions have the potential to play a regulatory role. Therefore, an important question to ask is whether a combination of inhibitors can be chosen such that the flux reduction is distributed evenly throughout the pathway. Such a regulatory structure would be optimally robust. This could be achieved through control at the transcriptional level, tuning SREBP2 activity, or at the enzyme level using pharmaceutical inhibition or miRNA induced degradation. Such an arrangement can be achieved in the models described above by using the combination of inhibitors with levels described in Supplementary section 9. The resulting flux profile is shown in Fig. 4B. Here, the inhibitor concentrations have been chosen so that, for each metabolite, the proportion of flux shunted down the off-pathway interaction to flux continuing down the main pathway is exactly the same. In addition to the superior robustness of such a distributed regulation, this profile of flux would also show greater specificity than a stepped, statin-like treatment as the upstream interactions would be significantly less flux suppressed, reducing the impact on off-pathway interactions. It is possible to speculate that a regulatory arrangement that targets a single interaction towards the end of the pathway would demonstrate even greater specificity to cholesterol as it would create a flux profile in which the flux remains at unregulated levels for most of the pathway, but steps down immediately before the point of cholesterol synthesis. However, such a regulatory structure would not enjoy the robustness of multiple interaction regulation. The improved specificity and robustness that comes with multiple interaction regulation suggests that future pharmacological therapies may be able to act more efficiently and with fewer side effects, if they are designed either to mimic the physiological response of IFNγ treatment or to create a distributed regulation of flux comparable to Fig. 4B.

Further interesting detail can be observed when we allow our simulations to run up to 24 h post infection or post IFNγ treatment. At 24 h we had access to only a limited number of QPCR measurements for a subset of the enzymes on the pathway, so to provide values for the unmeasured enzymes at 24 h, we assumed that the concentration did not change from the 12 h time point. From Supplementary Fig. 1, we can see that qualitative agreement between our simulated metabolite levels and the measured values improves in the 12–24 h window, further supporting the use of this model to explore the dynamics of the cholesterol biosynthesis pathway. Exploring the flux profiles at 24 h post IFNγ treatment and post mCMV infection (Supplementary Figs. 2 and 3, respectively), we can see that there is a difference in the relative timing of the flux suppression. In mCMV infection, pathway flux is suppressed at 12 h post infection and continues to be further suppressed over the 12–24 h interval. However, in response to IFNγ treatment, we see that the pathway starts to recover its pretreated profile of flux in the 12–24 h window. The simplest explanation for this difference would be that there is an inherent delay between mCMV infection and the induction of IFNγ signalling within cell populations. This would predict that running the mCMV infection experiments for a longer duration would yield enzyme expression measurements that would lead to a rise in pathway flux similar to that observed in IFNγ treatment. Although at 12 h post mCMV infection there was little evidence of the non-dominant interactions contributing to pathway regulation, it is also worth noting that, at 24 h post mCMV infection, the non-dominant interactions start to show a modest suppression of flux comparable to that seen in the 0–12 h window following IFNγ treatment. This further supports the idea that the two profiles are related.

The flux profiles shown in Figs. 2A and 3A and in Supplementary Figs. 2A and 3A exhibit interesting features. The broad trend is towards pathway suppression. However, there are points at which the flux exceeds unperturbed levels and where the flux downstream exceeds the flux into the pathway. This may be attributable, at least in part, to technical noise in the time course microarray measurements. However, we cannot rule out that this is biological. Points where the flux consuming a metabolite exceeds the flux producing it correspond to a reduction in the concentration of a metabolite (data not shown). In Figs. 2B and 3B, we see that there is a rise in flux at the 3K4M–4MZ interaction. In our simulations, this corresponded to a reduction in the concentration of 3K4M (data not shown).

Making quantitative predictions of pathway function requires a comprehensive set of high confidence parameters. Without such a parameter set, we are restricted to qualitative observations. Despite the importance of the cholesterol biosynthesis pathway to innate immunity and cardiovascular health and despite its value to industry as the target pathway of the statin class of drugs, we can see from the Table in Supplementary section 3 that the parameterization of this pathway is largely incomplete. This is a surprising result. The lack of a rigorous parameterization impedes studies of pathway behaviour and the development of more advanced scientific and therapeutic interventions. The incomplete nature of the parameterization led us to use estimates for the unknown parameters and this has implications for the flux profiles we obtained in simulations. The particular arrangement of dominant interactions and non-dominant interactions, for example, is parameter dependent and we can see from Supplementary section 3 that the two interactions identified as playing a dominant role in flux reduction (ACoA–HCoA and Squa-23Ox) are those with the lowest turnover (kcat) values. Accurate predictions of the precise arrangement of dominant and non-dominant interactions require a complete, high confidence parameter set. However, under the assumptions made in the model, it is a parameter independent result that all the interactions can contribute a modulation of flux.

The coordinate down-regulation of enzyme activity has been demonstrated to be an innate, immune, anti-viral response to mCMV infection mediated by IFNγ, rather than a pathogenic intervention [12] and so it would appear that flux reduction is also part of the innate immune response. The cholesterol biosynthesis pathway has been demonstrated to be critical to the optimal viral infection of a further range of viruses (Hepatitis C (HCV), HIV, Japanese Encephalitis (JEV), West Nile (WNV), Dengue (DENV), Measles viruses (MV) and African Swine Fever Viruses (ASFV) [7,8,12-16,19]). Since the flux suppression is in response to IFNγ, a general immunological cytokine, it is possible that such a flux suppression might be part of a general antiviral response.

Conclusion

Here, we have shown, with a simple dynamical model, that the coordinate down-regulation of the enzymes in the cholesterol biosynthesis pathway in macrophages, observed in response to mCMV infection and IFNγ treatment, leads to a graduated reduction in flux along the pathway that can be considered in terms of dominant interactions that significantly suppress the flux and non-dominant interactions that make modest, but not insignificant contributions to flux suppression. The model was built using a combination of Michaelis–Menten and mass action interactions and was validated using experimental measurements of metabolite concentrations at a range of time points. Our result suggests that all the interactions in the pathway could be candidate intervention points for future therapeutic strategies.

We were able to compare the impact on flux of coordinate regulation which acts to regulate multiple interactions to the impact of statin-like inhibitors which regulate a single enzyme. This allowed us to explore the advantages of multiple interaction regulation over single enzyme regulation. The graduated profile generated by coordinate regulation enjoys greater robustness and greater specificity and we highlight the importance of this for the prenylation and dolichylation non-sterol arms that fork from the cholesterol biosynthesis pathway. We suggest that such multiple interaction regulation could be exploited in future pharmaceutical therapies in order to improve their efficacy and specificity, perhaps mimicking the regulation of the pathway observed in immunological responses.

Author contribution

AM, SW, HG, KAR and PG compiled a description of the cholesterol biosynthesis pathway. SW, MLG, JH and PG modelled the pathway. SW, LL and MLG compiled parameters for the pathway. MB, GS and MRW completed the experiments. Analysis was completed by SW and PG. SW, KAR and PG wrote the manuscript with feedback from all authors.