Robust homeostasis of cellular cholesterol is a consequence of endogenous antithetic integral control

Although cholesterol is essential for cellular viability and proliferation, it is highly toxic in excess. The concentration of cellular cholesterol must therefore be maintained within tight tolerances, and is thought to be subject to a stringent form of homeostasis known as Robust Perfect Adaptation (RPA). While much is known about the cellular signalling interactions involved in cholesterol regulation, the specific chemical reaction network structures that might be responsible for the robust homeostatic regulation of cellular cholesterol have been entirely unclear until now. In particular, the molecular mechanisms responsible for sensing excess whole-cell cholesterol levels have not been identified previously, and no mathematical models to date have been able to capture an integral control implementation that could impose RPA on cellular cholesterol. Here we provide a detailed mathematical description of cholesterol regulation pathways in terms of biochemical reactions, based on an extensive review of experimental and clinical literature. We are able to decompose the associated chemical reaction network structures into several independent subnetworks, one of which is responsible for conferring RPA on several intracellular forms of cholesterol. Remarkably, our analysis reveals that RPA in the cholesterol concentration in the endoplasmic reticulum (ER) is almost certainly due to a well-characterised control strategy known as antithetic integral control which, in this case, involves the high-affinity binding of a multi-molecular transcription factor complex with cholesterol molecules that are excluded from the ER membrane. Our model provides a detailed framework for exploring the necessary biochemical conditions for robust homeostatic control of essential and tightly regulated cellular molecules such as cholesterol.


