A multimethod computational simulation approach for investigating mitochondrial dynamics and dysfunction in degenerative aging

Summary Research in biogerontology has largely focused on the complex relationship between mitochondrial dysfunction and biological aging. In particular, the mitochondrial free radical theory of aging (MFRTA) has been well accepted. However, this theory has been challenged by recent studies showing minimal increases in reactive oxygen species (ROS) as not entirely deleterious in nature, and even beneficial under the appropriate cellular circumstances. To assess these significant and nonintuitive observations in the context of a functional system, we have taken an in silico approach to expand the focus of the MFRTA by including other key mitochondrial stress response pathways, as they have been observed in the nematode Caenorhabditis elegans. These include the mitochondrial unfolded protein response (UPR mt), mitochondrial biogenesis and autophagy dynamics, the relevant DAF‐16 and SKN‐1 axes, and NAD +‐dependent deacetylase activities. To integrate these pathways, we have developed a multilevel hybrid‐modeling paradigm, containing agent‐based elements among stochastic system‐dynamics environments of logically derived ordinary differential equations, to simulate aging mitochondrial phenotypes within a population of energetically demanding cells. The simulation experiments resulted in accurate predictions of physiological parameters over time that accompany normal aging, such as the declines in both NAD + and ATP and an increase in ROS. Additionally, the in silico system was virtually perturbed using a variety of pharmacological (e.g., rapamycin, pterostilbene, paraquat) and genetic (e.g., skn‐1, daf‐16, sod‐2) schemes to quantitate the temporal alterations of specific mechanistic targets, supporting insights into molecular determinants of aging as well as cytoprotective agents that may improve neurological or muscular healthspan.

. Cellular and anatomical parameters of interest for 3-day-old C. elegans (N2); ranged/averaged literature values.
Individually optimized value, generated using AnyLogic calibration experiments to best fit or emulate reported kinetic curves.
Quantitated assumption based on qualitative evidence. Units converted to adhere to model standards.
Box S2: Equations and event statements pertaining to mitochondrion-level biochemical kinetics.
Differential equations governing the time-dependent amounts of key players (PQROS represents a function of paraquat-induced increases in ROS production rate, and PT represents a function of pterostilbene-mediated ROS mitigation, as outlined in the last section; when unperturbed by either compound, functions PT and PQROS = 1): *Rate-dependent discrete events for changes in mtDNA copies: Damage to mitochondrial DNA is assessed here in totality. All perturbations (mainly point mutations and deletions) are regarded as a change in mitochondria DNA that will lead to dysfunctional protein synthesis and deleterious faulty protein aggregation. Copy numbers vary from mitochondrion to mitochondrion (2-10 copies), but it is assumed here that the mean copy number is 6, and the distribution falls rapidly in either direction. It is also assumed here that all mitochondria begin with ~15% perturbed mitochondrial DNA (~1 copy per mitochondrion), and that all subsequent DNA damage is both rate-dependent and ROS-dependent [16]. The damage rate is also not always the same for every mitochondrion, and therefore is represented by a stochastic ROS-dependent rate equation (Eq. 1) with a mean number of events per week and a corresponding exponential distribution for the initial trigger and recurrence time. Additionally, once mtDNA has been altered, it is apparent that the cell works to remove it when unstressed [3]; however, following the stressor-mediated nuclear activation of ATFS-1, DmtDNA is more easily sustained [17] (Eq. 2).

*Dynamic expressions displayed in blue:
Concentrations of proteases, chaperone proteins, ETC assembly factors, and superoxide dismutase are determined by total cellular concentrations following nuclear synthesis (as outlined in the cell-level systemdynamics methods). Additionally, the inhibitory factor governing the repression of mtDNA-mediated OXPHOS protein production is determined by the total mitochondrial concentration of ATFS-1 (all mitochondria presumably contain the same concentration at any given time due to the nature of the global response over a cell).
Box S2 (continued): Equations and event statements pertaining to mitochondrion-level biochemical kinetics. Oxygen consumption (JO), as dictated by both OXPHOS activity and ROS production (pmol/min/worm): Sirtuin activity, as dictated by both NAD + breakdown and stilbenoid polyphenol exposure:

