Fine-tuning of AMPK–ULK1–mTORC1 regulatory triangle is crucial for autophagy oscillation

Autophagy is an intracellular digestive process, which has a crucial role in maintaining cellular homeostasis by self-eating the unnecessary and/or damaged components of the cell at various stress events. ULK1, one of the key elements of autophagy activator complex, together with the two sensors of nutrient and energy conditions, called mTORC1 and AMPK kinases, guarantee the precise function of cell response mechanism. We claim that the feedback loops of AMPK–mTORC1–ULK1 regulatory triangle determine an accurate dynamical characteristic of autophagic process upon cellular stress. By using both molecular and theoretical biological techniques, here we reveal that a delayed negative feedback loop between active AMPK and ULK1 is essential to manage a proper cellular answer after prolonged starvation or rapamycin addition. AMPK kinase quickly gets induced followed by AMPK-P-dependent ULK1 activation, whereas active ULK1 has a rapid negative effect on AMPK-P resulting in a delayed inhibition of ULK1. The AMPK-P → ULK1 ˧ AMPK-P negative feedback loop results in a periodic repeat of their activation and inactivation and an oscillatory activation of autophagy, as well. We demonstrate that the periodic induction of self-cannibalism is necessary for the proper dynamical behaviour of the control network when mTORC1 is inhibited with respect to various stress events. By computational simulations we also suggest various scenario to introduce “delay” on AMPK-P-dependent ULK1 activation (i.e. extra regulatory element in the wiring diagram or multi-phosphorylation of ULK1).


I. Describing the theoretical analysis
In this section, we briefly describe the mathematical approach used to study the autophagy induction. A system level view can be developed by bringing together the components and interactions reported in the literature. Such a network can be translated into a set of mathematical equations that describe how each component concentration/activity in the network changes with the time. The rate of change of a component is described by ordinary differential equation (ODE) based on biochemical reaction kinetics (see equation below). Each biochemical reaction is represented as a term on the right hand side of the ODE for a component participating in the reaction 1,2 . Each reaction in the network can be described either by using law of mass action or Michaelis-Menten kinetics [3][4][5] .
A generic differential equation describing the temporal changes of protein Xa is composed of two parts: production and consumption terms. The production can be given by protein synthesis and/or an activation term, while the consumption can be given by protein degradation and/or inactivation term. Usually synthesis, degradation, binding and dissociation reactions are described by mass action kinetics, whereas protein activity can be described either by mass action or Michaelis-Menten kinetics 3,6 . For example, if the activity of protein is controlled by covalent modification involving multi-site phosphorylations, Michaelis-Menten kinetics provides a good approximation for the process 7,8 . The value of parameters (rate constants, Michaelis constants) and initial conditions have to be specified in order to solve ODEs. The non-linear nature of biological processes makes it difficult to find the solution of ODEs analytically hence the equations has to be solved numerically. The equations can be solved using different numerical integration methods that are implemented as solvers in many of the freely available computer software.
Solving a set of non-linear ODEs gives the time evolution of the protein concentration/activity called time courses. Further, ODEs can be solved to obtain the input-output relationship called as signal response curves or as one parameter bifurcation diagram 2,3,9 . An input is the signal strength that is varied to obtain the steady state behaviour of the control system. This helps to capture the qualitative changes in the behaviour of the system. For example, the system behaviour can become abrupt and discontinuous when signal strength is increased from a low value to high value. A point at which such a qualitative change in the system occurs is defined as bifurcation point 2 .
In this work, the temporal profiles and signal response curves were computed numerical using XPP-AUT. All the simulations presented in the text are based on the following XPP codes. The rate constants (k) have the dimension of min -1 and Michaelis constants (J) are dimensionless. The proteins levels/activities are given in arbitrary units (a.u). The starting parameter set was able to refer to physiological conditions. The parameters values were perturbed to capture all the possible qualitative behaviours that the given network can exhibit.
The phosphorylation sites were also searched with the help of PhosphositePlus (https://www.phosphosite.org/homeAction.action). From this database, the consensus phosphorylation motif of AMPK was used, which was created with the help of several wellknown phosphorylated sequences of AMPK substrates. The same amino acids within the identified potential AMPK phosphorylation sequences were marked with red colour (see the Peptid column of Table 1

Introducing the theoretical analysis of phosphorylation site search on ULK1
The potential Ser and Thr phosphorylation sites of AMPK were identified on ULK1 sequence by Group-based Prediction System 5.0 (http://gps.biocuckoo.cn/). The threshold was high under prediction. The Supplementary Table 2. includes the position of Ser and Thr amino acids, the catalytic subunits of AMPK, the peptid sequence in the near of Ser and Thr amino acids and the score values. The score value was calculated by GPS algorithm to evaluate the potential of phosphorylation. The higher the value, the more potential the residue is phosphorylated (http://gps.biocuckoo.cn/).
The phosphorylation sites were verified by NetPhos 3.1 (http://www.cbs.dtu.dk/services/NetPhos). Those phosphorylation sites were collected and checked where the phosphorylation kinase was unknown (see the NetPhos column in the Table 2.). The score above 0.500 indicates positive predictions.
The phosphorylation sites were also searched with the help of PhosphositePlus (https://www.phosphosite.org/homeAction.action). From this database, the consensus phosphorylation motif of AMPK is used, which is created with the help of several well-known phosphorylated sequences of AMPK substrates. The same amino acids within the identified potential AMPK phosphorylation sequences were marked with red colour (see the Peptid column of Supplementary