Introduction
Cholesterol is an amphipathic lipid that plays an essential role in regulating the properties of cell membranes in mammalian cells.Specifically, cholesterol forms part of the assembly and function of lipid micro-environments on the cell surface, known as lipid rafts (Simons and Toomre, 2000;Ouweneel et al., 2020).Enriched with free cholesterol (FC) and glycosphingolipids, these dynamic micro-domains function to segregate and distribute lipids and proteins to the cell surface.Here, the lipid rafts act as a platform for cellular receptors through which cellular processes such as endocytosis, exocytosis, receptor trafficking and cell signalling can be realised.(Simons and Ikonen, 2000;Simons and Toomre, 2000).Signalling pathways that involve lipid rafts include immunoglobulin E (IgE) signalling during the allergic immune response (Sheets et al., 1999;Simons and Toomre, 2000), T-cell antigen receptor signalling (Janes et al., 2000;Simons and Toomre, 2000) and Hedgehog signalling (Incardona and Eaton, 2000;Simons and Toomre, 2000).
In contrast, the altered composition of lipid rafts (Kulkarni et al., 2022) is characteristic of neurodegenerative diseases such as Parkinson's, Alzheimer's, multiple sclerosis and lysosomal storage disease (Kulkarni et al., 2022), while pathogens such as HIV, for example, exploit lipid raft composition to enter and exit its host cell (Simons and Ikonen, 2000).Glycosphingolipids and cholesterol are central players in lipid raft biology (Sonnino and Prinetti, 2013;Grassi et al., 2020), and considering the impacts of disrupted lipid raft composition, the total cellular level of cholesterol and its distribution between membranes, and within a given membrane, must be precisely controlled (Simons and Ikonen, 2000).Furthermore, as a macronutrient, cholesterol plays an integral role in the synthesis of steroid hormones and vitamin D, while dysregulated cholesterol metabolism can lead to pathophysiological diseases such age-related macular degeneration and cardiovascular disease (Curcio et al., 2011;Mc Auley and Mooney, 2014;Scheepers et al., 2020).Recent experimental studies by Lu et al. (2022) have also highlighted the critical role that cholesterol homeostasis plays in terminal erythropoiesis, where disruption of intracellular cholesterol levels can interfere with terminal erythroid differentiation, ultimately leading to anaemia.
Evolution has endowed mammalian cells with intricate regulatory mechanisms to control cholesterol concentrations, ensuring adequate availability when required and shutting down cholesterol synthesis and uptake in conditions of cholesterol excess (Howe et al., 2016).The regulation of cholesterol takes place primarily in the plasma membrane (PM), endoplasmic reticulum (ER) and mitochondria, where the interplay between processes such as de novo biosynthesis, ingestion, efflux, storage and chemical modifications for specialised functions maintain steady cellular cholesterol levels (Tabas, 2002;Lange and Steck, 2016;Luo et al., 2020).Cholesterol regulation in the ER involves key proteins that respond to abundant cellular cholesterol levels by inhibiting cholesterol-promoting transcription pathways.Specifically, a family of transcription factors, collectively known as sterol regulatory element-binding proteins (SREBP), are sequestered by an ER-bound protein, SCAP, to prevent gene transcription of the cholesterol biosynthesis pathway as well as the exogenous uptake of cholesterol from the circulation and homeostasis of cholesterol is restored.Post translationally, oxygenated derivatives of excess cholesterol activate the reverse cholesterol removal pathways to remove excess cholesterol and maintain stable cholesterol levels (Howe et al., 2016).
But despite the growing literature on the many metabolic and signal transduction events involved in cholesterol regulation, and although some mathematical models of cholesterol metabolism have been formulated (Paalvast et al., 2015;Pool et al., 2018), there is currently no comprehensive mathematical description of the molecular network mechanisms that can tightly control intracellular cholesterol concentration within sub-cellular compartments, and maintain these concentrations within very tight concentration tolerances, in the face of significant variations in the delivery and intracellular uptake of cholesterol.
The homeostasis imposed on subcellular cholesterol concentrations corresponds to a stringent form of regulation known as Robust Perfect Adaptation (RPA) (Araujo and Liotta, 2018;Khammash, 2021;Araujo and Liotta, 2023b).RPA is a biological network property whereby one or more molecular concentrations within the network are maintained at a particular fixed "setpoint" at steady state, independently of specific network parameter choices (i.e., robustly), and in spite of altered inputs or disturbances to the system (Cannon, 1929;Kitano, 2007;Ang et al., 2010;Ferrell, 2016;Aoki et al., 2019;Jeynes-Smith and Araujo, 2021;Jeynes-Smith and Araujo, 2023).The mathematical principles underpinning RPA in complex biological networks are now well understood, and we now have access to a universal description of RPA at both the macroscale of biochemical reaction networks (Araujo and Liotta, 2018;Araujo and Liotta, 2023a), and at the network microscale of individual intermolecular interactions and chemical reactions (Araujo and Liotta, 2023b).In particular, it is now known that all RPA-capable networks, no matter how large or complex, are decomposable into two distinct and well-defined subnetwork modules-Opposer modules and Balancer modules (Araujo and Liotta, 2018;Araujo and Liotta, 2023a).Whereas Opposer modules have an overarching feedback structure, into which special types of interactions (known as "opposer nodes") are embedded, Balancer modules comprise a number of incoherent parallel pathways, incorporating specialised classes of interactions known as "balancer nodes" and "connector nodes".Particular network realisations of these two modular types can potentially be highly complex, and large and highly complicated RPA-capable networks can be constructed from the interconnections of these special RPA-conferring modular building blocks.We refer interested readers to (Araujo and Liotta, 2018;Khammash, 2021;Araujo and Liotta, 2023a) for a comprehensive overview of the macroscale topological properties of RPA-capable networks.
The microscale RPA problem, which considers the detailed chemical reaction networks (CRNs) that can implement RPA through potentially intricate intermolecular interactions, has also now been solved (Araujo and Liotta, 2023b).In particular, it is now apparent that networks of chemical reactions that can support RPA are subject to very stringent structural criteria, which allow topological features (opposer nodes, balancer nodes, connector nodes) of an associated RPA-module to be created by linear integral controllers.In this way, RPA-conferring intermolecular interactions at the network microscale can be decomposed into a collection of subsidiary RPA-problems, each obtained by a linear coordinate change, that collaborate to impose RPA on the concentrations of certain specific molecules in the CRN.
In the present study, we identify the molecular mechanisms responsible for conferring RPA on plasma membrane cholesterol.In particular, we identify (i) the topological characteristics of the RPAconferring network responsible for cholesterol homeostasis, (ii) the molecular mechanism by which cells sense a cholesterol concentration in excess of its association capacity with phospholipids in the membrane bilayers, and (iii) the integral that confers RPA on membrane-bound cellular cholesterol.In Section 2, we begin with an extensive review of the known molecular interactions and reactions involved in cholesterol regulation, including the critical transcriptional, translational and post-translational signalling events that regulate cholesterol within and between different subcellular compartments and organelles.In Section 3, we encode these detailed molecular mechanisms into a chemical reaction network (CRN), whose graphical structure and modularity can be analysed, and which can be used to identify molecular mechanisms directly responsible for RPA (Section 4).
2 The cell biology of cholesterol 2.1 Cholesterol regulation pathways