Mitochondrion-level agent-based computational methods
Box S3: Conditional and rate-dependent statements governing the rule-based categorization of mitochondrial stress, via the agent-based paradigm.
It is assumed that all mitochondria in the early stages of life are initially fully functioning and unstressed [18,19]. As life progresses, the mitochondria undergo various levels of stress [19,20], which have been roughly phenotyped and categorized for the purposes of the agent-based model. In our model, a mitochondrion can be (A.) functioning and optimal capacity (maximal ATP production), with little to no available stressors (ROS, accumulated proteins, damaged DNA), (B.) functioning below maximum capacity but still within the optimal range, in the presence of increased stressor concentrations, (C.) significantly stressed and malfunctioning with respect to energy production, or (D.) severely stressed with high amount of damaged genetic material, thus producing very little ATP and large amounts of superoxide. The transitions between these states are rate-dependent and governed by the corresponding stressors that account for these pathological micro-phenotypes [19,20], as seen in the equations below. The constants involved in facilitating these stress and recovery transitions have been optimized (StressConstant = 0.850 events %EPfaulty -1 nMROS -1 hr -1 ; RecoveryConstant = 1.775 events %EPfaulty nMROS hr -1 ) so that the stress-state results relate closely to true biological stressor levels in WT C. elegans [3,17]. Relative to these stress-states, the model accounts for mitochondrial autophagy (mitophagy), where the cell will recognize the existence of damaged mitochondria and effectively destroy such organelles in a rate-dependent manner following retrieval of such stress signals [21,22]. The stress signals emitted from the mitochondria vary in magnitude from state to state, as defined by the model-where insignificantly-stressed mitochondria will weakly activate the mitophagic mechanisms, significantly stressed mitochondria will most effectively activate mitophagy, and severely stressed mitochondria will suffer from hypothesized UPR mt -mediated overprotection and compromised stress signaling. Mitophagy here is assumed to be initiated by ineffective mitochondrial protein transport as well as SKN-1 or DAF-16 activation [22,23,34,55,56], and is therefore inversely correlated with these corresponding activities (outlined in the cell-level system-dynamics methods). Mitochondrial recycling/biogenesis is also a concomitant event that occurs within mitochondrial stress and quality control in the aging C. elegans [22], and it is here accounted for by (1.) the movement from high stress states to low stress states, and (2.) a timeout-triggered (every 1000-3000 min) agent-addition event. Such agent addition events occur more rapidly upon SKN-1 activation [22,55], where biogenesis is increased by 7-20% at maximum activity.
Stochastic rates for movement of increasing stress states, based on protein accumulation: Stochastic rates for movement of decreasing stress states, based on recovery of mitochondrial quality: Transition 5 : Mito-transport-, SKN-1-, and DAF-16-dependent [34] rates of mitochondrial loss via mitophagy (Segments "R" and "B" represent xenobiotic-induced perturbations of mitophagy rate, outlined in the last section. When unperturbed, R = 1 and B = 1):     *DAF16_activity is coupled with MnSOD production as a known antioxidant response element for this transcription factor. All mitochondrial proteins are also reliant on the efficacy of their transport systems, as delineated by the mathematical expression that utilizes the translocase abundance.

