Systems modelling predicts chronic inflammation and genomic instability prevent effective mitochondrial regulation during biological ageing

The regulation of mitochondrial turnover under conditions of stress occurs partly through the AMPK-NAD+-PGC1α-SIRT1 signalling pathway. This pathway can be affected by both genomic instability and chronic inflammation since these will result in an increased rate of NAD+ degradation through PARP1 and CD38 respectively. In this work we develop a computational model of this signalling pathway, calibrating and validating it against experimental data. The computational model is used to study mitochondrial turnover under conditions of stress and how it is affected by genomic instability, chronic inflammation and biological ageing in general. We report that the AMPK-NAD+-PGC1α-SIRT1 signalling pathway becomes less responsive with age and that this can prime for the accumulation of dysfunctional mitochondria.


Introduction
Genomic instability and mitochondrial dysfunction are two robust hallmarks of biological ageing (Fakouri et al., 2019;Lopez-Otin et al., 2013), and can influence each other during the ageing process (Fakouri et al., 2019). DNA repair is a key process that determines the extent of genomic instability that develops with age. One of the substrates needed for effective DNA repair is nicotinamide adenine dinucleotide (NAD + ), which is a substrate in poly-ADP ribosylation (PARylation) reactions mediated by Poly ADP-ribose polymerase (PARP) enzymes in the resolution of all types of DNA lesions (Fang et al., 2017;Wei and Yu, 2016). However, NAD + is also needed in energy metabolism and energy sensing (Canto et al., 2015). In fact, NAD + is part of a mito-nuclear signalling axis linking genomic stability with mitochondrial turnover and the energetic state of the cell (Canto et al., 2015).
The main molecular players in this pathway are AMP-activated protein kinase (AMPK), Peroxisome proliferator-activated receptor gamma coactivator 1-alpha (PGC1α), NAD + , Sirtuin 1 (SIRT1) and PARP1. Previous work on the interplay between DNA damage and mitochondrial function has found antagonistic effects between PARP1 activity responding to DNA damage levels and mitochondrial turnover responding to deacetylated PGC1α levels (Canto et al., 2015). This suggests that the regulatory circuit mediates dynamic trade-offs between DNA repair activity and other cellular functions.
The importance of mito-nuclear communication is highlighted by studies where NAD + supplementation has robustly increased both lifespan and healthspan in a variety of organisms (Balan et al., 2008;Belenky et al., 2007;Cerutti et al., 2014;Fang et al., 2016;Guan et al., 2017;Mouchiroud et al., 2013;Yaku et al., 2018). Conversely, low NAD + levels have been associated with ageing pathology and mitochondrial dysfunction (Gomes et al., 2013;Zhang et al., 2016;Zhu et al., 2015). Interestingly, a mito-nuclear communication loop mediated by p21 and reactive oxygen species (ROS) has also been established to drive cellular senescence (Passos et al., 2010). All of this suggests a key role of mito-nuclear communication in maintaining a balance between nuclear and mitochondrial maintenance with age (Fakouri et al., 2019).
Changes in the NAD + -mediated mito-nuclear communication axis have been reported with age (Mendelsohn and Larrick, 2017). Sirtuin activity, NAD + levels and AMPK responsiveness tend to decrease with age (Camacho-Pereira et al., 2016;Mendelsohn and Larrick, 2017;Salminen et al., 2016;Yaku et al., 2018) whilst DNA damage tends to increase with age (Freitas and de Magalhaes, 2011). NAD + depletion also contributes to cellular senescence and the associated secretory phenotype (Wiley and Campisi, 2021) as well as cell death (Sedlackova et al., 2020).
In this study, we investigated whether genomic instability and mitochondrial dysfunction may be linked through changes in mitonuclear communication using a systems modelling approach to explore how age-related alterations may affect the functionality of the NAD + -mediated mito-nuclear communication axis.

