Cell Level- Modeling of Aging and Rejuvenation

Andrey Koptyug ^{1, †, ‡, *}, Yurij Sukhovei ^2, †, Elena Kostolomova ³, Irina Unger ², Vladimir A. Kozlov ^{4, †}

1. Mid Sweden University, Akademigatan 1, Östersund, Sweden

2. Institute of Fundamental and Clinical Immunology, Tyumen branch, Kotovskogo str., 5, Tyumen, Russia

3. Tyumen State Medical University, Kotovskogo str., 5/2, Tyumen, Russia

4. Institute of Fundamental and Clinical Immunology, Yadrintcevskaya str., 14, Novosibirsk, Russia

‡ Current Affiliation: Associate Professor, SportsTech Research Center, Mid Sweden University, Östersund, Sweden

† These authors contributed equally to this work.

* Correspondence: Andrey Koptyug

Academic Editor: Raya Khanin

Special Issue: Utilizing Big Data to Elucidate Skin Aging [Big Data in Skin Aging]

Received: April 14, 2023 | Accepted: December 13, 2023 | Published: December 22, 2023

OBM Geriatrics 2023, Volume 7, Issue 4, doi:10.21926/obm.geriatr.2304263

Recommended citation: Koptyug A, Sukhovei Y, Kostolomova E, Unger I, Kozlov VA. Cell Level- Modeling of Aging and Rejuvenation. OBM Geriatrics 2023; 7(4): 263; doi:10.21926/obm.geriatr.2304263.

© 2023 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

Population aging is one of the achievements and, at the same time, burdens of modern society. According to the WHO, the proportion of people aged 60 years and above is steadily growing. Their number currently exceeds 1 billion and is expected to reach roughly 1.4 billion by 2030 [1]. According to a report by the UN, between the years 2019 and 2050, the share of people over 65 will grow by about 50% in the developed world and by up to 220% in the developing world [2]. The consequences of these changes will be both economic and social. In 2019, Japan was the world leader in the old age dependency ratio - calculated as the number of people older than 65 years of age divided by the number of people aged 15-64 - closely followed by the EU, USA, and China [3]. Research shows that this proportion is growing significantly for all countries in both the developed and developing world. In economic terms, a decreasing active labor force and a corresponding decrease in labor productivity are alarming and require policymakers to act (e.g. [4,5,6,7,8]). Obvious social issues related to population aging – such as an increasing burden on healthcare systems - are also growing. Such social issues are complex to separate from purely economic matters. Thus, many are propagating the necessity of “active and healthy aging” (e.g. [9,10,11]).

Skin aging is one of the factors strongly affecting self-esteem. It can often become a severe problem due to high societal pressure [12,13,14,15,16]. Decreased self-esteem and even potential depressive states, in many cases, lead to severe issues with overall health. Simultaneously, skin health (especially face and arms) is a good indicator of overall health, reflecting the aging process at other body system levels [17,18,19,20,21,22]. Many interventions aimed at supporting and improving skin health are connected to improvements in the immune status and general health. Thus, “measuring” aging is one of the critical aspects in both the assessment of the skin’s state and the dynamics of its changes (aging) and in the evaluation of the results of possible interventions (for example, aimed at rejuvenation). Different methods have been developed for quantitatively assessing skin aging and providing measurable parameters for objective comparison. These include image analysis, skin roughness and elasticity measurements, sebum and moisture content measurements, and cell profile studies, to name a few. Many cumulative predictors - for example, metabolic and epigenetic - are calculated from several different measurable parameters [23,24,25]. So-called biological age is also commonly used in such assessments [26,27,28,29,30,31,32,33,34,35]. Significant effort is also being put into developing anti-aging and rejuvenation strategies and methodologies [34,35,36,37,38,39,40,41,42,43]. Many of these are based upon top-level approaches developed for the body-organ-tissue levels (skin in particular). Many others are based on bottom-level methods (e.g., cell and sub-cell levels). Our research indicates that further progress in developing anti-aging and rejuvenating strategies and interventions can enormously benefit from cell-level modeling and data-centered approaches.

This project, related to the cell level-based approach and including the development of conceptual and computer models and experimental research, was initiated to explain several unique properties of natural substances believed to have anti-aging and rejuvenating properties (e.g. [33,34,35]). Empirical research on skin aging, skin rejuvenation, and skin aging quantification led to the conclusion that more attention should be given to cell-level studies simultaneously combining approaches common to cell biology, biophysics, and general formalisms from mathematics and modeling [21,22]. The present paper addresses issues related to defining aging and rejuvenation for cell populations. This is related to modeling that involves a vast number of cells as representative elements. Several key concepts - such as that of proliferation niche and homeostatic cell number stabilization in the cell population through coupled action of proliferation and apoptosis and regulation based on the maintenance of concentration of a single regulatory substance are introduced and described. We also present a computer simulation approach developed based on the conceptual modeling approach. Our initial conclusions are based on the results of modeling, simulations, and coupled experiments. We suggest possible ways forward, including application of big data methodologies and algorithms in modeling aging and rejuvenation.

2. Assessment of Cell Population Aging

2.1 Approaches to Quantifying Aging and Rejuvenation

One of the common approaches is studying aging and rejuvenation at the body-organ-tissue levels. Many issues related to learning and modeling top-level systems are their complexity, individual differences between humans, dynamics of the involved processes, and problems with registration and quantification of relatively slow changes. Most of the studies target changes in the facial skin and link them to the aging of the whole system [21,22,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49]. Corresponding studies often involve large groups of subjects, unique methods of data analysis aimed at minimizing influence from individual differences, and lengthy, sometimes life-long trials (e.g. [45,46,47,50]). Deriving conclusions from such data often depends on statistics, machine learning, and big data approaches [51,52,53,54,55,56,57].

