Systems biology offers valuable insights into aging by integrating experimental data with mathematical models and bioinformatics tools. Long-lived mutants of C. elegans, particularly clk-1, have provided extensive data on aging mechanisms. The clk-1 gene, which encodes a ubiquitin precursor, shows a pleiotropic phenotype characterized by slow behavior, high mitochondrial ROS levels, autophagy, and metabolic changes. However, the link between these changes and lifespan extension remains unclear. Using a Boolean network, we modeled genetic interactions and derived differential equations for a continuous approach. Our results highlight that aak-2 (AMPK) is crucial for clk-1 lifespan extension owing to its role in stress response regulation. We introduced a health index based on the attrition of neuromuscular behaviors to assess the health of various strains. Our findings suggest that while stress responses may enhance lifespan, overall health is determined by the extent of the damage.

Literature search to infer the clk-1 regulatory network

The clk-1 regulatory network was inferred through an extensive literature and database search. The primary sources of information included WormBase, WormBook, and PubMed. We systematically explored gene–gene interactions involving clk-1 and its associated pathways using WormBase as the main resource. This initial step yielded a preliminary network draft describing the expression changes in genes such as atfs-1, skn-1, taf-4, aak-2, and hsp-60.

Subsequently, we expanded our search to include the downstream and regulatory targets of these initially identified genes. This expansion was guided by experimentally validated interactions reported in the literature, which enabled the incorporation of additional key components, including mTOR, atgl-1, hif-1, and CREB. Through this iterative approach, the network was progressively refined to capture the interconnections between the diverse cellular processes.

To further contextualize the regulatory interactions, we performed a focused search for genes and molecular elements implicated in stress and mitochondrial stress responses, particularly those associated with mitochondrial mutations and reactive oxygen species (ROS) accumulation during aging. This step enabled the integration of pathways related to lipid metabolism, mitophagy, and the electron transport chain. The final curated network comprised 29 interconnected components, forming the foundation of the dynamic analysis presented in this study.

Construction of the boolean model for the clk-1 regulatory network

To build a comprehensive model for the clk-1 regulatory network, we thoroughly researched the existing literature to gather all documented interactions relevant to C. elegans clk-1 mutant. We were able to include the complete set of interactions to ensure that our regulatory network accurately reflects the current understanding of this organism model. In the network, each node symbolizes a distinct molecular component, such as genes, proteins, or small molecules, integral to the clk-1 network. The edges that connect these nodes denote the regulatory influence — positive or negative — one component exerts on another.

To translate experimental observations from the literature into Boolean regulatory rules, we first classified each interaction as either positive or negative, depending on whether the effector enhanced or repressed the expression or activity of its target. In Boolean formalism, positive interactions are represented by the identifier of the activating effector alone, whereas negative interactions are expressed using the negation operator NOT (or “!”).

When multiple effectors regulate a target, the formulation of the Boolean rule is determined by empirical evidence of combinatorial regulation. Specifically, when a target requires the concurrent action of two or more activators and the exclusion of specific inhibitors, the logical operators AND (“&”) and NOT are combined to capture these dependencies. Conversely, when activation can occur through alternative effectors acting independently, the OR (“|”) operator is used. Because multiple logical configurations can be consistent with the same set of interactions, several candidate rule sets were systematically tested across alternative scenarios, and the final rule set was selected as the one that best reproduced the experimentally observed conditions and behaviors.

Dynamic simulations were performed using a synchronous update scheme in which all nodes were updated simultaneously at each discrete time step according to their Boolean functions. This approach provides computational efficiency while preserving the attractor structure of the asynchronous formulation. The additional dynamic variability sometimes observed under asynchronous updates arises from the compositional structure of the logical functions; however, the synchronous implementation substantially reduces the computational cost. In particular, simulating a synchronous Boolean model comprising 29 nodes requires significantly fewer operations than its asynchronous counterpart, making it advantageous for exploratory analysis and parameter testing.