Model development
We developed a computational model of the AMPK-PGC1α-NAD + -SIRT1 axis in COPASI (Hoops et al., 2006) using coupled ordinary differential equations (ODEs) that simulate changes in molecule abundances over time. The time units of the model correspond to hours and the molecule and volume units of the model correspond to arbitrary units (AU).
The core set of interactions within the pathway are illustrated in Fig. 1. Upon AMPK activation via phosphorylation, there is rapid PGC1α phosphorylation and a slow PGC1α deacetylation through an increase in NAD + levels.
Interactions denoting (de)phosphorylation and (de)acetylation reactions are mechanistic, since they represent well-defined biochemical processes. However, the interaction representing an increase in NAD + levels caused by an increased AMPK activity is phenomenological, representing many underlying processes such as transcription factor activation, binding to a promoter, mRNA transcription and translation (Cantó et al., 2009).
To explore how the signalling pathway responds to relevant conditions requires further expansion as illustrated in Fig. 2. NAD + degradation is partitioned into PARP1-dependent and PARP1-independent reactions; NAD + synthesis is partitioned into AMPK-P-dependent and AMPK-P-independent reactions; AMPK dephosphorylation is partitioned into glucose-dependent and glucose-independent reactions; and so on. In addition, instead of using simple transitions to model phenomenological interactions, 'Delay' variables are introduced to capture the nonlinearity and timescale of more complex chains of interactions. A 'NegReg' dummy variable is used as a limiter of how much NAD + may accumulate in the cell to represent the negative feedback processes promoting increased conversion of NAD + into NADH or other species.
Transitions between light blue blocks show established mechanistic biochemical interactions such as (de)phosphorylation reactions. All other transitions are phenomenological, meaning they represent multiple underlying processes. Phenomenological transitions that were modelled using "Delay" variables are shown in dark blue. Phenomenological transitions modelled without using "Delay" variables are shown in pink, such as the effect of modelled stimuli. Endogenous biological proteins are shown in light blue. Note that 'composite' refers to a model species that is the product of two abundances (SIRT1 and NAD + ) so that the rate of the transition can occur non-linearly through a hill function but still respond to changes in the levels of two species (a Hill function can only respond to one species by design).
Although the model structure is informed byand devised to approximatecurrent biological knowledge, the transition kinetics that determine the dynamics of the model structure are based on parameters whose values are largely unknown. However, these can be quantified through the use of parameter estimation procedures. The calibration of the model with experimental data uses a parameter estimation procedure which explores numerous combinations of potential parameter values to select those which result in the closest fit of the model simulation to the experimental data. To fit all the necessary parameters, we gathered an extensive body of time course data from the literature. Importantly, this experimental data was generated in a single cell line of C2C12 myotubes. See supplementary Figs. S1 to S5 for details of model fitting to the calibration data.
The next task to ensure confidence in the model is to validate model predictions to novel experimental data. We extracted 22 datasets from 15 peer-reviewed publications (see supplementary Figs. S6 to S23). For the majority of network inputs the model predictions for perturbations as in the original experiments results in good qualitative and quantitative accordance with the validation data. We therefore concluded that the model was a good approximation of the underlying biological system. A summary of the depth and breadth of model validation can be seen in supplementary Tables 1. The equations and parameters modelling the AMPK-PGC1α-SIRT1-NAD + regulatory circuit can be found in supplementary Tables 2-5.

Increased PARP1 activity blunts the responsiveness of the AMPK-NAD + -PGC1α-SIRT1 pathway
We used the model to study how the AMPK-NAD + -PGC1α-SIRT1 pathway responds to energetic stress under different conditions of genomic instability. Genomic instability was simulated in the model as a 2.5 fold increase in basal PARP1 activity. A treatment with 0.5 mM AICAR (5-aminoimidazole-4-carboxamide ribonucleotide) was used as a proxy for energetic stress. AICAR is converted to the AMP analogue, ZMP, which is recognised by AMPK as an increase in AMP levels symptomatic of energetic stress. AICAR was introduced into the simulation after the system had equilibrated into a new steady state as a result of the increased genomic instability. As shown in Fig. 3, under conditions of genomic instability, the activation of the AMPK-NAD + -PGC1α-SIRT1 pathway was blunted. The extent of pathway activation as  reflected by the absolute levels reached by both NAD + and deacetylated PGC1α was reduced. As expected, this blunting of pathway activity occurred at the level of NAD + since the upstream activation of AMPK remained unaffected (lines are superimposed in the time course for this species).

