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 c i .

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 c i → k 1 S c i + 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 → p 1 S r + R . Similarly, S _h induces the production of HMGCR, S h → p 2 S h + H r , while the production of PCSK9 from S _p transcription is represented as S p → p 3 S p + P . Degradation rates of the SREBP transcription factors are expressed as S r → k 2 ∅ , S h → k 10 ∅ and S p → k 13 ∅ respectively (Sundqvist and Ericsson, 2003).

Receptor–mediated uptake and release of unesterified cholesterol via endosomes are captured in the reaction R + C L → k 3 R + C f . 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) for comprehensive modelling of this recycling process. Post–translational regulation of LDLR occurs via IDOL mediated degradation, simply captured as R → k 11 ∅ , while the binding of PCSK9 to LDLR facilitates the degradation of LDLR as well as its own degradation in the lysosomes. These degradation rates are expressed by the equations P + R → k 14 P , and P → k 16 ∅ , respectively.

C _f is transported to the PM, captured by the rate equation C f → k 4 C p . From here, fluxes between cholesterol in the PM, C _p , and the ER membrane cholesterol pool (C _e ), are captured as C p ⇌ k 6 k 5 C e .

A proportion of ER cholesterol is converted to cholesterol esters at a rate, C e → k 7 E , while a proportion of E can be converted back to unesterified ER cholesterol when required at a rate E → k 8 C e .

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 H R + H → k 9 H R + C e .

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 C e + H R → k 15 C e .

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 c i + 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 → k 12 ∅ .

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 δ = m − l − s , 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.

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 ₁ = 5 and s ₂ = 10 respectively) gives the rank of the full CRN (s = 15). Furthermore, the sum of the respective subnetwork deficiencies, δ ₁ = 1 and δ ₂ = 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 non-terminal 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.

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: f 1 = d S c i d t = μ − η S c i C , (1) f 2 = d C d t = θ C e − η S c i C , (2) f 3 = d S r d t = k 1 S c i − k 2 S r , (3) f 4 = d S h d t = k 1 S c i − k 10 S h , (4) f 5 = d R d t = p 1 S r − k 11 R − k 14 P R , (5) f 6 = d C f d t = k 3 C L R − k 4 C f , (6) f 7 = d C p d t = k 4 C f − k 5 C p + k 6 C e , (7) f 8 = d C e d t = k 5 C p − k 6 C e + k 9 H R H − k 12 C e − k 7 C e + k 8 E , (8) f 9 = d E d t = k 7 C e − k 8 E , (9) f 10 = d H R d t = p 2 S h − k 15 H R C e , (10) f 11 = d H d t = α − k 9 H R H , (11) f 12 = d S p d t = k 1 S c i − k 13 S p , (12) f 13 = d P d t = p 3 S p − k 16 P , (13)

Theorem 1 in (Araujo and Liotta, 2023b) makes clear that for any RPA–capable CRN with interacting molecules x ₁, …, x _n and corresponding mass–action rate equations f ₁, …, f _n , there always exists a linear combination of polynomials { r 1 , … , r n } ⊂ R [ x 1 , … , x n ] such that r 1 f 1 + ⋯ + r n f n = g x i , x j x i − c = ρ , (14) 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 ₁, …, f ₁₃⟩ 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 〈 f 1 , … , f 13 〉 ∩ R C e , C L = 〈 θ C e − μ 〉 . (15) 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: C p = ( k 6 + k 12 ) μ − θ α k 5 θ , C f = k 12 μ − θ α k 4 θ and E = k 7 μ k 8 θ , 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 d C d t − d S c i d t = θ C e − μ . (16) 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.