A hybrid Ornstein–Uhlenbeck–Branching framework unifies microbial and pediatric tumor evolution

Introduction: Pediatric tumors can relapse despite low mutation burdens, suggesting hybrid evolutionary dynamics shaped by stochastic variability and stabilizing forces. We develop a hybrid Ornstein–Uhlenbeck (OU)–Branching framework that couples mean-reverting stochastic trait dynamics with demographic birth–death processes to model lineage diversification under effective stabilizing constraints.

Methods: Using Escherichia coli long-term evolution experiment (LTEE) lineages (WT, priA, recG), we parameterized the equilibrium mean (μ), mean-reversion strength (θ), and diffusion scale (σ) on the log10 mutation-frequency axis via replicate-grouped likelihood inference. We performed forward simulations for predictive envelopes, uncertainty quantification, phase-plane dynamics, and OU–Branching lineage networks. We also ran illustrative in silico therapy simulations under fixed OU parameters with exposure-modulated birth/death rates.

Results: The fitted model recapitulated lineage-specific mutation dynamics and branching architectures. priA exhibited elevated stochastic dispersion and drift-prone behavior consistent with a high-plasticity regime, whereas recG showed constrained diversification and increased lineage turnover consistent with a collapse-prone regime. Illustrative therapy simulations generated oscillatory trait trajectories, suppression–rebound population dynamics, clonal pruning, and extinction-versus-persistence regimes.

Discussion: Although Y is directly observed as log10 mutation frequency in LTEE, in tumors Y can represent a longitudinally measurable phenotypic state (e.g., drug-tolerance scores from single-cell data, MRD/VAF-derived burden proxies, or pathway activity states). The balance between stabilizing strength (θ) and stochastic variability (σ) provides a quantitative axis governing plasticity and persistence, motivating future calibration to clinical longitudinal data for evolution-aware, patient-specific modeling.

Introduction

Pediatric cancers exhibit low mutation burdens and relatively simple karyotypes yet display aggressive behavior and complex evolutionary patterns (1–5). This paradox suggests that progression is driven less by continuous mutational accumulation and more by stochastic clonal selection and stabilizing forces operating within constrained phenotypic spaces (6–9). Understanding how variability, selection, and turnover interact to shape these dynamics remains a key challenge in precision oncology.

Traditional evolutionary models emphasize either deterministic selection or stochastic drift; however, pediatric cancers likely evolve through a hybrid regime in which stabilizing selection preserves essential traits while random fluctuations enable clonal diversification and adaptive escape (6, 7, 10, 11). Capturing this balance requires a framework that links phenotypic reversion toward equilibrium with probabilistic clonal branching.

In this study, a hybrid Ornstein–Uhlenbeck (OU) and branching-process model was developed to unify stochastic trait evolution with demographic birth–death dynamics (12, 13). The OU component describes mean-reverting phenotypic fluctuations around a lineage mean (μ) on the log₁₀ mutation-frequency scale, whereas the branching component captures random proliferation and extinction events that generate evolving lineage networks. Together, these coupled processes simulate how stabilizing selection and stochastic diversification jointly shape clonal evolution (Figure 1).

Figure 1

Figure 1. Linking OU traits to branching lineage dynamics. (A) Ornstein–Uhlenbeck (OU) dynamics illustrating mean reversion toward an equilibrium mean μ under stochastic diffusion σ and stabilizing strength θ. (B) Birth–death branching process showing stochastic clonal growth and extinction. (C) Schematic of E. coli LTEE lineages modeled under the OU–Branching framework. (D) Analogous schematic of pediatric tumor lineages showing diversification and extinction driven by OU-governed phenotypic evolution and branching proliferation.

As an experimental analog, the Escherichia coli long-term evolution experiment (LTEE) recapitulates general evolutionary phenomena relevant to tumors—mutation accumulation, clonal interference, and adaptive constraint—over measurable timescales (Supplementary Figures S1–S2) (14). A comparison of wild-type, priA (replication restart-deficient), and recG (replication fork repair-deficient) lineages quantifies lineage-specific mutation rates, phenotypic variance, and stabilizing selection strength (15–19). The priA mutant exhibits hypermutable, drift-prone behavior analogous to mutator subclones in relapsed pediatric cancers, whereas recG shows dynamics consistent with a collapse-prone regime under replication stress.

This framework links mutation-frequency dynamics, parameter inference, and lineage topology (Figures 1–6) to a unified view of evolutionary equilibria. Its extension to therapy simulations (Figure 7) provides an illustrative forward-simulation testbed for how cyclic exposure can reshape trait dynamics and clonal persistence. These simulations are hypothesis-generating and motivate future calibration on longitudinal clinical data.

Figure 2

Figure 2. Mutation-frequency trajectories and OU model fitting. Log-scale mutation-frequency trajectories in E. coli LTEE lineages (WT, priA, recG). Lines show replicate means ± SD. OU model fitting supports lineage-specific dynamics: priA exhibits an elevated equilibrium level (higher 10^μ) and greater stochastic dispersion, whereas WT and recG show more constrained trajectories. A likelihood-ratio test rejects a shared-parameter model (2ΔlogL = 33.18, df  = 6, p < 9.7 x 10^-6), supporting lineage-specific OU parameters (Table 2).

Overall, this study establishes a multiscale OU–Branching framework that links microbial evolutionary dynamics to tumor-relevant principles of constrained adaptation, stochastic diversification, and clonal turnover. By integrating stabilizing selection (θ) and stochastic variability (σ) with phenotype-coupled birth–death demography, the model yields interpretable regime behavior—stabilized, hypermutable/plastic, and collapse-prone—supported by likelihood-based inference, uncertainty quantification, and lineage-architecture simulations (Figures 1–6). Importantly, the same formalism is designed to transfer to pediatric tumors by defining $Y_t$ as any longitudinally measurable state variable, such as drug-tolerance or stress-response programs inferred from single-cell data, Minimal Residual Disease (MRD)–based burden proxies, or clonal fraction trajectories derived from time-series sequencing. In that setting, patient- or cohort-specific OU parameters (μ, θ, σ) can be estimated from repeated measurements using the same likelihood-based framework. Finally, the therapy simulations (Figure 7) are presented as an illustrative forward-simulation testbed showing how cyclic treatment can reshape trait dynamics and clonal persistence; this provides a practical roadmap for future calibration on clinical time-series data to support precision, evolution-aware pediatric oncology.

Methods

Estimation of lineage-specific parameters

Mutation-frequency data from E. coli LTEE populations (WT, priA, recG) (Supplementary Figures S1–S2; Supplementary Table S1) were preprocessed on a per-lineage basis to enable log-scale modeling. Let $x_{r,t}\ge 0$ denote the observed mutation frequency for replicate r at time t. In this microbial application, $Y_t$ is directly observed as $y_{r,t}$, the log₁₀ floored mutation frequency for replicate r at time t. The log₁₀ mutation frequency corresponds to that measured for replicate r. To avoid undefined log₁₀(0), zeros were replaced using a lineage-specific detection-floor constant