Increased NAD + degradation limits the accumulation of the NAD + signal
To explore what feature of NAD + metabolism dampens the activation of the signalling pathway we tested the ability of different treatments to rescue the pathway activation under conditions of genomic instability. Model simulations in Fig. 4 show that the preconditioning of cells with 0.5 mM nicotinamide riboside (NR) for 24 h does not rescue pathway activation by 0.5 mM AICAR as effectively as preconditioning of cells with 1 μM of the PARP1 inhibitor, PJ34.
A possible explanation is that NR supplementation does not reduce NAD + turnover as a result of high PARP1 activity, so much of the additional NAD + is removed, whilst PJ34 works through the inhibition of PARP1, slowing the rate of removal. Model simulations were confirmed by experimental measurements in MEF cells (Fig. 5) where genomic instability was introduced though the knockout of Rev1. The PDK4 data supports model findings. There is reduced activation by AICAR in the Rev1 knockouts and there is rescue, albeit slight, when supplemented with NR. For AICAR concentrations greater than 0.125 mM, there is also evidence of a potentiated AICAR stimulation in the WT cells when NR is present.
Note that this data does not mean that NR supplementation can never be as effective as inhibition of PARP1, but would require a high enough NR supplementation to overwhelm the increased PARP1 activity.
Increased degradation of NAD + is likely to dampen any attempt to stimulate production, even outside genomic instability. Fig. 6 shows a simple system of a signalling trigger with flux (k 3 ) acting on a molecule with synthesis (k 1 ) and degradation/use (k 2 ). A signalling event that introduces a given flux of a signalling molecule will result in a smaller change in the overall levels under higher turnover conditions ( Fig. 6b where blue shows low NAD + levels and yellow high). Importantly, changes in degradation through k2 resulted in pronounced changes in NAD + levels, changes in synthesis through k1 had little to no effect (Fig. 6c). Hence, any age-related mechanism that results in an increase in the rate of NAD + degradation, be it through increased PARP1 or CD38 activity, will have the same dampening effect on the AMPK-NAD + -PGC1α-SIRT1 pathway. This means that both genomic instability and chronic inflammation (through CD38) can interfere with the functionality of this pathway during the ageing process.
This perspective sheds some light on how nuclear maintenance may be prioritised over mitochondrial maintenance. An energetic stress signal activates the AMPK-NAD + -PGC1α-SIRT1 pathway by increasing SIRT1 activity through the increase in NAD + levels but not SIRT1 levels (Canto et al., 2009) -meaning that the reaction is substrate-limitedand so will also increase the substrate for PARP1-mediated reactions. Whereas a genomic stress signal in the form of increased PARP1 activity will reduce available NAD + for SIRT1-mediated reactions. Thus, the prioritisation lies in the asymmetrical influence on the parameters controlling NAD + as shown in Fig. 6a. One signal (energetic stress) will influence NAD + synthesis and the other (genomic stress) will influence NAD + utilisation. Under conditions of conflicting signals where both genomic and energetic stress may be present in an aged cell, genomic maintenance would be prioritised. Note that this mechanism of prioritisation can act in concert with other interactions between PARP1 and SIRT1 (Canto et al., 2015).

Reduced pathway activation promotes the accumulation of dysfunctional mitochondria
Considering that mitophagy can be triggered through the AMPK-NAD + -PGC1α-SIRT1 pathway under conditions of stress (Herzig and Shaw, 2018;Rabinovitch et al., 2017;Zhang et al., 2018), it is feasible that dampening its activity due to a higher rate of NAD + utilisation would lead to a reduced ability to maintain healthy mitochondrial populations and the potential accumulation of dysfunctional mitochondria.
This can be shown in principle by coupling this validated model of the AMPK-NAD + -PGC1α-SIRT1 signalling axis to a simple model of mitochondrial populations (Fig. 7). In this coupled model, newly formed (healthy) mitochondria can become damaged/old at a given rate, triggering AMPK activation through an increased AMP/ATP ratio in order to enhance mitophagy and remove the old mitochondria. Upon pathway activation, deacetylated PGC1α induces mitophagy and thus the disappearance of both types of mitochondrial populations at different rates (in the scale of those reported by (Dalle Pezze et al., 2014)). To model   curve 1), but increased dysfunctional mitochondria accumulation under genomic stress (2.5 a.u PARP1, comparing curve 3 to curve 1).
One explanation is that high PARP1 activity lowers basal deacetylated PGC1α levels below the value of 1 a.u as assumed in the simulations involving the mitochondrial module alone. Such simulations would suggest that increased NAD + degradation with age, caused by genomic instability, chronic inflammation or both, would interfere with the cell's ability to regulate mitochondria. This is especially true since the AMPK-NAD + -PGC1α-SIRT1 pathway also controls FOXO-induced mitophagy through the same mechanism (Canto et al., 2009) and so the pathway dampening would affect mitophagy induction through two main regulatory transcription factors: FOXO and PGC1α.
To make the simulations more realistic, the positive feedback loop driving a given rate of dysfunctional mitochondria accumulation over 80 years can be coupled to PARP1 so that the levels of this molecule gradually increase over the years instead of remaining fixed at a given value. Modelling such age-related changes to PARP1 levels still shows the same effect as that shown in Fig. 8 (supplementary Fig. 25). Perhaps unsurprisingly, simulating an age-related decrease in SIRT1 levels (supplementary Fig. 26) and AMPK levels (supplementary Fig. 27) in a similar manner also results in an increased accumulation of dysfunctional mitochondria due to the reduced ability of the latter to activate the AMPK-NAD + -PGC1α-SIRT1 signalling pathway to induce mitophagy. Interestingly, an age-related increase in AMPK dephosphorylation simulated by an age-related hyperglycaemia does not result in any changes to the accumulation of dysfunctional mitochondria since the constant stress signal causing an increased AMPK phosphorylation would balance the increased dephosphorylation (supplementary Fig. 28).
In accordance with the established beneficial effects of life-long NAD + supplementation, the simulation of this intervention by fixing NR to a value of 100 a.u (to model a life-long supplementation with 0.1 mM Nicotinamide Riboside at the cellular level) leads to a marked reduction in the accumulation of dysfunctional mitochondria over 80 years (supplementary Fig. 29). However, our model suggests that under the parameters we have calibrated to experimental data and then validated against additional data, the molecular network regulating NAD + responds more efficiently to interventions that reduce genomic stress than those supplementing the NAD + pool. Genomic stress could induce mitochondrial dysfunction by reducing PGC1α activity through NAD + depletion that would be difficult to counteract by NAD + supplementation.