Cell-level agent-based computational methods
Box S5: Conditional and rate-dependent statements governing the rule-based categorization of cell integrity, dependent upon the mitochondrial agent population within each cell agent.
Respiratory capacity is reportedly lost when deleterious mtDNA mutations hit a threshold value between 60% and 90% [17,21,29]. It is assumed that this threshold value can be translated directly as general loss of mitochondrial function, which have been outlined by the conditional eustress and distress states of the model. This distribution of threshold values (triangular distribution used in the model: 60%-62.5%-75%) is used to govern the loss of cellular integrity and the increased probability in cell death. Recovery of such cells in the model is possible, if their mitochondria can restore enough function to return below the lowest possible threshold value of 60%. This strict assumption is made to uphold the standards of true biological systems, which face a great deal of difficulty restoring homeostasis when mitochondrial populations are largely compromised [30]. It has been reported that humans begin to lose brain volume around the age of 40, at a rate of 5% every decade, mainly due to spontaneous neuronal degeneration [31]. If we assume the average human life expectancy at birth is 70 years and if we assume these degradation checkpoints are similarly conserved among eukaryotes like C. elegans (average N2 lifespan @ 25 O C: 20 days [32]), then neuronal death will begin to occur within the nematode around day 11.5 at a rate of 5% every 2.8 days, with a more rapid decline immediately before death. That being said, it has been reported that C. elegans do not lose substantial neuronal integrity [57], yet vacuolated compromised neuronal cells and neurodegenerative responses have been observed in these organisms [58]. In light of these findings, the model was designed to primarily determine the dynamics of the age-compromised state for energetically-demanding cells, while attempting to shed light on potential tissue degeneration thereafter. Once mitochondrial conditions force the cell into a compromised state, they have a limited time to recover or they are destroyed. This rate is defined in the model by an exponential distribution that has been characterized meet expected physiological neuronal lifespan.

Use of Rapamycin and Bafilomycin
One virtual pharmacological perturbation mode within the model addresses mitophagy flux. One perturbation deals with the established role of the mechanistic target of rapamycin (mTOR), which normally acts to inhibit selective mitophagy and reduce the rate of mitochondrial degradation. It has been identified that this protein is well conserved among eukaryotes, and such a homolog of this protein has been identified in C. elegans [33] which is indeed, as it is in mammals, inhibited upon rapamycin administration [34]. It has been well established that rapamycin acts through this mechanism to decrease heteroplasmic mtDNA accumulation in immortalized neurons [35] and to alleviate general mitochondrial dysfunction in whole rat brains (as assessed by oxidative damage biomarkers and ATP levels) [36]. Assuming this therapeutic effect is mechanistically conserved among eukaryotes, rapamycin is used here as a virtual xenobiotic to increase the rate of mitophagy. It is administered virtually within different dosing schemes at a concentration of 15 nM (which equates to 18.3 µg/L, within the therapeutic range for serum levels of 5-20 µg/L [35]). Once administered instantaneously, it virtually affects the cells by a maximal fraction (X R ) of 1.5 in relation to its EC 50 (estimated at ~10 nM [35]), and it begins to decrease in concentration with an aqueous half-life of about 10 h [37].
Box S6: Equations used to capture tissue concentrations of parameters of interest, where "i" is the current number of total mitochondria in the system and "n" is the current number of neurons. The system initially contains 65 neurons for simplicity. Bafilomycin A1, in contrast, facilitates an inhibition of autophagy in C. elegans by preventing the maturation and accumulation of autophagic vesicles [38]. Bafilomycin A1 is virtually administered here within different dosing schemes in a similar fashion at 10 nM, with a maximal fraction (X B ) of 0.5, an EC 50 of ~7 nM [38], and a low-end estimate of the culture half-life at 40 h (adapted from in vivo pharmacokinetic data for the structurally similar azithromycin [39]). The pharmacodynamic equations outlined below (Box S7) for rapamycin or bafilomycin administration are multiplied by the stochastic mitophagy transitions 7, 8 and 9 within the mitochondrial agent-based paradigm (Box S3), and subsequently integrated over time for different dosing schemes (note: there is no direct effect when concentrations of either chemical are zero).