$ϵ=0.5 \times min\{x_{r,t} :x_{r,t}>0\}$,

and we defined the floored frequency

$x_{r,t}^{*}=max(x_{r,t},ϵ)$.

The analyzed trait was

$y_{r,t}={\text{log}}_{10}(x_{r,t}^{*})$.

Rationale for OU dynamics

We model the log₁₀ mutation-frequency trajectory $y_{r,t}$ using an OU process because the observed time series exhibit bounded stochastic fluctuations around lineage-specific levels rather than unbounded random-walk drift, consistent with mean-reverting dynamics under stabilizing constraints. The OU model also yields a closed-form transition density under irregular sampling intervals, enabling likelihood-based inference and reproducible uncertainty quantification from replicate trajectories.

Lineage-specific Ornstein–Uhlenbeck (OU) parameters μ (equilibrium mean), θ (mean-reversion rate), and σ (diffusion scale) were estimated by maximizing the exact OU transition likelihood, using a replicate-grouped objective (sum of negative log-likelihood across replicate trajectories). Optimization used bounded L-BFGS-B with constraints θ > 0 and σ > 0, with multiple random restarts to mitigate local optima. Lineage-specific fits were compared against a shared-parameter null using a likelihood-ratio test (LRT; df = 6). Parameter stability and identifiability were evaluated by curvature near the optimum and by profile-likelihood surfaces in (θ, σ) space.

Simulation of Ornstein–Uhlenbeck and branching dynamics

All simulations were based on the OU stochastic differential equation:

$$d{Y}_t=\theta (\mu -Y_t)dt+\sigma d{W}_t$$

where $Y_t$ is the phenotypic trait, μ the equilibrium mean, θ the mean-reversion rate, σ the diffusion scale (volatility; variance is σ²), and $W_t$ a Wiener process (20–23). OU trajectories were simulated using the exact discrete-time solution over eight arbitrary time units (illustrative schematic for Figure 1; not used for parameter inference) with parameters θ = 0.8, σ = 0.3, and dt = 0.01. In this framework, θ is interpreted as an effective stabilizing strength acting on the modeled trait $Y_t$. Depending on the biological instantiation, θ may reflect direct stabilizing selection on $Y_t$ or an emergent mean-reverting constraint arising from population-genetic and cellular limits (e.g., repair capacity, trade-offs, clonal interference, or ecological feedback). Distinguishing these mechanisms would require additional orthogonal measurements beyond the scope of the present study.

Branching processes followed a stochastic birth–death model with per-capita rates b = 0.05 and d = 0.04, initialized with 40 individuals per replicate. Hybrid OU–Branching lineage schematics were implemented in Python (NumPy + NetworkX + Matplotlib), where each clone inherited its parental trait plus OU noise, and survival probability scaled with trait-dependent fitness. Extant nodes were colored by phenotype (viridis), and extinct nodes were shown in gray and identified at the snapshot time (t = 20; Figures 1, 6, 7).

Model validation and forward simulations

OU model validation used the exact transition density and forward simulations (10,000 replicates per lineage) under fitted parameters. Predictive envelopes were generated from the 2.5th to 97.5th percentiles at each time point, representing 95% predictive intervals. Observed mutation frequencies were overlaid with fitted mean trajectories. Parameter stability was corroborated by the curvature of the log-likelihood surface near the optimum (Figure 3) and the profile-likelihood geometry in (θ, σ) (Figure 4).

Figure 3

Figure 3. OU model fits and 95% envelopes for mutation dynamics. OU model fits for E. coli lineages (WT, priA, recG) with 95% predictive envelopes for mutation dynamics. Solid lines show fitted means; shaded areas indicate predictive intervals from forward simulation under fitted parameters. priA exhibits a higher mean mutation frequency and greater variability, whereas WT and recG remain within narrower envelopes. The OU-derived envelopes capture lineage-specific evolutionary constraints.

Figure 4

Figure 4. Profile-likelihood surfaces of OU parameters (θ, σ). (A-C) Profile-likelihood heatmaps for E. coli lineages (WT, priA, recG) showing ΔNLL relative to the lineage-specific minimum as a function of (θ, σ) with μ fixed at its MLE. Black contours denote fixed ΔNLL levels. Points mark the MLEs for each lineage. priA displays a broader, asymmetric basin consistent with strong θ–σ coupling, whereas WT and recG show more sharply localized optima.

Bootstrap uncertainty quantification of OU parameters

To quantify estimation uncertainty for lineage-specific OU parameters, we performed a nonparametric bootstrap over biological replicates within each lineage (WT, priA, recG). For each bootstrap iteration (B = 2000), replicate trajectories were sampled with replacement, and OU parameters (μ, θ, σ) were refit by maximum likelihood using the same replicate-grouped transition likelihood. Percentile-based 95% intervals were computed from the bootstrap distributions for μ, θ, and σ. To characterize parameter coupling, we summarized the correlation between log(θ) and log(σ) across bootstrap refits (Table 1).

Table 1

Table 1. OU parameter estimates with replicate-grouped bootstrap 95% intervals.

Profile-likelihood surface analysis

Profile-likelihood surfaces of θ and σ were computed in R (v4.3.0). For each lineage, the maximum-likelihood estimate (MLE) of $\mu$ was fixed, and log-spaced grids of θ and σ were explored within ±3 log units around the optimum.

Exact OU transition and NLL

Let $y_i$ denote the observed trait at time $t_i$ and $\text{Δ}{t}_i=t_i-t_{i-1}$. Under the OU transition,

$y_i | y_{i-1} ~ \mathcal{N}(m_i,v_i),$where

$m_i=\mu +(y_{i-1}-\mu )e^{-\theta \text{Δ}{t}_i}$,

$v_i=\frac{{\sigma }^2}{2\theta }(1-e^{-2\theta \text{Δ}{t}_i})$.

The negative log-likelihood is

$\text{NLL}(\mu ,\theta ,\sigma )=\frac{1}{2} {\displaystyle \sum }_i[\text{log}(2\pi v_i)+\frac{{(y_i-m_i)}^2}{v_i}]$.

Because the OU transition density is closed-form under irregular sampling intervals, likelihood-based estimation (MLE/NLL minimization) provides a standard, reproducible inference approach for longitudinal trajectories.

Contours at fixed ΔNLL levels were visualized using ggplot2 and viridis (Figure 4). Fit summaries (NLL, AIC, and ΔAIC) are reported for reference using AIC = 2k + 2×NLL (k = 3 parameters) (Table 2).

Table 2