Cholesterol flux
The ER contains an elaborate feedback system that is sensitive to the level of its membrane cholesterol by virtue of its small pool of cholesterol, which is set by the balance of PM cholesterol fluxes moving rapidly to and from the ER (Radhakrishnan et al., 2008;Lange and Steck, 2016).While little cholesterol moves from the plasma membrane to the ER at homeostatic cholesterol levels, the flux increases when plasma membrane cholesterol levels rise to the point where active cholesterol molecules, those above the association capacity of the bilayer phospholipids, escape from the surface of the bilayer and redistributes to the ER (Lange and Steck, 2020).
A proportion of cholesterol in the ER is converted to cholesterol esters (CE) by the ER enzyme acyl-coenzyme A (CoA):cholesterol acyltransferase (ACAT) (Chang et al., 2006;Chang et al., 2009).Newly synthesised CE accumulates within the rough ER and buds off as cytoplasmic lipid droplets, whereby it acts as a reservoir from which its hydrolysis can mobilise unesterified cholesterol for steroid hormone, oxysterols or bile acid production (Simons and Ikonen, 2000;Chang et al., 2009;Scheepers et al., 2020).
Reverse cholesterol transport begins with the formation of small high-density lipoprotein (HDL) particles by the liver and intestine.HDL acts as unesterified cholesterol acceptor that transports excess cholesterol within peripheral tissues to the plasma, a process mediated by ABCA1.Mature HDL particles then transport the cholesterol to the liver where it can be directly excreted into the bile or be metabolised into bile acids/salts before excretion (Feingold and Grunfeld, 2000;Ouimet et al., 2019).
Intracellular cholesterol is provided to the cellular organelles via biosynthesis and the uptake of circulating low-density lipoprotein complexes (LDL) via lipoprotein receptors (LDLR) in the plasma membrane.While the balance between these two pathways depends on cell type and the availability of LDL-derived cholesterol, the master regulator of gene expression to realise these pathways is a polytopic ER membrane protein called SCAP (Brown et al., 2018).

SREBP transcription factors
SCAP is inherently bound to sterol regulatory element-binding proteins (SREBP), a family of transcription factors anchored in the ER.Under abundant cholesterol conditions, the complexed SCAP molecule binds to cholesterol and the hydrophobic ER protein Insig-1 (Gong et al., 2006;Brown et al., 2018).The conformational changes that SCAP undergoes during this binding process prevent incorporation of the complex into COPII-coated vesicles for transport to the Golgi for processing (Sun et al., 2007).Consequently, the activation of gene transcription pathways further downstream is inhibited, biosynthesis and exogenous cholesterol uptake repressed, and cholesterol homeostasis is restored (Gong et al., 2006;Howe et al., 2016;Brown et al., 2018).
While in vitro studies confirmed the inhibition of SREBP processing under abundant cholesterol conditions, the mechanism by which the cell senses this abundance in the ER in vivo is not known (Lange and Steck, 1997).Motamed et al.(Motamed et al., 2011) showed the cholesterol-binding site in Scap is located in a membrane-associated loop that projects into the lumen of ER.Recognising that the study of Scap's interactions with membrane cholesterol is technically difficult (Radhakrishnan et al., 2004), Gay et al. (Gay et al., 2015) used another sensor molecule, the soluble bacterial toxin perfringolysin O (PFO), to confirm the inhibition of SREBP processing and gene expression as reported by (Radhakrishnan et al., 2008).Based on this evidence, it was hypothesised that SCAP may be binding to a pool of ER cholesterol that exceeds the sequestration capacity of the bilayer phospholipids and projects into the ER lumen, the exact location where the cholesterol-sensing domain of SCAP resides (Motamed et al., 2011;Gay et al., 2015).
Under cholesterol-depleted conditions, Insig-1 dissociates from the SREBP/Scap complex to be ubiquitylated and degraded (Gong et al., 2006).The Scap/SREBP complex is now free to be carried to the Golgi complex, where fusion with the Golgi membrane results in proteolysis of SREBPs to release its transcription factors (TF).This nuclear fragment is translocated to the nucleus to, amongst many other activities, activate genes for the expression of (i) enzyme 3-hydroxy-3-methylglutaryl-coenzyme A reductase (HMGCR), (ii) low-density lipoprotein receptor (LDLR) on the plasma membrane, and (iii) proprotein convertase subtilisin/kexin type-9 (PCSK9) (Lagace, 2014;Howe et al., 2016;Brown et al., 2018).

Cholesterol biosynthesis and down-regulation
The cholesterol biosynthesis process involves multiple complex pathways tightly controlled at various points, with the HMGCR enzyme catalysing the synthesis of mevalonate from 3-hydroxy-3-methylglutaryl-coenzyme A (HMG-CoA) in the first committed and rate-limiting step.This reaction forms the point of feedback control for the mevalonate pathway, and, as such, this enzyme is the primary target for medical intervention.Following a succession of intermediate changes, mevalonate is converted to lanosterol, the latter being converted to cholesterol after several sequential reactions.For the interested reader, Cerqueira et al. (Cerqueira et al., 2016) provides a detailed review of the biosynthesis of cholesterol.
The transcriptional down-regulation of cholesterol via the SREBP pathway is relatively slow, with the mRNA of target genes decreasing only after several hours.In contrast, 25-hydroxycholesterol (25HC) is an oxysterol produced by enzymatic conversion of ER cholesterol and is sensed by the membrane domain of HMG-CoA reductase to trigger its binding to Insig.Within minutes, this binding activates the enzyme's ubiquitylation, degradation, and subsequent inhibition of cholesterol biosynthesis (DeBose-Boyd, 2008;Sharpe and Brown, 2013;Odnoshivkina et al., 2022).