We used the Cytoscape software to visualize the biological network depicted in the graph. Through Cytoscape, we represented the clk-1 regulatory graph in a user-friendly format, facilitating a simplified exploration of its structure and dynamics. Furthermore, we enhanced the precision and predictive capability of the model using the R package BoolNet. This package was the computational framework for simulating the clk-1 network and identifying its stable states. These are crucial in representing network configurations that may endure over time, potentially correlating with specific phenotypes.

Continuous model formulation

To derive the continuous counterpart of the Boolean network, each logical rule was systematically converted into a Hill-type function that captured the cumulative regulatory influence of all relevant effectors on a specific target variable. This methodology preserves the logical structure inherent in Boolean formulation while introducing continuous dynamics through nonlinear activation functions. Each regulatory interaction is characterized by its own set of parameters, enabling precise control of the activation strength and decay rates. In the resulting system of ordinary differential equations (ODEs), the rate of change of each variable is defined by two components: a sigmoidal activation term representing synthesis or activation dynamics governed by the Hill function and a linear decay term accounting for degradation or deactivation.

We utilized the R package BoolNetPerturb to simulate the effects of various perturbations on the clk-1 network, gaining valuable insights into its resilience and adaptability. Furthermore, the deSolve R package was used to delve into the network's continuous dynamic behaviors, leading to a more comprehensive understanding of its temporal dynamics under different conditions.

C. elegans strains and culture

Nematodes were cultured in standard conditions on Nematode Growth Medium (NGM) plates and fed with the Escherichia coli strain OP50 at 20°C temperature. In the experiments, Bristol N2 was used as the wild-type strain, while the mutant strains used were aak-2(ok524) and clk-1(qm30). These strains were provided by the Caenorhabditis Genetics Center (CGC), which is funded by the NIH Office of Research Infrastructure Programs (P40 OD010440). To create a double mutant, aak-2(ok524) males were crossed with clk-1(qm30) hermaphrodites. F2 progeny were then separated on individual plates to identify the double mutants using the polymerase chain reaction (PCR) technique. To identify the aak-2 mutation, the following primer sequences were utilized: forward TTCCTGGCAACACCATAAGC and reverse CTCCAGAAAGTCTGGAGTTG. On the other hand, to detect the clk-1 mutation, the primer sequences used were forward TGTCGGTTCAGCACTTCTAG and reverse AGTATTGTCCGTGTCAGGAC. This method allowed for precise identification of the desired genetic configurations.

Measurement of neuromuscular behaviors

Our study assessed neuromuscular behaviors using the methodologies outlined by Hart, focusing on a trio of specific behaviors: pharyngeal pumping, swimming, and defecation. These behaviors were selected as indicators of neuromuscular activity and coordination, offering insights into the nematodes' physiological health and functional status under investigation.

The experimental protocol started with nematodes that had reached one day of adulthood. A cohort of 20 worms was selected to ensure age and developmental stage uniformity. This selection criterion was crucial for minimizing variability in neuromuscular function that might arise from differences in maturation. The behaviors of these worms were then systematically observed and recorded over nine days, with assessments conducted at 48-hour intervals. This longitudinal approach allowed us to capture possible changes in neuromuscular activity over time, providing a comprehensive overview of each behavior.

We measured pharyngeal pumping and swimming behaviors three times for each worm at every observation point. Pharyngeal pumping, a critical indicator of feeding behavior and metabolic activity, was quantified by counting the number of pharyngeal contractions within a predefined time frame (30 seconds). Swimming behavior, reflecting motor coordination and muscle function, was evaluated based on the worms' frequency of tail waving in an aquatic environment, within a predefined time frame (60 seconds). On the other hand, the defecation behavior of each worm was measured only once by recording the time between defecation cycles, characterized by periodic expulsion events, within a 10-minute time frame, to gain insight into the regulation of the gastrointestinal and neuromuscular systems.

Survival analysis

Our survival analysis protocol was designed following the procedures delineated by Park et al. to ensure the reliability and reproducibility of our results. This approach allowed for the precise monitoring of nematode longevity under controlled conditions. The experimental setup began with selecting nematodes at the L4 larval stage, an easily identifiable stage in their development. These selected individuals were then carefully transferred onto fresh Nematode Growth Medium (NGM) plates. This transfer was performed to standardize the starting conditions for each nematode, thereby reducing variability in survival outcomes that might arise from differences in environmental factors or developmental stages at the onset of the assay.