A typical top-level approach to age quantification uses chronological age (calculated from birth). However, it does not adequately represent the actual aging of a particular subject. So-called biological age, a cumulative parameter calculated based on the quantified results from several objective tests, more adequately describes the natural aging process of a particular person [21,22,23,24,25,26,27,28,29,30,31,32]. Researchers use different ways to calculate biological age and, sometimes, based on entirely different input parameters. This brings additional uncertainty and complicates the comparison of results from other studies. Nevertheless, comparing chronological and biological age values could be helpful. For example, experimental studies indicate that a difference between chronological and biological age could be successfully used for the assessment of the aging of human skin and could provide needed support for individualized decisions on possible treatment [21,22].

Cell-based studies provide certain advantages as they do not need complex ethical approvals, are easier to quantify, can be standardized, and can be replicated for objective validation. However, no reliable link allows a transfer of results from cell culture studies to tissue-organ-body levels. An additional issue is objectively defining age at the level of cell populations with large numbers of cells. Indeed, the age of individual cells could be clearly defined as the time passed from the moment the cell is born in the cell division process. The situation becomes complicated in the case of cell populations containing cells with newly born, young, aging, and old cells.

One can benefit from the analogies taken from modern mathematical physics in linking aging and rejuvenation at the cell level to that at the tissue-organ-body level. Indeed, the life of cells could be regarded as quantized. The cell division process starts from one cell and ends with two - and only two cells - if successful and only one if the operation fails. Apoptosis (non-pathological removal of old, senescent, and damaged cells) also deals with whole cells. Processes at the cell level are regulated by biochemical signaling. This means that when a single signaling molecule reaches a specific receptor, one act in a complicated chain of processes is initiated. In both cases, only integer numbers of cells and molecules are involved. From this point of view, basic process steps are quantized. Using the accepted estimate of 10¹⁴ cells in the human body, a representative tissue sample should contain about 10⁹ to 10¹⁰ cells/cm³. Similar numbers of cells are characteristic of reasonably sized cell cultures. For relatively large regulatory molecules, for example, with an atomic weight of 10⁴, one milligram of a substance contains 10¹⁶ individual molecules. Thus, models at the single cell level and a single dynamic process step involving cells should be expanded to models operating with significant numbers of acting elements - in the range of 10¹⁰ to 10¹⁶.

2.2 Age of a Cell and Age of a Cell Population

Clearly, the ages of all individual cells adequately represent the age of a population with large numbers of cells. At the most superficial model level, one can assume that we are dealing with cells of the exact origin and with the same functionality. We can also assume that newly divided cells have an age equal to zero, and all cells increase their age as time passes. In this case, we deal with a situation when elements of a chosen set differ only by one parameter (age), dynamically changing with passing time. Thus, we can adequately represent the cell population by its age profile, e.g., the distribution of the cells’ age (assuming that cells with the same age are indistinguishable). If one could experimentally acquire such a distribution or its parameters, it would allow for comparing changes happening to the cell population regarding aging and rejuvenation. The population is aging if the relative amount (proportion) of old and mature cells increases. The population becomes younger if the relative numbers of newly born and immature cells increases. Thus, the proposed approach should provide a method of objective comparison between two different cell populations or of the changes in the age profile of the same cell population with time or before and after specific interventions.

Figure 1 presents example cell age distributions reflecting changes: (a) when younger cells are added through intensified proliferation and older cells are simultaneously removed by apoptosis and (b) when individual cells are forced to live longer. The shown normalized distributions represent the relative numbers of cells with a certain age in a chosen population. The changes caused by adding newly born cells and removing older and senescent ones may be of internal and external origin to the selected cell population. For the most basic model, it is natural to concentrate attention on the internal mechanisms. In such cases, we disregard possible cell addition from outside or forced senescent cell removal and focus on cell division and apoptosis.

Figure 1 Exemplified cell age distributions reflecting changes when younger cells are ‘added’ through intensified proliferation, and simultaneously older cells are ’removed’ by apoptosis (a), and when individual cells are forced to live longer (b)—initial distribution- solid line; distributions after corresponding changes- dotted lines. The bottom bars represent the changing boundaries of the first and fourth quartiles.

Traditional low-level models related to cell division often use the approximation of the so-called stem cell niche dealing with stem cells most active in this process and their immediate microenvironment. Knowing that more adequate models should involve large numbers of cells capable of proliferation, we introduce the concept of proliferation niche. We define it as a space-restricted population of similar cells where two mechanisms, namely cell division and apoptosis, are simultaneously active [33,34,35]. In the case of adult tissue, one should also acknowledge continuous cell number maintenance (homeostasis), and cell number in the selected volume (cell concentration) is stable over time. To keep this balance in live tissues with the maintenance of cell numbers, two mechanisms should act simultaneously and synchronously: new cells are provided through proliferation; senescent, old, and damaged cells are removed through apoptosis [33,34,35]. As a result, the number of cells in a mature proliferation niche is stabilized (Figure 2). This statement is supported by experimental data, including a common observation of cuts and lesions healing when newly formed scar tissue only slightly exceeds the volume of the undamaged skin.

Figure 2 Exemplified homeostasis (cell number maintenance) in the proliferation niche, supported by the synchronized action of proliferation and apoptosis.

It should be noted that the boundaries of such proliferation niches could be somewhat arbitrary and are chosen only for modeling purposes. We conclude that the microenvironment of a single cell is inadequate for our current modeling purpose, and some cells with similar properties should be involved. Thus, the proliferation niche is defined as a chosen population of cells with similar functionality, spatially restricted, and with active proliferation and apoptosis. The situation is similar to defining the boundaries of the primary colors within an artist’s palette, where paints are added and mixed, generating multiple diffuse halftones.