Discussion
NAD + is a molecule involved in a myriad of biological processes (Fakouri et al., 2019). Its relevance to the ageing process is highlighted by a multitude of studies that report the modulation of lifespan and healthspan through changes in basal NAD + levels (Balan et al., 2008;Belenky et al., 2007;Cerutti et al., 2014;Guan et al., 2017;Mouchiroud et al., 2013). The lowering of basal NAD + levels with age has been associated with mitochondrial dysfunction (Camacho-Pereira et al., 2016), and here we show here that high NAD + turnover can prevent adaptation by increasing mitophagy. Interestingly, Dalle Pezze et al. (2014) reported a similar effect where pharmacological interventions aimed at reversing fibroblast senescence were less effective due to a higher mitochondrial turnover. This is not the only mechanism that can render a regulatory system less sensitive to signalling. We also demonstrate that reported changes to the AMPK-NAD + -PGC1α-SIRT1 pathway such as decreased AMPK levels and SIRT1 levels would result in a loss of pathway sensitivity. The computational model did not imply this for an age-related increase in AMPK dephosphorylation as the stress signal increases phosphorylation. However, this assumes that AMPK's ability to sense AMP/ADP and be subsequently phosphorylated remains unchanged with age, which may not always hold true (Hardman et al., 2014).

Conclusion
In this work we argue that the increased NAD + degradation reported with age is likely to lead to a dampening of the activation of the AMPK-NAD + -PGC1α-SIRT1 signalling pathway. As ageing involves increasing levels of genomic instability and activation of PARP1, supplementation of NAD+ or NR would require high doses to overcome the increased rate of degradation to reactivate SIRT1 and AMPK. Consequently, mitophagy and other downstream processes will decrease with age even as stress levels and damage to mitochondria rise, thus leading to the accumulation of dysfunctional mitochondria.

Computational
Model schematics were created in Cell Designer (version 4.4.2) consistent with SBGN PD standards.

Model simulation
A computational model of the AMPK-PGC1α-NAD + -SIRT1 axis was developed using the python programming language and the pycotools3 python package (Welsh et al., 2018). Pycotools3 (version 2.1.22) is a python package that operates as a front end to the COPASI program (Hoops et al., 2006) allowing users to call COPASI from python code to simulate or run parameter estimations on a given set of models and data. In particular, it allows models to be defined using antimony strings, a language used for defining systems biology models using the tellurium library (Medley et al., 2018).
Tellurium (version 2.2) is a systems biology python package for simulating biological systems that defines its own easy to read and edit, 'code like,' model definition called antimony. Antimony models can be simulated in Tellurium directly or converted to SBML for use with other software.
COPASI (version 4.27.217) is a suite of systems biology algorithms including algorithms for simulation, steady state finding and parameter estimation. It can be used via the command line as well as via a graphical user interface (GUI).
This model is based on coupled ordinary differential equations (ODEs) that simulate changes in species' abundances over time. The model is simulated deterministically using the LSODA algorithm with values of 'Relative Tolerance' and 'Absolute Tolerance' of 1e − 6 and 1e − 12 respectively for a 'Max internal steps' of 10,000.

Simulation of AICAR treatment
The simulation of a 0.5 mM AICAR treatment involves setting the initial abundance 'AICAR_treatment' variable in the model from 0 to 1. A unit increase in the value of this variable corresponds to a 0.5 mM increase in the dose of the AICAR treatment being simulated. Thus, the simulation of a 2 mM AICAR treatment corresponds to a setting of this variable to a value of 4.

