The Potential Roles of Transacylation in Intracellular Lipolysis and Related Qssa Approximations

Fatty acids (FAs) are crucial energy metabolites, signalling molecules, and membrane building blocks for a wide range of organisms. Adipose triglyceride lipase (ATGL) is the first and presumingly most crucial regulator of FA release from triacylglycerols (TGs) stored within cytosolic lipid droplets. However, besides the function of releasing FAs by hydrolysing TGs into diacylglycerols (DGs), ATGL also promotes the transacylation reaction of two DG molecules into one TG and one monoacylglycerol molecule. To date, it is unknown whether DG transacylation is a coincidental byproduct of ATGL-mediated lipolysis or whether it is physiologically relevant. Experimental evidence is scarce since both, hydrolysis and transacylation, rely on the same active site of ATGL and always occur in parallel in an ensemble of molecules. This paper illustrates the potential roles of transacylation. It shows that, depending on the kinetic parameters but also on the state of the hydrolytic machinery, transacylation can increase or decrease downstream products up to 80% respectively 30%. We provide an extensive asymptotic analysis including quasi-steady-state approximations (QSSA) with higher order correction terms and provide numerical simulation. We also argue that when assessing the validity of QSSAs one should include parameter sensitivity derivatives. Our results suggest that the transacylation function of ATGL is of biological relevance by providing feedback options and altogether stability to the lipolytic machinery in adipocytes.


Introduction
Fatty acids (FAs) are crucial for ATP production, synthesis of biological membranes, thermogenesis, and signal transduction (Glatz and Luiken 2015;Zechner et al. 2012). Animals and humans store FAs in the form of water-insoluble triacylglycerols (TGs, one species thereof is triolein, TO) within cytosolic lipid droplets (LDs) of specialised fat cells called adipocytes, but also in other cell types. From a biochemical perspective, a TG molecule is composed of three FAs esterified to a glycerol backbone. LD-associated TGs are degraded step-wise by an enzymatic cascade of lipid hydrolases commonly designated as "lipases" in a process called intracellular lipolysis.
Distinct lipases are responsible for each lipolytic step (Grabner et al. 2021;Young and Zechner 2013). The initial step in the degradation of TGs is catalysed by adipose triglyceride lipase (ATGL) releasing one FA and one diacylglycerol (DG, one species thereof is diolein, DO) molecule (Schreiber et al. 2019). In the second lipolytic step, hormone-sensitive lipase (HSL) hydrolyses DGs into FAs and monoacylglycerols (MGs; one species thereof is monoolein, MO) (Recazens et al. 2021). Finally, monoglyceride lipase (MGL) hydrolyses MGs into FAs and glycerol, both of which can be secreted from the adipocytes into the bloodstream to supply remote organs.
To ensure an adequate response to changes in the metabolic state of an organism, lipolysis is fine-tuned by a large number of molecular factors regulating the three lipases (Hofer et al. 2020). One of these factors is comparative gene identification-58 (CGI-58) which accelerates lipolysis by increasing ATGL activity . Despite literally thousands of studies on the physiological role of lipolysis and its complex regulation, a quantitative description of the intricate reaction processes of TG hydrolysis is, however, still missing.
Given the outstanding metabolic importance of lipolysis, it is remarkable that ATGL (and possibly HSL (Zhang et al. 2019)) not only catalyse the hydrolysis reaction but in parallel also a transacylation reaction, which transfers a FA from one DG onto a second DG to generate TG and MG (Jenkins et al. 2004;Zhang et al. 2019;Kulminskaya et al. 2021;Patel et al. 2022). In fact, measurement of ATGL activity is far from being trivial due its hydrolase and transacylase activity occurring simultaneously. Chemical inhibition of ATGL (Kulminskaya et al. 2021) as well as mutations in its active site (Patel et al. 2022) abolish both activities, implying that it is experimentally not possible to separate one activity from the other. For example, as soon as a TG molecule formed by DG transacylation is hydrolysed to DG, then the effective outcome (one FA and MG were formed) is identical to the outcome of a DG hydrolysis event, thereby masking DG transacylase activity.
The importance and potential metabolic function of transacylation is unknown. By its function, the transacylation reaction provides feedback towards the TG pool as well as an alternative in creating MG alongside DG hydrolysis. One may believe that transacylation plays a minor to no role in TG metabolism. However, in vitro experiments from our laboratory (Fig. 1) and others (Jenkins et al. 2004;Kulminskaya et al. 2021) showed remarkable ATGL-mediated DG transacylation.
The experimental procedure behind obtaining Fig. 1 uses a well-defined substrate of artificial LDs made of DO (a particular DG) and purified ATGL as sole lipase. More precisely, a truncated variant of mouse ATGL (mAT288) was used since ATGL is biochemically challenging enzyme and removing a functionally non-essential part of ATGL helps its purification. Figure 1 shows that when mAT288 was incubated with an emulsion of DO LDs, significant amounts of TO were formed both in absence or presence of the ATGL endogenous co-activator CGI-58. As mentioned above, it is impossible to directly distinguish DO transacylation and DO hydrolysis in Fig. 1, since TOs produced by transacylation are later hydrolysed by ATGL to DOs and FAs.
In this paper we study the potential roles of transacylation in relation to the hydrolytic pathway TG → DG → MG as sketched in Fig. 2. More precisely, we consider three activities: TG hydrolysis, DG hydrolysis and DG transacylation. We Fig. 1 Catalytic activities of ATGL towards DGs. Truncated, purified mouse ATGL covering amino acids 1-288 (mAT288) was incubated with a phospholipid-emulsified, 0.3 mM DO (a particular DG) substrate in the absence and presence of its physiological co-activator CGI-58. Incubation of the DO substrate with an equal volume of buffer served as negative control (blank). At the given time points, reactions were stopped and total lipids were extracted with CHCl 3 :MeOH = 2:1. Lipids were resolved by thin layer chromatography using the solvent CHCl 3 :acetone:acetic acid = 88:12:1 and visualised by charcoaling at 140 • C will not include the final step of lipolysis, in which MGL hydrolyses MG to glycerol and FA, since this step is purely downstream and does not interact with the question of understanding the role of DG transacylation.
TG hydrolysis is primarily catalyzed by ATGL showing robust TG hydrolase activity especially in presence of CGI-58 . Accordingly, ATGL ) as well as CGI-58 knockout mice (Radner et al. 2010) accumulate TGs to supraphysiological levels in various tissues. However, ATGL can also act as DG transacylase (Jenkins et al. 2004;Zhang et al. 2019;Kulminskaya et al. 2021) and DG hydrolase (Zimmermann et al. 2004), indicating that it is also involved in the second lipolytic step. Nevertheless, DG hydrolysis is believed to be primarily catalysed by HSL. Adding to the complexity of lipolysis, HSL also exhibits minor but empirically proven TG hydrolase activity (Fredrikson et al. 1981).
For the sake of simplicity and clarity, we will not detail the activities of ATGL and HSL, whose relative contributions are also not sufficiently well understood. Instead, we will combine the functions together into single reaction rates for TG hydrolysis, DG transacylation and DG hydrolysis, see Fig. 2. Specifically, we postulate Michaelis-Menten (MM) kinetics for TG and DG hydrolysis and mass action law kinetics (MAL) for DG transacylation.
However, both approaches have drawbacks. Firstly, choosing MM kinetics should not be considered as an ultimate modelling choice: Despite the fact that MM kinetic rates for TG and DG hydrolysis can be found in the literature (see e.g. Viertlmayr 2018; Kim et al. 2008 and the references therein) already the mere fact that transacylation occurs in parallel to hydrolysis implies that hydrolysis can at best be approximatively described by MM kinetics. In favour of our modelling choice, the data in Viertlmayr (2018) also show that MM kinetics for TG hydrolysis is still a decent approximation of the real, more complicated dynamics.
Secondly, using MAL as reaction rate model for an enzymatic transacylation reaction is also questionable. However, we compared parameter fittings for preliminary in vitro experimental data using mAT288 incubated with a phospholipid-emulsified TG substrate and found that using MM or Hill kinetics instead of MAL didn't improve the (already good) fitting quality obtain with the MAL model, and despite involving two resp. three fitting parameters instead of one. Hence, Occam's razor compels us to use a MAL transacylation rate, at least until the molecular mechanism of transacylation becomes better understood.
From a mathematical viewpoint, the stoichiometry of DG hydrolysis (involving one DG molecule) is different to DG transacylation (involving two DG molecules), which implies different nonlinearities of the respective reaction rates. Hence, the relative contributions of DG hydrolysis and DG transacylation depend on the DG concentration, which in return depends on the kinetics producing and processing DGs: We are faced with a complex nonlinear system behaviour. Finally, we also emphasise that we believe that the analysis of this paper can be extended to general superlinear and monotone increasing models for DG transacylation rates. Let where V 1 , K 1 and V 2 , K 2 are the corresponding maximum velocities and Michaelis constants of the TG and DG hydrolysis and σ denotes the MAL constant for the transacylation 2DG → TG + MG. Note that it is plausible to assume constant concentrations of ATGL and HSL. Hence, there is no explicit dependence of the parameter V 1 , σ and V 2 on ATGL or HSL.
The system of kinetic equations is as follows and we assume that We remark that p and f are downstream quantities, which are calculated from the solutions of the closed equations for s and q. Therefore, a big part of our analysis will focus on the first two equations.
The substrate TG will eventually be fully hydrolysed into products MG and FA through the processes in the reaction scheme in Fig. 2. It is fully hydrolysed even in the case when τ = 0. The cascade of MM enzyme reactions in Fig. 2 with τ = 0 has been previously studied (Storer and Cornish-Bowden 1974).
In either cases, the total concentrations of fatty acid and glycerol moieties are conserved in time, which leads to two conservations laws 3s + 2q + p + f = 3s 0 + 2q 0 and s + q + p = s 0 + q 0 (3) for all times. As a consequence of the conservation laws and s, q → 0 as time t → ∞, system (2) converges at t → ∞ to the following equilibrium state:

Outline
In Sect. 2, we non-dimensionalise model system (2) by introducing a suitable scaling. Moreover, we discuss the general behaviour and illustrate it by plotting a representative phase portrait. The results of this article are summarised as follows: In Sect. 3, we specify three perspectives of judging the role of transacylation and provide a comparative discussion in terms of numerical simulations. As expected from a nonlinear system like (2) the role of transacylation is complicated. It can both speed up and slow down significantly (up to 80% respectively 30%) the production of FAs and MGs. The complexity of the behaviour raises the question, which conditions might allow to use a simplified system as approximation of the full model (2). We will also emphasise that the dynamics depends not only on the parameters, but also on the current state of the hydrolytic machinery. It is possible, for instance, that DG transacylation initially slows down MG production only to altogether speed it up in the long-run of the evolution.
In Sect. 4, we exploit that q changes according to one gain and two loss terms, i.e.q = m 1 − m 2 − 2τ . We discuss in detail the involved time scales and provided Proposition 4.2, which states conditions for a QSSA to be reasonable justified. More precisely, we identify three non-dimensional parameters, which -when assumed sufficiently large -allow q to be approximated by a QSSA.
In Sect. 5, for each of these three parameters, we derive the corresponding first order QSSA models. We illustrate these approximations with numerical simulations including some higher order correction terms.
In the final Sect. 6, we provide a comprehensive discussion of the findings of the paper. We shall point out that DG transacylation can serve as failsafe mechanism, which produces MGs and FAs in the case of a lacking DG hydrolysis. We will show that this benefit comes at the price of slowing down the hydrolytic machinery in cases when DG transacylation and DG hydrolysis compete for substrate. We will also point out an example case, where the numerical plot of q(t) appears to the eye as already well approximated by its QSSA. Yet, the sensitivity derivative w.r.t. the DG transacylation rate constant κ, which is a qualitative measure of the change due to variations in κ, is predicted wrongly by the QSSA reduced system. Hence, as future works, we suggest to include parameter sensitivity derivatives into the discussion whether a QSSA is an admissible model reduction.

Nondimensionalisation and Qualitative Behaviour
We non-dimensionalise the enzymatic model system (2) by rescaling where  Segel and Slemrod (2017), Shoffner and Schnell (1989) to review some well-known and less well-known approaches and a comparison of the pros and cons of each. In the following we choose a nondimensionalisation specifically to our model system and our task of comparing the behaviour with and without transacylation. A natural choice for S ref in a model for TG breakdown is the initial TG concentration s 0 , i.e. S ref = s 0 > 0. We choose also P ref = F ref = s 0 as TG is ultimately processed into MG and FA. Since TG hydrolysis is the first hydrolytic step, which governs all the other activities, we normalise the time scale accordingly, i.e. T ref = s 0 /V 1 . This implies that the s equation rescales likė where the right hand side stated the rescaled non-dimensional s equation. The above choices of reference values lead to the following rescaled version of the q equation: A preferable choice for Q ref is less obvious. A first idea might be to also chose Q ref = s 0 like for the other concentrations. However, since DGs are processed by DG hydrolysis and DG transacylation alongside their production from DG hydrolysis, the actual DG concentration q is typically much smaller than s 0 . In view of this inflowoutflow dynamics of q, another possible choice of an intrinsic reference value Q ref could be a value for which inflow and outflow balance, i.e.q = 0 holds for a certain amount of substrate s, let's say initially when s = s 0 . However, choosing Q ref in this way would not apply to all the parameters we would like to analyse. In particular in the absence of transacylation, a value Q ref such thatq = 0 holds only exists when V 2 > V 1 , see Sect. 4. Even in the cases when such a reference value exists, it depends on σ , which makes a comparison between cases with and without transacylation, i.e. σ = 0 cumbersome. The same problem occurs when trying to scale the method by the principle of minimum simplification or pairwise comparision, see e.g. Shoffner and Schnell (1989).
In this paper, we therefore prefer to set Q ref = K 2 , which is independent of transacylation and allows to normalise the Michaelis constant of the DG hydrolysis rate, thus reducing one parameter. Following this choice, we introduce the dimensionless parameter L = s 0 /K 2 , which appears in all three terms in (5). Moreover, the dimensionless parameter V = V 2 /V 1 denotes the ratio of the two MM maximal velocities. Finally, the nondimensional transacylation term in (5) features the dimensionless parameter σ K 2 2 /V 1 , which compares DG transacylation to TG hydrolysis. We prefer to rewrite this parameter as where we introduce the dimensionless parameter κ = σ K 2 2 /V 2 . This choice has the advantage, that the dimensionless parameter κ directly compares DG transacylation to DG hydrolysis, thus the two downstream pathways for DG. Altogether, the rescaled model (2) becomes The model features three positive dimensionless parameters K , L, V > 0 and the non-negative parameter κ ≥ 0 whether transacylation is considered or not. Given the molecular biological background, one would expect that the initial TG concentration s 0 is at least not smaller than the Michaelis constants K 2 of DG hydrolysis, which suggests L ≥ 1 although we shall not assume this condition. There are also some data (see Kim et al. 2008) suggesting that V 2 should be somewhat larger than V 1 , so that one can maybe expect V ∼ 3, but again we shall not assume this. The scaled initial concentrations become and the non-dimensional form of the conservation laws (3) reads as which holds true for all times t ≥ 0. Moreover, p ∞ = 1 + q 0 /s 0 and r ∞ = 2 + q 0 /s 0 .