Table 2. Model comparison between a shared-parameter OU model and a lineage-specific OU model.

Formal quantitative trajectory separation metrics

To complement predictive envelopes, we quantified differences between lineage trajectories using model-implied predictive distributions. For each lineage, the OU predictive mean $m(t)$ and variance $v(t)$ were computed analytically from fitted parameters. Pairwise comparisons between lineages A and B were summarized across time points by:

i. mean absolute separation of predictive means $D_{mean}= mea{n}_t| m_A(t)- m_B(t)|$;

ii. the averaged 2-Wasserstein distance between the 1D Gaussian predictive distributions $\mathcal{N}(m_A(t), v_A(t))$ $\mathcal{N}(m_B(t), v_B(t))$, $W_2(\mathcal{N}(m_A,v_A), \mathcal{N}(m_B,v_B))= \sqrt{{(m_A-m_B)}^2+{(\sqrt{v_A}-\sqrt{v_B})}^2};$and

iii. a dominance probability $P(A>B)= mea{n}_t[\text{Φ}(\frac{m_A(t)-m_B(t)}{\sqrt{v_A(t)+v_B(t)}})]$. Uncertainty intervals for these trajectory metrics were obtained by recomputing the metrics across paired bootstrap draws from each lineage’s (μ, θ, σ) distribution (Table 3).

Table 3

Table 3. Pairwise trajectory separation metrics with bootstrap 95% intervals.

Phase-plane simulation of coupled dynamics

Phase-plane diagrams combined OU phenotypic evolution with phenotype-dependent birth–death dynamics. Deterministic vector fields visualized the skeleton dY/dt = θ(μ − Y) and dN/dt = (λ(Y) − δ(Y))N, where λ and δ depend on deviation from the lineage mean (Y − μ). Stochastic trajectories were simulated with exact discrete-time OU updates (Δt = 0.05) and Poisson birth/death events, and rendered with log-scaled N (Figure 5). For constant-rate simulations, we use b, d; for phenotype-coupled simulations, we use λ(Y), δ(Y) with b = λ₀ and d = δ₀ as baseline rates.

Figure 5

Figure 5. Phase-plane trajectories of the hybrid OU–Branching model. (A–C) Phase-plane plots coupling phenotypic evolution ($Y_t$ = log₁₀ mutation frequency) and population dynamics ($N_t$, log scale) for WT, priA, and recG, simulated using lineage-specific replicate-grouped OU parameters (Table 1). Gray arrows show the deterministic vector field defined by dY/dt = θ(μ − Y) and dN/dt = (λ(Y) − δ(Y))N. The gray dashed vertical line marks the lineage mean phenotype μ, and the orange dotted vertical line marks the initial phenotype $Y_0$; the demographic nullcline occurs where λ(Y) = δ(Y). Colored curves show representative stochastic OU–Branching trajectories, illustrating lineage-specific differences in constraint, stochastic dispersion, and demographic response. How to read this figure: first use the arrows to interpret the expected (deterministic) drift in (Y, N), then compare how stochastic trajectories fluctuate around μ and across the demographic nullcline to reveal differences in persistence and collapse tendencies across lineages.

Network construction of hybrid OU–Branching lineages

Lineage networks were generated using lineage-specific OU parameters (Table 1) and phenotype-coupled birth–death dynamics. At each time step, extant clones updated Y via the exact OU transition and produced offspring with division rate $\lambda (Y)= {\lambda }_0·exp(\alpha ·(Y-\mu ))$ and death rate $\delta (Y)= {\delta }_0·exp(-\beta ·(Y-\mu ))$, with rates clipped for numerical stability. Simulations were run for T = 40 time units, with extant and extinct clones identified at the snapshot (t = 20) for all lineages. Extinct clones were shown in gray, and extant clones were colored by phenotype using viridis; for visualization, the priA panel was downsampled to a maximum of 2000 extant nodes without altering the underlying simulation (Figure 6).

Figure 6

Figure 6. Hybrid OU–Branching lineage snapshots across E. coli lineages. (A–C) Lineage graphs for WT, priA, and recG were generated under phenotype-coupled birth–death dynamics using lineage-specific replicate-grouped OU parameters (Table 1) on the log₁₀ mutation-frequency scale. Extant clones are colored by phenotype $Y_t$ = log₁₀(mutation frequency) using a viridis scale, while extinct clones are shown in gray. Edge thickness is proportional to the instantaneous division rate λ(Y); darker edges denote extant branches, and lighter edges denote extinct branches. For visualization, the priA panel was downsampled to a maximum of 2000 extant nodes without altering the underlying simulation. How to read this figure: node color encodes phenotype, gray nodes indicate extinct lineages, and thicker/darker edges highlight branches with higher division rates and extant descendants.

Translational application of the hybrid OU–Branching framework

An OU–Branching simulator was developed in Python 3.11 to model therapy-modulated phenotypic and clonal dynamics under the priA profile. OU parameters were taken from replicate-grouped empirical fits on the log₁₀ mutation-frequency scale (Table 1; μ = −5.000, θ = 0.117, σ = 0.436) and coupled to a pharmacokinetic exposure function C(t) representing six 21-day treatment cycles.

Phenotypic evolution followed an OU stochastic differential equation, whereas population growth and extinction were modeled using a phenotype- and drug-dependent branching process (24). Division and death rates were defined as $\lambda (Y_t, C)= {\lambda }_0·exp(\alpha ·(Y_t-{\mu }_0)-k\lambda ·g(C))$ and $\delta (Y_t, C)= {\delta }_0·exp(-\beta ·(Y_t-{\mu }_0)+k\delta ·g(C))$, where $g(C)=C/(C+EC50)$ describes drug saturation (25). Here, μ₀ denotes the baseline (pre-therapy) optimum; in this illustrative implementation, we set μ₀ = μ from the replicate-grouped OU fit. In this formulation, increasing C suppresses division and increases death, capturing therapy-induced growth inhibition on a scale-invariant phenotype axis.

Each regimen comprised 256 stochastic replicates simulated for 200 days with Δt = 0.5 day, encompassing six dosing cycles and post-treatment periods. The ensemble distributions of $Y_t$ and $N_t$ describe therapy-driven phenotypic convergence, oscillatory tumor suppression, and extinction probabilities in the priA-mutant background, demonstrating the translational potential of the hybrid OU–Branching model for pediatric precision oncology (Figure 7).

Results

Mutation frequency trajectories and Ornstein–Uhlenbeck model fitting across bacterial lineages

Lineage-specific mutational dynamics were quantified in three E. coli strains—wild-type (WT), priA, and recG—using a hybrid OU–Branching framework. Mutation-frequency trajectories (Figure 2) showed marked divergence among lineages. The priA mutants accumulated mutations at a higher rate, consistent with replication restart defects leading to hypermutability. In contrast, recG displayed reduced mutation frequency; despite large diffusion (σ), its strong mean reversion (θ) constrained long-run dispersion, consistent with a collapse-prone regime.

