## Abstract

Recent studies on *Caenorhabditis elegans* reveal that gene manipulations can extend its lifespan several fold. However, how the genes work together to determine longevity is still an open question. Here we construct a gene regulatory network for worm ageing and quantify its underlying potential and flux landscape. We found ageing and rejuvenation states can emerge as basins of attraction at certain gene expression levels. The system state can switch from one attractor to another driven by the intrinsic or external perturbations through genetics or the environment. Furthermore, we simulated gene silencing experiments and found that the silencing of longevity-promoting or lifespan-limiting genes leads to ageing or rejuvenation domination, respectively. This indicates that the difference in depths between ageing and the rejuvenation attractor is highly correlated with worm longevity. We further uncovered some key genes and regulations which have a strong influence on landscape basin stability. A dynamic landscape model is proposed to describe the whole process of ageing: the ageing attractor dominates when senescence progresses. We also uncovered the oscillation dynamics, and a similar behaviour was observed in the long-lived creature *Turritopsis dohrnii*. Our landscape theory provides a global and physical approach to explore the underlying mechanisms of ageing.

## 1. Introduction

Ageing can be defined as a time progressive decline in intrinsic physiological functions, leading to an increased vulnerability to death. For many years, researchers have suggested that ageing is an inexorable and inevitable process. However, recent studies in several different organisms have shown that the rate of ageing can be significantly delayed by environmental and genetic interventions [1,2] and several transcription factors have been identified to change their expressions along with the senescence [3,4]. These discoveries have demonstrated that the ageing process is at least partially programmed, and that it can be controlled by altering the genetic programme. On the other hand, some types of intervention can decelerate the ageing process across the model organisms. This suggests that the interventions may also be effective in more complex creatures, such as humans.

The first gene that modulates lifespan, age-1, was discovered in *Caenorhabditis elegans* 30 years ago [5,6]. It led to an explosion of ageing studies in this organism. As a result, more genes and pathways that influence worm ageing were identified. The most widely studied pathway is the insulin/IGF-1 signalling (IIS) pathway [7], which has been suggested to regulate lifespan in yeast, worms, fruit flies and rodents [8]. Then several other conserved ageing-related pathways were uncovered, including the kinase target of rapamycin (TOR) [9], AMP kinase (AMPK) [10], sirtuins [11] and germline signalling [12]. These pathways regulate some important fundamental biological functions such as growth, fecundity, stress resistance and nutrient sensing. Additionally, more and more evidence has shown that crosstalk between pathways and the combinations of regulatory signals ultimately influence the systems-level phenotypes of ageing and ageing-associated functional decline [13,14]. Therefore, a systems-level approach is required to integrate various ageing-related factors to explore the fundamental mechanisms of ageing.

Here, we propose a two-step framework to address the issue. The first step is to construct the molecular network of *C. elegans* ageing from biological experiments, involving genes and their regulatory interactions related to ageing. The combinations of activities or expressions of the genes characterize different phenotypic behaviours of ageing. The signal transduction pathways and transcriptional regulations are quantified by the interactions to represent how genes influence each other. The second step is to model the network dynamics by nonlinear coupled differential equations, and then quantify the underlying global landscape and associated flux to explore the dynamics and mechanisms of the ageing process.

The epigenetic landscape concept was first introduced as a qualitative metaphor by Waddington in the 1940s [15], and has been quantified to study the biological processes in different fields, including cell differentiation [16], the cell cycle [17], the immune system [18] and complex diseases such as cancer [19]. Gene networks consist of regulatory elements with inherently stochastic underlying reactions due to the intrinsic and extrinsic sources. The stochastic dynamics of the gene networks can reveal the underlying molecular mechanisms of the biological processes [20]. Every point on the landscape represents a network state; each state is determined by the combination of different gene expressions; and the dynamics of the network becomes the evolution from one state to another. States with local minimum potentials are attractors corresponding to functional phenotypes. Ageing and rejuvenation states can all emerge as attractor basins. The global stability of the long-time dynamics of ageing can be quantified by the barrier height (BH) or transition rate between different attractors. Moreover, global sensitivity analysis based on the landscape topography and kinetics can help to identify the specific genes and regulations which have a significant impact on the ageing functions.

In the following section, we first discuss the assumptions and construction procedure of the ageing network, and then model the network to a set of nonlinear differential equations. We quantified the landscape topography and found interesting dynamical behaviours. Ageing can emerge as a phenotypic state. Furthermore, a rejuvenation state can coexist. We quantified the optimal path and switching rate from the rejuvenation state to the ageing state. We also quantified the optimal path from ageing back to the rejuvenation state, shedding light on possible medical strategies for rejuvenation. A dynamic landscape model was proposed to describe the full picture of the ageing process, because the landscape topography changes as ageing progresses. We further found oscillation dynamics between the rejuvenation state and the ageing state. Under certain conditions, an oscillation behaviour between rejuvenation and ageing has been observed in the long-lived creature *Turritopsis dohrnii* [21,22]. We explored the global stability and kinetic transition between the ageing and rejuvenation states. The key genes and regulations which control the *C. elegans* ageing behaviours are identified through global sensitivity analysis of the underlying landscape topography and kinetics. Our results are consistent with some gene silencing experiments.