Twenty-four hours post-transfer, marking the transition of the nematodes to one-day adults, the survival assays officially started, with this time point designated as t = 1. Starting the assays at the adult stage ensured that the survival data captured pertained exclusively to the adult phase of the nematode life cycle, facilitating the analysis on adult longevity.

Survival was assessed every 24 hours. Nematodes that deviated from the experimental parameters, such as those that crawled off the plate, exhibited bagging (a condition where eggs hatch inside the mother, leading to her death), or exploded, were censored at the time these events were observed. The censoring of these individuals was critical for the integrity of the survival analysis, as it accounted for non-standard mortality causes that do not directly relate to the natural lifespan of the organism. The nematodes were transferred to new NGM plates every 48 hours to maintain optimal conditions and minimize confounding factors such as overcrowding or depletion of food resources. This bi-daily transfer protocol was pivotal in ensuring that environmental degradation did not influence the observed survival outcomes over time. Additionally, this strategy prevented the inadvertent inclusion of new adults from potential reproduction, thus maintaining the purity of the adult cohort under study.

Health index estimation

In our study, we developed a nuanced approach to estimate the health index of nematodes, leveraging the variables derived from our behavioral analyses. This estimation process was designed to encompass individual and aggregated behavioral metrics, providing a holistic view of the nematodes' physiological state over time.

Aggregation of behavioral data into a health index

To construct the health index, we first relativize all the behavioral data by converting it to a relative scale. This standardization facilitated the comparison and aggregation of diverse behaviors on a standard metric, enhancing the interpretability of our analysis. Data from different behaviors were cumulatively added for the aggregation process to form a composite score. Notably, in the case of defecation—a behavior that increases its length over time—its relative metric was subtracted from the composite of the other behaviors. This adjustment was crucial for accurately reflecting the contribution of each behavior to the overall health index, acknowledging that not all behaviors equally signify healthiness.

Symbolic regression for health index estimation

The initial method used for estimating the health index was symbolic regression, executed through the software Eureqa (38). Symbolic regression diverges from traditional regression by searching for the mathematical expression that best fits the relationship within the data rather than fitting the data to a predefined model. Through this process, we first identified the optimal function that exhibited the highest fidelity to the behavioral data patterns, and then we calculated its definite integral over the observed period. This integral, representing the area under the curve generated by the function, served as a quantitative measure of the cumulative health index across the lifespan of the nematodes.

Direct area calculation from behavioral data

Additionally to the symbolic regression approach, we employed a direct method to estimate the health index by calculating the area under the data curve. This method drew lines through the behavioral data's first, second, or third quartiles for each measurement day. By connecting these points across the observation period, we delineated polygons whose boundaries were defined by the quartile lines and the temporal axis representing the consecutive days. The area enclosed by these polygons was then computed, offering an alternative measure of the health index that is directly derived from the raw behavioral data.

Statistical analysis

In our study, the statistical analysis assessed the effects of genetic variation and temporal dynamics on neuromuscular behaviors and survival across different nematode strains. For the variables of pharyngeal pumping frequency, defecation cycle speed, and waving frequency during swimming, we performed a two-way mixed Analysis of Variance (ANOVA) with nematode strain (fixed factor), time (as random factor), and their interaction as sources of variation. This analysis facilitated the discernment of strain-specific behavioral profiles and their progression over time, allowing us to pinpoint any significant deviations attributable to genetic differences or temporal changes. The test allows us to evaluate the interaction effect between strain and time, providing insights into whether the impact of genetic variation on neuromuscular behavior changes over the lifespan of the nematodes.

The Cox proportional hazards regression model was employed to analyze survival data. The Cox regression analysis was performed to verify the proportionality of the hazard assumption, employing both graphical diagnostics and statistical tests. This step ensured the validity of the model's findings. The survival curves generated from this model were then statistically compared, with the resulting p-values providing a basis for assessing the significance of differences observed between strains; all statistical analyses were conducted using the R programming language.