Analyzing the homeostasis in a proliferation niche applied to the cell population age immediately leads to several possible conclusions. Intensification of the cell proliferation and apoptosis (separately or jointly) in a cell population with cell number maintenance leads to rejuvenation of the system (more significant percentage of younger cells and smaller of older cells, Figure 1a). At the same time, enabling individual cells to live longer means that apoptosis should fall in intensity. In cell populations with cell number maintenance, this inevitably leads to the reduction of proliferation intensity. This leads to a decreasing proportion of younger cells and an increasing proportion of older ones. A corresponding shift in the first and fourth quartile boundaries in Figure 1b clearly illustrates this situation. It is further concluded that independent of the intrinsic or extrinsic nature of young cell addition and senescent and old cell removal, such intervention should lead to the rejuvenation of the cell population. Experimental observations on the effect of senescent cell removal support this suggestion [58,59,60,61,62,63,64].

This discussion means that the age of a cell population in vitro or selected tissue (in vivo or ex vivo) could be adequately represented by corresponding cell age distributions. Moreover, the effect of the interventions claiming rejuvenation and anti-aging products could be experimentally tested and analyzed based on the dynamics of cell age profiles. The proposed approach is cell-level-centered, and many tests could be performed in vitro or using ex vivo cell cultures. This would help avoid multiple risks and decrease the costs of experimental trials (e.g., [33]). This leaves a question about practical methods that could yield such distributions or other adequate quantitative measures related to cell age profiles.

2.3 Assessment of the Cell Population Age Profile Using Flow Cytometry

One of the methods that can be applied for monitoring cell population aging and rejuvenation is flow cytometry. It is a technology measuring and analyzing different cell characteristics as they flow through a collimated beam of light (e.g. [65,66]). The method allows quantifying various aspects of the cell life cycle using labeling with dyes or fluorescent markers attaching to specific cell receptors. The intensity of proliferation and apoptosis could also be quantified using corresponding assays [67,68]. Flow cytometry and specialized assays also discriminate among cell populations according to the cell cycle stages [69,70]. Our research has shown that using flow cytometry and corresponding assays to study apoptosis and cell distribution over the cell cycle stages can achieve practical results in assessing cell population aging and rejuvenation [33,34,35,70,71,72,73]. Additionally, basing main conclusions about the cell age profiles on the specific stages, such as synthesis and proliferation, is quite reasonable. The entire cell cycle is known to transit two of these main stages significantly faster than all others. Thus, a percentage of cells in these cell cycle stages (immediately preceding division and immediately following it) reflects how young or old this cell population is at a given time [71,72,73]. Using cells from cultures or extracted from live tissue, it is possible to assess and quantify the changes in the percentage of cells in the synthesis and proliferation phases after specific interventions or over time.

It is worth mentioning that flow cytometry yields statistically- analyzed data, as it uses samples with large numbers of cells (typically 10⁵ to 10⁶ cells), providing not only average and cumulative values but also corresponding distributions. Although this method does not offer accurate distributions of the cells over their age, it gives reliable data on the percentage of young cells in the studied cell population. Thus, monitoring changes in the portion of the cells in the synthesis and proliferation stages, together with the assessment of apoptosis intensity, could be successfully used to validate the cell-level models and find substances with desired anti-aging and rejuvenation potential. Substances with such properties should simultaneously intensify proliferation and apoptosis (supporting the proliferation niche homeostasis) and increase the proportion of the cells in the synthesis and proliferation phases [33,34,35,70].

2.4 Aging and Rejuvenation in the Simplified Dynamic Model of Homeostatic Proliferation Niche

As mentioned above, the homeostatic mechanism of maintaining cell numbers in a selected volume requires two feedback mechanisms: one- responsible for adding newly born cells (proliferation), and another- responsible for removing old, senescent, and damaged cells (apoptosis). These two feedback mechanisms should be synchronized to maintain the cell numbers in the selected tissue volume, as in mature organisms. Figure 3 illustrates the concept of aging/rejuvenation in the homeostatic proliferation niche by an intuitive model of a semi-elastic tank with two valves and water flowing through it (Figure 3a). Input water flow will represent positive feedback provided by proliferation (adding cells), and drainage will mean negative feedback provided by apoptosis (removing cells). The age of cells in this model is expressed by local temperature (hot – young cells, cold – old ones). It is clear that to prevent the tank from changing its shape, the water flow rate through both valves should be kept identical, and the valve throughputs should be synchronized.

Figure 3 Simplified water tank and valve model of the proliferation niche. Incoming flow represents proliferation, draining- apoptosis and temperature represents cell age. Warm water –young cells; cold water- old ones. Cooling the water while it flows through the tank means increasing age with time.

Hot water is injected through the input valve and cools gradually when flowing through the system. Suppose that both valves are initially set to 50% of their throughput capacity, which produces a specific temperature distribution in the tank (Figure 3b). Synchronous reduction of the flow to 25% leads to increased time water spending in the system, resulting in more intense cooling (Figure 3c). The temperature of the exiting water decreases and the water in the tank ‘on average’ cools. Synchronous opening of the valves to 75% of their capacity increases the flow rate; water spends less time in the system and cools less (Figure 3d), increasing the average temperature in the tank. Note that the average water temperature in the tank does not entirely represent water cooling and warming, depending on the flow regulation. The corresponding water temperature distribution through the tank volume is much more representative.