Phaseportrait and Global a Priori Bounds
The qualitative behaviour of the s − q system is illustrated by phaseportrait Fig. 3: Given initial data s 0 > 0 and assuming q 0 sufficiently small in comparison, the TG concentration s decays due to the predominant TG hydrolysis. In the same process, the DG concentration q builds up to a maximum value, before decaying alongside the breakdown of TG. We will discuss later the question under which conditions this coupled decay process can be approximated by a QSSA replacing the q equation. In a second regime for large initial values q 0 (compared to s 0 ), transacylation leads to a The large time behaviour of the s − q system is addressed by the following lemma:

Lemma 2.1 (Global bounds and exponential decay of TG/DG concentrations)
Consider solutions to (6) subject to initial data (7). Then, there exist positive constants s max , q max and C 1 , C 2 such that If additional q(0) is smaller than the value whereq(s max ) = 0 holds, then the constant q max can be chosen independently of L.
Proof The existence of constants s max and q max follows directly from the mass conservation laws (8) and the non-negativity of s and q, which imply the bounds s(t) ≤ 1 + 2q(0) 3 L =: s max and q(t) ≤ L + q(0). While s max is uniformly bounded in L ≥ 1 this is not true for this bound on q. However, the bound s max allows to estimatė Denote byq(s) the value of q, whereq = 0 holds (see Fig. 3). Then,q(s max ) constitutes an upper solution to q provided that q(0) ≤q(s max ). Hence, under this additional assumption, we deduce q(t) ≤q(s max ) =: q max , where q max is now independent of L (and also uniform in large V and κ, cf. Lemma 4.1 below).
Next, we estimate Let α and β be arbitrary numbers such that 0 < β < α ≤ 2β. For instance, we can choose α = 1 = 2β. Then, where we set C 2 = min α−β α(K +s max ) , LV 1+q max . It follows from the Gronwall inequality that s 0 and complete the proof.

The Qualitative Role of DG Transacylation
ATGL is the main TG hydrolase and thus believed to be the key regulator for the lipolytic cascade TG → DG → MG. This poses the question if the ability of AGTL to additionally facilitate DG transacylation is of physiological relevance or not?
The question of the qualitative role of DG transacylation to our model systems can be viewed from three different perspectives: (i) How does DG transacylation effect the processing of the substrate TG? (ii) How does DG transacylation effect the production of the product MG? (iii) How does DG transacylation effect the production of the product FA?
In the following, we present numerical simulations of the non-dimensionalised system (6) with MATLAB version 2021b and its built-in ode23s routines. In all plots, we consider initial data consisting only of TG, i.e. s 0 = 1 and q 0 = p 0 = f 0 = 0 and fix the parameters L = K = 1.