OU parameter estimates (Figure 2; Table 1) revealed significant heterogeneity across lineages. The equilibrium mean (μ) reflects the steady-state mutation frequency; on the original scale (10^μ), priA exceeded WT by ~60-fold (~1.8 orders of magnitude). The diffusion scale (σ) and stationary variance (${\sigma }^2$) also differentiated lineages, indicating more substantial random fluctuations in priA. The likelihood ratio test (2ΔlogL = 33.18, df = 6, p = 9.7 x 10^-6) confirmed that lineage-specific OU parameters provided a significantly better fit than a shared-parameter null model (Table 2). Collectively, these results demonstrate that bacterial replication–repair mutants evolve under distinct stabilizing selection regimes that govern mutation accumulation.

OU model fits and 95% envelopes for mutation-frequency dynamics in bacterial lineages

Model-based predictive accuracy was evaluated by comparing observed mutation-frequency data with OU model simulations using lineage-specific parameter estimates. The fitted trajectories (Figure 3) recapitulated empirical trends, with the priA mutant exhibiting higher mean mutation frequencies and broader variance envelopes than either WT or recG. The expanded 95% predictive interval for priA indicates increased stochastic drift and reduced phenotypic stability resulting from defective replication restart pathways.

In contrast, WT and recG lineages displayed tightly constrained envelopes centered around low equilibrium values, reflecting stronger stabilizing selection and minimal stochastic deviation. Across replicate simulations, 95% of observed data points fell within the OU-predicted intervals, validating the robustness of the lineage-specific maximum-likelihood parameterization. Collectively, these model fits reveal distinct selective regimes influencing mutation-frequency evolution across the three bacterial genetic backgrounds.

Uncertainty quantification of OU parameter estimates

To quantify uncertainty beyond point estimates and profile-likelihood geometry, we performed replicate-bootstrap refitting of the OU model (B = 2000). Across lineages, μ, θ, and σ showed lineage-specific uncertainty profiles and evident θ–σ coupling (Table 1). In particular, priA exhibited a low mean-reversion rate (θ ≈ 0.117) with strong θ–σ coupling (corr[logθ, logσ] ≈ 0.993), yielding the largest stationary variance (${\sigma }^2$/(2θ) median ≈ 0.817, 95% interval ≈ [0.753, 0.909]) on the log₁₀ scale.

Quantitative trajectory separation and sensitivity to parameter uncertainty

We complemented predictive envelopes with formal pairwise trajectory-separation metrics based on OU predictive distributions (Table 3), propagating uncertainty using replicate-grouped bootstrap parameter draws from Table 1. Across time points, priA showed strong separation from recG (D_mean ≈ 1.944 log10 units, 95% interval [1.912, 1.987], corresponding to ≈ 87.8-fold separation) and from WT (D_mean ≈ 1.045 log₁₀ units, 95% interval [1.018, 1.091], corresponding to ≈ 11.1-fold separation). The recG–WT comparison showed smaller divergence (D_mean ≈ 0.898 log₁₀ units, 95% interval [0.881, 0.917]). Bootstrap intervals were narrow relative to effect sizes, indicating robust qualitative regime separation.

Profile-likelihood analysis of OU parameters (θ, σ) for E. coli lineages

Two-dimensional profile-likelihood surfaces were computed for the mean-reversion rate (θ) and diffusion scale (σ) across E. coli lineages (Figures 4A–C). The WT and recG surfaces exhibited narrow elliptical contours centered on distinct optima, indicating well-constrained parameter estimates. In contrast, the priA mutant displayed a broad and skewed likelihood surface, reflecting strong coupling between θ and σ and greater stochastic uncertainty. The shallow curvature of the priA ΔNLL landscape suggests that multiple parameter combinations produce comparable likelihoods, consistent with higher intrinsic noise and reduced stabilizing selection in its evolutionary dynamics.

Fit summaries (NLL and AIC) are reported for reference (Table 2). Consistent with the profile-likelihood geometry, priA shows strong θ–σ coupling and broader uncertainty, whereas WT and recG exhibit more sharply constrained optima (Figure 4).

Coupled phenotypic and population-dynamics phase-planes

Phase-plane diagrams derived from stochastic simulations of the hybrid OU–Branching model captured the dynamic interplay between phenotypic evolution ($Y_t$) and population size ($N_t$, log scale). Each E. coli lineage exhibited a distinct coupling pattern between OU-driven phenotypic trajectories and branching-process demography.

The WT lineage exhibited smooth, monotonic trajectories, characterized by a gradual, well-regulated reduction in population size, accompanied by phenotypic reversion toward the lineage mean (μ). This behavior reflects strong stabilizing selection and minimal diffusion scale (σ), consistent with robust replication fidelity and steady population maintenance (Figure 5A). In contrast, the priA mutant showed large-amplitude oscillations in $Y_t$ and broad fluctuations in $N_t$ spanning several orders of magnitude, indicating mutator-driven stochastic drift and weakened phenotypic constraint (Figure 5B). The recG lineage remained confined near the equilibrium boundary, exhibiting modest noise and limited recovery following population collapse, which highlights impaired adaptive flexibility under replication-fork repair deficiency (Figure 5C).

The overlaid vector fields visualize the joint OU–Branching dynamics, with horizontal components representing phenotypic reversion toward μ and vertical components indicating exponential growth–decay cycles driven by the branching process. The blue trajectories trace reversion phases associated with population decline, most pronounced in WT, where stabilizing selection suppresses stochastic noise. Together, these results demonstrate that the hybrid OU–Branching framework accurately captures lineage-specific evolutionary equilibria, transient instability, and the contrasting balance between stabilizing selection and stochastic drift in DNA-repair-deficient lineages.

Hybrid OU–Branching lineage architectures across E. coli lineages

Hybrid OU–Branching simulations generated lineage networks that visualized the balance between stabilizing selection and stochastic diversification across DNA-repair backgrounds.

In the WT lineage, network connectivity remained compact, with few extant subclones and narrow phenotypic dispersion around the equilibrium mean, consistent with strong stabilizing selection and low diffusion-driven variability (Figure 6A). The priA lineage exhibited extensive branching, higher node density, and broader color gradients, reflecting increased phenotypic drift and elevated mutation rate (Figure 6B). In contrast, the recG lineage produced sparse and short-lived branches, consistent with impaired replication restart and reduced adaptive potential (Figure 6C). Despite high θ, recG lineage underwent rapid lineage extinction and minimal diversification, highlighting the interplay between deterministic constraint and stochastic collapse. Collectively, these results demonstrate that OU–Branching coupling produces lineage architectures that distinguish stable (WT), hypermutable (priA), and collapse-prone (recG) evolutionary modes.