Receptor-mediated LDL acquisition
LDLR on the plasma membrane facilitate endocytosis of LDL from the circulation.Endosomes release their content to lysosomes, where degradation of internalised LDL leads to the release of unesterified cholesterol and fatty acids (Brown and Goldstein, 1986).Recent studies by Trinh et al. (2020) confirmed that LDL-derived cholesterol from lysosomes are transported to the PM first to maintain optimal cholesterol levels, and subsequently from PM to the regulatory domain of the ER.
Transcriptional down-regulation of LDLR is facilitated by PCSK9, a soluble member of the proprotein convertase family of secretory serine endoproteases.Upon binding to LDLR on the cell surface, the complex is redirected for internalisation into lysosomes, where PCSK9 and LDLR are degraded.Consequently, LDLR recycling to the cell surface is prevented and LDL uptake from the circulation down-regulated (Lagace, 2014).
As an essential structural component of cell membranes, cholesterol modulates the properties of the membrane lipid bilayers by intercalating with phospholipids to maintain membrane fluidity, and rigidity (Lange and Steck, 2020).The association of cholesterol with diverse organelle phospholipids in specific ratios in cell membranes has been described by McConnell and Radhakrishnan (McConnell and Radhakrishnan, 2003) as stoichiometric complexes.Association equilibrium is achieved when cholesterol accumulates as phospholipid complexes up to a particular stoichiometry or set-point (Yeagle and Young, 1986;Lange et al., 2013).Conversely, uncomplexed sterol molecules, those above the association capacity of the bilayer phospholipids, have a relatively high tendency to escape from the surface of the bilayer.This tendency is also referred to as a constituent's activity or fugacity (Lange and Steck, 2008) and this fraction of cholesterol is called active cholesterol (Lange and Steck, 2016).These active molecules rapidly redistribute to the cytoplasmic membranes down its activity gradient via Aster proteins (Sandhu et al., 2018), where a rise in intracellular pools adjusts the synthesis and uptake of cholesterol to return the plasma membrane cholesterol to its physiological set-point (Lange and Steck, 2008;2016).

Liver X (LXR) transcription factors
27-hydroxycholesterol (27HC), another oxysterol derivative of ER cholesterol, is a transcriptional activator of a member of the nuclear receptor superfamily, the liver X receptor (LXR).This subfamily of transcription factors consists of two distinct members, LXRβ, being ubiquitously expressed, and LXRα, whose expression is restricted to tissues rich in lipid metabolism Lu et al. (2001).Being a Type-II receptor, LXR resides in the nucleus and is bound to its specific DNA response elements, even in the absence of ligand.27HC mediates signal transduction through direct binding to LXR in the nucleus at physiological concentrations (Janowski et al., 1999;Sever and Glass, 2013;Howe et al., 2016).Once activated, LXR upregulates ABCA1genes to facilitate cholesterol efflux via systemic HDL particles.Furthermore, activated LXR upregulates the inducible degrader of LDLR (Idol), to decrease cholesterol uptake from LDL (Zelcer et al., 2009;Howe et al., 2016).
We summarise all abbreviations used in the above descriptions in Table 1.