## 2. Results

### 2.1. Model network wiring and dynamics

Owing to the lack of full information on gene regulatory interactions, it is still challenging to construct a complete network for ageing. Here we focus on the construction of a core gene regulatory network to describe the basic properties of the *C. elegans* ageing process (figure 1).

We started building the network based on some widely studied pathways such as those of IIS, TOR and AMPK. How these pathways influence *C. elegans* longevity has already been discussed [7]. The IIS pathway plays diverse important roles in regulating the *C. elegans* development, energy metabolism and ageing. The IIS receptor DAF-2 is one of the earliest discovered genes whose mutations can increase the lifespan of *C. elegans* by more than twofold [7]. This lifespan extension requires the activity of transcription factor DAF-16, which induces the expressions of many downstream genes controlling ageing and promoting resistance to stresses [23]. A feedback regulation also exists from DAF-16 to DAF-2 through insulin-like peptides INS-7 [24]. Moreover, AMPK is an energy sensor that plays a central role in maintaining cellular energy homeostasis. DAF-16 synergistically works with AMPK subunits AAK-2 and AAKG-4 to form a positive feedback loop that significantly influences worm longevity [25,26]. In addition, TORC1 is a kinase-containing protein complex which integrates mitogen and nutrient signals into phosphorylation signalling events that regulate cell growth and proliferation. Inhibition of TORC1 extends the lifespan of *C. elegans*, which also leads to the activation of transcription factors DAF-16 and SKN-1 in the IIS pathway. SKN-1 upregulates TORC1 pathway gene expression in a negative feedback loop [27]. TORC1 also affects longevity by promoting mRNA translation through RSKS-1. The longevity phenotype of RSKS-1 requires PHA-4 and HIF-1, which regulate lifespan extension under dietary restriction conditions [9]. In addition, two microRNAs, miR-71 and miR-228, interact with transcription factors SKN-1 and PHA-4. They also mediate the longevity of the dietary-restricted *C. elegans* [28]. The detailed network construction procedures and supporting evidence from the associated network regulations are summarized in the electronic supplementary material, appendix.

After building the network wiring, we modelled the dynamics of the network through a set of differential equations. For each target gene, the sources that determine its activity or expression can be divided into three aspects: activation regulations, inhibition regulations and self-degradations, which can be summarized as: 2.1

where *x _{k}* represents the

*k*th gene expression and d

*x*/d

_{k}*t*represents the evolution of expression of gene

*k*(

*k*

*=*1,2, … ,11) in the network. A total of 11 equations are used to describe the ageing network dynamics. We assumed that the regulatory interactions among ageing genes can be quantitatively modelled by the activations and repressions described by the sigmoid-shaped Hill function [29]. Variables

*x*and

_{i}*x*, respectively, represent the gene expression levels and the corresponding concentrations of gene production that activate or inhibit the expression of gene

_{j}*k*. Parameter

*w*represents the relative strength of every regulation of gene

_{k}*k*, which determines the maximum value of the corresponding Hill function. If there is a total of

*N*regulations targeting gene

_{k}*k*, we assume that all

*w*= 1/

_{k}*N*; this assumption holds all the gene expression values between 0 and 1. The Hill coefficient, labelled as

_{k}*n*in the equation, defines the steepness or cooperativity of each regulation. Parameters

*s*and s

_{ik}*indicate the inflection points of the activation or inhibition regulation terms, which can be used to define the threshold of each regulatory interaction. In addition, self-degradation is quantified as the self-degradation rate*

_{jk}*µ*multiplied by the expression of gene

_{k}*k*. Finally, we sum all the activating and inhibiting Hill functions and self-degradation terms together to describe the evolution dynamics of each gene (see the electronic supplementary material, appendix for the setting of parameters).

Genes have their basal level of expression in nature. Some promoters are not very tightly regulated and show some degree of expressions before the addition of the inducers. However, ageing is such a complex process that most of the parameters are challenging to quantify. We do not have reliable estimates of the leaky expressions. We assume that constant production is relatively small, and use the parameter *w _{k}* to scale the gene expression values between 0 and 1. We simplify the complexity to focus on the relative change of expressions that can be used to characterize the ageing and rejuvenation attractors.

### 2.2. Potential landscape of the ageing network

