인간의 수명 추정 146~194(156세) 2026년 네이처

[문형철]

다른 노화 원인을 모두 없애도,

체세포 DNA 돌연변이(엔트로피성 정보 손실) 만으로 인간 중앙수명은

이론상 약 146–194년(대표값 156년)에 묶인다.

병목은

거의 교체되지 않는 뇌 신경세포와 심근세포다.

---

연구 질문

노화의 hallmarks 중 상당수는

원리상 되돌릴 수 있다(텔로미어, 노화세포, 일부 대사·염증 등).

반면

체세포 돌연변이는

비가역적 정보 손실에 가깝다.

질문:
가역적 노화를 전부 치료했을 때,

돌연변이만으로 남는 수명의 상한은 얼마인가?

이는

“인간이 영원히 살 수 있는가”가 아니라,

엔트로피적 상한을 수치로 잡으려는 모델링 논문이다.

---

방법 (점진적 모델)

인체를 신뢰성 이론의 직렬·병렬 시스템으로 본다.

핵심 장기(뇌, 심장, 간, 폐 등) 중 하나라도 임계 기능 이하로 떨어지면 사망.

단계를 나눠 복잡도를 올린다.

모델내용

I	노화 없음. 사망위험을 30세 스위스 코호트 수준으로 고정 → 중앙수명 1759년, 최대 약 2.9만 년
II	후분열 세포(신경·심근)에 돌연변이가 쌓여 세포가 죽고, 잔여 세포가 임계치 아래로 떨어지면 장기 부전
III	분열·재생 조직(간세포, 간 전구세포, 기도 기저세포 등)의 교체·고갈까지 포함

돌연변이 축적률은 조직별 시퀀싱 데이터(SNV, indel)에서 추정.
치사 돌연변이 확률 plethalp_{\text{lethal}}plethal​은 하한·상한을 두고 민감도 분석을 했다.
여러 장기 사망시각의 상관은 Fréchet–Hoeffding 경계로 묶어, 중앙수명을 146–194년 구간으로 제시한다.

가정(중요): 장기이식 없음, 돌연변이 축적률을 낮추는 개입 없음, 30세 배경사망률 고정.

---

핵심 결과

1.  
2. • 신경세포만 돌연변이로 죽게 하면 중앙수명 ~194년, 최대 ~557년
   • 심근세포만이면 중앙 ~208년, 최대 ~868년
   • 둘을 합치면 대표 중앙값 약 156년
3. • 간: 전구세포 지원이 있으면 10만 년 시뮬레이션에서 부전 거의 없음
   • 줄기세포 없이도 수만 년 단위 유지 가능
   • 기도 기저세포는 증식 고갈로 수천 년대에 한계
4.  
5.  

---

이론적 위치

• Szilard(무작위 염색체 hit) → Orgel(오류 파국) 계열을 정량 상한으로 옮긴 작업.
• “엔트로피적 상한”이라는 제목은, 돌연변이를 되돌리기 어려운 정보 손실로 본다는 뜻.
• 종 간 돌연변이율 ↔ 수명 반비례(Cagan 등)와도 맞닿아 있다. 다만 이 논문은 인간 내부, 조직 비대칭을 수치화한 점이 다르다.

---

한계 (저자·독자가 같이 봐야 할 점)

• 실험이 아니라 파라미터 기반 시뮬레이션.
• 치사 돌연변이 확률, 장기 부전 임계치, 장기 간 의존성에 민감.
• 기능 저하(시냅스, 흥분-수축 결합, 미토콘드리아, 후성유전)를 “세포 사망”으로 단순화.
• 클론 조혈, 암, 염증, mtDNA, 면역노화는 본 모델의 핵심 항이 아님.
• “156년”은 미래 예측이 아니라, 다른 노화를 모두 껐을 때의 가상 상한.

---

임상·통합의학으로 읽으면

사용자 관심(미토콘드리아 ROS, ISR, 대사, 재생, 암 프로토콜)과 연결하면 다음과 같다.

• 병목 장기는 뇌·심장. 신경·심근을 지키는 전략(대사 안정, ROS/철 항상성, 혈관·당대사, 수면·미주신경, 케톤/BHB 등)은 “수명 상한”이 아니라 현재 건강수명을 올리는 쪽에 더 직접적이다.
• 간·점막처럼 재생되는 조직은 돌연변이만으로는 잘 안 무너진다. 실제 노화·암·섬유화는 재생 실패 + 염증 + 클론 선택이 겹칠 때 문제가 된다.
• 돌연변이를 지금 임상에서 “지우는” 수단은 없다. 현실적 레버는
  • 돌연변이 속도를 줄이는 쪽(산화·복제 스트레스, 만성 염증, 흡연 등)
  • 손상 세포의 기능 유지·선택적 제거
  • 재생 가능한 구획의 줄기세포/전구세포 소진 지연
• “노화를 고치면 1759년”은 모델 I의 사고실험일 뿐이다. 이 논문의 메시지는 오히려
  가역적 hallmarks를 다 잡아도 정보 손실 항이 남는다는 점이다.

세 그림은 논문의 전체 논리 한 줄입니다.
Hallmarks 중 유전체 불안정성만 남긴다 → 장기별로 세포가 얼마나 버티는지 본다 → 인체를 직렬 시스템으로 합치면 중앙수명 156년.

---

Figure 1. 질문과 모델링 틀

a. 노화 hallmarks 바퀴에서 genomic instability(체세포 돌연변이) 만 남긴다. 나머지는 “치료됐다”고 가정.

b. 체세포 돌연변이 부담은 나이와 거의 직선으로 는다.

c. 생존곡선 세 갈래

• 분홍: 실제 노화 → 최고수명 122년에서 끝
• 주황: 돌연변이만 남긴 가상 노화 → 122년보다 오른쪽, 끝이 ?
• 청록: 노화 없음 → 거의 안 떨어짐

이 논문이 푸는 것이 주황 곡선의 ? 이다.

d. 인체 = 신뢰성 시스템. 뇌·폐·심장·간 중 하나라도 임계 이하로 가면 사망(직렬).

e. 장기 부전 정의: 세포 수가 임계선(점선) 아래로 떨어지는 시점 = organ death.

f. Incremental approach

• Model I: 사고사만. λ=λBG\lambda = \lambda_{BG}λ=λBG​ (상수). 노화 없음.
• Model II: 후분열 세포가 돌연변이로 죽음. λ=λBG+λmut\lambda = \lambda_{BG}+\lambda_{mut}λ=λBG​+λmut​. 교체 없음.
• Model III: 분열·줄기세포 보충. λ=λBG+λmut+replication\lambda = \lambda_{BG}+\lambda_{mut+replication}λ=λBG​+λmut+replication​.

---

Figure 2. 모델별 숫자 (그림에 찍힌 값)

Model I (a–c)
30세 사망위험으로 “동결”하면

• 중앙수명 1,759년
• 최대수명 29,221년
  c는 몇 살 위험을 고정하느냐에 따라 곡선이 수천~수만 년으로 밀리는 민감도.

Model II — 뇌·심장 (d–f)
교체 없는 세포. 잔여 세포 X(t)X(t)X(t)가 임계까지 감소.

뇌심장

임계 세포비율 XcritX_{crit}Xcrit​	60%	55%
임계 도달 tcritt_{crit}tcrit​	198년	212년
돌연변이만, 중앙 tmedt_{med}tmed​	194년	208년
돌연변이만, 최대 tmaxt_{max}tmax​	557년	868년
배경사망만	tmedt_{med}tmed​ 1,755 / tmaxt_{max}tmax​ 29,147

병목은 뇌·심장.

Model IIIA — 간, 간세포만 분열 (g–i)
임계 20%. tcritt_{crit}tcrit​ 37,664년.
간+배경: 중앙 37,821년, 최대 69,309년.
간은 돌연변이만으로는 사실상 안 무너진다.

Model IIIB — 간 + 전구세포 LPC (j–l)
10만 년 안에도 부전 없음. tcritt_{crit}tcrit​ unknown, tmed>100,000t_{med}>100{,}000tmed​>100,000.

Model IIIC — 폐 기도 기저세포 BBC (m–o)
증식 고갈. 임계 23%, tcritt_{crit}tcrit​ 4,423년.
폐+배경: 중앙 4,359년, 최대 7,132년.
간보다는 짧지만 뇌·심장보다는 한 자리 길다.

한 줄: 재생 안 되는 조직이 수명을 끊고, 재생되는 조직은 돌연변이 단독으로는 수명을 못 끊는다.

---

Figure 3. 장기를 직렬로 합친 최종 곡선

a. 배경사망 λBG\lambda_{BG}λBG​ + 뇌 + 폐 + 심장 + 간 → 하나라도 실패하면 사망.

b. 생존곡선

• 분홍 Normal aging: 중앙 79년, 10−3%10^{-3}\%10−3% 근처 107년대 (현실)
• 청록 No aging (Model I): 거의 안 떨어짐
• 보라 뇌(Model II), 빨강 심장(Model II), 노랑 간+LPC, 파랑 폐 BBC
• 검은 굵은 선 Combined: 장기 결합
  • 중앙 tmed=156t_{med}=156tmed​=156
  • 하한–상한 146–194 (회색 밴드, 장기 실패 시각의 상관 불확실성)
  • 최대 tmax=470t_{max}=470tmax​=470 (상한 쪽 557)

해석:
다른 hallmarks를 전부 꺼도, 돌연변이 + 뇌·심장 병목만으로 중앙수명은 약 156년(대략 현재의 2배)에 걸린다.
122년 현실 최고수명보다 길지만, 1,759년 비노화 모델보다는 훨씬 짧다.
따라서 돌연변이는 강한 상한이지 유일한 원인은 아니다.

---

그림만으로 남는 메시지

1. 엔트로피(비가역 돌연변이)만으로도 영생은 안 된다.
2. 상한의 실체는 전신이 아니라 후분열 뇌·심장.
3. 간·폐 재생 구획은 이 모델에서 수명 결정항이 아니다.
4. 156년은 예측이 아니라 “나머지를 다 고쳤을 때”의 가상 천장.

1. nature  
2. npj aging  
3. articles  
4. article

Somatic mutations impose an entropic upper bound on human lifespan

Download PDF

• Article
• Open access
• Published: 25 June 2026

Somatic mutations impose an entropic upper bound on human lifespan

• Evgeniy Efimov, 
• Vlad Fedotov, 
• Leonid Malaev, 
• Ekaterina E. Khrameeva & 
• Dmitrii Kriukov 

npj Aging volume 12, Article number: 116 (2026) Cite this article

Save article

•  
•  
•  

Abstract

Somatic mutations accumulate with age and can cause cell death, but their quantitative contribution to limiting human lifespan remains unclear. We developed an incremental modeling framework that progressively incorporates factors contributing to aging into a model of population survival dynamics, which we used to estimate lifespan limits if all aging hallmarks were eliminated except somatic mutations. Our analysis reveals fundamental asymmetry across organs: post-mitotic cells such as neurons and cardiomyocytes act as critical longevity bottlenecks, with somatic mutations reducing median lifespan from a theoretical non-aging baseline of 1759 years to 156 years. In contrast, proliferating tissues like liver maintain functionality for thousands of years through cellular replacement, effectively neutralizing mutation-driven decline. Multi-organ integration predicts median lifespans of 146–194 years—approximately twice current human longevity. This substantial yet incomplete reduction indicates that somatic mutations significantly drive aging but cannot alone account for observed mortality, implying comparable contributions from other hallmarks.

Similar content being viewed by others

Somatic mutation rates scale with lifespan across mammals

Article Open access13 April 2022

Targeting the “hallmarks of aging” to slow aging and treat age-related disease: fact or fiction?

Article Open access15 July 2022

Mitochondrial somatic mutation and selection throughout ageing

Article Open access15 February 2024

Explore related subjectsDiscover the latest articles and news in related subjects.

• Cardiology
• Cell biology

Introduction

Contemporary geroscience pursues the ultimate goal of extending human healthspan and lifespan, attracting substantial attention1,2,3. However, the very phenomenon of biological aging awaits a unified theory, and scientists worldwide engage in large-scale projects aiming to discover reliable biomarkers4,5, unravel aging mechanisms6, and develop computational models7,8,9. These findings are compiled into knowledge compendiums like the Hallmarks of Aging10, revealing aging as highly complex, interdependent, and multifactorial.

Some hallmarks appear reversible, and we can envision molecular pathways for drugs capable of correcting specific aging-related consequences beyond merely slowing down their accumulation. For example, telomerase overexpression‎ reverses telomere shortening11,12, while senolytics target senescent cells13,14, albeit actual longevity therapies are yet to appear. Other aging phenomena represent pure information loss (or entropy increase)15,16,17,18. One prominent example is somatic mutations—DNA alterations that occur throughout lifetime in non-germline cells19,20. The hypothesis that somatic mutations drive aging was first formalized by Leo Szilard21, who proposed that random hits to chromosomes accumulate over a lifetime and progressively impair cell function. Leslie Orgel later extended this framework by demonstrating that translational errors and somatic mutations form a coupled positive feedback loop: faulty protein synthesis machinery generates mutator polymerases, accelerating DNA mutation rates, which in turn further erode translational fidelity22. Although still considered part of genomic instability10,23 (Fig. 1a, b), the impact of somatic mutations on aging and life expectancy has been debated24,25. Mounting evidence suggests somatic mutations are not mere clock-like changes26, but drastically affect health by altering epigenetic landscapes27, inducing within-tissue genome mosaicism, generating pathogenic clones, and enabling age-related diseases28,29,30,31,32,33.

Fig. 1: Proposed incremental framework for modeling organismal aging.

Full size image

a Among the hallmarks of aging, we focus on genomic instability which is particularly challenging to reverse. b Somatic mutations accumulate throughout the human lifespan and correlate strongly with chronological age. c Normal survival curve reaches zero at the maximum documented human lifespan of 122 years. Here we ask: What would the maximum achievable lifespan be if somatic mutations were the sole driver of aging? It would presumably be more than if all drivers of aging are considered (“Normal aging”) and less than in case of “No aging” (no mortality increase over time). d The human body can be modeled as a complex reliability system60,61 composed of parallel and series elements; failure of any critical element leads to death. e Schematic representation of the hepatocyte population dynamics: starting at full capacity (100%), cell population declines over time, and organ failure occurs upon reaching a certain critical threshold which corresponds to the minimally viable population size. f Our incremental modeling framework: three successive stages of model refinement, each introducing an additional layer of biological complexity. Model I considers constant hazard rate (λ) due to background mortality fixed at some age. Model II adds mortality due to somatic mutations to overall mortality. Model III accounts for cell replication, which is supposed to compensate partially for cell loss.

Like any mutation, the somatic ones arise from infidelity in DNA damage repair, replication, or mitosis34. While natural selection strives to preserve DNA integrity, perfecting its maintenance is too costly in resources and time, alternatively spent on transcription and proliferation35, so species tolerate mutation rates that allow their survival36. The disposable soma theory37 implies that germline genetic information must be protected by all means, while somatic DNA is not inherited, so it needs less fidelity34,38 (experimentally confirmed by higher somatic versus germline mutation rates19,20,39). Accumulated somatic damage causes cell function to decline, leading to diseases and organismal aging40. Specific recurrent nucleotide alterations (mutational signatures) may confer cell survival benefits, causing clones harboring these signatures to proliferate faster and spread mutations widely in mitotically active tissues through clonal expansion events (e.g., clonal hematopoiesis), further accelerating aging29,32,41.