Chemical reaction network representation of cellular cholesterol regulation
Here we represent the biochemical processes reviewed in detail in Section 2 as a collection of chemical reactions, which gives rise to a CRN involving fourteen molecular species, comprising transcription factors, proteins, biochemical intermediates and macro-molecules (see Table 2).These key molecules participate in twenty-two chemical reactions (see Table 3).One of the species-the (extracellular) concentration of cholesterol-carrying lipoproteins (C L )serves as the distinguished input molecule for this system, since this determines the delivery of cholesterol to the cell, and is therefore the 'stimulus' to which intracellular cholesterol must adapt and exhibit homeostasis/RPA.
The biochemical processes captured by the individual reactions of the CRN are as follows: We capture the constitutive production rate of the SREBP/ SCAP/Insig1 complex in the ER by the reaction ∅→ μ S ci .
Under cholesterol-depleted conditions, S ci complexes are released to the Golgi for processing (dissociation and degradation of Insig-1 are omitted from our model).Proteolytic activity in the Golgi ensures the release of the SREBP transcription factors S r , S h and S p to the nucleus, expressed here as S ci → k1 S ci + S r + S h + S p .
The rate at which the transcription factor S r activates transcription of its target genes to produce LDLR protein is captured as S r → p1 S r + R. Similarly, S h induces the production of HMGCR, S h → p2 S h + H r , while the production of PCSK9 from S p transcription is represented as S p → p3 S p + P. Degradation rates of the SREBP transcription factors are expressed as S r → k2 ∅, S h → k10 ∅ and S p → k13 ∅ respectively (Sundqvist and Ericsson, 2003).
Receptor-mediated uptake and release of unesterified cholesterol via endosomes are captured in the reaction Here, the conservation of R represents the recycling of LDLR to the surface membrane.The interested reader is referred to (Tindall et al., 2009;Pool et al., 2018)  Returning to the cholesterol biosynthesis pathway, the rate of constitutive production of HMGCoA (H) is expressed as ∅ → α H.
Following a series of omitted intermediate steps, the final product in the biosynthesis pathway, C e , is produced via activation of HMGCoA by the HMGCR reductase enzyme, modelled as While cholesterol biosynthesis is tightly regulated via various feedback processes that lead to the degradation of HMGCR, we represent this cholesterol initiated degradation process as ER cholesterol molecules, C, exceeding the sequestration capacity of the phospholipid bilayer-that is, the active cholesterol molecules-are released into the ER lumen, modelled here by the reaction C e → θ C e + C.
The cholesterol molecules in the ER lumen are now sensed and bound by the S ci complex, resulting in the irreversible retention of a S ci C complex in the ER membrane, modelled as S ci + C→ η ∅.In this case, transport of the S ci complex to the Golgi for release of SREBP transcription factors is inhibited, resulting in a decrease in cellular cholesterol concentration.Finally, ER cholesterol efflux via systemic high-density lipoprotein (HDL) particles is represented as C e → k12 ∅.
In Figure 1, we provide a schematic overview of these interactions to summarise the overarching network structure associated with the regulation of cholesterol.We also provide a summary of symbols used to denote the various key molecules appearing in these interactions in Table 2, and a summary of the reactions themselves in Table 3.Now, a graph may be obtained from a set of chemical reactions by interpreting the multisets of species constituting either reactants or products [known collectively as "complexes" in chemical reaction network theory (CRNT)] as vertices, and the reactions themselves as directed edges.Connected components of the graph are referred to as "linkages classes" in CRNT (Feinberg, 2019).A "strong linkage class" corresponds to a strongly connected component (SCC) of the graph, being a maximal strongly-connected subgraph.A "terminal strong linkage class" is a strong linkage class in which no complex reacts to a complex in a different strong linkage class.Complexes belonging to terminal strong linkage classes are referred to as terminal complexes; all other complexes are non-terminal complexes.
For any such CRN graph structure, a key integer invariant known as the deficiency of the CRN may be computed, providing a quantitative measure of the linear independence of the CRN reactions given their distribution into linkage classes (Feinberg, 2019;Araujo and Liotta, 2023b).
The deficiency, δ, of a CRN is given by where m is the number of complexes, l is the number of linkage classes, and s is the dimension of the stoichiometric subspace (the span of the reaction vectors), also known as the rank of the CRN.In Figure 2 we organise the twenty-two reactions of the cholesterol homeostasis CRN into linkage classes.For this network, m = 23, l = 2 and s = 15, which gives a deficiency δ = 6.

Graph analysis of the CRN
By decomposing this CRN graph into algebraically independent subnetworks, the CRN reactions can be partitioned into subsets whose steady states may be determined independently from the rest of the CRN.In the context of RPA, such a decomposition into independent subnetworks allows us to distinguish the network reactions that contribute to the network's RPA capacity from reactions that play no role in the RPA capacity of the CRN (Araujo and Liotta, 2023b).Subnetworks are algebraically independent when their individual ranks sum to the rank of the parent network (Feinberg, 2019).Here we decompose the full CRN into two algebraically independent subnetworks, as shown in Figures 3, 4, where the sum of the ranks of the two subnetworks (s 1 = 5 and s 2 = 10 respectively) gives the rank of the full CRN (s = 15).Furthermore, the sum of the respective subnetwork deficiencies, δ 1 = 1 and δ 2 = 5 gives the full CRN deficiency of δ = 6.
Subnetwork 1 contains six non-terminal complexes C e , C p , C f , E, S ci and ∅.CRNs with a deficiency of precisely one can, by the Shinar-Feinberg Theorem (Shinar and Feinberg, 2010), achieve RPA if the CRN contains two distinct non-terminal complexes that differ in a single species.Since ∅ is a non-terminal complex in this algebraically-independent deficiency-one subnetwork, all nonterminal complexes comprising a single species (C e , C p , C f and E) are guaranteed to be RPA-capable.
We further decompose Subnetwork 1 into Subnetworks 1A and 1B, where the deficiency of Subnetwork 1A and Subnetwork 1B is δ 1A = 1 and δ 1B = 0 respectively.Subnetwork 1A contains the species involved in regulating the transcription of cholesterol-promoting processes via "entrapment" of the S ci C complex in the ER, including the non-terminal complex C e .Applying the Shinar-Feinberg Theorem (Shinar and Feinberg, 2010) to Subnetwork 1A, the two non-terminal complexes C e and ∅ differ in the single species C e only, identifying C e as the species capable of conferring RPA to the cellular cholesterol network.Subnetwork 1B contains the other three non-terminal complexes, C p , C f and E.
The reactions contained in Subnetwork 2 contribute to the overarching control of cellular cholesterol but do not have any influence on the RPA capacity of the CRN as a whole, and is therefore excluded from further analysis.
Having identified the controller reactions and the RPA-capable species of the cellular cholesterol network, we now proceed to determine the linear combination of rate equations that identify the integrator responsible for conferring RPA on the system, as well as the respective set points of each individual RPA species.