Hybrid OU–Branching dynamics under therapy

To illustrate how the OU–Branching formalism can serve as a forward-simulation testbed for therapy-modulated evolution, we performed in silico simulations using the priA profile parameters inferred from replicate-grouped OU fits (Table 1; μ = −5.000, θ = 0.117, σ = 0.436) (Figure 7). These simulations are illustrative and hypothesis-generating: OU parameters are held fixed, drug exposure enters phenomenologically through the birth–death rates, and the results are not calibrated to patient response data or intended as prescriptive scheduling recommendations.

Panel A shows phenotypic trajectories ($Y_t$, blue) and drug exposure C(t) (red) across six 21-day treatment cycles. Under this parameterization, phenotypes oscillate with dosing intervals and, on average, relax toward the long-run optimum μ, illustrating how cyclic exposure can interact with mean reversion (θ) and diffusion (σ) to shape temporal trait variability.

Panel B displays simulated population dynamics ($N_t$) on a logarithmic scale. Treated simulations (red) exhibit oscillatory suppression and rebound before approaching a lower plateau relative to untreated controls (gray dashed), reflecting repeated perturbations of the division and death rates $\lambda (Y_t, C)$ and $\delta (Y_t, C)$ that modulate net growth.

Panel C shows an example OU–Branching lineage architecture at t = 20 days, with clones colored by phenotype ($Y_t$) and edges scaled by the instantaneous division rate $\lambda (Y_t, C)$. The topology illustrates early diversification followed by pruning of peripheral branches under phenotype-coupled demography.

Panel D presents an illustrative cure-probability map P(cure) as a function of diffusion (σ) and selection strength (θ). Higher θ and lower σ favor extinction, whereas weaker mean reversion or higher stochastic diffusion promotes persistence/relapse.

Together, these simulations demonstrate—in an illustrative, non-predictive setting—how coupling OU-based phenotypic drift with branching process dynamics can generate therapy-associated stabilization, clonal pruning, and extinction-versus-persistence regimes that motivate future calibration on clinical longitudinal measurements.

Discussion

The hybrid Ornstein–Uhlenbeck (OU)–Branching framework offers a quantitative way to reason about therapy-associated dynamics and persistence by jointly modeling (i) therapy-shaped fluctuations of a latent or measured phenotypic state $Y_t$ and (ii) stochastic clonal survival through birth–death demography. Using priA-profile parameters inferred from replicate-grouped OU fits (μ = −5.000, θ = 0.117, σ = 0.436; Table 1) as a calibrated high-plasticity reference, we implemented an illustrative in silico cyclic-therapy regimen in which drug exposure modulates division and death rates while the trait dynamics remain mean-reverting (Figures 7A, B). In this setting, repeated dosing narrows phenotypic dispersion and induces oscillatory suppression–rebound dynamics in tumor burden as peripheral subclones are pruned during treatment and re-expand between cycles. Critically, $Y_t$ should be interpreted as an abstract phenotype axis in these simulations; in real tumors, it can be instantiated using longitudinal observables such as drug tolerance or differentiation-state scores from serial single-cell profiling, MRD trajectories, or clonal fraction proxies from time-series sequencing. This provides a direct path to patient-specific use: estimate (μ, θ, σ) from repeated measurements, couple them to therapy-dependent birth–death terms (and, where appropriate, therapy-dependent shifts in μ/θ/σ), and then explore how small changes in stabilizing strength or stochasticity alter extinction versus persistence regimes.

Linking microbial and tumor evolutionary dynamics

The E. coli LTEE provides a tractable analog for studying general principles of tumor evolution, enabling direct observation of clonal diversification, stabilizing selection, and extinction events over thousands of generations (14, 26–28). Analyses of priA and recG mutants revealed lineage-specific differences in mutation frequency, phenotypic stability, and adaptive potential. The priA lineage displayed high stochastic variance and weak stabilizing selection, consistent with a mutator phenotype that parallels hypermutable pediatric tumor subclones. Conversely, recG exhibited constrained diversification and frequent extinction, reminiscent of repair-deficient cancer lineages with limited adaptive flexibility. These distinct evolutionary regimes—stabilized, hypermutable, and collapse-prone—represent fundamental modes of tumor evolutionary behavior within a shared mathematical landscape.

Integrating stochasticity and selection through OU–Branching coupling

Classical models of cancer evolution often separate selection-driven deterministic growth from stochastic noise resulting from mutation or genetic drift (29–31). The hybrid OU–Branching model unifies these forces, allowing stochastic fluctuations to coexist with restoring selection toward a phenotypic optimum. The OU process captures trait-level stability and noise, whereas the branching component introduces demographic uncertainty via random birth–death processes (32–35). This integration explains how even low-mutation pediatric tumors can exhibit high intratumoral heterogeneity: stochastic phenotypic variation is continually generated and pruned by stabilizing selection, thereby maintaining evolutionary plasticity without requiring substantial mutational input.

Quantitative insights into lineage evolution

Parameter inference and simulation analyses provide quantitative evidence that stochastic variance (σ²) and mean-reversion strength (θ) determine lineage diversification and persistence. Lineages characterized by high σ² and low θ, such as priA, exhibit broad phenotypic dispersion and long-tailed lineage survival distributions. In contrast, lineages with strong mean reversion (high θ)—including WT and recG—exhibit constrained long-run dispersion (σ²/(2θ)) and reduced persistence of extreme phenotypic excursions. The model formalizes an evolutionary trade-off between plasticity and stability: excessive stochasticity promotes transient adaptation but risks instability, whereas strong reversion constrains exploration while preserving lineage integrity. This trade-off may underlie the tendency of some pediatric cancers to relapse following therapy, as treatment pressures perturb θ and σ², altering the effective “evolutionary temperature” of the tumor ecosystem. These regime distinctions are supported by trajectory-separation metrics (Table 3), with priA showing the strongest separation from recG and WT.

Translational implications for therapy modeling

The hybrid OU–Branching framework provides a unifying stochastic language for therapy-modulated evolution by coupling mean-reverting phenotypic dynamics with demographic birth–death branching. Here, we use therapy simulations as an illustrative forward-simulation testbed to show how the model behaves under cyclic exposure when OU parameters are held fixed, and treatment effects are introduced phenomenologically through exposure-modulated birth and death rates. Accordingly, these simulations are hypothesis-generating and are not calibrated to clinical response data, nor do they imply optimal scheduling or clinical efficacy.