Use of Paraquat
As an additional and different type of mechanistic perturbation, paraquat is virtually used here for its aging-related pro-oxidant properties within the mitochondria [40][41][42] and for its ability to induce the UPR mt [18]. Paraquat is virtually administered here in a similar fashion to rapamycin and bafilomycin, at a concentration of 100 µM (high) or 5 µM (low), with an environmental half-life that far exceeds the normal nematode lifespan [43], and so paraquat concentrations will simply be clamped within different dosing schemes. Paraquat perturbs the system using the same types of equations, however, since paraquat perturbs the system in two different areas (superoxide production rate and mitochondrial import of ATFS-1) in two distinctly different and opposing ways, the functions for each are separately categorized, as outlined below (Box S7). For the function pertaining to its pro-oxidant properties: the constant for maximal rate increase, X PQR , is set at 0.4, and the EC 50 for this effect is set at 25 µM [estimations; 44,45]. For the function pertaining to its UPR mt activation (presumably due to lack of ATFS-1 mitochondrial import efficiency): the constant for maximal impedance effect, X PQU , is set at 0.75, and the EC 50 for this effect is set at 50 µM [estimations ; 18]. These functions PQ ROS and PQ UPRmt (Box S7) are integrated over time for various dosing schemes and are then multiplied to the rates of mitochondrial superoxide production (Box S2) and mitochondrial ATFS-1 import (Box S4), respectively (note: there is no direct effect when exposure concentration of paraquat is zero).

Use of Pterostilbene
Plant-derived stilbenes have been extensively studied as anti-aging therapies, with special emphasis on resveratrol and pterostilbene (prominent constituents of blueberries, red grapes and red wine products). Pterostilbene has been recorded as a more potent neuromodulator than resveratrol in aging, even when administered at low doses [46], and therefore is used as the stilbene in this simulation study. Pterostilbene has direct radical scavenging properties and mitigates ROS through highly conserved oxidative stress response pathways as well, via sirtuin activity that acts upon transcription factors such as SKN-1 and DAF-16 in worms [47,52]. These effects ultimately manifest as reduced ROS content and increased bioenergetic efficiency in aging organisms, and specifically has reduced the levels of ROS in C. elegans by about 26% when constantly administered at a concentration of 100 µM [47]. Because of these observations, the pterostilbene function used in this model (same format as rapamycin) utilizes an EC 50 of 20-60 µM (conservative estimate) and a maximal ROS-mitigating (X PTR) and sirtuin-activating effect (X PTS ) of 0.26 and 0.5, respectively, which are applied to the non-peroxide-forming superoxide elimination rate and the Sirtuin_activity dynamic variable. Pterostilbene is administered here at a dose of 100 µM and assumes a culture half-life of 12 h due to its poor aqueous stability [48] and rapid pharmacokinetic clearance [49]. All of these factors are incorporated into the function PT (Box S7), which is applied to the rates (seen in Box S2) of mitochondrial superoxide elimination and peroxide elimination (note: there is no direct effect when exposure concentration of pterostilbene is zero).

Varied Dosing Schemes
Dosing schemes for rapamycin, bafilomycin, pterostilbene and paraquat were altered and/or combined to produce an array of exposure scenarios. For rapamycin, bafilomycin and pterostilbene, doses were administered at 0 min and every 2500 min thereafter; the resulting cellular exposure at any given time point was determined by their culture half-life kinetics. For paraquat, high or low doses were clamped to maintain a constant exposure value. For all of these xenobiotics, their administration was performed constantly, early (0-10000 min), mid-life (10000-15000 min), or late (15000-30000 min).

Genetic Perturbations
To identify the role of key genes as they relate to this system, several were virtually knocked down and the endpoints of interest were observed and compared to the control group. Such genes were those encoding for MnSOD, SKN-1, and DAF-16. To make these virtual perturbations, the rate of productive expression for these functional proteins were fractionalized for the duration of the simulation, to simulate such reduced genetic accessibility for the corresponding genes.