This simple model already helps to form a concept of the balance between the capacity and efficiency of the involved feedback branches in the homeostatic niche model. Figure 3 illustrates the case of regulation when both valves can be opened at least up to 75% of their maximum capacity (100%). The efficiency of the feedback mechanisms in all cases shown in Figure 3b-d changes without reaching its capacity limit. It can be proposed that the cooling of the water tank system (equivalent to cell population aging) is represented by two simultaneous trends: decreasing maximum capacity (potential) of the valve throughput and falling efficiency (currently set throughput). Modeling results conclude that aging systems have simultaneously lowered maximum capacity (potential) of the regulating mechanisms and, although supporting homeostasis, are working with lower and lower efficiency. The simplified water tank model of the old proliferation niche implies that both valves can no longer be opened up to 100% (for example, their current limit is only 60% of their initial, young state). Moreover, they are now both open only up to 25%. This means no interaction with the existing regulation mechanism can increase its efficiency above 60%. Increasing efficiency above this limit cannot be achieved in the present state and needs significant changes in the system itself, which in the present simplified model could mean valve and maybe even pump replacements.

It is possible to use this simple model to illustrate the consequences of malfunctioning feedback mechanisms caused by the loss of their synchronization (Figure 4). When the feedback mechanisms are synchronized, decreasing (increasing) efficiency in one of the mechanisms necessarily causes matching changes in the efficiency of another one, independent of which mechanism changed first (Figure 4a). When synchronization is lost, one of the valves cannot open (close) as much as is needed to maintain a stable tank volume (equivalent to cell numbers in the niche). A low relative efficiency of drainage leads to the uncontrollable expansion of the tank (equivalent of neoplasm formation, Figure 4b). Low relative efficiency of the input flow leads to shrinkage of the tank (dystrophy, Figure 4c). In both cases, the limits in which the system will keep its integrity (tank will not rupture or collapse) are determined not by the regulation mechanisms but by the properties of the tank and pumps supplying the water.

Figure 4 Simplified water tank model of the proliferation niche. Desynchronization of the input and output valves.

Returning to the analogy with quantum-level models, it is worth remembering that temperature is a macroscopic scale (continuous matter) cumulative parameter reflecting average kinetic energy, in the above example, of water molecules. In physics, certain approaches allow linking the properties at atomic and molecular levels to those at the macroscopic condensed matter levels. Adequate models of such systems must account for many elements and ‘quantized’ elementary events. This understanding has provided the initial ideas for conceiving the approach used in the described modeling. All cell-level events are quantized: a single cell can only divide forming two cells, and a single apoptosis act removes one cell. Moreover, cell life is governed by biochemical signaling: a single biomolecule reaches a cell membrane, triggering a particular event. Properties defined at cell-level modeling should be coupled with the properties at the level of significant cell numbers representing condensed-matter situations. A statistical approach can be effectively used in modeling the processes at the levels of large numbers of cells, giving many valuable results. However, the proper transition from the cell level to that of a human tissue, organ, or body is equivalent to a leap from single molecule properties towards the properties of continuous matter. Indeed, average human tissue contains about 10⁹ to 10¹⁰ cells in one cm³. One milligram of a bioactive substance with a relatively small molecular weight of 1 kilo Dalton (10³ au) has roughly 10¹⁷ molecules. In this paper, we concentrate on the conceptual model level and the modeling software (computer simulator) operating at modest element number levels not exceeding 10⁵. This still allows for valuable conclusions using basic statistical approaches. However, it is already clear that using cell-based models for making predictions at the level of human tissue, organ, and body needs the introduction of advanced approaches allowing operation with element numbers up to 10¹⁵ to 10²⁰.

3. Methods

3.1 Basics of Modeling Homeostasis in a Proliferation Niche

Earlier, we defined homeostasis in the proliferation niche as maintaining cell numbers in a selected volume. This led to a hypothesis that cell numbers in a selected volume (cell concentrations) and concentrations of some regulatory chemical factors could provide the desired link and maintenance of cell numbers in proliferation niche can be achieved by maintaining the concentration of such a regulating substance. However, this can only happen if this substance is a crucial regulator of both positive and negative feedback branches, e.g., proliferation and apoptosis. It should be noted that the modeling approach presented below results from lengthy discussions, analysis of the literature, formulating hypotheses and targeted experimental studies, inspired by the progress and the modeling results. As a result, it was possible to develop the search strategy for such a regulatory compound, separate it from the population of dividing cells, and prove the presence of cell number maintenance and its active role in it experimentally in vitro [33,34,35]. Consequently, the statements formulated below have already found specific substantiation.

Figure 5 presents a sequence of events leading to the regulation through the concentration of a critical chemical substance (regulating molecule) that became a basis for constructing a modeling program. Corresponding cell division only happens when a regulatory molecule reaches a corresponding receptor of the cell capable of proliferation. The cells that are not ready or incapable of the division would not react to the triggering molecule. Moreover, although many cells may be capable of division, only a certain number will divide, which would be defined by the availability of the regulator molecules. There may be insufficient regulator molecules for all currently ‘opened’ receptors to accept them. In addition, although there is a difference between the elementary proliferation act involving a stem cell and its daughters (a stem cell generates its copy, and its daughters divide, producing two identical cells), we can still model the elementary division act universally as we are constructing a model monitoring only cell age and numbers. In such a case, one cell is present before division and two cells after. In later versions of the model, it was possible to consider such differences and the potential loss of the proliferation capacity within consecutive generations of the daughters of the initial stem cell.