Mass action analysis
CRNs can also induce a set of rate equations under the mass-action assumption, whereby each reaction proceeds at a rate proportional to the concentration of each reactant.There now exist methods to analyse such polynomial dynamical systems to systematically extract the integral-control implementing features of RPA-capable CRNs (Araujo and Liotta, 2023b).For the cholesterol homeostasis CRN under consideration here, the mass-action rate equations are: Theorem 1 in (Araujo and Liotta, 2023b) makes clear that for any RPA-capable CRN with interacting molecules x 1 , . .., x n and corresponding mass-action rate equations f 1 , . .., f n , there always exists a linear combination of polynomials {r where ρ = g(x i , x j )(x i − c), in its lowest form, is the RPA polynomial of the CRN, being a special polynomial function in two variables.When this occurs, the CRN has the capacity for RPA in the variable x i , with setpoint x i = c, while x j is any variable that does not exhibit RPA.Moreover, if x i = c is a stable steady-state of the CRN's mass-action equations, the CRN exhibits RPA in the variable x i .In this case, the variables of the model are the fourteen species.Thus, we seek a projection of the ideal 〈f 1 , . .., f 13 〉 onto two species.We expect C e to be RPA-capable for this model, based on the analysis in the preceding section.Moreover, the input to the system, C L , is not RPA-capable by definition.We therefore compute the elimination ideal 〈f 1 , . . ., f 13 〉 ∩ R[C e , C L ] via computation of a suitable Gröbner basis.We provide a simple code to make this computation, which implements an efficient block elimination order, using the open-source software Singular (https://www.singular.uni-kl.de/) within the interactive Jupyter environment.Our code is available at https://github.com/RonelScheepers/RobustPerfectAdaptation.
Using this method, we compute the two-variable elimination ideal to be We note that the right-hand side of Eq 15 has the form of an RPA polynomial that is zero-order in the non-RPA-capable variable, C L , thereby confirming that the system is RPA-capable in the species C e , and also calculates that its steady-state setpoint is C e = μ/θ.Setpoints for the variables C p , C f and E can be computed similarly: k4θ and E k7μ k8θ , respectively.These calculations also confirm RPA capacity in all four variables suggested by the graphical analysis.Moreover, the linear combination of rate equations used to identify the RPA polynomial (15) within the steady-state ideal can be obtained directly using the lift function in Singular.In this case, we compute that This calculation reveals the linear coordinate change that identifies an "internal model" ( (Francis and Wonham, 1975;Francis and Wonham, 1976;Sontag, 2003;Bin et al., 2022;Araujo and Liotta, 2023b)), and demonstrates that the variable (C − S ci ) is the integral variable (or "actuator") for the integral controller.Other mathematical approaches tailored to the identification of linear controllers could also be used here to identify the integral variable and setpoint for the RPA-capable species (Cappelletti et al., 2020;Gupta and Khammash, 2022), given the simplicity of this particular controller.In addition, by computing the syzygies for the ideal, we confirm that this system has no linear syzygies.Thus, there are no mass conservation relationships embedded within the rate equations, and this linear coordinate change ( 16) is unique.We provide full details of these methods, along with an extended presentation of our analysis in our accompanying Supplementary Material.In Figure 5, we demonstrate via numerical simulation of the Eqs 1-13 that the species C e , C f , C p and E do indeed exhibit RPA for the indicated choice of parameters, and return to their expected setpoints for a range of step increases in the system input, C L .By contrast, none of the remaining nine species exhibits RPA, instead achieving a steady-state value that varies with C L .To confirm that these stable steady-state responses were achievable for a significant region of parameter space, we used Matlab's symbolic toolbox to compute the characteristic polynomial for the system's Jacobian matrix, and solved for the system's eigenvalues using 10,000 different random parameter sets, with individual parameters selected from a uniform distribution on the interval (1,50).This process was repeated for ten trials; for each trial, roughly 50% of parameter sets produced a stable system, with all calculated eigenvalues lying strictly in the left-half of the complex plane.This confirms that the CRN is not dependent on special fine-tuning of parameters to achieve stability, and can exhibit RPA in the noted variables for a wide range of different parameters.All Matlab code for this stability analysis is included at https://github.com/RonelScheepers/RobustPerfectAdaptation.

Discussion
The quest for the fundamental design principles that govern the robust implementation of important biological functions, particularly within the highly complex molecular interaction networks within single cells, is considered to be one of the most important grand challenges in the life sciences (Araujo and Liotta, 2023a;b).Robust Perfect Adaptation (RPA) is a keystone signalling phenomenon that is ubiquitously observed at all scales of biological organisation in contexts as diverse as the regulation of cellular signal transduction (Ferrell, 2016), sensory adaptation (Kaupp, 2010), single-celled chemotaxis (Alon et al., 1999), cellular stress responses (Eisner et al., 2018), in addition to playing a critical role in robust patterning during organism development (Ben-Zvi and Barkai, 2010).RPA is currently the only robust biological functionality for which there now exists a universal solution, at both the network macroscale (Araujo and Liotta, 2018) and the network microscale (Araujo and Liotta, 2023b).This universal solution makes clear that all RPA-capable networks are necessarily modular by design, and decomposable into subnetworks ('modules') from two and only two distinct classes-Opposer modules, with an overarching feedback structure; and Balancer modules, with an overarching feedforward structure.Moreover, these structures are characterised by well-defined topological features: opposer nodes (in Opposer modules), and balancer nodes and connector nodes (in Balancer modules).It is now known that, at the molecular level of biochemical reactions, all of these distinctive topological features (nodes) must be identifiable through linear coordinate changes in FIGURE 4 Subnetwork 2. Reactions in this Subnetwork contribute to the overarching controlling module but do not have any influence on the RPA capacity of the CRN as a whole.C e is the only species common to both Independent Subnetworks, and thus connects independent CRN subnetworks 1 and 2.
order for a chemical reaction network to exhibit RPA.In this way, control systems generated via intermolecular interactions and enzyme-catalysed reactions implement a special form of integral control that differs in important ways from modern engineering control systems that employ specially-designed components for the computation of integrals.
Although cellular cholesterol is thought to be subject to a very stringent form of homeostatic control, and thus RPA, the detailed   molecular mechanisms responsible for such robust control have been unknown until now, owing largely to the relative complexity of the signalling network underpinning cholesterol regulation in comparison with most known RPA-capable chemical reaction networks (CRNs).In this study, we develop a CRN for cholesterol regulation that encompasses a comprehensive review of known molecular interactions involved in the regulation of intracellular cholesterol concentrations.This has allowed us to undertake a detailed analysis of the signalling architectures that are involved in the tight regulation of cholesterol concentrations within cellular and sub-cellular membranes, from both a graphical analysis of the CRN structure, as well as an algebro-geometric analysis of the corresponding mass-action equations.
In common with the elegant recent mathematical model by Steck et al. (2021), our model captures the overarching negative feedback structure of the cellular cholesterol regulatory network.In fact, our analysis makes clear that the tight homeostatic control of membrane-bound cholesterol, in both the cellular plasma membrane and in the endoplasmic reticulum (ER), is the result of a single Opposer module, containing a single opposer node.We summarise the architecture of this network schematically in Figure 6, highlighting the overarching feedback architecture of the network.
Our graphical analysis suggests that it is the concentration of cholesterol in the ER membrane (C e ) that is subject to a molecular control mechanism, and thereby acts as a "sensor" molecule (see Subnetwork 1A).Once RPA is conferred to this specific subcellular cholesterol concentration, RPA is then also conferred to the cholesterol concentration at the cellular plasma membrane (C p ), to the cholesterol concentration released form the LDLR (C f ) and to the esterified cholesterol concentration (E), by virtue of the relationships of these molecules to C e in the CRN graph structure (see Subnetwork 1A).Molecular concentrations in the remaining independent subnetwork of the CRN (see Subnetwork 2) are not subject to RPA.
By computing a projection of the ideal generated by the mass-action equations onto two variables-the input to the network, C L , and the RPA-capable sensor molecule, C e -we demonstrate that the simple linear coordinate transformation that extracts an RPA polynomial from the system is unique, and corresponds to a known and well-characterised molecular control mechanism known as antithetic integral control (Briat et al., 2016;Briat et al., 2018).This simple RPA-conferring control mechanism has previously been identified in the form of sigma/anti-sigma factors in bacteria, and is now the subject of intense theoretical interest, having recently been implemented by synthetic integral feedback circuits within living cells (Aoki et al., 2019;Khammash, 2021).It is also now clear that antithetic integral control is an example from a special class of molecular control strategies known as maxRPA(Gupta and Khammash, 2022), which all involve setpoints determined by a minimum number of parameters.Our study now reveals that antithetic integral control is almost certainly the basis for the exquisitely controlled concentration of cholesterol in cellular and subcellular membranes.In this case, the counterpart to the interaction between sigma and anti-sigma factors is the sequestration of active cholesterol (C) by the SREBP/Scap/ Insig complex (S ci ).The linear combination of reaction rates involving these two variables gives rise to the single opposer node of the Opposer module, which we highlight in greater detail in Figure 7.Our analysis is also able to explicitly determine all setpoints for RPA-capable molecules as functions of network parameters.
The modelling and analysis approach presented here demonstrates how the fundamental network design principles required for robust homeostatic control can be realised within a relatively complex cellular regulatory network.Our study provides a detailed framework for exploring the necessary biochemical conditions for robust homeostasis and adaptation in other tightly regulated molecules, or sensory molecules, in complex networks throughout biology.Scheepers and Araujo 10.3389/fcell.2023.1244297

FIGURE 1
FIGURE 1 Schematic overview of molecular interactions associated with cholesterol regulation.ER cholesterol molecules (C e ) that exceed the sequestration capacity of membrane phospholipids protrude into the ER lumen (C).These are sensed and sequestered by constitutively produced SREBP/Scap/Insig complex, S ci (solid red line).Unbound S ci releases three transcription factors (S R , S H , S p ) into the nucleus to induce expression of LDLR (R), HMGCR (H R ) and PCSK9 (P).Cholesterol delivered from blood-borne lipoproteins (C L ) is taken up by LDLR, then engulfed and released via lysosomes to become free unesterified cholesterol (C f ) that supplements the PM cholesterol pool (C p ). Fluxes between PM cholesterol and ER membrane cholesterol (C e ) establish a balance, while esterification of C e into E further contribute to the overall regulation of cellular cholesterol concentration.PCSK9 binding to LDLR initiates degradation of both molecules via internalisation into lysosomes.Biosynthesis of cholesterol starts with the constitutive formation of HMG-CoA (H), to eventually be converted to C e , a process initialised by the catalytic enzyme H R .The rate-limiting step in endogenous production of cholesterol occurs when C e inhibits catalytic activity via degradation of H R .Black arrows indicate flux or upregulation, dashed arrows indicate de novo protein synthesis, blunt arrow heads indicate inhibition, round arrow head indicates catalytic target.Created with BioRender.com.

FIGURE 2
FIGURE 2CRN graph structure for the cholesterol homeostasis network.

FIGURE 3
FIGURE 3Subnetwork 1, with deficiency 1, confers RPA on the species C e , C p , C f and E. Subnetwork 1 can be decomposed into algebraically independent Subnetworks 1A and 1B respectively.Subnetwork 1A contains species C e , that regulates the transcription of cholesterol promoting processes, while Subnetwork 1B comprise the other three non-terminal complexes, C p , C f and E.

FIGURE 6
FIGURE 6Network schematic for the cellular cholesterol regulatory network.This network diagram captures the nature of the interactions among the thirteen species of the CRN, illustrating the overarching topology of the cellular cholesterol homeostatic machinery as a single Opposer module, with its characteristic feedback architecture.Grey arrows represent flux, or an activating/upregulating influence, while blunt arrow heads (flat bars) represent inhibition or a negative/downregulating influence.Solid circular line endings represent a catalytic reaction.The single opposer node, comprising the irreversible sequestration of cholesterol within the SREBP/Scap/Cholesterol complex (S ci C) is represented by the green arrows.

FIGURE 7
FIGURE 7Single opposer node for the cellular cholesterol regulatory network.The purple box highlights the sequestration process that implements a form of antithetic integral control, thereby constituting a single opposer node via a single linear coordinate change.These opposer interactions (indicated in red) confer RPA on the sensor molecule, C e , which in turn imparts RPA to several upstream molecules (not shown).

TABLE 1
List of molecular abbreviations.

TABLE 2
List of species and associated symbols.

TABLE 3
Full summary of all reactions comprising cellular cholesterol regulation network.
ptk12 ∅ Frontiers in Cell and Developmental Biology frontiersin.org