In our *C. elegans* ageing network, the system state is characterized by an 11–dimensional (11D) vector corresponding to the gene expression levels. The network dynamics can be determined through the evolution of the 11 nonlinear differential equations mentioned above (equation (2.1)). Based on the equations with additional noisy forces characterizing the intrinsic and environmental fluctuations, we can quantify the underlying potential landscape with the Langevin dynamics approach [16,30]. In addition, to visualize the landscape, we need to compress the 11D state space into a two-dimensional (2D) space which can still distinguish the system behaviours and preserve the neighbourhood relationships. The vertical coordinate of the three-dimensional (3D) landscape represents the potential landscape *U*, *U* = −ln(*P*). *P* represents the probability distribution of the state of the network.

The landscape of the *C. elegans* ageing network is shown in figure 2. We projected the 11D state vector into two core and conserved ageing-related transcription factors, DAF-16 and TORC1. DAF-16 is a transcription factor regulating the expression of genes that decelerate the rate of ageing, and the activity of TORC1 influences the transcription and translation of many downstream genes that block lifespan extension [2]. In this figure, we can clearly see that two attractors emerged which have the local minimum potentials and can be distinguished by the relative expressions of DAF-16 and TORC1 on the landscape. For the steady state with lower *U* on the left-hand side in figure 2, all genes with longevity-promoting effects, including DAF-16, SKN-1, AAK-2, AAKG-4, PHA-4 and miR-71, have relatively lower expression levels, and all genes with lifespan-limiting effects, including TORC1, DAF-2, RSKS-1, HIF-1 and miR-228, show higher expression levels than in the other steady state on the right-hand side. The expression or activity levels of all genes at these two steady states are summarized in the electronic supplementary material, table S2.

Based on the relative gene expression levels and their efficacy on the ageing process, we found that the two attractors reveal opposite system behaviours: one performs lifespan-limiting functions, and the other performs longevity-promoting functions. So we identified the two attractors as ageing and rejuvenation attractor states in figure 2. The ageing attractor leads to the decline in physiological functions, while the rejuvenation attractor tends to be against the deteriorations in physiological functions. Moreover, the potential of the ageing attractor is much lower than that of the rejuvenation attractor. This indicates that in this parameter space, from a statistical point of view, most of the worms are enduring a normal ageing process, which represents a progressive decline in physiological functions, but a small fraction of individuals show a distinctively opposite trend. More interestingly, the state of the system can switch from one attractor to another, driven by the interventions from intrinsic fluctuations or external forces of the system. The most probable paths of state switching are identified in the figure by black lines with directional arrows, which are quantified through a path integral method [31]. We also obtained the probabilistic flux [17], which is shown by magenta arrows on the contour plane. This curl flux force leads the optimal paths between the ageing state and the rejuvenation state to deviate from the steepest descent paths and this is irreversible.

### 2.3. Silencing of key genes

The gene silencing experiment is one of the most commonly used approaches to find genes and explore their functions in the *C. elegans* ageing process. In this section, we tested whether our model is capable of quantitatively explaining the gene silencing results, and identified the key genes that have the strongest influences on the *C. elegans* lifespan variations.

The parameter *μ* represents the self-degradation rate, which is initially set to 1 for every gene in the model. We can use the gene's increasing self-degradation rate *μ* to simulate the biological gene silencing operation. We computed the changes in the BH values and mean first passage times between different steady states to quantify the variations in the landscape stability. BH is defined as the difference between the minimum potential in one attractor and the potential of the saddle point from the current attractor to another corresponding attractor. BH_{AR} and BH_{RA} represent, respectively, the BH from ageing to rejuvenation and rejuvenation to ageing. The difference between BH_{AR} and BH_{RA} can be used to measure the relative stability of the ageing and rejuvenation attractors. The mean first passage time (MFPT) is defined as the average transition time from one attractor to another. For example, if the MFPT from ageing to rejuvenation is longer than the MFPT from rejuvenation to ageing, then *C. elegans* is more likely to stay in the ageing attractor state.

We gradually increased the self-degradation rate of each gene from 1.0 to 1.2 to see the corresponding variations in the landscape topography (figure 3*a*). The colour intensities of the heatmap denote the changing magnitude of BH_{AR}–BH_{RA}: the reddish areas indicate that the ageing attractor becomes harder to escape, and the greenish areas mean the rejuvenation attractor increasingly dominated in the system. We can see when accelerating the self-degradation rate (characterized by increasing the value of *μ*), the ageing attractor is more stable for the longevity-promoting genes (top half of figure 3*a*), and the rejuvenation attractor is more stable for the lifespan-limiting genes (bottom half). The results indicate that the silencing of longevity-promoting and lifespan-limiting genes shortens and extends the lifespan of *C. elegans*, respectively, consistent with previous experiments. The variations in the landscape stability provide a new perspective to quantitatively study how genes work together to influence the systems-level behaviour of worm ageing.