Figure 5 (1) Regulatory molecule R reaches the stem cell, forcing it to generate a copy with high proliferation potential (2). A regulatory molecule causes a cell with proliferation potential to divide (3) and generate two young cells (4). Stem cells can become dormant (5) and, at some point, can be woken up by a regulatory molecule (6). Regulatory molecules provide coupling between proliferation and apoptosis regulatory branches (7), causing apoptosis (8).

Synchronization of the cell division and apoptosis in the populations with non-pathologically dividing cells leads to the suggestion that the same regulation molecule should trigger and control both mechanisms, and its concentration has the role of a regulating factor.

The software simulator developed from the described conceptual model is based on the ‘reaction equations’ for the elementary acts of proliferation and apoptosis (see Supplementary Material S1 for details). Simulating software contains several cyclic computations over a few values defining numbers of cells with specific ages and numbers of regulatory molecules with their age. Each of the cycles represents certain time intervals with several consecutive actions and control elements. A detailed flow diagram illustrating the corresponding computing stages of the simulator is given in Supplementary Material S2. The simulator compares the number of cells that can divide to the number of available regulatory molecules. The corresponding number of cells that will divide and the number of regulatory molecules that will be consumed and synthesized is calculated using corresponding efficiency coefficients. The number of cells that the apoptosis should remove is calculated. Related data sets are updated, and critical parameter values are sent to the outputs and graphs before the next computation cycle goes forward. Newly born cells and newly synthesized regulatory molecules are assigned age equal to zero, and incorporating them into the data sets is straightforward. The simulator also generates particular subsets that calculate how many cells and regulatory molecules of a specific age will be removed based on the overall numbers and corresponding age-dependent probabilities. Connected modules using an exponential time decay function for the regulatory molecules and a weighting function favoring apoptotic removal of the cells with increasing age generate related subsets for properly updating main data sets.

4. Results and Discussion

4.1 Model Development and Preliminary Validation

One of the features incorporated into the later versions of the software is taking into account the so-called Hayflick limit. It is generally accepted that each cell has a maximum number of possible divisions, after which it only ages and loses the possibility to divide (e.g. [74,75,76]). It is assumed that this limit could differ for different cells, so corresponding values are used as initial parameters and varied in the simulation trials. In addition, specific other data and parameter values necessary for generating an adequate simulator were initially unavailable. Thus, certain flexibility was introduced into the software, and virtual experiments were carried out to find the parameter windows that allow the system to reach a steady state of cell numbers. Trials were performed using different initial conditions, searching for the parameter ranges leading to the system's stability state at a particular time.

Several virtual experiments for ‘what if’ cases allowed for setting up specific experiments to prove initial assumptions. However, some assumptions so far did not allow for straightforward validation. For example, simulations with the catalytic- and synthesis/consumption types of elementary reactions result in very similar system dynamics. Thus, it was impossible to suggest straightforward experiments to facilitate the choice between these two modes. On the other hand, several conclusions following the modeling were confirmed by direct experiments. For example, direct experiments with cell cultures confirm the dose saturation effect with the influence of adding regulatory substances to the system [33,34,35]. Even with relatively generic modeling, it becomes clear that with the intense proliferation following the cell division cascade, the main burden of building newly born cell numbers falls upon consecutive generations of daughter cells (Figure 6). Although it is known that stem cells in the mature organism divide quite rarely [71], cascades of the dividing daughter cells of different generations should be potent enough to support the balance between needed proliferation and apoptosis rates. Thus, the properties of stem cells initiating division cascades are quite important to account for. One such factor is the known ability of stem cells to become dormant, excluding themselves from active metabolism and closing corresponding characteristic receptors (e.g. [77,78,79,80,81,82]). Moreover, one of the possibilities is that the number of stem cells in mature tissues may be constant or only slowly decaying with age, but they are mainly dormant. They are woken up only when there is a need to intensify proliferation. We have hypothesized that within our model, the corresponding ‘wake-up call’ could be produced by an increasing concentration of the regulatory molecules, on the one hand, connected to the increased proliferation intensity, and on the other hand, demanding more robust support from more activated stem cells. Figure 6 schematically illustrates the cell division cascade considering the different aspects described above.

Figure 6 Cell division cascade representation used in the modeling process.

For modeling, we are simplifying the variety of cell types to only a few categories according to their proliferation potential. Modern cell biology distinguishes a number of such categories depending on different functional capacities and proliferation potential of the cells (totipotent/omnipotent, pluripotent, multipotent, oligopotent, unipotent, e.g. [83,84,85]). Some researchers also distinguish additional sub-categories depending on the stem cell origin, localization, specialization, etc. For the basic model, we use only three categories of cells about division: stem cells, multipotent stem cells, and somatic cells, which differ in their proliferation capacity, reducing the number of parameters involved. In the most basic approximation, both stem and multipotent cells can divide, although with different efficiency, and somatic cells cannot. This is one of the limitations inherent in the existing model. However, it already has a specific predictive capacity and can be used for further expansion, including more details.

As follows from earlier discussions, corresponding regulatory molecules should be naturally present in the populations of non-pathologically dividing cells. These molecules should simultaneously and synchronously intensify proliferation and apoptosis without decreasing cell vitality. This allowed the formulation of the criteria for searching for such molecules. Experiments with the supernatant extracted from the corresponding cell cultures confirm that molecules with such properties are present in specific chromatographically isolated fractions [33,34,35]. Further experiments confirmed that depletion of the corresponding fraction in the supernatant of the cell culture decreases proliferation intensity. At the same time, adding this fraction to the other culture increases the proliferation intensity and the concentration of the related molecules in the supernatant (autoinduction). The corresponding rejuvenating activity of the substances could be tested using cell cultures (e.g., human fibroblasts and mononuclear blood cells, both laboratory or ex-vivo ones). This allowed us to measure the lifetime of corresponding molecules in buffer solutions experimentally. The decay of activity measured by the intensification of the proliferation of test cell cultures is close to exponential, with a characteristic decay time of 36 hours [33,34,35].