Remark 3.1
We point our that normalising the parameters K and L doesn't really change the qualitative behaviour of our model system and in particular doesn't influence the discussion of the qualitative role of DG transacylation. More precisely, K modulates the transition of an approximately constant to an approximately linear Michaelis-Menten rate during the decay of the TGs. Changing K makes no significant qualitative changes to the decay of TGs also since the maximal feedback from DG transacylation is limited by 50%. On the other hand, the parameter L multiplies equally all the processes modelled on the right hand side of the q equation. Hence, changing L doesn't significantly influence the interplay between the processes of DG production and DG degradation.  (6) needed (i) to process 50% of the the initial substrate s 0 , (ii) to produce 50% of the final products p ∞ and (iii) to produce 50% of r ∞ . Bottom row: Relative contribution of DG transacylation 2κq 2 in the processing of DGs at the times as the above plots. The other parameters are L = K = 1 (Color figure online) Figure 4 plots in the top row the simulation times, which are needed-depending on values of the parameters V and κ-(i) to process 50% of the initial substrate TG s 0 , (ii) to produce 50% of the final amount of MG, p ∞ = s 0 = 1, and (iii) to produce 50% of the final amount of FA, f ∞ = 2s 0 = 2. The axes plot the decimal log of the parameter values (V , κ). Note that the times (ii) and (iii) to produce 50% of p ∞ and f ∞ will always be larger than time i) to process 50% of the initial substrate, which is clearly a requirement for both other times to be reached. The time to produce 50% of p ∞ can be much larger than the one of 50% of f ∞ , since TG hydrolysis can reach the FA target on its own, while in order to reach 50% of p ∞ , either DG transacylation or DG hydrolysis is required.
The interpretation of Fig. 4 is far from trivial. For a better understanding we provide the second row in Fig. 4, which plots the relative fraction of the DG transacylation rate to the total rate of DG processing, i.e. 2κq 2 /(2κq 2 + q 1+q ) in the same three cases, at the same times and under the same parameters as the top row. Note that thanks to our choice of non-dimensional parameters in Sect. 2, the ratio 2κq 2 /(2κq 2 + q 1+q ) only involves the parameter κ yet also depends via the solution q on V . All the images in the second row of Fig. 4 show that DG transacylation is dominant at a level of 90% on the left upper corner. On the other hand, DG hydrolysis is dominant at a level of 90% at the bottom and the extended right lower corner. In the in between area, the contour levels decay horizontally for small levels of V and diagonally for larger values of V . The reason for this kink in the contour levels is the change in the values of q in calculating the reaction rates 2κq 2 and q 1+q , respectively. Returning to the interpretation of the top row of Fig. 4, we see in (i) that increasing the value of κ in DG transacylation slows down the processing of the TG substrate. The slowing down effect is the strongest for values of V ∼ 1.
For larger value of V , i.e. in the right upper part, the contour lines of equal times follow the contour lines of equal DG fraction in the picture below. Accordingly, the increased TG processing times corresponds in the right upper part to the increase of the dominance of DG transacylation, which provides increasing feedback into the TG pool.
For small values of V , however, the contour lines in Fig. 4i do not correspond to the ones in the image below, which means that structure of the DG transacylation feedback changes: In this upper left corner, where DG transacylation is overall dominant and roughly constant, the contour lines in Fig. 4i follow the lines where the parameter product V κ is constant, i.e. log 10 κ = C − log 10 V for some constant C. In this part of the image, the increase of the product V κ determines the increase in feedback to the TG pool and the increase in the delay of the 50% TG processing time. Hence, we infer that in this part of the picture, changes in V and κ do not substantially change the values of q. Rather, the increase of the source term V κq 2 in the s equation is mainly due to the change in the parameter product V κ. This is different to the right upper area where the contour lines of Fig. 4i follow the contour lines of the image below! Figure 4ii shows for small values of κ vertical contour levels as the time needed to produce 50% of MG depends under this conditions entirely on DG hydrolysis and the parameter V . Such vertical contour levels for small values of κ are also shown in Fig. 4iii since DG hydrolysis also contributes to the production of FAs, yet with a much more moderate variation of the values compared to Fig. 4ii since TG hydrolysis is a constant source of FA production. In both images Fig. 4ii and iii, the influence of DG transacylation tilts the contour lines away from the middle for increasing values of κ. In the left upper half of Fig. 4ii, the contour lines follow again the lines log 10 κ = C −log 10 V for some constant C, which shows how DG transacylation helps in the case of lacking DG hydrolysis to produce MGs and thus lower the time of 50% p ∞ production. The left upper half of Fig. 4iii shows similar contour lines, since DG transacylation provides feedback into the TG pool, from which TG hydrolysis can release additional FAs. However, this DG feedback effects are somewhat more blurred compared Fig. 4ii due to the ongoing TG hydrolysis, which occurs anyway.
In the right upper half of Fig. 4ii, the contour lines turn to follow the contour lines in the image below and the production of MGs is actually slowed down by increasing κ and thus DG transacylation. At the first glance, this seems counter intuitive since DG transacylation provides an additional pathway to produce MGs. However, we recall that DG hydrolysis yields one MG out of one DG while DG transacylation results in only 1/2 MG. Hence, increased DG transacylation can hinder the twice more efficient production of MGs via DG hydrolysis, which causes the production of MGs to effectively slow down. The left upper half of Fig. 4iii shows again similar contour lines and a similar behaviour as DG transacylation hinders the production of FAs from DG hydrolysis.  (6) compared the same dynamics without DG transacylation, i.e. κ = 0. The left image shows the relative time change for different values of the parameter V and κ at the time when 50% of the the initial substrate s 0 is processed. The middle image shows that the relative time change in order to produce 50% of the final product p ∞ can be either negative or positive, representing cases where DG transacylation compensates for lacking DG hydrolysis or slows down MG production by substrate competition with DG hydrolysis. The right picture shows that also the relative time change in producing 50% of the final FAs f ∞ can be negative or positive, representing cases when the feedback of DG transacylation into the TG pool allows TG hydrolysis to release additional FAs and when DG transacylation shows substrate competition with DG hydrolysis and thus slows down the release of FAs. The other parameters are L = K = 1 (Color figure online) Figure 5 shows relative time changes compared to the dynamics without transacylation, i.e.

Relative Time Changes Compared to Dynamics Without DG Transacylation
where t κ>0 denotes the simulation times like in the cases (i), (ii), and (iii) of Fig. 4 for positive values of κ and t κ = 0 denotes the corresponding simulation time without DG transacylation (yet all other parameters being the same).
On the left, Fig. 5i plots the relative time change with respect to processing 50% of the initial substrate TG s 0 . For small values of κ, the image show small, positive relative time changes, which corresponds to the small, positive feedback of DG transacylation into the pool of TG delaying the time to have processed 50%. With increasing κ, the delay effect becomes more and more significant reaching up to 70% for values of V 1. The contour lines of Fig. 5i correspond to Fig. 4i, which is an easy consequence of the fact that for κ = 0, the processing time of 50% TG is due to TG hydrolysis alone and thus the same for all values of V .
The middle image of Fig. 5ii shows the complexity of the interplay DG transacylation vs. DG hydrolysis. For small values of V , increasing κ leads to DG transacylation becoming dominant and therefore speeding up MG production in Fig. 5ii up to 80% and higher. Here, DG transacylation acts as a by-pass process in order to produce MGs despite lacking DG hydrolysis. A similar speed is shown in Fig. 5iii, where DG transacylation feeds back into the TG pool, from which TG hydrolysis can release additional FAs. In contrast, for larger values of V 1, Fig. 5ii shows that increasing κ and DG transacylation results in a substrate competition to DG hydrolysis, which is nevertheless the more efficient process since it produces 1 MG out of 1 DG, while DG transacylation yields on 1/2 MG. Therefore, the 50% MG production times are slowed down by 40% and more. This is again similar in Fig. 5iii, where DG transacylation hinders the immediate production of a FA by DG hydrolysis (even if TG hydrolysis can later release 1/2 FA form the 1/2 TG produced by DG transacylation).

Effects of DG Transacylation Depend on the State of the Hydrolytic Machinery
We point out the effects of DG transacylation also depend on the current stage within the time evolution of the hydrolytic machinery. Figure 6 plots the relative time changes comparing system (6) with L = K = 1 and κ = 10 to the same system with κ = 0 as function of V and for different percentages x% of TG process and the corresponding percentages (100 − x)% of MG and FA production. Figure 6i plots the slowdown in the TG substrate processing times. For V 1, i.e. under conditions where DG hydrolysis is lacking, the slowdown of TG process is largest for small values of x, i.e. at the beginning of hydrolysis, when DG concentrations are largest and DG transacylation provides therefore the largest feedback into the TG pool. For V 1, i.e. under conditions where DG hydrolysis is present, the effects of DG transacylation and thus the slowdown get smaller and almost independent of x for large values of V . Figure 6ii shows for V 1 a strong 80% acceleration of MG production thanks to the MGs produced by DG transacylation, which serves under this conditions as a bypass of the lacking DG hydrolysis. In contrast, for values V 1, Fig. 6ii displays that MG production is slowed down for smallish x = 10, 25, . . ., i.e. at the beginning of hydrolysis. Yet, later in the evolution, when x = 75, 90, the early slow down is not only compensated but MG production is altogether accelerated due to DG transacylation. This example underlines the nonlinear behaviour of the effects of DG transacylation and how they depend on the current state of the hydrolytic machinery. This is important in view of biochemical experiments, in particular when comparing in vitro cell assays experiments to in vivo experiments, where effective DG concentrations are very difficult to determine. Figure 6iii shows for values V 1 and smallish x = 10, 25 a very strong speed up of 80% in the FA production time since DG transacylation creates additional TGs, which serve as a substrate for TG hydrolysis and therefore produce additional FAs. However, later in the evolution for x > 50%, there is no noticeable speed up anymore since only 50% of the final FAs come from TG hydrolysis (and in-between DG transacylation). The remaining x-50% of the final FAs can only be produced via DG hydrolysis, which is just as slow as it is.
On the other hand, Fig. 6iii shows for values V 1 a slow down effect on FA production from DG transacylation, as it competes with DG hydrolysis for substrate and (after a TG hydrolysis step) yields only half of the FA release compared to DG hydrolysis.