Moreover, the variation in the magnitude of landscape stability is different for different genes, and this is represented by the different colour intensities in figure 3*a*. We found that, in these 11 gene nodes, the silencing of DAF-16 and AAK-2 most significantly increases the value of BH_{AR}–BH_{RA}. This indicates that the ageing attractor is much more stable than the rejuvenation attractor. In *C. elegans*, DAF-16 encodes the FOXO transcription factor, which directly regulates the expressions of many ageing-related genes [32]. Mutations in DAF-16 suppress the lifespan extension caused by several genetic manipulations, including daf-2 [7]. AAK-2 is an energy sensor that couples with insulin-like signals to mediate the lifespan of *C. elegans*. Silencing of AAK-2 leads to a 37% decrease in median lifespan in a DAF-2 background and a 54% decrease in lifespan in a wild-type background [33,34]. At the bottom of figure 3*a*, the silencing of TORC1 shows a completely opposite efficacy, which effectively increases the depth of the rejuvenation attractor and decreases the depth of the ageing attractor. TORC1 regulates development, mRNA translation, autophagy and lipid storage in *C. elegans*. Inhibition of TORC1 extends the lifespan in different ways, including losing suppression control of the longevity-promoting genes SKN-1 and DAF-16 [27]. The well-studied lifespan-limiting gene DAF-2 is also shown to significantly decrease the stability of the ageing attractor by increasing its self-degradation rate [35].

We also plot the change in MFPT and BH against *μ* of DAF-16 and TORC1 separately in figure 3*b*; the changes in BH and MFPT show the same trend. The intuitive correspondence between landscape topography and the silencing of DAF-16 and TORC1 is shown in figure 3*c*.

These simulations reveal that our model results are consistent with gene silencing experiments. Furthermore, the model provides a bistable landscape topography to quantitatively explain the mechanism of how genes affect the ageing process, and how the balance of the potentials between the ageing and rejuvenation states determines lifespan.

### 2.4. Global sensitivity analysis of key regulations

We explored the ageing process as a systems-level phenomenon, so both genes and their interactions can influence the dynamics of ageing. However, due to the complexity of the ageing network, it is still very difficult to test the functions of all the regulations with a biological approach. In this section, we perturbed the regulation strengths to find key regulations that have important influences on the *C. elegans* ageing process. We used the relative BH between the ageing and rejuvenation attractors to measure the variations in the system stability.

The parameter *s* in the Hill function determines the inflection point of the sigmoidal function. It defines a threshold of regulator X at which the production Y reaches its half maximal value. When increasing the value of *s* for activated regulation X → Y, more X is needed to generate the same expression of Y. This indicates that the effective strength of this regulation becomes weaker. The function of *s* is similar in the case of inhibited regulation. For the gene expression value which is limited between 0 and 1 in the model, we defined the regulation strength as 1 − *s*. The initial values of all regulation strengths are set to 0.5, then we perturbed every regulation strength from 0.3 to 0.7 to see the variation in the landscape stability.

The result is shown in figure 4, where we used the change in BH_{AR}–BH_{RA} to measure the change in the landscape topography. The colour change from green to red indicates the gradual switching of the dominant attractor from rejuvenation to ageing. From this heatmap, we can clearly find three clusters of the 27 regulations divided by their different effects on the landscape stability. The top cluster of regulations shows that, when increasing their regulation strengths, the rejuvenation attractor will be more and more stable. The middle cluster has little impact on the landscape stability. The bottom cluster has an opposite trend to the top cluster. Strengthening of these regulations inhibits lifespan extension. We checked the regulations which show a significant impact on the landscape stability in the top and bottom clusters, and found some of their functions have already been confirmed by biological experiments. For example, the strengthening of TORC1 → RSKS-1 dramatically increases the stability of the ageing attractor in figure 4, which is consistent with the evidence that RSKS-1, as a key mRNA regulator downstream of TORC1, plays an important role in lifespan extension by inhibiting TORC1 in many organisms [9]. Moreover, the most widely studied path from DAF-2 to DAF-16 was also identified in our global sensitivity analysis [7]. In addition, we have identified several key regulations including AAK-2 TORC1 and AAK-2 HIF-1. Increasing their regulation strengths suppresses the ageing process. For SKN-1 → TORC1, increasing its regulation strength accelerates the rate of ageing. The functions of these regulations on worm ageing have not been clearly indicated before, and need further experimental verification. The predictions of the key regulations from our global sensitivity analysis are summarized in table 1.

In summary, we quantified regulation strengths and simulated landscape stability variations in response to perturbations in regulation strengths, which provides a quantitative way to uncover the key regulations that have a strong influence on *C. elegans* ageing. These predictions can be validated by the experiments.

### 2.5. Dynamic landscape of the ageing process