Using the priA profile parameters inferred from replicate-grouped OU fits (μ = −5.000, θ = 0.117, σ = 0.436; Table 1), the model generates cyclic oscillations of phenotype ($Y_t$) and population size ($N_t$) across six 21-day exposure cycles (Figures 7A, B). Under this parameterization, phenotypic trajectories fluctuate with dosing intervals and, on average, relax toward the long-run optimum μ, illustrating how the balance between mean reversion (θ) and stochastic diffusion (σ) can shape therapy-associated variability over time. In parallel, population trajectories exhibit a suppression–rebound pattern under repeated exposure, consistent with alternating periods of demographic contraction and regrowth in the branching process.

When the same fixed OU parameters are used to generate lineage networks at t = 20 days (Figure 7C), the resulting architecture shows early diversification followed by extinction of peripheral branches, reflecting stochastic pruning that can arise when demographic rates depend on the evolving phenotypic state.

Figure 7

Figure 7. Translational application of the hybrid OU–Branching framework (priA profile). (A) Phenotypic evolution under cyclic exposure showing the median trajectory of $Y_t$ = log₁₀(mutation frequency) (blue) alongside drug exposure C(t) (red) across six 21-day cycles. (B) Simulated population dynamics on a log scale comparing treated (red) versus untreated control (gray dashed) median trajectories. (C) Hybrid OU–Branching lineage architecture at the shared snapshot time t  = 20, with clones colored by phenotype and edges scaled by the instantaneous division rate λ($Y_t$, C(t)). (D) Illustrative outcome map showing how P(cure) varies with phenotypic diffusion σ (log₁₀ scale) and mean reversion θ. Panels (A–C) use the priA profile parameters inferred from replicate-grouped OU fits (Table 1). How to read this figure: exposure C(t) specifies the imposed schedule, while trajectories and lineage graphs are forward simulations under fixed OU parameters; the panels are intended to be hypothesis-generating and not predictive of clinical efficacy or optimal scheduling.

The outcome surface in Figure 7D is presented as an illustrative regime diagram: within the model, higher θ and lower σ shift simulations toward extinction, whereas weaker mean reversion or larger stochastic diffusion favors persistence. We emphasize that P(cure) here denotes a model-based extinction probability under the assumed dynamics, rather than a clinically validated cure probability.

In summary, the therapy simulations demonstrate—in a deliberately simplified, non-predictive setting—how coupling OU-based phenotypic drift with branching-process demography can produce oscillatory trajectories, clonal pruning, and extinction-versus-persistence regimes. A key next step is to calibrate $Y_t$ to longitudinal tumor measurements (e.g., MRD/VAF-derived burden proxies, inferred drug-tolerance or stress-response programs from serial single-cell profiling, or pathway activity scores) and to estimate patient- or cohort-specific parameters. Such calibration would enable principled sensitivity analyses and comparative testing of evolution-aware hypotheses across therapeutic contexts, without implying prescriptive optimization or clinical prediction in the present study.

Model limitations and future directions

Although the hybrid OU–Branching framework captures foundational evolutionary principles, several important limitations should be acknowledged. First, the current formulation reduces phenotypic evolution to a single latent trait evolving under effective stabilizing forces, whereas real tumors evolve on multidimensional, context-dependent fitness landscapes shaped by metabolism, epigenetic plasticity, microenvironmental heterogeneity, and therapy-induced state transitions. Extending the model to multi-trait (vector) OU dynamics, potentially with hierarchical coupling among traits and patient-specific covariates, would better reflect the polygenic and multi-axis nature of tumor adaptation.

Second, the present branching process is an effectively well-mixed, lineage-level approximation that does not explicitly encode spatial structure, local resource limitation, or niche partitioning. Spatial constraints, dispersal, and local competition can substantially alter clonal dynamics in tumors and generate distinct evolutionary modes (e.g., branching evolution, neutral expansion, or clonal sweeps) (36, 37). The well-mixed approximation is most reasonable when the modeled state variable and the data are lineage-aggregated (e.g., clonal fraction trajectories) and/or in regimes where mixing is enhanced, and population sizes are relatively small—such as minimal residual disease (MRD) or post-treatment residual disease in hematologic malignancies—but it is less appropriate for strongly structured solid tumors. Incorporating density dependence (carrying capacity), spatial embedding, and explicit niche structure would broaden biological realism and improve applicability across tumor types.

Relatedly, recent advances in spatial genomics and clone-mapping techniques (e.g., base-specific in situ sequencing) enable reconstruction of spatial clonal topologies and mutational landscapes in human tumors (38). These data provide a direct opportunity to align model behavior with spatially resolved lineage structure. Experimental validation remains critical: in pediatric leukemia or solid tumors, serial measurements such as time-series variant allele fractions, MRD-linked burden proxies, and (where feasible) single-cell lineage tracing could be used to estimate OU–Branching parameters (θ, σ²) under an explicit observation-noise model. Perturbation experiments targeting DNA repair, replication stress, or chromatin regulation could then test whether inferred parameter shifts are associated with increased lineage extinction versus persistence.

In future work, the model could be expanded to incorporate multiple traits, spatial structure, microenvironmental niches, and cell–cell interactions. Embedding OU–Branching dynamics within spatial frameworks such as the Cellular Potts Model or lattice-based spatial automata (e.g., BIO-LGCA) may help explore how migration, adhesion, and heterogeneity shape evolutionary trajectories (39–42). Finally, close collaboration with longitudinal oncology studies that collect repeated measurements across therapy will be essential to calibrate and validate the framework in real tumor systems.

For pediatric tumors, the OU–Branching state variable $Y_t$ can be instantiated as a longitudinally measurable cell-state or burden proxy, rather than mutation frequency per se. Practical choices include: (i) drug-tolerance or stress-response scores and differentiation-state programs derived from serial single-cell RNA-seq; (ii) MRD-linked burden metrics or log-transformed tumor fraction trajectories; and (iii) clonal fraction proxies from time-series sequencing (e.g., VAF- or clone-abundance–based summaries). With $Y_t$ defined and repeatedly measured, patient- or cohort-specific (μ, θ, σ) can be estimated using an observation model that accounts for measurement noise and sampling variability, enabling calibrated forward simulations that connect stabilization versus plasticity to lineage persistence and treatment sensitivity.

Conclusions

The hybrid OU–Branching framework provides a quantitative bridge between evolutionary biology and translational oncology by linking constrained (mean-reverting) phenotypic dynamics with stochastic lineage diversification and extinction. By uniting microbial and cancer evolution within a shared stochastic–stabilizing paradigm, it offers a principled way to interpret how pediatric tumors may retain adaptive flexibility despite low mutation burdens, while remaining subject to effective stabilizing forces.

More broadly, the key contribution of this work is a generalizable inferential formalism that couples continuous trait evolution and discrete branching within a single likelihood-based framework. This unified perspective enables coherent estimation of latent trajectories and lineage architecture from longitudinal data, and it provides a transparent foundation for future calibration to clinical time-series measurements to support evolution-aware hypothesis testing in pediatric oncology.

Data availability statement