While there are approaches that might eventually mitigate somatic mutations (including elimination of high-burden cells34, whole-genome editing34,42, or harnessing revertant mosaicism43,44), they remain extremely far from clinical application. Consequently, hypothetical anti-aging therapies might only reverse conditionally reversible hallmarks while being fundamentally unable to affect entropic phenomena like somatic mutations. To date, no attempts have been made to estimate how long a future human, having undergone a treatment against reversible aging, could possibly live (Fig. 1c). What would be the maximum effect of radical life extension45,46 excluding extravagant interventions like massive cell replacement47,48 or transferring consciousness onto a silicon chip49,50,51? In other words, what is all this strife for?

This work aims to estimate how long humans would live if we cured most aging processes except somatic mutations (Fig. 1c). Given accumulated sequencing data, relative simplicity of mutation mechanisms52, and, on the other hand, virtual irreversibility of somatic mutation burden, focusing on this phenomenon provides a plausible assessment of human lifespan’s upper bound.

Previous models eval‎uating the relationship between somatic mutation accumulation and lifespan found that mutation rates are inversely proportional to mammalian species lifespan53,54—a pattern envisioned in the Peto’s paradox, which implies that larger, longer-lived species such as whales have evolved lower per-cell mutation rates and more robust genome maintenance mechanisms to compensate for their greater cell numbers and longevity55,56. While intriguing, this relashionship oversimplifies the dynamics of mutation-driven aging without providing frameworks for theoretical investigations. Other attempts to construct dynamical systems for assessing cell population dynamics focus mainly on stem cells57 and are better suited for studying cancer development58,59 rather than normal aging.

Estimating maximum possible lifespan while theories of aging remain immature is challenging and requires its own novel methodology. Comprehensive estimation would require complex models of somatic mutation dynamics across all tissues and organs plus a host of public data on mutation rates and stem cell differentiation across ages, tissues, and individuals. However, guided by reliability theory, we instead represent the human body as a system of conditionally-independent components60,61 (Fig. 1d), configured in parallel (redundantly) or in series (non-redundantly). This model predicts that if any unique component (e.g., brain, heart, liver) or all instances of parallel components (both lungs or kidneys) reach critical malfunction—due to cell death, in the simplest scenario (Fig. 1e)—the whole system dies, as do the real organisms unless organ function is restored artificially. We can thus interrogate each essential organ individually, focusing on those which are well-studied in terms of somatic mutation burden.

Another technique routinely used in mathematical and physical reasoning involves eval‎uating lower and upper bounds of given parameters. For terminally differentiated, non-dividing cells like neurons, stochastic death from deleterious mutations would provide an upper lifespan bound, considering brain failure equivalent to organism death. When modeling other tissues, incorporating cell division and differentiation would markedly improve the estimates.

In our work, we propose an incremental approach to modeling aging, which involves incorporating aging factors progressively to increase model complexity when necessary (Fig. 1f). We developed dynamical system models describing particular critical organs that deteriorates under somatic mutation burden, with model complexity being adjusted based on tissue replication capabilities. These models apply reliability theory to biological systems and employ upper/lower bound estimates for variables lacking solid experimental evidence. By refining and generalizing the models—first per organ, then within interconnected critical blocks (Fig. 1d)—we enhance their capacity to reflect biological reality and estimate lifespans. Applying this framework, we leverage data from brain, heart, liver, and lungs, offering an in-depth view of organ-specific cell population dynamics under somatic mutation pressure, allowing us to estimate the limits to human lifespan when aging is driven purely by somatic mutagenesis. We thus present a robust framework for eval‎uating the impact of somatic mutations on aging (Fig. 1c–e) and propose a novel approach to modeling aging, which could potentially enable ranking of fundamental aging mechanisms by their contribution to organismal aging, and determining which should be prioritized in our efforts to counteract them.

Results

Designing an incremental approach for modeling aging under somatic mutations

Before modeling aging, we must define it precisely. Despite numerous theories, no consensus exists62,63, except that aging ultimately manifests as a demographic phenomenon through population-level survival curves. Here, we adopt a working definition of aging as the age-dependent increase in mortality risk (also known as hazard rate). Modeling large-scale population aging while accounting for all mechanistic biological factors presents a formidable challenge. Existing models capture only specific molecular aspects (e.g., senescent cell accumulation64) or omit detailed mechanisms in favor of phenomenological perspectives, treating organisms as systems that accumulate abstract damage17,65,66.

To combine biological factors with population-level survival, we designed an incremental strategy. At each step, we introduce an elementary survival-contributing factor, progressively refining the model. Our baseline (Model I) describes a non-aging population where all individuals have constant mortality rates. Next, we view every individual as a system of critical elements (organs), a failure of each causing instant death. Organs fail when a certain percentage of their cells die. We focus solely on cell death from somatic mutations, assuming all other death causes are incorporated into background mortality from Model I and do not increase with age. Model II considers background mortality plus cell loss from somatic mutagenesis. Model III adds replication for actively proliferating tissues, in three variants. Model IIIA includes somatic cells with replication limited by proliferative potential (Hayflick limit67,68). Model IIIB includes somatic cells plus adult stem cells which replenish somatic pools via differentiation, possess unlimited proliferative potential, but die from mutations as well. Model IIIC describes stem-like cells with more prolonged but still limited proliferation that leads to both self-replenishment and differentiation. Accounting for variability in maximum organ capacity and cell death/replication rates yields varying cell population trajectories per organism, incorporating cell-level uncertainty when calculating population-level survival curves.

Finally, we integrate individual aging organs into a single system with the reliability theory framework, yielding combined survival trajectories that unite background and organ-specific mortality due to somatic mutations.

Model I: Baseline model of a non-aging population yields an enormous theoretical maximum lifespan

The first step of our framework establishes a baseline model to describe mortality dynamics in the complete absence of aging. Technically, a non-aging population is assumed to have constant mortality rate throughout lifespan, allowing to calculate its survival for Model I (Fig. 2a) as:

(1)

where λ is constant mortality rate from demographic data fixed at a specific age, t is time elapsed since birth, and S(t) is population survival function. The Gompertz law69 of exponentially increasing mortality degenerates into constant mortality. But how to determine this constant λ?

Fig. 2: Results of the incremental modeling framework for predicting maximum human lifespan under somatic mutation accumulation.

Full size image

a Model I: baseline constant-risk model which accounts only for background (BG) mortality λBG. b Hazard rate (HR) estimated for different ages in the Swiss population. According to Model I, maximum and median lifespans are defined as ages at which the survival probability drops to 10−5 (or 10−3%) and to 0.5 (or 50%), respectively. Both lifespans are calculated using HR frozen at every 10 years, except for the point of lowest mortality (HR freeze at 8 years instead of 10). The lifespans calculated for the 30-years mortality are highlighted with black frames and larger font size. c Simulated survival curves for Model I across different ages of HR freeze. Age numbers in the legend correspond to the numbers along the curves. The curve calculated for the 30-years mortality is highlighted with a black frame around its age number. d Model II: adds organ-specific somatic mutation accumulation which causes cell death at rate μ as an additional source of mortality risk λmut. e Cell population dynamics simulated according to Model II for brain (left) and heart (right). Organ failure occurs when cell population falls below a critical threshold Xcrit. f Survival curves simulated according to Model II for brain (left) and heart (right) under three scenarios: BG-only, organ-only, and combined (mortality from organ-specific accumulation of somatic mutations and BG mortality). g Model IIIA: further includes cell replication via function ρ(P), where P denotes the average replication potential of the whole cell population. h Somatic cell population dynamics simulated under Model IIIA for hepatocytes. i Survival curves simulated for Model IIIA. j Model IIIB: extends IIIA by adding adult stem cells (including stem cell proliferation ρy with unlimited potential) which replenish the somatic pool at rate σ and are also prone to death from somatic mutations at rate μy. k Somatic hepatocyte population dynamics simulated according to Model IIIB in the presence of liver progenitor cells (LPCs); the death threshold is not reached within 100,000 years of simulation. l Survival curves simulated for Model IIIB. m Model IIIC: stem-like cell type which exhibits mutation-driven mortality (at rate μb), differentiation, and limited proliferation (at rate ρb = f(P)). n Lung cell population dynamics simulated under Model IIIC for bronchial basal cells. o Survival curves simulated for Model IIIC. In e, h, k, and n, shaded regions represent percentile bands (from the 2.5th to the 97.5th in 5% increments) for the simulated X(t) trajectories; bright solid line corresponds to the median trajectory; red dotted line is the critical threshold of cell population required for organ survival (Xcrit); the time at which the median cell population trajectory crosses this critical threshold is annotated as tcrit. In f, i, l, and o, horizontal dotted gray lines correspond to survival probabilities of 0.5 and 10−5; median and maximum lifespans are annotated for organ-only and combined survival curves (and for the BG curve in f) as tmed and , respectively. In i, l, and o, survival curves for models calculated using BG and combined mortalities overlap each other and are represented next to each other for visualization purposes only.