The potential landscape above can be used to describe the behaviour of *C. elegans* ageing. Perturbations of the network will change the relative depth between the ageing and rejuvenation attractors, which determines the rate of ageing. However, this model is only a time slice of the natural ageing process. To understand the whole process of ageing, we need to explore a dynamic landscape which continuously changes with ageing.

Several transcription factors have been identified whose expressions increase or decrease during the natural ageing process. The associated studies show that these regulators significantly influence the lifespan of *C. elegans* [3,4]. Why do their expressions change synchronously with ageing? What is the origin of initiation and promotion for the ageing process? Ageing theories such as accumulation of damage, hyperfunction of the development pathways and others have been suggested to explain these questions [1,44]. It seems none of them are fully satisfactory at present. Nonetheless, the variations in regulatory interactions probably play an important role in ageing because they directly control gene expressions.

To visualize the dynamic ageing model, we assumed that the regulation strength of DAF-16 TORC1 gradually decreases with time. We use the equation d*S*/d*t* = −*λS* to model the continuous change of regulation strength parameter *S* in our model. *λ =* 10^{−7} represents the rate of change. We assume the change of parameter *S* is synchronous with the rate of ageing, which is much slower than the dynamics of gene expression in our model. By this means, the system has enough time to relax to the ageing and rejuvenation attractors at each time point of the whole ageing process.

In figure 5, the changes in regulation strength lead to a continuous transformation of the landscape topography. When the worm is young, the rejuvenation attractor dominates the landscape, and the physiological functions are perfect and adaptable. As time progresses, with the changes in regulation strength, the ageing attractor becomes more and more stable, which leads to a progressive decline of the system functions and an increased vulnerability to death. There are two time scales in the ageing process. One is the changes in the gene regulation strength and the other is the switching times between the rejuvenation and ageing attractors. This provides a dynamic landscape as shown in figure 5. The ageing process can be described by the landscape changes, which are determined by the change in the regulation strength, which is likely to be caused by the accumulation of damage. The top of the landscape before significant regulation changes with only the rejuvenation state can be defined as the birth (young) state, and the landscape at the bottom after significant regulation changes (deterioration) from damage accumulations can be defined as the death (old) state. A slower rate of changes in regulation strength leads to a longer average lifespan from birth to death (attractors), and a faster rate of change leads to a shorter average lifespan.

Many experiments have supported that the ageing process can be significantly delayed [1,2,44,45]. The rejuvenation and ageing states can coexist only at certain stages of the lifespan (see figure 5 at certain regulation strength), while usually the switching rate is faster than the ageing rate (change rate of regulation strength). On the one hand, the slower switching rate is due to the higher barrier between the rejuvenation and ageing attractors. This makes it harder to communicate between these two attractors and can lead to one or a few distinct paths. On the other hand, the faster switching rate is due to the lower barrier between the rejuvenation and ageing attractors. Then it becomes easier to communicate between these two attractors. This can lead to multiple paths. The lifespan is expected to have a narrower or wider distribution depending on how many possible states or routes are available.

From daily experience and experiments, ageing seems to be a one-way process [46]. We believe that, in certain gene disruptions and stress experiments, the perturbations do not significantly change the process of regulation deterioration. Or, in other words, there are no significant changes in the BH from the ageing attractor at the bottom to the birth attractor at the top along the changes in regulations shown in figure 5, which is in general high by the genetic constructions. The reason for rescaling the lifespan distribution without significant reversal of the ageing process in some experiments can be that the barrier from birth to death and the relative stability between birth and death might be changing while the reverse barrier from death to birth is not changing significantly or becomes even larger. This effectively makes it a one-way process from the birth to the death attractor without the reversal. This naturally leads to a lifespan distribution scaling during the ageing process. Therefore, the lifespan distribution does not change its shape significantly because it is effectively a single direction downhill process without the reversal. This is our landscape perspective.

A further question is can we reverse the ageing process to make the worm rejuvenate from old to young? From the dynamic landscape model, if certain key genes and regulations do not deteriorate significantly, then the bottom level of figure 5 will rise up. As a result, there can be coexistence of rejuvenation and ageing attractors or even a dominant rejuvenation attractor. This provides the opportunity for the emergence of a stable rejuvenation state, leading to the possibility of ageing reversal. It is, therefore, important to find those key genes and regulations critical for the ageing and ageing reversal process, for the purpose of control and medical applications.