All data analyzed in this study are included in the article and its Supplementary Material. Code used for analyses and figure generation is publicly available in a GitHub repository: https://github.com/shkim9391/hybrid-OU-branching-pediatric-evolution.

Ethics statement

This study did not involve human participants, human data, or vertebrate animals. The work analyzed microbial (E. coli) mutation-frequency time series and performed in silico simulations; therefore, ethical approval was not required.

Author contributions

SK: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declared that financial support was not received for this work and/or its publication.

Acknowledgments

The author thanks Fisher College for institutional support, including research time and resources that made this work possible.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2026.1727973/full#supplementary-material

References

1. Thatikonda V, Islam SMA, Autry RJ, Jones BC, Gröbner SN, Warsow G, et al. Comprehensive analysis of mutational signatures reveals distinct patterns and molecular processes across 27 pediatric cancers. Nat Cancer. (2023) 4:276–89. doi: 10.1038/s43018-022-00509-4

PubMed Abstract | Crossref Full Text | Google Scholar

2. Gröbner SN, Worst BC, Weischenfeldt J, Buchhalter I, Kleinheinz K, Rudneva VA, et al. The landscape of genomic alterations across childhood cancers. Nature. (2018) 555:321–7. doi: 10.1038/nature25480

PubMed Abstract | Crossref Full Text | Google Scholar

3. Brady SW, Roberts KG, Gu Z, Shi L, Pounds S, Pei D, et al. The genomic landscape of pediatric acute lymphoblastic leukemia. Nat Genet. (2022) 54:1376–89. doi: 10.1038/s41588-022-01159-z

PubMed Abstract | Crossref Full Text | Google Scholar

4. Sweet-Cordero EA and Biegel JA. The genomic landscape of pediatric cancers: implications for diagnosis and treatment. Science. (2019) 363:1170–5. doi: 10.1126/science.aaw3535

PubMed Abstract | Crossref Full Text | Google Scholar

5. Greenhalgh R, Yang W, Brady SW, Flasch DA, Liu Y, Szlachta KA, et al. The landscape of structural variation in pediatric cancer. bioRxiv. (2025) 2025.04.24.650284. doi: 10.1101/2025.04.24.650284

PubMed Abstract | Crossref Full Text | Google Scholar

6. Rozhok AI and DeGregori J. Toward an evolutionary model of cancer: causes, consequences, and therapeutic implications. Proc Natl Acad Sci U.S.A. (2015) 112:8914–21. doi: 10.1073/pnas.1501713112

PubMed Abstract | Crossref Full Text | Google Scholar

7. Lipinski KA, Barber LJ, Davies MN, Ashenden M, Sottoriva A, and Gerlinger M. Cancer evolution and the limits of predictability in precision cancer medicine. Trends Cancer. (2016) 2:49–63. doi: 10.1016/j.trecan.2015.11.003

PubMed Abstract | Crossref Full Text | Google Scholar

8. Huang S and Soto AM. The end of the genetic paradigm of cancer. PloS Biol. (2025) 23:e3003052. doi: 10.1371/journal.pbio.3003052

PubMed Abstract | Crossref Full Text | Google Scholar

9. Whiting FJH, Househam J, Baker A-M, Sottoriva A, and Graham TA. Phenotypic noise and plasticity in cancer evolution. Trends Cell Biol. (2024) 34:451–64. doi: 10.1016/j.tcb.2023.10.002

PubMed Abstract | Crossref Full Text | Google Scholar

10. Davies PC and Agus DB. Stochasticity and determinism in cancer creation. Converg Sci Phys Oncol. (2016) 1:26003. doi: 10.1088/2057-1739/1/2/026003

PubMed Abstract | Crossref Full Text | Google Scholar

11. Ben-David U, Rameen B, and Golub TR. Genomic evolution of cancer models: perils and opportunities. Nat Rev Cancer. (2019) 19:97–109. doi: 10.1038/s41568-018-0095-3

PubMed Abstract | Crossref Full Text | Google Scholar

12. Bartoszek K, Pienaar J, Mostad P, Andersson S, and Hansen TF. A phylogenetic comparative method for studying multivariate adaptation. J Theor Biol. (2012) 314:204–15. doi: 10.1016/j.jtbi.2012.08.005

PubMed Abstract | Crossref Full Text | Google Scholar

13. Bartoszek K. Using the Ornstein–Uhlenbeck process to model the evolution of interacting populations. J Theor Biol. (2017) 429:35–45. doi: 10.1016/j.jtbi.2017.06.011

PubMed Abstract | Crossref Full Text | Google Scholar

14. Lenski RE. Experimental evolution and the dynamics of adaptation and genome evolution in microbial populations. ISME J. (2017) 11:2181–94. doi: 10.1038/ismej.2017.69

PubMed Abstract | Crossref Full Text | Google Scholar

15. Kim S-H, Pytlos MJ, and Sinden RR. Replication restart: a pathway for (CTG)·(CAG) repeat deletion in Escherichia coli. Mutat Res. (2006) 595:5–22. doi: 10.1016/j.mrfmmm.2005.07.010

PubMed Abstract | Crossref Full Text | Google Scholar

16. Sandler SJ and Marians KJ. Role of PriA in replication fork reactivation in Escherichia coli. J Bacteriol. (2000) 182:9–13. doi: 10.1128/jb.182.1.9-13.2000

PubMed Abstract | Crossref Full Text | Google Scholar

17. Windgassen TA, Leroux M, Satyshur KA, Sandler SJ, and Keck JL. Structure-specific DNA replication-fork recognition directs helicase and replication restart activities of the PriA helicase. Proc Natl Acad Sci U.S.A. (2018) 115:E9075–84. doi: 10.1073/pnas.1809842115

PubMed Abstract | Crossref Full Text | Google Scholar

18. Romero ZJ, Chen SH, Armstrong T, Wood EA, van Oijen A, Robinson A, et al. Resolving toxic DNA repair intermediates in every Escherichia coli cell. Nucleic Acids Res. (2020) 48:8445–60. doi: 10.1093/nar/gkaa579

PubMed Abstract | Crossref Full Text | Google Scholar

19. McKenzie AM, Henry C, Myers KS, Place MM, and Keck JL. Identification of genetic interactions with priB links the PriA/PriB DNA replication restart pathway to double-strand DNA break repair in Escherichia coli. G3 (Bethesda). (2022) 12:jkac295. doi: 10.1093/g3journal/jkac295

PubMed Abstract | Crossref Full Text | Google Scholar

20. Hansen TF. Stabilizing selection and the comparative analysis of adaptation. Evolution. (1997) 51:1341–51. doi: 10.1111/j.1558-5646.1997.tb01457.x

PubMed Abstract | Crossref Full Text | Google Scholar