From empirical life tables (http://www.lifetable.de), we can eval‎uate mortality rates at different ages and solve equation (1) for t (see “Methods”). We can thus obtain the lifespan of a population with “frozen” mortality λ as the value of t at any S(t) (percentage remaining from the initial population), assuming its initial size as N = 8 billion individuals—approximately current human population70. We define median survival as time until half the initial population survived and maximum survival as time until 10−3% (1 in 100,000 individuals from the starting cohort) survived, yielding:

(2)

Naturally, mortality rates at different ages and survival curves vary substantially by population. We used life tables to estimate median and maximum lifespans across 16 populations in our dataset (Supplementary Fig. 1) and focused on Swiss data (Fig. 2b, c) for downstream modeling, as Switzerland exhibits the lowest extrinsic hazard rate at age 30 among countries with the longest post-1998 observation records.

Importantly, non-aging populations are not immortal: death causes continue acting on individuals, inevitably leading to eventual extinction, which is reflected by survival curves for mortality rates fixed at every ten years: all of them eventually fall to nearly zero (Fig. 2c).

All 16 studied populations display J-shaped curves of mortality rate, with their minima lying between 5–15 years depending on the country, followed by steep exponential rise appearing as linear growth in the logarithmic scale, with some populations (Switzerland, Singapore, Denmark, Israel, US, etc.) exhibiting intermediate mortality plateaus between ages 20–30 (Supplementary Fig. 1).

We calculated median and maximum lifespans for several age points of mortality “freeze” (Fig. 2b, Supplementary Fig. 1). For the Swiss population, mortality fixed at age 20 yields median lifespan of 1992 years and maximum of 33,088 years. Though extraordinary, these numbers follow directly from our framework. Conversely, fixing mortality at the 110-year-old supercentenarian levels yields median remaining lifespan of only 1 year—consistent with intuition. The maximum remaining lifespan is 12 years, reflecting probabilistic survival of at least one individual among 100,000 under such elevated mortality risks, which aligns with the reflections of the late Mikhail Blagosklonny who noted that “A mere application of standard medical care to centenarians, as rigorously as to younger adults, would probably extend lifespan beyond 122, even without the need of a scientific breakthrough”71.

The childhood minimum mortality point (8 years in the Swiss dataset) corresponds to highest predicted lifespans, both median (10,395 years) and maximum (172,665 years)—a striking illustration of how age-dependent differences in mortality profoundly affect projected lifespan in the non-aging model.

The choice of baseline mortality risk links to the question: “At what age does aging begin?” Perspectives on this issue vary widely62,63—from aging starting during embryogenesis72 to the arguments that it begins only when body growth ceases, consistent with disposable soma theory37,73. Mortality risk curves (Fig. 2b, Supplementary Fig. 1) reveal that the exponential mortality increase typically emerges around age 30—a threshold coinciding with the commonly recognized onset of middle age. We adopt the age 30 mortality rate as our baseline, which corresponds to median and maximum remaining lifespans of 1759 and 29,221 years. This choice reflects the mortality risk of a typical early adult characterized by good health, low incidence of age-related pathologies, and full societal engagement. At this age, individuals are generally active, employed, and routinely exposed to diverse extrinsic hazards including car accidents, drug abuse, occupational hazards, infectious diseases, etc., thus embodying composite baseline mortality risks without the influence of aging.

Consequently, all further analyses assume that mortality rate is fixed at the levels of a modern 30-year-old person and referred to as background mortality (λ30 y.o. = λBG). Building upon this assumption, we incrementally incorporate additional factors related to somatic mutagenesis and cell proliferation that contribute to organismal survival, progressively enhancing model realism and explanatory capacity.

Model II. Somatic mutations in post-mitotic tissues drive substantial mortality increase but cannot fully account for observed aging

Building upon the baseline non-aging model, we next introduced somatic mutagenesis in non-dividing cells, focusing on neurons and cardiomyocytes as representative post-mitotic cell types (Fig. 2d-f). Cell loss dynamics driven by deleterious mutations can be described by an ordinary differential equation (Model II):

(3)

where X is the current number of cells at risk in an organ, dX is the number of cells lost over time dt, and μ = μ0 ⋅ plethal is death probability per cell, which equals to the product of: μ0—somatic mutation accumulation rate measured as number of mutations per cell per year, and plethal—the probability of a single mutation to be lethal for a cell. Both μ0 and plethal can be estimated separately for every mutation type—we focused on single nucleotide variations (SNVs) and insertions-deletions (indels)—and then summed up for a combined estimate of μ:

To obtain μ0 for each cell type, tissue-specific mutation accumulation rates were estimated via mixed-effects linear modeling of single-cell somatic mutation burdens against donor age74,75,76,77, fitted separately for SNVs and indels (Supplementary Fig. 2, Table 4). To derive plethal—the probability that a given somatic mutation kills the cell—we first defined tissue-specific panels of essential genes from in vivo CRISPR knockout/activation screens, and then scored each observed mutation by how likely it is to disable such a gene, combining: whether it falls within an essential gene’s functional regions (gene body, splice sites, promoters, and enhancers), how damaging it is predicted to be based on a median estimate of multiple deleteriousness predictors (provided in Ensembl VEP and FAVOR databases78,79), and how likely the gene is to be haploinsufficient. Summing these scores across all observed mutations and dividing by their number yields the fraction that are lethal, allowing to compute the combined annual lethal mutation rate μ for each cell type (Table 1), yielding μ ≈ 2.44 ⋅ 10−3 for cortical neurons and μ ≈ 2.61 ⋅ 10−3 for cardiomyocytes. These values were sanity-checked against a lower bound based on the known CDS size and loss-of-function fraction of a single essential gene (POLR2A), and an upper bound derived from the maximum possible neuron loss rate across aging, statistically compatible with published cortical neuron count data80. We assumed that organ failure occurs once the number of functional cells falls below a critical fraction determined from empirical physiological thresholds (see “Methods”).

Table 1 Estimated mutation accumulation rates (μ0), per-mutation lethality probabilities (plethal), and resulting annual cell death rates (μ) for each analyzed cell type

Full size table

Equation (3) for Model II corresponds to an exponential decay solution, X(t) = Ke−μt, with decay rate governed by the μ coefficient, and starting from organ cell capacity .

To translate cell-level dynamics into organism-level survival, we determined when each organ reaches its critical failure threshold based on the exponential cell loss from equation (3). Accounting for variability in both decay rate μ and initial organ capacity K, we computed the times to reach organ failure (Fig. 2e). We then combined this somatic mutation-driven mortality with baseline mortality from Model I by summing their respective hazard rates λmut (organ-specific) and λBG, assuming independence between mortality sources that is equivalent to multiplying their survival functions (see “Methods”):

(4)

Importantly, the survival curve in this Model (and in all others in this work) is not obtained by comparing a mean cell-count trajectory to the failure threshold. Instead, each individual in the synthetic cohort receives an analytically computed failure time drawn from the joint distribution of organ capacity and mutation rate (K, μ), so the resulting population survival curve correctly captures individual-level variability, including early-failing individuals with smaller organs or higher mutation rates (see “Methods”).

Combining baseline mortality (fixed at age 30) with somatic mutation-driven cell death yields a population survival curve with median lifespan of 194 years and maximum of 557 years for neuron-driven aging, and 208 and 868 years, respectively, for cardiomyocyte-driven aging (Fig. 2f). Compared to the non-aging baseline (median: 1759 years, maximum: 29,221 years), this represents a substantial reduction driven by a single elementary aging mechanism. However, the resulting maximum lifespans remain markedly longer than the observed maximum of 122 years, while the resulting median lifespans are approximately three times as long as median life expectancies in developed countries, indicating that somatic mutations in post-mitotic tissues, while contributing substantially, cannot alone explain observed human mortality patterns.

Model III: Cellular replication in proliferating tissues provides strong protection against somatic mutation-driven aging

The previous model applies only to post-mitotic tissues that cannot replace lost cells. However, most human tissues comprise actively proliferating cells capable of compensating for cell loss through replication. To determine whether cellular replication can mitigate somatic mutation-driven aging, we extended our model to include both dividing somatic cells with finite replicative potential and a pool of stem cells that self-renew and differentiate (see “Methods”):

(5)

where μ and μy represent cell death rates for somatic and stem populations, while ρ and ρy represent their respective replication rates.

Somatic cells possess limited proliferative potential P (Hayflick limit), which depletes with each division at rate ψ, but can be replenished at rate ψy through stem cell differentiation, which adds cells to the somatic pool at rate σ. Importantly, this proliferative ceiling is a constitutive physiological property of somatic cells present from embryonic development—its age-associated progressive worsening (driven by immune decline, mitochondrial dysfunction, and other processes assumed to be eliminated in our framework), rather than its mere existence, is what constitutes an aging hallmark (see “Discussion”).

All mentioned rates, except the death rates, are functions that depend on current somatic/stem cell populations (see “Methods”).

Although the resulting system of differential equations lacks a closed-form analytical solution, it can be solved numerically. Using empirical data, we estimated the required parameters (see “Methods”; Supplementary Figs. 2, 7, 8, 9, 11, and 12) for two representative proliferating tissues: liver (hepatocytes and liver progenitor cells, LPCs) and respiratory epithelium (bronchial basal cells, BBCs). We implemented three variants of this model.

Model IIIA: hepatocytes, finite replicative capacity, no stem cells

Simulation results for Model IIIA (Y = 0; Fig. 2g-i) reveal that even without the LPC support, hepatocyte replicative capacity provides robust protection against somatic mutation-driven cell loss. Median time to reach the critical threshold extends to 37,664 years, far exceeding failure times in post-mitotic organs (198 years for neurons, 212 years for cardiomyocytes). Consequently, liver failure contributes negligibly to population-level mortality: the median and maximum lifespans remain virtually unchanged at 1755 and 29,147 years (Fig. 2h, i), matching the non-aging baseline.

Model IIIB: hepatocytes + liver progenitor cells

Model IIIB (calculated using full equation (5); Fig. 2j-l) incorporates LPC-mediated hepatocyte replenishment and yields even more striking results. None of the simulated trajectories reached organ failure within 100,000 years, yielding 100% survival (Fig. 2k, l). These findings indicate that when supported by stem cells, liver exhibits exceptional resilience to somatic mutation-driven damage and does not constitute a lifespan-limiting bottleneck. Hence, liver mortality can be safely neglected within our modeling framework.

Model IIIC: respiratory basal cells

Unlike hepatocytes in Model IIIB, which divide symmetrically and are continuously replenished by a separate upstream stem cell pool (LPCs), respiratory basal epithelial cells play the stem-like role themselves: they undergo the full spectrum of stem cell division fates—symmetric self-renewal, symmetric differentiation, and asymmetric division—with no external reservoir to restore their numbers or proliferative potential (which drops with each division, unlike that of LPCs), so that each division of any type irreversibly draws down their remaining replicative capacity.

Model IIIC for these basal cells (see “Methods”; Fig. 2m-o) reveals somewhat distinct aging patterns compared to Models IIIA and IIIB. Despite finite replicative capacity and continuous loss through differentiation, basal cells initially maintain stable population dynamics: cell population trajectories remain nearly horizontal for thousands of years before declining linearly as cells exhaust their proliferative potential and succumb to accumulated mutations. Median cell population trajectory crosses the critical threshold at 4423 years (Fig. 2n). This delayed decline manifests in the organism-level survival curve as prolonged stability at 100% survival, followed by a rapid drop, yielding median lifespan of 4359 years and maximum of 7617 years (Fig. 2o). When combined with baseline mortality, the resulting population survival curve closely resembles the non-aging baseline throughout most of the lifespan, with median survival remaining at 1755 years. However, unlike Models IIIA and IIIB, maximum lifespan becomes reduced from the baseline to 7132 years. These results indicate that while respiratory epithelium does not limit median lifespan, it may eventually impose an upper bound on maximum survival, suggesting that even tissues with substantial regenerative capacity can become lifespan-limiting factors at extreme ages.

Collectively, these findings demonstrate fundamental asymmetry in how somatic mutations affect aging across tissue types. While somatic mutations drive substantial mortality increases in post-mitotic tissues, incorporating cellular replication effectively neutralizes this effect in studied proliferating tissues (Supplementary Figs. 4 and 5). Across all three replication-based models (IIIA, IIIB, and IIIC), organ failure times extend far beyond the lifespans observed in normal aging and projected for non-aging, rendering the individual contribution of somatic mutations negligible in actively proliferating tissues.

Predicted lifespans are robust to uncertainty in lethal mutation probability

To assess the robustness of our predictions to uncertainty in plethal, which is by far the most difficult parameter to estimate empirically, we conducted a sensitivity analysis, in which  and  were independently sampled from log-uniform distributions spanning the full biologically plausible range for each organ (n = 300 parameter sets; see “Methods”). Across all post-mitotic organs, predicted lifespan decreases monotonically with increasing plethal on a log-log scale, with no threshold effects or discontinuities (Supplementary Fig. 13). Our empirically derived default values sit near the upper end of the sampled range, meaning our main estimates represent conservative (shorter) predictions within the plausible space. The uncertainty bands on survival curves span roughly one to two orders of magnitude in time, yet the qualitative conclusion (post-mitotic organs fail well within the simulation window, proliferating tissues persist) holds across the entire parameter space (Supplementary Fig. 14). Most strikingly, the liver + LPC model never reaches either the median or maximum survival threshold under any sampled parameter combination, confirm‎ing that the resilience of proliferating tissues to somatic mutation pressure is not sensitive to plethal.

Multi-organ model reveals theoretical lifespan limits under pure somatic mutation-driven aging

We now integrate these organ-specific models into a multi-organ framework inspired by the reliability theory60. An organism is represented as a series system of organs, surviving only if all critical organs remain functional, i.e., failure of any single organ causes death (Fig. 3a). This series configuration also includes an additional element representing background mortality λBG, as defined in Model I. The series arrangement reflects biological independence between the modeled cell populations across organs; a parallel sub-structure would only be appropriate for anatomically separable organ units with independent cell pools capable of mutually compensating for each other’s failure. For instance, in contrast to lung endothelial cells, which are confined to their respective lobes, respiratory basal cells that we use in our modeling form a single continuous epithelial pool that migrates freely along the airways81, making a lobe-wise parallel decomposition of lungs biologically unmotivated in our multi-organ model. Assuming independence among aging processes in individual organs, the organismal survival curve can be expressed as the product of each organ’s survival functions and baseline survival (Fig. 3a).

Fig. 3: Multi-organ integration reveals theoretical lifespan limits under somatic mutation-driven aging.

Full size image

a The organism is modeled as a series reliability system in which four critical organs (brain neurons, heart cardiomyocytes, liver hepatocytes, and respiratory basal epithelium) age independently through somatic mutation-driven cell loss, alongside baseline mortality (λBG) fixed at age 30 levels. Failure of any single organ causes organismal death. Under the independence assumption, overall survival is the product of individual organ survival functions. b Population survival curves under different inter-organ dependency assumptions. The independence model is represented by solid black line. Fréchet bounds establish theoretical limits shown as dashed gray lines: from perfect positive dependence (upper bound), assuming that all organs age in synchrony, to perfect negative dependence (lower bound), assuming anti-correlated aging. Fréchet bounds establish theoretical limits shown as dashed gray lines: from perfect positive dependence (upper bound), assuming that all organs age in synchrony, to perfect negative dependence (lower bound), assuming anti-correlated aging. These extreme scenarios are biologically unlikely but serve as mathematical guarantees—the true survival curve must lie within the shaded region regardless of the actual inter-organ dependency structure, which remains empirically undetermined. The shaded region with diagonal hatches, therefore, represents the full range of mathematically possible outcomes, not a range of biologically probable ones. Pink line indicates normal aging (for the Swiss population). Gray horizontal dotted lines correspond to survival probabilities of 0.5 and 10−5. Median and maximum lifespans are annotated as tmed and , respectively, for the multi-organ and normal aging models. Survival curves for models IIIB and IIIC overlap each other and are represented next to each other for visualization purposes only.

The resulting composite survival curve describes a hypothetical human population aging solely through somatic mutations in four critical cell types (neurons, cardiomyocytes, hepatocytes, and respiratory basal epithelial cells), with background mortality fixed at 30-year-old levels. According to this model (Fig. 3b), estimated median lifespan is 156 years, corresponding to 197% of the empirically observed median human lifespan (in the Swiss cohort) of 79 years. Maximum lifespan reaches 470 years, or 385% of the empirically observed maximum of 122 years.

The independence assumption represents a strong simplification of biological reality, because organs do not age in isolation. To quantify uncertainty arising from potential inter-organ dependencies, we applied Fréchet inequalities82 (see “Methods”). These inequalities establish upper and lower bounds for the joint survival probability. The upper bound assumes perfect positive dependence—all organs age in perfect synchrony, as if their health states were perfectly correlated. In this scenario, organismal survival is limited by whichever organ fails first. The lower bound assumes perfect negative dependence—one organ’s integrity maintained at another’s expense. While neither extreme is biologically realistic, they provide conservative boundaries that quantify maximal possible deviations from independence.

Under these bounds, median lifespan ranges from 146 years (perfect negative dependence) to 194 years (perfect positive dependence) (Fig. 3b). Maximum lifespan shows greater sensitivity to inter-organ dependency, spanning 210–557 years across the same range of bounding assumptions. This broad interval for maximum lifespan reflects compounding effects of dependency assumptions when extrapolating to extreme population tails.

Discussion

Motivation for this study stems from a deceptively simple yet profoundly consequential question: What if we were to eliminate all reversible mechanisms of aging, leaving only irreversible, entropy-driven processes such as somatic mutations? How long could a human then live? To address this, we developed an incremental modeling approach—a methodological strategy common in mathematics and physics but still rare in geroscience—that progressively incorporates biological aging factors upon a baseline non-aging model. We employed reliability theory60, which conceptualizes organs as critical system components with finite failure tolerance, and conducted extensive literature review to estimate key biological parameters with rigorous validation against empirical bounds.

Thus, this work introduces a framework that explicitly links mechanistic drivers of aging to population-level manifestations, with potential to rank canonical hallmarks by their relative impact through shifts in survival curves. Importantly, our estimates assume: (1) no organ/tissue transplantation; (2) no interventions reducing mutation accumulation; and (3) baseline all-cause mortality fixed at contemporary 30-year-old levels.

Our results suggest that even if all hallmarks of aging were eliminated except for somatic mutations, median human lifespan would reach only 146–194 years—roughly twice the current 79 years. This approximate 2-fold gap in median survival indicates that other hallmarks of aging10 likely contribute at least as much to limiting human lifespan as somatic mutations do. Maximum lifespan would range from 210 to 557 years depending on inter-organ dependency structure (470 years under independence). While these estimates substantially exceed current longevity, they fall well short of the optimistic projections of achieving “longevity escape velocity”83.

A key insight is that post-mitotic cell types (neurons, cardiomyocytes) act as critical longevity bottlenecks under somatic mutation pressure, with maximum failure times of 557 years (brain) and 868 years (heart)—far shorter than the 29,147-year non-aging baseline. This reflects inexorable mutation accumulation in irreplaceable cells, driving median lifespan from 1755 to 194 years for brain aging alone. In stark contrast, proliferating tissues show remarkable resilience. For liver with stem cell support (Model IIIB), no simulated trajectories reached failure within 100,000 years. Even without stem cells (Model IIIA), hepatocytes maintained functionality beyond 37,000 years, as dead cells are continuously replaced by undamaged progeny. The exception is airway basal cells (Model IIIC), which fail at ~ 4359 years median time due to proliferative exhaustion, suggesting that even replication-competent tissues may still constrain extreme longevity.

Interestingly, applying Model II to both hepatocytes and airway basal cells yields median lifespans approximately equal to BG mortality even without any replication (Supplementary Fig. 3).

This asymmetry between tissue types supports the notion that aging is not uniform across the organism and provides rational basis for prioritizing therapeutic interventions. Somatic mutations alone are unlikely to bottleneck longevity of proliferating tissues like liver and respiratory epithelium, as regenerative capacity neutralizes their impact (provided protective mechanisms like immune surveillance remain fully functional). Conversely, post-mitotic organs such as brain and heart emerge as critical lifespan bottlenecks. Especially, neuronal irreplaceability imposes severe constraints even if other organs are rejuvenated or replaced. In this context, translational cell therapy aimed at neuronal replacement48 may represent our only viable path to 150+ year lifespans, even assuming all other hallmarks are mitigated.

A systemic rejuvenation approach that restores tissue renewal capacity to youthful levels would, in our framework, correspond to proliferating tissues operating under the most favorable Model III parameters—and our results show that under such conditions, somatic mutations contribute negligibly to mortality. The residual bottleneck would then be entirely determined by post-mitotic organs, where mutation accumulation is irreversible and not directly dependent of systemic signaling. Slowing this accumulation requires orthogonal interventions that target mutagenesis itself rather than the downstream signaling environment.

Notably, both neurons and, especially, cardiomyocytes feature elevated physiological ROS production, which is also reflected in their oxidative damage-driven mutational signatures76,84,85. In the idealized non-aging scenario, physiological ROS production would remain, while its pathological age-progressive amplification would not. Conversely, in real aging, pathological processes such as Alzheimer’s disease85, multiple sclerosis76, and ischemic heart disease77 further accelerate somatic mutation accumulation beyond the baseline rate, increasing cellular susceptibility and organ dysfunction. Our model is not intended to capture this feedback loop, and including it would pull post-mitotic organ failure earlier than our estimates suggest.

Among the limitations, our model cannot be falsified in the conventional sense, as no organism ages exclusively through somatic mutations. However, it permits a critical sanity check: predicted lifespans substantially exceed the observed ones (median: 79 years; maximum: 122 years) by ~2-fold for median and 2–5-fold for maximum lifespan. This outcome validates our modeling approach: had the predictions fallen below observed values, it would indicate either fundamental flaws in model structure or systematic biases in parameter estimation (e.g., overestimated mutation accumulation rates). Instead, the model suggests that somatic mutations in the four examined organs markedly increase mortality from the non-aging levels, but still cannot fully explain observed human mortality.

Our framework incorporates only four organs and assumes independence—a simplification of real physiological interdependencies. The Fréchet bounds (median: 146–194 years; maximum: 210–557 years) provide mathematical limits on how inter-organ dependencies could shift our predictions, but the true dependency structure remains unknown.

Additionally, our model considers only lethal mutations that cause immediate cell death, neglecting potential pathways through which mutations affect tissue function. First, sublethal mutations may impair cellular function without triggering death, gradually degrading tissue maintenance capacity. Second, we assume malignant transformation is effectively cleared by immune surveillance, thereby excluding cancer as an age-dependent failure mode. We note, however, that age-independent cancer risk, which reflects the baseline probability of immune surveillance failure present even in young healthy individuals, is already embedded in the background all-cause mortality rate, which we freeze at the level of a 30-year-old throughout our models. Third, and perhaps most important for proliferating tissues, we do not model clonal expansion—the preferential proliferation of mutant cell lineages that outcompete other cells. Fourth, we do not model non-malignant but functionally suboptimal expanding clones, that is, lineages that proliferate actively while being metabolically or functionally compromised, such as the mtDNA-deficient intestinal crypt stem cells86, where mitochondrial clones expand via crypt fission despite respiratory chain dysfunction. Such clones would progressively erode tissue function in proliferating tissues over extended timescales through a mechanism distinct from both cell death and malignant transformation.

The relevance of these limitations varies substantially across the cell types we model. For neurons and cardiomyocytes, clonal expansion is mechanistically impossible: both are post-mitotic, so mutations are strictly cell-autonomous and no clonal dynamics are possible. For hepatocytes in a healthy liver, clonal expansion is spatially constrained by lobule architecture: clonal patches grow slowly and then quiesce, bounded by lobule boundaries, with pronounced expansion requiring chronic injury or cirrhosis87,88. For respiratory basal cells in non-smokers, clonal dynamics follow neutral drift among equipotent progenitors89, with non-neutral selection strongly associated with extrinsic carcinogenic insults such as tobacco smoke90,91 that are outside our model’s scope. In stark contrast to that is clonal hematopoiesis of indeterminate potential (CHIP): hematopoietic stem cells are both a spatially unrestricted, high-competition niche92 and a source of circulating immune cells that drive systemic inflammation across organs93—a regime none of our modeled cell types operate in. Where clonal selection is relevant for proliferating tissues, its incorporation would only accelerate functional decline relative to our predicted lifespans, consistent with the upper bound character of our estimates.

Reliability theory in our framework captures organ-level failure dynamics but does not account for systemic signaling and inter-organ crosstalk that characterize the body as an integrated organism. Crucially, including each of these omitted factors and other hallmarks of aging would only pull predicted lifespan further downward—our framework is therefore intended to produce an upper bound on longevity constrained by somatic mutagenesis, not a projection of achievable lifespan.

Accordingly, our outputs should be interpreted as estimates grounded in experimentally measured mutation rates, empirically derived cell division kinetics, functionality thresholds, and parsimonious assumptions regarding failure dynamics—analogous to how sample means estimate population parameters.

In perspective, this work establishes an extensible, quantitative foundation for modeling how individual aging mechanisms constrain lifespan. By demonstrating that somatic mutations alone would permit median lifespans of 146–194 years—substantially above current values yet well below the non-aging baseline—we provide an empirical reference against which other hallmarks can be compared. Ongoing efforts profiling somatic mutations across tissues94 and species will enable incorporating additional organs and cross-species validation, which may reveal further lifespan-limiting bottlenecks. Parameter estimation would significantly benefit from in vivo experiments on replication rates, mutation lethality, and proliferation limits. Notably, cross-species extension of our framework would require species-specific parameterization—at minimum μ0 and plethal, and other additional parameters for the proliferating tissues—and Peto’s paradox, the long-standing puzzle of why whales with ~107 more cells and ~102 longer lives than mice do not show proportionally higher cancer incidence55, is naturally accommodated within it. The paradox dissolves once species-specific genome maintenance is taken into account: bowhead whales, for instance, accumulate mutations far more slowly than other mammals, owing likely in part to ~100-fold elevated expression‎ of CIRBP, a protein enhancing double-strand DNA break repair fidelity56. A whale entering our model with its own empirical μ0 would yield estimated lifespans consistent with its observed longevity.

Our Model III incorporates the Hayflick limit for cell proliferation as a constitutive biological property of somatic cells, one that is present from birth and independent of age-progressive deterioration, and therefore should not be conflated with an aging hallmark. Indeed, being post-mitotic is just as physiological for certain cell types as having a limited but non-zero proliferative potential: neurons, cardiomyocytes, and respiratory ciliated cells (the differentiated progeny of BBCs modeled in Model IIIC) can all be described by the same equations as Model III with remaining proliferative potential (P) set to 0—which immediately reduces to Model II. Moreover, many rapidly proliferating cell types exhaust their replicative limit within days: for example, intestinal transit-amplifying cells divide every ~12 hours and complete only 4–5 cycles before giving rise to post-mitotic daughters95. The appearance of senescent cells per se is also a physiological phenomenon, observed to support normal mammalian embryogenesis, tissue regeneration, and wound healing96,97,98,99. What would constitute an aging hallmark is instead the progressive exhaustion of the replicative limit (H) of stem cells contributing to stem niche exhaustion100,101 or the marked accumulation of senescent cells. Both appear to be driven by other processes: immune decline impairs senescent cell clearance100,102,103,104, and mitochondrial dysfunction accelerates telomere attrition105,106,107—all of which are assumed to be eliminated in our framework. In the idealized setting of Model III, senescent cells are therefore promptly cleared and stem cell niches remain intact, so the Hayflick limit acts purely as a fixed biological ceiling on replicative capacity in non-stem cells rather than a progressive aging driver.

Integrating all hallmarks (mitochondrial dysfunction, telomere attrition, proteostasis loss, epigenetic alterations, etc.) requires substantially greater complexity, encompassing extracellular matrix and intercellular networks108. However, our work provides a critical first step towards dissecting aging into quantifiable mechanistic components. Ultimately, incorporating other major aging factors could pave the way to a comprehensive, mechanistic theory of aging. Achieving it would require coordinated efforts across multiple research groups—a challenging but increasingly attainable prospect in the near future.

Methods

Model I: constant hazard rate

Survival function S(t) is defined as the probability that an individual survives beyond the age of t years—in other words, it is the probability to die at age τ which is later than t, therefore:

(6)

where F(t) = P(τ ≤ t) is equivalent to the fraction of population that has died by age t.

Hazard function λ(t) is the instantaneous rate of death at age t, conditional on having survived to this age t. In other words, it is the probability that an individual dies within time interval Δt, given that they survived until t, divided by this interval Δt:

(7)

Let t1 < t2 < ⋯ < tk be the times of observed death events in a certain population. Empirical survival probability for this population according to the Kaplan-Meier estimator109 is:

where di is the number of deaths at each ti and ni is the number of individuals at risk (still alive) just before ti.

With the help of a common formula from survival analysis110, we can obtain hazard rate as:

(8)

where  is equivalent to probability density function (PDF) of survival , and is essentially the probability of dying between time points t and t + Δt. The derivative  was computed from  using numpy.gradient from the numpy111 Python package.

If we freeze hazard rate λ(t) at its value for a certain age t and denote it as λ, then we can derive the survival function S(t) for Model I from equations (6), (7), and (8) as:

(9)

Within this constant mortality framework, we define median lifespan tmedian as the time at which half the initial population remains, and maximum lifespan tmaximum as the time at which the survival probability falls to 10−5 (or 10−3%)—corresponding to 1 in 100,000 individuals, or approximately 10 surviving trajectories in our Monte Carlo simulations of 106 individuals. These times are found by solving S(t) = 1/2 and S(t) = 10−5 for t, yielding:

(10)

Survival function for the Switzerland population derived from the empirical life-tables (see below) was utilized as background (BG) survival function for all further models, hereafter referred to as SBG(t). Switzerland was selected as the reference population because, among all countries in our dataset with the longest post-1998 observation records, it yields the highest non-aging median lifespan at age 30—reflecting the lowest extrinsic hazard rate among well-characterized developed-world populations (see Table 2 and Supplementary Fig. 1). The hazard rate at age 30 calculated for this population is accordingly referred to as the background mortality rate λBG.

Table 2 Populations included in the analysis of background mortality, with the ranges of calendar years for which data were available in the Human Life-Table Database after filtering by year 1998 onwards (populations ordered alphabetically)

Full size table

Empirical hazard rates calculation

To calculate hazard rates for each age point, we estimated empirical survival probabilities using life-table data obtained from the Human Life-Table Database (HLD; https://www.lifetable.de). Our analysis includes 16 geographically diverse populations (Australia, Brazil, Canada, Chile, Denmark, Hong Kong, India, Israel, Japan, Malaysia, South Korea, Russia, Singapore, Spain, Switzerland, and the United States).

Raw life table data were obtained from the HLD contained files with age-specific central death rates m(x) and person-years lived L(x) across multiple calendar years, for each population. Only records from 1998 onward were retained, reflecting contemporary mortality conditions and excluding historical data that may not be representative of the present-day hazard levels (Table 2). To construct a single representative hazard profile per population, we first averaged the central death rate at each integer age a across all retained calendar years, weighting each year’s estimate by the corresponding number of person-years lived L(x):

yielding approximately 110 age-specific hazard estimates per population. Isolated single-age anomalies were suppressed by applying a centered 3-point rolling median over consecutive ages (function Series.rolling from the pandas Python package). The resulting cleaned discrete profile was then log-transformed and interpolated onto a uniform grid of 2000 equally-spaced ages (from the earliest to the latest observed age) using shape-preserving piecewise cubic Hermite interpolation (PCHIP; function PchipInterpolator from scipy.interpolate), which passes exactly through each data point without introducing additional smoothing. A final boxcar moving average with a window of 4.6 years was applied log-space using a uniform convolution kernel (numpy.convolve) to suppress residual high-frequency noise while preserving biologically meaningful features of the hazard curve, including the childhood mortality decline and the Gompertz inflection at older ages; the bandwidth was chosen empirically as the smallest window eliminating visible spikes in the smoothed profiles. The smoothed log-hazard values were then exponentiated to recover the final continuous hazard curve  used in all subsequent analyses.

Model II: cell death

We begin our incremental construction of the aging model by incorporating the phenomenon of cell death driven by somatic mutations. It can be modeled as a simple exponential decay process that acts independently on each cell of a given organ. When cell population falls below a critical threshold, the entire organ and, therefore, organism dies. This scenario can be captured by combining the baseline mortality introduced in Model I with an additional, somatic mutation-specific mortality component, which we introduce below. The model is based on the following assumptions:

•  
•  
•  
•  
•  

Model II applies to post-mitotic tissues such as neurons and cardiomyocytes, which do not undergo cell division. In the absence of replication, there is no branching process, no clonal expansion, and no selection acting on somatic variants. Each cell independently accumulates mutations at a constant rate, making the linear constant-rate (exponential decay) formulation the correct description rather than an approximation. Heavy-tailed failure dynamics driven by rare early clonal events is a concern relevant to proliferating tissues that does not arise in this setting.

We model cell population X = X(t) change over time as:

(11)

where μ = μ0 ⋅ plethal is per-cell hazard rate equal to the product of μ0—the accumulation rate of somatic mutations measured as number of mutations (SNVs and indels) per cell per year, and plethal—the probability of a single mutation to be lethal for a cell. The calculation of μ0 and plethal is described in greater detail in sections “Estimating accumulation rates μ0 of somatic mutation burden” and “Estimating probability of lethal somatic mutations”. We set the initial number of the cells to X(t = 0) = X0 = K, with K denoting total cell capacity for a given organ.

Individuals are allowed to vary in organ capacity K (some people have larger or smaller organs) and the per-cell hazard rate μ (some accumulate damage faster or slower), both of which are treated as random variables across the population.

The closed-form solution for this model is

(12)

Deriving population survival function for Model II

To derive the population-level survival function from individual cell dynamics, we account for inter-individual variation in organ capacity and mutation rates. While individual cells follow exponential decay (X(t) = Ke−μt), translating this to population mortality requires integrating over heterogeneity in initial conditions, which makes the calculations rather complex.

We start by denoting a critical cell population threshold as Xcrit. When X(t) ≤ Xcrit, the organ can no longer maintain function and the organism dies. Solving X(t = tcrit) = Xcrit gives the time-to-organ-failure for an individual with parameters (K, μ):

(13)

However, instead of asking “when does cell count hit the threshold?”, we can reformulate the problem by asking: “How large must an organ have been at birth to be still alive at age t?”

An individual with organ capacity K and per-cell mutation hazard μ dies when X(t) = Xcrit. Inverting this relationship, we obtain the minimum initial capacity required to survive until t:

(14)

Thus, an individual i (with a fixed pair of K and μ) still survives despite their mutation burden at time t if and only if their initial capacity exceeded this threshold:

(15)

To obtain population survival, we marginalize over the inter-individual distributions (probability density functions, PDFs) of organ capacity ϕ(K) and lethal mutation rate ξ(μ). Population survival is therefore an integral over the joint distribution of K and μ (assuming K and μ are independent across individuals):

(16)

(17)

Because death occurs for K values below Xcriteμt, we can alter the integration limits to include only K values above it:

(18)

For a fixed μ, the inner integral becomes the probability of survival, and it is simply the probability that the initial capacity exceeds the capacity threshold:

(19)

where Φ(…) denotes the cumulative distribution function (CDF) of organ capacity.

To account for the variability in μ—some people’s cells die faster, so they need even larger organs to survive to the same age,—we calculate the expectation (i.e., average) of this probability over the distribution of mutation rates by:

(20)

This expression‎ has a natural interpretation: for each mutation hazard μ, only individuals with sufficiently large initial organ capacity survive to time t; the population survival is the average of this fraction across all μ values.

Assuming independence between somatic mutation-driven population mortality and background (BG) hazards, the full survival function is the following:

(21)

(22)

This formulation captures how population survival emerges from the interplay between initial capacity heterogeneity (variation in K), aging heterogeneity (variation in μ), and the exponential accumulation of cellular damage.

Since analytical solution for S(t) is intractable, we employed Monte Carlo simulation with 10,000 iterations. For each iteration, we sampled {Ki, μi} from log-normal prior distributions ϕ(K) and ξ(μ) (parameterized further below in sections “Estimating accumulation rates μ0 of somatic mutation burden” and “Estimating other model parameters for each organ”), computed the survival function, and averaged across samples. Maximum and median lifespans were obtained by searching for the first entry where S(t) ≤ 10−5 and S ≤ 1/2, respectively.

We note that this procedure integrates over the full joint distribution of (K, μ) across the synthetic population, directly yielding the distribution of individual failure times rather than comparing a mean trajectory to the threshold. Within each individual, the deterministic ODE approximation is justified by the law of large numbers: for organs with K ~ 109–1011 cells, demographic fluctuations around the mean scale as –10−4 and are negligible for the purposes of threshold crossing.

Model IIIA: cell death and replication of somatic cells

We next add the ability to replicate into our framework in order to model proliferating cell types. To do it, we rely on the following assumptions:

•  
•  
•  
•  
•  
•  

To create this model, we employ a system of two equations:

(23)

Modeling replication rate depending on current population and proliferative potential

To model ρX, we choose the simplest law of sigmoidal population growth, which models proliferation as a function of current population fraction x = X/K via a simple logistic gate function g(x) = (1 − X/K) = (1 − x). The replication rate ρ is proportional to this gate ρ ~ g(x) with a coefficient r, that is: . Thus, when X → K (cell count approaches organ capacity), then x → 1 and ρ → 0 (no replication). As X falls down close to zero, so that x → 0, ρ increases up to r. Hence, r represents maximum rate at which the cells of a given type can proliferate (expressed in number of divisions per year). If we consider the X(t) curve starting from X0 = 1 (population growth from 1 cell), then its inflection point would be where the rate of change of the growth rate switches from speeding up to slowing down (i.e., X(t) acceleration is zero)—equivalently, where the absolute growth dX/dt is largest. Because we assumed , then the inflection point would lie at xinfl = 0.5 (half of organ capacity).

To introduce proliferative potential, we multiply  by the remaining fraction of replicative potential P/H which decreases from 1 to 0, depending on how many cellular divisions have passed. P/H = 1 for new cells at t = 0, later becoming P = 0 and turning proliferating cells into the senescent ones. Therefore, the division rate—in other words, the fraction of cells that divide and increase current population due to replication across dt—becomes (Supplementary Fig. 8):

(24)

Deriving the rate of proliferative potential change

To describe how average proliferative potential changes over time dP/dt = − ψP (where ψ represents the rate of potential loss), we introduce total potential Ψ—that is, the number of divisions left across the whole current cell population: . Since P is average potential, all cells have the same amount of it. Thus, Ψ = P ⋅ X.

When an i-th cell dies, the system loses its P remaining divisions, in average: . Then, net loss of Ψ due to cell death over dt is dΨdeath/dt = − P ⋅ μX. When an i-th cell divides, the system “loses” one mother cell with its P divisions and “gains” two daughter cells, each with P − 1 divisions left. The net change in Ψ per cell division is . Over dt, the number of divided cells is ρX. Summing up all divisions that occur over this period, net change of Ψ becomes . Combining Ψ loss from both death and division, we get

Since Ψ = P ⋅ X, then, according to the product rule for derivatives,

We can rearrange it to solve for dP/dt:

Dividing the whole equation by X and expanding the result, we obtain

Remembering that dP/dt = − ψP, we can now derive ψ (which is the proportion of potential lost on average across all cells over a given year) as follows:

(25)

Since  and , the r/H coefficient becomes . The  part is dimensionless, therefore , yielding loss of proliferative potential averaged across dt = 1 year, as we desired.

Full system of equations for Model IIIA

Combining the derivations above, the resulting system for Model IIIA (somatic only) becomes:

(26)

Model IIIB: Death and replication of somatic and stem cells

After birth, most human tissues harbor a population of progenitor cells which replenish the somatic cell pool, preventing a tissue from the Hayflick-determined lifespan. Historically, these cells have been called various synonyms, depending on the tissue of origin: adult stem cells, somatic stem cells, resident stem cells, progenitor cells, etc. For convenience, in this work, we refer to all of them simply as stem cells (SCs) because we do not consider the germline for our modeling. SC proliferation is a continuum112: in some tissues, SCs proliferate constantly, in order to support the fast turnover of their somatic descendants, the most prolific examples of which are hematopoietic stem cells producing blood and immune lineages93, intestinal crypt cells sustaining intestinal epithelium lining113, epidermal stem cells responsible for skin renewal114. Other SCs mostly lie quiescent and become activated only after a serious injury occurs, so there is urgent need to replace the dead somatic cells and restore organ function: for example, muscle satellite cells, adipose stem cells, etc.112. Some other tissues are virtually devoid of SCs and are, supposedly, not replenished (we employ Model II for those).

Modeling SCs that belong to rapid turnover niches (such as hematopoietic stem cells) appears challenging in our framework, because these are especially prone to clonal expansion, meaning that somatic mutations can hardly be assumed independent, and that an organ can deteriorate functionally, while still being able to sustain its cell counts. Moreover, the plasticity of these SCs allows them to choose their fate dynamically, such that the proportions between self-renewing division and differentiation into various somatic lineages (multipotency) might change dynamically as well, making the models ever more cumbersome. For the quiescent SCs, on the other hand, their data on somatic mutagenesis and proliferation parameters are markedly more scarce. As a middle ground, we chose liver progenitor cells (LPCs, also known as oval cells) which are known to be bipotent and differentiate into either hepatocyte or cholangiocyte populations of adult human liver: cholangiocytes are epithelial cells that line out the bile ducts, while hepatocytes are the main players responsible for blood detoxification, lipid formation, gluconeogenesis, bile secretion, protein synthesis and storage, etc.115,116,117. With so many jobs, hepatocytes are under constant stress, especially considering their key role in toxin removal. However, they can apparently proliferate quite a handful of times, so that the LPCs are normally quiescent throughout human life, unless severe liver damage is done: e.g., the removal of a sizable bulk of liver cells ( >50% and probably even >70% hepatocyte loss) or large amounts of toxins pouring into the bloodstream115,118,119.

The magnificent proliferative capacity of hepatocytes has recently been attributed to the heterogeneity within their ranks: in mouse livers, there are hepatocyte subtypes which express high levels of telomerase (TERT), thus maintaining their telomere length without senescence120. The members of this TERThigh species represent approximately 30% of murine hepatocyte population, reside further away from the ducts, and are likely the ones responsible for restoring the common hepatocyte pool by differentiating into their more abundant, TERTlow siblings, which either reached their replicative senescence or succumbed to external damage. In adult humans, however, TERT expression‎ has not been reported121,122. Other putatively interesting hepatocyte subtypes have been described in humans and mice116,123,124,125, but incorporating all of them into our model is currently unfeasible, because there is no data on replication rates and mutagenesis with such subtype-level resolution. Taking all that into account, we treat hepatocytes as a single population, relying on the average estimates of their parameters.

Our assumptions for this model are:

•  
•  
•  
•  
•  
•  
•  
•  

To incorporate SCs into our modeling, we need to add a few terms and a separate equation:

(27)

where

Modeling rates related to stem cell death and division

As for somatic cells, μyY is calculated using baseline mutation accumulation rate μ0∣y and the probability that a given mutation becomes lethal for a stem cell:

(28)

To model both σ and ρy, we can employ a formula similar to somatic replication: ry ⋅ g(x), where ry is maximum stem cell replication rate, and a sigmoid gate g(x) is responsible for the non-linear population growth between homeostasis and recovery. Given that LPCs engage in replication only after a significant portion of somatic cells is gone (x = X/K is low), we can employ a more flexible, Richards-type logistic gate (in one of its reparameterizations)126: g(x) = (1−x)m. To find the inflection point (the point of zero acceleration), we can solve the d2X/dt2 = 0 equation and obtain:

Thus, if we want to set the inflection point to xinfl = 1/3 (maximum population growth at 1/3 of organ capacity), then we can use . Evidently, when m = 1, the Richards-type logistic becomes the simple logistic g(x) = 1 − x, with xinfl = 1/2 (maximum population growth at half of organ capacity), which we use for all other cell types except for LPCs.

However, LPCs proliferate not only after large hepatocyte loss (when x = X/K is low), but also in response to toxin-induced damage to the LPC pool (when y = Y/Q is low)115,119, so the proliferation law for LPCs has to reflect both hepatocyte density and LPC niche occupancy. We can modify our Richards-type gate as a function g(x, s) of both x and s in various ways, the simplest of which are multiplication (logical AND-like) and addition (logical OR-like):

Because the addition lets either hepatocyte loss OR low stem occupancy trigger high proliferation (Supplementary Fig. 10, we chose the additive gate to represent LPC proliferation in our model:

(29)

where mx and ms are shape parameters that describe how sharply each gate falls off as its corresponding population fraction (x = X/K or y = Y/Q) increases. We will further refer to g(x, y) as g for simplicity.

The resulting product of ry ⋅ g is total stem division rate . To find which proportion of this rate yields actual SCs (ρy) and not the differentiated ones (σ), we should adjust total rate  by the appropriate balance (denoted as γ) between the mentioned division types an SC can undergo. Whichever is the fraction of asymmetric divisions, the stem pool neither loses nor gains any cells from it (0 ⋅ fxy). The number of SCs produced (or lost) due to their division is therefore a balance between the fraction of self-renewing divisions (+1 ⋅ fyy) and the fraction of differentiation events (−1 ⋅ fxx), yielding:

(30)

Somatic cells can arise either as 2 daughters of a symmetric differentiation (+2 ⋅ fxx), or as 1 daughter of an asymmetric one (+1 ⋅ fxy), while the effect of symmetric self-renewal on somatic population is zero (0 ⋅ fyy). Therefore, net somatic gain is γsoma = (2fxx + 1fxy + 0fyy), and

(31)

Deriving the rates of proliferative potential change considering stem cells

Finally, to obtain a formula for proliferative potential replenished by the fresh, stem-derived somatic cells, we once again resort to the concept of total potential Ψ = P ⋅ X (number of divisions left across the whole population). For somatic population only, Ψ change per division was either due to cell death (by ), or due to cell division . This time, we have the third player—stem-driven replenishment. When an i-th cell differentiates from an SC, it increases total potential by . Over dt, the number of newly differentiated cells is σY. Thus, total potential change due to differentiation over this time is . Combining the three equations yields

As in the somatic-only model, Ψ = P ⋅ X, allowing us to derive:

Inserting the combined equation for dΨ/dt yields:

As before, we divide the whole equation by X, expand the result, and obtain

(32)

Remembering that dP/dt = − ψP + ψyY and that − ψP = − 2ρ, we can now derive ψy (the proportion of potential gained from newly differentiated cells per year) by replacing σ with its definition and dividing  by Y:

(33)

The coefficient 1/X can be called a dilution coefficient: it adjusts the input of new cells into the existing potential by dividing it over the current population. Thus, if X is high, then the effect of stem-derived cells on average potential is tiny. On the contrary, when X drops, every new differentiation markedly impacts the population’s proliferative potential.

Similarly with the (H − P) multiplier, which can be thought of as the “room for improvement” or as the magnitude of potential deficit—it is the difference between the maximum possible potential (H) of a new cell arriving from the stem pool and the current average potential (P) of somatic population. When somatic cells are young, P is close to H, meaning that (H − P) → 0. The arrival of newly differentiated cells does not make a huge impact on average potential, because current population is already near-perfect. As somatic cells grow old, P → 0, meaning that (H − P) → H. The potential deficit is enormous: every new cell with P = H provides a massive boost to the average potential.

Thus, the (H − P)/X factor creates the situation such that the restorative effect is at its strongest when the tissue is most exhausted. Curiously, this factor arises by itself from our simple principles; recall that we did not assume such dilution or deficit properties for our equations from the start.

Full system of equations for Model IIIB

Summarizing equations (28), (29), (30), (31), (32), and (33), the resulting system for Model IIIB (somatic & stem) becomes:

(34)

Model IIIC: death of differentiating stem cells with limited proliferative potential

Among the cell types whose single-cell profiles of somatic mutagenesis were obtained, there is one peculiar species: of human bronchial basal cells (hereafter referred to as BBCs, also called airway basal cells and HBECs). These multipotent cells give rise to the many cell types of lung airway epithelium, playing the key role in its homeostasis and regeneration, which makes them a valid member of the somatic stem cell family. However, their proliferation was also shown to be finite in culture: earlier studies at atmospheric 21% O2, together with a ROCK (Rho-associated coiled-coil kinase) inhibitor resulted in BBC senescence at 80 population doublings127, while milder, 2% O2 conditions allowed them to survive 200+ doublings128. Even though the in vivo data are scarce, the fact that BBCs require exceptionally mild conditions to continue proliferating points us towards assuming that their division in vivo is also limited.

Recently, it was shown that there are actually two morphologically identical subpopulations of airway basal cells (in mouse trachea), one of which is basal stem cells, while the other is their luminal progenitors129. Only the progenitors can differentiate further, thus serving as a transient cell type between the stem cells and the secretory or ciliated cells. However, whether the same two subtypes exist in human airway epithelium and, if yes, then what their respective mutation rates are, remains unresolved. Therefore, we treat airway basal cells as a single population, relying on the average estimates of their parameters (as with hepatocytes).

Our assumptions for this model are:

•  
•  
•  
•  
•  
•  

To model BBC behavior, we need to re-adjust our system:

(35)

As with other cell types,

(36)

Modeling BBC replication rate

The rate of division ρb can be modeled with the similar rb ⋅ gb as for other stem cells, adjusted by the current fraction of divisions left P/H. As there is little information on BBC-specific patterns of population growth, we assume a simple sigmoid (1 − B/K), where K is total BBC capacity of human airways. Combining maximum division rate rb, sigmoid growth based on current population size (1 − B/K), and the proliferative potential coefficient P/H, we obtain total division rate:

Since only a fraction of this rate yields BBCs, while the rest goes away to differentiation, we need to adjust  by the fraction of BBC-producing divisions among all divisions, denoted as γb. As with model 3, the BBC pool neither loses nor gains any cells from asymmetric divisions (0 ⋅ fxb), and the net gain of BBCs due to their division is a balance between self-renewal (+1 ⋅ fss) and differentiation (−1 ⋅ fxx), yielding:

Combining total division rate  and this fraction of BBC-producing divisions γb, we obtain BBC-producing division rate:

(37)

As for dP/dt, we resort once again to the concept of total system potential. The BBC system has no external source of new potential—it is only consumed, even though BBCs are multipotent (stem-like). The degradation of average potential P should be proportional to the total rate of division events occurring in the population, as each division consumes BBC potential.

Deriving the rate of proliferative potential change in BBCs

Due to cell death, total potential Ψ changes over time by dΨdeath/dt = − P ⋅ μbB. The change of total potential due to divisions can be found through multiplying the net change of potential per division  by the total number of divisions over time . Due to any division, the BBC pool loses 1 mother cell with potential P and gains daughter cells, each with potential P − 1. If both daughters are BBCs, then ; if only one daughter is a BBC, then the amount of potential added to (or, in other words, remained within) the BBC pool is . Symmetric differentiation only consumes mother potential, without adding any daughters to the total pool . Total potential change due to different division types is

Thus, combining BBC death and division,

Similarly to the previous models, Ψ = P ⋅ B, allowing us to derive:

Inserting the combined equation for dΨ/dt yields:

As before, we divide the whole equation by B, expand the result, expand γb, and obtain

Remembering that dP/dt = − ψbP, we can now derive ψb (rate of potential loss per year):

(38)

Full system of equations for Model IIIC

Summarizing the equations for dB/dt and dP/dt, the resulting system for Model IV becomes:

(39)

For clearer understanding of the last equation for dP/dt, we can expand the second bracket and replace :

There is no fxx-containing term, because differentiated cells make no change to the average potential of BBCs.

Processing data on somatic mutations

The landscape of somatic mutations is challenging to profile, as they occur in each cell independently and therefore can hardly be detected in bulk tissue samples. Traditional genomic techniques were able to detect only large-scale chromosomal alterations (e.g., aneuploidy, megabase copy number variations), while the estimation of single nucleotide variants (SNVs) and insertions/deletions (indels) remained unfeasible130,131. To overcome these limitations, techniques like single-cell cloning and laser capture microdissection (LCM) were developed. However, their applicability is constrained by tissue type—particularly for non-proliferating tissues. The emergence of single-cell whole-genome sequencing (scWGS) now provides the most robust method for measuring the somatic mutation burden within individual cells, including rare, non-clonal variants94. ScWGS almost always requires whole-genome amplification to produce enough DNA for sequencing, which can be done in several ways. Among these, the most widely used are MDA132 (multiple displacement amplification) and PTA133 (primary template-directed amplification). Although MDA was among the first methods to produce large DNA fragments that yielded good coverage with a relatively straightforward protocol, it suffers from allelic dropout (one allele may fail to amplify), biased genome coverage, and amplification artifacts (artificial SNVs/indels), especially due to late amplification errors and single-strand dropout85,133. PTA allows more uniform, accurate, sensitive, and reproducible profiling of somatic mutations than any other existing scWGS approach74,85,133. We therefore prioritized PTA-based scWGS datasets for our analysis. For tissues where PTA data were unavailable, we utilized MDA-based datasets as an alternative.

For Model II (only death, no replication), we focused on brain neurons and heart cardiomyocytes. Models IIIA (somatic only) and IIIB (somatic & stem) were employed for liver hepatocytes and liver progenitor cells (LPCs). Model IIIC was run on human bronchial basal cells (BBCs) retrieved from proximal lung airways. Model II was additionally run on hepatocytes and BBCs to check their survival dynamics in the absence of any replication (Supplementary Fig. 3).

For neurons, we obtained Variant Call Format (VCF) files, metadata, and/or other data tables from the studies by Luquette et al.74, Ganz et al.75, and Motyer et al.76. From these, we selected only neuronal samples sequenced with PTA. For heart, we downloaded supplementary data from the PTA-based pre-print study on left ventricle cardiomyocytes by Choudhury et al.77. Data on liver hepatocytes and LPCs came from the study by Brazhnik et al.134, while BBCs were profiled by Huang et al.135. For both liver and lung, we retrieved the respective VCFs (genome build hg38) from SomaMutDB—a database of somatic mutations136, while sample metadata were downloaded from the supplementary sections of these articles. From all datasets, we selected samples exclusively from healthy controls. For lung tissue, we further restricted to non-smokers. We excluded outliers with exceptionally high or low mutational burden, as well as samples with extended post-mortem intervals (PMI ≥ 50 hours). The resulting outlier sample list is: “4638-Neuron-4”, “CT2_PL_Neu2”, and “CT2_PL_Neu3” (Brain); “1039_A1”, “1039_A5”, and “5828_E6” (Heart); “N1274_2”, “N1274_8A”, “N1276_4”, “N1407_2”, “N1410_1”, “N1410_6”, and “N1415_3” (Lung).

A summary of all analyzed somatic mutation datasets and the resulting sample counts is provided in Table 3.

Table 3 Datasets of somatic mutations included in the analysis

Full size table

The hg19-based datasets were converted into the hg38 assembly using Python liftover (v1.3.0) package (https://github.com/jeremymcrae/liftover).

Estimating accumulation rates μ 0 of somatic mutation burden

For datasets from Luquette et al.74 and Ganz et al.75, we calculated mutation burden per cell from raw VCF files by counting mutation calls and dividing by the reported sensitivity (48.7% and 46.2% for SNV and indel calling in Luquette et al.; 48% and 41% for SNV and indel calling in Ganz et al.). For the remaining datasets, we obtained mutation burden values directly from supplementary materials, which were already sensitivity-adjusted.

Accumulation rates per year of mutation burden per cell were estimated with a mixed-effect linear model accounting for donor-specific affects. Mixed-effects model fitting was performed for each organ separately, and for SNVs and indels separately, using the statsmodels137 Python package and its function mixedlm(mut_burden ~ age + (1 ∣ donor_id)), with REML (restricted maximum likelihood)-based estimation turned on. Accumulation rates were retrieved as slope coefficients returned by this function. Slope confidence intervals (CIs) were calculated from regression values assuming normal distribution and using the scipy138 Python package function norm.ppf, at the 95% confidence level.

Unfortunately, indel calls were unavailable from Choudhury et al.77 and Motyer et al.76. Additionally, Luquette et al.74, Ganz et al.75, and Brazhnik et al.134 did not report per-cell indel burdens, though their datasets (retrieved from supplementary files or SomaMutDB) contained raw indel calls. To estimate indel accumulation rates for brain and liver cells, we scaled raw indel burdens by the ratio of adjusted-to-raw SNV burdens, assuming similar detection sensitivities for SNVs and indels within each sample, then applied mixed-effects modeling. Since cardiomyocytes share post-mitotic, long-lived characteristics with neurons, we hypothesized similar SNV-to-indel accumulation ratios. In the absence of direct cardiomyocyte indel measurements, we inferred these rates by dividing cardiomyocyte SNV accumulation parameters by the neuronal SNV-to-indel ratio.

SNV burdens for LPCs did not show an upward trend with age (slope = −16.176, 95% CI = [−83.627, 51.275], P-value = 0.64, intercept = 1021.476), likely due to the limited sample size (10 samples from 3 donors, all below 18 years). In the original article, Brazhnik et al.134 explored LPC-derived organoids to compare their SNV accumulation rates with hepatocytes.

Rather than employing the LPC-derived organoid data from Brazhnik et al.134, we maintained an in vivo framework by finding the hepatocyte-to-LPC ratio of total SNVs accumulated per cell in young (≤36 years) individuals (similarly to the analysis by Brazhnik et al.134 demonstrated in the paper’s Fig. 1C), obtaining values of 1.565 (SNVs) and 1.671 (indels). LPC mutation accumulation rates were thus estimated by dividing the hepatocyte-modeled coefficients by these ratios.

We note that our empirical accumulation rates μ0 reflect mutation burdens accrued under conditions that include age-progressive oxidative stress, which is itself an aging hallmark assumed eliminated in our framework. Aging-independent mutation accumulation rate would therefore likely be somewhat lower than our empirical μ0, but such value is very difficult to estimate.

Results of the accumulation rates estimation are demonstrated in Table 4.

Table 4 Mutation accumulation rates per tissue according to the mixed-effects modeling

Full size table

Estimating the probability of lethal somatic mutations

To find the lethal mutations rates (μ), we adjusted mutation accumulation rates obtained in the previous section (μ0) by the probability of a mutation to be lethal for a cell of a given type (plethal). Estimating this probability is no trivial task, because, to our knowledge, there are no data which could yield these probabilities directly. Therefore, before proceeding with the estimations, we first eval‎uated the upper and lower bounds as sanity check—any estimated probability should be within this range.

Establishing lower/upper bounds for lethal mutation probability

We estimated the lower bound as the probability that a mutation occurs in the coding sequence (CDS) of a single gene whose disruption is lethal for any human cell type. The largest subunit of human RNA polymerase II (POLR2A gene) can serve as a representative example of such gene. The CDS of POLR2A (the principal transcript: RefSeq NM_000937.5) is 5,913 bp long, encoding a protein of 1970 amino acids.

For each 3 base long codon, there are 9 possible SNVs. Among 549 total SNVs (61 codons × 9 mutations), 23 are nonsense changes that generate a premature stop codon. Fraction of nonsense SNVs is thus fnonsense = 23/(61 ⋅ 9) = 0.0419 = 4.19%.

Missense changes make up the largest share of possible SNVs (392 out of 549, which is 0.714 = 71.4%), but only a subset of them will actually knock POLR2A out of function. Deep mutational scanning (DMS)—where every possible single-amino-acid substitution is assayed in pooled functional screens—and other experiments find that on the order of 20–30% of missense variants measurably impair protein activity in enzymes and other essential domains139,140.

Given that POLR2A comprises 29 exons (as per RefSeq GRCh38.p14 assembly), it also has 28 intervening introns. Every intron is flanked by splice sites: the canonical GT at the 5’ end (donor) and AG at the 3’ end (acceptor). An SNV hitting either of those bases typically abolishes proper splicing of the neighboring exon, and is therefore treated as a separate mutational target of Lsplice = 28 × 4 = 112 bp, distinct from the CDS.

The human genome comprises ~3.1 ⋅ 109 base pairs (as per RefSeq GRCh38.p14 assembly), yielding ~6.2 ⋅ 109 base pairs for a diploid cell. Since loss-of-function variants in POLR2A cause haploinsufficiency141,142, a single hit on either copy is sufficient to impair cell viability, justifying the use of the diploid genome size as the denominator. Combining the fractions of nonsense (stop-gain) mutations and deleterious missense variants (taking 20% as the lower bound), and adding the splice site target separately, we estimate the lower bound for the probability of lethal mutations as follows:

To set the upper bound, we resorted to the SNV burdens obtained in the section on “Estimating accumulation rates μ0 of somatic mutation burden”. By age 50, a cell of any analyzed type acquires more than N50 yo≥1000 SNVs (Supplementary Fig. 2). If we conservatively assume that at least one of these is lethal, then the annual rate of lethal SNV accumulation is μ = 1/50years = 0.02 lethal SNVs/year. The implied annual cell death rate would then be μ = 2%/year. Leveraging this μ as a risk variable (usually denoted as λ) in an exponential decay model (as our Model II), we would obtain the cumulative cell loss over a 60-year period of 1 − e−μt = 1 − e−0.02⋅60 ≈ 70%, which contradicts the latest empirical measurements of human cerebral cortex neurons stating that there is no observable neuron loss in cognitively healthy subjects between 25 and 87 years80. This figure can still be used as an upper bound, but we might be able to tighten it more.

To derive a data-grounded upper bound on μ, we used measurements of frontal cortex (FC) neuron counts from the mentioned paper80 (nsubjects = 43; , SD(NFC neurons) = 0.71 ⋅ 109). The reported Pearson’s correlation between FC neuron count and age is r = − 0.10 (P = 0.54), which is non-significant. Rather than treating this as evidence of zero loss, we constructed a 95% confidence interval around the observed r = − 0.10—that is, the range of correlations that are compatible with the observed data—in the sense that observing a sample value of r = − 0.10 would not be a surprising outcome if the true age–count correlation lays anywhere within it. The lower (most negative) end of this interval is then our worst-case estimate: if the true correlation were any more negative than this, the observed r would have been an unlikely event at the 5% significance level. This lower bound is thus the most pessimistic neuron loss rate the neuron count data cannot rule out, which is exactly what we need for an upper bound on μ ().

In other words, we asked: given that we observe r = − 0.10 in a sample of 43 subjects, which is the most negative r that we could still observe if we measured neuron counts across the general population? This lower bound onr could be then transformed into an upper bound onμ, which would then be used to derive an upper bound on plethal.

The issue is that correlation coefficients are not supposed to have a symmetric, normally distributed uncertainty, so one cannot simply write r ± 1.96 ⋅ SE to get a confidence interval. The Fisher z-transform z = arctanh(r) is a standard reparametrization for this issue: in the transformed space, the uncertainty is symmetric, normal, and depends only on sample size: , making a standard confidence interval straightforward to compute. We therefore obtain z = arctanh(−0.10) = −0.100 with . The 95% confidence interval on z is then z ± 1.96 ⋅ SE(z), and we take its lower (most negative) bound, corresponding to the largest plausible loss rate:

Transforming back to a correlation coefficient: {r}_{95 % }=\tanh (-0.410)=-0.389. This is the most negative true correlation between neuron count and age that the data cannot rule out at the 95% confidence level. To convert this worst-case correlation into an annual cell loss rate, we note that a correlation coefficient is just a normalized regression slope: the actual slope of neuron count versus age equals r ⋅ SD(N)/SD(age), telling us how many neurons are lost per year on average. Dividing by the mean neuron count  then gives the fractional loss rate—the proportion of cells lost per year—which is directly comparable to :

For FC neurons this gives  year−1. Because finding comparable data on cell loss is challenging to perform for every cell type analyzed in our paper, we treated this neuron-derived  as universal across all our cell types. Nevertheless, cell type-specific mutation accumulation rates have been obtained in the section on “Estimating accumulation rates μ0 of somatic mutation burden”, allowing us to calculate the cell type-specific upper bounds on plethal. To perform that, we note that the total annual cell death rate μ is constrained by  across all mutation types jointly:

Since we seek an upper bound rather than the exact probabilities per cell type, we assume a shared lethality probability across SNVs and indels, giving:

where  and  are the annual SNV and indel accumulation rates for each cell type, respectively, obtained from the mixed-effects modeling described above (Table 4). The resulting upper bounds are reported in Table 5.

Table 5 Upper and lower bounds on lethality probabilities across cell types used in our modeling

Full size table

The range between lower bound 1.94 ⋅ 10−7 and upper bounds ~10−4 is admittedly wide, but it provides us with the prior boundaries necessary to perform the sanity check of our estimations.

Essential genes annotation

Moving on from the single POLR2A gene, we next sought other genes that could be essential for cell survival in each analyzed tissue. The concept of gene essentiality is a rather challenging topic, with many caveats and unsolved issues143,144,145,146,147. Large-scale studies of genes essential for survival, development, or regeneration are typically performed using CRISPR-based inhibition/activation screens. The in vivo human screens of gene essentiality remain ethically unfeasible. Recently, such screens were performed for neurons which had been differentiated from induced pluripotent stem cells (iPSCs)148,149. However, the iPSC-derived neurons may not fully replicate gene expression‎ patterns and chromatin features of mature human neurons150. Therefore, we relied on in vivo mouse screens to determine our panels of essential genes for each cell type, admitting that human essential genes might not be fully captured in the knockout (KO) mouse models143.

Neuron-specific151,152, cardiomyocyte-specific153,154, and liver-specific155,156 essential genes were obtained from the respective papers of CRISPR-based screening. Genes passing the p-value threshold of 0.1 after the Benjamini-Hochberg multiple testing correction (as provided in the original papers) were selected. Mouse gene symbols were converted into their human homologs using the gseapy157 Python package. Since healthy lung-specific essentiality screens remain unavailable to date, we defined lung-essential genes as those reported as essential in > 50% of human lung-derived cell lines in the OGEEv3 database158. Each tissue-specific signature of essential genes was united with core-essential human genes (genes reported to be essential across 80+% of all human cell lines) from the OGEEv3 database158. The UpSet plot of the resulting signatures and their overlaps is created using the UpSetPlot159 Python package and displayed in Supplementary Fig. 6.

Annotation by deleteriousness

Mutations do not occur randomly across the genome, and every cell type accrues them in a different way; therefore, we could not assume that the mutations accumulated at rate μ0 fall in essential and non-essential genomic regions with equal probabilities. Instead, we relied on the empirically observed mutation distributions across the genome.

The probability of a mutation to be lethal (for each j-th mutation class) was modeled as a frequency of encountering a lethal mutation among all observed mutations:

We incorporate two mutation classes in our model: SNVs and indels. NSNV s and Nindels denote the total numbers of SNVs and indels observed per tissue in the same datasets that we used to estimate accumulation rates μ0 of somatic mutation burden. We estimated the expected numbers of lethal mutations (Nlethal SNV s and Nlethal indels) by computing lethality scores ci for each mutation—the probability that that individual mutation is lethal—and summing them across all observed mutations:

The lethality score c is a compound metric which consists of three key factors: 1) whether a mutation falls into a region associated with an essential gene, 2) how deleterious this mutation is, and 3) the probability of this mutation to be haploinsufficient:

The essentiality coefficient kess takes values 1 or 0 (essential-related or unrelated) and is inferred by intersecting the VEP- and FAVOR-merged annotation (including genes, their upstream/downstream regions, and CAGE promoters/enhancers) of a given mutation with the tissue-specific list of essential genes (see section on “Essential genes annotation”).

Importantly, the quantity ∑ici does not aggregate per-mutation death probabilities into a cell-level death probability—it estimates the expected number of lethal mutations, Nlethal, among all empirically observed mutations in a given tissue dataset. To make this more intuitive: we could assume a limiting case where pdel = phaplo = 1 for all mutations (that is, every mutation in an essential region is assumed to be unconditionally lethal), then ∑ici reduces to a plain count of mutations falling within essential genomic regions.

Estimating the effect of a given mutation on its associated gene is another challenging issue, especially for the non-coding regions. To do that, we annotated each mutation using the Ensembl Variant Effect Predictor (VEP)78 and the Functional Annotation of Variants—Online Resource (FAVOR)79 databases. Then, we calculated pdel in two steps: first, we aggregated scores from several deleteriousness predictors available in these databases by taking their median. Second, since raw predictor scores do not directly represent the probability of deleteriousness, we applied score calibration.

Pejaver et al.160 show posterior probability curves of multiple pathogenicity predictors, from which we inferred that the predictors with the strongest evidence for pathogenicity can be calibrated by a power law pdel ≈ 0.9 ⋅ (score)6, with many weaker predictors falling even further below this curve. Because their study did not include all predictors we used in our work, and therefore we could not calibrate all predictors separately, we applied the same calibration to our scores pdel = 0.9 ⋅ (median score)6.

Raw CADD scores were transformed using a sigmoid model , with A and B parameters estimated by Benevenuta et al.161. PHRED-scaled scores from FAVOR (protein function annotation PC and conservation annotation PC) were rescaled to the [0, 1] range via . Haploinsufficiency probabilities (phaplo) were taken directly from the Ensembl VEP annotation.

Because there are no data on indel burden in cardiomyocytes, we estimated plethal for cardiomyocyte indels using the same plethal∣SNV s/plethal∣indels ratio as in brain neurons (similarly to how we estimated indel accumulation rate for cardiomyocytes from the brain data in the “Estimating accumulation rates μ0 of somatic mutation burden” section).

After obtaining plethal SNV s and plethal indels for every tissue, we calculated μ as:

The resulting plethal, μ0, and μ for every cell type are summarized in Table 1. The 95% confidence intervals (CI) on μ are propagated from the regression-derived uncertainty in μ0, with CI bounds computed assuming normality of the slope estimates.

Estimating other model parameters for each organ

Aside from the mutation accumulation rate μ0 and the probability of a mutation to be lethal plethal, we needed to estimate other parameters to be able to solve our ODEs: organ capacity K, critical thresholds of cell population necessary for organ survival Xcrit, maximum replication rates r, ry, and rb (for hepatocytes, LPCs, and BBCs, respectively), Hayflick limit H, liver stem capacity Q, Hill-type sigmoid parameters θ and n (for the liver model), and all fractions of different division types that the LPCs and BBCs can undergo.

Brain modeling parameters

Modeling brain as a whole single organ is difficult, because its regions are rather heterogeneous, with varying energy demands and neuron/non-neuron ratios. As most neurons in the analyzed datasets came from the frontal cortex, we focused on this region to provide the parameters for our model. While there have lots of different estimates of neuron counts in the human brain throughout history, the latest research reports that frontal cortex contains (3.5 ± 0.7) ⋅ 109 neurons (40% of brain mass, 34% of brain neurons), with no significant difference between males and females and with no significant decline across the human lifespan80.

Defining the critical threshold of cortex neurons necessary for human survival is rather vague: data from functional hemispherectomy (removal of one cerebral hemisphere) in children and adults162,163 and frontal lobotomy history164 can be used to argue that the frontal cortex is certainly important165, but not required for vital autonomic control, suggesting a relatively low survival floor for the X/K ratio. However, the same works show that people can survive biologically with far less cortex, but cognition (in the everyday, independent-function sense) fails much earlier—therefore, we switch from “vital for survival” to “vital for cognition” (that is, not demented; capable of independent living). Still, no study gives a neuron-count “cut-point” for being cognitively normal. We inferred the threshold from human data on how much neuronal loss (or closely tracked proxies like cortical thickness/atrophy) accompanies clinically manifest dementia in Alzheimer’s disease (AD) and frontotemporal dementia (FTD). According to those, once ~30–50% neuronal loss/marked cortical thinning accrues in the association cortex (including frontal cortex), cognitive deficits are the rule166,167,168; therefore, we place the cognitive “floor” in the middle of this range of neuronal loss, at Xcrit ∣ brain = 0.6K.

Heart modeling parameters

According to the stereological and 14C birth-dating studies, the most cited cardiomyocyte (CM) numbers are as follows. The whole heart CM estimates vary widely169: some researchers find increase from ~109 at birth to ~4 ⋅ 109 in adults (≤20 years)170, while others report CM numbers already high at 1 month of age, remaining constant through life at (3.2 ± 0.75) ⋅ 109171. We treat (3.2 ± 0.75) ⋅ 109 as a whole-heart estimate with reasonable uncertainty for adults.

The number of cardiomyocytes is reported to remain constant during the human lifespan171, and CM turnover occurs mostly during the first two decades of life; after that, turnover rate declines exponentially in adults from ~1%/year at the age of 20 years to ~0.5%/year in the elderly171,172. Even though some heart regeneration appears possible172, stem cell activity is negligibly minimal, because it is evidently not sufficient to restore heart function after severe injury due to myocardial infarction. Loss of ≥40% left ventricle mass (that is, when ~0.6K remains) was reported to necessitate cardiogenic shock and very high mortality173. Even with recent advances in medicine, long-term patient survival is strongly compromised in such severe conditions174. Admittedly our model of gradual CM loss should not be matched directly to the cases of acute myocardial infarction, but in case of radical life extension, on the other hand, the remaining cardiac output must be able to sustain human body for a lot longer than the current human lifespan. To account for a possibility of long-term adaptation to the gradually lowering cardiac output, we set critical CM threshold slightly lower than 0.6K, at Xcrit ∣ heart = 0.55K.

The resulting K and Xcrit for neurons and cardiomyocytes are summarized in Table 6.

Table 6 Variables included in the proposed Model II for human brain neurons and heart cardiomyocytes

Full size table

Liver modeling parametersLiver somatic and stem capacities

We calculated liver hepatocyte capacity via multiplying total liver mass by hepatocellularity (hepatocyte density, in cells per gram). For liver mass, we took estimates from a large autopsy-based study performed by Bell et al.175. Regarding hepatocellularity, we used the protein ratio-based measurements by Sohlenius–Sternbeck176 estimating (139 ± 25) ⋅ 106 hepatocytes/g. LPC definitions vary in terms of namings, morphology, and cellular markers (liver progenitor cells, hepatic stem cells, EpCAM+ adult hepatic progenitors, “oval cells”, biliary epithelial cell–derived progenitors, mesenchymal liver stem cells, etc.), introducing notable uncertainty to their population counts. For a mean estimate of Q (LPC capacity), we used the most recent measurements from spatial and single-cell transcriptomics data177 reporting LPCs as 2.95 ± 1.91% of epithelial liver population. For the lower bound, murine marker-based measurements178 reporting LPCs as 0.5% of epithelial liver population can be used. Epithelial liver population comprises LPCs, hepatocytes, and cholangiocytes, therefore we also retrieved hepatocyte (60 ± 10%)179,180 and cholangiocyte (4 ± 1%)181 fractions of total liver cells to compute Q.

Monte-Carlo sampling (500,000 samples) was performed to propagate uncertainty from the respective studies. All parameters were modeled as normal distributions with truncation applied to maintain biological plausibility. The analysis computed total epithelial cells as hepatocytes multiplied by an epithelial multiplier derived from cholangiocyte and hepatocyte fractions, then calculated LPCs as a fraction of this epithelial population. As a result, we obtained the following liver cell capacities:

•  
•  

Liver critical threshold

Before liver surgery, the future liver remnant (FLR) measure (as determined by computer tomography volumetry) is used as a clinical proxy for how much liver remains after planned resection (partial hepatectomy, also called PHx)—and hence, how many functioning hepatocytes remain. Multiple reviews and clinical practice sources converge on FLR of at least 30% for healthy livers182,183. Some meta-reviews even report successful, survivable resections of >70% liver volume184. Therefore, we estimate the critical threshold for hepatocyte population required for liver survival as Xcrit ∣ liver = 20 ± 5% of hepatocyte capacity = (0.20 ± 0.05) ⋅ K.

Hepatocytes proliferation limit

Albeit capable of sustained proliferation in vivo, human hepatocytes do not divide as readily in vitro, rapidly losing their proliferative ability and function in standard culture and typically surviving only a few divisions before senescence or loss of phenotype185. The approaches to culturing murine hepatocytes are also difficult to translate to the human ones185, suggesting that the TERThigh hepatocyte subspecies featured in murine livers120 might not be as prominent in humans. Excluding the studies of hepatocyte de-differentiation, we can rely on a study which reached at least 40 hepatocyte population doublings185 for the lower bound of our proliferation limit: Hlower = 40. From the upper side, we set the bound at Hlower = 200 based on the adult human hepatocytes immortalized via telomerase overexpression‎ and inhibition of tumor suppressors186. Due to the lack of a better measure, we estimated H for hepatocytes via geometric mean:  divisions.

Hepatocytes replication rate

In healthy human liver (i.e., steady state), diploid hepatocytes replicate on average at rlower = 0.71 year−1187. To estimate r, we leveraged datasets of human liver regeneration after partial hepatectomy with varying resection percentages (that is, varying FLR). During such a short timeframe, we assume that cell mortality is negligible (μ → 0) and proliferative potential loss is also negligible (P/H → 1). Because we also know from literature that LPCs engage in liver regeneration meaningfully only at 2/3 liver loss and more, the dX/dt equation from our Model III becomes very simple:

We grouped 9 datasets from 6 different PHx studies based on their X0/K ratio (X/K immediately after resection) as follows. Lower bound (LB) group: X0/K = [0.70, 1.00] (less than 1/3 resection, regeneration is not at its fullest)188,189,190; upper bound (UB) group: X0/K = [0.00, 0.35] (large resection encompassing 2/3 PHx and slightly less than that; LPC division contributes to regeneration, so this group cannot be used to give a mean r estimate)191,192; and a mean estimate (MEAN) group: X0/K = (0.35, 0.70) 189,190,192,193. We fit X(t) curves using the equation using the maximum likelihood approach for each dataset propagating the reported CI intervals, and then inferred a combined curve and r for each group (Supplementary Fig. 7), which resulted in r values and its CI range displayed in Table 7.

Table 7 Hepatocyte replication rates inferred from datasets of liver regeneration after hepatectomy

Full size table

Liver progenitor cell replication rate

As mentioned above, the nature of liver stem cells is relatively elusive, and different approaches to defining and quantifying them exist, which makes it challenging to estimate their parameters in our model: maximum stem proliferation rate ry, sigmoid gate parameters θX, θY, nx, and ny, and the three division probabilities fyy, fxx, and fxy. Based on the numerous evidence of LPCs being activated only upon >70% of liver loss115,118,119, we assumed both θX = 1/3 and θY = 1/3 (which reflect X/K and Y/Q ratios, respectively, at which half of LPCs are activated and proliferate).

For ry estimation, we relied on the upper/lower bounds approach: as upper bound, we inferred ry ∣ upper from primary-culture population doubling times using:

Calculating ry ∣ in vitro from two studies of adult (Td = 62h)194 and fetal (Td = 46h)195 human LPCs yields ry ≈ 0.268 day−1 and ry ≈ 0.362 day−1, respectively. Another, later study cultured human hepatocytes-derived LPCs (HepLPCs), biliary epithelium cells-derived LPCs (BecLPCs) and the actual, resident LPCs-derived LPCs (reLPCs) and demonstrated notably shorter doubling times (~25.8, ~28.5, and ~21.6 h, respectively)196, which converts to much higher proliferation rates of 0.645, 0.584, and 0.770 day−1, respectively. Evidently, even though LPCs do not appear to proliferate in homeostasis (ry ∣ lower → 0), they are still capable of a much higher activity in a wide range of conditions. We provide an estimate for ry via a geometric mean of adult LPC-derived upper bounds and an approximation for rs ∣ lower ≈ 10−4 day−1 ≈ 0.037 year−1:

This value is still lower than the lowest in vitro-inferred upper bound (0.268 day−1), and is not pegged down to zero, as it would be if we considered ry only at homeostasis.

Liver progenitor cell division probabilities

The proportions between division types change relatively to the current deficits of both stem (Y/Q) and somatic (X/K) cell populations. In homeostasis, asymmetric divisions (fxs) predominate, and the symmetric self-renewal (fyy) and differentiation (fxx) cancel each other out, according to the neutral-drift models197. After a severe injury causing widespread cell loss, stem cells self-renew rapidly in order to later switch to symmetric differentiation and quickly replenish the organ’s functional capacity198,199. Therefore, all division type probabilities should be modeled as functions of both y = Y/Q and x = X/K. At any y or x, the division fractions must be non-negative and sum to one: fxx + fyy + fxy = 1.

A convenient choice to model such relationships is to form utilities for each division outcome that are linear (or affine) functions of y and x, and then convert these utilities to probabilities by a softmax200. The softmax approach enforces positivity and the sum-to-one automatically. Because the softmax depends only on utility differences, one of the utilities can be fixed as a reference. We choose Uxy = 0 as the baseline (equivalently we could subtract the same constant from all utilities). This zero baseline allows us to avoid an extra set of parameters that would otherwise be redundant and ensures identifiability: without fixing one of the utilities, the parameters are under-determined. Thus, the system of utilities can be parametrized as:

From these utilities, we can derive the probabilities of each division type as:

(40)

As one can infer, , which supports the sum-to-one rule.

The softmax lets us transform observed probabilities into log-odds (utilities) very simply. Divide each observed probability by the baseline fxy and take the natural logarithm:

For each state (y, x) and triplet (fyy, fxx, fxy) observed at this state, we thus have two linear equations as the ones above.

Each parameter has a simple interpretation:

•  
•  
•  
•  

Because utilities enter the probabilities exponentially, a change of 1 unit in a utility multiplies the corresponding odds by e1 ≈ 2.7—so magnitudes are interpretable in odds terms.

Given several anchor states (y, x) with the respective (fyy, fxx, fxy) triplets, we can solve the softmax system algebraically. The problem with LPCs, however, is that the precise division fractions are hard to find in the literature. One study198 quantified symmetric self-renewals and asymmetric divisions of Lgr5+ stem cells in murine livers injured via toxin administration (Lgr5+ supposedly marks damage-induced LPCs)115. The injured mice were treated with either FXR (farnesoid X receptor) agonist to force LPC symmetric self-renewal, or PPARα (peroxisome proliferator-activated receptor α) agonist to promote asymmetric divisions. After FXR treatment, the probabilities were demonstrated as fasymmetric = fxy ≈ 15%, fsymmetric ≈ 85%. After PPARα treatment: fasymmetric = fxy ≈ 80%, fsymmetric ≈ 20%.

Because PPARα activation drives LPCs to quiescence and maintenance, we assumed that the PPARα-induced fractions can be used to approximate division probabilities at homeostasis. We decomposed the fraction of symmetric divisions into fyy ∣ homeo = fxx ∣ homeo = 1/2 ⋅ fsymmetric ∣ PPARα ≈ 20/2 = 10%, so that self-renewal and differentiation cancel each other at homeostasis, leaving most divisions asymmetric fxy ∣ homeo = 80%.

On the contrary, the FXR-induced triplet corresponds to the state of hyperactive self-renewal following injury. In physiological injury response without exogenous FXR activation, differentiation would not be impeded, and a fraction of symmetric divisions would continue leading to it: we assume it to be the same as homeostasis, fxx ∣ stem wave = fxx ∣ homeo = 10%; then fyy ∣ stem wave = fsymmetric ∣ FXR − fxx ∣ stem wave = 85 − 10 = 75%. As this study did not trace the later differentiation of FXR-induced Lgr5+ daughters, we can only assume that symmetric differentiation rises to the same levels as does symmetric self-renewal in order to supply functional hepatocytes to the injured organ from the newly made LPCs: fsymmetric∣stem wave = fsymmetric ∣ differentiation wave. Hence, swapping the two kinds of symmetric divisions between each other, fxx ∣ differentiation wave = fyy ∣ stem wave = 75% and fyy ∣ differentiation wave = fxx ∣ stem wave = 10%. As the exact (y, x) pairs for stem and differentiation waves are not known, we hypothesized these pairs based on our addition-based sigmoid rule for stem cell replication (29) and (30), according to which LPC proliferation rate ρy is high at y = 1/3, at x = 1/3, and continues to be high at y = 2/3 if x is still 1/3.

Thus, the three calibration points become as shown in Table 8.

Table 8 LPC division probabilities inferred from ref. 198

Full size table

Solving the softmax system at these calibration points with the lstsq function of numpy Python package, we obtained:

These utilities are then fed in the softmax system (40) to compute LPC division probabilities at each x and y of our modeling (Supplementary Fig. 9).

Liver modeling summary

The resulting values for all liver-related variables included in our modeling are summarized in Table 9.

Table 9 Variables included in the proposed Models IIIA and IIIB for human hepatocytes and LPCs

Full size table

Lung modeling parametersBasal cell capacity

The data on somatic mutagenesis in lungs came from the proximal bronchial basal cells; however, it is unreasonable to estimate basal cell capacity based only on the proximal bronchi: the airway tree is continuous, and epithelial cells can proliferate along its whole length. Therefore, to estimate basal cell capacity K in human airway epithelium, we employed the measurements of tracheal geometry and histology data of the whole human tracheobronchial tree. The inputs for this estimation are: CT-derived tracheal lumen volume V = 32.0 ± 8.3 cm3 and tracheal length L = 102.8 ± 9.9 mm 201, diameter reduction factors k = 0.79 and 0.94 (depending on generation number) to generate bronchial generation-wise diameters via a symmetric Weibel model202,203, a length-to-diameter ratio LtD = 1.46 ± 0.15 to convert diameters into bronchial generation lengths202, total conducting airway surface area Ameasured = 2471 ± 320 cm2204, total epithelial cell count Nepithelial = 10.5 ⋅ 109204, number of conducting airway generations G = 25 203, basal cell fractions at large and terminal airways205, and sex-specific total lung volumes (TLVs)206.

Each bronchial generation has ng = 2g airways. This is an idealization: asymmetric branching will change counts and therefore area allocation, but 2g is the standard Weibel-type starting point. Trachea is normally oval-like, but we approximate it by computing the equivalent circular tracheal diameter D0 from its volume V and length L as

Diameter change by each bronchial generation g depends on the per-generation reduction factor k, which changes from 0.79 before ~15th generation to 0.94 after it202. We used the following piecewise model, which ensures Dg is continuous at the transition (no artificial jump):

Given a generation-specific length-to-diameter ratio LtD = Lg/Dg (dimensionless), the lateral surface area of one bronchial branch in generation g would be

Each generation’s area is . These areas are relative predictions based on idealized assumptions (perfect symmetry, cylinders, no irregularities, etc.). Morphometry204 reports a measured total conducting airway surface area Ameasured and a total epithelial cell count Nepithelial, so we enforce the empiric total area via the scaling factor

which preserves the model’s relative distribution of areas while ensuring the sum of modeled areas equals the empiric total.

We treat Nepithelial as a reference at TLVcombined sexes. For each sex group with TLV equal to TLVsex group, we scale the epithelial count linearly:

Basal cells are primarily located in the proximal conducting airways (trachea, bronchi, and bronchioles)207. In trachea, basal cells make up roughly 40 ± 5% of the epithelium (at least in mice)129. Normal human bronchial epithelium contains 6–31% basal cells, with proportion varying along the proximal-distal axis205. The basal cell-containing conducting epithelium extends distally to terminal bronchioles of about 0.5 mm diameter, and only the respiratory bronchioles are lined by a simple cuboidal epithelium lacking basal cells208. We used linear interpolation to calculate the basal fractions for G = 25 bronchial generations (between the 31 and 6% basal fraction endpoints): let pg be the fraction of epithelial cells at generation g that are basal. For each g-th generation,

Each generation contains Kg = Ag ⋅ ρsex group ⋅ pg basal cells. Finally, the total basal cell capacity K is

Uncertainty is propagated by Monte-Carlo sampling (500,000 draws): all uncertain inputs (tracheal V, L, LtD, Nepithelial, TLVs, and class basal fractions pg) are drawn from truncated normal distributions to avoid nonphysical negative draws from the normal distributions (implemented with truncnorm from the scipy Python package).

As a result, we obtained the following airway basal cell capacities:

•  

Lung critical threshold

Estimating the critical threshold for the B/K ratio is difficult because there is no direct human experiment that reports a single numeric threshold. From the transplantation studies, we can assume that at least 60 million basal cells are required to restore lung function81, placing the lower bound at . From lung resection studies showing that the net loss of at least 10% lung volume is perfectly tolerable209,210, we can assume that the same proportion of basal cell loss is also tolerable; therefore Bcrit ∣ upper/K = 0.90 = 90%. We estimated Bcrit as a geometric mean of these bounds, making .

Basal cell replication rate

To estimate maximum replication rate for BBCs, we again resorted to the upper/lower bound logic. A steady-state study of murine tracheal basal cells by Watson et al.129 found that these cells divide, on average, once per 11 ± 4.4 days. From it, we can infer the lower bound for BBC replication rate (assuming tracheal and bronchial basal cells have similar cycling dynamics) of:

In two studies by Peters-Hall et al.127,128, BBCs were cultured in different conditions, exploring their maximum number of population doublings (PD): at 21% or 2% O2, with or without Rho-associated protein kinase inhibition (ROCKi), in normal and cystic fibrosis (CF)-derived BBCs. We assumed that under extremely mild, stress-free conditions (2% O2, ROCKi, healthy donor), BBCs show virtually unrestricted proliferative potential (with no differentiation bias). Hence, the settings explored by Peters-Hall et al. can approximate the maximum replication rate.

Both studies feature population doubling curves PD(t). To estimate rb ∣ upper, recall that, under unrestricted growth:

Thus, solving for rb:

After fitting a linear regression on data points (corresponding to 2% O2, ROCKi, healthy donors) provided in Fig. 2A of Ref. 127 and Fig. 1A of Ref. 128, we obtained slopes of ~0.83 and ~0.41 BBC population doublings (PD) per day (dPD/dt), respectively (Supplementary Fig. 11). Estimating rb ∣ upper from each of those yields  and . To find rb ∣ mean and its 95% CI, we calculated geometric mean between the one lower and two upper bounds and propagated the uncertainty from rb ∣ lower measurement error and rb ∣ upper fitting errors (via Monte-Carlo sampling with N = 10,000 samples), yielding:

Basal cell proliferation limit

In standard culturing conditions, airway basal cells proliferate for 30-40 PDs127,211, but the 21% oxygen concentration in these conditions is not physiological: basal cells do not come in direct contact with ambient atmospheric oxygen in the living airways. To mitigate the effect of oxidative damage on BBC replication limits, the ROCKi studies used 2% O2 to reach > 220 PDs127,128, but they did so without reporting the exact O2 concentration of the in vivo BBC niche. Other studies incured hypoxia response in basal cells at < 1%212 and, to some extent, at 5% O2 (considered normoxic for many other tissues)213, suggesting that the 2% conditions are sub-physiological. To estimate BBC replication limit, we set Hlower = 40 from the 21% O2 conditions, Hupper = 220 from the 2% O2 conditions, and perform linear interpolation between these two endpoints, yielding Hmean ≈ 172.6 ≈ 170 divisions at 7% O2 (putative oxygen concentration at the BBC level; slightly above the still hypoxic 5% O2).

Basal cell division probabilities

In contrast to the liver progenitor division probabilities, we cannot derive BBC division probabilities as functions of both b = B/K and some x = X/(differentiated cell capacity), because there are no data on somatic mutagenesis in the differentiated daughters of basal cells yet, and, therefore, we cannot consider them in our model. Hence, we derive the division fractions only as functions of b. We also have to assume that fbb = fxx at any b (assuming differentiation and self-renewal cancel out in the long run), because there is no x dimension which we relied on to model the fbb and fxx dynamics in case of liver progenitors.

The aforementioned study of murine tracheal basal cells at homeostasis129 reports the following division fractions: fxb∣homeo = 0.94 ± 0.03 = 94 ± 3%, fbb ∣ homeo = fxx ∣ homeo = 0.3 = 3%. In another study, influenza infection was shown to shift basal cell division to fsymmetric ∣ injury = 0.48 = 48%, thus fasymmetric ∣ injury = fxb = 1 − 0.48 = 0.52 = 52% 214. From our fbb = fxx assumption, fbb ∣ injury = fxx ∣ injury = 0.48/2 = 0.24 = 24%. This 8-days long influenza infection resulted in the decrease of ciliated area per micron from ~4.2 to ~0.5 and in the decrease of Foxj+ (ciliated airway cell marker) cells per micron from ~0.024 to ~0.005 in infected murine airways214, suggesting that the population of differentiated ciliated cells dropped to levels somewhere between 0.5/4.2 ≈ 0.12 = 12% and 0.005/0.024 ≈ 0.21 = 21%. Assuming that the population of basal cells experienced a similar drop, we set the B/K anchor at injury to 25% (or 1/4).

We leveraged these mouse-based values as the best existing estimate of airway basal cell division fractions at homeostasis and during regeneration and constructed the softmax system, resembling our approach to deriving the LPC division probabilities. Let b = B/K. Similarly to the LPCs, the system of utilities can be parametrized as:

Because of our fbb = fxx assumption, a1 = a2 = a and q1 = q2 = q. With a > 0, decreasing b increases (1 − b) and thus increases Ubb(b) = Uxx(b).

Solving this softmax system for the anchors at injury and homeostasis as before, we obtain a ≈ 3.56, q ≈ − 3.45, and therefore:

Thus, we compute BBC division probabilities at each b of our modeling as (Supplementary Fig. 12):

Lung modeling summary

The resulting values for all lung-related variables included in our modeling are summarized in Table 10.

Table 10 Variables included in the proposed Model IV for human airway basal cells

Full size table

Fréchet–Hoeffding bounds

To quantify the uncertainty in the joint survival function of multiple organs due to unknown dependence among organ failure times, we employed the Fréchet-Hoeffding bounds—sharp extremal limits on any joint distribution with given marginals82.

For a joint survival function  with marginal survival functions , the Fréchet-Hoeffding bounds are:

(41)

Note that these bounds apply to the joint survival at a common time point t (i.e., t1 = t2 = … = tn = t), which is relevant when estimating organismal survival as a time until any organ fails. The upper bound corresponds to comonotonicity (perfect positive dependence), where all organs fail simultaneously; the lower bound corresponds to the most antagonistic dependence compatible with the given marginals. Notably, while the upper Fréchet-Hoeffding bound is always attainable (via comonotonicity), the lower bound is sharp only for n = 2. For n≥3 this bound remains valid but may not be achievable by any joint survival distribution with the given marginals—i.e., it is not necessarily tight. Nevertheless, it is the best possible universal lower bound that depends solely on the marginal survival functions215.

In our case, the number of marginals is 4, which means that the inequalities change to:

(42)

As a reference model under the assumption of mutual independence, we use:

(43)

Finally, to account for the background (non-organ-specific) mortality, we multiply the resulting joint survival function and its bounds by an exponential baseline survival term , yielding the adjusted survival estimates used in lifespan prediction.

Monte-Carlo simulation for empirical survival probability attribution

Because Model III variants provide intractable survival probability, we estimate it using Kaplan-Meier. Failure dynamics were governed by a mechanistic ordinary differential equation (ODE) system describing the interactions between state variables.

To account for uncertainty in biological parameters (e.g., mutation rates, carrying capacities, regenerative rates), we performed Monte-Carlo simulations under a log-normal parameterization (justified by multiplicative biological noise). For each simulation:

1. 1.
2. 2.
3. 3.

The survival function  was estimated non-parametrically using the Kaplan-Meier (KM) estimator, with right-censored observations at . Confidence intervals for the KM curve were computed using the Greenwood’s formula216. Additionally, to model background (non-organ-specific) mortality, the mechanistic survival estimate was multiplicatively combined with an exponential random mortality term , where λBG was estimated from human demographic data.

Note that our simulation of the synthetic population uses samples from the parameters that are estimated for the mean distribution.

Sensitivity analysis of the lethal mutation probability p lethal

To assess the robustness of model predictions in regards to uncertainty in the probability that a somatic mutation is lethal for a cell (plethal), we conducted a sensitivity analysis over the joint space of two key parameters: the probability that a single-nucleotide variant (SNV) is lethal () and the analogous probability for insertion-deletion mutations (). Both parameters were varied independently across a biologically plausible range spanning several orders of magnitude.

Parameter sampling

For each organ, n = 300 parameter sets were generated by independently sampling  and  from a log-uniform distribution over the lower/upper bound interval , where  and , and the latter is varied across organs (Table 5). Sampling on the log scale ensures uniform coverage across orders of magnitude, avoiding the underrepresentation of small values that would arise from linear sampling. Formally, , which is then exponentiated to recover p. A fixed random seed was used to ensure reproducibility.

For Model IIIB (liver hepatocytes and a liver progenitor cell compartment, LPC), an additional pair of lethal probabilities— and —was drawn independently from a separate log-uniform interval, allowing the LPC compartment to be explored over its own.

Mutation rate propagation

For each sampled pair (, ), the effective mean mutation rate μi and its associated standard deviation were computed analytically from the organ-specific SNV and indel rate estimates:

This exact propagation avoids Monte Carlo estimation of μ itself and preserves the dependence of uncertainty in μ on the sampled lethal probabilities.

Survival rate estimation

For each of the n sampled parameter sets, organ-specific survival trajectories were computed using 5000 Monte Carlo replicates. For Model II, each replicate drew a mutation rate from a log-normal distribution parameterised by (μmean, μstd), then the respective survival was eval‎uated analytically using the log-normal cumulative distribution function applied to the organ lethality threshold-crossing condition. For Model III, full ODE simulations were run in parallel, with event times extracted and interpolated onto a common time grid via the Kaplan–Meier estimator. In both cases, organ-specific survival was multiplied by a background (non-organ) survival term, , to obtain the composite survival function.

Hazard rate estimationModel II

For the brain and heart models, we derived a semi-analytical expression‎ for the organ survival probability:

where Φ(…) is the cumulative distribution function (CDF) of the organ capacity K.

Using the identity , we obtain the corresponding hazard rate:

(44)

where ϕ(…) is the probability density function (PDF) of K.

To estimate the expectation over the distribution ξ(μ), we used the Monte-Carlo method with 10,000 samples.

Model III

For the liver and lungs, where a closed-form solution for the survival probability is not available, we employed the nonparametric Nelson-Aalen estimator for the cumulative hazard function Λ(t):

(45)

where ti are the observed event times, di is the number of failures at time ti and ni is the amount of subjects at risk just prior to ti.

To estimate the hazard rate h(t) from the cumulative hazard, we differentiated  numerically and applied a spline smoothing function from scipy.interpolate.UnivariateSpline with the power k = 5. This helps to reduce noise and obtain a smooth hazard rate curve.

Data availability

No additional experimental data were generated during the current study. The data which are necessary to reproduce the results of this paper, including model parameters and demographic data, can be found in the GitHub repository at https://github.com/ComputationalAgingLab/somatic-mutations-simulation.

Code availability

Code scripts and notebooks, including the model and parameters calculation, which are necessary to reproduce the results of this paper, are provided in the GitHub repository at https://github.com/ComputationalAgingLab/somatic-mutations-simulation.

References

1.  
2.  
3.  
4.  
5.  
6.  
7.  
8.  
9.  
10.  
11.  
12.  
13.  
14.  
15.  
16.  
17.  
18.  
19.  
20.  
21.  
22.  
23.  
24.  
25.  
26.  
27.  
28.  
29.  
30.  
31.  
32.  
33.