Model Reduction and QSSA Prerequisites
The most important equations in system (6) is the closed subsystems for TG and DG concentrations: Here, we have scaled (w.l.o.g.) the initial TG concentration to one and introduce the function d(q) as a short hand for the sum of the DG hydrolysis and DG transacylation reaction rates in Eq. (11).
In this section, we investigate under which conditions system (10)-(11) can be reduced to a simpler model system by means of a quasi-steady-state-approximation (QSSA). QSSAs exploit that some quantities change slowly in time compared to typical time scales of the overall dynamics. Then, a zero order QSSA sets the time derivative of these quantities to zero, thus treating it quasi like for a steady state. From the viewpoint of asymptotic analysis, QSSAs are singularly perturbed ODE problems, see e.g. Holmes (1995), Kuehn (2015) and the many references therein. Hence, approximating slowly changing quantities leads to initial layers (except in the special case where the initial data of the slowly changing quantities already matches their QSSA values), which can be resolved, for instance, by matched asymptotics. More accurate approximations can be obtained by calculating higher order correction terms. Note that (10) can't qualify for a QSSA since the positive feedback provided by DG transacylation is at most half of TG hydrolysis. Even if DG hydrolysis could be neglected, then the feedback DG transacylation and subsequent TG hydrolysis provides is only like TG → DG →1/2 TG → 1/2 DG …, see Fig. 2.
However, Eq. (11) features a gain and two loss terms (TG hydrolysis vs. DG hydrolysis and DG transacylation), where the loss terms are monotone increasing in q. Hence, the dynamics of (11) tends to balance gain and loss terms. If this happens sufficiently fast, a QSSA is expected to serve as a simplified approximative model replacing (11) by an algebraic equation (and possible higher order correction terms).
As a first example, we divide Eq. (11) by L and assume that L is a large parameter. Then, the small parameter ε = L −1 suggests a QSSA εq 0 whenever L is sufficiently large. In the following, we will see that also the parameters V and κ yield QSSAs provided they are sufficiently large. In these two cases, one rescales q → εq in (11) with either ε = V −1 or ε = κ −1/2 , which lead again to a left hand side εq 0.
In all three cases, system (10)-(11) can be approximated by a zero order QSSA model system by setting εq = 0 and replacing the differential equation (11) by the algebraic equation where the functionq =q(s) is the unique positive solution to (12) (recall the monotonicity inq) as long as it exists. In particular,q is not defined when κ = 0 and m 1 (s) V ≥ 1 which contradicts the bound m 2 (q) < 1 for allq ∈ R + . Note, that the explicit solution to the third order polynomial (12) by Cardano's formula is of little practical use and that we will make use ofq either implicitly or approximatively.
The following Lemma 4.1 states a parameter condition which implyq ≤ 1/2 and provides a useful approximation formula ofq in this case. Note that in practical examples, we expect QSSA levelsq to be rather small, in particular if V and/or κ are large enough.
has a unique positive solutionq for every input value I > 0. The explicit solution by Cardano's formula is too cumbersome for practical use. Instead, we derive an approximation formula for sufficiently smallq, i.e. inputs I . In particular, since d(q) is monotone increasing, a sufficient condition forq ≤ 1 2 is obtained from estimating m 1 V ≤ 1 V for all s ≥ 0 and K > 0: Next, for sufficiently smallq (to be determined), we can rewrite d(q) = 2κq 2 +q −q 2 +q 3 1 +q = I The last term on the left hand side δ(q) :=q 3 1+q will be of lower order for sufficiently smallq. A formal asymptotic expansionq =q 0 + δq 1 yields In order to make sure that the zero order approximationq 0 is real in the cases 0 ≤ 2κ < 1, we enforce 1 > 4I (1 − 2κ) by estimating I ≤ 1 V independently of s and by imposing the sufficient constraint Next, we calculate thatq 0 ≤ 1 2 implies δ(q 0 ) =q 3 0 1+q 0 ≤ 1 12 . Thus, (15) yields δ(q) = δ(q 0 ) + O(δ) and thus δ(q) is of order O(10 −1 ). On the other hand, straightforward calculation shows thatq 0 ≤ 1 2 is equivalent to I ≤ (1 + 2κ)/4. Since m 1 (s) ≤ 1, a sufficient condition which ensuresq 0 ≤ 1 2 independently of s and K is given by assumption (13) above. Finally, it is straightforward to check that (13) always implies constraint (16), which completes the proof.

Conditions for QSSA
A QSSA which approximates εq 0 is expected to feature an initial layer whenever the initial data q(0) do not equal the initial values of the approximationq(s(0)), see e.g. Holmes (1995). Indeed, we expectq ∼ O(ε −1 ) during the initial layer except for initial data q(0), which are already O(ε)-close toq(s(0)).
Here, it is natural to consider initial data q(0) = 0, i.e. that there is initially no DG, which implies that there is an initial layer behaviour. (Other initial data which are not O(ε)-close toq(s(0)) can be treated the same.) For a QSSA to be a useful and valid approximation, the initial layer must be short compared to the dynamics which is of interest.
In the following, we derive necessary conditions on the parameters as well as bounds on the shortness of the initial layer. The proximity of q(t) toq(s(t)) is -in our opinion -best studied after rewriting (11) as perturbation π around the QSSA variableq(s(t)), i.e.
Note that the squared bracket on the r.h.s. is always strictly positive (since q = π +q ≥ 0). Hence, the equality follows, which is a useful observation. Figure 7 compares the evolution of q(t) from (10) andq(s(t)) for different values of V : Starting at the initial data without DG, i.e. q(0) = 0 resp. π(0) = −q(s(0)), TG hydrolysis yields s(t) to continuously decreases from its starting value s(0) = 1 (not plotted) and q(t) to grow towards a maximal value q(t m ) at a time t m > 0 wheṅ q(t m ) = 0 holds. Thus (18) implies that at t m > 0 holds: q(t m ) =q(s(t m )) =:q m and π(t m ) = 0.
At time t m of maximal DG concentration, the rates of DG hydrolysis and DG transacylation balance exactly TG hydrolysis. For times larger than t m , the DG concentration q(t) decays. Hence, (18) implies π(t) > 0 for t > t m , which can be interpreted that the decay of the rates of DG hydrolysis and transacylation lags behind the decay of the TG hydrolysis rate caused by the decay of s(t), see also phaseportrait Fig. 3.
Altogether, Fig. 7 suggests that a QSSA can be used provided that the initial layer interval [0, t m ] is sufficiently short and that after t m the perturbation π(t) is sufficiently small. Figure 7 illustrates that we can expect a QSSA for q to hold for sufficiently large values of V . Analog plots can be made in favour if a QSSA for sufficiently large values of L and κ.
We would like to estimate the duration of the initial layer in the following. First, we derive the evolution equation of π(t). By taking the derivate of (12), we calculate withq short forq(s)q where we use the formulas: By recallingq =q +π , we obtain the evolution equation for π : subject to the initial data π(0) = −q(s(0)) = −q(1) < 0. The time scales of Eq. (21) are key to understanding necessary conditions for a QSSA. The first term in (21) acts as a stabilising terms and aims to relax π(t) towards zero (since the square bracket is positive). On the other hand, the second term is always positive (asṡ < 0). Hence, it helps q to reach the maximal levelq(s(t m )) faster, but afterwards becomes a destabilising factor causing π > 0 as long asq < 0 resp.ṡ < 0. Roughly speaking, the second term quantifies the adaptation of the QSSAq due to the decay of s.
The eigenvalue −LV d (q) < 0 characterises the linear stability of the QSSAq(s) provided that the decay of s is negligible.