Several additional factors were also considered and implemented in the simulation software. To predict the results of therapeutic applications for further comparison with experimental data, possibilities to add or remove specific amounts of different cells and regulator molecules during the simulation were introduced. We have also incorporated an equivalent of the mitogen action upon the cell population. According to the developed model, mitogen should force a certain number of cells with proliferation potential to divide immediately (e.g. [86,87]) without any involvement of regulating molecules. As it is experimentally known, mitogen application destroys the cell number maintenance mechanism: mitogens enhance only proliferation without accelerating apoptosis. In the simulations, we have also introduced the possibility of somatic cell induction into the multipotent state. There is a lot of experimental data pointing out that such installation could happen under the influence of so-called Yamanaka factors (e.g. [88,89]), and corresponding interventions may have high rejuvenating potential [90,91,92,93].

The simulation starts with generating initial age distributions (initial arrays) of all elements and continuously runs until it is forced to stop. We have also added one extra output parameter, representing the energy the system consumes. All elementary reactions (acts of division, apoptosis, regulator molecule consumption, and generation and stem cell wakeup) were assigned relative energy consumption shares that were summed up at any particular simulation moment. These shares are somewhat arbitrary, as are many other unknown parameters, and so far, are mainly based on reasonable assumptions and available data. The present simulation software version was developed based on the above concepts and corrected according to the experimental data, and it incorporates a few possible elementary reaction equation choices. Supplementary Material S1 includes a list of main postulates integrated into the model and resulting simulation software code, specifically mentioning which are axiomatic and which are supported by experimental data. Supplementary Material S2 describes the modeling program and its functional flow diagram. Supplementary Material S3 provides the description of the controls, parameter displays, and graphic displays realized in the simulator and shows the view of the virtual control panel of the simulator.

4.2 Main Results of the Simulations: Youth, Maturity, and Old Age of The Simulated System

Several trials were required before it was possible to outline the window of initial condition parameters leading to the desired homeostatic stabilization of the system. It was concluded that the most critical is the interplay of the corresponding reaction rates and characteristic decay time of the regulator molecules (corresponding rates are defined as a fixed number of time steps between the initiation of the process and the availability of its result). After establishing this balance, it was found that homeostatic conditions are reached for a wide range of initial conditions and critical parameters. Moreover, it became clear that the competition between the rates of production and consumption of the regulator molecules and the balance between the rates of cell proliferation and apoptosis (represented by corresponding delays of the elementary acts), in a way, play a role of the ‘measuring mechanism’ for homeostatic feedback.

Figure 7, Figure 8 and Figure 9 illustrate some of the significant simulation results. All of these simulations were carried out with the same initial conditions. The gradually decreasing probabilities with predefined overall element numbers represented all initial age distributions. For clarity of the discussion, all output graphs are represented by the trend curves, excluding the influence of stochastic contributions introduced in the latest version of the simulation program. To illustrate the difference, one of the curves in Figure 8 shows the trend line over the corresponding output graph produced by the simulator with the noise due to the presence of stochastic contributions.

Figure 7 Output graphs showing the changes in critical parameters characterizing the ‘full life cycle’ of the homeostatically-stabilized cell population from ‘early infancy’ to ‘death’ without external interventions and with no ability for somatic cells to divide. Horizontal axis: time (number of simulator cycles). Vertical dotted lines indicate corresponding characteristic age intervals: infancy (up to 200 simulation cycles), youth (200 to 680 cycles), maturity (680 to 1950 cycles), old age, and decline (from 1950 cycles).

Figure 8 Output graphs showing the changes in critical parameters characterizing the ‘life’ of the homeostatically- stabilized cell population without external interventions and with a specific capacity of somatic cells to be promoted into a multipotent state and to divide (corresponding graph shows the number of somatic cell divisions). Horizontal axis: time (number of simulator cycles). Vertical dotted lines indicate corresponding characteristic time intervals: infancy (up to 150 cycles), youth (150 to 310 cycles), maturity (310 to 680 cycles), intense aging and decline (from 680 cycles), and ‘end of life’ (at 2100 cycles).

Figure 9 System reaction to (a) ‘stem cell therapy’ (addition of the young stem cells) and (b) addition of the new regulator molecules into the system at different states. Arrows mark corresponding intervention moments. Horizontal axis: time (number of simulator cycles).

Figure 7 presents the dynamic changes in the number of cells (including the number of currently active stem cells), the number of dying cells (due to apoptosis), the number of regulator molecules, estimated energy consumption, and the number of somatic cells. The corresponding graphs represent the ‘full lifetime’ of the modeled system from ‘early infancy’ to ‘death.’ The current simulation was performed under conditions that did not allow the induction of somatic cells into the multipotent state. The simulations also produce graphs showing dormant stem cell numbers, numbers of ‘dying’ regulator molecules, and other parameters and corresponding age distributions not represented in the current figure but quite valuable for the analysis of the system dynamics (for reference- see Figure S3.1 in the Supplementary Material S3).

Figures 7 and 8 illustrate typical simulation results in two different cases when regulator molecules are consumed in each reaction. Figures 7 shows the case when stem cells can age, and apoptosis is active for all cell types, including stem cells; cells have a limited number of divisions (Hayflick limit), and somatic cells cannot be induced to the multipotent state. Although we focused our conceptual models on mature systems with established homeostasis, corresponding results represent all stages of system development from a young age through maturity to old age and death.