21. Butler MA and King AA. Phylogenetic comparative analysis: a modeling approach for adaptive evolution. Am Nat. (2004) 164:683–95. doi: 10.1086/426002

PubMed Abstract | Crossref Full Text | Google Scholar

22. Blomberg SP, Rathnayake SI, and Moreau CM. Beyond Brownian motion and the Ornstein–Uhlenbeck models of trait evolution. Am Nat. (2020) 195:145–65. doi: 10.1086/706339

PubMed Abstract | Crossref Full Text | Google Scholar

23. Ho LST and Ane C. Intrinsic inference difficulties for trait evolution with Ornstein–Uhlenbeck models. Methods Ecol Evol. (2014) 5:1133–46. doi: 10.1111/2041-210X.12285

Crossref Full Text | Google Scholar

24. Waclaw B, Bozic I, Pittman ME, Hruban RH, Vogelstein B, and Nowak MA. A spatial model predicts that dispersal and cell turnover limit intratumour heterogeneity. Nature. (2015) 525:261–4. doi: 10.1038/nature14971

PubMed Abstract | Crossref Full Text | Google Scholar

25. Simeoni M, Magni P, Cammia C, De Nicolao G, Croci V, Pesenti E, et al. Predictive pharmacokinetic–pharmacodynamic modeling of tumor growth kinetics in xenograft models after anticancer treatment. Cancer Res. (2004) 64:1094–101. doi: 10.1158/0008-5472.can-03-2524

PubMed Abstract | Crossref Full Text | Google Scholar

26. Barrick JE and Lenski RE. Genome dynamics during experimental evolution. Nat Rev Genet. (2014) 14:827–39. doi: 10.1038/nrg3564

PubMed Abstract | Crossref Full Text | Google Scholar

27. Maddamsetti R and Grant NA. Divergent evolution of mutation rates and biases in the long-term evolution experiment with Escherichia coli. Genome Biol Evol. (2020) 12:1591–605. doi: 10.1093/gbe/evaa178

PubMed Abstract | Crossref Full Text | Google Scholar

28. León D and González A. Modeling evolution in a long-term evolution experiment with Escherichia coli. arXiv. (2018) arXiv:1804.02660. doi: 10.48550/arXiv.1804.02660

Crossref Full Text | Google Scholar

29. Beerenwinkel N, Schwarz RF, Gerstung M, and Markowetz F. Cancer evolution: mathematical models and computational inference. Syst Biol. (2015) 64:e1–25. doi: 10.1093/sysbio/syu081

PubMed Abstract | Crossref Full Text | Google Scholar

30. Yin A, Moes DJAR, van Hasselt JGC, Swen JJ, and Guchelaar H-J. A review of mathematical models for tumor dynamics and treatment resistance evolution of solid tumors. CPT Pharmacometrics Syst Pharmacol. (2019) 8:720–37. doi: 10.1002/psp4.12450

PubMed Abstract | Crossref Full Text | Google Scholar

31. Novozhilov AS, Karev GP, and Koonin EV. Biological applications of the theory of birth-and-death processes. Brief Bioinform. (2006) 7:70–85. doi: 10.1093/bib/bbl006

PubMed Abstract | Crossref Full Text | Google Scholar

32. Cooper N, Thomas GH, Venditti C, Meade A, and Freckleton RP. A cautionary note on the use of Ornstein–Uhlenbeck models in macroevolutionary studies. Biol J Linn Soc Lond. (2015) 118:64–77. doi: 10.1111/bij.12701

PubMed Abstract | Crossref Full Text | Google Scholar

33. Noble R, Burri D, Le Sueur C, Lemant J, Viossat Y, Kather JN, et al. Spatial structure governs the mode of tumour evolution. Nat Ecol Evol. (2022) 6:207–17. doi: 10.1038/s41559-021-01615-9

PubMed Abstract | Crossref Full Text | Google Scholar

34. Carmona-Fontaine C, Bucci V, Akkari L, Deforet M, Joyce JA, and Xavier JB. Emergence of spatial structure in the tumor microenvironment due to the Warburg effect. Proc Natl Acad Sci U.S.A. (2013) 110:19402–7. doi: 10.1073/pnas.1311939110

PubMed Abstract | Crossref Full Text | Google Scholar

35. Chkhaidze K, Heide T, Werner B, Williams MJ, Huang W, Caravagna G, et al. Spatially constrained tumour growth affects the patterns of clonal selection and neutral drift in cancer genomic data. PloS Comput Biol. (2019) 15:e1007243. doi: 10.1371/journal.pcbi.1007243

PubMed Abstract | Crossref Full Text | Google Scholar

36. Fu X, Zhao Y, Lopez JI, Rowan A, Au L, Fendler A, et al. Spatial patterns of tumour growth impact clonal diversification in clear cell renal cell carcinoma. Nat Ecol Evol. (2022) 6:88–102. doi: 10.1038/s41559-021-01586-x

PubMed Abstract | Crossref Full Text | Google Scholar

37. Lomakin A, Svedlund J, Strell C, Gataric M, Shmatko A, Rukhovich G, et al. Spatial genomics maps the structure, nature and evolution of cancer clones. Nature. (2022) 611:594–602. doi: 10.1038/s41586-022-05425-2

PubMed Abstract | Crossref Full Text | Google Scholar

38. Deutsch A, Nava-Sedeño JM, Syga S, and Hatzikirou H. BIO-LGCA: a cellular automaton modelling class for analysing collective cell migration. PloS Comput Biol. (2021) 17:e1009066. doi: 10.1371/journal.pcbi.1009066

PubMed Abstract | Crossref Full Text | Google Scholar

39. Alber MS, Kiskowski MA, Glazier JA, and Jiang Y. On Cellular Automaton Approaches to Modeling Biological Cells. In: Rosenthal J and Gilliam DS, editors. Mathematical Systems Theory in Biology, Communications, Computation, and Finance. The IMA Volumes in Mathematics and its Applications, vol. 134 . Springer, New York, NY (2003). doi: 10.1007/978-0-387-21696-6_1

Crossref Full Text | Google Scholar

40. Ghaemi M and Shahrokhi A. Combination of the cellular Potts model and lattice gas cellular automata for simulating avascular cancer growth. arXiv. (2006) 4173:297–303. doi: 10.1007/11861201_35

Crossref Full Text | Google Scholar

41. Van Liedekerke P, Palm MM, Jagiella N, and Drasdo D. Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results. Comput Part Mech. (2015) 2:401–44. doi: 10.1007/s40571-015-0082-3

Crossref Full Text | Google Scholar

42. Metzcar J, Wang Y, Heiland R, and Macklin P. A review of cell-based computational modeling in cancer biology. JCO Clin Cancer Inform. (2019) 3:1–13. doi: 10.1200/CCI.18.00069

PubMed Abstract | Crossref Full Text | Google Scholar