Remark 4.1
The fast timescale of classical QSSA for a single substrate -single enzyme reaction is often calculated from the eigenvalue LV d (q) (Murray 2002, p. 179), i.e. T = 1/ (LV d (q)). The same time scale is obtained when using Rice's scaling by inverse rates to our system, see e.g. Shoffner and Schnell (1989). In the following, we will derive a better estimation of the timescale, see (24) below, which is more specific to our system.

Time Scale Analysis
First, we consider (21) subject to the initial data π(0) = −q 0 = −q(s(0)), q(0) = 0 andṡ(0) = −m 1 (s(0)) and obtaiṅ By recalling (12) where the last equality holds due toq(t m ) = 0, i.e. m 1 = V d(q) and q(t m ) =q m . Thirdly, we consider the halfway-up timepoint when π = −q (s) 2 holds in [0, t m ]. Evaluatingπ yieldṡ For the first square bracket, (12) implies the elementary bounds while for the second square bracket, we obtain where the lower bound follows fromq > 0 implying 2V κq 2 < m 1 − V m 2 . Note also that 1 2+q 1 2 for small values ofq, which we expect for large values of V and κ. Hence, we observe thatπ We therefore take this average as a typical value ofπ throughout the initial layer and introduce as a typical time scale for the initial layer the change in q (that is q(t m ) − q(0) =q(s(t m ))) divided by the above typical value ofπ : In the following, we derive for T the more convenient upper bound T fast in (29) depending only onq m =q(t m ) = q(t m ) and the parameter K . First, sinceqd (q) = 4κq 2 +q (1+q) 2 , we obtain from (12) the bounds and moreover, θ ∈ [ 1 2 , 2] forq ≤ 1. Next, we rewrite T as and identify the terms T 1 and T 2 in (22) and T 3 in (23) using the notationsq 0 := q(s(0)) =q(1) andq m =q(s(t m )): where we have used for T 3 in (23) thatq m =q(t m ) = q(t m ) and m 1 − V κq 2 = ψm 1 for some ψ ∈ [ 1 2 , 1].
In particular, we will estimate (27) and (28) as sinceq m ≤q 0 , θ ≤ 2 and ψ −1 ≤ 2. We then obtain the upper bound The time scales for the initial layer T and its more practical upper bound T fast have to be compared a time scale representing the decay of the substrate s. Since m 1 2 ≤ −ṡ ≤ m 1 holds during the initial layer (due toq > 0 implying 2V κq 2 < m 1 − V m 2 ), we choose the Michaelis Menten modelṡ = − 3 4 m 1 (s) to provide a reference time of the decay of s(t). Straightforward integration with s(0) = 1 yields As comparison for the times scale of the initial layer, we introduce the time T 90% when s(T 90% ) = 9 10 , which means that at time T 90% about 90% of the TG substrate is still to be processed. We calculate A QSSA relies on the separation of times scales. Here we impose that during the fast time scale of the initial layer, at most 10% of the substrate are processed: We state the following Proposition:  Fig. 7, condition (32) for L = 1 = K requiresq m 0.068. This is clearly satisfied for V = 10 whereq m ∼ 0.028, but not for V = 2, whereq m ∼ 0.073. In the latter case, we calculate T 90% ∼ 0.26 while (29) yields that T fast 0.29 and we see that T fast ≤ T 90% is not satisfied. Nevertheless, when looking at the case V = 2 in Fig. 7, one would think that the QSSAq constitutes already a very nice approximation of q after the initial layer. One could suggest to relax the condition T fast ≤ T 90% , to for instance T fast ≤ T 84% and calculate analog to above the conditionq 0.11, which is satisfied for V = 2. However, we will elaborate this example further in the discussion and show that one has to be careful with being generous in accepting a QSSA.

Remark 4.3
We point out that T 1 = 2(K + 1)q m L is crucial for deriving condition (31). Without T 1 , we obtain the contradiction This is not surprising since T 2 and T 3 are contributions, which are cause by the adaptation ofq due to the decay of s(t) and it would be indeed surprising if this feedback alone would enable a sufficiently short initial layer.

Remark 4.4
From the three large parameters, which allow for QSSAs, the case L large is analog to the classical argument which justifies the Michaelis-Menten QSSA from a large ratio of substrate to enzyme, see e.g. Murray (2002), Shoffner and Schnell (1989).
The cases V and κ large are related, but also different in the sense that a large downstream rate (either V 2 V 1 or σ K 2 2 V 1 , see (6)) implies that the maximal values of q are small and thus quickly reached after a short initial layer, altogether allowing to replace q by its QSSAq.

Asymptotic Expansions
In this section, we perform asymptotic expansions of the non-dimensionalised full model system (6), which we recall for the convenience of the reader We will consider the three cases which allow a QSSA: sufficiently large L, V and κ. In the case of large L, the left hand side of q equation in (35) features the small parameter ε = 1 L . By setting ε = 0, the q equation reduces to the algebraic equation (12). The corresponding zero order QSSA model system reads as We point out that the QSSAq(s(t)) from (36) is not the same asq(s(t)) with s(t) solving the full system (35), which we have studied in the previous Sect. 4.1. However, the feedback fromq(s(t)) to the evolution ofs in (36) is rather weak (recall m 1 2 ≤ −ṡ ≤ m 1 ) and was therefore neglected for the sake of clarity in the derivation the conditions of Proposition 4.2, which are anyway not sharp.
In the following, we improve the the approximative quality of the zero order QSSA system (36) by deriving higher order corrections. We will first consider the case of large L and later the cases of large V and κ. We will see that all these cases lead to somewhat different approximative systems and explain the difference in terms of the interpretation of the parameters. In the case of large V , we will even provided second order corrections, since it has been suggested that V ∼ 3 is a realistic value, which would make second order terms necessary.