One of the most striking results the modeling produces is distinguishing the stages in the system’s aging. These stages have quite characteristic parameter dynamics and can be easily indicated in the graphs presented in Figures 7 and 8. Vertical dotted lines mark corresponding ‘age zones’ in the charts.

Youth and infancy stage: all stem cells are active (none are asleep); numbers of somatic cells and regulator molecules grow, pass the maximum, and start to fall; numbers of multipotent cells grow and stabilize. Simultaneously, apoptosis intensifies (the number of dying cells increases). Consumed energy reaches the maximum and falls, reflecting active proliferation and apoptosis. Two dotted lines within the youth stage (more clearly seen in Figure 7) represent the early (infancy) and later youth stages. Death of somatic and multipotent cells starts with a delay: all cells in the initial population (initial distributions) are of relatively small age (relatively ‘young’).

Maturity is characterized by stable numbers of somatic and multipotent cells, regulator molecules and apoptosis intensity (dying cell numbers), and energy consumption. Stem cells continuously become dormant, decreasing the number of active stem cells. In this period, the overall number of cells stabilizes; the number of regulatory molecules in the system supports adequate numbers of the multipotent cell divisions needed for compensation of cell removal by apoptosis. The system enters homeostasis; all major parameters stabilize, including the number of falling asleep and awoken stem cells.

Old age is characterized by stem cells starting to lose their proliferation potential, reaching the Hayflick limit of their division numbers. Consequently, the number of multipotent cells decreases, causing, in turn, a decreasing number of newly born somatic cells. The number of regulator molecules is still stable, and apoptosis continuously removes the cells, causing their numbers to fall. Energy consumption also falls as the intensity of the processes in the system decreases.

Death: after about 22000 simulation cycles, the number of somatic and stem cells becomes zero, and energy consumption falls almost to zero. Some regulator molecules and multipotent cells remain, but after a few more simulation cycles, all critical parameter values become zero.

Figure 8 illustrates the results of a typical simulation with the same parameters as described above and illustrated by Figure 7 - except somatic cells can be induced into the multipotent state and can divide with nonzero probability. Qualitative changes due to this difference are not significant. The modeling results also clearly illustrate a few infancy-youth stages of the system described earlier. The time scale for the number of somatic cell divisions graph is extended to demonstrate the increased impact of allowing division for somatic cells. After roughly the 2250^th simulation cycle somatic and stem cells stop to divide. Shortly after this, all changes cease and the system ends its life.

4.3 Main Results of the Virtual Experiments

The many controls incorporated into the latest simulator versions allow for virtual experiments with adding or removing regulator molecules, adding external stem cells, eliminating cells, or adding a mitogen. Adding mitogen forces cell division without affecting apoptosis, and this type of division is not followed by the synthesis or consumption of the regulator (corresponding experimental indications to this exist). Adding external stem cells (with the current version of the program, all of them are active and have zero age) mimics stem cell therapy. Adding a certain amount of regulatory molecules according to the present conceptual model should represent the rejuvenation of the system without destroying homeostasis. Virtual experimentation allows for validation of the model by comparing corresponding modeling results with data from expressly set experiments or from available literature. So far, all known data support the conclusions following from the modeling.

Figure 9 illustrates the simulated results of the ‘virtual therapy’. Figure 9a: the system is undergoing an intense virtual stem cell therapy (addition of a fixed amount of active stem cells in several consecutive interventions) without adding regulator molecules within the mature state of the system. This immediately intensifies the cell division cascade and causes an increase in multipotent and somatic cell numbers. Simultaneously, the amount of the available regulator molecules drops as they are consumed in the intensified proliferation. However, homeostasis in the system is active, cells continue to divide, and the amount of regulator slowly grows to the pre-therapy state.

Figure 9b: the system is undergoing ‘rejuvenation therapy’ by adding regulator molecules at different system maturity stages (the bottom image separately shows the youth period and onset of the maturity period). Adding regulator molecules to the ‘young’ system produces relatively weak and only short-lived changes (first three arrows). At that stage, proliferation and ‘own production’ of the regulator molecules are quite active, and almost all capacity for the intensification of both proliferation and apoptosis is already utilized (interventions up to 4,000 simulation cycles, Figure 9b). Thus, additional regulator molecules are excessive, will not be used, and die within a relatively short time. A small rejuvenation impact (increased share of young cells) can be seen in the corresponding cell age distributions, but this ceases within a relatively short time.

System reaction in the early maturity state (after ca 10000 simulation cycles, Figure 9b) is quite different. Intensification of the proliferation (and matching apoptosis) lasts significantly longer. This happens because at the system's mature state, intensifying proliferation and increased amount of regulator molecules wakes dormant stem cells, which initiate new proliferation cascades. As long as the amount of regulation molecules is high, awakened stem cells do not go to sleep, and support intensified proliferation. With the time, the system returns to the steady state condition. This is due to degradation effects (in the model, this is reflected by the limited probabilities related to cell division, regulator production, and limited regulator lifetime).

Corresponding ‘therapy interventions’ with only additional regulator molecules in the old system state are much less effective and - at some point close to the system expiration - are becoming completely useless, mainly because almost all stem cells have reached the division limit, and such intervention cannot lead to new proliferation cascades. At the same time, the addition of external stem cells with unused proliferation potential can prolong the system’s lifetime, especially when accompanied by the addition of regulatory molecules. However, such interventions cannot extend the system’s life forever, and at some point, they also become ineffective, and the system expires, having exhausted its proliferation potential.

4.4 Prospects for Further Model Development and Validation