From the dynamic landscape model, we quantified the optimal path from the young attractor state to the old attractor state as shown in figure 5. Furthermore, we also found that an optimal path from the old attractor state to the young attractor state exists and can be quantified (figure 5). In this dynamic landscape picture, the process of ageing reversal can coexist with the normal ageing process. However, the ageing path and rejuvenation path are irreversible and different due to the presence of the flux breaking the detailed balance. Furthermore, the chance or the probability of the two processes can be very different. We calculated the probability of the optimal paths from young to old as *P*_{Y→O} and from old to young as *P*_{O→Y} from a path integral formulation based on the ageing network wiring in our model. The ratio of the probability becomes *P*_{Y→O}/*P*_{O→Y} = 34.61. This means that rejuvenation as an ageing reversal process is much harder and has a significantly lower chance of occurring than the normal ageing process. We believe this difference can be a reason why rejuvenation or ageing reversal has not been observed in experiments. Through genetic and environmental manipulations, the underlying landscape of ageing can be altered. In this way, one hopes to increase the chance to realize ageing reversal.

In [47], the authors reviewed some experiments showing that the rejuvenation process can emerge at the cell or tissue level. They also suggested that reprogramming from differentiated adult cells to induced pluripotent stem cells (iPSCs) can be partially seen as resetting of the ageing clock, because the iPSCs resemble embryonic stem cells (ESCs), and ESCs can generate all cell types and support complete cell fate development. In this respect, ageing is like a stem cell developmental process while reversal of ageing is like a reprogramming process. Finding a reversal of ageing from old to young has some resemblance to finding a reprogramming path from a differentiated cell back to a pluripotent stem cell. Although the cell differentiation and reprogramming are not equivalent to the ageing and rejuvenation process, the analogy here can provide new insights to understand and explore these questions.

### 2.6. Oscillation dynamics between ageing and rejuvenation

When we vary the regulation strengths, in the range of (0, 1), an oscillation dynamics can emerge in a certain parameter space (electronic supplementary material, appendix and table S3). The transitions between the oscillation and bistable dynamics are determined by the regulations: SKN-1 DAF-2, SKN-1 → miR-71 and miR-71 DAF-2, which are direct or indirect interactions from SKN-1 to DAF-2. The changes in landscape topography are shown in figure 6*a*–*c*; RS represents the regulation strength of SKN-1 DAF-2. The landscape shows oscillation dynamics with a Mexican hat shape when RS is 0.75 (figure 6*a*). The two relatively deeper regions coloured with blue on the oscillation ring correspond to the ageing and rejuvenation states. The system state rotates along the oscillation ring valley around the centre hill of the Mexican hat. When the regulation strength RS is increased to 0.8 (figure 6*b*), we can see that the barrier between the two attractors becomes higher, which means the switching between ageing and rejuvenation is harder. The system has then switched from oscillation to bistable dynamics. When RS = 0.85 (figure 6*c*), the landscape clearly reflects a bistable phase with the ageing attractor and rejuvenation attractor as discussed before.

It is very interesting to try to explain the oscillation dynamics for ageing and rejuvenation. If the oscillation dynamics exists in real life, it means a creature can have the normal ageing process. But instead of dying eventually, it will gradually rejuvenate from old back to young. This cyclic process can be repeated as naturally as breathing. Although this behaviour has not been found in *C. elegans* at present, it has been observed in the long-lived creature *T. dohrnii*, which is the first metazoan known that can repetitively rejuvenate from the medusae stage back to the polyp stage, thus escaping death and achieving potential immortality [21,22]. If the assumption of oscillation dynamics can be generalized to more complex biological systems, it will give new perspectives to understand and control the ageing process.

In addition, we calculated the entropy production rate (EPR) for the phase transition from oscillation to bistability by increasing the regulation strength of SKN-1 DAF-2 (figure 6*d*). The EPR represents the total entropy variations. Decreasing the EPR means the system costs less energy to maintain (see electronic supplementary material, appendix for detailed definition of EPR). From figure 6*d*, we can see that the EPR remains at a high level when the system shows oscillation dynamics. When the system switches from oscillation to bistability, the EPR sharply decreases and then stays at a low level. This implies that the oscillation phase costs much more energy to maintain than the bistable attractors. We also calculated the flux integral and period of the oscillation dynamics (figure 6*e*). The flux integral shows the opposite trend to the period; this is because a larger flux indicates more energy or a stronger driving force in the oscillation, which leads to a faster speed and shorter period of oscillations.

From the above discussion, we can see that, to achieve oscillation between ageing and rejuvenation or immortality, more energy has to be pumped into the system. The oscillation speed between the ageing and rejuvenation states is accelerated by the energy pump. The energy pump to the system breaks the detailed balance and leads the system to be in a non-equilibrium state. In fact, the energy pump manifests itself as a driving force of the dynamics in addition to the gradient of the landscape in the form of the curl flux. A larger flux is correlated with a larger deviation from equilibrium and bigger energy costs or bigger EPR. While the landscape helps to stabilize the ageing and rejuvenation states, it is the curl flux that drives the dynamics and accelerates the switching between ageing and rejuvenation.

## 3. Discussion

In this study, we constructed the gene regulatory network of *C. elegans* ageing by integrating evidence from biological experiments, and modelled the system dynamics through 11 nonlinear coupled differential equations. We quantified the potential landscape of worm ageing by the Langevin dynamics method. The model shows bistable dynamics and the two steady states can be characterized as ageing and rejuvenation by the relative activities of genes with longevity-promoting and lifespan-limiting functions. While the ageing attractor leads to a decline in physiological functions, the rejuvenation attractor tends to be against the deteriorations of physiological functions. We proposed that the balance of depths between the ageing and rejuvenation attractors determines the speed of *C. elegans* senescence. Moreover, this balance can be disturbed by internal fluctuations or external interventions of the system.

The BH and MFPT between the ageing and rejuvenation attractors can be used to quantify the landscape topography. The increase in the self-degradation rate *µ* will increase the relative stability of the ageing attractor for longevity-promoting genes on the one hand, and increase the relative depth of the rejuvenation attractor for lifespan-limiting genes on the other hand. These model behaviours can give quantitative explanations of some gene silencing experiments. We also found the key genes and regulations for ageing and rejuvenation by perturbing the regulation strengths to see their impacts on the landscape stability. To explore the whole picture of the ageing process, we explored the dynamic landscape, which describes the continuous transformations of the landscape topography along the natural ageing process. Importantly, we found that oscillation dynamics can emerge between the rejuvenation state and the ageing state, and similar behaviour has been observed in the long-lived creature *T. dohrnii*.

The mechanism of ageing is a topic under intense discussions. We believe that the ageing process is at least partially programmed. This is because modifications of ageing-related genes can greatly elongate the lifespan and several transcription factors have been identified to change their expressions along with senescence. These studies support that the genes and the regulatory interactions among them play important roles in worm ageing.

It is still challenging to construct a complete ageing network due to the lack of information on gene regulatory interactions. Therefore, we have focused on the most widely studied genes and regulations to build a core gene regulatory network of *C. elegans* ageing. The IIS, TOR and AMPK are the most well-studied pathways related to ageing. Manipulations of the genes in these pathways have already been found to significantly influence the ageing process. Transcription factors such as DAF-16 and TORC1 have many downstream target genes, like elt-3, etc. We do not consider them all, particularly those without feedback in the network. Furthermore, there are some genes that have been identified to be differentially expressed between young and old worms. However, the underlying regulatory interactions of these transcription factors or genes with others in the ageing network are not very clear. Therefore, it is still challenging to use these results directly for network construction. Further detailed investigations are needed to explore the more complete regulatory interactions among the genes in the ageing process.

Genes performing longevity-promoting and lifespan-limiting functions are coloured as green and red in figure 1. There are three exceptional regulations, SKN-1 → TORC1, PHA-4 → miR-228 and miR-71 PHA-4, that either activate other genes in their own group or repress genes in the other group in our network. The expressions of miR-228 and PHA-4 are used together with other genes to distinguish the rejuvenation and ageing states, which can be found in the electronic supplementary material, table S2. Their expression values are also changed in our perturbation studies. We perturbed the regulation strengths, which showed that the regulation SKN-1 → TORC1 significantly influences the landscape topography but the other two regulations do not have obvious effects (figure 4). Additionally, we separately removed the three regulations from the network; the result is consistent with the global sensitivity analysis of the regulation strength (electronic supplementary material, figure S1). This shows the exceptional regulation SKN-1 → TORC1 is important for the ageing process while the other two are less relevant even though their expression values are changed in ageing and rejuvenation upon perturbations (the underlying landscape topography in terms of the ageing and rejuvenation attractors is not changed significantly). While mutual repressions often lead to a toggle switch, the combination of the activation–repression regulations can limit cycle oscillations.

In the landscape picture, the attractors of ageing and rejuvenation as well as the oscillations are all the result of the regulatory interactions in the ageing network. Otherwise there would not be such distinct attractors or ageing and rejuvenation states. Without interactions, every state will have the same negligible weight and is not observable. Accumulated damage can cause changes in the genes and gene regulations in the network. In other words, damage leads to changes in the network wiring and, therefore, changes in the underlying landscape topography. Therefore, we proposed a dynamic landscape picture to describe the whole process of ageing. With the deterioration of the key regulations, the landscape changes from a dominant rejuvenation attractor at the beginning to a dominant ageing attractor at the end, which represents the natural ageing process from young to old. In this picture, the process of ageing reversal can coexist with the normal ageing process. We quantified the optimal paths for these two processes. The ageing and rejuvenation paths are irreversible and different from each other due to the presence of the flux breaking the detailed balance. Furthermore, the probability of a rejuvenation or ageing reversal path is much lower than the probability of a normal ageing path. This means that the rejuvenation or ageing reversal process has much less chance of being realized than the normal ageing process. This could be a reason why ageing reversal still has not been seen in *C. elegans*. Although ageing reversal is difficult to realize in complex organisms, some experiments [47] have already shown that rejuvenation can emerge at the cell or tissue level.

The variations in network wiring can also limit cycle oscillations in certain parameter spaces, as in our *T. dohrnii* example. This seems to be the only example of a living creature being immortal with oscillations between rejuvenation and ageing, which can be quantified in our Mexican hat-shaped landscape and flux flow picture [17,48]. The landscape attracts the system to the oscillation ring while the flux drives the oscillation along the ring. In this way, the rejuvenation and ageing states can stably oscillate among each other, reaching immortality.

In this work, we have provided a new systems-level approach based on landscape theory to integrate various ageing-related factors to explore the fundamental mechanism, and to predict key genes and interactions in the *C. elegans* ageing process. This approach is general and can be used for studies of other ageing systems.

## 4. Material and methods

The ageing process in real life is affected by the intrinsic or external fluctuations of the system. So a Langevin equation is appropriate to describe the corresponding stochastic time evolution of gene expression dynamics in the model. These stochastic differential equations are as follows: 4.1

where **x** is a vector of the concentrations of genes or gene products. **F**(**x**) is the driving force of the gene regulations formulated in equation (2.1). *η*(**x**, *t*) is the noise term that has a Gaussian probability distribution with correlation function 〈*η*(**x**, *t*),*η*(**x**, *t′*)〉 = 2*Dδ*(*t* − *t*′), where *D* is the diffusion coefficient characterizing the strength of the fluctuations and *δ* is the Dirac delta function. We can quantify the global probability distribution *P* for the state space. For a simple example of a 2D system, we can divide the 2D state space into 200 × 200 grids. Following a long-time Langevin dynamic trajectory by solving the stochastic differential equations in equation (4.1) iteratively, we obtained the frequency of the state of the system in each grid. Then the frequencies can be converted to a probability distribution, and the underlying landscape can be quantified from *U* = −ln(*P*).

The individual stochastic trajectory is unpredictable. We will follow the evolution of the probability distribution to describe the system behaviours of the ageing process. This evolution of probability distribution obeys the Fokker–Planck equation [17,48,49], 4.2and 4.3

Equation (4.2) denotes the evolution of probability, the variation of local probability *P*(**x**,*t*) is equal to the divergence of flux **J**(**x**, *t*) in and out of the region. The divergence of steady-state flux should be zero, i.e. ∇ . **J**_{ss}(**x**) = 0 due to the steady-state condition. However, the flux itself is not necessarily equal to zero. In fact, it can rotate as a curl. The probability flux can be calculated from equation (4.3), which leads the optimal paths not to follow the one with the steepest gradient. The flux drives the system state to rotate along the oscillating ring. The flux can be used to evaluate the degree of detailed balance breaking and quantify the non-equilibriumness. It is also the origin of the energy dissipation quantified by the EPR. The dynamics of the non-equilibrium system is determined by both the landscape (gradient) and the flux [17,48,49], **F**(**x**) = −*D*∇*U*(**x**) + **J**_{ss}(**x**)/*P*_{ss}(**x**), where the landscape *U* is defined as the negative of the logarithm of the steady-state probability: *U* = −ln(*P*_{ss}).

The optimal paths indicate the most probable paths when the system switches from one state to another. In this work, the optimal paths, as shown in figure 2, represent how the worm's state switches between the ageing and rejuvenation states by the inductions, or the intrinsic or environmental fluctuations, of the system. The quantification of the optimal paths is critical in uncovering the underlying mechanism of this biological process, which is obtained through the path integral approach [16,31].

We formulate the probability of the system dynamics to switch from **x**_{initial} at time 0 to **x** _{final} at time *t* with the path integral as:
where the integral over *D***x** indicates the sum over all the possible trajectories starting at state **x** _{initial} at time 0 and ending at state **x** _{final} at time *t*. **F**(**x**) represents the driving forces originating from the regulatory interactions. **D**(**x**) represents the strength of the fluctuations characterized by the diffusion coefficient matrix. *S*(**x**) and *L*(**x**(*t*)) refer to the action and Lagrangian of the associated path. This quantifies the weight of the path.

An optimal path is the path with the largest weight. As the contribution of each path is exponentially weighted by the negative of the corresponding action, we can identify the optimal paths through minimizing the action and ignore the subleading contributions.

## Authors' contributions

J.W. designed the research; L.Z. and J.W. performed the research, analysed the data and wrote the paper.

## Competing interests

We declare no competing interests.

## Funding

We are grateful for support from NSFC grant no. 91430217 and MOST, China, grant nos. 2016YFA0203200 and 2013YQ170585. J.W. is grateful for support from grant no. NSF-PHY 76066.

## Acknowledgment

L.Z. thanks Dr Han Yan and Kun Zhang for discussions and comments.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3577184.

- Received May 30, 2016.
- Accepted November 7, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.