First Order Asymptotic Expansion for L Large
The parameter L = s 0 K 2 represents to ratio of initial TG substrate to the MM constant of DG hydrolysis. One could think that s 0 should be significantly larger than a realistic value K 2 , which would be in favour of considering L large. However, this is not so clear since ATGL and HSL are assumed bound to the surface of lipid droplets, the TG substrate which is actually available might need to be also near the surface. So effectively, the available s 0 might not be larger than K 2 .
With ε = 1 L being a small parameter, we rewrite the q equation of (35) as Inserting the asymptotic expansion q = q 0 + εq 1 + O(ε 2 ) into (37) and comparing coefficients yields in order ε 0 yields m 1 (s) = V d(q 0 ), which implies and the zero order approximation is given by the solution of the algebraic QSSA equation (12). In order ε 1 , we obtainq 0 = −V d (q 0 )q 1 , which yields where we have calculatedq from differentiating m 1 (s) = V d(q). In order to obtain an expression of q 1 in terms of s only, we insert the s equation of (35) into the right hand side of (39), we have withq =q(s(t)) Hence, we have as first order QSSA We are now able to derive the first order approximative system. Inserting Eqs. (38) and (39) into the s equation of full system (35) yieldṡ and, thereforeṡ Next, we insert Eqs. (38) and (39) into the p equation of (35): where q 1 is given by (40). Finally, inserting Eqs. (38) and (39) into the f equation of (35) yieldsḟ where q 1 is given by (40). Altogether, we obtain as first order QSSA model system Figure 8 plots comparisons of q solving the full system (35) with the QSSAq and the first order approximationq + εq 1 for increasing values of L = 1, 2, 4, 10 while the other parameters are normalised. The case L = 4 illustrates particular nicely the improvement due to the first order approximation.

Remark 5.1
In the first order reduced model system (45), the s equation is a scalar equation. Its solution then allows to integrate up p and f . However, system (45) depends heavily onq(s(t)), which is the awkward solution of the third order polynomial (12), which make the practical use of system (42)-(44) questionable.

First and Higher Order Asymptotic Expansion for V Large
We consider the case V large and V 2 V 1 holds. In this parametric regime, DG hydrolysis, that is the loss term m 2 in the q equation in system (35), is much faster than the inflow term m 1 from TG hydrolysis. Hence, q is expected to be O(ε) and ought to be rescaled accordingly, i.e. q = ε Q with ε = 1/V . Hence, we obtain from non-dimensionalised system (35) the following rescaled model Note that with q = O(ε), DG hydrolysis is in zero order approximation a linear term (since there is little substrate) and that DG transacylation is a term of order O(ε). As a consequence, the asymptotic expansion in the case V large is simpler and more explicit than in the previous case of L large. We insert the ansatz Q = Q 0 + ε Q 1 + O(ε 2 ) into the Q equation in (46) and compare the coefficients to obtain Q 0 (t) = m 1 (s(t)), Inserting (47) into (46) yields the first order QSSA system Note, that solutions to (48) conserves total glycerol and total FAs only up to the zero order terms. In first order, the conservation laws do not hold due to the term εm 1 m 1 (s) in the equations forṗ andḟ , which stems form the first term in Q 1 and is hence proportional toṡ and a results of the asymptotic expansion. The s equation in (48) is scalar like in system (45) and has the added benefit of not requiring to solve forq.
We again illustrate the asymptotic analysis by numerical simulations, which plot the asymptotic expansion of the QSSA quantity q = ε Q and ε = 1/V . Unfortunately, as of now we are lacking experimental data, which would allow to identify V from experiments. However, an educated guess suggests that values V ∼ 3 could be expected. Hence, we provide the asymptotic expansion of q up to third order terms. Straightforward (and somewhat lengthy) calculations show that the following asymptotic expansion of q: Remark 5.2 Note that in the asymptotics of large V , the zero order approximation q 0 = 0, which follows from expanding the rescaled system (46), but could have already been see from the q equation in (35): There, for large V follows already that q = O 1 V holds. Figure 9 plots q(t) as part of the solution of (35) and shows in comparison the approximations by q 1 , q 1 + q 2 and q 1 + q 2 + q 3 given by (49) for parameters L = K = κ = 1 and where s(t) is part of the solution of (35). We see that already for V = 3, the solution q and the third order approximation q 1 + q 2 + q 3 are in very good agreement after the initial layer, while this is not really true for the first and second order approximation. Even better approximation is obtained for larger values of L.

First Order Asymptotic Expansion of q: Ä Large
We consider the case κ large. This is a parametric regime where σ K 2 2 V 2 holds and DG transacylation, which is a loss term in the q equation in (35) is a fast process. Like in the case of V large, q is again expected to be small. Given the quadratic nonlinearity of transacylation, it turns out that the correct rescaling is q = ε Q with ε = 1/ √ κ. By accordingly rescaling the non-dimensionalised system (35), we obtain Note that the this rescaling DG transacylation is dominant and or order O(1) while DG hydrolysis is a O(ε) term and linear up to order O(ε 2 ) (since the substrate concentration q is small). Hence, the asymptotic expansion in this case is again different to the cases V or L large. We apply the asymptotic expansion Q = Q 0 + ε Q 1 + O(ε 2 ). Insertion into (50) and comparison of the coefficients yields where we expressedṡ by means of the s equation in (50). Inserting (51) into (46) yields the first order approximative system We illustrate the asymptotic analysis by numerical simulations, which plot the asymptotic expansion of the QSSA quantity q = ε Q and ε = 1/ √ κ. Rescaling (51) back yields the following asymptotic expansion of q up to second order terms, where (like in the case V large) the zero order approximation q 0 = 0: Figure 10 plots q(t) as part of the solution of (35) and shows in comparison the approximations by q 1 and q 1 + q 2 given by (53) for parameters L = K = V = 1 and where s(t) is part of the solution of (35). We see very good agreement only vor very large values like κ = 100, which reflects that the scaling factor is √ κ. Already with κ = 25, the second order QSSA q 1 + q 2 is not longer a very good approximation.