Some of the conclusions from modeling are experimentally proven at the cell culture level and in vivo, and thorough experimental research on many related issues is underway. Critical cell culture experiments were carried out with human dermal fibroblasts (both commercial laboratory cultures and ex-vivo ones produced from skin samples taken from the gluteal area of healthy volunteers) and from the pooled ex-vivo mononuclear blood cells (extracted from venous blood of healthy volunteers). Both cell types indicate skin aging [21,22] and rejuvenation. Search criteria based on the formulated hypotheses and following the modeling results allowed the extraction of substances with supposed activity from the supernatants derived from the cultures of actively dividing cells [33,34,35]. Experiments with healing skin lesions also confirm that extracted substances non-pathologically intensify cell proliferation and tissue regeneration, with minimal scar formation in many cases. At the same time, details of the corresponding mechanism of the regulation compound interaction with the cells are not fully understood. Existing models and simulating software are extensively used to suggest specific experiments that can further prove the validity of the proposed approach. As far as we can see, this is a primary way forward for the corresponding experimental validation of both conceptual and computational models.

At the same time, significant limitations of the cell and regulator molecule numbers that the existing simulating software can effectively handle are critical to expanding the modeling towards the level of human tissues, organs, and the body. In addition, it became clear that further progress within currently used approximations is significantly limited, as the model does not include spatial effects. The present version of the model cannot consider the diffusion of the regulating substance, mass transfer and migration of the cells. Due to these reasons, the current model is poorly applicable to the circulatory system, where the critical immune cells are born in one place and transported by blood and lymph flow to other sites.

Although the current general modeling approach provides many valuable conclusions, the complexity of the needed procedure far exceeds the possibilities of the presented model with a single degree of freedom (time) and its current opportunity to handle only 10⁵ to 10⁶ elements (cells and regulator molecules). Thus, significant progress is connected to both further rectifications of the existing model and involving different modeling approaches. In this respect, the most promising could be multi-dimensional agent-based models (e.g. [94,95]) in two or three dimensions, possibly incorporating spatial constraints and mass transfer. Another good approach can be connected to solving discretized differential equations, moving from time-continuous functions towards quantized elementary reactions.

5. Conclusions

A cell-based approach is used to model critical processes related to aging and rejuvenation. Simulator software was compiled using a conceptual model of homeostasis in the proliferation niche describing cell number maintenance under simultaneous and synchronized action of cell proliferation and apoptosis. An adjusted definition of cell population age is proposed to allow adequate computer modeling and quantitative comparison. It is shown that the cell age profile (age distribution over cell age) is sufficient for generic modeling and analysis of the issues related to aging and rejuvenation, and its critical parameters can be measured experimentally using flow cytometry.

The corresponding project involved intense modeling and validation of the main results using literature data and specific experiments. So far, no critical deficiencies of the chosen approach or significant discrepancies with known experimental data have been found. The following results can be formulated:

1. Aging of a homeostatic system can be defined as simultaneous ‘degradation of regulation mechanism capacity’ and ‘loss of its efficiency.’ Corresponding external ‘rejuvenating interventions’ aimed at increasing the ‘efficiency’ can only give results within the existing ‘system capacity.’

2. Cell number maintenance in such a system is governed by two simultaneously acting synchronized mechanisms: proliferation and apoptosis. Maintaining the concentration of a certain regulating substance could be a possible control of homeostasis.

3. Such a system has a specific ‘stability window’ over the initial conditions and critical parameter values and can return to its homeostatic state after external disturbances.

4. Simulation results allow interpreting system development stages as ‘youth,’ ‘maturity,’ ‘old age,’ and ‘death’ with an adequate representation of critical parameter dynamics.

5. The developed simulator allows for modeling results of ‘external interventions’ using stem cells, regulator molecules, Yamanaka factors, and mitogens.

6. Simulations show relatively low rejuvenating effect from therapeutic interventions using stem cells and regulator molecules at a young age, much higher efficiency and prolonged action at a mature age, and small or negligible effects at ancient age.

7. Simulations show mutual enhancement of the rejuvenation efficiency with simultaneous application of stem cell and regulator substance therapy.

8. Simulations show a ‘dose saturation effect’ of both stem cell and regulator substance interventions (limiting concentration depends on the system's current status).

The scope of the present paper does not allow for proper discussion of the validation of modeling concepts and simulation results, so a few separate publications are in the pipeline.

The corresponding simulator can be further developed by incorporating additional parameters and effects. However, several apparent deficiencies limit further progress in modeling reaching a tissue-organ-body level. Along with increasing numbers of elements (cells, regulator molecules) adequate for representing tissues and organs, the incorporation of spatial effects and more potent modeling approaches is needed. We hope to attract wider attention from specialists in advanced modeling and big data to the studies of aging and rejuvenation in general and cell-level modeling in particular.

Author Contributions

A. Koptyug, Yu. Sukhovei and V. Kozlov; conceptualization, and methodology; Yu. Sukhovei, I. Unger and E. Kostolomova; cell culture preparation and cell culture experiments; Yu. Sukhovei, I. Unger and A. Koptyug; data handling; Yu. Sukhovei, A. Koptyug and V. Kozlov; formal analysis; writing, review and editing. All authors have read and agreed to the published version of the manuscript.

Competing Interests

The authors have declared that no competing interests exist.

Additional Materials

The following Supplementary Materials are uploaded together with the main text of this paper.

1. Supplementary Material S1: Simplified reaction equations used for modeling elementary acts of proliferation and apoptosis; main statements used during the development of simulator software code.

2. Supplementary Material S2: Description of the modeling program, and flow diagram of the simulator software code.

3. Supplementary Material S3: Virtual control panel of the simulator program, brief description of the output parameters and graphs.