Simulation of NR supplementation
The simulation of NR supplementation in the model involves setting the initial abundance of the NR-NMN variable from 0 to X where X is the NR concentration in μM that aims to be simulated. Hence, the simulation of a 500 μM NR supplementation involves setting this variable to a value of 500.

Simulation of glucose restriction
The simulation of glucose restriction involves changing the value of the Glucose_source to the desired value of a glucose steady state in mM. For example, the simulation of a glucose restriction protocol from 25 mM to 5 mM glucose involves changing the value of Glucose_source from 25 to 5 (or vice-versa for the simulation of a glucose addition protocol). Note that the default value of the latter parameter is 25 since the standard cell culture protocol of C2C12 myotubes, from which most data is utilised in the development of the model, involves 25 mM glucose.

Simulation of a PJ34 treatment
The simulation of PARP inhibition through PJ34 is modelled by setting the initial abundance of the 'PARP' variable from 1 to 0.

Simulation of genomic instability
The simulation of genomic instability feeding through increased PARP1 activity is modelled by changing the initial abundance of the 'PARP' level to the fold-change measured experimentally. This means that if an experimental model of genomic instability displays a 3-fold increase in PARP1 activity over controls, the initial abundance of the 'PARP' variable should be changed from 1 to 3 in the model.

Simulation of SIRT1 knockout
The simulation of SIRT1 knockout involves reducing the initial abundance of the 'SIRT1' variable. The value of the 'SIRT1' variable in the model should be altered from 1 to 0.

Parameter estimation
The parameter estimation procedure was carried out using multiple rounds where parameter values shown to be well behaved as indicated by profile likelihood's (Welsh et al., 2018), were fixed for the next round. The parameter estimation involved the use of the 'Particle Swarm' algorithm with an 'Iteration Limit' setting of 7000 and a 'Swarm size' of 350 (with 'Std. Deviation' left at standard setting of 1e − 6 ). Full source code used for parameter estimation and creation of figures can be found at (https://github.com/npc101-ncl/NAD-SIRT1-model).

Data storage
SBML versions of both models are stored in BioModels under the same ID: MODEL2206300001.

Cell culture preparation
WT and Rev1− /− mice embryonic fibroblasts (MEF) cells were provided by Dr. Niels de Wind from the department of human genetics at Leiden University Medical Center (LUMC). MEF cells were cultured in Dulbecco's Modified Eagle Medium (DMEM) plus 10 % (v/v) fetal calf serum (Invitrogen) and 100 U/ml penicillin-streptomycin (Gibco, Life Technologies) in a humidified atmosphere of 5 % CO2 at 37 • C. All the experiments are carried out two weeks after the initiation of the cell culture or as stated. One day before treatment, cells were plated in 6 well plates (Costar, Corning) at density of 4 × 10^4 cells/ml (determined by trypan blue staining (1:1 ratio) after trypsinization and counting under haemocytometer). The next day, cells are treated with and without 3 mM NR (NIAGEN, ChromaDex, USA) for 18 h and treatment was terminated by washing with 1XPBS (D8537, Sigma-Aldrich) twice. Then cells are again treated with different concentrations of AICAR (Toronto Research Chemicals, Canada, 0.125-1 mM) for 12 h and cells were washed with 1XPBS and collected in lysis buffer, and RNA was extracted as per protocol with Qiagen RNeasy kit (ID:74004)..

Real Time-PCR
Total RNA from Rev1+/+ and Rev1− /− MEF cells was extracted and RNA concentration and purity were determined using an (λ = 260/280 nm; Epoch Spectrophotometer, Bio-tek). The first strand of cDNA was synthesised from total RNA using the high-capacity cDNA reverse transcription kit (Applied Biosystems, USA). Primer sequences were designed using the NCBI primer design tool and all primers were purchased from oligos, Denmark. Sequence: F GAGCTGTTCTCCCGCTACAG R CGGTCAGGCAGGATGTCAAT A 20 μl reaction mixture containing 10 μl SYBR-Green Master mix (Amplicon), 100 nM PDK4 gene-specific primers and 1000 ng/μl of cDNA template was loaded into each well of a 96-well plate. All PCR reactions were performed in triplicate. PCR was performed under these thermal conditions: 95 • C for 12 min, followed by 40 cycles of 95 • C for 15 s, 60 • C for 20 s, and 72 • C for 20s. A dissociation melt curve analysis was performed to verify the specificity of the PCR products. Actin primers were used to amplify the endogenous control product. The relative mRNA expression analysed using the 2− ΔΔCT method.