Conclusion and Outlook
Intracellular lipolysis is a central metabolic pathway involved in a variety of cellular processes including energy production, signal transduction, and lipid remodelling. The main lipolytic enzymes (i.e. ATGL, HSL, and MGL) have been functionally characterised in great detail. Over the last decades, gain-and loss-of-function studies in different model organisms corroborated their crucial roles in neutral lipid catabolism (Grabner et al. 2021;Schreiber et al. 2019). Nevertheless, lipolysis is far from being completely understood. We are convinced that the classical view of lipolysis as a linear process involving the consecutive action of ATGL, HSL, and MGL is too simplistic. On the one hand, published data suggest that lipolysis is an enzymatically redundant process (ATGL not only hydrolyses TGs but also DGs (Zimmermann et al. 2004), vice versa HSL not only hydrolyses DGs but also TGs and MGs (Fredrikson et al. 1981)). On the other hand, three independent studies convincingly showed that ATGL also exhibits DG transacylation activity at least in vitro (Jenkins et al. 2004;Zhang et al. 2019;Kulminskaya et al. 2021), indicating that intracellular lipolysis is a non-linear process. The in vivo relevance of DG transacylation is still elusive though. Kinetic parameters of model system (6) are currently not available. In particular, the parameters κ and V for DG transacylation and DG hydrolysis are unknown since these two functions can not be distinguished by any current experimental protocol. In a work in progress, we are developing mathematical approaches to identify kinetic parameters from experimental data using minimisation of tracking functional as well as Bayesian estimators. First preliminary results of parameter identification (yet for more specific mathematical models than system (6)) suggest values of κ and V , which put the system dynamics somewhere near the middle of the parametric plots Figs. 4 and 5. A plausible guess suggests V somewhat larger that one and DG transacylation on average maybe 20% of DG hydrolysis, but with a strong dependence on the substrate and an increase up to 50% in some experiments.
For these parameter values, the effects of DG transacylation as shown in Fig. 5 are rather neutral: A slight delay in substrate processing time but negligible effects on the time of MG and FA production compared to the system without DG transacylation. This might be a reason why some researchers believe that DG transacylation does not play a significant role in vivo: Not because it is not active, but because it behaves rather neutrally.
However, if the hydrolytic machinery is shifted away from these conditions, DG transacylation will play a more significant role. First, if DG hydrolysis is lacking, DG transacylation serves as a failsafe mechanism by producing MGs and even indirectly FAs when TGs formed are hydrolysed further to DGs. However, this failsafe mechanism comes at a price: If DG hydrolysis is present, a highly effective DG transacylation will compete with DG hydrolysis. This will cause downstream delays by forcing the lipolysis cascade TG → DG → MG rather into a semi-open loop of TG → DG → 1/2 TG → 1/2 DG → 1/4 TG etc, see e.g. Fig 5 and its illustration of the downstream delays.
It seems highly remarkable that our current preliminary guess on parameter values places the kinetics of the lipolytic machinery in the middle of these two options: failsafe mechanism versus DG substrate competition. One can speculate that the lipolytic machinery evolved to operate at a flexible state rather than at the maybe more efficient state without DG transacylation. We interpret our observations that there is biological purpose behind both DG hydrolysis and DG transacylation constituting the lipolysis machinery.
As a second biological purpose, we can speculate that various effects of DG transacylation have the ability to homogenise the process of lipid hydrolysis over the surface of lipid droplets. We recall that the ratio of DG transacylation to DG hydrolysis depends on the DG concentration levels and that DG transacylation becomes dominate for sufficiently large concentrations of DGs. The enzyme localisation on the surface of LDs can be temporally and spatially heterogeneous (Egan et al. 1992;Granneman et al. 2009;Sztalryd et al. 2003). Local absence of HSL would lead to local DG accumulation as DGs can't leave the LD. It can by speculated that DG transacylation may serve as a regulator of spatial heterogeneities of DG concentrations on the surface of the lipid droplet by processing local surpluses (and without having to wait for diffusion to spread DGs). Hence, DG transacylation would contribute to an overall more homogeneous and controlled downstream processing, which is certainly an important physiological factor. However, ATGL-mediated DG transacylation averaged over the entire surface of the lipid droplet cannot fully compensate for complete HSL deficiency in vivo, since HSL knockout mice accumulate DGs in adipose and non-adipose tissues (Haemmerle et al. 2002).
There are many potential directions of future research. From a mathematical perspective, a rigorous proof of the convergence of the fast-slow ODE systems to the limiting QSSA is work in progress. The proof requires a sufficiently sharp perturbation argument for nonlinear, non-autonomous ODE systems. Unfortunately, estimations of the fundamental matrix of non-autonomous systems like in Teschl (2012, Chapter 3) are currently insufficient, which we strongly consider to be a technical artefact rather than the results of the system behaviour. As a by-pass, we are currently trying to find a suitable change of variables, which might allow to improve the estimates on the non-autonomous semigroup.
From a modelling and interdisciplinary perspective, a central goal is the identification of the lipolytic kinetic parameters from experimental data. In particular, we aim to understand the dynamics of TG breakdown in adipocytes via cell-based assays and mice experiments. However, we will have to work with a limited number of data points both in adipocyte cell assays and even less in mice experiments. We expect that reduced QSSA models like presented in Sect. 5 will be important in order to fit parameter from in vivo data.
There is, however, an additional mathematical question concerning the use of QSSAs. As part of the process of mathematical modelling, a sensitivity analysis of the parameter is an important step. In particular, we are studying parameter sensitivity derivatives and how they evolve along solutions of the system. Parameter sensitivity derivatives provided very useful insights into the qualitative behaviour of the system and also into how challenging the inverse problem of identifying parameters from measurements will be.
In the spirit of this paper, which is to discuss the role of DG transacylation, the following system considers the sensitivity derivatives with respect to the parameter κ: subject to homogeneous initial data ds dκ t=0 = dq dκ t=0 = dp dκ t=0 = d f dκ t=0 = 0.
Note that, for instance, if ds dκ is positive, then the dynamics observed for increased values of κ will show also increased values of s.
Hence, differentiation shows dp dκ = − ds dκ , which implies that dp/dκ is negative for all times, initially decreases until eventually approaches 0, see Fig. 11. Figure 11 plots the parameter sensitivity derivatives ds/dκ and dp/dκ of the QSSA model for the parameters L = K = 1 and κ = 16, which illustrates this analysis for different values of V = 0.1, 1, 2, 10.
In comparison, Fig. 12 plots the evolution of the κ-sensitivity derivatives for the full system (54) for parameters L = K = 1, κ = 16 and the same values of V = 0.1, 1, 2, 10. We observe that the qualitative behaviour of ds dκ and ds dκ corresponds well. However, in the cases V = 0.1, 1, 2, the product sensitivity derivative dp dκ of the full system (54) is positive for some initial time interval. This is in contrast to our analysis of the QSSA system, where dp dκ = − ds dκ ≤ 0 holds. For the full system (54), the differentiated conservation law (8), i.e. ds dκ + dq dκ + dp dκ = 0 allows for dp dκ > 0 as long as dq dκ is sufficiently negative, see Fig. 12. Fig. 11 The evolution of the κ-sensitivity derivative ds/dκ and dp/dκ of the QSSA system for different values of V and L = K = 1, κ = 16 (Color figure online) Fig. 12 The evolution of the κ-sensitivity derivative of ds/dκ, dq/dκ, dq/dκ and dr/dκ of the full system (54) for different values of V and L = K = 1, κ = 16. Note that in contrast to the QSSA system plotted in Fig. 11, the derivative dq/dκ turns positive for short times (Color figure online) Hence, there are parameters for which the full system predicts an increase of MG production due to increased DG transacylation when the QSSA system predicts a decrease of MG production due to increased DG transacylation. In particular, this contradiction is observed for parameter value V ≤ 2, where for V = 2 Fig. 7 plots (for the same parameters L = K = 1 and κ = 16) a seemingly very good approximation of q byq.
This example of disagreement between full and QSSA model suggests that an analysis of validity of a QSSA should also include the sensitivity derivatives. Note that already for the classical Michaelis-Menten QSSA analysis, it is known that there are parameter regions where a QSSA can't be validated, see Eilertsen et al. (2022). In our current paper, Proposition 4.2 identifies parameter conditions, which are necessary for a QSSA, but we do not provide an analysis which shows that these conditions are sufficient. Such a theorem is significantly beyond the scope of the this paper and left for future works.